Welford 算法 | 优雅地计算海量数据的均值与方差

当内存装不下数据时

在机器学习特征工程或数据分析中，我们经常遇到这样的场景：手头有成百上千个独立的特征文件（CSV、Parquet 或 Numpy 格式），总量达到了几百 GB 甚至 TB 级别。现在需要计算这些特征的全局统计量（平均值、方差、标准差）来进行归一化（Standardization）。然而，开发机内存只有 16GB。如果尝试简单的 pandas.read_csv() 或 numpy.concatenate() 把所有数据读入内存，程序会瞬间 OOM（Out of Memory）崩溃。面对据量 >> 内存的场景，我们需要一种流式（Streaming）或增量（Incremental）的计算方案——Welford 算法。

教科书公式的陷阱

在统计学教科书中，我们学到的方差计算公式通常是这样的：
σ2=∑i=1nxi2−(∑i=1nxi)2nn\sigma^2 = \frac{\sum_{i=1}^n x_i^2 - \frac{(\sum_{i=1}^n x_i)^2}{n}}{n}σ2=n∑i=1n​xi2​−n(∑i=1n​xi​)2​​

或者它的变形：
σ2=1n∑i=1nxi2−xˉ2\sigma^2 = \frac{1}{n} \sum_{i=1}^n x_i^2 - \bar{x}^2σ2=n1​i=1∑n​xi2​−xˉ2

这种公式非常适合手算，但在计算机工程实现中，它有两个致命缺陷：

1. 内存不友好：你需要维护两个巨大的和（∑x\sum x∑x 和 ∑x2\sum x^2∑x2），虽然这可以通过累加解决，但无法避免第二个问题。
2. 数值稳定性极差（Catastrophic Cancellation）：这是最严重的问题。假设你的数据数值很大（例如 10910^9109），那么 ∑x2\sum x^2∑x2 会变成一个天文数字。当公式中计算 ∑x2−nxˉ2\sum x^2 - n\bar{x}^2∑x2−nxˉ2 时，实际上是在做两个非常巨大的数字相减。在浮点数运算中，大数相减会丢失大量的有效位数，导致结果精度极低，甚至算出负的方差（这在数学上是不可能的，但在计算机里常发生）。

我们需要一种算法，它既不需要存储历史数据，又能在计算过程中保持数值较小，这就是 Welford 在线算法。

Welford 算法的核心思想

B. P. Welford 于 1962 年提出了一种迭代算法。它的核心思想是：每读入一个新的数据点，就利用旧的统计量，修正出新的统计量。
我们只需要维护三个变量：

• nnn：当前的样本计数。
• μ\muμ：当前的平均值（Mean）。
• M2M_2M2​：当前的平方差聚合值（用于计算方差）。

迭代公式每当流式数据中进来一个新的数值 xxx，更新步骤如下：

• 计数加一：n←n+1n \leftarrow n + 1n←n+1
• 更新均值：δ=x−μold,μnew←μold+δn\delta = x - \mu_{\text{old}}, \quad \mu_{\text{new}} \leftarrow \mu_{\text{old}} + \frac{\delta}{n}δ=x−μold​,μnew​←μold​+nδ​
• 更新 M2M_2M2​（最精妙的一步）：δ2=x−μnew,M2←M2+δ×δ2\delta_2 = x - \mu_{\text{new}}, \quad M_2 \leftarrow M_2 + \delta \times \delta_2δ2​=x−μnew​,M2​←M2​+δ×δ2​
• 最终的标准差计算：std=M2n−1\text{std} = \sqrt{\frac{M_2}{n-1}}std=n−1M2​​​

为什么它更优秀？

在这个公式中，我们始终在计算 x−μx - \mux−μ（数据点与均值的差值）。无论原始数据 xxx 有多大，这个差值通常都比较小。计算小数字的平方和，永远比计算大数字的平方和要精确得多。

代码实践

import numpy as np classWelfordStats:"""Welford's Online Algorithm for calculating Mean and Variance. 适合流式数据、大数据集的增量计算。 """def__init__(self): self.n =0 self.mean =0.0 self.M2 =0.0defupdate(self, x):"""处理单个数值""" self.n +=1 delta = x - self.mean self.mean += delta / self.n delta2 = x - self.mean self.M2 += delta * delta2 defupdate_batch(self, chunk):""" 处理一批数据 (Numpy Array) 在实际文件读取中，分块读取效率远高于逐行读取。 """ chunk = np.asarray(chunk)# 针对这一批数据进行循环更新for x in chunk: self.update(x)@propertydefvariance(self):"""样本方差 (Sample Variance, 分母 n-1)"""if self.n <2:return0.0return self.M2 /(self.n -1)@propertydefstd(self):"""样本标准差"""return np.sqrt(self.variance)def__str__(self):returnf"Count: {self.n}, Mean: {self.mean:.6f}, Std: {self.std:.6f}"

回到文章开头的需求，假设有一堆 CSV 文件，我们可以这样计算全局分布：

import pandas as pd import glob # 1. 初始化统计器 stats = WelfordStats()# 2. 获取文件列表 file_list = glob.glob("features/*.csv")# 3. 流式读取并计算for file_path in file_list:print(f"Processing {file_path}...")# 分块读取 (chunksize)，避免一次性读入整个文件# 假设特征列名为 'feature_value'for chunk in pd.read_csv(file_path, chunksize=10000): data_chunk = chunk['feature_value'].values # 喂给 Welford 算法 stats.update_batch(data_chunk)# 4. 输出最终结果print("-"*30)print("全局统计结果：")print(stats)

推导过程

公式 M2,n=M2,n−1+(xn−μn−1)(xn−μn)M_{2,n} = M_{2,n-1} + (x_n - \mu_{n-1})(x_n - \mu_n)M2,n​=M2,n−1​+(xn​−μn−1​)(xn​−μn​) 是成立的核心在于利用均值的递推关系进行代数变换。我们要计算的是前 nnn 个数据的离差平方和：
M2,n=∑i=1n(xi−μn)2M_{2,n} = \sum_{i=1}^n (x_i - \mu_n)^2M2,n​=i=1∑n​(xi​−μn​)2

将其拆解为前 n−1n-1n−1 项和第 nnn 项，并引入旧均值 μn−1\mu_{n-1}μn−1​ 进行配方，经过一系列消项和化简（利用 ∑(xi−μ)=0\sum(x_i - \mu)=0∑(xi​−μ)=0 的性质），最终可以证明：
∑i=1n(xi−μn)2=∑i=1n−1(xi−μn−1)2+(xn−μn−1)(xn−μn)\sum_{i=1}^n (x_i - \mu_n)^2 = \sum_{i=1}^{n-1} (x_i - \mu_{n-1})^2 + (x_n - \mu_{n-1})(x_n - \mu_n)i=1∑n​(xi​−μn​)2=i=1∑n−1​(xi​−μn−1​)2+(xn​−μn−1​)(xn​−μn​)

这正是 Welford 算法中 M2M_2M2​ 的更新逻辑：新的 M2M_2M2​ = 旧的 M2M_2M2​ + (当前值 - 旧均值) ×\times× (当前值 - 新均值)

总结

Welford 算法是处理大规模数据统计的“最佳实践”。

• 空间复杂度：O(1)O(1)O(1)。无论数据多少，只占 3 个变量的内存。
• 时间复杂度：O(N)O(N)O(N)。只需遍历一次数据。
• 数值稳定性：极高，避免了大数吃小数的问题。