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

谢谢！

简短概述：如何在Cython中实现线性代数？

- BVJ

---

在BLAS方面，你看过tokyo了吗？ - ali_m

是的，但它非常有限，并且似乎没有得到维护（我认为上次提交是两年前）。 - BVJ

尽管Tokyo还不完整，但我认为其一般方法仍然非常有启示性。您可以使用当前的scipy开发版本，该版本似乎具有BLAS和LAPACK的Cython包装器，或者在此期间，您可以编写一个特定的.pxd头文件来公开所需的特定例程。 - ali_m

这些数组的维度是多少？如果足够小，最快的方法可能是在C中完整地写出代数表达式。 - hpaulj

请注意，Cython接口最近在2015 Scipy conf上进行了展示：https://www.youtube.com/watch?v=R4yB-8tB0J0 - Caleb Hattingh

Answer 1

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

如果在重复调用的BLAS函数移植到Cython后没有体验到任何加速，那么无论是numpy中执行此操作的Python开销不重要（例如，如果计算本身是最昂贵的部分），还是您做错了一些事情（不必要的内存复制等）。

我想说，与GSL相比，这种方法应该更快且更易于维护，前提是您使用优化的BLAS（OpenBLAS、ATLAS、MKL等）编译scipy。

- rth

---

嗨，谢谢。似乎最后一部分是问题所在。我安装了ATLAS，现在结果好多了。做矩阵乘法确实需要做很多工作，因此scipy BLAS Cython API可以节省很多工作量。顺便说一句，使用Tokyo的无验证函数也可以得到非常相似的结果，因此对于那些只需要简单操作的人来说，使用它肯定更容易。 - BVJ

Answer 2

作为“einsum”补丁的一部分，我编写了一个Python模拟函数的“c”代码。虽然我在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 - 你可能能够使用它们，或者像使用 numpy 数组一样使用 cython 数组。这份文档举例说明，你可以通过 [None,:] 添加一个维度。那是否意味着你可以通过广播来执行常见的 numpy 外积操作？

- hpaulj

---

谢谢。正如我之前所写的，我可以像你一样重新用C或Cython编写操作，但我真的相信那句话“不要重复发明轮子，除非你真的想了解轮子”。另外，关于memoryviews的观点，据我所知，你不能直接相乘memoryviews，所以你最终还是会循环。我认为通过混合cpython数组和东京实现，我可能会有所发现……如果我得到任何有趣的东西，我会发布出来。 - BVJ