c++ - Eigen 雅可比 SVD

Question

我正在尝试使用 Eigen 计算 SVD(奇异值分解)。 C 是一个 27x18 矩阵，秩为 15。

JacobiSVD<MatrixXd> svd( C, ComputeFullV | ComputeFullU );
cout << svd.computeU() << endl;
cout << svd.computeV() << endl;

MatrixXd Cp = svd.matrixU() * svd.singularValues().asDiagonal() * svd.matrixV().transpose();
MatrixXd diff = Cp - C;
PRINT_MAT( "diff", diff );

PRINT_MAT 只是一个cout。令人惊讶的是，我看到 diff 的某些值是非常大的数字，例如 4.0733184565807887e+250。

我会不会做错了什么？

Answer 1

如果您查看矩阵元素的大小，您会注意到 svd.matrixU() 是 18x18，svd.singularValues() 是 18，并且svd.matrixV() 是 27x27。当您编写 svd.matrixU() * svd.singularValues().asDiagonal() 时，结果是一个无法乘以 svd.matrixV() 的 18x18 矩阵。您已经定义了禁用边界检查的 -DNDEBUG。您看到的随机数是分配之前内存中的内容。您可以使用以下代码解决此问题:

MatrixXd res(C.rows(), C.cols());
res.setZero();
res.topLeftCorner(C.rows(), C.rows()) = (svd.matrixU() * svd.singularValues().asDiagonal());

MatrixXd Cp = res * svd.matrixV().transpose();
MatrixXd diff = Cp - C;
cout << "diff:\n" << diff.array().abs().sum();

正如 ggael 指出的那样，您可以要求只计算薄矩阵，这看起来像:

#include <Eigen/Core>
#include <Eigen/SVD>
#include <iostream>

using namespace Eigen;
using std::cout;

int main()
{
    MatrixXd C;
    C.setRandom(27,18);
    JacobiSVD<MatrixXd> svd( C, ComputeThinU | ComputeThinV);
    MatrixXd Cp = svd.matrixU() * svd.singularValues().asDiagonal() * svd.matrixV().transpose();
    MatrixXd diff = Cp - C;
    cout << "diff:\n" << diff.array().abs().sum() << "\n";
    return 0;
}