Re: [eigen] about band/skyline matrices

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

oook now i've thought a bit more. i understand now the importance of
addressing rows and cols, sorry :) wasn't fully awake.

About BandedBase, i'm also not sure anymore of what i said because the
efficient implementation of product and of assignment depends so much
on the storage that not much can be factored.

So yes now i mostly agree with your proposal,

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

with only one exception: I think that the default should always make
the diagonals be contiguously stored. So the default should always be
row-major, even if that's contrary to the rest of Eigen that defaults
to column-major. Although that's surprising (which is bad) it's not
worse than to have the default order change between BandMatrix and
Tridiagonal; and it's better for performance.

The alternative would be to default to col-major but do the
"transposed" storage scheme, but then we need to document better how
BLAS interoperability works

Cheers,
Benoit

2009/7/15 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
>