Re: [eigen] [patch] LDLt decomposition with rank-deficient matrices |

[ Thread Index | Date Index | More lists.tuxfamily.org/eigen Archives ]

*To*: eigen@xxxxxxxxxxxxxxxxxxx*Subject*: Re: [eigen] [patch] LDLt decomposition with rank-deficient matrices*From*: Ben Goodrich <bgokgm@xxxxxxxxxxxxxx>*Date*: Tue, 15 Jun 2010 15:25:21 -0400*Dkim-signature*: v=1; a=rsa-sha256; c=relaxed/relaxed; d=googlemail.com; s=gamma; h=domainkey-signature:mime-version:received:received:in-reply-to :references:date:message-id:subject:from:to:content-type; bh=ibJgyVrmqA/37feqmoK17WDJ/JlkRCmGvGTMYzvBXSU=; b=W5/sO2pmo2V9fG3dHambb6yAildcR01sOdXg/XA1bDFGrfpCSJPjxHKNwEChC5wRlC Uyfwn1SHpJLFtKmd4o1FpQS3IckFS5w/zIPy6t9A8T0bzjktGc78dkO+VvIRqdzbVr3Q nmUYjpIAhoJ3cXUTXSNU355MtKjDlL5u1gbpw=*Domainkey-signature*: a=rsa-sha1; c=nofws; d=googlemail.com; s=gamma; h=mime-version:in-reply-to:references:date:message-id:subject:from:to :content-type; b=CV6kvv+yHMiVZJvlwu3SqWLSbClqdZzevyBxYAvtQokKVkc6lOChL6e+LZ4fhPIpVN AWMuYWl1oVoeVndVgyVBSQbUupSm05t+1p4IdYyzF3Tjhq/OIGA5PFo3r+IURYKN9e3W XGtFfJGLuqlwCroP2QOT5ES2T1w1HhZP7rBUQ=

Hi, On Sat, Jun 5, 2010 at 11:34 PM, Ben Goodrich <bgokgm@xxxxxxxxxxxxxx> wrote: > Hi, > > On Sun, Feb 28, 2010 at 1:44 AM, Benoit Jacob <jacob.benoit.1@xxxxxxxxx> wrote: >> If I understand correctly your other email, you'd be happy with >> instead a non-pivoting LDLt. That sounds more reasonable than >> "unpivoting" a pivoting LDLt. >> >> So: OK to add a non-pivoting option to LDLt. >> >> Benoit > > Here is a unified patch to reimplement the non-pivoting option. Sorry > for the more than 3 month delay; had to write words in dissertation > rather than code :( . Has anyone tried my patch or has comments on the things below? I rebased the patch in the attachment to make sure it would still apply and to improve the documentation a little more. -- Ben > Some thoughts: > > 0) If you changed your mind and don't want this option anymore, that's > fine. I can apply the patch locally. That said, having a pivot option > makes it possible to be compatible with eigen2 and is probably useful > to someone besides just me. And it is easy to do. > > 1) Gael made some big changes to LDLT a couple of days ago that I > think I understand but might have misunderstood something when I > reimplemented the non-pivoting option on top of them. > > 2) For example, in the positive semi-definite but (numerically) > rank-deficient case, I had to put exact zeros on the diagonal of D to > get the reconstruction right. This contravenes Benoit's statement > earlier in this thread that Eigen does not "finish" decompositions of > rank-deficient matrices. But I didn't immediately see any other good > way to do it, in light of Gael's recent changes to LDLT. > > 3) I implemented the option like > > if(pivot) > { > .... > } > else > { > ... // mostly copy-and-paste from above > } > > which resulted in some code duplication but avoids the need to > repeatedly evaluate if(pivot) inside the loop. I guess that is > slightly faster and clearer, but if you would rather not have the code > duplication, I can do it the other way. > > 4) I made some changes to the documentation explaining the pivot > option, but just noticed (and haven't really read yet) today's long > documentation thread. > > 5) I copied-and-pasted a block inside /test/cholesky.cpp and exercised > the pivot=false option. It seems to work when you do ./check.sh > cholesky. I did some other tests locally with singular matrices, but > /test/cholesky.cpp does not seem to have any tests with singular > matrices, so maybe some should be added? > > Thanks, > Ben >

# HG changeset patch # User Ben Goodrich <bgokGM@xxxxxxxxx> # Date 1276629596 14400 # Branch pivot # Node ID c6bf0801d6e0f89fff5d7168ce97ce00e70bf58c # Parent ee38d0fe25eca4b019f2544ab05c1aa09c213d67 LDLT: Introduce and test pivot option, which defaults to true, but when false decomposes the matrix without pivoting, similar to eigen2 behavior. diff -r ee38d0fe25ec -r c6bf0801d6e0 Eigen/src/Cholesky/LDLT.h --- a/Eigen/src/Cholesky/LDLT.h Sat Jun 05 22:58:36 2010 -0400 +++ b/Eigen/src/Cholesky/LDLT.h Tue Jun 15 15:19:56 2010 -0400 @@ -44,7 +44,7 @@ * The decomposition uses pivoting (by default) to ensure stability, so that L will have * zeros in the bottom right rank(A) - n submatrix. Avoiding the square root * on D also stabilizes the computation. If the pivoting option is not utilized, \f$ A \f$ - * must be positive semi-definite by construction, which is NOT checked. In that positive semi-definite + * must be positive semi-definite by construction, which is NOT checked. In the positive semi-definite * case, numerical problems may arise if \f$ A \f$ is rank-deficient and pivoting is not utilized. * * Remember that Cholesky decompositions are not rank-revealing. Also, do not use a Cholesky @@ -98,6 +98,11 @@ m_isInitialized(false) {} + /** \brief Default Constructor with decomposition + * + * Performs the decomposition of \a matrix with pivoting by default, which + * can be changed by setting \a pivot to false. + */ LDLT(const MatrixType& matrix, const bool pivot=true) : m_matrix(matrix.rows(), matrix.cols()), m_transpositions(matrix.rows()), @@ -302,7 +307,11 @@ if(sign) *sign = 1; // this is by assumption rather than calculated } - + + // partition the matrix: + // A00 | - | - + // lu = A10 | A11 | - + // A20 | A21 | A22 Index rs = size - k - 1; Block<MatrixType,Dynamic,1> A21(mat,k+1,k,rs,1); Block<MatrixType,1,Dynamic> A10(mat,k,0,1,k); @@ -364,6 +373,7 @@ }; /** Compute / recompute the LDLT decomposition A = L D L^* = U^* D U of \a matrix + * with pivoting by default. */ template<typename MatrixType, int _UpLo> LDLT<MatrixType,_UpLo>& LDLT<MatrixType,_UpLo>::compute(const MatrixType& a, const bool pivot)

**Follow-Ups**:**Re: [eigen] [patch] LDLt decomposition with rank-deficient matrices***From:*Benoit Jacob

**References**:**Re: [eigen] [patch] LDLt decomposition with rank-deficient matrices***From:*Ben Goodrich

**Messages sorted by:**[ date | thread ]- Prev by Date:
**Re: [eigen] Matrix product between some MatrixType and triangularView<UnitLower>** - Next by Date:
**Re: [eigen] Matrix product between some MatrixType and triangularView<UnitLower>** - Previous by thread:
**Re: [eigen] [patch] LDLt decomposition with rank-deficient matrices** - Next by thread:
**Re: [eigen] [patch] LDLt decomposition with rank-deficient matrices**

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