使用内部数组定义提高 Cython Lapack 性能？

Question

我正在尝试加快我正在构建的模型中的一些矩阵求逆，特别是通过在 Cython 中实现一些线性代数例程。我有代码工作，但正在尝试优化。目前我的 Cython 看起来像这样：

import cython
import numpy as np
cimport numpy as cnp
cimport scipy.linalg.cython_lapack as cython_lapack

cdef int ZERO = 0
ctypedef cnp.float64_t REAL_t

def inverse(mat,identity,pivots,k):
    cdef int K = k
    cdef REAL_t* mat_pointer = <REAL_t *>cnp.PyArray_DATA(mat)
    cdef REAL_t* iden_pointer = <REAL_t *>cnp.PyArray_DATA(identity)
    cdef int* piv_pointer = <int *>cnp.PyArray_DATA(pivots)
    cython_lapack.dgesv(&K,&K,mat_pointer,&K,piv_pointer,iden_pointer,&K,&ZERO)

一旦我编译（如 linalg.pyx），我可以运行这很好，但是因为 dgesv 修改了一些输入变量，我必须使用包装函数和做一些复制：

import linalg
import numpy as np

identity = np.eye(K) # K is the dimensionality of the square matrix
def lapack_inverse(a):
    b = identity.copy()
    linalg2.inverse(a,b,np.zeros(K,dtype=np.intc),K)
    return b

对于小 K，我们看到相当可靠的速度改进 (2-3x)：

In [29]: K=10
In [30]: x = np.random.random((K,K))
In [31]: identity = np.eye(K)
In [32]: print np.allclose(np.linalg.inv(x),lapack_inverse(x))
True
In [33]: %timeit np.linalg.inv(x)
10000 loops, best of 3: 28 µs per loop
In [34]: %timeit lapack_inverse(x)
100000 loops, best of 3: 10.3 µs per loop

但对于较大的 K，returns 明显减小（例如对于 K=200）：

In [8]: %timeit np.linalg.inv(x)
100 loops, best of 3: 10.1 ms per loop
In [9]: %timeit lapack_inverse(x)
100 loops, best of 3: 8.1 ms per loop

我想我可以通过在 C 中定义标识和主元数组来获得更好的性能，而不是像我目前正在做的那样进行复制和传递。换句话说，我想修改 Cython 函数，让它只需要我传递要反转的数组。

我的问题是 (a)：这会导致速度提高（即使是很小的提高）吗？并且，假设它会，(b) 我如何着手直接在 Cython 中定义枢轴和恒等数组？

Answer 1

在 Cython 中有很多分配内存的方法——基本上你可以创建任何 Python 对象或使用标准的 C 分配技术。 This previous question 列出了很多（答案中给出了时间）并且几乎肯定比我的答案更有用。不过，它处理的用例略有不同。

在您的情况下，您必须将单位矩阵分配为一个 numpy 数组，因为答案已写入其中并且您想从函数中 return （你显然可以分配它，否则复制，但这似乎有点毫无意义）。我怀疑最好是调用 eye 而不是每次都复制它。

您可以随心所欲地分配枢轴。在下面的示例中，我使用 C 分配机制作为分配大块内存的最直接方式（因此希望最快）。这样做的缺点是您必须使用 free.

解除分配它

cimport numpy as cnp
import numpy as np
#cimport scipy.linalg.cython_lapack as cython_lapack
from libc.stdlib cimport malloc, free

ctypedef cnp.float64_t REAL_t

def inverse(mat):
    # you should probably check that mat.shape[0]==mat.shape[1]
    # and that mat is actually a float64 array
    cdef int k = mat.shape[0]
    cdef int info = 0

    cdef REAL_t* mat_pointer = <REAL_t*>cnp.PyArray_DATA(mat)
    # I suspect you should be doing a check that mat_pointer has been assigned
    cdef REAL_t* iden_pointer

    cdef int* piv_pointer = <int*>malloc(sizeof(int)*k)
                # this is uninitialised (the contents are arbitrary)
                # but that's OK because it's used as an output
    try:
        identity = np.eye(k)

        iden_pointer = <REAL_t*>cnp.PyArray_DATA(identity)

        cython_lapack.dgesv(&k,&k,mat_pointer,&K,
                     piv_pointer,iden_pointer,&K,&info)
        # you should check info to ensure it's worked
        return identity
    finally:
        free(piv_pointer) # the "try ... finally" ensures that this is freed

有几个地方我注意到您应该添加检查以确保一切都有效...这会减慢速度但会更安全。您可能还可以添加一些 Cython 类型定义以获得额外的速度，但我还没有真正研究这个例子。

最后 - 想一想为什么要使用逆函数。如果您只是将另一个矩阵乘以逆矩阵，那么直接对该矩阵使用 dgsev 会更好（更快更准确）。很少有充分的理由直接计算逆。

使用内部数组定义提高 Cython Lapack 性能？

Improving Cython Lapack performance with internal array definitions?

python

linear-algebra

cython

lapack