Eigen C++: Performance of sparse-matrix manipulations

Can anybody explain the following behavior of Eigen sparse matrices? I have been looking into aliasing and lazy evaluation, but I don't seem to be able to improve on the issue. Tech specs: I am using the latest Eigen stable release on Ubuntu 16.10 with g++ compiler and no optimization flags.

Suppose I define a simple identity in the following way:

SparseMatrix<double> spIdent(N,N);
spIdent.reserve(N);
spIdent.setIdentity();

then perform these operations with it

spIdent-spIdent;
spIdent*spIdent;
spIdent - spIdent*spIdent;

and measure the computation times for all three. What I get is this

0 Computation time: 2.6e-05
1 Computation time: 2e-06 
2 Computation time: 1.10706

Meaning that either operation is fast, but the combination is super slow. The noalias() method is only defined for dense matrices, plus in my dense example it did not make much of a difference. Any enlightenment?

MCVE:

#include <iostream>
#include <ctime>
#include "../Eigen/Sparse"

using namespace std;
using namespace Eigen;

int main() {

unsigned int N=2000000;

SparseMatrix<double> spIdent(N,N);
spIdent.reserve(N);
spIdent.setIdentity();

clock_t start=clock();
spIdent*spIdent;
cout << "0 Computation time: " << float(clock() - start)/1e6 << '\n';

start=clock();
spIdent-spIdent;
cout << "1 Computation time: " << float(clock() - start)/1e6 << '\n';

start=clock();
spIdent - (spIdent*spIdent);
cout << "2 Computation time: " << float(clock() - start)/1e6 << '\n';

return 0;

}

Solution

It's not so much that it gets optimized out as much as the lazy evaluation is, well, very lazy. Look at the product. The code that's called is (at least in whatever version of Eigen that was included on this machine):

template<typename Derived>
template<typename OtherDerived>
inline const typename SparseSparseProductReturnType<Derived,OtherDerived>::Type
SparseMatrixBase<Derived>::operator*(const SparseMatrixBase<OtherDerived> &other) const
{
  return typename SparseSparseProductReturnType<Derived,OtherDerived>::Type(derived(), other.derived());
}

which returns an expression of the product (i.e. lazy). Nothing is done with this expression, so the cost is zero. The same goes for the difference. Now, when doing a-a*a the a*a is an expression. It then meets the operator-. This sees the expression on the right hand side. This expression is then evaluated to a temporary (i.e. costs time) in order to use it in operator-. Why evaluate to a temporary? Read this for their logic (search for "The second circumstance").

operator- is a CwiseBinaryOp with the product expression as the right hand side. The first thing CwiseBinaryOp does is assign the right hand side to a member:

EIGEN_STRONG_INLINE CwiseBinaryOp(const Lhs& aLhs, const Rhs& aRhs, const BinaryOp& func = BinaryOp())
      : m_lhs(aLhs), m_rhs(aRhs), m_functor(func)

(m_rhs(aRhs)) which in turn calls a SparseMatrix constructor:

/** Constructs a sparse matrix from the sparse expression \a other */
template<typename OtherDerived>
inline SparseMatrix(const SparseMatrixBase<OtherDerived>& other)
  : m_outerSize(0), m_innerSize(0), m_outerIndex(0), m_innerNonZeros(0)
{
  ...
  *this = other.derived();
}

which in turn calls operator= which will (somebody correct me if I'm wrong) always trigger an evaluation, in this case, to a temporary.