CODEWITHSUNDEEP
c++eigen3
mkuse
mkuse

Reputation: 2488

Eigen Jacobi SVD

I am trying to compute SVD (Singular Value Decomposition) with Eigen. C is a 27x18 matrix with rank 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 is just a cout. Surprisingly, I see some of the values of diff as very large numbers, something like 4.0733184565807887e+250.

Could I be doing something wrong?

Upvotes: 2

Views: 23030

Answers (1)

Avi Ginsburg
Avi Ginsburg

Reputation: 10596

If you look at the sizes of the matrix elements, you'll notice that svd.matrixU() is 18x18, svd.singularValues() is 18, and svd.matrixV() is 27x27. When you write svd.matrixU() * svd.singularValues().asDiagonal() the result is a 18x18 matrix which cannot multiply svd.matrixV(). You've defined -DNDEBUG which disables bounds checking. The random numbers you saw are what was in memory before allocation. You can get around this using the following code:

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();

As ggael pointed out, you can ask that only the thin matrices be computed, which would look like:

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

Upvotes: 4

Related Questions