Why does SVD result of Armadillo differ from NumPy?

Question

In my Python code, I was computing SVD of some data using numpy.linalg.svd:

from numpy import linalg
(_, _, v) = linalg.svd(m)

V matrix returned by this was:

[[ 0.4512937  -0.81992002 -0.35222884]
 [-0.22254721  0.27882908 -0.93419863]
 [ 0.86417981  0.4999855  -0.05663711]]

While porting my code to C++, I switched to using Armadillo for computing SVD:

#include <armadillo>

arma::fmat M; // Input data
arma::fmat U;
arma::fvec S;
arma::fmat V;
arma::svd(U, S, V, M);

The resulting V for the same data is:

  0.4513  -0.2225  -0.8642
 -0.8199   0.2788  -0.5000
 -0.3522  -0.9342   0.0566

We can see that the transpose of V from Armadillo matches V from NumPy. Except that is, for the last column of V from Armadillo. Those values have the opposite sign of the values in last row of NumPy result.

What is happening here? Why do the SVD results from two popular libraries differ like this? And which of the two is the correct result?

Answer 1

Both are correct... The rows of the v you got from numpy are the eigenvectors of M.dot(M.T) (the transpose would be a conjugate transpose in the complex case). Eigenvectors are in the general case defined only up to a multiplicative constant, so you could multiply any row of v by a different number, and it will still be an eigenvector matrix.

There is the additional constraint on v that it be a unitary matrix, which loosely translates to its rows being orthonormal. This reduces your available choices for every eigenvector to only 2: the normalized eigenvector pointing in either direction. But you still get to multiply any row by -1 and still have a valid v.

If you want to test it for your matrix, which I have loaded as a:

>>> u, d, v = np.linalg.svd(a)
>>> D = np.zeros_like(a)
>>> idx = np.arange(a.shape[1])
>>> D[idx, idx] = d
>>> np.allclose(a, u.dot(D).dot(v))
True
>>> v[2] *= -1
>>> np.allclose(a, u.dot(D).dot(v))
True

---

Actually, you can only multiply the rows of v by -1 in the real domain, but in the complex case you can multiply them by any complex number of absolute value 1:

>>> vv = v.astype(np.complex)
>>> vv[0] *= (1+1.j)/np.sqrt(2)
>>> np.allclose(a, u.dot(D).dot(v))
True