Re: [eigen] SVD with singular matrices
• To: eigen@xxxxxxxxxxxxxxxxxxx
• Subject: Re: [eigen] SVD with singular matrices
• From: Keir Mierle <mierle@xxxxxxxxx>
• Date: Wed, 29 Sep 2010 17:50:48 -0700
• Domainkey-signature: a=rsa-sha1; c=nofws; d=gmail.com; s=gamma; h=mime-version:in-reply-to:references:date:message-id:subject:from:to :content-type; b=F9zvRZ7B3+i48Ps3yZRvTtfSwd5f8H38uQKiOTLGc/yCbwD8WXrO6YYTALOIaFROM6 cJ7uf2RNAA+Dxbl042q9dqjMQwh1Nu7MJfxVuAeO8ugfmeqJ6MzvxHJs1dPk6xBXXAtS 6QMb6ximlW7Oqj2UGuBTSsArcFHa4R/sE45kE=

On Wed, Sep 29, 2010 at 8:10 AM, Benoit Jacob wrote:
Here's a diff for the SVD unit test, exposing the problem.

The sorting of eigenvalues isn't the problem, the bidiagonalization is.

The most productive thing I can say is: let's declare SVD
not-for-exactly-singular-matrices for now, and ... yeah , yeah, make
the new SVD happen :-P .. in the meanwhile you have JacobiSVD.

This is a serious problem. Almost all multiview geometry hackers regularly solve matrices that are exactly singular (i.e. start out as 3x5, get extended with two rows of zeros to be 5x5). This is because many problems in multiview geometry produce homogenous linear equations. I use the SVD now in libmv and it seems to work, but now I'm worried.

Keir

Benoit

2010/9/29 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
>>>
>>
>>
>>
>>
>>
>

 Mail converted by MHonArc 2.6.19+ http://listengine.tuxfamily.org/