Pourquoi PCA de données au moyen de SVD des données?

Cette question concerne un moyen efficace de calculer les principaux composants.

De nombreux textes sur l'ACP linéaire préconisent l'utilisation de la décomposition en valeurs singulières des données casewise . Autrement dit, si nous avons des données et que nous voulons remplacer les variables (ses colonnes ) par des composants principaux, nous faisons SVD: , des valeurs singulières (racines carrées des valeurs propres) occupant la diagonale principale de , vecteurs propres droits sont la matrice de rotation orthogonale des variables-axes en composantes-axes, les vecteurs propres gauche sont comme , seulement pour les cas. Nous pouvons alors calculer les valeurs des composants comme . $\bf X$ $\bf X=USV'$ $\bf S$ $\bf V$ $\bf U$ $\bf V$ $\bf C=XV=US$
Une autre façon de réaliser l'ACP des variables consiste à décomposer la matrice carrée $\bf R=X'X$ (c'est-à-dire que $\bf R$ peut être des corrélations ou des covariances, etc., entre les variables). La décomposition peut être une décomposition propre ou une décomposition en valeurs singulières: avec une matrice semi-définie positive symétrique carrée, ils donneront le même résultat $\bf R=VLV'$ avec des valeurs propres comme la diagonale de $\bf L$ et $\bf V$ comme décrit précédemment. Les valeurs des composants seront $\bf C=XV$ .

Maintenant, ma question: si les données $\bf X$ sont une grande matrice et que le nombre de cas est (ce qui est souvent un cas) beaucoup plus grand que le nombre de variables, alors la voie (1) devrait être beaucoup plus lente que la voie (2), parce que way (1) applique un algorithme assez cher (tel que SVD) à une grande matrice; il calcule et stocke une énorme matrice $\bf U$ dont nous n'avons vraiment pas besoin dans notre cas (le PCA des variables). Si c'est le cas , pourquoi tant de texbooks semblent-ils préconiser ou simplement mentionner la seule façon (1)? Peut-être qu'il est efficace et je manque quelque chose?

pca algorithms svd matrix-decomposition ttnphns
la source

En général, nous nous intéressons uniquement aux quelques composantes principales qui expliquent la majeure partie de la variance. Il est possible de faire une SVD réduite; par exemple , si

est de dimension

où

puis de » la fonction calcule uniquement la première

gauche et droite vecteurs singuliers par défaut.

X

$X$

N \times p

$N \times p$

p << N

$p << N$ Rsvd

p

$p$

M. Berk

@ M.Berk:

cependant le même dans les deux approches: ils donnent des résultats équivalents (égaux à des changements de signe). De plus, par exemple, R ne calcule

que sur demande.

p

$p$

C

$\mathbf C$

cbeleites prend en charge Monica

Avez-vous des références pour le chemin (1)? Je ne suis conscient que de la mise en œuvre de l'ACP via SVD sur la matrice de covariance (c'est-à-dire de la manière 2), car cela évite certains problèmes numériques et s'adapte évidemment à la dimensionnalité, pas à la taille de l'ensemble de données. Way (1) J'appellerais SVD, pas PCA du tout. Je ne l'ai vu que dans un contexte SVD pur, où l'on n'effectuerait pas une décomposition complète en fait.

Anony-Mousse -Restate Monica

@ Anony-Mousse, juste une Joliffe, Principal component analysis, 2nd ed.chose à mentionner, en fait, Joliffe décrit les deux façons, mais dans le chapitre central sur PCA, il dit à propos de la voie 1, pour autant que je me souvienne.

ttnphns

@ Anony-Mousse, la voie 1 est pour moi importante du point de vue théorique car elle montre clairement comment l'ACP est directement liée à l' analyse de correspondance simple .

ttnphns

Réponses:

Voici mon 2ct sur le sujet

La conférence de chimiométrie où j'ai d'abord appris la PCA a utilisé la solution (2), mais elle n'était pas orientée vers le numérique, et ma conférence de numérique n'était qu'une introduction et n'a pas discuté de SVD pour autant que je m'en souvienne.
Si je comprends bien Holmes: Fast SVD for Large-Scale Matrices , votre idée a été utilisée pour obtenir un SVD rapide sur le plan des matrices longues.
Cela signifierait qu'une bonne implémentation SVD peut suivre en interne (2) si elle rencontre des matrices appropriées (je ne sais pas s'il existe encore de meilleures possibilités). Cela signifierait que pour une implémentation de haut niveau, il est préférable d'utiliser le SVD (1) et de laisser au BLAS le soin de choisir l'algorithme à utiliser en interne.

Vérification pratique rapide: le svd d'OpenBLAS ne semble pas faire cette distinction, sur une matrice de 5e4 x 100, svd (X, nu = 0)prend 3,5 s médian, tout svd (crossprod (X), nu = 0)en 54 ms (appelé depuis R avec microbenchmark).
La quadrature des valeurs propres est bien sûr rapide, et jusqu'à ce que les résultats des deux appels soient équivalents.

timing  <- microbenchmark (svd (X, nu = 0), svd (crossprod (X), nu = 0), times = 10)
timing
# Unit: milliseconds
#                      expr        min         lq    median         uq        max neval
#            svd(X, nu = 0) 3383.77710 3422.68455 3507.2597 3542.91083 3724.24130    10
# svd(crossprod(X), nu = 0)   48.49297   50.16464   53.6881   56.28776   59.21218    10

mise à jour: Jetez un œil à Wu, W .; Massart, D. & de Jong, S .: Les algorithmes PCA du noyau pour les données étendues. Partie I: Théorie et algorithmes, Chemometrics and Intelligent Laboratory Systems, 36, 165 - 172 (1997). DOI: http://dx.doi.org/10.1016/S0169-7439(97)00010-5

Cet article discute des propriétés numériques et informatiques de 4 algorithmes différents pour l'ACP: SVD, décomposition propre (EVD), NIPALS et POWER.

Ils sont liés comme suit:

computes on      extract all PCs at once       sequential extraction    
X                SVD                           NIPALS    
X'X              EVD                           POWER

Le contexte de l'article est large , et ils fonctionnent sur (noyau PCA) - c'est juste la situation opposée à celle que vous demandez. Donc, pour répondre à votre question sur le comportement à longue matrice, vous devez échanger la signification de "noyau" et "classique". $\mathbf X^{(30 \times 500)}$ $\mathbf{XX'}$

comparaison des performances

Sans surprise, EVD et SVD changent de place selon que les algorithmes classiques ou noyau sont utilisés. Dans le contexte de cette question, cela signifie que l'un ou l'autre peut être meilleur selon la forme de la matrice.

Mais à partir de leur discussion sur les SVD et EVD "classiques", il est clair que la décomposition de est un moyen très courant de calculer l'APC. Cependant, ils ne spécifient pas quel algorithme SVD est utilisé autre que celui qu'ils utilisent la fonction de Matlab . $\mathbf{X'X}$ svd ()

    > sessionInfo ()
    R version 3.0.2 (2013-09-25)
    Platform: x86_64-pc-linux-gnu (64-bit)

    locale:
     [1] LC_CTYPE=de_DE.UTF-8       LC_NUMERIC=C               LC_TIME=de_DE.UTF-8        LC_COLLATE=de_DE.UTF-8     LC_MONETARY=de_DE.UTF-8   
     [6] LC_MESSAGES=de_DE.UTF-8    LC_PAPER=de_DE.UTF-8       LC_NAME=C                  LC_ADDRESS=C               LC_TELEPHONE=C            
    [11] LC_MEASUREMENT=de_DE.UTF-8 LC_IDENTIFICATION=C       

    attached base packages:
    [1] stats     graphics  grDevices utils     datasets  methods   base     

    other attached packages:
    [1] microbenchmark_1.3-0

loaded via a namespace (and not attached):
[1] tools_3.0.2

$ dpkg --list libopenblas*
[...]
ii  libopenblas-base              0.1alpha2.2-3                 Optimized BLAS (linear algebra) library based on GotoBLAS2
ii  libopenblas-dev               0.1alpha2.2-3                 Optimized BLAS (linear algebra) library based on GotoBLAS2

cbeleites soutient Monica
la source

Ainsi, vos tests (3,5 secondes contre 54 ms) soutiennent ma ligne que la "voie 1" est considérablement plus lente. Droite?

ttnphns

@ttnphns: oui. Mais comme svd est fourni par le BLAS, cela pourrait être différent avec un BLAS différent. Je m'attendais à ce qu'un bon BLAS optimisé fasse quelque chose comme ça. Cependant, cela ne semble pas être le cas avec OpenBLAS. Je suis trop paresseux pour vérifier d'autres BLAS, mais peut-être que quelques personnes pourraient vérifier leurs autres BLAS afin que nous découvrions ceux qui sont optimisés pour ce cas et ceux qui ne le sont pas. (J'ai envoyé un e-mail au développeur OpenBLAS et je lui ai envoyé un lien vers cette question, alors peut-être qu'il pourrait ajouter des informations, par exemple des raisons pour ne pas basculer l'algorithme vers les svd (X'X)matrices longues.)

cbeleites prend en charge Monica

X^{'}

$X'$

X

$X$

n < p

$n<p$

X^{'} X

$X'X$

u_{n + 1} = X^{'} X u_{n} / | | X^{'} X u_{n} | |

$u_{n+1} = X'Xu_n / ||X'Xu_n||$

v_{1}

$v_1$

X

$X$

X^{'} X

$X'X$

X^{'} \times (X u_{n})

$X' \times (X u_n)$

X X^{T}

$\mathbf X \mathbf X^T$

Je parlais de votre mise à jour, où Nipals est impliqué. Je confirme que Nipals n'est pas impliqué dans le SVD de Lapack. À propos de votre expérience de référence, quelque chose comme ça microbenchmark(X <- matrix(rnorm(5e6), ncol=100), Y <- t(X), svd(X), svd(Y), control=list(order="inorder"), times = 5)peut aussi être intéressant.

Elvis

La SVD est plus lente mais est souvent considérée comme la méthode préférée en raison de sa précision numérique plus élevée.

$\mathbf X$ $\frac{1}{n-1}\mathbf X^\top \mathbf X$ $\mathbf{XX}^\top$ $n\ll p$ ).

Voici ce qui est écrit dans l' aide de lapca() fonction de MATLAB :

Algorithme de composant principal pca utilise pour effectuer l'analyse du composant principal [...]:

'svd' - Par défaut. Décomposition en valeurs singulières (SVD) de X.

$n$ $p$

La dernière phrase met en évidence le compromis crucial entre vitesse et précision qui est en jeu ici.

$1000\times 100$

X = randn([1000 100]);

tic; svd(X); toc         %// Elapsed time is 0.004075 seconds.
tic; svd(X'); toc        %// Elapsed time is 0.011194 seconds.
tic; eig(X'*X); toc      %// Elapsed time is 0.001620 seconds.
tic; eig(X*X'); toc;     %// Elapsed time is 0.126723 seconds.

$n \ll p$ $\mathbf{XX}^\top$ (quatrième ligne) sera le moyen le plus rapide à la place. La SVD de la matrice de données elle-même sera plus lente dans les deux cas.

$\mathbf X$ $X$ $XX^⊤$

X = (\begin{matrix} 1 & 1 & 1 \\ ϵ & 0 & 0 \\ 0 & ϵ & 0 \\ 0 & 0 & ϵ \end{matrix}),

$\mathbf X = \begin{pmatrix}1&1&1\\\epsilon & 0 & 0\\ 0 & \epsilon & 0 \\ 0 & 0 & \epsilon\end{pmatrix},$

3 + ϵ^{2}

$3+\epsilon^2$

ϵ^{2}

$\epsilon^2$

ϵ^{2}

$\epsilon^2$

ϵ = 10^{- 5}

$\epsilon = 10^{-5}$

eps = 1e-5;
X = [1 1 1; eye(3)*eps];
display(['Squared sing. values of X: ' num2str(sort(svd(X),'descend').^2')])
display(['Eigenvalues of X''*X:       ' num2str(sort(eig(X'*X),'descend')')])

obtenir des résultats identiques:

Squared sing. values of X: 3       1e-10       1e-10
Eigenvalues of X'*X:       3       1e-10       1e-10

$\epsilon = 10^{-10}$

Squared sing. values of X: 3       1e-20       1e-20
Eigenvalues of X'*X:       3           0 -3.3307e-16

$\mathbf X$ $\mathbf X$ a des colonnes presque colinéaires (c.-à-d. De très petites valeurs singulières), le calcul de la matrice de covariance puis le calcul de sa composition par eigend résultera en perte de précision par rapport à la SVD directe.

Je dois ajouter que l'on est souvent heureux d'ignorer cette perte de précision [minuscule] potentielle et d'utiliser plutôt la méthode plus rapide.

amibe dit réintégrer Monica
la source

Aussi pour le SVD, vous pouvez écrire (ou simplement utiliser) l'algorithme pour

X^{T}

$\mathbf X^T$ qui a le comportement d'exécution opposé (car la décomposition a une symétrie par rapport à la transposition

X

$\mathbf X$ ) et sinon le même comportement numérique. La stabilité numérique est un point important (+1) - bien que je suppose que cela dépend beaucoup des données, que ce soit important. Par exemple, je suis généralement trop peu de cas et de mesures bruyantes - donc la stabilité de mes modèles n'est généralement pas limitée par la stabilité numérique des SVD sous-jacentes, mais par les lots de culture de patients / cellules ou le bruit instrumental.

cbeleites prend en charge Monica

Merci pour la réponse et pour un examen approfondi des avantages et des inconvénients.

ttnphns

amibe, est-ce que vous trouvez le temps de montrer un exemple concret où la stabilité numérique souffre de l' eig()approche? (Les lecteurs en bénéficieront: il y a un compromis entre vitesse et stabilité. Comment décider dans une situation concrète?)

ttnphns

@ttnphns J'ai réécrit toute la réponse, en donnant un exemple concret. Regarde.

amibe dit Réintégrer Monica

@amoeba, merci beaucoup d'être revenu et de donner un exemple! J'ai essayé les deux exemples epsilon dans SPSS et j'ai obtenu des résultats comme les vôtres, sauf pour la toute dernière ligne: au lieu de 3 0 -3.3307e-16eigen dans spss, je suis retourné 3 0 0. Il semble que la fonction ait une certaine valeur de tolérance intégrée et fixe au-delà de laquelle elle se réinitialise. Dans cet exemple, la fonction est apparue comme pour pirater le nœud de l'instabilité numérique en mettant à zéro les deux petites valeurs propres, le «0» et le «-16».

ttnphns