performance tuning on Eigen sparse matrix

Question

I've implemented something using Eigen's SparseMatrix, basically it's something like,

SparseMatrix W;
...
W.row(i) += X.row(j);  // X is another SparseMatrix, both W and X are row major.
...

and I did some perf-profiling on the code via google-pprof, and I think the above code is problematic, see figure below,

fig 1

then fig 2

finally fig 3

looks like the operator+= brings in much memory-copy stuff.

I don't know much about the internals of SparseMatrix operations, but is there any recommended way to optimize the above code?

Does W.row(i) and X.row(j) have the same sparsity patterns?? If yes, then this could indeed be optimized. — ggael
@ggael, yes, they're supposed to be, at least the each row of W should cover whatever sparse rows in X, which means, in terms of sparsity patterns, each single row of X should be a subset of any row in W. — avocado
@ggael, if W is a dense matrix, the above W.row(i) += X.row(j) could still be optimized further? Also I wonder if eigen's SparseMatrix elem-wise addition is less performant than ordinary c++ impl using stl hashmap or vector? — avocado
If W is a dense matrix (e.g., MatrixXd) then this should be fine. If W is a sparse matrix but you are sure that W.row(i) += X.row(j) can be done in-place without any reallocation/copies, then this could be optimized within Eigen itself by (1) extending the API to let Eigen knows and (2) writing the respective in-place evaluation code. Meanwhile, you could write your own inplace evaluation. I'll post an exemple. — ggael

ggael ggael · Accepted Answer · 2017-10-09T20:21:15

If the sparsity of X is a subset of the sparsity of W, then you can wrote your own function doing the addition in-place:

namespace Eigen {
template<typename Dst, typename Src>
void inplace_sparse_add(Dst &dst, const Src &src)
{
  EIGEN_STATIC_ASSERT( ((internal::evaluator<Dst>::Flags&RowMajorBit) == (internal::evaluator<Src>::Flags&RowMajorBit)),
                      THE_STORAGE_ORDER_OF_BOTH_SIDES_MUST_MATCH);

  using internal::evaluator;
  evaluator<Dst> dst_eval(dst);
  evaluator<Src> src_eval(src);

  assert(dst.rows()==src.rows() && dst.cols()==src.cols());
  for (Index j=0; j<src.outerSize(); ++j)
  {
    typename evaluator<Dst>::InnerIterator dst_it(dst_eval, j);
    typename evaluator<Src>::InnerIterator src_it(src_eval, j);
    while(src_it)
    {
      while(dst_it && dst_it.index()!=src_it.index())
        ++dst_it;
      assert(dst_it);
      dst_it.valueRef() += src_it.value();
      ++src_it;
    }
  }
}
}

Here is a usage example:

int main()
{
  int n = 10;
  MatrixXd R = MatrixXd::Random(n,n);
  SparseMatrix<double, RowMajor> A = R.sparseView(0.25,1), B = 0.5*R.sparseView(0.65,1);

  cout << A.toDense() << "\n\n" << B.toDense() << "\n\n";

  inplace_sparse_add(A, B);

  cout << A.toDense() << "\n\n";

  auto Ai = A.row(2);
  inplace_sparse_add(Ai, B.row(2));

  cout << A.toDense() << "\n\n";
}

performance tuning on Eigen sparse matrix

1 Answers