簡(jiǎn)介GP-EnKF是一份基于Python實(shí)現(xiàn)的在線高斯過程回歸算法代碼源自Fusion 2018論文所述方法面向需要處理流式數(shù)據(jù)實(shí)時(shí)預(yù)測(cè)與不確定性估計(jì)的研究者、工程師與算法學(xué)習(xí)者。該方案將高斯過程回歸與集合卡爾曼濾波EnKF相結(jié)合通過狀態(tài)集合的預(yù)測(cè)與更新步驟交替迭代緩解了傳統(tǒng)高斯過程回歸在數(shù)據(jù)規(guī)模增大時(shí)計(jì)算復(fù)雜度快速上升的問題適用于環(huán)境科學(xué)、控制工程、信號(hào)處理等動(dòng)態(tài)系統(tǒng)在線監(jiān)測(cè)場(chǎng)景。資源包為zip壓縮格式整體約22KB內(nèi)含可運(yùn)行的Python核心腳本可直觀理解高斯過程先驗(yàn)設(shè)定、EnKF狀態(tài)集合構(gòu)建與觀測(cè)更新融合的完整流程并便于對(duì)照論文復(fù)現(xiàn)在線學(xué)習(xí)效果再遷移到自身流式數(shù)據(jù)任務(wù)中開展預(yù)測(cè)與不確定性分析。壓縮包內(nèi)具體文件構(gòu)成暫未顯示。目前已有305人瀏覽學(xué)習(xí)適合具備Python與概率模型基礎(chǔ)、希望快速掌握高斯過程與濾波融合在線回歸方法的讀者。1. 在線高斯過程回歸的高成本困局GP-EnKF把O(n3)變成在線更新批量高斯過程回歸每次加入新數(shù)據(jù)都要重算 n×n 核矩陣的逆復(fù)雜度隨數(shù)據(jù)量三次方增長(zhǎng)。數(shù)據(jù)量一上三千單次更新就能把機(jī)器卡到懷疑人生在線數(shù)據(jù)流場(chǎng)景下根本跑不動(dòng)。GP-EnKF 的思路是先讓歸納點(diǎn)把訓(xùn)練數(shù)據(jù)壓縮成 m 個(gè)偽樣本再用集合卡爾曼濾波器在數(shù)據(jù)到達(dá)時(shí)同步更新歸納點(diǎn)和超參數(shù)本身——既保留 GP 的不確定性估計(jì)能力又把單步更新降到只與歸納點(diǎn)數(shù)量相關(guān)的量級(jí)。這篇筆記拆的是 Fusion 2018 論文的配套 Python 代碼從原理講到復(fù)現(xiàn)參數(shù)再落到避坑。適合正在做在線預(yù)測(cè)、流式數(shù)據(jù)建模、以及需要預(yù)測(cè)方差而不是只看點(diǎn)估計(jì)的從業(yè)者。2. 從批量GP到EnKF歸納點(diǎn)與狀態(tài)估計(jì)的核心原理2.1 批量GP的O(n3)瓶頸與在線化矛盾高斯過程回歸的預(yù)測(cè)式寫出來很漂亮均值是 m(x)k(x,X)[K(X,X)σ2I]?1y方差是 v(x)k(x,x)?k(x,X)[K(X,X)σ2I]?1k(X,x)。但漂亮背后有個(gè)殘酷的現(xiàn)實(shí)——每次新觀測(cè)到達(dá)都要對(duì) K(X,X)σ2I 做一次 Cholesky 分解或求逆。n 從兩千漲到四千計(jì)算量直接翻八倍這在流式數(shù)據(jù)場(chǎng)景里不可接受。行業(yè)里通常有兩條替代路線。一條是稀疏近似Sparse GP用 m 個(gè)歸納點(diǎn)替代 n 個(gè)訓(xùn)練點(diǎn)把復(fù)雜度降到 O(nm2)另一條是遞歸濾波把超參數(shù)當(dāng)常數(shù)用 Kalman 類方法更新后驗(yàn)。但兩條路線各有各的坑稀疏 GP 把歸納點(diǎn)固定在初始化位置的話數(shù)據(jù)分布一漂移預(yù)測(cè)立刻崩遞歸濾波則要手動(dòng)推導(dǎo)協(xié)方差傳播公式觀測(cè)模型稍微非線性就推不動(dòng)。這兩條路線的共同盲區(qū)是歸納點(diǎn)放哪、超參數(shù)取多少在線場(chǎng)景下其實(shí)是動(dòng)態(tài)量。數(shù)據(jù)分布會(huì)漂移最優(yōu)長(zhǎng)度尺度會(huì)變化把這些東西當(dāng)常數(shù)處理等于假設(shè)世界不變。GP-EnKF 的出發(fā)點(diǎn)就是把這個(gè)盲區(qū)正面解決——把歸納點(diǎn)位置、對(duì)應(yīng)函數(shù)值 u、核函數(shù)超參數(shù)全部塞進(jìn)一個(gè)狀態(tài)向量用 EnKF 做聯(lián)合估計(jì)。數(shù)據(jù)到達(dá)時(shí)更新的不是某一個(gè)值而是整個(gè)狀態(tài)的分布。2.2 歸納點(diǎn)把n個(gè)數(shù)據(jù)壓縮成m個(gè)偽樣本歸納點(diǎn)的思想可以追溯到 Sparse GP。假設(shè)有 n 個(gè)訓(xùn)練點(diǎn)我們選出 m 個(gè)偽輸入 Z{z?,...,z_m}再用這 m 個(gè)點(diǎn)上的函數(shù)值 u 來近似完整的 GP 后驗(yàn)。關(guān)鍵推導(dǎo)是如果 u 的先驗(yàn)是 GP(0, K(Z,Z))那么給定 u 時(shí)任意測(cè)試點(diǎn) x 的預(yù)測(cè)服從高斯分布均值是 K(x,Z)K(Z,Z)?1u方差是 k(x,x)?K(x,Z)K(Z,Z)?1K(Z,x)。注意這里完全沒有 n 參與計(jì)算計(jì)算量只取決于 m。m 怎么取我一般看輸入維度定一維問題 5–15 個(gè)點(diǎn)就夠二維至少 20維度再高建議從 30 起步。取太少擬合不了非線性取太多就失去了在線更新的意義。歸納點(diǎn)初始位置可以用 k-means 聚類中心也可以直接在輸入范圍內(nèi)均勻撒點(diǎn)。GP-EnKF 的好處是后面這些點(diǎn)會(huì)自己動(dòng)——數(shù)據(jù)來了它自動(dòng)調(diào)整位置這是靜態(tài)歸納點(diǎn)方案給不了的。實(shí)際工程里我會(huì)額外注意歸納點(diǎn)的排序約束。一維輸入時(shí) Z 必須保持單調(diào)否則核矩陣 K(Z,Z) 的相鄰行可能會(huì)因?yàn)閮蓚€(gè)歸納點(diǎn)距離過近而近乎線性相關(guān)直接導(dǎo)致矩陣奇異。這個(gè)問題在后面的避坑章節(jié)還會(huì)詳細(xì)展開。2.3 EnKF憑什么能做GP狀態(tài)估計(jì)EnKF 是集合卡爾曼濾波器的縮寫核心是用一組 ensemble 粒子比如 50 個(gè)狀態(tài)向量近似后驗(yàn)分布用粒子的樣本協(xié)方差代替解析協(xié)方差。它比粒子濾波簡(jiǎn)單——不需要重要性重采樣不存在權(quán)重退化問題又比標(biāo)準(zhǔn) Kalman 對(duì)非線性觀測(cè)模型的寬容度高得多。GP 的觀測(cè)模型恰好是非線性的預(yù)測(cè)觀測(cè) \hat{y}K(x,Z)K(Z,Z)?1u這個(gè)映射關(guān)于 u 是線性的但關(guān)于 Z 是非線性的。EnKF 的觀測(cè)擾動(dòng)形式在這里特別自然對(duì)每個(gè) ensemble 成員把 y_obs 加上 N(0,σ_n2) 的隨機(jī)擾動(dòng)然后計(jì)算 Kalman 增益并更新狀態(tài)。更新的核心公式是s_i^a s_i^f K (y_obs ε_(tái)i ? \hat{y}_i)其中 K 由 ensemble 樣本協(xié)方差估計(jì)K P H? (H P H? R)?1這里 P 是狀態(tài)向量的樣本協(xié)方差H 是觀測(cè)對(duì)狀態(tài)的敏感度用 ensemble 統(tǒng)計(jì)量近似R 是觀測(cè)噪聲方差。還有一個(gè)容易被忽略的優(yōu)勢(shì)EnKF 更新后得到的是 ensemble不是一個(gè)點(diǎn)估計(jì)。預(yù)測(cè)的不確定性直接由 ensemble 的離散度給出不用額外推導(dǎo)協(xié)方差傳播公式。這意味著工程上我們不需要維護(hù)復(fù)雜的解析協(xié)方差方差是統(tǒng)計(jì)出來的不是推出來的——這個(gè)特性讓 GP-EnKF 在新場(chǎng)景里落地特別快。3. GP-EnKF的Python實(shí)現(xiàn)狀態(tài)向量、預(yù)測(cè)步與更新步3.1 狀態(tài)向量設(shè)計(jì)與ensemble初始化這里給出核心代碼完整實(shí)現(xiàn)按 Fusion 2018 論文的思路來。先定義 RBF 核函數(shù)注意一維輸入專用版本避免多維廣播的細(xì)節(jié)干擾import numpy as np def rbf_1d(X1, X2, ell, sigma_f): 一維輸入的RBF核矩陣 X1, X2: 一維數(shù)組或類數(shù)組 ell: 長(zhǎng)度尺度, sigma_f: 信號(hào)標(biāo)準(zhǔn)差 返回形狀 (len(X1), len(X2)) 的核矩陣 X1 np.atleast_1d(X1).reshape(-1, 1) X2 np.atleast_1d(X2).reshape(-1, 1) dist2 (X1 - X2.T) ** 2 # 廣播成 (n, m) 的平方距離矩陣 return sigma_f**2 * np.exp(-0.5 * dist2 / ell**2)核函數(shù)里 dist2 的廣播是關(guān)鍵——X1 和 X2 先都 reshape 成列向量然后相減得到二維距離矩陣。這個(gè)寫法比雙重循環(huán)快一個(gè)數(shù)量級(jí)而且代碼量少。注意ell和sigma_f是標(biāo)量參數(shù)。接下來定義 GP-EnKF 類。狀態(tài)向量設(shè)計(jì)為 s[u?...u_m, log ?, log σ_f, z?...z_m]一維輸入時(shí)維度是 2m2。歸納點(diǎn)函數(shù)值 u 是核心狀態(tài)量超參數(shù)取 log 空間是為了保證更新后不會(huì)變成負(fù)值class GPEnKF: def __init__(self, n_ens50, n_inducing10, sigma_n0.05, ell_init1.0, sigma_f_init1.0, process_noise0.01): self.n_ens n_ens # ensemble成員數(shù) self.m n_inducing # 歸納點(diǎn)數(shù) self.sigma_n sigma_n # 觀測(cè)噪聲標(biāo)準(zhǔn)差 self.process_noise process_noise # 過程噪聲尺度 # 歸納點(diǎn)初始位置輸入范圍內(nèi)均勻撒點(diǎn), 所有成員共享 Z0 np.linspace(-2, 2, n_inducing) self.Z np.tile(Z0, (n_ens, 1)) # 形狀 (n_ens, m) # 歸納點(diǎn)函數(shù)值從N(0,1)采樣, 每個(gè)成員獨(dú)立 self.u np.random.randn(n_ens, n_inducing) # (n_ens, m) # 超參數(shù)在log空間采樣, 加少量擾動(dòng)讓ensemble有初始散布 self.log_ell np.log(ell_init) 0.1 * np.random.randn(n_ens) self.log_sf np.log(sigma_f_init) 0.1 * np.random.randn(n_ens)初始化里np.tile讓所有 ensemble 成員共享初始?xì)w納點(diǎn)位置但 u 和超參數(shù)各自有獨(dú)立擾動(dòng)。這個(gè)設(shè)計(jì)保證初始 ensemble 有足夠的多樣性避免一開始就塌縮成一個(gè)點(diǎn)。sigma_n是最敏感的參數(shù)它決定了更新步的信噪比——設(shè)太小模型會(huì)狂追噪聲設(shè)太大預(yù)測(cè)會(huì)過度平滑。3.2 預(yù)測(cè)步隨機(jī)游走與過程噪聲EnKF 的預(yù)測(cè)步在 GP 場(chǎng)景里沒有物理模型驅(qū)動(dòng)所以用隨機(jī)游走近似。邏輯是狀態(tài)量在兩次觀測(cè)之間有小幅隨機(jī)漂移漂移幅度由 process_noise 控制def predict_step(self): 狀態(tài)演化: 歸納值隨機(jī)游走, 超參數(shù)微擾 # 歸納點(diǎn)函數(shù)值按隨機(jī)游走演化, 幅度正比于過程噪聲 self.u self.process_noise * np.random.randn(*self.u.shape) # 歸納點(diǎn)位置慢速漂移, 只允許小步移動(dòng) self.Z 0.001 * np.random.randn(*self.Z.shape) # 超參數(shù)在log空間小步擾動(dòng), 保證非負(fù)性 self.log_ell 0.005 * np.random.randn(self.n_ens) self.log_sf 0.005 * np.random.randn(self.n_ens)這里三個(gè)隨機(jī)游走的幅度是有講究的。歸納點(diǎn)函數(shù)值 u 的擾動(dòng)幅度process_noise設(shè)為 0.01代表狀態(tài)在相鄰兩步之間的先驗(yàn)不確定性歸納點(diǎn)位置擾動(dòng)是 0.001比 u 小一個(gè)量級(jí)防止 Z 漂移太快導(dǎo)致核矩陣形狀劇變超參數(shù)擾動(dòng) 0.005 控制在 log 空間相當(dāng)于每次最多變化 0.5%。有朋友會(huì)問為什么不直接把 process_noise 設(shè)成 0那樣狀態(tài)就完全確定更新步會(huì)退化成確定性映射ensemble 方差持續(xù)縮小最后徹底塌縮。過程噪聲的本質(zhì)是給 ensemble 持續(xù)注入不確定性讓濾波器保持可被新數(shù)據(jù)修正的狀態(tài)。這個(gè)參數(shù)在非平穩(wěn)數(shù)據(jù)上尤其重要——它本質(zhì)上告訴了濾波器世界在變你要跟得上。3.3 更新步EnKF分析公式與代碼更新步是整套實(shí)現(xiàn)的核心。流程分四段先算每個(gè)成員的預(yù)測(cè)觀測(cè)再組裝狀態(tài)矩陣并估計(jì)樣本協(xié)方差然后算 Kalman 增益做協(xié)方差膨脹最后更新每個(gè)成員的狀態(tài)def update_step(self, x_obs, y_obs, inflation1.05): EnKF分析步: 用當(dāng)前觀測(cè)更新每個(gè)ensemble成員 x_obs: 當(dāng)前觀測(cè)的輸入(標(biāo)量) y_obs: 當(dāng)前觀測(cè)的目標(biāo)值(標(biāo)量) inflation: 協(xié)方差膨脹因子, 防止方差塌縮 n self.n_ens m self.m ell np.exp(self.log_ell) sf np.exp(self.log_sf) # 1. 預(yù)測(cè)觀測(cè): 對(duì)每個(gè)成員計(jì)算 \hat{y}_i K(x,Z)K(Z,Z)^{-1}u H_ens np.zeros(n) for i in range(n): Kzz rbf_1d(self.Z[i], self.Z[i], ell[i], sf[i]) 1e-6 * np.eye(m) Kxz rbf_1d(np.array([x_obs]), self.Z[i], ell[i], sf[i]) # 用solve替代inv, 數(shù)值更穩(wěn)定 H_ens[i] Kxz np.linalg.solve(Kzz, self.u[i]) # 2. 組裝狀態(tài)矩陣并估計(jì)統(tǒng)計(jì)量 state np.hstack([self.u, self.log_ell[:, None], self.log_sf[:, None], self.Z]) state_mean state.mean(axis0) H_mean H_ens.mean() # 協(xié)方差膨脹: 把ensemble圍繞均值拉開, 抵消更新步的方差收縮 state state_mean np.sqrt(inflation) * (state - state_mean) # 3. Kalman增益: 標(biāo)量觀測(cè)時(shí)退化為向量形式 # PH_T Cov(state, \hat{y}), HPH_R Var(\hat{y}) sigma_n^2 PH_T ((state - state_mean).T (H_ens - H_mean)) / (n - 1) HPH_R np.sum((H_ens - H_mean)**2) / (n - 1) self.sigma_n**2 K PH_T / HPH_R # 形狀 (state_dim,) # 4. 觀測(cè)擾動(dòng) 更新 y_perturbed y_obs self.sigma_n * np.random.randn(n) for i in range(n): innovation y_perturbed[i] - H_ens[i] state[i] K * innovation # 拆回狀態(tài)分量 self.u state[:, :m] self.log_ell state[:, m] self.log_sf state[:, m1] self.Z state[:, m2:]這段代碼里值得注意幾個(gè)工程細(xì)節(jié)。np.linalg.solve(Kzz, self.u[i])替代np.linalg.inv(Kzz) self.u[i]前者用 LU 分解避免顯式求逆數(shù)值穩(wěn)定性好得多。Kzz對(duì)角線上加的1e-6是 jitter專門對(duì)付歸納點(diǎn)距離過近導(dǎo)致的近奇異矩陣。協(xié)方差膨脹放在統(tǒng)計(jì)量計(jì)算之前這是標(biāo)準(zhǔn) EnKF 流程——膨脹作用于預(yù)測(cè) ensemble而不是更新步之后膨脹后再計(jì)算 PH_T 和 HPH_R增益本身就包含了對(duì)塌縮的修正。觀測(cè)擾動(dòng)y_obs sigma_n * np.random.randn(n)是 EnKF 的隨機(jī)擾動(dòng)形式它保證了更新后的 ensemble 方差不會(huì)系統(tǒng)性偏小。如果你希望實(shí)現(xiàn)完全確定性的更新可以用平方根版本的 EnKFETKF但代碼復(fù)雜度會(huì)明顯上升一般場(chǎng)景沒有這個(gè)必要。3.4 在線預(yù)測(cè)從ensemble到后驗(yàn)分布預(yù)測(cè)時(shí)把每個(gè) ensemble 成員的歸納點(diǎn)信息代入 GP 預(yù)測(cè)式得到一組預(yù)測(cè)值再統(tǒng)計(jì)均值和方差def predict(self, x_query): 預(yù)測(cè)均值與方差 x_query: 查詢點(diǎn)(標(biāo)量) 返回: (均值, 方差), 方差包含觀測(cè)噪聲項(xiàng) ell np.exp(self.log_ell) sf np.exp(self.log_sf) preds np.zeros(self.n_ens) for i in range(self.n_ens): Kzz rbf_1d(self.Z[i], self.Z[i], ell[i], sf[i]) 1e-6 * np.eye(self.m) Kxz rbf_1d(np.array([x_query]), self.Z[i], ell[i], sf[i]) preds[i] Kxz np.linalg.solve(Kzz, self.u[i]) mean preds.mean() var preds.var() self.sigma_n**2 # ensemble方差 觀測(cè)噪聲 return mean, var預(yù)測(cè)方差由兩部分構(gòu)成ensemble 方差代表了模型對(duì)函數(shù)值的不確定性sigma_n**2是觀測(cè)噪聲。這個(gè)加法很重要——如果不加置信區(qū)間會(huì)系統(tǒng)性偏窄做不確定性量化時(shí)覆蓋率會(huì)明顯低于理論值。在線學(xué)習(xí)主循環(huán)很簡(jiǎn)潔。數(shù)據(jù)流持續(xù)進(jìn)入每步先 predict_step 再 update_step每隔若干步做一次評(píng)估# 在線學(xué)習(xí)循環(huán)示例: 300個(gè)數(shù)據(jù)點(diǎn), 每步更新一次 np.random.seed(42) X_stream np.sort(np.random.uniform(-5, 5, 300)) y_stream np.sin(X_stream) 0.05 * np.random.randn(300) model GPEnKF(n_ens50, n_inducing10, sigma_n0.05) rmse_list [] for t in range(300): model.predict_step() model.update_step(X_stream[t], y_stream[t]) # 每10步評(píng)估一次在固定測(cè)試點(diǎn)上的預(yù)測(cè)精度 if t 20 and t % 10 0: m, v model.predict(np.array([1.2])) rmse_list.append((m - np.sin(1.2))**2) print(Test RMSE:, np.sqrt(np.mean(rmse_list)))這個(gè)循環(huán)是 GP-EnKF 最基本的用法。300 個(gè)數(shù)據(jù)點(diǎn)全程在線更新沒有重新訓(xùn)練單步開銷取決于 n_ens 和 m 的乘積與累計(jì)數(shù)據(jù)量無關(guān)。實(shí)際場(chǎng)景里如果數(shù)據(jù)到達(dá)是批量突發(fā)一次來 50 條可以把循環(huán)改造成 mini-batch 形式——對(duì)一批數(shù)據(jù)逐條調(diào)用 update_step或者把批量觀測(cè)向量化后者需要把 H_ens 擴(kuò)展為矩陣形式。4. 參數(shù)設(shè)置與Fusion 2018復(fù)現(xiàn)要點(diǎn)4.1 影響精度的四個(gè)參數(shù)表參數(shù)設(shè)置直接影響收斂速度、預(yù)測(cè)精度和數(shù)值穩(wěn)定性。我把 Fusion 2018 論文里涉及的關(guān)鍵參數(shù)整理成一張表然后逐個(gè)說明選擇依據(jù)參數(shù)含義推薦范圍設(shè)置過小的后果設(shè)置過大的后果n_ensensemble成員數(shù)30–100樣本協(xié)方差噪聲大估計(jì)不穩(wěn)定計(jì)算量線性增長(zhǎng)收益遞減n_inducing歸納點(diǎn)數(shù)5–30一維擬合不了非線性結(jié)構(gòu)失去在線計(jì)算優(yōu)勢(shì)sigma_n觀測(cè)噪聲標(biāo)準(zhǔn)差0.01–0.1或數(shù)據(jù)噪聲的估計(jì)值模型狂追噪聲預(yù)測(cè)方差偏小過度平滑細(xì)節(jié)丟失process_noise過程噪聲0.001–0.05狀態(tài)演化過慢非平穩(wěn)數(shù)據(jù)滯后狀態(tài)抖動(dòng)大預(yù)測(cè)方差虛高inflation協(xié)方差膨脹因子1.0–1.1ensemble提前塌縮方差人為放大置信區(qū)間失真n_ens 是精度和速度的主要權(quán)衡項(xiàng)。50 是多數(shù)場(chǎng)景的甜點(diǎn)——樣本協(xié)方差已經(jīng)有足夠統(tǒng)計(jì)精度單步更新在普通筆記本上毫秒級(jí)完成。如果你的數(shù)據(jù)噪聲特別小協(xié)方差矩陣的條件數(shù)不好建議把 n_ens 提到 80 以上。n_inducing 的選擇邏輯不同一維平滑函數(shù) 5 個(gè)就夠帶多個(gè)波峰的函數(shù)要 10–15二維輸入至少 20 起步。Fusion 2018 論文的實(shí)驗(yàn)里一維基準(zhǔn)用了 10 個(gè)歸納點(diǎn)我在復(fù)現(xiàn)時(shí)發(fā)現(xiàn)這個(gè)值在大多數(shù)平滑函數(shù)上足夠但遇到劇烈振蕩的函數(shù)比如頻率超過 3 的正弦疊加需要追加到 15。sigma_n 是最容易翻車的參數(shù)。很多人在初始化時(shí)設(shè) 0.05但這個(gè)值必須和數(shù)據(jù)的真實(shí)噪聲水平匹配。一個(gè)可行的估計(jì)方式拿前 20 個(gè)數(shù)據(jù)算相鄰點(diǎn)差分的標(biāo)準(zhǔn)差再除以√2得到噪聲的粗略估計(jì)。初始化階段寧可從大到小調(diào)不要一開始就設(shè)成 0.001 這種值。4.2 初始化策略超參數(shù)先從數(shù)據(jù)里猜超參數(shù)初始化對(duì) GP-EnKF 的收斂速度影響巨大。隨機(jī)初始化不是好主意——長(zhǎng)度尺度差一個(gè)數(shù)量級(jí)核矩陣的形狀會(huì)完全不同EnKF 要花很多步才能把超參數(shù)拉回正軌。我一般按照下面的流程做初始化# 用前20個(gè)數(shù)據(jù)點(diǎn)估算超參數(shù)初始值 init_X X_stream[:20] init_y y_stream[:20] # 長(zhǎng)度尺度: 輸入范圍的1/4左右 ell_init (init_X.max() - init_X.min()) / 4.0 # 信號(hào)方差: 目標(biāo)值的方差 sigma_f_init np.sqrt(np.var(init_y)) # 觀測(cè)噪聲: 相鄰點(diǎn)差分標(biāo)準(zhǔn)差 / sqrt(2) diff_std np.std(np.diff(init_y)) sigma_n_init diff_std / np.sqrt(2.0) print(fell_init{ell_init:.3f}, sf_init{sigma_f_init:.3f}, sn_init{sigma_n_init:.3f})長(zhǎng)度尺度取輸入范圍的 1/4是為了保證初始核矩陣覆蓋大部分?jǐn)?shù)據(jù)點(diǎn)的相互作用。如果把長(zhǎng)度尺度設(shè)成輸入范圍的幾倍核函數(shù)會(huì)過于平滑前幾步的預(yù)測(cè)偏差會(huì)被 EnKF 放大。信號(hào)方差直接取目標(biāo)值方差這是一個(gè)無偏估計(jì)——GP 先驗(yàn)的邊際方差就是 σ_f2。觀測(cè)噪聲用相鄰點(diǎn)差分估計(jì)是時(shí)間序列里常用的小技巧假設(shè)相鄰點(diǎn)函數(shù)值接近差分主要由噪聲主導(dǎo)。4.3 訓(xùn)練與評(píng)估流程完整的訓(xùn)練評(píng)估流程按下面的順序走每一步都有明確的檢查點(diǎn)# 1. 加載數(shù)據(jù)并劃分warm-start段和正式評(píng)估段 n_warm 50 n_eval 250 # 2. 用warm-start段做超參數(shù)初始化估計(jì) init_X, init_y X_stream[:n_warm], y_stream[:n_warm] # ... 按4.2節(jié)代碼計(jì)算ell_init等 # 3. 初始化模型 model GPEnKF(n_ens50, n_inducing10, sigma_nsigma_n_init, ell_initell_init, sigma_f_initsigma_f_init) # 4. 正式在線學(xué)習(xí), 全程記錄預(yù)測(cè)誤差和置信區(qū)間覆蓋率 test_points np.linspace(-5, 5, 20) true_test np.sin(test_points) mean_pred np.zeros(20) std_pred np.zeros(20) for t in range(n_warm, n_warm n_eval): model.predict_step() model.update_step(X_stream[t], y_stream[t]) # 每20步做一次全測(cè)試點(diǎn)預(yù)測(cè) if (t - n_warm) % 20 0: for j, xq in enumerate(test_points): mean_pred[j], var_pred model.predict(np.array([xq])) std_pred[j] np.sqrt(var_pred) # 5. 計(jì)算RMSE和95%區(qū)間覆蓋率 rmse np.sqrt(np.mean((mean_pred - true_test)**2)) coverage np.mean((true_test mean_pred - 1.96*std_pred) (true_test mean_pred 1.96*std_pred)) print(fRMSE: {rmse:.4f}, 95% interval coverage: {coverage:.2%})RMSE 衡量預(yù)測(cè)均值的精度覆蓋率衡量不確定性量化的質(zhì)量。一個(gè)健康的實(shí)現(xiàn)在覆蓋率上應(yīng)該落在 90%–98% 之間。如果覆蓋率低于 85%說明預(yù)測(cè)方差系統(tǒng)性偏小優(yōu)先檢查 sigma_n 是否設(shè)置過小以及 inflation 是否被關(guān)掉了。覆蓋率超過 99% 則說明方差偏大模型太保守適合處理高噪聲場(chǎng)景但預(yù)測(cè)均值精度可能受損。這里有一個(gè)評(píng)估上的常見誤區(qū)覆蓋率不能用訓(xùn)練數(shù)據(jù)算必須在模型從未見過的測(cè)試點(diǎn)上算。在線場(chǎng)景下測(cè)試點(diǎn)必須在數(shù)據(jù)流開始前就劃定不能在跑完后再挑表現(xiàn)好的點(diǎn)來算——那等于拿著答案找答案。5. GP-EnKF避坑指南五個(gè)常規(guī)翻車點(diǎn)與排查手段5.1 協(xié)方差奇異與Cholesky失敗現(xiàn)象運(yùn)行過程中突然報(bào)LinAlgError: Matrix is not positive definite或者numpy.linalg.solve拋奇異矩陣錯(cuò)誤程序直接中斷。多半發(fā)生在更新步計(jì)算Kzz時(shí)。原因兩個(gè)或更多歸納點(diǎn)位置距離過近。RBF 核矩陣的列會(huì)因距離近而近乎線性相關(guān)加上浮點(diǎn)精度限制矩陣條件數(shù)爆炸。常見于更新步跑了幾百輪之后歸納點(diǎn)在 EnKF 的驅(qū)動(dòng)下慢慢擠到一起或者初始 Z 設(shè)置得過密。解決兩條防線。第一在 Kzz 對(duì)角線上加 jitter代碼里已經(jīng)寫的是 1e-6 * np.eye(m)如果問題復(fù)現(xiàn)就把 jitter 提到1e-5或1e-4。第二在每次 update_step 結(jié)束后強(qiáng)制檢查歸納點(diǎn)間距小于閾值就重新均勻散布。# 歸納點(diǎn)間距檢查與修復(fù) min_gap 1e-3 for i in range(self.n_ens): Z_i np.sort(self.Z[i]) # 排序, 保證單調(diào)性 gaps np.diff(Z_i) if gaps.min() min_gap: # 重新在[min, max]范圍內(nèi)均勻散布 self.Z[i] np.linspace(Z_i.min(), Z_i.max(), self.m)追 min_gap 閾值時(shí)先從小往大加不要一上來設(shè) 1e-2否則會(huì)頻繁觸發(fā)修復(fù)影響歸納點(diǎn)自由度。這個(gè)修復(fù)邏輯要在每次 update_step 之后、下一次 predict_step 之前執(zhí)行。5.2 ensemble塌縮與方差過小現(xiàn)象跑了一段時(shí)間后ensemble 成員幾乎完全一致self.u各行相差極小預(yù)測(cè)方差趨近于 0置信區(qū)間窄成一條線。數(shù)值上np.std(self.u, axis0)的最大值小于 1e-4。原因EnKF 更新步本質(zhì)上是線性收縮——所有成員都朝觀測(cè)靠攏方差系統(tǒng)性減小。如果過程噪聲設(shè)得極小比如 0.001 以下預(yù)測(cè)步注入的不確定性遠(yuǎn)小于更新步的收縮量幾十步后 ensemble 就塌縮成一個(gè)點(diǎn)。這是 EnKF 的已知問題不是代碼 bug。解決分三步排查。確認(rèn)process_noise至少為 0.01不要低于這個(gè)值。打開協(xié)方差膨脹把inflation從 1.0 提到 1.05–1.1。如果還不行檢查觀測(cè)噪聲sigma_n是否設(shè)置過小——觀測(cè)噪聲越小更新步的收縮越猛烈對(duì)膨脹的需求越大。# 每次更新后監(jiān)控ensemble離散度, 快速發(fā)現(xiàn)塌縮 spread np.mean(np.std(model.u, axis0)) if spread 1e-4: print(Warning: ensemble collapsed, spread , spread)記住一個(gè)判斷準(zhǔn)則預(yù)測(cè)標(biāo)準(zhǔn)差應(yīng)該和預(yù)測(cè)殘差在同一個(gè)量級(jí)。如果標(biāo)準(zhǔn)差比殘差小一個(gè)量級(jí)塌縮已經(jīng)發(fā)生立即調(diào)大 inflation 或 process_noise。5.3 非平穩(wěn)數(shù)據(jù)滯后與預(yù)測(cè)偏差現(xiàn)象數(shù)據(jù)分布中段漂移比如函數(shù)形狀從低頻變成高頻之后預(yù)測(cè)均值跟不上去殘差系統(tǒng)性增大RMSE 逐步惡化。收斂但滯后滯后長(zhǎng)度跟漂移幅度成正比。原因過程噪聲是隨機(jī)游走模型它假設(shè)狀態(tài)在單位步長(zhǎng)內(nèi)的變化幅度有限。如果漂移速度遠(yuǎn)超process_noiseEnKF 的增益系數(shù)會(huì)低估真實(shí)變化預(yù)測(cè)步注入的不確定性不足以覆蓋漂移。這和溫度計(jì)測(cè)體溫一樣——溫度計(jì)的熱慣性太大體溫已經(jīng)升高它還在慢慢爬。解決把process_noise從 0.01 提到 0.05 再觀察。如果滯后明顯改善但預(yù)測(cè)方差同步增大說明之前的過程噪聲確實(shí)太小。另一個(gè)辦法是引入遺忘因子只讓最近一段時(shí)間的觀測(cè)參與狀態(tài)估計(jì)實(shí)施方式是每 N 步把歸納點(diǎn)的后驗(yàn)方差初始化為當(dāng)前方差的 2–3 倍模擬重新開始。# 每N步增強(qiáng)一次過程噪聲, 應(yīng)對(duì)突發(fā)漂移 if t % 100 0: model.process_noise * 1.5 model.process_noise min(model.process_noise, 0.05)這個(gè)策略對(duì)突發(fā)式漂移有效但不要濫用——持續(xù)放大過程噪聲會(huì)讓穩(wěn)態(tài)預(yù)測(cè)方差虛高正常時(shí)期的表現(xiàn)會(huì)變差。5.4 歸納點(diǎn)退化與覆蓋不足現(xiàn)象數(shù)據(jù)分布在 [?5, 5]但歸納點(diǎn)慢慢集中到 [?2, 2] 的區(qū)間內(nèi)測(cè)試點(diǎn) x4 處的預(yù)測(cè)方差比其他位置大好幾倍。歸納點(diǎn)位置在更新步的牽拉下喪失了全局覆蓋。原因EnKF 更新步對(duì)歸納點(diǎn)的修改是數(shù)據(jù)驅(qū)動(dòng)的——靠近觀測(cè)位置的歸納點(diǎn)其 Kxz 值大受到的影響強(qiáng)遠(yuǎn)處歸納點(diǎn)的 Kxz 值指數(shù)級(jí)衰減幾乎不參與更新。長(zhǎng)此以往遠(yuǎn)處歸納點(diǎn)失去數(shù)據(jù)支撐逐漸漂移或被噪聲主導(dǎo)實(shí)際有效覆蓋收縮。解決定期檢查歸納點(diǎn)的覆蓋范圍覆蓋不足就強(qiáng)制重新散布。常見做法是每 50 步把歸納點(diǎn)按當(dāng)前數(shù)據(jù)分布重排一次# 每50步重新分配歸納點(diǎn)位置, 保持覆蓋 if t % 50 0: data_min, data_max X_stream[t-50:t].min(), X_stream[t-50:t].max() # 重新在最近50步的數(shù)據(jù)范圍內(nèi)生成歸納點(diǎn) model.Z np.random.uniform(data_min, data_max, size(model.n_ens, model.m))注意重新散布?xì)w納點(diǎn)時(shí)u 值不能直接丟棄——應(yīng)該用原本的 u 在舊 Z 上的后驗(yàn)對(duì)新的 Z 做插值。簡(jiǎn)化做法是用 GP 預(yù)測(cè)式重新計(jì)算u_new K(Z_new, Z_old) K(Z_old, Z_old)?1 u。完整做一次插值計(jì)算量不大但能避免歸納點(diǎn)重排引起的預(yù)測(cè)跳變。5.5 超參數(shù)發(fā)散與核寬度失衡現(xiàn)象self.log_ell在運(yùn)行穩(wěn)定期持續(xù)向一個(gè)方向漂移最終長(zhǎng)度尺度變成 0.01 或 100 這種極端值。長(zhǎng)度尺度接近 0 時(shí)核函數(shù)幾乎無平滑能力預(yù)測(cè)跟著噪聲走接近 100 時(shí)核函數(shù)完全平滑預(yù)測(cè)退化成直線。原因EnKF 對(duì) log 超參數(shù)的更新依賴 PH_T 里超參數(shù)與預(yù)測(cè)觀測(cè)的協(xié)方差項(xiàng)。如果觀測(cè)對(duì)超參數(shù)不敏感數(shù)據(jù)量太少或歸納點(diǎn)位置不佳這個(gè)協(xié)方差估計(jì)噪聲很大導(dǎo)致超參數(shù)被觀測(cè)噪聲牽著做隨機(jī)游走長(zhǎng)時(shí)間無約束漂移。解決給超參數(shù)加軟約束的偏好項(xiàng)——在預(yù)測(cè)步里把 log 超參數(shù)往初始值方向拉回一點(diǎn)幅度與偏離距離成正比# 帶約束的超參數(shù)演化 log_ell_target np.log(self.ell_init) # 初始值作為目標(biāo) self.log_ell 0.005 * np.random.randn(self.n_ens) self.log_ell - 0.002 * (self.log_ell - log_ell_target) # 回歸力回歸力系數(shù) 0.002 的含義是偏離初始值 1 個(gè) log 單位每步會(huì)被拉回 0.2%。這個(gè)強(qiáng)度足夠防止長(zhǎng)時(shí)間漂移又不會(huì)壓制真實(shí)變化。如果你用的是論文原版代碼確認(rèn)它是否包含這個(gè)約束——多數(shù)復(fù)現(xiàn)版本沒有需要自己加。6. 驗(yàn)證你的實(shí)現(xiàn)合成數(shù)據(jù)基準(zhǔn)測(cè)試與不確定性檢查拿到代碼先別急著上真實(shí)數(shù)據(jù)用已知真值的合成數(shù)據(jù)跑一遍驗(yàn)證。我的做法是從一個(gè)帶真實(shí)超參數(shù)的 GP 里采樣一條時(shí)間序列然后對(duì)比 GP-EnKF 和標(biāo)準(zhǔn)批量 GP 的預(yù)測(cè)結(jié)果。批量 GP 是標(biāo)準(zhǔn)答案如果兩者差異在 20% 以內(nèi)基本可以確認(rèn)實(shí)現(xiàn)正確。from scipy.linalg import cholesky # 從真值GP采樣: 長(zhǎng)度尺度1.0, 信號(hào)方差1.0 X_all np.linspace(-5, 5, 200) K_true rbf_1d(X_all, X_all, 1.0, 1.0) 1e-6 * np.eye(200) L cholesky(K_true, lowerTrue) y_all L np.random.randn(200) 0.05 * np.random.randn(200) # GP-EnKF在線學(xué)習(xí) model GPEnKF(n_ens50, n_inducing10, sigma_n0.05, ell_init1.0, sigma_f_init1.0) pred_mean np.zeros(200) pred_std np.zeros(200) for t in range(200): model.predict_step() model.update_step(X_all[t], y_all[t]) pred_mean[t], pred_var model.predict(np.array([X_all[t]])) pred_std[t] np.sqrt(pred_var) # 指標(biāo)1: 預(yù)測(cè)RMSE rmse np.sqrt(np.mean((pred_mean - y_all)**2)) # 指標(biāo)2: 95%區(qū)間覆蓋率 coverage np.mean((y_all pred_mean - 1.96*pred_std) (y_all pred_mean 1.96*pred_std)) print(fRMSE: {rmse:.4f}, Coverage: {coverage:.2%})兩個(gè)指標(biāo)各有側(cè)重。RMSE 衡量跟蹤精度覆蓋率衡量不確定性校準(zhǔn)度。單個(gè)指標(biāo)過關(guān)不算數(shù)必須兩個(gè)同時(shí)達(dá)標(biāo)——RMSE 很低的實(shí)現(xiàn)在覆蓋率上可能只有 50%說明模型過度自信預(yù)測(cè)方差嚴(yán)重偏小覆蓋率接近 100% 但 RMSE 偏高說明方差虛高模型太保守。還有一個(gè)更敏感的檢查項(xiàng)逐點(diǎn)殘差的標(biāo)準(zhǔn)差應(yīng)該接近預(yù)測(cè)標(biāo)準(zhǔn)差的中位數(shù)。如果殘差標(biāo)準(zhǔn)差是預(yù)測(cè)標(biāo)準(zhǔn)差的兩倍以上說明不確定性被系統(tǒng)性低估反之則被高估。這個(gè)比值是判斷 EnKF 參數(shù)是否匹配數(shù)據(jù)的快速方法。從那以后我每次在新數(shù)據(jù)集上跑 GP-EnKF都會(huì)先做一遍這個(gè)合成驗(yàn)證然后再看真實(shí)數(shù)據(jù)?；ㄊ昼娕芡昊鶞?zhǔn)能省掉后面一整天排查參數(shù)的時(shí)間。這套驗(yàn)證流程也建議你在下載代碼包后第一時(shí)間跑一遍——確認(rèn)實(shí)現(xiàn)沒有問題再上自己的數(shù)據(jù)比直接沖進(jìn)去調(diào)參靠譜得多。希望幫到你。本文還有配套的精品資源點(diǎn)擊獲取