Re: [eigen] Another LDLt issue
• 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=

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

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