2009/3/30 Bill Greene <w.h.greene@xxxxxxxxx>:
> Benoit,
>
>>I see. So we can always provide general selfadjoint LDLt "at the
>>user's own risk". The generic case is that it will "work" since
>>division by exactly 0 would be incredibly bad luck, but it will be
>>unstable because one will be dividing by small numbers.
>
> I wasn't proposing continuing the factorization if a diagonal term
> is (or becomes) very small.
Ah , OK.
> I just propose that the test be
> if(abs(diag) < eps) return. This should insure that the algorithm
> doesn't compute a bad result, shouldn't it?
Yes, it should. But if the user has to test that the factorization
succeeded and in the converse case use another factorization, then I
am afraid that he will prefer doing the other factorization right away
! Especially if that other choice is going to be more stable and not
much slower (i'm thinking partial-pivoting block LU -- we could also
implement block Cholesky but that's yet more work and we don't have
enough arms).
> My own experience
> is that this will succeed often enough to make it useful.
Indeed, I was saying that's the generic case.
>> >But the problem is that pivoting won't help here if all the remaining
>>
>> >diagonal coefficients are zero. Indeed, a meaningful pivoting LDLt is
>>
>> >of the form:
>
> Are you talking about a case where an LDLt factorization exists but the
> matrix is singular?
No, I was still thinking about this 2x2 example,
0 1
1 0
that matrix is nonsingular and selfadjoint, yet its diagonal is 0.
> To be honest, before starting this current thread, I didn't even realize
> such a case could
> exist. What do people use this factorization for-- calculating matrix rank?
The basic use case is solving a system of equations,
> I think most
> people are interested in the non-singular case but, of course, would like to
> know if the
> factorization algorithm fails.
The main issue is that in computer / floating-point terms, all
matrices are "non-singular" since the determinant is virtually never
exactly 0, the real issue is numerical stability when the matrix is
"close" to being singular.
>Anyway,
> routine dsytrf does an LDLt factorization for a symmetric but not
> necessarily positive definite
> matrix.
Our LDLt does the same,
> They *have* chosen to include diagonal pivoting in the
> implementation; obviously a lot
> more work.
Yes, that's exactly what our LDLt does (didn't know it was called
diagonal pivoting but that makes a lot of sense), see my previous
email, transforming a matrix X into PXPt has exactly the effect of
permuting along the diagonal direction.
Cheers,
Benoit