Re: [eigen] [patch] LDLt decomposition with rank-deficient matrices

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

2010/2/24 FMDSPAM <fmdspam@xxxxxxxxx>:
> Thank you Benoit,
>
> (Sorry for my wired sentences in the previous mail. please forget the noise,
> here is my second attempt)
>
> Am 24.02.2010 14:26, schrieb Benoit Jacob:
>>
>> 2010/2/24 FMDSPAM<fmdspam@xxxxxxxxx>:
>>
>>>
>>> Hi Experts,
>>>
>>> please, for my personal clarification, let me ask:
>>> (I'm not asking for a brief lesson in lin alg. but I have to happily live
>>> with those answers :-) )
>>>
>>> Def.:
>>>    D_labil := "how many zeros in vectorD()"
>>>    D_instable : "how many entries in vectorD() are smaller or equal 0.0"
>>>
>>> a) LDLt for non singular matrix is rank-revealing (of course, i hope),
>>>
>>
>> This is a tautology: if you assume your matrices to be nonsingular,
>> then you know their rank. Nonsingular means rank N.
>>
>>
>
> Yes, of course. It was listed for the sake of completeness.
>>>
>>> n_labil = 0
>>>
>>
>> Yes, for (sufficiently) nonsingular matrices, there won't be too small
>> entries in the vectorD.
>>
>>
>>>
>>> b) LDLt for positive definite matrix =>  D_instable = 0
>>>
>>
>> Yes.
>>
>>
>>>
>>> c) LDLt for singular matrix is not rank-revealing, but: n_deficit>  0
>>>
>>
>> What is n_deficit ? Number of zeros in vectorD?
>
> Yes exactly, sorry.
>>
>>  I think this kind of
>> works, yes, as you're only interested in knowing if it's 0 or>0. As
>> the LDLt decomposition still gives you the determinant, it in
>> particular tells you very vaguely if the matrix is singular or not,
>> but don't use it for that, it shouldn't be very precise.
>>
>>
>
> A better alternative exist, right? (see my (ab-)use case below).

If you want to measure the rank (deficit), you need a rank-revealing
decomposition. For example:
 - full-pivoting LU
 - column-pivoting QR (yes it's enough)
 - full-pivoting QR
 - SVD

>>>
>>> d) LDLt for indefinite matrix =>  D_instable>  0  but counting D(i)<= 0.0
>>> is
>>> meaning less?
>>>
>>
>> Yes. The "zeros" in the vectorD will actually be small positive
>> values. There's no magic cutoff where you'd decide that a small value
>> can be considered zero.
>>
>>
>
> OK. Maybe from a practical point of view: D(i,i) < D.max*1e-8 could
> considered as "small value"? But that may depend on concrete use cases.

I don't know, this could be trickier depending on how big the
off-diagonal elements of the L matrix can get. All I can do is point
you to the above Higham article, I'm not an expert myself.

> My (ab-)use case or "how semi-education cause numerical nightmares":
>
> I have infiltrated a very ancient but nevertheless successfully just in use
> Fortran FE-code with my own taucs-LLt-solver binding (with a very welcome
> performance boost).
> The  taucs-LLt-solver itself use blockwise a traditional Lapack-solving
> routine.
>
> In that Lapack solving routine, I do observe and count "nearly" negative or
> fully diagonal elements (A(i,i) < +eps).
> If I detect a "negative" diagonal elements, I increment my a counter (say
> "++D_instabile") and brutally set the diagonal element to A(i,i) = +1.0.. the
> rest of the Lapack routine thereon is used without further modifications.
> That count of "negative" diagonal elements is used as an indicator for an
> outer iteration loop to find the instability limit (D_instabile = 0 <->
> D_instabile=1).
>
> This algorithm is in use, but tend to be unstable for bigger number of
> degrees of freedom.
> Now, thank you, I have several new staring points.

Good luck :)

Benoit

>
> Frank
>
>
>