Valeurs propres et SVD : calculer, reconstruire et approximer

Décomposer l’action d’une matrice selon quelques directions permet de comprendre comment elle étire les vecteurs, quelles informations la transformation efface et quelle erreur reste si l’on conserve moins de directions. La vue d’ensemble de l’algèbre linéaire présente les vecteurs, le rang, la projection orthogonale et les moindres carrés. Partons d’un calcul de valeurs propres pour une matrice 2×2, puis utilisons une matrice 3×2 pour effectuer une SVD, reconstruire la matrice et l’approximer par une matrice de rang inférieur. Toutes les matrices sont réelles, les vecteurs sont des colonnes et $T$ désigne la transposée.

Une droite invariante sous la transformation

Pour une matrice carrée $B$, un vecteur non nul $v$ qui vérifie $Bv=\lambda v$ est un vecteur propre ; $\lambda$ est sa valeur propre. Les notes du MIT 18.06SC sur les valeurs propres partent de cette relation : la droite engendrée par $v$ reste invariante. Si $\lambda>0$, le sens est conservé ; si $\lambda<0$, il est inversé ; si $\lambda=0$, le vecteur est envoyé sur zéro. Le vecteur nul est exclu, car il vérifierait l’égalité pour toute valeur de $\lambda$.

Prenons

$$B=\begin{bmatrix}2&1\\1&2\end{bmatrix}.$$

Pour que $(B-\lambda I)v=0$ admette une solution non nulle, $B-\lambda I$ doit être singulière. Résolvons d’abord l’équation caractéristique :

$$\det(B-\lambda I)=(2-\lambda)^2-1 =\lambda^2-4\lambda+3=(\lambda-3)(\lambda-1)=0.$$

Nous obtenons $\lambda_1=3$ et $\lambda_2=1$. Calculons ensuite les noyaux correspondants :

• Pour $\lambda=3$, $B-3I=\begin{bmatrix}-1&1\\1&-1\end{bmatrix}$, donc les coordonnées vérifient $v_2=v_1$.
• Pour $\lambda=1$, $B-I=\begin{bmatrix}1&1\\1&1\end{bmatrix}$, donc les coordonnées vérifient $v_2=-v_1$.

Après normalisation, choisissons

$$q_1=\frac{1}{\sqrt2}\begin{bmatrix}1\\1\end{bmatrix},\qquad q_2=\frac{1}{\sqrt2}\begin{bmatrix}1\\-1\end{bmatrix}.$$

La multiplication confirme $Bq_1=3q_1$ et $Bq_2=q_2$. Tout vecteur $x=c_1q_1+c_2q_2$ vérifie donc $Bx=3c_1q_1+c_2q_2$ : la première coordonnée est multipliée par trois et la seconde conserve sa longueur.

Une matrice réelle $n\times n$ s’écrit sous la forme $B=P\Lambda P^{-1}$ avec $P,\Lambda$ réelles si et seulement si elle possède $n$ vecteurs propres réels linéairement indépendants. Une matrice réelle peut avoir des valeurs propres complexes ; l’interprétation par une droite invariante s’applique aux vecteurs propres réels. Par exemple, $J=\begin{bmatrix}1&1\\0&1\end{bmatrix}$ possède la valeur propre double 1, mais le noyau de $J-I$ se réduit à la droite engendrée par $(1,0)^T$. Ses vecteurs propres ne peuvent pas former une base de dimension deux.

La symétrie permet une décomposition orthogonale

Le théorème spectral présenté dans les notes du MIT sur les matrices symétriques affirme qu’une matrice réelle symétrique possède des valeurs propres réelles et une base orthonormée de vecteurs propres. Pour deux valeurs propres distinctes $\lambda_i\ne\lambda_j$, l’orthogonalité se déduit directement :

$$\lambda_j q_i^Tq_j=q_i^TBq_j=(Bq_i)^Tq_j =\lambda_i q_i^Tq_j, \qquad q_i^Tq_j=0.$$

Dans un sous-espace propre associé à une valeur propre multiple, deux vecteurs choisis arbitrairement ne sont pas nécessairement orthogonaux. On peut toutefois y choisir une base orthonormée. Avec $Q=[q_1\ q_2]$, nous avons $Q^TQ=I$ et $Q^{-1}=Q^T$ ; l’exemple s’écrit

$$B=Q\begin{bmatrix}3&0\\0&1\end{bmatrix}Q^T =3q_1q_1^T+q_2q_2^T.$$

Chaque terme à droite projette sur un axe, puis multiplie par la valeur propre correspondante. La symétrie n’impose pas des valeurs propres positives ; celles de cet exemple le sont toutes les deux.

SVD rectangulaire et dimensions des facteurs

Une matrice 3×2 envoie un vecteur de dimension deux dans un espace de dimension trois. La relation $Av=\lambda v$ ne permet donc pas de comparer directement l’entrée et la sortie. Les notes du MIT sur la SVD établissent que toute matrice réelle admet une décomposition en valeurs singulières, ou SVD. Deux ensembles de directions orthonormées relient les deux espaces :

$$A=U\Sigma V^T,\qquad Av_i=\sigma_i u_i.$$

Le vecteur singulier droit $v_i$ appartient à l’espace d’entrée ; le vecteur singulier gauche $u_i$ appartient à l’espace de sortie. Les valeurs singulières sont ordonnées selon $\sigma_1\ge\cdots\ge\sigma_p\ge0$, où $p=\min(m,n)$. La multiplication par $V^T$ extrait les coordonnées d’entrée, $\Sigma$ les étire et $U$ les exprime dans l’espace de sortie.

Forme, $A\in\mathbb R^{m\times n}$	$U$	$\Sigma$	$V^T$
SVD complète	$m\times m$	$m\times n$	$n\times n$
SVD économique, $p=\min(m,n)$	$m\times p$	$p\times p$	$p\times n$
Termes non nuls seulement, $r=\operatorname{rank}(A)$	$m\times r$	$r\times r$	$r\times n$

Dans la forme complète, $U,V$ sont des matrices carrées orthogonales. Dans les deux autres formes, leurs colonnes sont orthonormées, mais les matrices peuvent être rectangulaires. Une SVD économique peut encore contenir des valeurs singulières nulles : $p$ n’est pas nécessairement le rang $r$.

Le tutoriel d’algèbre linéaire de SciPy nomme les objets renvoyés U, s, Vh. s est un tableau unidimensionnel de valeurs singulières ; pour une matrice réelle, Vh représente $V^T$, et non $V$. La documentation de scipy.linalg.svd précise les dimensions économiques obtenues avec full_matrices=False.

Reconstruire une matrice avec deux composantes singulières

Prenons

$$A=\begin{bmatrix}2&0\\0&2\\-2&-2\end{bmatrix},\qquad A^TA=\begin{bmatrix}8&4\\4&8\end{bmatrix}=4B.$$

Puisque $A^TA v_i=\sigma_i^2v_i$, les vecteurs propres déjà calculés conviennent : $v_1=q_1$ et $v_2=q_2$, avec les valeurs propres 12 et 4. Ainsi,

$$\sigma_1=\sqrt{12}=2\sqrt3,\qquad \sigma_2=2.$$

Pour chaque valeur singulière non nulle, calculons $u_i=Av_i/\sigma_i$ :

$$\begin{aligned} Av_1&=\begin{bmatrix}\sqrt2\\\sqrt2\\-2\sqrt2\end{bmatrix},& u_1&=\frac{1}{\sqrt6}\begin{bmatrix}1\\1\\-2\end{bmatrix};\\[6pt] Av_2&=\begin{bmatrix}\sqrt2\\-\sqrt2\\0\end{bmatrix},& u_2&=\frac{1}{\sqrt2}\begin{bmatrix}1\\-1\\0\end{bmatrix}. \end{aligned}$$

Les deux vecteurs ont une norme égale à 1 et $u_1^Tu_2=0$. En posant $U_p=[u_1\ u_2]$, $\Sigma_p=\operatorname{diag}(2\sqrt3,2)$ et $V_p=[v_1\ v_2]$, nous obtenons une SVD économique de dimensions $(3\times2)(2\times2)(2\times2)$. Pour la forme complète, ajoutons $u_3=(1,1,1)^T/\sqrt3$ et une ligne de zéros au bas de $\Sigma$ ; $A^Tu_3=0$.

Développons le produit en produits extérieurs. Chaque composante non nulle est de rang 1 :

$$\begin{aligned} A&=\sigma_1u_1v_1^T+\sigma_2u_2v_2^T\\ &=\begin{bmatrix}1&1\\1&1\\-2&-2\end{bmatrix} +\begin{bmatrix}1&-1\\-1&1\\0&0\end{bmatrix} =\begin{bmatrix}2&0\\0&2\\-2&-2\end{bmatrix}. \end{aligned}$$

Changer simultanément les signes de $u_i,v_i$ ne change pas leur produit extérieur. Pour une matrice réelle symétrique, les valeurs singulières sont les valeurs absolues des valeurs propres. Le signe d’une valeur propre négative est porté par les orientations relatives des vecteurs singuliers gauche et droit. Les valeurs singulières ne coïncident avec les valeurs propres elles-mêmes que dans le cas semi-défini positif.

Troncature et erreur d’approximation

Pour $1\le k<r$, conservons seulement les $k$ plus grandes composantes :

$$A_k=\sum_{i=1}^{k}\sigma_i u_i v_i^T.$$

La section 5 des notes de Stanford CS168 sur l’approximation de faible rang donne le résultat d’optimalité : parmi les matrices de rang au plus $k$, $A_k$ minimise l’erreur en norme de Frobenius. Cette norme est définie par $\|M\|_F^2=\sum_{a,b}M_{ab}^2$. Les produits extérieurs distincts sont orthogonaux pour le produit scalaire correspondant, car

$$\langle u_iv_i^T,u_jv_j^T\rangle_F =(u_i^Tu_j)(v_i^Tv_j).$$

Le carré de la norme du résidu est donc la somme des carrés des valeurs singulières écartées :

$$\|A-A_k\|_F=\sqrt{\sum_{i=k+1}^{r}\sigma_i^2},\qquad \|A-A_k\|_2=\sigma_{k+1}.$$

La norme matricielle $2$ mesure l’amplification maximale de la longueur d’un vecteur d’entrée unitaire. Le résidu atteint la valeur indiquée dans la direction $v_{k+1}$. Si $k\ge r$, l’erreur de reconstruction est nulle.

Dans notre exemple, conserver la première composante donne

$$A_1=\begin{bmatrix}1&1\\1&1\\-2&-2\end{bmatrix},\qquad A-A_1=\begin{bmatrix}1&-1\\-1&1\\0&0\end{bmatrix}.$$

Ainsi, $\|A-A_1\|_F=2$, $\|A\|_F=\sqrt{12+4}=4$ et l’erreur relative en norme de Frobenius vaut $1/2$. La fraction conservée du carré de la norme est $12/(12+4)=75\%$. Conserver 75 % du carré de la norme ne signifie pas une erreur relative de 25 % en norme : le carré de l’erreur représente 25 %, et sa racine carrée donne l’erreur relative en norme.

PCA : des carrés des valeurs singulières aux variances

La documentation de scikit-learn sur la PCA décrit les composantes orthogonales qui expliquent le plus de variance dans les données centrées. Le centrage ne met pas automatiquement les caractéristiques à la même échelle.

Chaque ligne de $X\in\mathbb R^{N\times d}$ représente une observation ; les moyennes des colonnes ont été soustraites et $N>1$. Avec la convention de covariance empirique utilisant le dénominateur $N-1$, comme pour la variance expliquée dans l’API PCA, remplaçons $X$ par sa SVD complète $U\Sigma V^T$ :

$$C=\frac{X^TX}{N-1} =V\frac{\Sigma^T\Sigma}{N-1}V^T.$$

Chaque $v_i$ est donc une direction principale, de variance $\sigma_i^2/(N-1)$. Pour les $k$ premières directions, la matrice des scores est $Z=XV_k=U_k\Sigma_k$ et la reconstruction vaut $ZV_k^T=X_k$. Si les données originales n’étaient pas centrées, il faut ajouter la moyenne pour retrouver leurs coordonnées d’origine.

Les colonnes de notre matrice $A$ ont une somme nulle ; elle peut donc servir directement de $X$. Avec $N=3$, la covariance est $C=\begin{bmatrix}4&2\\2&4\end{bmatrix}$ et les variances des composantes valent 6 et 2. Les scores de la première composante sont $(\sqrt2,\sqrt2,-2\sqrt2)^T$. Leur multiplication par $v_1^T$ reconstruit $A_1$, avec une proportion de variance expliquée de $6/(6+2)=75\%$.

Clustering et réduction de dimension traite du choix des échelles, des données utilisées pour l’ajustement et de l’évaluation en aval. Les 75 % mesurent la variance conservée dans ces coordonnées ; ils ne permettent pas de conclure qu’une tâche de classification conserve assez d’information.

Moindres carrés : valeurs singulières nulles et petites

Pour $\min_x\|Ax-b\|_2^2$, les notes du MIT sur la pseudo-inverse construisent la pseudo-inverse à partir de la SVD en prenant l’inverse des valeurs singulières non nulles. La solution des moindres carrés de norme euclidienne minimale est

$$x_*=A^+b=\sum_{i=1}^{r}\frac{u_i^Tb}{\sigma_i}v_i.$$

Ce calcul projette $b$ sur l’espace des colonnes, puis annule l’étirement dans chaque direction. Toute autre solution des moindres carrés s’écrit $x_*+z$, avec $Az=0$. Les directions du noyau étant orthogonales à $x_*$, celui-ci possède la plus petite norme.

Dans notre exemple, pour $b=(1,0,0)^T$, les coefficients selon les deux directions singulières valent $1/(6\sqrt2)$ et $1/(2\sqrt2)$. Nous obtenons

$$x_*=\begin{bmatrix}1/3\\-1/6\end{bmatrix},\quad Ax_*=\begin{bmatrix}2/3\\-1/3\\-1/3\end{bmatrix},\quad b-Ax_*=\frac13\begin{bmatrix}1\\1\\1\end{bmatrix}.$$

Le résidu appartient au noyau gauche : $A^T(b-Ax_*)=0$. Comme $A$ est de rang colonne plein, les coefficients sont uniques. En la remplaçant par $A_1$, nous avons $A_1v_2=0$ : ajouter $tv_2$ à une solution ne change pas les prédictions. L’identifiabilité des coefficients exige que la matrice n’efface complètement aucune direction d’entrée.

Une petite valeur singulière non nulle pose un autre problème : la division par cette valeur amplifie les perturbations des données. Pour $A$ fixée, $\delta b=\varepsilon u_i$ entraîne $\delta x=(\varepsilon/\sigma_i)v_i$. Par exemple, avec $D=\operatorname{diag}(1,0.001)$, une perturbation de $10^{-6}$ sur la deuxième coordonnée de sortie devient une perturbation de coefficient de $0.001$.

Pour une matrice de rang colonne plein, le nombre de conditionnement spectral vaut $\kappa_2(A)=\|A\|_2\|A^+\|_2=\sigma_1/\sigma_n$. Les notes du MIT sur les méthodes numériques, cours 5, relient ces normes aux valeurs propres extrêmes de $A^TA$ et montrent que $\kappa_2(A^TA)=\kappa_2(A)^2$ : les valeurs propres sont les carrés des valeurs singulières. Dans notre exemple, ces nombres valent respectivement $\sqrt3$ et 3 ; celui de $D$ vaut 1000. L’analyse numérique distingue la sensibilité du problème de la stabilité de l’algorithme ; les moindres carrés ordinaires et la régularisation approfondissent le choix du solveur et l’interprétation des coefficients.

Tronquer une SVD pour réduire la dimension contrôle l’erreur de reconstruction de la matrice. Tronquer une pseudo-inverse pour résoudre un problème écarte des directions de coefficients et modifie l’estimation. Le seuil dépend de l’erreur des données et de l’usage prévu ; une petite valeur singulière non nulle n’est pas un zéro mathématique.

Vérifier les composantes avec Python

Ce code utilise uniquement la bibliothèque standard. Il vérifie les vecteurs propres et les composantes singulières calculés ci-dessus, puis calcule la reconstruction, l’erreur de rang un et les coefficients des moindres carrés. Les fonctions auxiliaires prennent des listes : les deux vecteurs d’un produit scalaire doivent avoir la même longueur, et chaque ligne de la matrice doit avoir la longueur du vecteur d’entrée. Des dimensions incompatibles déclenchent ValueError.

from math import sqrt, isclose


def dot(x, y):
    if len(x) != len(y):
        raise ValueError("Vector lengths must match")
    return sum(a * b for a, b in zip(x, y))


def mv(matrix, vector):
    return [dot(row, vector) for row in matrix]


B = [[2, 1], [1, 2]]
vectors = [[1 / sqrt(2), 1 / sqrt(2)],
           [1 / sqrt(2), -1 / sqrt(2)]]
for eigenvalue, v in zip([3, 1], vectors):
    assert all(isclose(a, eigenvalue * b, abs_tol=1e-12)
               for a, b in zip(mv(B, v), v))

A = [[2, 0], [0, 2], [-2, -2]]
s = [sqrt(12), 2.0]
left = [[a / sigma for a in mv(A, v)]
        for sigma, v in zip(s, vectors)]
assert isclose(dot(left[0], left[1]), 0, abs_tol=1e-12)
assert all(isclose(dot(u, u), 1) for u in left)
components = [[[sigma * u[i] * v[j] for j in range(2)]
               for i in range(3)]
              for sigma, u, v in zip(s, left, vectors)]
reconstructed = [[sum(c[i][j] for c in components) for j in range(2)]
                 for i in range(3)]
assert all(isclose(reconstructed[i][j], A[i][j], abs_tol=1e-12)
           for i in range(3) for j in range(2))
A1 = components[0]
error = sqrt(sum((A[i][j] - A1[i][j]) ** 2
                 for i in range(3) for j in range(2)))
b = [1, 0, 0]
x = [sum(v[j] * dot(u, b) / sigma
         for sigma, u, v in zip(s, left, vectors)) for j in range(2)]
print("eigenvalues:", [3, 1])
print("singular values:", [round(t, 6) for t in s])
print("reconstructed:", [[round(t, 6) for t in row] for row in reconstructed])
print("rank-one:", [[round(t, 6) for t in row] for row in A1])
print(f"Frobenius error: {error:.6f}")
print(f"retained variance: {s[0] ** 2 / sum(t * t for t in s):.6f}")
print("least-squares x:", [round(t, 6) for t in x])

Sortie :

eigenvalues: [3, 1]
singular values: [3.464102, 2.0]
reconstructed: [[2.0, 0.0], [0.0, 2.0], [-2.0, -2.0]]
rank-one: [[1.0, 1.0], [1.0, 1.0], [-2.0, -2.0]]
Frobenius error: 2.000000
retained variance: 0.750000
least-squares x: [0.333333, -0.166667]

Mes lectures · Les notes ouvertes sont mémorisées dans ce navigateur.