Compile Scipy function with Cython

Question

I am running a simulation in Python 3.4 - that involves a lot of dot products between a sparse array (in csr format) and a dense vector. I am using Scipy for the sparse matrix, numpy for everything else.

Using Cython gave me a massive boost (~x6 speed increase), after making sure that I cdef everything properly and after minimizing Python interaction (bt going through the html file that Cython gives me and modifying my code).

Now, I profile the code and 50% of the simulation time is spent on the line with the dot product. I am wondering if it is possible to somehow accelerate this line, say by complining this one dot function in Cython?

I know I could write my own implementation for (csr sprase 2d matrix) dot (dense vector), but I am trying to avoid that.

Edit: I have included a minimal example of the code. I am sorry, I can't see how to make it smaller. It is a textbook exercise in statistical mechanics. Place marbles in pots until one of the pots exceeds capacity. Then, start a cascade which propagates according to a (here sparse) matrix. I am using batch sampling.

Please focus on the line towards the end.

from __future__ import division
import  numpy           as np
import  cython
cimport numpy           as np
cimport cpython.array
import  scipy.sparse    as sps


@cython.cdivision(True)
@cython.nonecheck(False)
@cython.boundscheck(False)
@cython.wraparound(False)
def simulate(long[:] capacity_vec,
             int random_array_size,
             long n,
             int seed,
             int[:] A_col,
             int[:] A_row,
             long[:] A_data):


    #### Initialise ####################################################################################################

    # Initialise states
    cdef int[:] damage  = np.random.randint(0, int(np.min(capacity_vec)/2), n).astype(np.int32)
    cdef int[:] dr_list = np.random.choice(n, 1000).astype(np.int32)
    cdef int[:] states  = np.zeros(n).astype(np.int32)
    cdef int[:] states_ = np.zeros(n).astype(np.int32)
    cdef int[:] change  = np.zeros(n).astype(np.int32)

    # Initialise counters
    cdef int k, violations, violations_, counter= 0, dr_id=0, increment_index = 0


    # Build Sparse Adjecency Matrix
    cA_sps = sps.csr_matrix( (A_data, (A_row, A_col) ), shape=(n,n) ).astype(np.int32)


    while counter < 1000:

        #### Place damage until a cascade starts #######################################################################
        while damage[increment_index] <= capacity_vec[increment_index]:# Check for violations

            increment_index         = dr_list[dr_id]                   # Where do we place the marble?

            damage[increment_index] = damage[increment_index] + 1      # place the marble

            dr_id                   = dr_id + 1                        # another random number used

            if dr_id == random_array_size - 1:                         # Check if we run out of random numbers

                dr_list = np.random.choice(n, random_array_size).astype(np.int32) # if so, pick new increment_index

                dr_id   = 0                                            # and reset the counter


        #### Initialise cascade ########################################################################################
        violations, violations_  = 1, 0
        states[increment_index]  = 1


        #### Propagate cascade #########################################################################################
        while violations > violations_:                                # check for fixed point, propagate cascade
            for k in range(n): change[k] = states[k] - states_[k]
            ### THIS LINE IS THE PROBLEM. It takes up half of all simulation time.
            damage      = damage + cA_sps.dot(change).astype(np.int32) # spread violations  

            states_     = states.copy()                                # store previous states

            # Determine previous and current violations
            violations, violations_ = 0 , violations

            for k in range(n):

                states_[k]  = 0

                if damage[k] > capacity_vec[k]:

                    violations = violations + 1

                    states[k]  =  1                                    # deactivate any node that has a violation


        for k in range(n): damage[k] = 0
        counter  = counter + 1                                         # progress cascade id after storing

Show us the code you are asking about, especially the one dot function. — J_H
csr matrix product already uses cython. It takes different routes depending whether the other is vector or 2d and dense or sparse. So tracing the code will take some work. — hpaulj
I will edit the code down to a simple example - preserving the dot product - and edit my post. I didn't do that because I think it will still be too long. — Dionysios Georgiadis
I don't think you can beat dot-multiplication of scipy, it is already very fast. I would look somewhere else for a potential speed-up. — ead
I am fully confident that I won't be smarter than the Scipy folks. The reason I suggest to implement my own dot product is that I could do that in Cython. Currently, It is my impression that there is a lot of overhead in Cython interfacing with Python. I would avoid that with a cython_dot_product function. No? — Dionysios Georgiadis

Dawid Laszuk Dawid Laszuk · Accepted Answer · 2017-11-21T07:36:53

I'd discourage from writing your own matrix multiplication. SciPy is done by smart people who know what they are doing, and unless your confident in numerical computing, just don't. Most of SciPy code is already compiled.

However, what you might look at is code for sparse.csr_matrix.dot. Getting into definition directly here and then here, you'll see that there are few checks done in Scipy. If you know what exact form you want, you could write your own method (modify your SciPy copy) and compute your product directly. Not sure how much that would help, though.

If you want to build Scipy yourself it is as easy as checking out whole project from GitHug and then running

python setup.py build
python setup.py install

For more direct instructions check build documentation.

Compile Scipy function with Cython

1 Answers