I am trying to program a function in cython to simulate monte-carlo. The function includes several small operations with linear algebra, such as point products and matrix inversions. Since this function is called a hundred thousand times, numpy overheads get most of the cost. Three years ago, someone asked this question: call point products and linear algebra operations in Cython?
I tried using the recommendations of both answers, but the first scipy.linalg.blas still goes through the python shell, and I am not getting any improvements. Secondly, using the gsl wrapper is also quite slow and tends to freeze my system when the sizes of the vectors are very large. I also found a Ceygen package that looked very promising, but it seems that the installation file broke in the latest Cython update. On the other hand, I saw that scipy is working on a cython shell for lapack, but it looks like it is not yet available ( scipy-cython-lapack)
In the end, I can also code my own C routines for these operations, but it looks like they reinvent the wheel.
So, to summarize: is there a new way for this kind of operations in Cython? (Therefore, I do not think this is a duplicate) Or did you find a better way to deal with a problem that I have not yet seen?
Example of required code: (This is just an example, of course, it can still be improved, but just to give an idea)
cimport numpy as np
import numpy as np
cpdef double risk(np.ndarray[double, ndim=2, mode='c'] X,
np.ndarray[double, ndim=1, mode='c'] v1,
np.ndarray[double, ndim=1, mode='c'] v2):
cdef np.ndarray[double, ndim=2, mode='c'] tmp, sumX
cdef double ret
tmp = np.exp(X)
sumX = np.tile(np.sum(tmp, 1).reshape(-1, 1), (1, tmp.shape[0]))
tmp = tmp / sumX
ret = np.inner(v1, np.dot(X, v2))
return ret
Thank!!
tl; dr: practical linear algebra in cython?