25

我正在尝试使用 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，没有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 调用它只会导致大部分时间花在 numpy 调用的优化 C 代码中。但是，就我而言，我对小矩阵上的这些线性代数运算有很多调用——在这种情况下，重复从 Cython 回到 numpy 再回到 Cython 所引入的开销将远远超过实际计算运算所花费的时间布拉斯。因此，对于这些简单的操作，我想将所有内容都保留在 C/Cython 级别，而不是通过 python。

我不想通过 GSL，因为这增加了另一个依赖项，而且还不清楚 GSL 是否得到积极维护。由于我假设代码的用户已经安装了 scipy/numpy，我可以放心地假设他们拥有与这些库一起使用的所有相关 C 代码，所以我只想能够利用该代码并调用它来自赛通。

编辑：我在 Cython ( https://github.com/tokyo/tokyo ) 中找到了一个包含 BLAS 的库，它很接近但不是我想要的。我想直接调用 numpy/scipy C 函数（我假设用户已经安装了这些函数。）

4

3 回答 3

25

调用与 Scipy 捆绑在一起的 BLAS “相当”简单，这是调用 DGEMM 来计算矩阵乘法的一个示例：https ://gist.github.com/pv/5437087 请 注意，BLAS 和 LAPACK 期望所有数组都是 Fortran 连续的（模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 XXXXCython 访问（虽然它们是私有的，所以不要这样做）。

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

我认为目前没有任何健全且可移植的方式来访问执行繁重工作的函数gamln，（尽管您可以窥探 UFunc 对象，但这不是一个健全的解决方案:)，所以目前可能最好从 scipy.special 中获取源代码的相关部分并将其与您的项目捆绑在一起，或者使用例如 GSL。

于 2013-04-22T18:16:24.380 回答
1

如果您接受使用 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 中不存在的对称矩阵乘积。

于 2014-03-15T20:11:10.450 回答
0

因为我刚刚遇到了同样的问题，并且写了一些额外的函数，所以我会把它们包括在这里，以防其他人发现它们有用。我编写了一些矩阵乘法，还调用了 LAPACK 函数来进行矩阵求逆、行列式和 Cholesky 分解。但是你应该考虑尝试在任何循环之外做线性代数的东西，如果你有的话，就像我在这里做的那样。顺便说一句，如果你有建议，这里的行列式函数不太有效。另外，请注意，我不会检查输入是否符合要求。

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
于 2018-04-30T01:46:49.217 回答