python - Cython 中用于循环调用的快速基本线性代数

Question

我正在尝试在 cython 中为蒙特卡洛模拟编写一个函数。该函数涉及多个小型线性代数运算，如点积和矩阵求逆。由于该函数被调用了数十万次，因此 numpy 开销占了很大一部分成本。 三年前有人问过这个问题:calling dot products and linear algebra operations in Cython? 我尝试使用这两个答案中的建议，但第一个 scipy.linalg.blas 仍然通过 python 包装器，我并没有真正得到任何改进。第二，使用 gsl 包装器也相当慢，并且当向量的维度非常大时往往会卡住我的系统。我还找到了 Ceygen 包，它看起来很有前途，但安装文件似乎在上次 Cython 更新中损坏了。 另一方面，我看到 scipy 正在为 lapack 开发一个 cython 包装器，但它看起来仍然不可用(scipy-cython-lapack) 最后，我还可以为这些操作编写自己的 C 例程，但这似乎是在重新发明轮子。

总结一下:在Cython中有没有一种新的方式来处理这种操作？ (因此我不认为这是重复的)或者您是否找到了一种更好的方法来处理我还没有见过的这类问题？

强制性代码示例: (这只是一个例子，当然还可以改进，但只是给出思路)

 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

谢谢!!

tl;dr:如何使用 cython 进行线性代数？

Answer 1

答案you link to仍然是从 Cython 调用 BLAS 函数的好方法。它并不是真正的 python 包装器，Python 只是用于获取指向该函数的 C 指针，这可以在初始化时完成。所以你基本上应该获得类似 C 的速度。我可能是错的，但我认为即将发布的 Scipy 0.16 版本将提供一个方便的 BLAS Cython API，基于这种方法，它不会改变性能方面的事情。

如果在将重复调用的 BLAS 函数移植到 Cython 后您没有体验到任何加速，那么在 numpy 中执行此操作的 python 开销无关紧要(例如，如果计算本身是最昂贵的部分)或者您是做错事(不必要的内存副本等)

我会说这种方法应该比 GSL 更快、更容易维护，当然前提是您使用优化的 BLAS(OpenBLAS、ATLAS、MKL 等)编译 scipy。