python - 加速 scipy 中集成的函数评估

Question

我正在尝试将代码从 Matlab 移植到 SciPy。这是我迄今为止编写的代码的简化版本:https://gist.github.com/atmo/01b6e007be9ef90e402c 。然而，Python 版本比 Matlab 慢得多。我在要点中包含了分析结果，它们表明 python 几乎 90% 的时间都花在评估函数 f 上。除了用 C 或 Cython 重写之外，还有什么方法可以加快其评估速度？

Answer 1

正如我在评论中提到的，如果您考虑到 quad(以及复杂的函数 f)，您可以消除大约一半的调用em>矩阵是对称的。

通过重写这个复杂的函数，可以进一步提高速度，仍然是在纯Python中。我做的大部分都是出于同情。

最后，我尝试使用 np.vectorize 对 quad 的调用进行矢量化。

from scipy.integrate import quad
from scipy.special import jn as besselj
from scipy import exp, zeros, linspace
from scipy.linalg import norm
import numpy as np

def complicated_func(lmbd, a, n, k):
    u,v,w = 5, 3, 2
    x = a*lmbd
    fac = exp(2*x)
    comm = (2*w + x)
    part1 = ((v**2 + 4*w*(w + 2*x) + 2*x*(x - 1))*fac**5
        + 2*u*fac**4
        + (-v**2 - 4*(w*(3*w + 4*x + 1) + x*(x-2)) + 1)*fac**3
        + (-8*(w + x) + 2)*fac**2
        + (2*comm*(comm + 1) - 1)*fac)
    return part1/lmbd *besselj(n+1, lmbd) * besselj(k+1, lmbd)

def perform_quad(n, k, a):
    return quad(complicated_func, 0, np.inf, args=(a,n,k))[0]

def improved_main():
    sz = 20
    amatrix = np.zeros((sz,sz))
    ls = -np.linspace(1, 10, 20)/2
    inds = np.tril_indices(sz)
    myv3 = np.vectorize(perform_quad)
    res = myv3(inds[0], inds[1], ls.reshape(-1,1))
    results = np.empty(res.shape[0])
    for rowind, row in enumerate(res):
        amatrix[inds] = row
        symm_matrix = amatrix + amatrix.T - np.diag(amatrix.diagonal())
        results[rowind] = norm(symm_matrix)
    return results

计时结果显示速度提高了 5 倍(如果我只运行一次，请原谅我，因为它需要足够长的时间):

In [11]: %timeit -n1 -r1 improved_main()
1 loops, best of 1: 6.92 s per loop

In [12]: %timeit -n1 -r1 main()
1 loops, best of 1: 35.9 s per loop

如果您立即将 v 替换为其平方，还会获得微量增益，因为这是它唯一一次在该复杂函数中使用:作为其平方。

对 besselj 的调用也存在大量重复，但我不知道如何避免这种情况，因为 quad 将确定 lmbd，因此您无法轻松地预先计算这些值，然后执行查找。

如果您分析 improved_main，您会发现对 complicated_func 的调用量几乎减少了 2 倍(仍需要计算对角线) )。所有其他速度提升都可以归因于 np.vectorize 以及对 complicated_func 的改进。

我的系统上没有Matlab，所以如果你改进那里的复杂功能，我无法对它的速度增益做出任何声明。