|Re: [eigen] LDLt and LLt fixes|
[ Thread Index |
| More lists.tuxfamily.org/eigen Archives
- To: eigen@xxxxxxxxxxxxxxxxxxx
- Subject: Re: [eigen] LDLt and LLt fixes
- From: Benoit Jacob <jacob.benoit.1@xxxxxxxxx>
- Date: Mon, 30 Mar 2009 22:22:21 -0400
- 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=G20HwoZwvPrPqZNC66dhJi/kgNtrkqdZx2xXSuSwmnM=; b=WOtR6FipMc3X0YoWpDoPsve6hpg6cbuZapv+QSklddXfvB7NGoIOtAj8viTdE8jMfU MfCG9o4vb3ReXiv7iKJan+2M6v/HGCcBYh4iDcvC4rg8pOVl4PsTvze8fVEF6h1hZweW sHPH+Azwv67U1crD0378PSwP9sJbazoTJsuho=
- 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=jpj3AGubKmozqC2sYX6YiZaAmqIdXabo98Yn6dh72e7Cylll+q+Z/OnPfth8rthoOV p5L0FMOTMYUOyWr5vfZvZsptyKJqhdQqbTgTr4E0D21plekblemsF+oz8dpf1eJWsDwj ewOzhKGDNT06a5ruYoNak7t9NsFYykRN8sIWI=
Thanks for the testcase. First, I just want to clarify that our LDLt
does diagonal pivoting which is all the pivoting that a Cholesky dec
In our unit tests (which succeed) we have random positive and negative
semidef matrices, so I'm confident that our code handles them well.
Your matrix is:
which has trace 1 and determinant -4, hence is indefinite (it has one
positive and one negative eigenvalue). So it is not surprising that
our LDLt fails on it: we only claim to support semidefinite matrices
(pos or neg).
What is more surprising is how other libraries can do, I'd like to
know, but I wont be able to investigate this as I have a ton other
things to do in Eigen first !
2009/3/30 Bill Greene <w.h.greene@xxxxxxxxx>:
> I think I accidentally sent you a test case that actually *requires*
> pivoting because the (0,0) term in the
> matrix is zero. So please ignore it if you like.
> On Mon, Mar 30, 2009 at 9:12 PM, Bill Greene <w.h.greene@xxxxxxxxx> wrote:
>> Wow, thanks for the quick work!
>> I just downloaded the latest code from svn and tried a quick test..
>> Surprisingly, it didn't work. After some
>> debugging and investigation, I think this may be an unusual case. My other
>> test cases work correctly--
>> exactly as I expect. Here is the code for my failed case
>> void bathe2DOF()
>> const int n = 2;
>> Eigen::Matrix2d K, M;
>> K << 6., -2., -2., 4.;
>> M << 2., 0., 0., 1.;
>> const double omega2 = 3.0;
>> Eigen::Matrix2d A = K - omega2*M;
>> Eigen::LDLT<Eigen::Matrix2d> factA(A);
>> std::cout << "Is invertible: " << factA.isInvertible() << std::endl;
>> Eigen::Vector2d b, x;
>> b << 0.0, 10.0;
>> factA.solve(b, &x);
>> prtMat(x, "x");
>> The correct value for x, by the way, is [-5.0, 0].
>> The lapack routines dsytrf and dsytrs compute this solution. I promised to
>> provide more details of the lapack
>> solution-- and I intend to-- but I'm somewhat puzzled by the factored
>> matrix returned by dsytrf. Even though
>> the solution from dsytrs is correct, the factored matrix is not what I
>> expected-- or perhaps its just that the
>> documentation is so cryptic I don't know what I'm seeing. At any rate,
>> I'll provide more lapack details as
>> soon as I understand them.
>> On Mon, Mar 30, 2009 at 5:55 PM, Benoit Jacob <jacob.benoit.1@xxxxxxxxx>
>>> 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...