|[eigen] LDLt and LLt fixes|
[ Thread Index |
| More lists.tuxfamily.org/eigen Archives
- To: eigen@xxxxxxxxxxxxxxxxxxx
- Subject: [eigen] LDLt and LLt fixes
- From: Benoit Jacob <jacob.benoit.1@xxxxxxxxx>
- Date: Mon, 30 Mar 2009 23:55:46 +0200
- Dkim-signature: v=1; a=rsa-sha256; c=relaxed/relaxed; d=gmail.com; s=gamma; h=domainkey-signature:mime-version:received:date:message-id:subject :from:to:content-type:content-transfer-encoding; bh=QxCJvPbWZ3o6Yuun1f2I7oDWHeOSW5tgVEO63HiUok4=; b=SKe7moKrhJmVEsGI+ZMoaZjOs7WEMPiituk+sEznbkjcEFqcMsg7sbM8mTxnwGR8FN HzeQLqhor9GU98Mo57VpS3ytx2ZrH41v3SCHW6Ia4N9f0d543D0ylvVJOl5zrjvQJeLc 0HtaZrBh/AQsljPZ4qzrCXLD8JUjBiqXI7YVM=
- Domainkey-signature: a=rsa-sha1; c=nofws; d=gmail.com; s=gamma; h=mime-version:date:message-id:subject:from:to:content-type :content-transfer-encoding; b=LqfExtfuQX/5F4qh24GBZsRF5LynNLVLboF6McAGw/jr66QSgwhGYhwJPehCAO+SGW XI4YjKTt+sk68Nb2/sc8/iFpkVWvldaLraLT5xkvxi5I9l5dwWui4bfpeakrAIKwL4pD AKOQffZyG9Ug2R13DauhNDLf2VoZwES33gqUY=
the thread with Bill inspired me... so now (trunk revision 947097),
LDLt supports the negative semidefinite case.
There's a couple of things to discuss.
First, the floating-point comparisons in both LLt and LDLt had issues. LLt did
if(x < eps)
and LDLt did
if(Djj <= 0)
and replacing these by proper fuzzy comparisons (experimenting with a
few different ones) resulted in a big improvement in the accuracy of
some methods, like LDLt::rank() which used to be untested and had
typically 3% imprecision, and now is exact in 95% of cases (However
for an exact rank computation, it's still not nearly as reliable as
the full-pivoting LU which can be explained by the different between
bi-directional pivoting and the diagonal pivoting that's all one can
do with Cholesky).
The other issue is that LLt and LDLt have claims in their dox that
they use only a triangular half of the matrix and don't need the other
half to be initialized, well that's not currently the case, one can
check that by uncommenting a line I added to the cholesky.cpp
unit-test. Help fixing this welcome...