Re: [eigen] [RFC] Full pivoting LDLT

[ Thread Index | Date Index | More Archives ]

On Wed, Feb 4, 2009 at 10:10 AM, Jitse Niesen <jitse@xxxxxxxxxxxxxxxxx> wrote:
> On Tue, 3 Feb 2009, Keir Mierle wrote:
>> There are a couple of issues.
>> 1. Pivoting is not optional. Any objections?
>> 2. If the matrix is only semidefinite, then matrixL() gives wrong results
>> (1's in diagonal where they shouldn't be).
>> Minus objections, I'll commit this.
> I don't understand issue 2. I thought the L in the LDL^T decomposition has
> per definition 1s on the diagonal?
> It looks like you assume that the matrix is positive semi-definite. However,
> I believe that the LDL^T decomposition is also useful for symmetric
> indefinite matrices.

If A is symmetric indefinite, then a very different algorithm is
required because D is block-diagonal as 1x1 and 2x2 blocks. I did not
implement this because (for my use cases) A is always (semi)definite,
and it would substantially slow the code.

> Solving Ax = b with A symmetric is twice as fast with
> LDL^T than with LU (if I remember correctly).
> Why do you store P in two places: m_p and m_transpositions? Is using
> m_transpositions faster than m_p?

I left that in for the non-solve-in-place case. If you are solving in
place, you need to apply the transpositions via m_transpositions. If
you are solving into a new matrix, then m_p should be used instead
(see LU.h's use of m_p in solve()). I haven't implemented solve()
using m_p yet though (feel free!). m_p is a direct mapping from one
permutation to another; transpositions is the series of swaps to apply
to get to the other permutation.

> I was a bit surprised that you use both the triangles in m_matrix during the
> computation, since one half is just a rescaling of the other half. So I
> tried using only the lower triangle, but my timings (using g++ 4.1) are
> inconclusive on whether this is actually faster. I attach the diff for your
> entertainment, but I don't have any good reason for it.

Hmm, that is interesting that it's slower. I'll look at this tonight.


> Cheers,
> Jitse

Mail converted by MHonArc 2.6.19+