艾根·雅各比 SVD

Eigen Jacobi SVD

本文关键字：SVD 艾根      更新时间：2023-10-16

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

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 .

我会做错什么吗？

如果您查看矩阵元素的大小，您会注意到svd.matrixU()是 18x18，svd.singularValues()是 18，svd.matrixV()是 27x27。当你写svd.matrixU() * svd.singularValues().asDiagonal()结果是一个 18x18 矩阵，它不能乘以svd.matrixV()。您已经定义了禁用边界检查的 -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;
}
相关文章：