J'observe un comportement très étrange dans le résultat SVD de données aléatoires, que je peux reproduire à la fois dans Matlab et R. Il ressemble à un problème numérique dans la bibliothèque LAPACK; est-ce?
Je tire échantillons de la gaussienne dimensionnelle avec une moyenne nulle et une covariance d'identité: . Je les assemble dans une matrice de données X . (Je peux éventuellement centrer X ou non, cela n'influence pas ce qui suit.) Ensuite, j'effectue une décomposition en valeurs singulières (SVD) pour obtenir X = U S V ⊤ . Prenons deux éléments particuliers de U , par exemple U 11 et U , et demandez quelle est la corrélation entre les entredifférents tirages de . Je m'attendrais à ce que si le nombre de tirages soit raisonnablement grand, alors toutes ces corrélations devraient être autour de zéro (c'est-à-dire que les corrélations de population devraient être nulles, et les corrélations d'échantillon seraient petites).
Cependant, j'observe des corrélations étrangement fortes (autour de ) entre , , et , et uniquement entre ces éléments. Toutes les autres paires d'éléments ont des corrélations autour de zéro, comme prévu. Voici à quoi ressemble la matrice de corrélation pour les éléments "supérieurs" de ( premiers éléments de la première colonne, puis les premiers éléments de la deuxième colonne):U 10 10
Remarquez des valeurs étrangement élevées dans les coins supérieurs gauches de chaque quadrant.
C'est ce commentaire de @ whuber qui a attiré mon attention sur cet effet. @whuber a fait valoir que les PC1 et PC2 ne sont pas indépendants et a présenté cette forte corrélation comme preuve de cela. Cependant, j'ai l'impression qu'il a accidentellement découvert un bogue numérique dans la bibliothèque LAPACK. Qu'est-ce qui se passe ici?
Voici le code R de @ whuber:
stat <- function(x) {u <- svd(x)$u; c(u[1,1], u[2, 2])};
Sigma <- matrix(c(1,0,0,1), 2);
sim <- t(replicate(1e3, stat(MASS::mvrnorm(10, c(0,0), Sigma))));
cor.test(sim[,1], sim[,2]);
Voici mon code Matlab:
clear all
rng(7)
n = 1000; %// Number of variables
k = 2; %// Number of observations
Nrep = 1000; %// Number of iterations (draws)
for rep = 1:Nrep
X = randn(n,k);
%// X = bsxfun(@minus, X, mean(X));
[U,S,V] = svd(X,0);
t(rep,:) = [U(1:10,1)' U(1:10,2)'];
end
figure
imagesc(corr(t), [-.5 .5])
axis square
hold on
plot(xlim, [10.5 10.5], 'k')
plot([10.5 10.5], ylim, 'k')
la source
Réponses:
Ce n'est pas un bug.
Comme nous l'avons exploré (en détail) dans les commentaires, deux choses se produisent. La première est que les colonnes deU sont contraintes pour répondre aux exigences SVD: chacune doit avoir une longueur unitaire et être orthogonale à toutes les autres. Affichage d'un U comme une variable aléatoire créée à partir d' une matrice aléatoire X via un algorithme de SVD particulier, nous notons donc que ces k ( k + 1 ) / 2 contraintes fonctionnellement indépendantes créer des dépendances statistiques entre les colonnes de U .
Ces dépendances pourraient être révélées plus ou moins en étudiant les corrélations entre les composantes deU , mais un deuxième phénomène émerge : la solution SVD n'est pas unique. Au minimum, chaque colonne de U peut être négativement indépendante, donnant au moins 2k solutions distinctes avec k colonnes. De fortes corrélations (excédant 1 / 2 ) peut être induite en changeant les signes des colonnes de manière appropriée. (Une façon de le faire est donnée dans mon premier commentaire à la réponse d'Amoeba dans ce fil: je force tous les uje je,i=1,…,k pour avoir le même signe, les rendant tous négatifs ou tous positifs avec une probabilité égale.) En revanche, toutes les corrélations peuvent être faites pour disparaître en choisissant les signes au hasard, indépendamment, avec des probabilités égales. (Je donne un exemple ci-dessous dans la section "Modifier".)
Avec des soins, nous pouvons discerner partiellement ces deux phénomènes lors de la lecture des matrices de diagramme de dispersion des composants deU . Certaines caractéristiques - telles que l'apparition de points répartis de manière presque uniforme dans des régions circulaires bien définies - démentent un manque d'indépendance. D'autres, tels que les diagrammes de dispersion montrant des corrélations claires non nulles, dépendent évidemment des choix effectués dans l'algorithme - mais de tels choix ne sont possibles qu'en raison du manque d'indépendance en premier lieu.
Le test ultime d'un algorithme de décomposition comme SVD (ou Cholesky, LR, LU, etc.) est de savoir s'il fait ce qu'il prétend. Dans ce cas, il suffit de vérifier que lorsque SVD renvoie le triple des matrices(U,D,V) , que X est récupéré, jusqu'à l'erreur en virgule flottante anticipée, par le produit UD V′ ; que les colonnes de U et de V sont orthonormées; et que ré est diagonal, ses éléments diagonaux sont non négatifs et sont disposés en ordre décroissant. J'ai appliqué de tels tests à l'
svd
algorithme dansR
et je ne l'ai jamais trouvé erroné. Bien que ce ne soit pas une assurance qu'il est parfaitement correct, une telle expérience - qui, je crois, est partagée par un grand nombre de personnes - suggère que tout bogue nécessiterait une sorte d'entrée extraordinaire pour être manifeste.Ce qui suit est une analyse plus détaillée des points spécifiques soulevés dans la question.
En utilisantk augmente, les corrélations entre les coefficients de U s'affaiblissent, mais elles sont toujours non nulles. Si vous deviez simplement effectuer une simulation plus large, vous constateriez qu'ils sont importants. (Lorsque k = 3 , 50000 itérations devraient suffire.) Contrairement à l'affirmation de la question, les corrélations ne "disparaissent pas complètement".
R
lasvd
procédure de, vous pouvez d'abord vérifier qu'à mesure queDeuxièmement, une meilleure façon d'étudier ce phénomène est de revenir à la question fondamentale de l' indépendance des coefficients. Bien que les corrélations tendent à être proches de zéro dans la plupart des cas, le manque d'indépendance est clairement évident. Ceci est rendu plus évident en étudiant la distribution à plusieurs variables complète des coefficients deU . La nature de la distribution émerge même dans de petites simulations dans lesquelles les corrélations non nulles ne peuvent pas (encore) être détectées. Par exemple, examinez une matrice de nuage de points des coefficients. Pour rendre cela possible, j'ai défini la taille de chaque jeu de données simulé à 4 et j'ai gardé k = 2 , dessinant ainsi 1000 réalisations de la matrice 4 × 2 U , créant une matrice 1000 × 8 . Voici sa matrice de nuage de points complète, avec les variables répertoriées par leur position dans U :
Le fait de parcourir la première colonne révèle un manque d'indépendance intéressant entreu11 et les autres uje j : regardez par exemple comment le quadrant supérieur du nuage de points avec u21 est presque vide; ou examiner le nuage elliptique incliné vers le haut décrivant la relation ( u11, u22) et le nuage incliné vers le bas pour la paire ( u21, u12) . Un examen attentif révèle un manque d'indépendance évident entre presque tous ces coefficients: très peu d'entre eux semblent indépendants à distance, même si la plupart d'entre eux présentent une corrélation proche de zéro.
(NB: la plupart des nuages circulaires sont des projections d'une hypersphère créée par la condition de normalisation forçant la somme des carrés de tous les composants de chaque colonne à l'unité.)
Les matrices de nuages de points aveck = 3 et k = 4 présentent des schémas similaires: ces phénomènes ne se limitent pas à k = 2 , ni ne dépendent de la taille de chaque jeu de données simulé: ils deviennent plus difficiles à générer et à examiner.
Les explications de ces modèles vont à l'algorithme utilisé pour obtenirU dans la décomposition en valeurs singulières, mais nous savons que de tels modèles de non-indépendance doivent exister par les propriétés très définissantes de U : puisque chaque colonne successive est (géométriquement) orthogonale à la précédente celles-ci, ces conditions d'orthogonalité imposent des dépendances fonctionnelles parmi les coefficients, ce qui se traduit par des dépendances statistiques parmi les variables aléatoires correspondantes.
Éditer
En réponse aux commentaires, il peut être intéressant de noter dans quelle mesure ces phénomènes de dépendance reflètent l'algorithme sous-jacent (pour calculer une SVD) et dans quelle mesure ils sont inhérents à la nature du processus.
Les modèles spécifiques de corrélations entre les coefficients dépendent beaucoup des choix arbitraires faits par l'algorithme SVD, car la solution n'est pas unique: les colonnes deU peuvent toujours être indépendamment multipliées par - 1 ou 1 . Il n'y a aucun moyen intrinsèque de choisir le signe. Ainsi, lorsque deux algorithmes SVD font des choix de signe différents (arbitraires ou peut-être même aléatoires), ils peuvent entraîner des motifs différents de diagrammes de dispersion des valeurs ( uje j, uje′j′) . Si vous souhaitez voir cela, remplacez la
stat
fonction dans le code ci-dessous parCe premier réordonne au hasard les observationsui , j et ui , j′ ).
x
, effectue la SVD, puis applique l'ordre inverse àu
pour correspondre à la séquence d'observation d'origine. Comme l'effet est de former des mélanges de versions réfléchies et tournées des diagrammes de dispersion d'origine, les diagrammes de dispersion dans la matrice seront beaucoup plus uniformes. Toutes les corrélations de l'échantillon seront extrêmement proches de zéro (par construction: les corrélations sous-jacentes sont exactement nulles). Néanmoins, le manque d'indépendance restera évident (dans les formes circulaires uniformes qui apparaissent, notamment entreLe manque de données dans certains quadrants de certains des diagrammes de dispersion d'origine (illustrés dans la figure ci-dessus) provient de la façon dont l'
R
algorithme SVD sélectionne les signes pour les colonnes.Rien ne change dans les conclusions. Parce que la deuxième colonne deU est orthogonale à la première, elle (considérée comme une variable aléatoire multivariée) dépend de la première (également considérée comme une variable aléatoire multivariée). Vous ne pouvez pas avoir tous les composants d'une colonne indépendants de tous les composants de l'autre; tout ce que vous pouvez faire est de regarder les données de manière à masquer les dépendances - mais la dépendance persistera.
Voici unk > 2 et dessiner une partie de la matrice de nuage de points.
R
code mis à jour pour gérer les casla source
svd(X,0)
parsvds(X)
dans mon code Matlab fait disparaître l'effet! Pour autant que je sache, ces deux fonctions utilisent des algorithmes SVD différents (les deux sont des routines LAPACK, mais apparemment différentes). Je ne sais pas si R a une fonction similaire à celle de Matlabsvds
, mais je me demande si vous allez toujours maintenir que c'est un "vrai" effet et non un problème numérique.U
vous décidez au hasard si chacune de ses colonnes doit rester telle quelle ou si elle doit changer son signe, les corrélations dont vous parlez ne disparaîtront-elles pas?stat
S
soient dans un ordre particulier; c'est une question de commodité. D'autres routines le garantissent (par exemple MATLABsvds
) mais ce n'est pas une exigence générale. @amoeba: En regardantsvds
(qui semble exempt de ce comportement problématique) le calcul est basé sur le calcul des valeurs propres en premier (donc il n'utilise pas les routines standarddgesdd
/dgesvd
LAPACK - je soupçonne fortement qu'il utilisedsyevr
/dsyevx
au premier abord).Cette réponse présente une réplication des résultats de @ whuber dans Matlab, ainsi qu'une démonstration directe que les corrélations sont un "artefact" de la façon dont l'implémentation SVD choisit de signer pour les composants.
Compte tenu de la longue chaîne de commentaires potentiellement déroutants, je tiens à souligner pour les futurs lecteurs que je suis entièrement d'accord avec ce qui suit:
Ma question était: pourquoi voyons-nous des corrélations élevées de même pour un grand nombre de tirages aléatoires N r e p = 1000 ?∼ 0,2 Nr e p= 1000
Voici une réplication de l'exemple de @ whuber avec , k = 2 et N r e p = 1000 dans Matlab:n = 4 k = 2 Nr e p= 1000
À gauche se trouve la matrice de corrélation, à droite - des diagrammes de dispersion similaires à ceux de @ whuber. L'accord entre nos simulations semble parfait.
Maintenant, suite à une suggestion ingénieuse de @ttnphns, j'attribue des signes aléatoires aux colonnes de , c'est-à-dire après cette ligne:U
J'ajoute les deux lignes suivantes:
Voici le résultat:
Toutes les corrélations disparaissent, exactement comme je m'y attendais depuis le début !
PS. Félicitations à @whuber pour avoir dépassé la réputation de 100k aujourd'hui!
la source
stat
stat <- function(x) { u <- svd(x)$u; as.vector(sign(runif(1) - 1/2)*u %*% diag(sign(diag(u)))) }
svds
svd
R
la source