Décomposition QR

Tous le code Python proposé est accessible dans un notebook Jupyter interactif.

Notations

On adopte les notations suivantes :

$T_{n}^{+ +} (R)$ désigne l’ensemble des matrices triangulaires supérieures de $M_{n} (R)$ à coefficients diagonaux strictement positifs.
$⟨ \cdot, \cdot ⟩$ désigne le produit scalaire usuel sur $M_{n, 1} (R)$ i.e. $\forall (X, Y) \in M_{n, 1} (R)^{2}$ , $⟨ X, Y ⟩ = X^{⊤} Y$ .

Théorème

Soit $A \in G L_{n} (R)$ . Alors il existe $Q \in O_{n} (R)$ et $R \in T_{n}^{+ +} (R)$ telles que $A = Q R$ .

Intérêt

Dès lors que l’on a une décomposition QR d’une matrice $A$ , il devient trivial de résoudre un sytème $A X = Y$ . En effet ce système équivaut à $R X = Q^{- 1} Y = Q^{⊤} Y$ . Le calcul de $Q^{⊤}$ est immédiat et ce dernier système est triangulaire : il se résout donc très simplement.

Rappels sur le procédé d’orthonormalisation de Gram-Schmidt

Si $(f_{1}, \dots, f_{n})$ est une base d’un espace euclidien $E$ de dimension $n$ , on peut la transformer en une base orthonormée $(e_{1}, \dots, e_{n})$ en appliquant le procédé suivant :

on norme le premier vecteur $f_{1}$ en un vecteur $e_{1}$ ;
une fois construits les vecteurs $e_{1}, \dots, e_{k}$ , on norme le projeté orthogonal de $f_{k + 1}$ sur $vect (e_{1}, \dots, e_{k})^{⊥}$ pour obtenir le vecteur $e_{k + 1}$ .

Du point de vue du calcul, cela revient à construire par récurrence les vecteurs $f_{1}, \dots, f_{n}$ via la relation de récurrence

\forall k \in [[0, n - 1]], e_{k + 1} = \frac{{\tilde{e}}_{k + 1}}{‖ {\tilde{e}}_{k + 1} ‖} où {\tilde{e}}_{k + 1} = f_{k + 1} - \sum_{j = 1}^{k} ⟨ f_{k + 1}, e_{j} ⟩ e_{j}

Remarque. On convient que la somme est nulle lorsque $k = 0$ .

En plus du fait que $(e_{1}, \dots, e_{n})$ est une base orthonormée de $E$ , on vérifie aisément que

\forall k \in [[1, n]], (vect (e_{1}, \dots, e_{k}) = vect (f_{1}, \dots, f_{k})) \land (⟨ e_{k}, f_{k} ⟩ > 0)

La décomposition QR via le procédé de Gram-Schmidt

Notons $f_{1}, \dots, f_{n}$ les colonnes de la matrice $A \in G L_{n} (R)$ . Comme la matrice $A$ est inversible, la famille $(f_{1}, \dots, f_{n})$ est une base de $M_{n, 1} (R)$ . On lui applique le procédé d’orthonormalisation de Gram-Schmidt pour obtenir une base orthonormée (pour le produit scalaire canonique) $(e_{1}, \dots, e_{n})$ de $M_{n, 1} (R)$ .

Notons $Q$ la matrice dont les colonnes sont $e_{1}, \dots, e_{n}$ . Comme $(e_{1}, \dots, e_{n})$ est une base orthonormée, $Q$ est une matrice orthogonale. En notant $R$ la matrice de la base $(f_{1}, \dots, f_{n})$ dans la base $(e_{1}, \dots, e_{n})$ , on a alors bien $A = Q R$ .

De plus, on rappelle que

\forall k \in [[1, n]], (vect (e_{1}, \dots, e_{k}) = vect (f_{1}, \dots, f_{k})) \land (⟨ e_{k}, f_{k} ⟩ > 0)

donc la matrice $R$ est triangulaire supérieure et ses coefficients diagonaux sont strictement positifs. En effet, comme $(e_{1}, \dots, e_{n})$ est une base orthonormée, $R_{i, j} = ⟨ f_{j}, e_{i} ⟩$ pour tout $(i, j) \in {[[1, n]]}^{2}$ et notamment $R_{k, k} = ⟨ f_{k}, e_{k} ⟩ > 0$ pour tout $k \in [[1, n]]$ .

Algorithme en Python

import numpy as np

def QR_gram_schmidt(A):
    Q = np.zeros_like(A)
    R = np.zeros_like(A)
    for i in range(A.shape[1]):
        dp = Q[:, :i].T@A[:, i]
        R[:i, i] = dp
        Q[:, i] = A[:, i]-Q[:, :i]@dp
        Q[:, i] = Q[:, i]/np.linalg.norm(Q[:, i])
        R[i, i] = Q[:, i].T@A[:, i]
    return Q, R

Méthode de Householder

Dans ce qui suit, $E$ désigne un espace euclidien de dimension $n$ et $B$ une base orthonormée de $E$ .

Matrices de Householder

Soient $n$ un vecteur non nul de $E$ et $s$ la réflexion par rapport à l’hyperplan $vect (n)^{⊥}$ . Le projeté orthogonal d’un vecteur $x \in E$ sur $vect (n)$ est $\frac{⟨ n, x ⟩}{| n |^{2}} n$ . On en déduit que

s (x) = x - 2 \frac{⟨ n, x ⟩}{‖ n ‖^{2}} n

Si on note $N$ sa matrice de $n$ dans la base $B$ , alors la matrice de $s$ dans la base $B$ est donc

I_{n} - 2 \frac{N N^{⊤}}{N^{⊤} N}

Une telle matrice est appelée une matrice de Householder.

Réflexion échangeant deux vecteurs de même norme

Soient $E$ un espace euclidien de dimension $n$ ainsi que $u$ et $v$ deux vecteurs distincts et de même norme. Il existe alors une (unique) réflexion $s$ échangeant $u$ et $v$ . Il s’agit de la réflexion par rapport à l’hyperplan $vect (v - u)^{⊥}$ .

Notamment si $S$ , $U$ et $V$ sont respectivement les matrices de $s$ , $u$ et $v$ dans cette base $B$ , alors

S = I_{n} - 2 \frac{(U - V) (U - V)^{⊤}}{(U - V)^{⊤} (U - V)}

Application à la décomposition QR

On procède à un raisonnement par récurrence. Il est clair que toute matrice de $G L_{1} (R)$ admet une décomposition $Q R$ . Supposons que ce soit vrai pour toute matrice de $G L_{n - 1} (R)$ , où $n$ est un entier supérieur ou égal à $2$ .

Soit $A \in G L_{n} (R)$ . Notons $U$ sa première colonne (nécessairement non nulle). Comme les vecteurs $U$ et $V = (| U |, 0, \dots, 0)^{⊤}$ ont même norme, il existe une matrice de réflexion $S$ telle que $S U = V$ . La première colonne de $S A$ est alors $V$ . Notons $\tilde{A}$ la matrice $S A$ dont on a supprimé les premières lignes et premières colonnes. Autrement dit

S A = (\begin{matrix} ‖ U ‖ & \begin{matrix} ⋆ & \dots & ⋆ \end{matrix} \\ \begin{matrix} 0 \\ ⋮ \\ 0 \end{matrix} & \tilde{A} \end{matrix})

Comme $S$ et $A$ sont inversibles, $det (S A) = | U | det (\tilde{A}) \neq 0$ donc $\tilde{A} \in G L_{n - 1} (R)$ . Il existe donc une matrice $\tilde{Q} \in O_{n - 1} (R)$ et $\tilde{R} \in T_{n - 1}^{+ +} (R)$ telle que $\tilde{A} = \tilde{Q} \tilde{R}$ . Puisque $S = S^{- 1}$ ( $S$ est une matrice de réflexion), on a alors $A = Q R$ avec

Q = S (\begin{matrix} 1 & \begin{matrix} 0 & \dots & 0 \end{matrix} \\ \begin{matrix} 0 \\ ⋮ \\ 0 \end{matrix} & \tilde{Q} \end{matrix}) \in O_{n} (R) et R = (\begin{matrix} ‖ U ‖ & \begin{matrix} ⋆ & \dots & ⋆ \end{matrix} \\ \begin{matrix} 0 \\ ⋮ \\ 0 \end{matrix} & \tilde{R} \end{matrix}) \in T_{n}^{+ +} (R)

On en déduit également un algorithme :

Entrée Une matrice $A \in G L_{n} (R)$

Initialisation
- $R \leftarrow A$
- $Q \leftarrow I_{n}$
Pour $k$ variant de $1$ à $n$ :
- On extrait la première colonne $U \in M_{n - k + 1, 1} (R)$ de la matrice extraite $(R_{i, j})_{k \leq i, j \leq n}$
- Si $U_{1} \neq | U |$
  - $V \leftarrow (| U |, 0, \dots, 0)$ .
  - $S \leftarrow I_{n - k + 1} - 2 \frac{(U - V) (U - V)^{⊤}}{(U - V)^{⊤} (U - V)}$
  - $Q_{1} \leftarrow (\begin{matrix} I_{k - 1} & 0 \\ 0 & S \end{matrix})$
  - $Q \leftarrow Q Q_{1}$
  - $R \leftarrow Q_{1} R$

Renvoyer $Q$ et $R$

On définit d’abord une fonction calculant la matrice de Householder associée à un vecteur $U$ .

import numpy as np

def matrice_householder(U):
    return np.eye(U.shape[0]) - 2 * (U@U.T) / (U.T@U)

On implémente alors un algorithme itératif calculant la décomposition QR d’une matrice inversible $A$ .

def QR_householder(A):
    n = A.shape[0]
    R = A.copy()
    Q = np.eye(n)
    for k in range(n):
        U = R[k:, k].copy()
        norme_U = np.linalg.norm(U)
        if U[0] != norme_U:
            U[0] -= norme_U
            S = matrice_householder(U)
            Q[:, k:] = Q[:, k:]@S
            R[k:, :] = S@R[k:, :]
    return Q, R

On peut également proposer une fonction récursive.

def QR_householder_recu(A):
    if A.shape[0] == 1:
        return (np.matrix([[1]]), A) if A[0, 0] > 0 else (np.matrix([[-1]]), -A)
    U = A[:, 0].copy()
    norme_U = np.linalg.norm(U)
    if U[0] != norme_U:
        U[0] -= norme_U
        Q = matrice_householder(U)
        R = Q@A
    else:
        Q = np.eye(A.shape[0])
        R = A
    Q1, R1 = QR_householder_recu(R[1:, 1:])
    Q[:, 1:] = Q[:, 1:]@Q1
    R[1:, 1:] = R1
    return Q, R