在 Cython 中调用点积和线性代数运算? [英] calling dot products and linear algebra operations in Cython?

问题描述

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

I'm trying to use dot products, matrix inversion and other basic linear algebra operations that are available in numpy from Cython. Functions like numpy.linalg.inv (inversion), numpy.dot (dot product), X.t (transpose of matrix/array). There's a large overhead to calling numpy.* from Cython functions and the rest of the function is written in Cython, so I'd like to avoid this.

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

If I assume users have numpy installed, is there a way to do something like:

#include "numpy/npy_math.h"

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

as an extern, and call these functions? Or alternatively call BLAS directly (or whatever it is that numpy calls for these core operations)?

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

To give an example, imagine you have a function in Cython that does many things and in the end needs to make a computation involving dot products and matrix inverses:

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 模块的开销仍然会使整体速度更快.

how can this be done? If there's a library that implements these in Cython already, I can also use that, but have not found anything. Even if those procedures are less optimized than BLAS directly, not having the overhead of calling numpy Python module from Cython will still make things overall faster.

我想调用的示例函数:

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

• dot product (np.dot)
• matrix inversion (np.linalg.inv)
• matrix multiplication
• taking transpose (equivalent of x.T in numpy)
• gammaln function (like scipy.gammaln equivalent, which should be available in C)

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

I realize as it says on numpy mailing list (https://groups.google.com/forum/?fromgroups=#!topic/cython-users/XZjMVSIQnTE) that if you call these functions on large matrices, there is no point in doing it from Cython, since calling it from numpy will just result in the majority of the time spent in the optimized C code that numpy calls. However, in my case, I have many calls to these linear algebra operations on small matrices -- in that case, the overhead introduced by repeatedly going from Cython back to numpy and back to Cython will far outweigh the time spent actually computing the operation from BLAS. Therefore, I'd like to keep everything at the C/Cython level for these simple operations and not go through python.

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

I'd prefer not to go through GSL, since that adds another dependency and since it's unclear if GSL is actively maintained. Since I'm assuming users of the code already have scipy/numpy installed, I can safely assume that they have all the associated C code that goes along with these libraries, so I just want to be able to tap into that code and call it from Cython.

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

edit: I found a library that wraps BLAS in Cython (https://github.com/tokyo/tokyo) which is close but not what I'm looking for. I'd like to call the numpy/scipy C functions directly (I'm assuming the user has these installed.)

推荐答案

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

Calling BLAS bundled with Scipy is "fairly" straightforward, here's one example for calling DGEMM to compute matrix multiplication: https://gist.github.com/pv/5437087 Note that BLAS and LAPACK expect all arrays to be Fortran-contiguous (modulo the lda/b/c parameters), hence order="F" and double[::1,:] which are required for correct functioning.

通过在单位矩阵上应用 LAPACK 函数 dgesv 可以类似地计算逆运算.有关签名，请参阅此处.所有这一切都需要下降到相当低级的编码，你需要自己分配临时工作数组等等——但是这些可以封装到你自己的便利函数中，或者只是重用 tokyo 将 lib_* 函数替换为以上述方式从 Scipy 获取的函数指针.

Computing inverses can be similarly done by applying the LAPACK function dgesv on the identity matrix. For the signature, see here. All this requires dropping down to rather low-level coding, you need to allocate temporary work arrays yourself etc etc. --- however these can be encapsulated into your own convenience functions, or just reuse the code from tokyo by replacing the lib_* functions with function pointers obtained from Scipy in the above way.

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

If you use Cython's memoryview syntax (double[::1,:]) you transpose is the same x.T as usual. Alternatively, you can compute the transpose by writing a function of your own that swaps elements of the array across the diagonal. Numpy doesn't actually contain this operation, x.T only changes the strides of the array and doesn't move the data around.

可能可以重写 tokyo 模块以使用 Scipy 导出的 BLAS/LAPACK 并将其捆绑在 scipy.linalg 中，这样您就可以来自 scipy.linalg.blas cimport dgemm.拉取请求被接受，如果有人想深入研究的话.

It would probably be possible to rewrite the tokyo module to use the BLAS/LAPACK exported by Scipy and bundle it in scipy.linalg, so that you could just do from scipy.linalg.blas cimport dgemm. Pull requests are accepted if someone wants to get down to it.

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

As you can see, it all boils down to passing function pointers around. As alluded to above, Cython does in fact provide its own protocol for exchanging function pointers. For an example, consider from scipy.spatial import qhull; print(qhull.__pyx_capi__) --- those functions could be accessed via from scipy.spatial.qhull cimport XXXX in Cython (they're private though, so don't do that).

但是，目前 scipy.special 不提供这个 C-API.然而实际上提供它会很简单，因为 scipy.special 中的接口模块是用 Cython 编写的.

However, at the present, scipy.special does not offer this C-API. It would however in fact be quite simple to provide it, given that the interface module in scipy.special is written in Cython.

我认为目前没有任何理智和可移植的方式来访问为 gamln 做繁重工作的函数，(虽然你可以窥探 UFunc 对象，但这不是一个理智的解决方案:)，所以目前最好从 scipy.special 获取源代码的相关部分并将其与您的项目捆绑在一起，或者使用例如GSL.

I don't think there is at the moment any sane and portable way to access the function doing the heavy lifting for gamln, (although you could snoop around the UFunc object, but that's not a sane solution :), so at the moment it's probably best to just grab the relevant part of source code from scipy.special and bundle it with your project, or use e.g. GSL.

这篇关于在 Cython 中调用点积和线性代数运算?的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！