Re: [eigen] Another LDLt issue |
[ Thread Index |
Date Index
| More lists.tuxfamily.org/eigen Archives
]
- To: eigen@xxxxxxxxxxxxxxxxxxx
- Subject: Re: [eigen] Another LDLt issue
- From: Benoit Jacob <jacob.benoit.1@xxxxxxxxx>
- Date: Mon, 30 Mar 2009 20:15:04 +0200
- Dkim-signature: v=1; a=rsa-sha256; c=relaxed/relaxed; d=gmail.com; s=gamma; h=domainkey-signature:mime-version:received:in-reply-to:references :date:message-id:subject:from:to:content-type :content-transfer-encoding; bh=JNo0QY1aNHfu6ueahoKCUQ+DjandzV3emrXFnMSte6k=; b=rNyQbGQ1Qausgip3xcEN0EE07ncmvQmAYrpOfijzCd8hEuFjrkiaSu6H3t/79VBAcu JEngH7/M9XYwvf8YV8dOlMbgtbTzAyk1RUMw123hZXukil4CQG/bxY5+4KA7iLVfeYys u8PnoxMdPKjhB0IrIEXEYcdnBbfQrPZaDVGBk=
- Domainkey-signature: a=rsa-sha1; c=nofws; d=gmail.com; s=gamma; h=mime-version:in-reply-to:references:date:message-id:subject:from:to :content-type:content-transfer-encoding; b=QwEggsvvBWnfpqlub6JbaU5r21HjOWgjTzxJLwb/2GwgAzNDTxd/8V0PRrIIGxDOyJ n8+XdT8WX0guMfUXDRXuOzq9gfSvspiMzrPXMSH/n/uySlAVIsbfE5W9Nxq529mZ3TYK LCEgU52nGc8QZ2hFsUr0FC+SFmNo1bAH1wzqo=
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