Re: [eigen] Sparse usage feedback

[Gael Guennebaud <gael.guennebaud@xxxxxxxxx>]

On Mon, Oct 18, 2010 at 1:20 PM,  <bernard.hugueney@xxxxxxxxxx> wrote:


> I meant that you can have access to all the relevent info with your
>
> current API,
>
> but not through an "itarator API" because you get them by calling member
>
> functions
>
> it.value();
>
> it.row();   // row index
>
> it.col();   // col index (here it is equal to k)
>
> it.index(); // inner index, here it is equal to it.row()
>
>
>
> which means that you cannot replace your loops with calls to std
>
> algorithms.
>
> I.e, I'd like to rewrite your previous example as :
>
> SparseMatrix m;
>
>
>
> for(int j=0;j<m.outerSize();++j)
>
>  std::for_each(m.innerVector(j).begin(), m.innerVector(j).end(), …);
>
>
>
> But it means that the iterators must return the relevent info through the
>
> idiomatic operator*() instead of custom member functions (.value(),…).
>
>
>
> I wrote iteratorS because as you only have one operator*(), you'd have to
>
> give more than one iterator, or you'd make your user pay for unused data.
>
> (i.e you could give rowsCols(j).begin()/.end() with operator*() returning
>
> and std::pair<IndexType, IndexType>
>
> and values(j).begin()/.end() returning a Scalar.
>


ok, but in many cases you need both the values and the coordinates in
row/col or inner/outer. For instance, in Eigen itself all usages of
the InnerIterator make use of the trio value/outer/inner.  So we would
need a ton of different iterators to satisfy all needs. Also
std::pair<> is a very bad thing as it makes code unreadable. What is
first() ? what is second() ?

Actually, iterators on values and "inner indices" are very easy to do
since they are just Scalar* and Index*. For instance something like:

for(int j=0;j<m.outerSize();++j)
  std::for_each(m.innerVector(j)._valuePtr(),
m.innerVector(j)._valuePtr()+m.innerVector(j).nonZeros(), …);

should already work for sparse objects with storage.

Also in case I move from the current iterator mechanism to the
begin/end paradigm (while keeping the rest intact, i.e., with the full
set of attributes like value(), row(), col()), then it would be
straightforward to write custom adaptators to STL compatible
iterators. To give an idea:

template<typename EigenIterator>
class value_adaptator {
  EigenIterator m_it;
public:
  value_adaptator(const EigenIterator& it) : m_it(it) {}
  operator++() { m_it++; }
  Scalar& operator* () {return m_it.value();}
.....
};

> Failing to do that means giving up on the std algorithms (or having to
>
> rewrite them
>
> as does GIL :
>
> http://stlab.adobe.com/gil/html/giltutorial.html#STLEquivalentsSec


Actually, Eigen already allows you to perform a lot of operations, and
I'd be interested to know what you are missing to have to fallback to
the STL algorithms. Actually, I plan to add a method returning the non
zeros as a dense Eigen vector. We could also to the same with the
indices. I guess that that alone would already fill many of your needs
?

> Btw, I'll also investigate whether it would be prohibitively expensive to
>
> allow
>
>  std::for_each(m.begin(), m.end(), …);
>
> to walk the whole matrix.(IMVHO, for sparse matrices, you don't usually
>
> fit many
>
> matrices in L1 cache so the extract checks will be lost in the data fecth
>
> latency).


I think the overhead would be zero most of the time.


>>> 3°) This is not possible to construct a Sparse from a Dense
>
>>>
>
>>
>
>> This is already possible:
>
>>
>
>> sparse = dense.sparseView();
>
>> sparse = dense.sparseView(ref,epsilon); // skip coefficients smaller
>
>> than ref*epsilon
>
>
>
> Thanks. I missed (and I think it should deserve some ref in
>
> http://eigen.tuxfamily.org/dox-devel/TutorialSparse.html
>
>
>
> However, I think that Sparse classes should have constructors taking dense
>
> matrices
>
> calling sparseview() because if you can construct an object of class A
>
> from
>
> an object of classe B, the idiomatic way to do it in C++ is to have a A(B
>
> const&)
>
> constructor.

Yes then we could also have an explicit constructor pruning elements
which are exactly 0.

I think this was the plan, so patch welcome :)


Note that the sparseView approach is still needed because it offers
many more features:

* control of the epsilon
* it does not create any temporary object, this is just a "view" expression..

For instance:

DenseMatrix D;
SparseMatrix S;

S * D; // valid, use a sparse x dense algorithm and returns a dense object

S * D.sparseView(); // valid, use a sparse x sparse algorithm, returns
a sparse object, D is used directly without intermediate temporary

S * SparseMatrix(D); // create a temporary sparse matrix object as a
copy of the non zeros of D, and returns a sparse object

>>> 4°) I think some scalar operation traits are missing
>
>>>
>
>>> (i.e. trying to multiply an int matrix with a double fails to find
>
>>>
>
>>> ‘struct Eigen::ei_scalar_product_traits<double, int>’ )
>
>>
>
>> This is because this is explicitly not allowed. You have to explicitly
>
>> cast the int matrix to double. Our opinion is that such an automagic
>
>> cast would be dangerous.
>
>
>
> When sizeof(int)*8 <= std::numeric_limits<double>::digits I believe that
>
> such
>
> a cast is not dangerous (and I have 32 bits for ints and 53 bits for
>
> double).


Yes, in this case this is safe regarding the accuracy, but I think
that when we hit such a mix of integer/doubles, there are 3
possibilities:

1 - the user made a mistake (typo), and this is not what he wanted to
do (might be hard to debug if we allow this)

2 - this is indeed what the user wanted to do, but if we silently do
the (very) *expensive* cast for him he won't notice it while maybe, it
would be better (faster) to directly use double, or make a copy as
double because the cast happens many times of the same data. So the
idea is that by enforcing the user to write the cast we somehow
enforce him to think a bit about his algorithm and the appropriate use
of scalar types.

3 - justified use case. Even in this case I think it is good thing to
not hide the cast to help readers to better understand the code and
what really happen.


> However, I think maybe a catch all template could be used to give an hint
>
> to the
>
> problem with an helpful message in the source code of this template
>
> (designed to trigger an error)
>
> such as /* no implicit cast provided, use the .cast<> template member
>
> function
>
> to mix those types*/).


I thought this was already the case, I'll check.

>
>>actually I've already started a rewrite of them with a
>
>> slightly simpler/safer API. The main change is that I'm removing these
>
>> Backend template parameters.
>
>
>
> I am very interested in this !
>
> Where is it availbale ?
>

currently in my local hard drive, I'll commit/push it soon for discussions.


>
> Currently, I'm giving the backend choice at runtime and deriving
>
> the implementations with dynamic dispatching.
>
>    template<typename Backend> struct LUFact : LUFactBase
>
>
>
> and I convert the returned matrices types to a Backend-independant
>
> sparse type.
>
>
>
> As I wrote above, I think that sparses matrices (contrary to dense
>
> matrices)
>
> are never small enough to mandate the all same optimizations as dense
>
> matrices.
>
> For instance, I can't think of a relevent use case forbidding dynamic
>
> dispatching
>
> (maybe using covariant return type if various backend return different
>
> concrete
>
> implementations for the result matrices of the LU decompositions).


Yes you are right, for sparse objects dynamic dispatching is fine.
Typically, the shape (selfadjoint-ness, triangular-ness, etc.) of the
matrix could also be runtime flags... This puzzled me for a while. The
problem is that I would also like that the API is as close as possible
to the dense API, such that it is possible to write codes which work
for both worlds... So I kept the triangular/selfadjoint-View
mechanisms.

For the decompositions, the plan is to have one "low level" class per
algorithm, with a common API (compute()/solve()) and specific
fine-tunning features. For instance the new CholmodDecomposition has
options to control some basic aspects such as whether you want a
LLt/LDLt, simplicial/supernodal, etc. plus a full access to
cholmod_common object.

Then we could imagine to have a unified solver class built on top of
these low level classes with runtime dispatching as you did with a
runtime manager to register/remove solvers on the fly, etc. I see it
as a separated module.

gael