Re: [eigen] first benchmark of large and sparse matrix

["Gael Guennebaud" <gael.guennebaud@xxxxxxxxx>]

On Mon, Jun 23, 2008 at 7:38 PM, Schleimer, Ben <bensch128@xxxxxxxxx> wrote:
> Hi guys,
>   I didn't know gmm++ matrices were so slow either.

well, so far I've only tested  the sum which is not a critical
operation.  I guess they mainly focused on the optimization of the
matrix product and triangular solver which are the main bricks to
build linear solvers.

> I wonder if we could take the gmm++ algos and use templated eigen sparse matrices with them...

that was more or less the plan, but I'm not sure anymore because of
licensing issues....

>
> btw, how are you designing sparse templated matrices?
> Do you have any papers to point to?

well, nothing is fixed yet ! you might find some info on the wiki.  So
far,  I found that the most common storage scheme is the "Compressed
Column Storage" strategy:
http://www.netlib.org/linalg/html_templates/node90.html
so I started to play with it such that we are able to use optimized
backends for the high level algorithms. To manipulate such matrices,
the idea is to replace the for loops (random access) by Iterators: we
can specialize such Iterators for each type of expression, and we can
even nest the Iterators as we nest the expressions.
CwiseBinaryOp::InnerIterator is a good example of what I mean (in
src/Sparse/CoreIterators.h).

One of the main challenge with such a storage is to dynamically set a
matrix.  GMM++ does not allow any write operation, so you have to use
another kind of sparse matrix, i.e., a row of std::map and eventually
you do a copy to build your final and efficient sparse matrix. MTL4
uses a similar principle: they preallocate 5 entries (or more if you
want more) per column and they manage overflows using a std::map.
Currently I implemented something different which allows to directly
fill such a matrix under the constraint that the coefficients must be
set by increasing indices. I guess this is enough in most cases but we
definitely need something more flexible, at least to be able to
implement any kind of matrix product.

of course any idea and suggestion are more than welcome !

cheers,
gael


>
> Cheers
> Ben
>
>
> --- On Sun, 6/22/08, Benoît Jacob <jacob@xxxxxxxxxxxxxxx> wrote:
>
>> From: Benoît Jacob <jacob@xxxxxxxxxxxxxxx>
>> Subject: Re: [eigen] first benchmark of large and sparse matrix
>> To: eigen@xxxxxxxxxxxxxxxxxxx
>> Cc: "Gael Guennebaud" <gael.guennebaud@xxxxxxxxx>
>> Date: Sunday, June 22, 2008, 11:05 PM
>> Man, that is good...!!!
>>
>> It's normal that the sparse matrix is a bit slower than
>> the dense one with 10%
>> nonzero coeffs; also I suppose that the memory usage is
>> already significantly
>> lower. So your numbers look very good!
>>
>> I wasn't expecting GMM++ to perform so poorly; it
>> sounds to me as if their
>> sparse matrices are not intented for doing matrix
>> arithmetic, but instead are
>> just intended to apply solving algorithms; still, being
>> able to form sums is
>> at least useful to construct the matrices in the first
>> place! (And I
>> definitely personnally see use cases for full support for
>> arithmetic on
>> sparse matrices, as this models a mathematical object known
>> as a 'convolution
>> algebra' and is thus very interesting).
>>
>> So it is sounding like you're solving the problem of
>> sparse matrix storage; if
>> there exists a possibility of providing an abstraction
>> layer for sparse
>> algorithms supporting multiple backends (as discussed on
>> IRC, implementing
>> sparse algos ourselves is really too ambitious, so we want
>> to reuse existing
>> code for that) then we would get a very, very good
>> Eigen-Sparse module (or
>> auxiliary library, depending on how you see it).
>>
>> FYI this is of high interest for Krita and for Step and
>> probably has countless
>> interesting applications; a while ago Benjamin did play
>> with some research
>> papers on vector graphics using sparse matrices to
>> interpolate intelligently
>> between vector shapes, which was very interesting (still
>> have the app he
>> coded using gmm++ and eigen1, can send it to you)... so no
>> need to convince
>> me about the usefulness of all this.
>>
>> Cheers,
>> Benoit
>>
>> P.S. I'm adding vectorization to Redux
>>
>> On Monday 23 June 2008 02:09:15 Gael Guennebaud wrote:
>> > Hi,
>> >
>> > this WE I started looking at sparse matrices and I
>> started to write
>> > some experimental piece of code to see how they could
>> be integrated
>> > into the current base code. The current sparse matrix
>> class implement
>> > the "Compressed Column Storage" (CCS) scheme
>> and  does not require any
>> > change in Core: an InnerIterator template class is
>> specialized for
>> > each matrix and expression type. These iterators allow
>> to efficiently
>> > skip the zero coefficients of a row while being
>> compatible with the
>> > expression template mechanism.
>> >
>> > To check whether it was worth it, I benchmarked the
>> sum of 3 matrices
>> > of 4000x4000 with 10% of non zero coefficients
>> (random) :
>> >
>> > for (int k=0; k<10; ++k)
>> >    m4 = m1 + m2 + m3;
>> >
>> > the results:
>> >
>> > * dense matrix (with vectorization and all eigen
>> optimizations) : 0.7sec
>> >
>> > * dense matrix without xpr template, i.e.:
>> >     m4 = (m1 + m2).eval() + m3
>> >     => 1.6 sec
>> >
>> > * experimental sparse matrix with generic xpr
>> template: 1.6 sec
>> >
>> > * experimental sparse matrix without xpr template (but
>> with a
>> > specialized operator* for the sum): 2.1 sec
>> >
>> > * GMM++: 19sec !! (+ very high memory consumption)
>> >
>> > so yeah it was worth it ! Even with sparse matrix the
>> ET is pretty
>> > useful. Moreover the current implementation is far to
>> be optimized
>> > (e.g. it uses std::vector<>) so we can expect
>> even better performance.
>> >
>> > I believe the poor perf of GMM++ is partly due to the
>> use of "mapped"
>> > sparse matrix (one std::map per row) for the result
>> because their
>> > implementation of CCS matrix is read-only.
>> >
>> > cheers,
>> > gael.
>
>
>