Les maths derrière l'ACP

Cette page explique le fonctionnement de l'ACP dans le cas d'un système multidimensionel avec une distribution gaussienne. Nous allons expliquer comment calculer l'ACP d'un nuage de points. Précisément, nous allons calculer :

le centre (la moyenne) du nuage de points
la matrice de covariance (ou l'écart-type)
l'angle de rotation (ou la matrice de rotation)

Le script Matlab utilisé sur cette page peut être téléchargé depuis la section Téléchargement de cette page.

Nuage de points

Considérons le nuage de points suivant :

Ce nuage de $ N $ points a été tiré au hasard avec les propriétés suivantes :

Le centre du nuage est donné par :

$$ C = \begin{bmatrix} 8 && 2 \end{bmatrix} $$

La covariance du nuage de points est données par:

$$ \Sigma= \begin{bmatrix} 3 && 0\\ 0 && 1 \end{bmatrix} $$

L'angle de rotation du nuage est : $ \alpha = \dfrac{\pi}{3} $

Consultez cette page pour obtenir davanatage de détail sur les covariances. L'objectif de l'ACP est de calculer les paramètres ci-dessous à partir des coordonnées des points :

Dans la suite, nous allons considérer que les points auront les coordonnées suivantes :

$$ X = \begin{bmatrix} x_1 && y_1 \\ x_2 && y_2 \\ x_3 && y_3 \\ ... && ...\\ x_n && y_n \end{bmatrix} $$

Le script Matlab ci-dessous permet de générer le nuage de points :

% Number of sample 
N=1000;

% Standard deviation along X and Y
StdX=3;
StdY=1;

% Offset along X and Y
offsetX = 8;
offsetY = 2;

% Rotation of sample (correlation)
alpha=pi/3;

%% Create and display random sample
X = randn(N,2)
X = X*[StdX,0;0,StdY]';
X = X*[cos(alpha) , -sin(alpha) ; sin(alpha) , cos(alpha) ]';
X = X+[offsetX offsetY];

Centre

Calculer le centre du nuage de points est certainement la partie la plus facile. Le centre est la moyenne des coordonnées $ x $ et $ y $ :

$$ \hat{C} = \begin{bmatrix} x_c && y_c \end{bmatrix} = \begin{bmatrix} \sum\limits_{i=1}^n x_i && \sum\limits_{i=1}^n y_i \end{bmatrix} $$

Le code Matlab ci-dessous calcule le centre utilisé dans la suite du calcul de l'ACP :

>> mean(X)

ans =

    7.8775    1.9065

Le centre est ici proche du résultat escompté : $ \begin{bmatrix} 8 && 2 \end{bmatrix} $.

Covariance

L'étape suivante consiste à calculer la matrice de covariance du nuage de points centré. Avant de calculer la matrice de covariance, nous devons centrer les points. Définissons $ X_c $ les coordonnées des points après translation (centrage en 0,0).

$$ X_c = \begin{bmatrix} x_1-x_c && y_1-y_c \\ x_2-x_c && y_2-y_c \\ x_3-x_c && y_3-y_c \\ ... && ...\\ x_n-x_c && y_n-y_c \end{bmatrix} $$

La matrice de covariance est obtenue en calculant le produit suivant :

$$ \Sigma = \dfrac{1}{N} \times X_c^\top \times X_c $$

Ou encore à partir des équations suivantes :

$$ \Sigma = \begin{bmatrix} \dfrac{1}{N} \sum\limits_{i=1}^N (x_i-x_c)^2 && \dfrac{1}{N} \sum\limits_{i=1}^N (x_i - x_c)(y_i - y_c) \\ \dfrac{1}{N} \sum\limits_{i=1}^N (x_i - x_c)(y_i - y_c) && \dfrac{1}{N}\sum\limits_{i=1}^N (y_i-y_c)^2 \end{bmatrix} $$

Avec Matlab, la covariance peut être directement calculée avec la fonction cov() :

>> cov (X-mean(X))

ans =

    3.0458    3.7081
    3.7081    7.5771

Décomposition en valeurs singulières

La matrice de covariance contient à la fois la matrice de rotation et les écart-types. Ces données peuvent être extraites grâce à la décomposition en valeurs singulières (SVD). Le résultat de la décomposition en valeurs singulière est donné ci-dessous :

$$ \Sigma=R.(S.S).R^T $$

avec :

$R$ est la matrice de rotation (vecteurs propres);
$S$ est une matrice d'homothétie (le carré des valeurs propres).

Pour comprendre cette décomposition, je vous invite à lire cette section. Je ne vais pas expliquer le calcul de la SVD, car je l'ai déjà détaillé sur cette page: Décomposition en valeurs singulières (SVD) d’une matrice 2×2.

Avec Matlab, la meilleure option est d'utiliser la fonction svd() :

>> sigma = cov (X-mean(X));
>> [U,S,D] = svd(sigma)

U =

   -0.4892   -0.8722
   -0.8722    0.4892

S =

    9.6569         0
         0    0.9660

D =

   -0.4892   -0.8722
   -0.8722    0.4892

$ U $ est la matrice de rotation et $S$ contient les écarts-types.

Écarts-types

Lorsque la décomposition en valeurs singulières est réalisée, il devient relativement facile d'obtenir les écarts-types le long de chaque axe en calculant :

$$ \Sigma = \begin{bmatrix} \sigma_x && 0 \\ 0 && \sigma_y \end{bmatrix} = \sqrt(S) $$

$ \sigma_x $ et $ \sigma_y $ sont les écarts-types le long de chaque axe. Ils peuvent aussi être vus comment des facteurs ou des coefficients d'homothétie pour chaque axes.

Vérifions avec Matlab si nos écarts-types sont corrects :

>> sqrt(S)

ans =

    3.1076         0
         0    0.9829

Comme les valeurs attendues sont 3 et 1, le résultat est suffisament proche pour être validé.

Rotation

La paramètre suivant que nous allons calculer est l'angle de rotation, la matrice de rotation ayant déjà été calculée : c'est la matrice $ U $ de la décomposition en valeurs singulières. Nous pouvons déduire l'angle de rotation à partir de cette matrice.

La matrice de rotation est en réalité une matrice de rotation autour de l'axe Z. Cette matrice est décrite par l'expression suivante :

$$ U = = \begin{bmatrix} U_{11} && U_{12} \\ U_{21} && U_{22} \end{bmatrix} = R_z(\alpha) = \begin{bmatrix} cos(\gamma) && -sin(\gamma \\ sin(\gamma) && cos(\gamma) \end{bmatrix} $$

L'angle $ \alpha$ peut être calculé à partir des coefficients de la matrice. Comme la majorité des langages de programmation possède une fonction atan2(), le meilleur moyen d'obtenir l'angle est de calculer :

$$ \alpha = atan2(U_{12}, U_{11}) $$

Vérifions avec Matlab :

>> alpha = atan2(U(1,2), U(1,1))

alpha =

   -2.0847

L'angle attendu est $ \dfrac{\pi}{3} \approx 1.0472 $. Comme l'ACP est symétrique il y a deux options pour l'angle : $ \dfrac{\pi}{3} $ et $ \dfrac{\pi}{3} - \pi \approx -2.0944$. la seconde solution est celle calculée par Matlab.

Vecteurs de l'ACP

Il peut parfois être utile de calculer les vecteurs représentant l'ACP comme sur cet exemple :

Pour obtenir le premier vecteur, il suffit de tourner et de mettre à l'échelle un vecteur unité le long de l'axe X :

$$ \vec{u} = U \times \sqrt{S} \times \begin{bmatrix} 1 \\ 0 \end{bmatrix} $$

Pour obtenir le second vecteur, il suffit de tourner et de mettre à l'échelle un vecteur unité le long de l'axe Y :

$$ \vec{v} = U \times \sqrt{S} \times \begin{bmatrix} 0 \\ 1 \end{bmatrix} $$

Avec Matlab, nous obtenons les vecteurs suivants :

>> u=U*sqrt(S)*[1 ; 0]

u =

   -1.4885
   -2.6012

>> v=U*sqrt(S)*[0 ; 1]

v =

   -0.8952
    0.5123

Le vecteur $ \vec{u} $ est affiché en rouge sur la figure ci-dessus, $ \vec{v} $ est affiché en vert. La superposition des vecteurs avec le nuage de points confirme nos résultats.

Téléchargement

Ci-dessous, vous pouvez télécharger le script utilisé pour créer les figures de cette page et la vérification des équations :

principal_component_analysis.m

Voir aussi

Dernière mise à jour : 24/11/2022