[eigen] discussion on special matrices

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

Hi,

After a long break away from coding I am trying to get back to the
important discussion on special matrices. I know that's long overdue
and hope Mauro didn't lose all hope...

So Gael, I studied your page
http://eigen.tuxfamily.org/index.php?title=SpecialMatrix
and I think I agree with most of your ideas but sometimes I can't be
sure because it's hard to tell what you meant. For example you say
"remove Part" but then "move solveTriangular* to Part::solve* " so I
guess that you meant TriangularView::solve*.

Anyway, let me make a proposal following closely yours with a few
differences, to see if we agree. I use email for now but once we agree
on a solid base, we come back to the wiki.

*** CHANGES TO THE PART CLASS ***

A. replace if() by ei_assert() in all coefficient accessors in Part.
So, to schematize,

Part<UpperTriangular>::coeff(int i, int j) { return i<=j ? m_matrix(i,j) : 0; }

becomes

Part<UpperTriangular>::coeff(int i, int j) { ei_assert(i<=j);
m_matrix(i,j) : 0; }

In other words, Part becomes completely trivial (useless) from the
point of view of coefficient access. The advantage is that it does not
have a cost.

Notice that we already have this behavior in coeffRef().


B. Part no longer inherits MatrixBase

C. We keep the operator=

Part::operator=(const MatrixBase&)

that allows to do " m.part<UpperTriangular>() = other_xpr;"

Yes we'll discuss later the ugliness of the triangular/rectangular
mismatch in the API here... I don't see an obvious solution. Maybe we
can name it "copyFrom" instead of "operator=" ?

D. We add a few methods like setZero() as that is widely used, trivial
wrappers around C.

E. We add the triangular solving here as Part::solve().

F. If we implement compact triangular storage as a separate class
CompactTriangular, then it seems indeed to make sense to move some
methods (like solve) to a TriangularBase, however it's nontrivial,
because for high performance, the coefficient accesses need to be done
differently in a CompactTriangular matrix. Indeed, since every coeff
index by row and column takes 7 arithmetic operations, we don't want
to implement the solving with this kind of access! Without
vectorization, i'd be confident that we can factor coeff accesses
using pointers to rows/columns depending on storage orders. So that's
the kind of coefficient access that we should implement in a
TriangularBase.

But to say everything, between the difficulty of implementing
efficient triangular solving and the fact that each access requires 7
ops, I'm not convinced about CompactTriangular.

Another thing: if we make Part inherit TriangularBase, then we are
restricting Part to triangular parts, so that means more classes, and
some meta kungfu in MatrixBase::part() to determine the return type.

***** DIAGONAL MATRICES *****

I agree 100% with your proposal. The removal of the MatrixBase
inheritance is completely aligned with my chakras of simplifying the
whole scheme and making all these special matrix classes just simple,
dumb, specialized and efficient.

**** OTHERS ****

I assume that for the others too, you meant no MatrixBase inheritance
(that was my opinion when we first discussed it so I may be wishful
here, correct me if I am).

For Band matrices, we really need non-square ones. The most important
example is bidiagonalization, that is useful for the SVD. It should
work for all rectangular matrices (no reason at all to restrict to
square).

On the other hand, it's also very important to have the special case
of self-adjoint band matrices. For tridiagonal matrices, in most cases
they're selfadjoint (used for selfadjoint eigensolving) and then we
want to optimize the storage. Then of course they're square.


***** SECTION "How to glue everything together" ************

About the Common base class, I don't understand. You suggest to move
operator+ there ; first of all I was not convinced that such an
operator would be useful at all, I was more seeing these special
matrix classes as holding intermediate data to implement and/or call
solving/decomposition algorithms, and I don't see the use case for
performing a full addition of band matrices....? Then, assuming we
really want operator+, it would be implemented in a completely
different way between different special matrix classes, so could not
share a common implementation in a base class.

About providing specializations of assignment, i completely agree.

For the same reason as above, I wouldn't bother at all about CwiseBinaryOp.

Cheers,
Benoit