Nous introduisons une décomposition matricielle qui nous permettra plus tard d'implémenter les calculs des décompositions en valeurs propres et en valeurs singulières.
Nous cherchons des vecteurs orthonormaux qui forment une base pour les espaces couverts par les colonnes successives d'une matrice $\mathbf{A} \in \mathbb{R}^{M \times N}$. Notons :
\begin{align*}
<\mathbf{a_1}> &: \text{l'espace (la ligne) couvert par la première colonne de $\mathbf{A}$} \\
<\mathbf{a_1}, \mathbf{a_2}> &: \text{l'espace (le plan) couvert par les deux premières colonne de $\mathbf{A}$} \\
Sans perte de généralité, considérons la situation où $M>N$. Nous appelons alors cette expression, $\mathbf{A}=\mathbf{\hat{Q}}\mathbf{\hat{R}}$ (avec $\mathbf{\hat{Q}} \in \mathbb{R}^{M \times N}$ et $\mathbf{\hat{R}} \in \mathbb{R}^{N \times N}$), la forme réduite de la factorisation QR.
Nous pouvons aussi considérer la forme dite complète de la factorisation QR, avec $\mathbf{Q} \in \mathbb{R}^{M \times M}$ et $\mathbf{R} \in \mathbb{R}^{M \times N}$. Dans ce cas, les $M-N$ dernières lignes de $\mathbf{R}$ sont nulles et il suffit de compléter $\mathbf{\hat{Q}}$ en lui ajoutant $M-N$ colonnes orthogonales.
# Orthogonalisation de Gram-Schmidt
Nous présentons l'algorithme de Gram-Schmidt pour la décomposition QR.
Le principe est de construire $\mathbf{q_j}$ orthogonal à $\mathbf{q_1}, \mathbf{q_2}, \dots, \mathbf{q_{j-1}}$ en retirant à la colonne $\mathbf{a_j}$ ses composantes selon les directions de $\mathbf{q_1}, \mathbf{q_2}, \dots, \mathbf{q_{j-1}}$ :
## Dérivation de l'algorithme de Gram-Schmidt modifié
Nous présentons une autre manière de voir ce résultat qui nous mènera à un algorithme équivalent mais numériquement plus stable.
Nous pouvons considérer que $\mathbf{q_j}$ est obtenu en projetant $\mathbf{a_j}$ sur l'espace orthogonal à $<\mathbf{a_1}, \mathbf{a_2}, \dots, \mathbf{a_{j-1}}>$, cet espace étant par construction égal à $<\mathbf{q_1}, \mathbf{q_2}, \dots, \mathbf{q_{j-1}}>$ :
De façon équivalente, $\mathbf{P_j}$ peut s'exprimer comme le produit de projecteurs sur les espaces orthogonaux à $\mathbf{q_1}$, puis $\mathbf{q_2}$, etc. jusqu'à $\mathbf{q_{j-1}}$ :
Ainsi, nous observons que dès que $\mathbf{q_1}$ est connu ($\mathbf{q_1} = \mathbf{v_1} / \|\mathbf{v_1}\|$, $\mathbf{v_1} = \mathbf{a_1}$), nous pouvons appliquer $\mathbf{P_{\bot \mathbf{q_{1}}}}$ à $\mathbf{a_2}, \mathbf{a_3}, \dots, \mathbf{a_N}$. Puis, quand $\mathbf{q_2}$ est connu, nous pouvons appliquer $\mathbf{P_{\bot \mathbf{q_{2}}}}$ à $\mathbf{P_{\bot \mathbf{q_{1}}}}\mathbf{a_3}, \mathbf{P_{\bot \mathbf{q_{1}}}}\mathbf{a_4}, \dots, \mathbf{P_{\bot \mathbf{q_{1}}}}\mathbf{a_N}$. Etc. Nous en déduisons l'algorithme de Gram-Schmidt modifié.
Nous comparons notre implémentation (voir la fonction `mgs` dans le code source associé à ce module) avec le résultat d'un petit exemple issu de [Wikipedia](https://en.wikipedia.org/wiki/Gram–Schmidt_process#Example).
```{r}
mat <- matrix(c(3,2,1,2),nrow=2,ncol=2,byrow = TRUE)
## Comparaison de la stabilité numérique des deux versions de l'algorithme de Gram-Schmidt
Observons comment de petites erreurs de calculs se propagent dans le cas de l'algorithme de Gram-Schmidt (GS). Mettons qu'une petite erreur ait été commise sur le calcul de $\mathbf{q_2}$ qui s'en trouve de ce fait n'être pas parfaitement perpendiculaire à $\mathbf{q_1}$ :
\[
\mathbf{q_1}^T\mathbf{q_2} = \epsilon > 0
\]
Quelle est l'effet de cette erreur sur la valeur de $\mathbf{v_3}$ ?
Il n'y a pas d'erreur pour l'orthogonalité entre $\mathbf{q_3}$ et $\mathbf{q_2}$. Etudions maintenant l'erreur commise sur l'orthogonalité entre $\mathbf{q_3}$ et $\mathbf{q_1}$.
Ainsi, comme pour la version classique de Gram-Schmidt, une erreur de l'ordre de grandeur de $\epsilon$ est commise sur l'orthogonalité entre $\mathbf{q_1}$ et $\mathbf{q_3}$. Par contre, il n'y a pas de propagation de l'erreur pour l'orthogonalité entre $\mathbf{q_2}$ et $\mathbf{q_3}$. C'est pourquoi l'algorithme de Gram-Schmidt modifié est numériquement plus stable que l'algorithme de Gram-Schmidt classique.
# Autres algorithmes
Notons qu'une autre approche très souvent employée pour le calcul de la décomposition QR est la méthode dite de Householder. Elle assure le maintien d'une quasi parfaite orthogonalité des colonnes de $\mathbf{Q}$. Il existe également l'approche dite des rotations de Givens, plus facilement parallélisable et permettant la mise à jour de la décomposition suite à l'ajout d'une ligne à la matrice $\mathbf{X}$.
Par contre, Gram-Schmidt permet de calculer la forme réduite de la décomposition QR, alors que les algorithmes de Givens et Householder calculent la forme complète. Ainsi, Gram-Schmidt n'est pas tout à fait aussi stable mais peut parfois être plus efficace.