Re: [eigen] Efficient way to store symmetric matrix for use with Cholesky? |

[ Thread Index | Date Index | More lists.tuxfamily.org/eigen Archives ]

*To*: eigen@xxxxxxxxxxxxxxxxxxx*Subject*: Re: [eigen] Efficient way to store symmetric matrix for use with Cholesky?*From*: Gael Guennebaud <gael.guennebaud@xxxxxxxxx>*Date*: Wed, 9 Jun 2010 14:00:21 +0200*Dkim-signature*: v=1; a=rsa-sha256; c=relaxed/relaxed; d=gmail.com; s=gamma; h=domainkey-signature:received:mime-version:received:in-reply-to :references:from:date:message-id:subject:to:content-type; bh=JCMketq8BNqBfSODtay0YxX8u6+YHj5hfNQhKlTcNw4=; b=MIC5uCGySaYnyhc34YRCBA2FmrbZW+SiCTo4lgGUsS9O2vRHw4DTHG/7XcEO15pf1+ JN7lnHUpmZTbyXWh0FCKEIwe1MLUMHWkzKno1TKslYLp1EiLWleKh+yRUaa+y15/xfp2 pNRV83LdUlNMloB+gIMZzaK5twqFIq+EhSF8E=*Domainkey-signature*: a=rsa-sha1; c=nofws; d=gmail.com; s=gamma; h=mime-version:in-reply-to:references:from:date:message-id:subject:to :content-type; b=YHg4iUT0eKCvzxrlJ2NsOQ8QJtkVMlAtdKHJKqrIx09dPwz1CYsIbA3TRlcskw+zUh f/q4hjMISAkytRaaioXsQ5LdS6wjieUc2MF/DXyeuah/93KmBkzr/ESa7rT3W+fQiazx SPjY4uRhxU2UrmiqZxUX4Y7SXgYkRSgLx2HFM=

arf, sorry you are right, my unit test was stupid enough to not catch that issue.

Anyway, now that's fixed. Thank you for the report!

gael

On Wed, Jun 9, 2010 at 1:00 PM, SHIROKOBROD Oleg <Oleg.Shirokobrod@xxxxxxxxxx> wrote:

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

I.

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);

I.

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 OfGael 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 +++

___________________________________________________________________________

This e-mail is confidential and is for the addressee only. Please refer to

www.oxinst.com/email-statement for regulatory information.

**References**:**[eigen] Efficient way to store symmetric matrix for use with Cholesky?***From:*Manoj Rajagopalan

**Re: [eigen] Efficient way to store symmetric matrix for use with Cholesky?***From:*Gael Guennebaud

**RE: [eigen] Efficient way to store symmetric matrix for use with Cholesky?***From:*SHIROKOBROD Oleg

**Messages sorted by:**[ date | thread ]- Prev by Date:
**[eigen] beta1 on june 26** - Next by Date:
**[eigen] efficient SoA with Eigen** - Previous by thread:
**RE: [eigen] Efficient way to store symmetric matrix for use with Cholesky?** - Next by thread:
**[eigen] stable norm**

Mail converted by MHonArc 2.6.19+ | http://listengine.tuxfamily.org/ |