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 12:27:17 -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=XaPjRwhL8GZvmUEcEv5vqPII7XNdBi+T6TB6yqsAcGk=; b=UpMbUGmdvVfM326kJCsWRRAptRJqwbamSuVm5HsHh5PU08EyXDbhhAK1kMyjhuiAPt g7/x2JVQ75FiMWWpAWbWoOFrhRNOutnM9GuFxA+NhSqRQ5ulAq1GCl6pb52KQ1ySn+J9 nfRpZ4VUSEOlv9EOX2REsqnWW1w2FnVBqtvUU=*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=VRmuO0QtuZkefsj2mZADTSCF11E5xarZv5Hs2lr+IxLS7GoAF0dUyRPvKf9KBt7lNS VQEhyJ37Pla1colRjbiIWeViYi+648P5J6jfc9LKZWNB9bwtE7r4fRSyhBX7KwNQxx+W uH3nVVdreZnDMtcsvAMH0kNVxF4coCDsba5W4=

Yes, as your example shows clearly, you can't guarantee an LDLt factorization

for symmetric, non-positive definite matrices without pivoting. However, I think

in many practical cases-- at least the ones I have experience with-- you can

do the basic LDLt factorization. In the mechanical example I mentioned previously,

the matrix to be factored is computed as the difference of two SPD matrices,

A = K - alpha*M. You can certainly pick values of alpha that make A singular and

the LDLt factorization doesn't exist. But as long as you're not "too close" to this

singular case, LDLt with no pivoting seems to work quite well.

Bill

On Mon, Mar 30, 2009 at 9:21 AM, Benoit Jacob <jacob.benoit.1@xxxxxxxxx> wrote:

Oops. In fact, invertibility won't help, and I'm not sure at all that

LDLt (even pivoting) exists for all self-adjoint matrices. I mean, I

read here,

http://en.wikipedia.org/wiki/Cholesky_decomposition#Avoiding_taking_square_roots

that

"this factorization can be used for any square, symmetrical matrix"

but I don't see anything there backing this assertion, it seems to me

that it's wrong, and at least the formulas that they give fail for the

matrix

0 1

1 0.

Cheers,

Benoit

>So a generalization of LDLt handling non-positive matrices would have

>to assume invertibility. It would work for self-adjoint invertible

>matrices. Why not. Are you still interested in that?

>The next question would then be: implement as separate dec or as a

>variant of LDLt controlled by a template parameter?

2009/3/30 Benoit Jacob <jacob.benoit.1@xxxxxxxxx>:

> I had a closer look.

>

> What happens is that the "positive" assumption allows to restrict the

> pivot lookup to the diagonal. So we get full pivoting at the price of

> partial pivoting.

>

> What allows us to do that is that the norm of a positive matrix is

> bounded by a multiple of the norm of its diagonal. In particular, if

> the diagonal is zero, then the whole matrix is zero.

>

> This is really characteristic of positive matrices, it doesn't work

> for general self-adjoint matrices, e.g.

> 0 1

> 1 0

>

> So a generalization of LDLt handling non-positive matrices would have

> to assume invertibility. It would work for self-adjoint invertible

> matrices. Why not. Are you still interested in that?

>

> The next question would then be: implement as separate dec or as a

> variant of LDLt controlled by a template parameter?

>

> By the way, Keir, at line 162:

> m_matrix.block(j, j, size-j, size-j).fill(0); // Zero unreliable data.

> Unless this is backed by experiment, I kind of disagree with this

> line, I don't see how zeroing this noise can be useful, and I can

> think of disadvantages of zeroing: it makes the reconstruction of the

> original matrix more inaccurate, not to mention that it generates more

> code.

>

> Cheers,

> Benoit

>

>

> 2009/3/30 Benoit Jacob <jacob.benoit.1@xxxxxxxxx>:

>> This is interesting. Indeed it seems that this works for all

>> self-adjoint matrices. I'll try to have a look at it.

>>

>> Cheers,

>> Benoit

>>

>> 2009/3/29 <w.h.greene@xxxxxxxxx>:

>>> I just finished reading a long recent thread on some issues with LDLt.

>>>

>>> In some numerical experiments I was doing today, I've come across a

>>> related issue. It appears that the current implementation of LDLt requires

>>> that the matrix be positive definite. This is not strictly necessary for

>>> the LDLt factorization to succeed. It is necessary only that the diagonal

>>> term not be zero during factorization. A negative diagonal term is not,

>>> by itself, a problem.

>>>

>>> From reading the previous posts, it appears that a main reason for the

>>> LDLt implementation compared with LLt was to avoid the performance

>>> penalty of a square root. In fact, I think the main benefit of LDLt is that

>>> it can factor both negative- and positive-definite symmetric matrices.

>>> (As an aside, these routinely occur in dynamic analysis of mechanical

>>> systems)

>>>

>>> I suggest the test in the LDLt factorization routine be changed to

>>> check that abs(diag_term) > eps.

>>>

>>> Bill Greene

>>

>

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

**Messages sorted by:**[ date | thread ]- Prev by Date:
**[eigen] Re: discussion on special matrices** - 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/ |