Re: [eigen] patch for tridiagonal matrices in eigen

["Mauro Iazzi" <mauro.iazzi@xxxxxxxxx>]

Hi, sorry for replying so late.
I agree with most of Benoit' and Gael's ideas, but I have some
comments relative to my use case.
Also, I'm starting to understand this issue better, thanks to gael's
explanation, but I may have missed something. Tell me if this is the
case.

2008/12/30 Benoit Jacob <jacob.benoit.1@xxxxxxxxx>:
>
> 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.
>

It's a perfect criterion. However, when something else than an
expression is used, memory may be wasted. This is quite a concern with
me: a single matrix storage may be the few hundred megs that exceed my
memory and start using swap. I have no obvious solutions.

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

I think this is a good idea, since someone could implement a custom
storage class in this self contained way.
However, product by a diagonal matrix can be optimized for both left
and right multiplication. Moreover, I would not rely on template
specialization to choose the right path (like template<...>
ReturnType<> operator+(const MatrixType &a, const DiagonalType &b); )
but this is because I find this things quite unpredictable.

As Gael, I think it could still derive from MatrixBase. It seems to me
that a lot of things should be reimplemented in the opposite case.

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

I have to add that the sum shows similar problems. and while the
product expression between a diagonal and a square matrix can be
written without if()s, the sum expression cannot (unless I am missing
some obvious way). So also just for diagonals the problem remains.

In general return expressions should follow the general rule for Bad
Ideas you stated before. I think it's worth listing a few cases,
assuming there is a DiagonalType and TridiagonalType

DiagonalType D;
TridiagonalType T;
MatrixType A;

A = D + A;
A = T * A;
A = D + T;

all of those can be optimized if written directly. If one uses an
expression there are if()s and if one does not it creates an
unnecessary temporary (which can be lethal in my case). Last case has
the strongest possible optimization. Things like

A = (D + A) * A;

make this problem even more so. Seems to me that the problem is one
cannot optimize unless he knows where to put the result (e.g. when
there is an assigment).

One possible solution would be to allow the class to specify if it can
do better than the normal algorithm. Or, even better to use a cost
system like the one that is already in place for expressions vs.
cache. So the return type of (D+A) could override the assignment to
the fast version. On multiplication it would create a temporary as
normal. (D+T) could return a tridiagonal expression etc...

I guess this is just a more complicated version of what Gael suggests.
So even slower compilation times, since a lot of paths would need to
be specified. I can't really say this is a good solution. But it seems
one that saves memory.

A less complicated and less user-friendly method would be to just use
expressions and allow for specialization of eval() (and if possible of
assignment). In this way one can choose what he wants to optimize. A
relatively sane default for these expressions' flags can make this
solution quite painless in the general case.

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

I'd be really happy about this. I recently had a sparse matrix defined
in a way that makes it very easy to compute the product without even
saving the coefficients. If it was possible to define an expression
for it cheaply it would be a big gain.

I hope this made sense. Cheers,

mauro

---