Cython with variable-length arrays

Question

Dynamically growing arrays are a type of array. They are very useful when you don't know the exact size of the array at design time. First you need to define an initial number of elements. (Wikipedia)

I have written a Python solution and converted it to Cython. Cython can be used to improve the speed of nested for loops in Python. Where my Cython code is slightly faster. My Cython solution is obviously not the fastest. I am trying to perform a nested for loop similar to the one in my Python code as fast as possible in Cython.

It would help to have some experience in C, which I don't. The main problem that I ran into is that Cython has different scoping rules to Python. Since C and Python have different scoping rules. In other words, we cannot create a new vector in the loop and assign it to the same name.

My solution works but is too slow. Can anyone improve Cython code above by using a more C like approach.

Python

import numpy as np

my_list = [1,2,3]
n = 10
a = 0.5

Estimate_1_list = []
Estimate_2_list = []

for l in my_list:

    # Resizable matrices
    a_mat = np.zeros((l,n+1),float)
    b_mat = np.zeros((l,n+1),float)
    
    for i  in range(n):
        t = i*a
        
        for j in range(l):
            
            # Fill matrices
            a_mat[j,i+1] = a_mat[j,i+1] + np.random.random()
            
            b_mat[j,i+1] = a_mat[j,i+1]/(2*t+3)
    
    # Append values of interest to use at different values of matrix size
    Estimate_1_list.append(np.mean(a_mat[:,n]))
    Estimate_2_list.append(np.std(a_mat[:,n]))   
results = [Estimate_1_list,Estimate_2_list]

Cython

import cython
# Load cython extension
%load_ext Cython

%%cython
import numpy as np

def my_function(list my_list, int n, int a ):
cdef list Estimate_1_list = []
cdef list Estimate_2_list = []
cdef int l,i,t,j
for l in my_list:

    # Resizable matrices (could I use memory view?)
    a_mat = np.zeros((l,n+1),float)
    b_mat = np.zeros((l,n+1),float)

    for i  in range(n):
        t = i*a

        for j in range(l):

            # Fill matrices
            a_mat[j,i+1] = a_mat[j,i+1] + np.random.random()

            b_mat[j,i+1] = a_mat[j,i+1]/(2*t+3)

    # Append values of interest to use at different values of matrix size
    Estimate_1_list.append(np.mean(a_mat[:,n]))
    Estimate_2_list.append(np.std(a_mat[:,n]))  
    
# Return results 
results = [Estimate_1_list,Estimate_2_list]
return results

Tests

# Test cython to show that the function is running
my_list = [1,2,3]
n = 10
a = 0.5
my_function(my_list, n, a)

[[0.13545224609230933, 0.6603542545719762, 0.6632002117071227],
 [0.0, 0.19967544614685195, 0.22125180486616808]]

Answer 1

In addition to Reinderien's answer, may I suggest:

Yes, Use Memoryviews to speed up access

In addition to the code you have, I would also type the a_mat and b_mat matrixes as double[:,::1] following the Typed Memoryviews guide. (the "1" means contiguous and is allows for slightly faster access). You are right in that you can not cdef declare variables in a loop, but you can cdef declare variables at the top of the function, and reassign within the loop.

    ...
    cdef list Estimate_1_list = []
    cdef list Estimate_2_list = []
    cdef int l,i,t,j
    cdef double[:, ::1] a_mat # declare here
    cdef double[:, ::1] b_mat
    for l in my_list:
        
        # Resizable matrices (could I use memory view?)
        a_mat = np.zeros((l,n+1),float) # assign in loop
        b_mat = np.zeros((l,n+1),float)
    ...

Turn off Boundscheck And Wrap Around check

Once you are sure that code is bug free. You can locally turn off bounds and wraparound checks using the decorators or with statements for additional speed up. (Compiler Directives, Local Directives)

Answer 2

Zeros

This:

a_mat = np.zeros

is not the right call for your purpose. You want np.empty instead, because you don't actually care what the initial values are since you do a comprehensive initialization loop right after.

Further to that: since you're adding random() to every element of a_mat, simply initialize a_mat to a single call of random() with the correct shape, rather than having to do element-wise addition.

Vectorization

You have an outer dimension (l through my_list), a second dimension (l) and a third dimension (n + 1). The second dimension is variable; the first and third are constant. This means that you can more efficiently represent this if you rearrange your dimensions so that the fixed ones are on the interior. In other words, if you had

len(my_list) = 3
l = 1, 2, 3
n = 10

then you can actually represent this as a single three-dimensional matrix of dimensions 6 * 3 * 10, where 6 == len(my_list) * (len(my_list) - 1). I think it's possible to do all of this without a single for loop, which is the ideal for vectorized performance.