python - 在 Cython 中调用点积和线性代数运算？

Question

我正在尝试使用 Cython 的 numpy 中提供的点积、矩阵求逆和其他基本线性代数运算。 numpy.linalg.inv(反转)、numpy.dot(点积)、X.t(矩阵/数组的转置)等函数。从 Cython 函数调用 numpy.* 的开销很大，而函数的其余部分是用 Cython 编写的，所以我想避免这种情况。

如果我假设用户安装了 numpy，有没有办法做类似的事情:

#include "numpy/npy_math.h"

作为一个extern，并调用这些函数？或者直接调用 BLAS(或者 numpy 调用这些核心操作的任何内容)？

举个例子，假设你在 Cython 中有一个函数可以做很多事情，最终需要进行涉及点积和矩阵求逆的计算:

cdef myfunc(...):
  # ... do many things faster than Python could
  # ...
  # compute one value using dot products and inv
  # without using 
  #   import numpy as np 
  #   np.*
  val = gammaln(sum(v)) - sum(gammaln(v)) + dot((v - 1).T, log(x).T)

如何做到这一点？如果已经有一个在 Cython 中实现这些的库，我也可以使用它，但没有找到任何东西。即使这些过程不如直接优化 BLAS，没有从 Cython 调用 numpy Python 模块的开销仍然会使事情总体上更快。

我想调用的示例函数:

• 点积(np.dot)
• 矩阵求逆(np.linalg.inv)
• 矩阵乘法
• 转置(相当于numpy中的x.T)
• gammaln 函数(类似于 scipy.gammaln 等价物，应该在 C 中可用)

我意识到，正如它在 numpy 邮件列表 (https://groups.google.com/forum/?fromgroups=#!topic/cython-users/XZjMVSIQnTE) 上所说，如果您在大型矩阵上调用这些函数，那么从 Cython 执行它是没有意义的，因为从 numpy 调用它只会导致大多数在 numpy 调用的优化 C 代码中花费的时间。然而，就我而言，我在小矩阵上多次调用这些线性代数运算——在这种情况下，反复从 Cython 回到 numpy 再回到 Cython 所带来的开销将远远超过实际计算 BLAS 操作所花费的时间。因此，对于这些简单的操作，我想将所有内容都保留在 C/Cython 级别，而不是通过 python。

我不希望通过 GSL，因为这会增加另一个依赖项，而且尚不清楚 GSL 是否得到积极维护。由于我假设代码的用户已经安装了 scipy/numpy，我可以安全地假设他们拥有与这些库一起使用的所有相关 C 代码，所以我只想能够利用该代码并调用它来自 Cython。

编辑:我找到了一个在 Cython ( https://github.com/tokyo/tokyo ) 中包装 BLAS 的库，它很接近但不是我想要的。我想直接调用 numpy/scipy C 函数(我假设用户已经安装了这些函数。)

Answer 1

调用与 Scipy 捆绑的 BLAS “相当”简单，这是调用 DGEMM 计算矩阵乘法的一个示例:https://gist.github.com/pv/5437087请注意，BLAS 和 LAPACK 期望所有数组都是 Fortran 连续的(以 lda/b/c 参数为模)，因此 order="F"和 double[::1,:]这是正常运行所必需的。

通过应用 LAPACK 函数 dgesv 可以类似地计算逆。在单位矩阵上。签名见here .所有这些都需要降级到相当低级的编码，您需要自己分配临时工作数组等 --- 但是这些可以封装到您自己的便利函数中，或者只是重用 tokyo 中的代码通过替换 lib_*通过上述方式从 Scipy 获得函数指针的函数。

如果您使用 Cython 的 memoryview 语法 (double[::1,:])，则转置是相同的 x.T照常。或者，您可以通过编写自己的函数来计算转置，该函数在对角线上交换数组的元素。 Numpy 实际上不包含此操作，x.T只会改变数组的步幅，不会移动数据。

可能会重写 tokyo模块使用 Scipy 导出的 BLAS/LAPACK 并将其捆绑到 scipy.linalg , 这样你就可以做 from scipy.linalg.blas cimport dgemm . Pull requests如果有人想深入了解，就会被接受。

---

如您所见，这一切都归结为传递函数指针。正如上面提到的，Cython 实际上确实提供了自己的协议(protocol)来交换函数指针。例如，考虑 from scipy.spatial import qhull; print(qhull.__pyx_capi__) --- 这些功能可以通过 from scipy.spatial.qhull cimport XXXX 访问在 Cython 中(虽然它们是私有(private)的，所以不要这样做)。

但是，目前，scipy.special不提供此 C-API。然而实际上提供它非常简单，因为 scipy.special 中的接口(interface)模块是用 Cython 编写的。

我认为目前没有任何合理且可移植的方式来访问为 gamln 执行繁重工作的功能。 , (尽管您可以窥探 UFunc 对象，但这不是一个理智的解决方案 :)，因此目前最好从 scipy.special 中获取源代码的相关部分并将其与您的项目捆绑在一起，或者使用例如GSL。