numpy.random.normal同种子跨系统生成值微差原因咨询

问题背景 #

因需要对接验证其他遗留程序输出，故采用旧版np.random.seed接口⁽¹⁾，在两套Linux系统上对比np.random.normal的生成结果，测试代码如下：

import numpy as np

np.random.seed(0)
x = np.random.normal(scale=1e-3, size=10**5)
np.save('test.npy', x)

将系统A生成的test.npy拷贝至系统B后对比两个数组，运行结果如下：

>>> other = np.load('test.npy')
>>> (x != other).sum(), len(x)
(29, 100000)
>>> mask = x != other
>>> np.abs(x[mask] - other[mask])
array([5.42101086e-20, 1.35525272e-20, 2.71050543e-20, 5.42101086e-20,
       1.08420217e-19, 1.08420217e-19, 2.16840434e-19, 2.16840434e-19,
       1.35525272e-20, 1.08420217e-19, 1.08420217e-19, 5.42101086e-20,
       2.71050543e-20, 1.08420217e-19, 2.16840434e-19, 5.42101086e-20,
       2.71050543e-20, 2.16840434e-19, 2.16840434e-19, 2.71050543e-20,
       2.71050543e-20, 1.08420217e-19, 1.08420217e-19, 1.08420217e-19,
       5.42101086e-20, 1.08420217e-19, 1.08420217e-19, 5.42101086e-20,
       2.71050543e-20])
>>> x[mask]
array([ 4.52489093e-04,  9.78961454e-05, -1.47113076e-04, -3.67859222e-04,
       -5.33279620e-04,  8.40794952e-04, -7.75987295e-04,  1.34205479e-03,
        6.34459482e-05,  5.07109360e-04, -7.68363366e-04,  3.33350262e-04,
       -2.19367067e-04,  6.11402140e-04, -1.30486526e-03, -4.42699624e-04,
        1.45463287e-04, -1.22491651e-03,  1.05226781e-03, -2.43032730e-04,
       -2.40551279e-04,  4.95396595e-04, -7.25454745e-04, -8.50779215e-04,
       -2.66274662e-04,  7.28854386e-04,  8.38515107e-04,  3.36152654e-04,
       -1.26550328e-04])

10万个生成元素中有29个存在量级为1e-20~1e-19的极微小数值差异。已完成以下校验排除外部因素：

• 两套系统Python版本均为3.9.4，NumPy版本测试了1.20.2与1.23.0，版本完全一致；
• 调用np.random.normal前后，通过np.random.get_state()获取的RNG状态在两套系统上完全一致；
• test.npy文件经多次拷贝、MD5校验，确认传输过程无损坏，差异来自随机数生成环节。

测试环境信息 #

• 系统A（生成test.npy的设备）：内核版本为3.10.0-1160.31.1.el7.x86_64（CentOS 7系列），同OS版本、不同CPU的设备A2生成结果与A完全一致，初步排除CPU硬件影响，怀疑差异与OS版本相关，系统信息查询命令输出如下：

$ uname -a
Linux SystemA 3.10.0-1160.31.1.el7.x86_64 #1 SMP Thu Jun 10 13:32:12 UTC 2021 x86_64 x86_64 x86_64 GNU/Linux

• 系统B（对比加载的设备）：内核版本为5.4.0-113-generic（Ubuntu 20.04系列），系统信息查询命令输出如下：

$ uname -a
Linux SystemB 5.4.0-113-generic #127-Ubuntu SMP Wed May 18 14:30:56 UTC 2022 x86_64 x86_64 x86_64 GNU/Linux

补充测试结论⁽¹⁾ #

• 使用旧版接口写法rs = RandomState(MT19937(SeedSequence(0)))时，跨系统差异仍然存在；
• 改用新版随机数接口np.random.default_rng(seed=0)（基于PCG64位生成器）时，跨系统生成结果完全一致，差异消失。

根本原因 #

这个差异的核心来源是旧版Legacy RandomState（MT19937位生成器）的正态分布采样实现，对系统数学库的浮点运算实现存在依赖，没有做跨平台的结果一致性保证：

1. 旧版np.random.normal基于Box-Muller变换实现正态分布采样，计算过程中会调用三角函数、对数、平方根等基础浮点运算，这些运算的实现来自系统自带的libm数学库。
2. 两套系统搭载的glibc版本存在代差：系统A（CentOS 7）搭载glibc 2.17，系统B（Ubuntu 20.04）搭载glibc 2.31。两个版本的libm对超越函数（sin、cos、log等）的浮点实现精度、舍入逻辑存在极微小差异，这种差异刚好落在双精度浮点数的最低有效位上。观测到的最大差值2.168e-19，恰好是双精度浮点数在1e-3数值量级下的1ULP（最小精度单位），和误差表现完全吻合。
3. RNG状态跨系统完全一致但输出不同，是因为MT19937生成的原始均匀随机位流确实完全一致，差异出在把均匀随机数转换为正态分布随机数的数学运算环节，和位生成器本身的逻辑无关。
4. 新版np.random.default_rng使用的PCG64位生成器搭配的分布采样实现，没有依赖系统libm的超越函数实现，而是内置了跨平台逐位一致的运算逻辑做分布转换，因此相同种子下跨系统生成结果完全一致，不会出现这类最低位差异。

补充说明：这种1ULP量级的浮点差异不属于计算错误，IEEE 754双精度浮点标准本身允许不同数学库实现的超越函数存在1ULP以内的输出差异，没有任何通用浮点计算环境会承诺不同系统libm的计算结果逐位相等。如果需要跨环境逐位复现随机数结果，必须使用新版numpy.random接口。

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