Cython中可能存在的优化：NumPy数组

Question

以下是我的Cython代码，用于从多元正态分布中抽取样本。我使用循环，因为每次密度都不同。（conLSigma是Cholesky分解的结果）

由于每次循环都要进行逆矩阵和Cholesky分解，所以这需要很长时间。虽然比纯Python代码快，但我想知道是否有更好的方法来提高速度。

from __future__ import division

import numpy as np 

cimport numpy as np 

ctypedef np.float64_t dtype_t

cimport cython
@cython.boundscheck(False)
@cython.wraparound(False)

def drawMetro(np.ndarray[dtype_t, ndim = 2] beta,
              np.ndarray[dtype_t, ndim = 3] H,
              np.ndarray[dtype_t, ndim = 2] Sigma,
              float s):

    cdef int ncons = betas.shape[0]
    cdef int nX = betas.shape[1]
    cdef int con

    cdef np.ndarray betas_cand = np.zeros([ncons, nX], dtype = np.float64)
    cdef np.ndarray conLSigma = np.zeros([nX, nX], dtype = np.float64)

    for con in xrange(ncons):
        conLSigma = np.linalg.cholesky(np.linalg.inv(H[con] + Sigma))
        betas_cand[con] = betas[con] + s * np.dot(conLSigma, np.random.standard_normal(size = nX))

    return(betas_cand)

- joon

Answer 1

Cholesky分解创建一个下三角矩阵。这意味着在np.dot中进行的乘法运算接近一半都不需要执行。如果您更改该行

betas_cand[con] = betas[con] + s * np.dot(conLSigma, np.random.standard_normal(size = nX))

转换为

tmp = np.random.standard_normal(size = nX)
for i in xrange(nX):
    for j in xrange(i+1):
        betas_cand[con,i] += s * conLSigma[i,j] * tmp[j]

然而，您还需要进行更改。

cdef np.ndarray betas_cand = np.zeros([ncons, nX], dtype = np.float64)

转换为

cdef np.ndarray betas_cand = np.array(betas)

当然，你可以为乘法使用切片，但我不确定它是否比我建议的方法更快。无论如何，希望你能理解。我认为除了这些，你几乎没有其他办法来加速。

- Justin Peel

---

谢谢你的建议。我在想，是否可以直接调用C库（如BLAS），而不是使用np.dot等函数，这样会更有帮助。 - joon

Answer 2

先计算Cholesky分解，然后通过回代法求解下三角矩阵的逆矩阵。这比使用linalg.cholesky(linalg.inv(S))更快。

- Axel