Eigenvectors computed with numpy's eigh and svd do not match

Question

Consider singular value decomposition M=USV*. Then the eigenvalue decomposition of M* M gives M* M= V (S* S) V*=VS* U* USV*. I wish to verify this equality with numpy by showing that the eigenvectors returned by eigh function are the same as those returned by svd function:

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))

The last assertion fails. Why?

Calendula

You can add a positive number to all the diagonal elements, i.e. make M2=M+a*I, where a is large enough to make M2 positive semidefnite. Then SVD and eigh should agree better. – Heading

Answer 1

Just play with small numbers to debug your problem.

Start with A=np.random.randn(3,2) instead of your much larger matrix with size (50,20)

In my random case, I find that

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

and for v2:

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

they only differ for a sign, and obviously, even if normalized to have unit module, the vector can differ for a sign.

Now if you try the trick

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

for your original big matrix, it fails... again, this is OK. What happens is that some vectors have been multiplied by -1, some others haven't.

A correct way to check for equality between the vectors is:

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

and indeed, to get a feeling of how this works you can print this quantity:

(V1*V2).sum(0)

that indeed is either +1 or -1 depending on the vector:

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

EDIT: This will happen in most cases, especially if starting from a random matrix. Notice however that this test will likely fail if one or more eigenvalues has an eigenspace of dimension larger than 1, as pointed out by @Sven Marnach in his comment below:

There might be other differences than just vectors multiplied by -1. If any of the eigenvalues has a multi-dimensional eigenspace, you might get an arbitrary orthonormal basis of that eigenspace, and to such bases might be rotated against each other by an arbitraty unitarian matrix

Galliard

@Calendula OK, I'm lost :) But I trust in your judgment, and so I'll remove my comments not to confuse future readers. Cheers! – Inhuman

There might be other differences than just vectors multiplied by -1. If any of the eigenvalues has a multi-dimensional eigenspace, you might get an arbitrary orthonormal basis of that eigenspace, and to such bases might be rotated against each other by an arbitraty unitarian matrix. – Tisatisane

@SvenMarnach, this is a very valid point. I will edit the post to point out this caveat – Galliard