title: Scipy sparse matrix exponentiation: a**16 is slower than a*a*a*a*a*a*a*a*a*a*a*a*a*a*a*a?
body: <p>I'm doing a simple sparse matrix exponentiation, <code>a**16</code>, using scipy-0.17. (Note, <em>not</em> element-wise multiplication). However, on my machines (running Debian stable and Ubuntu LTS) this is ten times slower than using a for loop or doing something silly like <code>a*a*a*a*a*a*a*a*a*a*a*a*a*a*a*a</code>. This doesn't make sense, so I assume I'm doing something wrong, but what?</p>
  
  <pre><code>import scipy.sparse
  from time import time
  
  a=scipy.sparse.rand(2049,2049,.002)
  
  print ("Trying exponentiation (a**16)")
  t=time()
  x=a**16
  print (repr(x))
  print ("Exponentiation took %f seconds\n" % (time()-t))
  
  print ("Trying expansion (a*a*a*...*a*a)")
  t=time()
  y=a*a*a*a*a*a*a*a*a*a*a*a*a*a*a*a
  print (repr(y))
  print ("Expansion took %f seconds\n" % (time()-t))
  
  print ("Trying a for loop (z=z*a)")
  t=time()
  z=scipy.sparse.eye(2049)
  for i in range(16):
      z=z*a
  print (repr(z))
  print ("Looping took %f seconds\n" % (time()-t))
  
  # Sanity check, all approximately the same answer, right? 
  assert (abs(x-z)&gt;=1e-9).nnz==0
  assert (abs(x-y)&gt;=1e-9).nnz==0
  </code></pre>