隱馬爾可夫模型hidden Markov modelHMM這個名字聽起來像是概率論教材里最不好惹的那一章可如果你寫過語音識別、做過中文分詞、跑過基因序列標(biāo)注回頭一看會發(fā)現(xiàn)早就跟它打過照面了。它的核心設(shè)定樸素得有點反直覺系統(tǒng)在背后按一條你根本看不見的狀態(tài)鏈悄悄演化每個時刻只吐給你一個觀測值你要靠這串觀測值去反推狀態(tài)是什么、參數(shù)該怎么估。這種看得見結(jié)果、看不見原因的結(jié)構(gòu)在現(xiàn)實里比比皆是——用戶點沒點廣告你能看到他心里想買什么你看不到設(shè)備報出來的溫度你能讀到內(nèi)部是正常運轉(zhuǎn)還是已經(jīng)開始劣化你不知道。我寫這篇的動機(jī)很簡單。市面上講 HMM 的材料要么停在公式推導(dǎo)層面三個問題列出來、公式貼上去就結(jié)束了讀者抄完公式不知道代碼該長什么樣要么直接甩一句調(diào) hmmlearn 就行等到真出問題——概率全部下溢成 0、模型收斂到一個明顯不合理的結(jié)果、狀態(tài)數(shù)怎么選——完全無從下手。我打算把這兩頭接上從為什么這個算法必須這么設(shè)計講到我的實際代碼長什么樣中間那些只有踩過坑才知道的細(xì)節(jié)一個都不省。適合誰看如果你已經(jīng)能看懂基本的條件概率和矩陣乘法想真正把 HMM 用起來這篇是給你寫的。如果你只是想快速回顧一下公式也可以只看第 1 節(jié)和第 7 節(jié)一個幫你建立直覺一個幫你判斷這東西今天還值不值得用。1. 從看不見的狀態(tài)說起HMM 要解決的到底是什么問題1.1 兩個骰子的賭局先建立對隱藏狀態(tài)的物理直覺拋開公式先想象一個場景。桌上有兩個骰子一個是正常骰子六個面概率均等另一個被做過手腳擲出 6 的概率是 0.5其余五個面平分剩下的 0.5?，F(xiàn)在有人每次隨機(jī)挑一個骰子擲一次記下點數(shù)然后換不換骰子你也不知道——但你只能看到一長串點數(shù)6 6 1 6 6 6 2 6 ...。問題來了這串點數(shù)背后每一輪到底用的是哪個骰子你腦子里會自然冒出兩個判斷。第一連續(xù)看到好幾個 6那大概率是那個動了手腳的骰子。第二骰子的切換本身也有慣性一個人不太可能每一輪都換著擲上一輪用了哪個這一輪繼續(xù)用的概率應(yīng)該更高。這兩個判斷恰好就是 HMM 的兩根支柱發(fā)射概率某個狀態(tài)下觀測值出現(xiàn)的可能性和轉(zhuǎn)移概率狀態(tài)與狀態(tài)之間的切換傾向。再換個更貼近應(yīng)用的例子。假設(shè)你有個朋友住在很遠(yuǎn)的地方每天只能在社交平臺上看到他發(fā)了什么——散步購物打掃房間。你想知道他那邊天氣怎么樣但天氣預(yù)報你拿不到。這里晴/雨是隱藏狀態(tài)散步/購物/打掃是你能看到的觀測。晴天的日子里他更可能去散步雨天的日子里他更可能待在家打掃。你手上有一串他的活動記錄想反推這幾天的天氣序列。這就是 HMM 最典型的用法觀測序列已知狀態(tài)序列未知。我說這兩個例子的意思是想強(qiáng)調(diào)一點HMM 不是先有數(shù)學(xué)再有應(yīng)用而是先有隱藏原因驅(qū)動可見結(jié)果這類問題然后才被人用概率語言寫下來的。很多人學(xué) HMM 卡住就是卡在直接看公式腦子里沒有一個具體場景在跑。你先把骰子或者天氣這個畫面刻進(jìn)腦子后面所有的 α、β、γ 都只是給這個畫面配上的符號而已。1.2 五個要素各自管什么π、A、B 的分工不能混一個標(biāo)準(zhǔn)的 HMM 由五樣?xùn)|西確定通常記成一個元組 λ (N, M, π, A, B)。N 是隱藏狀態(tài)的個數(shù)M 是觀測值的個數(shù)離散情形下就是詞典大小。剩下三個才是真正需要估的參數(shù)也是新手最容易混的地方。參數(shù)名稱形狀物理含義約束π初始狀態(tài)分布長度 N 的向量序列第一個時刻處于各狀態(tài)的概率非負(fù)和為 1A狀態(tài)轉(zhuǎn)移矩陣N × N 矩陣從狀態(tài) i 跳到狀態(tài) j 的概率 a_ij每行非負(fù)行和為 1B發(fā)射概率矩陣N × M 矩陣狀態(tài) i 下觀測到符號 k 的概率 b_i(k)每行非負(fù)行和為 1這張表我建議你對著代碼多看幾遍因為九成的實現(xiàn) bug 都出在這三個矩陣的行列方向上。A 是行是當(dāng)前狀態(tài)、列是下一狀態(tài)也就是說A[i][j]表示 i → jB 是行是狀態(tài)、列是觀測符號B[i][k]表示狀態(tài) i 生成符號 k。一旦你把 A 寫成了轉(zhuǎn)置前向算法照樣能跑出數(shù)字但那個數(shù)字是錯的而且不會報錯這才是最要命的。另外注意行和為 1 這個約束。π 是一個概率分布A 的每一行、B 的每一行都是一個概率分布。你在做參數(shù)更新比如 Baum-Welch或者手工調(diào)參的時候如果破壞了行歸一化整個模型就失去了概率解釋算出來的似然再也無法比較大小。我在代碼里習(xí)慣在每次參數(shù)更新后強(qiáng)制加一步歸一化寧可多一次除法也不留下臟參數(shù)。狀態(tài)數(shù) N 的選擇是個獨立話題第 5 節(jié)會專門講。這里先記住一個原則N 是你對系統(tǒng)內(nèi)部有幾種模式的先驗假設(shè)它不是一個能靠數(shù)據(jù)自動確定的量只能靠實驗去試。1.3 三大問題與算法對應(yīng)別把它們記成三個孤立知識點教科書里 HMM 一定會列出三個基本問題很多人把它們當(dāng)作三個獨立的考點背下來其實它們是同一條流水線上的三道工序。問題一句話描述經(jīng)典算法時間復(fù)雜度評估問題給定 λ 和觀測序列 O算 P(O | λ)前向算法或前向-后向O(N2T)解碼問題給定 λ 和 O找最可能的狀態(tài)序列 Q*Viterbi 算法O(N2T)學(xué)習(xí)問題只給 O反推最合適的 λBaum-WelchEM 的一種特例O(N2T) 每輪它們的關(guān)系是這樣的學(xué)習(xí)問題內(nèi)部要反復(fù)調(diào)用評估問題來算似然和新舊參數(shù)而解碼問題在訓(xùn)練完之后才用得上。所以真正動手實現(xiàn)的時候順序一定是先寫前向、再寫 Viterbi、最后寫 Baum-Welch。反過來寫會非常痛苦因為 Baum-Welch 的正確性完全依賴前向和后向的中間結(jié)果而這兩個東西寫錯了訓(xùn)練出來的參數(shù)雖然能跑但毫無意義。還有個容易被忽略的點這三個問題的復(fù)雜度都是 O(N2T)而不是 O(N^T)。這個差距是決定性的也是第 2 節(jié)要展開的核心。簡單說如果復(fù)雜度是 N 的 T 次方那 HMM 根本不可能被用在真實數(shù)據(jù)上一條一百幀的語音就足以讓任何計算機(jī)算到天荒地老。2. 前向算法為什么不能暴力枚舉一次復(fù)雜度上的生死抉擇2.1 先算算暴力枚舉到底有多離譜要算 P(O | λ)最直覺的做法是把所有可能的狀態(tài)序列都列出來對每一條序列算這條序列產(chǎn)生這串觀測的概率最后加總。公式寫出來很干凈P(O | λ) Σ_Q π_{q1} b_{q1}(o_1) a_{q1q2} b_{q2}(o_2) ... a_{q_{T-1}q_T} b_{qT}(o_T)問題出在那個 Σ_Q 上。Q 是長度為 T 的狀態(tài)序列每個位置有 N 種取值所以一共有 N^T 條序列。每條序列要算 2T 次乘法?？偟挠嬎懔渴?O(T · N^T)。我們代入具體數(shù)字感受一下。假設(shè) N 3一個很小的模型T 100一百幀語音這在語音識別里算極短的。3 的 100 次方大約是 5×10^47。就算你的機(jī)器每秒能算十億次乘法也需要 10^31 秒。宇宙年齡才 4×10^17 秒。這個差距不是靠優(yōu)化常數(shù)或者換硬件能彌補(bǔ)的它是量級上的不可能。提示這個指數(shù)爆炸是所有序列模型都要面對的第一道墻。前向算法、Viterbi、CTC 的前向后向本質(zhì)上都是在用動態(tài)規(guī)劃把這堵墻推倒把指數(shù)復(fù)雜度壓成多項式復(fù)雜度。2.2 前向變量的遞推結(jié)構(gòu)把公共子問題緩存下來關(guān)鍵洞察是不同路徑之間共享了大量前綴。比如兩條狀態(tài)序列1→2→1→3和1→2→1→2它們前三個位置完全一樣產(chǎn)生前三個觀測的概率也一樣只是最后一個位置分岔了。暴力枚舉會把這段前綴算兩遍而動態(tài)規(guī)劃只算一遍。于是定義前向變量α_t(i) P(o_1, o_2, ..., o_t, q_t i | λ)讀法是到第 t 個時刻為止觀測序列是 o_1 到 o_t并且第 t 個時刻恰好處于狀態(tài) i 的聯(lián)合概率。注意這是聯(lián)合概率不是條件概率這一點很多資料講得不清楚導(dǎo)致后面跟 β 相乘的時候有人想不通為什么可以直接乘。遞推關(guān)系來自一個很自然的分解α_t(j) [ Σ_i α_{t-1}(i) · a_ij ] · b_j(o_t)翻譯成人話要想到達(dá)第 t 刻在狀態(tài) j 且看到了 o_t這個局面可以從上一刻的任意狀態(tài) i 出發(fā)跳過來概率是 α_{t-1}(i) · a_ij把所有 i 的可能性加起來再乘上在狀態(tài) j 下吐出觀測 o_t的概率。初始條件是 α_1(i) π_i · b_i(o_1)終止時 P(O | λ) Σ_i α_T(i)。復(fù)雜度一下就降到了 O(N2T)。還是 N 3、T 100計算量大概是 900 次乘法眨眼就完事了。從 10^47 到 900這就是動態(tài)規(guī)劃在序列建模里的分量。2.3 數(shù)值下溢浮點數(shù)會在你毫無察覺的時候把結(jié)果抹成 0上面那段推導(dǎo)數(shù)學(xué)上完全正確但直接照著寫代碼跑到 T 稍微大一點就會出事而且是那種結(jié)果直接變成 0.0的靜默錯誤。原因是你一直在做連乘。α_1(i) 是個小于 1 的數(shù)每往前推一步就要乘一堆小于 1 的概率數(shù)值單調(diào)遞減。double 類型能表示的最小正規(guī)格化數(shù)大約在 1e-308 量級一旦乘積掉到這個線以下就成了非規(guī)格化數(shù)再掉下去直接變成 0。粗略估算如果每一步平均乘一個 0.1 量級的因子那么大約 T 308 步之后α 就全線歸零了。后果不只是結(jié)果不準(zhǔn)而是整個后續(xù)計算全部失效。似然變成 0 意味著對數(shù)似然變成負(fù)無窮Baum-Welch 里 α · β 全是 0γ 做歸一化的時候是 0/0得到 NaNNaN 一旦進(jìn)入?yún)?shù)矩陣后面所有迭代都是 NaN。我見過不止一個項目在這上面浪費了一整天以為是數(shù)據(jù)有問題實際就是沒做數(shù)值處理。注意只要你的序列長度可能超過幾十就一定要處理下溢。這不是優(yōu)化項是必選項。2.4 縮放因子與 log 域兩條路線的取舍主流解法有兩條。第一條是縮放因子法scaling也是教科書上最常推的做法。思路是每一步算完 α_t(i) 之后除以當(dāng)前時刻所有狀態(tài)的和 c_tc_t 1 / Σ_i α_t(i)α?_t(i) α_t(i) · c_t因為做了歸一化α?_t(i) 的數(shù)值始終在 0 到 1 之間而且每時刻所有狀態(tài)的和恰好是 1永遠(yuǎn)不會下溢。同時這個過程可逆原始的似然滿足 P(O | λ) 1 / Π_t c_t取對數(shù)就是 log P(O | λ) -Σ_t log c_t。所以縮放不但解決了下溢還免費給了你一個數(shù)值穩(wěn)定的對數(shù)似然。第二條是全程對數(shù)域log domain。把所有概率取對數(shù)乘法變加法加法用 log-sum-exp 處理。這條路的好處是不會下溢負(fù)數(shù)域可以到 -700 以下壞處是要自己實現(xiàn) log-sum-exp而且不能再用 numpy 的矩陣乘法直接操作——log(exp(a) exp(b))每次都要手工展開代碼會變啰嗦。方案優(yōu)點缺點適用場景縮放因子法可復(fù)用矩陣乘法代碼簡潔順帶得到對數(shù)似然概率本身接近 0 時縮放因子會很大絕大多數(shù)離散 HMM 實現(xiàn)對數(shù)域數(shù)值最穩(wěn)無下溢風(fēng)險需要 log-sum-exp向量化困難序列超長、或要做混合精度直接計算代碼最少必下溢只適合 T 20 的演示我個人的習(xí)慣是縮放法為主。Viterbi 那部分則直接全程對數(shù)域因為它只涉及加法和 max不需要 log-sum-exp天然適合對數(shù)域。這兩條路線在同一份代碼里并存其實非常自然后文第 6 節(jié)會給出完整實現(xiàn)。3. Viterbi 解碼從觀測序列反推最可能的狀態(tài)軌跡3.1 動態(tài)規(guī)劃表里每一項代表什么評估問題問的是這串觀測整體出現(xiàn)的概率有多大解碼問題問的是最可能是哪條狀態(tài)鏈產(chǎn)生了它。這兩個問題的目標(biāo)不同算法結(jié)構(gòu)卻驚人地相似只是把求和換成了取最大。Viterbi 定義了一個新變量δ_t(j) max over q_1..q_{t-1} of P(q_1, ..., q_{t-1}, q_t j, o_1, ..., o_t | λ)注意這里取的是路徑上的最大值而不是所有路徑的總和。遞推式寫成δ_t(j) [ max_i δ_{t-1}(i) · a_ij ] · b_j(o_t)直覺上很好理解我要找一條最優(yōu)路徑走到第 t 刻在狀態(tài) j那么它的前綴一定是走到第 t-1 刻某個狀態(tài) i的最優(yōu)路徑。這是最優(yōu)子結(jié)構(gòu)也是動態(tài)規(guī)劃能成立的前提。這里有一個必須點破的細(xì)節(jié)Viterbi 求出來的路徑不等于每個時刻單獨取最可能狀態(tài)拼起來的那條路徑。前者是所有時刻聯(lián)合起來最可能的整條路徑后者是逐點最大化邊緣概率。這兩個結(jié)果經(jīng)常會不一樣。當(dāng)年我第一次接觸這個區(qū)別時也覺得別扭但仔細(xì)想想就通了逐點最優(yōu)拼起來的那條路徑很可能包含了一個概率極低的轉(zhuǎn)移比如從狀態(tài) A 一步跳到狀態(tài) D而這種跳轉(zhuǎn)在實際模型里幾乎不會發(fā)生。聯(lián)合最優(yōu)則會繞開它。3.2 回溯指針別忘了記下從哪來Viterbi 的實現(xiàn)必須維護(hù)兩張表一張存 δ 值一張存指針 ψ。ψ_t(j) 記錄的是在第 t 刻到達(dá)狀態(tài) j 的那條最優(yōu)路徑上一刻是從哪個狀態(tài)來的。為什么必須存指針因為動態(tài)規(guī)劃是自底向上填表的填到最后一刻你只知道終點在哪不知道路徑怎么走。沒有指針的話你得從頭重新推一遍那就退化成指數(shù)搜索了。存了指針最后從 δ_T 的最大值位置開始一步步順著 ψ 往回跳就能在 O(T) 時間里還原完整路徑。我見過有人寫 Viterbi 只在最后對 δ_T 取了個 argmax 就當(dāng)成答案輸出那其實只輸出了最后一個時刻的狀態(tài)前面全是錯的。這個 bug 特別隱蔽因為輸出長度對不上通常會被立刻發(fā)現(xiàn)但如果他恰好在最后加了個 padding 或者只關(guān)心最后一個狀態(tài)就完全看不出來了。3.3 一段可以直接跑的 Viterbi 實現(xiàn)import numpy as np def viterbi(obs, pi, A, B): obs: 觀測序列元素為 0..M-1 的整數(shù) pi : (N,) 初始狀態(tài)分布 A : (N, N) 轉(zhuǎn)移矩陣A[i, j] P(i - j) B : (N, M) 發(fā)射矩陣B[i, k] P(觀測 k | 狀態(tài) i) 返回: (最優(yōu)狀態(tài)路徑, 該路徑的對數(shù)概率) T len(obs) N len(pi) EPS 1e-300 log_pi np.log(pi EPS) log_A np.log(A EPS) log_B np.log(B EPS) delta np.full((T, N), -np.inf) psi np.zeros((T, N), dtypeint) # 初始化 delta[0] log_pi log_B[:, obs[0]] # 遞推 for t in range(1, T): for j in range(N): # 上一刻所有狀態(tài)轉(zhuǎn)移到 j 的候選值 candidates delta[t - 1] log_A[:, j] best_i int(np.argmax(candidates)) psi[t, j] best_i delta[t, j] candidates[best_i] log_B[j, obs[t]] # 回溯 path np.zeros(T, dtypeint) path[T - 1] int(np.argmax(delta[T - 1])) for t in range(T - 2, -1, -1): path[t] psi[t 1, path[t 1]] return path, delta[T - 1].max()代碼里有幾個地方值得單獨說一下。第一加EPS再取對數(shù)。如果 A 或 B 里有精確的 0np.log(0) 會返回 -inf 并拋出運行時警告。加上一個極小的正數(shù)就把 -inf 變成了一個很大的負(fù)數(shù)約 -690這在后續(xù)做加法時行為正確不會污染結(jié)果。這也是工程實踐和教科書推導(dǎo)的差別數(shù)學(xué)上 0 就是 0工程上你得給它一個地板值。第二delta初始化為 -inf 而不是 0。如果初始化成 0那些理論上不可能的狀態(tài)在第一天就會有個 0 的對數(shù)概率比真實可能狀態(tài)的負(fù)數(shù)值還大argmax 會直接選錯。用一個明確的負(fù)無窮表示不可達(dá)語義才干凈。第三內(nèi)層循環(huán)里的delta[t-1] log_A[:, j]是向量化寫法比再套一層 Python 循環(huán)快得多。別看 N 小的時候差別不明顯當(dāng) N 到幾十、T 到幾千的時候這個改寫能帶來十倍以上的速度差異。3.4 Viterbi 與前向算法只差一個運算符把兩個算法并排放在一起看你會發(fā)現(xiàn)它們的骨架幾乎完全一致。對比項前向算法Viterbi 算法中間量語義前綴聯(lián)合概率之和前綴路徑的最大概率遞推核心Σ_i α_{t-1}(i)·a_ijmax_i δ_{t-1}(i)·a_ij是否需要指針不需要需要 ψ 用于回溯數(shù)值處理縮放因子或?qū)?shù)域天然適合對數(shù)域輸出標(biāo)量 P(O|λ)狀態(tài)序列 路徑得分這個對比給我的啟發(fā)是很多看起來不同的算法骨子里是同一個動態(tài)規(guī)劃框架只是聚合算子不同——求和得到的是整體可能性取最大得到的是最優(yōu)解。同樣的思路你在 CTC 的前向后向、條件隨機(jī)場的維特比解碼里都會再遇到一次。理解了這一層再看那些算法就不會覺得是一個個孤立的黑盒了。順帶提一個工程上的細(xì)節(jié)Viterbi 的輸出得分不是概率而是對數(shù)域上的路徑得分。你沒法拿它和 P(O|λ) 直接比大小也不應(yīng)該對它做 exp 之后再解釋成這條路徑的概率。它是最優(yōu)路徑的得分僅此而已。我見過有人拿這個分?jǐn)?shù)去做閾值判斷結(jié)果閾值怎么調(diào)都不對就是因為混淆了這兩個量。4. Baum-Welch 訓(xùn)練沒標(biāo)注數(shù)據(jù)時參數(shù)是怎么學(xué)出來的4.1 EM 的套路先猜狀態(tài)再更新參數(shù)來回拉鋸前面兩節(jié)都假設(shè) λ (π, A, B) 已知?？涩F(xiàn)實里最常見的情況是你只有一堆觀測序列一個標(biāo)注都沒有。這時候要估參數(shù)靠的就是 Baum-Welch 算法。它的本質(zhì)是 EM期望最大化在 HMM 上的具體形式。EM 的思路可以粗暴地概括成一句話如果我大概知道狀態(tài)序列參數(shù)就好估如果我大概知道參數(shù)狀態(tài)序列就好猜。那就先隨便給個參數(shù)用參數(shù)猜狀態(tài)再用猜出來的狀態(tài)更新參數(shù)反復(fù)來回。具體到 HMM這個猜是通過計算兩個后驗量來完成的。γ_t(i) P(q_t i | O, λ)第 t 刻處于狀態(tài) i 的后驗概率ξ_t(i, j) P(q_t i, q_{t1} j | O, λ)第 t 刻在 i、下一刻在 j 的后驗概率這兩個量都可以從前向變量和后向變量直接算出來。后向變量 β_t(i) P(o_{t1}, ..., o_T | q_t i, λ)含義是給定第 t 刻在狀態(tài) i后面那串觀測出現(xiàn)的概率。有了它倆γ_t(i) α_t(i)·β_t(i) / P(O|λ) ξ_t(i, j) α_t(i)·a_ij·b_j(o_{t1})·β_{t1}(j) / P(O|λ)這里就能看出為什么前向和后向必須都用縮放后的版本分子分母里的縮放因子會精確地相互抵消而 P(O|λ) 用縮放因子算出來的對數(shù)似然對應(yīng)上整套推導(dǎo)是一致的。4.2 重估公式背后其實是很樸素的數(shù)數(shù)拿到 γ 和 ξ 之后參數(shù)更新公式長這樣π?_i γ_1(i)a_ij Σ_{t1}^{T-1} ξ_t(i, j) / Σ_{t1}^{T-1} γ_t(i)b?_i(k) Σ_{t: o_t k} γ_t(i) / Σ_{t1}^{T} γ_t(i)別看符號多翻譯成人話非常直觀。a_ij 的分子是整個序列里從 i 跳到 j 的期望次數(shù)分母是整個序列里待在 i 的期望次數(shù)兩者一除就是待在 i 的時候往 j 走的比例。這和統(tǒng)計學(xué)里最原始的頻率估計完全是一個路子只不過因為狀態(tài)看不見我們用的是期望次數(shù)而不是實際次數(shù)。b?_i(k) 同理分子是期望在狀態(tài) i 并且觀測到符號 k 的次數(shù)分母是期望待在狀態(tài) i 的總次數(shù)。π?_i 直接取第一時刻的狀態(tài)后驗。我在紙上推第一遍的時候總覺得用估計出來的 γ 去更新參數(shù)這件事有點循環(huán)論證的味道。后來想明白了EM 每一步都保證似然不下降這是它有收斂性保證的原因所以這個循環(huán)是朝上走的不會越走越差。代價是它可能停在局部最優(yōu)這就引出了下一小節(jié)。4.3 初始化決定你能不能跳出局部最優(yōu)Baum-Welch 是局部優(yōu)化算法初始值給得不好結(jié)果就會很差。幾種常見做法隨機(jī)初始化最省事但結(jié)果不穩(wěn)定。同一個模型跑兩次可能一次收斂到似然 -320另一次 -410。如果你在做實驗對比這會讓你懷疑人生。我的做法是固定隨機(jī)種子并且至少跑三到五次取最好的那個。均勻初始化即 π 和 A 全填 1/N、B 全填 1/M然后在上面加一點小擾動。不加擾動的話所有狀態(tài)完全對稱梯度方向可能也一樣導(dǎo)致多個狀態(tài)收斂到同一個解等于白設(shè)了那么多個狀態(tài)。這是個很容易忽略的坑你以為設(shè)了 5 個狀態(tài)訓(xùn)練完發(fā)現(xiàn) 5 行 B 幾乎一模一樣?；谙闰炛R的初始化在我做過的項目里效果最好。比如做設(shè)備狀態(tài)監(jiān)測我知道大概有正常、亞健康、異常三類那我就可以讓 π 和 A 偏向正常狀態(tài)持續(xù)時間長、異常狀態(tài)一旦進(jìn)入不太容易出來這種結(jié)構(gòu)。把領(lǐng)域知識編碼進(jìn)初始值往往比多跑幾十輪迭代更有用。還有一個技巧值得單獨提參數(shù)綁定與結(jié)構(gòu)化約束。如果你的應(yīng)用里能確信某些轉(zhuǎn)移概率應(yīng)該為 0比如狀態(tài)只能沿著 1→2→3 單向演進(jìn)那就把對應(yīng)的初始化值設(shè)成 0 并加地板值保護(hù)。這相當(dāng)于給模型加了個強(qiáng)先驗?zāi)茱@著減少要估的參數(shù)數(shù)量在小數(shù)據(jù)上尤其有用。4.4 收斂判據(jù)什么時候該停手標(biāo)準(zhǔn)的停止條件是似然增量小于閾值比如連續(xù)兩輪的對數(shù)似然變化小于 1e-4 就停。但我實際用下來發(fā)現(xiàn)光看似然不夠還得盯兩個東西。第一是迭代次數(shù)上限。EM 在接近收斂的時候會變得非常慢一輪只漲 0.0001 的似然也能磨幾十輪。設(shè)一個比如 100 輪的上限能省下大量時間。第二是過擬合。如果訓(xùn)練數(shù)據(jù)很少而狀態(tài)數(shù)很多似然會一直漲但那是在記憶訓(xùn)練數(shù)據(jù)。判斷方法很簡單留一小部分?jǐn)?shù)據(jù)做驗證集看驗證集似然什么時候開始掉。掉頭的位置才是真正該停的地方。注意對數(shù)似然單調(diào)不降是 EM 在完全數(shù)據(jù)下的理論性質(zhì)。如果你發(fā)現(xiàn)似然在下降那說明實現(xiàn)里有 bug最常見的是縮放因子沒處理好、γ 的歸一化用了錯誤的軸、或者前后向的初始化邊界寫錯了。5. 把模型落到工程里觀測表示、平滑與狀態(tài)數(shù)5.1 連續(xù)觀測離散化的代價前面所有討論都假設(shè)觀測是離散符號。但現(xiàn)實數(shù)據(jù)大多是連續(xù)的語音的 MFCC 特征、傳感器的溫度讀數(shù)、金融的收益率。最省事的做法是先做向量量化VQ用一個聚類算法把連續(xù)向量映射成若干個離散碼字然后套用離散 HMM。這么做的問題很明顯量化誤差是硬損失。兩個相差很小的向量可能被分到不同碼字相差很大的也可能因為落在聚類邊界上被分到同一個。而且碼本大小 M 又成了一個要調(diào)的超參數(shù)M 太小丟信息M 太大則 B 矩陣參數(shù)量暴漲N × M在小數(shù)據(jù)上直接過擬合。我一般的判斷標(biāo)準(zhǔn)是如果特征維度低比如 1 到 3 維、數(shù)據(jù)量中等VQ 離散 HMM 夠用且實現(xiàn)簡單如果特征維度高十幾維以上或者對精度有要求就別省這個事直接上連續(xù)觀測模型。5.2 GMM-HMM語音識別時代的經(jīng)典配置連續(xù)觀測下最常見的選擇是每個狀態(tài)配一個高斯混合模型也就是所謂的 GMM-HMM。此時發(fā)射概率不再是查表而是b_i(o) Σ_{m1}^{K} w_{im} · N(o | μ_{im}, Σ_{im})每個狀態(tài) i 有 K 個高斯分量每個分量有自己的權(quán)重、均值和協(xié)方差。參數(shù)量比離散情形大得多Baum-Welch 的更新公式也要相應(yīng)擴(kuò)展到每個觀測點對每個高斯分量各有一個后驗責(zé)任度。這套結(jié)構(gòu)在上世紀(jì)九十年代到 2010 年代初是語音識別的主流直到被深度神經(jīng)網(wǎng)絡(luò)取代。如果你想找一份完整、嚴(yán)謹(jǐn)?shù)?GMM-HMM 推導(dǎo)那部分內(nèi)容量足夠單獨寫一篇這里只提兩個我踩過的坑一是協(xié)方差矩陣的對角化假設(shè)。為了減少參數(shù)量實踐中幾乎都用對角協(xié)方差。這意味著你假設(shè)特征各維之間獨立。如果特征之間有強(qiáng)相關(guān)比如你直接塞了原始 MFCC 而沒有做 DCT 去相關(guān)這個假設(shè)會顯著傷害性能。二是方差地板。某些高斯分量的方差會在訓(xùn)練中越縮越小最后接近 0導(dǎo)致那個分量的似然在遇到稍微偏離均值的點時就變成天文數(shù)字把整個后驗壓成獨熱向量。做法是給方差設(shè)一個下限比如特征全局方差的 1% 或者一個固定的小常數(shù)。5.3 零概率與平滑地板值該給多少離散 HMM 里如果某個觀測符號在訓(xùn)練數(shù)據(jù)中從未在某個狀態(tài)下出現(xiàn)過那 b_i(k) 的估計值就是 0。這會帶來兩個后果一是取對數(shù)變 -inf二是如果測試時真出現(xiàn)了這個符號整條路徑的得分直接變成負(fù)無窮模型完全無法處理。標(biāo)準(zhǔn)解法是拉普拉斯平滑分子加 1、分母加 Mb?_i(k) (count_i(k) 1) / (Σ_k count_i(k) M)如果 M 很大比如碼本有幾千個加 1 已經(jīng)不夠了更常用的是加一個更小的偽計數(shù)比如 0.01 或者 0.1。這個數(shù)怎么選我的經(jīng)驗是先給 0.1 跑一輪看未登錄符號的處理是否正常再根據(jù)驗證集表現(xiàn)微調(diào)。對 A 矩陣其實也需要類似處理只是問題沒那么嚴(yán)重因為狀態(tài)數(shù)通常比觀測數(shù)少A 里出現(xiàn)零的概率低得多。但結(jié)構(gòu)化約束人為把某些 a_ij 設(shè)為 0的時候就要小心了這些人為的 0 不應(yīng)該被平滑掉需要在更新后重新置零。5.4 狀態(tài)數(shù)怎么選一個沒有標(biāo)準(zhǔn)答案的問題這是 HMM 實踐里最玄學(xué)的一環(huán)。狀態(tài)數(shù)不是從數(shù)據(jù)里學(xué)出來的是你自己定的。我的做法一般是這樣數(shù)據(jù)規(guī)模建議狀態(tài)數(shù)范圍理由幾百幀2 - 4參數(shù)少避免過擬合幾千幀4 - 10能建模一定復(fù)雜度幾萬幀以上10 - 30數(shù)據(jù)足以支撐參數(shù)估計有領(lǐng)域知識按物理含義定比如正常/預(yù)警/故障就是 3核心判斷方法是畫似然-狀態(tài)數(shù)曲線。狀態(tài)數(shù)從 2 加到 10看訓(xùn)練集似然和驗證集似然各怎么變。如果驗證集似然在某個點之后不再上升甚至下降那個點就是上限。如果訓(xùn)練集似然一直漲、驗證集不漲說明已經(jīng)過擬合了該選更小的狀態(tài)數(shù)。還有個更直觀的診斷方法訓(xùn)練完之后把 B 矩陣每一行打印出來或者畫成熱力圖看狀態(tài)之間有沒有明顯的差異。如果好幾個狀態(tài)的發(fā)射分布幾乎一樣說明狀態(tài)數(shù)給多了或者初始化沒擾動好那幾個狀態(tài)是冗余的。這個檢查我基本每次都會做比看數(shù)字直觀得多。6. 實戰(zhàn)從零寫一個能用的 HMM 工具類6.1 數(shù)據(jù)結(jié)構(gòu)與接口設(shè)計前面把三塊算法分別講完了這里把它們拼成一個完整的類。接口設(shè)計上我堅持兩點一是參數(shù)全部用 numpy 數(shù)組傳形狀固定為 π(N,)、A(N,N)、B(N,M)二是所有對外方法都返回對數(shù)域的數(shù)值避免調(diào)用者拿到一個已經(jīng)下溢成 0 的概率還不知道。import numpy as np EPS 1e-300 class DiscreteHMM: def __init__(self, n_states, n_obs, seed0): rng np.random.default_rng(seed) self.N n_states self.M n_obs # 均勻初始化 小擾動避免狀態(tài)對稱 self.pi np.full(n_states, 1.0 / n_states) self.A np.full((n_states, n_states), 1.0 / n_states) self.B np.full((n_states, n_obs), 1.0 / n_obs) self.A rng.normal(0, 0.01, self.A.shape) self.B rng.normal(0, 0.01, self.B.shape) self._normalize() def _normalize(self): self.pi np.clip(self.pi, EPS, None) self.pi / self.pi.sum() self.A np.clip(self.A, EPS, None) self.A / self.A.sum(axis1, keepdimsTrue) self.B np.clip(self.B, EPS, None) self.B / self.B.sum(axis1, keepdimsTrue)_normalize這一步是整套代碼的安全網(wǎng)。不管前面的更新公式寫了什么只要最后過一遍這個函數(shù)參數(shù)就一定滿足概率約束。我建議在每個對外修改參數(shù)的方法末尾都調(diào)它一次成本極低收益極高。6.2 前向、后向與 Baum-Welch 的實現(xiàn)要點def forward(self, obs): T len(obs) alpha np.zeros((T, self.N)) c np.zeros(T) alpha[0] self.pi * self.B[:, obs[0]] c[0] 1.0 / (alpha[0].sum() EPS) alpha[0] * c[0] for t in range(1, T): # alpha[t-1] 與 A 的轉(zhuǎn)置做乘法等價于 sum_i alpha[t-1,i] * A[i,j] alpha[t] (alpha[t - 1] self.A) * self.B[:, obs[t]] c[t] 1.0 / (alpha[t].sum() EPS) alpha[t] * c[t] log_prob -np.sum(np.log(c EPS)) return alpha, c, log_prob def backward(self, obs, c): T len(obs) beta np.zeros((T, self.N)) beta[T - 1] c[T - 1] for t in range(T - 2, -1, -1): beta[t] (self.A (self.B[:, obs[t 1]] * beta[t 1])) * c[t] return beta前向里的alpha[t-1] self.A是關(guān)鍵。alpha[t-1]是長度 N 的行向量self.A是 N×N矩陣乘法的結(jié)果第 j 個分量正好是 Σ_i alpha[t-1,i]·A[i,j]一步到位把內(nèi)層循環(huán)消掉了。這個改寫帶來的速度提升非常明顯值得記住。后向里的beta[T-1] c[T-1]容易寫錯。標(biāo)準(zhǔn)推導(dǎo)里如果前向做了縮放后向的初始化也必須用同樣的縮放因子才能保證 α_t · β_t 的乘積在縮放意義下自洽。很多人這里直接初始化成 1結(jié)果 γ 的歸一化之后數(shù)值看著差不多對但訓(xùn)練出來的參數(shù)總有點偏。def baum_welch_step(self, obs): T len(obs) alpha, c, log_prob self.forward(obs) beta self.backward(obs, c) gamma alpha * beta gamma / (gamma.sum(axis1, keepdimsTrue) EPS) # xi: (T-1, N, N) xi np.zeros((T - 1, self.N, self.N)) for t in range(T - 1): xi[t] (alpha[t][:, None] * self.A * self.B[:, obs[t 1]][None, :] * beta[t 1][None, :]) xi[t] / (xi[t].sum() EPS) new_pi gamma[0].copy() denom_A gamma[:-1].sum(axis0)[:, None] EPS new_A xi.sum(axis0) / denom_A denom_B gamma.sum(axis0)[:, None] EPS new_B np.zeros_like(self.B) obs_arr np.asarray(obs) for k in range(self.M): mask (obs_arr k) if mask.any(): new_B[:, k] gamma[mask].sum(axis0) new_B / denom_B return new_pi, new_A, new_B, log_probxi那段用的是廣播把 N×N 的中間結(jié)果一次算出來。如果寫成兩層循環(huán)代碼可讀性其實更高但速度慢很多。我建議先用循環(huán)版本驗證正確性確認(rèn)無誤后再換成廣播版本兩個版本的結(jié)果應(yīng)該逐元素一致這也是一種很好的自測手段。6.3 用天氣-活動這組玩具數(shù)據(jù)做端到端驗證現(xiàn)在造一組數(shù)據(jù)來跑通。隱藏狀態(tài)是天氣觀測是活動。import numpy as np rng np.random.default_rng(42) # 真實參數(shù)只用于生成數(shù)據(jù)模型不知道 pi_true np.array([0.6, 0.4]) # 晴、雨 A_true np.array([[0.7, 0.3], [0.4, 0.6]]) B_true np.array([[0.6, 0.3, 0.1], # 晴: 散步/購物/打掃 [0.1, 0.4, 0.5]]) # 雨: 散步/購物/打掃 def sample_seq(T, pi, A, B): states, obs [], [] s rng.choice(len(pi), ppi) for _ in range(T): states.append(s) obs.append(rng.choice(B.shape[1], pB[s])) s rng.choice(A.shape[0], pA[s]) return np.array(states), np.array(obs) _, obs sample_seq(3000, pi_true, A_true, B_true) model DiscreteHMM(n_states2, n_obs3, seed7) prev None for it in range(50): new_pi, new_A, new_B, logp model.baum_welch_step(obs) model.pi, model.A, model.B new_pi, new_A, new_B model._normalize() if prev is not None and abs(logp - prev) 1e-4: print(f第 {it} 輪收斂, log P {logp:.2f}) break prev logp print(學(xué)到的 A:\n, np.round(model.A, 3)) print(學(xué)到的 B:\n, np.round(model.B, 3))跑出來的 A 和 B 不會和真實值完全一樣但結(jié)構(gòu)應(yīng)該對得上一個狀態(tài)的發(fā)射分布偏向散步另一個偏向打掃。這是無監(jiān)督學(xué)習(xí)的固有特性——狀態(tài)的編號是不可辨識的學(xué)出來的狀態(tài) 0 可能對應(yīng)真實的晴天也可能對應(yīng)雨天取決于初始化。所以比較結(jié)果的時候不要直接比矩陣要看兩個狀態(tài)的發(fā)射分布是否有明確區(qū)分。如果你希望狀態(tài)編號穩(wěn)定可以在訓(xùn)練完之后按某個規(guī)則重排比如按 B 的第一列大小排序。這個操作在需要跨實驗對比的時候很有用。6.4 和現(xiàn)成庫對拍確認(rèn)自己沒寫錯從零實現(xiàn)最大的風(fēng)險是跑出來了但數(shù)字不對。我強(qiáng)烈建議用一個成熟庫做交叉驗證。Python 生態(tài)里常用的是 hmmlearn。這里有個很多人踩過的坑hmmlearn 在較新的版本里改了接口語義。老的MultinomialHMM在 0.3.0 之后表示的是多項分布輸入是計數(shù)對分類觀測要用CategoricalHMM。如果你按老教程寫MultinomialHMM會發(fā)現(xiàn)怎么都對不上甚至直接報錯。from hmmlearn.hmm import CategoricalHMM ref CategoricalHMM(n_components2, n_iter50, random_state7) ref.fit(obs.reshape(-1, 1)) print(參考實現(xiàn) log 似然:, ref.score(obs.reshape(-1, 1))) print(參考實現(xiàn)的發(fā)射矩陣:\n, np.round(np.exp(ref.emissionprob_), 3))對拍的時候注意兩點。一是似然值不能直接比大小因為隨機(jī)初始化和迭代輪數(shù)不同兩個實現(xiàn)會收斂到不同的局部最優(yōu)似然差幾個點很正常。二是要看結(jié)構(gòu)把兩個實現(xiàn)學(xué)出來的 B 矩陣按行排序后比較如果分布形態(tài)一致比如都是一個狀態(tài)偏散步、一個偏打掃說明你的實現(xiàn)是對的。真正需要警惕的是量級誤差。如果你的實現(xiàn)算出來似然是 -20000參考實現(xiàn)是 -300那基本可以確定你的縮放或者歸一化寫錯了這時候就該回去逐行檢查forward和backward里 c 的用法。7. HMM 的邊界在哪里它輸給了誰又在哪里還在贏7.1 觀測獨立假設(shè)是它的硬傷HMM 有兩個核心假設(shè)馬爾可夫性當(dāng)前狀態(tài)只依賴上一狀態(tài)和觀測獨立性當(dāng)前觀測只依賴當(dāng)前狀態(tài)。第二個假設(shè)在現(xiàn)實中經(jīng)常站不住。舉個具體的例子。做中文詞性標(biāo)注的時候一個詞的詞性顯然受到它前后若干個詞的影響而 HMM 的生成式結(jié)構(gòu)讓它在建模這種觀測之間的相互依賴時非常吃力。更麻煩的是HMM 是生成式模型它要對 P(O) 建模也就是要建模任意觀測序列出現(xiàn)的概率。對于一個有幾十萬詞的詞表這件事既沒必要也做不好。判別式模型比如條件隨機(jī)場 CRF直接對 P(Q|O) 建模繞開了對 P(O) 的建模還天然支持任意豐富的觀測特征——前后詞、詞的后綴、是否包含數(shù)字等等。這就是為什么在序列標(biāo)注任務(wù)上CRF 和后來的神經(jīng)網(wǎng)絡(luò)模型能全面壓過 HMM。7.2 與主流序列模型的定位對比模型建模對象觀測特征訓(xùn)練方式典型強(qiáng)項HMMP(O, Q) 生成式僅當(dāng)前觀測無監(jiān)督 EM 或有監(jiān)督計數(shù)數(shù)據(jù)少、缺標(biāo)注、需要可解釋狀態(tài)CRFP(Q | O) 判別式任意上下文特征有監(jiān)督序列標(biāo)注精度CTCP(Q | O) 邊緣化對齊神經(jīng)網(wǎng)絡(luò)特征有監(jiān)督無需幀級對齊端到端語音識別序列到序列模型P(Q | O)神經(jīng)網(wǎng)絡(luò)特征有監(jiān)督長距離依賴、生成任務(wù)這張表不是要評判誰高誰低而是想說清楚一件事它們解決的約束條件不同。CRF 和神經(jīng)網(wǎng)絡(luò)模型精度更高但它們都需要大量標(biāo)注數(shù)據(jù)。如果你的場景里根本拿不到標(biāo)注HMM 這種能靠 EM 從無標(biāo)注數(shù)據(jù)里學(xué)出東西的能力就是別人替代不了的。7.3 今天還值得用 HMM 的幾類場景說了這么多它的短板但我在實際工作里并沒有把它徹底淘汰掉反而有幾類場景還在持續(xù)使用。第一類是數(shù)據(jù)量小又有強(qiáng)結(jié)構(gòu)先驗的任務(wù)。比如工業(yè)設(shè)備的退化階段識別我大概知道健康狀態(tài)是緩慢變化的、故障狀態(tài)一旦進(jìn)入就不容易退出這種結(jié)構(gòu)可以直接寫進(jìn) A 矩陣的初始值里。幾百條序列就夠訓(xùn)練出一個可用的模型換成神經(jīng)網(wǎng)絡(luò)早就過擬合了。第二類是需要可解釋性的場合。HMM 學(xué)出來的狀態(tài)是有具體含義的——你打開 B 矩陣一看就知道狀態(tài) 2 是個高異常傾向的狀態(tài)。而在很多神經(jīng)網(wǎng)絡(luò)模型里你要額外做一大堆可解釋性工作還未必看得清楚。在需要給非技術(shù)同事解釋為什么系統(tǒng)判斷這個設(shè)備異常的時候HMM 的優(yōu)勢非常實在。第三類是作為更強(qiáng)模型的組件。比如用 HMM 做初步的序列切分把長序列分段之后交給下游模型處理或者用它生成偽標(biāo)簽再用這些偽標(biāo)簽去訓(xùn)練神經(jīng)網(wǎng)絡(luò)。這種用小模型撬動大模型的用法我在數(shù)據(jù)標(biāo)注資源緊張的項目里用過幾次效果比直接硬上神經(jīng)網(wǎng)絡(luò)好。第四類是快速原型驗證。想驗證一個序列建模的方向有沒有戲用 HMM 搭個基線只要幾十行代碼跑幾分鐘就出結(jié)果。對比之下搭一個神經(jīng)網(wǎng)絡(luò)基線可能要一兩天。先用最便宜的手段確認(rèn)問題本身可解再決定要不要投入更多資源這個順序我建議一直保持。最后分享一個我在實際使用中的體會HMM 的調(diào)參重點從來不在算法本身而在狀態(tài)數(shù)的選擇和初始值的設(shè)置。算法是確定的代碼寫對了就是對的但狀態(tài)數(shù)和初始值這兩個東西直接決定了你能不能用它解決問題。我見過太多人在前向后向公式上反復(fù)糾結(jié)卻從來沒打印過 B 矩陣看看到底學(xué)出了什么。多做幾次訓(xùn)完就打開參數(shù)看一眼的動作比多讀幾遍推導(dǎo)有用得多。