使用 SVD 方法时 Matlab 和 Eigen 结果之间的差异

Question

我正在尝试将之前用 Matlab 编写的直接线性变换算法重构为 C++。 在比较对矩阵 A 执行 SVD 所得到的矩阵时，我注意到生成的右奇异矩阵 V 的某些元素存在很大的值差异。

确切地说，我感兴趣的是获取与最小奇异值相对应的向量（应该在矩阵 V 的最后一列中找到，就像 Matlab 中的情况一样），但是特征矩阵 V 最右列中的值与Matlab的V有很大不同。我还注意到这些值之间存在符号差异，但我知道这来自 SVD 对于 U 和 V 矩阵的不同符号值有多种解决方案。

这是我正在执行 SVD 的矩阵 A：

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代码：

[U,S, V] = svd(A);

disp('Singular Values S:');
disp(S);


disp('Matrix V:');
disp(V);

disp(U*S*V');

C++代码：

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;

正如您在 C++ 代码中看到的，我尝试设置计算阈值，但这也没有产生任何不同的结果。

以下是 C++ 的结果：

这是 Matlab 中的结果：

我发现有趣的是，重建矩阵A等于输入矩阵A（除了一些零被垃圾值替换），但V矩阵不同。此外，对于方阵，C++ 和 Matlab 的结果都匹配。

我做错了什么，还是因为与奇异值矩阵相乘时的精度误差？

Answer 1

您运行的 MATLAB 版本是什么？无论如何，我在 MATLAB 2023b 和 Octave 8.3、Windows 10 x64 中测试了以下代码。它们之间产生相同的结果，但与你的答案不同。

Z= [ 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];

format bank;
    
A = double(-Z);

disp('------------------------------------------------------------------------------------------------------------');

[U, S, V] = svd(A);

disp(U*S*V');

disp('------------------------------------------------------------------------------------------------------------');

disp(V);