蒙特卡洛方法：算不动的问题就扔骰子

前言

有些问题Deterministic的路子根本走不通，硬解是天文数字，但是拿随机数去"扔"，反而能又快又好的逼近答案。这就是蒙特卡洛方法，高性能计算领域里最气质独特的一支。

从算圆周率说起

最经典的入门例子，往单位正方形里随机撒点，数一数落在四分之一圆里的比例，命中数 $m$ 与总样本数 $N$ 之比就是面积的估计，于是 $\hat{\pi} = \dfrac{4m}{N}$：

import random

def pi_by_mc(n):
    hit = 0
    for _ in range(n):
        x = random.random()
        y = random.random()
        if x*x + y*y <= 1.0:  # 落在四分之一圆内
            hit += 1
    return 4.0 * hit / n  # 比例乘4就是pi的估计

print(pi_by_mc(1000000))  # 大约 3.141x，每次跑结果都不一样

一百万个点能算出3.14，看起来平平无奇，但这个套路厉害的地方在于：它对问题的维度完全不敏感。算一维的定积分是这个套路，算三百维空间里粒子的运动轨迹还是这个套路，代码几乎不用改。而传统的数值积分碰到高维直接指数爆炸，维度诅咒不是开玩笑的。

这方法是哪来的

蒙特卡洛这个名字来自摩纳哥那个著名赌场。上世纪四十年代曼哈顿计划期间，Los Alamos实验室要研究中子在核材料里的扩散输运问题，这玩意解析解不存在，网格法维度太高算不动。乌拉姆（Ulam）和冯诺依曼（von Neumann）这些人想出了用随机数模拟单个中子的一生：它走多远、撞不撞、往哪个方向弹，全按概率分布抽样，模拟几百万个中子的轨迹，统计出来的结果就是答案。

当时这些模拟跑在ENIAC上，是人类第一批真正的大规模科学计算。后来赌场的名字就成了方法的名字，也算赛博考古里的浪漫了。

它的命门：收敛速度

蒙特卡洛的误差 $\epsilon \propto \dfrac{1}{\sqrt{N}}$：

$$ \epsilon \sim N^{-1/2} \; \Rightarrow \; \text{样本数扩大100倍，精度才能提高10倍} $$

这就是它出名的地方，也是它被诟病的地方。想要小数点后多一位精度，代价是实打实的一百倍计算量。所以蒙特卡洛天生适合并行：每个样本之间完全独立，一亿个样本撒到一万个核上各算各的最后汇总，并行效率接近线性，这也是它算高性能计算经典算法的原因。

针对收敛慢的问题有一套方差缩减技术，最著名的是重要性采样：别傻乎乎均匀撒点，在贡献大的区域多撒点，再按概率把权重除回去，同样的样本量精度能高出好几倍。

如今它无处不在

• 金融工程的期权定价、风险价值计算
• 物理化学的分子模拟
• 机器学习里的各种采样：强化学习里agent在环境里摸爬滚打估计回报，本质上就是蒙特卡洛
• 去年AlphaGo四比一战胜李世石震动全球，它下棋搜索用的Monte Carlo Tree Search，名字里就写着Monte Carlo

某种意义上说，当今AI的半个江山，都是当年Los Alamos那批人扔下的骰子的回声。

最后

蒙特卡洛的随机数质量是生命线，拿质量差的线性同余发生器糊弄，高维问题里样本会呈现出肉眼可见的规律性条纹（相关性问题），结论直接歪掉。正经的科学计算请认准梅森旋转（Python的random默认就是它）起步，讲究点的直接上Sobol、Halton这类低差异序列，用确定性低差异序列代替随机数的玩法叫拟蒙特卡洛（QMC），同样样本量精度高出一截，配合蒙特卡洛食用风味更佳。