|Re: [eigen] [patch] LDLt decomposition with rank-deficient matrices|
[ Thread 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=
On Sat, Jun 5, 2010 at 11:34 PM, Ben Goodrich <bgokgm@xxxxxxxxxxxxxx> wrote:
> 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.
> 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
> ... // 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?
# 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 @@
+ /** \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()),
@@ -302,7 +307,11 @@
*sign = 1; // this is by assumption rather than calculated
+ // partition the matrix:
+ // A00 | - | -
+ // lu = A10 | A11 | -
+ // A20 | A21 | A22
Index rs = size - k - 1;
@@ -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)