Speed-up cython code
<p>I have code that is working in python and want to use cython to speed up the calculation. The function that I've copied is in a .pyx file and gets called from my python code. V, C, train, I_k are 2-d numpy arrays and lambda_u, user, hidden are ints. 
I don't have any experience in using C or cython. What is an efficient 
way to make this code faster.
Using <code>cython -a</code> for compiling shows me that the code is flawed but how can I improve it. Using <code>for i in prange (user_size, nogil=True):</code>
results in <code>Constructing Python slice object not allowed without gil</code>.</p>

<p>How has the code to be modified to harvest the power of cython?</p>

<pre><code> @cython.boundscheck(False)
 @cython.wraparound(False)
 def u_update(V, C, train, I_k, lambda_u, user, hidden):
    cdef int user_size = user
    cdef int hidden_dim = hidden
    cdef np.ndarray U = np.empty((hidden_dim,user_size), float)
    cdef int m = C.shape[1]

    for i in range(user_size):
        C_i = np.zeros((m, m), dtype=float)
        for j in range(m):
            C_i[j,j]=C[i,j]

        U[:,i] = np.dot(np.linalg.inv(np.dot(V, np.dot(C_i,V.T)) + lambda_u*I_k), np.dot(V, np.dot(C_i,train[i,:].T)))

    return U
</code></pre>

<p>You are trying to use <code>cython</code> by diving into the deep end of pool.  You should start with something small, such as some of the numpy examples.  Or even try to improve on <code>np.diag</code>.</p>

<pre><code>    i = 0
    C_i = np.zeros((m, m), dtype=float)
    for j in range(m):
        C_i[j,j]=C[i,j]
</code></pre>

<p>v.</p>

<pre><code>    C_i = diag(C[i,:])
</code></pre>

<p>Can you improve the speed of this simple expression?  <code>diag</code> is not compiled, but it does perform an efficient indexed assignment. </p>

<pre><code> res[:n-k].flat[i::n+1] = v
</code></pre>

<p>But the real problem for <code>cython</code> is this expression:</p>

<pre><code>U[:,i] = np.dot(np.linalg.inv(np.dot(V, np.dot(C_i,V.T)) + lambda_u*I_k), np.dot(V, np.dot(C_i,train[i,:].T)))
</code></pre>

<p><code>np.dot</code> is compiled.  <code>cython</code> won't turn that in to <code>c</code> code, nor will it consolidate all 5 <code>dots</code> into one expression.  It also won't touch the <code>inv</code>.  So at best <code>cython</code> will speed up the iteration wrapper, but it will still call this Python expression <code>m</code> times.</p>

<p>My guess is that this expression can be cleaned up.  Replacing the inner <code>dots</code> with <code>einsum</code> can probably eliminate the need for <code>C_i</code>. The <code>inv</code> might make 'vectorizing' the whole thing difficult.  But I'd have to study it more.  </p>

<p>But if you want to stick with the <code>cython</code> route, you need to transform that <code>U</code> expression into simple iterative code, without calls to numpy functions like <code>dot</code> and <code>inv</code>.</p>

<p>===================</p>

<p>I believe the following are equivalent:</p>

<pre><code>np.dot(C_i,V.T)
C[i,:,None]*V.T
</code></pre>

<p>In:</p>

<pre><code>np.dot(C_i,train[i,:].T) 
</code></pre>

<p>if <code>train</code> is 2d, then <code>train[i,:]</code> is 1d, and the <code>.T</code> does nothing.</p>

<pre><code>In [289]: np.dot(np.diag([1,2,3]),np.arange(3))
Out[289]: array([0, 2, 6])
In [290]: np.array([1,2,3])*np.arange(3)
Out[290]: array([0, 2, 6])
</code></pre>

<p>If I got that right, you don't need <code>C_i</code>.</p>

<p>======================</p>

<p>Furthermore, these calculations can be moved outside the loop, with expressions like (not tested)</p>

<pre><code>CV1 = C[:,:,None]*V.T   # a 3d array
CV2 = C * train.T  

for i in range(user_size):
    U[:,i] = np.dot(np.linalg.inv(np.dot(V, CV1[i,...]) + lambda_u*I_k), np.dot(V, CV2[i,...]))
</code></pre>

<p>A further step is to move both <code>np.dot(V,CV...)</code> out of the loop.  That may require <code>np.matmul</code> (@) or <code>np.einsum</code>.   Then we will have</p>

<pre><code>for i...
    I = np.linalg.inv(VCV1[i,...])  
    U[:,i] = np.dot(I+ lambda_u), VCV2[i,])
</code></pre>

<p>or even</p>

<pre><code>for i...
     I[...i] = np.linalg.inv(...) # if inv can't be vectorized
U = np.einsum(..., I+lambda_u, VCV2)
</code></pre>

<p>This is a rough sketch, and details will need to be worked out.</p>

python numpy cython