Decoding Hidden Markov Models Faster Than Viterbi
Via Online Matrix-Vector (max, +)-Multiplication
Gabriele Farina

Massimo Cairo

Dipartimento di Elettronica, Informazione e Bioingegneria
Politecnico di Milano
I-20133, Milan, Italy
gabriele2.farina@mail.polimi.it

Department of Mathematics
University of Trento
Trento, Italy
massimo.cairo@unitn.it

arXiv:1512.00077v2 [cs.LG] 11 Dec 2015

Romeo Rizzi

Department of Computer Science
University of Verona
Verona, Italy
romeo.rizzi@univr.it
Abstract
In this paper, we present a novel algorithm for the
maximum a posteriori decoding (MAPD) of timehomogeneous Hidden Markov Models (HMM), improving the worst-case running time of the classical Viterbi algorithm by a logarithmic factor. In
our approach, we interpret the Viterbi algorithm as
a repeated computation of matrix-vector (max, +)multiplications. On time-homogeneous HMMs, this
computation is online: a matrix, known in advance, has
to be multiplied with several vectors revealed one at a
time. Our main contribution is an algorithm solving
this version of matrix-vector (max, +)-multiplication
in subquadratic time, by performing a polynomial preprocessing of the matrix. Employing this fast multiplication algorithm, we solve the MAPD problem in
O(mn2 / log n) time for any time-homogeneous HMM
of size n and observation sequence of length m, with
an extra polynomial preprocessing cost negligible for
m > n. To the best of our knowledge, this is the first
algorithm for the MAPD problem requiring subquadratic
time per observation, under the only assumption – usually verified in practice – that the transition probability
matrix does not change with time.

Introduction
Hidden Markov Models (HMMs) are simple probabilistic
models originally introduced (Viterbi 1967) to decode convolutional codes. Due to their universal and fundamental nature, these models have successfully been applied in
several fields, with many important applications, such as
gene prediction (Haussler and Eeckman 1996), speech, gesture and optical character recognition (Gales 1998; Huang,
Ariki, and Jack 1990; Starner, Weaver, and Pentland 1998;
Agazzi and Kuo 1993), and part-of-speech tagging (Kupiec 1992). Their applications to bioinformatics began in
the early 1990 and soon exploded to the point that curCopyright c 2016, Association for the Advancement of Artificial
Intelligence (www.aaai.org). All rights reserved.

rently they hold a recognized place in that field (Yoon 2009;
Mäkinen et al. 2015).
A HMM describes a stochastic process generating a sequence of observations y1 , y2 , . . . , yn . Internally, a sequence of hidden states x1 , x2 , . . . , xn is generated according to a Markov chain. At each time instant t = 1, 2, . . . , n,
a symbol yt is observed according to a probability distribution depending on xt . We consider only time-homogeneous
HMMs, i.e. models whose parameters do not depend on
the time t. While this assumption covers the majority
of applications, some notable exceptions involving timeinhomogeneous models are known (Lafferty, McCallum,
and Pereira 2001).
Maximum a posteriori decoding (MAPD). Since the states
of the model are hidden, i.e. only the generated symbols can
be observed, a natural problem associated with HMMs is
the MAPD problem: given a HMM M and an observed sequence of symbols Y of length m, find any state path X
through M maximizing the joint probability of X and Y .
We call any such X a most probable state path explaining the
observation Y . Traditionally, the MAPD problem is solved
by the Viterbi algorithm (Viterbi 1967), in O(mn2 ) time and
O(mn) memory for any model of size n and observation sequence of length m.
Over the years, much effort has been put into lowering the
cost of the Viterbi algorithm, both in terms of memory and of
running time. (Grice, Hughey, and Speck 1997) showed that
a checkpointing technique√can be employed to reduce the
memory complexity to O( m · n); refinements of this idea
(embedded checkpointing) deliver a family of time-memory
tradeoffs, culminating into an O(n log m) memory solution
with a slightly increased running time O(mn2 log m).
At the same time, several works reducing the time complexity of the algorithm in the average-case were developed (Šrámek 2007; Churbanov and Winters-Hilt 2008;
Felzenszwalb, Huttenlocher, and Kleinberg 2004; Esposito
and Radicioni 2009; Kaji et al. 2010). Many of these works
make assumptions on the structure of T , and may lose the
optimality or degenerate to the worst case Θ(mn2 ) opera-

tions when these assumptions are not fulfilled. In (Lifshits
et al. 2009), the authors show a method to speed up the decoding of HMMs by a O(log m) factor, by precomputing all
possible observation sequences of length log m, in a fashion similar to the Four Russians method. This requires the
number of such sequences to be “small”. In the same work,
the authors also show that it is possible to compress the observation sequence to achieve speedups proportional to the
compression ratios. However, this latter method seems to
require the observation sequence to be available in advance.
To the best of our knowledge no algorithm achieving a
worst-case running time better than O(mn2 ) is known under
the only assumption of time-homogeneousness.
Approach. We give an algorithm solving the MAPD problem for time-homogeneous HMMs with time complexity
asymptotically lower than O(mn2 ), in the worst case. We
regard the MAPD problem as an iterated computation of a
matrix-vector multiplication. For time-homogeneous models, the matrix is known in advance and does not change with
time. However, the sequence of vectors to be multiplied cannot be foreseen, as each vector depends on the result of the
previous computation; this rules out the possibility to batch
the vectors into a matrix and defer the computation. We call
this version of the problem, in which a fixed matrix has to be
multiplied with several vectors revealed one at a time, “the
online matrix-vector multiplication (OMV M UL) problem”.
Consider the problem of multiplying a n × n matrix with
a column vector of size n. Without further assumptions,
the trivial O(n2 ) time algorithm is optimal, since all the
n2 elements of the matrix have to be read at least once.
However, under the assumption that the matrix is known in
advance and can be preprocessed, this trivial lower bound
ceases to hold. Algorithms faster than the trivial quadratic
one are known for the OMV M UL problem over finite semirings (Williams 2007), as well as over real numbers with standard (+, ·)-multiplication, if the matrix has only a constant
number of distinct values (Liberty and Zucker 2009).
However, none of the above algorithm can be applied to time-homogeneous HMMs, as their decoding relies on online real matrix-vector (max, +)-multiplication
(ORMV (max, +)-M UL). In the specific case of real
(max, +)-multiplication, subcubic algorithms have been
known for years (Dobosiewicz 1990; Chan 2008; 2015) for
the matrix-matrix multiplication problem, with important
applications to graph theory and boolean matrix multiplication, among others. However, we are not aware of any
algorithm solving the ORMV (max, +)-M UL problem in
subquadratic time. Note that the ORMV (max, +)-M UL
can be used to compute the OMV M UL over the Boolean
semiring: for this problem, it has been conjectured (Henzinger et al. 2015) that no “truly polynomially subquadratic”
algorithm1 exists for the ORMV (max, +)-M UL problem.
We reduce the ORMV (max, +)-M UL problem to a
multi-dimensional geometric dominance problem, following an approach similar to that of (Bremner et al. 2006;
Chan 2008). Then, the geometric problem is solved by a
1
That is, running in time O(n2−ε ) for some ε > 0 after a polynomial preprocessing or the matrix.

divide-and-conquer algorithm, which can be regarded as a
transposition of the algorithm of (Chan 2008) to the online
setting. Our technique yields a worst-case O(mn2 / log n)
algorithm, called GDFV, solving the MAPD problem after a
polynomial preprocessing of the model.
Contributions. Our key contributions are as follows: (i)
we extend the geometric dominance reporting problem introduced in (Chan 2008) to the online setting; (ii) we solve
the ORMV (max, +)-M UL problem in O(n2 / log n) time
after a polynomial preprocessing of the n × n matrix; (iii)
we show an algorithm solving the MAPD problem on timehomogeneous HMMs in O(mn2 / log n) time in the worstcase, after a polynomial preprocessing of the model.
Finally, we experimentally evaluate the performance of
our algorithms, with encouraging results. Currently the
problem sets in which we outperform Viterbi are limited,
but we hope that the approach we propose will open the way
to further improvements on this problem in future works.

Preliminaries
Notation
The i-th component of a vector v is denoted by v[i]; similarly, M[i, j] denotes the entry of row i and column j, in
matrix M. Indices will always be considered as starting
from 1. Given two vectors a and b of dimension n, such
that a[i] ≤ b[i] for every coordinate index i = 1, . . . , n, we
write a  b and say that b dominates a, or, equivalently,
that (a, b) is a dominating pair.
Given a matrix or vector M with non-negative entries, we
write log M to mean the matrix or vector that is obtained
from M by applying the logarithm on every component. We
will almost always work with the extended set R∗ = R ∪
{−∞}, so that we can write log 0 = −∞. We assume that
−∞ + x = x + (−∞) = −∞ and x ≥ −∞ for all x ∈ R∗ .

Hidden Markov Models (HMMs)
We formally introduce the concept of time-homogeneous
Hidden Markov Models.
Definition 1. A time-homogeneous HMM is a tuple M =
(S, A, Π, T , E), composed of:
• a set S = {s1 , . . . , sn } of n hidden states; n is called the
size of the model,
• an output alphabet A = {a1 , . . . , a|A| },
• a probability distribution vector Π = {π1 , . . . , πn } over
the initial states,
• a matrix T = {ts (s0 )}s,s0 ∈ S of transition probabilities
between states,
• a matrix E = {es (a)}as ∈∈ SA of emission probabilities.
Matrices T and E are stochastic, i.e., the entries of every
row sum up to 1.
For notational convenience, we relabel the states of a
HMM with natural numbers, i.e. we let S = {1, . . . , n}.
As stated in the introduction, HMMs define generative
processes over the alphabet A. The initial state x0 ∈ S
is chosen according to the distribution Π; then, at each step,
a symbol a is generated according to the probability distribution ex (a), where x is the current state; a new state x0 is

chosen according to the probability distribution induced by
tx (x0 ), and the process repeats. The probability of a state
path X = (x1 , . . . , xm ) joint to an observation sequence
Y = (y1 , . . . , ym ) is computed as:
! m
!
m−1
Y
Y
Pr(X, Y ) = πx1
txi (xi+1 )
exi (yi ) .
i=1

i=1

The Viterbi algorithm
The Viterbi algorithm consists of two phases: in the first
phase, a simple dynamic programming approach is used to
determine the probability of the most probable state path
ending in each possible state. In the second phase, the data
stored in the dynamic programming table is used to reconstruct a most probable state path.
Definition 2. Assume given a HMM M = (S, A, Π, T , E)
and an observed sequence A = (a1 , . . . , am ). For every
s ∈ S and i = 1, . . . , m, denote by qi (s) the probability
of any most probable path ending in state s explaining the
observation Ai−1 = (a1 , . . . , ai−1 ).
By definition of qi (s), any most probable path explaining
A has probability maxs∈S {es,am · qm (s)}. The qi (s) values can be computed inductively. Indeed, q1 (s) = πs for all
s ∈ S, while for every i > 1 and s ∈ S it holds:
qi (s) = max
{qi−1 (s0 ) · ts0 (s) · es0 (ai−1 )} .
(1)
0
s ∈S

In order to compute all the n values qi (s), for any fixed i >
1, Θ(n2 ) comparisons have to be performed. This phase is
in fact the bottleneck of the algorithm.
The second phase of the algorithm uses the qi (s) values to reconstruct an optimal path in O(mn) time. We
will not deal with this second and faster part, and only
mention that most of the previously developed solutions
for it, including the memory saving ones (Šrámek 2007;
Churbanov and Winters-Hilt 2008), are still applicable once
the first part has been carried out based on our approach.

Online matrix-vector multiplication problem
We formally introduce the online matrix-vector (max, +)multiplication problem briefly discussed in the introduction.
Definition 3 (Matrix-vector (max, +)-multiplication).
Given a m×n matrix A and a n-dimensional column vector
b over R∗ , their (max, +)-multiplication A ∗ b is a column
vector of dimension m whose i-th component is:
n
(A ∗ b)[i] = max {A[i, j] + b[j]} .
j=1

Problem 1 (ORMV (max, +)-M UL problem). Given a
m × n matrix A over R∗ , perform a polynomial-time preprocess of it, so that the product A ∗ b can be computed for
any input vector b of size n.
Lemma 1 is a simple result that bridges the gap between
the MAPD problem and the ORMV (max, +)-M UL problem, showing that any fast solution to Problem 1 can be
turned into a fast solution for the MAPD problem.
Lemma 1. Any algorithm for Problem 1 computing A ∗ b
in T (m, n) time after U (m, n) preprocessing time can be
used to solve the MAPD problem for any time-homogeneous

HMM of size n and observation sequence of length m in
O(m · T (n, n)) time after a U (n, n) time preprocessing.
Proof. Denote by t T the transpose of the state transition
matrix T of M. Furthermore, for every i = 1, . . . , m, introduce the pair of vectors
qi = (qi (1), . . . , qi (n)), ei = (e1 (ai ), . . . , en (ai )).
The value qi (s) corresponding to instant i > 0 and state
s ∈ S can be computed as follows, in logaritmic scale:
log qi (s) = log max
{qi−1 (s0 ) · ts0 (s) · es0 (ai−1 )}
0
s ∈S

= max
{log(ts0 (s)) + log(qi−1 (s0 ) es0 (ai−1 ))}
0
s ∈S

= ((log t T ) ∗ (log qi−1 + log ei−1 ))[s].
Notice that the n × n matrix log t T depends only on the
model and is time-invariant; therefore, we can compute
qi (s) for all s = 1, . . . , n in batch with an instance of
ORMV (max, +)-M UL:
log qi = (log t T ) ∗ (log qi−1 + log ei−1 ).
Once the qi (s) values have been computed, we can use
the second part of the Viterbi algorithm out of the box. The
time required for the multiplication – the bottleneck of the
algorithm – is O(T (n, n)) by hypothesis, hence the final algorithm has time complexity O(m·T (n, n)). The time complexity of the preprocessing is U (n, n).

From multiplication to geometric dominance
Consider the following geometric problem, which apparently has no relation with the (max, +)-multiplication, nor
with the Viterbi algorithm.
Problem 2 (Online geometric dominance reporting).
Let B be a set of d-dimensional vectors2 . Given a vector
p ∈ B, we define its domination set as
δB (p) = {b ∈ B : b  p}.
Preprocess B so that at a later time the set δB (p) can be
computed for any input vector p.
Lemma 2. Any algorithm solving Problem 2 in
O(T (d, |B|) + |δB (p)|) time after a preprocessing
time O(U (d, |B|)) can be turned into an algorithm solving
the ORMV (max, +)-M UL problem for any m × t matrix
in O(m + t · T (t, m)) time, with preprocessing time
O(mt2 + t · U (t, m)).
Proof. This constructive proof comes in different stages.
First, the result is obtained under the simplifying assumptions that (i) neither A nor b have any −∞ entry, and (ii)
the maximum sum on any row, A[·, j] + b[j], is achieved
by exactly one value of the column index j. Later, these
conditions will be dropped.
Observe that (A ∗ b)[i] = A[i, j ∗ ] + b[j ∗ ] iff
A[i, j ∗ ] + b[j ∗ ] ≥ A[i, j] + b[j]
(2)
for every column index j. Under assumption (i), this inequality can be rewritten as
A[i, j] − A[i, j ∗ ] ≤ b[j ∗ ] − b[j].
2

The coordinates of the vectors can range over any chosen
totally-ordered set, as long as any two coordinate values can be
compared in constant time.

Defining the values
ai,j ∗ (j) = A[i, j] − A[i, j ∗ ], bj ∗ (j) = b[j ∗ ] − b[j]
for all the feasible values of i, j, j ∗ , we obtain
(A ∗ b)[i] = A[i, j ∗ ] + b[j ∗ ] ⇐⇒ ai,j ∗ (j) ≤ bj ∗ (j) ∀j.
(3)
Notice that the last expression is actually a statement
of geometric dominance, between the two t-dimensional
vectors Ãi,j ∗ = (ai,j ∗ (1), . . . , ai,j ∗ (t)) and b̃j ∗ =
(bj ∗ (1), . . . , bj ∗ (t)). This immediately leads to the following algorithm.
Algorithm 1 m × t matrix-vector (max, +)-multiplication
1: procedure PREPROCESS(A)
2:
for j ∗ = 1, . . . , t do
3:
for i = 1, . . . , m do
4:
Ãi,j ∗ ← (ai,j ∗ (1), . . . , ai,j ∗ (t))
∗
5:
Bj ← {Ã1,j ∗ , . . . , Ãm,j ∗ }
6:
Preprocess Bj ∗
1: procedure MULTIPLY(b)
. Returns A ∗ b
2:
for j ∗ = 1, . . . , t do
3:
b̃j ← (bj ∗ (1), . . . , bj ∗ (t))
4:
δ ← δBj ∗ (b̃j ∗ ), as defined in Problem 2
5:
for all Ãi,j ∗ ∈ δ do
6:
m[i] ← A[i, j ∗ ] + b[j ∗ ]
7:
return m

Line 6 in procedure PREPROCESS requires U (t, m) time,
making the total preprocessing cost of MULTIPLY O(mt2 +
t · U (t, m)). As for the running time, Line 4 of MULTI PLY takes O(T (t, m) + |δ|) time, and is executed t times.
Under assumption (ii), there is only one column j ∗ that satisfies Equation 3 for each row i, hence the total number of
elements appearing in δ is exactly m. As a consequence,
Line 6 of MULTIPLY is executed m times and the total time
complexity of MULTIPLY is O(m + t · T (m, t)).
We now relax assumptions (i) and (ii) by applying a transformation of the input. Instead of working on R∗ , we work
on triples over N × R × N. The matrix A is transformed
as follows: each element x > −∞ is replaced by h0, x, 0i,
while each occurence of −∞ is replaced by h−1, 0, 0i. The
input vector b is transformed similary, but the third coordinate is used to hold the index of the replaced element.
Namely, element b[j] = x is replaced by h0, x, ji if x >
−∞ and by h−1, 0, ji otherwise. One can informally regard
the first two entries of each triple ha, b, ·i as a shorthand for
the value a · ∞ + b. Any two triples are compared according
to lexicographical order, while addition and subtraction are
performed element-wise.
The crucial observation is that, if A[i, j ∗ ] + b[j ∗ ] >
A[i, j] + b[j] before the transformation, then the same holds
also after the transformation. Hence, we can solve the transformed problem to obtain the solution of the original problem. Our algorithm can be applied as it is to the transformed
problem: indeed, the inequality in Equation 2 can be rearranged without any further assumption; moreover, there can
be no two distinct columns j and j 0 achieving the maximum,
as the two triples A[i, j] + b[j] and A[i, j 0 ] + b[j 0 ] differ at

least on the third element. Once the output vector is obtained, replace each triple hk, x, ji with x if k = 0 and with
−∞ otherwise. Conveniently, the third element j holds the
index of the column achieving the maximum.
Lemma 3 shows that every fast algorithm for the
OMV (max, +)-M UL of narrow rectangular matrices can
be turned into a OMV (max, +)-M UL algorithm for square
matrices.
Lemma 3. Any algorithm computing the OMV (max, +)M UL of a m × t matrix in T (m, t) time and U (m, t) preprocessing time, can be used to multiply any m × n matrix,
n ≥ t, in O(n/t · (T (m, t) + m)) time and O(n/t · U (m, t))
preprocessing time.
Proof. Assume without loss of generality that n is an integer
multiple of t (otherwise, add columns to A and elements to
b with value −∞ until the condition is met). The idea is to
split A and b into n/t blocks, each of size m × t and t × 1
respectively:
A = (A1 | · · · |An/t ), b = (b1 | · · · |bn/t ).
n/t

Observe that (A ∗ b)[i] = max`=1 {(A` ∗ b` )[i]} . This immediately leads to the following algorithm. First, we preprocess each block A` in U (m, t) time with the given algorithm, so that the product A` ∗ b` can be later computed in
T (m, t) time. As soon as the vector b is received, compute
m` = A` ∗ b` for all ` = 1, . . . , n/t, and finally the output
n/t
vector by (A ∗ b)[i] = max`=1 {m` [i]}.
The time analysis is straightforward. The computation of
each m` takes T (m, t) time. There are n/t such computations and merging the results takes O(m·n/t) time, yielding
a total time of O(n/t·(T (m, t)+m)). The total preprocessing time is O(n/t · U (m, t)).

Before continuing with the next section, where the main
result of this paper will be discussed, we state the following
theorem, whose proof is available in the Appendix of this
paper.
Theorem 1. Problem 2 can be solved in O(d log |B| +
|δB (p)|) time per query, after a O(|B|d+1 ) time and space
preprocessing.

An O(n2 / log n) algorithm (GDFV)

Lemma 4. There exists an algorithm solving Problem 2 in
O(d cdε |B|ε + |δB (p)|) time for every ε ∈ (0, 1], where
d
cε := 1/(2ε − 1). The preprocessing requires O(c0ε |B|1+ε )
time and memory for every ε ∈ (0, log2 3/2], where c0ε :=
1/(21+ε − 2).
Proof. We develop a simple divide-and-conquer algorithm;
we assume without loss of generality that |B| is a power of
two.
Overview. If d = 0, return δB (p) = B. If B contains only
one vector b, check if b  p and return either {b} or the
empty set accordingly. In all the other cases, split B into two
sets B − and B + of size |B|/2, according to the median d-th
coordinate γ of the vectors in B, so that b− [d] ≤ γ ≤ b+ [d]
for all b− ∈ B − and b+ ∈ B + . Now consider the d-th coordinate of the query vector p: if it is strictly less than γ, then

all the vectors in B + do not occur in the solution. Hence,
solve the problem recursively on B − . Otherwise, both the
sets B + and B − need to be considered; however, the d-th
coordinate for the vectors in B − is known to be ≤ p[d] and
can be dropped. Hence, solve the problem recursively on
B + and p and on (B − )0 and p0 , where the apostrophe denotes the discard of the last coordinate, and merge the solutions. The recursive step is summarized by the following
recurrence:

δ − (p)
if p[d] < γ,
δB (p) = B
0
−
0
+
δB (p) ∪ δ(B ) (p ) otherwise.

In order to make the algorithm faster, we exploit the fact
that B is known in advance. At preprocessing time, we
build a tree that guides the execution of the algorithm, where
each node u corresponds to a subproblem Bu over a du dimensional space. The root corresponds to the original set
B. If |Bu | ≥ 2 and du ≥ 1, then the node u stores the
median value γ and has three children, corresponding to the
subproblems Bu+ , Bu− and (Bu− )0 . Otherwise, u is a leaf storing the content of Bu . We analyze the cost of building the
tree later on. For now, notice that the size of the tree is at
most polynomial in |B|: the height is at most log |B|, as
the value |Bu | halves at each level, so the nodes are at most
O(3log2 |B| ) = O(|B|log2 3 ) = O(|B|1.59 ).
Time analysis. Our algorithm starts from the root node,
and visits recursively the nodes in the tree that are needed
to solve the problem. When we reach a leaf u with du =
0, we output Bu (which is not empty) in O(|Bu |) time. If
instead du > 0 and |Bu | = 1, we pay O(d) time to check if
b  p, and O(1) to output b if needed. On internal nodes,
we only pay constant extra time as the median coordinate γ
is known from the tree. The cost of producing the output is
O(|δB (p)|), and is measured separately. Hence, the running
time is O(Td (|B|)+|δB (p)|) where Td (n) satisfies the linear
recurrence relation

Td (n/2)
Td (n) = 1 + max
Td−1 (n/2) + Td (n/2)
(4)
= Td−1 (n/2) + Td (n/2) + 1,
with base cases Td (1) = d and T0 (n) = 0. (The time required to handle this last case is included in O(|δB (p)|)).
We show by induction that
Td (n) ≤ T d (n) := d cdε nε ,
for any chosen ε ∈ (0, 1], where cε := 1/(2ε − 1). Notice
that cε ≥ 1 for ε ∈ (0, 1], thus the statement is true for the
base cases as T d (n) ≥ d. Assuming the inductive hypothesis, we obtain for n ≥ 2 and d ≥ 1:
Td (n) ≤ T d (n/2) + T d−1 (n/2) + 1
= d cdε (n/2)ε + (d − 1) cd−1
(n/2)ε + 1
ε

≤ d cdε (n/2)ε + d cd−1
(n/2)ε
ε

1 + c−1
1 + 2ε − 1
ε
=
T
(n)
= T d (n),
d
2ε
2ε
completing the induction. Thus, the time complexity of the
algorithm is O(d cdε |B|ε + |δB (p)|).
Preprocessing. The tree is built starting from the root.
Finding the median d-th coordinate γ, computing Bu+ , Bu−
= d cdε nε

and (Bu− )0 , and storing the data in the node u, all require
O(|Bu |) time and memory. Hence, the time and memory
cost to build the tree is O(Ud (|B|)), where Ud (n) satisfies
the recurrence
Ud (n) = 2Ud (n/2) + Ud−1 (n/2) + n
with base cases U0 (n) = n and Ud (1) = 1. We show by
induction that
d
Ud (n) ≤ U d (n) := 3c0ε n1+ε − 2n
for any chosen ε ∈ (0, log2 3/2], where c0ε := 1/(21+ε −
2). Notice that c0ε ≥ 1 for ε ∈ (0, log2 3/2]. Hence, the
d
statement is true for the base cases, as U d (n) ≥ 3c0ε n1+ε −
2n ≥ 3n − 2n ≥ n. Assuming the inductive hypothesis, we
obtain for n ≥ 2 and d ≥ 1:
Ud (n) ≤ 2 U d (n/2) + U d−1 (n/2) + n
d

= 2 · (3 c0ε (n/2)1+ε − n)
+ (3 c0ε
d

d−1

(n/2)1+ε − n) + n
d−1

= 2 · 3 c0ε (n/2)1+ε + 3 c0ε
(n/2)1+ε − 2n
−1
2 + c0ε
d
= 3 c0ε n1+ε ·
− 2n
21+ε
d
= 3 c0ε n1+ε − 2n = U d (n).
completing the induction. Hence, the time and memory cost
d
of the preprocessing phase is O(Ud (|B|)) = O(c0ε nε ).

We remark that the time complexity and the recurrence
of the simple algorithm given in the above proof differ from
those of Chan by necessity, as the online setting requires
each vector to be treated separately. On the other hand, his
result follows directly from ours:
Theorem 2 ((Chan 2008), Lemma 2.1). Given n red/blue
points in Rd∗ we can report all K dominating pairs in
O(kεd n1+ε + K) time for any ε ∈ (0, 1), where kε :=
2ε /(2ε − 1).
Proof. After reading the set B of blue points, preprocess
them as described in the proof of Lemma 4. Then, for each
red point p, perform a query to find all the dominators of p,
i.e. δB (p); simply flush out the union of all the dominating
pairs obtained. By Lemma 4, the cost of the preprocessing is O((21+ε − 2)−d n1+ε ) = O(kεd n1+ε ). On the other
hand, each of the n queries takes time O(d(2ε − 1)−d nε ),
that is3 O(kεd nε ), excluding the output; the overhead due to
the actual output of the pairs is O(K). The final cost of the
algorithm is therefore O(cdε n1+ε + K) as desired.
Finally, Lemmas 1, 4, and 3 combine into the following.
Theorem 3. There exists an algorithm solving the MAPD
problem in O(mn2 / log n) time after a polynomial preprocessing of the model, for any HMM of size n and observation
sequence of length m.
Proof. It is enough to show how to solve the
ORMV (max, +)-M UL problem for any n × n matrix in O(n2 / log n) time. To this end, apply Lemma 3 with
m = n and t = α log2 n, where α ∈ (0, 1/2) ⊆ R. This
3

Indeed, using the binomial expansion formula we have:

d

d−1
1
1
kεd = 1 + ε
≥d
= Ω(d cdε ).
2 −1
2ε − 1

gives a running time of O(n · T (α log2 n, n) + n2 / log n),
where T (d, n) is the cost of computing the domination set
of a d-dimensional vector over a fixed set of n points (see
Problem 2). Substituting the bounds of Lemma 4 into the
time complexity, yields a polynomial preprocessing cost,
and the following time bound for each multiplication:
O((2ε − 1)−α log2 n αn1+ε log2 n + n2 / log n) =
ε

= O(n1+ε−α log2 (2 −1) log2 n + n2 / log n)
for all ε ∈ (0, 1). Setting ε = 2α, the exponent of the first
term becomes 1 + α(2 − log2 (4α − 1)) < 2 for all 0 <
α < 1/2. Therefore, the time complexity of the algorithm is
O(n2 / log n).
We call the resulting algorithm geometric dominance
faster Viterbi (GDFV).

6
5
4
3
2
1
6
5
4
3
2
1

6
5
4
3
2
1

t =2

Trivial
50

150

200

250

300

350

400

t =3

Trivial
400

600

800

1000

1200

1400

1600

1800

2000

t =4

Trivial
5000

Experimental evaluation
Methodology. All the algorithms are implemented in the
C++11 language, compiled using the clang compiler and
run on the OSX 10.10.3 operating system. The main memory is a 8GB 1600MHz DDR3 RAM, and the processor is
an Intel Core i7-4850HQ CPU, with 6MB shared L3 cache.
All the matrices and vectors used for the experiments have
entries sampled from a uniform distribution over (0, 1] ⊆ R.
Results. How does our proposed ORMV (max, +)-M UL
algorithm for narrow matrices compare to the trivial one?
We analyze the throughput of Algorithm 1, based on the geometric subroutines exposed in the proof of Lemma 4, comparing it with the trivial multiplication approach. For every
chosen pair (n, t), we run 25 tests, each of which consists of
an online multiplication of a n × t matrix with 10 000 vectors. The results of our tests are summarized in Figure 1,
where we see that our algorithm can be up to 4 times faster
than the trivial one. This is mainly due to the fact that the
number of accesses to the tree is much less than n · t, and
to the lower number of comparisons needed to find the answer. See the Appendix for a more thorough analysis of the
average number of accesses to the decision trees.
How does the complete GDFV algorithm compare with
the Viterbi algorithm? We experimentally evaluate the first
phase of the GDFV algorithm, i.e. the computation of the
qi (s) values defined in Equation 1. This is the most expensive task in the decoding of HMMs. We implement the
algorithm as described in the proof of Theorem 3, using
α = 0.25, that is splitting the n × n transition probability
matrix T of the model in approximately n/2 blocks when
n ≤ 4000. We summarize the results in Figure 2, where we
see that our algorithm is roughly twice as fast as the Viterbi
algorithm, in line with expectations. However, we note that
the amount of memory required by our algorithm makes it
impractical for larger values of α. Indeed, we have verified that when the memory pressure becomes high other factors slow down the implemented algorithm, such as cache
and page misses, or, for bigger allocations, the hard drive
latency.

100

10000
Number of rows (n)

15000

20000

Figure 1: Relative throughput of Algorithm 1 when compared to the trivial quadratic algorithm on matrices of size
n×t for t = 2, 3, 4. A higher throughput implies faster computation. Legend: dotted black lines denote min and max
values, red solid lines denote median values, black dashed
lines denote mean values (µ), and gray shading denotes the
range µ ± σ, where σ is the standard deviation.
2.0
1.8
1.6
1.4
1.2

Viterbi

1.0
20

30

40

50
60
Model size (n)

70

80

90

Figure 2: Relative throughput of the GDFV algorithm when
compared to the Viterbi algorithm. A higher throughput implies faster computation. Legend: as in Figure 1.

Conclusion and future works
In this paper, we give the first algorithm for the maximum a
posteriori decoding (MAPD) of time-homogeneous Hidden
Markov Models requiring asymptotically less than O(mn2 )
operations in the worst case. To this end, we first introduce
an online geometric dominance reporting problem, and propose a simple divide-and-conquer solution, generalizing the
classical result by (Chan 2008). At an intermediate step, we
also give the first algorithm solving the online matrix-vector
(max, +)-multiplication problem over R∗ in subquadratic
time after a polynomial preprocessing of the matrix. Finally,
we apply the faster multiplication to the MAPD problem.
Furthermore, we think that our proposal paves the way to
several unexplored questions which we intend to explore in
future works:
• cut larger polylogarithmic factors, by splitting cases in
Equation 4 in a different manner, as in (Chan 2015);
• study and implement a more succinct version of the deci-

sion tree, in order to mitigate the memory footprint;
• analyze the relationship of our work with other existing
heuristics, such as CarpeDiem (Esposito and Radicioni
2009);
• combine our speed-up to the one delivered by the approach in (Lifshits et al. 2009). Notice that this would
require further assumptions on the observation sequence;
• we note that the decision trees built by our algorithm
could be implemented at a hardware level, resulting
in specialized chips performing asymptotically less that
O(n2 ) operations per observed symbol, in the worst case.
• investigate “truly polynomially subquadratic” solutions
for the MAPD problem, at the expense of an exponential preprocessing of the model. As a concrete example,
we present the following theorem, which is a corollary of
Theorem 1 (a formal proof can be found in the Appendix):
Theorem 4. The MAPD problem on
√ time–homogeneous
3/2
HMMs
can
be
solved
in
O(mn
log n) time with an


√
O n1+ n/ log n time and space preprocessing.

Gales, M. J. 1998. Maximum likelihood linear transformations for hmm-based speech recognition. Computer speech
& language 12(2):75–98.

Acknowledgements

Kaji, N.; Fujiwara, Y.; Yoshinaga, N.; and Kitsuregawa, M.
2010. Efficient staggered decoding for sequence labeling.
In Proceedings of the 48th Annual Meeting of the Association for Computational Linguistics, 485–494. Association
for Computational Linguistics.

Massimo Cairo was supported by the Department of Computer Science, University of Verona under PhD grant “Computational Mathematics and Biology”.
We would like to thank Marco Elver, Nicola Gatti, Zu
Kim, and Luigi Laura for their valuable suggestions.

References
Agazzi, O. E., and Kuo, S.-s. 1993. Hidden markov model
based optical character recognition in the presence of deterministic transformations. Pattern recognition 26(12):1813–
1826.
Bremner, D.; Chan, T. M.; Demaine, E. D.; Erickson, J.;
Hurtado, F.; Iacono, J.; Langerman, S.; and Taslakian, P.
2006. Necklaces, convolutions, and x+ y. In Algorithms–
ESA 2006. Springer. 160–171.
Chan, T. M. 2008. All-pairs shortest paths with real weights
in o (n 3/log n) time. Algorithmica 50(2):236–243.
Chan, T. M. 2015. Speeding up the four russians algorithm by about one more logarithmic factor. In Proceedings
of the Twenty-Sixth Annual ACM-SIAM Symposium on Discrete Algorithms, 212–217. SIAM.
Churbanov, A., and Winters-Hilt, S. 2008. Implementing
em and viterbi algorithms for hidden markov model in linear
memory. BMC bioinformatics 9(1):224.
Dobosiewicz, W. 1990. A more efficient algorithm for the
min-plus multiplication. International journal of computer
mathematics 32(1-2):49–60.
Esposito, R., and Radicioni, D. P. 2009. Carpediem: Optimizing the viterbi algorithm and applications to supervised
sequential learning. The Journal of Machine Learning Research 10:1851–1880.
Felzenszwalb, P. F.; Huttenlocher, D. P.; and Kleinberg,
J. M. 2004. Fast algorithms for large-state-space hmms
with applications to web usage analysis. Advances in NIPS
16:409–416.

Grice, J.; Hughey, R.; and Speck, D. 1997. Reduced
space sequence alignment. Computer applications in the
biosciences : CABIOS 13(1):45–53.
Haussler, D. K. D., and Eeckman, M. G. R. F. H. 1996. A
generalized hidden markov model for the recognition of human genes in dna. In Proc. Int. Conf. on Intelligent Systems
for Molecular Biology, St. Louis, 134–142.
Henzinger, M.; Krinninger, S.; Nanongkai, D.; and Saranurak, T. 2015. Unifying and strengthening hardness for
dynamic problems via the online matrix-vector multiplication conjecture. In Proceedings of the Forty-Seventh Annual
ACM on Symposium on Theory of Computing, STOC ’15,
21–30. New York, NY, USA: ACM.
Huang, X. D.; Ariki, Y.; and Jack, M. A. 1990. Hidden
Markov models for speech recognition, volume 2004. Edinburgh university press Edinburgh.

Kupiec, J. 1992. Robust part-of-speech tagging using
a hidden markov model. Computer Speech & Language
6(3):225–242.
Lafferty, J. D.; McCallum, A.; and Pereira, F. C. N. 2001.
Conditional random fields: Probabilistic models for segmenting and labeling sequence data. In Proceedings of the
Eighteenth International Conference on Machine Learning,
ICML ’01, 282–289. San Francisco, CA, USA: Morgan
Kaufmann Publishers Inc.
Liberty, E., and Zucker, S. W. 2009. The mailman algorithm: A note on matrix–vector multiplication. Information
Processing Letters 109(3):179–182.
Lifshits, Y.; Mozes, S.; Weimann, O.; and Ziv-Ukelson, M.
2009. Speeding up hmm decoding and training by exploiting
sequence repetitions. Algorithmica 54(3):379–399.
Mäkinen, V.; Belazzougui, D.; Cunial, F.; and Tomescu,
A. I. 2015. Genome-Scale Algorithm Design. Cambridge
University Press.
Šrámek, R. 2007. The on-line viterbi algorithm. KAI FMFI
UK, Bratislava, máj.
Starner, T.; Weaver, J.; and Pentland, A. 1998. Real-time
american sign language recognition using desk and wearable computer based video. Pattern Analysis and Machine
Intelligence, IEEE Transactions on 20(12):1371–1375.
Viterbi, A. J. 1967. Error bounds for convolutional codes
and an asymptotically optimum decoding algorithm. Information Theory, IEEE Transactions on 13(2):260–269.
Williams, R. 2007. Matrix-vector multiplication in subquadratic time:(some preprocessing required). In SODA,
volume 7, 995–1001.

Yoon, B.-J. 2009. Hidden markov models and their applications in biological sequence analysis. Current genomics
10(6):402.

Appendix
Missing proofs
We prove Theorem 1 and Theorem 4.
Theorem 1. Problem 2 can be solved in O(d log |B| +
|δB (p)|) time per query, after a O(|B|d+1 ) time and space
preprocessing.
Proof. For every coordinate k = 1, . . . , d, define the sequence bk1 , . . . , bk|B| containing all the vectors in B ordered
by their k–th coordinate. Namely:
bk1 [k] ≤ bk2 [k] ≤ · · · ≤ bk|B| [k].
Given an input vector p, let rk ∈ {0, . . . , |B|} be the
last position in which p can be inserted in the sequence
bk1 , . . . , bk|B| , while maintaining its ordering according to the
k–th coordinate. That is, rk is the only index that satisfies:
bk1 [k], . . . , bkrk [k] ≤ p[k] < bkrk +1 [k], . . . , bk|B| [k].
Now, observe that the solution to a query can be obtained as:
δB (p) = {b ∈ B : b[k] ≤ p[k] ∀k}
\
= {bk1 , · · · , bkrk }.
k

The algorithm works as follows. At preprocessing time, we
order the vectors according to each coordinate k = 1, . . . , d.
Then, we precompute
and store in a lookup table the solution
T
δB (p) = k {bk1 , · · · , bkrk } for each of the |B + 1|d possible choices of the values r1 , . . . , rd . This takes O(|B|d+1 )
time and space, and dominates the cost of the preprocessing
phase.
Once the input vector p is received, rk is computed by binary search for each coordinate k = 1, . . . , d. Then, the
solution for the values r1 , . . . , rd just obtained is looked
up in the table. The time for lookup is dominated by the
O(d log |B|) time of the binary searches, yielding a total running time of O(d log |B| + |δB (p)|) to compute and return
the output.
Theorem 4. The MAPD problem on
√ time–homogeneous
3/2
HMMs
can
be
solved
in
O(mn
log n) time with an


√
O n1+ n/ log n time and space preprocessing.

Proof. We apply consecutively Theorem 1, Lemma 2,
Lemma 3 and Lemma 1. First, by Theorem 1, we solve
the online dominance reporting problem in O(d log |B|)
time. Second, we apply Lemma 2 to this algorithm, and
obtain a solution for the online matrix–vector (max, +)multiplication for a rectangular m × t matrix running in
O(m + t2 log |m|) time. Third, by Lemma 3, this algorithm can be used to solve the multiplication for n × n
2
matrices in O(n/t · (n +
pt log n)) time for any choice
of t. Choosing t =
n/ log n, the running time is
√
3/2
O(n n log n) =
O(n
log1/2 n) with a preprocess√
ing cost of O(n n/ log n+1 ) time and space. Finally, by
Lemma 1, we employ this algorithm to solve the
√ MAPD
problem in O(mn3/2 log1/2 n) time, with O(n1+ n/ log n )
time and space required for the preprocessing.

Experimental evaluation
7

Min, Max
Mean (µ)
Median
µ ±σ

6
5
4
3
2

Trivial

1
50

100

150

200
250
Number of rows (n)

300

350

400

Figure 3: Relative throughput of Algorithm 1 when compared to the trivial quadratic algorithm on matrices of size
n × 2. A higher throughput implies faster computation.
7

Min, Max
Mean (µ)
Median
µ ±σ

6
5
4
3
2

Trivial

1
400

600

800

1000
1200
1400
Number of rows (n)

1600

1800

2000

Figure 4: Relative throughput of Algorithm 1 when compared to the trivial quadratic algorithm on matrices of size
n × 3. A higher throughput implies faster computation.
7
6
5
4

Min, Max
Mean (µ)
Median
µ ±σ

3
2

Trivial

1
5000

10000
Number of rows (n)

15000

20000

Figure 5: Relative throughput of Algorithm 1 when compared to the trivial quadratic algorithm on matrices of size
n × 4. A higher throughput implies faster computation.

