維特比演算法和隱馬爾可夫模型的解碼

一、概述

  維特比演算法是安德魯.維特比(Andrew Viterbi)於1967年為解決通訊領域中的解碼問題而提出的，它同樣廣泛用於解決自然語言處理中的解碼問題，隱馬爾可夫模型的解碼是其中典型的代表。無論是通訊中的解碼問題還是自然語言處理中的解碼問題，本質上都是要在一個籬笆網路中尋找得到一條最優路徑。
  所謂籬笆網路，指的是單向無環圖，呈層級連線，各層節點數可以不同。如圖是一個籬笆網路，連線上的數字是節點間概念上的距離（如間距、代價、概率等），現要找到一條從起始點到終點的最優路徑。

  在實際問題中，節點數和層數往往是大量的，因而採取遍歷所有的路徑計算其距離進行比較的方式是不可行的。維特比演算法正是通過動態規劃的方式高效求得這條最優路徑。

二、維特比演算法

1.演算法原理

  該問題具有這樣一個特性，對於最優（如最短距離）的路徑，任意一段子路徑一定是該段兩端點間所有可達路徑中最優的，如若不然，將該段中更優的子路徑接到兩端點便構成了另一個整體最優路徑，這是矛盾的。或者說，最優路徑中，從起始點到由近及遠的任一點的子路徑，一定是該段所有可達路徑中最優的。也即，整體最優，區域性一定最優。

  該特性也就是說，對每個節點而言，如果最優路徑經過這一點，則一定是經過從起始點到這點的這段最優路徑。那麼，只要從頭開始，由區域性向整體推進，漸次地找到起始點到當前點的最優路徑，算至終點便得到了整體最優路徑。這樣的方式叫做動態規劃，是維特比演算法的基本思想。
維特比演算法求解籬笆網路最短路徑的過程為：
  從第一層開始，對每一層的每一個節點，計算出起始點到它的最短距離，並記錄下相應最短路徑下它的前一個節點，逐層遞推，算至終止點時便得到了整體最短距離，再依照節點記錄下的前置節點進行回溯，就得到了最短路徑的序列。對第\(i\)層第\(j\)個節點\(P_{i,j}\)，假設起始點到它的最短距離為\(D\left( P_{i,j} \right)\)，相應最短路徑下它的前一個節點為\(Pre\left( P_{i,j} \right)\)，則

\[D\left( P_{i,j} \right)=\min_{1\leq k \leq N}{\left[ D_{i-1,k}+d(P_{i-1,k},P_{i,j}) \right]} \]

也就是，對前一層的所有節點，計算每一個節點的記錄的最短距離D與它到當前節點的距離d的和，取其中最小的那個值（其中， \(d(A,B)\)表示A,B兩節點間的距離。）.

\[Pre\left( P_{i,j} \right)=P_{i-1,j\ast}=arg\min_{1\leq k \leq N}{\left[ D_{i-1,k}+d(P_{i-1,k},P_{i,j}) \right]} \]

也就是，滿足距離最短的那條路徑上在前一層的節點.

2.示例

  在如下網路中，連線上的數字是節點間的距離，求S點到E點的最短距離和與之對應的路徑。

第一層：
對節點\(P_{1,1}\)，
起始點到它只有一條路徑，最短距離\(D(P_{1,1})=2\)，前一個節點\(Pre(P_{1,1}) = S\)；

對節點\(P_{1,2}\)，
起始點到它只有一條路徑，最短距離 \(D(P_{1,2}) = 1\)，前一個節點\(Pre(P_{1,2}) = S\)；

對節點\(P_{1,3}\)，
起始點到它只有一條路徑，最短距離\(D(P_{1,3}) = 3\)，前一個節點\(Pre(P_{1,3}) = S\) 。

第二層：
對節點 \(P_{2,1}\)，
\(D(P_{1,1}) + d(P_{1,1}，P_{2,1}) = 2 + 3 = 5\)，
\(D(P_{1,2}) + d(P_{1,2}，P_{2,1}) = 1 + 2 = 3\) ，
\(D(P_{1,3}) + d(P_{1,3}，P_{2,1}) = 3 + 9 = 12\)，
最短距離 \(D(P_{2,1}) = min\left\{ 5,3,12 \right\} = 3\)，前一個節點\(Pre(P_{2,1}) = P_{1,2}\) ;

對節點\(P_{2,2}\) ，
\(D(P_{1,1}) + d(P_{1,1}，P_{2,2}) = 2 + 6 = 8\)，
\(D(P_{1,2}) + d(P_{1,2}，P_{2,2}) = 1 + 5 = 6\)，
\(D(P_{1,3}) + d(P_{1,3}，P_{2,2}) = 3 + 2 = 5\)，
最短距離\(D(P_{2,2}) = min\left\{ 8,6,5 \right\} = 5\)，前一個節點\(Pre(P_{2,2}) = P_{1,3}\) ;

對節點\(P_{2,3}\)，
\(D(P_{1,1}) + d(P_{1,1}，P_{2,3}) = 2 + 4 = 6\) ，
\(D(P_{1,2}) + d(P_{1,2}，P_{2,3}) = 1 + 7 = 8\)，
\(D(P_{1,3}) + d(P_{1,3}，P_{2,3}) = 3 + 6 = 9\) ，
最短距離\(D(P_{2,3}) = min\left\{ 6,8,9 \right\} = 6\)，前一個節點\(Pre(P_{2,3}) = P_{1,1}\) ;

第三層：
對節點 \(P_{3,1}\)，
\(D(P_{2,1}) + d(P_{2,1}，P_{3,1}) = 3 + 9 = 12\)，
\(D(P_{2,2}) + d(P_{2,2}，P_{3,1}) = 5 + 2 = 7\)，
\(D(P_{2,3}) + d(P_{2,3}，P_{3,1}) = 6 + 6 = 12\)，
最短距離\(D(P_{3,1}) = min\left\{ 12,7,12 \right\} = 7\)，前一個節點\(Pre(P_{3,1}) = P_{2,2}\);

對節點\(P_{3,2}\)，
\(D(P_{2,1}) + d(P_{2,1}，P_{3,2}) = 3 + 3 = 6\)，
\(D(P_{2,2}) + d(P_{2,2}，P_{3,2}) = 5 + 6 = 11\)，
\(D(P_{2,3}) + d(P_{2,3}，P_{3,2}) = 6 + 3 = 9\)，
最短距離\(D(P_{3,2}) = min\left\{ 6,11,9 \right\} = 6\)，前一個節點\(Pre(P_{3,2}) = P_{2,1}\);

對節點\(P_{3,3}\)，
\(D(P_{2,1}) + d(P_{2,1}，P_{3,3}) = 3 + 8 = 11\)，
\(D(P_{2,2}) + d(P_{2,2}，P_{3,3}) = 5 + 7 = 12\)，
\(D(P_{2,3}) + d(P_{2,3}，P_{3,3}) = 6 + 4 = 10\)，
最短距離\(D(P_{3,3}) = min\left\{ 11,12,10 \right\} = 10\)，前一個節點\(Pre(P_{3,3}) = P_{2,3}\);

第四層（終點）：
對節點 \(E\)，
\(D(P_{3,1}) + d(P_{3,1}，E) = 7 + 3 = 10\)，
\(D(P_{3,2}) + d(P_{3,2}，E) = 6 + 7 = 13\)，
\(D(P_{3,3}) + d(P_{3,3}，E) = 10 + 6 = 16\)，
最短距離\(D(E) = min\left\{ 10,13,16 \right\} = 10\)，前一個節點\(Pre(E) = P_{3,1}\);

又\(Pre(E) = P_{3,1}\)，\(Pre(P_{3,1}) = P_{2,2}\)，\(Pre(P_{2,2}) = P_{1,3}\)，\(Pre(P_{1,3}) = S\)
故最短距離為10，與之對應的路徑為（\(S\)，\(P_{1,3}\)，\(P_{2,2}\)，\(P_{3,1}\)，\(E\)）.

三、隱馬爾可夫模型的解碼

1.問題描述

  隱馬爾可夫模型（HMM）的解碼問題指，給定模型和輸出序列，如何找出最有可能產生這個輸出的狀態序列。自然語言處理中，也即如何通過觀測訊號確定最有可能對應的實際語義。在狀態序列上，每個狀態位是狀態集合中的元素之一，因此該問題等價於在狀態集合中的節點構成的有向網路（籬笆網路）中找出一條概率最大的路徑（最優路徑），如圖。該問題可以通過維特比演算法得到高效的解決。

2.演算法敘述

  假設 \(P(s_{t,j})\)表示從起始時刻到\(s_{t,j}\)的最優路徑的概率，\(Pre(s_{t,j})\)表示從起始時刻到 \(s_{t,j}\)的最優路徑上前一個節點，則隱馬爾可夫模型的維特比解碼演算法為：

輸入：隱馬爾可夫模型 \(\lambda=(\pi,A,B)\)和觀測 \(O=(o_1,o_2,...,o_T)\)；
輸出：最優狀態序列\(S^{\ast}=(s_{1}^{\ast},s_{2}^{\ast},...,s_{T}^{\ast})\).
(1)初始化
    \(P(s_{1,j})=\pi_{j}b_{j}(o_1)\)，
    \(Pre(s_{1,j})=None\)，\(j=1,2,...,N\)

(2)遞推
 對 \(t=2,3,...,T\)

\[P(s_{t,j})=\max_{1\leq k \leq N}{\left[ P(s_{t-1,k})a_{kj} \right]b_{j}(o_t)} \]

\[Pre(s_{t,j})=arg\max_{1\leq k \leq N}{\left[ P(s_{t-1,k})a_{kj} \right]} \]

，\(j=1,2,...,N\).

(3)遞推終止
 最大概率$$P^{\ast}=\max_{1\leq j \leq N}{P(s_{T,j})}$$
 最優路徑上的最後一個狀態$$s_{T}^{\ast}=arg\max_{1\leq j \leq N}{\left[ P(s_{T,j}) \right]}$$

(4)回溯路徑，確定最優狀態序列
    \(S^{\ast}=\left( s_{1}^{\ast},s_{2}^{\ast},...,s_{T-1}^{\ast},s_{T}^{\ast} \right)\)
     \(=\left( Pre(s_{2}^{\ast}),Pre(s_{3}^{\ast}), ...,Pre(s_{T}^{\ast}),s_{T}^{\ast}\right)\)

3.示例

（參考自《統計學習方法》）
狀態集合 \(Q=\left\{ q_1, q_2, q_3 \right\}\)，觀測集合\(V=\left\{ 0,1 \right\}\)，模型 \(\lambda=\left( \pi,A,B \right)\) ，

\(A=\begin{bmatrix} 0.5 & 0.2 & 0.3 \\ 0.3 & 0.5 & 0.2 \\ 0.2 & 0.3 & 0.5 \\ \end{bmatrix}\) , \(B=\begin{bmatrix} 0.5 & 0.5 \\ 0.4 & 0.6 \\ 0.7 & 0.3 \end{bmatrix}\),\(\pi=\left( 0.2, 0.4, 0.4 \right)^{T}\)

已知觀測序列\(O=\left( 0, 1, 0 \right)\)，求最優狀態序列。

解：
(1)在t=1時（初始化），對每一個狀態，求觀測為0的最大概率
 \(P(s_{1,1})=0.2\times0.5=0.1\)，\(Pre(s_{1,1})=None\)
 \(P(s_{1,2})=0.4\times0.4=0.16\)，\(Pre(s_{1,2})=None\)
 \(P(s_{1,3})=0.4\times0.7=0.28\)，\(Pre(s_{1,3})=None\)

(2)在t=2時，對每一個狀態，求觀測為1的
 最大概率$$P(s_{2,j})=\max_{1 \leq k \leq 3}{\left[ P(s_{1,k})a_{kj} \right]b_{j}(1)}$$
 當前最優的前一個狀態$$Pre(s_{2,j})=arg\max_{1 \leq k \leq 3}{\left[ P(s_{1,k})a_{kj} \right]}$$，\(j=1,2,3.\)
\(P(s_{2,1})=max\left\{ 0.1\times0.5\times0.5, 0.16\times0.3\times0.5, 0.28\times0.2\times0.5 \right\}\)\(=0.028\)

\(Pre(s_{2,1})=s_{1,3}=q_3\)

\(P(s_{2,2})=max\left\{ 0.1\times0.2\times0.6,0.16\times0.5\times0.6,0.28\times0.3\times0.6 \right\}\)\(=0.0504\)

\(Pre(s_{2,2})=s_{1,3}=q_3\)

\(P(s_{2,3})=max\left\{ 0.1\times0.3\times0.3, 0.16\times0.2\times0.3,0.28\times0.5\times0.3 \right\}\)\(=0.042\)

\(Pre(s_{2,3})=s_{1,3}=q_3\)

(3)在t=3時，對每一個狀態，求觀測為0的
 最大概率 $$P(s_{3,j})=\max_{1 \leq k \leq 3}{\left[ P(s_{2,k})a_{kj} \right]b_{j}(0)}$$
 當前最優的前一個狀態$$Pre(s_{3,j})=arg\max_{1 \leq k \leq 3}\left[ P(s_{2,k})a_{kj} \right]$$，\(j=1,2,3.\)

\(P(s_{3,1})=max\left\{ 0.028\times0.5\times0.5, 0.0504\times0.3\times0.5,0.042\times0.2\times0.5 \right\}\)\(=0.00756\)

\(Pre(s_{3,1})=s_{2,2}=q_2\)

\(P(s_{3,2})=max\left\{ 0.028\times0.2\times0.4, 0.0504\times0.5\times0.4, 0.042\times0.3\times0.4 \right\}\)\(=0.01008\)

\(Pre(s_{3,2})=s_{2,2}=q_2\)

\(P(s_{3,3})=max\left\{ 0.028\times0.3\times0.7,0.0504\times0.2\times0.7,0.042\times0.5\times0.7 \right\}\)\(=0.0147\)

\(Pre(s_{3,3})=s_{2,3}=q_3\)

(4)得到結果.
 最大概率
  \(P^{\ast}=\max_{1 \leq j \leq 3}P\left( s_{3,j} \right)\)
   \(=max\left\{ 0.00756,0.01008,0.0147 \right\}\)
   \(=0.0147\)
 最優狀態序列
  \(S^{\ast}=\left( Pre(s_{2}^{\ast})t,Pre(s_{3}^{\ast}),s_{3}^{\ast} \right)\)
   \(=\left( s_{1,3},s_{2,3},s_{3,3} \right)\)
   \(=\left( q_3,q_3,q_3 \right)\)

4.python實現

對上述HMM解碼示例的python實現程式為

import numpy as np

def viterbi(pi, A, B, Q, V, obs_seq):
    '''
    :param pi:HMM初始狀態概率向量，list型別
    :param A:HMM狀態轉移概率矩陣，list型別
    :param B:HMM觀測生成概率矩陣，list型別
    :param Q:狀態集合，list型別
    :param V:觀測集合，list型別
    :param obs_seq:觀測序列，list型別
    :return:最優狀態序列的概率sta_pro，float型別；最優狀態序列sta_seq，list型別
    '''

    # HMM模型引數轉換為array型別
    pi = np.array(pi)
    A = np.array(A)
    B = np.array(B)

    # 1.定義動態計算結果儲存矩陣
    rowNum = len(Q)  # 行數，狀態數
    colNum = len(obs_seq)  # 列數，生成的觀測數，即時刻數

    # 儲存節點當前最大概率的矩陣
    probaMatrix = np.zeros((rowNum,colNum))

    # 儲存當前最優路徑下的前一個節點的矩陣
    preNodeMatrix = np.zeros((rowNum,colNum))

    # 2.初始化（第1時刻）
    probaMatrix[:,0] = pi*np.transpose(B[:,obs_seq[0]])
    preNodeMatrix[:,0] = [-1]*rowNum  # 第1時刻節點的前一個節點不存在，置為-1

    # 3.遞推，第2時刻至最後
    for t in range(1, colNum):
        list_pre_max_proba = []  # 節點最大前置概率列表
        list_pre_node = []  # 節點當前最優路徑中前一個節點列表
        for j in range(rowNum):
            pre_proba_list = list(np.array(probaMatrix[:,t-1])*np.transpose(A[:,j]))  # 前置概率列表，前一時刻的節點最大概率與到當前節點轉移概率的乘積
            '''
            注：因為計算機的二進位制機制對小數的表達是有限的，所以對小數作運算將產生一定的誤差。
            在使用函式獲取pre_proba_list中的最大值和對應的索引時，為有效降低這種誤差，將資料放大後再進行操作。
            '''
            pre_proba_list = [x*pow(10,5) for x in pre_proba_list]  # 放大100000倍
            prePro = max(pre_proba_list)/pow(10,5)  # 最大前置概率
            preNodeIndexNum = pre_proba_list.index(max(pre_proba_list))  # 前置節點的索引號
            list_pre_max_proba.append(prePro)  # 最大前置概率加入列表
            list_pre_node.append(preNodeIndexNum)  # 前置節點的索引號加入列表

        probaMatrix[:,t] = np.array(list_pre_max_proba)*np.transpose(B[:,obs_seq[t]])  # 最大前置概率乘上觀測概率，即為當前最大概率
        preNodeMatrix[:,t] = list_pre_node  # 將該列前置節點索引號加入矩陣

    # 此時，得到了完整的probaMatrix和preNodeMatrix，對這兩個矩陣進行操作便可得到需要的結果
    # 4.得到最大概率
    maxPro = np.max(probaMatrix[:, colNum-1])  # 全域性最大概率（即最後一列的最大值）

    # 5.得到最優狀態序列的狀態索引號列表
    lastStateIndexNum = np.argmax(probaMatrix[:, colNum-1])  # 最優狀態序列中最後一個狀態的索引號
    stateIndexList = []  # 定義最優狀態的索引號列表
    stateIndexList.append(lastStateIndexNum)

    # 回溯，完成狀態索引號列表
    currentIndex = lastStateIndexNum;
    for t in range(colNum-1, 0, -1):
        fls = preNodeMatrix[:, t].tolist()  # 矩陣中的數值是浮點型
        ls = list(map(int, fls))  # 轉為整型
        currentIndex = ls[currentIndex]
        stateIndexList.append(currentIndex)

    stateIndexList.reverse()  # 反轉列表

    # 6.由索引號序列得到最優狀態序列
    stateSeq = [Q[i] for i in stateIndexList]

    return maxPro,stateSeq

if __name__=='__main__':
    # 狀態集合
    Q = ["q1", "q2", "q3"]

    # 觀測集合
    V = [0, 1]

    # 初始狀態概率向量
    pi = [0.2, 0.4, 0.4]

    # 狀態轉移概率矩陣
    A = [[0.5, 0.2, 0.3],
         [0.3, 0.5, 0.2],
         [0.2, 0.3, 0.5]]

    # 觀測概率矩陣
    B = [[0.5, 0.5],
         [0.4, 0.6],
         [0.7, 0.3]]

    # 觀測序列
    obs_seq = [0, 1, 0]

    maxPro, stateSeq = viterbi(pi, A, B, Q, V, obs_seq)

    print("最大概率為：", maxPro)
    print("最優狀態序列為：", stateSeq)

End.

參考

1. 吳軍. 數學之美(第二版). 人民郵電出版社.
2. 李航. 統計學習方法. 清華大學出版社.
3. https://www.cnblogs.com/zhibei/p/9391014.html