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

Question

我正在尝试使用Cython的numpy中提供的点积，矩阵求逆和其他基本线性代数运算。诸如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）
• γ函数（与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。

编辑
：我发现一个库，将BLAS包装在Cython（https://github.com/tokyo/tokyo）中，这很接近，但不是我要找的东西。我想直接调用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应用于单位矩阵，可以类似地计算逆。有关签名，请参见此处。所有这些都需要下拉到相当底层的编码，您需要自己分配临时工作数组等。—但是，这些可以封装到您自己的便捷函数中，或者只是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 XXXXCython访问（尽管它们是私有的，所以不要这样做）。

但是，目前scipy.special不提供此C-API。但是，鉴于scipy.special中的接口模块是用Cython编写的，实际上提供它非常简单。

我认为目前尚无任何明智且可移植的方式来访问功能，以进行繁重的工作gamln（尽管您可以窥探UFunc对象，但这不是理智的解决方案：），所以目前可能最好只是从scipy.special中获取源代码的相关部分，并将其与您的项目捆绑在一起，或者使用例如GSL。