Fast Numpy Loops
<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>