如何加速多变量复杂函数积分？Python计算效率优化需求

1. 降低不必要的积分精度要求 #

原代码设置的epsrel=1e-20和epsabs=1e-20属于极端过高的精度要求，数值积分根本无法达到该精度，反而会迫使积分器反复迭代浪费算力。根据实际需求调整到合理范围（如1e-6或1e-8），同时缩小无效积分区间（比如M1的后验分布集中在4-16，无需从0积分到30），能大幅减少计算时间。

修改示例：

# 替换精度设置并缩小M1积分区间
options={'epsrel':1e-6, 'epsabs':1e-6}
K, errorK = integrate.nquad(G, ranges=[[4, 16],[0, 10],[0, 10],[0, 10]], opts=[options, options, options, options])

2. 向量化G函数，消除循环与嵌套函数开销 #

原G函数内部嵌套多个小函数，且用列表推导循环计算每个样本的密度值，完全可以用numpy向量化操作替代，利用广播机制一次性完成所有样本计算，大幅提升效率。

修改后的G函数：

def G(M1, Mn, S1, Sn):
    n = N
    i = np.array(I[1:])
    x = np.array(XX[1:])
    
    # 向量化计算loc和scale
    loc = M1*(n - i)/(n-1) + Mn*(i - 1)/(n-1)
    scale = S1*(n - i)/(n-1) + Sn*(i - 1)/(n-1)
    
    # 向量化计算正态密度
    y = (x - loc)/scale
    fnorm = np.exp(-y**2/2) / np.sqrt(2*np.pi)
    F_vals = fnorm / scale
    
    # 直接计算乘积
    return np.prod(F_vals)

3. 缓存重复的积分计算 #

原代码中pdf_m1和F1执行相同的三重积分，每次调用都重新计算，浪费大量算力。可以用functools.lru_cache缓存结果，避免重复计算。

修改示例：

from functools import lru_cache

# 缓存要求参数为可哈希类型，将M1转为float类型
@lru_cache(maxsize=None)
def F1(M1):
    return integrate.tplquad(pdf, 0, 10, 0, 10, 0, 10, args=(float(M1),), epsabs=1e-6, epsrel=1e-6)[0]

# 复用F1的缓存结果生成pdf_m1
def pdf_m1(M1_list):
    return [F1(float(m)) for m in M1_list]

4. 替换高维积分方法：蒙特卡洛积分替代自适应积分 #

对于4维及以上的积分，scipy的自适应积分（nquad/tplquad）效率极低，蒙特卡洛积分更适合高维场景，误差随样本数的平方根降低，计算速度提升几个数量级。

示例：自行实现蒙特卡洛计算4维积分K

def monte_carlo_G(num_samples=500_000):
    # 在积分区间内生成随机样本
    M1_samples = np.random.uniform(4, 16, num_samples)
    Mn_samples = np.random.uniform(0, 10, num_samples)
    S1_samples = np.random.uniform(0, 10, num_samples)
    Sn_samples = np.random.uniform(0, 10, num_samples)
    
    # 批量计算G的值
    vals = np.array([G(M1, Mn, S1, Sn) for M1, Mn, S1, Sn in zip(M1_samples, Mn_samples, S1_samples, Sn_samples)])
    
    # 蒙特卡洛积分公式：区间体积 × 样本均值
    volume = (16-4)*(10-0)*(10-0)*(10-0)
    return volume * np.mean(vals), volume * np.std(vals)/np.sqrt(num_samples)

# 调用，样本数可根据精度需求调整
K, errorK = monte_carlo_G(500_000)

5. 使用Numba JIT编译加速函数计算 #

用numba对G函数进行JIT编译，能将numpy操作的速度提升至接近C语言的水平，进一步减少计算耗时。

示例：

from numba import jit

@jit(nopython=True)
def G_numba(M1, Mn, S1, Sn):
    n = 20
    i = np.array([1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20])
    x = np.array([11.3, 14.8, 7.6, 10.5, 12.7, 3.9, 11.2, 5.4, 8.5, 5.0, 4.4, 7.3, 2.9, 5.7, 6.2, 7.3, 3.3, 4.2, 5.5, 4.2])
    
    loc = M1*(n - i)/(n-1) + Mn*(i - 1)/(n-1)
    scale = S1*(n - i)/(n-1) + Sn*(i - 1)/(n-1)
    
    y = (x - loc)/scale
    fnorm = np.exp(-y**2/2) / np.sqrt(2*np.pi)
    F_vals = fnorm / scale
    
    return np.prod(F_vals)

6. 并行化积分计算 #

scipy的积分器支持并行计算，通过workers参数启用多核心运算，充分利用CPU资源缩短积分时间。

修改示例：

# nquad启用并行，workers=-1表示使用所有可用核心
K, errorK = integrate.nquad(G, ranges=[[4,16],[0,10],[0,10],[0,10]], opts=options, workers=-1)

# tplquad同样启用并行
@lru_cache(maxsize=None)
def F1(M1):
    return integrate.tplquad(pdf, 0, 10, 0, 10, 0, 10, args=(float(M1),), epsabs=1e-6, epsrel=1e-6, workers=-1)[0]

内容的提问来源于stack exchange，提问作者Alexey_Proz