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

[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).
> > 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.
> > e) Do all this also hold for LLt decompositions?
> >      
> Oops, for LLt decomposition it's even worse. First you're taking sqrt
> of small values, so you're losing half of the precision. Second, our
> LLt is non-pivoting so it's totally unstable on singular matrices.
>
>    
Strike! You sunk my ship :-) .
All that is reasonable to me.

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.

Frank