from numpy import *
# This line configures matplotlib to show figures embedded in the notebook,
%matplotlib inline
import matplotlib.pyplot as plt
import pylab as pl # import MATLAB-like API
On génère un nuage de $N$ points dans $\mathbb{R^2}$, $X \in \mathbb{R}^{N \times 2}$, avec $N = 400$. Chaque ligne de $X$ représente un point dans $\mathbb{R}^2$.
N = 400
X = random.randn(N,2) # distribution normale (0,1) avec tirages indépendants
T = array([[-1,1],[4,2]]) # matrice quelconque 2x2
X = dot(X,T) # on transforme X linéairement X*T
Prétraitement des données $X$. On recentre les colonnes de $X$ avec des moyennes en zéro: $X_c = (X - \bar X)$, avec $\bar X$ la moyenne de $X$ prise sur chaque colonne.
Xc = zeros(X.shape)
for k in range(0,2):
Xc[:,k] = X[:,k] - mean(X[:,k])
Représentation des données originales $X_c$.
pl.plot([-30, 30],[0, 0],'k')
pl.plot([0, 0],[-20, 20],'k')
pl.plot(Xc[:,0],Xc[:,1],'b.')
pl.plot(Xc[1,0],Xc[1,1],'ro') # X[1,:] tagué rouge
pl.plot(Xc[8,0],Xc[8,1],'co') # X[8,:] tagué cyan
pl.axis('equal')
pl.axis([-10,10,-10,10])
[-10, 10, -10, 10]
Décomposition en valeurs singulières: $X = U \Sigma V^*$.
U, s, Vs = linalg.svd(Xc, full_matrices=True)
La matrice Vs est la matrice adjointe de V: Vs = $V^*$, et s est un vecteur, avec $\Sigma$ = diag(s).
s, Vs
(array([ 91.34598705, 28.05125181]), array([[-0.90078622, -0.43426281], [-0.43426281, 0.90078622]]))
Le nuage de points se décompose $X_c = U \Sigma V^*$. Le long des composantes principales $P = X_c V = U\Sigma$. $P$ est $X_c$ dans la base donnée par les colonnes de $V$. Pour le voir, on prend l'adjoint de $X_c$ pour exprimer avoir un ensemble de $N$ vecteurs colonne. Les colonnes de $V^{-1} X_c^*$ sont alors les colonne de $X_c^*$ exprimées dans la base des vecteurs colonnes de $V$. Mais $V^{-1} X_c^* = (X_c (V^{-1})^*)^* = (X_c V)^* = P^*$, d'où les vecteurs lignes de $P$ sont les points de $X_c$ exprimés dans la base de $V$.
P = dot(Xc,Vs.T)
pl.plot(P[:,0],P[:,1],'b.')
pl.plot(array([-30,30]),array([0, 0]),'k')
pl.plot(array([0, 0]),array([-30,30]),'k')
pl.plot(P[1,0],P[1,1],'ro') # X[1,:] tagué rouge
pl.plot(P[8,0],P[8,1],'co') # X[8,:] tagué cyan
pl.axis('equal')
pl.axis([-10,10,-10,10])
[-10, 10, -10, 10]
Les valeurs singulières $\sigma_1$ et $\sigma_2$ sont les racines carrées des valeurs propres de la matrice 2x2 $C = X_c^* X_c$. Les coefficients $(i,j)$ de la matrice $C$ sont les produits scalaires $(X_i-\bar X)^*(X_j - \bar X_j)$, avec $X_i$ la colonne $i$ de $X$ et $\bar X_i$ la moyenne de la colonne $i$. Donc $C_{i,j}/(N-1)$ est une matrice de covariance. Les valeurs singulières de $X_c$, qui sont les racines carrées des valeurs propres de la matrice $C$, sont donc les écarts types à un facteur de normalisation $\sqrt{N-1}$ près.
len1 = s[0]*2/sqrt(N-1) # 2 écarts types le long du semi-axe u1
len2 = s[1]*2/sqrt(N-1) # 2 écarts types le long du semi-axe u2
(On a selectionné les longueurs égales à deux fois les écarts types parce que l'intervalle [-len1, +len1] contient environs 96% des points.)
pl.plot(P[:,0],P[:,1],'b.')
pl.plot(array([-30,30]),array([0, 0]),'k')
pl.plot(array([0, 0]),array([-30,30]),'k')
pl.plot([0,len1],[0,0],'g',linewidth=4)
pl.plot([0,0],[0,len2],'g',linewidth=4)
pl.plot(P[1,0],P[1,1],'ro') # X[1,:] tagué rouge
pl.plot(P[8,0],P[8,1],'co') # X[8,:] tagué cyan
pl.axis('equal')
pl.axis([-10,10,-10,10])
pl.xlabel('premiere composante')
pl.ylabel('deuxieme composante')
<matplotlib.text.Text at 0x107c6ced0>
On définit un cercle unité, qu'on transformera en ellipse avec semi-axes principaux donnés par les vecteurs propres de $C$, c.-à-d. les colonnes de $V$, normalisés par les facteurs len1 et len2.
theta = linspace(0,2*pi,101)
circle = zeros((101,2))
circle[:,0] = cos(theta)
circle[:,1] = sin(theta)
ellipse = dot(circle*s*2/sqrt(N-1),Vs)
pl.plot([-30, 30],[0, 0],'k') # axe x
pl.plot([0, 0],[-20, 20],'k') # axe y
pl.plot(Xc[:,0],Xc[:,1],'b.') # nuage de points
pl.plot(ellipse[:,0],ellipse[:,1],'g',linewidth=4) # ellipse avec semi-axes orientés par V et de longueurs 2*sigma/sqrt(N-1)
pl.plot([0,Vs[0,0]*len1],[0,Vs[0,1]*len1],'g',linewidth=4) # semi-axes orientés par V1 et de longueur 2*sigma1/sqrt(N-1)
pl.plot([0,Vs[1,0]*len2],[0,Vs[1,1]*len2],'g',linewidth=4) # semi-axes orientés par V2 et de longueur 2*sigma2/sqrt(N-1)
pl.plot(Xc[1,0],Xc[1,1],'ro') # X[1,:] tagué rouge
pl.plot(Xc[8,0],Xc[8,1],'co') # X[8,:] tagué cyan
pl.axis('equal') # ratio d'aspect des axes 1:1 important pour l'interprétation géométrique
pl.axis([-10,10,-10,10])
pl.xlabel('X1')
pl.ylabel('X2')
<matplotlib.text.Text at 0x107e7d110>
pl.plot(P[:,0],P[:,1],'b.')
pl.plot(array([-30,30]),array([0, 0]),'k')
pl.plot(array([0, 0]),array([-30,30]),'k')
pl.plot(circle[:,0]*len1,circle[:,1]*len2,'g',linewidth=4) # ellipse avec semi-axes le long des composantes principales
pl.plot(P[1,0],P[1,1],'ro') # X[1,:] tagué rouge
pl.plot(P[8,0],P[8,1],'co') # X[8,:] tagué cyan
pl.axis('equal')
pl.axis([-10,10,-10,10])
pl.xlabel('premiere composante')
pl.ylabel('deuxieme composante')
<matplotlib.text.Text at 0x107c65210>