Ziggurat算法：高效生成正态分布随机数的原理与实现

1. 从“慢”到“快”：为什么我们需要Ziggurat算法？

如果你用过MATLAB里的randn函数，或者在任何需要生成正态分布随机数的场景里写过代码，你可能觉得这很简单，一行命令就搞定了。但如果你真的去深究过这行命令背后的实现，或者在一个对性能极其敏感的场景（比如高频交易策略回测、大规模蒙特卡洛模拟，或者实时物理引擎）里调用过百万、千万次randn，你很快就会意识到，这个“简单”的操作，可能正是你整个程序里最耗时的瓶颈之一。

传统的正态随机数生成方法，比如经典的Box-Muller变换或者反函数法，在数学上很优雅，但计算上并不“便宜”。Box-Muller需要计算三角函数（sin,cos）和对数（log），这些都是相对昂贵的浮点运算。反函数法则需要数值求解误差函数erf的反函数，计算量更大。当你的循环需要生成数以亿计的随机数时，这些开销累积起来会非常可观。

这就引出了我们今天要拆解的核心：Ziggurat算法。它不是MATLAB的专属，而是一种被广泛应用在众多科学计算库和编程语言（如NumPy的早期版本、GNU Scientific Library等）中的高效正态随机数生成器。它的目标非常明确：在保证统计性质正确的前提下，用尽可能少的计算资源，尤其是尽可能避免那些昂贵的超越函数（如log,sin,exp），来快速产生服从标准正态分布 N(0,1) 的随机数。

我第一次在MATLAB的底层优化文档里注意到它时，感觉这个名字很神秘——“金字形神塔”。后来才明白，这个算法的核心思想，正是用一系列堆叠的矩形（就像古代美索不达米亚的金字形神塔）来覆盖正态分布曲线的右半部分，通过一种巧妙的拒绝采样和分层查找策略，使得绝大多数情况下，只需要几次简单的比较和一次均匀随机数生成，就能得到一个合格的正态随机数。复杂计算被“摊销”到了前期一次性的初始化过程中。理解了它，你不仅能优化自己的代码，更能体会到算法设计中对“计算成本”的极致权衡。

2. Ziggurat算法的核心：用“堆盒子”来替代复杂计算

要理解Ziggurat为什么快，我们得先看看它要解决的问题长什么样。标准正态分布的概率密度函数（PDF）是 $f(x) = \frac{1}{\sqrt{2\pi}} e^{-x^2/2}$。我们的目标是生成服从这个分布的随机数。

最朴素的想法是“反函数法”：先产生一个[0,1)间的均匀随机数U，然后计算 $X = F^{-1}(U)$，其中F是正态分布的累积分布函数（CDF）。但 $F^{-1}$（即probit函数）没有解析解，需要数值近似，计算慢。

Ziggurat算法采用了截然不同的思路：用一系列容易抽样的区域（矩形）去覆盖难以直接抽样的PDF曲线。它主要针对PDF的右半部分（x>=0），因为分布是对称的，生成一个右半部分的数，再随机赋予正负号即可。

2.1 “神塔”的构建：分层与覆盖

想象一下正态分布钟形曲线的右半部分。Ziggurat算法用n个水平放置的矩形去覆盖它，这些矩形一层层向上堆叠，最底层（第0层）最宽，最高层（第n-1层）最窄，形成一个类似金字塔或阶梯状的结构——这就是“Ziggurat”名字的由来。

具体构建过程是这样的：

1. 确定层数n：通常取n=128或256。层数越多，每个矩形越薄，拒绝采样的概率越低，效率越高，但初始化表格也越大。这是一个时空权衡。
2. 确定矩形高度：让所有矩形具有相同的面积A。这是一个关键设定。因为PDF曲线下的总面积（x>=0部分）是0.5，所以 $n * A$ 应该略大于0.5。实际上，我们会寻找一个A，使得最顶层的矩形刚好触碰到PDF曲线，或者以某种方式终止。
3. 计算每层的边界：设第i个矩形（i从0到n-1）覆盖的x轴范围是 $[x_i, x_{i+1}]$，其中 $x_0 = 0$。矩形的高度是 $y_i = f(x_i)$（即PDF在左边界处的值）。由于面积相等，有 $A = (x_{i+1} - x_i) * y_i$。同时，矩形右上角 $(x_{i+1}, y_i)$ 应该落在PDF曲线上方或恰好落在曲线上。通过迭代可以求解出这一系列 $x_i$ 和 $y_i$。

最终，我们得到两个预先计算好的表格：x[i]和y[i]（其中y[i] = f(x[i])），以及面积A。这个初始化过程只需要做一次，之后生成随机数时直接查表。

2.2 抽样“三步走”：查表、比较、决定

当我们需要生成一个随机数时，算法步骤如下：

1. 随机选层：生成一个随机整数 $i$，均匀分布在 $[0, n-1]$。
2. 随机选横坐标：在该层对应的x轴区间 $[x_i, x_{i+1}]$ 内，均匀生成一个随机数 $u_1$。那么点 $(u_1, y_i)$ 就落在了第i个矩形内部。
3. 接受/拒绝判断：
   • 快速接受：如果 $u_1 < x_i$（注意，对于i=0，$x_0=0$，这个条件不成立），那么这个点一定位于PDF曲线下方（因为矩形左半边完全在曲线下）。此时，直接返回 $u_1$ 作为采样值（别忘了最后随机赋予正负号）。这是最快的情况。
   • 常规判断：否则，点位于矩形的右半部分。我们需要计算PDF曲线在该点的真实值 $f(u_1)$，并与一个在 $[0, y_i]$ 区间内均匀生成的随机数 $u_2 * y_i$ 比较。
     • 如果 $u_2 * y_i < f(u_1)$，则点落在曲线下方，接受 $u_1$ 作为采样值。
     • 否则，拒绝这次采样，回退到步骤1重试（或者进入更复杂的“尾部”处理，见下文）。

这里的关键优化在于：绝大多数情况下，我们都会落入“快速接受”区域。因为对于靠下的层（i较小），$x_i$ 很小，矩形左半边很窄，“快速接受”区域小。但对于靠上的层（i较大），$x_i$ 很大，矩形左半边非常宽，几乎整个矩形都属于“快速接受”区。而随机选层是均匀的，这意味着我们有很大概率选到上层，从而以极低的成本（一次整数随机、一次浮点随机、一次比较）就完成采样。只有当下层被选中，且点落在右半边时，才需要计算相对昂贵的 $f(u_1) = e^{-u_1^2/2}$（可以预先计算 $\frac{1}{\sqrt{2\pi}}$ 的系数）。

2.3 处理“塔尖”和“塔底”

上面的流程还有两个特例需要处理：

• 最底层（i=0）：这一层比较特殊，因为它覆盖了从0到 $x_1$ 的区域，并且曲线在0点最高。这一层没有“快速接受”区。算法通常将最底层进一步细分，或者采用其他更直接的采样方法（比如简单的拒绝采样），因为这块区域面积占比小，对整体效率影响不大。
• 最顶层（i=n-1）与尾部：最顶层的矩形无法完全覆盖PDF曲线无限延伸的“尾部”（$x > x_n$ 的部分）。因此，当随机选层恰好选中最顶层时，我们需要一个专门的“尾部采样”算法。一个经典方法是：令 $x = -\ln(U_1)/x_n$, $y = -\ln(U_2)$，如果 $2y > x^2$，则返回 $x + x_n$，否则拒绝重试。这实际上是在对指数分布的尾部进行采样。虽然涉及对数运算，但触发尾部采样的概率极低（由层数n控制，例如n=128时概率小于0.5%），因此其计算开销被“摊销”到了海量的快速采样中。

通过这种精妙的分层和查表设计，Ziggurat算法将昂贵的指数函数计算概率降到了很低，大部分计算成本只是生成均匀随机数和整数比较，因此速度比Box-Muller等传统方法快得多。

3. 在MATLAB中验证与实现一个简易Ziggurat

MATLAB内置的randn函数在历史上确实使用过Ziggurat算法（具体版本可能随更新而变化）。我们可以通过一些简单的测试来感知其效率，并尝试自己实现一个简化版来加深理解。

3.1 对比MATLAB内置randn的性能

虽然我们无法看到MATLAB的源码，但可以进行一个简单的速度测试：

num_samples = 1e7; % 测试内置randn tic; samples_builtin = randn(num_samples, 1); time_builtin = toc; fprintf('内置 randn 生成 %d 个样本耗时: %.4f 秒\n', num_samples, time_builtin); % 测试基于Box-Muller的方法（作为对比） tic; u1 = rand(num_samples/2, 1); u2 = rand(num_samples/2, 1); % Box-Muller变换 z0 = sqrt(-2 * log(u1)) .* cos(2 * pi * u2); z1 = sqrt(-2 * log(u1)) .* sin(2 * pi * u2); samples_bm = [z0; z1]; time_bm = toc; fprintf('Box-Muller 生成 %d 个样本耗时: %.4f 秒\n', num_samples, time_bm); fprintf('速度比 (Box-Muller / randn): %.2f\n', time_bm / time_builtin);

在我的测试环境中，内置randn通常比纯MATLAB实现的Box-Muller快2到5倍。这其中的优势部分来自于MATLAB底层是编译优化的C/C++代码，但算法本身的效率（Ziggurat vs Box-Muller）是根本原因。

3.2 实现一个简化版Ziggurat生成器

为了真正搞懂，我们来写一个高度简化、用于教学的Ziggurat生成器。这个版本层数较少（比如8层），省略一些优化，但完整呈现算法骨架。

classdef SimpleZiggurat < handle % 一个简化的Ziggurat正态随机数生成器教学示例 properties (Constant, Access = private) % 使用一个较小的层数便于理解 N = 8; % 预先计算好的表 (对于N=8，这些值需要预先算好) % x[i]: 第i层的左边界 % y[i]: f(x[i])，即PDF在左边界的值 X = [0, 0.258, 0.522, 0.800, 1.099, 1.427, 1.797, 2.230, 3.949]; % x[N]是尾部起点 Y = [0.3989, 0.385, 0.350, 0.290, 0.226, 0.160, 0.099, 0.044, 0.000]; % y[N] = f(x[N]) ~ 0 % 每层的矩形面积 A A = 0.138; end properties (Access = private) % 用于生成均匀随机数的状态 rng_state; end methods function obj = SimpleZiggurat(seed) % 构造函数，初始化随机数种子 if nargin < 1 seed = sum(100*clock); end obj.rng_state = seed; end function [x, state] = randn(obj, num) % 生成num个标准正态随机数 x = zeros(num, 1); for k = 1:num [x(k), obj.rng_state] = obj.sampleOne(obj.rng_state); end state = obj.rng_state; end function [sample, new_state] = sampleOne(obj, state) % 核心：生成一个样本 while true % 1. 随机选择一层 (0 到 N-1) [u1, state] = obj.uni_rand(state); i = floor(obj.N * u1); % i in [0, N-1] % 2. 在该层随机生成横坐标 u [u, state] = obj.uni_rand(state); u_scaled = obj.X(i+1) + u * (obj.X(i+2) - obj.X(i+1)); % 3. 快速接受测试 (对于i>0且u_scaled落在左半区) if i > 0 && abs(u_scaled) < obj.X(i+1) sample = u_scaled; break; end % 4. 常规测试或尾部测试 if i == obj.N - 1 % 尾部采样 (简化版，使用拒绝采样) [sample, state] = obj.sampleTail(state); break; else % 常规拒绝采样：比较 f(u_scaled) 和 一个均匀随机数 [v, state] = obj.uni_rand(state); f_u = obj.normpdf(u_scaled); if v * obj.Y(i+1) < f_u sample = u_scaled; break; end % 否则拒绝，继续while循环 end end % 5. 随机赋予正负号 [sign_bit, state] = obj.uni_rand(state); if sign_bit < 0.5 sample = -sample; end new_state = state; end function [sample, new_state] = sampleTail(obj, state) % 尾部采样简化实现 (非标准Ziggurat尾部算法，仅为演示) % 实际Ziggurat使用更高效的 -ln(U)/x_N 方法 while true [u1, state] = obj.uni_rand(state); [u2, state] = obj.uni_rand(state); % 在尾部区域进行简单的拒绝采样 x_tail = obj.X(end) + u1 * 5; % 假设尾部延伸到 x_N+5 f_x = obj.normpdf(x_tail); if u2 * obj.Y(end-1) < f_x % 与倒数第二层的高度比较 sample = x_tail; new_state = state; return; end end end function y = normpdf(~, x) % 标准正态分布概率密度函数 y = exp(-0.5 * x.^2) / sqrt(2*pi); end function [u, new_state] = uni_rand(obj, state) % 一个简单的线性同余生成器 (LCG) 用于演示 % 注意：这不是高质量的RNG，仅用于教学 a = 1664525; c = 1013904223; m = 2^32; state = mod(a * state + c, m); u = double(state) / m; new_state = state; end end end

这个简化版本有很多不完美的地方：层数太少导致拒绝率高，均匀随机数生成器(LCG)质量差，尾部采样实现低效等等。但它清晰地展示了Ziggurat算法的框架：查表（X, Y）、分层随机选择、快速接受判断、常规拒绝判断、尾部特殊处理。你可以运行它，并和randn比较结果的基本统计特性（均值、方差、直方图），会发现虽然慢，但确实能产生正态分布的数据。

注意：生产级别的Ziggurat实现复杂得多。例如，MATLAB或NumPy中的实现会使用更大的查找表（256层）、更高质量的基础均匀随机数生成器（如Mersenne Twister）、经过精心优化的尾部分布采样算法，并且所有判断和计算都针对CPU流水线和缓存做了优化。我们的简化版只是为了理解原理。

4. Ziggurat的优缺点与适用场景

没有一种算法是万能的，Ziggurat也不例外。了解它的边界，才能更好地决定何时使用它，或者何时选择其他方法。

4.1 优势：速度与确定性

• 极高的速度：这是Ziggurat最大的优点。通过用查表和比较代替绝大多数超越函数计算，它的平均采样时间非常短。对于需要海量正态随机数的应用（如蒙特卡洛模拟），这能带来显著的性能提升。
• 可重复性与确定性：和大多数基于查找表的算法一样，Ziggurat的行为是确定性的。给定相同的随机数种子和初始化表格，它总是产生相同的序列。这对于调试和可重复的科学研究非常重要。
• 适用于向量化：算法的核心流程（生成随机整数索引、生成均匀随机数、比较）很容易进行向量化操作，在现代CPU的SIMD指令集下可以获得极高的吞吐量。

4.2 劣势与潜在陷阱

• 初始化开销与内存占用：Ziggurat需要预先计算并存储x[i]和y[i]表格。对于层数n=256，这只需要几千字节的内存，在现代计算机上微不足道。但这也意味着它有一个小小的“启动成本”。如果你只需要生成几十个随机数，这个开销可能显得不划算。
• 周期性与随机性质量：Ziggurat的随机性质量完全依赖于底层均匀随机数生成器（URNG）的质量。它本身只是一个“分布变换器”。如果你使用的URNG周期短或有统计缺陷，那么Ziggurat输出的序列也会有同样的问题。在生产环境中，必须为其配备一个强大的URNG（如Mersenne Twister、PCG或硬件RDRAND）。
• 尾部精度：由于尾部采样是独立且触发概率极低的特殊路径，理论上尾部分布的精度可能和主体部分略有差异。但对于绝大多数实际应用（即使是6-sigma的极端事件），这种差异可以忽略不计。
• 不适用于非常规正态分布：标准的Ziggurat生成的是N(0,1)。如果你需要N(μ, σ²)，只需进行线性变换：X = μ + σ * Z，其中Z来自Ziggurat。这不会影响性能。但如果你需要生成截断正态分布、多元正态分布或其他复杂变体，Ziggurat不能直接使用，可能需要结合其他算法。

4.3 与其他算法的对比选型

如何选择？对于绝大多数应用，直接使用你所用的科学计算库（MATLAB的randn、NumPy的numpy.random.randn、C++的std::normal_distribution）的默认实现就是最佳选择。这些库的维护者已经为你选择了在当前平台上经过充分优化和测试的算法，很可能就是某种变体的Ziggurat或高度优化的Box-Muller/Marsaglia Polar。

只有当你在以下情况时，才需要考虑自己实现或选择特定算法：

1. 极致性能追求：你在编写底层库，且性能分析表明随机数生成是瓶颈，你可能需要针对你的CPU架构（如利用AVX指令集）手写一个向量化的Ziggurat。
2. 特殊环境限制：在内存极其有限（如嵌入式系统）或没有浮点运算单元的环境下，你可能需要寻找不使用查找表或超越函数的近似方法。
3. 可重复性研究：你需要一个算法，其行为在不同编程语言、不同硬件上完全一致（比特级可重复）。这时你可能需要指定一个特定的、公开参考实现的算法。

5. 超越基础：Ziggurat的变体与现代优化

经典的Ziggurat算法由Marsaglia和Tsang在2000年提出。此后，研究者们提出了各种改进和变体，旨在进一步提高速度、降低拒绝率或简化实现。

5.1 改进的拒绝判断与“阶梯”优化

在经典Ziggurat的“常规判断”中，我们需要计算 $f(u_1)$ 并与 $u_2 * y_i$ 比较。计算 $f(u_1)$ 需要求指数 $e^{-u_1^2/2}$。一个常见的优化是注意到，我们其实只需要比较 $u_2 * y_i$ 和 $f(u_1)$ 的大小，而不需要 $f(u_1)$ 的精确值。因此，可以比较两者的对数： $ \ln(u_2) + \ln(y_i) < -\frac{u_1^2}{2} + \text{const} $ 这样就把指数比较转化成了对数和平方比较。虽然log运算也不便宜，但有时比exp快，或者可以结合其他技巧。

更进一步的优化是使用“阶梯状”区域而非严格的矩形。例如，矩形-楔形混合法（Rectangle-Wedge Tail）或者Doornik改进算法。这些方法通过使用更贴合PDF曲线的区域形状（如梯形），来进一步降低拒绝率，虽然查表可能稍复杂，但减少了昂贵函数调用的概率。

5.2 向量化实现与SIMD

现代CPU支持单指令多数据流（SIMD）指令集（如SSE、AVX、NEON）。Ziggurat算法非常适合向量化，因为它的核心流程是对大量独立数据执行相同操作：

1. 向量化生成一批均匀随机整数（层索引）。
2. 向量化生成一批均匀随机浮点数（横坐标u1和判断值u2）。
3. 向量化执行比较操作（u1 < x[i]?）。
4. 对于未能“快速接受”的少数样本，回退到标量路径进行常规或尾部处理。

这种“向量化主路径+标量处理异常”的模式，在现代JIT编译器（如MATLAB的JIT、Julia的LLVM）或手写汇编中能带来巨大性能提升。这也是为什么像NumPy这样的库，在生成大量随机数时速度惊人的原因之一。

5.3 与现代随机数生成器的结合

Ziggurat只是一个分布变换器。它的“原料”是均匀分布的随机数。因此，其最终的速度和质量上限，很大程度上取决于底层均匀随机数生成器（URNG）。

• 高质量URNG：如Mersenne Twister (MT19937)，周期极长（$2^{19937}-1$），统计性质良好，是许多库的默认选择。但它状态空间大（2.5KB），速度不是最快。
• 高性能URNG：如PCG系列、xoshiro/xoroshiro，它们在保证良好统计性质的同时，速度更快，状态空间更小，更适用于并行和向量化。
• 硬件RNG：如Intel CPU上的RDRAND/RDSEED指令，提供基于物理熵源的随机数，适用于密码学安全场景，但速度较慢。

一个高性能的随机数库通常会分层设计：底层是一个快速的、非密码学的URNG（如PCG），中间层是Ziggurat这样的分布变换器，共同为上层应用提供高速的randn()调用。

5.4 在MATLAB环境下的实践建议

对于MATLAB用户，我的建议是：

1. 信任并优先使用内置randn：MathWorks的工程师已经为你的平台做了深度优化。在R2015a之后，MATLAB默认使用Mersenne Twister作为基础生成器，并结合了高效的变换算法（很可能是优化版的Ziggurat或类似算法）。在性能测试中，它几乎总是最优选择。
2. 向量化调用：避免在循环中调用randn(1)。总是尝试生成向量或矩阵。例如，用randn(1e6, 1)代替一百万次循环中的randn()调用，性能有数量级差异。
3. 管理随机数状态：使用rng函数来控制随机数种子，确保实验的可重复性。例如，在脚本开头执行rng(0, 'twister')。
4. 需要更高质量或不同分布时：查看Statistics and Machine Learning Toolbox中的函数，如normrnd（参数化正态分布）、mvnrnd（多元正态）等。对于极端性能需求，可以考虑用MEX文件调用C/C++实现的高性能库（如Intel MKL的VSL）。

理解Ziggurat算法，最终不是为了让你在MATLAB中自己重写一个randn，而是让你明白，在你轻松调用那些高级函数时，底层经历了怎样精妙的设计和权衡。这种理解能帮助你在遇到性能瓶颈时，知道可能的优化方向，也能让你在需要将算法移植到其他环境（如嵌入式C代码）时，做出更明智的选择。随机数生成是科学计算的基石之一，而Ziggurat无疑是这块基石上一颗高效而优雅的明珠。