cython 中可能的优化：numpy 数组答案

【问题标题】：Possible optimization in cython: numpy arraycython 中可能的优化：numpy 数组
【发布时间】：2010-11-21 05:22:55
【问题描述】：

以下是我从多元正态分布中绘制的 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)

【问题讨论】：

标签： python numpy cython

【解决方案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)

您当然可以使用切片进行乘法运算，但我不确定它是否会比我建议的方式更快。无论如何，希望你能明白这一点。我不认为你可以做很多其他事情来加快速度。

【讨论】：

感谢您的提示。我想知道如果不使用 np.dot 的东西，直接调用诸如 BLAS 之类的 C 库会有所帮助。

【解决方案2】：

先计算cholesky分解，然后通过反向替换来反转下三角矩阵。这应该比 linalg.cholesky(linalg.inv(S)) 更快。

【讨论】：