These Doctora

Méthodes itératives pour la résolution d’équations
matricielles
El Mostafa Sadek
To cite this version:

El Mostafa Sadek. Méthodes itératives pour la résolution d’équations matricielles. Modélisation et
simulation. Université du Littoral Côte d’Opale; Université Cadi Ayyad (Marrakech, Maroc). Faculté
des sciences et techniques Guéliz, 2015. Français. �NNT : 2015DUNK0434�. �tel-01545833�
HAL Id: tel-01545833

https://theses.hal.science/tel-01545833v1
Submitted on 23 Jun 2017
HAL is a multi-disciplinary open access L’archive ouverte pluridisciplinaire HAL, est

archive for the deposit and dissemination of sci- destinée au dépôt et à la diffusion de documents
entific research documents, whether they are pub- scientifiques de niveau recherche, publiés ou non,
lished or not. The documents may come from émanant des établissements d’enseignement et de
teaching and research institutions in France or recherche français ou étrangers, des laboratoires
abroad, or from public or private research centers. publics ou privés.
UNIVERSITÉ CADI AYYAD UNIVERSITÉ DU LITTORAL.
FACULTÉ DES SCIENCES CÔTE D’OPALE, LMPA.
ET TECHNIQUES - CALAIS -
MARRAKECH FRANCE.
THÈSE DE DOCTORAT
Préparée dans le cadre d’une cotutelle entre

l’Université de Cadi Ayyad et l’Université du Littoral
Spécialité : Mathématiques Appliquées et Informatique
CED : Sciences de l’Ingénieur
Méthodes itératives pour la résolution

d’équations matricielles.
Thèse présentée par : SADEK EL MOSTAFA

Soutenue publiquement, le 23 Mai 2015 devant le jury :
Rapporteurs :
Abouir Jilali Professeur, FST Université Hassan II, Mohammedia, Maroc

El Hajji Said Professeur, Université Mohamed V, Rabat, Maroc
Guedda Mohammed Professeur, Université de Picardie Jules Verne, Amiens France
Examinateur :
Sadok Hassane Professeur, Université du Littoral Côte d’Opale, France
Directeurs de Thèse :
Bentbib Abdeslem Hafid Professeur, FST, Université Cadi Ayaad, Marrakech

Jbilou Khalide Professeur, Université du Littoral Côte d’Opale, France
1
Cette thèse a été préparée aux :
Laboratoire de Mathématiques Appliquées et Informatique

Université Cadi Ayyad
Faculté des Sciences et Techniques.
B.P. 549, Av. Abdelkarim Elkhattabi
Guéliz Marrakech - Maroc
% 05 24 43 34 04 / 05 24 43 31 63
Fax : 05 24 43 31 70
Site : www.fstg-marrakech.ac.ma
Laboratoire de mathématiques pures et appliquées

Centre Universitaire de la Mi-Voix
Maison de la Recherche Blaise Pascal
50, Rue Ferdinand Buisson
CS 80699
Calais Cedex, France
% 03 21 46 55 90
fax : 03 21 46 55 86
Site : www-lmpa.univ-littoral.fr
EL MOSTAFA SADEK Thèse de Doctorat

FICHE PRESENTATIVE DE LA
THESE
- Nom et Prénom de l’auteur : SADEK EL MOSTAFA

- Intitulé du travail : Méthodes itératives pour la résolution d’équations matricielles.
Encadrants :
- Abdeslem Hafid BENTBIB, Professeur de l’enseignement supérieur
- Laboratoire : LAMAI, Laboratoire de Mathématiques Appliquées et Informatique.
- Khalide JBILOU, Professeur de l’enseignement supérieur

- Laboratoire : LMPA, Laboratoire de Mathématiques Pures et Appliquées.
- Lieux de réalisation des travaux :

• Faculté des Sciences et Techniques, Laboratoire de Mathématiques Appliquées et In-
formatique, Marrakech.
• Université du Littoral Côte d’Opale, Laboratoire de Mathématiques Pures et Appli-
quées, Calais France.
- Période de réalisation du travail de thèse : 4 ans
- Rapporteurs de thèse :
Abouir Jilali Professeur, FST Université Hassan II, Mohammedia, Maroc
El Hajji Said Professeur, Université Mohamed V, Rabat, Maroc
Guedda Mohammed Professeur, Université de Picardie Jules Verne, Amiens, France.
- Cadres de coopération : Cotutelle de thèse avec l’université de Littorale Côte

d’Opale, France
- Ce travail a donné lieu aux résultats suivants :
2
3
Liste des publications
* S. Agoujil, Abdeslem Hafid Bentbib, Khalide Jbilou and EL Mostafa Sadek. A Mi-
nimal residal norm method for large-scale Sylvester matrix equations. Electronic Tran-
sactions on Numerical Analysis. Volume 43, pp. 45–59, 2014.
* Abdeslem Hafid Bentbib, Khalide Jbilou and EL Mostafa Sadek. On some Krylov
subspace based methods for large-scale nonsymmetric algebraic Riccati problems, soumis.
* Abdeslem Hafid Bentbib, Khalide Jbilou and EL Mostafa Sadek. A minimal residual
method for large scale Riccati matrix equations, en préparation.
Liste des communications
* Méthode itérative pour la résolution de l’équation de Riccati, au congrès international

de la Société Marocaine de Mathématiques Appliquées (SM2A), organisé à Marrakech
en 2012.
* A minimal residual method for large scale Sylvester matrix equations, au congrès inter-
national : Numerical Analysis and Scientific Computition with Applications (NASCA13),
organisé à Calais France en 2013.
* A new projection method for solving large-scale nonsymmetric algebric Riccati equa-
tions , au congrès international JANO à EST ESSAOUIRA en 2013.
* A Minimal residal norm method for large-scale Sylvester matrix equations, au congrès
international : Modélisation et Calcul Scientifique pour L’Ingénierie Mathématiques MO-
CASIMà Marrakech en 2014.

Table des matières
Table des Matières 4
Liste des Tableaux 7
Liste des Figures 8
Notations 9
Introduction Générale 10
1 Outils de développement théoriques 13

1.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13
1.2 Matrices particulières . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13
1.3 Norme de Frobenius . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14
1.4 Formule de Sherman-Morrison Woodbury . . . . . . . . . . . . . . . . . . 14
1.5 Inverse généralisée d’une matrice . . . . . . . . . . . . . . . . . . . . . . . 14
1.6 La décomposition en valeurs singulières (SVD) . . . . . . . . . . . . . . . 15
1.7 Approximation de rang inférieur . . . . . . . . . . . . . . . . . . . . . . . 16
1.8 Le produit de Kronecker et le -produit . . . . . . . . . . . . . . . . . . . 16
1.9 Méthodes de Krylov . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18
1.9.1 Sous-espace de Krylov par blocs . . . . . . . . . . . . . . . . . . . 18
1.9.2 Sous-espace de Krylov étendu par blocs . . . . . . . . . . . . . . . 19
1.9.3 Sous-espace de Krylov étendu global . . . . . . . . . . . . . . . . . 22
1.9.4 Algorithme d’Arnoldi étendu global . . . . . . . . . . . . . . . . . 23
1.10 Contrôlabilité et l’observabilité . . . . . . . . . . . . . . . . . . . . . . . . 25
2 Méthodes de type globale-minimisation pour les équations de Lyapu-

nov de grandes dimensions 27
2.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27
2.2 Solution exacte de l’équation de Lyapunov . . . . . . . . . . . . . . . . . . 29
2.3 Méthode de Galerkin . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30
2.4 Méthode de minimisation du résidu . . . . . . . . . . . . . . . . . . . . . . 33
2.5 Méthodes itératives pour résoudre le problème réduit . . . . . . . . . . . . 34
2.5.1 LSQR global . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34
2.5.2 La méthode du gradient conjugué globale préconditionné . . . . . 38
2.6 Exemples numériques . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41
2.6.1 Exemple 1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42
2.6.2 Exemple 2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43
4
TABLE DES MATIÈRES 5
2.7 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44
3 L’Extended-Bloc Arnoldi algorithme avec minimisation de résidu pour

les équations de Sylvester creuse et de grandes taille 45
3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45
3.2 Méthode de type Galerkin . . . . . . . . . . . . . . . . . . . . . . . . . . . 47
3.4 Méthodes pour résoudre le problème de minimisation réduit . . . . . . . . 50
3.4.1 La méthode directe basée sur la formulation de Kronecker . . . . . 50
3.4.2 Méthode de Hu et Reichel . . . . . . . . . . . . . . . . . . . . . . . 51
3.4.3 La Méthode LSQR globale . . . . . . . . . . . . . . . . . . . . . . 54
3.4.4 La méthode du gradient conjugué globale préconditionné . . . . . 56
3.5 Forme factorisée de l’approximation de la solution . . . . . . . . . . . . . 58
3.6 L’algorithme GA-LRSE . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59
3.7 L’algorithme MR-LR-Sylvester . . . . . . . . . . . . . . . . . . . . . . . . 60
3.8 Résolution de l’équation de Stein non symétrique . . . . . . . . . . . . . . 60
3.8.1 Méthode de Galerkin . . . . . . . . . . . . . . . . . . . . . . . . . . 62
3.8.2 Méthode de minimisation du résidu . . . . . . . . . . . . . . . . . 65
3.9.1 Tests numérique pour l’équation de Sylvester . . . . . . . . . . . . 67
3.9.2 Tests numérique pour l’équation de Stein non symétrique . . . . . 72
3.10 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 74
4 Méthode de minimisation du résidu pour la résolution de l’équation

de Riccati continue de grande taille 76
4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76
4.2 Méthode de Galerkin . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 78
4.4 Résolution du problème de minimisation réduit . . . . . . . . . . . . . . . 81
4.5 Méthode de MINRES globale. . . . . . . . . . . . . . . . . . . . . . . . . . 84
4.6 Forme factorisée de la solution approchée . . . . . . . . . . . . . . . . . . 88
4.7 Algorithme GA-CAREs . . . . . . . . . . . . . . . . . . . . . . . . . . . . 89
4.8 Algorithme MR-CAREs . . . . . . . . . . . . . . . . . . . . . . . . . . . . 89
4.9.1 Exemple 1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 91
4.9.2 Exemple 2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 93
4.9.3 Exemple 3 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 93
4.10 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 93
5 Équation matricielle de Riccati non symétrique et application à l’équa-

tion de transport 95
5.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 95
5.2 La méthode de Newton-Krylov par blocs . . . . . . . . . . . . . . . . . . . 98
5.3 Solution approchée de rang inférieur de NAREs . . . . . . . . . . . . . . . 100
5.4 Applications à la théorie de transport . . . . . . . . . . . . . . . . . . . . 106
5.5 Résultats numériques . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 108
5.6 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 110

TABLE DES MATIÈRES 6
Bibliographie 112
Résumé 120
Abstract 122

Liste des tableaux
3.1 Résultats pour l’exemple 1.6 . . . . . . . . . . . . . . . . . . . . . . . . . . 71

3.2 Résultats de l’exemple 1.7. . . . . . . . . . . . . . . . . . . . . . . . . . . . 72
3.3 Résultats de l’exemple 2.4 . . . . . . . . . . . . . . . . . . . . . . . . . . . 74
4.1 Résultats de l’exemple 2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . 93

4.2 Résultats de l’exemple 3. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 93
5.1 Résultats pour les exemples de NAREs dans la théorie de transport avec
c = 0.5 et α = 0.5. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 109
5.2 Résultats pour les exemples de NAREs dans la théorie de transport avec
c = 0.9999 et α = 10−8 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 109
5.3 Résultats pour les exemples 4, comparaisons avec la méthode SDA. . . . . 110
7
Table des figures
2.1 Résultats de l’exemple 1 : MR(PCGG) . . . . . . . . . . . . . . . . . . . 42

2.2 Résultats de l’exemple 1 : MR(GL-LSQR) . . . . . . . . . . . . . . . . . 43
2.3 Résultats de l’exemple 2 : MR(PCGG) . . . . . . . . . . . . . . . . . . . . 43
2.4 Résultats de l’exemple 2 : MR(GL-LSQR) . . . . . . . . . . . . . . . . . . 44
3.1 Résultat pour l’exemple 1.1. . . . . . . . . . . . . . . . . . . . . . . . . . . 68

4.1 Valeurs singulières de la solution exacte. . . . . . . . . . . . . . . . . . . . 88

4.3 Résultats de l’exemple 1.2. . . . . . . . . . . . . . . . . . . . . . . . . . . . 92
5.1 Les valeurs singulières de la solution non négative minimale de l’équation

NAREs . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 101
8
Table des notations
R le corps des nombres réels.

C le corps des nombres complexes.
Rm×n l’espace des matrices à m lignes et n colonnes, à coefficients dans R.
In la matrice identité de taille n × n.
I In s’il n’y a pas de confusion.
0n la matrice nulle d’ordre n.
ek la k ième colonne de I.
A−1 l’inverse de la matrice A.
sp(A) le spectre de la matrice A.
AT la matrice transposée de A.
A−T (A−1 )T = (AT )−1 .
det(A) le déterminant de la matrice A.
λ(A) valeur propre de la matrice A.
xT le transposé du vecteur x.
⊗ le produit de Kronecker.
le produit de diamond.
Km (A, V ) le sous espace de Krylov par blocs.
e (A, V )
Km le sous espace de Krylov étendu par blocs.
Kem (A, V ) le sous espace de Krylov étendu global.
SepF (A1 , A2 ) la séparation entre A1 et A2 , SepF (A1 , A2 ) = minkXkF =1 kA1 X − XA2 k.
<, > le produit scalaire usuel.
<, >F le produit scalaire de Frobenius.
X ⊥F Y < X, Y >F = 0.
kxk la norme euclidienne du vecteur x.
kAkF la norme de Frobenius de la matrice A.
δi,j le symbole de Kronecker.
fin de la démonstration.
9
Introduction Générale
Les équations matricielles interviennent dans de nombreux domaines des mathéma-

tiques, des sciences physiques, dans différents domaines des sciences de l’ingénieur, trai-
tement du signal, restauration d’images, filtrage, réduction de modèle en théorie du
contrôle, contrôle optimal, discrétisation d’équations aux dérivées partielles, et la ré-
solution d’équations différentielles ordinaires ou aux dérivées partielles, la théorie de
transport, voir par exemple les références suivantes [21, 39, 56, 62].
En analyse numérique, pour résoudre les équations algébriques matricielles, il existe plu-
sieurs méthodes possibles : certaines sont directes et d’autres itératives. Les méthodes
directes sont très efficaces : elles donnent la solution exacte (aux erreurs d’arrondi près)
du système linéaire considéré (ou d’une équation matricielle). Pour des problèmes de
très grande taille, cependant ces méthodes directes peuvent devenir très chères, elles
ont l’inconvénient de nécessiter une assez grande place mémoire car elles nécessitent le
stockage de toute la matrice en mémoire vive. Les méthodes itératives de type Kry-
lov, employées depuis une trentaine d’années, sont généralement plus efficaces dans les
problèmes de grande taille. La plupart des méthodes itératives manipulent le système
linéaire au travers de produits "matrice-vecteur", ce qui réduit la place mémoire.
Les équations algébriques matricielles ont été obtenues à partir de la discrétisation

d’équations aux dérivés partielles, elles sont en général de grande taille. Les méthodes
itératives basées sur une technique de projection sur un sous-espaces de Krylov de dimen-
sion m (ou sur des sous-espaces de Krylov étendus) sont plus efficaces et rapides (la vi-
tesse de convergence). Ces méthodes peuvent être classées en deux catégories principales :
les méthodes de projections basées sur la condition de Petrov-Galerkin [62, 98, 100, 102]
et les méthodes du résidu minimal MR (Minimal Residual) [74, 82, 94, 95, 99, 100] (par
exemple GMRES, MINRES).
La résolution d’équations matricielles telles que l’équation de Sylvester [13, 18, 21, 39,
56, 57, 74], Lyapunov [15, 60, 61, 68, 98, 102] ou celle de Riccati [20, 55, 64, 65] a connu
un développement considérable ces vingt dernières années. Ceci est du à leurs applica-
tions dans de nombreux domaines : le contrôle optimal, le traitement de signal et la
10
TABLE DES FIGURES 11
restauration d’images, réduction de modèle en théorie du contrôle, calcul d’un contrôle

optimal en contrôle linéaire quadratique. Le cas d’équations de Riccati non symétriques
jouent aussi un rôle important dans de nombreuses applications dont principalement :
la théorie de transport [25, 38, 71, 72], théorie des jeux [1].
Pour des matrices creuses et de grandes tailles, peu de méthodes efficaces existent. Les
méthodes itératives de type de projection sur des sous-espaces de Krylov sont généra-
lement plus efficaces et rapides. Nous proposons dans ce travail de thèse d’étudier de
nouvelles méthodes itératives de type de projection sur des sous espaces de Krylov étendu
(extended Krylov subspace) [36] pour résoudre les équations matricielles suivantes : Lya-
punov, Sylvester, Stein, Riccati symétrique et Riccati non symétrique.
Cette thèse se décompose de la façon suivante :
Le premier Chapitre sera consacré aux rappels de plusieurs définitions, théorèmes et des
rappels des sous-espaces de type Krylov qui seront utilisés tout au long de cette thèse.
Le deuxième Chapitre portera sur les équations algébriques de Lyapunov AX + XAT +

BB T = 0, où A matrice de grande taille et B une matrice de rang inférieur. Nous
proposons deux méthodes itératives de projection sur des sous espaces de Krylov étendu
global. La première basée sur la condition d’orthogonalité de Galerkin et la deuxième
avec la condition de minimisation de résidu. Nous donnerons des résultats de majoration
de l’erreur et l’efficacité des méthodes proposées par des tests numériques.
Dans le Chapitre 3, nous nous intéresserons à la résolution numérique des équations de

Sylvester AX + XB = EF T , où les matrices A et B sont de grandes tailles et creuses
et les matrices E et F sont de rang inférieur. Nous proposons une nouvelle méthode de
projection sur des sous espaces de Krylov étendu par blocs. Cette méthode est basée
sur la condition de minimisation du résidu MR (Minimal Residual) [74, 82, 94, 99, 100].
Le problème réduit obtenu est résolu soit par des méthodes directes (méthode de Hu et
Reichel [57]) ou bien des méthodes itératives : Gl-LSQR, PGCG. Les approximations de
la solution sont données sous forme factorisée, nous permettant d’économiser de la place
mémoire. Des exemples numériques sont présentés montrant l’efficacité de la méthode
proposée.
Le Chapitre 4 porte sur la résolution d’équations matricielles non linéaires de Riccati

symétriques AT X + XA − XBB T + CC T = 0, qui intervient dans plusieurs problèmes,
notamment dans les problèmes de calcul de contrôle linéaire quadratique. Dans ce cha-
pitre, nous appliquons la méthode de minimisation du résidu MR (Minimal Residual) à
l’équation de Riccati. Cette méthode est basée sur des sous espaces de Krylov étendu
par blocs, l’algorithme d’Arnoldi étendu et la méthode MINRES [94] pour résoudre le
problème de minimisation projetée. Les solutions approchées (symétrique) sont données

TABLE DES FIGURES 12
T , où Z une matrice de rang inférieur. Enfin, ce chapitre se

sous la forme Xm = Zm Zm m
termine par des exemples numériques pour valider notre approche.
Le dernier chapitre de cette thèse, concerne le cas de l’équation de Riccati non symé-
trique NARE (Nonsymmetric Algebraic Riccati Equation) XCX − XD − AX + B = 0.
Nous proposons deux nouvelles méthodes itératives pour résoudre NARE. La première
est une méthode de type Newton-Krylov, basée sur la méthode de Newton et des sous-
espaces de Krylov par blocs : à chaque itération de la méthode de Newton, on résout
une équation de Sylvester de grande taille par l’utilisation de l’algorithme d’Arnoldi
par blocs. La deuxième méthode est une méthode itérative de projection sur des sous
espaces de Krylov étendus par blocs avec la condition d’orthogonalité de Galerkin. Nous
allons également présenter des résultats théoriques de la majoration de l’erreur. La per-
formance des approches proposées au regard d’autres méthodes est confirmée par des
tests comparatifs.

Chapitre 1
Outils de développement
théoriques
1.1 Introduction
Dans ce chapitre, nous introduisons les notations et les définitions qui seront utilisées
tout au long de cette thèse. Nous donnons quelques définitions importantes par la suite :
matrices particulières, la décomposition en valeurs singulières d’une matrice , puis l’in-
verse généralisée d’une matrice. Nous rappelons ensuite le produit de Kronecker et le
-produit [12, 30]. Enfin, nous donnons un rappel sur des sous espaces de type Krylov
[36, 64, 65, 100, 101] et des algorithmes pour construire une base orthonormée des sous
espaces de Krylov, voir [55, 56, 101, 102] pour plus de détails.
1.2 Matrices particulières
Soit A = (aij ) une matrice carrée de taille n × n.
Définition 1.2.1. La matrice A est dite :

– symétrique si AT = A,
– diagonale si aij = 0 pour i 6= j,
– triangulaire inférieure si aij = 0 pour i < j,
– triangulaire supérieure si aij = 0 pour i > j,
– orthogonale si AT A = AAT = In ,
– Hessenberg supérieure si aij = 0 pour i > j + 1,
– définie positive si (Ax, x) > 0 pour tout x ∈ Rn /{0},
– non négative si tous les éléments aij de A sont positifs et aij 6= 0,
– Une matrice A dont tous les termes extra-diagonaux sont négatifs ou nuls est appelée
Z-matrice.
– A est une M -matrice si A est une Z-matrice inversible et non négative
13
1.3 Norme de Frobenius 14
1.3 Norme de Frobenius
Pour deux matrices Y et Z ∈ Rn×s , le produit scalaire de Frobenius de Y et Z est défini

par
hY, ZiF = tr(Y T Z),
où tr(Y T Z) désigne la trace de la matrice Y T Z. La norme associée est celle de Frobenius

noté par k . kF , et définie par
n X
X m
k A kF = tr(AT A) = (ai,j )2 , où A = [ai,j ] ∈ Rn×m .
i j
1.4 Formule de Sherman-Morrison Woodbury
Proposition 1.4.1. Soient B une matrice inversible d’ordre n et u, v ∈ Rn . Alors la

matrice B + uv T est inversible si et seulement si 1 + v T B −1 u 6= 0. Dans ce cas, nous
avons la formule suivante :
(B −1 u)v T
(B + uv T )−1 = (I − )B −1 (1.1)
1 + v T B −1 u
La formule généralisée, connue sous le nom de formule de Sherman-Morrison-Woodbury,

est définie dans la proposition suivante.
Proposition 1.4.2. Soient A ∈ Rn×n et C ∈ Rm×m deux matrices inversibles. Soient

U ∈ Rn×m et V ∈ Rm×n . Si la matrice (C −1 +V A−1 U ) est inversible, alors (A+U C −1 V )
est inversible, et vérifie la formule de Sherman-Morrison-Woodbury donnée comme suit :
(A + U C −1 V )−1 = A−1 − A−1 U (C −1 + V A−1 U )−1 V A−1
1.5 Inverse généralisée d’une matrice
Définition 1.5.1. La pseudo-inverse, A+ , d’une matrice A de taille m×n, à coefficients

réels ou complexes est l’unique matrice de taille n × m satisfaisant les quatre critères
suivants :
1. AA+ A = A,
2. A+ AA+ = A+ ,
3. (AA+ )T = AA+ (AA+ est une matrice symétrique),
4. (A+ A)T = A+ A (A+ A est également symétrique).

1.6 La décomposition en valeurs singulières (SVD) 15
Remarque 1.5.2. Pour les matrices dont les composantes sont des nombres complèxes
au lieu des nombres réels, AT sera remplacée par A∗ , où A∗ désigne l’adjointe de la
matrice A. La pseudo-inverse appelée aussi inverse au sens de Moore-Penrose existe et
unique pour toute matrice réelles ou complexe.
Nous donnons la proposition suivante
Proposition 1.5.3 ([43]). Le pseudo-inverse d’une matrice A ∈ Rm×n de rang maximal

est définie par :
– A+ = AT (AAT )−1 si rang(A) = m,
– A+ = (AT A)−1 AT Si rang(A) = n.
1.6 La décomposition en valeurs singulières (SVD)
En mathématiques, le procédé d’algèbre linéaire de décomposition en valeurs singulières

(ou la décomposition SVD : Singular Value Decomposition) d’une matrice est un outil
important de factorisation de matrices rectangulaires réelles ou complexes, elle est l’une
des méthodes de factorisation la plus générale et la plus utile. La décomposition en
valeurs singulières est utilisée, entre autres, dans le calcul de pseudo-inverse, dans la
résolution de problèmes aux moindres carrés et dans les méthodes de régularisation.
On l’emploie également en traitement d’image, en traitement du signal, et pour les
approximations de rang inférieur d’une matrice.
Théorème 1.6.1 (Golub et Van Loan, 1990). Soit A une matrice réelle de dimension
m × n. Il existe deux matrices orthogonales U et V respectivement de dimension m × m
et n × n tel que
A = U ΣV T , (1.2)
où Σ = diag[σ1 , σ2 , ..., σm ] ∈ Rm×n . Les σi sont appellés les valeurs singulières de A et

elles vérifient σ1 ≥ σ2 ≥ ... ≥ σr > σr+1 = ... = σm = 0, r ≤ min{m, n}. Les vecteurs
ui (colonnes de U) et vi (colonnes de V) sont appellés les vecteurs singuliers de gauche
respectivement de droite.
Remarque 1.6.2. Le rang de la matrice A est égal au nombre de valeurs singulières

non nulles que possède la matrice A.
La décomposition SVD d’une matrice A de taille m × n nécessite 4m2 n + 8mn2 + 9n3
opérations élémentaires. Avec MATLAB on obtient la décomposition singulière avec la
commande[U,S,V]=svd(A).
Nous donnons le théorème suivant
Théorème 1.6.3 ([43]). Soit A ∈ Rm×n une matrice de rang r, et soit sa décomposition
en valeurs singulières A = U ΣV T , où Σ = diag[σ1 , σ2 , ..., σr , 0, ..., 0] ∈ Rm×n . Alors le
1.7 Approximation de rang inférieur 16
pseudo-inverse de A est donné par :
A+ = V Σ+ U T ,
où Σ+ = diag[1/σ1 , 1/σ2 , ..., 1/σr , 0, ..., 0] ∈ Rn×m .
1.7 Approximation de rang inférieur
Soit A ∈ Rm×n une matrice. Déterminer une matrice X de même taille que A mais de
rang inférieur k de telle sorte que kA − XkF soit minimale, est un problème classique.
La solution de ce problème est donnée par le théorème d’Eckart-Young suivent
Théorème 1.7.1 (Eckart et Young). Soit A une matrice de rang r, et soit sa décom-
position en valeurs singulières
A = U ΣV T ,
où Σ = diag[σ1 , σ2 , ..., σr , 0, ..., 0] ∈ Rm×n , et soit k un entier inférieur ou égal à r. Si

on note Σk = diag[σ1 , σ2 , ..., σk ], alors nous avons
 1
r 2
X
min kA − XkF = kA − Ak kF =  σi2 (A) ,
rang(X)≤k
i=k+1
où !
Σk 0
Ak = U V T.
0 0
Ce théorème établit une relation entre le rang k de l’approximation X de A et les

valeurs singulières σk+1 , ..., σr de A. Par conséquent, si les valeurs singulières décroissent
rapidement vers zéro alors nous pouvons espérer déterminer une approximation de A
possédant un rang très faible.
1.8 Le produit de Kronecker et le -produit
Le produit de Kronecker est un outil important dans l’algèbre linéaire et les équations
matricielles. Il permet de transformer une équation matricielle en un système linéaire,
pour plus de détails voir par exemple le livre [96]. Le -produit est défini dans [12, 30],
c’est un outil important pour résoudre les équations matricielles par les méthodes de
type Krylov global, voir [12, 30].
Définition 1.8.1. Soient A et B deux matrices de taille m × n et s × t respectivement.

Le produit de Kronecker A ⊗ B est la matrice de taille ms × tn définie par :

1.8 Le produit de Kronecker et le -produit 17
 
a1,1 B ··· ··· a1,n B
..
 
 a2,1 B . a2,n B
 

A⊗B =
 .. .. ..



 . . . 

am,1 B am,2 B · · · am,n B
Pour Y = [yi,j ] ∈ Rn×s ; on note par vec(Y ) le vecteur de Rns défini par vec(Y ) =
[y(., 1)T , (y., 2)T , ..., y(., s)T ] où y(., j), j = 1, ..., s est la jème colonne de Y .
On a les propriétés suivantes :
Proposition 1.8.2 ([12, 30]). Soient A et B deux matrices de taille m × n et s × t

réspctivement. Nous avons les proprétés suivantes :
1. (A ⊗ B)T = AT ⊗ B T .
2. (A ⊗ B)(C ⊗ D) = (AC ⊗ BD).
3. Si A et B sont des matrices non singulières de dimension n × n et p × p, respecti-
vement, alors (A ⊗ B)−1 = A−1 ⊗ B −1 .
4. Si A ∈ Rn×n et B ∈ Rp×p , alors
det(A ⊗ B) = det(A)p det(B)n et tr(A ⊗ B) = tr(A)tr(B).
5. vec(ABC) = (C T ⊗ A)vec(B).
6. vec(A)T vec(B) = trace(AT B).
On rappelle le produit diamond noté par et qui est défini de la façon suivante :
Définition 1.8.3 ([12, 30]). Soient A = [A1 , ..., Ap ] et B = [B1 , ..., Bl ] deux matrices
de taille n × ps et n × ls, respectivement, où Ai et Bj (i = 1, ..., p; j = 1, ..., l) sont des
matrices de Rn×s . Alors la matrice AT B ∈ Rp×l est définie par :
 
hA1 , B1 iF hA1 , B2 iF · · · hA1 , Bl iF

.. 
 hA2 , B1 iF .
 
T

A B =
 .. ..



 . . 

hAp , B1 iF ··· · · · hAp , Bl iF
Remarque 1.8.4.
1. Si s = 1 alors AT B = AT B.
2. Si s = 1, p = 1 et l = 1, alors on pose A = u ∈ Rn et B = v ∈ Rn , on aura
AT B = uT v ∈ R.
3. La matrice A = [A1 , ..., Ap ] est F −orthonormale si et seulement si AT A = Ip .
4. Si X ∈ Rn×s , alors X T X =k X k2F .

1.9 Méthodes de Krylov 18
On rappelle aussi les propriétés suivantes satisfaites par le −produit.
Proposition 1.8.5 ([12, 30]). Soient A, B, C ∈ Rn×ps , D ∈ Rn×n , L ∈ Rp×p et α ∈ R.

Alors on a :
1. (A + B)T C = AT C + B T C.
2. AT (B + C) = AT B + AT C.
3. (αA)T C = α(AT C).
4. AT B)T = B T A.
5. (DA)T B = AT (DT B).
6. AT (B(L ⊗ Is )) = (AT B)L.
7. k AT B kF ≤k A kF k B kF .
1.9 Méthodes de Krylov
Les méthodes de projection sur les sous-espace de Krylov sont très utilisé dans dans des
problèmes de grande dimension. Ces méthodes se basent sur des techniques de projection
sur un sous-espace, appelé sous-espace de Krylov. Les méthodes itératives de type Krylov
a connu un développement considérable dans les dernières années, pour résoudre les
équations matricielles de grande taille.
Dans cette section, on va rappeler le sous-espace de Krylov par blocs, sous-espace de

Krylov étendu par blocs ou global. On va rappeler aussi des algorithmes de type Arnoldi,
pour construire une base orthogonale ou orthonormée d’un sous espace de type Krylov.
1.9.1 Sous-espace de Krylov par blocs
Le sous espace de Krylov par blocs associé à la paire (A, V ) est défini par
Km (A, V ) = Image({V, AV, A2 V, · · · , Am−1 V }).
L’algorithme d’Arnoldi par blocs ci-dessus consiste à construire une base orthonormale
de l’espace Km (A, V ).

Algorithm 1 Algorithme d’Arnoldi par blocs

– Calculer la décomposition QR de V : [V1 , R] =qr(V )
– Pour j = 1 : m
– W = A Vj
– Pour i = 1, . . . , j
– Hi,j = ViT W
– W = W − Vi Hi,j
– Fin Pour i
– Calculer la décomposition QR de W : [Vj+1 , Hj+1,j ] =qr(W )
– Fin Pour j
On note Vm la matrice réelle de taille n × mr définie par
Vm = [V1 , . . . , Vm ]
et H̃m la matrice réelle de taille (m + 1)r × mr de type Hessenberg supérieure par blocs,
dont les blocs non nuls Hi,j de taille r × r sont produits par l’algorithme précédent. On
vérifie les identités suivantes
T
A Vm = Vm Hm + Vm+1 Hm+1,m Em
= Vm+1 H̃m ,
où Em = [0r , . . . , 0r , Ir ] est la matrice de taille mr × r composée des r dernières colonnes

de la matrice identité Imr et Hm est la matrice de taille mr × mr obtenue en supprimant
les r dernières lignes de H˜m .
1.9.2 Sous-espace de Krylov étendu par blocs
Pour la résolution des systèmes linéaires, il existe plusieurs méthodes possibles ; cer-
taines sont directes et d’autres itératives. Dans des problèmes de très grande taille, ces
méthodes peuvent devenir très chères, et sont limitées par la taille du problème et bien
que les dernières versions de Matlab prennent en charge des problèmes de taille allant
jusqu’à O(104 ), le temps CPU ( Central Processing Unit) ainsi que la taille mémoire
peuvent poser problème. Les méthodes efficaces sont les méthodes de sous-espace Kry-
lov et consistent à construire une suite de sous espaces emboîtés de taille croissante qui
vont contenir à la limite la solution recherchée. Elles sont basées sur une technique de
projection sur un sous espace de Krylov, de dimension m (m ≤ n) plus petite que la
taille n du problème. Parmi les nombreux ouvrages de référence traitant des méthodes
de Krylov, on pourra citer le livre de Yousef Saad [100]. Ces méthodes diffèrent l’une de
l’autre par le type de projection qui est appliquée et le choix des sous-espaces de Krylov.

Dans cette section, on va rappeler le sous-espace de Krylov étendu par blocs.
Supposons que nous calculons des approximations des expressions de la forme w :=

f (A)v, où f est une fonction non linéaire définie sur A.
Dans [36], Druskin et Knizhnerman ont montré que f ne peut pas être approchée avec
précision par un polynôme de degré m − 1. Pour cette raison, ils ont proposé l’utilisation
de la méthode de sous-espace de Krylov étendu, ce qui permet d’approcher f par une
fonction rationnelle.
Le sous-espace de projection que l’on considère peut être considéré comme une somme
de deux sous-espaces de Krylov par blocs. Le premier lié à Km (A, V ) et le deuxième à
Km (A−1 , V ). Plus précisément, en supposant que la matrice A est inversible, le sous-
espace de Krylov étendu par blocs associé à (A, V ) est donné par [36, 102]
e
Km (A, V ) = Image{V, A−1 V, AV, A−2 V, A2 V, · · · , Am−1 V, A−m+1 V }
= Km (A, V ) + Km (A−1 , A−1 V ).
nous avons,
e (A, V ) ⊂ Ke
Km e e
m+1 (A, V ) et Km (A, V ) ⊂ AKm (A, V ).
Les méthodes de projection sont basées sur des algorithme qui permettent de construire
une base orthonormée de l’espace de Krylov. parmi ces algorithmes, le plus connu est
l’algorithme d’Arnoldi qui est basée sur d’orthonormalisation de Gram-Schmidt.
e (A, V ) en ap-
Pour construire une base orthonormée de l’espace de Krylov étendu Km
plique l’algorithme d’Arnoldi étendu par blocs ci-dessous ([55, 102] ).
Après m étapes de l’algorithme d’Arnoldi étendu par blocs, on obtient la matrice
Vm = [V1 , V2 , . . . , Vm ] où Vi ∈ Rn×2r ,
et les matrices de Hessenberg par blocs :
 
H1,1 · · · ··· ··· H1,m
. ..
 
H2,1 . .
 

 . 

 .. .. .. 
 0 . . . 
Hm = 
 
.. .. .. .. .. 

 . . . . .


..
 
 .. .. 

 . . . Hm,m 

0 ··· ··· 0 Hm+1,m

Algorithm 2 Algorithme d’Arnoldi étendu par blocs (EBA)

1. Entrées : les matrices A ∈ Rn×n et V ∈ Rn×r et m entier.
2. Calculer la décomposition QR de [V, A−1 V ], c’est-à-dire, [V, A−1 V ] = V1 Λ ;
3. On pose V0 = [ ] ;
4. Pour j = 1, 2, ..., m faire
(1) (2)
(a) On prend Vj = Vj (:, 1 : r) et Vj = Vj (:, r + 1 : 2r).
h i
(1) (2)
(b) Vj = [Vj−1 , Vj ] ; V̂j+1 = A Vj , A−1 Vj ;
(c) Orthogonalisation de V̂j+1 avec les Vi , i = 1, 2, ..., j, pour obtenir Vj+1 , c’est-
à-dire,
i. Pour i = 1, 2, . . . , j faire
ii. Hi,j = ViT V̂j+1 ;
iii. V̂j+1 = V̂j+1 − Vi Hi,j ;
iv. Fin Pour
(d) Calculer la décomposition QR de V̂j+1 , c’est-à-dire, V̂j+1 = Vj+1 Hj+1,j ;
5. Fin Pour
et  
H1,1 · · · · · · ··· H1,m
 . .. 
H2,1 . . .
 
 
.. .. ..
 
Hm =  . .
 
 0 . 

 .. .. .. .. .. 

 . . . . . 

0 0 Hm,m−1 Hm,m
La matrice Hm de taille 2mr × 2mr est obtenue en supprimant la dernière ligne par
blocs de Hm . Par construction, la matrice Vm est orthonormale, ce qui signifie que les
matrices V1 , V2 , · · · , Vm constituent une base orthonormale de l’espace de Krylov étendu
e (A, V ). Soit T ∈ R2mr×2mr la restriction de la matrice A sur le sous-espace de Krylov
Km m
e (A, V ), c’est-à-dire T
étendu Km T
m = Vm A Vm . Dans [102], Simoncini a démontré que
la matrice Tm est également matrice de Hessenberg supérieure par blocs, et on peut
calculer la matrice Tm à partir Hm .
On établit les identités suivantes de manière immédiate
Proposition 1.9.1. [55] Soit T̄m = VTm+1 AVm , alors nous avons
AVm = Vm+1 T̄m

= Vm Tm + Vm+1 Tm+1,m ETm .
où Ti,j ∈ R2r×2r est le bloc (i, j) de la matrice Tm et Em est la matrice de taille

mr × r constituée des r dernières colonnes de la matrice identité Imr (c.à.d Em =
[02(m−1)r×2r , I2r ]T ).

1.9.3 Sous-espace de Krylov étendu global
Le sous-espace de Krylov étendu global associé à (A, V ) est donne par [36, 56, 102]
g
Km (A, V ) = vect{V, A−1 V, AV, A−2 V, A2 V, · · · , Am−1 V, A−m+1 V }
= Kgm (A, V ) + Kgm (A−1 , A−1 V ).
g (A, V ),
Pour construire une base une base orthonormale {V1 , V2 , ..., Vm } de l’espace Km
on a besoin de définir La factorisation QR global.
La factorisation QR global
Soit Z = [Z1 , ..., Zk ] ∈ Rn×ks ( avec Zi ∈ Rn×s , i = 1, ..., k). La décomposition glo-
bale de QR de Z est basée sur la procédure de Gram-Schmidt globale, elle permet
de calculer une matrice F -orthonormale Q = [Q1 , ..., Qk ] de taille n × ks telle que
vect{Q1 , ..., Qk } = vect{Z1 , ..., Zk } avec hQi , Qi iF = 1 et hQi , Qj iF = 0 si i 6= j.
L’algorithme est décrit ci-dessous :
Algorithme 3 : L’algorithme de Gram-Schmidt modifié global (GGSG)
Entrée : Z = [Z1 , ..., Zk ] une matrice réelle de dimension n × sk.

Sortie : Q = [Q1 , ..., Qk ] ∈ Rn×sk F-orthonormale et R ∈ Rk×k triangulaire supérieure.
1. R = (ri,j ) = 0.
2. r1,1 =k Z1 kF .
3. Q1 = Z1 /r1,1 .
4. Pour i = 2, ..., k, faire :
5. Q = Zi ,
6. Pour j = 1, ..., i − 1, faire :
7. rj,i = hQ, Zj iF ,
8. Q = Q − rj,i Qj ,
9. Fin Pour j.
10. ri,i =k Q kF ,
11. Qi = Q/ri,i ,
12. Fin Pour i.

Proposition 1.9.2. [55] Soit Z = [Z1 , ..., Zk ] une matrice de taille n × ks avec Zi ∈
Rn×s , pour i = 1, ..., k. Donc en utilisant l’algorithme 3, la matrice Z peut être factorisée
sous la forme suivante
Z = Q(R ⊗ Is ),
avec Q = [Q1 , ..., Qk ] une matrice F -orthonormale de taille n × ks et qui satisfait QT

Q = Ik et R une matrice triangulaire supérieure de taille k × k.
On note que QT Z = QT (Q(R ⊗ Is )), donc en utilisant la Proposition 1.8.5, on aura

QT Z = R.
1.9.4 Algorithme d’Arnoldi étendu global
Dans ce paragraphe, on va décrire la procédure d’Arnoldi étendu global pour construire

g (A, V ) avec v ∈ Rn×s .
une base F -orthonormale {v1 , ..., v2m } de Km i
Les deux blocs v1 et v2 sont calculés de la manière suivante :
V
v1 = , et ω2,2 v2 = A−1 V − ω1,2 v1 , (1.3)
ω1,1
où les paramètres ω1,1 et ω2,2 sont tels que la matrice [v1 , v2 ] de dimension n × 2s est
F -orthonormale. Donc ω1,1 =k V kF , ω1,2 = tr(V1T (A−1 V )) et ω2,2 =k U kF , avec
U = A−1 V − ω1,2 v1 .
Pour calculer les blocs v2j+1 et v2j+2 , pour j = 1, ..., m − 1, on utilise les formules
suivantes : ( P2j
h2j+1,2j−1 v2j+1 = Av2j−1 − i=1 hi,2j−1 vi ,
P2j+1 (1.4)
h2j+2,2j v2j+2 = A−1 v2j − i=1 hi,2j vi .
En utilisant les conditions d’orthogonalités v2j+1 ⊥F v1 , ..., v2j et v2j+2 ⊥F v1 , ..., v2j+1 , on
obtient :
j j
hi,2j−1 = tr(viT Av2j−1 ) = tr(viT Ui,1 ) et hi,2j = tr(viT A−1 v2j ) = tr(viT Ui,2 ),
où
i−1 i−1
j j
A−1 v2j −
X X
Ui,1 = Av2j−1 − hk,2j−1 vk et Ui,2 hk,2j vk .
k=1 k=1
Les paramètres h2j+1,2j−1 et h2j+2,2j sont tels que k v2j+1 kF = 1 et k v2j+2 kF = 1

respectivement. Donc,
j j
h2j+1,2j−1 =k U2j+1,1 kF et h2j+2,2j =k U2j+1,2 kF .
Maintenant, afin de décrire une nouvelle méthode pour calculer les vecteurs colonnes de
la matrice Vm = [vi ]i=1,...,2m de dimension n×2ms et les entrées non nulles de la matrice
bloc de Hessenberg supérieure Hm = [hi,j ]j=1,...,2m

i=1,...,2m , on introduit
!
(1) (2) h2i−1,2j−1 h2i−1,2j
Vi = [v2i−1 , v2i ] = [Vi , Vi ], et Hi,j = , (1.5)
h2i,2j−1 h2i,2j
qui sont, respectivement, le ième vecteur bloc de dimension n × 2s de la matrice Vm ,

et le (i, j)-ème bloc de taille 2 × 2 de la matrice bloc de Hessenberg supérieure Hm . En
utilisant le produit de Kronecker et les formules (1.4) et (1.5), on obtient
j
(1) (2)
Vj+1 (Hj+1,j ⊗ Is ) = [AVj , A−1 Vj ] −
X
Vi (Hi,j ⊗ Is ).
i=1
Comme les blocs V1 , ..., Vm sont orthogonaux par rapport au −produit (i.e., ViT Vj = 02
pour i 6= j), alors en utilisant les propriétés du −produit, on trouve
(1) (2)
Hi,j = ViT [AVj , A−1 Vj ] = ViT Uij , i = 1, ..., j,
où Uij =
Pi−1
k=1 Vk (Hk,j ⊗ Is ). Pour les blocs triangulaires supérieurs Hj+1,j , on rappel
que les deux matrices Hj+1,j et Vj+1 sont calculées telles que Vj+1 Vj+1 = I2 . Donc on
j
obtient Vj+1 et Hj+1,j en calculant la factorisation QR globale de Uj+1 .
Soit Ω = [ωi,j ] une matrice triangulaire supérieure de dimension 2×2 telle que ses entrées
non nulles sont définies en (1.3), et on peut vérifier facilement que Ω, et V1 le premier
bloc de la matrice Vm peuvent être obtenus facilement par la décomposition QR globale
de [V, A−1 V ].
Finalement, on peut résumer les résultats précédents dans l’algorithme suivant :
Algorithme 4 : Algorithme d’Arnoldi étendu global (EGA)
Entrée : A ∈ Rn×n , V ∈ Rn×s et m .

Sortie : Vm = [V1 , ..., Vm ] de dimension n×2ms et Hm ∈ R2m×2m matrice de Hessenberg
supérieure.
1. On calcule la factorisation QR globale de [V, A−1 V ], i.e , [V, A−1 V ] = V1 (Ω ⊗ Is ) ;
2. On pose V0 = [ ] ;
3. Pour j = 1, ..., m, faire :
(1) (2)
4. On prend Vj = Vj (:, 1 : s) ; Vj = Vj (:, s + 1 : 2s) ;
(1) (2)
5. Vj = [Vj−1 , Vj ]; U = [AVj , A−1 Vj ] ;
6. Pour i = 1, ..., j, faire :
7. Hi,j = ViT U ;
8. U = U − Vi (Hi,j ⊗ Is ) ;
9. Fin(i).
1.10 Contrôlabilité et l’observabilité 25
10. On calcule la décomposition QR globale de U , i.e., U = Vj+1 (Hj+1,j ⊗ Is ) ;

11. Fin(j).
Si les matrices triangulaires supérieures Hj+1,j ne sont pas de rang maximal, l’algorithme
1.3 calcule une matrice F −orthonormale Vm avec Vm = [V1 , ..., Vm ] et Vi ∈ Rn×2s , (i =
1, ..., m) et Hm une matrice bloc de Hessenberg supérieure de dimension 2(m + 1) × 2m
avec Hm = [Hi,j ] et Hi,j ∈ R2×2 .
On introduit une matrice de dimension 2m × 2m donnée par Tm = VTm (AVm ) = [Ti,j ],
où Ti,j = ViT (AVj ) ∈ R2×2 , pour i, j = 1, ..., m.
Soit Tm = VTm+1 (AVm ) et Em
T = [0
2×2(m−1) , I2 ] une matrice qui correspond aux
deux dernières colonnes de la matrice identité 2m × 2m, en utilisant le même argument
utilisé dans le cas d’Arnoldi extended par bloc, on peut montrer qu’après m itérations
de l’algorithme 3, on aura :
T
AVm = Vm+1 (Tm ⊗ Is ) = Vm (Tm ⊗ Is ) + Vm+1 (Tm+1,m Em ⊗ Is ).
La proposition suivante, donnée en [55], permet de calculer Tm sans utiliser le produit

matrice-vecteur avec A.
Proposition 1.9.3. ([55]) Soit Tm = [t:,1 , ..., t:,2m ] et Hm = [h:,1 , ..., h:,2m ] où t:,i , h:,i ∈
R2(m+1) sont les ièmes colonnes de deux matrices blocs de Hessenberg supérieures Tm et
Hm de dimension 2(m + 1) × 2m. Donc les colonnes impaires sont telles que
t:,2j−1 = h:,2j−1 , pour j = 1, ..., m,
tandis que les colonnes paires satisfont
1 2(m+1)
t:,2 = (ω1,1 e1 − ω1,2 h:,1 ),
ω2,2
1 2(m+1)
t:,2j+2 = (e − t:,1:2j+1 h1:2j+1,2j ) pour j = 1, ..., m − 1,
h2j+2,2j 2j
(k)
où ei est la ième colonne de la matrice identité Ik et ω1,1 , ω1,2 sont déjà définis en
(1.3).
1.10 Contrôlabilité et l’observabilité
Dans cette partie, nous rappelons les résultats relatifs à la contrôlabilité et à l’observa-
bilité de la paire (A, B) et de la paire (C, A) définies par un système dynamique linéaire
à coefficients constants du type :
(
ẋ = Ax + Bu,
(1.6)
y = Cx,
1.10 Contrôlabilité et l’observabilité 26
Ces notions sont largement développées dans [111]. Nous rappelons la notion de matrices
stables dans la définition suivante :
Définition 1.10.1. Une matrice A est dite stable, si
Sp(A) ⊂ C− .
Définition 1.10.2. La paire (A, B) est dite stabilisable s’il existe une matrice F telle
que A + BF est stable.
Remarque 1.10.3. Un système dynamique linéaire à coefficients constants, donné par

l’équation (1.6) est dit asymptotiquement stable, si A est stable.
Les définitions de la stabilisation et de la détectabilité du système linéaire (1.6) sont

données dans [31].
Définition 1.10.4. La paire (C, A) est dite détectable si (AT , C T ) est stabilisable.

Chapitre 2
Méthodes de type
globale-minimisation pour les
équations de Lyapunov de grandes
dimensions
2.1 Introduction
les équations matricielles de Lyapunov interviennent dans plusieurs domaines comme

dans la théorie du contrôle [42], la théorie de la communication, l’analyse de la stabi-
lité des systèmes dynamiques [87], la stabilité des équations différentielles ordinaires,
la réduction de modèles [4, 5, 53], et la restauration d’images [18]. Par exemple, nous
P
considérons le système linéaire dynamique invariant dans le temps :
(
ẋ(t) = Ax(t) + Bu(t)
Σ := (2.1)
y(t) = Cx(t)
où A est une matrice creuse de grande taille n × n.
Le problème de réduction consiste à approcher Σ par
( dx̂(t)
dt = Âx̂(t) + B̂u(t)
Σ̂ :=
ŷ(t) = Ĉ x̂(t),
où Â ∈ Rk×k , B̂ ∈ Rk×s et Ĉ ∈ Rk avec k << n.
Parmi les méthodes de réduction des modèles de grande dimension, la plus connue est
la troncation équilibrée. Cette méthode est basée sur le calcul des grammiens de com-
mandabilité et d’observabilité. Les grammiens de commandabilité Wc et d’observabilité
27
2.1 Introduction 28
Wo sont obtenus par résolution des équations de Lyapunov suivantes :
AWc + Wc AT + BB T = 0,
AT Wo + Wo A + C T C = 0.
Il a été démontré que [42] :
k Σ − Σ̂ k∞ ≤ 2(σk+1 + ... + σn )
où les σi sont Les valeurs singulières de Hankel associées à Σ donnée par : σi =

p
λi (Wc Wo ) arrangées en ordre décroissant, et la norme k Σ k∞ est la plus grande
valeur singulière de la fonction de transfert associée à Σ et définie par k Σ k∞ =
supω∈R (σmax (G(jω)), où G(.) représente la fonction de transfert associée à Σ et jω
varie sur la totalité de l’axe imaginaire. Pour plus de détails voir [42, 68]. L’équation
matricielle de Lyapunov est de la forme suivante
AX + XAT + BB T = 0, (2.2)
où A, B et X sont des matrices respectivement de taille n × n, n × s, n × n avec s << n.
L’équation matricielle de Lyapunov peut être formulée comme un système linéaire de

taille n2 × n2 . En utilisant le produit de Kronecker, nous obtenons
(In ⊗ A + A ⊗ In )vec(X) + vec(BB T ) = 0. (2.3)
Donc l’équation de Lyapunov possède une solution unique X si et seulement si λi (A) +

λj (A) 6= 0 pour tout i, j = 1, 2, · · · , n
Pour les équations matricielles de lyapunov de petite taille, il existe plusieurs méthodes
dans la littérature, comme par exemple, l’algorithme de Bartels-Stewart [11], Hessenberg-
Schur [45].
Pour la résolution de l’équation de Lyapunov (2.2), nous distinguons deux catégories.

La première se sont Les méthodes directes, comme par exemple, l’algorithme de Bartels-
Stewart [11], Hessenberg-Schur [45]. Ces méthodes ne sont pas utilisables lorsque la
matrice A est de grande taille ou creuse.
La deuxième catégorie est constituée des méthodes itératives pour les problèmes de
grande taille. La plupart de ces méthodes sont les méthodes itératives de projections
sur des sous-espaces de type Krylov, basées sur les algorithmes de type Arnoldi [36,
55, 68, 100, 102] pour construire une base orthonormée de l’espace de Krylov. Une
caractéristique commune de toutes ces méthodes de projection est que les approximations
Xm de la solution exacte de l’équation de Lyapunov sont données sous forme factorisée,
T avec Z une matrice rang inférieur. De plus il y a la possibilité d’obtenir
Xm = Zm Zm

2.2 Solution exacte de l’équation de Lyapunov 29
une bonne approximation Xm de rang inférieur (low rank approximate solution), il est
soutenu par des résultats théoriques récents montrant la décroissance rapide des valeurs
propres vers zéro de la solution exacte X, voir [61, 68, 82, 98, 102]. Cette approche est
très importante lorsqu’il s’agit des problèmes de grande taille, puisque stocker la matrice
Xm complète en général nécessiterait une quantité significative de l’espace mémoire.
Nous proposons dans cette partie une nouvelle méthode pour résoudre l’équation de
Lyapunov de grande taille déterminer l’approximation de la solution de l’équation de
2.2 Solution exacte de l’équation de Lyapunov
Dans cette section, nous allons donner l’expression exacte de la solution X de l’équation
de Lyapunov. Supposons que l’équation de Lyapnuov admette une solution unique X.
Il a été montré dans [54, 105] que l’expression de la solution est la suivante
p X
X p
X= γi,j Ai−1 BB T (AT )j−1 ,
j=1 i=1
où p le degré du polynôme minimal de A pour la matrice B (P (A)B = 0) et la matrice

T = (γi,j )1≤i,j≤p est la solution d’une équation matricielle de Lyapnuov de petite taille
(Voir [103]).
Jbilou et Riquet ont proposé, dans [68], une méthode permettant de calculer les coeffi-
cients γi,j à partir d’une base F -orthonormée de l’espace de Krylov
Kp (A, B) = span{B, AB, A2 B, ..., Ap−1 B},
construit par l’algorithme d’Arnoldi global. Ils ont aussi montré que la solution X est
de la forme
X = Vp (T ⊗ Is )VpT , (2.4)
où Vp est la matrice F-orthogonale construite par l’algorithme d’Arnoldi global associé

Kp (A, B).
Comme la solution exacte est donnée par l’expression (2.4), l’approximation Xm de la

solution exacte est définie par :
T
X = Vm (T ⊗ Is )Vm . (2.5)

2.3 Méthode de Galerkin 30
2.3 Méthode de Galerkin
Dans cette section, nous présentons une méthode de projection pour résoudre l’équation
de Lyapunov de grande taille. En utilisant des projections dans des sous espaces de
Krylov étendus (extended) globale avec la condition d’orthogonalité de Galerkin.
Nous allons chercher des approximations Xm de la solution exacte X de l’équation de

Lyapunov sous la forme
GA T
Xm = Vm (Ym ⊗ Is )Vm , (2.6)
où Vm est la matrice F-orthogonale construite en appliquant m itérations de l’algorithme

e (A, B). Donc on a les relations suivantes
d’Arnoldi étendu global à Km
T
Tm = Vm (AVm ).
La condition d’orthogonalité de Petrov-Galerkin permet de donner l’équation de Lyapu-

nov projetée suivante :
Tm YmGA + YmGA TTm + kBk2F e2m 2m T

1 (e1 ) = 0. (2.7)
La solution de l’équation projeté au-dessus peut être obtenu par une méthode directe
comme la méthode Hessenberg-Schur [45].
Pour réduire le calcul dans le test d’arrêter les itérations de la méthode proposée, nous
donnons une borne supérieure (approximation supérieure) de la norme de Frobenius du
résidu R(Xm ).
Théorème 2.3.1. Soient Ym la solution exacte de l’équation de Lyapunov projetée (2.7)

GA ⊗I )V T la solution approchée de l’équation de Lyapunov 2.2. La norme
et Xm = Vm (Ym s m
de Frobenius vérifie l’inégalité
q
GA
k R(Xm ) kF ≤ 2(m + 1)α := rm , (2.8)
où α =k Tm+1,m Ym kF
Démonstration. A l’étape m, nous avons :
GA T GA T T
R(Xm ) = AVm (Ym ⊗ Is )Vm + Vm (Ym ⊗ Is )Vm A + BB T
GA T GA TT
= Vm+1 (Tm ⊗ Is )(Ym ⊗ Is )Vm + Vm (Ym ⊗ Is )(Tm ⊗ Is )Vm+1 + k B k2F V1 V1T .
En utilisant la relation
(2m+2)
V1 = Vm+1 [e1 ⊗ Is ],

et le fait que
Vm = [V1 , V2 , V3 , ..., Vm ] (2.9)

 
I2s
 

 I2s 


= [V1 , V2 , V3 , ..., Vm+1 ]  .. 
(2.10)
 . 

 

 I2s 

02s
 
[I2 ⊗ Is ]
 

 [I2 ⊗ Is ] 


= [V1 , V2 , V3 , ..., Vm+1 ]  .. 
(2.11)
 . 

 

 [I2 ⊗ Is ] 

[02 ⊗ Is ]
" ! #
I2m
= [V1 , V2 , V3 , ..., Vm+1 ] ⊗ Is (2.12)
02×2m
!
I2m
= Vm+1 (I m ⊗ Is ), où I m = . (2.13)
02×2m
Nous obtenons,
T
h
GA GA
R(Xm ) = Vm+1 (Tm Ym ⊗ Is )(I m ⊗ Is ) + (I m ⊗ Is )(Ym Tm ⊗ Is )
i
(2m+2) (2m+2) T T
+(k B k2F e1 (e1 ) ⊗ Is ) Vm+1 .
T (2m+2) (2m+2)T
hh i i
GA GA T
= Vm+1 Tm Ym I m + I m Ym Tm + k B k2F e1 e1 ⊗ Is Vm+1 .
" # " #T " # " #T
Tm GA I2m I2m GA Tm
= Vm+1  Ym + Ym
Tm+1,m 02×2m 02×2m Tm+1,m
" (2m) #" (2m) #T  
e1 e1  ⊗ Is  V T
+kB k2F m+1 .
02×1 02×1
  
(2m) (2m)T
T Y + Ym TTm + k B k2F e1 e1 T
Ym Tm+1,m
= Vm+1  m m  ⊗ Is  V T .
m+1
Tm+1,m Ym 0
(2m) (2m)T
Comme Y solution de l’équation Tm Ym + Ym TTm + k B k2F e1 e1 = 0,
" ! #
0 T
Ym Tm+1,m T
Rm (Xm ) = Vm+1 ⊗ Is Vm+1 .
Tm+1,m Ym 0

Par conséquent
" ! #
0 T
Ym Em Tm+1,m
T
kRm (Xm )kF ≤ Vm+1 × Vm+1 TY
⊗ Is
F Tm+1,m Em m 0 F
Comme {V1 , V2 , ..., Vm+1 } est une base F -orthonormée de l’espace de Krylov étendu
√
Km+1 (A, B) et Vm+1 = [V1 , V2 , ..., Vm+1 ], nous avons kVm+1 kF = m + 1. D’autre part
" ! # !
0 Ym em T
1 Tm+1,m 0 T
Ym Em Tm+1,m
Vm+1 ⊗ Is =
Tm+1,m eTm Ym 0 F
Tm+1,m (e2m T
1 ) Ym 0 F
√ T
= 2 Tm+1,m Em Ym
F
Par conséquent q
kRm (Xm )kF ≤ 2(m + 1)α,
TY
où α = Tm+1,m Em .
m
F
Pour économiser de la mémoire il est possible d’écrire Xm sous forme factorisée. Soit la
décomposition en valeurs singulières de la matrice Ym = P ΣP T , où Σ = diag(σ1 , σ2 , ..., σ2m )
et σ1 ≥ σ2 ≥ ... ≥ σ2m . Nous fixons dtol et définissons Pl la matrice constituée des k
premières colonnes de l correspondant aux l valeurs singulières supérieures ou égales à
dtol.
1/2
En posant Σl = diag(σ1 , σ2 , ..., σl ) et Zm = Vm Pl (Σl ⊗ Is ), nous obtenons l’approxi-
mation Xm = T
Zm Zm (qui est la meilleure approximation de rang l de la matrice Xm ,
on pourra se référer à [44, 56] pour un exposé complet sur ce sujet).
L’algorithme de la méthode d’Arnoldi étendu global pour résoudre l’équation de Lyapu-

nov, se résume comme suit :
Algorithme 2.1 : pour résoudre l’équation de Lyapunov (GA-Lyap)
1. On choisit une tolérance > 0 et un nombre maximal d’itérations mmax

2. Pour m = 1, 2, 3, ..., mmax
3. Construire les matrices Vm et Tm par l’algorithme d’Arnoldi étendu global appliqué
au couple (A, B).
4. Résoudre l’équation de Lyapunov projetée suivante
Tm YmGA + YmGA TTm + kBk2F e2m 2m T

1 (e1 ) = 0.
GA en utilisant (2.8).
5. Calculer la norme de Frobenius de rm
– Si kRm k < , Aller à 8
– Sinon m = m + 1

2.4 Méthode de minimisation du résidu 33
– Fin si
6. Fin pour m
7. Calculer la décomposition en valeurs singulières (SVD) de la matrice YmGA , c’est à
dire YmGA = VΣVT , où Σ = diag(σ1 , σ2 , ..., σ2mr ) et σ1 ≥ σ2 ≥ ... ≥ σ2mr ;
trouver l telle que σl+1 ≤ toltrn < σl et soit Σl = diag(σ1 , σ2 , ..., σl ) ;
1/2
Zm = Vm Vl (Σl ⊗ Is ).
T.
8. La solution approchée Xm est donnée par Xm = Zm Zm
2.4 Méthode de minimisation du résidu
Dans ce paragraphe, nous proposons une nouvelle méthode de projection pour résoudre
l’équation de Lyapunov de grande taille. En utilisant des projections sur des sous espaces
de Krylov étendus (extended) globale avec la condition de minimisation de la norme de
Frobenius du résidu associé à l’équation de Lyapunov.
Les approximations Xm de la solution exacte X de l’équation de Lyapunov sont de la

forme
Xm = Vm (Ym ⊗ Is )VTm , (2.14)
où Vm est la matrice F-orthogonale construite en appliquant m itérations de l’algorithme

d’Arnoldi étendu global à Keg (A, B).
M R qui minimisent la norme du résidu
Nous allons chercher des approximations Xm
R(Xm ) = AXm + Xm AT + BB T , où Xm = Vm (Ym ⊗ Is )VTm , c’est-à dire :
MR
Xm = arg min kAXm + Xm AT + BB T k (2.15)
Xm =Vm (Ym ⊗Is )VT
m
M R = V (Y M R ⊗ I )VT , nous avons le résultat suivant

Comme Xm m m s m
Théorème 2.4.1. Soit Vm la matrice F -orthogonale construite en appliquant m itéra-

tions de l’algorithme d’Arnoldi étendu globale à (A, B). Alors le problème de minimisa-
tion :
MR
Xm = arg min AXm + Xm AT + BB T
m
F
peut s’écrire sous la forme

" #T " # 
(2m) (2m)T

MR I2m I2m T kBk2F e1 e1 0 
Ym = arg min Tm Ym + Ym Tm + .
Ym 02×2m 02×2m 0 0
F
(2.16)

2.5 Méthodes itératives pour résoudre le problème réduit 34
Démonstration. A l’étape m, nous avons :
min AXm + Xm AT + BB T
m
F
T
= min AVm (Ym ⊗ Is ) Vm A + kBk2F V1 V1T
T T
+ Vm (Ym ⊗ Is ) Vm
Ym F
 " #T  " # !
I2m I2m T

= min Vm+1  Tm Ym ⊗ Is  ⊗ Is  + ⊗ Is Ym Tm ⊗ Is
Ym 02×2m 02×2m
(2m+2)T
i
(2m+2)
+ kBk2F e1 e1 T
⊗ Is Vm+1 .
F
 " #T " # " (2m) #" (2m) #T  
I2m I2m T e1 e1
= min Vm+1 Tm Ym + Ym Tm + kBk2F  ⊗ Is  V
Ym 02×2m 02×2m 02×1 02×1
Comme Vm+1 independent de Ym donc nous avons :

" #T " # 
(2m) (2m)T

Ym = arg min Tm Ym + Ym Tm +  .
Ym 02×2m 02×2m 0 0
F
Pour résoudre le problème de minimisation projeté (2.16), il existe plusieurs stratégies

possibles, nous proposons dans la section suivante deux méthodes itératives pour déter-
miner la solution approchée de (2.16) .
2.5 Méthodes itératives pour résoudre le problème réduit
2.5.1 LSQR global
Dans cette section, nous proposons une méthode itérative de type LSQR global (GL-
LSQR) pour résoudre le problème de minimisation projeté (2.16). La méthode LSQR
classique proposée par Paige et Saunders (1982) est une méthode itérative pour résoudre
les systèmes linéaires ou les problèmes aux moindres carrés. Pour plus de précision voir
[95]. L’algorithme LSQR est basé sur l’algorithme de Golub-Kahan [43], en utilisant
la procédure de bidiagonalisation de Lanczos dans [95]. La version globale de cette
méthode est donnée dans [106], il y a aussi la version bloc voir [75]. Dans ce qui suit,
nous allons proposer une méthode globale-LSQR (GLSQR) pour résoudre les problèmes
de minimisation de type minY kLm (Y ) − CkF , où Lm est un opérateur linéaire. Comme
notre objectif est de résoudre le problème de minimisation suivant :
" #T " # 
(2m) (2m)T

Ym = arg min Tm Ym + Ym Tm + ,
Ym 02×2m 02×2m 0 0
F

Nous proposons le problème sous la forme
min kLm (Y) − CkF , (2.17)

Y
où !
I
Lm (Y) = T̄m Y I 0 + Y T̄Tm (2.18)
0
et  
(2m) (2m)T
kBk2F e1 e1 0 
C = − .
0 0
L’opérateur adjoint de Lm par rapport au produit scalaire de Frobenius est donné par
!
T I
Lm (Y ) = T̄Tm Y + I 0 Y T̄B
m. (2.19)
0
Procédure de bidiagonalisation globale de Lanczos (Golub-Kahan) :

Cette procédure itérative est initialisée par
e1 = C
β1 U β1 = kCkF

T
α1 V1 = Lm U
e1 α1 = Lm T Ue1 .
F
et les relations de récurrence,

 ei+1 = Lm (Vei ) − αi U
βi+1 U ei

 αi+1 Vei+1 = Lm T U
ei − βi Vei
les scalaires αi > 0 et βi > 0 sont choisi de telle sorte que U

ei = Vei = 1, pour tout
F F
i = 1, 2, . . . ..
On définit les matrices suivantes
e := [U
U e1 , U
e2 , ..., U
ek ];
k
e := [Ve1 , Ve2 , ..., Ve ]

V k k
et  
α1
 
 β α2 
 2 

.. 
T̄k := 

β3 . .

 
 .. 

 . αk 

βk+1

Par construction les matrices U

ei et Vei sont F -orthonormées, ce qui signifie qu’ils sont or-
thonormées par rapport au produit scalaire de Frobenius : hUi , Uj iF = δi,j et hVi , Vj iF =
δi,j où δi,j est le symbole de Kronecker. L’utilisation de la procédure de bidiagonalisation
globale de Lanczos permet de trouver à chaque itération k des approximations Y k de la
solution exacte YmM R du problème de minimisation réduit (2.16). Nous pouvons écrire
les relations de récurrence précédentes sous forme matricielle comme suit
[Lm (Ve1 ), Lm (Ve2 ), ..., Lm (Vek )] = U

e
k+1 (T̄k ⊗ Is ). (2.20)
La méthode consiste à chercher la solution approchée Y k ayant la forme suivante
k
X
Yk = z (i) Vei .
i=1
Nous appliquons l’opérateur Lm à Y k , nous obtenons
k
!
X
k
Lm (Y ) = Lm z (i) Vei
i=1
k
X
= z (i) Lm (Vei )
i=1
= [Lm (Ve1 ), ..., Lm (Vek )] (zk ⊗ Is ) , où zk = (z (1) , z (2) , · · · , z (k) )T
= U k+1 (T̄k ⊗ Is )(zk ⊗ Is )
e
= U k+1 (T̄k zk ⊗ Is ).
e
et
C − Lm (Y k ) = e1 − U
β1 U e
k+1 (T̄k zk ⊗ Is )
F F
= U k+1 (β1 e1 ⊗ Is ) − Uk+1 (T̄k zk ⊗ Is )
e e
F
h i
= Uk+1 (β1 e1 ⊗ Is ) − (T̄k zk ⊗ Is )
F
= (β1 e1 − T̄k zk ) ⊗ Is = kβ1 e1 − Tk zk kF .
Donc le problème minimize C − Lm (Y k ) est alors équivalent à résoudre le problème

F
suivant
minimize β1 e1 − T̄k zk . (2.21)
2
Alors, la résolution du problème de minimisation (2.16), par la méthode LSQR globale

(Gl-LSQR) permet de donner une solution approchée de la solution exacte YmM R du
problème (2.16).
La méthode LSQR globale (Gl-LSQR) est résumé dans l’algorithme suivant :

Algorithme 2.2 : Algorithme LSQR globale (Gl-LSQR)
1. Initialisation : X0 = 0
β1 = kN kF , U
e1 = N/β1 ,
α1 = (Lm )T (U
e1 ) ,
F
Ve1 = (Lm )T (Ue1 )/α1 ,
Wf1 = Ve1 , Φ1 = β1 , ρ1 = α1 .
2. Pour i = 1, 2, ..., kmax :

fi = Lm (Vei ) − αi U
W ei , βi+1 = W
fi ,
U
ei+1 = W
fi /βi+1 ,
Li = (Lm )T (U
ei+1 ) − βi+1 Vei ,
αi+1 = kLi kF ,
Vei+1 q
= Li /αi+1 ,
ρi = 2 )
(ρ1 2 + βi+1
ci = ρ1 /ρ1
si = βi+1 /ρ1
θi+1 = si αi+1
ρi+1 = ci αi+1
Φi = ci Φi
Φi+1 = si Φi
Y i = Y i−1 + (Φi /ρi )W
fi
Wi+1 = Vi+1 − (θi+1 /ρi )W
fi .
Si Φi+1 est petit alors on s’arrête.
3. Fin pour.
L’algorithme précédent permet de calculer une solution approchée Y k avec 1 ≤ k ≤

M R approche
kmax, du problème de minimisation (2.20). Afin de savoir si une solution Xm
suffisamment bien la solution exacte X, nous devons être en mesure de calculer de façon
M R . Ensuite, nous donnons l’expression
simple et peu couteuse le résidu Rm associé à Xm
de kRm (Xm )k selon la norme du résidu associé à l’algorithme LSQR globale (Gl-LSQR)
appliquée au problème de minimisation réduit.
M R = V (Y M R ⊗ I )VT la solution approchée de l’équation
Théorème 2.5.1. Soit Xm m m s m
matricielle de Lyapunov obtenu après m itérations de l’algorithme d’Arnoldi étendu glo-
M R la solution exacte du problème
bal et Ym
" #T " # 
(2m) (2m)T

Ym = arg min Tm Ym + Ym Tm +  .
Ym 02×2m 02×2m 0 0
F
Alors, nous avons :

MR
√
kRm (Xm )kF ≤ m + 1Φ

où Φ =| Φkmax | est donné dans l’algorithme 1 (Gl-LSQR).
Démonstration. nous avons
MR MR MR T
kR(Xm )kF = kAXm + Xm A + BB T kF , où Xm
MR MR
= Vm (Ym ⊗ Is )VTm
h i
MR
= Vm+1 Lm (Ym ) − C ⊗ Is VTm+1
h i F
MR
≤ kVm+1 kF × − C ⊗ Is VTm+1
Lm (Ym )
F
√ MR
≤ m + 1 Lm (Ym )−C
F
√
≤ m + 1 min kLm (Y) − CkF
Y
√
≤ m + 1Φ.
2.5.2 La méthode du gradient conjugué globale préconditionné
Dans ce paragraphe, nous considérons la méthode du gradient conjugué globale précon-

ditionnés (PGCG) pour résoudre le problème de minimisation réduit (2.16). La méthode
du gradient conjugué préconditionné (PCG) classique appliquée à l’équation normale.
L’équation normale associée à (2.16) est donnée par
L∗m (Lm Y )) = L∗m (C), (2.22)
où Lm et L∗m les opérateurs donnés par (2.18) et (2.19) respectivement.
En conséquence,
min kLm (Y ) − CkF ⇐⇒ L∗m (Lm (Y ) − C) = 0.

Y
La matrice Tm sous la forme !

Tm
Tm = ,
hm
où hm représentent les 2r dernières lignes de la matrice Tm . Par conséquent, l’équation

normale (2.22) peut se écrire comme suit
T T
Tm Tm Y + Y Tm Tm + TTm Y TTm + Tm Y Tm − C1 = 0, (2.23)
où C1 = LTm (C). Considérons la décomposition en valeurs singulières (SVD) de Tm :
T
Tm = U Σ U , (2.24)
alors nous obtenons la décomposition suivante
T
Tm Tm = QDQT , (2.25)
T
où Q = U , et D = Σ Σ.
Soient Ye = QT Y Q et Ce = QT C1 Q, donc l’équation normale (2.22) peut aussi s’écrire
e m Ye T
DYe + Ye D + T e T Ye T
em + T e T − Ce = 0, (2.26)
m m
e m = QT Tm Q et Ye = QT Y Q. Cette expression propose que l’on peut utiliser l’opé-

où T
rateur matriciel suivant comme un préconditionneur
M(Ye ) = DA Ye + Ye DB . (2.27)
On peut voir que l’expression (2.22) correspondant à l’équation normale de l’opérateur

matriciel suivant !
Q
Lem (Ye ) = Tem Ye QT 0 + Ye TemT (2.28)
0
où Tem = Tm Q. Par conséquent, l’algorithme du gradient conjugué préconditionné glo-

bale est obtenu en appliquant le préconditionneur (2.27) à l’équation normale associée
à l’opérateur matriciel défini par (2.28). Ceci est résumé dans l’algorithme suivant :

Algorithme 2.3 : Algorithme du gradient conjugué globale préconditionné(PGCG)
1. On choisit une tolérance tol > 0, et un nombre maximal d’itérations jmax ;

e0 = C − L
Choisir Ye0 et on fixe R e∗ (R
em (Ye0 ) ; S0 = L e 0 ), Z0 = M−1 (S0 ) ; P0 = S0 .
m
2. Pour j = 0, 1, 2, ..., jmax

(a) Wj = Lem (Pj )
(b) αj = hSj , Zj iF /|Wj |2F
(c) Yej+1 = Yej + αj Pj
(d) R e j − αj Wj
e j+1 = R
(e) Si kR
e j+1 kF < tol. Arrêt
Sinon
(f) Sj+1 = Le∗m (R
e j+1 )
(g) Zj+1 = M−1 (Sj+1 )

(h) βj = hSj+1 , Zj+1 iF / hSj , Zj iF
(i) Pj+1 = Zj+1 + βj Pj .
3. Fin Pour.
On note que l’utilisation du préconditionneur M nécessite à chaque itération la résolu-

tion d’une équation de Sylvester. Comme la matrice D de ces équations matricielles de
Lyapunov est une matrice diagonale, ceci réduit les coûts de calcul.
M R ) à chaque itération m nécessite
Notons que le calcul de la norme du résidu R(Xm
seulement la connaissance de YmM R sans avoir à construire la solution approximative
M R . Cette norme du résidu R(X M R ) peut être calculée par
Xm m
 " #T " # 
(2m) (2m)T

MR I2m I2m T kBk2F e1 e1 0 
rm := R(Xm ) = Tm Ym + Ym Tm + .
0 0 0 0
F
(2.29)
M R n’est calculée que si la convergence est réalisée. Dans
La solution approximative Xm
la section suivante, nous proposons des méthodes (directes ou itératives) pour calculer
la solution YmM R du problème de minimisation réduit (2.15).
Pour économiser la place de la mémoire, il est possible d’écrire Xm sous forme fac-
torisée. Soit Ym = P ΣP T , la décomposition en valeurs singulières de la matrice Ym ,
où Σ ∈ R2m×2m désigne la matrice des valeurs singulières de Ym rangées dans l’ordre
décroissant, P est une matrice unitaire. Nous fixons un seuil dtol et définissons Pk la
matrice constituée des k premières colonnes de P correspondant aux k valeurs singu-
lières supérieures ou égales à dtol. Nous obtenons la décomposition en valeur singulière
tronquée
Ym ≈ Pk Σk PkT
2.6 Exemples numériques 41
1/2
où Σk = diag[σ1 , . . . , σk ]. Soit Zm = Vm Pk (Σk ⊗ Is ), nous obtenons alors l’approxi-
mation de Xm sous la forme d’une expression factorisée
T
Xm ≈ Zm Zm . (2.30)
Cette factorisation est très importante pour les problèmes de grande dimension, quand
on n’a pas besoin de calculer l’approximation de la solution Xm , mais on a besoin de la
stocker à chaque itération.
La méthode de minimisation du résidu MR pour résoudre l’équation de Lyapunov est

résumé dans l’algorithme suivant
Algorithme 2.4 : pour résoudre l’équation de Lyapunov (MR-LE)
1. On choisit une tolérance tol > 0 et un nombre maximum d’itérations itermax .

2. Pour m = 1, 2, 3, ..., itermax
3. Construire les matrices Vm et Tm par l’algorithme EGA appliqué à (A, B).
4. Résoudre le problème de minimisation
MR
Ym = arg min kLm (Ym ) + CkF .
Ym
5. Si kRm kF ≤ tol, Arrêt,

T.
6. La solution approchée Xm est donnée par Xm = Zm Zm
2.6 Exemples numériques
Dans cette section, nous présentons quelques exemples numériques de l’équation algé-
brique de Lyapunov de grande taille. Nous allons comparer l’approche de Galerkin (GA)
avec notre approche de minimisation de résidu (MR). En utilisant : la méthode LSQR
globale (Gl-LSQR) et la méthode du gradient conjugué globale préconditionné (PGCG)
pour résoudre le problème de minimisation réduit. Les algorithmes ont été codés avec
Matlab 2009. A chaque itération m de l’algorithme d’Arnoldi étendu global, l’équation
matricielle de Lyapunov projetée de taille 2m × 2m a été résolue par la méthode directe
de Bartels-Stewart [11]. Pour résoudre le problème de minimisation réduit, on s’arrête
si la norme de Frobenius du résidu est inférieure à tolp = 10−10 ou quand le nombre
maximum d’itérations kmax = 500 est atteint. Pour l’algorithme minimisation de résidu
(MR) et l’algorithme de Galerkin on s’arrête lorsque la norme de Frobenius du résidu est
inférieure à une certaine tolérance tol = 10− 12 ou bien le nombre maximum d’itérations
itermax = 50. Dans tous les tests, les coefficients de la matrice B ont été générés par
des valeurs aléatoires uniformément distribuées dans l’intervalle [0, 1]. La matrice A est
obtenue par la discrétisation par différences finies de l’opérateur suivant
∂u ∂u
Lu = ∆u − f1 (x, y) + f2 (x, y) + g(x, y), (2.31)
∂x ∂y
dans le domaine [0, 1] × [0, 1] avec les condition de Dirichlet homogène. La taille de la
matrice construite est n = n20 , où n0 est le nombre de points de la grille dans chaque
direction. La discrétisation de l’opérateur Lu donne des matrices qu’on pourra trouver
dans la bibliothèque Lyapack [89] et désignés par
A = fdm2d_matrix(n0,’f_1(x,y)’,’f_2(x,y)’,’g(x,y)’).
Les figures qui apparaissent dans la suite, nous donnent les courbes représentant la
norme de Frobenius du résidu dans une échelle logarithmique en fonction du nombre
d’itérations m.
2.6.1 Exemple 1
Dans ce premier exemple, la matrice A est donnée par :
A = f dm(n0, x2 + 2y, ey+x , 5).
Pour r = 3 et n = 6400, nous avons :
4
GA
2 MR
0
Log10 of resid. norms
−2
−4
−6
−8
−10
−12
−14
0 10 20 30 40 50
Iterations
Figure 2.1: Résultats de l’exemple 1 : MR(PCGG)

4
GA
2 MR

−2
−4
−6
−8
−10
−12
−14
0 10 20 30 40 50
Iterations
Figure 2.2: Résultats de l’exemple 1 : MR(GL-LSQR)
2.6.2 Exemple 2
Dans ce deuxième exemple, la matrice A est donnée par :
A = f dm(n0, x2 + y 2 , sin(y + x), 100).
Pour r = 2 et n = 10000, nous avons :
4
GA
2 MR
0
−2
−4
−6
−8
−10
−12
−14
0 10 20 30 40 50
Iterations
Figure 2.3: Résultats de l’exemple 2 : MR(PCGG)

2.7 Conclusion 44
4
GA
2 MR

−2
−4
−6
−8
−10
−12
−14
0 10 20 30 40 50
Iterations
Figure 2.4: Résultats de l’exemple 2 : MR(GL-LSQR)
2.7 Conclusion
Dans ce chapitre, nous avons proposé deux méthodes itératives pour résoudre les équa-
tions matricielles de Lyapunov de grande taille. Ces méthodes sont basées sur la projec-
tion sur des sous espaces de Krylov étendus global. La première est basée sur la condition
d’orthogonalité de Galerkin et la deuxième est basée sur la condition de minimisation
de résidu MR. Le problème de minimisation réduit de la deuxième méthode MR est
résolu par : la méthode LSQR global ou bien la méthode PCGG. Nous avons donné des
exemples numériques pour montrer l’efficacité des deux méthodes proposées.

Chapitre 3
L’Extended-Bloc Arnoldi
algorithme avec minimisation de
résidu pour les équations de
Sylvester creuse et de grandes
taille
Dans ce chapitre, nous présentons une nouvelle méthode pour résoudre des équations
matricielles de Sylvester de grande taille. La méthode proposée est une méthode itérative
basée sur la projection dans des sous-espaces de Krylov étendu (extended) par blocs et
la minimisation de résidu (MR minimal residual). nous résolvons le problème réduit
soit par des méthodes directes (méthode de Hu et Reichel [57]) ou bien des méthodes
itératives (globale LSQR, PGCG).
3.1 Introduction
Dans ce chapitre, nous considérons l’équation matricielle de Sylvester
AX + XB + EF T = 0 (3.1)
où A et B sont des matrices de taille n × n et s × s respectivement, E ∈ Rm×n et

F ∈ Rm×n sont des matrices de rang maximale, X ∈ Rn×s est la matrice inconnue à
déterminer.
Les équations matricielles de Sylvester jouent un rôle important dans de nombreuses

théories : automatique, traitement du signal, restauration d’images, filtrage, réduction
de modèle en théorie du contrôle, calcul d’un contrôle optimal en contrôle linéaire quadra-
tique, discrétisation d’équations aux dérivées partielles, voir par exemple les références
45
3.1 Introduction 46
suivantes [21, 39, 56, 62].
Pour résoudre l’équation matricielle de Sylvester, on distinguera deux cas :

Dans le cas où les matrices A et B sont de petite ou moyenne taille, il existe plu-
sieurs méthodes directes dont par exemple la méthode Bartels-Stewart [11] et celle de
Hessenberg-Schur [45]. Ces méthodes sont basées sur la réduction de Hessenberg de la
plus grande des deux matrices A et B et la décomposition de Schur de la plus petite.
Les dernières années, plusieurs méthodes de projection basées sur les sous-espaces de
Krylov ont été proposées, citons par exemple les références suivantes [36, 39, 56, 60, 62,
68, 82, 98, 102]. L’idée principale de ces méthodes est d’utiliser des sous-espaces de Kry-
lov [57, 98] ou Krylov étendu par blocs ou global [55, 56, 64] ou sous-espace de Krylov
rationnel [82] puis appliquer un algorithme de type Arnoldi pour construire une base
orthonormale de l’espace de projection. L’équation matricielle de Sylvester de grande
taille est alors projetée sur ces sous-espaces de Krylov. La solution de l’équation matri-
cielle peut s’exprimer en fonction de la solution de l’équation projetée et des matrices
orthogonales construites à partir les sous-espaces de Krylov. On peut citer aussi la mé-
thode ADI (Alternating Directional Implicit) [21, 23, 78]. La méthode ADI (ou Lr-ADI)
permet d’obtenir une convergence plus rapide sous certaines conditions de A et B.
L’équation matricielle de Sylvester (3.1) peut s’écrire comme un système linéaire :
(Is ⊗ A + B T ⊗ In )vec(X) = −vec(EF T ), (3.2)
où Ip désigne la matrice identité de taille p × p.
Cette écriture qui ramène l’équation de Sylvester à la résolution d’un système linéaire de
taille ns × ns, ne présente en réalité que peu d’avantage en terme de résolution, qu’elle
soit directe ou itérative. Cette formulation permet d’établir le résultat suivant :
Proposition 3.1.1. L’équation AX + XB + EF T = 0 admet une solution unique si

et seulement si les matrices A et −B n’ont aucune valeur propre commune, c’est-à-dire :
σ(A) ∩ σ(−B) = ∅.
où σ(A) et σ(−B) représentent respectivement le spectre de A et −B
Démonstration. En effet, le système linéaire (3.2) possède une unique solution si et

seulement si la matrice In ⊗ A + B T ⊗ Im est inversible. Comme les valeurs propres
de In ⊗ A + B T ⊗ Im sont (λi − µj ) où λi et µj sont les valeurs propres de A et B T
respectivement. La valeur 0 ne sera pas dans le spectre de In ⊗ A + B T ⊗ Im si et
seulement si
σ(A) ∩ σ(−B) = ∅.

3.2 Méthode de type Galerkin 47
Cette proposition donne une condition nécessaire et suffisante d’existence et d’unicité de

la solution de l’équation de Sylvester (3.1). Des hypothèses plus fortes nous permettent
de donner une écriture explicite de la solution de l’équation de Sylvester.
Définition 3.1.2. Une matrice carrée est dite stable si toutes ses valeurs propres ap-
partiennent au sous ensemble R∗− + iR de C.
Proposition 3.1.3. Si les matrices A et B sont stables, alors l’équation (3.2) possède
une unique solution X, de plus
Z +∞
X= etA EF T etB dt.
0
La suite de ce chapitre est organisé de la façon suivante : dans un premier temps, nous
allons rappeler la méthode de projection de Petrov-Galerkin sur le sous-espace de Krylov
étendu (extended). Le paragraphe suivant est consacré à notre méthode : la méthode
de la minimisation de résidu MR. Puis, nous exposons des méthodes itératives ou di-
rectes pour résoudre le problème de minimisation réduit. La Section 5 est consacrée à la
forme factorisée de l’approximation de la solution approchée (low rank approximate so-
lutions). Dans les Sections 6 et 7, nous donnons les algorithmes pour résoudre l’équation
matricielle de Sylvester par l’approche de Galerkin et la méthode de la minimisation de
résidu MR. Dans l’avant dernier paragraphe, nous appliquons notre approche MR, et la
méthode de Galerkin à l’équation matricielle de Stein non Symétrique.
Finalement, dans le dernier paragraphe, nous donnons quelques résultats numériques

qui nous permettront de comparer les méthodes de minimisation de résidu MR et la
méthode de Galerkin ou la méthode Lr-ADI.
3.2 Méthode de type Galerkin
Dans cette section nous allons rappeler l’approche de Galerkin pour résoudre l’équation
de Sylvester de grande taille (Pour plus de détails voir [56]). Cette approche est basée
sur les sous-espaces de Krylov étendus (extended) par blocs et la condition de Galerkin.
Considérons l’équation matricielle de Sylvester
AX + XB + EF T = 0, (3.3)
où A ∈ Rn×n , B ∈ Rs×s et E, F ∈ Rn×r avec r << n. Les matrices E et F étant

supposées avoir un rang r petit par rapport aux tailles de A et B. Nous allons chercher
GA de la solution X sous la forme

à construire des approximations Xm
GA
Xm = Vm YmGA WTm , (3.4)
où Vm , Wm sont les matrices orthonormales construites en appliquant simultanément

m itérations de l’algorithme d’Arnoldi étendu par blocs aux paires de matrices (A, E)
et (B T , F ) respectivement. Dans ce cas, on a
TA = VTm AVm , TB = WTm BWm .
La matrice YmGA est déterminée en imposant la condition d’orthogonalité
GA GA
Rm := AXm + Xm B − EF T ⊥ Lm (A, B),
où Lm (A, B) est le sous espace de R2mr×2mr constitué des matrices de la forme Vm Y WTm .
Donc, nous avons
VTm Rm Wm = 0.
On montre facilement que YmGA ∈ Rmr×mr est la solution de l’équation de Sylvester de

taille réduite
T
TA GA GA B e eT
m Ym + Ym Tm + E F = 0, (3.5)
e = VT E et Fe = WT F . La solution Y GA de l’équation de Sylvester de taille réduite
où E m m m
(3.5) peut être obtenue par une méthode directe, par exemple la méthode Bartels-Stewart
[11]. En pratique l’équation de Sylvester projetée (3.5) est résolue par la fonction " lyap"
de Matlab. Pour tester la convergence, il est nécessaire de calculer la norme du résidu
Rm , pour cela le résultat suivant donne une expression du résidu Rm en fonction des
matrices de petite taille.
Théorème 3.2.1. [56] Soit YmGA la solution exacte de l’équation projetée (3.5). Alors
la norme de Frobenius du résidu Rm est donnée par
q
k Rm kF = A
k Ym Em (Tm+1,m B
)T k2F + k Tm+1,m ETm Ym k2F (3.6)
L’intérêt pratique de ce résultat sur le résidu permet de tester la convergence dans

GA = V Y GA WT à chaque itération.
l’algorithme sans avoir à calculer Xm m m m
Dans les dernières années, plusieurs méthodes de projection sur les sous-espaces de Kry-
lov ont été proposées pour trouver les approximations Xm de la solution exacte X de
l’équation Lyapunov, Riccati ou Sylvester ; voir par exemple, [39, 60, 62, 68, 98, 102]. La
plupart de ces méthodes de projection utilisent la condition de Galerkin pour trouver la

solution approchée Xm .
Dans cette section, nous nous intéressons à la résolution de l’équation de Sylvester par
la méthode de projection sur les sous-espaces de Krylov étendus avec la condition de
minimisation du résidu MR (minimal residual). Donc nous allons chercher à construire
M R de la solution X sous la forme
des approximations Xm
MR
Xm = Vm YmM R WTm , (3.7)
M R minimise le résidu R (X ), où X
de telle sorte que Xm T
m m m = Vm Ym Wm . Pour cela,
nous nous intéressons à la résolution du problème de minimisation suivant :
MR
Xm = arg min AX + XB + EF T , (3.8)
X=Vm Ym WT
m
F
où Vm = [V1 , V2 , · · · , Vm ], Wm = [W1 , W2 , · · · , Wm ] ∈ Rn×mr sont les matrices orthonor-

males construites en appliquant simultanément m itérations de l’algorithme d’Arnoldi
étendu par blocs aux paires de matrices (A, E) et (B T , F ) respectivement. Dans ce cas,
nous avons :
AVm = Vm+1 T̄A
m (3.9)
et
B T Wm = Wm+1 T̄B
m. (3.10)
Nous avons alors le résultat suivant
Théorème 3.3.1. Soient Vm et Wm sont les matrices orthonormales construites en

appliquant simultanément m itérations de l’algorithme d’Arnoldi étendu par blocs aux
paires de matrices (A, E) et (B T , F ) respectivement. Alors résoudre le problème de mi-
nimisation
MR
Xm = arg min AXm + Xm B + EF T
Xm =Vm Ym WT
m
F
revient à résoudre le problème suivant

! !
I RE RFT 0
YmM R = arg min T̄A
m Ym I 0 + Ym T̄BT
m + (3.11)
Ym 0 0 0 F
où E = V1 RE et F = W1 RF sont les décomposition QR de E et F respectivement.

3.4 Méthodes pour résoudre le problème de minimisation réduit 50
Démonstration. Soient E = V1 RE et F = W1 RF sont les décomposition QR de E et F

respectivement. D’après les relations (3.9) et (3.10) nous avons
min AX + XB + EF T
X=Vm Ym WT
m
F
= min AVm Ym WTm + Vm Ym WTm B + V1 RE RFT W1T

Ym F
= min Vm+1 T̄A T

m Ym Wm + Vm Ym T̄BT T
m Wm+1 + V1 RE RFT W1T
Ym F
" ! !#
I RE RFT 0
= min Vm+1 T̄A
m Ym I 0 + Ym T̄BT
m + WTm+1
Ym 0 0 0 F
" ! !#
I RE RFT 0
= min T̄A
m Ym I 0 + Ym T̄BT
m + .
Ym 0 0 0 F
M R ) à chaque itération m nécessite

Notons que le calcul de la norme du résidu R(Xm
seulement la connaissance de YmM R sans avoir à construire la solution approximative
M R . Cette norme du résidu R(X M R ) peut être calculée par
Xm m
" ! !#
I RE RFT 0
rm := T̄A MR
m Ym I 0 + YmM R T̄BT
m + . (3.12)
0 0 0 F
M R est calculée que si la convergence est réalisée. Dans la section

La solution approchée Xm
suivante, nous proposons des méthodes (directes ou itératives) pour calculer la solution
YmM R du problème de minimisation réduit (3.11).
3.4 Méthodes pour résoudre le problème de minimisation

réduit
3.4.1 La méthode directe basée sur la formulation de Kronecker
Dans ce paragraphe, nous nous intéressons à la résolution du problème (3.11), en utili-

sant les propriétés du produit de Kronecker, le problème (3.11) peut s’écrire comme un
problème aux moindres carrés vectorielle :
" ! !#
I I
min ⊗ T̄A
m + T̄B
m ⊗ y+c , (3.13)
y 0 0 F
!
RE RFT 0
où c = vec et y = vec(Y ). Ce problème aux moindres-carrés est alors
0 0
résolu par des méthodes directes ou itératives. Par exemple :

posons " ! !#
I I
H= ⊗ T̄A
m + T̄B
m ⊗ .
0 0
Si H est de rang maximal, nous avons y = −(H T H)−1 H T c = −H † c est la solution aux
moindres carrés, où H † désigne le pseudo-inverse de la matrice H, et le résidu est donner
par rm = Hy + c = (I − H(H T H)−1 H T )c.
Pour un petit nombre d’itérations m, cette méthode donne une bonne approximation.
Toutefois, lorsque m augmente, la méthode est très lente. Dans les prochaines sous-
sections, nous exposons d’autre méthodes efficaces pour résoudre le problème de mini-
misation (3.11).
3.4.2 Méthode de Hu et Reichel
Nous exposons dans cette partie la méthode de Hu et Reichel (voir [57]) pour résoudre le
problème réduit (3.11). Plus précisément, nous appliquons l’approche donnée dans [74]
et la stratégie donnée dans [102].
Écrivons les matrices TA et TB sous la forme :

! !
TA TB
TA = et TB =
hA hB
où hA et hB représentent les 2r dernières lignes des matrices TA et TB respectivement.

Soient TA = U TA U T et TB = V TB V T la décomposition de Schur de TA et TB respecti-
vement.
Alors le problème de minimisation est équivalent à

" ! !#
I RE RFT 0
min T̄A
m Ym I 0 + Ym T̄BT
m +
Ym ∈R2mr×2mr 0 0 0 F
" # ! " #T !
TA I TB RE RFT 0
= min  Ym I 0 + Ym +  ,
Ym ∈R2mr×2mr hA 0 hB 0 0
! F
TA Y + Y TBT + RE RFT 0
= min
Ym ∈R2mr×2mr hA Y Y hTB F
ce qui est équivalent à

2 2
min TA Y + Y TTB + RE RFT + khA Y k2F + Y hTB
Ym ∈R2mr×2mr F F
2 2
= min T
U TA U Y + Y (V TB V ) + RE RF T T
+ khA Y k2F + Y hTB
Ym ∈R2mr×2mr F F
2 2 2
= min TA Ye + Ye TBT + U T RE RFT V + hA U Ye + Ye V T hTB , où Ye = U T Y V
em ∈R2mr×2mr
Y F F F

On utilisant le produit de Kronecker nous obtenons

    2
I ⊗ TA + TB ⊗ I V ect(U T RE RFT V )
   
min  I ⊗ hA U  V ect(Ye ) +  0  (3.14)
Y ∈R2mr×2mr
   
hB V ⊗ I 0 2
Nous pouvons écrire le problème (3.14) de la façon suivante

" # " # 2
R d
min y+ (3.15)
Y S 0 2
" #
I ⊗ hA U
où R = I ⊗ TA + TB ⊗ I, S = , y = V ect(Ye ) et d = V ect(U T RE RFT V ).
hB V ⊗ I
L’équation normale associée est
(RT R + S T S)y + RT d = 0. (3.16)
Si la matrice R est inversible, alors nous avons
(I + (SR−1 )T SR−1 )ye + d = 0, où y = R−1 ye.
Maintenant, en utilisant la formule de Sherman-Morrison-Woodbury, nous obtenons
ye = −d + (SR−1 )T (I + SR−1 (SR−1 )T )−1 SR−1 d.
Par conséquent, à partir de ye, nous obtenons facilement Ym la solution du problème

(3.11).
Pour calculer SR−1 on utilise la stratégie suivante :

! !
−1 I ⊗ hA U −1 (I ⊗ hA U )(I ⊗ TA + TB ⊗ I)−1
SR = (I ⊗ TA + TB ⊗ I) =
hB V ⊗ I (hB V ⊗ I)(I ⊗ TA + TB ⊗ I)−1
Posons que ! !
Z1 (I ⊗ hA U )(I ⊗ TA + TB ⊗ I)−1
=
Z2 (hB V ⊗ I)(I ⊗ TA + TB ⊗ I)−1
Alors
Z1 = (I ⊗ hA U )(I ⊗ TA + TB ⊗ I)−1 et Z2 = (hB V ⊗ I)(I ⊗ TA + TB ⊗ I)−1
Donc
(I ⊗ TA + TB ⊗ I)T Z1T = (I ⊗ (hA U )T ).

Pour calculer Z1T , il suffit de résoudre le système suivant (I ⊗ TA + TB ⊗ I)T z = (I ⊗

(hA U )T )ei pour i = 1, ..., 4mr2 . Nous remarquons que I ⊗ (hA U )T est une matrice
diagonale par blocs, avec tous les éléments diagonaux égaux à ceux de la matrice (hA U )T .
Par conséquent, chaque colonne hl := (I ⊗ (hA U )T )el , l = 1, ..., 4mr2 , satisfait à la
relation V ec((hA U )T ei eTj ) = hi+2r(j−1) , i = 1, ..., 2r, j = 1, ..., 2mr.
Donc le système
(I ⊗ TA + TB ⊗ I)T z = (I ⊗ (hA U )T )el
est équivalent à l’équation de Sylvester
TA Z + ZTB = (hA U )T ei eTj
avec z = V ect(Z), pour i = 1, ..., 2r, j = 1, ..., 2mr, avec l = i + 2r(j − 1). Donc pour
calculer Z1T il suffit de résoudre les 2mr2 équations de Sylvester suivantes :
TA Z + ZTB = (hA U )T ei eTj .
De la même manière, pour le deuxième bloc de Z2 . Donc pour calculer Z2T il suffit de
résoudre les 2mr2 équations de Sylvester suivantes :
TA Z + ZTB = (eTj hB V )ei .
La méthode de Hu et Reichel est résumé dans l’algorithme suivant :
Algorithm 3 Algorithme de Hu-Reichel pour résoudre le problème (3.14)
1. Calculer la décomposition de Schur de TA et TB , TA = U TA U T et TB = U TB U T

2. Calculer d = V ect(U T RE RFT V ) :
3. Calculer SR−1 par la résolution de 2mr équations de Sylvester
4. Résoudre ye = −d + (SR−1 )T (I + SR−1 (SR−1 )T )−1 SR−1 d.
5. Récupérer YmM R
On peut également résoudre le problème réduit (3.14) par la méthode LSQR globale ou
la méthode du gradient conjugué globale préconditionné.
Remarque 3.4.1. La méthode aux moindres carrés par Kronecker et la méthode de Hu et

Reichel nécessitent O(m3 r3 ) multiplications qui est le même ordre que pour la méthode
de Galerkin avec l’utilisation de l’algorithme Bartels-Stewart pour résoudre l’équation
de Sylvester projetée (3.2). Lorsque m augmente, l’approche des moindres carrés et la
Méthode Hu-Reichel deviennent très lentes.
Dans les deux sections suivantes, nous exposons des méthodes itératives pour résoudre
le problème de minimisation réduit.
3.4.3 La Méthode LSQR globale
A chaque étape du processus m, nous déterminons des approximations de la solution

YmM R du problème de minimisation réduit (3.11).
Le problème de minimisation réduit (3.11) s’écrit de la manière suivante
min kLm (Y ) − CkF , (3.17)

Y
où !
I
Lm (Y ) = T̄A
mY I 0 + Y T̄BT
m (3.18)
0
et !
RE RFT 0
C=− .
0 0
Nous appliquons la méthode la méthode LSQR Globale pour résoudre le problème de

minimisation (3.17) avec l’opérateur adjoint de Lm par rapport au produit scalaire de
Frobenius donné par
!
T I
Lm (Y ) = (T̄A T
m) Y + I 0 Y T̄B
m. (3.19)
0
L’algorithme global LSQR permet de calculer une solution approchée Y k avec 1 ≤

MR
k ≤ kmax, du problème de minimisation (3.17 ). Afin de savoir si une solution Xm
approche suffisamment la solution exacte X, nous devons être en mesure de calculer
M R . Ensuite, nous donnons
de façon simple et peu couteuse le résidu Rm associé à Xm
l’expression de kRm (Xm )k selon la norme du résidu associé à l’algorithme LSQR globale
(Gl-LSQR) appliquée au problème de minimisation réduit.
M R = V Y M R WT la solution approchée de l’équation matri-
Théorème 3.4.2. Soit Xm m m m
cielle de Sylvester obtenu après m itérations de l’algorithme d’Arnoldi étendu par blocs,
où " # " #
h i I RE RFT 0
mY I 0 + Y T̄BT
m + . (3.20)
Y 0 0 0 F
Alors,
MR √
kRm (Xm )kF = mΦ
où Φ =| Φkmax | est donné dans l’algorithme 3 (Gl-LSQR).

Démonstration. nous avons
MR MR MR
kRm (Xm )kF = kAXm + Xm B + EF T kF où Xm
MR
= Vm YmM R WTm
" # " #
h i I RE RFT 0
= T̄A MR
m Ym I 0 + YmM R T̄BT
m +
0 0 0 F
= k(β1 e1 − T̄k zk ) ⊗ Im kF
√
= mk(β1 e1 − T̄k zk )k2
√
= mΦ.
Le Théorème 3.4.2 joue un rôle très important dans la pratique ; il permet de calculer
la norme du résidu sans calculer la solution approchée pour tester la convergence. Des
hypothèses plus fortes nous permettent de donner une écriture de la solution de l’équa-
tion de Sylvester et une majoration en norme de l’erreur X − Xm .
Si les matrices A et B sont stables (Re(λi (A)) < 0 et Re(λi (B)) < 0), alors l’ équation
de Sylvester (1.1) a une solution unique donnée par la représentation intégrale suivante
(voir [80]) Z ∞
X= etA EF T etB dt.
0
La norme logarithmique "2-norme" de la matrice stable A est définies par
1
µ2 (A) = λmax (A + AT ) < 0.
2
La norme logarithmique donne une majoration de l’exponentielle de matrices. Il est

connu que (voir [90]).
k etA k2 ≤ eµ2 (A)t ; t ≥ 0. (3.21)
Dans le théorème suivant et avec les mêmes notations que dans le théorème précédent,
nous donnons une majoration en norme de l’erreur X − Xm
MR =
Théorème 3.4.3. Supposons que A et B soient des matrices stables et soit Xm
Vm YmM R WTm la solution approchée de l’équation de Sylvester obtenu après m itérations
de l’algorithme d’Arnoldi étendu par blocs. Alors nous avons
MR −rm
k X − Xm k2 ≤ , (3.22)
µ2 (A) + µ2 (B)
où rm est la norme du résidu donné par (3.12).
M R ) correspond au résidu
Démonstration. Rappelons que R(Xm
MR MR
AXm + Xm B + EF T = R(Xm
MR
)

de l’équation matricielle de Sylvester :
AX + XB + EF T = 0
La soustraction terme à terme permet de retrouver l’équation suivante :
MR MR MR
A(Xm − X) + (Xm − X)B = R(Xm ). (3.23)
Comme A et B sont supposés être stable, nous trouvons

Z ∞
MR
Xm −X = etA R(Xm
M R tB
)e dt.
0
Par conséquent, en utilisant la norme logarithmique, on obtient

Z ∞
MR MR
kX− Xm k≤k R(Xm ) k2 et(µ2 (A)+µ2 (B)) dt.
0
M R ) k ≤k R(X M R ) k et (µ (A) + µ (B)) < 0, nous obtenons :

Comme k R(Xm 2 m F 2 2
MR −rm
k X − Xm k2 ≤ ,
µ2 (A) + µ2 (B)
M R ) k ) est la norme du résidu donné par (3.12).

où rm (rm =k R(Xm F
La convergence de l’algorithme LSQR globale (Gl-LSQR) peut être lente. Ensuite nous
devons utiliser cette méthode avec un préconditionneur. Pour cela, il existe plusieurs
stratégies que nous pouvons utiliser, on pourrait adapter par exemple les techniques
utilisées dans [7].
3.4.4 La méthode du gradient conjugué globale préconditionné
Dans cette section, nous considérons la méthode du gradient conjugué globale précondi-
tionnés (PGCG) pour résoudre le problème de minimisation réduit (3.11). La méthode
du gradient conjugué préconditionné (PCG) classique appliquée à l’équation normale.
L’équation normale associé à (3.11) est donnée par
L∗m (Lm Y )) = L∗m (C), (3.24)
où Lm et L∗m les opérateur donné par (3.18) et (3.19) respectivement.

A B
Soient les matrices Tm et Tm sous la forme
! !
A TA
m TB
m
Tm = et TB = ,
hA
m hB
m
A B
où hA et hB représentent les 2r dernières lignes de la matrice Tm et Tm , respectivement.
Par conséquent, l’équation normale (3.24) peut s’écrire comme suit
AT A BT B T T
Tm Tm Y + Y Tm Tm + TA B A B
m Y Tm + Tm Y Tm − C1 = 0, (3.25)
A
où C1 = LTm (C). Considérons la décomposition en valeurs singulières (SVD) de Tm et
B
Tm :
A T B T
Tm = U A ΣA V A ; Tm = U B ΣB V B , (3.26)
alors nous obtenons la décomposition suivante
AT A BT B
Tm Tm = QA DA QTA , Tm Tm = QB DB QTB , (3.27)
T
où QA = V A , QB = V B et DA = ΣA ΣA .
Soient Ye = QTA Y QB et Ce = QTA C1 QB , donc l’équation normale (3.25) peut aussi s’écrire
e A Ye T
DA Ye + Ye DB + T e B + (T
e A )T Ye T
e B )T − Ce = 0, (3.28)
m m m m
e A = QT TA QA , T
où T e B = QT TB QB et Ye = QT Y QB . Cette expression propose que
m A m m B m A
l’on peut utiliser l’opérateur matriciel suivant comme un préconditionneur
P(Ye ) = DA Ye + Ye DB . (3.29)
On peut voir que l’expression (3.28) correspondant à l’équation normale de l’opérateur

matriciel suivant
!
QA
Lem (Ye ) = TemA Ye QTB 0 + Ye (TemB )T (3.30)
0
A B
où TemA = Tm QA et TemB = Tm QB . Par conséquent, l’algorithme du gradient conjugué
préconditionné globale est obtenu en appliquant le préconditionneur (3.29) à l’équation
normale associée à l’opérateur matriciel défini par (3.30). Ceci est résumé dans l’algo-
rithme suivant :

3.5 Forme factorisée de l’approximation de la solution 58
Algorithme 3 : Algorithme du gradient conjugué globale préconditionné(PGCG)
1. On choisit une tolérance tol > 0, et un nombre maximal d’itérations jmax ;

e0 = C − L
Choisir Ye0 et on fixe R e∗ (R
em (Ye0 ) ; S0 = L e 0 ), Z0 = P −1 (S0 ) ; P0 = S0 .
m
2. Pour j = 0, 1, 2, ..., jmax

(a) Wj = Lem (Pj )
(b) αj = hSj , Zj iF /|Wj |2F
(c) Yej+1 = Yej + αj Pj
(d) R e j − αj Wj
e j+1 = R
(e) Si kR
e j+1 kF < tol. Arrêt
Sinon
(f) Sj+1 = Le∗m (R
e j+1 )
(g) Zj+1 = P −1 (Sj+1 )

(h) βj = hSj+1 , Zj+1 iF / hSj , Zj iF
(i) Pj+1 = Zj+1 + βj Pj .
3. Fin Pour.
On note que l’utilisation du préconditionneur P nécessite à chaque itération la résolution

d’une équation de Sylvester. Comme les matrices DA et DB de ces équations matricielles
de Sylvester sont des matrices diagonales, ceci réduit les coûts de calcul.
3.5 Forme factorisée de l’approximation de la solution
On pourra décomposer la solution Xm comme produit de deux matrices de rang inférieur,

c’est-à-dire Xm = Z1 Z2T où Z1 et Z2 sont des matrices de petite rang (inférieure à
2mr). Pour cela, nous considérons la décomposition en valeurs singulières de la matrice
Ym ∈ R2mr×2mr .
Ym = U1 Σ U2T
où Σ est la matrice diagonale des valeurs singulières de Ym rangées dans l’ordre décrois-
sant, U1 et U2 étant deux matrices unitaires. Nous fixons une certaine tolérance dtol et
définissons U1,l et U2,l les matrices constituées des l premières colonnes de U1 et de U2
correspondants aux l valeurs singulières supérieures ou égales à dtol. Nous obtenons la
décomposition en valeur singulière tronquée
Ym ≈ U1,l Σl U2,l T

3.6 L’algorithme GA-LRSE 59
1/2 1/2
où Σl = diag[σ1 , . . . , σl ]. Soit Z1,m = Vm U1,l Σl , et Z2,m = Wm U2,l Σl , nous obte-
nons alors l’approximation de Xm sous la forme d’une expression factorisée
T
Xm ≈ Z1,m Z2,m . (3.31)
3.6 L’algorithme GA-LRSE
Dans cette section, nous présentons l’algorithme pour résoudre l’équation de Sylvester de
grande taille (GA-LRSE) par la méthode de projection de Galerkin sur les sous-espaces
de Krylov étendus (extended) par blocs.
Algorithme 3 : pour résoudre l’équation de Sylvester (GA-Lr-SE)

2. Pour m = 1, 2, 3, ..., mmax
3. Construire les matrices Vm et TA
m par l’algorithme 1 appliqué au couple (A, E).
4. Construire les matrices Wm et TB T

m par l’algorithme 1 appliqué au couple (B , F ).
5. Résoudre l’équation de Sylvester projetée suivante
TA GA GA B T e eT
m Ym + Ym (Tm ) + E F = 0
6. Calculer la norme de Frobenius du résidu :

q
k Rm kF = B
k YmGA Em (Tm+1,m A
)T k2F + k Tm+1,m ETm YmGA k2F

– Sinon m = m + 1
– Fin si
7. Fin pour m
dire YmGA = VΣW T , où Σ = diag(σ1 , σ2 , ..., σ2mr ) et σ1 ≥ σ2 ≥ ... ≥ σ2mr ;
A =V VΣ 1/2
B 1/2
Zm m l l , etZm = Wm Wl Σl
A ZB . T
9. La solution approchée Xm est donnée par Xm ≈ Zm m

3.7 L’algorithme MR-LR-Sylvester 60
3.7 L’algorithme MR-LR-Sylvester
Dans cette section, nous présentons l’algorithme pour résoudre l’équation de Sylvester
de grande taille (MR-Lr-SE) par la méthode de minimisation du résidu sur les sous-
espaces de Krylov étendus (extended) par blocs.
Algorithme 4 : pour résoudre l’équation de Sylvester (MR-Lr-SE)
1. On choisit une tolérance tol > 0 et un nombre maximum d’itérations itermax .

2. Pour m = 1, 2, 3, ..., itermax
A
3. Construire les matrices Vm et Tm par l’algorithme (EBA) appliqué au couple
(A, E).
B
4. Construire les matrices Wm et Tm par l’algorithme (EBA) appliqué au couple
(B T , F ).
5. Résoudre le problème de minimisation
" # " #
h i I T RE RFT 0
m Ym I 0 + Ym T̄B
m + .
Ym 0 0 0 F
6. Si kRm kF ≤ tol, Arrêt,

T .
7. La solution Xm est donnée par Xm ≈ Z1,m Z2,m
3.8 Résolution de l’équation de Stein non symétrique
Dans ce paragraphe, nous considérons l’équation de Stein non symétrique
AXB − X + EF T = 0, (3.32)
où A et B sont des matrices carrées de taille n × n et s × s, respectivement. Les matrices

E et F sont de dimension n × r et r × s respectivement. Les matrices A et B sont
supposées creuses de très grande taille r << n, s.
L’équation matricielle (3.32) (également appelée équation généralisée de Stein ou équa-

tion de Sylvester discrète) joue un rôle important dans la théorie du contrôle, la restau-
ration d’images, filtrage de systèmes dynamiques non linéaire à temps discret de grande
taille, à chaque étape de la méthode de Newton pour résoudre les équations algébriques
de Riccati discrète et d’autres problèmes. Pour plus d’informations, voir les références
suivantes [17, 18, 32, 33, 42, 66, 80].

3.8 Résolution de l’équation de Stein non symétrique 61
Les méthodes directes pour résoudre l’équation matricielle de Stein non symétrique sont
proposés dans [8, 11, 45, 98]. Ces méthodes sont intéressantes lorsque les matrices A
et B sont de petite taille. L’équation matricielle (3.32) peut être formulée comme un
système linéaire de taille ns × ns :
(B T ⊗ A − Is ⊗ In )vec(X) = −vec(EF T ), (3.33)
On peut utiliser les méthodes de type Krylov pour résoudre le système linéaire (3.33),
par exemple l’algorithme de GMRES [99]. Lorsque les matrices A et B sont de grande
taille et creuses, cette approche ne peuvent être appliquées directement. Dans ce cas, il
existe des méthodes pour résoudre l’équation matricielle de Stein non symétrique, par
exemple la méthode low-rank Stein block-Arnoldi. Cette méthode basée sur la projection
de l’équation (3.32) dans les sous-espaces de Krylov par blocs [63, 66].
L’équation matricielle de Stein non symétrique admet solution unique si et seulement si

λi (A)λj (B) 6= 1 pour tout i = 1, 2, ..., n et j = 1, 2, ..., s, où λi (A) est la i-ème valeur
propre de la matrice A. En particulier, si ρ(A) < 1 et ρ(B) < 1, où ρ(A) désigne le rayon
spectral de la matrice A, alors l’équation (3.33) admet une solution unique X. De plus,
elle peut s’exprimer sous la forme d’une série matricielle suivante :
∞
X
X= Ai E F T B i .
i=0
Nous nous référons à [80] pour les démonstrations de ces résultats. Nous supposons que
la conditions d’existence et d’unicité sont vérifiées pour chacune des équations de Stein
rencontrées dans la suite de cette section. La méthode proposée par Jbilou (low-rank
Stein block-Arnoldi) permet de construire des approximations Xm de petit rang de la
forme
T
Xm = Vm Ym Wm
où Vm , Wm sont les matrices orthonormales construites en appliquant simultanément m

itérations de l’algorithme d’Arnoldi par blocs aux paires de matrices (A, E) et (B T , F )
respectivement, et Ym la solution de l’équation de Stein projetée suivante :
HA Ym HBT − Ym + EeFe T = 0, (3.34)
T E et F
où Ee = Vm e = W T F.
m
Dans cette section, nous exposons ensuite : la même méthode de Jbilou mais dans les
sous-espaces de Krylov étendus (extended) par blocs, et la méthode de la minimisation
du résidu (MR) pour résoudre l’équation de Stein non symétrique.

3.8.1 Méthode de Galerkin
Dans cette section nous allons rappeler l’approche de Galerkin pour résoudre l’équation
de Stein non symétrique de grande taille. Cette approche a été introduite par Jbilou
dans [63]. On considère l’équation de Stein non symétrique (3.32). Nous supposons que
les matrices E et F sont de rang maximal, nous allons chercher à construire une suite
GA )
des approximations de petit rang (Xm m∈N∗ de la solution X sous la forme
GA
Xm = Vm YmGA WTm , (3.35)
où Vm , Wm sont les matrices orthogonales construites en appliquant simultanément m

itérations de l’algorithme d’Arnoldi étendu par blocs aux paires de matrices (A, E) et
(B T , F ) respectivement. Dans ce cas, on a
TA = VTm AVm , et TB = WTm BWm .
Pour trouver la matrice YmGA ∈ R2mr×2mr nous utilisons la condition d’orthogonalité

suivante
RGA GA GA T
m := AXm B − Xm + EF ⊥F Lm (A, B),
où Lm (A, B) est le sous espace des matrices de la forme Vm Ym WTm . La condition d’or-
thogonalité ci-dessus permet de donner l’équation suivante :
VTm RGA
m Wm = 0.
Comme les matrices Vm et Wm sont orthogonales, donc il est plus facile de monter que
YmGA ∈ R2mr×2mr est la solution de l’équation de Stein de taille réduite
TA YmGA TB − YmGA + E
e Fe T = 0, (3.36)
e = VT E et Fe = WT F .
où E m m
En supposant que λi (TA ).λj (TB ) 6= 1 pour toutes i = 1, 2, ..., 2mr et j = 1, 2, ..., 2ms,
pour assurer l’existence et l’unicité de la solution YmGA . La résolution de l’équation de
Stein de taille réduite (3.36) peut être obtenue par une méthode directe, par exemple la
méthode Hessenberg-Schur [33]. Pour tester la convergence, Il est nécessaire de calculer
la norme de Frobenius du résidu RGA
m , pour cela le résultat suivant donne une expression
du résidu RGA
m en fonction des matrices de petite taille.
Théorème 3.8.1. [63] Soient YmGA la solution exacte de l’équation de Stein réduite
GA = V Y GA WT la solution approchée de l’équation (3.32). Alors la norme
(3.36) et Xm m m m
de Frobenius du résidu Rm est donnée par
q
kRG GA
m (Xm )kF =
2 + β2 + γ2 ,
αm m m (3.37)

où les réels αm , βm et γm sont définis par
αm = TA B
m Ym Em (Tm+1,m )
T A
, βm = Tm+1,m ETm (TB
m)
T
,
F F
et
A
γm = Tm+1,m ETm Ym Em (Tm+1,m
B
)T .
F

RGA
m = AVm Ym WTm B − Vm Ym WTm + EF T
F
= Vm+1 T̄A BT T T T T
m Ym T̄m Wm+1 − Vm Ym Wm + V1 RE RF W1
F
!
I
= Vm+1 T̄A BT T
m Ym T̄m Wm+1 − Vm+1 Ym I 0 WTm+1 + V1 RE RFT W1T
0 F
" ! !#
I RE RFT 0
= Vm+1 T̄A BT
m Ym T̄m − Ym I 0 + WTm+1
0 0 0 F
Comme les matrices Vm+1 et Wm+1 sont orthogonales, et

" # " #
A TA
m B TB
m
Tm = A
, et Tm =
Tm+1,m ETm B
Tm+1,m ETm
Donc, nous avons

" # " #T ! !
TA
m TB
m Ym 0 RE RFT 0
RGA
m = 
A
Ym − + 
Tm+1,m ETm B
Tm+1,m ETm 0 0 0 0
F
" ! !#
TA B T
m Ym (Tm ) TA B
m Ym Em (Tm+1,m )
T −Ym + RE RFT 0
= A
+
Tm+1,m ETm (TB
m)
T A
Tm+1,m ETm Ym Em (Tm+1,m
B )T 0 0 F
!
TA B T
m Ym (Tm ) − Ym + RE RF
T TA B
m Ym Em (Tm+1,m )
T
= A
Tm+1,m ETm (TB
m)
T A
B )T F
Puisque YmGA la solution exacte de l’équation de Stein (3.36) donc nous avons
!
0 TA B
m Ym Em (Tm+1,m )
T
RGA
m = A
.
Tm+1,m ETm (TB
m)
T A
B )T F
Par conséquent
2 2 2 2
RGA
m = TA B
m Ym Em (Tm+1,m )
T A
+ Tm+1,m ETm (TB
m)
T A
+ Tm+1,m ETm Ym Em (Tm+1,m
B
)T .
F F F F
h i
Nous avons finalement le résultat (3.37) et nous rappelons que ETm = 02r×2(m−1)r I2r .

Ce résultat est très important puisqu’il nous permet de calculer la norme de Frobenius
GA ) sans avoir a effectuer le produit des matrices de grande taille.
du résidu Rm (Xm
La méthode ainsi obtenue sera notée GA-N-Stein, et algorithme correspondant peut être
formulé comme suit :
Algorithme 3 : pour résoudre l’équation de Stein non symétrique (GA-Lr-SE)

2. Pour m = 1, 2, 3, ..., mmax
3. Construire les matrices Vm et TA
m par l’algorithme 2 (Ch 1) appliqué au couple
(A, E).
4. Construire les matrices Wm et TB
m par l’algorithme 2 (Ch 1) appliqué au couple
(B T , F ).
5. Résoudre l’équation de Stein projetée suivante
TA GA B T GA e eT
m Ym (Tm ) − Ym + E F = 0
6. Calculer la norme de Frobenius du résidu :

q
kRG GA
m (Xm )kF =
2 + β2 + γ2 ,
αm m m
où
αm = TA B
m Ym Em (Tm+1,m )
T A
, βm = Tm+1,m ETm (TB
m)
T
,
F F
et
A
γm = Tm+1,m ETm Ym Em (Tm+1,m
B
)T .
F

– Sinon m = m + 1
– Fin si
7. Fin pour m
dire YmGA = VΣW T , où Σ = diag(σ1 , σ2 , ..., σ2mr ) et σ1 ≥ σ2 ≥ ... ≥ σ2mr ;
A =V VΣ B1/2 1/2
Zm m l l , et Zm = Wm Wl Σl
A ZB . T
9. La solution approchée Xm est donnée par Xm ≈ Zm m
L’équation de Stein projetée est résolu par les méthode directes.

3.8.2 Méthode de minimisation du résidu
Dans ce paragraphe, nous appliquons la méthode de minimisation de résidu (MR) à

l’équation de Stein non symétrique. A chaque itération m de l’algorithme d’Arnoldi
M R de l’équation de Stein non symétrique
étendu par blocs, on a une approximation Xm
M R la solution du problème de minimisation suivant :
(3.32), où Xm
X M R = arg min AXB − X + EF T , (3.38)

X=Vm Ym WT
m
F
où Vm , Wm ∈ Rn×mr sont les matrices orthogonales construites en appliquant simulta-

nément m itérations de l’algorithme d’Arnoldi étendu par blocs aux paires de matrices
(A, E) et (B T , F ) respectivement. Donc, nous avons les relations suivantes :
AVm = Vm+1 T̄A

m (3.39)
et
B T Wm = Wm+1 T̄B
m. (3.40)
Nous avons le résultat suivant
Théorème 3.8.2. Le problème de minimisation
MR
Xm = arg min AXm B − Xm + EF T
Xm =Vm Ym WT
m
F
est équivalent au problème de minimisation réduit suivant

" " # " ##
I h i RE RFT 0
YmM R = arg min T̄A B T
m Ym (T̄m ) − Ym I 0 + . (3.41)
Ym 0 0 0 F
où E = V1 RE et F = W1 RF sont les décomposition QR de E et F respectivement.

min AXB − X + EF T
X=Vm Ym WT
m
F
= min AVm Ym WTm B − Vm Ym WTm + V1 RE RFT W1T

Ym F
!
I
= min Vm+1 T̄A BT T
m Ym T̄m Wm+1 − Vm+1 Ym I 0 WTm+1 + V1 RE RFT W1T
Ym 0 F
" ! !#
I RE RFT 0
= min Vm+1 T̄A BT
m Ym T̄m − Ym I 0 + WTm+1
Ym 0 0 0 F
" ! !#
I RE RFT 0
= min T̄A BT
m Ym T̄m − Ym I 0 + .
Ym 0 0 0 F

Le problème de minimisation réduit (3.41) est résolu par les méthodes directes ou itéra-
tives comme le cas du problème de minimisation associé à l’équation de Sylvester.
On peut énoncer l’algorithme MR pour résoudre l’équation de Stein non symétrique

comme suit
Algorithm 4 Algorithme de MR pour résoudre l’équation de Stein non symétrique

• On choisit une tolérance tol > 0 et un nombre maximum d’itérations itermax .
• Pour m = 1, 2, 3, ..., itermax
A
• Construire les Vm , Tm , exécuter l’algorithme EBA pour (A, E).
B
• Construire les Wm , Tm , exécuter l’algorithme EBA pour (B T F ).
• Résoudre le problème de minimisation
! !
I RE RFT 0
YmM R = arg min T̄A B T
m Ym (T̄m ) − Ym I 0 + .
Ym 0 0 0
F
• Si kRm kF ≤ tol, Arrêt,

T .
• La solution Xm est donnée par Xm ≈ Z1,m Z2,m
Dans cette section, nous présentons quelques exemples numériques de l’équation matri-
cielle de Sylvester et l’équation matricielle de Stein non Symétrique. Nous allons com-
parer l’approche de Galerkin (GA) avec notre approche : minimisation de résidu (MR),
pour des problèmes de grande taille. En utilisant : la méthode LSQR globale (Gl-LSQR),
la méthode du gradient conjugué globale préconditionné (PGCG), direct (le produit de
Kronecker), et la méthode Hu et Reichel (HR), pour la résolution le problème de mi-
nimisation réduit. Les algorithmes ont été codés avec Matlab 2009. A chaque itération
m de l’algorithme d’Arnoldi étendu par blocs, équation matricielle de Sylvester projetée
de taille 2mr × 2mr a été résolue par la méthode Bartels-Stewart [11] pour l’approche
de Galerkin (GA). Le problème de minimisation réduit est résolue par LSQR global
(Gl-LSQR) et la méthode du gradient conjugué globale préconditionné (PGCG). On
s’arrête si la norme de Frobenius du résidu est inférieure à toll = 10−12 ou quand le
nombre maximum d’itérations kmax = 1000 est atteint. Pour l’algorithme minimisation
de résidu (MR) et l’algorithme de Galerkin s’arrête lorsque la norme de Frobenius du
résidu est inférieure à une certaine tolérance tol (qui sera précisée pour chaque figure) ou
bien le nombre maximum d’itérations itermax = 50. Dans tous les tests, les coefficients

des matrices E et F ont été générés par des valeurs aléatoires uniformément distribuées
dans l’intervalle [0, 1]. Les matrices A et B sont obtenues par à partir de la discrétisation
par différences finies de l’opérateur suivant
∂u ∂u
Lu = ∆u − f1 (x, y) + f2 (x, y) + g(x, y), (3.42)
∂x ∂y
matrice construite est n = n20 , où n0 est le nombre de points de la grille dans chaque
direction. La discrétisation de l’opérateur Lu donne des matrices dans la bibliothèque
Lyapack [89] et désignés par
On utilisera aussi des matrices de la collection Matrix-Market et aussi la matrice ’flow_meter’

de la collection Oberwolfach (imtek.de/simulation/benchmark).
Les figures qui apparaissent dans la suite, nous donnons les courbes représentant la
norme de Frobenius du résidu dans une échelle logarithmique en fonction du nombre
d’itérations m. La courbe de la méthode de minimisation de résidu (MR) est représentée
par ligne pointillée (verte) et celle de Galerkin (GA) par linge continue (en bleu). On
donne également des tableaux pour comparer les méthodes proposées. Pour cela, on a
fixé la tolérance à tol = 10−7 pour les deux méthodes, et nous donnons : la norme de
Frobenius du résidu, le nombre d’itérations ainsi que le temps CPU en secondes.
3.9.1 Tests numérique pour l’équation de Sylvester
Dans ce paragraphe, nous présentons quelque exemples numériques de l’équation de

Sylvester de grande talle.
Exemple 1.1 Dans ce premier exemple, les matrices A et B sont obtenus à partir de la
discrétisation de l’opérateur Lu et LB respectivement
∂u ∂u
LA = ∆u − exy − sin(xy) − (y 2 − x2 ) u,
∂x ∂y
∂u ∂u q 2
LB = ∆u − 100 ex − 10 xy − x + y 2 u,
∂x ∂y
sur le carré unité [0, 1] × [0, 1] avec des conditions aux bord de Dirichlet homogènes. Les
matrices A et B sont de taille n = n20 et s = s20 respectivement, où n0 et s0 sont les
nombres de points de la grille intérieure dans la direction x et y respectivement.
Pour ce exemple, nous avons choisi n = 4900, s = 3600 et r = 2. on obtient la figure 1.1
ci-dessus.

8
GA
MR
6

2
−2
−4
−6
−8
0 10 20 30 40 50 60
Iterations
Figure 3.1: Résultat pour l’exemple 1.1.
Notons que le problème de minimisation réduit (associé à la méthode de minimisation

du résidu MR) est résolu par la méthode itérative LSQR globale (Gl-LSQR).
Exemple 1.2 Dans cet exemple, nous avons utilisé les matrices A = pde2961 et B =
thermal de la collection Harwell-Boeing (http ://math.nist.gov/MatrixMarket). Les ma-
trices A et B sont de taille n = 2961 et s = 3456 respectivement. Pour ce test, nous
avons choisi r = 2 et le problème de minimisation réduit a été résolu par la méthode
direct des moindres carrés basé sur le produit de Kronecker (dans la section 3.4.1).
8
GA
MR
6
2
−2
−4
−6
−8
−10
0 5 10 15 20 25 30 35 40
Iterations
Exemple 1.3 Dans cet exemple, nous avons utilisé les matrices A et B sont obtenus à
partir de la discrétisation de l’opérateur Lu. avec des dimensions n = 8100 et s = 3456
respectivement. Dans la figure 3.3, nous utilisons les matrices

A =fdm2d_matrix(90,’cos(x + y)’,’sin(y 2 )’,’100’),
et B = mmread(0 thermal.mtx0 ); r = 5 et tol = 10−10 . Le problème de minimisation de

taille réduit est résolu par la méthode de Hu et Reichel.
6
GA
4 MR
2
−2
−4
−6
−8
−10
−12
0 5 10 15 20 25 30 35
Iterations
Exemple 1.4 Dans La figure 3.4. les matrices A et B sont obtenus à partir de la
discrétisation de l’opérateur Lu et LB respectivement, avec des dimensions n = 4900
et s = 3600, respectivement. Dans la figure 3.4, nous utilisons les matrices A = add32 ;
n = 4960 et B =fdm2d_matrix(p0 ,0 x + y 0 ,0 x × y 0 ,0 1000 ); r = 4 et tol = 10−10 . Le
problème de minimisation projeté résolu par la méthode Gl-LSQR.
4
GA
MR
2
0
−2
−4
−6
−8
−10
−12
0 5 10 15 20 25 30
Iterations

Exemple 1.5. Dans cet exemple, nous considérons l’équation particulier de Sylvester
suivante :
AX + XA = EF T (3.43)
Cette équation intervienne dans le calcul du grammien croisé d’un système dynamique
linéaire invariant par rapport au temps. Pour plus de détails sur cette équation et la
connexion avec les problèmes de réduction de modèle [5]. On choisi n = s = 6400, r = 3
et tol = 10−9 . Le problème de minimisation projeté résolu par la méthode PGCG.
8
GA
6 MR
4
−2
−4
−6
−8
−10
0 10 20 30 40 50
Iterations
Les exemples numériques confirment l’efficacité de la méthode proposée par rapport à

l’approche de Galerkin.
Exemple 1.6. Dans cet exemple, nous avons comparé les performances de la méthode
de minimisation du résidu MR et l’approche Galerkin. Le problème de minimisation
réduit a été résolu par la méthode des moindres carrés direct, par LSQR globale (Gl-
LSQR), par le gradient conjugué globale préconditionné (GPCG) ou par la méthode de
Hu-Reichel (HR). Dans le tableau 3.1 suivant, nous donnons les résultats obtenus par
différentes méthodes de MR et la méthode de Galerkin.
Les matrices A et B sont issues de la collection Matrix-Market collection et aussi la
matrice ’flow_meter’ de la collection Oberwolfach (imtek.de/simulation/benchmark),
ou bien provenant de la discrétisation de l’opérateur Lu (3.42). Pour ce test aussi, nous
avons choisi r = 2.
D’après les résultats du tableau 3.1 que les méthode MR(GPCG) et MR(direct) sont
plus rapides que les autre méthodes, et on observe aussi que la méthode MR(GPCG)
donne les meilleurs résultats en terme de temps de calcul.

Les matrices Méthode Temps N.d’itér k Rm kF
A = thermal, MR(PGCG) 26s 10 1.1 × 10−8

B = add32 MR (direct) 26s 9 2.1 × 10−8
n = 3456, s = 4960 MR (Gl-LSQR) 28s 10 1.3 × 10−8
MR(HR) 27s 9 3.5 × 10−8
GA 25s 9 1.4 × 10−8
2
A = f dm(cos(xy), ey x , 100) MR(PGCG) 45s 15 1.9 × 10−8
B = add32 MR (direct) 56s 15 2.3 × 10−8
n = 90000, s = 4960 MR (Gl-LSQR) 49s 17 1.1 × 10−8
MR(HR) 88s 15 3.1 × 10−8
GA 83s 21 8.6 × 10−5
A = f dm(sin(xy), exy , 10) MR(PGCG) 56s 18 2.5 × 10−8

B = thermal MR(direct) 58s 16 3.1 × 10−8
n = 122500, s = 3456 MR(Gl-LSQR) −− 50 − − −−
MR(HR) 110s 16 2.0 × 10−8
GA 82s 25 8.4 × 10−5
A = f low MR (PGCG) 49s 16 2.0 × 10−8

B = f dm(sin(xy), xy, 1000) MR (direct) 700s 35 3.1 × 10−8
n = 9669, s = 40000 MR(Gl-LSQR) −− 50 − − −−
MR(HR) 150s 35 2.4 × 10−8
GA 95s 50 2.5 × 10−5
A = f dm(xy, y 2 , 1) MR(PGCG) 706s 18 2.1 × 10−8

B = f dm(xy, cos(xy), 10) MR(direct) 1500s 42 1.5 × 10−8
n = 122500, s = 48400 MR(Gl-LSQR) −− −− − − −−
MR(HR) −− −− − − −−
GA 805s 50 4.2 × 10−4
Table 3.1: Résultats pour l’exemple 1.6
Maintenant, on va comparer la méthode MR(GPCG) avec la méthode ADI (Lr_ADI)

[21, 23].
Exemple 1.7
Dans cet exemple, nous comparons notre méthode MR(GPCG) avec deux autres mé-
thodes qui font référence dans la littérature : l’approche de GA et la méthode de rang
inférieur ADI (Lr_ADI)[21, 23]. Pour cette expérience, les matrices A et B sont les
mêmes que celles données dans [Exemple 1, [23]]. Elles sont obtenues à partir la dis-
crétisation de l’opérateur Lu (donnée par la relation (3.42)) à 5 points sur un maillage
régulier. Les matrices A et B sont de taille n = 6400 et s = 3600 respectivement. Pour
cette expérience, nous avons choisi r = 4 et la tolérance est tol = 10−11 .

Méthode Res. Norm. Temps

MR(PGCG) 1.2 × 10−12 2.8s
GA 6.4 × 10−12 3.2s
Lr_ADI 2.5 × 10−12 5.2s
Table 3.2: Résultats de l’exemple 1.7.
3.9.2 Tests numérique pour l’équation de Stein non symétrique
Dans ce paragraphe, nous présentons quelque exemples numériques de l’équation de

Stein non symétrique de grande talle. Le problème de minimisation réduit associer à la
méthode MR sera résolu soit par la méthode direct de Kronecker (direct) ou bien par
la méthode du gradient conjugué globale préconditionné (PGCG). L’équation de Stein
projetée (3.41) à été résolue en utilisant la fonction "dlyap" de Matlab. Les matrices A
et B sont obtenus à partir de la discrétisation de l’opérateur Lu (3.42), ou bien de la
collection "Matrix-Market". Les matrices E et F sont des valeurs aléatoires uniformément
distribuées dans l’intervalle [0, 1].
Exemple 2.1 Dans cet exemple, nous avons choisi
A =fdm2d_matrix(90,’cos(x + y)’,’sin(y 2 )’,’100’),
et B =0 thermal.mtx0 , donc les matrices A et B sont de taille n = 8100 et s = 3456

respectivement. On choisit aussi r = 3 et tol = 10−9 .
On obtient, pour les deux méthodes : MR(direct) et GA, les graphes suivants :
6
GA
MR
4
2
−2
−4
−6
−8
−10
0 5 10 15 20 25 30 35 40
Iterations
Figure 3.6: Résultats de l’exemple 2.1

Exemple 2.2 On prend dans cet exemple, A = f low et B = thermal, donc les matrices
A et B sont de taille n = 9669 et s = 3456 respectivement. On prend aussi r = 3 et
tol = 10−9 , et on obtient, pour la méthode MR(PGCG) et l’approche de Galrkin (GA),
les courbes suivantes :
8
GA
MR
6
4
−2
−4
−6
−8
0 10 20 30 40 50
Iterations
Exemple 2.3 Dans cet exemple, on augmente la taille de l’équation de Stein non Sy-
métrique n = 10000 et s = 8100. Nous considérons A = fdm2d_matrix(100,’cos(x +
y)’,’sin(y 2 )’,’100’), B = fdm2d_matrix(90,’x + y’,’xy’,’200’). Pour r = 4 et tol = 10−7 ,
nous obtenons les courbes suivantes :
4
GA
MR
2
−2
−4
−6
−8
0 5 10 15 20 25 30 35 40 45
Iterations

3.10 Conclusion 74
Exemple 2.4. Dans cet exemple, on compare les méthodes trois méthodes MR(PGCG),
MR(direct) et GA pour différentes matrices A et B. Pour les tests dans le tableau (3.3),
nous avons choisi tol = 10−7 et r = 2. Dans le tableau (5.3), nous donnons les résultats
obtenus par les trois méthodes.
Les matrices Méthode Temps N.d’itér k Rm kF
A = thermal, MR(PGCG) 55s 4 5.05 × 10−9

B = add32 MR (direct) 58s 6 1.22 × 10−9
n = 3456, s = 4960 GA 57s 7 3.53 × 10−8
A = f dm(cos(x + y), siny 2 , 100) MR(PGCG) 5s 7 9.32 × 10−8

B = f dm(x + y, yx, 200) MR (direct) 4s 3 3.03 × 10−8
n = 10000, s = 6400 GA 18s 42 4.55 × 10−8
A = f dm(sin(xy), exy , 10) MR(PGCG) 5s 3 7.29 × 10−8

B = pde2961 MR(direct) 128s 18 9.10 × 10−8
n = 14400, s = 3456 GA 31s 50 9.81 × 10−7
A = f low MR (PGCG) 14s 14 6.65 × 10−8

B = f dm(cos(x + y), sin(y 2 ), 1000) MR (direct) 178s 20 4.06 × 10−8
n = 9669, s = 12100 GA 44s 50 4.78 × 10−5
A = f dm(exy , sin(xy), y 2 − x2 ) MR(PGCG) 8s 9 4.59 × 10−8

B = f dm(ex , cos(xy), y 2 ) MR(direct) 230s 26 1.78 × 10−8
n = s = 10000 GA 60s 50 2.24 × 10−2
Table 3.3: Résultats de l’exemple 2.4
Nous remarquons que la méthode MR(PGCG) permet de donner une bonne approxima-
tion, avec de meilleurs temps de calcul.
3.10 Conclusion
Dans ce chapitre, nous avons présenté une nouvelle méthode itérative pour résoudre
l’équation de Sylvester ou Stein non symétrique de grande taille. La méthode proposée
est basée sur la projection dans des sous-espaces de Krylov étendu (extended) par blocs
avec la condition de minimisation du résidu, elle se ramène à la résolution d’un problème
de minimisation réduit. Ce problème de minimisation réduit est résolu par des méthodes
itératives ou directes. Pour accélérer la convergence de la résolution du problème de
minimisation d’ordre inférieur, nous utilisons un préconditionnement pour la méthode de
gradient conjugué global associé à l’équation normale. les approximations de la solution
exacte sont données sous la forme factorisée : produit de deux matrices de rang inférieur

3.10 Conclusion 75
(low rank approximation), qui permet d’économiser de la mémoire pour les problèmes
de grande dimension. L’avantage de la méthode de minimisation de résidu par rapport
à l’approche de Galerkin est la stabilité numérique et également les meilleurs temps de
calcul. Cela est démontré par les tests numériques effectués.

Chapitre 4
Méthode de minimisation du résidu

pour la résolution de l’équation de
Riccati continue de grande taille
4.1 Introduction
Les équations de Riccati algébriques jouent un rôle fondamental dans de nombreux

problèmes de la théorie du contrôle linéaire : contrôle linéaire avec coût quadratique,
contrôle de type H∞ ou H2 , réduction de modèle, des équations différentielles etc, pour
plus de détails voir [2, 3, 20, 33, 55, 64, 65, 80]. Le calcul de la solution de l’équation
algébrique de Riccati fréquemment utilisée dans différents domaines comme par exemple
les problèmes de calcul de contrôle linéaire quadratique (QLQ) suivant
Z +∞
min y(t)T y(t) + u(t)T u(t) dt (QLQ) (4.1)
0
où y(t) ∈ Rn et le contrôle u(t) ∈ Rn vérifient l’équation différentielle suivante :
(
x0 (t) = Ay(t) + Bu(t)
(4.2)
y(t) = Cx(t); x(0) = x0,
avec A ∈ Rn×n , B ∈ Rn×s , C ∈ Rr×n , s << n et r << n.
Nous supposons que les matrices B et C sont de rang maximal s et r respectivement.

Sous les conditions énoncées dans [109] : si le couple (A, B) est c-stabilisable (c’est à
dire qu’il existe une matrice K telle que A − BK soit stable) et (C, A) est c-détectable
(c’est à dire que (AT , C T ) est c-stabilisable) ; alors le problème (4.2) est minimisée par
u(t) = −B T Xx(t) (4.3)
76
4.1 Introduction 77
où X est la solution symétrique semi-définie positive de l’équation de Riccati matricielle

continue
AT X + XA − XBB T X + C T C = 0. (4.4)
Considérons la matrice Hamiltonienne H de dimension 2n × 2n associée à l’équation

(4.4)
" #
A BB T
H=
CT C −AT
La solution de l’équation de Riccati est liée au calcul de sous espaces invariants de
la matrice Hamiltonienne H associée. Si les matrices Y, Z, W dans Rn×n , avec Z non
singulière, satisfont : ! !
Y Y
H = W,
Z Z
alors X = Y Z −1 est l’unique solution de l’équation algébrique de Riccati (4.4).

De nombreuses méthodes numériques ont été proposées pour résoudre l’équation de Ric-
cati (4.4), on peut citer dans un premier temps, et pour les problèmes de petite taille
l’approche par les valeurs propres consiste à calculer les sous-espaces invariants de la ma-
trice Hamiltonienne H [16, 80, 91], comme les méthodes du type Lanczos symplectiques
[16, 80, 107], ou la méthode de Newton [6, 14, 15, 80]. Pour les problèmes de grande
taille, on a les méthodes de projection sur le sous-espaces de Krylov ou sur le sous-espace
de Krylov par blocs [62, 65], ou bien sur le sous-espace de Krylov étendu (extended) par
blocs [55]. Par exemple dans l’article de Heyouni et Jbilou [55], on cherche à construire
GA de la solution X de l’équation de Riccati (4.4) sous la forme
une solution approchée Xm
GA
Xm = Vm YmGA VTm .
dont la condition de Galerkin est définie comme suit :
VTm (AT XmGA + Xm

GA GA
A − Xm BB T Xm
GA
+ C T C)Vm = 0
M R de la solution X de
Notre approche consiste à chercher la solution approchée Xm
l’équation de Riccati (4.4) sous la forme
MR
Xm = Vm YmM R VTm
mais avec la condition de minimisation suivante
X M R = arg min AX + XAT − XBB T X + CC T (4.5)

X=Vm Ym VT
m
F
où Vm est une matrice orthonormée.

Nous proposons dans ce travail une nouvelle méthode itérative basée sur la projection
sur un sous-espace de Krylov étendu (extended) par blocs et la minimisation de résidu
MR (Minimal Residual). Nous montrons, à travers des exemples numériques en fin de

ce chapitre, que la méthode proposée donne de bons résultats numériques comparés à la

méthode de Galerkin [55]. Le reste du chapitre est organisé comme suit : Dans la Section
3, nous allons rappeler la méthode de Galerkin pour résoudre l’équation de Riccati de
grande taille, cette méthode à été proposée par Heyouni et Jbilou [55]. Puis, la Section
4 introduit notre approche. La section 6 est consacrée à la méthode de MINRES globale
pour résoudre une équation matricielle linéaire de type Lyapunov généralisée de petite
taille. La Section 7, concernant la forme factorisée de la solution approchée, Puis, la
Section 8 et 9 l’algorithme générale pour résoudre l’équation de Riccati symétrique de
grande taille par l’approche de Galerkin et l’approche de minimisation de résidu, et on
termine ce chapitre par les exemples numériques et une petite conclusion.
4.2 Méthode de Galerkin
Dans cette section, on va rappeler la méthode de projection sur le sous-espace de Krylov

étendu par blocs avec la condition d’orthogonalité de Galerkin pour résoudre l’équation
de Riccati de grande taille. Cette méthode est proposée par Heyouni et Jbilou en 2009
[55].
On considère l’équation de Riccati (4.4)
AT X + XA − XBB T X + C T C = 0
GA )
On cherche à construire une suite d’approximations de petit rang (Xm m∈N∗ de la
solution X sous la forme
GA
Xm = Vm YmGA VTm (4.6)
où Vm est la matrice orthonormale construite en appliquant l’algorithme d’Arnoldi

étendu par blocs à la paire (AT , C T ). La matrice YmGA ∈ R2mr×2mr est déterminée
en imposant la condition d’orthogonalité de Galerkin
Rm (Xm ) := AT Xm
GA GA
+ Xm GA
A − Xm BB T Xm
GA
+ C T C ⊥F Lm (AT , C T ).
où Lm (AT , C T ) est le sous espace vectoriel de Rn×n constitués des matrices de la forme
X = Vm Y VTm .
c’est à dire
VTm Rm Vm = 0.

Nous pouvons facilement montrer que YmGA ∈ R2mr×2mr est la solution de l’équation de
Riccati de taille réduite
TA YmGA + YmGA TTA − YmGA B

eBe T Y GA + C
m
eCe T = 0, (4.7)
où TA = VTm AT Vm , Ce = VTm C T et B
e = VT Bm . La solution Y GA de l’équation de
m m
Riccati de taille réduite peut être obtenue par une méthode directe comme celles qui
sont décrites dans [16, 27, 88]. Pour tester la convergence, il est nécessaire de calculer la
norme du résidu Rm , pour cela le résultat suivant donne une expression de la norme de
Frobenius du résidu Rm en fonction des matrices de petite taille.
Théorème 4.2.1. Soit YmGA la solution exacte de l’équation réduite (4.7) obtenue par
l’algorithme d’Arnoldi étendu par blocs. On a alors la norme de Frobenius du résidu
Rm := Rm (Xm ) est donnée par
√ A
k Rm kF = 2 k Tm+1,m ETm YmGA kF (4.8)
Démonstration. La démonstration est donnée dans ([64]).
L’intérêt pratique de ce résultat sur le résidu est de nous permettre de mettre en place
un test d’arrêt sur l’algorithme sans avoir a calculer le produit Vm YmGA Vm
T à chaque
itération.
GA k, où X est
Nous donnons ci-dessous une majoration de la norme de l’erreur kX − Xm
la solution exacte de l’équation de Riccati (4.4).
Théorème 4.2.2. Soient Ỹm la matrice de taille 2s × 2ms correspondante aux 2s der-
nières lignes de la matrice de YmGA , γm = kTm+1,m Ỹm k, η = kBB T k, Am = A −
BB T Xm
GA et supposons que σ T 2
m = sep(Am , −Am ) > 0. Alors si 4γm η/δm < 1, nous
avons
GA 2γm
kX − Xm k≤ p
2 − 4γ η
. (4.9)
δm + δm m
Démonstration. La preuve est donnée dans ([64]).
La plupart des méthodes proposées ces dernières années pour résoudre l’équation de
Riccati continue de grande taille sont des méthodes de projection sur les sous-espaces de
Krylov, par exemple : méthode de projection sur le sous espaces de Krylov par blocs [62]
ou sur sous espaces de Krylov étendu par blocs [55]. Pour d’autres méthodes on peut
citer [14, 22, 60, 64, 66, 98, 99]. Ces méthodes génèrent des approximations de petit rang
de la solution de l’équation (4.4). Ces méthodes qui utilisent la condition de Galerkin,

présentent l’avantage de converger rapidement. Nous proposons une nouvelle méthode

de type projection sur le sous espace de Krylov étendu avec la condition de minimisation
du résidu à la place de la condition de Galerkin. L’objectif de cette méthode MR est de
e (AT , C T ).
minimiser la norme du résidu sur Km
Comme les approximations de la solution exacte sont de la forme :
Xm = Vm Ym VTm , (4.10)
M R à l’étape m, liée Ke (AT , C T ). Cette

on va chercher à construire une approximation Xm m
M R est de la forme
approximation Xm
X M R = Vm YmM R VTm ,
M R est la solution du problème de minimisation suivant

et Xm
MR
Xm = arg min AT X + XA − XBB T X + C T C , (4.11)
X∈Lm (AT ,C T ) F
où Lm (AT , C T ) est l’ensemble des matrices de la forme X = Vm Y VTm , avec Vm =

[V1 , V2 , · · · , Vm ] est la matrice orthogonale construite en appliquant l’algorithme d’Ar-
noldi étendu par blocs à (AT , C T ) et m un entier naturel. Dans ce cas nous avons la
relation suivante
AVm = Vm+1 T̄A

m
Le théorème suivant permet de transformer le problème de minimisation grande taille

(4.11) en un problème de minimisation de petite taille.
Théorème 4.3.1. Soit Vm une matrice orthonormale construite en appliquant l’algo-

rithme d’Arnoldi étendu par blocs, aux paire (AT , C T ). Alors le problème de minimisation
MR
Xm = arg min AT Xm + Xm A − XBB T X + C T C
Xm =Vm Ym VT
m
F
est équivalent au problème de minimisation réduit suivant

" # " # " #
h i I I h i T
RC RC 0
YmM R = arg min T̄m Ym I 0 + Ym T̄Tm − T
Ym B̃ B̃ Ym I 0 +
Ym 0 0 0 0 F
(4.12)
où V1 RC est la décomposition QR de C T et B̃ = VT B.
Démonstration. Soient V1 RC la décomposition QR de C T , et Vm la matrice orthogo-

nale construite en appliquant l’algorithme d’Arnoldi étendu par blocs à (AT , C T ). On a

4.4 Résolution du problème de minimisation réduit 81
VTm AT = Vm+1 T̄A

m et la matrice Vm+1 est orthogonale. Ce qui nous donne,
min AT X + XA − XBB T X + CC T
X=Vm Ym VT
m
F
= min AT Vm Ym VTm + Vm Ym VTm A − Vm Ym B̃ B̃ T Vm Ym VTm + V1 RC RC

T
W1T
Ym F

I I T
RC RC 0
= min Vm+1 T̄m Ym [ I 0 ]+ 0
Ym T̄Tm − 0
T
Ym B̃ B̃ Ym [ I 0 ]+ 0 0
VTm+1
Ym F
" " # " # " ##
I I h i T
RC RC 0
= min T̄m Ym I 0 + Ym T̄Tm − Ym B̃ B̃ T Ym I 0 + .
Ym 0 0 0 0 F
où B̃ = VT B.
4.4 Résolution du problème de minimisation réduit
A chaque étape m de l’algorithme d’Arnoldi étendu par blocs, nous déterminons une
approximation de la solution YmM R du problème (4.11) que nous posons sous la forme,
min kF (Y )kF (4.13)

Y
où F est l’opérateur non linéaire
F (Y ) = AY B + B T Y AT + B T Y EY B + F.
Plus précisément
" #
h i T
RC RC 0
A = Tm , E = −B
eBeT , B = I 0 et F = .
0 0
Le problème de minimisation (4.13) est équivalent au problème suivant
min f (Y ), (4.14)
Y
où
f (Y ) := min kF (Y )k2F =< F (Y ), F (Y ) >F .
Y
Il est claire que f est différentiable, et la différentielle de f en Y est :
f (Y + tH) − f (Y )
dfY (H) = lim , ∀ H ∈ Rn×n . (4.15)
t−→0 t
On a
F (Y ) = AY B + B T Y AT + B T Y EY B + F

Donc pour tout réel t on a
F (Y +tH) = AY B+tAHB+B T Y AT +tB T HAT +B T Y E(Y +tH)B+tB T HE(Y +tH)B+F
= AY B+B T Y AT +F+t(AHB+B T HAT )+B T Y EY B+tB T Y EHB+tB T HEY B+t2 B T HEHB
= F (Y ) + t(AHB + B T HAT + B T Y EHB + B T HEY B) + t2 (B T HEHB).
Pour simplifier le calcul on pose
Φ(H) = AHB + B T HAT + B T Y EHB + B T HEY B.
Donc
f (Y + tH) = < F (Y + tH), F (Y + tH) >F

= < F (Y ) + tΦ(H) + t2 B T HEHB, F (Y ) + tΦ(H) + t2 B T HEHB >F
= < F (Y ), F (Y ) >F +2t < F (Y ), Φ(H) >F +t2 (2 < F (Y ), B T HEHB >F
+ < Φ(H), Φ(H) >F ) + 2t3 (< B T HEHB, Φ(H) >F ) + t4 k B T HEHB kF
Alors
dfY (H) = 2 < F (Y ), Φ(H) >F

= 2 < F (Y ), AHB + B T HAT + B T Y EHB + B T HEY B >F
= 2 < AT F (Y )B T + BF (Y )A + E T Y T BF (Y )B T + BF (Y )B T Y T E T , H >F
Comme Y est symétrique, nous avons
dfY (H) = 2 < AT F (Y )B T + BF (Y )A + E T Y BF (Y )B T + BF (Y )B T Y E T , H >F
La condition nécessaire de minimisation est que dfY = 0 ce qui permet de donner

l’équation suivante :
∀H : 2 < AT F (Y )B T + BF (Y )A + E T Y BF (Y )B T + BF (Y )B T Y E T , H >F = 0.
Ce qui est équivalent à :
Ψ(Y ) := AT F (Y )B T + BF (Y )A + E T Y BF (Y )B T + BF (Y )B T Y E T = 0 (4.16)
Nous allons résoudre numériquement l’équation matricielle non linéaire (4.16) par la
méthode de Newton, donc la dérivée de Fréchet de l’application non linéaire Ψ(Y ) en Y

est donnée par :
Ψ0Y (Z) = AT FY0 (Z)B T + BFY0 (Z)A + E T ZBF (Y )B T + E T Y BFY0 (Z)B T

+ BFY0 (Z)B T Y E T + BF (Y )B T ZE T
= (AT + E T Y B)FY0 (Z)B T + BFY0 (Z)(A + B T Y E T ) + E T ZBF (Y )B T
+BF (Y )B T ZE T .
où FY0 (Z) = B T Z(AT + EY B) + (A + B T Y E)ZB.

On remplace FY0 (Z) dans Ψ0Y (Z) ce qui donne
Ψ0Y (Z) = (AT + E T Y B)B T Z(AT + EY B)B T + (AT + E T Y B)(A + B T Y E)ZBB T

+ BB T Z(AT + EY B)(A + B T Y E T ) + B(A + B T Y E)ZB(A + B T Y E T )
+E T ZBF (Y )B T + BF (Y )B T ZE T
Comme BB T = I et la matrice E symétrique donc nous avons
Ψ0Y (Z) = (AT B T + E T Y )Z(AT B T + E T Y ) + (AT B T + E T Y )T Z(AT B T + E T Y )T

+(AT + E T Y B)(A + B T Y E)Z + Z(AT + EY B)(A + B T Y E T )
+EZ(BAY + Y AT B T + Y EY + BF B T ) + (BAY + Y AT B T + Y EY + BF B T )ZE.
Pour simplifier
M(Z) := Ψ0Y (Z) = A1 ZA1 + AT1 ZAT1 + A2 Z + ZA2 + EZA3 + A3 ZE (4.17)
où A1 = AT B T + EY, A2 = (AT + EY B)(A + B T Y E),
A3 = BAY + Y AT B T + Y EY + BF B T .
On va construire une suite de matrices Yk d’approximations de la solution exacte Y de

l’équation Ψ(Y ) = 0. On choisit un premier terme Y0 = 0 et, en supposant Yk connu, on
définit Yk+1 comme suit
Yk+1 = Yk + Nk
où la matrice Nk est solution de l’équation linéaire matricielle
Ψ0Yk (Nk ) = −Ψ(Yk ) (4.18)
Remarque 4.4.1. L’opérateur linéaire M défini par (4.17) est symétrique par rapport
au produit scalaire de Frobenius, c’est-à-dire :
∀, V et W, on a hM(V ), W iF = hV, M(W )iF .

4.5 Méthode de MINRES globale. 84
Pour résoudre l’équation matricielle linéaire symétrique (4.18), nous proposons la mé-
thode MINRES globale.
4.5 Méthode de MINRES globale.
La méthode MINRES (MINimum RESidual method) de Paige et Saunders [ [94], SIAM

J. Numer. Anal., 1975], est une méthode itérative très souvent utilisée pour résoudre des
systèmes linéaires creux Ax = b dans le cas A est symétrique, en utilisant le processus
Lanczos et la propriété de minimisation de la norme du résidu.
Dans cette section, nous proposons la méthode de MINRES globale (Gl-MINRES) pour
résoudre les équations matricielles linéaires symétriques comme :
M(Z) = C (4.19)
où M : Rn×n −→ Rn×p l’opérateur matriciel linéaire symétrique et C est une matrice de

taille n × p.
Soit V une matrice réelle de taille n × p . Soit k un entier naturel, on définit l’opérateur
Mk par la relation de récurrence suivante
M0 = I et Mk (V ) = M(Mk−1 (V )).
On définit le sous-espace de Krylov d’ordre k, associé à M et V par :
Kk (M, V ) = sev{V, M(V ), ..., Mk−1 (V )}
engendré par les matrices V, M(V ), ..., Mk−1 (V ), où k est un entier naturel non nul.
L’ensemble Kk (M, V ) ainsi défini est un sous-espace vectoriel de Rn×n .
L’algorithme de Lanczos global [70], dont les étapes sont données ci-dessous, produit
une base orthonormale {V1 , ..., Vk } de l’espace de Krylov Kk (M, V ) au sens de la norme
de Frobenius, c’est à dire telle que hVi , Vj iF = δi,j , 1 ≤ i, j ≤ k.

Algorithm 5 L’algorithme de Lanczos symétrique global.

1. Initialisation : V1 = V /kV kF
2. Calculer Ṽ1 = M(V1 )
3. α1 = hV1 , Ṽ1 iF
4. Ṽ1 = Ṽ1 − α1 V1
5. β1 = kṼ1 kF
6. V2 = Ṽ1 /β1
7. Itérations :Pour j = 2 : k faire
8. Ṽj = M(Vj ) − βj Vj−1
9. αj = hVj , Ṽj iF
10. Ṽj = Ṽj − αj Vj
11. βj+1 = kṼj kF
12. Vj+1 = Ṽj /βj+1
13. Fin pour j
On définit la matrice Vk = [V1 , ..., Vk ] ∈ Rn×pk et T̃k , la matrice tridiagonale de taille (k+
1) × k dont les entrées non nulles sont calculées par l’algorithme de Lanczos symétrique
global.
 
α1 β1 0 ... 0
.. .. 
 
.

 β1 α2 β2 . 
 
 .. .. 
0 β2 . . 0 
T̃k =  . (4.20)
 
 . .. .. .. 
 . . . . βk−1 

 ..
 
.. .. 
 .

. . αk  
... ... 0 βk
La matrice Tk de taille k × k est obtenue en supprimant la dernière ligne de T̃k . Par

construction, la matrice par blocs Vk est F-orthonormale, ce qui signifie que les matrices
V1 , ..., Vk constituent un système orthonormal au sens du produit scalaire de Frobenius,
c’est à dire hVi , Vj iF = δi,j , 1 ≤ i, j ≤ k. On établit les identités suivantes de manière
immédiate
Proposition 4.5.1. [19] Après k itérations de l’algorithme Lanczos symétrique global,

nous avons les relations suivantes :
[M(V1 ), ..., M(Vk )] = Vk (Tk ⊗ Ip ) + Ek+1 , (4.21)

où
Ek+1 = hk+1,k [0n×p , ..., 0n×p , Vk+1 ],
et
[M(V1 ), ..., M(Vk )] = Vk+1 (T̃k ⊗ Ip ). (4.22)
La méthode MINRES globale consiste à construire itérativement une suite (Zk )k∈N∗
d’approximations de la solution Z ∗ de l’équation (4.19) de la façon suivante
On commence par choisir un premier terme X0 ∈ Rn×n et en notant R0 le résidu

correspondant : R0 = C − M(Z0 ), on construit les itérés Zk de telle manière à avoir
Zk = Z0 + Z̃k où Z̃k ∈ Kk (M, R0 ) (4.23)
Rk = C − M(Zk )⊥F Kk (M, M(R0 )). (4.24)
Par construction, le résidu Rk = C − M(Zk ) est le projeté F -orthogonal de R0 sur

le sous-espace Kk (M, M(R0 )) = sev{M(R0 ), ..., Mk (R0 )} engendré par les matrices
M(R0 ), ..., Mk (R0 ). Par conséquent, la matrice xk est solution du problème de minimi-
sation
min kC − M(Z)kF (4.25)

Z∈Z0 +KK (M,R0 )
Le résultat suivant [19] établit que le problème de moindres carrés (4.25) est équivalent
à un problème de dimension réduite
Proposition 4.5.2. L’approximation Zk construite par la méthode MINRES globale est

donnée par
Zk = Z0 + Vk (zk ⊗ In ),
où zk est solution du problème aux moindres carrés
min kkR0 kF e1 − T̃k zk2 (4.26)

z∈Rk
où e1 est le premier vecteur de la base canonique de Rk+1 .
Pour résoudre le problème (4.26), nous considérons la décomposition QR de la matrice

T̃k
T̃k : R̃k = Qk T˜k ,
où R˜k ∈ R(k+1)×k est triangulaire supérieure et Qk ∈ R(k+1)×(k+1) est orthogonale.

Posons gk = kR0 kF Qk e1 . En notant R1 la matrice k × k obtenue en éliminant la dernière
ligne de R̃k , le vecteur yk est calculé en résolvant le système triangulaire R1 yk = gk .

A chaque itération, le résidu Rk doit être calculé, ce qui peut s’avérer coûteux. La
proposition suivante permet de calculer kRk kF sans évaluer M(Zk ).
Proposition 4.5.3 ([19]). À l’étape K, le résidu Rk = C − M(Zk ) obtenu par la

méthode MINRES globale vérifie les deux identités suivantes
Rk = γk+1 Vk+1 (QT ek+1 ⊗ Ip ) (4.27)
et
kRk kF = |γk+1 |, (4.28)
où γk+1 est la dernière composante du vecteur gk = kR0 kF Qk e1 et ek+1 est le dernier

vecteur de la base canonique de Rk+1 : ek+1 = (0, . . . , 0, 1)T .
Dans l’algorithme MINRES globale ci-dessous. En utilisant la technique de redémar-

rage après un nombre choisi d’itérations. A chaque redémarrage, on choisit la dernière
approximation calculée comme terme initial pour l’algorithme de MINRES global.
On peut maintenant écrire l’algorithme de la méthode MINRES global avec redémar-

rage :
Algorithm 6 Algorithme MINRES global avec redémarrage

1. Initialisation : on choisit Z0 , une tolérance et on pose iter = 0
2. Calculer R0 = C − M(Z0 ), β = kR0 k, et V1 = R0 /β
3. Calculer Vk et Tek par l’algorithme de Lanczos global appliqué à (M, V1 ).
4. Calculer zk réalisant min kkR0 kF e1 − T̃k zk2
z∈Rk
5. Calculer Zk = Z0 + Vk (zk ⊗ Ip )
6. Calculer la norme du résidu kRk kF en utilisant la proposition (4.28)
7. Si kRk kF <
8. Arrêt
9. Sinon
10. Z0 = Zk , R0 = Rk , β = kR0 kF , V1 = R0 /β, iter=iter+1, Aller à 2 :
11. Fin si
la généralisation au cas M non symétrique est la méthode GMRES global ([19]).

4.6 Forme factorisée de la solution approchée 88
4.6 Forme factorisée de la solution approchée
Dans cette section, on va approximer la solution approchée Xm par le produit de deux

matrices de rang inférieurs (low-rank approximation), pour économiser la mémoire dans
les test numériques de grande taille.
À travers les exemples numériques on peut montrer que la solution exacte de l’équation
de Riccati continue est de petite rang par une approximation, comme dans l’exemple
suivant. Dans cet exemple, on considère le système dynamique
(
ẋ(t) = A x(t) + B u(t),
y(t) = C x(t),
où la matrice A résulte de la discrétisation par différences finies de l’équation aux dérivées

partielles de l’opérateur
∂u ∂u
L(u) = ∆u − exy − sin(xy) − (y 2 − x2 ) u,
∂x ∂y
dans le domaine Ω = [0 1] × [0 1] avec des conditions aux bord de Dirichlet homogènes.
La taille de la matrice A est donnée par n = n20 , où n0 est le nombre de points de la grille
dans chaque direction. Pour ce test, nous prenons n0 = 30 et r = 4. Les coefficients des
matrices B ∈ Rn×r et C ∈ Rr×r sont des valeurs aléatoires uniformément distribuées
dans l’intervalle [0, 1].
−2
Singular values of the solution
−4
−6
−8
−10
−12
−14
−16
−18
−20
0 100 200 300 400 500 600 700 800 900
Figure 4.1: Valeurs singulières de la solution exacte.
La figure 4.1 représente les valeurs singulières de la solution exacte dans le cas où n = 900
et r = 4. On constate que les valeurs singulières décroissent rapidement vers zéro. Cela
implique que le rang numérique de la solution exacte est petit. Dans ce cas particulier,

4.7 Algorithme GA-CAREs 89
on a rang(X) ' 25. Donc la solution approchée Xm peut être écrite comme produit de
deux matrices de petit rang.
Nous considérons la décomposition en valeurs singuliers de la matrice Ym ∈ R2mr×2mr .
Ym = U Σ U T
où Σ est la matrice diagonale des valeurs singulières de Ym rangées dans l’ordre décrois-
sant, U est une matrice unitaire. Nous fixons une certaine tolérance dtol et définissons Ul
la matrice constituée des l premières colonnes de U correspondant aux l valeurs singu-
lières supérieures ou égales à dtol. Nous obtenons la décomposition en valeur singulière
tronquée
Ym ≈ Ul Σl Ul T
1/2
où Σl = diag[σ1 , . . . , σl ]. Soit Zm = Vm Ul Σl , nous obtenons alors l’approximation de
Xm sous la forme d’une expression factorisée
T
Xm ≈ Zm Zm . (4.29)
4.7 Algorithme GA-CAREs
Dans ce paragraphe, nous rappelons l’algorithme d’Arnoldi étendu par blocs en utilisant
l’approche de Galerkin (GA-CARE) pour résoudre l’équation de Riccati continue de
grande taille, donc l’algorithme GA-CARE est comme suit :
4.8 Algorithme MR-CAREs
la méthode de minimisation de résidu (MR-CAREs) pour la résolution de l’équation de

Riccati continue de grande taille est résumé dans l’algorithme suivant

Algorithm 7 Algorithme GA-CAREs

1. On choisit une tolérance > 0, un nombre maximum d’itérations mmax et une
tolérance toltrunc .
2. Pour m = 1, 2, 3, ..., mmax faire
3. Appliquer l’algorithme 1 à (AT , C T ) pour construire les matrices Vm et Tm .
4. Résoudre l’équation de Riccati
T
Tm Ym + Ym TTm − Ym B̃m B̃m T
Ym + C̃m C̃m = 0;
où B̃m = Vm B et C̃m = Vm C.
5. Calculer la norme de Frobenius du résidu
√
rm = 2kTm+1,m Ỹm kF ,
où Ỹm est la matrice de taille 2s × 2ms correspondant à les 2s dernières lignes de

la matrice de YmGA .
6. Si rm < , Aller à 10.
7. sinon m = m + 1.
8. Fin Si
9. Fin Pour
10. Calculer la décomposition en valeurs et vecteurs singuliers SVD de Ym , c’est-à-dire
Ym = U Σ U T , où Σ = diag[σ1 , σ2 , · · · , σ2ms ] et σ1 ≥ σ2 ≥ · · · ≥ σ2ms .
Déterminer l telle que σl ≥ dtol > σl+1 , soit Σl = diag[σ1 , σ2 , · · · , σl ].
1/2
Calculer Zm = Vm Ul Σl
T.
11. L’approximation de Xm donnée par Xm ≈ Zm Zm
Dans cette section, nous présentons quelques exemples numériques de l’équation de Ric-
cati continue de grande taille. Les algorithmes ont été codés en Matlab R2009a. Les
équations de Riccati de taille réduite (4.7) ont été résolues en utilisant la fonction "care"
de Matlab. Le nombre d’itérations de l’algorithme de Newton a été limité à itermax = 40
et le nombre d’itérations de l’algorithme MINRES globale a été limité à max = 1000
avec le nombre d’itérations de redémarrage fixé à kredém = 12. La résolution d’équations
de type Lyapunov généralisée ont été effectuées jusqu’à ce que la norme du résidu soit
inférieure à tolM IN RES = 10−7 . Le critère d’arrêt pour l’algorithme de Newton est défini
par le paramètre
∆m =k Ym+1 − Ym k / k Ym k< tolN ewton = 10−9 .
Les résultats suivants donnent une comparaison en résidu entre notre méthode MR-
CAREs et la méthode GA-CARES [55]. Les courbes dans la suite donnent la norme du

Algorithm 8 Algorithme MR-CAREs

1. On choisit une tolérance > 0, un nombre maximum d’itérations mmax et une
tolérance toltrunc .
2. Pour m = 1, 2, 3, ..., mmax faire
3. Appliquer l’algorithme 1 à (AT , C T ) pour construire les matrices Vm et Tm .
4. Résoudre l’équation (4.16)
Ψ(YmM R ) = 0
5. Calculer la norme de Frobenius du résidu rm = kF (YmM R )kF .

6. Si rm < , Aller à 10.
7. sinon m = m + 1.
8. Fin Si
9. Fin Pour
10. Calculer la décomposition en valeurs et vecteurs singuliers SVD de Ym , c’est-à-dire
Ym = U Σ U T , où Σ = diag[σ1 , σ2 , · · · , σ2ms ] et σ1 ≥ σ2 ≥ · · · ≥ σ2ms .
Déterminer l telle que σl ≥ dtol > σl+1 , soit Σl = diag[σ1 , σ2 , · · · , σl ].
1/2
Calculer Zm = Vm Ul Σl
T.
11. L’approximation de Xm donnée par Xm ≈ Zm Zm
résidu k Rm kF , où Rm = AT X + XA − XBB T + C T C dans une échelle logarithmique

à chaque itération m.
4.9.1 Exemple 1
Dans tous les tests effectués dans sous section, la matrice A provenant de la discrétisation
de l’opérateur
∂u ∂u
Lu = ∆u − f1 (x, y) + f2 (x, y) + g(x, y), (4.30)
∂x ∂y
matrice construite est n = n20 , où n0 est le nombre de points de la grille dans chaque di-
rection. La discrétisation de l’opérateur Lu construites des matrices dans la bibliothèque
Lyapack [89] et désignés par
dans le domaine Ω = [0 1] × [0 1] avec des conditions aux bord de Dirichlet homogènes,

où les fonctions f _1, f _2 et g sera précisées dans chaque exemple. Pour tous les test,
nous prenons r = 2.
Exemple 1.1
Dans cet exemple, nous considérons n = 25000 et

A= fdm2d_matrix(50,’x + y’,’exp(x ∗ y)’,’10’)
8
A.G
M.R
6
2
log10(||R(m)|| )
F
−2
−4
−6
−8
0 5 10 15 20 25 30 35 40
les iteration m
Exemple 1.2
Dans cet exemple, nous considérons n = 64000 et la matrice A
A= fdm2d_matrix(80,’cos(x + y)’,’sin(y 2 )’,’y 2 − x2 ’)
8
A.G
M.R
6
2
log10(||R(m)|| )
F
−2
−4
−6
−8
0 5 10 15 20 25 30 35 40
les iteration m
Figure 4.3: Résultats de l’exemple 1.2.

4.10 Conclusion 93
4.9.2 Exemple 2
Nous considérons l’exemple donné par Jbilou [64, 65] tel que la matrice A est donnée
par,  
a 1 − d ··· 0 1
 .. 

 1+d a . ··· 0


.. .. ..
 
A = − . . .
 
 0 1+d .
 .. .. 

 0 . . 1−d 

1 0 0 1+d a
Pour n = 4000, C = In×s , B = rand(n, 10), et d = 0.2, nous obtenons le résultat suivant,
Méthode Itér Temps k Rm kF

MR 5 13 s 9.00 × 10−9
GA 24 24 s 4.17 × 10−4
Table 4.1: Résultats de l’exemple 2
Notons que l’algorithme GA s’arrête à l’itération m = 25 : c’est un problème d’existence

et d’unicité de la solution de l’équation de Riccati projetée.
4.9.3 Exemple 3
Soit A la matrice opposée de la matrice précédente. Pour d = 0, 5, nous avons
Méthode N.Itérations Temps k Rm kF

MR 7 11 s 2.17 × 10−10
GA 11 6s 1.39 × 103
Table 4.2: Résultats de l’exemple 3.
Notons que l’algorithme GA s’arrête à l’itération m = 12 : c’est un problème d’existence

et d’unicité de la solution de l’équation de Riccati projetée.
4.10 Conclusion
Dans ce chapitre, nous avons proposé une nouvelle méthode pour la résolution numérique
de l’équation de Riccati continue de grande taille. Basée sur la méthode de minimisation
de résidu MR (Minimal Residual) et la méthode de projection dans un sous-espaces de
Krylov étendu par blocs, elle se ramène à chaque itérations à la résolution numérique
d’un problème de minimisation réduit. Ce dernier problème est résolu à l’aide de la

4.10 Conclusion 94
méthode de Newton et la méthode de MINRES globale. Les tests numériques montrent

l’efficacité de notre méthode.

Chapitre 5
Équation matricielle de Riccati non

symétrique et application à
l’équation de transport
5.1 Introduction
Les équations algébriques de Riccati non symétriques NAREs (Nonsymmetric Algebraic

Riccati Equation) jouent un rôle fondamental dans de nombreux problèmes de la théorie
des jeux différentiels, des modèles de Markov, et dans la théorie de transport [1, 24, 26,
32, 38, 41, 46, 71, 72, 110]. L’équation de Riccati non symétrique est de la forme suivante
R(X) = XCX − XD − AX + B = 0 (NARE) (5.1)
où A ∈ Rn×n , D ∈ Rp×p , B ∈ Rn×p , et C ∈ Rp×n . On peut associer l’équation de Riccati

non symétrique (5.1) à la matrice M qui se présente sous la forme suivante
!
D −C
M= ∈ R(n+p)×(n+p) . (5.2)
−B A
On suppose que M est une M-matrice inversible. Nous rappelons que M est une M-
matrice si elle peut être écrite comme M = aI − N où N matrice avec des éléments
non-négatifs et a ≥ ρ(N ), où ρ(N ) est le rayon spectral de la matrice N . Cela assure
l’existence de la solution non négative minimale X ∗ de l’équation de Riccati non symé-
trique (5.1), pour plus de détails voir [24, 25, 41, 46, 71, 72].
Soit H la matrice définie par
H = JM
95
5.1 Introduction 96
!
In 0
où J = , avec In la matrice identité de taille n. Donc nous avons
0 −Ip
!
D −C
H= . (5.3)
B −A
La solution de l’équation (5.1) est liée au calcul de sous espaces invariants de la matrice
H associée. En effet, si X est une solution de l’équation de Riccati non symétrique (5.1),
alors on a
! !
In In
H = (A − XB)
X X
Si les matrices Y , Z et W dans Rn×n , avec Z non singulière, satisfont :

! !
Y Y
H = W,
Z Z
alors X = Y Z −1 est une solution de l’équation algébrique de Riccati non symétrique,

pour plus de détails sur les propriétés de la matrice H et l’équation de Riccati non sy-
métrique voir [24–26, 41, 46, 110].
Le théorème de l’existence et l’unicité des solutions de Riccati non symétrique est liée à
la matrice M voir [24, 41, 46, 71, 72]. De plus, la structure particulière de la matrice M
assure l’existence de la solution minimale non négative X ∗ tel que X ∗ ≥ 0 et X ≥ X ∗
pour toute solution X de l’équation de Riccati non symétrique (5.1), pour plus de détails
voir [24, 41, 46]. Cela est indiqué dans le théorème suivant
Théorème 5.1.1. [46] Supposons que M soit une M-matrice. Alors l’équation NARE
(5.1) admet une unique solution minimale non négative X ∗ . De plus, si M est inversible,
alors A − X ∗ C et D − CX ∗ sont inversibles et M-matrices.
Pour les équations NAREs de petites ou moyennes tailles, il y a plusieurs méthodes

numériques pour trouver la solution non-négative minimale X ∗ de l’équation de Ric-
cati non symétrique. Par exemple, les méthodes exactes de Schur ont été étudiées dans
[41, 48, 49]. La méthode de Newton a également été étudiée dans [26, 41, 46, 48], elle
nécessite à chaque étape de calculer la solution d’une équation matricielle de Sylvester.
La méthode peut être coûteuse lorsqu’on résout l’équation de Sylvester par les méthodes
directs. D’autres méthodes comme l’algorithme SDA (Structure-preserving Doubling Al-
gorithm) [50] et la méthode ALI (Alternately Linearized Implicit) [10] ont été proposées
ces dernières années. En général, les méthodes itératives du point-fixe [25, 41, 46, 84–86]
sont moins coûteuses que la méthode de Newton, la méthode de Schur ou les méthodes

5.1 Introduction 97
SDA. Certaines techniques d’accélération de la convergence basées sur des méthodes

d’extrapolation vectorielles [69] ont été proposées récemment dans [38] pour accélérer la
convergence de certaines de ces méthodes itératives de point fixe, comme celles intro-
duites dans [38, 40, 47].
Dans ce chapitre, nous considérons l’équation matricielle de Riccati non symétrique

(NAREs) de grande taille avec le second membre s’écrit comme le produit de deux ma-
trices de petit rang. Ce type d’équations NAREs intervient notamment dans la théorie du
transport et ses applications. Nous nous proposons dans ce chapitre, deux nouvelles mé-
thodes pour la résolution numérique de l’équation de Riccati non symétrique (NAREs)
de grande taille. La première est une méthode itérative de projection sur des sous es-
paces de Krylov étendu EBA (Extended Block Arnoldi) [36, 56, 102]. Pour obtenir des
solutions approchées de rang inférieur "Low-rank", nous allons traiter aussi le cas parti-
culier correspondant à l’équation matricielle de Riccati non symétrique NAREs dans la
théorie de transport [38, 72].
La deuxième méthode, consiste à combiner la méthode de Newton-Kleinman et la mé-

thode d’Arnoldi par blocs. A chaque étape la méthode de Newton-Kleinman nécessi-
tera la résolution d’une équation matricielle de Sylvester de grande taille. Cette équa-
tion ne pouvait pas être résolu par les méthodes directes comme l’algorithme Bartels-
Stewart [11]. On peut utiliser les méthodes de de type Krylov, comme celles définies
dans [29, 56, 62, 65], pour résoudre l’équation de Sylvester de grande taille.
Le reste de ce chapitre est organisé comme suit : Dans la Section 2, nous allons combi-
ner la méthode de Newton-Kleinman avec les méthodes de sous-espaces de Krylov par
bloc ou bien des sous-espaces de Krylov étendus par blocs. La Section 3 est consacrée
à la résolution de l’équation NAREs de grande taille par une méthode de projection
à l’aide du sous espace de Krylov étendu et la condition d’orthogonalité de de Galer-
kin. Les solutions approchées de l’équation NAREs sont données sous la forme factorisée
(low rank approximation). L’application aux équations algébriques de Riccati non symé-
triques dans la théorie de transport sera considérée dans la Section 4. Enfin, la dernière
partie est consacrée à quelques expériences numériques.
Dans ce chapitre, nous utilisons les notations suivantes : Si K et L sont deux matrices
de taille n × m, alors K ≥ L si Li,j ≥ Ki,j pour tout i, j. La matrice K est dite non
négative si Li,j ≥ 0. Enfin, la séparation, par rapport à la norme de Frobenius, entre
deux matrices K et L est donnée par
sepF (K, L) = min kKX − XLkF .

kXkF =1

5.2 La méthode de Newton-Krylov par blocs 98
5.2 La méthode de Newton-Krylov par blocs
La méthode de Newton-Kleinman a été introduite par D.L. Kleinman en 1968 [79], pour
résoudre l’équation de Riccati continue, basée sur une linéarisation préliminaire par la
méthode de Newton, qui nous ramène à chaque étape à la résolution d’une équation
matricielle de Lyapunov avec un second membre décomposé (low rank). Nous nous atta-
chons dans ce paragraphe à résoudre l’équation de Riccati non symétrique. La méthode
de Newton nécessite à chaque itération la résolution d’une grande équation de Sylvester.
Cette équation est résolue par la méthode de Galerkin basée sur l’algorithme d’Arnoldi
par blocs. On suppose dans cette partie que
B, C 6= 0, (I ⊗ A + DT ⊗ I)V ect(B) > 0 (5.4)
et (I ⊗ A + DT ⊗ I) inversible et M -matrice. Ces conditions sont suffisantes pour assurer

la convergence de la méthode de Newton avec l’itération initiale de Newton X0 = 0.
Soit
R(X) = −XC1 C2T X + XD + AX − EF T . (5.5)
Nous appliquons la méthode de Newton-Kleinman pour résoudre l’équation
R(X) = 0.
La dérivée de Fréchet de l’application non linéaire R(X) en Xk est donnée par
R0Xk (Z) = (A − Xk C1 C2T )Z + Z(D − C1 C2T Xk ). (5.6)
En choisissant comme premier terme X0 = 0n×n , on construit la suite (Xk )k∈N des
approximations de la solution de (5.5) définie par la relation de récurrence
R0Xk (Xk+1 − Xk ) = −R(Xk ). (5.7)
On montre facilement que les itérations Xk+1 de Newton peuvent être vues comme
solutions de l’équation matricielle de Sylvester
Ak Xk+1 + Xk+1 Dk + Lk MkT = 0 (5.8)
où Ak = A − Xk C1 C2T , Dk = D − C1 C2T Xk , Lk = [Xk C1 , −E] et Mk = [XkT C2 , F ].
La convergence de la méthode de Newton pour la résolution de l’équation de Riccati

non symétrique est donnée par le théorème suivant ;
Théorème 5.2.1 (Guo et Higham,(2007) [46]). Si X ∗ la solution non négative minimale

de l’équation (5.14), alors les itérés de la méthode de Newton {Xi } sont bien défini, avec

5.2 La méthode de Newton-Krylov par blocs 99
X0 = 0, et nous avons
∀k ≥ 1, X0 ≤ X1 ≤ Xk ≤ X ∗ , et lim (Xk ) = X ∗ .
k→∞
De plus, si la matrice M est inversible la convergence est quadratique.
L’algorithme suivant décrit le processus de la méthode de Newton-Kleinman pour ré-

soudre l’équation de Riccati non symétrique
Algorithm 9 [La Méthode de Newton-Kleinman pour résoudre NAREs]

– On choisit une estimation initiale X0 , une tolérance tol et on fixe kmax.
– Pour k = 0, . . . , kmax
– Calculer Xk+1 la solution de l’équation de Sylvester de grande taille suivante
Ak Xk+1 + Xk+1 Dk + Lk MkT = 0 (5.9)
– Si k R(Xk+1 ) kF < tol, Arrêt

– Fin pour.
Dans chaque itération de l’algorithme de Newton-Kleinman, on résout l’équation ma-

tricielle de Sylvester avec le second membre donné comme le produit de deux matrices
de petit rang. Pour les petites et moyennes dimensions, on peut utiliser des méthodes
directes comme la méthode Bartels-Stewart [11]. Dans le cas où n est grand, les mé-
thodes de projection sur des espaces de type Krylov, plusieurs méthodes numériques ont
été proposées ces dernières années ; voir [29, 55, 56, 64, 65, 100, 102]. Ici, nous avons
utilisé des sous espaces de Krylov étendus par blocs (ou le sous espaces de Krylov par
bloc) pour résoudre l’équation matricielle de Sylvester de grande taille (5.8). Le procédé
est défini comme suit : On applique d’abord l’algorithme d’Arnoldi étendu par blocs
(ou l’algorithme d’Arnoldi par blocs) à (Ak , Lk ) et (DkT , Mk ) à l’étape m, il permet de
générer, les matrices orthonormés Vm,k et Wm,k dont les colonnes sont des bases ortho-
normales du sous espace de Krylov étendu par blocs de Km (Ak , Lk ) et Km (DkT , Mk ),
m de la solution X
respectivement. Les approximations Xk+1 k+1 sont données par :
m T
Xk+1 = Vm,k Zm Wm,k , (5.10)
et Zm est obtenu à partir de la condition d’orthogonalité de Galerkin :
T m
Vm,k R(Xk+1 ) Wm,k = 0, (5.11)
m ) = A X m + X m D + L M T est le résidu de X m
où R(Xk+1 k k+1 k+1 k k k k+1 correspondant à
l’équation de Sylvester (5.8). Par conséquent, nous obtenons la matrice Zm est la solution
de l’équation de Sylvester de petit taille suivante :
T
T(A) (D)
m Zm + Zm Tm + Lm MTm = 0, (5.12)
5.3 Solution approchée de rang inférieur de NAREs 100
(A) T A V (D) T DT W T T
où Tm = Vm,k k m,k , Tm = Wm,k k m,k , Lm = Vm,k Lk et Mm = Wm,k Mk .
L’équation de Sylvester de taille réduite (5.12) sera résolu par la méthode de Bartels-
m ) k sans avoir calculer la
Stewart [11]. On peut calculer la norme du résidu k R(Xk+1 F
m , voir [56, 62] pour plus de détails. Ici l’approximation X m
solution approchée Xk+1 k+1
peut être exprimée comme un produit de deux matrices de rang inférieur et cela permet
d’économiser de la place mémoire.
Notons que dans le cas de l’utilisation de l’algorithme d’Arnoldi étendu par blocs pour
résoudre l’équation matricielle de Sylvester de grande taille (5.9), nous avons besoin de
calculer produits de la forme Ak−1 Y et Dk−T Y où Ak = A−Xk C1 C2T et Dk = D−C1 C2T Xk
avec Xk sous la forme suivante Xk = Z1 Z2T . Puisque A et D sont creuses, les matrices
Ak et Dk sont pas plus creuses, puis le calcul de les produits A−1 −T
k Y et Dk Y devient
très coûteux. La meilleure méthode pour éviter ce problème est d’utiliser la formule de
Sherman-Morrison-Woodbury donnée par
(A + U V T )−1 Y = A−1 Y − A−1 U (I + V T A−1 U )V T A−1 Y, (5.13)
où U et V sont des matrices de taille n × r. Notons que, si nous utilisons la méthode

d’Arnodi par bloc [29, 62] pour résoudre l’équation matricielle de Sylvester (5.9), alors on
a seulement les produits de la forme Ak Y et DkT Y qui sont nécessaires. Généralement,
cette méthode nécessite plusieurs itérations, par rapport à la méthode d’Arnodi étendu
par bloc, pour obtenir de bonnes approximations mais pourrait être moins cher si les
produits inverse "matrice-matrice" sont dominants.
5.3 Solution approchée de rang inférieur de NAREs
Dans cette section, nous considérons les équations algébriques de Riccati non symétriques
avec un second membre donné comme produit de deux matrices E et F de petit rang.
Dans ce cas l’équation de Riccati non symétrique est de la forme suivante :
AX + XD − XC1 C2T X − EF T = 0, (LrNARE) (5.14)
où les matrices A et D sont de grande taille et inversibles. E ∈ Rn×s , F ∈ Rp×s , C1 ∈

Rp×s et C2 ∈ Rn×s avec s n, p. Ces équations matricielles présentent de nombreuses
applications telles que dans la théorie de transport et autres. A titre d’exemple, nous
considérons l’équation de Riccati non symétrique de la théorie de transport avec n =
p = 900 et c = α = 0.5 (voir la section suivante). Sur la figure 5.1, nous avons tracé la
décomposition en valeurs singulières (SVD) de la solution non négative minimale exacte
obtenue par la méthode de Schur.

5
10
0
10
−5
10
−10
10
−15
10
−20
10
0 100 200 300 400 500 600 700 800 900
Figure 5.1: Les valeurs singulières de la solution non négative minimale de l’équation
NAREs
Comme le montre cette figure, les valeurs singulières décroisent rapidement vers zéro,
donc le rang numérique de la solution exacte (non négative minimale) est petit, ce qui
donne une solution de rang inférieur "low-rank". Ici on a rang(X ∗ ) = 20. Cette remarque
permet de rechercher des méthodes qui donne des approximations de petit rang "low-
rank" de la solution de l’équation LrNARE.
Dans la suite, nous proposons une méthode itérative de projection sur des sous espaces
de Krylov étendu par blocs pour résoudre l’équation LrNARE, basée sur l’algorithme
d’Arnoldi étendu par blocs et la condition d’orthogonalité de Galerkin. Cette méthode
est décrite ci-dessous comme suit :
Nous appliquons l’algorithme d’Arnoldi étendu par blocs à (A, E) et (DT , F ), nous obte-
nons les matrices orthogonales Vm+1 et Wm+1 et les matrices de Hessenberg supérieure
par blocs HA D
m et Hm . Nous définissons aussi les matrices de Hessenberg supérieures par
blocs suivantes
TA T A T
m = Vm AVm , Tm = Vm+1 A Vm ,
et
TD T T D T T
m = Wm D Wm , Tm = Wm+1 D Wm .
Comme nous avons expliqué précédemment que les matrices Hessenberg par blocs TA
m
et TD A D
m sont obtenus à partir Hm et Hm sans avoir calculer le produit "matrice-vecteur".
Ensuite, nous considérons des solutions approchées de petit rang de la forme
Xm = Vm Ym WTm , (5.15)

avec Vm = [V1 , . . . , Vm ], Wm = [W1 , . . . , Wm ] et Ym ∈ R2ms×2ms . La matrice Ym est

obtenue à partir de la condition de Galerkin suivante
VTm R(Xm )Wm = 0, (5.16)
où R(Xm ) est le résidu correspondant à l’approximation Xm et défini par :
R(Xm ) = AXm + Xm D − Xm C1 C2T Xm − EF T .
Nous utilisons la condition d’orthogonalité de Galerkin (5.16) et le fait que les matrices
Vm et Wm sont orthogonales, nous obtenons une équation projetée de type Riccati non
symétrique
TA D T T e eT
m Ym + Ym (Tm ) − Ym C1,m C2,m Ym − Em Fm = 0, (5.17)
où C1,m = WTm C1 , C2,m = VTm C2 , E

em = VT E et Fem = WT F .
m m
Notons que si Mm la matrice associée a l’équation projetée (5.17), c’est à dire

!
(TD
m)
T T
−C1,m C2,m
Mm = .
em Fe T
−E TmA
m
Soit la matrice orthogonale Qm définie par :

!
Wm 0
Qm = .
0 Vm
Alors nous avons la relation entre la matrice M associée à l’équation (5.14) et la matrice
Mm :
Mm = QTm M Qm . (5.18)
Nous supposons que l’équation projetée de Riccati non symétrique (5.17) à une solution
non négative minimale unique qui pourrait être obtenue par des méthodes classiques.
Nous avons le résultat suivant qui nous permet de calculer la norme du résidu sans
calculer la solution approchée Xm .
Théorème 5.3.1. Soit Ym la solution exacte de l’équation NARE réduite (5.17), obtenue
à l’itération m de l’algorithme d’Arnoldi étendu par blocs . Alors la norme de Frobenius
du résidu Rm = R(Xm ) est donnée par :
T
k Rm k2F =k Ym Em (TD 2 A T 2
m+1,m ) kF + k Tm+1,m Em Ym kF , (5.19)
où TA D A D
m+1,m et Tm+1,m sont les dernier bloc des matrices Tm et Tm respectivement.

Démonstration. Soit Vm+1 = [Vm , Vm+1 ] et Wm+1 = [Wm , Wm+1 ] les matrices ortho-
normales construites en appliquant simultanément m itérations de l’algorithme EBA
aux paires de matrices (A, E) et (DT , F ) respectivement, donc on a TA T
m = Vm AVm et
TD T T
m = Wm D Wm . Le résidu R(Xm ) peut être formulé comme un produit de matrices
R(Xm ) = AXm + Xm D − Xm C1 C2T Xm − EF T

" #
TA D T e eT e eT
m Ym + Ym (Tm ) − Ym C1 C2 Ym − E F
D
Ym Em (Tm+1,m )T
= Vm+1 A
Wm+1
Tm+1,m ETm Ym 02r
où Ce1 = WTm C1 , Ce2 = VTm C2 , E

e = VT E et Fe = WT F .
m m
Comme Ym est la solution de l’équation projetée
TA D T e eT e eT
m Ym + Ym (Tm ) − Ym C1 C2 Ym − E F = 0,
nous avons
!
0 D
Ym Em (Tm+1,m )T
Rm = R(Xm ) = Vm+1 A
WTm+1 . (5.20)
Tm+1,m ETm Ym 0
Puisque Vm+1 et Wm+1 sont orthonormales, nous obtenons :
BT A
k Rm k2F =k Ym Em Tm+1,m k2F + k Tm+1,m ETm Ym k2F .
Le théorème précédent joue un rôle très important dans la pratique, il permet de calculer
la norme du résidu sans calculer l’approximation de la solution et de minimiser le temps
d’exécution.
pour économiser de la mémoire, la solution approchée Xm = Vm Ym WTm pourrait être

donné comme un produit de deux matrices de petit rang. En effet, considérons la dé-
composition en valeurs singulières de la matrice Ym , c’est-à-dire
e Σ Ve T ,
Ym = U
où Σ est la matrice diagonale des valeurs singulières de Ym rangées dans l’ordre dé-
croissant. Soient U
el et Vel les matrices constituées des l premières colonnes de U
e et
de Ve correspondant aux l valeurs singulières supérieures ou égales certaine tolérance
el Σl Ve T où
dtol. Nous obtenons la décomposition de valeur singulière tronquée Ym ≈ U l
(1) 1/2 (2) 1/2
Σl = diag[σ1 , . . . , σl ]. Soit Zm = Vm U
el Σ , et Zm = Wm Vel Σ , il en résulte que
l l
(1) (2) T
Xm ≈ Zm (Zm ) . (5.21)

Nous donnons maintenant un résultat de perturbation
Théorème 5.3.2. Soit Xm l’approximation de rang inférieur de l’équation LrNARE

(5.14). Nous avons :
(A − Km )Xm + Xm (D − Jm ) − Xm C1 C2T Xm − EF T = 0, (5.22)
A
où Km = Vm+1 Tm+1,m VmT et Jm = Wm (Tm+1,m
D )T Wm+1
T .
Démonstration. En multipliant l’équation de Riccati non symétrique de dimension ré-

duite (1.5) à gauche et à droite par Vm et WTm respectivement, et en utilisant les les
relations suivante
AVm = Vm TA A T
m + Vm+1 Tm+1,m Em ,
et
DT Wm = Wm TD D T
m + Wm+1 Tm+1,m Em .
nous obtenons
A
[AVm − Vm+1 Tm+1,m ETm ]Ym WTm + Vm Ym [DT Wm − Wm+1 Tm+1,m
D
ETm ]T
−Vm Ym WTm BC T Vm Ym WTm + EF T = 0.
Comme Vm et Wm sont deux matrices orthonormées, cela nous permet de donner les
relations suivantes
A
Vm+1 Tm+1,m ETm Ym WTm = Vm+1 Tm+1,m
A
ETm VTm Xm
et
D
Vm Ym Em (Tm+1,m )T Wm+1
T D
= Xm Wm Em (Tm+1,m )T Wm+1
T
d’autre part, on a également
Vm Em = Vm et Wm Em = Wm .
Par conséquent
(A − Km )Xm + Xm (D − Jm ) − Xm BC T Xm + EF T = 0,
A
où Km = Vm+1 Tm+1,m VmT et Jm = Wm (Tm+1,m
D )T Wm+1
T .

Lorsque D = AT , C1 = C2 et E = F l’équation (5.14) sera une équation de Riccati

symétrique, et le résultat du théorème 5.3.2 coïncide dans ce cas avec le résultat de la
perturbation donnée dans l’article de Jbilou [65].
Ensuite, nous donnons un résultat montrant que l’erreur X − Xm est une solution exacte
de l’équation NARE perturbée.
Théorème 5.3.3. Soit Xm la solution approchée obtenue par m itération de l’algorithme

d’Arnoldi étendu par blocs, et X la solution de l’équation LrNARE (5.14). Alors, l’erreur
X − Xm est une solution de l’équation algébrique de Riccati non symétrique perturbée
suivante
Am (X − Xm ) + (X − Xm )Dm − (X − Xm )C1 C2T (X − Xm ) + Km Xm + Xm Jm = 0 (5.23)
où Am = A − Xm C1 C2T et Dm = D − C1 C2T Xm , Km = Vm+1 TA T

m+1,m Vm et Jm =
Wm (TD T T
m+1,m ) Wm+1 .
Démonstration. D’après le théorème (5.3.2), nous avons :
(A − Km )Xm + Xm (D − Jm ) − Xm C1 C2T Xm − EF T = 0. (5.24)
La soustraction (5.24) et (5.14), permet de donner l’équation suivante
A(X − Xm ) + (X − Xm )D − XC1 C2T X + Xm C1 C2T Xm + Km Xm + Xm Jm = 0. (5.25)
Enfin, en développant (5.25), nous obtenons (5.23) se qui termine la démonstration.
Nous donnons maintenant un résultat de majoration de la norme de l’erreur X − Xm au

pas m de l’algorithme d’Arnoldi étendu par blocs.
Théorème 5.3.4. Soit Xm l’approximation de faible rang de la solution exacte X de

γm η
l’équation LrNARE. Si δm = sepF (Am , −Dm ) > 0 et 2 < 1/4, il existe une solution
δm
X de (5.14) satisfaisant
2γm
k X − Xm kF ≤ p
2 − 4γ η
.
δm + δm m
où γm =k Rm kF et η =k C1 C2T kF .
Démonstration. Tout d’abord, remarquer que

!
0 Ym Em (TD
m+1,m )
T
Km Xm + Xm Jm = Vm+1 WTm+1 , (5.26)
TA T
m+1,m Em Ym 0

5.4 Applications à la théorie de transport 106
par conséquent, à partir de (5.20), nous obtenons
Km Xm + Xm Jm = Rm . (5.27)
Maintenant nous appliquons [ le Théorème 2.1 dans [104]] à l’équation algébrique de

Riccati non symétrique (5.23), se qui termine la démonstration.
Nous remarquons que δm = sepF (Am , −Dm ) est également défini par
δm = inf k T(P ) kF ,
kP k=1
où T(P ) = Am P + P Dm .
Les étapes de la méthode que nous proposons (pour résoudre l’équation de Riccati non
symétrique de grande taille) sont résumées dans l’algorithme suivant
Algorithm 10 [La méthode d’Arnoldi étendu par blocs pour résoudre NAREs]
– Entrées : les matrices A, D, C1 , C2 , E et F . Le nombre maximum d’itérations mmax ,
et les tolérances toler et dtol.
(1) (2)
– Sorties : les matrices Zm et (Zm ) de telle sorte que la solution approchée est donnée
par
(1) (2) T
Xm ≈ Zm (Zm ) .
– Pour m = 1, . . . , mmax
– Appliquer l’algorithme d’Arnoldi étendu par blocs à (A, E) et (DT , F ) pour obtenir
les matrices orthogonales Vm , Wm et les matrices de Hessenberg par blocs TA
m et
TD
m.
– Résoudre l’équation NAREs de taille réduite
TA D T T e eT
m Ym + Ym (Tm ) − Ym C1,m C2,m Ym − Em Fm = 0,
– Si k Rm kF < toler, Arrêt

– Fin pour.
Dans la section suivante, nous allons appliquer les méthodes itératives proposées pour
la résolution d’un particulier de l’équation algébrique de Riccati utilisée dans la théorie
de transport.
5.4 Applications à la théorie de transport
Dans cette section, on considère le cas particulier de l’équation algébrique de Riccati

non symétrique, largement utilisée dans la théorie de transport, voir [24, 40, 73, 92]. Le
5.4 Applications à la théorie de transport 107
problème de transport consiste à résoudre une équation intégrodifférentielle. Après la

discrétisation de cette équation intégrodifférentiel, le problème peut être exprimé sous
la forme de l’équation matricielle de Riccati non symétrique suivante
(∆ − eq T ) X + X(Γ − qeT ) − Xqq T X − eeT = 0, (5.28)
où les matrices et les vecteurs ci-dessus dépendent des paramètres c et α satisfaisant

0 < c ≤ 1, 0 ≤ α < 1, pour plus de détails voir [24, 38, 40, 73, 92]. Cette équation
particulière est de la forme LrNARE donnée à (5.14) avec les matrices suivantes
A = ∆ − eq T , D = Γ − qeT , C1 = C2 = q, et E = F = e. (5.29)
Les matrices ∆ et Γ dans l’équation NARE (5.28) sont données par
∆ = diag(δ1 , . . . , δn ); Γ = diag(γ1 , . . . , γn ), (5.30)
où
1 1
δi = , et γi = , i = 1, . . . , n. (5.31)
cxi (1 − α) cxi (1 + α)
Les vecteurs e et q sont donnés par
wi
e = (1, . . . , 1)T , q = (q1 , . . . , qn )T avec qi = , i = 1, . . . , n. (5.32)
2xi
Nous appliquons ensuite la méthode d’Arnoldi étendue à l’équation NARE (5.28) pour
obtenir des solutions approximatives bas de classement. Nous remarquons que l’appli-
cation de la méthode ci-dessus, nous utilisons des opérations « matrice-vecteur» de la
forme A−1 v et D−T v. Comme les matrices A et D impliquées dans l’équation algébrique
de Riccati non symétrique (5.28) sont la somme d’une matrice diagonale et d’une ma-
trice de petit rang, nous pouvons utiliser la formule Sherman-Morrison-Woodbury pour
calculer A−1 u et D−1 u, où u ∈ Rn .
Donc nous avons :
−1 ∆−1 e q T ∆−1 u
A−1 u = (∆ − eq T ) u = ∆−1 u + , (5.33)
1 − q T ∆−1 e
et
−1 Γ−1 q eT Γ−1 u
D−1 u = (Γ − qeT )(∆ − eq T ) u = Γ−1 u + . (5.34)
1 − eT Γ−1 q
Cette formule permet de réduire le temps de calcul et permet également de réduire la
place en mémoire.

5.5 Résultats numériques 108
5.5 Résultats numériques
Dans cette section, nous donnons quelques exemples numériques de l’équation de Riccati
non symétrique NAREs pour montrer l’efficacité des méthodes proposées. Les différentes
expériences numériques ont été réalisées sur un ordinateur Portable d’Intel Core i5 à
1,6 GHz et 4G de RAM. Les algorithmes ont été codés dans Matlab2009. Nous avons
comparé la performance de la méthode d’Arnoldi étendu par blocs et la méthode Newton-
Block-Arnoldi avec la méthode RRE (Reduced Rank Extrapolation) proposée par El-
Moallem et Sadok pour résoudre l’équation de Riccati non symétrique, pour plus de
détails voir [38]. Pour l’algorithme d’Arnoldi étendu par blocs, le critère d’arrêt était
kR(Xm )kF / kE F T kF < 10−11 ,
où la norme de kR(Xm )k résiduelle a été calculée en utilisant le théorème 5.3.1. Pour la

méthode de Newton-BA et la méthode RRE, les itérations étaient arrêtée lorsque
k Xk+1 − Xk kF /kXk kF < 10−11 .
Nous considérons le problème venant de la théorie de transport comme expliqué dans la

section 5. Nous avons utilisé différentes valeurs des paramètres α et c, et également diffé-
rent tailles du problème de transport. Le tableau (??) et le tableau (3.2) rapportent sur
les résultats obtenus par la méthode d’Arnoldi étendu par blocs pour résoudre NAREs
(algorithme 2) notée (EBA-NARE), la méthode d’Arnoldi Newton-Bloc (Newton-BA)
résumée dans l’algorithme 3 pour laquelle l’équation de Sylvester, apparaissant dans
chaque étape de l’itération de Newton, a été résolu de manière itérative par l’algorithme
d’Arnoldi par blocs et la méthode RRE [38]. Pour la méthode EBA-NARE, l’équation
de Riccati non symétrique projetée NARE (5.17) a été résolue en utilisant le programme
sda_affine_mnare de Bini, Iannazzo et Meini pour plus de détails voir[24].
Exemple 1.
Dans ce premier exemple, nous avons choisi c = 0.5 et α = 0.5. Nous avons d’abord
calculé les approximations XN et XE données par le Newton-BA et par l’EBA-NARE
pour le taille n = 2000 et obtenir le suivant norme d’erreur k XN − XE kF = 1.2 · 10−10
cette résultat signifie que la solution approchée XE de la méthode EBA-NARE obtenue
se rapproche de la solution non négative minimale XN obtenue par la méthode de New-
ton.
Nous considérons maintenant des problèmes avec les tailles suivantes n = 4000, n =
10000, n = 20000 et n = 36000. Dans le tableau (5.1), une liste des normes résiduelles re-
latives obtenues pour chaque méthode et le temps de calcul correspondant (en secondes).
5.5 Résultats numériques 109
Pour toutes les expériences, les itérations de la méthode de Newton ne dépassent pas 5
itérations. Le nombre maximum d’itérations intérieures était itermax = 50 et ces itéra-
tions intérieures ont également été arrêtée lorsque le résidu correspondant est inférieure
à tol = 10−10 . Nous remarquons que les approximations Xm produite par l’algorithme
EBA-NARE sont données sous la forme factorisée (5.21). Ces approximations sont nu-
mériquement de petit rang, par exemple, lorsque n = 4000, le rang de l’approximation
obtenue par la méthode EBA-NARE était égale à 22, et que le rang correspondant à
l’approximation obtenue par l’approche de Newton-BA était 23.
Table 5.1: Résultats pour les exemples de NAREs dans la théorie de transport avec
c = 0.5 et α = 0.5.
Méthode EBA-NARE Méthode Newton-BA Méthode RRE

n Res. Norm CPU time Res. Norm CPU time Res. Norm CPU time
4000 2.7 · 10−12 1.1s 3.5 · 10−12 7.5s 2.5 · 10−12 6.7s
10000 1.5 · 10−12 2.5s 4.5 · 10−12 34s 1.7 · 10−12 34.5s
20000 5.3 · 10−12 6.5s 6.5 · 10−12 415s 6.5 · 10−12 790s
36000 9.2 · 6.2−12 194s 5.2 · 10−12 1430s −−− > 3000s
Exemple 2.
Dans le deuxième exemple, nous considérons également l’équation matricielle de Riccati

non symétrique dans la théorie de transport avec c = 0.9999 et α = 10−8 . Cet exemple
correspond à un cas difficile à résoudre par les méthodes itératives (le problème lié à
la singularité de la matrice Jacobienne) pour plus de détails voir [38, 49, 52]. Dans le
tableau 5.2, nous avons présenté les normes des résidus et les temps de calcul CPU obtenu
pour les trois méthodes : (EBA-NARE), (Newton-BA) et méthode RRE. Un maximum
de kmax = 30 itérations a été autorisé aux itérations Newton extérieures. ici aussi, nous
avons utilisé les mêmes critères d’arrêt que dans l’exemple 1. Comme on le voit à partir
du Tableau 5.2, les temps de calcul pour la méthode de Newton sont beaucoup plus
élevés que ceux obtenus dans le tableau 5.1. Une possibilité pour éviter la difficulté « le
problème de singularité de la matrice Jacobienne » de la méthode Newton dans ce cas
est l’utilisation des techniques de décalage « shift technique», en transformant l’équation
de départ avec une technique de «shift», qui permet d’éviter le problème de singularité
pour la matrice Jacobienne. Pour plus de détails voir [38, 49, 52].
Table 5.2: Résultats pour les exemples de NAREs dans la théorie de transport avec
c = 0.9999 et α = 10−8 .
Méthode EBA-NARE Méthode Newton-BA Méthode RRE

n Res.Norm CPU time Res.Norm CPU time Res.Norm CPU time
4000 1.7 · 10−12 1.6s 3.5 · 10−11 35s 1.5 · 10−12 13.5s
10000 2.5 · 10−12 2.7s 4.5 · 10−11 125s 2.3 · 10−12 72s
20000 7.6 · 10−12 7.4s 1.5 · 10−11 1950s 4.4 · 10−12 2140s
36000 7.3 · 6.2−12 211s −−− > 3000s − − − > 3000s
120000 5.2 · 6.2−12 763s −−− > 3000s − − − > 3000s

5.6 Conclusion 110
Exemple 3. Dans ce troisième exemple, nous considérons l’équation de Riccati non

symétrique (LrNARE) donné en (5.12). Les matrices A et D sont obtenues à partir de
la discrétisation par différences finie centrées des opérateurs
∂u ∂u
LAu = ∆u − exy + sin(xy) + 100x, (5.35)
∂x ∂y
et
∂u ∂u
LDu = ∆u − x2 y + x cos(xy) + 10, (5.36)
∂x ∂y
respectivement, sur le carré unité [0, 1] × [0, 1] avec des conditions aux bord de Dirichlet
homogènes. Le nombre de points de la grille intérieure dans chaque direction a été n0
pour l’opérateur LAu et p0 pour l’opérateur LDu . La dimension des matrices A et C
sont n = n20 et p = p20 respectivement. Différents choix de valeurs ont été faits pour n0 et
m0 . Les matrices E, F , C1 et C2 sont des matrices aléatoires uniformément distribuées
dans l’intervalle [0, 1] et s = 5. Nous utilisé le même critère d’arrêt comme pour les
exemples premiers.
Le tableau 5.3 donne les résultats obtenus par les trois méthodes EBA-NARE, Newton-
BA et la méthode SDA [24, 50]. Pour les trois méthodes, nous avons reporté la norme
de Frobenius du résidu et le temps CPU (en secondes) de calcul. Comme indiqué sur le
tableau 5.3, les méthodes proposées donne les meilleurs résultats en terme de temps de
calcul.
Table 5.3: Résultats pour les exemples 4, comparaisons avec la méthode SDA.
Méthode EBA-NARE Méthode Newton-BA Méthode SDA

n, p Res.Norm CPU time Res.Norm CPU time Res.Norm CPU time
n = 400,
p = 225. 3.7 · 10−12 0.6s 6.5 · 10−12 1.8s 2.5 · 10−12 4.7s
n = 40000,
p = 22500 2.4 · 6.2−12 82.5s 4.5 · 10−12 124s −−−−−
n = 122500
p = 90000 6.3 · 6.2−12 240s 6.5 · 10−12 415s − − −−
5.6 Conclusion
Nous avons présenté dans ce chapitre, deux nouvelles méthodes itératives pour calculer la
solution approchée de rang inférieur d’équations algébriques de Riccati non symétriques
de grande taille. La première basée sur des sous-espace de Krylov étendus par blocs et la
condition d’orthogonalité de Galerkin. Dans la deuxième méthode, nous avons présenté
la méthode de Newton-bloc Arnoldi, basée sur la méthode de Newton et des sous-espaces
de Krylov par blocs : à chaque itération de la méthode de Newton, on résout une équation
de Sylvester de grande taille par l’utilisation de l’algorithme d’Arnoldi par blocs. Nous
5.6 Conclusion 111
avons également présenté des résultats de majoration de l’erreur. Nous avons appliqué
ces méthodes itératives à l’équation de Riccati non symétrique bien connue dans la
théorie du transport. Les résultats numériques montrent que les approches proposées
sont efficaces pour les problèmes de grande taille.

Bibliographie
[1] H. Abou-Kandil, G. Freiling, V. Ionescu, and G. Jank, Matrix Riccati Equations in

Control and Systems Theory, Systems and Control : Foundations and Applications,
Birkh auser, Boston, 2003.
[2] B. D. O. Anderson, and J. B. Moore, Linear Optimal Control, Prentice-Hall, En-

glewood Cliffs, 1971.
[3] B. D. O. Anderson, and J. B. Moore, Optimal Control-Linear Quadratic Methods,

Prentice-Hall, Englewood Cliffs, NJ, 1990.
[4] A. C. Antoulas, D.C. Sorensen, Projection methods for balanced model reduction,
Technical Report, Rice University, Houston, TX, 2001.
[5] A. C. Antoulas, Approximation of Large-Scale Dynamical Systems, SIAM, 2005.
[6] W. F. Arnold and A. J. Laub, Generalized eigenproblem algorithms and software

for algebraic Riccati equations. Proc. IEEE, 72, pp : 1746–1754, 1984.
[7] J. Baglama, L. Reichel, D. Richmond, An augmented LSQR method, Numer. Al-

gorithms, 64 (2013), pp. 263–293.
[8] A. Y. Barraud, A numerical algorithm to solve AT XA − X = Q, IEEE Trans.

Autom. Contr., AC-22(1977), pp. 883–885.
[9] Z-Z. Bai, Y-H. Gao, L-Z. Lu, Fast iterative schemes for nonsymmetric algebraic
Riccati equations arising from transport theory, SIAM J. Sci. Comp., 30(2008), pp.
804–818.
[10] Z-Z. Bai, X-X. Guo, S-F. Xu, Alternately linearized implicit iteration methods for
the minimal nonnegative solutions of the nonsymmetric algebraic Riccati equations,
Num. Lin. Alg. with App, 13(2006), pp. 655–674.
[11] R.H. Bartels, G.W. Stewart, Solution of the matrix equation AX + XB = C,

Algorithm 432, Comm. ACM 15(1972), pp. 820–826.
112
BIBLIOGRAPHIE 113
[12] M. Bellalij, K. Jbilou, H. Sadok, New convergence results on the global GMRES
method for diagonalizable matrices, Journal of Computational and Applied Mathe-
matics 219 (2008), pp. 350–358.
[13] P. Benner, Factorized solution of Sylvester equations with applications in control,

in : B. De Moor, B. Motmans, J. Willems, P. Van Dooren, V. Blondel (Eds.),
Proceedings of the Sixteenth International Symposium on Mathematical Theory of
Network and Systems, MTNS, Leuven, Belgium, 2004.
[14] P. Benner and R. Byers. An exact line search method for solving generalized conti-
nous algebraic Riccati equations. IEEE Trans. Automat. Control, 43, pp : 101–107,
(1998).
[15] P. Benner, J. Li, and T. Penzl. Numerical solution of large Lyapunov equations,
Riccati equations and linear-quadratic optimal control problems. Numer. Linear
Alg. Appl., 15, pp. 755–777, 2008.
[16] P. Benner and H. Faβbender. An implicitly restarted symplectic Lanczos method

for the hamiltonian eigenvalue problem. Linear Alg. Appl., 263, pp. 75–111, 1997.
[17] A. Bouhamidi , K. Jbilou, Stien implicit Runge-Kutta methods with high stage order
for large-scale ordinary differential equations. Applid Numerical Mathematics, 61,
pp. 149–159, 2011.
[18] A. Bouhamidi , K. Jbilou, Sylvester Tikhonov-regularization methods in image

restoration, J. Comput. Appl. Math., 206(1), pp. 86–98, 2007.
[19] A. Bouhamidi , K. Jbilou, A note on the numerical approximate solutions for ge-
neralized Sylvester matrix equations with applications, Applied Mathematics and
Computation 206 (2008), pp. 687–694.
[20] S. Bittanti, A. Laub, and J. C. Willems, The Riccati equation, Springer-Verlag,

Berlin, 1991.
[21] P. Benner, R.C. Li, N. Truhar, On the ADI method for Sylvester equations, J. of
Comput. Appl. Math., (233)2009, pp. 1035–1045.
[22] P. Benner, H. Mena, and J. Saak, On the parameter selection problem in the
Newton-ADI iteration for large-scale Riccati equations, Electron. Trans. Numer.
Anal., 29 (2008), pp. 136–149.
http ://etna.math.kent.edu/vol.29.2007–2008/pp136–149.dir/
[23] P. Benner and P. Kürschner, Computing real low-rank solutions of Sylvester equa-
tions by the factored ADI method, Comput. Math. Appl., 67 (2014), pp. 1656–1672.
[24] D.A. Bini, B. Iannazzo, B. Meini, Numerical Solution of Algebraic Riccati Equa-
tions, SIAM, Philadelphia, PA, 2012.

BIBLIOGRAPHIE 114
[25] D.A. Bini, B. Iannazzo, G. Latouche, B. Meini, On the solution of algebraic Riccati
equations arising in fluid queues, Linear Algebra Appl., 413(2006), pp. 474–494.
[26] D.A. Bini, B. Iannazzo and F. Poloni, A fast Newtons method for a nonsymmetric
algebraic Riccati equation, SIAM J. Matrix Anal. Appl., 30(2008), pp : 276–290.
[27] R. Byers, Solving the algebraic Riccati equation with the matrix sign function.
Linear Algebra and its Applications, 85, pp : 267–279, 1987.
[28] J-P Chehab, M. Raydan, An implicit preconditioning strategy for large-scale ge-
neralized Sylvester equations, Applied Mathematics and Computation, 217(2011),
pp : 8793–8803.
[29] A. Bouhamidi, M. Hached, M. Heyouni and K. Jbilou, A preconditioned block

Arnoldi for large Sylvester matrix equations, Numer. Lin. Alg. Appl., 2(2013), pp :
208–219.
[30] R. Bouyouli, K. Jbilou, R. Sadaka, H. Sadok, Convergence properties of some block

Krylov subspace methods for multiple linear systems, Journal of Computational
and Applied Mathematics, 196 (2006), pp. 498–511.
[31] Callier, F.M and Desoer, Linear System Theory, SpringerVerlag (1991).
[32] B. N. Datta, Numerical Methods for Linear Control Systems Design and Analysis,
Elsevier Academic Press, 2003.
[33] B.N. Datta, Numerical Methods for Linear Control Systems, Elsevier Academic
press, 2004.
[34] B.N. Datta, Krylov-subspace methods for large scale matrix problems in control,
Generation of Computer Systems, 19(2003), pp. 125–126.
[35] B.N. Datta, K. Datta, Theoretical and computational aspects of some linear algebra
problems in control theory, in :C.I. Byrnes, A. Lindquist (EDs), Computational and
Combinatorial Methods in Systems Theory, Elsevier, Amsterdam (1986), pp 201–
212.
[36] V. Druskin, L. Knizhnerman, Extended Krylov subspaces : approximation of the

matrix square root and related functions, SIAM J. Matrix Anal. Appl. 19 (3) (1998),
pp. 755–771.
[37] E. J. C. Doyle, K. Glover, P. P. Khargonekar, and B. A. Francis, State space so-

lutions to standard H2 and H∞ control problems, IEEE Trans. Automat. Control,
34 (1989), pp. 831–846.
[38] R. El-Moallem and H. Sadok, Vector extrapolation methodsapplied to algebraic

Riccati equations arising in transport theory, Elect. Trans. Numer. Anal., 40(2013),
pp. 489–506.
BIBLIOGRAPHIE 115
[39] A. El Guennouni, K. Jbilou, A.J. Riquet, Block Krylov subspace methods for solving
large Sylvester equations, Numer. Alg. 29 (2002), pp.75–96.
[40] B. D. Ganapol, An investigation of a simple transport model, Transport Theory

Statist. Phys., 21(1992), pp. 1–37.
[41] C-H. Guo, Nonsymmetric algebraic Riccati equations and Wiener-Hopf factoriza-
tion for M-matrices, SIAM J. Matrix Anal. Appl., 23(1)(2001), pp. 225–242.
[42] K. Glover, D.J. N. Limebeer, J.C. Doyle, E.M. Kasenally and M.G. Safonov, A
characterisation of all solutions to the four block general distance problem, SIAM
J. Control Optim., 29(1991), pp. 283–324.
[43] G.H. Golub, and W. Kahan, Calculating the singular values and pseudoinverse of
a matrix, SIAM J. Numer. Anal. 2 (1965), pp. 205–224.
[44] G.H. Golub and C.F. Van Loan, Matrix Computations, (Johns Hopkins Studies in
Mathematical Sciences)(3rd Edition). The Johns Hopkins University Press,1996.
[45] G.H. Golub, S. Nash, C. Van Loan, A Hessenberg Schur method for the problem
AX + XB = C , IEEE Trans. Automat. Control 24 (1979), pp. 909–913.
[46] C. H. Guo, N. J. Higham, Iterative solution of a nonsymmetric algebraic Riccati

equation. SIAM, J. Matrix Anal. Appl. 29(2)(2007), pp : 396–412.
[47] C. H. Guo, Nonsymmetric algebraic Riccati equations and Wiener-Hopf factoriza-

tion for M -matrices, SIAM J. Matrix Anal. Appl. 23 (1) (2001), pp 225–242.
[48] C-H. Guo and A. J. Laub, On the iterative solution of a class of nonsymmetric
algebraic Riccati equations, SIAM J. Matrix Anal. Appl., 22(2000), pp. 376–391.
[49] C.H. Guo, B. Iannazzo, and B. Meini, On the doubling algorithm for a (shifted)
nonsymmetric algebraic Riccati equation, SIAM J. Matrix Anal. Appl., 29(4)(2007),
pp : 1083–1100.
[50] X.-X. Guo, W.-W. Lin, and S.-F. Xu, A structure-preserving doubling algorithm for
nonsymmetric algebraic Riccati equation, Numer. Math., 103(2006), pp. 393–412.
[51] Guang-Da Hu and Qiao Zhu, Bounds of modulus of eigenvalues based on Stein
equation. Kybernetika, 46(4), pp. 655–664, 2010.
[52] C. He, B. Meini, NH. Rhee, A shifted cyclic reduction algorithm for quasi- birth-
death problems, SIAM Journal on Matrix Analysis and Applications, 23(3)(2001),
pp. 673–691.
[53] S.J. Hammarling, Numerical solution of the stable, nonnegative definite Lyapunov
equation, IMA J. Numer. Anal. 2 (1982), pp. 303–323.

BIBLIOGRAPHIE 116
[54] R.E. Hartwig, Resultants and the solution of AX − XB = −C, SIAM J. Appl.
Math. 23 (1) (1972), pp. 104–117.
[55] M. Heyouni, K. Jbilou, An extended Block Arnoldi algorithm for large-scale solu-
tions of the continuous-time algebraic Riccati equation, ETNA 33 (2009), pp. 53–62.
http ://etna.mcs.kent.edu/vol.33.2008-2009/pp53-62.dir
[56] M. Heyouni, Extended Arnoldi methods for large low-rank Sylvester matrix equa-
tions, Appl. Num. Math., 60(11)(2010), pp. 1171–1182.
[57] D.Y. Hu, L. Reichel, Krylov-subspace methods for the Sylvester equation, Linear
Algebra Appl. 172 (1992), pp. 283–313.
[58] M. Herty, R. Pinnau, M. Sead, On Optimal Control Problems in Radiative Transfer,

Optimization Methods and Software, 22(2007), pp. 917–936.
[59] Carl Jagels and Lothar Reichel, Recursion relations for the extended Krylov sub-
space method. Linear Algebra Appl., 434(7)(2011), pp. 1716–1732.
[60] I. M. Jaimoukha and E. M. Kasenally, Krylov subspace methods for solving large
Lyapunov equations, SIAM J. Numer. Anal.,31(1994), pp. 227–251.
[61] K. Jbilou, ADI preconditioned Krylov methods for large Lyapunov matrix equa-
tions, Linear Algebra and its Applications. 432 (2010) pp. 2473–2485.
[62] K. Jbilou, Low-rank approximate solution to large Sylvester matrix equations, App.
Math. Comput., 177(2006), pp. 365–376.
[63] K. Jbilou, Approximate solutions to large Stein matrix equations, World Academy
of Science, Engineering and Technology Vol :70(2012), pp : 19–25.
[64] K. Jbilou, Block Krylov subspace methods for large continuous-time algebraic Ric-
cati equations, Num. Alg., 34(2003), pp : 339–353.
[65] K. Jbilou, An Arnoldi based algorithm for large algebraic Riccati equations, Appl.
Math. Lett., 19 (2006), pp. 437–444.
[66] K. Jbilou and A. Messaoudi, A Computational Method for Symmetric Stein Ma-
trix Equations, Numerical Linear Algebra in Signals, Systems and Control. Volume
80(2011), pp 295–311.
[67] K. Jbilou, A. Messaoudi, H. Sadok, Global FOM and GMRES algorithms for matrix
equations, Appl. Numer. Math. 31(1999), pp. 49–63.
[68] K. Jbilou and A. J. Riquet, Projection methods for large Lyapunov matrix equa-
tions, Lin. Alg. Appl, 415(2006), pp. 344–358.
[69] K. Jbilou and H. Sadok, Vector extrapolation methods. Applications and numerical
comparison, J. Comput. Appl. Math., 122(2000), pp : 149–165.
BIBLIOGRAPHIE 117
[70] K. Jbilou, H. Sadok, Global Lanczos-based methods with applications, Technical

Report LMA 42, Univérsité du Littoral, Calais, France, 1997.
[71] Jonq Juang and Wen-Wei Lin, nonsymmetric algebric Riccati equations and
Hamiltonian-Like matrices. SIAM J. MATRIX ANAL. APPL. Vol. 20 (1998), pp.
228–243.
[72] J. Juang, Existence of algebraic matrix Riccati equations arising in transport theory,
Lin. Alg. Appl., 230(1995), pp : 89–100.
[73] J. Juang and D. Chen, Iterative solution for a certain class of algebraic matrix
Riccati equations arising in transport theory, Transp. Theor. Statis. Phys., 22(1),
pp : 65–80.
[74] Z. Jia and Y. Sun, A QR decomposition based solver for the least squares problems
from the minimal residual method for the Sylvester equation, J. Comput. Math.,
25(2007), pp. 531–542.
[75] S. Karimi, B. Zali, The block preconditioned LSQR and GL-LSQR algorithms for
the block partitioned matrices, Applied Mathematics and Computation, 227(2014),
pp. 811–820.
[76] A. Klein and P.J.C. Spreij, On Stein’s equation, Vandermonde matrices and Fisher’s
information matrix of time series processes. Part I : The autoregressive moving
average process. Lin. Alg. Appl., 329(13), pp. 9–47, 2001.
[77] André Klein and Peter Spreij. On the solution of Stein’s equation and Fisher’s
information matrix of an ARMAX process. Linear Algebra and its Applications,
396(2005), pp :1–34.
[78] D. Kressner and C. Tobler, Low-rank tensor Krylov subspace methods for parame-
trized linear systems, SIAM J. Matrix Anal. Appl., 32 (2011), pp. 1288–1316.
[79] D.L. Kleinman, On an iterative technique for Riccati equation computations, IEEC
Trans. Autom. Contr., 13(1968), pp : 114–115.
[80] P. Lancaster, and L. Rodman, Algebraic Riccati Equations, Clarendon Press, Ox-
ford, (1995).
[81] P. Lancaster and M. Tismenetsky, The Theory of Matrices, Academic Press, Or-
lando, 1985.
[82] Y. Lin and V. Simoncini, Minimal residual methods for large scale Lyapunov equa-
tions, App. Num. Math., 72(2013), pp : 52–71.
[83] P. Lancaster, L. Rodman, The Algebraic Riccati Equations, Clarendon Press, Ox-
ford, 1995.

BIBLIOGRAPHIE 118
[84] Y. Lin, I. Bao and Y. Wei, A modified Newton method for solving non-symmetric
algebraic Riccati equations arising in transport theory, IMA J. Numer. Anal.,
29(2008), pp : 215–224.
[85] L.-Z. Lu, Newton iterations for a non-symmetric algebraic Riccati equation, Numer.
Lin. Alg. Appl., 12(2005), pp : 191–200.
[86] L.-Z. Lu, Solution form and simple iteration of a nonsymmetric algebraic Riccati
equation arising in transport theory, SIAM J. Matrix Anal. Appl., 26(2005), pp :
679–685.
[87] J. Lasalle, S. Lefschetz, Stability of Lyapunov Direct Methods, Academic Press,

New York, 1961.
[88] A. J. Laub, A Schur method for solving algebraic Riccati equations, IEEE Trans.
Automat. Control, 24 (1979), pp. 913–921.
[89] T. Penzl, LYAPACK : A MATLAB toolbox for large Lyapunov and Riccati equa-
tions, model reduction problems, and linear-quadratic optimal control problems,
software available at
https ://www.tu-chemnitz.de/sfb393/lyapack/.
[90] C. Moler and C. Van loan, Ninteen dubious ways to compute the exponential of a
matrix, SIAM Rev., 20(1978), pp. 801–836.
[91] V. Mehrmann. The Autonomous Linear Quadratic Problem, Theory and Numerical
Solution. Lecture Notes in Control and Information Sciences, Vol. 63, Springer,
Heidelberg, (1995).
[92] V. Mehrmann and H. Xu, Explicit solutions for a Riccati equations from transport
theory, SIAM J. Matrix Anal., 30(2008), pp : 1339–1357.
[93] V. Mehrmann, H. Xu, Explicit solutions for a Riccati equation from transport
theory, SIAM J. Matrix Anal. Appl. 30 (4)(2008), pp : 1339–1357.
[94] C. C. Paige and M. A. Saunders, Solution of sparse indefinite systems of linear

equations. SIAM J. Numer. Anal., 12(4)(1975), pp : 617–629.
[95] C.C. Paige, A. Saunders, LSQR : an algorithm for sparse linear equations and sparse
least squares, ACM Trans. Math. Software 8 (1982), pp. 43–71.
[96] Frédéric Rotella, Pierre Borne, Théorie et pratique du calcul matriciel, 1995.
[97] L. C. G. Rogers, Fluid models in queueing theory and Wiener-Hopf factorization of

Markov chains, Ann. Appl. Probab., 4(2)(1994), pp : 390–413.
[98] Y. Saad, Numerical solution of large Lyapunov equations, in Signal Processing,

Scattering, Operator Theory and Numerical Methods, M.A. Kaashoek, J.H. Van
Shuppen and A.C. Ran, eds., Birkhaser, Boston, 1990, pp. 503–511.
BIBLIOGRAPHIE 119
[99] Y. Saad and M. H. Schultz, GMRES : A generalized minimal residual algorithm for
solving nonsymmetric linear systems, SIAM J. Sci. Statis. Comput., 7 (1986), pp.
856–869.
[100] Y. Saad, Iterative Methods for Sparse Linear Systems, Society for Industrial and
Applied Mathematics 2éme edition. SIAM, Philadelpha, PA, 2003.
[101] M. Sadkane, Block Arnoldi and Davidson methods for unsymmetric large eigen-
value problems, Numer. Math., 64(1993), pp. 687–706.
[102] V. Simoncini, A new iterative method for solving large-scale Lyapunov matrix
equations, SIAM J. Sci. Comput, 29(2007), pp. 1268–1288.
[103] E. de Souza, S.P. Bhattacharyya, Contollability, observability and solution of AX −

XB = C, Linear Algebra Appl. 39 (1981), pp. 167–188.
[104] G.W. Stewart and J.G. Sun, Matrix Perturbation Theory, Academic Press, New
York, (1990).
[105] R.A. Smith, Matrix calculations for Liapunov quadratic forms, J. Different. Eqs.
2(1966), pp. 208–217.
[106] F. Toutounian and S. Karimi, Global least squares method (Gl-LSQR) for sol-
ving general linear systems with several right-hand sides, Applied Mathematics
and Computation 178 (2006), pp.452–460.
[107] P. Van Dooren, A generalized eigenvalue approach for solving Riccati equations.
SIAM J. Sci. Statist. Comput., 2, pp :121–135, (1981).
[108] P. Van Dooren, Gramian based model reduction of large-scale dynamical systems,
in Numerical Analysis, Chapman and Hall, pp. 231–247, CRC Press London, 2000.
[109] W.M. Wonham, On a matrix Riccati equation of stochastic control. SIAM J.

Contr., 6, pp : 681–697, 1968.
[110] D. Williams, A potential-theoretical note on the quadratic Wiener-Hopf equation

for Q-matrices, in Seminar on Probability XVI, Lecture Notes in Mathematics 920,
pp. 91–94, Springer-Verlag, Berlin, 1982.
[111] Zhou, K., Doyle, J.C. and Glover, Robust Optimal Control. Prentice Hall, New-
Jersey (1995).

Résumé
Nous nous intéressons dans cette thèse, à l’étude des méthodes itératives pour la réso-
lution d’équations matricielles de grande taille : Lyapunov, Sylvester, Riccati et Riccati
non symétrique.
L’objectif est de chercher des méthodes itératives plus efficaces et plus rapides pour ré-
soudre les équations matricielles de grande taille. Nous proposons des méthodes itératives
de type projection sur des sous espaces de Krylov par blocs
Km (A, V ) = Image{V, AV, . . . , Am−1 V },
ou des sous espaces de Krylov étendus par blocs
e
Km (A, V ) = Image{V, A−1 V, AV, A−2 V, A2 V, · · · , Am−1 V, A−m+1 V }
Ces méthodes sont généralement plus efficaces et rapides pour les problèmes de grande
dimension.
Nous avons traité d’abord la résolution numérique des équations matricielles linéaires :
Lyapunov, Sylvester et Stein. Nous avons proposé une nouvelle méthode itérative basée
sur la minimisation de résidu MR et la projection sur des sous espaces de Krylov étendus
e (A, V ). L’algorithme d’Arnoldi étendu par blocs permet de donner un
par blocs Km
problème de minimisation projeté de petit taille. Le problème de minimisation de taille
réduit est résolu par différentes méthodes directes ou itératives. Nous avons présenté
aussi la méthode de minimisation de résidu basée sur l’approche global à la place de
l’approche bloc. Nous projetons sur des sous espaces de Krylov étendus Global
e
Km (A, V ) = sev{V, A−1 V, AV, A−2 V, A2 V, · · · , Am−1 V, A−m+1 V }.
Nous nous sommes intéresses en deuxième lieu à des équations matricielles non linéaires,
et tout particulièrement l’équation matricielle de Riccati dans le cas continu et dans
le cas non symétrique appliquée dans les problèmes de transport. Nous avons utilisé la
méthode de Newton et l’algorithme MINRES pour résoudre le problème de minimisation
projeté. Enfin, nous avons proposé deux nouvelles méthodes itératives pour résoudre les
équation de Riccati non symétriques de grande taille : la première basée sur l’algorithme
d’Arnoldi étendu par bloc et la condition d’orthogonalité de Galerkin, la deuxième est
120
BIBLIOGRAPHIE 121
de type Newton-Krylov, basée sur la méthode de Newton et la résolution d’une équation

de Sylvester de grande taille par une méthode de type Krylov par blocs. Pour toutes ces
méthodes, les approximations sont données sous la forme factorisée, ce qui nous permet
d’économiser la place mémoire en programmation. Nous avons donné des exemples nu-
mériques qui montrent bien l’efficacité des méthodes proposées dans le cas de grandes
tailles.
Mots clés : Sous espaces de Krylov étendu, méthodes blocs et globales, Approximation
de rang inférieur, Lyapunov, Sylvester, Riccati continu, Riccati non symétrique, Méthode
de Newton, Théorie de transport, Condition de Galerkin, Méthode de minimisation de
résidu.

Abstract
In this thesis, we focus in the studying of some iterative methods for solving large ma-
trix equations such as Lyapunov, Sylvester, Riccati and nonsymmetric algebraic Riccati
equation. We look for the most efficient and faster iterative methods for solving large
matrix equations. We propose iterative methods such as projection on block Krylov
subspaces
Km (A, V ) = Range{V, AV, . . . , Am−1 V },
or block extended Krylov subspaces
e
Km (A, V ) = Range{V, A−1 V, AV, A−2 V, A2 V, · · · , Am−1 V, A−m+1 V }.
These methods are generally most efficient and faster for large problems.
We first treat the numerical solution of the following linear matrix equations : Lyapunov,
Sylvester and Stein matrix equations. We have proposed a new iterative method based
e (A, V ).
on Minimal Residual MR and projection on block extended Krylov subspaces Km
The extended block Arnoldi algorithm gives a projected minimization problem of small
size. The reduced size of the minimization problem is solved by direct or iterative me-
thods. We also introduced the Minimal Residual method based on the global approach
instead of the block approach. We projected on the global extended Krylov subspace
e
Km (A, V ) = Span{V, A−1 V, AV, A−2 V, A2 V, · · · , Am−1 V, A−m+1 V }.
Secondly, we focus on nonlinear matrix equations, especially the matrix Riccati equation
in the continuous case and the nonsymmetric case applied in transportation problems.
We used the Newton method and MINRES algorithm to solve the projected minimization
problem. Finally, we proposed two new iterative methods for solving large nonsymme-
tric Riccati equation : the first based on the algorithm of extended block Arnoldi and
Galerkin condition, the second type is Newton-Krylov, based on Newton’s method and
the resolution of the large matrix Sylvester equation by using block Krylov method.
122
BIBLIOGRAPHIE 123
For all these methods, approximations are given in low rank form, wich allow us to save
memory space. We have given numerical examples that show the effectiveness of the
methods proposed in the case of large sizes.
Keywords : Extended Krylov subspaces, Low-rank approximation, Lyapunov equation,

Matrix Sylvester equation, Riccati equation, Nonsymmetric Riccati equation, Transport
theory, Newton method, minimal residual method MR, Galerkin condition.

These Doctora

Transféré par

Droits d'auteur :

Formats disponibles

These Doctora

Transféré par

Informations du document

Copyright

Formats disponibles

Partager ce document

Partager ou intégrer le document

Options de partage

Avez-vous trouvé ce document utile ?

Ce contenu est-il inapproprié ?

Droits d'auteur :

Formats disponibles

These Doctora

Transféré par

Droits d'auteur :

Formats disponibles

Méthodes itératives pour la résolution d’équations

To cite this version:

HAL Id: tel-01545833

HAL is a multi-disciplinary open access L’archive ouverte pluridisciplinaire HAL, est

Préparée dans le cadre d’une cotutelle entre

Spécialité : Mathématiques Appliquées et Informatique

CED : Sciences de l’Ingénieur

Méthodes itératives pour la résolution

Thèse présentée par : SADEK EL MOSTAFA

Abouir Jilali Professeur, FST Université Hassan II, Mohammedia, Maroc

Sadok Hassane Professeur, Université du Littoral Côte d’Opale, France

Bentbib Abdeslem Hafid Professeur, FST, Université Cadi Ayaad, Marrakech

Cette thèse a été préparée aux :

Laboratoire de Mathématiques Appliquées et Informatique

Laboratoire de mathématiques pures et appliquées

EL MOSTAFA SADEK Thèse de Doctorat

- Nom et Prénom de l’auteur : SADEK EL MOSTAFA

- Khalide JBILOU, Professeur de l’enseignement supérieur

- Lieux de réalisation des travaux :

- Période de réalisation du travail de thèse : 4 ans

- Cadres de coopération : Cotutelle de thèse avec l’université de Littorale Côte

Liste des publications

Liste des communications

* Méthode itérative pour la résolution de l’équation de Riccati, au congrès international

EL MOSTAFA SADEK Thèse de Doctorat

Table des Matières 4

Liste des Tableaux 7

Liste des Figures 8

1 Outils de développement théoriques 13

2 Méthodes de type globale-minimisation pour les équations de Lyapu-

3 L’Extended-Bloc Arnoldi algorithme avec minimisation de résidu pour

4 Méthode de minimisation du résidu pour la résolution de l’équation

5 Équation matricielle de Riccati non symétrique et application à l’équa-

EL MOSTAFA SADEK Thèse de Doctorat

EL MOSTAFA SADEK Thèse de Doctorat

3.1 Résultats pour l’exemple 1.6 . . . . . . . . . . . . . . . . . . . . . . . . . . 71

4.1 Résultats de l’exemple 2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . 93

2.1 Résultats de l’exemple 1 : MR(PCGG) . . . . . . . . . . . . . . . . . . . 42

3.1 Résultat pour l’exemple 1.1. . . . . . . . . . . . . . . . . . . . . . . . . . . 68

4.1 Valeurs singulières de la solution exacte. . . . . . . . . . . . . . . . . . . . 88

5.1 Les valeurs singulières de la solution non négative minimale de l’équation

R le corps des nombres réels.

Les équations matricielles interviennent dans de nombreux domaines des mathéma-

Les équations algébriques matricielles ont été obtenues à partir de la discrétisation

restauration d’images, réduction de modèle en théorie du contrôle, calcul d’un contrôle

Cette thèse se décompose de la façon suivante :

Le deuxième Chapitre portera sur les équations algébriques de Lyapunov AX + XAT +

Dans le Chapitre 3, nous nous intéresserons à la résolution numérique des équations de

Le Chapitre 4 porte sur la résolution d’équations matricielles non linéaires de Riccati

EL MOSTAFA SADEK Thèse de Doctorat

T , où Z une matrice de rang inférieur. Enfin, ce chapitre se

EL MOSTAFA SADEK Thèse de Doctorat

1.2 Matrices particulières

Soit A = (aij ) une matrice carrée de taille n × n.

Définition 1.2.1. La matrice A est dite :

1.3 Norme de Frobenius

Pour deux matrices Y et Z ∈ Rn×s , le produit scalaire de Frobenius de Y et Z est défini

où tr(Y T Z) désigne la trace de la matrice Y T Z. La norme associée est celle de Frobenius

1.4 Formule de Sherman-Morrison Woodbury

Proposition 1.4.1. Soient B une matrice inversible d’ordre n et u, v ∈ Rn . Alors la

1.8 Le produit de Kronecker et le -produit

On rappelle aussi les propriétés suivantes satisfaites par le −produit.

avec Q = [Q1 , ..., Qk ] une matrice F -orthonormale de taille n × ks et qui satisfait QT

On note que QT Z = QT (Q(R ⊗ Is )), donc en utilisant la Proposition 1.8.5, on aura