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 ... +1 > >>> 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. Phil > > 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 >>> http://sb9.astro.ulb.ac.be/~pourbaix * mailto:pourb...@astro.ulb.ac.be >>> >>> --------------------------------------------------------------------- >>> To unsubscribe, e-mail: dev-unsubscr...@commons.apache.org >>> For additional commands, e-mail: dev-h...@commons.apache.org >>> >> >> --------------------------------------------------------------------- >> To unsubscribe, e-mail: dev-unsubscr...@commons.apache.org >> For additional commands, e-mail: dev-h...@commons.apache.org >> >> > > > --------------------------------------------------------------------- > To unsubscribe, e-mail: dev-unsubscr...@commons.apache.org > For additional commands, e-mail: dev-h...@commons.apache.org > --------------------------------------------------------------------- To unsubscribe, e-mail: dev-unsubscr...@commons.apache.org For additional commands, e-mail: dev-h...@commons.apache.org
