I'm trying to refactor a Direct Linear Transform algorithm that was previously written in Matlab into C++. While comparing the resulting matrices from performing a SVD on my matrix A, I noticed a big value difference in some elements of the resulting right singular matrix V.
Precisely I'm interested in taking the vector that corresponds to the smallest singular value (which should be found in the last column of matrix V as is the case in Matlab), but the values in rightmost columns of the Eigen matrix V differ greatly from Matlab's V. I also noticed a sign difference between the values but I know that comes from SVD having multiple solutions for different sign values of U and V matrices.
Here is the matrix A on which I'm performing the SVD:
Matrix -A
1 1 1 1 -0 -0 -0 -0 -1 -1 -1 -1
-0 -0 -0 -0 1 1 1 1 -1 -1 -1 -1
2 2 2 1 -0 -0 -0 -0 -4 -4 -4 -2
-0 -0 -0 -0 2 2 2 1 -4 -4 -4 -2
3 3 3 1 -0 -0 -0 -0 -9 -9 -9 -3
-0 -0 -0 -0 3 3 3 1 -9 -9 -9 -3
4 4 4 1 -0 -0 -0 -0 -16 -16 -16 -4
-0 -0 -0 -0 4 4 4 1 -16 -16 -16 -4
5 5 5 1 -0 -0 -0 -0 -25 -25 -25 -5
-0 -0 -0 -0 5 5 5 1 -25 -25 -25 -5
6 6 6 1 -0 -0 -0 -0 -36 -36 -36 -6
-0 -0 -0 -0 6 6 6 1 -36 -36 -36 -6
7 7 7 1 -0 -0 -0 -0 -49 -49 -49 -7
-0 -0 -0 -0 7 7 7 1 -49 -49 -49 -7
8 8 8 1 -0 -0 -0 -0 -64 -64 -64 -8
-0 -0 -0 -0 8 8 8 1 -64 -64 -64 -8
9 9 9 1 -0 -0 -0 -0 -81 -81 -81 -9
-0 -0 -0 -0 9 9 9 1 -81 -81 -81 -9
10 10 10 1 -0 -0 -0 -0 -100 -100 -100 -10
-0 -0 -0 -0 10 10 10 1 -100 -100 -100 -10
Matlab code:
[U,S, V] = svd(A);
disp('Singular Values S:');
disp(S);
disp('Matrix V:');
disp(V);
disp(U*S*V');
C++ code:
Eigen::JacobiSVD<Eigen::MatrixXd> svd;
//svd.setThreshold(1e-10);
svd.compute(A, Eigen::ComputeFullU | Eigen::ComputeFullV);
Eigen::VectorXd singularValues = svd.singularValues();
Eigen::MatrixXd singularValueMatrix = singularValues.asDiagonal();
int dimensionDiff = A.cols() - A.rows();
if (dimensionDiff <= 0) { //Formatting the singular value matrix so that it matches Matlab's
singularValueMatrix.conservativeResize(singularValueMatrix.rows() + -dimensionDiff, singularValueMatrix.cols() );
for (size_t i = 1; i <= -dimensionDiff; i++)
{
for (size_t j = 0; j < singularValueMatrix.cols(); j++)
{
singularValueMatrix(singularValueMatrix.rows()-i,j) = 0;
}
}
}
else {
singularValueMatrix.conservativeResize(singularValueMatrix.rows(), singularValueMatrix.cols() + dimensionDiff);
for (size_t i = 0; i < singularValueMatrix.rows()-1; i++)
{
for (size_t j = 1; j < dimensionDiff; j++)
{
singularValueMatrix(i, singularValueMatrix.cols() - j) = 0;
}
}
}
Eigen::MatrixXd U = svd.matrixU();
Eigen::MatrixXd V = svd.matrixV();
//Displaying the results
std::cout << "Singular Values:\n" << singularValueMatrix << "\n\n";
std::cout << "Left Singular Vectors (U):\n" << U << "\n\n";
std::cout << "Right Singular Vectors (V):\n" << V << "\n\n";
std::cout << "Reconstructed Matrix:\n" << U * singularValueMatrix * V.transpose() << std::endl;
As you can see in the C++ code, I tried setting the calculation threshold but that also didn't yield any different results.
Here are the results in C++: C++ SVD Results
And here are the results in Matlab: Matlab SVD Results
What I find interesting is that the reconstructed matrix A is equal to the input matrix A (except for some zeroes being replaced by garbage values), but the V matrices are different. Also, for square matrices the results from both C++ and Matlab match.
Am I doing something wrong, or is it because of the precision error while multiplying with the singular value matrix?
What version of MATLAB are you running? Anyway, I tested following code in both MATLAB 2023b and Octave 8.3, Windows 10 x64. And they produce close enough but not identical V between themselves.
Octave V:
MATLAB 2023b V: