commons-dev mailing list archives

Site index · List index
Message view « Date » · « Thread »
Top « Date » · « Thread »
From Phil Steitz <>
Subject Re: [math] Work going on on SVD
Date Wed, 03 Feb 2010 11:33:46 GMT
Luc Maisonobe wrote:
> Luc Maisonobe a écrit :
>> Dimitri Pourbaix a écrit :
>>> Hi,
>> Hi Dimitri,
>>> In its present form (CM2.0), SingularValueDecomposition suffers some
>>> problems when the matrix is (numerically) singular.  Luc proposed a
>>> way to improve the situation by limiting the singular values to non
>>> zero ones.  Whereas the result is OK if someone is interested in getting
>>> the singular values, it is not from a purely mathematical point of view.
>>> Indeed, the resulting U matrix does no longer hold the right dimension.
>>> Instead of a number of columns equal to the number of columns of the
>>> original matrix, U now has as many columns as non-zero singular values.
>>> The product U*S*V^T yields the original matrix but the wrong size might
>>> put some users into trouble (that is true for the size of S as well).
>> I agree with this analysis. The current behavior is interesting for some
>> use cases but not for all of them. So I suggest we provide both cases.
>> This could be done either by renaming the current implementation as
>> TruncatedSVD and having your implementation named
>> SingularValueDecomposition, or by using some constructor parameters to
>> select the desired behavior.
> What do we do here ? If Dimitri can provide a new mathematically
> compliant SVD I would suggest to have it use the name
> SingularValueDecomposition and to rename the existing class
> TruncatedSVD.

 This class would take arguments to distinguish between
> Compact decomposition (automatic cutoff at zero value) and truncated
> decomposition (cutoff by singular values indices).

I like the optional constructor arguments approach better - no need
for a new impl class and no need to change the interface.

>>> I propose to compute U ... without taking advantage of V.  That means
>>> calling EigenDecomposition a second time but should work even in case
>>> of singular matrices.  That is the solution I am working on.  However,
>>> doing so, I notice that EigenDecomposition also suffers major problems
>>> in case of singular matrices.  A 3x3 singular matrix where 0 is an
>>> eigen value with multiplicity 2 ... yields only 2 distinct eigen
>>> vectors.  The vectors associated to the null eigen value are equal!!
>> Yes, this is JIRA issue MATH-333. The problem is that in the current
>> implementation, the vector is computed by
>> EigenDecomposition.findEigenVector which takes one eigenvalue as its
>> argument. so eigenvalues with order greater than 1 simply result in the
>> same computation repeated several times ...
>> Perhaps one way to solve this is to reproduce what is done in DLARRV
>> from LAPACK. The routine uses the index of the eigenvalue and not only
>> its value.
> Does someone had a look at DLARRV and could explain how it works (the
> vector part only, the singular values part is already well understood I
> think) ?
> I think we should wait for this issue to be solved before publishing
> 2.1. It is currently targeted at 2.2 but I would vote to solve it
> earlier. It is really a big drawback of the current implementation. Not
> being able to decompose identity is rather sad ...

>>> So, before I can improve SVD, I have to improve EigenDecomposition!
>>> By the way, going through SVD, EigenDecomposition, I noticed that
>>> BidiagonalTransformer and TridiagonalTransformer both use the
>>> Householder vector computation deeply imbedded in their code.  In
>>> order to make both classes easier to read (and to debug), I wonder
>>> if it might be useful to introduce a class Householder which would
>>> take care of the computation of the vector in a unique place.
>> This would be a nice improvement and would probably be useful for other
>> linear algebra algorithms later on. The rationale for the current
>> implementation was to avoid copying data back and forth between a low
>> level elementary Householder class called n times and a high level
>> transformer like bi-diagonal or tri-diagonal transformer. If we can find
>> a way to have an elementary transform  processing in-place the array
>> provided to it by the  high level transformer, this would be fine.
> Perhaps this could wait for 2.2 or 3.0.

+1 for waiting on this unless Dimitri or whoever ends up patching
MATH-333 finds that pulling out the Householder decomp makes it
easier to complete and test the fix.

> Luc
>> Once again, we see we lack a way to have partial view of matrices slices.
>> Luc
>>> Regards,
>>>  Dim.
>>> ----------------------------------------------------------------------------
>>> Dimitri Pourbaix                         *
>>> Institut d'Astronomie et d'Astrophysique *      Don't worry, be happy
>>> CP 226, office 2.N4.211, building NO     *         and CARPE DIEM.
>>> Universite Libre de Bruxelles            *
>>> Boulevard du Triomphe                    *      Tel : +32-2-650.35.71
>>>  B-1050 Bruxelles                        *      Fax : +32-2-650.42.26
>>>     *
>>> ---------------------------------------------------------------------
>>> To unsubscribe, e-mail:
>>> For additional commands, e-mail:
>> ---------------------------------------------------------------------
>> To unsubscribe, e-mail:
>> For additional commands, e-mail:
> ---------------------------------------------------------------------
> To unsubscribe, e-mail:
> For additional commands, e-mail:

To unsubscribe, e-mail:
For additional commands, e-mail:

View raw message