[eigen] Re: Special Matrices: first get diagonal matrices done

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

Hi,

Finally, this is done, in this fork:

http://bitbucket.org/bjacob/eigen2_special_matrices/

I imported permissions from eigen/eigen2 so you can commit directly to it.

The main reasons why I don't push directly to eigen/eigen2 are that
it's quite invasive, you may have comments, and also I was too lazy to
update the Sparse module accordingly, so I hope that Gael can help me
here; what needs to be done seems mainly to update
SparseDiagonalProduct.h, in the same way as the file
Core/DiagonalProduct.h: decouple it completely from the other products
by implementing a separate operator* for DiagonalBase.

So let's explain the new class hierarchy. I did some renaming:
DiagonalMatrixBase--->DiagonalBase
DiagonalMatrixWrapper--->DiagonalWrapper


DiagonalBase is still the common parent of DiagonalMatrix and of
DiagonalWrapper.

What's new is 1) DiagonalBase doesn't inherit MatrixBase anymore, and
2) there are new "abstract" base classes encompassing both
DiagonalBase and MatrixBase. These new classes are MultiplierBase and
AnyMatrixBase. See file MatrixBase.h.

A AnyMatrixBase is anything that can be copied into a MatrixBase.

A MultiplierBase is a AnyMatrixBase that furthermore can be multiplied
with a MatrixBase, producing a new MatrixBase.


AnyMatrixBase and MultiplierBase just serve to overload operator= ,
operator*, etc. The only method in them is derived().

DiagonalBase() is very minimalistic, for example I didn't implement
rows(), cols() etc. Instead, it has a method diagonal() that gives the
actual vector of diagonal coefficients. The idea was that at best we
would implement a very limited subset of the MatrixBase API in
DiagonalBase. That would force the user to look up the docs frequently
for availability of methods in DiagonalBase. So instead, I thought
that it would be easiest to force the user to address explicitly the
vector of diagonal coeffs, diagonal(), and from that point on the user
is in known territory with an actual MatrixBase. That said we can
always add rows() etc if you think that's useful.

The next steps for me are to code BandedMatrix (should now be a piece
of cake) and start the SVD...
Gael--->if you can take care of adapting the Sparse module, that would
be wonderful

on another note, I've been thinking too that the whole marked() stuff
for triangular and selfadjoint matrices could go away, in fact I don't
believe anymore that triangularness should be part of the Flags.
Instead it would be simpler and IMO just as good to have separate
methods like selfAdjointProduct().

Cheers,
Benoit

2009/5/11 Benoit Jacob <jacob.benoit.1@xxxxxxxxx>:
> Hi,
>
> Before getting started with BandedMatrix and other kinds of special
> matrices, I think the perfect starting point is to first nail the
> situation with diagonal matrices.
>
> Currently we have:
> * class Diagonal nests a matrix xpr and is a vector xpr
> * class DiagonalMatrixWrapper nests a vector xpr and is a matrix xpr
> * class DiagonalMatrix stores a vector and is a matrix xpr
>
> I think that we agree that we want instead DiagonalMatrix to NOT be a
> matrix xpr (ie NOT inherit MatrixBase) and have a coeff() method that
> just asserts row==col.
>
> DiagonalMatrixWrapper is still as much needed as before for situations like
>
> mat * vec.asDiagonal();
>
> where asDiagonal() is actually just being used as a thin wrapper
> around the vector allowing to select the matrix product strategy.
>
> What is open for discussion is whether DiagonalMatrixWrapper should
> still inherit MatrixBase? And if not, should there be an AsDense kind
> of class constructing a MatrixBase expression out of a
> DiagonalMatrixWrapper without evaluating?
>
> Let me explain why I think that the answer to both questions is _no_.
> The solution i'm going to advocate will be to keep the same class
> names, only make DiagonalMatrixBase no longer inherit MatrixBase.
>
> First of all we could have a DiagonalMatrixBase::toMatrix() that
> returns a plain matrix (see at the end for naming discussion), so that
> the user who wants something that behaves like a matrix, can get a
> plain matrix if he wants to:
>
> cout << vec.asDiagonal().toMatrix() << endl;
>
> The argument for having instead a asDense() returning a MatrixBase, is
> to avoid the cost of evaluating the matrix, but then remember that
> this expression would have a very expensive coeff() method as it would
> have to do a if(row==col). So as soon as you're accessing a large part
> of the matrix, it's not good anyway (performance wise) to have an
> asDense expression here.
>
> I also think that there isn't much of a use case for doing directly
> "vec.asDiagonal()(i,j)" -- it's OK, i think, to force the user to
> replace that by "i==j?vec(i):0" making him realize where there is room
> for optimization.
>
> Then let's look at the main use case for asDiagonal():
>
> mat * vec.asDiagonal()
>
> Currently, since DiagonalMatrixWrapper is a MatrixBase, we have to
> catch that in Product.h, adding complexity in there.
>
> If instead DiagonalMatrixWrapper were NOT a MatrixBase, we could catch
> that in a completely separate operator*. Less code overall and
> especially in Product.h, shorter compilations times, "you don't pay
> for what you don't use", same result.
>
> Is there anything I didn't think of?
>
> Yes that means that one could no longer do:
> mat * mat2.marked<DiagonalBits>()
> and hope for that product to be optimized, but frankly !! if you know
> that a matrix happens to be diagonal, then you can probably avoid
> storing it as a full matrix in the first place! So I don't think
> there's a real use case for that; and even if there were one could
> still do:
> mat * mat2.diagonal().asDiagonal().
>
> To summarize my plan:
> *no more changes to Diagonal (the ex-DiagonalCoeffs)
> *make DiagonalMatrixBase no longer inherit MatrixBase. It doesn't even
> need a coeff() method.
> *DiagonalMatrix still inherits DiagonalMatrixBase
> *DiagonalMatrixWrapper still inherits DiagonalMatrixBase
> *DiagonalMatrixBase is the class against which we implement the
> operator* catching diagonal products, and by the way we can also have
> operator= etc.
> *DiagonalMatrixBase::toMatrix() returns a plain Matrix.
>
> (the last method toMatrix(), i couldn't find a nicer name: eval()
> suggests it returns a DiagonalMatrix, toDense() doesn't make it clear
> it evaluates, ...)
>
> Notice that if you are really in favor of a diagonal matrix
> expression, we can still add this to this scheme using multiple
> inheritance on MatrixBase and DiagonalMatrixBase...
>
> Cheers,
> Benoit
>