万倍加速——高性能正态积分的终极方案附代码
发表更新数据科学 / 量化交易

万倍加速——高性能正态积分的终极方案附代码

在毫秒必争的衍生品定价与大规模风险回测中，传统的正态分布积分运算往往是代码毒瘤和性能瓶颈。本文将介绍一种完全向量化的数值算法，助你利用现代 CPU、GPU 甚至TPU 的并行能力突破卡点。

本文适用于性能敏感型场景，而非极致精度敏感型场景

导言

在量化金融的工程实践中，标准正态分布的累积分布函数是衍生品定价、风险度量以及各种随机微分方程数值解的核心组件。在 Python 生态中，绝大多数 Research 阶段的代码会直接调用 scipy.stats.norm.cdf 或 scipy.stats.multivariate_normal。

然而，这种“标准做法”存在严重的性能瓶颈：

  1. 调用开销与 GIL 限制：scipy 的底层虽然是 C/Fortran，但每次调用的 Python wrapper 开销不可忽视。特别是在很多非纯矩阵运算的场景中，循环内的 Python 函数调用会成为致命拖累。
  2. 缺乏 SIMD 向量化支持：标准的库函数往往难以被编译器自动优化为 AVX-512 或 Neon 指令集。当我们需要在大规模张量上并发计算时，黑盒函数阻断了流水线优化的可能。
  3. 框架移植性差：现代高性能交易系统或回测框架往往基于 C++、Rust 或 CUDA，甚至是在 Python 内部使用 Numba/JAX 进行 JIT 编译。引入一个沉重的 scipy 依赖导致在某些 JIT 环境（如 Numba nopython 模式）下根本无法通过。

因此，我们需要一种 基于纯初等数学运算（加减乘除、乘方、指数）的近似算法。它们不依赖复杂的外部库，天然适合向量化（SIMD），且极易移植到 GPU Kernel 或 FPGA 逻辑中，同时在精度上能够满足需求。

一元正态分布

对于单变量标准正态分布 \(N(0,1)\)，最经典的数值算法莫过于利用多项式拟合来逼近误差函数，从而规避掉昂贵的积分运算，仅保留一个指数运算 \(e^{-x^2/2}\) 和若干次乘加运算。

由于正态分布关于 \(y\) 轴对称，\(\Phi(-x) = 1 - \Phi(x)\)，我们只需要计算 \(x \ge 0\) 的情况。

对于 \(x \ge 0\)，累积分布函数的近似计算公式如下：

$$ \Phi(x) \approx 1 - Z(x) \left( b_1 t + b_2 t^2 + b_3 t^3 + b_4 t^4 + b_5 t^5 \right) $$

其中：

\(Z(x) = \frac{1}{\sqrt{2\pi}} e^{-x^2/2}\) 是标准正态分布的概率密度函数（PDF）。

\(t = \frac{1}{1 + p x}\) 是辅助变量。

为了达到金融计算所需的精度（最大绝对误差 \(|\epsilon(x)| < 7.5 \times 10^{-8}\)），我们需要使用以下特定的拟合参数：

\(p = 0.2316419\)

\(b_1 = 0.319381530\)

\(b_2 = -0.356563782\)

\(b_3 = 1.781477937\)

\(b_4 = -1.821255978\)

\(b_5 = 1.330274429\)

整个过程仅包含 1 次指数运算、1 次倒数运算和 5 次乘法累加，计算复杂度极低。\(10^{-7}\) 的误差对于绝大多数期权定价（甚至包括希腊字母的差分计算）已经足够。这种算法在底层可以被编译器轻松优化为 FMA（Fused Multiply-Add）指令，是在 C++ 或 Numba 中手写 NormCDF 的首选方案。

二元正态分布

对于标准二元正态分布函数 \(\Phi_2(x, y; \rho)\)，其计算公式为：

$$ \Phi_2(x, y; \rho) = \frac{1}{2\pi\sqrt{1-\rho^2}} \int_{-\infty}^{x} \int_{-\infty}^{y} \exp\left( -\frac{u^2 - 2\rho u v + v^2}{2(1-\rho^2)} \right) du \, dv $$

其中 \(\rho\) 是相关系数，满足 \(-1 < \rho < 1\)，

要解决二重积分 \(\iint\) 的计算困难，最优雅的路径是利用 Plackett 恒等式 (1954) 。这是一个非常深刻的数学发现：二元正态分布对相关系数 \(\rho\) 的偏导数，竟然恰好等于其概率密度函数。

$$\frac{\partial \Phi_2(x,y;\rho)}{\partial \rho} = \phi_2(x,y;\rho)$$

利用微积分基本定理，我们可以将 \(\Phi_2\) 拆解为“起点”加上“变化量的积分”。

  • 起点 (\(\rho=0\))：当资产不相关时，二元分布退化为两个一元分布的乘积，即 \(\Phi(x)\Phi(y)\)。
  • 过程 (\(\rho\))：从 0 到 \(\rho\) 的累积影响。

由此得到 Drezner-Wesolowsky 的核心积分公式：

$$\Phi_2(x,y;\rho) = \Phi(x)\Phi(y) + \frac{1}{2\pi} \int_0^{\rho} \frac{1}{\sqrt{1-t^2}} \exp\left( -\frac{x^2 - 2txy + y^2}{2(1-t^2)} \right) dt$$

  1. 降维：原本关于 \(dxdy\) 的二维积分，被神奇地转换成了关于相关系数 \(t\) 的一维积分。

  2. 平滑性：被积函数 \(g(t; x, y)\) 在区间 \([0, \rho]\) 内是非常平滑的解析函数。这意味着我们不需要大量的采样点就能获得极高的精度。

  3. 变量角色转换：在积分内部，\(x\) 和 \(y\) 变成了固定的“参数”，不再参与积分运算。

公式虽然变简单了，但工程上还有一个大坑：积分上限是 \(\rho\)。对于量化场景，我们有 \(N\) 个期权，每个期权的 \(\rho_i\) 都不一样。这意味着每个样本的积分路径长度不同。为了向量化，我们必须将所有参差不齐的积分区间 \([0, \rho_i]\) 映射到标准区间 \([-1, 1]\) 上。

令积分变量 \(t\) 与标准变量 \(u\) 的关系为：

$$t(u) = \frac{\rho}{2}(u+1), \quad u \in [-1, 1]$$

此时微分元 \(dt\) 变为：

$$dt = \frac{\rho}{2} du$$

代入上述换元，积分项变为：

$$\int_0^\rho g(t) dt = \frac{\rho}{2} \int_{-1}^1 g\left( \frac{\rho}{2}(u+1) \right) du$$

根据 Gauss-Legendre 理论，对于 \([-1, 1]\) 上的积分，可以近似为 \(K\) 个固定节点 \(u_j\) 和权重 \(w_j\) 的加权和：

$$\text{Integral} \approx \frac{\rho}{2} \sum_{j=1}^{K} w_j \cdot g\left( \frac{\rho}{2}(u_j+1); x, y \right)$$

无论 \(\rho\) 是 0.1 还是 0.9，我们永远只需要计算那固定的 \(K\) 个点（例如 \(K=10\)）。\(u_j\) 和 \(w_j\) 是数学常数，可以预计算。

代码与测试

以下是复制即可使用的jax实现，能够在cpu/gpu/tpu上高效运行。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
import jax
import jax.numpy as jnp

jax.config.update("jax_enable_x64", True)

GL_NODES = jnp.array([
    -0.9739065285171717, -0.8650633666889845, -0.6794095682990244, 
    -0.4333953941292472, -0.1488743389816312,  0.1488743389816312,
    0.4333953941292472,  0.6794095682990244,  0.8650633666889845, 
    0.9739065285171717
])

GL_WEIGHTS = jnp.array([
    0.0666713443086881, 0.1494513491505806, 0.2190863625159820,
    0.2692667193099963, 0.2955242247147529, 0.2955242247147529,
    0.2692667193099963, 0.2190863625159820, 0.1494513491505806,
    0.0666713443086881
])

@jax.jit
def fast_norm_cdf(x):
    p  = 0.2316419
    b1 = 0.319381530
    b2 = -0.356563782
    b3 = 1.781477937
    b4 = -1.821255978
    b5 = 1.330274429
    
    abs_x = jnp.abs(x)
    z = (1.0 / jnp.sqrt(2 * jnp.pi)) * jnp.exp(-0.5 * abs_x**2)
    t = 1.0 / (1.0 + p * abs_x)
    
    poly = t * (b1 + t * (b2 + t * (b3 + t * (b4 + t * b5))))
    y = 1.0 - z * poly
    
    return jnp.where(x >= 0, y, 1.0 - y)

@jax.jit
def fast_bivariate_norm_cdf(x, y, rho):
    base_prob = fast_norm_cdf(x) * fast_norm_cdf(y)
    rho_clipped = jnp.clip(rho, -0.9999, 0.9999)
    
    x_col = x[:, None]
    y_col = y[:, None]
    rho_col = rho_clipped[:, None]
    nodes_row = GL_NODES[None, :]
    weights_row = GL_WEIGHTS[None, :]
    
    t_grid = rho_col / 2.0 * (nodes_row + 1)
    
    one_minus_t2 = 1.0 - t_grid**2
    numerator = x_col**2 - 2 * t_grid * x_col * y_col + y_col**2
    denominator = 2 * one_minus_t2
    
    log_integrand = -jnp.log(2 * jnp.pi) - 0.5 * jnp.log(one_minus_t2) - (numerator / denominator)
    integrand = jnp.exp(log_integrand)
    
    integral_sum = jnp.sum(weights_row * integrand, axis=1)
    
    return base_prob + (rho_clipped / 2.0) * integral_sum

GPU Benchmark（Tesla T4，CPU: Xeon 7B13）

Samples (N) Task Method Time (s) Speedup
10,000 Univariate JAX 0.15163 0.0×
10,000 Bivariate JAX 0.45608 12.1×
100,000 Univariate JAX 0.26846 0.0×
100,000 Bivariate JAX 0.46016 153.1× (Est.)
1,000,000 Univariate JAX 0.15443 0.4×
1,000,000 Bivariate JAX 0.27804 1139.2× (Est.)
10,000,000 Univariate JAX 0.08185 5.3×
10,000,000 Bivariate JAX 0.18332 12730.8× (Est.)
100,000,000 Univariate JAX 0.12218 36.6×
100,000,000 Bivariate JAX 0.22598 101728.1× (Est.)
  • 对于二元正态分布积分，随着样本规模增大，GPU 并行度被充分释放，速度提升呈指数级放大，在 \(10^8\) 规模下达到 \(10^5\times\) 级别加速（估算）。
  • 一元 CDF 在小批量下受 kernel launch 与 JIT 开销影响明显，而在超大 batch 场景中开始体现 GPU 吞吐优势。

benchmark

万倍加速——高性能正态积分的终极方案附代码

https://heth.ink/FastNorm/

作者

YK

发布于

2026-01-22

更新于

2026-01-22

许可协议

目录

分类

最新文章

×