Re: [eigen] about band/skyline matrices

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

Hi,

it's great to have you thinking about these storage schemes. Naively I
would have coded it as n separate vectors. That's what I call below
the naive scheme.

On my side, my brain is still stuck in an infinite loop with the
design of MatrixBase wrt special matrices :/

2009/7/15 Gael Guennebaud <gael.guennebaud@xxxxxxxxx>:
>
> Hi all,
>
> I've started to implement a generic BandMatrix class with storage which aims
> to efficiently and compactly represents matrices where only a few sub- and
> super-diagonals near the main diagonal are non zero. See:
> http://bitbucket.org/bjacob/eigen2_special_matrices/src/tip/Eigen/src/Core/BandMatrix.h
>
> My question is about the storage scheme. Currently I've simply packed the
> diagonal in a single vector, e.g.:
>
> [11 12 13 00]
> [21 22 23 24]
> [00 32 33 34]
> [00 00 43 44]
>
> is stored:
>
> size = 4
> nb supers = 2
> nb subs = 1
> data = [11 22 33 44 12 23 34 13 24 21 32 43]
>
> The main advantages are:
>  - the diagonal elements are sequentially stored in memory
>  - minimal memory consumption

This isn't much different that the naive storage scheme: for fixed
size it amounts to the same; for large enough dynamic size (for a
given number of diagonals) it is almost the same.

>  - it is easy to extract any diagonals

But not as easy as with the naive scheme, right? You still need to do
an offset computation to know where in your vector the k-th diagonal
is stored.

> Drawbacks:
>  - it is not possible to extract a row or a column
>  - => it is not possible to implement, e.g., an efficient matrix product

I don't know if that's a big problem! Look at how you'd implement
bidiagonal*matrix product: you wouldn't need to address rows of the
bidiagonal matrix, right ?

It seemed to me that by far the most important thing is to be able to
efficiently address each diagonal.

>  - limited to square matrices

OK so that wouldn't be the unique solution, but we can always have a
BandedBase and have your solution as one of several derived classes,
like what we have with DiagonalBase.

> So another possibility would be to implement the BLAS storage scheme where
> the key idea is to "move" the columns up and down such that the diagonals
> become rows, and then you prune the empty rows. The previous matrix would be
> stored in a column major matrix like this:
>
> [__ __ 13 24]
> [__ 12 23 34]
> [11 22 33 44]
> [21 32 43 __]
>
>
> The main advantages are:
>  - it is easy to extract the non zero part of any column,
>  - the columns are sequentially stored in memory => products can be
> efficiently implemented
>  - it is still possible to extract any row using something similar that
> DiagonalCoeffs
>  - the diagonals are still easy to extract
>  - compatible with BLAS/LAPACK
> Drawbacks:
>  - a little waste of memory (negligible if the number of sub and supers
> diagonal is small compared to the size of the matrix)

safe assumption ;)

>  - the diagonals are not sequentially stored in memory that is not very good
> as a representation for tridiagonal matrices.
>
> Note that the last drawback can be overcame using a row major matrix. In
> that case the storage become compatible with the BLAS triadiagonal storage
> scheme but not with the band storage scheme which impose a colmajor matrix.

That's weird...!
Indeed it seems very important to me that the diagonals are sequentially stored.

>
> So here is what I suggest:
>
> - Implement BandMatrix using the BLAS/LAPACK storage scheme with the
> additional possiblity to use a rowmajor matrix for the underlying storage
> - Add a TridiagonalMatrix<Scalar,Size,Options> class inheriting with
> simplified template parameters and simplified ctor, e.g.:
>
> template<typename Scalar, int Size, int Supers, int Subs, int Options =
> ColMajor>
> class BandMatrix;
>
> Options = ColMajor or RowMajor | SelfAdjoint
>
> template<typename Scalar, int Size, int Options>
> class TridiagonalMatrix : public BandMatrix<Scalar,Size,1,1,Options |
> RowMajor>;
>
> Another possiblity would be to automatically fall back to a RowMajor storage
> if Supers<=1 and Subs<=1 and the storage order is left unspecified by the
> user, e.g:
>
> BandMatrix<float,Dynamic,1,1>               => creates a BLAS compatible
> tridiagonal matrix
> BandMatrix<float,Dynamic,1,1,ColMajor> => creates a BLAS compatible band
> matrix

I'm still not convinced that either storage scheme discussed here has
an advantage over the naive scheme that i described above. I still
don't see the use case for addressing the rows or cols of a banded
matrix.

But the fact that we're already discussing 3 different schemes, and
that BLAS interoperability might matter to people, suggests that
perhaps we should be flexible and support different schemes as you
suggest.

I think that the best way to be flexible is to have a BandedBase class
with the API to address the diagonals, and then implement the various
storage schemes as derived classes of it. Like what we have in
DiagonalBase.

This way, you can start implementing your favorite storage scheme and
if I like another one i can add it too. One could also have separate
classes for special cases, like Bidiagonal, TridiagonalSelfadjoint,
etc.

By the way, let me insist that in the tridiagonal matrices case,
what's most important is selfadjoint tridiagonal matrices. So i don't
see much use case for a general tridiagonalmatrix class, but of course
SelfadjointTridiagonal is very important.

I still don't have an opinion whether the selfadjoint variants should
be completely different classes or a template parameter as you
suggest.

Cheers,
Benoit