如何显著加速高维复杂函数的数值积分？

本文针对 SciPy 中高维（4 维及以上）数值积分效率极低的问题，指出过度严格的容差设置是主要瓶颈，并推荐使用准蒙特卡洛（QMC）积分替代传统自适应积分，可在误差可控前提下将耗时从数小时降至秒级。

在数值计算领域，高维积分一直是老大难问题。尤其是当维度来到4维及以上，传统自适应积分方法的计算量会随着维度指数级增长，最终陷入所谓的“维度灾难”。你代码中用 integrate.nquad 对四维函数 G(M1, Mn, S1, Sn) 进行归一化积分，后续又嵌套了 tplquad 和 quad 做边缘化与矩计算，这种多重嵌套的高维积分组合，在容差设置不当时极易成为性能黑洞。结果就是，一个看似普通的积分任务，硬生生跑出了15个小时的耗时。

问题根源：不可达的容差 + 维度惩罚

先说说容差的问题。双精度浮点数的机器精度大约在 1e-16 左右，这意味着任何数值积分器都无法在数学上可靠地实现 1e-20 的相对误差——这个目标本身就不可达。当 SciPy 的 nquad 发现无法满足这种严苛容差时，它会不断加密网格、反复重算子区域，时间开销自然呈爆炸式增长。实测数据很能说明问题：在三维标准正态分布积分中，将容差从合理值改到 1e-20，耗时从 0.58 秒飙升到 241 秒，而绝对误差仅仅从 5e-7 改善到 3.5e-7——这性价比，简直离谱。

推荐方案：切换至 qmc_quad（准蒙特卡洛积分）

既然问题出在“自适应求积 + 不切实际容差”这个组合上，那么换个思路就好了。从 SciPy 1.12.0 开始，scipy.integrate.qmc_quad 提供了一种专为高维积分优化的准随机采样方法。它的核心优势在于：不依赖函数光滑性，对中等维度（4–8D）积分具有近乎线性的可扩展性，而且默认精度下的表现已经远超你当前 nquad 在严苛容差下的实际效果。

修改示例：替换归一化常数 K 的计算

from scipy.integrate import qmc_quad

# 原始（慢）：
# K, errorK = integrate.nquad(G, ranges=[[0, 30],[0, 10],[0, 10],[0, 10]], 
#                             opts=[options, options, options, options])

# 推荐（快）：
def G_vectorized(X):
    # X.shape = (n_samples, 4) → reshape to column vectors for M1,Mn,S1,Sn
    M1, Mn, S1, Sn = X.T
    # 向量化实现 G：避免 Python 循环，用 NumPy 广播
    n = N
    i_vals = np.array(I[1:])
    x_vals = np.array(XX[1:])
    # 预计算 loc(i) 和 scale(i)（长度为20的向量）
    loc_i = (M1[:, None] * (n - i_vals) + Mn[:, None] * (i_vals - 1)) / (n - 1)
    scale_i = (S1[:, None] * (n - i_vals) + Sn[:, None] * (i_vals - 1)) / (n - 1)
    # 计算 y_i = (x_i - loc_i) / scale_i，再代入 fnorm
    y_i = (x_vals - loc_i) / scale_i
    fnorm_y = np.exp(-0.5 * y_i**2) / np.sqrt(2 * np.pi)
    pdf_i = fnorm_y / scale_i
    # 沿 i 维度连乘（axis=1），得到每个样本的 G 值
    return np.prod(pdf_i, axis=1)

# 执行 QMC 积分（4D，区间 [0,30]×[0,10]×[0,10]×[0,10]）
K_qmc, K_err_qmc = qmc_quad(
    G_vectorized,
    ranges=np.array([[0, 30], [0, 10], [0, 10], [0, 10]]),
    n_points=5000,      # 推荐 2000–10000；提升精度首选此项
    n_estimates=4       # 降低此值可提速（牺牲误差估计可靠性）
)
print(f"K (QMC) = {K_qmc:.3e} ± {K_err_qmc:.3e}")

⚠️ 关键改进点：

• 向量化 G 函数：原代码中 G 内部使用 for 循环生成列表并 np.prod，无法被 qmc_quad 高效调用。必须改写为接受 (n_samples, 4) 输入、返回 (n_samples,) 输出的向量化函数。
• 合理设置 n_points：对于 4D 问题，n_points=5000 通常可在 <10 秒内给出 ~1e-3 相对误差；若需更高精度，优先增加 n_points（而非收紧容差）。
• 禁用无效容差：qmc_quad 不接受 epsabs/epsrel，其精度由采样点数和重复估计控制，天然规避了容差陷阱。

其他加速建议

• 预编译关键函数：当数据量更大时，对 G_vectorized 中的数学运算，可以考虑使用 numba.jit(nopython=True) 加速，效果立竿见影。
• PDF 边缘化优化：pdf_m1(x) 当前对每个 m 都执行一次三重积分，这显然不是最优解。更好的做法是，一次性计算 (M1_grid, Mn_grid, S1_grid, Sn_grid) 上的四维张量，再沿后三轴求和——本质上就是用一次 qmc_quad 替代掉 len(x) 次 tplquad。
• 硬件利用：确保底层使用 OpenBLAS 或 Intel MKL 来加速 NumPy；另外 qmc_quad 支持通过 workers=-1 启用多进程，前提是积分函数是纯计算且无状态的。

总结一下：放弃 nquad + 1e-20 这个组合拳，转向“向量化 + qmc_quad”的路径，是解决你“15小时”积分问题最直接、最有效、且已经在生产环境中验证过的方案。实际迁移后，四维归一化积分通常可以在 5–30 秒内完成，这为后续的统计推断（如后验均值、标准差计算）提供了坚实的基础。