在 Cython 中调用点积和线性代数运算？
Posted
技术标签:
【中文标题】在 Cython 中调用点积和线性代数运算？【英文标题】：calling dot products and linear algebra operations in Cython? 【发布时间】：2013-04-13 09:45:18 【问题描述】：我正在尝试使用 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 函数（我假设用户已经安装了这些函数。）
【问题讨论】：
你做了cimport numpy，对吧？
@LevLevitsky：是的，但我认为这不会改变重复调用 np.dot 和其他需要返回 Python 的函数的间接成本。
cimport 的想法是避免回到 Python... 但我只在文档中看到过，从未尝试过。
@LevLevitsky：不，这不正确。
我猜this 是我的意思。
【参考方案1】：
调用与 Scipy 捆绑的 BLAS “相当”简单，这是调用 DGEMM 来计算矩阵乘法的一个示例：https://gist.github.com/pv/5437087 请注意，BLAS 和 LAPACK 期望所有数组都是 Fortran 连续的（以 lda/b/c 参数为模)，因此 order="F" 和 double[::1,:] 是正确运行所必需的。
通过在单位矩阵上应用 LAPACK 函数dgesv，可以类似地计算逆。签名见here。所有这些都需要降级到相当低级的编码，您需要自己分配临时工作数组等 --- 但是这些可以封装到您自己的便利函数中，或者只是通过替换 @ 来重用 tokyo 中的代码987654328@函数通过上述方式从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 实际上确实提供了自己的协议来交换函数指针。例如，考虑 from scipy.spatial import qhull; print(qhull.__pyx_capi__) --- 这些函数可以通过 Cython 中的 from scipy.spatial.qhull cimport XXXX 访问（虽然它们是私有的，所以不要这样做）。
但是，目前scipy.special 不提供此 C-API。然而实际上提供它非常简单，因为 scipy.special 中的接口模块是用 Cython 编写的。
我认为目前没有任何合理且可移植的方式来访问为gamln 执行繁重工作的函数，（尽管您可以窥探 UFunc 对象，但这不是一个合理的解决方案 :)，所以目前最好从 scipy.special 中获取源代码的相关部分并将其与您的项目捆绑在一起，或者使用例如GSL。
【讨论】：
numpy 中的x.T 或 np.inv(x) 怎么样？我试图找到它要求转置和反转的 C 函数，但找不到。您能否解释一下代码中的 order="F" 调用以及将 numpy 数组传递给这些 C 函数时数组顺序的重要性？
谢谢。您能否详细说明“可能可以重写 tokyo 模块以使用 Scipy 导出的 BLAS/LAPACK 并将其捆绑在 scipy.linalg 中，这样您就可以从 scipy.linalg.blas cimport dgemm 中执行”请，我还是看不出这和东京有什么不同？如果我能更多地了解该解决方案将如何进行，我很乐意尝试这样做。东京不是直接叫blas吗？
东京正在直接呼叫 BLAS。 Scipy 也是如此（通过 f2py）；更好的函数（例如inv(x)）是用 Python 编写的。将类似 Tokyo 的 Cython BLAS 界面与 Scipy 捆绑会很有用，因为 tokyo 没有广泛预安装，并且必须处理与您的项目一起分发 BLAS 是一件麻烦事。或者，通过函数指针获取 BLAS 例程将是避免分发 BLAS 的另一种方法，根据您的问题，这似乎是您想要的一件事。这将需要修改东京。 （顺便说一句，我没有详细看过东京，所以我不知道它有多好。）
请注意，您实际上不需要在 Python 中设置 order='F' - 您可以使用 C 有序数组，将维度设置为 shape[1],shape[0] 在调用 dgemm 并将 transpose 参数设置为 't'。
@PatrickMineault：你能举个例子吗？如果可能，我想避免设置 order="F"【参考方案2】：
如果您接受使用 GSL，也许最简单的方法是使用这个 GSL->cython 接口 https://github.com/twiecki/CythonGSL 并从那里调用 BLAS（参见示例 https://github.com/twiecki/CythonGSL/blob/master/examples/blas2.pyx）。它还应该注意 Fortran 与 C 的排序。 没有多少新的 GSL 功能，但您可以放心地假设它已得到积极维护。 CythonGSL 比 tokyo 更完整；例如，它具有 numpy 中不存在的对称矩阵乘积。
【讨论】：
【参考方案3】：因为我刚刚遇到了同样的问题，并且写了一些额外的函数，所以我会把它们包括在这里，以防其他人发现它们有用。我编写了一些矩阵乘法，还调用了 LAPACK 函数来进行矩阵求逆、行列式和 Cholesky 分解。但是你应该考虑尝试在任何循环之外做线性代数的东西，如果你有的话，就像我做的here。顺便说一句，如果您有建议，这里的行列式函数不太有效。另外，请注意，我不会检查输入是否符合要求。
from scipy.linalg.cython_lapack cimport dgetri, dgetrf, dpotrf
cpdef void double[:, ::1] inv_c(double[:, ::1] A, double[:, ::1] B,
double[:, ::1] work, double[::1] ipiv):
'''invert float type square matrix A
Parameters
----------
A : memoryview (numpy array)
n x n array to invert
B : memoryview (numpy array)
n x n array to use within the function, function
will modify this matrix in place to become the inverse of A
work : memoryview (numpy array)
n x n array to use within the function
ipiv : memoryview (numpy array)
length n array to use within function
'''
cdef int n = A.shape[0], info, lwork
B[...] = A
dgetrf(&n, &n, &B[0, 0], &n, &ipiv[0], &info)
dgetri(&n, &B[0,0], &n, &ipiv[0], &work[0,0], &lwork, &info)
cpdef double det_c(double[:, ::1] A, double[:, ::1] work, double[::1] ipiv):
'''obtain determinant of float type square matrix A
Notes
-----
As is, this function is not yet computing the sign of the determinant
correctly, help!
Parameters
----------
A : memoryview (numpy array)
n x n array to compute determinant of
work : memoryview (numpy array)
n x n array to use within function
ipiv : memoryview (numpy array)
length n vector use within function
Returns
-------
detval : float
determinant of matrix A
'''
cdef int n = A.shape[0], info
work[...] = A
dgetrf(&n, &n, &work[0,0], &n, &ipiv[0], &info)
cdef double detval = 1.
cdef int j
for j in range(n):
if j != ipiv[j]:
detval = -detval*work[j, j]
else:
detval = detval*work[j, j]
return detval
cdef void chol_c(double[:, ::1] A, double[:, ::1] B):
'''cholesky factorization of real symmetric positive definite float matrix A
Parameters
----------
A : memoryview (numpy array)
n x n matrix to compute cholesky decomposition
B : memoryview (numpy array)
n x n matrix to use within function, will be modified
in place to become cholesky decomposition of A. works
similar to np.linalg.cholesky
'''
cdef int n = A.shape[0], info
cdef char uplo = 'U'
B[...] = A
dpotrf(&uplo, &n, &B[0,0], &n, &info)
cdef int i, j
for i in range(n):
for j in range(n):
if j > i:
B[i, j] = 0
cpdef void dotmm_c(double[:, :] A, double[:, :] B, double[:, :] out):
'''matrix multiply matrices A (n x m) and B (m x l)
Parameters
----------
A : memoryview (numpy array)
n x m left matrix
B : memoryview (numpy array)
m x r right matrix
out : memoryview (numpy array)
n x r output matrix
'''
cdef Py_ssize_t i, j, k
cdef double s
cdef Py_ssize_t n = A.shape[0], m = A.shape[1]
cdef Py_ssize_t l = B.shape[0], r = B.shape[1]
for i in range(n):
for j in range(r):
s = 0
for k in range(m):
s += A[i, k]*B[k, j]
out[i, j] = s
【讨论】：
以上是关于在 Cython 中调用点积和线性代数运算？的主要内容，如果未能解决你的问题，请参考以下文章