r6.03.01 Résolution de systèmes non réguliers par une méthode de décomposition en valeurs singulières

Résumé :

Ce document est consacré à la résolution des systèmes d’équations linéaires non réguliers. Les matrices prises en compte peuvent être carrées non inversibles ou rectangulaires.

Après avoir rappelé le cadre théorique des solutions au sens des moindres carrés, nous concentrons l’exposé sur la méthode par décomposition en valeurs singulières qui fournit, d’une part, un outil de diagnostic du degré de régularité du système, et, d’autre part, une famille d’algorithmes de résolution à la fois plus généraux et plus stables que ceux dérivant de l’approche par les équations normales.

Enfin, nous détaillons l’algorithme mis en œuvre dans le Code_Aster qui résorbe la fonctionnalité équivalente de la librairie Nag (F04JDF pour la version 12 et F04JDE pour la version 15) utilisée pour la modélisation du comportement métallurgique des aciers [R4.04.01].

Solution d’un système linéaire rectangulaire

Dans ce paragraphe nous allons définir une notion de solution pour le système linéaire [éq 1-1] qui jouit des propriétés d’ existence et d’ unicité. La démarche procède en deux temps:

• D’abord, par une approche du type moindres carrés nous construisons un problème d’optimisation différentiable et convexe (section 2.1) qui admet toujours au moins une solution (section 2.2). La situation 2) du paragraphe 1 est alors éliminée,

• Puis, analysant la propriété d’unicité (section 2.3) pour constater qu’elle n’est pas toujours garantie nous imposerons une contrainte supplémentaire (section 2.4) à la solution caractérisée dans la section 2.1 de façon à rétablir l’unicité.

Formalisme des moindres carrés

L’unique solution d’un système linéaire \(\text{Ax}=b\) de matrice carrée et régulière réalise le minimum de la quantité \(\parallel \text{Ay}-b\parallel\) quand \(y\) décrit \({ℝ}^{n}\) . Cette propriété nous ouvre la voie qui conduit à une notion de solution pour un système linéaire général du type [éq1-1] qui lui confère les mêmes propriétés que celles du cas particulier du système régulier. Nous dirons donc d’un point \(x\) de \({ℝ}^{n}\) qu’il est solution du système [éq1-1] s’il est solution du problème d’optimisation :

\(\parallel \text{Ax}-b\parallel =\underset{y\in {ℝ}^{n}}{\text{Min}}\parallel \text{Ay}-b\parallel\) éq 2.1-1

Cette approche est naturelle car elle définit une solution dont le résidu \(r=\text{Ax}-b\) est nul dans le cas où le second membre est élément de \(\text{Im}A\) et est de norme minimale dans le cas contraire, ce qui constitue le mieux qu’on puisse attendre.

Pour analyser le problème [éq 2.1-1], il est commode de lui substituer le problème d’optimisation sans contraintes équivalent suivant:

\(\text{trouver}x\in {ℝ}^{n}\text{tel que}J(x)=\underset{y\in {ℝ}^{n}}{\text{MinJ}}(y)\) éq 2.1-2

où \(J(.)\) est la fonctionnelle définie par:

\(J:y\in {ℝ}^{n}\to J(y)=\frac{1}{2}\parallel \text{Ax}-b\parallel\)

L’intérêt du problème [éq 2.1-2] tient au fait que la fonctionnelle \(J(.)\) vérifie les propriétés suivantes:

• \(J(.)\) est deux fois continument différentiable:

\(\text{DJ}(x):h\in {ℝ}^{n}\to \text{DJ}(x)h=({A}^{T}\text{Ax}-{A}^{T}\text{b,h})\in ℝ\) éq 2.1-3

\({D}^{2}J(x):(\text{h,k})\in {ℝ}^{n}\times {ℝ}^{n}\to \text{DJ}(x)(\text{h,k})=({A}^{T}\text{Ah},k)\in ℝ\) éq 2.1-4

• \(J(.)\) est quadratique et convexe.

Ainsi, le problème [éq 2.1-2] s’inscrit dans le cadre de l’optimisation différentiable et convexe de sorte que nous disposons des résultats suivants ([bib1] p. 156 et 146):

1. De la convexité: tout optimum local est en fait un optimum global, c’est à dire une solution de [éq2.1-2],

2. De la différentiabilité: tout optimum local vérifie l’équation d’Euler \(\text{DJ}(x)=0\) sur Än qui, compte‑tenu de [éq 2.1-3], conduit à la caractérisation par les équations dites normales :

\({A}^{T}\text{Ax}={A}^{T}b\) éq 2.1-5

Existence d’optimums

Dans [bib1] p. 171 on trouve une démonstration de l’existence d’au moins une solution aux équations normales [éq 2.1-5]. Cette démonstration s’appuie sur des arguments destinés à la prise en compte du cas de la dimension infinie (théorème de projection sur un convexe fermé d’un espace de Hilbert).

Notre cas étant nettement plus simple, nous donnons une démonstration du résultat qui n’utilise que des arguments algébriques simples qui, de plus, nous serons utiles dans le paragraphe 3. Montrer que, pour tout \(b\) élément de \({ℝ}^{m}\) , les équations normales [éq 2.1-5] admettent une solution équivaut à l’établissement de l’inclusion \(\text{Im}{A}^{T}\subset \text{Im}{A}^{T}A\) . Or, pour toute matrice réelle \(M\) d’ordre \(m\times n\) nous avons \(\text{Im}{M}^{T}={(\text{Ker}M)}^{\perp }\) ([bib3] p. 28). Aussi, l’inclusion à établir qui équivaut à \({(\text{Ker}A)}^{\perp }\subset {(\text{Ker}{A}^{T}A)}^{\perp }\) qui est elle même équivalente à \(\text{Ker}{A}^{T}A\subset \text{Ker}A\) . Soit donc \(x\in \text{Ker}{A}^{T}A\) ; alors \(\text{Ax}\in \text{Ker}{A}^{T}\) , c’est-à-dire \(\text{Ax}\in {(\text{Im}A)}^{\perp }\) . Comme \(\text{Ax}\) est aussi élément de \(\text{Im}A\) , il ne peut être que nul ce qui signifie que \(x\in \text{Ker}A\) et achève la démonstration.

A ce stade du propos, nous pouvons dire que tout système du type [éq 1-1] admet au moins une solution au sens de [éq 2.1-3] et toutes ces solutions sont caractérisées comme solution au sens de Cramer des équations normales de [éq 2.1-5]. La situation 2) du paragraphe 1 est éliminée.

Reste à éliminer la situation 3), c’est-à-dire à garantir l’unicité.

Unicité de l’optimum et rang du système

Il est clair que les équations normales [éq 2.1-5], caractérisant les optimums que nous cherchons, admettent une unique solution sous la condition nécessaire et suffisante que \({A}^{T}A\) soit régulière. Comme \({A}^{T}A\) est toujours semi-définie positive, son inversibilité équivaut à sa définie positivité, de sorte que, compte tenu de l’expression [éq 2.1-4] de la dérivée seconde de la fonctionnelle \(J(.)\) , nous retrouvons le théorème bien connu d’unicité de l’optimum du problème [éq2.1-2] pour une fonctionnelle convexe deux fois continument différentiable ([bib1] th 7.4-3 et 7.4-4).

En toute généralité, rien n’empêche la matrice \({A}^{T}A\) d’être singulière, la solution du système [éq1‑1] au sens de [éq 2.1-2] n’est donc pas toujours unique. Nous disposons néanmoins d’un critère pour détecter cette situation. A la section 2.2 nous avons établi que \(\text{Im}{A}^{T}\subset \text{Im}{A}^{T}A\) et comme l’inclusion réciproque est trivialement vraie, nous pouvons conclure à l’identité \(\text{Im}{A}^{T}=\text{Im}{A}^{T}A\) . L’introduction du rang \(\text{rg}(A)\) de la matrice \(A\) , la dimension de son espace image, nous permet alors de dire qu’une condition nécessaire et suffisante pour que \({A}^{T}A\) soit inversible est que \(\text{rg}({A}^{T}A)=n\) ce qui équivaut à \(\text{rg}(A)=n\) car \(\text{rg}({A}^{T}A)=\text{rg}({A}^{T})=\text{rg}(A)\) .

L’intérêt de ce critère tient au fait qu’il limite l’analyse à la seule matrice \(A\) sans qu’il soit nécessaire de former explicitement \({A}^{T}A\) . Ce critère nous montre aussi que les équations normales associées à un système linéaire strictement sous-contraint admettent toujours une infinité de solutions. En effet, le rang d’une matrice est aussi égal au nombre de colonnes indépendantes qu’elle possède; aussi, pour que ce rang atteigne la valeur \(n\) il est nécessaire que les colonnes de la matrice soit d’ordre au moins \(n\) .

Solution au sens des moindres carrés

Nous venons de constater que l’ensemble des points qui minimisent le résidu du système [éq 1-1] n’est pas nécessairement réduit à un seul point. Pour rétablir l’unicité nous affinons la notion de solution du système [éq 1-1] de la section 2.1 en définissant la solution au sens des moindres carrés comme l’élément de norme minimale de l’ensemble des points qui minimisent le résidu. Cette solution \(x\) est alors caractérisée par:

\(x\in {S}^{\text{def}}=\left\lbrace y\in {ℝ}^{n};{A}^{T}\text{Ax}={A}^{T}b\right\rbrace \text{et}\parallel x\parallel =\underset{y\in S}{\text{Inf}}\parallel y\parallel\)

Cette caractérisation n’est pas satisfaisante sur le plan pratique car elle demande la résolution d’un problème d’optimisation sous contraintes. Nous allons lui substituer une autre caractérisation plus adaptée au sens où elle conduira (voir la section 4.2) à une procédure de calcul nettement plus simple.

L’ensemble \(S\) ci-dessus est le translaté de noyau de \({A}^{T}A\) par l’un quelconque des vecteurs solutions des équations normales [éq 2.1-5]. Aussi, la condition supplémentaire de minimisation de la norme s’interprète comme une simple projection: la solution au sens des moindres carrés du système [éq 1-1] n’est rien d’autre que la projection de l’origine de Än sur l’ensemble des solutions des équations normales. Aussi, nous pouvons la caractériser comme le point d’intersection entre l’ensemble \(S\) et l’orthogonal du noyau de \({A}^{T}A\) .

La définition d’une solution au système [éq 1-1] peut alors être résumée comme suit:

\(x\text{est solution de}\text{Ax}=b\iff \lbrace \begin{array}{}{A}^{T}\text{Ax}={A}^{T}b\\ x\in {(\text{Ker}{A}^{T}A)}^{\perp }\end{array}\) éq 2.4-1

La première condition fait de \(x\) un vecteur de résidu minimal tandis que la deuxième sélectionne, parmi les vecteurs de résidu minimal, celui de norme minimale.

La définition [éq 2.4-1] est une généralisation classique de la notion de solution d’un système équi‑contraint régulier et confère à tout système du type [éq 1-1] une solution et une seule.

Valeurs singulières

Dans ce paragraphe nous présentons quelques résultats utiles pour la conception d’une méthode de résolution opérationnelle du système [éq 1-1]. Ces résultats dérivent de la notion de valeurs singulières (section 3.1) et permettent de construire une base du noyau et une base de l’image de la matrice du système (section 3.2) à partir desquelles il est possible de donner un sens, adapté au calcul de la solution au sens des moindres carrés, à l’inverse d’une matrice quelconque (section 3.3).

Décomposition en valeurs singulières

Commençons par rappeler la définition des valeurs singulières. On appelle valeurs singulières d’une matrice réelle \(A\) d’ordre \(m\times n\) les racines carrées des valeurs propres de la matrice carrée \({A}^{T}A\) d’ordre \(n\) qui, rappelons-le, est semi-définie positive.

La notion de diagonalisation des matrices carrées (lorsqu’elles sont diagonalisables) se généralise aux matrices rectangulaires (sans restriction) par le concept de décomposition (ou factorisation) en valeurs singulières.

Pour toutes matrices réelles \(A\) d’ordre \(m\times n\) , il existe deux matrices carrées unitaires \(Q\) et \(P\) d’ordre respectif \(m\) et \(n\) telles que:

\({A=Q\Sigma P}^{T}\) éq 3.1-1

où \(\Sigma\) est une matrice d’ordre \(m\times n\) dont la structure est schématisée ci-dessous:

\(\Sigma =∣\begin{array}{cccc}{\mu}_{1}& & & \\ & {\mu}_{2}& & \\ & & \mathrm{...}& \\ & & & {\mu}_{n}\end{array}∣\begin{array}{cccc}& & & \\ & 0& & \\ & & & \\ & & & \end{array}\mid \text{si}m\le n\)

\(\Sigma =∣\begin{array}{c}\underline{\begin{array}{cccc}{\mu}_{1}& & & \\ & {\mu}_{2}& & \\ & & \mathrm{...}& \\ & & & {\mu}_{n}\end{array}}\\ 0\end{array}∣\text{si}m>n\)

Les \({\mu}_{i}\) sont les valeurs singulières de \(A\) que nous supposons ordonnées par ordre décroissant :

\({\mu}_{1}\ge {\mu}_{2}\ge \text{...}\ge {\mu}_{n}\)

On peut trouver une démonstration de ce résultat dans [bib1] p.10 pour le cas équi-contraint et dans [bib3] p.73 pour le cas sur-contraint, le cas sous-contraint s’en déduit alors par transposition.

La factorisation SVD [éq 3.1-1] de \(A\) donne \({A}^{T}A={P\Sigma }^{T}{\mathrm{SP}}^{T}\) et \({\text{AA}}^{T}={Q\Sigma \Sigma }^{T}{Q}^{T}\) de sorte que, \({\Sigma}^{T}\Sigma\) et \({\Sigma \Sigma }^{T}\) étant des matrices carrées diagonales, la matrice \(P\) est constituée des vecteurs propres orthonormalisés de la matrice \({A}^{T}A\) tandis que la matrice \(Q\) est constituée des vecteurs propres orthonormalisés de la matrice \({\text{AA}}^{T}\) .

Rang, image et noyau

Le paragraphe [§2] a montré le rôle fondamental que joue le rang de la matrice \(A\) et le noyau de la matrice \({A}^{T}A\) pour la résolution d’un système linéaire non régulier du type [éq 1-1]. Nous allons voir maintenant comment la factorisation [éq 3.1-1] peut être utilisée pour déterminer ce rang ainsi qu’une base de \(\text{Ker}{A}^{T}A\) .

Soit \(r\) l’indice de la plus petite valeur singulière non nulle. La factorisation [éq 3.1-1] s’écrit aussi \({Q}^{T}\text{AP}=\Sigma\) où la prise en compte des valeurs singulières nulles permet de préciser la décomposition en bloc de \(\Sigma\) :

où \({\Sigma}_{r}=\text{Diag}({\mu}_{1},{\mu}_{2},\text{...},{\mu}_{r})\) est la matrice diagonale d’ordre \(r\) des valeurs singulières non nulles dans l’ordre croissant.

Puisque les matrices \(Q\) et \(P\) sont régulières, les matrices \(A\) et \(\Sigma\) sont équivalentes de sorte que leur noyau et image respectifs coïncident. Nous en déduisons donc que:

• Le rang de \(A\) coïncide avec le nombre de valeurs singulières non nulles:

\(\text{rg}A=r\)

• Les vecteurs colonnes de \(P\) d’indice \(r+1\) à \(n\) forment une base de \(\text{Ker}A\)

• Les vecteurs colonnes de \(Q\) correspondant aux valeurs singulières non nulles forment une base de \(\text{Im}A\)

D’autre part, à la section 2.3 nous avons vu que \(\text{Im}{A}^{T}=\text{Im}{A}^{T}A\) . L’identité \(\text{Im}{M}^{T}={(\text{Ker}M)}^{\perp }\) nous donne alors \(\text{Ker}A=\text{Ker}{A}^{T}A\) de sorte que la deuxième condition de la définition [éq 2.4-1] est simplement réalisée par tout vecteur qui s’exprime comme une combinaison linéaire des vecteurs colonnes de \(P\) correspondant aux valeurs singulières non nulles.

Pseudo-inverse et solution au sens des moindres carrés

Une autre application de la décomposition en valeurs singulières consiste en la notion de pseudo‑inverse (ou inverse Moore-Penrose) qui généralise la notion habituelle d’inverse d’une matrice carrée régulière aux matrices rectangulaires d’une part, et aux matrices carrées singulières d’autre part.

Tout d’abord, la pseudo-inverse d’une matrice \(\Sigma\) de la décomposition en valeurs singulières [éq3.1-1] est définie par:

où \({\Sigma}_{r}^{-1}=\text{Diag}(\frac{1}{{\mu}_{1}},\frac{1}{{\mu}_{2}},\text{...},\frac{1}{{\mu}_{r}})\) est l’inverse au sens habituel de \({S}_{r}\) .

Ceci étant, nous utilisons la décomposition [éq 3.1-1] de la matrice \(A\) pour définir sa pseudo‑inverse \({A}^{+}\) par:

\({A}^{+}=P{\Sigma}^{+}{Q}^{T}\) éq 3.3-1

De même, de la décomposition [éq 3.1-1] de la matrice \(A\) nous tirons \({A}^{T}A=P{\Sigma}^{T}\Sigma {P}^{T}\) , de sorte que la pseudo-inverse \({({A}^{T}A)}^{+}\) de la matrice \({A}^{T}A\) est définie par:

\({({A}^{T}A)}^{+}=P{\Sigma}^{+}{({\Sigma}^{T})}^{+}{P}^{T}\) éq 3.3-2

Nous sommes maintenant en mesure de fournir une interprétation simple de la solution au sens des moindres carrés définie par [éq 2.4-1].

La restriction à \({(\text{Ker}{A}^{T}A)}^{\perp }\) de l’application linéaire associée à la matrice \({A}^{T}A\) définit un isomorphisme de \({(\text{Ker}{A}^{T}A)}^{\perp }\) sur \(\text{Im}{A}^{T}A\) . Comme, d’une part \({(\text{Ker}{A}^{T}A)}^{\perp }={(\text{Ker}A)}^{\perp }=\text{Im}{A}^{T}\) , et, d’autre part, \(\text{Im}{A}^{T}A=\text{Im}{A}^{T}\) , cette restriction est en fait un automorphisme de \({(\text{Ker}{A}^{T}A)}^{\perp }\) .

Dans la base de \({(\text{Ker}{A}^{T}A)}^{\perp }\) constituée par les \(r\) premières colonnes de la matrice \(P\) , cet automorphisme est représenté par la matrice \({\Sigma}_{r}^{2}\) . Aussi, son automorphisme réciproque \(y\) est représenté par la matrice \({\Sigma}_{r}^{-2}\) . L’extension à \({ℝ}^{n}\) de cet automorphisme est alors représentée, dans la base associée à la matrice \(P\) , par la matrice \({({\Sigma}^{T}\Sigma )}^{+}={\Sigma}^{T}{({\Sigma}^{T})}^{+}\) , et donc, dans la base canonique, par la matrice \({({A}^{T}A)}^{+}\) .

Il suit que:

• Nous retrouvons le fait que, pour tout \(b\) élément de \({ℝ}^{n}\) , il existe un unique vecteur \(x\in {(\text{Ker}{A}^{T}A)}^{\perp }\) solution de \({A}^{T}\text{Ax}={A}^{T}b\) , soit l’existence et l’unicité de la solution au sens des moindres carrés [éq 2.4-1) du système \(\text{Ax}=b\) ,

• Cette unique solution est donnée par:

\(x={({A}^{T}A)}^{+}{A}^{T}b\) éq 3.3-3

La pseudo-inverse d’une matrice est définie à partir de la décomposition SVD de cette matrice. Comme la décomposition SVD n’est pas unique, la matrice pseudo-inverse n’est pas unique. Par contre, du point de vue des applications linéaires associées aux matrices, l’application pseudo‑inverse est unique. Toutes les matrices pseudo-inverses associées aux différentes décompositions SVD d’une matrice donnée ne sont alors que des représentantes matricielles particulières qui expriment cette application pseudo-inverse relativement aux bases induites par les matrices orthogonales des décompositions SVD. Aussi, l’expression [éq 3.3-3] a un sens: elle définit un vecteur dont \(\mathrm{x}\) représente les composantes par rapport à la base d’arrivée (matrice \(P\) ) de la décomposition SVD.

Résolution d’un système linéaire rectangulaire

Les deux méthodes de résolution du système [éq 1-1] que nous présentons aux sections 4.1 (méthode des équations normales) et 4.2 (décomposition en valeurs singulières) visent à la résolution des équations normales [éq 2.1-5]. Ces deux méthodes se distinguent non seulement par le choix des algorithmes qu’elles mettent en œuvre (inversion contre pseudo-inversion), mais aussi par leur degré de généralité et par leurs propriétés numériques qui sont comparées à la section 4.3

Méthode des équations normales

La résolution du système \(\text{Ax}=b\) par la méthode des équations normales consiste à calculer la solution au sens des moindres carrés [éq 2.4-1] de façon “directe”, c’est-à-dire en utilisant directement la relation \(x={({A}^{T}A)}^{-1}{A}^{T}b\) . Pour cela, il s’agit d’abord de calculer \({A}^{T}A\) et \({A}^{T}b\) , puis de résoudre le système obtenu soit par méthode itérative soit par une factorisation de \({A}^{T}A\) .

Nous pouvons d’ores et déjà remarquer que cette méthode est limitée aux matrices \({A}^{T}A\) régulières, ce qui limite son domaine d’application au système [éq 1-1] dont la matrice est de plein rang (voir la section 2.3). En particulier, la méthode des équations normales ne peut traiter ni les systèmes strictement sous-contraints, ni les systèmes équi-contraints singuliers (voir la section 2.4).

Méthode par décomposition en valeurs singulières

Nous avons vu à la section 3.3 que la solution du système \(\text{Ax}=b\) au sens des moindres carrés définie par [éq 2.4-1] peut être caractérisée par la relation \(\mathrm{x}={({\mathrm{A}}^{\mathrm{T}}\mathrm{A})}^{+}{\mathrm{A}}^{\mathrm{T}}\mathrm{b}\) [éq 3.3-3]. La méthode de résolution du système basée sur cette propriété est dite méthode par décomposition en valeurs singulières car elle construit la pseudo-inverse [éq 3.3-2] de \({A}^{T}A\) via la décomposition SVD [éq3.1‑1] de la matrice \(A\) .

Comme toute matrice peut être décomposée en valeurs singulières, il suit que tout système du type [éq1-1] peut être résolu au sens de [éq 2.4-1] par cette méthode qui présente donc, au moins, l’avantage de la généralité par rapport à la méthode des équations normales.

Ce n’est pas tout. La méthode par décomposition en valeurs singulières, contrairement à la méthode des équations normales, ne demande pas la construction explicite de la matrice \({A}^{T}A\) et du vecteur \({A}^{T}b\) (nous verrons à la section 4.3 l’intérêt sur le plan numérique de cette propriété). En effet, il est aisé de vérifier que la matrice \(\Sigma\) des valeurs singulières de la factorisation [éq 3.1-1] satisfait à l’identité \({\Sigma}^{+}{({\Sigma}^{T})}^{+}{\Sigma}^{T}={\Sigma}^{+}\) , de sorte que, prémultipliant \({A}^{T}\) par la pseudo-inverse de \({A}^{T}A\) et tenant compte de la factorisation [éq 3.1-1], nous obtenons \({({A}^{T}A)}^{+}{A}^{T}=P{\Sigma}^{+}{({\Sigma}^{T})}^{+}{\Sigma}^{T}{P}^{T}P{\Sigma}^{T}{Q}^{T}\) , qui, par orthogonalité de \(P\) nous donne \({({A}^{T}A)}^{+}{A}^{T}={A}^{+}\) . Dès lors, les caractérisations [éq 3.3-3] et [éq 2.4-1] de la solution cherchée sont équivalentes à la caractérisation:

\(x\text{est solution de}\text{Ax}=b\iff x={A}^{+}b\) éq 4.2-1

Comparaison de la méthode des équations normales à la méthode de décomposition en valeurs singulières

Dans les deux sections précédentes, nous venons de constater, qu’algébriquement, la méthode de décomposition en valeurs singulières est plus générale et plus simple que la méthode des équations normales. Nous allons maintenant constater, en suivant [bib3] p. 336, qu’elle lui est aussi supérieure sur le plan numérique. Cette supériorité s’exprime d’une part, en terme de stabilité non seulement de la résolution, mais aussi de la construction du problème et, d’autre part, à un niveau moins critique, en terme d’adaptation au traitement des matrices creuses.

Conditionnement

Le conditionnement d’une matrice \(A\) d’ordre \(m\times n\) est défini comme le rapport de ses valeurs singulières extrêmes et non nulles:

\(\text{cond}(A)=\frac{{\mu}_{1}}{{\mu}_{r}}\)

où \(r\) est le rang de la matrice \(A\) .

Les résultats présentés dans [bib3] p.184, utilisant les équations normales comme un outil d’analyse et non comme un outil de calcul, montrent que la perturbation de la solution du problème d’optimisation [éq 2.1-1] due aux erreurs d’arrondis peut être proportionnelle à \(\text{cond}{(A)}^{2}\) . Mais les résultats classiques de l’analyse de stabilité de la solution d’un système linéaire par rapport à ces mêmes erreurs montrent une proportionnalité au nombre de conditionnement de la matrice. Si bien que dans le cas d’une résolution directe des équations normales, nous obtenons une erreur toujours proportionnelle à \(\text{cond}({A}^{T}A)=\text{cond}{(A)}^{2}\) , ce qui est moins bon que \(\text{cond}(A)\) .

La méthode de résolution par décomposition en valeurs singulières n’utilise que des transformations orthogonales (voir le paragraphe 4), si bien qu’elle ne modifie pas le conditionnement initial du problème [éq 1-1] et est donc, de ce point de vue, plus attrayante que la méthode des équations normales.

Perte de précision

Nous venons de voir que les erreurs d’arrondis conduisent à une dégradation de la solution plus sensible lorsqu’elle est calculée via les équations normales plutôt que par une décomposition en valeurs singulières. L’exemple suivant, tiré de [bib 2], montre que la construction même du système [éq2.1-5] des équations normales est perturbée par les erreurs d’arrondi.

Soit donc la matrice suivante:

\(A=\left[\begin{array}{cc}1& 1\\ \varepsilon & 0\\ 0& \varepsilon \end{array}\right]\) conduisant à \({A}^{T}A=\left[\begin{array}{cc}1+{\varepsilon}^{2}& 1\\ 1& 1+{\varepsilon}^{2}\end{array}\right]\)

dont les valeurs singulières sont \({\mu}_{1}=\sqrt{2+{\varepsilon}^{2}}\) et \({\mu}_{2}=\mid \varepsilon \mid\) , de sorte que le rang de \(A\) est \(2\) dès que \(\varepsilon \ne 0\) . Si \(\varepsilon\) vérifie \({\varepsilon}^{2}<{\varepsilon}_{\text{mach}}<\varepsilon\) où \({\varepsilon}_{\text{mach}}\) est la précision machine, alors tous les coefficients de \({A}^{T}A\) seront calculés à la valeur \(1\) et les valeurs singulières calculées seront, au mieux, \({\mu}_{1}=\sqrt{2}\) et \({\mu}_{2}=0\) . Il suit que le rang numérique, calculé par les équations normales, sera \(1\) , alors que celui calculé par une décomposition SVD de la matrice \(A\) serait égal à \(2\)

Structure creuse

A un niveau moindre, la construction de la matrice des équations normales induit un remplissage du système associé que la méthode utilisant la décomposition en valeurs singulières évite.

Conclusion

Le tableau suivant résume la discussion des sous-sections précédentes:

	Généralité	Conditionnement	Perte de précision à la construction du problème	remplissage
Équations normales	systèmes de plein rang	\(\text{cond}{(A)}^{2}\)	possible	oui
SVD	tout système	\(\text{cond}(A)\)	impossible	non

Algorithme SVD pour la résolution d’un système linéaire équi ou sous-contraint

Dans ce paragraphe, nous détaillons la méthode de résolution des systèmes non-réguliers mise en œuvre dans le Code_Aster . Cette méthode s’applique au système sous-contraint ou équi-contraint singulier et fournit la solution au sens des moindres carrés [éq 2.4-1].

Le calcul d’une décomposition SVD de \(A\) est équivalent au calcul du spectre de la matrice normale associée \({A}^{T}A\) . Aussi, il ne peut être obtenu qu’à la convergence d’un procédé itératif.

La section 5.1 expose le principe de l’algorithme et montre en particulier comment l’application de deux transformations orthogonales permet de réduire le problème à la simple recherche de la décomposition SVD d’une matrice bi-diagonale supérieure. Les sections 5.2 et 5.3 sont consacrées à l’algorithmique de ces réductions . La section 5.4 présente l’algorithme de la décomposition SVD de la matrice bi-diagonale.

Les algorithmes seront décrits avec la convention de notation dans laquelle:

• \(R(i,j,\theta )\) désigne la rotation de Givens du plan \((i,j)\) et d’angle \(\theta\) ,

• \({A}^{(k)}\) désigne l’itéré d’indice \(k\) d’une itération matricielle et \({A}^{(k,l)}\) l’itéré \(l\) d’une itération interne à l’itéré \({A}^{(k)}\) .

Réduction du problème et principe de l’algorithme

Dans cette section, nous présentons l’algorithme de résolution d’un sytème linéaire équi ou sous‑contraint par la méthode SVD.

Nous réduisons le problème à la recherche de la décomposition SVD d’une matrice bidiagonale comme dans [bib2] mais nous effectuons la réduction d’une autre façon que celle proposée dans [bib2]: nous commençons par réduire la matrice à une forme triangulaire supérieure, puis, nous réduisons cette triangulaire à une forme bidiagonale supérieure. Ces deux réductions sont effectuées par transformations orthogonales.

Les opérations de calcul de la décomposition SVD s’enchaînent comme suit:

Réduction à la forme triangulaire supérieure

\(\begin{array}{ccc}A=[U\text{}0]{P}_{1}^{T}& & \text{si}m<n\\ A={\mathrm{UP}}_{1}^{T}& & \text{si}m=n\end{array}\) éq 5.1-1

où \({P}_{1}\) est une matrice orthogonale d’ordre \(n\) et \(U\) une matrice triangulaire supérieure d’ordre \(m\) .

Réduction à la forme bi-diagonale supérieure

\(U={Q}_{2}{\text{BP}}_{2}^{T}\) éq 5.1-2

où \({Q}_{2}\) et \({P}_{2}\) sont deux matrices orthogonales d’ordre \(m\) et \(B\) une matrice bidiagonale supérieure d’ordre \(m\) .

Décomposition SVD de la bi-diagonale supérieure

\(B={Q}_{3}\Sigma {P}_{3}^{T}\) éq 5.1-3

où \({Q}_{3}\) et \({P}_{3}\) sont deux matrices orthogonales d’ordre \(m\) et \(\Sigma\) une matrice diagonale d’ordre \(m\) de la forme:

Combinant les relations [éq 5.1-1], [éq 5.1-2] et [éq 5.1-3], nous obtenons une décomposition SVD de la matrice \(A\) :

éq 5.1-4

La solution au sens des moindres carrés [éq 2.4-1] du système [éq 1-1] est alors obtenue par l’application de la pseudo-inverse [éq 3.3-1] de \(A\) déduite de la décomposition en valeurs singulières [éq 5.1-4]. Nous obtenons donc:

\(x={P}_{1}\left[\begin{array}{c}{P}_{2}{P}_{3}{S}^{+}{Q}_{3}^{T}{Q}_{2}^{T}b\\ 0\end{array}\right]\) éq 5.1-5

L’algorithme proposé consiste donc en l’enchaînement des factorisations [éq 5.1-1], [éq 5.1-2] et [éq5.1-3] préalablement à l’application de la relation [éq 5.1-5].

Réduction à la forme triangulaire supérieure

A partir d’une matrice \(A\) d’ordre \(m\times n\) pour \(m\le n\) on détermine une matrice triangulaire supérieure \(U\) d’ordre \(m\) et une matrice orthogonale \(P\) d’ordre \(n\) telles que:

\(\begin{array}{ccc}A=\left[U\text{}0\right]{P}^{T}& & \text{si}m<n\\ & & \\ A={\mathrm{UP}}^{T}& & \text{si}m=n\end{array}\)

Algorithmiquement, la factorisation utilise une méthode d’élimination qui s’interprète, comme la construction d’une suite de matrices \({A}^{(k)}\) par:

\(\lbrace \begin{array}{ccc}{A}^{(m)}=A& & \\ & & \\ {A}^{(k-1)}={A}^{(k)}{P}^{(k)}& & \text{pour}k=m,m-1,\text{...},1\end{array}\)

où chaque matrice courante \({A}^{(k)}\) présente la structure schématisée ci-dessous:

Les coefficients \({a}_{i,j}^{(k)}\) des matrices \({A}^{(k)}\) de l’itération vérifient donc:

\({a}_{i,j}^{(k)}=0\text{si}\lbrace \begin{array}{cc}k+1\le i\le m& \text{et}m+1\le j\le n\\ k+1\le i\le m& \text{et}1\le j\le k\\ k+1\le j\le i\le m& \end{array}\) éq 5.2-1

de sorte qu’à l’issue de la récurrence, nous aurons:

\(U={A}^{(1)}\text{et}P=\underset{k=m,m-1,\text{...},1}{\Pi {P}^{(k)}}\)

Le problème se réduit donc à la préservation de la structure [éq 5.2-1] lors du passage de \({A}^{(k)}\) à \({A}^{(k-1)}\) par une transformation \({P}^{(k)}\) qui doit être orthogonale. Le problème de l’orthogonalité est réglé en choisissant la transformation comme un produit de rotations de Givens et le problème de la préservation de la structure est résolu en effectuant ce produit dans un ordre qui ne détruit pas les zéros créés.

Tenant compte de la structure rectangulaire de la matrice \({A}^{(k)}\) , nous construisons l’itéré en \({A}^{(k-1)}\) deux phases:

• La phase \(1\) annule successivement les coefficients \({a}_{k,j}^{(k-1)}\) correspondant aux colonnes \(j=k-1,k-2,\text{...},1,\) , ce qui se traduit par:

\(\begin{array}{cc}{A}^{(k-1,k-1)}={A}^{(k)}& \\ {A}^{(k-1,j-1)}={A}^{(k-1,j)}R{(k,j,{\theta}_{j}^{(k)})}^{T}& \text{pour}j=k-1,k-2,\text{...},1\end{array}\)

• La phase \(2\) annule successivement les coefficients \({a}_{k,j}^{(k-1)}\) correspondant aux colonnes \(j=n,n-1,\text{...},m+1,\) , ce qui se traduit par la récurrence:

\(\begin{array}{cc}{A}^{(k-1,n)}={A}^{(k-1,0)}& \\ {A}^{(k-1,j-1)}={A}^{(k-1,j)}R{(k,j,{\theta}_{j}^{(k)})}^{T}& \text{pour}j=n,n-1,\text{...},k+1\end{array}\)

L’angle \({\theta}_{j}^{(k)}\) de la rotation de Givens du plan \((k,j)\) est choisi pour annuler le coefficient en position \((k,j)\) de \({a}^{(k-1,j)}\) . L’application de chaque rotation ne modifie donc que les colonnes \(k\) et \(j\) ce qui ne détruit pas les coefficients nuls produits par les étapes précédentes. Nous constatons que la colonne \(k\) joue un rôle particulier (celui de pivot) car elle seule est systématiquement modifiée par chaque rotation alors que les autres colonnes ne sont modifiées que par la rotation qui annule leur coefficient à la ligne \(k\) .

A l’issue de ces récurrences, nous avons \({A}^{(k-1)}={A}^{(k-1,k)}\) . La matrice \({P}^{(k)}\) est alors donnée par:

\({P}^{(k)}=\underset{j=m+1}{\overset{j=n}{\Pi}}R(k,j,{\theta}_{j}^{(k)})\underset{j=1}{\overset{j=k-1}{\Pi}}R(k,j,{\theta}_{j}^{(k)})\)

de sorte que la matrice \(P\) vaut:

\(P=\underset{k=1}{\overset{k=m}{\Pi}}(\underset{j=m+1}{\overset{j=n}{\Pi}}R(k,j,{\theta}_{j}^{(k)})\underset{j=1}{\overset{j=k-1}{\Pi}}R(k,j,{\theta}_{j}^{(k)}))\)

Réduction à la forme bi-diagonale supérieure

Réduire une matrice carrée triangulaire supérieure \(A\) d’ordre \(m\) à la forme bidiagonale supérieure consiste à trouver deux matrices orthogonales \(P\) et \(Q\) et une matrice bidiagonale supérieure \(B\) , toutes trois d’ordre \(m\) , telles que:

\(A={\text{QBP}}^{T}\)

Algorithmiquement, la factorisation procède comme celle de la section précédente en utilisant une méthode d’élimination qui s’interprète algébriquement comme la construction d’une suite de matrices \({A}^{(k)}\) par:

\(\lbrace \begin{array}{cc}{A}^{(1)}=A& \\ {A}^{(k+1)}={Q}^{{(k)}^{T}}{A}^{(k)}{P}^{(k)}& \text{pour}k=1,2,\text{...},m-2\end{array}\)

où chaque matrice courante \({A}^{(k)}\) présente la structure diagonale par bloc suivante:

• Le bloc diagonal supérieur (indices de ligne et de colonne variant de \(1\) à \(k-1\) ) est une matrice bi-diagonale supérieure d’ordre \(k-1\) ,

• Le bloc diagonal inférieur (indices de ligne et de colonne variant de \(k\) à \(m\) ) est une matrice triangulaire supérieure d’ordre \(m-k\) .

Les coefficients \({a}_{i,j}^{(k)}\) des matrices \({A}^{(k)}\) de l’itération vérifient donc:

\({a}_{i,j}^{(k)}=0\text{si}\lbrace \begin{array}{cc}1\le i\le k-1& \text{et}i+2\le j\le m\\ 1\le i\le k-1& \text{et}i<j\\ k\le i\le m& \text{et}1\le j<i\end{array}\) éq 5.3-1

de sorte qu’à l’issue de la récurrence, nous aurons:

\(B={A}^{(m+1)}\text{},\text{}Q=\underset{k=1}{\overset{k=m}{\Pi}}{Q}^{(k)}\text{et}P=\underset{k=1}{\overset{k=m}{\Pi}}{P}^{(k)}\)

Comme pour la factorisation de la section précédente, le problème se réduit à la préservation de la structure [éq 5.3-1] lors du passage de \({A}^{(k)}\) à \({A}^{(k+1)}\) . L’orthogonalité des transformations \({Q}^{(k)}\) et \({P}^{(k)}\) est obtenue en les construisant comme produit de rotations de Givens et le problème de la préservation de la structure est résolu en effectuant ces produits dans un ordre qui ne détruit pas les zéros créés par les étapes précédentes.

L’algorithme annule donc successivement les coefficient \((k,j+1)\) pour \(j=m-1,m-2,\text{...},k+2\) par l’application à droite d’une rotation de Givens du plan \((j,j+1)\) . Cette rotation ne modifie que les colonnes \(j\) et \(j+1\) , ce qui crée un coefficient parasite en position \((j+1,j)\) . Ce coefficient parasite est alors éliminé par l’application à gauche de la transposée d’une rotation de Givens dans le plan \((j,j+1)\) .

Le procédé de passage de \({A}^{(k)}\) à \({A}^{(k+1)}\) peut être alors formalisé par:

\(\begin{array}{}{A}^{(k+1,m-1)}={A}^{(k)}\\ \lbrace \begin{array}{ccc}{A}^{(k+1,j=1/2)}={A}^{(k+1,j)}R(J,J+1,{\theta}_{j}^{(k)})& & \\ {A}^{(k+1,j=1)}=R{(J,J+1,{\theta}_{j=1/2}^{(k)})}^{T}{A}^{(k+1,j=1/2)}& & \text{pour}j=m-1,m-2,\text{...},k+1\end{array}\\ {A}^{(k+1)}={A}^{(k+1,k)}\end{array}\)

où les angles \({\theta}_{j}^{(k)}\) et \({\theta}_{j=1/2}^{(k)}\) sont choisis pour annuler respectivement le coefficient en position \((k,j+1)\) de \({A}^{(k+1,j)}\) et le coefficient en position \((j+1,j)\) de \({A}^{(k+1,j=1/2)}\) .

La structure des matrices \({A}^{(k+1,j)}\) est illustrée dans la figure suivante:

A l’issue de cette récurrence, les matrices et \({P}^{(k)}\) et \({Q}^{(k)}\) sont données par:

\({P}^{(k)}=\underset{j=m-1}{\overset{j=k+1}{\Pi}}R(j,j+1,{\theta}_{j}^{(k)})\text{et}{Q}^{(k)}=\underset{j=m-1}{\overset{j=k+1}{\Pi}}R(k,j,{\theta}_{j=1/2}^{(k)})\)

de sorte que les matrices \(P\) et \(Q\) valent:

\(P=\underset{k=1}{\overset{k=m-2}{\Pi}}\underset{j=m-1}{\overset{j=k+1}{\Pi}}R(j,j+1,{\theta}_{j}^{(k)})\text{et}Q=\underset{k=1}{\overset{k=m-2}{\Pi}}\underset{j=m-1}{\overset{j=k+1}{\Pi}}R(j,j+1,{\theta}_{j=1/2}^{(k)})\)

Décomposition SVD d’une bidiagonale supérieure

Nous présentons un algorithme de construction de la décomposition SVD d’une matrice bidiagonale supérieure \(A\) d’ordre \(m\) . L’algorithme construit donc deux matrices orthogonales \(Q\) et \(P\) et une matrice diagonale \(D\) telles que:

\(A={\text{QDP}}^{T}\)

L’algorithme est tiré de [bib2].

Principe de l’algorithme

L’algorithme de calcul de la décomposition SVD diagonalise itérativement la matrice \(A\) au moyen de la récurrence:

\(\lbrace \begin{array}{cc}{A}^{(1)}=A& \\ {A}^{(k+1)}={Q}^{{(k)}^{T}}{A}^{(k)}{P}^{(k)}& \text{pour}k=1,2,\text{...}\end{array}\) éq 5.4.1-1

où les matrices \({Q}^{(k)}\) et \({P}^{(k)}\) sont orthogonales et les matrices \({A}^{(k)}\) sont bi-diagonales supérieures.

A la convergence, nous aurons:

\(D={A}^{(\infty )},P=\underset{k=1}{\overset{k=\infty }{\Pi}}{P}^{(k)}\text{et}Q=\underset{k=1}{\overset{k=\infty }{\Pi}}{Q}^{(k)}\)

L’idée de l’itération consiste à:

• Choisir \({P}^{(k)}\) pour faire converger l’algorithme \(\text{QR}\) appliqué à la diagonalisation de la matrice (dite normale) \({A}^{T}A\) sans la former explicitement. En effet, la matrice \(P\) de la décomposition SVD de \(A\) n’est rien d’autre que la matrice des vecteurs propres de

\({A}^{T}A\) ,

• Choisir \({Q}^{(k)}\) pour préserver la structure bidiagonale supérieure des itérés successifs.

Comme dans le cas des factorisations présentées aux sections 5.2 et 5.3, les matrices \({Q}^{(k)}\) et \({P}^{(k)}\) sont construites comme produit de rotations de Givens. Le passage de \({A}^{(k)}\) à \({A}^{(k+1)}\) est alors réalisé par:

\({A}^{(k+1)}={\left[{Q}^{(k,2)}{Q}^{(k,3)}\text{...}{Q}^{(k,m)}\right]}^{T}{A}^{(k)}=\left[{P}^{(k,2)}{P}^{(k,3)}\text{...}{P}^{(k,m)}\right]\) éq 5.4.1-2

où les \({Q}^{(k,i)}\) et \({P}^{(k,i)}\) sont deux rotations du plan \((i-1,i)\) d’angle respectif \({\theta}_{j}\) et \({\phi }_{i}\) :

\({Q}^{(k,i)}=R(i-1,i,{\theta}_{i}^{(k)})\text{et}{P}^{(k,i)}=R(i-1,i,{\phi }_{i}^{(k)})\)

Les rotations sont alternativement appliquées à droite puis à gauche de façon à ce que \({A}^{(k+1)}\) conserve la structure bidiagonale supérieure de \({A}^{(k)}\) . Pour ce faire:

L’angle \({\phi }_{2}^{(k)}\) est, pour le moment, choisi arbitrairement; l’application de la rotation \({P}^{(k,2)}\) crée alors un coefficient en position \((2,1)\) ,

L’angle \({\theta}_{2}^{(k)}\) est choisi pour que l’application de la rotation \({Q}^{(k,2)}\) annule le coefficient en position \((2,1)\) , ce qui crée un coefficient non nul en position \((1,3)\) ,

L’angle \({\phi }_{3}^{(k)}\) est choisi pour que l’application de la rotation \({P}^{(k,3)}\) annule le coefficient en position \((1,3)\) , ce qui crée un coefficient non nul en position \((3,2)\) ,

\(\begin{array}{}.\\ .\\ .\end{array}\)

L’angle \({\theta}_{m-1}^{(k)}\) est choisi pour que l’application de la rotation \({Q}^{(k,m-1)}\) annule le coefficient en position \((m-1,m-2)\) , ce qui crée un coefficient non nul en position \((m-2,m)\) ,

L’angle \({\phi }_{m}^{(k)}\) est choisi pour que l’application de la rotation \({P}^{(k,m)}\) annule le coefficient en position \((m-2,m)\) , ce qui crée un coefficient non nul en position \((m,m-1)\) ,

L’angle \({\theta}_{m}^{(k)}\) est choisi pour que l’application de la rotation \({Q}^{(k,m)}\) annule le coefficient en position \((m,m-1)\) , et la matrice \({A}^{(k+1)}\) soit bidiagonale supérieure.

Pour toute valeur de l’angle, ce procédé assure le maintient de la structure bidiagonale supérieure aux itérés [éq 5.4.1-2]. Nous allons voir maintenant comment il est possible de choisir cet angle pour faire converger l’itération [éq 5.4.1-1].

Diagonalisation implicite de la matrice normale

L’algorithme de la sous-section précédente laisse indéterminé l’angle \({\phi }_{2}^{(k)}\) de la première rotation de \({P}^{(k)}\) . Nous allons lever cette indétermination de façon à faire de la matrice \({P}^{(k)}\) la matrice orthogonale d’un pas \(\text{QR}\) , avec décalage spectral, appliqué à la diagonalisation de la matrice normale \(M={A}^{T}A\) .

A l’itération SVD [éq 5.4.1-1] de la matrice \(A\) , nous associons une itération sur la matrice normale \(M={A}^{T}A\) :

\({M}^{(k+1)}={A}^{{(k+1)}^{T}}{A}^{(k+1)}={P}^{{(k)}^{T}}{M}^{(k)}{P}^{(k)}\)

Itération \(\text{QR}\) pour la diagonalisation de la matrice normale

La transformation \(\text{QR}\) , avec décalage spectral \({\sigma}_{k}\) , appliquée à \({M}^{(k)}\) s’écrit:

Factoriser \({M}^{(k)}-{\sigma}_{k}I\) sous la forme \({M}^{(k)}-{\sigma}_{k}I={P}_{\sigma}{R}_{\sigma}\)

Construire \({M}_{\sigma}^{(k+1)}\) par \({M}_{\sigma}^{(k+1)}={P}_{\sigma}{R}_{\sigma}+{\sigma}_{k}I\)

où \({P}_{\sigma}\) et \({R}_{s}\) sont deux matrices respectivement orthogonale et triangulaire supérieure. Les matrices \({M}^{(k)}\) et \({M}_{\sigma}^{(k+1)}\) sont donc tridiagonales et semblables:

\({M}_{\sigma}^{(k+1)}={P}_{\sigma}^{T}{M}^{(k)}{P}_{\sigma}\)

Du point de vue pratique, la matrice \({P}_{\sigma}\) se présente comme un produit de rotations de Givens:

\({P}_{\sigma}^{T}=R(n-1,n,{\psi}_{n})R(n-2,n-1,{\psi}_{n-1})\text{...}R(1,2,{\psi}_{2})\)

Les angles \({\psi}_{k}\) sont choisis pour que l’application à gauche de \(R(k-1,k,{\psi}_{k})\) à la matrice \(\left[\underset{1=k-1,k-2,\text{...},2}{\Pi}R(1-1,1,{\psi}_{1})\right]({M}^{(k)}-{\sigma}_{k}I)\) annule le coefficient de position \((k,k-1)\) dans la matrice résultat.

Francis a montré que le passage de \({M}^{(k)}\) à \({M}_{\sigma}^{(k+1)}\) ne nécessite pas la formation explicite de la matrice \({M}^{(k)}-{\sigma}_{k}I\) : le décalage peut être effectué implicitement. Le théorème s’énonce comme suit:

Théorème (Francis) : Soit \(X\) une matrice orthogonale dont la première colonne coïncide avec celle de \({P}_{\sigma}\) . Sous les hypothèses:

1. \({M}^{(k+1)}={X}^{T}{M}^{(k)}X\)

2. \({M}^{(k+1)}\) est tridiagonale,

3. Les éléments sous-diagonaux de \({M}^{(k)}\) sont tous non nuls (irréductibilité de \({M}^{(k)}\) ),

on a

\({M}^{(k+1)}={\text{DM}}_{\sigma}^{(k+1)}D\)

où \(D\) est une matrice diagonale de coefficients diagonaux tous égaux à \(\pm 1\) .

Application à l’algorithme SVD

Dès lors, le choix de l’angle \({\phi }_{2}^{(k)}\) de la première rotation de l’itération SVD [éq 5.4.1-1] consiste en \({\phi }_{2}^{(k)}\text{=-}{\psi}_{2}\) . Ainsi \(R(1,2,{\phi }_{2}^{(k)})=R{(1,2,{\psi}_{2})}^{T}\) de sorte que la première colonne de \({P}^{(k)}\) coïncide avec la première colonne de \({P}_{\sigma}\) . Donc, si tous les éléments sous-diagonaux de \({M}^{(k)}\) sont non nuls, alors, les matrices \({P}^{(k)}\) et \({P}_{\sigma}\) s’identifient (à un facteur multiplicatif \(\pm 1\) près des colonnes) et l’itération SVD [éq 5.4.1-1] est équivalente à l’application de la transformation \(\text{QR}\) , avec décalage, à la matrice \({M}^{(k)}\) .

Choix du décalage spectral

Le décalage est habituellement choisi comme la valeur propre du mineur inférieur d’ordre deux de \({A}^{(k)}\) la plus proche de \({a}_{m,m}^{(k)}\) . Ce choix assure une convergence globale qui, généralement, est cubique.

Applicabilité du théorème de Francis et phénomène de décomposition

L’utilisation du théorème de Francis suppose la non nullité de tous les coefficients sous-diagonaux des matrices \({M}^{(k)}\) , qui n’est en rien garantie. De plus, dans le cadre de la méthode \(\text{QR}\) de diagonalisation d’une matrice tridiagonale symétrique, l’apparition de coefficients sous-diagonaux nuls est:

• Souhaitable: les coefficients nuls découplent les blocs diagonaux qu’ils encadrent, ce qui ramène la diagonalisation de la matrice complète à la diagonalisation de ses blocs diagonaux (ce phénomène est souvent appelé “décomposition”),

• Inévitable: la convergence de l’algorithme vers une valeur propre s’interprète algébriquement comme l’apparition d’un bloc diagonal précédent d’ordre 1.

Dans la sous-section suivante, nous allons voir quel traitement il convient d’adopter en présence d’une décomposition.

Analyse de décomposition

L’analyse de décomposition porte sur chaque itéré pris indépendamment des autres, aussi, nous n’utiliserons pas l’indice supérieur \((k)\) .

Soit \({d}_{1},{d}_{2},\text{...},{d}_{m}\) et \({e}_{2},{e}_{3},\text{...},{e}_{m}\) les éléments respectivement diagonaux et sur-diagonaux de \(A\) . Les éléments sous-diagonaux de la matrice normale \(M={A}^{T}A\) sont alors donnés par:

\({m}_{i+1,i}={d}_{i}{e}_{l+1}\text{pour}i=1,2,\text{...},m-1\)

Supposons, pour simplifier, que seul le coefficient \({m}_{l-1,l}\) est nul. La matrice \(M\) présente alors une structure de deux blocs diagonaux dont la réunion des spectres respectifs donne le spectre de \(M\) . Cette décomposition a lieu soit pour \({e}_{1}=0\) soit pour \({e}_{1}\ne 0\text{et}{d}_{l-1}=0\) .

Le cas \({e}_{1}=0\) ne pose aucune difficulté. La matrice \(A\) possède alors une structure diagonale de deux blocs qui fournissent chacun une partie complémentaire de la décomposition SVD de \(A\) . Chaque bloc étant bi-diagonal supérieur, sans coefficients sur-diagonaux nuls, sa décomposition SVD est calculable par l’itération [éq 5.4.1-1].

Le cas \({e}_{1}\ne 0\text{et}{d}_{l-1}=0\) est plus délicat. En effet, l’itération [éq 5.4.1-1] ne peut être appliquée ni à la matrice \(A\) , pour ne pas violer les hypothèses du théorème de Francis, ni à aucune sous matrices de \(A\) , pour assurer la structure bi-diagonale aux itérés. Ce problème est contourné par la post-multiplication de \(A\) par une série de rotations de Givens dans les plans successifs \((l-1,l),(l-1,l+1),\text{...},(l-1,m)\) :

• La rotation du plan \((l-1,l)\) annule le coefficient \((l-1,l)\) et crée un coefficient en position \((l-1,l+1)\) ,

• La rotation du plan \((l-1,l+1)\) annule le coefficient \((l-1,l+1)\) et crée un coefficient en position \((l-1,l+2)\) ,

• La rotation du plan \((l-1,l+1)\) annule le coefficient \((l-1,l+1)\) et crée un coefficient en position \((l-1,l+2)\) ,

\(\begin{array}{}.\\ .\\ .\end{array}\)

• La rotation du plan \((l-1,m)\) annule le coefficient \((l-1,m)\) et ne crée pas de coefficient.

De sorte que la matrice produite par ce procédé présente la même structure que celle correspondant au cas \({e}_{1}=0\) .

Organisation de l’algorithme

L’algorithme isole successivement chaque valeur singulière, aussi, il existe un indice \({k}_{p}\) tel que l’itéré \({A}^{({k}_{p})}\) se décompose en deux blocs diagonaux:

où \({B}^{({k}_{p})}\) est une matrice bidiagonale supérieure d’ordre \(p\) et \({D}^{({k}_{p})}\) est une matrice diagonale d’ordre \(m-p+1\) rassemblant sur sa diagonale les valeurs singulières trouvées.

A partir de cet itéré, l’algorithme applique l’itération de la sous-section 5.4.1 à la sous-matrice \({B}^{({k}_{p})}\) jusqu’à l’annulation du coefficient en position \((p-1,p)\) , signal de la convergence de la \({p}^{\text{ième}}\) valeur singulière. Chaque pas de l’itération interne, ainsi définie, s’organise comme suit:

• Analyse de décomposition,

• Dans le cas où la décomposition trouvée correspond à un élément diagonal nul, une série de rotations supplémentaires est appliquée pour retrouver la structure de décomposition générée par un élément sur-diagonal nul. Ces rotations sont construites suivant la méthode présentée à la sous-section 5.4.3,

• Si le coefficient en position \((p-1,p)\) ne produit pas de décomposition alors la sous‑matrice \({B}^{({k}_{p})}\) est l’objet d’un pas de l’itération [éq 5.4.1-1] où, conformément à l’analyse de la sous-section 5.4.2, l’angle de la première rotation est choisi pour que ce pas soit équivalent à l’application d’une transformation \(\text{QR}\) , avec décalage spectral implicite, sur la matrice normale associée.

La convergence complète de l’itération est alors obtenue à l’indice \({k}_{m}\) pour lequel la sous-matrice \({D}^{({k}_{p})}\) est d’ordre \(m\) .

Bien entendu, dans la pratique, un coefficient est considéré comme nul dès qu’il est inférieur, en valeur absolue, à une certaine tolérance. La tolérance généralement utilisée pour les problèmes de valeurs singulières est choisie comme le produit de la précision machine par \({\parallel A\parallel }_{1}\) . Remarquons que dans le cas d’une décomposition produite par un élément diagonal nul à la tolérance choisie, l’application de la série de rotations supplémentaires décrites à la sous-section 5.4.3 crée une sous colonne d’éléments non nuls sous ce coefficient. Ces éléments ne sont pas gênants car ils sont tous nuls à la tolérance choisie.

Bibliographie

1. 16. 7. CIARLET : “Introduction à l’analyse numérique matricielle et à l’optimisation” _ MASSON (1985).

2. 7. 8. GOLUB, C. REINSCH “Singular value decomposition and least squares solutions” in “Handbook for Automatic Computation - Linear Algebra, Vol. 2” J.H. WILKHINSON, C.REINSH Editors _ SPINGER VERLAG (1971).

3. 16. LASCAUX, R. THEODOR : “ Analyse numérique matricielle appliquée à l’art de l’ingénieur”, tomes 1 et 2 _ MASSON (1986).

Description des versions du document

Version Aster	Auteur(s) Organisme(s)	Description des modifications
3	B.QUINNEZ, R.MICHEL EDF-R&D/MMN	Texte initial