Re: [eigen] [RFC] Full pivoting LDLT

[Keir Mierle <mierle@xxxxxxxxx>]

On Thu, Feb 5, 2009 at 10:31 AM, Jitse Niesen <jitse@xxxxxxxxxxxxxxxxx> wrote:
> On Wed, 4 Feb 2009, Keir Mierle wrote:
>
>> If A is symmetric indefinite, then a very different algorithm is
>> required because D is block-diagonal as 1x1 and 2x2 blocks.
>
> I don't understand why D would have 2x2 blocks. Could you please explain
> this, or tell me where you read this?
>
>> 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.
>
> Fair enough, thank you for your explanation.
>
> I'm not sure transpositions are the fastest way to do this.

Neither am I :)

I did this because I didn't want to make a sizeof(b.rows()) copy when
solving in place. It is probably faster to do that in some cases
though. I'm open to measurements. It's an easy change to make.

Actually, this brings up something else I was thinking of: We should
make efficient permutation matrix class, which has a specialization
for matrix multiply which permutes the matrix. This way we could
return a meaningful matrixP() from the pivoting decompositions, and it
would remove code duplication (both LU and LDLT do pivoting now, so
will QR when I get a chance).

> For instance,
> suppose you have a 3x3 matrix, and your pivots are (1,1), (2,2) and (0,0).
> Thus m_transpositions is (1,2,2) so you permute the vector b on the
> right-hand side by doing two swaps, thus six reads and writes. However, it's
> possible to do it in four:
> tmp = b(0); b(0) = b(1); b(1) = b(2); b(2) = tmp;
>
> On the other hand, the triangular solve is O(n^2) while the permutation is
> O(n) so it probably doesn't make sense to optimize the permutation.

Again, for the in-place solve I didn't want to allocate a O(n)
temporary. However, it is probably worthwhile. We should measure this.

Keir

>
> Cheers,
> Jitse
>
>
>