La PCA comme moindres carrés

Le modèle linéaire du cours suppose la matrice $\mathbf{H}$ connue : ses colonnes sont choisies, puissances, sinusoïdes, indicatrices. Ce projet pose la question suivante : que se passe-t-il quand $\mathbf{H}$ est inconnue, quand les données doivent découvrir elles-mêmes le sous-espace qui les explique ? La réponse est l'analyse en composantes principales (principal component analysis, PCA), l'outil standard de la réduction de dimension, qui n'est donc pas une méthode nouvelle : ce sont les moindres carrés du chapitre 2, dont la base est devenue une inconnue de plus.

Le modèle s'écrit avec les seuls symboles du cours : $M$ observations en dimension $d$, rangées en colonnes d'une matrice $\mathbf{X}$, vivant près d'un sous-espace de dimension $k\ll d$ :

$$\mathbf{X}=\mathbf{H}\mathbf{\Theta }+\mathbf{W},$$

où

• $\mathbf{H}$ ($d\times k$) porte le sous-espace, désormais inconnu,
• $\mathbf{\Theta }=[{\mathit{\theta }}^{(1)}\cdots {\mathit{\theta }}^{(M)}]$ rassemble en colonnes les paramètres de chaque observation,
• $\mathbf{W}$ est le bruit.

1. Lien avec le cours

Le projet prolonge directement le modèle linéaire : la solution des moindres carrés sert deux fois par itération, le projecteur orthogonal $\mathbf{P}=\mathbf{H}({\mathbf{H}}^T\mathbf{H}{)}^{-1}{\mathbf{H}}^T$ devient l'objet estimé, et la matrice de covariance du chapitre 1 fait le pont avec la formulation classique de la PCA.

2. Structure du code

Le cadre commun des projets s'applique : le modèle d'abord, la loss ensuite, puis l'algorithme à la main.

class LowRankModel:
    """M observations in dimension d, living near a k-dimensional subspace."""

    def __init__(self, d=10, k=2, sigma=0.1):
        ...

    def rvs(self, m, rng):
        """Simulate (H_true, Theta_true, X) with X = H @ Theta + W."""
        ...


def loss(H, Theta, X):
    """Frobenius misfit ||X - H @ Theta||^2."""
    ...


def als(X, k, n_iter, rng):
    """Alternating least squares: solve Theta given H, then H given Theta."""
    ...

Tests imposés. Sans bruit et pour des données exactement de rang $k$, la loss finale de als est numériquement nulle et le sous-espace est retrouvé (les projecteurs coïncident) ; la loss de als décroît à chaque demi-étape, sans exception, chaque demi-étape étant un moindres carrés exact ; et les erreurs de reconstruction de als, de la SVD et de sklearn.decomposition.PCA coïncident sur les mêmes données centrées.

3. Travail demandé

1. Base connue. Pour $M$ observations partageant la même $\mathbf{H}$ connue, montrer que l'estimation de chaque ${\mathit{\theta }}^{(i)}$ est la solution des moindres carrés du cours, appliquée colonne par colonne, et qu'elle s'écrit d'un bloc $\hat{\mathbf{\Theta }}=({\mathbf{H}}^T\mathbf{H}{)}^{-1}{\mathbf{H}}^T\mathbf{X}$.
2. Base inconnue. Poser le problème joint $\arg \underset{\mathbf{H},\mathbf{\Theta }}{min}\Vert \mathbf{X}-\mathbf{H}\mathbf{\Theta }{\Vert }^2$. Montrer que pour toute matrice inversible $\mathbf{A}$, le couple $(\mathbf{H}\mathbf{A},{\mathbf{A}}^{-1}\mathbf{\Theta })$ donne la même loss : la base n'est pas identifiable, seul le sous-espace l'est, et le projecteur $\mathbf{P}$ est l'objet bien défini. Convenir d'une représentation, colonnes orthonormées.
3. Les moindres carrés alternés, à la main. Implémenter als : à $\mathbf{H}$ fixée, résoudre $\mathbf{\Theta }$ (la formule de l'étape 1) ; à $\mathbf{\Theta }$ fixée, résoudre $\mathbf{H}$ (la même formule, transposée). Vérifier la décroissance monotone de la loss, choisir un critère d'arrêt, et étudier l'effet de l'initialisation aléatoire : la loss finale change-t-elle ? le sous-espace ?
4. Le bon outil. Centrer les données, soustraire la moyenne empirique, c'est estimer $\mathit{\mu }$, puis obtenir le sous-espace optimal d'un coup par la SVD de $\mathbf{X}$ centrée, les $k$ premiers vecteurs singuliers gauches (théorème d'Eckart-Young, admis), et confronter à sklearn.decomposition.PCA. Comparer à als : erreurs de reconstruction, distance entre projecteurs $\Vert {\mathbf{P}}_{\mathrm{ALS}}-{\mathbf{P}}_{\mathrm{SVD}}\Vert$.
5. Interprétation. Tracer l'erreur de reconstruction en fonction de $k$, de $1$ à $d$ : la courbe en coude, l'outil de choix de la dimension. Pour $k=2$, tracer le nuage des colonnes de $\hat{\mathbf{\Theta }}$, la représentation réduite des données, et interpréter les colonnes de $\hat{\mathbf{H}}$, les directions que les données ont élues.
6. Ouverture : la PCA probabiliste. Le geste bayésien du cours s'applique aussi ici : la pPCA (Tipping et Bishop) modélise $\mathbf{x}=\mathbf{H}\mathit{\theta }+\mathit{\mu }+\mathbf{w}$ avec $\mathit{\theta }\sim \mathcal{N}(\mathbf{0},\mathbf{I})$ et $\mathbf{w}\sim \mathcal{N}(\mathbf{0},{\sigma }^2\mathbf{I})$. Son MLE retrouve le sous-espace principal, avec en prime ${\hat{\sigma }}^2$ = moyenne des valeurs propres écartées de la covariance empirique (résultat admis). Le vérifier par Monte-Carlo sur des données simulées où $\sigma$ est connue. Remarquer enfin que l'alternance de als, une affectation puis une mise à jour, est aussi la structure de K-means : la réduction de dimension et le partitionnement sont deux argmin alternés.

4. Livrables

Le notebook reproductible des modalités communes, avec la démonstration de non-identifiabilité rédigée, la courbe de décroissance de la loss de als, la comparaison chiffrée ALS/SVD/sklearn, la courbe en coude et le nuage réduit commentés, et la vérification pPCA de l'ouverture.