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)

Answer 1

先计算出Cholesky分解，然后再通过回代的方法来求解下三角矩阵的逆，这样应该比直接用linalg.cholesky(linalg.inv(S))要快。

Answer 2

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)

当然，你可以使用切片来进行乘法运算，但我不确定这样做是否会比我建议的方法更快。总之，希望你能明白这个思路。我觉得要进一步加速这个过程，可能没有太多其他的方法了。