2009/3/30 Bill Greene <w.h.greene@xxxxxxxxx>:
> Benoit,Ah , OK.
>>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.
Yes, it should. But if the user has to test that the factorization
> 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?
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
Indeed, I was saying that's the generic case.
> My own experience
> is that this will succeed often enough to make it useful.
No, I was still thinking about this 2x2 example,
>> >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?
that matrix is nonsingular and selfadjoint, yet its diagonal is 0.
The basic use case is solving a system of equations,
> 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 main issue is that in computer / floating-point terms, all
> I think most
> people are interested in the non-singular case but, of course, would like to
> know if the
> factorization algorithm fails.
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.
Our LDLt does the same,
> routine dsytrf does an LDLt factorization for a symmetric but not
> necessarily positive definite
Yes, that's exactly what our LDLt does (didn't know it was called
> They *have* chosen to include diagonal pivoting in the
> implementation; obviously a lot
> more work.
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.