title: C++ Eigen Sparse Matrix multiplication much slower than python scipy.sparse
body: <p><strong>Edit:</strong> The huge difference in performance is due to a bug in the test, when set up properly <code>Eigen</code> is 2 to 3 times faster.</p>
  
  <p>I noticed that sparse matrix multiplication using  <code>C++</code> <a href="http://eigen.tuxfamily.org/dox/group__TutorialSparse.html" rel="nofollow">Eigen</a> library is much slower than using <code>Python</code> <a href="http://docs.scipy.org/doc/scipy/reference/sparse.html" rel="nofollow">scipy.sparse</a> library. I achieve in <code>scipy.sparse</code> in <code>~0.03</code> seconds what I achieve in <code>Eigen</code> in <code>~25</code> seconds. Maybe I doing something wrong in Eigen?</p>
  
  <p>Here Python code:</p>
  
  <pre><code>from scipy import sparse
  from time import time
  import random as rn
  
  N_VALUES = 200000
  N_ROWS = 400000
  N_COLS = 400000
  
  rows_a = rn.sample(range(N_COLS), N_VALUES)
  cols_a = rn.sample(range(N_ROWS), N_VALUES)
  values_a = [rn.uniform(0,1) for _ in xrange(N_VALUES)]
  
  rows_b = rn.sample(range(N_COLS), N_VALUES)
  cols_b = rn.sample(range(N_ROWS), N_VALUES)
  values_b = [rn.uniform(0,1) for _ in xrange(N_VALUES)]
  
  big_a = sparse.coo_matrix((values_a, (cols_a, rows_a)), shape=(N_ROWS, N_COLS))
  big_b = sparse.coo_matrix((values_b, (cols_b, rows_b)), shape=(N_ROWS, N_COLS))
  
  big_a = big_a.tocsr()
  big_b = big_a.tocsr()
  
  start = time()
  
  AB = big_a * big_b;
  
  end = time()
  
  print 'time taken : {}'.format(end - start)
  </code></pre>
  
  <p>C++ code:</p>
  
  <pre><code>#include &lt;iostream&gt;
  #include &lt;cstdlib&gt;
  #include &lt;vector&gt;
  #include &lt;algorithm&gt;
  #include &lt;Eigen/Dense&gt;
  #include &lt;Eigen/Sparse&gt;
  
  using namespace Eigen;
  
  std::vector&lt;long&gt; gen_random_sample(long min, long max, long sample_size);
  double get_random_double(double min, double max);
  std::vector&lt;double&gt; get_vector_of_rn_doubles(int length, double min, double max);
  
  int main()
  {
  
    long N_COLS = 400000;
    long N_ROWS = 400000;
    long N_VALUES = 200000;
  
    SparseMatrix&lt;double&gt; big_A(N_ROWS, N_COLS);
    std::vector&lt;long&gt; cols_a = gen_random_sample(0, N_COLS, N_VALUES);
    std::vector&lt;long&gt; rows_a = gen_random_sample(0, N_COLS, N_VALUES);
    std::vector&lt;double&gt; values_a = get_vector_of_rn_doubles(N_VALUES, 0, 1);
  
    for (int i = 0; i &lt; N_VALUES; i++)
      big_A.insert(cols_a[i], cols_a[i]) = values_a[i];
    // big_A.makeCompressed(); // slows things down
  
    SparseMatrix&lt;double&gt; big_B(N_ROWS, N_COLS);
    std::vector&lt;long&gt; cols_b = gen_random_sample(0, N_COLS, N_VALUES);
    std::vector&lt;long&gt; rows_b = gen_random_sample(0, N_COLS, N_VALUES);
    std::vector&lt;double&gt; values_b = get_vector_of_rn_doubles(N_VALUES, 0, 1);
  
    for (int i = 0; i &lt; N_VALUES; i++)
      big_B.insert(cols_b[i], cols_b[i]) = values_b[i];
    // big_B.makeCompressed();
  
    SparseMatrix&lt;double&gt; big_AB(N_ROWS, N_COLS);
  
    clock_t begin = clock();
  
    big_AB = (big_A * big_B); //.pruned();
  
    clock_t end = clock();
    double elapsed_secs = double(end - begin) / CLOCKS_PER_SEC;
    std::cout &lt;&lt; "Time taken : " &lt;&lt; elapsed_secs &lt;&lt; std::endl;
  
  }
  
  std::vector&lt;long&gt; gen_random_sample(long min, long max, long sample_size)
  {
    std::vector&lt;long&gt; my_vector(sample_size); // THE BUG, is right std::vector&lt;long&gt; my_vector
  
    for (long i = min; i != max; i++)
      {
        my_vector.push_back(i);
      }
  
    std::random_shuffle(my_vector.begin(), my_vector.end());
  
    std::vector&lt;long&gt; new_vec = std::vector&lt;long&gt;(my_vector.begin(), my_vector.begin() + sample_size);
  
      return new_vec;
  }
  
  double get_random_double(double min, double max)
  {
     std::uniform_real_distribution&lt;double&gt; unif(min, max);
     std::default_random_engine re;
     double a_random_double = unif(re);
  }
  
  std::vector&lt;double&gt; get_vector_of_rn_doubles(int length, double min, double max)
  {
    std::vector&lt;double&gt; my_vector(length);
    for (int i=0; i &lt; length; i++)
      {
        my_vector[i] = get_random_double(min, max);
      }
    return my_vector;
  }
  </code></pre>
  
  <p>I compiled with: <code>g++ -std=c++11 -I/usr/include/eigen3 time_eigen.cpp -o my_exec -O2 -DNDEBUG</code>.</p>
  
  <p>Am I missing a way to do sparse multiplication fast using Eigen?</p>