\n", "\n", "

\n", "Following Section 5.3.3.2 of [Müller, FMP, Springer 2015], we describe in this notebook the Viterbi algorithm, which efficiently solves the uncovering problem of hidden Markov models (HMMs). For a detailed introduction of HMMs, we refer to the famous tutorial paper by Rabiner.\n", "\n", "

- \n",
"
- \n",
"Lawrence R. Rabiner:
**A Tutorial on Hidden Markov Models and Selected Applications in Speech Recognition.**Proceedings of the IEEE, 77 (1989), pp. 257–286. \n", "

\n", " Bibtex \n", " \n",
"

\n",
"**Note:** Due to Python conventions, the indexing in the implementation starts with index

"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Observation sequence: O = [0 2 0 2 2 1]\n",
"Optimal state sequence: S = [0 0 0 2 2 1]\n",
"D =\n",
"[[ 0.4200 0.1008 0.0564 0.0135 0.0033 0.0000]\n",
" [ 0.0200 0.0000 0.0010 0.0000 0.0000 0.0006]\n",
" [ 0.0000 0.0336 0.0000 0.0045 0.0022 0.0003]]\n",
"E =\n",
"[[0 0 0 0 0]\n",
" [0 0 0 0 2]\n",
" [0 2 0 2 2]]\n"
]
}
],
"source": [
"import numpy as np\n",
"from numba import jit\n",
"\n",
"@jit(nopython=True)\n",
"def viterbi(A, C, B, O):\n",
" \"\"\"Viterbi algorithm for solving the uncovering problem\n",
"\n",
" Notebook: C5/C5S3_Viterbi.ipynb\n",
"\n",
" Args:\n",
" A (np.ndarray): State transition probability matrix of dimension I x I\n",
" C (np.ndarray): Initial state distribution of dimension I\n",
" B (np.ndarray): Output probability matrix of dimension I x K\n",
" O (np.ndarray): Observation sequence of length N\n",
"\n",
" Returns:\n",
" S_opt (np.ndarray): Optimal state sequence of length N\n",
" D (np.ndarray): Accumulated probability matrix\n",
" E (np.ndarray): Backtracking matrix\n",
" \"\"\"\n",
" I = A.shape[0] # Number of states\n",
" N = len(O) # Length of observation sequence\n",
"\n",
" # Initialize D and E matrices\n",
" D = np.zeros((I, N))\n",
" E = np.zeros((I, N-1)).astype(np.int32)\n",
" D[:, 0] = np.multiply(C, B[:, O[0]])\n",
"\n",
" # Compute D and E in a nested loop\n",
" for n in range(1, N):\n",
" for i in range(I):\n",
" temp_product = np.multiply(A[:, i], D[:, n-1])\n",
" D[i, n] = np.max(temp_product) * B[i, O[n]]\n",
" E[i, n-1] = np.argmax(temp_product)\n",
"\n",
" # Backtracking\n",
" S_opt = np.zeros(N).astype(np.int32)\n",
" S_opt[-1] = np.argmax(D[:, -1])\n",
" for n in range(N-2, -1, -1):\n",
" S_opt[n] = E[int(S_opt[n+1]), n]\n",
"\n",
" return S_opt, D, E\n",
"\n",
"# Define model parameters\n",
"A = np.array([[0.8, 0.1, 0.1], \n",
" [0.2, 0.7, 0.1], \n",
" [0.1, 0.3, 0.6]])\n",
"\n",
"C = np.array([0.6, 0.2, 0.2])\n",
"\n",
"B = np.array([[0.7, 0.0, 0.3], \n",
" [0.1, 0.9, 0.0], \n",
" [0.0, 0.2, 0.8]])\n",
"\n",
"\n",
"O = np.array([0, 2, 0, 2, 2, 1]).astype(np.int32)\n",
"#O = np.array([1]).astype(np.int32)\n",
"#O = np.array([1, 2, 0, 2, 2, 1]).astype(np.int32)\n",
"\n",
"# Apply Viterbi algorithm\n",
"S_opt, D, E = viterbi(A, C, B, O)\n",
"#\n",
"print('Observation sequence: O = ', O)\n",
"print('Optimal state sequence: S = ', S_opt)\n",
"np.set_printoptions(formatter={'float': \"{: 7.4f}\".format})\n",
"print('D =', D, sep='\\n')\n",
"np.set_printoptions(formatter={'float': \"{: 7.0f}\".format})\n",
"print('E =', E, sep='\\n')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Log-Domain Implementation of Viterbi Algorithm\n",
"\n",
"In each iteration of the Viterbi algorithm, the accumulated probability values of $\\mathbf{D}$ at step $n-1$ are multiplied with two probability values from $A$ and $B$. More precisely, we have\n",
"\n",
"$$\n",
" \\mathbf{D}(i,n) = \\max_{j\\in[1:I]}\\big(a_{ij} \\cdot\\mathbf{D}(j,n-1) \\big) \\cdot b_{i_{k_n}}.\n",
"$$\n",
"Since all probability values lie the interval $[0,1]$, the product of such values decreases exponentially with the number $n$ of iterations. As a result, for input sequences $O=(o_{1},o_{2},\\ldots o_{N})$ with large $N$, the values in $\\mathbf{D}$ typically become extremely small, which may finally lead to a numerical underflow. A well-known trick when dealing with products of probability values is to work in the **log-domain**. To this end, one applies a logarithm to all probability values and **replaces multiplication by summation**. Since the logarithm is a strictly monotonous function, ordering relations are preserved in the log-domain, which allows for transferring operations such as **maximization** or **minimization**.\n",
"\n",
"In the following code cell, we provide a **log variant** of the Viterbi algorithm. This variant yields exactly the same optimal state sequence $S^\\ast$ and the same backtracking matrix $E$ as the original algorithm, as well as the accumulated log probability matrix $\\log(\\mathbf{D})$ (`D_log`) (where the logarithm is applied for each entry). We again test the implementation on our toy example. Furthermore, as a sanity check, we apply the exponential function to the computed log-matrix `D_log`, which should yield the matrix `D` as computed by the original Viterbi algorithm above. "
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Observation sequence: O = [0 2 0 2 2 1]\n",
"Optimal state sequence: S = [0 0 0 2 2 1]\n",
"D_log =\n",
"[[ -0.87 -2.29 -2.87 -4.30 -5.73 -714.35]\n",
" [ -3.91 -711.57 -6.90 -713.57 -715.00 -7.44]\n",
" [-710.01 -3.39 -712.30 -5.40 -6.13 -8.25]]\n",
"exp(D_log) =\n",
"[[ 0.4200 0.1008 0.0564 0.0135 0.0033 0.0000]\n",
" [ 0.0200 0.0000 0.0010 0.0000 0.0000 0.0006]\n",
" [ 0.0000 0.0336 0.0000 0.0045 0.0022 0.0003]]\n",
"E =\n",
"[[0 0 0 0 0]\n",
" [0 0 0 0 2]\n",
" [0 2 0 2 2]]\n"
]
}
],
"source": [
"@jit(nopython=True)\n",
"def viterbi_log(A, C, B, O):\n",
" \"\"\"Viterbi algorithm (log variant) for solving the uncovering problem\n",
"\n",
" Notebook: C5/C5S3_Viterbi.ipynb\n",
"\n",
" Args:\n",
" A (np.ndarray): State transition probability matrix of dimension I x I\n",
" C (np.ndarray): Initial state distribution of dimension I\n",
" B (np.ndarray): Output probability matrix of dimension I x K\n",
" O (np.ndarray): Observation sequence of length N\n",
"\n",
" Returns:\n",
" S_opt (np.ndarray): Optimal state sequence of length N\n",
" D_log (np.ndarray): Accumulated log probability matrix\n",
" E (np.ndarray): Backtracking matrix\n",
" \"\"\"\n",
" I = A.shape[0] # Number of states\n",
" N = len(O) # Length of observation sequence\n",
" tiny = np.finfo(0.).tiny\n",
" A_log = np.log(A + tiny)\n",
" C_log = np.log(C + tiny)\n",
" B_log = np.log(B + tiny)\n",
"\n",
" # Initialize D and E matrices\n",
" D_log = np.zeros((I, N))\n",
" E = np.zeros((I, N-1)).astype(np.int32)\n",
" D_log[:, 0] = C_log + B_log[:, O[0]]\n",
"\n",
" # Compute D and E in a nested loop\n",
" for n in range(1, N):\n",
" for i in range(I):\n",
" temp_sum = A_log[:, i] + D_log[:, n-1]\n",
" D_log[i, n] = np.max(temp_sum) + B_log[i, O[n]]\n",
" E[i, n-1] = np.argmax(temp_sum)\n",
"\n",
" # Backtracking\n",
" S_opt = np.zeros(N).astype(np.int32)\n",
" S_opt[-1] = np.argmax(D_log[:, -1])\n",
" for n in range(N-2, -1, -1):\n",
" S_opt[n] = E[int(S_opt[n+1]), n]\n",
"\n",
" return S_opt, D_log, E\n",
"\n",
"# Apply Viterbi algorithm (log variant)\n",
"S_opt, D_log, E = viterbi_log(A, C, B, O)\n",
"\n",
"print('Observation sequence: O = ', O)\n",
"print('Optimal state sequence: S = ', S_opt)\n",
"np.set_printoptions(formatter={'float': \"{: 7.2f}\".format})\n",
"print('D_log =', D_log, sep='\\n')\n",
"np.set_printoptions(formatter={'float': \"{: 7.4f}\".format})\n",
"print('exp(D_log) =', np.exp(D_log), sep='\\n')\n",
"np.set_printoptions(formatter={'float': \"{: 7.0f}\".format})\n",
"print('E =', E, sep='\\n')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
" "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
"\n",
"`0`

. Also, there is an index shift when applying the algorithm to our toy example, which becomes `O = [0, 2, 0, 2, 2, 1]`

.\n",
"