Algorithme de Bartels-Stewart

En algèbre linéaire numérique, l'algorithme de Bartels-Stewart est un algorithme utilisé pour résoudre numériquement l'équation de Sylvester $AX-XB=C$ .

Historique[modifier | modifier le code]

Développé par Richard Harold Bartels et Gilbert Wright Stewart en 1972^[1], l'algorithme était la première méthode numériquement stable qui pouvait être systématiquement appliquée pour résoudre de telles équations. L'algorithme opère en utilisant les décompositions de Schur réelles de $A$ et $B$ pour transformer l'équation $AX-XB=C$ en un système triangulaire qui peut ensuite être résolu en utilisant la substitution avant ou arrière. En 1979, Gene H. Golub, Charles F. Van Loan et Stephen Nash ont introduit une version améliorée de l'algorithme^[2] connue sous le nom d'algorithme de Hessenberg-Schur. La méthode de Bartels et Stewartt reste une approche standard pour résoudre les équations de Sylvester lorsque $X$ est de taille petite ou moyenne.

L'algorithme[modifier | modifier le code]

Soient $X,C\in \mathbb {R} ^{m\times n}$ , et supposons que les valeurs propres de $A$ sont distinctes des valeurs propres de $B$ . L'équation matricielle $AX-XB=C$ a alors une solution unique. L'algorithme de Bartels-Stewart calcule $X$ en appliquant les étapes suivantes^[1] :

Calculer les décompositions de Schur réelles :
$R=U^{T}AU,$

$S=V^{T}B^{T}V.$

Les matrices $R$ et $S$ sont des matrices triangulaires par blocs, avec des blocs diagonaux de taille 1 ou 2.
Calculer $F=U^{T}CV.$
Résoudre le système simplifié
$RY-YS^{T}=F$ ,

où $Y=U^{T}XV$ . Pour cela, on peut utiliser la substitution avant par blocs. Concrètement, si $s_{k-1,k}=0$ , alors on a :
$(R-s_{k,k}I)y_{k}=f_{k}+\sum _{j=k+1}^{n}s_{k,j}y_{j},$

où $y_{k}$ est la $k$ -ième colonne de $Y$ . Quand $s_{k-1,k}\neq 0$ , les colonnes $y_{k-1}$ et $y_{k}$ doivent être concaténées et résolues simultanément.
Calculer $X=UYV^{T}.$

Coût du calcul[modifier | modifier le code]

En utilisant l'algorithme QR, les décompositions réelles de Schur à l'étape 1 nécessitent environ $10(m^{3}+n^{3})$ opérations flottantes, de sorte que le coût global du calcul global est : $10(m^{3}+n^{3})+2,5(mn^{2}+nm^{2})$ ^[2].

Simplifications et cas particuliers[modifier | modifier le code]

Dans le cas particulier où $B=-A^{T}$ et $C$ est une matrice symétrique, la solution $X$ est également symétrique. Cette symétrie peut être exploitée pour calculer $Y$ plus efficacement à l'étape 3 de l'algorithme^[1].

L'algorithme de Hessenberg-Schur[modifier | modifier le code]

L'algorithme de Hessenberg–Schur^[2] remplace la décomposition $R=U^{T}AU$ de l'étape 1 par la décomposition

H=Q^{T}AQ

,

où $H$ est une matrice de Hessenberg supérieure. Cela conduit à un système de la forme

HY-YS^{T}=F

qui peut être résolu en utilisant la substitution directe. L'avantage de cette approche est que $H=Q^{T}AQ$ peut être calculée en utilisant les réflexions de Householder au prix de $(5/3)m^{3}$ flops, par rapport aux $10m^{3}$ flops nécessaires pour calculer la décomposition réelle de Schur de $A$ .

Logiciel et implémentation[modifier | modifier le code]

Les sous-programmes requis pour la variante Hessenberg-Schur de l'algorithme Bartels-Stewart sont implémentés dans la bibliothèque SLICOT. Ceux-ci sont utilisés dans la boîte à outils du système de contrôle MATLAB.

Approches alternatives[modifier | modifier le code]

Pour les systèmes de grande taille, le coût en ${\mathcal {O}}(m^{3}+n^{3})$ le coût de l'algorithme de Bartels-Stewart peut être prohibitif. Quand $A$ et $B$ sont creuses ou structurées, de sorte que les résolutions linéaires et les multiplications vectorielles matricielles sont efficaces, les algorithmes itératifs peuvent être plus rapides. Ceux-ci incluent des méthodes basées sur la projection, qui utilisent des itérations Méthode itérative de sous-espace de Krylov, des méthodes basées sur l'itération implicite de direction alternée (ADI) et des hybridations qui impliquent à la fois la projection et l'ADI^[3]. Des méthodes itératives peuvent également être utilisées pour construire directement des approximations de rang inférieur à $X$ lors de la résolution $AX-XB=C$ .

Notes et références[modifier | modifier le code]

(en) Cet article est partiellement ou en totalité issu de l’article de Wikipédia en anglais intitulé « Bartels–Stewart algorithm » (voir la liste des auteurs).

↑ ^{a b et c} R. H. Bartels et G. W. Stewart, « Algorithm 432 : Solution of the matrix equation AX + XB = C [F4] », Communications of the ACM, vol. 15, n^o 9,‎ 1972, p. 820–826 (DOI 10.1145/361573.361582).
↑ ^{a b et c} Gene H. Golub, Stephen Nash et Charles F. Van Loan, « A Hessenberg–Schur method for the problem AX + XB= C », IEEE Transactions on Automatic Control, vol. 24, n^o 6,‎ 1979, p. 909–913 (DOI 10.1109/TAC.1979.1102170, hdl 1813/7472).
↑ Valeria Simoncini, « Computational Methods for Linear Matrix Equations », SIAM Review, vol. 58, n^o 3,‎ 2016, p. 377–441 (DOI 10.1137/130912839, hdl 11585/586011, S2CID 17271167).

Portail des mathématiques

[:0-1] {a b et c} R. H. Bartels et G. W. Stewart, « Algorithm 432 : Solution of the matrix equation AX + XB = C [F4] », Communications of the ACM, vol. 15, n^o 9,‎ 1972, p. 820–826 (DOI 10.1145/361573.361582).

[:1-2] {a b et c} Gene H. Golub, Stephen Nash et Charles F. Van Loan, « A Hessenberg–Schur method for the problem AX + XB= C », IEEE Transactions on Automatic Control, vol. 24, n^o 6,‎ 1979, p. 909–913 (DOI 10.1109/TAC.1979.1102170, hdl 1813/7472).

[3] Valeria Simoncini, « Computational Methods for Linear Matrix Equations », SIAM Review, vol. 58, n^o 3,‎ 2016, p. 377–441 (DOI 10.1137/130912839, hdl 11585/586011, S2CID 17271167).

[1]

[2]

[3]