Pourquoi ne puis-je pas obtenir un SVD valide de X via la décomposition de valeurs propres de XX 'et X'X?

9

J'essaie de faire SVD à la main:

m<-matrix(c(1,0,1,2,1,1,1,0,0),byrow=TRUE,nrow=3)

U=eigen(m%*%t(m))$vector
V=eigen(t(m)%*%m)$vector
D=sqrt(diag(eigen(m%*%t(m))$values))

U1=svd(m)$u
V1=svd(m)$v
D1=diag(svd(m)$d)

U1%*%D1%*%t(V1)
U%*%D%*%t(V)

Mais la dernière ligne ne revient pas m. Pourquoi? Cela semble avoir quelque chose à voir avec les signes de ces vecteurs propres ... Ou ai-je mal compris la procédure?

statisticien raté
la source
On me dit à plusieurs reprises que le signe n'a pas d'importance dans les SVD ... comme ça
statisticien raté
@Amoeba Merci d'avoir clarifié cela. Je me concentrais sur la question anglaise plutôt que sur le code. Statisticien défaillant: voyez ce qui D=diag(c(-1,1,1)*sqrt(eigen(m%*%t(m))$values))fait et gardez à l'esprit que la racine carrée (ainsi que tout vecteur propre normalisé) n'est définie que jusqu'à signer. Pour plus d'informations, passez mà m <- matrix(-2,1,1)et incluez ,1,1)à la fin de chacun des appels à diag. Ceci est un exemple qui crée le même problème - mais c'est tellement simple que la nature du problème deviendra complètement évidente. 1×1
whuber
1
Je l'ai. Je vous remercie! Avez-vous une règle générale pour déterminer le vecteur c (-1, 1, 1)? Ou comment relier les signes des deux décompositions?
failedstatistician
1
Notez que le truc de @ whuber avec c(-1,1,1)fonctionne, mais Ddéfini comme cela ne vous donne pas de valeurs singulières. Les valeurs singulières doivent toutes être positives par définition. La question de savoir comment relier les signes de Uet Vest bonne, et je n'ai pas de réponse. Pourquoi ne faites-vous pas simplement un SVD? :-)
amibe

Réponses:

13

Analyse du problème

La SVD d'une matrice n'est jamais unique. Soit la matrice dimensions et que sa SVD soitn × kAn×k

A=UDV

pour un la matrice avec des colonnes orthonormales, une diagonale la matrice avec des entrées non-négatives, et un la matrice avec des colonnes orthonormales.U p × p D k × p Vn×pUp×pDk×pV

Maintenant, choisissez arbitrairement n'importe quelle diagonale matrice ayant s sur la diagonale, de sorte que soit l' identité . alorsS ± 1 S 2 = I p × p I pp×pS±1S2=Ip×pIp

A=UDV=UIDIV=U(S2)D(S2)V=(US)(SDS)(VS)

est également un SVD de parce que montre que a des colonnes orthonormées et un calcul similaire montre que a des colonnes orthonormées. De plus, puisque et sont diagonaux, ils font la navette, d'où montre que toujours des entrées non négatives.( U S ) ( U S ) = S U U S = S I p S = S S = S 2 = I p U S V S S D S D S = D S 2 = D DA

(US)(US)=SUUS=SIpS=SS=S2=Ip
USVSSD
SDS=DS2=D
D

La méthode implémentée dans le code pour trouver un SVD trouve un qui diagonalise et, de la même façon, un qui diagonale Il procède au calcul de en fonction des valeurs propres trouvées dans . Le problème est que cela ne garantit pas une correspondance cohérente des colonnes de avec les colonnes de .A A = ( U D V ) ( U D V ) = U D V V D U = U D 2 U V A A = V D 2 V . D D 2 U VU

AA=(UDV)(UDV)=UDVVDU=UD2U
V
AA=VD2V.
DD2UV

Une solution

Au lieu de cela, après avoir trouvé un tel et un tel , utilisez-les pour calculerVUV

UAV=U(UDV)V=(UU)D(VV)=D

directement et efficacement. Les valeurs diagonales de ce ne sont pas nécessairement positives. D (C'est parce qu'il n'y a rien dans le processus de diagonalisation de ou qui garantira que, puisque ces deux processus ont été effectués séparément.) Rendez-les positifs en choisissant les entrées le long de la diagonale de pour égaler les signes des entrées de , de sorte que ait toutes les valeurs positives. Compensez cela en multipliant à droite par :AAAASDSDUS

A=UDV=(US)(SD)V.

C'est un SVD.

Exemple

Soit avec . Un SVD estn=p=k=1A=(2)

(2)=(1)(2)(1)

avec , et .U=(1)D=(2)V=(1)

Si vous diagonalisezAA=(4)U=(1)D=(4)=(2)AA=(4)V=(1)

UDV=(1)(2)(1)=(2)A.
D=UAV=(1)(2)(1)=(2).
S=(1)UUS=(1)(1)=(1)DSD=(1)(2)=(2)
A=(1)(2)(1),

Code

Voici le code modifié. Sa sortie confirme

  1. La méthode recrée mcorrectement.
  2. U et sont vraiment encore orthonormés.V
  3. Mais le résultat n'est pas le même SVD renvoyé par svd. (Les deux sont également valables.)
m <- matrix(c(1,0,1,2,1,1,1,0,0),byrow=TRUE,nrow=3)

U <- eigen(tcrossprod(m))$vector
V <- eigen(crossprod(m))$vector
D <- diag(zapsmall(diag(t(U) %*% m %*% V)))
s <- diag(sign(diag(D)))  # Find the signs of the eigenvalues
U <- U %*% s              # Adjust the columns of U
D <- s %*% D              # Fix up D.  (D <- abs(D) would be more efficient.)

U1=svd(m)$u
V1=svd(m)$v
D1=diag(svd(m)$d,n,n)

zapsmall(U1 %*% D1 %*% t(V1)) # SVD
zapsmall(U %*% D %*% t(V))    # Hand-rolled SVD
zapsmall(crossprod(U))        # Check that U is orthonormal
zapsmall(tcrossprod(V))       # Check that V' is orthonormal
whuber
la source
1
+1. C’est très clair. J'ajouterais seulement qu'en pratique, il suffit de calculer l'un Uou l' autre Vet puis d'obtenir une autre matrice en multipliant avec A. De cette façon, on n'effectue qu'une (au lieu de deux) compositions de composition, et les signes sortiront bien.
amoeba
2
@Amoeba C'est vrai: dans l'esprit du calcul manuel d'un SVD, qui est évidemment un exercice éducatif, aucune attention n'est accordée ici à l'efficacité.
whuber
2
Merci pour ton aide! Je pense que je comprends ce problème (enfin).
failedstatistician
3
@Federico Merci pour ce rappel. Vous avez tout à fait raison - j'ai implicitement supposé que toutes les valeurs propres sont distinctes, car en effet, cela va sûrement être le cas dans les applications statistiques et on sort de l'habitude de considérer les ambiguïtés avec les espaces propres "dégénérés".
whuber
3
Vous avez raison, ce n'est qu'un cas de bord, et un cas délicat en effet. Dans un certain sens, il est une autre manifestation du même problème que vous décrivez dans votre réponse, que cette méthode ne garantit pas un « correspondant » entre les colonnes de et . Le calcul de la SVD à partir des compositions des eigendes est toujours un excellent exemple d'apprentissage. UV
Federico Poloni
5

Comme je l'ai souligné dans un commentaire à la réponse de @ whuber, cette méthode pour calculer le SVD ne fonctionne pas pour chaque matrice . Le problème ne se limite pas aux panneaux.

Le problème est qu'il peut y avoir des valeurs propres répétées, et dans ce cas, la composition propre de et n'est pas unique et tous les choix de et ne peuvent pas être utilisés pour récupérer le facteur diagonal de la SVD. Par exemple, si vous prenez une matrice orthogonale non diagonale (disons ), alors . Parmi tous les choix possibles pour la matrice de vecteur propre de , retournera , donc dans ce cas n'est pas diagonal.A A ' U V A = [ 3 / 5 4 / 5 - 4 / 5 3 / 5 ] A A ' = A ' A = I I U = V = I U ' A V = AAAAAUVA=[3/54/54/53/5]AA=AA=IIeigenU=V=IUAV=A

Intuitivement, ceci est une autre manifestation du même problème que @whuber souligne, qu'il doit y avoir une "correspondance" entre les colonnes de et , et le calcul de deux compositions à eigend séparément ne le garantit pas.VUV

Si toutes les valeurs singulières de sont distinctes, alors la composition par eigend est unique (jusqu'à l'échelle / signes) et la méthode fonctionne. Remarque: ce n'est toujours pas une bonne idée de l'utiliser en code de production sur un ordinateur avec une arithmétique à virgule flottante, car lorsque vous les produits et le résultat calculé peut être perturbé par une quantité de l'ordre de , où est la précision de la machine. Si les grandeurs des valeurs singulières diffèrent considérablement (de plus de , grosso modo), cela nuit à la précision numérique des plus petites.A A A A A 2 uAAAAAA2uu2×1016108

Le calcul de la SVD à partir des deux compositions de type eigend est un excellent exemple d'apprentissage, mais dans les applications réelles, utilisez toujours la svdfonction de R pour calculer la décomposition en valeurs singulières.

Federico Poloni
la source
1
Ce commentaire est un bon conseil. Veuillez noter, cependant, que ce fil ne se préoccupe pas de la bonne façon de calculer la SVD (et je pense que personne ne contesterait votre recommandation). L'OP accepte implicitement que cela svdfonctionne. En effet, ils l'utilisent comme un standard par rapport auquel comparer un calcul manuel, dont le but est de vérifier la compréhension, de ne le remplacer svden aucune façon.
whuber
@whuber Observation correcte; J'ai reformulé le dernier paragraphe.
Federico Poloni