Re: [eigen] LDLt and LLt fixes |

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

*To*: eigen@xxxxxxxxxxxxxxxxxxx*Subject*: Re: [eigen] LDLt and LLt fixes*From*: Keir Mierle <mierle@xxxxxxxxx>*Date*: Mon, 30 Mar 2009 19:27:30 -0700*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; bh=7fGzY+48YUZbNu/jzTBF+Z3CowtFN6BAX+IttU++2to=; b=xLWEZ5RZfP8oscJcqaVagABjwIJSL/BWA0cHuxsPc3VcLsFYYvX/8XhcZlAhx3qh0C vyXh8l5OJkuceyoNz5zbY9s7QutEhpbmDlN0BrSfVwPFCs6GexuJ56RjO/Eatkg95nNq uvDYWODTTdeGe/nw3VMkYdpQeXKO4yGG1cDvk=*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; b=guT9IokpIlRApQGXUEH6FGLX9IL/XQr7HSVO6FxIebfF1gOqgopVRW5tMxQ3KPOhtf h37/r+4bayRuXf4KPWGUq+jRQnko/lvtqVhWI2by3H8tBKivpdsHvHVXiNIqNo5MQtlA hBJ/zNKgFopGU1doAQ4er1rwGfyLNAU9z0SZk=

On Mon, Mar 30, 2009 at 7:22 PM, Benoit Jacob <jacob.benoit.1@xxxxxxxxx> wrote:

Note that the "Cholesky" of an indefinite matrix results in a block-diagonal D with 1x1 and 2x2 blocks. This is known as the block Cholesky. If we want to support that, it should be in an entirely different routine because it will ruin the performance of normal LDLt due to all the special casing.

Also, in the case of a semidefinite matrix, a pivoting cholesky implementation will produce a PLDLPt decomposition such that PL has a triangular part that is zeros at the bottom.

See page 2 of http://eprints.ma.man.ac.uk/1193/01/high90c.pdf

Keir

Bill,

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

can do.

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:

0 -2

-2 1

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 !

Note that the "Cholesky" of an indefinite matrix results in a block-diagonal D with 1x1 and 2x2 blocks. This is known as the block Cholesky. If we want to support that, it should be in an entirely different routine because it will ruin the performance of normal LDLt due to all the special casing.

Also, in the case of a semidefinite matrix, a pivoting cholesky implementation will produce a PLDLPt decomposition such that PL has a triangular part that is zeros at the bottom.

See page 2 of http://eprints.ma.man.ac.uk/1193/01/high90c.pdf

Keir

Benoit

2009/3/30 Bill Greene <w.h.greene@xxxxxxxxx>:

> Benoit,

>

> 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.

>

> Bill

>

> On Mon, Mar 30, 2009 at 9:12 PM, Bill Greene <w.h.greene@xxxxxxxxx> wrote:

>>

>> Benoit,

>>

>> 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.

>>

>> Bill

>>

>> On Mon, Mar 30, 2009 at 5:55 PM, Benoit Jacob <jacob.benoit.1@xxxxxxxxx>

>> wrote:

>>>

>>> Hi,

>>>

>>> 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...

>>>

>>> Cheers,

>>> Benoit

>>>

>>>

>>

>

>

**Follow-Ups**:**Re: [eigen] LDLt and LLt fixes***From:*Benoit Jacob

**References**:**[eigen] LDLt and LLt fixes***From:*Benoit Jacob

**Re: [eigen] LDLt and LLt fixes***From:*Bill Greene

**Re: [eigen] LDLt and LLt fixes***From:*Bill Greene

**Re: [eigen] LDLt and LLt fixes***From:*Benoit Jacob

**Messages sorted by:**[ date | thread ]- Prev by Date:
**Re: [eigen] LDLt and LLt fixes** - Next by Date:
**Re: [eigen] sse asin implementation** - Previous by thread:
**Re: [eigen] LDLt and LLt fixes** - Next by thread:
**Re: [eigen] LDLt and LLt fixes**

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