1

I'm currently trying to implement some prototype Python code into C++ and am running into an issue where I am getting different results for svd calculate between the 2 when using the same exact array

Code:

python:

                        (U, S, V) = np.linalg.svd(H)
                        print("h\n",H)
                        print("u\n",U)
                        print("v\n",V)
                        print("s\n",S)
                        rotation_matrix = np.dot(U, V)

prints

h
 [[ 1.19586781e+00 -1.36504900e+00  3.04707238e+00]
 [-3.24276981e-01  4.25640964e-01 -6.78455372e-02] 
 [ 4.58970250e-02 -7.33566042e-02 -2.96605698e-03]]
u
 [[-0.99546325 -0.09501679  0.0049729 ]
 [ 0.09441242 -0.97994807  0.17546529]
 [-0.01179897  0.17513875  0.98447306]]
v
 [[-0.34290622  0.39295764 -0.85322893]
 [ 0.49311955 -0.6977843  -0.51954806]
 [-0.79953014 -0.59890012  0.04549948]]
s
 [3.5624894  0.43029207 0.00721429]

C++ code:

    std::cout << "H\n" << HTest << std::endl;
    Eigen::JacobiSVD<Eigen::MatrixXd> svd;

    svd.compute(HTest, Eigen::ComputeThinV | Eigen::ComputeThinU);
    std::cout << "h is" << std::endl << HTest << std::endl;
    std::cout << "Its singular values are:" << std::endl << svd.singularValues() << std::endl;
    std::cout << "Its left singular vectors are the columns of the thin U matrix:" << std::endl << -1*svd.matrixU() << std::endl;
    std::cout << "Its right singular vectors are the columns of the thin V matrix:" << std::endl << -1*svd.matrixV() << std::endl;

prints:

h is
    1.19587    -1.36505     3.04707
  -0.324277    0.425641  -0.0678455
   0.045897  -0.0733566 -0.00296606
Its singular values are:
   3.56249
  0.430292
0.00721429
Its left singular vectors are the columns of the thin U matrix:
 -0.995463 -0.0950168  0.0049729
 0.0944124  -0.979948   0.175465
 -0.011799   0.175139   0.984473
Its right singular vectors are the columns of the thin V matrix:
-0.342906   0.49312  -0.79953
 0.392958 -0.697784   -0.5989
-0.853229 -0.519548 0.0454995

So H,U,S are eqivalent between the 2, but V is not. What could cause this?

Jesper Juhl
  • 30,449
  • 3
  • 47
  • 70
cheesemas46
  • 139
  • 12
  • https://stackoverflow.com/questions/588004/is-floating-point-math-broken – Jesper Juhl Nov 04 '22 at 19:32
  • @JesperJuhl Do you think the difference between .39 and .49 is from floating point math? Theres differences in the V matrixx of .1 +, I dont think this a floating point issue........ – cheesemas46 Nov 04 '22 at 19:34
  • 4
    Just in case, you haven't notice, one V is just the transpose of the other. Is your question is "why is it so?". Or, was it that you did not notice that it is the same matrix, transpose? – chrslg Nov 04 '22 at 19:47
  • @chrslg ohhh I see, I missed that. Why is that so tho? – cheesemas46 Nov 04 '22 at 19:51
  • 2
    I'd say, "because" :-). I don't think there is a good reason. Just 2 implementations. In maths lessons, you've probably learned SVD decomposition with formula M=U.S.Vᵀ. So C++ library probably stick to this formula, and gives U, S, V such as M=U.S.Vᵀ. Where as linalg documentation says that it returns U,S,V such as `M=(U*S)@V`. So one call V what the other calls Vᵀ. Hard to say which one is right. As long as they do what their doc say they do... – chrslg Nov 04 '22 at 19:57

1 Answers1

1

I didn't notice that the V's were just transposes of each other. user chrslg has a good explanation for why this is so Ill just copy it here:

"I'd say, "because" :-). I don't think there is a good reason. Just 2 implementations. In maths lessons, you've probably learned SVD decomposition with formula M=U.S.Vᵀ. So C++ library probably stick to this formula, and gives U, S, V such as M=U.S.Vᵀ. Where as linalg documentation says that it returns U,S,V such as M=(U*S)@V. So one call V what the other calls Vᵀ. Hard to say which one is right. As long as they do what their doc say they do"

cheesemas46
  • 139
  • 12
  • In regards to the why: Eigen uses column-major order (Fortran order). Numpy commonly uses row-major order (C-order). So accessing individual singular vectors in numpy is more efficient if they are in one row, and in Eigen when they are in a column – Homer512 Nov 06 '22 at 19:57