title: Parallelization of outer product on multiple vectors using cython
body: <p>Suppose to have a numpy array m of shape (N, M) and I want to compute</p>
  <pre><code>res = np.zeros((M, M))
  for i in range(N):
     res += np.outer(m[i], m[i])
  </code></pre>
  <p>This loop can be made more efficient using <code>einsum</code>, i.e.</p>
  <pre><code>res = np.sum(np.einsum('ij,ik-&gt;ijk', m , m, axis=0)
  </code></pre>
  <p>but this requires the storage of a N x M x M matrix, which can be (<em>and in my case is</em>) very demanding.</p>
  <p>I thought to build this function in cython, using parallelzation</p>
  <pre><code>import numpy as np
  cimport numpy as np
  from cython.parallel import prange
  
  def get_s(double[:,:] m):
      cdef Py_ssize_t i = 0
      cdef int n = m.shape[1]
      
      res = 0.
      for i in prange(n, nogil=True):
         res += np.outer(m[i], m[i])
      return res  
  
  </code></pre>
  <p>The idea of this code is that
  The run of this code produces a lot of errors since I am using python objects, disallowed operations, and I don't know how to properly initialize <code>res</code>.</p>