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, 23 Feb 2010 21:29:10 -0500*Dkim-signature*: v=1; a=rsa-sha256; c=relaxed/relaxed; d=googlemail.com; s=gamma; h=domainkey-signature:mime-version:received:in-reply-to:references :date:message-id:subject:from:to:content-type; bh=/G/YcLSu4X3ZLZhc9WhS5WssHzSMNr8mcDkU6mxn1SQ=; b=Q9e/ByQvg9gpfQfJaeIjYCY6wYEnSYLJ+2VfOAjWzFljPRoMLzEQm7rjL4MgIvMtYG ygSJtojkrapmVNU+fe19t+qtELmcKEk2MVFYYTZRppE0HsU+HxTmAanLtucikp/ClwQz W3AgDAG4wo5U7207kxOL5k0RwHEqA53aFMGkg=*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=AHpbs90GvIEJHKdLFHNj+KjPjrMxSaHdNK0IlW5Ygl9iU7W59rb/fGhaX8jHBRJYn0 ubeh1PvPY98hAHQwE70lZwH70p+oB35I8fvOs2vfhAEGmhfCEd+DeJ6Fp6jQEi6WUZUZ eUn/qVXtacC7EKh2gQPtWFULu0UGRGmFgjE/4=

Hi, I have found another robustness problem with ldlt() in the current development branch. Some combination of rank-deficiency, constant diagonal, pivoting, and double precision can make it impossible to roundtrip a matrix through the LDLt decomposition. I am attaching another test case that yields the output below. Thank you, Ben A = LDL' has a unit diagonal (apart from numerical noise) 1 1 1 1 1 This is a numerical disaster (compare the second cells). Here is the true D 1 0.915092 0 0 0 Here is the calculated D 1 0.0785762 2.01553e-16 7.97973e-17 1.39591e-15 This is good (because the diagonal of A is forced to have 1.0 in all its cells). Here is the true D 1 0.915092 0 0 0 Here is the calculated D 1 0.915092 -6.06069e-16 -8.0296e-16 -8.4893e-16

#include "eigen/Eigen/Cholesky" #include<iostream> using namespace Eigen; using namespace std; extern "C" { int main() { LDLT<MatrixXd> A_factored; VectorXd D(5); D << 1.0, 0.915092047236159, 0.0, 0.0, 0.0; MatrixXd L(5,5); L << 1.0, 0.0, 0.0, 0.0, 0.0, 0.2913896922745232, 1.0, 0.0, 0.0, 0.0, 0.3506763450086252, 0.9789801109802591, 1.0, 0.0, 0.0, -0.0551313984496821, 1.0437742813828543, 0.0, 1.0, 0.0, -0.2269670280948324, -1.0180827163703894, 0.0, 0.0, 1.0; MatrixXd A = L * D.asDiagonal() * L.transpose(); A_factored = A.ldlt(); cout << "A = LDL' has a unit diagonal (apart from numerical noise)" << endl << A.diagonal() << endl; cout << "This is a numerical disaster (compare the second cells)." << endl; cout << "Here is the true D" << endl << D << endl << "Here is the calculated D" << endl << A_factored.vectorD() << endl; for(int i=0; i < 5; i++) A(i,i) = 1.0; A_factored = A.ldlt(); cout << "This is good (because the diagonal of A is forced to have 1.0 in all its cells)." << endl; cout << "Here is the true D" << endl << D << endl << "Here is the calculated D" << endl << A_factored.vectorD() << endl; return 0; } }

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

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

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

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

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

**Messages sorted by:**[ date | thread ]- Prev by Date:
**Re: [eigen] [patch] LDLt decomposition with rank-deficient matrices** - Next by Date:
**Re: [eigen] a branch for SMP (openmp) experimentations** - 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/ |