Hi

In a recent thread there was an error in how a matrix is reconstructed from its SVD decomposition. I apologize if this is just an old and settled issue and I am just adding noise, but I got bitten by numpy's unfamiliar output myself a long time ago and I see others get confused as well. So what is the issue?

Let the svd decomposition of X be USV:
U, S, V = np.linalg.svd(X),
Then U has X's left singular vectors in its columns, S contains the singular values, and V has X's right singular vectors in its *rows*. The reconstruction of X is (matrix notation) will be: U*diag(S)*V  (not U*diag(S)*V.T).
All other high-level software (Matlab, Scilab, Mathematica, R, etc), outputs the right singular vectors columnwise, that is V = [v1, v2, v3, ... vn], where vn is a column (eigenvector) in V, thus the reconstruction would be U*diag(S)*V.T. Also, as far as I know, most linear algebra textbooks operate with eigenvectors consistent as column vectors in explanations of the SVD. I think numpy's svd should do so too.

I know lapack's dgesdd returns V.T (or conjugate), and specify that in its documentation so this is a true interface of the library, but I still think its wrong and its just too confusing for any beginner who usually has experience in other software, such as matlab, prior to numpy. Also, ,as I am typing here, I realize that changing the output would break lots of stuff, and pass silently through many tests as the shape of V is similar (if full_matrices=0). Oh well, I guess that proposal is off the table? Perhaps some *stronger* hints in the documentation are needed.

Arnar

PS:
In the docs at http://www.scipy.org/NumPy_for_Matlab_Users , the svd equivalents have wrong notation, this is not helping :-). I didnt manage to change it, perhaps some other may be so kind?