Gli autovettori calcolati con numpy eigh e svd non corrispondono a

Question

Considera la decomposizione del valore singolare M = USV *. Quindi la scomposizione dell'autovalore di M * M fornisce M * M = V (S * S) V * = VS * U * USV *. Vorrei verificare questa uguaglianza con NumPy mostrando che gli autovettori restituiti dalla eigh funzione sono gli stessi di quelli restituiti da svd funzione:

import numpy as np 
np.random.seed(42) 
# create mean centered data 
A=np.random.randn(50,20) 
M= A-np.array(A.mean(0),ndmin=2) 

# svd 
U1,S1,V1=np.linalg.svd(M) 
S1=np.square(S1) 
V1=V1.T 

# eig 
S2,V2=np.linalg.eigh(np.dot(M.T,M)) 
indx=np.argsort(S2)[::-1] 
S2=S2[indx] 
V2=V2[:,indx] 

# both Vs are in orthonormal form 
assert np.all(np.isclose(np.linalg.norm(V1,axis=1), np.ones(V1.shape[0]))) 
assert np.all(np.isclose(np.linalg.norm(V1,axis=0), np.ones(V1.shape[1]))) 
assert np.all(np.isclose(np.linalg.norm(V2,axis=1), np.ones(V2.shape[0]))) 
assert np.all(np.isclose(np.linalg.norm(V2,axis=0), np.ones(V2.shape[1]))) 

assert np.all(np.isclose(S1,S2)) 
assert np.all(np.isclose(V1,V2)) 

L'ultima asserzione fallisce. Perché?

Answer 1

Basta giocare con numeri piccoli per eseguire il debug del problema.

Inizia con A=np.random.randn(3,2) al posto del tuo molto più grande matrice con la dimensione (50,20)

Nel mio caso a caso, trovo che

v1 = array([[-0.33872745, 0.94088454], 
    [-0.94088454, -0.33872745]]) 

e per v2:

v2 = array([[ 0.33872745, -0.94088454], 
    [ 0.94088454, 0.33872745]]) 

si differenziano solo per un segno, e ovviamente, anche se normalizzato per avere un modulo unità, il vettore può differire per un segno.

Ora, se si cerca il trucco

assert np.all(np.isclose(V1,-1*V2)) 

per il vostro grande matrice originaria, non riesce ... ancora una volta, questo è OK. Quello che succede è che alcuni vettori sono stati moltiplicati per -1, altri no.

Un modo corretto per verificare la presenza di parità tra i vettori è:

assert allclose(abs((V1*V2).sum(0)),1.) 

e anzi, per avere un'idea di come questo si lavora in grado di stampare questa quantità:

(V1*V2).sum(0) 

che in effetti è sia +1 o -1 seconda del vettore:

array([ 1., -1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 
    1., -1., 1., 1., 1., -1., -1.]) 

MODIFICA: Questo accadrà nella maggior parte dei casi, specialmente se si parte da una matrice casuale. Si noti, tuttavia, che questo test sarà probabilmente esito negativo se uno o più autovalori ha un'autospazio di dimensione più grande di 1, come sottolineato da @Sven Marnach nel suo commento qui sotto:

Ci potrebbero essere altre differenze non solo vettori moltiplicato per -1. Se uno qualsiasi degli autovalori ha un autospazio multi-dimensionale, si potrebbe ottenere una base ortonormale arbitraria di quella autospazio, e per tali basi potrebbero essere ruotati l'uno contro l'altro da una matrice unitaria arbitraty