Re: [eigen] SVD with singular matrices

[Benoit Jacob <jacob.benoit.1@xxxxxxxxx>]

OK, i've attached a compilable variant that also prints the
reconstructed matrix.

It turns out that the reconstructed matrix is
  1 0
  1 0

I.e. we seem to have a problem with row/column permutations. Which
could be related to the sorting of singvals.

I'm trying to understand why our unit test is not catching this. It's
true that we're not testing matrices with exactly 0 singvals.

Benoit

2010/9/29  <hamelin.philippe@xxxxxxx>:
> This simple test fails:
>
> bool testSingularMatrixSVD()
> {
>    MatrixXd a(2,2), avalidation(2,2);
>    unsigned int rows = 2;
>    unsigned int cols = 2;
>
>    a <<    1,     1,
>            0,     0;
>
>    Eigen::SVD<MatrixXd> svd(a);
>    MatrixXd sigma = MatrixXd::Zero(rows,cols);
>    MatrixXd matU  = MatrixXd::Zero(rows,rows);
>    sigma.diagonal() = svd.singularValues();
>    matU = svd.matrixU();
>    avalidation = matU * sigma * svd.matrixV().transpose();
>
>    cout << "testSingularMatrixSVD() : ";
>    if((a).isApprox(avalidation, 1e-12))
>    {
>        cout << "Success." << endl;
>        return true;
>    }
>    else
>    {
>        cout << "Fail." << endl;
>        return false;
>    }
>
> }
>
>
> -----Message d'origine-----
> De : Listengine [mailto:listengine@xxxxxxxxxxxxxxxxx] De la part de Benoit Jacob
> Envoyé : 29 septembre 2010 10:32
> À : eigen@xxxxxxxxxxxxxxxxxxx
> Objet : Re: [eigen] SVD with singular matrices
>
> 2010/9/29  <hamelin.philippe@xxxxxxx>:
>> Hello,
>>
>> when updating from eigen2 to eigen3, I found that my pseudo-inverse was not
>> working correctly for singular matrices. However, just replacing SVD with
>> JacobiSVD makes it work. Before looking further, is there any limitation
>> with SVD with singular matrices such as:
>>
>> [1 1;0 0]
>
> JacobiSVD always works, and is always precise. But it's slow for large matrices.
>
> SVD uses the Golub-Kahan bidiagonalization approach. So it's faster,
> but potentially inaccurate for certain singular matrices.
>
> I would be surprised though if the above 2x2 matrix was enough to
> expose problems in it. That would be a bug. Can you make a compilable
> test case?
>
> Thanks,
> Benoit
>
>>
>> Thank you,
>>
>> ------------------------------------
>> Philippe Hamelin, ing. jr, M. Ing
>> Chercheur / Researcher
>>
>> T: 450-652-8499 x2198
>> F: 450-652-1316
>>
>> Expertise robotique et civil
>> Institut de recherche d'Hydro-Québec (IREQ)
>> 1740, boul. Lionel-Boulet
>> Varennes (QC) J3X 1S1, Canada
>>
>
>
>
>
>

#include<Eigen/Eigen>
#include<iostream>
using namespace Eigen;
using namespace std;

int main()
{
   MatrixXd a(2,2), avalidation(2,2);
   unsigned int rows = 2;
   unsigned int cols = 2;

   a <<    1,     1,
           0,     0;

   Eigen::SVD<MatrixXd> svd(a);
   MatrixXd sigma = MatrixXd::Zero(rows,cols);
   MatrixXd matU  = MatrixXd::Zero(rows,rows);
   sigma.diagonal() = svd.singularValues();
   matU = svd.matrixU();
   avalidation = matU * sigma * svd.matrixV().transpose();

   cout << "relative error: " << (a-avalidation).norm() / a.norm() << endl;
   cout << "reconstructed matrix:\n" << avalidation << endl;
}