Re: [eigen] patch for tridiagonal matrices in eigen

["Gael Guennebaud" <gael.guennebaud@xxxxxxxxx>]

Hi Benoit,

I mostly agree with your proposal. Two things:

1) I would suggest the following declaration:

template<typename Scalar, int Size, bool AsXpr = false> DiagonalMatrix;

the main reason: this is easier to use as a "diagonal matrix with
storage" class, and this does not prevent to use it as an expression
via the additional template parameter. Moreover, such an interface is
more suited for a TridiagonalMatrix class or for any other matrices
with skyline storage. hm... this last remark let though that all these
special matrices could be handled by a unique and generic
SkylineMatrix class: template<typename Scalar, int Size, int Bands,
bool AsXpr> SkylineMatrix; In practice, DiagonalMatrix could inherit
SkylineMatrix<Scalar,Size,1>....

2) In order to guarantee a high consistency, I'm not 100% sure that
DiagonalMatrix should not inherit MatrixBase. It is very important
than any function which takes as input an Eigen matrix expression, can
also takes as input a DiagonalMatrix. Also, inheriting MatrixBase
could avoid a couple of code duplication. For instance, trace() is
simply .diagonal().sum(), so once we have overloaded both diagonal()
(to return a vector) no need to reimplement neither sum() nor trace...
 However this would require to make sure that all calls to a
MatrixBase<> object are prefixed by ".derived()" that is currently not
the case. A robust and general solution to that problem is to go back
to pairs of publicfunctionname() _implementationfunctionname() where
the former would be implemented only in MatrixBase as:
publicfunctionname() { return derived()._implementationfunctionname(); }
and MatrixBase() would provide a default implementation of
_implementationfunctionname(); Perhaps, one or two macro would be good
here. As usual, the drawback is compilation time.

Cheers,
Gael.

On Tue, Dec 30, 2008 at 12:23 PM, Benoit Jacob <jacob.benoit.1@xxxxxxxxx> wrote:
> Hi Mauro and Gael,
>
> Finally I started actually thinking about this discussion (until now
> my attitude was "noooo resist entering into yet another technical
> discussion you don't have time")
>
> First let's repeat one comment about the DiagonalMatrix class: Its
> coeff() method does an if(). So it's slow (except in the fixed-size
> cases where the if() can go away when loops are unrolled, but that's
> not the general case). In 99% of practical use cases, one doesn't want
> and one doesn't need to pay that price. For example, the diagonal
> product is rewritten in such a way that this coeff() method is
> bypassed.
>
> That makes me think that the DiagonalMatrix expression was a bad idea
> (my bad, it dates back to early eigen2 days) in the first place. Right
> now I think it's completely useless, read on to see what I propose to
> have instead.
>
> First let's look at what happened since this DiagonalMatrix class was written:
> 1) Gael wrote the Diagonal product code that actually completely
> bypasses the DiagonalMatrix class (i.e. the coeff() method). In Gael's
> DiagonalProduct, the DiagonalMatrix is just used as a trivial wrapper
> around the diagonal vector (which is directly used).
> 2) Mauro is now puzzled as to how to implement TridiagonalMatrix
> because the wrong design of DiagonalMatrix is even worse when applied
> to tridiagonal matrices.
>
> Before continuing, let's state a rule to determine when an expression
> class (here, DiagonalMatrix) is a bad idea (of course it's a
> sufficient condition, not a necessary condition):
>
> An expression class MyExpression is a bad idea whenever the
> MyExpression::coeff() method is slow and in practice it can be
> bypassed when implementing actual algorithms.
>
> For exemple, as we said above, in the case of DiagonalMatrix, the slow
> coeff() method is bypassed when implementing the diagonal matrix
> product.
>
> What do I propose to do then?
>
> If MyExpression represents a special class of matrices for which the
> storage can be optimized and certain algorithms can be optimized, then
> instead of an expression class MyExpression, I think that what we want
> is a special matrix class MySpecialMatrix that does NOT inherit
> MatrixBase, so it's not at all an expression. MySpecialMatrix should
> have only the storage that it needs, and should be where we implement
> the special optimized algorithms.
>
> For example the DiagonalMatrix class, rethought as a "special matrix",
> would look like:
>
> template<typename DiagonalVector>
> class DiagonalMatrix
> {
>  public:
>    template<typename OtherDerived>
>    DiagonalProductReturnType<blabla> operator*(const
> MatrixBase<OtherDerived>& other);
>
>  protected:
>    DiagonalVector m_diagonal;
> };
>
> Here it remains to be seen what DiagonalProductReturnType<blabla> is.
> I think that we can make it a plain expression (i.e. inheriting
> MatrixBase) because there the coeff() method can be implemented
> without an if(). Not sure though, how that plays with vectorization
> (depending on the matrix storage order). For TridiagonalProduct
> though, it's not possible to make it return an expression (again
> because the coeff() would do an if()). Probably we want some meta
> logic here, like we have in Product.h, to distinguish cases where we
> want an expression and cases where we don't.
>
> You probably are going to ask: am I talking about a solution for the
> diagonal-matrix-expression or for the diagonal-matrix-storage? The
> answer is that the DiagonalMatrix above works for both! If you take
> DiagonalMatrix<VectorXd> you have a diagonal-matrix-storage. On the
> other hand, the diagonal-matrix-expression returned by
> MatrixBase<Derived>::asDiagonal() would be a DiagonalMatrix<const
> Derived&>. The user is not expected to construct that himself, he just
> calls asDiagonal().
>
> Another aspect of my proposal is that Product.h does no longer have to
> dispatch all the special products, they are now implemented in the
> special matrix classes. So it's more extensible and potentially
> lighter on the compiler.
>
> What do you think?
>
> Cheers,
> Benoit
>
> ---
>
>

---