title: Why is numpy's einsum faster than numpy's built in functions?
body: <p>Lets start with three arrays of <code>dtype=np.double</code>. Timings are performed on a intel CPU using numpy 1.7.1 compiled with <code>icc</code> and linked to intel's <code>mkl</code>. A AMD cpu with numpy 1.6.1 compiled with <code>gcc</code> without <code>mkl</code> was also used to verify the timings. Please note the timings scale nearly linearly with system size and are not due to the small overhead incurred in the numpy functions <code>if</code> statements these difference will show up in microseconds not milliseconds:</p>
  
  <pre><code>arr_1D=np.arange(500,dtype=np.double)
  large_arr_1D=np.arange(100000,dtype=np.double)
  arr_2D=np.arange(500**2,dtype=np.double).reshape(500,500)
  arr_3D=np.arange(500**3,dtype=np.double).reshape(500,500,500)
  </code></pre>
  
  <p>First lets look at the <code>np.sum</code> function:</p>
  
  <pre><code>np.all(np.sum(arr_3D)==np.einsum('ijk-&gt;',arr_3D))
  True
  
  %timeit np.sum(arr_3D)
  10 loops, best of 3: 142 ms per loop
  
  %timeit np.einsum('ijk-&gt;', arr_3D)
  10 loops, best of 3: 70.2 ms per loop
  </code></pre>
  
  <p>Powers:</p>
  
  <pre><code>np.allclose(arr_3D*arr_3D*arr_3D,np.einsum('ijk,ijk,ijk-&gt;ijk',arr_3D,arr_3D,arr_3D))
  True
  
  %timeit arr_3D*arr_3D*arr_3D
  1 loops, best of 3: 1.32 s per loop
  
  %timeit np.einsum('ijk,ijk,ijk-&gt;ijk', arr_3D, arr_3D, arr_3D)
  1 loops, best of 3: 694 ms per loop
  </code></pre>
  
  <p>Outer product:</p>
  
  <pre><code>np.all(np.outer(arr_1D,arr_1D)==np.einsum('i,k-&gt;ik',arr_1D,arr_1D))
  True
  
  %timeit np.outer(arr_1D, arr_1D)
  1000 loops, best of 3: 411 us per loop
  
  %timeit np.einsum('i,k-&gt;ik', arr_1D, arr_1D)
  1000 loops, best of 3: 245 us per loop
  </code></pre>
  
  <p>All of the above are twice as fast with <code>np.einsum</code>. These should be apples to apples comparisons as everything is specifically of <code>dtype=np.double</code>. I would expect the speed up in an operation like this:</p>
  
  <pre><code>np.allclose(np.sum(arr_2D*arr_3D),np.einsum('ij,oij-&gt;',arr_2D,arr_3D))
  True
  
  %timeit np.sum(arr_2D*arr_3D)
  1 loops, best of 3: 813 ms per loop
  
  %timeit np.einsum('ij,oij-&gt;', arr_2D, arr_3D)
  10 loops, best of 3: 85.1 ms per loop
  </code></pre>
  
  <p>Einsum seems to be at least twice as fast for <code>np.inner</code>, <code>np.outer</code>, <code>np.kron</code>, and <code>np.sum</code> regardless of <code>axes</code> selection. The primary exception being <code>np.dot</code>  as it calls DGEMM from a BLAS library. So why is <code>np.einsum</code> faster that other numpy functions that are equivalent?</p>
  
  <p>The DGEMM case for completeness:</p>
  
  <pre><code>np.allclose(np.dot(arr_2D,arr_2D),np.einsum('ij,jk',arr_2D,arr_2D))
  True
  
  %timeit np.einsum('ij,jk',arr_2D,arr_2D)
  10 loops, best of 3: 56.1 ms per loop
  
  %timeit np.dot(arr_2D,arr_2D)
  100 loops, best of 3: 5.17 ms per loop
  </code></pre>
  
  <hr>
  
  <p>The leading theory is from @sebergs comment that <code>np.einsum</code> can make use of <a href="http://en.wikipedia.org/wiki/SSE2">SSE2</a>, but numpy's ufuncs will not until numpy 1.8 (see the <a href="https://github.com/numpy/numpy/blob/master/doc/release/1.8.0-notes.rst">change log</a>). I believe this is the correct answer, but have <em>not</em> been able to confirm it. Some limited proof can be found by changing the dtype of input array and observing speed difference and the fact that not everyone observes the same trends in timings.</p>