RE: [eigen] Efficient way to store symmetric matrix for use with Cholesk

Hi Gael,

It seems that there is still a bug in LDLT. I discussed this case with Benoit

>Hi,

>a quick look at the code suggests that SelfAdjoitView::ldlt() is unimplemented, unless i missed it. So either we should implement it, or remove the declaration.

>Right now, you can do: matrix.ldlt() but not selfadjointview.ldlt().

>I guess the real problem is that we used to have a little bug in ldlt that makes it use more than just one triangular half of the matrix.

>See comment in LDLT.h. Maybe the first thing to do, is to fix this bug.

>Benoit

Here is a test case

typedef Eigen::Matrix<double, Eigen::Dynamic, 1> Vector;

typedef Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic> Matrix;

Matrix M(3, 3);

Vector values(3);

// Fill the Vector velues

vlues =

1.0

3000.0

9000000.0

M.selfadjointView<Eigen::Upper>().rankUpdate(values);

40.0 63000.0 129150000.0

0.0 129150000.0 297675000000.0

0.0 0.0 731699325000000.0

II.

M += values * values.transpose()

40.0 63000.0 129150000.0

63000.0 129150000.0 297675000000.0

129150000.0 297675000000.0 731699325000000.0

Eigen::LDLT<Matrix, Eigen::Upper> ldltOfM(M);

Vector rightHandSide(3);

// Fill the Vector rightHandSide

rightHandSide =

45.73845

93646.978500000027

215689932.67500001

ldltOfM.solveInPlace(rightHandSide);

result rightHandSide

1.1434612500000001

0.00072510242740998852

2.9477946105116333e-007

II.

result rightHandSide

0.00067356578946468231

0.00073205916040101719

-3.1608187134546198e-009

As you see in the first case strictly lower triangular part of M is zero apart the second case with full selfajoint matrix filled.

Correct result is for the second case. So LDLT still uses by the bug lower triangular part of a matrix.

Thanks,

Oleg

From: Listengine [mailto:listengine@xxxxxxxxxxxxxxxxx] On Behalf Of Gael Guennebaud
Sent: Tuesday, June 08, 2010 22:26
To: eigen@xxxxxxxxxxxxxxxxxxx
Subject: Re: [eigen] Efficient way to store symmetric matrix for use with Cholesky?

On Tue, Jun 8, 2010 at 7:02 PM, Manoj Rajagopalan <rmanoj@xxxxxxxxx> wrote:

Hi,

Right now I am using a full square matrix to store a symmetric matrix and
am calling LDLT on it. I believe a bug was recently fixed that gets the
Cholesky module accessing only one triangular part of the matrix, as it
should. Is there some way to just store one triangular part and call LDLT on
that instead of storing the entire matrix?

Hi,

indeed, now the LDLT dec works as the LLT in the sense you can specify which triangular part has to considered:

mat.selfadjointView<Upper>().ldlt()...

or

LDLT<MatrixType,Upper> ldlt(mat);
.....

Of course you still need a full dense matrix, but note that one half is free to store another triangular or selfadjoint matrix.

Now it seems you are looking for a compact storage for triangular matrices. This is currently not available in Eigen but if someone is willing to do it the idea is to split the triangular matrix into two triangles and one square and transpose one of the triangle such that it can be put in the remaining empty space as it is better explained there: http://netlib.org/lapack/lapack-3.2.html#_9_6_rectangular_full_packed_formate The main advantage of this scheme is that's quite easy to implement matrix products, solvers and decompositions on top of it reusing dense algorithms. Since the benefits are relatively small, this is a low priority for us. Nevertheless, I think this is a quite nice project for someone who want to contribute because it mainly consists in relatively high level coding.

cheers,

gael

Thanks,
Manoj

+++ Virus-scanned by MailControl for Oxford Instruments +++