|Re: [eigen] SimplicialCholesky - subclass or patch?|
[ Thread Index |
| More lists.tuxfamily.org/eigen Archives
- To: eigen <eigen@xxxxxxxxxxxxxxxxxxx>
- Subject: Re: [eigen] SimplicialCholesky - subclass or patch?
- From: Gael Guennebaud <gael.guennebaud@xxxxxxxxx>
- Date: Thu, 30 Jun 2011 08:32:16 +0200
- Dkim-signature: v=1; a=rsa-sha256; c=relaxed/relaxed; d=gmail.com; s=gamma; h=mime-version:in-reply-to:references:from:date:message-id:subject:to :content-type:content-transfer-encoding; bh=2MZAuCjpDCwwLUKP2aSIVGA7hSNgjTQby7tLi5nj6m8=; b=xTGkpGrPDL+0Z+MFUJojj7F2SnALXzjrpdU1mljare6SFwfIYJBEANsICTNgoUxxNz upFO4qWNwhQxUQ2/qxJnSwDDpO3lH4qKAZYVj3tx6NtSEHGKLIs3rGlmi+YDkazxKMbG QpG15jW/KPxzuuMVNGKfD2sFPFXzKsUHnZTgk=
hm, good point. That would suggest to have two classes, SimplicialLDLt
and SimplicialLLt that could be based on a common
On Wed, Jun 29, 2011 at 12:19 AM, Douglas Bates <bates@xxxxxxxxxxxxx> wrote:
> 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.