Target audience: Advanced
Estimated reading time: 7'
Many machine learning algorithms utilize dynamic programming, and the Hidden Markov Model (HMM) is a prime example of this. Scala's tail recursion optimization offers a superb alternative to the conventional iterative approach for implementing the three fundamental forms of the HMM.
What you will learn: How to write a fast implementation of Viterbi algorithm in Scala
Notes:
- Scala version: 2.12.4
- Source code available at github.com/patnicolas/ScalaforMachineLearning/supervised/hmm
- This post assumes a basic understanding of Markov Model and Scala programming language [ref 1].
Theory
Markov chain
A Markov process is a stochastic process characterized by its dependence on time, where the future's state is independent of past states if the present state is known.
The intricacies of Markov process theory vary significantly based on the nature of the time space, whether it is discrete or continuous, and the nature of the state space, which can be discrete (countable and with all subsets measurable) or a broader topological space.
The following diagram describes a simple 2-state Markov chain [ref 2]:
Fig. 1 State transition diagram for a 2-state Markov chain
The state transition probabilities matrix is defined as: \[\begin{bmatrix} 1-\alpha & \alpha \\ \beta & 1-\beta \end{bmatrix}\]
Hidden Markov model
Hidden Markov Models (HMMs) are versatile tools used across various domains to infer sequences of data that are not directly observable. Their applications include:
- Neuroscience
- Speech analysis and recognition
- Gene sequencing prediction
- Time series analysis in general
- Logistics and supply chain optimization
- Handwriting recognition
- Computational finance and quantitative analysis
A HMM algorithm uses 3 key components [ref 3].
- A set of observations
- A sequence of hidden states
- A model that maximizes the joint probability of the observations and hidden states, known as the Lambda model
HMM are usually visualized as lattice of states with transition and observations as illustrated in fig 2.
Fig. 2 Diagram of lattice of states, transition and observations
There are 3 use cases (or canonical forms) of the HMM [ref 4]:
- Evaluation: Evaluate the probability of a given sequence of observations, given a model.
- Training: Identify (or learn) a model given a set of observations.
- Decoding Estimate the state sequence with the highest probability to generate a given as set of observations and a model.
Viterbi algorithm
The terms "Viterbi path" and "Viterbi algorithm" are now commonly used to describe the application of dynamic programming algorithms in solving maximization problems that deal with probabilities.
Given a sequence of states {qt} and sequence of observations {oj}, the probability δt(i) for any sequence to have the highest probability path for the first T observations is defined for the state Si.
Definition of delta function (eq 1):\[\delta _{t}(i)= \max_{q_{i}:0,T-1}\ p(q_{0:T-1} =S_{i},\ O_{0:T-1} \ | \lambda )\]
Initialization of delta (eq 2): \[\delta _{0}(i)=\pi_{i}b_{i}(O_{0})\ \ \ \psi_{t}(i)=0 \ \ i:0,T-1\]
Recursive computation of delta (eq 3):\[\delta_{t}(i)=\max_{i}( \delta_{t-1}(i)a_{ij}b_{j}O_{t})\ \ \ \ \psi_{t}=arg \max_{i}(\delta_{t-1}(i).a_{ij})\]
Computation of the optimum sequence of states Q (eq 4):\[q^{*}_{t}=\psi_{t+1}(q^{*}_{t+1})\ \ \ \ q^{*}_{T}=arg \max_{i}\delta_{T}(i)\]
Scala implementation
State transition
- delta: sequence to have the highest probability path for the first i observations is defined for a specific test δt(i)
- psi Matrix that contains the indices used in the maximization of the probabilities of pre-defined states
- Qstar: the optimum sequence q* of states Q0:T-1
final class HMMState(hmmModel: HMMModel, maxIters: Int) {// Matrix of elements (t, i) that defines the highest probability of a single path // of t observations reaching state S(i) - Initialized val delta = Array.fill(hmmModel.getT)(Array.fill(hmmModel.getN)(0.0))
// Auxiliary matrix of indices that maximize the probability
//of a given sequence of states - Initialized
val psi = Array.fill(hmmModel.getT)(Array.fill(hmmModel.getN)(0))
// Singleton to compute the sequence Q* of states with
// the highest probability given a sequence of observations
object QStar {
private val qStar = Array.fill(lambda.getT)(0)
// Update Q* the optimum sequence of state using backtracking..
def update(t: Int, index: Int): Unit
...
}
def getModel: HMMModel = hmmModel
}
The class HMMModel is defined by the 3 components of the Markov processes
The algorithm is conveniently illustrated by the following diagram.
- State transition matrix A
- Observation emission matrix B
- Initial state probabilities, pi
and the number of observations, numObs.
case class HMMModel (
A: Array[Array[Double]],
B: Array[Array[Double]],
pi: Array[Double],
numObs: Int
) {
def numStates: Int = A.nRows
def numSymbols: Int = B.nCols
}
The number of states and symbols are directly extracted from the state-transition A, and observations B matrices.
Viterbi path
Fig. 3 Illustration of state-transition and observation matrices
First, let's create the key member and method for the Lambda model for HMM. The model is defined as a tuple of the transition probability matrix A, emission probability matrix B and the initial probability π
Then let's define the basic components for implementing the Viterbi algorithm. The class ViterbiPath encapsulates the member variables critical to the recursion.
The Viterbi algorithm is fully defined with
The Viterbi algorithm is fully defined with
- lambda: Lambda model as described in the previous section (line 2)
- obsIdx: Index of observed states
class ViterbiPath(state: HMMState, obsIdx: Array[Int]) { val hmmModel = state.getModel() // Access the model through its state
def getPath: HMMPrediction = { val path = recurse(1) HMMPrediction(path, qStar())
}
}
The predicted output from the Viterbi algorithm consists of
- List of predicted states
- Likelihood of the sequence of observed states
case class HMMPrediction(likelihood: Double, states: Array[Int])
Tail recursion
The recursive method, recurse that implements the formula or steps defined earlier. The method relies on the tail recursion. Tail recursion or tail elimination algorithm avoid the creation of a new stack frame for each invocation, improving the performance of the entire recursion [ref 5].
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26 | @scala.annotation.tailrec
private def recurse(t: Int): Double = {
// Initialization of the delta value for the first observation
if( t == 0)
initial // then for the subsequent observations ... else { // Update the maximum delta value and its state index for observation t Range(0, hmmModel.getN).foreach( updateMaxDelta(t, _) )
// If we reached the last observation... exit by backtracking the
// computation from the last observation
if( t == obs.size-1) {
val idxMaxDelta = Range(0, hmmModel.getN) .map(i => (i, state.delta(t, i)))
.maxBy(_._2)
// Update the Q* value with the index that maximize /the delta.A
state.QStar.update(t+1, idxMaxDelta._1)
idxMaxDelta._2
}
else
recurse(t+1)
}
}
|
Once initialized (line 5), the maximum value of delta, maxDelta, is computed through the method updateMaxDelta after the first iteration (line 10). Once the step t reaches the maximum number of observation labels, last index in the sequence of observations obs.size-1) (line 15), the optimum sequence of state q* / state.QStar is updated (line 21). The index of the column of the transition matrix A, idxMaxDelta corresponding to the maximum of delta is computed (lines 15-18). The last step is to update the matrix QStar (line 21).
The method updateMaxDelta computes the maximum value of delta of the states for this observation at time t and the state index j.
def updateMaxDelta(t: Int, j: Int): Unit = {
val (psiValue, index) = (0 until hmmModel.numStates)
.map(i => (i, delta(t - 1, i) * hmmModel.A(i, j)))
.maxBy(_._2)
psi(t)(j) = psiValue
delta += (t, j, index)
}
References
- [1] Scala Programming Language
- [2] Pattern Recognition and Machine Learning; Chap 13 Sequential Data/Hidden Markov Models - C. Bishop - Springer 2009
- [3] Scala for Machine Learning; Chap 7 Sequential Data Models - P. Nicolas - Pack Publishing 2017
- [4] Machine Learning: A Probabilistic Perspective; Chap 17 Markov and Hidden Markov models - K. Murphy - The MIT Press 2012
- [5] Wikipedia: Tail Recursion
---------------------------
Patrick Nicolas has over 25 years of experience in software and data engineering, architecture design and end-to-end deployment and support with extensive knowledge in machine learning.
He has been director of data engineering at Aideo Technologies since 2017 and he is the author of "Scala for Machine Learning" Packt Publishing ISBN 978-1-78712-238-3
He has been director of data engineering at Aideo Technologies since 2017 and he is the author of "Scala for Machine Learning" Packt Publishing ISBN 978-1-78712-238-3
