在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来重用tokyo中的代码. c15>具有通过上述方式从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中，这样您就可以执行from 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__) ---这些功能可以通过Cython中的from scipy.spatial.qhull cimport XXXX访问(尽管它们是私有的，所以不要这样做).

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屋！