Which SVD algorithm implementation is used in MATLAB?

Hi,

I am comparing singular value decomposition function [U,S,V] = svd(A) to some C implementations of the algorithm. However, I am getting somewhat different results: for example, columns of the output matrix U are mixed-up, or some output values have different sign, etc.

Does anyone know which SVD algorithm implementation is used in MATLAB?

[Walter Roberson]

Did you compare against MKL and against LAPACK ?

https://www.intel.com/content/www/us/en/developer/articles/technical/solve-svd-problem-for-sparse-matrix-with-intel-math-kernel-library.html

http://www.netlib.org/lapack/explore-html/d1/d7e/group__double_g_esing_ga84fdf22a62b12ff364621e4713ce02f2.html

[Danijel Domazet]

I didn't compare to MKD nor LAPACK. MKD is copyrighted, and LAPACK I can't use since it written in Fortran and I need a plain C version.

[Danijel Domazet]

Walter, I have tried LAPACK, but no success. Please check here:

https://www.mathworks.com/matlabcentral/answers/1655690-matlab-and-lapack-implementation-of-the-svd-algorithm-not-the-same-output.

[Walter Roberson]

I am not sure how it is relevant that MKL is copyrighted? LAPACK is copyrighted too, with a modified BSD license.

MKL is probably closed source, but you did not indicate that you needed source access; for everything you have mentioned so far you only need to be able to link against it.

Answer 1

[Steven Lord]

Keep in mind that the SVD of a matrix is not unique. Quoting from Wikipedia: "Non-degenerate singular values always have unique left- and right-singular vectors, up to multiplication by a unit-phase factor (for the real case up to a sign). Consequently, if all singular values of a square matrix M are non-degenerate and non-zero, then its singular value decomposition is unique, up to multiplication of a column of U by a unit-phase factor and simultaneous multiplication of the corresponding column of V by the same unit-phase factor."

A = magic(5);
[U, S, V] = svd(A);
format longg
check1 = U*S*V' - A

check1 = 5×5

1.0e+00 * 3.5527136788005e-15 -1.4210854715202e-14 -1.93178806284777e-14 2.48689957516035e-14 1.59872115546023e-14 0 -1.77635683940025e-14 -1.4210854715202e-14 2.48689957516035e-14 1.77635683940025e-14 -3.5527136788005e-15 6.21724893790088e-15 1.77635683940025e-14 7.105427357601e-15 7.105427357601e-15 -3.5527136788005e-15 3.19744231092045e-14 1.77635683940025e-14 -1.77635683940025e-14 -8.43769498715119e-15 5.32907051820075e-15 3.5527136788005e-14 1.4210854715202e-14 -2.17603712826531e-14 -1.95399252334028e-14

U2 = -U;
V2 = -V;
check2 = U2*S*V2' - A

check2 = 5×5

1.0e+00 * 3.5527136788005e-15 -1.4210854715202e-14 -1.93178806284777e-14 2.48689957516035e-14 1.59872115546023e-14 0 -1.77635683940025e-14 -1.4210854715202e-14 2.48689957516035e-14 1.77635683940025e-14 -3.5527136788005e-15 6.21724893790088e-15 1.77635683940025e-14 7.105427357601e-15 7.105427357601e-15 -3.5527136788005e-15 3.19744231092045e-14 1.77635683940025e-14 -1.77635683940025e-14 -8.43769498715119e-15 5.32907051820075e-15 3.5527136788005e-14 1.4210854715202e-14 -2.17603712826531e-14 -1.95399252334028e-14

The elements of check1 and check2 both look sufficiently small as to be close enough to 0 to say both (U, S, V) and (U2, S, V2) satisfy the SVD contract. It's possible the columns of the output U aren't "mixed up" they're just a different (correct) answer (when combined with the corresponding V) than you expected.

[Danijel Domazet]

Thanks. Let me investigate some more.

[Danijel Domazet]

Tried to use LAPACK library, please check:

https://www.mathworks.com/matlabcentral/answers/1655690-matlab-and-lapack-implementation-of-the-svd-algorithm-not-the-same-output