Re: [eigen] Another LDLt issue

[Benoit Jacob <jacob.benoit.1@xxxxxxxxx>]

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