Eigen Jacobi 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()时，结果是一个18×18的矩阵，它不能乘以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; 
}