Décomposition QR

Soient x un vecteur colonne arbitraire de dimension m et α = ± ||x||, où || || désigne la norme euclidienne. Pour des raisons de stabilité du calcul, α doit de plus être du signe opposé au premier élément de x.

Soit e₁ le vecteur (1, 0, ..., 0)^T, et définissons, si x n'est pas colinéaire à e₁ :

${\displaystyle \mathbf {u} =\mathbf {x} -\alpha \mathbf {e} _{1},}$

${\displaystyle \mathbf {v} ={\mathbf {u} \over ||\mathbf {u} ||},}$

${\displaystyle Q_{1}=I-2\mathbf {v} \mathbf {v} ^{T}.}$

Q₁ est la matrice de Householder ou matrice orthogonale élémentaire et

${\displaystyle Q_{1}x=(\alpha \ ,0,\cdots ,0)^{T}.}$

(Si x est colinéaire à e₁, on a le même résultat en prenant pour Q la matrice identité.)

Nous pouvons utiliser ces propriétés pour transformer une matrice A de dimension m×n en une matrice triangulaire supérieure. Tout d'abord, on multiplie A par la matrice de Householder Q₁ en ayant pris le soin de choisir pour x la première colonne de A. Le résultat est une matrice QA avec des zéros dans la première colonne excepté du premier élément qui vaudra α.

${\displaystyle Q_{1}A={\begin{pmatrix}\alpha _{1}&\star &\dots &\star \\0&&&\\\vdots &&A'&\\0&&&\end{pmatrix}}}$

Ceci doit être réitéré pour A' qui va être multipliée par Q’₂ (Q’₂ est plus petite que Q₁). Si toutefois, vous souhaitez utiliser Q₁A plutôt que A', vous devez remplir la matrice de Householder avec des 1 dans le coin supérieur gauche :

${\displaystyle Q_{k}={\begin{pmatrix}I_{k-1}&0\\0&Q_{k}'\end{pmatrix}}}$

Après t itérations, t = min(m – 1, n),

${\displaystyle R=Q_{t}\cdots Q_{2}Q_{1}A}$

est une matrice triangulaire supérieure. Si Q = Q^T
₁ Q^T
₂ ... Q^T
_t alors A = QR est la décomposition QR de A. De plus, par construction les matrices Q_k sont non seulement orthogonales mais aussi symétriques, donc Q = Q₁ Q₂ ... Q_t.

Calculons la décomposition QR de

${\displaystyle A={\begin{pmatrix}12&-51&4\\6&167&-68\\-4&24&-41\end{pmatrix}}}$

On choisit donc le vecteur a₁ = (12, 6, -4)^T. On a donc ||a₁|| = √12² + 6² + (–4)² = 14. Ce qui nous conduit à écrire ||a₁|| e₁ = (14, 0, 0)^T.

Le calcul nous amène à u = 2(–1, 3, –2)^T et ${\displaystyle v={1 \over 14^{1 \over 2}}(-1,3,-2)^{T}}$. La première matrice de Householder vaut

${\displaystyle {\begin{aligned}Q_{1}&=I-{2 \over 14}{\begin{pmatrix}-1\\3\\-2\end{pmatrix}}{\begin{pmatrix}-1&3&-2\end{pmatrix}}\\&=I-{1 \over 7}{\begin{pmatrix}1&-3&2\\-3&9&-6\\2&-6&4\end{pmatrix}}={\begin{pmatrix}6/7&3/7&-2/7\\3/7&-2/7&6/7\\-2/7&6/7&3/7\\\end{pmatrix}}.\end{aligned}}}$

Observons que

${\displaystyle Q_{1}A={\begin{pmatrix}14&21&-14\\0&-49&-14\\0&168&-77\end{pmatrix}}.}$

Nous avons maintenant sous la diagonale uniquement des zéros dans la 1^re colonne.

Pour réitérer le processus, on prend la sous matrice principale

${\displaystyle A'=M_{11}={\begin{pmatrix}-49&-14\\168&-77\end{pmatrix}}}$

Par la même méthode, on obtient

${\displaystyle \alpha _{2}={\sqrt {(-49)^{2}+168^{2}}}=175,\quad u_{2}=(-224,168)^{T},\quad v_{2}=(-4/5,3/5)^{T},\quad Q'_{2}={\begin{pmatrix}-7/25&24/25\\24/25&7/25\end{pmatrix}}.}$

La 2^e matrice de Householder est donc

${\displaystyle Q_{2}={\begin{pmatrix}1&0&0\\0&-7/25&24/25\\0&24/25&7/25\end{pmatrix}}.}$

Finalement, on obtient

${\displaystyle Q=Q_{1}Q_{2}={\begin{pmatrix}6/7&-69/175&58/175\\3/7&158/175&-6/175\\-2/7&6/35&33/35\end{pmatrix}}}$

${\displaystyle R=Q_{2}Q_{1}A=Q^{T}A={\begin{pmatrix}14&21&-14\\0&-175&70\\0&0&35\end{pmatrix}}.}$

La matrice Q est orthogonale et R est triangulaire supérieure, par conséquent, on obtient la décomposition A = QR.

Le coût de cette méthode pour une matrice n×n est en : 4/3n³ Ce coût est relativement élevé (la méthode de Cholesky, pour les matrices symétriques définies positives est en 1/3n³). Cependant, la méthode de Householder présente l'avantage considérable d'être beaucoup plus stable numériquement, en limitant les divisions par des nombres petits. La méthode de Givens, malgré un coût encore supérieur à celui-ci, offrira encore davantage de stabilité.

On considère le procédé de Gram-Schmidt appliqué aux colonnes de la matrice ${\displaystyle A=[\mathbf {a} _{1},\cdots ,\mathbf {a} _{n}]}$, muni du produit scalaire ${\displaystyle \langle \mathbf {v} ,\mathbf {w} \rangle =\mathbf {v} ^{\top }\mathbf {w} }$ (ou ${\displaystyle \langle \mathbf {v} ,\mathbf {w} \rangle =\mathbf {v} ^{*}\mathbf {w} }$ pour le cas complexe). L'algorithme présenté ci-dessous convient à une matrice de rang ${\displaystyle n}$. Pour des matrices de rang inférieur, il est à adapter à chaque fois que le vecteur ${\displaystyle \mathbf {u} _{i}}$obtenu est nul.

On définit la projection :

${\displaystyle \Pi _{\mathbf {e} }\mathbf {a} ={\frac {\langle \mathbf {e} ,\mathbf {a} \rangle }{\langle \mathbf {e} ,\mathbf {e} \rangle }}\mathbf {e} }$

puis les vecteurs :

${\displaystyle {\begin{aligned}\mathbf {u} _{1}&=\mathbf {a} _{1},&\mathbf {e} _{1}&={\mathbf {u} _{1} \over \|\mathbf {u} _{1}\|}\\\mathbf {u} _{2}&=\mathbf {a} _{2}-\Pi _{\mathbf {e} _{1}}\,\mathbf {a} _{2},&\mathbf {e} _{2}&={\mathbf {u} _{2} \over \|\mathbf {u} _{2}\|}\\\mathbf {u} _{3}&=\mathbf {a} _{3}-\Pi _{\mathbf {e} _{1}}\,\mathbf {a} _{3}-\Pi _{\mathbf {e} _{2}}\,\mathbf {a} _{3},&\mathbf {e} _{3}&={\mathbf {u} _{3} \over \|\mathbf {u} _{3}\|}\\&\vdots &&\vdots \\\mathbf {u} _{k}&=\mathbf {a} _{k}-\sum _{j=1}^{k-1}\Pi _{\mathbf {e} _{j}}\,\mathbf {a} _{k},&\mathbf {e} _{k}&={\mathbf {u} _{k} \over \|\mathbf {u} _{k}\|}.\end{aligned}}}$

On réarrange ensuite les équations de sorte que les a_i soient à gauche, en utilisant le fait que les e_i sont des vecteurs unitaires :

${\displaystyle {\begin{aligned}\mathbf {a} _{1}&=\langle \mathbf {e} _{1},\mathbf {a} _{1}\rangle \mathbf {e} _{1}\\\mathbf {a} _{2}&=\langle \mathbf {e} _{1},\mathbf {a} _{2}\rangle \mathbf {e} _{1}+\langle \mathbf {e} _{2},\mathbf {a} _{2}\rangle \mathbf {e} _{2}\\\mathbf {a} _{3}&=\langle \mathbf {e} _{1},\mathbf {a} _{3}\rangle \mathbf {e} _{1}+\langle \mathbf {e} _{2},\mathbf {a} _{3}\rangle \mathbf {e} _{2}+\langle \mathbf {e} _{3},\mathbf {a} _{3}\rangle \mathbf {e} _{3}\\&\vdots \\\mathbf {a} _{k}&=\sum _{j=1}^{k}\langle \mathbf {e} _{j},\mathbf {a} _{k}\rangle \mathbf {e} _{j}\end{aligned}}}$

où ${\displaystyle \langle \mathbf {e} _{i},\mathbf {a} _{i}\rangle =\|\mathbf {u} _{i}\|}$. Ceci s'écrit matriciellement :

${\displaystyle A=QR}$

avec

${\displaystyle Q=\left[\mathbf {e} _{1},\cdots ,\mathbf {e} _{n}\right]\qquad {\text{et}}\qquad R={\begin{pmatrix}\langle \mathbf {e} _{1},\mathbf {a} _{1}\rangle &\langle \mathbf {e} _{1},\mathbf {a} _{2}\rangle &\langle \mathbf {e} _{1},\mathbf {a} _{3}\rangle &\ldots \\0&\langle \mathbf {e} _{2},\mathbf {a} _{2}\rangle &\langle \mathbf {e} _{2},\mathbf {a} _{3}\rangle &\ldots \\0&0&\langle \mathbf {e} _{3},\mathbf {a} _{3}\rangle &\ldots \\\vdots &\vdots &\vdots &\ddots \end{pmatrix}}.}$

On reprend la matrice de l'exemple

${\displaystyle A={\begin{pmatrix}12&-51&4\\6&167&-68\\-4&24&-41\end{pmatrix}}.}$

Rappelons qu'une matrice orthogonale Q vérifie

${\displaystyle {\begin{matrix}Q^{T}\,Q=I.\end{matrix}}}$

On peut alors calculer Q par les moyens de Gram-Schmidt comme suit :

${\displaystyle U={\begin{pmatrix}\mathbf {u} _{1}&\mathbf {u} _{2}&\mathbf {u} _{3}\end{pmatrix}}={\begin{pmatrix}12&-69&-58/5\\6&158&6/5\\-4&30&-33\end{pmatrix}};}$

${\displaystyle Q={\begin{pmatrix}{\frac {\mathbf {u} _{1}}{\|\mathbf {u} _{1}\|}}&{\frac {\mathbf {u} _{2}}{\|\mathbf {u} _{2}\|}}&{\frac {\mathbf {u} _{3}}{\|\mathbf {u} _{3}\|}}\end{pmatrix}}={\begin{pmatrix}6/7&-69/175&-58/175\\3/7&158/175&6/175\\-2/7&6/35&-33/35\end{pmatrix}}.}$

Dans ce cas, on a :

${\displaystyle {\begin{matrix}Q^{T}A=Q^{T}Q\,R=R;\end{matrix}}}$

${\textstyle {\begin{matrix}R=Q^{T}A=\end{matrix}}{\begin{pmatrix}14&21&-14\\0&175&-70\\0&0&35\end{pmatrix}}.}$

La décomposition RQ transforme une matrice A en produit d'une matrice triangulaire supérieure R et une matrice orthogonale Q. La seule différence avec la décomposition QR est l'ordre de ces matrices.

La décomposition QR est l'application du procédé de Gram-Schmidt sur les colonnes de A, en partant de la première colonne ; la décomposition RQ est l'application du procédé de Gram-Schmidt sur les lignes de A, en partant de la dernière ligne.

Dans cette méthode, la matrice Q utilise des rotations de Givens (en). Chaque rotation annule un élément de la partie triangulaire inférieure stricte de la matrice, construisant la matrice R, tandis que la concaténation des rotations engendre la matrice Q.

Dans la pratique, les rotations de Givens ne sont pas effectivement assurées par la construction d'une matrice pleine et une multiplication matricielle. Une procédure de rotation de Givens est utilisé à la place qui est l'équivalent de la multiplication par une matrice de Givens creuse, sans efforts supplémentaires de la manipulation des éléments non nuls. La procédure de rotation de Givens est utile dans des situations où seul un nombre relativement restreint hors éléments diagonaux doivent être remis à zéro, et est plus facilement parallélisée que les transformations de Householder.

Reprenons le même exemple

${\displaystyle A={\begin{pmatrix}12&-51&4\\6&167&-68\\-4&24&-41\end{pmatrix}}.}$

On doit d'abord construire une matrice de rotation qui annulera l'élément le plus bas de la colonne de gauche, a₃₁ = –4, qu'on construit par une méthode de rotation de Givens. On appelle cette matrice G₁. On va d'abord faire une rotation du vecteur (6,-4), pour le ramener sur l'axe X. Ce vecteur forme un angle θ = arctan –(–4)/6. La matrice G₁ est donc donnée par :

${\displaystyle G_{1}={\begin{pmatrix}1&0&0\\0&\cos(\theta )&-\sin(\theta )\\0&\sin(\theta )&\cos(\theta )\end{pmatrix}}}$

${\displaystyle \approx {\begin{pmatrix}1&0&0\\0&0{,}83205&-0{,}55470\\0&0{,}55470&0{,}83205\end{pmatrix}}.}$

Le produit G₁A annule le coefficient a₃₁ :

${\displaystyle G_{1}A\approx {\begin{pmatrix}12&-51&4\\7{,}21110&125{,}6396&-33{,}83671\\0&112{,}6041&-71{,}83368\end{pmatrix}}.}$

Par suite, on construit des matrices de Givens G₂ et G₃, qui vont respectivement annuler a₂₁ et a₃₂, engendrant la matrice R. La matrice orthogonale Q^T est formée de la concaténation de toutes les matrices de Givens créées Q^T = G₃G₂G₁.