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

Question

我正在尝试使用点积,矩阵求逆和其他基本线性代数运算,这些运算可以从Cython中获得.函数如numpy.linalg.inv(反转),numpy.dot(点积),X.t(矩阵/数组的转置).numpy.*从Cython函数调用有很大的开销,其余的函数是用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,没有numpy从Cython 调用Python模块的开销仍然会使整体事情变得更快.

我要调用的示例函数:

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

我在numpy邮件列表(https://groups.google.com/forum/?fromgroups=#!topic/cython-users/XZjMVSIQnTE)上说,如果你在大型矩阵上调用这些函数,那么就没有意义了.从Cython执行此操作,因为从numpy调用它只会导致在优化的C代码中花费大部分时间来调用numpy.但是,在我的情况下,我在小矩阵上有很多调用这些线性代数运算 - 在这种情况下,从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连续的(modulo)将LDA/b/C参数),因此,order="F"和double[::1,:]其中所需要的正确功能.

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

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

可能有可能重写tokyo模块以使用Scipy导出的BLAS/LAPACK并将其捆绑scipy.linalg,以便您可以这样做from scipy.linalg.blas cimport dgemm.如果有人想要提取请求,则接受拉取请求.

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

但是,目前scipy.special还没有提供这种C-API.然而,实际上提供它是非常简单的,因为scipy.special中的接口模块是用Cython编写的.

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