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中有new这种操作的方法吗？（因此我不认为这是重复的）或者你有没有找到更好的方法来处理我还没有见过的这类问题？

强制性代码示例：（这只是一个例子，当然还是可以改进的，只是给个思路）

 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

作为 einsum 补丁的一部分，我编写了函数 c 代码的 Python 模拟。当我在 Python 中解析字符串时，我使用 cython 来编写积和计算。

pyx 在 https://github.com/hpaulj/numpy-einsum/blob/master/sop.pyx

您甚至可以通过 Python 调用来调用编译 einsum 模块。有很多SO题比较np.dot和np.einsum，不同的计算比较好

http://docs.cython.org/src/userguide/memoryviews.html

另请参阅 cython memoryviews - 您也许可以使用这些或 cython 数组，就像使用 numpy 数组一样。例如，该文档提到您可以使用 [None,:] 添加维度。这是否意味着您可以通过广播进行常见的 numpy 外积？

Answer 2

答案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。

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

Fast basic linear algebra in Cython for recurrent calls

python

numpy

cython

scipy