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

Question

我试着使用点积，矩阵求逆和其他基本的线性代数运算，这些都可以从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调用numpyPython模块的开销仍然会使整体速度更快。

我要调用的示例函数：

• 点积（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函数（我假设用户已经安装了这些函数）

Answer 1

由于我刚刚遇到了同样的问题，并且编写了一些附加函数，所以我将在这里包含它们，以防其他人发现它们有用。我编写了一些矩阵乘法，并调用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

Answer 2

如果您接受使用GSL，最简单的方法可能是使用这个GSL->；cython接口https://github.com/twiecki/CythonGSL，然后从那里调用BLAS（请参见示例https://github.com/twiecki/CythonGSL/blob/master/examples/blas2.pyx）。它还应该负责Fortran与C的排序。 GSL的新特性并不多，但您可以放心地假设它是主动维护的。与东京相比，CythonGSL更完整；例如，它具有numpy中没有的对称矩阵产品。

Answer 3

调用与Scipy绑定的BLAS是“相当”简单的，这里有一个调用DGEMM来计算矩阵乘法的例子：https://gist.github.com/pv/5437087注意BLAS和LAPACK希望所有数组都是Fortran连续的（lda/b/c参数的模），因此order="F"和double[::1,:]是正确工作所必需的。

通过对恒等矩阵应用LAPACK函数dgesv，可以类似地计算逆矩阵。有关签名，请参见here。所有这些都需要降低到较低级别的编码，您需要自己分配临时工作数组等等——但是这些可以封装到您自己的方便函数中，或者只需使用从tokyo获得的函数指针替换lib_*函数，就可以重用来自tokyo的代码。

如果使用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。