title: Fast Numpy Loops
body: <p>How do you optimize this code (<strong><em>without</em></strong> vectorizing, as this leads up to using the semantics of the calculation, which is quite often far from being non-trivial):</p>
  
  <pre><code>slow_lib.py:
  import numpy as np
  
  def foo():
      size = 200
      np.random.seed(1000031212)
      bar = np.random.rand(size, size)
      moo = np.zeros((size,size), dtype = np.float)
      for i in range(0,size):
          for j in range(0,size):
              val = bar[j]
              moo += np.outer(val, val)
  </code></pre>
  
  <p>The point is that such kind loops correspond quite often to operations where you have double sums over some vector operation.</p>
  
  <p>This is quite slow:</p>
  
  <pre><code>&gt;&gt;t = timeit.timeit('foo()', 'from slow_lib import foo', number = 10)
  &gt;&gt;print ("took: "+str(t))
  took: 41.165681839
  </code></pre>
  
  <p>Ok, so then let's cynothize it and add type annotations likes there is no tomorrow:</p>
  
  <pre><code>c_slow_lib.pyx:
  import numpy as np
  cimport numpy as np
  import cython
  @cython.boundscheck(False)
  @cython.wraparound(False)
  
  def foo():
      cdef int size = 200
      cdef int i,j
      np.random.seed(1000031212)
      cdef np.ndarray[np.double_t, ndim=2] bar = np.random.rand(size, size)
      cdef np.ndarray[np.double_t, ndim=2] moo = np.zeros((size,size), dtype = np.float)
      cdef np.ndarray[np.double_t, ndim=1] val
      for i in xrange(0,size):
          for j in xrange(0,size):
              val = bar[j]
              moo += np.outer(val, val)
  
  
  &gt;&gt;t = timeit.timeit('foo()', 'from c_slow_lib import foo', number = 10)
  &gt;&gt;print ("took: "+str(t))
  took: 42.3104710579
  </code></pre>
  
  <p>... ehr... what? Numba to the rescue!</p>
  
  <pre><code>numba_slow_lib.py:
  import numpy as np
  from numba import jit
  
  size = 200
  np.random.seed(1000031212)
  
  bar = np.random.rand(size, size)
  
  @jit
  def foo():
      bar = np.random.rand(size, size)
      moo = np.zeros((size,size), dtype = np.float)
      for i in range(0,size):
          for j in range(0,size):
              val = bar[j]
              moo += np.outer(val, val)
  
  &gt;&gt;t = timeit.timeit('foo()', 'from numba_slow_lib import foo', number = 10)
  &gt;&gt;print("took: "+str(t))
  took: 40.7327859402
  </code></pre>
  
  <p>So is there really no way to speed this up? The point is:</p>
  
  <ul>
  <li>if I convert the inner loop into a vectorized version (building a larger matrix representing the inner loop and then calling np.outer on the larger matrix) I get <em>much</em> faster code.</li>
  <li><strong>if I implement something similar in Matlab (R2016a) this performs quite well due to JIT.</strong></li>
  </ul>