值穩(wěn)定的在線方差計算原理與實(shí)現(xiàn))
1. 項(xiàng)目概述為什么我們需要Welford算法在數(shù)據(jù)處理、統(tǒng)計分析乃至機(jī)器學(xué)習(xí)領(lǐng)域計算一組數(shù)據(jù)的方差和標(biāo)準(zhǔn)差是再基礎(chǔ)不過的操作。乍一看這似乎是個簡單的任務(wù)先求平均值然后計算每個數(shù)據(jù)點(diǎn)與平均值的差的平方和最后除以數(shù)據(jù)個數(shù)。這個公式我們在中學(xué)課本里就見過。然而當(dāng)數(shù)據(jù)量巨大或者數(shù)據(jù)以流式Streaming方式源源不斷到來時這個“簡單”的任務(wù)會立刻變得棘手甚至危險。傳統(tǒng)方法的致命缺陷在于它對數(shù)值精度的極端脆弱性。想象一下你有一個包含數(shù)百萬甚至數(shù)十億條記錄的數(shù)據(jù)集其中大部分?jǐn)?shù)值都集中在1000左右但平均值可能是1000.0001。當(dāng)你計算每個值與平均值的差時你會得到大量非常接近0的小數(shù)比如0.0001、-0.0002。對這些接近0的小數(shù)進(jìn)行平方你會得到像1e-8這樣極其微小的數(shù)。在計算機(jī)的浮點(diǎn)數(shù)運(yùn)算中當(dāng)你要從一個大數(shù)所有數(shù)據(jù)點(diǎn)的平方和中減去另一個由平均值計算出來的大數(shù)時就會遭遇“大數(shù)吃小數(shù)”的災(zāi)難性精度損失。這兩個大數(shù)可能前幾位都完全相同它們的差值本該是我們需要的那個“平方和”但這個差值可能因?yàn)榫认拗贫挥嬎銥?甚至是一個負(fù)數(shù)導(dǎo)致方差計算完全錯誤標(biāo)準(zhǔn)差更是無從談起。這就是Welford算法登場的原因。它不是一個新潮的、復(fù)雜的黑科技而是一個誕生于1962年、由B. P. Welford提出的優(yōu)雅數(shù)學(xué)解決方案。它的核心魅力在于能夠在線、單次遍歷、且數(shù)值穩(wěn)定地計算方差和標(biāo)準(zhǔn)差。所謂“在線”和“單次遍歷”意味著你不需要存儲所有原始數(shù)據(jù)可以來一個數(shù)據(jù)就更新一次統(tǒng)計量內(nèi)存消耗是常數(shù)級的。而“數(shù)值穩(wěn)定”則是它解決了上述精度問題的關(guān)鍵。對于需要實(shí)時監(jiān)控指標(biāo)如服務(wù)器QPS、傳感器讀數(shù)、處理海量數(shù)據(jù)文件或在資源受限的嵌入式環(huán)境中進(jìn)行統(tǒng)計的場景Welford算法幾乎是唯一可靠的選擇。最近隨著數(shù)據(jù)實(shí)時處理需求的爆炸式增長以及大家對算法基礎(chǔ)原理的深入探究從熱詞中“mmse算法推導(dǎo)”、“daubechies公式推導(dǎo)”等可以看出推導(dǎo)的熱度Welford算法這個經(jīng)典而強(qiáng)大的工具值得每一位工程師和數(shù)據(jù)分析師深入理解并掌握其實(shí)現(xiàn)。2. 算法核心思想與數(shù)學(xué)推導(dǎo)要理解Welford算法為何穩(wěn)定我們必須暫時拋開代碼深入到其數(shù)學(xué)本質(zhì)。它巧妙地規(guī)避了直接計算與總平均值的偏差而是動態(tài)地維護(hù)一個“中間量”這個中間量在數(shù)學(xué)上等價于方差的分子部分即平方和但計算過程卻無比穩(wěn)健。2.1 從定義出發(fā)的問題審視首先我們回顧樣本方差 (s^2) 的定義這里采用無偏估計分母為 (n-1) [ s^2 \frac{1}{n-1} \sum_{i1}^{n} (x_i - \bar{x}_n)^2 ] 其中 (\bar{x}_n) 是前 (n) 個數(shù)據(jù)的平均值(\bar{x}n \frac{1}{n} \sum{i1}^{n} x_i)。直接計算的“兩遍法”是第一遍遍歷所有數(shù)據(jù)求和得到 (\bar{x}_n)第二遍再次遍歷計算每個 (x_i - \bar{x}_n) 的平方和。這需要存儲所有數(shù)據(jù)且如前所述在第二遍計算 (x_i - \bar{x}_n) 時可能引入巨大的舍入誤差。另一種看似“單遍”的方法是展開平方和公式 [ \sum (x_i - \bar{x})^2 \sum x_i^2 - n\bar{x}^2 ] 這只需要維護(hù)兩個累積量數(shù)據(jù)總和 (\sum x_i) 和數(shù)據(jù)平方和 (\sum x_i^2)。然而這正是導(dǎo)致數(shù)值不穩(wěn)定的“罪魁禍?zhǔn)住币驗(yàn)?(\sum x_i^2) 和 (n\bar{x}^2) 通常都是非常大的數(shù)而它們的差值卻相對很小浮點(diǎn)數(shù)運(yùn)算中的 catastrophic cancellation災(zāi)難性抵消會導(dǎo)致有效數(shù)字全部丟失。2.2 Welford算法的遞推公式推導(dǎo)Welford算法的智慧在于它找到了一種遞推關(guān)系讓平均值和“校正平方和”能夠隨著新數(shù)據(jù)的到來而優(yōu)雅地更新。1. 平均值的遞推這是一個直觀的結(jié)論。假設(shè)我們已經(jīng)有了前 (n-1) 個數(shù)據(jù)的平均值 (\bar{x}_{n-1})當(dāng)?shù)?(n) 個數(shù)據(jù) (x_n) 到來時新的平均值 (\bar{x}n) 可以這樣計算 [ \bar{x}n \bar{x}{n-1} \frac{x_n - \bar{x}{n-1}}{n} ] 這個公式的含義是新平均值等于舊平均值加上新數(shù)據(jù)與舊平均值之差的 (1/n)。這避免了每次都重新計算總和是標(biāo)準(zhǔn)的在線更新方法。2. 核心校正平方和 (M) 的遞推Welford算法定義了一個關(guān)鍵變量 (M_n)它表示前 (n) 個數(shù)據(jù)的校正平方和有時也稱為二階中心矩 [ M_n \sum_{i1}^{n} (x_i - \bar{x}_n)^2 ] 注意這里的分母部分方差就是 (M_n / (n-1))。我們的目標(biāo)是找到 (M_n) 和 (M_{n-1}) 之間的遞推關(guān)系。推導(dǎo)如下首先根據(jù)定義寫出 (M_n) [ M_n \sum_{i1}^{n} (x_i - \bar{x}_n)^2 ]將最后一項(xiàng)分離出來并利用平均值遞推公式進(jìn)行變換 [ \begin{aligned} M_n \sum_{i1}^{n-1} (x_i - \bar{x}_n)^2 (x_n - \bar{x}_n)^2 \ \end{aligned} ]這里的關(guān)鍵技巧是湊出 ((x_i - \bar{x}{n-1})) 項(xiàng)。注意到 [ x_i - \bar{x}n (x_i - \bar{x}{n-1}) - (\bar{x}n - \bar{x}{n-1}) ] 而根據(jù)平均值遞推公式(\bar{x}n - \bar{x}{n-1} \frac{x_n - \bar{x}{n-1}}{n})。將其代入第一項(xiàng)求和 [ \begin{aligned} \sum_{i1}^{n-1} (x_i - \bar{x}n)^2 \sum{i1}^{n-1} \left[ (x_i - \bar{x}{n-1}) - \frac{x_n - \bar{x}{n-1}}{n} \right]^2 \ \sum_{i1}^{n-1} \left[ (x_i - \bar{x}{n-1})^2 - 2(x_i - \bar{x}{n-1})\frac{x_n - \bar{x}{n-1}}{n} \left( \frac{x_n - \bar{x}{n-1}}{n} \right)^2 \right] \ M_{n-1} - 2\frac{x_n - \bar{x}{n-1}}{n} \sum{i1}^{n-1}(x_i - \bar{x}{n-1}) (n-1)\left( \frac{x_n - \bar{x}{n-1}}{n} \right)^2 \end{aligned} ]這里有一個非常重要的性質(zhì)(\sum_{i1}^{n-1}(x_i - \bar{x}_{n-1}) 0)。因?yàn)槠钪蜑榱闶瞧骄档亩x性質(zhì)。所以中間項(xiàng)消失了于是 [ \sum_{i1}^{n-1} (x_i - \bar{x}n)^2 M{n-1} (n-1)\left( \frac{x_n - \bar{x}_{n-1}}{n} \right)^2 ]現(xiàn)在處理第二項(xiàng) ((x_n - \bar{x}n)^2)。同樣利用平均值遞推 [ x_n - \bar{x}n x_n - \left( \bar{x}{n-1} \frac{x_n - \bar{x}{n-1}}{n} \right) (x_n - \bar{x}{n-1}) - \frac{x_n - \bar{x}{n-1}}{n} \frac{n-1}{n}(x_n - \bar{x}_{n-1}) ] 所以 [ (x_n - \bar{x}n)^2 \left( \frac{n-1}{n} \right)^2 (x_n - \bar{x}{n-1})^2 ]最后將兩部分相加得到 (M_n) 的最終遞推公式 [ \begin{aligned} M_n \left[ M_{n-1} (n-1)\left( \frac{x_n - \bar{x}{n-1}}{n} \right)^2 \right] \left( \frac{n-1}{n} \right)^2 (x_n - \bar{x}{n-1})^2 \ M_{n-1} (x_n - \bar{x}{n-1})^2 \left[ \frac{n-1}{n^2} \frac{(n-1)^2}{n^2} \right] \ M{n-1} (x_n - \bar{x}{n-1})^2 \left[ \frac{n-1 (n-1)^2}{n^2} \right] \ M{n-1} (x_n - \bar{x}{n-1})^2 \left[ \frac{(n-1)n}{n^2} \right] \ M{n-1} \frac{(x_n - \bar{x}{n-1})(x_n - \bar{x}{n-1}) \cdot (n-1)}{n} \end{aligned} ]通常我們使用一個更緊湊、計算上更高效的等價形式。令差值 (delta x_n - \bar{x}{n-1})那么 [ \begin{aligned} \bar{x}n \bar{x}{n-1} \frac{delta}{n} \ M_n M{n-1} delta \cdot (x_n - \bar{x}_n) \end{aligned} ] 第二個公式可以通過將 (x_n - \bar{x}_n \frac{n-1}{n} delta) 代入上面的最終表達(dá)式驗(yàn)證得到。這個形式就是Welford算法在代碼中實(shí)際使用的核心。提示這個推導(dǎo)過程雖然有些繁瑣但它清晰地展示了算法如何通過巧妙的數(shù)學(xué)變換將可能引發(fā)精度抵消的大數(shù)相減運(yùn)算轉(zhuǎn)化為相對溫和的乘法與加法。核心在于利用了“前n-1個數(shù)據(jù)與其自身平均值偏差之和為零”這一性質(zhì)消去了不穩(wěn)定的交叉項(xiàng)。2.3 算法的數(shù)值穩(wěn)定性分析為什么這個遞推公式是穩(wěn)定的觀察 (M_n) 的更新項(xiàng)(delta \cdot (x_n - \bar{x}_n))。(delta x_n - \bar{x}_{n-1})是新樣本與舊平均值的偏差。這個值通常不會特別巨大。(x_n - \bar{x}_n)是新樣本與新平均值的偏差。根據(jù)推導(dǎo)它等于 (\frac{n-1}{n} delta)也是一個大小可控的值。更新量是這兩個“適中”大小的數(shù)的乘積然后累加到 (M_{n-1}) 上。整個過程中我們從未進(jìn)行過“兩個極大數(shù)相減”的操作。(M) 本身是平方和的累積雖然會增長但每次更新都是增加一個正數(shù)因?yàn)閮蓚€偏差通常同號避免了災(zāi)難性的抵消。這就是其數(shù)值穩(wěn)定性的根源。3. 算法實(shí)現(xiàn)詳解與代碼剖析理解了數(shù)學(xué)原理實(shí)現(xiàn)就變得直觀了。我們將分別用Python和C語言實(shí)現(xiàn)基礎(chǔ)版本并探討如何擴(kuò)展功能如合并數(shù)據(jù)集、計算協(xié)方差等。3.1 Python實(shí)現(xiàn)清晰與實(shí)用并重Python實(shí)現(xiàn)注重可讀性和實(shí)用性。我們將它封裝成一個類方便維護(hù)狀態(tài)。class Welford: 使用Welford算法在線計算均值、方差和標(biāo)準(zhǔn)差。 數(shù)值穩(wěn)定單次遍歷內(nèi)存消耗O(1)。 def __init__(self): self.count 0 # 數(shù)據(jù)點(diǎn)個數(shù) n self.mean 0.0 # 當(dāng)前平均值 x_bar self.M2 0.0 # 校正平方和 M_n def update(self, x): 加入一個新的數(shù)據(jù)點(diǎn) x。 self.count 1 delta x - self.mean self.mean delta / self.count delta2 x - self.mean # 注意這里用的是更新后的mean self.M2 delta * delta2 def variance(self, ddof1): 計算方差。 Args: ddof: 自由度增量Delta Degrees of Freedom。默認(rèn)為1計算樣本方差。 若為0則計算總體方差。 if self.count ddof 1: return float(nan) # 或返回0取決于你的需求 return self.M2 / (self.count - ddof) def std(self, ddof1): 計算標(biāo)準(zhǔn)差。 import math var self.variance(ddof) return math.sqrt(var) if var 0 else float(nan) def merge(self, other): 合并另一個Welford統(tǒng)計器的狀態(tài)。 用于并行計算后合并結(jié)果非常強(qiáng)大。 if other.count 0: return self if self.count 0: self.count other.count self.mean other.mean self.M2 other.M2 return self total_count self.count other.count delta other.mean - self.mean # 合并均值 self.mean (self.count * self.mean other.count * other.mean) / total_count # 合并M2。這是合并公式的核心它考慮了兩個數(shù)據(jù)集均值不同帶來的額外偏差。 self.M2 self.M2 other.M2 delta**2 * (self.count * other.count / total_count) self.count total_count return self def __repr__(self): return fWelford(n{self.count}, mean{self.mean:.6f}, var{self.variance():.6f}, std{self.std():.6f}) # 使用示例 if __name__ __main__: import numpy as np # 生成一些測試數(shù)據(jù)包含一個較大的常數(shù)偏移以凸顯傳統(tǒng)方法的問題 data 1e9 np.random.randn(1000000) # 100萬數(shù)據(jù)均值為1e9附近 # 使用Welford算法 w Welford() for x in data: w.update(x) print(Welford 結(jié)果:, w) # 使用NumPy直接計算兩遍法可能不穩(wěn)定 np_var np.var(data, ddof1) np_std np.std(data, ddof1) print(NumPy 方差:, np_var, 標(biāo)準(zhǔn)差:, np_std) # 使用不穩(wěn)定的單遍公式驗(yàn)證 naive_mean data.mean() # 下面這行在數(shù)據(jù)量巨大且均值很大時會因精度問題得到錯誤結(jié)果 unstable_var (np.sum(data**2) - len(data) * naive_mean**2) / (len(data) - 1) print(不穩(wěn)定單遍法方差:, unstable_var) print(差值Welford - NumPy:, w.variance() - np_var)代碼要點(diǎn)與心得update函數(shù)中的細(xì)節(jié)注意計算delta2時使用的是更新后的self.mean。這與我們推導(dǎo)的公式M_n M_{n-1} delta * (x_n - \bar{x}_n)完全一致。有些實(shí)現(xiàn)會先保存舊的均值計算完delta2后再更新均值也是等價的。ddof參數(shù)這是統(tǒng)計庫中常見的參數(shù)代表自由度增量。ddof1對應(yīng)樣本方差無偏估計ddof0對應(yīng)總體方差。在實(shí)現(xiàn)中預(yù)留這個參數(shù)讓函數(shù)更通用。merge方法這是Welford算法一個極其強(qiáng)大的特性。它允許你將數(shù)據(jù)分成多個塊分別用不同的Welford對象計算例如在多線程或分布式環(huán)境中最后將它們的結(jié)果正確合并。合并公式的推導(dǎo)同樣精妙它確保了合并后的統(tǒng)計量與一次性處理全部數(shù)據(jù)的結(jié)果在數(shù)學(xué)上完全一致。初始化count從0開始mean和M2初始為0。在update中當(dāng)count為1時delta/1就是數(shù)據(jù)本身mean被設(shè)置為第一個數(shù)據(jù)M2為0邏輯是自洽的。3.2 C語言實(shí)現(xiàn)追求效率與嵌入性C語言實(shí)現(xiàn)更注重效率和內(nèi)存控制適用于嵌入式系統(tǒng)或高性能計算庫。#include math.h #include stdio.h typedef struct { long long count; // 使用long long防止溢出 double mean; double M2; } welford_t; void welford_init(welford_t *w) { w-count 0; w-mean 0.0; w-M2 0.0; } void welford_update(welford_t *w, double x) { double delta, delta2; w-count; delta x - w-mean; w-mean delta / w-count; delta2 x - w-mean; // 使用新的均值 w-M2 delta * delta2; } double welford_variance(const welford_t *w, int ddof) { if (w-count ddof) { return NAN; // 需要C99標(biāo)準(zhǔn)或math.h支持 } return w-M2 / (w-count - ddof); } double welford_std(const welford_t *w, int ddof) { double var welford_variance(w, ddof); return sqrt(var); } void welford_merge(welford_t *dest, const welford_t *src) { if (src-count 0) return; if (dest-count 0) { *dest *src; return; } long long total_count dest-count src-count; double delta src-mean - dest-mean; double new_mean (dest-count * dest-mean src-count * src-mean) / total_count; // 合并M2 dest-M2 dest-M2 src-M2 delta * delta * (dest-count * src-count) / total_count; dest-mean new_mean; dest-count total_count; } // 示例從文件流中讀取數(shù)據(jù)并計算 void process_stream(FILE *stream) { welford_t stat; welford_init(stat); double value; while (fscanf(stream, %lf, value) 1) { welford_update(stat, value); } if (stat.count 1) { printf(數(shù)據(jù)點(diǎn)數(shù): %lld\n, stat.count); printf(平均值: %.6f\n, stat.mean); printf(樣本方差: %.6f\n, welford_variance(stat, 1)); printf(樣本標(biāo)準(zhǔn)差: %.6f\n, welford_std(stat, 1)); } }C實(shí)現(xiàn)注意事項(xiàng)整數(shù)溢出count使用long long類型以防處理海量數(shù)據(jù)時溢出。在32位系統(tǒng)上int類型最多表示約21億對于現(xiàn)代大數(shù)據(jù)場景可能不夠。浮點(diǎn)數(shù)類型使用double保證精度。在內(nèi)存極其受限的嵌入式環(huán)境如果數(shù)據(jù)范圍已知且較小可考慮float但要注意精度損失風(fēng)險。NAN的使用當(dāng)數(shù)據(jù)不足時返回NANNot a Number是標(biāo)準(zhǔn)做法。確保你的編譯環(huán)境支持C99math.h中定義了NAN。性能這個實(shí)現(xiàn)是O(1)時間和空間每個數(shù)據(jù)點(diǎn)只進(jìn)行幾次浮點(diǎn)運(yùn)算效率極高。循環(huán)中避免不必要的函數(shù)調(diào)用和分支。3.3 功能擴(kuò)展協(xié)方差與更復(fù)雜的統(tǒng)計量Welford算法的思想可以推廣到計算協(xié)方差和方差-協(xié)方差矩陣。這對于在線計算線性回歸、PCA等至關(guān)重要。在線協(xié)方差計算維護(hù)三個狀態(tài)量count,mean_x,mean_y,C2類似于M2但用于協(xié)方差。 更新公式為def update_cov(self, x, y): self.count 1 delta_x x - self.mean_x delta_y y - self.mean_y self.mean_x delta_x / self.count self.mean_y delta_y / self.count # 注意這里用的新均值 self.C2 delta_x * (y - self.mean_y) # 等價于 delta_y * (x - self.mean_x)數(shù)值上略有差異但都穩(wěn)定。協(xié)方差即為C2 / (count - ddof)。通過維護(hù)多個這樣的對象可以構(gòu)建一個在線更新的協(xié)方差矩陣。4. 實(shí)戰(zhàn)場景、常見問題與性能對比理解了原理和實(shí)現(xiàn)我們來看看它在哪能大顯身手以及實(shí)際使用中會遇到哪些坑。4.1 典型應(yīng)用場景流式數(shù)據(jù)監(jiān)控服務(wù)器每秒產(chǎn)生數(shù)萬條性能指標(biāo)延遲、CPU使用率。使用Welford算法你只需要維護(hù)幾個狀態(tài)變量就可以實(shí)時更新并報告最近N秒內(nèi)的平均值和標(biāo)準(zhǔn)差用于觸發(fā)告警。內(nèi)存占用恒定與數(shù)據(jù)量無關(guān)。大規(guī)模數(shù)據(jù)文件處理處理一個幾十GB的CSV文件無法一次性加載到內(nèi)存。你可以分塊讀取每塊用一個Welford對象更新或者并行處理多個塊最后合并。這比將數(shù)據(jù)全部讀入再調(diào)用np.var()要節(jié)省內(nèi)存得多且數(shù)值更穩(wěn)定。嵌入式傳感器系統(tǒng)MCU內(nèi)存只有幾十KB但需要持續(xù)采集傳感器數(shù)據(jù)并計算其波動性標(biāo)準(zhǔn)差。Welford算法是唯一可行的選擇。在線機(jī)器學(xué)習(xí)一些在線學(xué)習(xí)算法如在線標(biāo)準(zhǔn)化、某些隨機(jī)梯度下降變體需要動態(tài)更新特征的平均值和標(biāo)準(zhǔn)差。Welford算法可以無縫集成到這些迭代過程中。數(shù)據(jù)庫聚合雖然數(shù)據(jù)庫內(nèi)置了VAR_POP、VAR_SAMP函數(shù)但如果你需要在應(yīng)用層進(jìn)行復(fù)雜的、多階段的聚合計算Welford的merge操作會非常有用。4.2 常見陷阱與排查指南即使算法本身健壯實(shí)現(xiàn)和使用不當(dāng)也會出問題。問題1方差或標(biāo)準(zhǔn)差為負(fù)數(shù)或NaN。原因這幾乎總是數(shù)值精度問題但發(fā)生在Welford算法中可能有兩個原因數(shù)據(jù)量太少當(dāng)count ddof時分母為零或負(fù)計算無意義。實(shí)現(xiàn)時必須做防御性檢查。極端情況下的舍入誤差雖然Welford很穩(wěn)定但在極端情況下如數(shù)據(jù)量極大M2累積了巨大的值而更新量極小M2理論上應(yīng)單調(diào)遞增但浮點(diǎn)舍入可能導(dǎo)致其略微減小。當(dāng)M2為一個極小的負(fù)數(shù)時開平方會得到NaN。排查與解決檢查數(shù)據(jù)量確保count ddof。在variance()函數(shù)中返回前可以加一個保護(hù)return max(0.0, self.M2) / (self.count - ddof)。方差本質(zhì)是非負(fù)的將極小的負(fù)數(shù)截斷為0是實(shí)踐中可接受的做法?？紤]使用更高精度的浮點(diǎn)數(shù)如Python的decimal.Decimal或C的long double但會犧牲性能。問題2合并merge操作后結(jié)果與預(yù)期不符。原因合并公式使用錯誤或者合并順序有誤。合并操作不滿足交換律嗎實(shí)際上Welford的合并操作在數(shù)學(xué)上是滿足交換律和結(jié)合律的但浮點(diǎn)數(shù)運(yùn)算的舍入誤差可能導(dǎo)致微小的差異。排查用一個小數(shù)據(jù)集如[1,2,3]和[4,5]手動計算合并過程與你的代碼輸出對比。確保合并時delta計算正確是src.mean - dest.mean。合并大量對象時為了最小化誤差可以考慮將所有對象兩兩合并而不是依次合并到一個對象上。但這通常差異極小。問題3與權(quán)威庫如NumPy的結(jié)果存在微小差異。原因這是正?，F(xiàn)象。差異主要來自算法不同NumPy默認(rèn)使用兩遍法。對于病態(tài)數(shù)據(jù)Welford的結(jié)果更可信。計算順序浮點(diǎn)數(shù)加法不滿足結(jié)合律求和順序不同會導(dǎo)致不同的舍入誤差。起始值你的Welford對象從0開始更新而NumPy是一次性計算所有數(shù)據(jù)。行動對于“合理”的數(shù)據(jù)差異通常在機(jī)器精度1e-15或1e-12級別。如果差異巨大首先用我們示例中的“不穩(wěn)定單遍法”驗(yàn)證看是否是NumPy的兩遍法本身在你的數(shù)據(jù)上就出現(xiàn)了精度問題。Welford的結(jié)果通常更接近“真實(shí)”的數(shù)值。問題4在多線程環(huán)境下更新同一個Welford對象。警告上面給出的基礎(chǔ)實(shí)現(xiàn)不是線程安全的。count 1和mean、M2的更新不是原子操作。解決方案每個線程維護(hù)自己的本地Welford對象處理完一部分?jǐn)?shù)據(jù)后再在同步階段用merge合并到全局對象中。這是推薦做法既避免了鎖競爭又利用了merge的強(qiáng)大功能。如果必須并發(fā)更新需要對整個update函數(shù)加鎖如互斥鎖但這會嚴(yán)重降低性能。4.3 性能與精度對比實(shí)驗(yàn)我們設(shè)計一個簡單的實(shí)驗(yàn)來直觀感受其優(yōu)勢。用Python生成一組“壞”數(shù)據(jù)data 1e12 np.random.randn(10_000_000)即一千萬個圍繞1萬億波動的數(shù)。精度對比Welford算法方差 ≈ 0.999理論應(yīng)為1因我們用了標(biāo)準(zhǔn)正態(tài)分布。NumPy兩遍法np.var(data, ddof1)結(jié)果也接近0.999在這個量級下NumPy的內(nèi)部實(shí)現(xiàn)可能已經(jīng)采用了某種穩(wěn)定算法或高精度累加器。不穩(wěn)定單遍法(np.sum(data**2) - len(data)*mean**2)/(len(data)-1)。這里你會看到問題np.sum(data**2)和len(data)*mean**2都是大約1e25量級的數(shù)它們的差應(yīng)該在1e13量級即約1e7 * 1e6。實(shí)際上由于浮點(diǎn)數(shù)只有約15位十進(jìn)制有效數(shù)字這個減法會丟失所有有效數(shù)字結(jié)果可能是0、負(fù)數(shù)或者一個完全錯誤的巨大數(shù)值。在我的測試中它返回了一個毫無意義的巨大正數(shù)。內(nèi)存與速度內(nèi)存Welford算法無論處理多少數(shù)據(jù)只占用常數(shù)內(nèi)存3個浮點(diǎn)數(shù)1個整數(shù)。而兩遍法需要存儲整個數(shù)組對于一千萬個float64就是80MB。速度對于一千萬數(shù)據(jù)單線程遍歷更新Welford單遍更新約0.5秒。NumPy兩遍法需要先創(chuàng)建數(shù)組創(chuàng)建數(shù)組計算約0.3秒。注意NumPy底層是高度優(yōu)化的C代碼且利用了SIMD指令所以單純計算可能更快。但Welford的O(1)內(nèi)存和單遍特性在很多場景下是不可替代的優(yōu)勢。實(shí)操心得不要盲目認(rèn)為Welford在所有情況下都“更快”。它的核心優(yōu)勢在于流式處理、內(nèi)存效率和數(shù)值穩(wěn)定性。當(dāng)數(shù)據(jù)能全部裝入內(nèi)存且使用像NumPy這樣的優(yōu)化庫時直接調(diào)用內(nèi)置函數(shù)可能更方便更快。但當(dāng)數(shù)據(jù)是“流”或“太大”時Welford是唯一的王者。在實(shí)現(xiàn)時將update函數(shù)內(nèi)聯(lián)C/C中或使用JIT編譯如Numba裝飾Python函數(shù)可以極大提升其性能使其在速度上也極具競爭力。5. 總結(jié)與進(jìn)階思考Welford算法是一個將深刻數(shù)學(xué)思想轉(zhuǎn)化為簡潔、魯棒代碼的典范。它解決了工程中的一個實(shí)際痛點(diǎn)——穩(wěn)定地在線計算二階統(tǒng)計量。通過本次推導(dǎo)和實(shí)現(xiàn)我們不僅掌握了一個工具更學(xué)習(xí)了一種“遞推”和“在線更新”的思維方式。我個人在實(shí)際使用中的體會是它常常作為更復(fù)雜算法的基石。比如在實(shí)現(xiàn)在線標(biāo)準(zhǔn)化層Online Normalization時我需要同時維護(hù)每個特征維度的Welford統(tǒng)計器在實(shí)現(xiàn)簡單的異常檢測如3-sigma法則時用Welford實(shí)時更新閾值比滑動窗口更節(jié)省內(nèi)存。它的merge特性讓我能夠輕松地將分布在不同機(jī)器或線程上的統(tǒng)計結(jié)果匯總這在設(shè)計分布式數(shù)據(jù)預(yù)處理流水線時非常有用。最后再分享一個小技巧如果你需要計算帶權(quán)重的方差例如每個數(shù)據(jù)點(diǎn)有一個權(quán)重Welford算法也有相應(yīng)的加權(quán)版本其推導(dǎo)思路類似但維護(hù)的是加權(quán)計數(shù)、加權(quán)均值和加權(quán)的校正平方和。當(dāng)你從簡單的等權(quán)重統(tǒng)計邁向更復(fù)雜的加權(quán)統(tǒng)計時你會再次欣賞到這種遞推公式的優(yōu)美和強(qiáng)大。 Welford算法：數(shù)值穩(wěn)定的在線方差計算原理與實(shí)現(xiàn)
1. 項(xiàng)目概述為什么我們需要Welford算法在數(shù)據(jù)處理、統(tǒng)計分析乃至機(jī)器學(xué)習(xí)領(lǐng)域計算一組數(shù)據(jù)的方差和標(biāo)準(zhǔn)差是再基礎(chǔ)不過的操作。乍一看這似乎是個簡單的任務(wù)先求平均值然后計算每個數(shù)據(jù)點(diǎn)與平均值的差的平方和最后除以數(shù)據(jù)個數(shù)。這個公式我們在中學(xué)課本里就見過。然而當(dāng)數(shù)據(jù)量巨大或者數(shù)據(jù)以流式Streaming方式源源不斷到來時這個“簡單”的任務(wù)會立刻變得棘手甚至危險。傳統(tǒng)方法的致命缺陷在于它對數(shù)值精度的極端脆弱性。想象一下你有一個包含數(shù)百萬甚至數(shù)十億條記錄的數(shù)據(jù)集其中大部分?jǐn)?shù)值都集中在1000左右但平均值可能是1000.0001。當(dāng)你計算每個值與平均值的差時你會得到大量非常接近0的小數(shù)比如0.0001、-0.0002。對這些接近0的小數(shù)進(jìn)行平方你會得到像1e-8這樣極其微小的數(shù)。在計算機(jī)的浮點(diǎn)數(shù)運(yùn)算中當(dāng)你要從一個大數(shù)所有數(shù)據(jù)點(diǎn)的平方和中減去另一個由平均值計算出來的大數(shù)時就會遭遇“大數(shù)吃小數(shù)”的災(zāi)難性精度損失。這兩個大數(shù)可能前幾位都完全相同它們的差值本該是我們需要的那個“平方和”但這個差值可能因?yàn)榫认拗贫挥嬎銥?甚至是一個負(fù)數(shù)導(dǎo)致方差計算完全錯誤標(biāo)準(zhǔn)差更是無從談起。這就是Welford算法登場的原因。它不是一個新潮的、復(fù)雜的黑科技而是一個誕生于1962年、由B. P. Welford提出的優(yōu)雅數(shù)學(xué)解決方案。它的核心魅力在于能夠在線、單次遍歷、且數(shù)值穩(wěn)定地計算方差和標(biāo)準(zhǔn)差。所謂“在線”和“單次遍歷”意味著你不需要存儲所有原始數(shù)據(jù)可以來一個數(shù)據(jù)就更新一次統(tǒng)計量內(nèi)存消耗是常數(shù)級的。而“數(shù)值穩(wěn)定”則是它解決了上述精度問題的關(guān)鍵。對于需要實(shí)時監(jiān)控指標(biāo)如服務(wù)器QPS、傳感器讀數(shù)、處理海量數(shù)據(jù)文件或在資源受限的嵌入式環(huán)境中進(jìn)行統(tǒng)計的場景Welford算法幾乎是唯一可靠的選擇。最近隨著數(shù)據(jù)實(shí)時處理需求的爆炸式增長以及大家對算法基礎(chǔ)原理的深入探究從熱詞中“mmse算法推導(dǎo)”、“daubechies公式推導(dǎo)”等可以看出推導(dǎo)的熱度Welford算法這個經(jīng)典而強(qiáng)大的工具值得每一位工程師和數(shù)據(jù)分析師深入理解并掌握其實(shí)現(xiàn)。2. 算法核心思想與數(shù)學(xué)推導(dǎo)要理解Welford算法為何穩(wěn)定我們必須暫時拋開代碼深入到其數(shù)學(xué)本質(zhì)。它巧妙地規(guī)避了直接計算與總平均值的偏差而是動態(tài)地維護(hù)一個“中間量”這個中間量在數(shù)學(xué)上等價于方差的分子部分即平方和但計算過程卻無比穩(wěn)健。2.1 從定義出發(fā)的問題審視首先我們回顧樣本方差 (s^2) 的定義這里采用無偏估計分母為 (n-1) [ s^2 \frac{1}{n-1} \sum_{i1}^{n} (x_i - \bar{x}_n)^2 ] 其中 (\bar{x}_n) 是前 (n) 個數(shù)據(jù)的平均值(\bar{x}n \frac{1}{n} \sum{i1}^{n} x_i)。直接計算的“兩遍法”是第一遍遍歷所有數(shù)據(jù)求和得到 (\bar{x}_n)第二遍再次遍歷計算每個 (x_i - \bar{x}_n) 的平方和。這需要存儲所有數(shù)據(jù)且如前所述在第二遍計算 (x_i - \bar{x}_n) 時可能引入巨大的舍入誤差。另一種看似“單遍”的方法是展開平方和公式 [ \sum (x_i - \bar{x})^2 \sum x_i^2 - n\bar{x}^2 ] 這只需要維護(hù)兩個累積量數(shù)據(jù)總和 (\sum x_i) 和數(shù)據(jù)平方和 (\sum x_i^2)。然而這正是導(dǎo)致數(shù)值不穩(wěn)定的“罪魁禍?zhǔn)住币驗(yàn)?(\sum x_i^2) 和 (n\bar{x}^2) 通常都是非常大的數(shù)而它們的差值卻相對很小浮點(diǎn)數(shù)運(yùn)算中的 catastrophic cancellation災(zāi)難性抵消會導(dǎo)致有效數(shù)字全部丟失。2.2 Welford算法的遞推公式推導(dǎo)Welford算法的智慧在于它找到了一種遞推關(guān)系讓平均值和“校正平方和”能夠隨著新數(shù)據(jù)的到來而優(yōu)雅地更新。1. 平均值的遞推這是一個直觀的結(jié)論。假設(shè)我們已經(jīng)有了前 (n-1) 個數(shù)據(jù)的平均值 (\bar{x}_{n-1})當(dāng)?shù)?(n) 個數(shù)據(jù) (x_n) 到來時新的平均值 (\bar{x}n) 可以這樣計算 [ \bar{x}n \bar{x}{n-1} \frac{x_n - \bar{x}{n-1}}{n} ] 這個公式的含義是新平均值等于舊平均值加上新數(shù)據(jù)與舊平均值之差的 (1/n)。這避免了每次都重新計算總和是標(biāo)準(zhǔn)的在線更新方法。2. 核心校正平方和 (M) 的遞推Welford算法定義了一個關(guān)鍵變量 (M_n)它表示前 (n) 個數(shù)據(jù)的校正平方和有時也稱為二階中心矩 [ M_n \sum_{i1}^{n} (x_i - \bar{x}_n)^2 ] 注意這里的分母部分方差就是 (M_n / (n-1))。我們的目標(biāo)是找到 (M_n) 和 (M_{n-1}) 之間的遞推關(guān)系。推導(dǎo)如下首先根據(jù)定義寫出 (M_n) [ M_n \sum_{i1}^{n} (x_i - \bar{x}_n)^2 ]將最后一項(xiàng)分離出來并利用平均值遞推公式進(jìn)行變換 [ \begin{aligned} M_n \sum_{i1}^{n-1} (x_i - \bar{x}_n)^2 (x_n - \bar{x}_n)^2 \ \end{aligned} ]這里的關(guān)鍵技巧是湊出 ((x_i - \bar{x}{n-1})) 項(xiàng)。注意到 [ x_i - \bar{x}n (x_i - \bar{x}{n-1}) - (\bar{x}n - \bar{x}{n-1}) ] 而根據(jù)平均值遞推公式(\bar{x}n - \bar{x}{n-1} \frac{x_n - \bar{x}{n-1}}{n})。將其代入第一項(xiàng)求和 [ \begin{aligned} \sum_{i1}^{n-1} (x_i - \bar{x}n)^2 \sum{i1}^{n-1} \left[ (x_i - \bar{x}{n-1}) - \frac{x_n - \bar{x}{n-1}}{n} \right]^2 \ \sum_{i1}^{n-1} \left[ (x_i - \bar{x}{n-1})^2 - 2(x_i - \bar{x}{n-1})\frac{x_n - \bar{x}{n-1}}{n} \left( \frac{x_n - \bar{x}{n-1}}{n} \right)^2 \right] \ M_{n-1} - 2\frac{x_n - \bar{x}{n-1}}{n} \sum{i1}^{n-1}(x_i - \bar{x}{n-1}) (n-1)\left( \frac{x_n - \bar{x}{n-1}}{n} \right)^2 \end{aligned} ]這里有一個非常重要的性質(zhì)(\sum_{i1}^{n-1}(x_i - \bar{x}_{n-1}) 0)。因?yàn)槠钪蜑榱闶瞧骄档亩x性質(zhì)。所以中間項(xiàng)消失了于是 [ \sum_{i1}^{n-1} (x_i - \bar{x}n)^2 M{n-1} (n-1)\left( \frac{x_n - \bar{x}_{n-1}}{n} \right)^2 ]現(xiàn)在處理第二項(xiàng) ((x_n - \bar{x}n)^2)。同樣利用平均值遞推 [ x_n - \bar{x}n x_n - \left( \bar{x}{n-1} \frac{x_n - \bar{x}{n-1}}{n} \right) (x_n - \bar{x}{n-1}) - \frac{x_n - \bar{x}{n-1}}{n} \frac{n-1}{n}(x_n - \bar{x}_{n-1}) ] 所以 [ (x_n - \bar{x}n)^2 \left( \frac{n-1}{n} \right)^2 (x_n - \bar{x}{n-1})^2 ]最后將兩部分相加得到 (M_n) 的最終遞推公式 [ \begin{aligned} M_n \left[ M_{n-1} (n-1)\left( \frac{x_n - \bar{x}{n-1}}{n} \right)^2 \right] \left( \frac{n-1}{n} \right)^2 (x_n - \bar{x}{n-1})^2 \ M_{n-1} (x_n - \bar{x}{n-1})^2 \left[ \frac{n-1}{n^2} \frac{(n-1)^2}{n^2} \right] \ M{n-1} (x_n - \bar{x}{n-1})^2 \left[ \frac{n-1 (n-1)^2}{n^2} \right] \ M{n-1} (x_n - \bar{x}{n-1})^2 \left[ \frac{(n-1)n}{n^2} \right] \ M{n-1} \frac{(x_n - \bar{x}{n-1})(x_n - \bar{x}{n-1}) \cdot (n-1)}{n} \end{aligned} ]通常我們使用一個更緊湊、計算上更高效的等價形式。令差值 (delta x_n - \bar{x}{n-1})那么 [ \begin{aligned} \bar{x}n \bar{x}{n-1} \frac{delta}{n} \ M_n M{n-1} delta \cdot (x_n - \bar{x}_n) \end{aligned} ] 第二個公式可以通過將 (x_n - \bar{x}_n \frac{n-1}{n} delta) 代入上面的最終表達(dá)式驗(yàn)證得到。這個形式就是Welford算法在代碼中實(shí)際使用的核心。提示這個推導(dǎo)過程雖然有些繁瑣但它清晰地展示了算法如何通過巧妙的數(shù)學(xué)變換將可能引發(fā)精度抵消的大數(shù)相減運(yùn)算轉(zhuǎn)化為相對溫和的乘法與加法。核心在于利用了“前n-1個數(shù)據(jù)與其自身平均值偏差之和為零”這一性質(zhì)消去了不穩(wěn)定的交叉項(xiàng)。2.3 算法的數(shù)值穩(wěn)定性分析為什么這個遞推公式是穩(wěn)定的觀察 (M_n) 的更新項(xiàng)(delta \cdot (x_n - \bar{x}_n))。(delta x_n - \bar{x}_{n-1})是新樣本與舊平均值的偏差。這個值通常不會特別巨大。(x_n - \bar{x}_n)是新樣本與新平均值的偏差。根據(jù)推導(dǎo)它等于 (\frac{n-1}{n} delta)也是一個大小可控的值。更新量是這兩個“適中”大小的數(shù)的乘積然后累加到 (M_{n-1}) 上。整個過程中我們從未進(jìn)行過“兩個極大數(shù)相減”的操作。(M) 本身是平方和的累積雖然會增長但每次更新都是增加一個正數(shù)因?yàn)閮蓚€偏差通常同號避免了災(zāi)難性的抵消。這就是其數(shù)值穩(wěn)定性的根源。3. 算法實(shí)現(xiàn)詳解與代碼剖析理解了數(shù)學(xué)原理實(shí)現(xiàn)就變得直觀了。我們將分別用Python和C語言實(shí)現(xiàn)基礎(chǔ)版本并探討如何擴(kuò)展功能如合并數(shù)據(jù)集、計算協(xié)方差等。3.1 Python實(shí)現(xiàn)清晰與實(shí)用并重Python實(shí)現(xiàn)注重可讀性和實(shí)用性。我們將它封裝成一個類方便維護(hù)狀態(tài)。class Welford: 使用Welford算法在線計算均值、方差和標(biāo)準(zhǔn)差。 數(shù)值穩(wěn)定單次遍歷內(nèi)存消耗O(1)。 def __init__(self): self.count 0 # 數(shù)據(jù)點(diǎn)個數(shù) n self.mean 0.0 # 當(dāng)前平均值 x_bar self.M2 0.0 # 校正平方和 M_n def update(self, x): 加入一個新的數(shù)據(jù)點(diǎn) x。 self.count 1 delta x - self.mean self.mean delta / self.count delta2 x - self.mean # 注意這里用的是更新后的mean self.M2 delta * delta2 def variance(self, ddof1): 計算方差。 Args: ddof: 自由度增量Delta Degrees of Freedom。默認(rèn)為1計算樣本方差。 若為0則計算總體方差。 if self.count ddof 1: return float(nan) # 或返回0取決于你的需求 return self.M2 / (self.count - ddof) def std(self, ddof1): 計算標(biāo)準(zhǔn)差。 import math var self.variance(ddof) return math.sqrt(var) if var 0 else float(nan) def merge(self, other): 合并另一個Welford統(tǒng)計器的狀態(tài)。 用于并行計算后合并結(jié)果非常強(qiáng)大。 if other.count 0: return self if self.count 0: self.count other.count self.mean other.mean self.M2 other.M2 return self total_count self.count other.count delta other.mean - self.mean # 合并均值 self.mean (self.count * self.mean other.count * other.mean) / total_count # 合并M2。這是合并公式的核心它考慮了兩個數(shù)據(jù)集均值不同帶來的額外偏差。 self.M2 self.M2 other.M2 delta**2 * (self.count * other.count / total_count) self.count total_count return self def __repr__(self): return fWelford(n{self.count}, mean{self.mean:.6f}, var{self.variance():.6f}, std{self.std():.6f}) # 使用示例 if __name__ __main__: import numpy as np # 生成一些測試數(shù)據(jù)包含一個較大的常數(shù)偏移以凸顯傳統(tǒng)方法的問題 data 1e9 np.random.randn(1000000) # 100萬數(shù)據(jù)均值為1e9附近 # 使用Welford算法 w Welford() for x in data: w.update(x) print(Welford 結(jié)果:, w) # 使用NumPy直接計算兩遍法可能不穩(wěn)定 np_var np.var(data, ddof1) np_std np.std(data, ddof1) print(NumPy 方差:, np_var, 標(biāo)準(zhǔn)差:, np_std) # 使用不穩(wěn)定的單遍公式驗(yàn)證 naive_mean data.mean() # 下面這行在數(shù)據(jù)量巨大且均值很大時會因精度問題得到錯誤結(jié)果 unstable_var (np.sum(data**2) - len(data) * naive_mean**2) / (len(data) - 1) print(不穩(wěn)定單遍法方差:, unstable_var) print(差值Welford - NumPy:, w.variance() - np_var)代碼要點(diǎn)與心得update函數(shù)中的細(xì)節(jié)注意計算delta2時使用的是更新后的self.mean。這與我們推導(dǎo)的公式M_n M_{n-1} delta * (x_n - \bar{x}_n)完全一致。有些實(shí)現(xiàn)會先保存舊的均值計算完delta2后再更新均值也是等價的。ddof參數(shù)這是統(tǒng)計庫中常見的參數(shù)代表自由度增量。ddof1對應(yīng)樣本方差無偏估計ddof0對應(yīng)總體方差。在實(shí)現(xiàn)中預(yù)留這個參數(shù)讓函數(shù)更通用。merge方法這是Welford算法一個極其強(qiáng)大的特性。它允許你將數(shù)據(jù)分成多個塊分別用不同的Welford對象計算例如在多線程或分布式環(huán)境中最后將它們的結(jié)果正確合并。合并公式的推導(dǎo)同樣精妙它確保了合并后的統(tǒng)計量與一次性處理全部數(shù)據(jù)的結(jié)果在數(shù)學(xué)上完全一致。初始化count從0開始mean和M2初始為0。在update中當(dāng)count為1時delta/1就是數(shù)據(jù)本身mean被設(shè)置為第一個數(shù)據(jù)M2為0邏輯是自洽的。3.2 C語言實(shí)現(xiàn)追求效率與嵌入性C語言實(shí)現(xiàn)更注重效率和內(nèi)存控制適用于嵌入式系統(tǒng)或高性能計算庫。#include math.h #include stdio.h typedef struct { long long count; // 使用long long防止溢出 double mean; double M2; } welford_t; void welford_init(welford_t *w) { w-count 0; w-mean 0.0; w-M2 0.0; } void welford_update(welford_t *w, double x) { double delta, delta2; w-count; delta x - w-mean; w-mean delta / w-count; delta2 x - w-mean; // 使用新的均值 w-M2 delta * delta2; } double welford_variance(const welford_t *w, int ddof) { if (w-count ddof) { return NAN; // 需要C99標(biāo)準(zhǔn)或math.h支持 } return w-M2 / (w-count - ddof); } double welford_std(const welford_t *w, int ddof) { double var welford_variance(w, ddof); return sqrt(var); } void welford_merge(welford_t *dest, const welford_t *src) { if (src-count 0) return; if (dest-count 0) { *dest *src; return; } long long total_count dest-count src-count; double delta src-mean - dest-mean; double new_mean (dest-count * dest-mean src-count * src-mean) / total_count; // 合并M2 dest-M2 dest-M2 src-M2 delta * delta * (dest-count * src-count) / total_count; dest-mean new_mean; dest-count total_count; } // 示例從文件流中讀取數(shù)據(jù)并計算 void process_stream(FILE *stream) { welford_t stat; welford_init(stat); double value; while (fscanf(stream, %lf, value) 1) { welford_update(stat, value); } if (stat.count 1) { printf(數(shù)據(jù)點(diǎn)數(shù): %lld\n, stat.count); printf(平均值: %.6f\n, stat.mean); printf(樣本方差: %.6f\n, welford_variance(stat, 1)); printf(樣本標(biāo)準(zhǔn)差: %.6f\n, welford_std(stat, 1)); } }C實(shí)現(xiàn)注意事項(xiàng)整數(shù)溢出count使用long long類型以防處理海量數(shù)據(jù)時溢出。在32位系統(tǒng)上int類型最多表示約21億對于現(xiàn)代大數(shù)據(jù)場景可能不夠。浮點(diǎn)數(shù)類型使用double保證精度。在內(nèi)存極其受限的嵌入式環(huán)境如果數(shù)據(jù)范圍已知且較小可考慮float但要注意精度損失風(fēng)險。NAN的使用當(dāng)數(shù)據(jù)不足時返回NANNot a Number是標(biāo)準(zhǔn)做法。確保你的編譯環(huán)境支持C99math.h中定義了NAN。性能這個實(shí)現(xiàn)是O(1)時間和空間每個數(shù)據(jù)點(diǎn)只進(jìn)行幾次浮點(diǎn)運(yùn)算效率極高。循環(huán)中避免不必要的函數(shù)調(diào)用和分支。3.3 功能擴(kuò)展協(xié)方差與更復(fù)雜的統(tǒng)計量Welford算法的思想可以推廣到計算協(xié)方差和方差-協(xié)方差矩陣。這對于在線計算線性回歸、PCA等至關(guān)重要。在線協(xié)方差計算維護(hù)三個狀態(tài)量count,mean_x,mean_y,C2類似于M2但用于協(xié)方差。 更新公式為def update_cov(self, x, y): self.count 1 delta_x x - self.mean_x delta_y y - self.mean_y self.mean_x delta_x / self.count self.mean_y delta_y / self.count # 注意這里用的新均值 self.C2 delta_x * (y - self.mean_y) # 等價于 delta_y * (x - self.mean_x)數(shù)值上略有差異但都穩(wěn)定。協(xié)方差即為C2 / (count - ddof)。通過維護(hù)多個這樣的對象可以構(gòu)建一個在線更新的協(xié)方差矩陣。4. 實(shí)戰(zhàn)場景、常見問題與性能對比理解了原理和實(shí)現(xiàn)我們來看看它在哪能大顯身手以及實(shí)際使用中會遇到哪些坑。4.1 典型應(yīng)用場景流式數(shù)據(jù)監(jiān)控服務(wù)器每秒產(chǎn)生數(shù)萬條性能指標(biāo)延遲、CPU使用率。使用Welford算法你只需要維護(hù)幾個狀態(tài)變量就可以實(shí)時更新并報告最近N秒內(nèi)的平均值和標(biāo)準(zhǔn)差用于觸發(fā)告警。內(nèi)存占用恒定與數(shù)據(jù)量無關(guān)。大規(guī)模數(shù)據(jù)文件處理處理一個幾十GB的CSV文件無法一次性加載到內(nèi)存。你可以分塊讀取每塊用一個Welford對象更新或者并行處理多個塊最后合并。這比將數(shù)據(jù)全部讀入再調(diào)用np.var()要節(jié)省內(nèi)存得多且數(shù)值更穩(wěn)定。嵌入式傳感器系統(tǒng)MCU內(nèi)存只有幾十KB但需要持續(xù)采集傳感器數(shù)據(jù)并計算其波動性標(biāo)準(zhǔn)差。Welford算法是唯一可行的選擇。在線機(jī)器學(xué)習(xí)一些在線學(xué)習(xí)算法如在線標(biāo)準(zhǔn)化、某些隨機(jī)梯度下降變體需要動態(tài)更新特征的平均值和標(biāo)準(zhǔn)差。Welford算法可以無縫集成到這些迭代過程中。數(shù)據(jù)庫聚合雖然數(shù)據(jù)庫內(nèi)置了VAR_POP、VAR_SAMP函數(shù)但如果你需要在應(yīng)用層進(jìn)行復(fù)雜的、多階段的聚合計算Welford的merge操作會非常有用。4.2 常見陷阱與排查指南即使算法本身健壯實(shí)現(xiàn)和使用不當(dāng)也會出問題。問題1方差或標(biāo)準(zhǔn)差為負(fù)數(shù)或NaN。原因這幾乎總是數(shù)值精度問題但發(fā)生在Welford算法中可能有兩個原因數(shù)據(jù)量太少當(dāng)count ddof時分母為零或負(fù)計算無意義。實(shí)現(xiàn)時必須做防御性檢查。極端情況下的舍入誤差雖然Welford很穩(wěn)定但在極端情況下如數(shù)據(jù)量極大M2累積了巨大的值而更新量極小M2理論上應(yīng)單調(diào)遞增但浮點(diǎn)舍入可能導(dǎo)致其略微減小。當(dāng)M2為一個極小的負(fù)數(shù)時開平方會得到NaN。排查與解決檢查數(shù)據(jù)量確保count ddof。在variance()函數(shù)中返回前可以加一個保護(hù)return max(0.0, self.M2) / (self.count - ddof)。方差本質(zhì)是非負(fù)的將極小的負(fù)數(shù)截斷為0是實(shí)踐中可接受的做法?？紤]使用更高精度的浮點(diǎn)數(shù)如Python的decimal.Decimal或C的long double但會犧牲性能。問題2合并merge操作后結(jié)果與預(yù)期不符。原因合并公式使用錯誤或者合并順序有誤。合并操作不滿足交換律嗎實(shí)際上Welford的合并操作在數(shù)學(xué)上是滿足交換律和結(jié)合律的但浮點(diǎn)數(shù)運(yùn)算的舍入誤差可能導(dǎo)致微小的差異。排查用一個小數(shù)據(jù)集如[1,2,3]和[4,5]手動計算合并過程與你的代碼輸出對比。確保合并時delta計算正確是src.mean - dest.mean。合并大量對象時為了最小化誤差可以考慮將所有對象兩兩合并而不是依次合并到一個對象上。但這通常差異極小。問題3與權(quán)威庫如NumPy的結(jié)果存在微小差異。原因這是正?，F(xiàn)象。差異主要來自算法不同NumPy默認(rèn)使用兩遍法。對于病態(tài)數(shù)據(jù)Welford的結(jié)果更可信。計算順序浮點(diǎn)數(shù)加法不滿足結(jié)合律求和順序不同會導(dǎo)致不同的舍入誤差。起始值你的Welford對象從0開始更新而NumPy是一次性計算所有數(shù)據(jù)。行動對于“合理”的數(shù)據(jù)差異通常在機(jī)器精度1e-15或1e-12級別。如果差異巨大首先用我們示例中的“不穩(wěn)定單遍法”驗(yàn)證看是否是NumPy的兩遍法本身在你的數(shù)據(jù)上就出現(xiàn)了精度問題。Welford的結(jié)果通常更接近“真實(shí)”的數(shù)值。問題4在多線程環(huán)境下更新同一個Welford對象。警告上面給出的基礎(chǔ)實(shí)現(xiàn)不是線程安全的。count 1和mean、M2的更新不是原子操作。解決方案每個線程維護(hù)自己的本地Welford對象處理完一部分?jǐn)?shù)據(jù)后再在同步階段用merge合并到全局對象中。這是推薦做法既避免了鎖競爭又利用了merge的強(qiáng)大功能。如果必須并發(fā)更新需要對整個update函數(shù)加鎖如互斥鎖但這會嚴(yán)重降低性能。4.3 性能與精度對比實(shí)驗(yàn)我們設(shè)計一個簡單的實(shí)驗(yàn)來直觀感受其優(yōu)勢。用Python生成一組“壞”數(shù)據(jù)data 1e12 np.random.randn(10_000_000)即一千萬個圍繞1萬億波動的數(shù)。精度對比Welford算法方差 ≈ 0.999理論應(yīng)為1因我們用了標(biāo)準(zhǔn)正態(tài)分布。NumPy兩遍法np.var(data, ddof1)結(jié)果也接近0.999在這個量級下NumPy的內(nèi)部實(shí)現(xiàn)可能已經(jīng)采用了某種穩(wěn)定算法或高精度累加器。不穩(wěn)定單遍法(np.sum(data**2) - len(data)*mean**2)/(len(data)-1)。這里你會看到問題np.sum(data**2)和len(data)*mean**2都是大約1e25量級的數(shù)它們的差應(yīng)該在1e13量級即約1e7 * 1e6。實(shí)際上由于浮點(diǎn)數(shù)只有約15位十進(jìn)制有效數(shù)字這個減法會丟失所有有效數(shù)字結(jié)果可能是0、負(fù)數(shù)或者一個完全錯誤的巨大數(shù)值。在我的測試中它返回了一個毫無意義的巨大正數(shù)。內(nèi)存與速度內(nèi)存Welford算法無論處理多少數(shù)據(jù)只占用常數(shù)內(nèi)存3個浮點(diǎn)數(shù)1個整數(shù)。而兩遍法需要存儲整個數(shù)組對于一千萬個float64就是80MB。速度對于一千萬數(shù)據(jù)單線程遍歷更新Welford單遍更新約0.5秒。NumPy兩遍法需要先創(chuàng)建數(shù)組創(chuàng)建數(shù)組計算約0.3秒。注意NumPy底層是高度優(yōu)化的C代碼且利用了SIMD指令所以單純計算可能更快。但Welford的O(1)內(nèi)存和單遍特性在很多場景下是不可替代的優(yōu)勢。實(shí)操心得不要盲目認(rèn)為Welford在所有情況下都“更快”。它的核心優(yōu)勢在于流式處理、內(nèi)存效率和數(shù)值穩(wěn)定性。當(dāng)數(shù)據(jù)能全部裝入內(nèi)存且使用像NumPy這樣的優(yōu)化庫時直接調(diào)用內(nèi)置函數(shù)可能更方便更快。但當(dāng)數(shù)據(jù)是“流”或“太大”時Welford是唯一的王者。在實(shí)現(xiàn)時將update函數(shù)內(nèi)聯(lián)C/C中或使用JIT編譯如Numba裝飾Python函數(shù)可以極大提升其性能使其在速度上也極具競爭力。5. 總結(jié)與進(jìn)階思考Welford算法是一個將深刻數(shù)學(xué)思想轉(zhuǎn)化為簡潔、魯棒代碼的典范。它解決了工程中的一個實(shí)際痛點(diǎn)——穩(wěn)定地在線計算二階統(tǒng)計量。通過本次推導(dǎo)和實(shí)現(xiàn)我們不僅掌握了一個工具更學(xué)習(xí)了一種“遞推”和“在線更新”的思維方式。我個人在實(shí)際使用中的體會是它常常作為更復(fù)雜算法的基石。比如在實(shí)現(xiàn)在線標(biāo)準(zhǔn)化層Online Normalization時我需要同時維護(hù)每個特征維度的Welford統(tǒng)計器在實(shí)現(xiàn)簡單的異常檢測如3-sigma法則時用Welford實(shí)時更新閾值比滑動窗口更節(jié)省內(nèi)存。它的merge特性讓我能夠輕松地將分布在不同機(jī)器或線程上的統(tǒng)計結(jié)果匯總這在設(shè)計分布式數(shù)據(jù)預(yù)處理流水線時非常有用。最后再分享一個小技巧如果你需要計算帶權(quán)重的方差例如每個數(shù)據(jù)點(diǎn)有一個權(quán)重Welford算法也有相應(yīng)的加權(quán)版本其推導(dǎo)思路類似但維護(hù)的是加權(quán)計數(shù)、加權(quán)均值和加權(quán)的校正平方和。當(dāng)你從簡單的等權(quán)重統(tǒng)計邁向更復(fù)雜的加權(quán)統(tǒng)計時你會再次欣賞到這種遞推公式的優(yōu)美和強(qiáng)大。