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

[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).

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