Re: [eigen] Proposal: Matrices requiring compact storage of triangular part

[Benoit Jacob <jacob.benoit.1@xxxxxxxxx>]

2010/7/19 Gael Guennebaud <gael.guennebaud@xxxxxxxxx>:
> Hi Manoj,
>
> thanks for this detailed email, but as Benoit said this is a bit too
> much for now, and most of it can be easily implemented without any
> modifications in Eigen. Actually, I think there are only two welcome
> changes that are affecting Core:
> - rename the current TriangularBase to TriangularViewBase
> - add support for symmetric complex matrix (not only selfadjoint).

OK for both points.

Benoit

>
> For the rest, this can be easily implemented in a separated module
> providing compact triangular representations and related algorithms.
> Also, regarding the design, I'd rather see something like:
>
> CompactTriangularBase
>    <- CompactTriangularMatrix
>    <- CompactSelfAdjointMatrix
>    <- CompactSymmetricMatrix
>
> rather than a selfadjoint/symmetric view on a CompactTriangularMatrix.
>
> Also, since CompactSelfAdjointMatrix and CompactSymmetricMatrix are
> very very similar, they could either be factorized in a unique class
> with a template option (but I cannot find a good name for either
> symmetric or selfadjoint), or they could inherit a common base class
> factoring most of the code.
>
> LAPACK original triangular layout has zero interest IMO, so let's
> forget about it.
>
> that was just some quick thoughts ;)
>
> gael
>
>
>
> On Wed, Jul 14, 2010 at 11:00 PM, Manoj Rajagopalan <rmanoj@xxxxxxxxx> wrote:
>> Hi Eigen developers,
>>
>>   Gael's recent reference to the RFPF storage layout for triangular (portions
>> of) matrices got me experimenting with implementing a SymmetricMatrix class
>> and related products. The question of sharing the code with Eigen came up.
>> Here are my thoughts on slight changes to Eigen3 that might be required
>> before the integration can be performed.
>>
>> 1. Define a Triangular_RFPF class (for lack of a better name). I don't know if
>> this can be folded into the existing TriangularBase hierarchy (which assumes
>> rectangular/dense 'models' being viewed). If not, then I suggest renaming the
>> existing TriangularBase to TriangularViewBase and the introducing a
>> TriangularBase as the CRTP base-class for all types of triangle-based views.
>>
>> 2. Define a selfadjointView() method that returns a hermitian view of the
>> stored triangle. Since the Triangular_RFPF class itself will take a Lower|
>> Upper mode indicator, this method won't be templated like it is for
>> MatrixBase<> derivatives.
>>
>> 3. Also define a symmetricView() method that returns a symmetric view of the
>> stored triangle.
>>
>> 4. The two views will be identical for real Scalar types. For complex the two
>> are mutually-exclusive. So we will need to introduce a Symmetric=0x20 into
>> the mode enum in util/Constants.h. For real Scalar cases, Flags &
>> (SelfAdjoint|Symmetric) should be true. In other words, both bits should be
>> set. For the complex case only one should.
>>
>> 5. Some kind of ei_meta_if<> will need to be used to return identical wrapper
>> classes for symmetricView() and selfadjointView() if Scalar is real, so that
>> existing product implementations etc. can work as is. Then, only symmetric
>> products will have to be implemented for complex Scalars.
>>
>> 6. New product wrappers will have to be written when Triangular_RFPF or its
>> (selfadjoint|symmetric)View() are involved but this will involve working with
>> blocks of the RFPF storage matrix as discussed in the original paper so this
>> part should be easy.
>>
>> 7. New triangular solver wrappers will have to be written but this should
>> again be easy - the paper discusses the lower-triangular case.
>>
>> 8. The SelfAdjointView<> wrapper serves as a validator for LLT and LDLT. So
>> symmetricView() for real Scalars should have a SymmetricViewReturnType
>> typedef that maps to SelfAdjointView<> for real scalars. This will ensure
>> minimal changes to Eigen3. Since there are no specializations of L(D)LT
>> implementations for classes other than SelfAdjointView, we are good.
>>
>> 9. Since TriangularView<>, Triangular_RFPF<> and a hypothetical vanilla
>> TriangularMatrix class (with LAPACK layout) share a common abstraction, this
>> begs the question of whether they should inherit from siblings at some level
>> since they differ either in ownership of the data or in the layout but are
>> all triangle-accessors. In fact, layout and ownership are orthogonal (you can
>> have a view with RFPF layout instead of the vanilla LAPACK layout), there is
>> scope for more generalization but this might be too much work.
>>
>> 10. Complex symmetric matrices arise in certain differential equations
>> situations and there are efficient solvers for certain classes that use
>> complex orthogonal (as opposed to unitary) transformations. Having complex
>> symmetric classes would be desirable at some point.
>>
>> 11. I think the complex case can be safely postponed to a later release of
>> Eigen depending on use-community's interest. Then SymmetricProduct wrappers
>> can be deferred. Common cases will work with the SelfAdjointProduct wrappers
>> and Eigen releases won't be stalled :-)
>>
>>
>> Thanks,
>> Manoj
>>
>>
>>
>
>
>