Re: [eigen] Another LDLt issue |

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

*To*: eigen@xxxxxxxxxxxxxxxxxxx*Subject*: Re: [eigen] Another LDLt issue*From*: Bill Greene <w.h.greene@xxxxxxxxx>*Date*: Mon, 30 Mar 2009 14:45:02 -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; bh=NhLZl0fI0ird3fEZSQbPBB9ulBVJQ04l0Trsoo0JboU=; b=vU01fYClvqwihUs/+zO5kMAq0QLiYZI0PdnolhQvPGR/WlvHlMT+Jcs8v3lP/FzyO+ +Hr5Sv6IiQE8D3B8Ohh28USeRQv/qyENTYpgHMT9MKJDr4+ne3qctyIAbKVKeQCva1Ok ysMzaQE0btz3cEJKkGPebw0CECKos5+pcDpPU=*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=kSFIUsnRetMqa8TSgXoPED4bzSCchYktkOqppNVgba/9lyD/acY8ysNRZ7y+tMS1YO kuTiDPwmga4Bgl8pJjVr/cHjZ3gpk55FWtmqt4xs0ZVVavhoZ5mLIzj9ZREIOLz04/M4 W9qMjL/EGE0QkRPtZ1fk/GFfkdDutcZ529oPE=

Benoit,

Just a couple of notes.

>No, I was still thinking about this 2x2 example,

>0 1

1 0

>that matrix is nonsingular and selfadjoint, yet its diagonal is 0.

Yes, but pivoting *will* allow this matrix to be factored. When I say

"pivoting" here I mean changing the order of the variables. Perhaps

you mean partial pivoting?

>Our LDLt does the same.

No, I don't think so. dsytrf succeeds for the negative definite case

and Eigen::LDLt simply stops. It isn't the pivoting that allows dsytrf

to proceed-- at least in my example-- it just doesn't consider a negative

diagonal term to be an error. I can produce a small test case if that

would be helpful?

Bill

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

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

> Benoit,Ah , OK.

>

>>I see. So we can always provide general selfadjoint LDLt "at the

>>user's own risk". The generic case is that it will "work" since

>>division by exactly 0 would be incredibly bad luck, but it will be

>>unstable because one will be dividing by small numbers.

>

> I wasn't proposing continuing the factorization if a diagonal term

> is (or becomes) very small.

Yes, it should. But if the user has to test that the factorization

> I just propose that the test be

> if(abs(diag) < eps) return. This should insure that the algorithm

> doesn't compute a bad result, shouldn't it?

succeeded and in the converse case use another factorization, then I

am afraid that he will prefer doing the other factorization right away

! Especially if that other choice is going to be more stable and not

much slower (i'm thinking partial-pivoting block LU -- we could also

implement block Cholesky but that's yet more work and we don't have

enough arms).

Indeed, I was saying that's the generic case.

> My own experience

> is that this will succeed often enough to make it useful.

No, I was still thinking about this 2x2 example,

>> >But the problem is that pivoting won't help here if all the remaining

>>

>> >diagonal coefficients are zero. Indeed, a meaningful pivoting LDLt is

>>

>> >of the form:

>

> Are you talking about a case where an LDLt factorization exists but the

> matrix is singular?

0 1

1 0

that matrix is nonsingular and selfadjoint, yet its diagonal is 0.

The basic use case is solving a system of equations,

> To be honest, before starting this current thread, I didn't even realize

> such a case could

> exist. What do people use this factorization for-- calculating matrix rank?

The main issue is that in computer / floating-point terms, all

> I think most

> people are interested in the non-singular case but, of course, would like to

> know if the

> factorization algorithm fails.

matrices are "non-singular" since the determinant is virtually never

exactly 0, the real issue is numerical stability when the matrix is

"close" to being singular.

Our LDLt does the same,

>Anyway,

> routine dsytrf does an LDLt factorization for a symmetric but not

> necessarily positive definite

> matrix.

Yes, that's exactly what our LDLt does (didn't know it was called

> They *have* chosen to include diagonal pivoting in the

> implementation; obviously a lot

> more work.

diagonal pivoting but that makes a lot of sense), see my previous

email, transforming a matrix X into PXPt has exactly the effect of

permuting along the diagonal direction.

Cheers,

Benoit

**Follow-Ups**:**Re: [eigen] Another LDLt issue***From:*Benoit Jacob

**References**:**[eigen] Another LDLt issue***From:*w . h . greene

**Re: [eigen] Another LDLt issue***From:*Benoit Jacob

**Re: [eigen] Another LDLt issue***From:*Benoit Jacob

**Re: [eigen] Another LDLt issue***From:*Benoit Jacob

**Re: [eigen] Another LDLt issue***From:*Bill Greene

**Re: [eigen] Another LDLt issue***From:*Benoit Jacob

**Re: [eigen] Another LDLt issue***From:*Bill Greene

**Re: [eigen] Another LDLt issue***From:*Benoit Jacob

**Messages sorted by:**[ date | thread ]- Prev by Date:
**Re: [eigen] Another LDLt issue** - Next by Date:
**Re: [eigen] Another LDLt issue** - Previous by thread:
**Re: [eigen] Another LDLt issue** - Next by thread:
**Re: [eigen] Another LDLt issue**

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