Table Of ContentSUBMITTEDTOIEEETRANSACTIONSONINFORMATIONTHEORY,SEPTEMBER2008 1
Learning Hidden Markov Models using
Non-Negative Matrix Factorization
George Cybenko, Fellow, IEEE, and Valentino Crespi, Member, IEEE
Abstract—TheBaum-Welshalgorithmtogetherwithitsderiva- Collectively,matricesA(k)andΓcompletelydefinetheHMM
tives and variations has been the main technique for learning and we say that a model for the HMM is λ = ({A(k) | 1 ≤
8 Hidden Markov Models (HMM) from observational data. We k ≤M},Γ).
0 present an HMM learning algorithm based on the non-negative
We presentan algorithmfor learning an HMM from single
0 matrixfactorization (NMF)of higherorderMarkovian statistics
2 that is structurally different from the Baum-Welsh and its asso- or multiple observation sequences. The traditional approach
ciated approaches. The described algorithm supports estimation forlearninganHMMistheBaum-WelshAlgorithm[1]which
p
of the number of recurrent states of an HMM and iterates the has been extended in a variety of ways by others [11], [12],
e
non-negative matrix factorization (NMF) algorithm to improve
S [13]. More recently, a different type of HMM approximation
thelearnedHMMparameters.Numericalexamplesareprovided
algorithm has been proposed by Finesso et al. [14], based
4 as well.
on stochastic realization techniques. Our approach is related
2
IndexTerms—HiddenMarkovModels,machinelearning,non-
to the Finesso approach which is significantly different from
negative matrix factorization.
] Baum-Welsh based methods. In particular,our algorithm is in
G
the spirit of the Myhill-Nerode [15] construction for building
L I. INTRODUCTION automata from languages. In that construction, states are
.
s Hidden Markov Models (HMM) have been successfully defined as equivalence classes of pasts which produce the
c used to modelstochastic systems arising in a varietyof appli- same future outputs. In an HMM, the “future” of a state is
[
cationsrangingfrombiologytoengineeringtofinance[1],[2], a distribution over future observations.
1 [3], [4], [5], [6]. Following accepted notation for representing At a conceptual level, our algorithm operates as follows.
v the parameters and structure of HMM’s (see [7], [8], [9], [1], We firstestimatethematrixofanobservationsequence’shigh
6
[10] for example), we will use the following terminologyand order statistics. This matrix has a natural non-negativematrix
8
0 definitions: factorization (NMF) [16] which can be interpreted in terms
4 1) N is the numberof states ofthe Markovchain underly- of the probabilitydistributionof future observationsgiventhe
9. ing the HMM. The state space is S ={S1,...,SN} and currentstateoftheunderlyingMarkovChain.Onceestimated,
0 the system’s state process at time t is denoted by x ; these probabilitydistributionscan be usedto directlyestimate
t
8 2) M is the number of distinct observables or symbols the transition probabilities of the HMM.
0 generatedbythe HMM. Thesetof possibleobservables The estimated HMM parameters can be used, in turn, to
:
v is V ={v ,...,v } and the observationprocessat time compute the NMF matrix factors as well as the underlying
1 M
i t is denoted by y . We denote by yt2 the subprocess higher order correlation matrix from data generated by the
X t t1
y y ...y ; estimated HMM. We present a simple example in which an
r t1 t1+1 t2
a 3) The joint probabilities NMF factorization is exact but does not correspond to any
HMM.Thisis afactthatcan beestablishedbycomparingthe
a (k)=P(x =S ,y =v |x =S );
ij t+1 j t+1 k t i factors computed by the NMF with the factors computed by
are the time-invariant probabilities of transitioning to the estimated HMM parameters. This kind of comparison is
state S at time t+1 and emitting observationv given not possible with other approaches [14].
j k
that at time t the system was in state S . Observation It is important to point out that the optimal non-negative
i
v is emitted duringthe transition fromstate S to state matrix factorization of a positive matrix is known to be NP-
k i
S . We use A(k) = (a (k)) to denote the matrix of Hard in the general case [17], so in practice one computes
j ij
state transition probabilities that emit the same symbol only locally optimal factorizations. As we will show through
v . Note that A = A(k) is the stochastic matrix examples, the repeated iteration of the factorization and tran-
k k
representing the MarkPov chain state process xt. sition probability estimation steps improves the factorizations
4) The initial state distribution, at time t = 1, is given by and overall model estimation. Details are provided below.
Γ = {γ ,...,γ } where γ = P(x = S ) ≥ 0 and
1 N i 1 i
γ =1.
i i A. Preliminaries and Notation
P
G.CybenkoiswiththeThayerSchoolofEngineering,DartmouthCollege, The only inputto our algorithm is an observationsequence
Hanover, NH03755USAe-mail: [email protected].
of length T of the HMM, namely:
V. Crespi is with the Department of Computer Science, California State
University atLosAngeles, LA,90032USAe-mail: [email protected].
Manuscript submittedSeptember 2008 O1:T =O1O2...OT
SUBMITTEDTOIEEETRANSACTIONSONINFORMATIONTHEORY,SEPTEMBER2008 2
where O ∈V is the HMM output at observation time t. distributiononsubstringsoflengths,P(·|S ,s,λ).Itisknown
t i
Wedonotassumethattheobservationtimet=1coincides that the distributions P(·|S ,s,λ) for p+ s ≥ 2N −1 are
i
with the process’ initial state so that the initial distribution of complete characterizations of the ergodic states of the HMM
states is not necessarily governed by Γ. In fact, at present, with respect to the observables of the HMM [20], [14].
ouralgorithmis capableof learningonly the ergodicpartition We now focus attention on substrings that precede state
of an HMM, namely the set of states that are recurrent. occupancyin the HMM’s underlyingMarkov chain. Over the
Consequently, our model of an HMM refers only to the course of a long observation sequence such as O , there
1:T
transition probability component λ = {A(k)} that identifies is some probability, P(S |U,p,λ) that the HMM is in state
k i
this ergodic partition (see [18], [19] for some background on i given that we have just observed the length p substring
this concept). U =v v ...v . These probabilitiescan be computedfrom
j1 j2 jp
Given O , we construct two summary statistics repre- the A(k)’s according to
1:T
sented as matrices Rp,s and Fp,s for positive integers p and π′ p A(j )e
s. Rp,s is simply a histogram of contiguous prefix-suffix P(S |U,p,λ)= r=1 r si0 , (2)
i0 QP(U|p,λ)
combinations whose rows are indexed by observations subse-
quences of length p and columns are indexed by observation whereπisthestationarydistributionoftheunderlyingMarkov
subsequences of length s. chain process and P(U|p,λ)=π′ p A(j )1.
r=1 r
If there are M symbols in the observation alphabet, then Note that formulae (1) and (2) aQre closely related to com-
Rp,s is an Mp by Ms matrix whose (i,j)th entry is the putations arising in the classical Viterbi algorithm [1].
number of times the prefix substring corresponding to i is Let U,V be two strings of observations of length p and s
immediatelyfollowedbythesubstringcorrespondingtoj.The respectively.LetU andV beidentifiedwithintegersuandvas
correspondence between strings and integers is lexicographic already explained before so that P(V|U,λ) = Fp,s. Assume
u,v
inourexamplesbelowalthoughanyothercorrespondencewill V was emitted after time t and U immediately preceded V.
do as well. We call U the prefix string and V the suffix string. Then by
The matrix Fp,s is simply Rp,s normalized to be row applying elementary properties of probability we can write:
stochastic. Specifically, if G = (g ) where g = Rp,s
then Fip,j,s = Rip,,js/gi for gi 6= 0 anid Fip,j,s = i0 forPgij=i,0j. Fup,,vs ∼ P(ytt++1s =V | ytt−p+1 =U,λ)
RowsofRp,s,andcorrespondinglyFp,s,arezeroiftheprefix N
corresponding to the row label is not observed in the data. = P(ytt++1s =V,xt =Sk|ytt−p+1 =U,λ)
Zero rows of these matrices can be deleted reducing the size Xk=1
ofthematriceswithoutaffectingthealgorithmdescribebelow. N
Accordingly, Fp,s is constructed to be row stochastic. = P(V|Sk,s,λ)P(Sk|U,p,λ) .
Entry Fp,s is essentially an estimate of P(V|U) the prob- Xk=1
u,v
ability of observing observation sequence V of length s, Consequently we can express the distribution Fup,,:s ∼
indexed by v, following observation sequence U of length p, P(·|U,λ) as a mixture
indexed by u.
N
Note thatwhile R, F and G haveexponentiallymanyrows Fp,s ∼ P(·|S ,s,λ)P(S |U,p,λ) . (3)
u,: k k
and columns with respect to p and s, the actual number of kX=1
nonzero entries in these matrices are bounded above by T so
If the underlying state process x is ergodic then in the limit
that, stored as sparse matrices, they require no more storage t
asT →∞relation(3) becomesan equalityalmostsurely. As
thantheoriginalobservationsequence.NotethatBaum-Welsh
a result of the above observations, for sufficiently large p and
methods require storing and repeatedly accessing the original
s, the matrix Fp,s has the following properties:
observation sequence.
• rank(Fp,s) ≤ N, where N is the minimal number of
AsimplebutkeyobservationaboutstatesofanHMMisthat
states representing the HMM, λ;
each state of an HMM induces a probability distribution on
• Each row of Fp,s is a convex combination (mixture) of
symbol subsequences of any length s. Specifically, suppose
the N generators, P(·|S ,s,λ), for i=1,2,...,N.;
an HMM, λ, is in state S (having not yet emitted an i
i0 • Letting D be the N × Ms nonnegative matrix whose
observationinthatstate)andconsiderthesymbolsubsequence
rows are the distributions P(·|S ,s,λ), i.e., D =
V =v v ...v . Then k k,:
j1 j2 js P(·|S ,s,λ), for k = 1,2,...,N, we can rewrite (3)
k
P(V|S ,s,λ)=P(yt+s =v v ...v |x =S ) as
i0 t+1 j1 j2 js t i0
is independent of t under the ergodic assumption and can be Fup,,:s ∼[P(S1|U,p,λ)P(S2|U,p,λ) ··· P(SN|U,p,λ)]∗D.
computed from the A(k)’s according to
Consequently, if we let C = (c ) be the Mp × N
u,k
s nonnegative matrix with c = P(S |U,p,λ) we can
P(V|Si0,s,λ)=e′Si0 A(jr)1 , (1) write Fp,s ∼ C ∗ D. Obsue,rkve that Ckand D are both
rY=1 (row) stochastic.
where ei denotes the (0,1)-vector whose only nonzero entry • The factorization depends on the model λ. Moreover
is in position i and 1 = [1 1 ... 1]′. Call this probability factors C and D can be computed directly from λ using
SUBMITTEDTOIEEETRANSACTIONSONINFORMATIONTHEORY,SEPTEMBER2008 3
(1) and (2). Consequently,the size of the smallest model approximateNMFwhichconsistsofdeterminingC and
compatible with the data is equal to prank(Fp,s), the D ofdimensionsMp×N andN×Ms respectivelythat
positive rank of Fp,s. (The positive rank, prank(A), of minimize D (Fp,s||C∗D), where
ID
an m×n nonnegative matrix A is the smallest integer
K
ij
N such that A factors in the productof two nonnegative DID(K||W)= (Kijlog −Kij +Wij)
W
matrices of dimensions m×N and N ×n respectively.) Xij ij
Itisknownthatrank(A)≤prank(A)≤min{m,n}and is the I-divergence function [7] (observe that if
that the computation of prank(A) is NP-hard [17]. So it 1′K1 = 1′W1 = 1 then D (K||W) =
ID
would appear that in general it is NP-hard to estimate K logK /W sotheI-divergencefunctionisa
N given Fp,s even in ideal conditions (T → ∞) since gPeni,ejraliiz,jationoif,jthe Ki,jullback-Leiblerdistance between
rank(Fp,s)≤N.However,itisnotobvioushowdifficult probabilitydistributions).Thisoptimizationproblemcan
it is to estimate when prank(Fp,s) < rank(Fp,s) in be solved through iterative methods [16], [22] that re-
the case Fp,s was built from a typical realization of an quireinitialmatricesC ,D andcanonlybeguaranteed
0 0
HMM. We discuss an example at the end of this paper toconvergetolocaloptima.Afterexecutingthisstep,we
that illustrates the open problems and challenges. have a locally optimal estimate of the true distributions
To summarize this discussion, note that the matrix Fp,s is P(·|S ,s,λ).
i
based on the distribution of length p prefixes and corre- 4) Estimate matrices A(k), k = 1,2,...,N, from D. Let
sponding length s suffixes and completely characterizes an usconsiderA(1)=(a (1)),theothermatricesareesti-
i,j
HMM providing p+s ≥ 2N −1. Its positive rank is equal matedinasimilarmanner.LetV(s−1) =v v ···v
j1 j2 js−1
to the minimal number of states in the underlying Markov be a generic sequence of s−1 observations. Then by
chain. Moreover,an appropriatelyconstructedfactorizationof marginalization we can write
Fp,s exposesthe state transition and emission probabilitiesof
M
the HMM. The algorithm presented below extracts the state P(V(s−1)|S ,s−1,λ)= P(V(s−1)v |S ,s,λ) .
i k i
transition matrices, {A(k)} from this factorization. In turn,
k Xk=1
as shown above, the A(k)’s can be used to construct the
probability distributions over suffixes that generate Fp,s and Consequently,the conditionaldistributionsoversuffixes
of length s − 1, P(·|S ,s − 1,λ), can be estimated
so can be used to compute a new factorization. This iteration i
from D by adding columns of D appropriately. Let
is essentially the basis for our algorithm.
H be the matrix thus obtained from D so that H ≈
In the machine learning context, we have access only to a i,:
P(·|S ,s−1,λ). Those conditionaldistributions satisfy
finite amount of observation data (T bounded). Consequently i
rank(Fp,s) will be generally higher than N. This requires a the following equality for any V(s−1):
decision about the HMM’s order, N, not unlike that arising N
in principal component analysis (PCA) [21] to estimate the P(v1V(s−1)|Si,s,λ)= ai,j(1)P(V(s−1)|Sj,s−1,λ).
number of components. Xj=1
II. THE ALGORITHM Therefore P(v1 ·|Si,s,λ) = Nj=1ai,j(1)P(·|Sj,s−
Based on the above discussion, our algorithm is outlined 1,λ) so we can obtain the unkPnown values ai,j(1) by
solving the following systems of linear equations:
below.Numericalexampleswithdiscussionsfollowtheformal
description. Di,1:Ms−1 =Ai,:(1)∗H, i=1,2,...,N
1) Compute Fp,s and G from the input observation data,
whereA (1)=[a (1)a (1)···a (1)].Compactly
O , defined above. i,: i,1 i,2 i,N
1:T
2) Estimatethenumberofstates,N,byanalyzingtheeither D:,1:Ms−1 = A(1) ∗H. As in step 2, because of the
Fp,s ordiag(G)∗Fp,s, bothcomputedin Step1. Inthe finiteness of T and working with bounded arithmetic
cases in which prank(Fp,s) = rank(Fp,s) (e.g. when precision we need to content ourselves with a solution
that minimizes some distance (for example, the L
rank(F)≤2) one typicalway to obtain this estimate is 1
to compute the SVD (singular value decomposition) of norm)between Di,1:Ms−1 and Ai,:(1)∗H, for all i. We
have formulated these problems as linear programming
theaforementionedmatricesandthenobservetherateof
problems using the L norm.
decrease of the singular values. For T sufficiently large 1
a significant gap between the Nth and the (N +1)th 5) OutputestimatedHMMλ′ ={A(k)|k =1,2,...,N}.
largest singular value becomes appreciable. Note that This algorithm can be iterated using the estimated λ′ and
since prank ≥ rank, an estimate based on the singular formula (1) to compute new matrices C0′ and D0′, and then
values is a lower bound for the order of the HMM. restartingfromstep3abovewithmatricesC′ andD′ asinitial
0 0
3) Estimate distributions P(·|S ,s,λ), for i=1,2,...,N. factors in the approximate NMF. In particular:
i
This step is achieved through the Nonnegative Matrix 6) Compute D′ =[P(j|S ,s,λ′)] using formula (1).
0 i i,j
Factorization(NMF) of Fp,s. This yieldsFp,s ≈C∗D 7) ComputeC′ bysolvingthelinearprogrammingproblem
0
with D ≈P(·|S ,s,λ) as observed before. Fp,s =C′ ∗D′, for a row stochastic C′.
i,: i 0 0 0
Note that because of the finiteness of T in general 8) Set C :=C′ and D :=D′.
0 0 0 0
prank(Fp,s) > N. So it is necessary to solve the 9) goto 3).
SUBMITTEDTOIEEETRANSACTIONSONINFORMATIONTHEORY,SEPTEMBER2008 4
Another possibility for step 7) above is to compute C′ 2) EstimateN =prank(F).AnalyzesingularvaluesofF:
0
using formulae (2) and (3) and then use the resulting C′ and
D′ as initial guesses for the NMF algorithm. We have0tried 0.88 0.48 0.033 0.011 .
0
thisvariantbutitdoesnotproducesignificantlydifferentfinal (cid:2) (cid:3)
This suggests rank(F)=prank(F)=2.
results.
3) Estimate distributions P(·|S ,3,λ ) and P(·|S ,3,λ )
1 1 2 1
by solving argmin D (F||C∗D):
C,D ID
III. NUMERICAL EXAMPLES
0.02 0.98
WecallanHMM“Deterministic”(DHMM)ifforeachstate
1 0
there exists at most one outgoing transition labeled with the C = ,
0 1
sameobservable.WedemonstrateourmethodonaDHMM,on
0.34 0.66
an HMM that can be transformed into an equivalent DHMM
and also on an HMM for which such a transformation does 0 0 0 0.0 0.25 0.24 0.0 0.5
D =
not exist. We finally discuss an example that illustrates the (cid:20) 0.13 0.13 0 0.25 0 0 0.25 0.24 (cid:21)
situation when rank < prank.
4) Estimate matrices A(0) and A(1):
Itisimportanttonotethatthesignificantmetricforlearning
anHMMisnottheextenttowhichthetransitionprobabilities 2.2e−18 0 0.0077 0.99
are accurately learned but the extent to which the observation A˜(0)= ,A˜(1)= .
(cid:20) 6.9e−18 0.51 (cid:21) (cid:20) 0.49 0 (cid:21)
statistics are learned. This is a consequence of the fact that
After a second iteration of the algorithm the reconstructed
HMM’s with different transition probabilities and different
matrices become:
numbersofstatescanproduceobservationssequenceswiththe
samestatistics sothatlearningaspecific transitionprobability 0 0 0.0077 0.99
Aˆ(0)= Aˆ(1)= .
characterization is not a well-posed problem. (cid:20) 5.6e−17 0.51 (cid:21) (cid:20) 0.49 0 (cid:21)
In our examples we measure the accuracy of our estimates
Thereconstructedmodelisessentiallyidenticaltotheoriginal
by computing the I-divergence rate of the finite dimensional
oneexceptforstatereordering.Thisresultisparticularlystrik-
distributions associated with the observation process of the
ing when compared with existing techniques. For example,
original model from those associated with the observation
Crutchfieldetal[23]demonstratedtheirCausal-StateSplitting
process of the estimated model. Formally, each HMM λ
Reconstruction(CSSR)ǫ-machinereconstructionalgorithmon
induces a family of finite dimensional distributions
the same Even Process. Their result is essentially equivalent
N to ours in accuracy but it required the processing of 104
Pn(y1n)= πiP(y1n|x1 =Si,λ) observations instead of the 103 that we use.
Xi=1
on sequences of observations of length n, where π is the
B. An HMM that has an equivalent DHMM
stationary distribution of the underlying state process. Let λ
andλ′betwoHMM’swithPnandQntheirrespectiveinduced Consider the model λ2 =({A(0),A(1)},Γ=[0 1]) with
finite dimensional distributions. The I-divergence rate of λ
0.67 0.33 0 0
from λ′ is defined as A(0)= and A(1)= .
(cid:20) 0 0 (cid:21) (cid:20) 1 0 (cid:21)
1
DID = lim DID(Pn||Qn) We simulated this process and produced a sequence of T =
n→∞n
10000 observations. Then we ran our algorithm with p = 2
when the limit exists.
and s=3:
1) Build F2,3 from data O:
A. A DHMM Example
0.3 0.15 0.22 0 0.22 0.11 0 0
Consider the stochastic process described by model λ =
1 0.45 0.22 0.33 0 0 0 0 0
({A(0),A(1)},Γ=[0 1]) with F =
0.29 0.16 0.22 0 0.23 0.1 0 0
0.5 0 0 0.5 0 0 0 0 0 0 0 0
A(0)= and A(1)= .
(cid:20) 0 0 (cid:21) (cid:20) 1 0 (cid:21)
2) Estimate N. Analyze singular values of F:
Thisissometimesreferredtoasthe“EvenProcess”[23],[24].
0.86 0.24 0.012 0
We simulated this process and produced a sequence of T =
1000observations.Thenweranouralgorithmwithp=2and (cid:2) (cid:3)
to estimate N. This suggests again N =2.
s=3:
3) Estimate distributions P(·|S ,3,λ ) and P(·|S ,3,λ ).
1 2 2 2
1) Build F2,3 from data O: Solve argmin D (F||C∗D):
C,D ID
0.14 0.13 0 0.26 0 0 0.22 0.26
0.6 0.4
F = 0 0 0 0 0.25 0.24 0 0.5 1 0
0.13 0.14 0 0.23 0 0 0.26 0.24 C = ,
0.59 0.41
0.08 0.07 0 0.17 0.08 0.08 0.17 0.33
0 0
SUBMITTEDTOIEEETRANSACTIONSONINFORMATIONTHEORY,SEPTEMBER2008 5
0.45 0.22 0.33 0 0 0.0 0 0
D =(cid:20) 0.07 0.05 0.06 0 0.55 0.27 0 0 (cid:21) 2.5x 10−6 Divergence rate of the original HMM from the estimated HMM
4) Reconstruct matrices A(0) and A(1):
2
0.6 0.4 0 3.2e−21
A˜(0)=(cid:20) 0.1 0.072 (cid:21),A˜(1)=(cid:20) 0.83 0 (cid:21) . Q)/nn1.5
P||n
After a second iteration of the algorithm the reconstructed D(ID 1
model becomes λˆ :
2
0.5
0.6 0.4 0 4.2e−21
Aˆ(0)= ,Aˆ(1)= .
(cid:20) 0.1 0.072 (cid:21) (cid:20) 0.83 0 (cid:21) 0
0 5 10 15
Length of observation sequences
These computed transition probabilities are different enough
from the transition probabilities of the original HMM used 10x 10−3 Divergence rate of the original HMM from the estimated HMM
to generate the data but the statistics of the observation
9
sequencesare very close. Figure 1 shows the accuracyof this
8
estimateintermsoftheI-divergencerateoftheoriginalmodel
from the estimated one. We computed DID(Pn||Qn)/n for Q)/nn 7
n = 1,2,...,15, with Pn (Qn) being the finite dimensional P||n 6
probability distributions over sequences of observations of D(ID 5
length n, emitted by model λ (λˆ ) in stationary conditions.
2 2
4
We can observe that this quantity, the divergence rate of
P from Q , stabilizes to a very small value (smaller than 3
n n
2.5·10−6) as expected. 2
0 5 10 15
In fact, this example is equivalent to a DHMM model as Length of observation sequences
the reader can readily check independently.
Fig. 1. Accuracy of HMM’s λˆ2 (top) and λˆ3 (bottom). Here Pn is the
distribution over sequences of observations of fixed length n induced by
theoriginal modelwhereas Qn referstotheestimatedmodel.Thesequence
C. An HMM that has no equivalent finite state DHMM
DID(Pn||Qn)/niscalculated forincreasing values ofn.
Consider the model λ =({A(0),A(1)},Γ=[1 0 0]) with
3
0.5 0.5 0 0 0 0 where ⊗ is the Kronecker product. We can verify that
A(0)= 0 0.5 0 and A(1)= 0 0 0.5 . rank(F) = 3 whereas prank(F) = 4 [22]. Moreover
0.5 0.5 0 0 0 0 F =CD exactly with
We simulated this process and produced a sequence of T = 0.5 0 0.5 0
10000 observations. Then we ran our algorithm with p = 4 0.5 0.5 0 0
C =
and s = 5. After a second iteration of the algorithm the 0 0 0.5 0.5
reconstructed model becomes λˆ3: 0 0.5 0 0.5
0.29 0.43 0.28 0.0 0 0 and
Aˆ(0)= 0 0.29 0.18 ,Aˆ(1)= 0.44 0 0.081 .
1 0 0 0
0.29 0.71 0 0.0 0 0
1 0 0 1 0
D = ⊗[1 1 1 1 1 1 1 1] .
As before, Figure 1 (bottom) shows the accuracy of this 8 0 1 0 0
estimateintermsoftheI-divergencerateoftheoriginalmodel 0 0 0 1
from the estimated one.
Assume that F was obtained from a typical sequence of
Observe that this HMM cannot be transformed into an
observations emitted by an HMM λ with 4 states so that
equivalent deterministic HMM [25].
F2,5 =F =CD. Then it must be that C =[P(S |i,2,λ)]
j i,j
and D = [P(k|S ,5,λ)] . Consider the following model
j j,k
D. Discussion of Rank vs Prank λ={A(1),A(2)} with
We first provide an example of a stochastic matrix whose 0.5 0 0.5 0 0 0 0 0
prank differsfromits rankbutthatmatrix doesnotrepresent 0 0 0 0 0.5 0 0.5 0
A(0)= ,A(1)= .
the statistics of any HMM. 0 0.5 0 0.5 0 0 0 0
0 0 0 0 0 0.5 0 0.5
1 1 0 0
1 1 0 1 0 One can verify that λ is the only four-state model such that
F = · ⊗[1 1 1 1 1 1 1 1] ,
16 0 1 0 1 D = [P(k|Sj,5,λ)]j,k. In fact observe that the system of
0 0 1 1 equations defining λ in stage 4 of the algorithm admits, in
SUBMITTEDTOIEEETRANSACTIONSONINFORMATIONTHEORY,SEPTEMBER2008 6
this case, only one solution. Nevertheless, using formula (2): for learning HMM’s and has the advantage of a estimating
the HMM order from spectral properties of the high order
1 1 1 1
[P(S |i,2,λ)] =(1/4)∗ ⊗ 6=C . correlation statistics of the observation sequence. The algo-
j i,j (cid:20) 1 1 (cid:21) (cid:20) 1 1 (cid:21)
rithm effectively compresses data by summarizing it into a
Consequently no HMM can generate F. statistical matrix. Options for recursively computing the steps
1) An example of prank > rank for an exact HMM model: ofthealgorithmtoachieveon-linealgorithmswillbeexplored.
The following four-state model λ is an example of an HMM Additionally, sparse matrix algorithms can be explored for
whose induced F2,5 matrix has rank 3 but positive rank space and time efficiency when the underlying matrices are
prank =4: large and sparse.
0.5 0 0 0 0 0.5 0 0
0 0 0 0 0 0 1 0 ACKNOWLEDGMENTS
A(0)= , A(1)= .
0 0 0 0 0 0 0 1 This work results from research programs at the Thayer
0 0 0 0 1 0 0 0 SchoolofEngineering,DartmouthCollege,andattheDepart-
ment of Computer Science at the California State University,
To verify the claim we computed factors C = [P(S |i,λ)]
j i,j
Los Angeles, supported by the AFOSR Grant FA9550-07-
and D = [P(j|S ,5,λ)] , for N = 4, using formulae (1)
i i,j
1-0421. The Dartmouth effort was also supported by U.S.
and (2) and then obtained F2,5 = C ∗D. Then we verified
DepartmentofHomelandSecurityGrantAwardNumber2006-
numerically that rank(F2,5)=3. Finally, we applied Lemma
CS-001- 000001 and DARPA Contract HR001-06-1-0033.
2.4 in [22] to confirm that prank(F2,5) = 4. We also
The authors also sincerely thank Dr. Robert Savell for
verified the character of this model by directly applying our
algorithm to it in order to obtain F2,5 empirically (for T = pointing out important related work to us.
10000). An analysis of the singular values of F2,5, namely
[0.8530 0.4825 0.1799 0.0114], demonstrates the difficulty REFERENCES
of this case. The fourth singular value is nonzero due to
[1] L. R. Rabiner, “A Tutorial on Hidden Markov Models and Selected
the finiteness of T. Consequently it is difficult to determine Applications inSpeechRecognition,” Proceeding oftheIEEE,vol.77,
whether N =3 or N =4. no.2,pp.257–286,1989.
[2] D.X.Song, D.Wagner, and X.Tian, “TimingAnalysis ofKeystrokes
andTimingAttacks onSSH,”inUSENIXSecurity Symposium,2001.
IV. OPENQUESTIONS AND FUTURE WORK [3] A. Apostolico and G. Bejerano, “Optimal Amnesic Probabilistic Au-
tomata or How to Learn and Classify Proteins in Linear Time and
A crucial issue is the estimation of N, the size of the Space,”JournalofComputationalBiology,vol.7,no.3/4,pp.381–393,
smallest HMM that generates the stream of data. Under ideal 2000.
[4] D. Ron, Y. Singer, and N. Tishby, “Learning Probabilistic Automata
conditions, (T → ∞), we have seen that N = prank(Fp,s).
with variable memory length,” in Proceedings of the workshop on
However,filtering out“noise” from the empiricalmatrix Fp,s computational learningtheory, 1994.
inordertohaveanaccurateestimateofthepositiverankisan [5] F. Perronnin, J.-L. Dugelay, and K. Rose, “A Probabilistic Model
of Face Mapping with Local Transformations and its Application to
open challenge. Observed that a spectral analysis of F may,
PersonRecognition,” PatternAnalysis andMachine Intelligence, IEEE
in general, produce only a lower bound to N. Transactions on,vol.27,no.7,pp.1157–1171, July2005.
A second important issue in our methodologyconcerns the [6] A.Krogh,M.Brown,I.S.Mian,K.Sjolander,andD.Haussler,“Hidden
Markov Models in Computational Biology: Applications to Protein
computation of the approximate NMF. Existing methods are
Modeling,” UCSC,Tech.Rep.UCSC-CRL-93-32, 1993.
suboptimaldue to the presence of local optima. This problem [7] L. Finesso and P. Spreij, “Nonnegative Matrix Factorixation and I-
affectsthe accuracyof the producedestimate ateach iteration divergence Alternating Minimization,” Linear Algebra and its Appli-
cations, vol.416,pp.270–287,2006.
of our algorithm. Consequently it is important to investigate
[8] E.Vidal,F.Thollard,C.delaHiguera,F.Casacuberta,andR.Carrasco,
convergencepropertieswhenstages3−5of thealgorithmare “Probabilistic Finite-State Machines-Part i,” PAMI, vol. 27, no. 7, pp.
iteratedwithnewinitialfactorsC′,D′ toseedtheapproximate 1013–1025, July2005.
NMF, usingC′ andD′ ascompu0ted0accordingto steps6−8, [9] ——,“ProbabilisticFinite-StateMachines-Partii,”PAMI,vol.27,no.7,
0 0 pp.1026–1039, July2005.
from model λ′ that was estimated in the preceding step. [10] Y. Ephraim, “Hidden Markov Processes,” IEEETransactions onInfor-
A third question concerns with properties of Fp,s as s → mationTheory,vol.48,no.6,pp.1518–1569, 2002.
[11] E.B.Fox,E.B.Sudderth, M.I.Jordan,andA.S.Willsky,“AnHDP-
∞.Inotherwords,cantheAsymptoticEquipartitionProperty HMMforSystemswithState Persistence,” inProceedings ofthe25th
beappliedtodistributionsP(·|Si,s,λ)sothatthedistribution InternationalConferenceonMachineLearning,Helsinky,Finland,2008.
on the “typical” finite suffixes is uniform and the rest of the [12] D.BleiandM.Jordan,“VariationalMethodsfortheDirichletProcess,”
inProceedingsofthe21stInternationalConferenceonMachineLearn-
distribution is zero?
ing,2004.
[13] Y. Teh, M. Jordan, M. Beal, and D. Blei, “Hierarchical Dirichlet Pro-
cesses,” Technical Report 653, Department Of Statistics, UC Berkeley,
V. CONCLUSION
2003.Advances inNeuralInformation ProcessingSystems, Cambridge
We have presented a new algorithm for learning an HMM (MA),vol.17,2005.
[14] L. Finesso and P. Spreij, “Approximation of Stationary Processes by
from observations of the HMM’s output. The algorithm is
HiddenMarkovModels,”arXiv:math/0606591v2, February2008.
structurally different from traditional Baum-Welsh based ap- [15] J. E. Hopcroft and J. D. Ullman, Introduction to Automata Theory,
proaches[1],[11],[12],[13].Itisrelatedtobutdifferentfrom Languages and Computation. Addison-Wesley Publishing, Reading
Massachusetts, 1979.
recent approaches in stochastic systems realization [14]. We
[16] D. D. Lee and H. S. Seung, “Algorithms for Non-negative Matrix
believethismethodopensanewlineofalgorithmdevelopment Factorization,” inNIPS,2000,pp.556–562.
SUBMITTEDTOIEEETRANSACTIONSONINFORMATIONTHEORY,SEPTEMBER2008 7
[17] S. Vavasis, “On the Complexity of Nonnegative Matrix Factorization,”
arXiv:0708.4149v2, February2008.
[18] A. Berman and R. J. Plemmons, Nonnegative Matrices in the Mathe-
matical Sciences. SIAM,1994.
[19] P. Billingsley, Probability and Measure, 3rd Edition. Wiley-
Interscience, 1995.
[20] J. Carlyle, Stochastic Finite-State system theory, ser. Systems Theory.
McGrawHill, NewYork,1969,ch.10.
[21] P. Comon, “Independent component analysis: a new concept?” Signal
Processing,vol.36,pp.287–314,1994.
[22] J. Cohen and U. Rothblum, “Nonnegative Ranks, Decompositions and
Factorizations ofNonnegative Matrices,”LinearAlgebraanditsAppli-
cations, vol.190,pp.149–168,1993.
[23] C. R. Shalizi, K. L. Shalizi, and J. P. Crutchfield, “Pattern discovery
in time series, part i: Theory, algorithm, analysis, and convergence,”
JournalofMachine LearningResearch,(inpreparation), 2002.
[24] B.Weiss,“Subshiftsoffinitetypeandsophicsystems,”Monatsheftefur
Mathematik, vol.77,pp.462–474, 1973.
[25] J. P. Crutchfield and C. R. Shalizi, “Thermodynamic Depth of Causal
States: Objective Complexity via Minimal Representations,” Physical
ReviewE,vol.59,pp.275–283,1999.
GeorgeCybenkoistheDorothyandWalterGramm
Professor of Engineering at the Thayer School of
Engineering at Dartmouth College. Prior to join-
ing Dartmouth, he was Professor of Electrical and
PLACE Computer Engineering at the University of Illinois
PHOTO at Urbana-Champaign. His current research inter-
HERE ests are in machine learning, signal processing and
computersecurity.CybenkowasfoundingEditor-in-
ChiefofIEEEComputinginScienceandEngineer-
ing and IEEE Security & Privacy. He is presently
FirstVice-President oftheIEEEComputerSociety.
Heearned his B.Sc. (Toronto)and Ph.D(Princeton) degrees inmathematics
andisaFellow oftheIEEE.
Valentino Crespi received his Laurea Degree and
hisPh.D.DegreeinComputerSciencefromtheUni-
versityofMilan, Italy, inJuly1992andJuly1997,
respectively. FromSeptember 1998toAugust2000
PLACE hewasanAssistantProfessorofComputerScience
PHOTO attheEasternMediterraneanUniversity,Famagusta,
HERE North Cyprus and from September 2000to August
2003 he worked at Dartmouth College, Hanover,
NH,asaResearchFaculty.SinceSeptember2003he
hasbeenanAssistantProfessorofComputerScience
attheCalifornia StateUniversity, LosAngeles.His
research interests include Distributed Computing, Tracking Systems, UAV
Surveillance, Sensor Networks, Information and Communication Theory,
Complexity Theory and Combinatorial Optimization. At Dartmouth College
he developed the TASK project and participated to the Process Query
Systemsproject,directed byProf.GeorgeCybenko.AtCSULAhehasbeen
teaching lower division, upper division and master courses on Algorithms,
Data Structures, Java Programming, Compilers, Theory ofComputation and
Computational Learning of Languages and Stochastic Processes. During his
professionalactivityDrCrespihaspublishedanumberofpapersinprestigious
journals and conferences of Applied Mathematics, Computer Science and
Engineering. Moreover Dr Crespi is currently a member of the ACM and
oftheIEEE.
