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; 
}