[eigen] SimplicialCholesky - subclass or patch?

[Douglas Bates <bates@xxxxxxxxxxxxx>]

I forgot to cc: the list on this reply.


---------- Forwarded message ----------
From: Douglas Bates <bates@xxxxxxxxxxxxx>
Date: Tue, Jun 28, 2011 at 5:19 PM
Subject: Re: [eigen] SimplicialCholesky - subclass or patch?
To: Tom Makin <t.makin@xxxxxxxxxx>
Cc: gael.guennebaud@xxxxxxxxx


On Tue, Jun 28, 2011 at 2:59 PM, Tom Makin <t.makin@xxxxxxxxxx> wrote:
> I would very much support any additions to Eigen that allow direct access to
> the L matrix in a sparse Cholesky decomposition. I am currently working on a
> surrogate modelling application and have recently switched over to Eigen in
> all my code. The main reason I need access to L is to update the
> factorisation without having to re-factorise the whole thing when I add a
> couple of extra rows. I am comfortable doing this for the dense case but I
> am yet to get to grips with the sparse version because I can't get at L
> directly. I have been able to do a limited amount of sparse updating using
> cholmod, but I would like to try and implement it in way that is independent
> of the sparse solver.
>
> Douglas: Would you be willing to send me your SimplicalCholesky subclass?

I spoke too soon.  I had not yet tested the code that I wrote and,
naturally, it blew up on me.

The problem I am encountering is that the return type of a matrixL
method will be determined by the MatrixType and the current value of
m_LDLt. If m_LDLt is true the return type should be
TriangularView<MatrixType, UnitLower>, otherwise
TriangularView<MatrixType, Lower>.  Because m_LDLt can be set or unset
for an instance of the class, I don't think there is a way to
determine at compile time which it should be.  The only way I can
think of to give a definite return type is always to return the Lower
form and, when LDLt is true, add the identity to the strict lower
triangle, which wastes some storage, and will have less efficient
solve methods than are otherwise possible.

Am I missing something here?


> On 28/06/2011 19:00, Douglas Bates wrote:
>>
>> I realize that the Sparse module is not stable in Eigen 3.0.1 and the
>> SimplicialCholesky class is in the unsupported directory.
>> Nevertheless, these are exactly what I need for a central calculation
>> in my research.  (If you want to know the details I include a draft of
>> a paper describing the particular calculations being done.  You can
>> skip to sections 2.3, 2.4 and 2.6 for the relevant steps).
>>
>> Current versions of our code are based on SparseSuite with a rather
>> clunky interface from R.  I would much prefer to use Eigen but, as
>> shown in the enclosed description, I need to access the triangular
>> matrix L from the LLt decomposition and I need a method like
>> logAbsDeterminant for the SimplicialCholesky classes.  I can do this
>> by subclassing the SimplicialCholesky templated class and have done
>> so.  A better long term solution would be to produce and submit
>> patches for SimplicialCholesky.h but before doing so I would like to
>> check that it would be worthwhile.  If that code is not considered
>> viable in the long term then I will stick with subclassing rather than
>> modifying that class.  Could someone from the core group indicate if
>> you would want to receive patches for that file?  I can't guarantee
>> that I will get all the myriad details of the Traits right the first
>> time around but I am willing to learn.
>
>