Re: [eigen] Instability in LLT and LDLT methods. |

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

*To*: eigen@xxxxxxxxxxxxxxxxxxx*Subject*: Re: [eigen] Instability in LLT and LDLT methods.*From*: Benoit Jacob <jacob.benoit.1@xxxxxxxxx>*Date*: Wed, 28 Jan 2009 14:35:28 +0100*Dkim-signature*: v=1; a=rsa-sha256; c=relaxed/relaxed; d=gmail.com; s=gamma; h=domainkey-signature:mime-version:received:in-reply-to:references :date:message-id:subject:from:to:content-type :content-transfer-encoding; bh=H7PbI0cpdOGa5f4rCkWjjbNr/jCiY0gh41wETuLsM8A=; b=DNJtwWqYYlVFvWcvdq9fOHl4P6yN5/cKLYyBE6L2wvmjkGCH8r4ySXlg4NNDHGwx6t f2ammZHmSkY4oZNG+EfMiPnDWmMq1cn4zxFt0fkg6ObVHkKJOA/qf9QN2qM4Hkqny4zw sWMUAoFu+IwIozj8YF5g0SoiebZfPLmoK5FCc=*Domainkey-signature*: a=rsa-sha1; c=nofws; d=gmail.com; s=gamma; h=mime-version:in-reply-to:references:date:message-id:subject:from:to :content-type:content-transfer-encoding; b=jdvIwVwtHQIj/NnxS31z78O0+Pz53gxkv97EmIW1oAjsPKFPHEXK5TCzZj1bY+Hho/ oMCNVDSilY6Srtio3/zF8K0cm0h01w5/ZBUOcATL4mzJhZAk+5hLilqA5SkMB3C9Kd5M siEPQewtf4J17VO3D34r9XOZCezPKeCGpRFng=

2009/1/28 Jitse Niesen <jitse@xxxxxxxxxxxxxxxxx>: > On Tue, 27 Jan 2009, Keir Mierle wrote: > >> Probably it's better to do full pivoting. Apparently cholesky is >> stable for semidefinite matrices when full pivoting is used: >> >> http://eprints.ma.man.ac.uk/1101/01/covered/MIMS_ep2008_56.pdf > > Cholesky (LL^T) is stable without any pivoting for positive definite > matrices. So I think the possibility of Cholesky without pivoting should > remain. Here is the thing that's got me thinking since yesterday and why I didn't reply so far. - on the one hand, since pivoting is essential for the stability of LU, and LLt is a close variant of LU, I would assume that pivoting is also important for the stability of LLt. - on the other hand, something very special happens with LLt: the matrix L is _unique_. Since it is well-determined mathematically, one would expect it to be computable in a stable way. Now, if we do pivoting, it's no longer a LLt decomposition but rather a (PL)(PL)t decomposition with P a permutation matrix. So implicitly if the L of LLt is computable in a stable way (and I believe it is), it has to be doable without pivoting. This argument doesn't apply to LU as, since the L and U aren't unique, one can't talk of "the L of LU" being computable. Just to say that I'm puzzled by LLt and I don't know whether or not pivoting is needed here (by contrast at least for LU and SVD pivoting is clearly needed). > Why do you need an epsilon? Perhaps a better question is, why does the > routine test whether the matrix is positive definite? As you say, it's not > easy to write such a test, and it could be said that it's the user's > responsibility that the matrix is positive definite - though of course a > test is nice when the code is run in debug mode. I agree. If an algorithm works only for positive definite matrices, it has to be the responsibility of the user. Moreover we should not check for that. Just like the compiler doesn't check for x!=0 when we do 1/x. Also, here by "positive definite" we have to mean "the ratio of the biggest to the lowest eigenvalue is not too big" or if you prefer "the lowest eigenvalue is not negligible when compared to the biggest eigenvalue" i.e. the notion of positive definite has to be "fuzzy". > A related issue I've been wondering about is: why does LU decomposition use > complete pivoting? I believe that's quite a bit more expensive than partial > pivoting It's about 2x more expensive currently, as the pivot lookup part is not yet vectorized. If it were, it would still be 1.5x more expensive. > When solving a linear system with > an invertible matrix, partial pivoting performs as well as complete > pivoting, doesn't it? LAPACK uses partial pivoting, as does GSL. No, partial pivoting fails even for invertible matrices of size 200x200, and fails consistently at sizes like 1000x1000. I think that LAPACK and GSL do a hybrid of partial and full pivoting that I have yet to understand. In any case, anything less than full pivoting will fail in certain cases. I believe that LAPACK's implementation will only fail on specially tailored matrices, but still, if I understand how their hybrid pivoting works, I can construct an invertible matrix making it fail, even without specially weird coefficients. So I'm not comfortable with this, all this insecurity for, what, at 50% speedup. Plus, a user on the forum reports that our current implementation is already as fast as LAPACK for size 1500x1500. Cheers, Benoit

**Follow-Ups**:**Re: [eigen] Instability in LLT and LDLT methods.***From:*Jitse Niesen

**References**:**[eigen] Instability in LLT and LDLT methods.***From:*Keir Mierle

**Re: [eigen] Instability in LLT and LDLT methods.***From:*Gael Guennebaud

**Re: [eigen] Instability in LLT and LDLT methods.***From:*Keir Mierle

**Re: [eigen] Instability in LLT and LDLT methods.***From:*Jitse Niesen

**Messages sorted by:**[ date | thread ]- Prev by Date:
**Re: [eigen] internal compiler error in cachefriendlyproduct.h** - Next by Date:
**Re: [eigen] internal compiler error in cachefriendlyproduct.h** - Previous by thread:
**Re: [eigen] Instability in LLT and LDLT methods.** - Next by thread:
**Re: [eigen] Instability in LLT and LDLT methods.**

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