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](http://docs.scipy.org/doc/numpy/reference/generated/numpy.linalg.svd.html):

    
    
    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](http://arma.sourceforge.net/) 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: 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](http://en.wikipedia.org/wiki/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
    

