1-Q-25
Speech Parameter Sequence Modeling with
Latent Trajectory Hidden Markov Model ∗˓ Hirokazu Kameoka
NTT Communication Science Laboratories / The University of Tokyo

1Introduction

The weakness of hidden Markov models (HMMs)
is that they have diﬃculty in modeling and capturing the local dynamics of feature sequences dueto the piecewise stationarity assumption and theconditional independence assumption on feature sequences. Traditionally, in speech recognition systems, this limitation has been circumvented by appending dynamic (delta and delta-delta) components to the feature vectors. HMM-based speechsynthesis systems [1] also use the joint vector of
static and dynamic features as an observed vector
in the training process. In the synthesis process,on the other hand, a sequence of static features isgenerated according to the output probabilities of
the trained HMM given an input sentence by taking
account of the explicit constraint between the staticand dynamic features [2]. Although the HMM-based
speech synthesis framework has many attractive features, one drawback is that the criteria used fortraining and synthesis are inconsistent. While thejoint likelihood of static and dynamic features ismaximized during the training process, the likelihood of only the static features is maximized during
the synthesis process. This implies that the model
parameters are not trained in such a way that thegenerated parameter sequences become optimal. To
address this problem, Zen [3] introduced a variantof HMM called the “trajectory HMM,” which wasobtained by incorporating the explicit relationshipbetween static and dynamic features into the traditional HMM. This has made it possible to provide
a uniﬁed framework for the training and synthesisof speech parameter sequences, however, it causes
diﬃculty as regards parameter inference. Since theconditional independence assumption on the feature vectors is lost, eﬃcient algorithms for training
and decoding regular HMMs such as the Viterbi algorithm and the Forward-Backward algorithm are
no longer applicable to the trajectory HMM. Thus,
some approximations and brute-force methods areusually necessary to obtain training and decoding
algorithms [3, 4].
In this paper, we propose formulating a new model
called the “latent trajectory HMM.” In contrast
with the conventional trajectory HMM, the present
model splits the generative process of an observed
feature sequence into two processes, one for a sequence of the joint vectors of static and dynamicfeatures given HMM states and the other for an observed feature sequence given the sequence of thejoint vectors. By treating the joint vector of staticand dynamic features as a latent variable to be
marginalized out, we obtain a probability density
function of an observed feature sequence with a dif∗ferent form from the likelihood function of the trajectory HMM. As described below, this new formulation naturally allows the combined use of powerful inference techniques such as the expectationmaximization (EM) algorithm, Viterbi algorithm
and Forward-Backward algorithm for training anddecoding, while still retaining the spirit of the original trajectory HMM.
This work is not only directed towards speech synthesis applications but also towards several diﬀerent
applications such as voice conversion and acousticto-articulatory mapping, in which trajectory modeling has proven to be eﬀective [5–7]. Another interesting application we have in mind is audio source
separation. Recently, we proposed methods forsingle- and multi-channel audio source separation
based on factorial HMMs [8–11], where the spectrogram of a mixture signal is modeled as the sum of
the outputs emitted from multiple HMMs, each representing the spectrogram of an underlying source.
One promising way to improve this approach would
be to incorporate the dynamics of source spectra.
This can be accomplished by plugging the present
model into the factorial HMM formulation. The
present formulation will play a key role in making
this possible.

2Trajectory Hidden Markov Model

We start by brieﬂy reviewing the original formulation of the trajectory HMM [3]. Let us use ct todenote a D-dimensional static feature vector and deﬁne the joint vector of ct and its velocity and accelT
2 T T
3Deration components ot := [cTt , Δct , Δ ct ] ∈ R
as the observed vector at time t. We write the sequences of the static features and the observed vecT TTT Ttors as c = [cT
1 , . . . , cT ] and o = [o1 , . . . , oT ] ,respectively. Thus, the dimensions of c and o become DT and 3DT . The relationship between c ando can be described explicitly using a constant 3DT
by DT matrix W aso = W c,(1)where W is a sparse matrix that appends ﬁrst andsecond order time derivatives to the static feature
vector sequence.
Within the traditional HMM framework, a sequence of observed vectors, o, is simply assumed tobe generated from an HMM. Here, if we assume theemission probability density to be a single Gaussian
distribution, the probability density function of ogiven a state sequence s = [s1 , . . . , sT ] and an HMM
parameter set λ = {μ, U , π}, with μ = {μi }1≤i≤I ,U = {Ui }1≤i≤I , and π = {πi,j }1≤i≤I,1≤j≤J , isજࡏτϥδΣΫτϦӅΕϚϧίϑϞσϧʹΑΔԻ੠ಛ௃ྔ‫ྻܥ‬ϞσϦϯά, ُԬ߂࿨ (NTT CS ‫ݚ‬/౦େ)

日本音響学会講演論文集- 263 -2015年9月

given asp(o|s, λ) = N (o; μs , Us ) =TN (ot ; μst , Ust ),t=1(2)

where N (x; μ, Σ) denotes a Gaussian distribution
with mean μ and covariance Σ:N (x; μ, Σ) ∝1|Σ|112T

e− 2 (x−μ) Σ−1(x−μ).

(3)μs and Us denote the mean sequence and a block
diagonal matrix whose diagonal elements are givenby the sequence of the covariance matrices of theemission densities over time:T T
μs = [μTs 1 , . . . , μ sT ] ,Us = diag(Us1 , . . . , UsT ).

(4)
(5)

In HMM-based speech synthesis systems, the parameter set is typically trained by solving the maximum likelihood estimation problem
p(o|s, λ)p(s),(6)λ̂ = argmax log
λs

where p(s) is given by the product of the state transition probabilities. At the synthesis stage, given astate sequence s and with the trained parameter λ,a static feature sequence c is generated according toĉ = argmax p(c|s, λ),c− 12 (cT W T Us−1 W c−2cT W T Us−1 μs )∝e= N (c; c̄s , Vs ).

3(8)
(9)
(10)

By completing the square in the exponent of (9) withrespect to c, we obtain c̄s and Vs asc̄s = (W T Us−1 W )−1 W T Us−1 μs ,Vs = (WTUs−1 W )−1 .

(11)
(12)

Thus, the solution to (7) is c̄s . Geometrically, (10)can be viewed as a cutting plane of the densityp(o|s, λ) at o = W c.
As shown above, the traditional HMM-based
framework uses diﬀerent criteria for training andsynthesis: While p(o|s, λ) is used for training,p(c|s, λ) is used for synthesis. This implies that λ̂
is not necessarily optimal for generating optimal c.
To address this inconsistency between the training
and synthesis criteria, Zen [3] proposed introducing
a framework called the “trajectory HMM”, which
also uses (10) as the training criterion. Instead ofsolving (6), the parameter set λ is trained by solving{λ̂, ŝ} = argmax log p(c|s, λ)p(s),λ,s

Unlike the regular HMM, the conditional independence assumption on observed vectors does not holdin the trajectory HMM: While the regular HMM assumes that each observed vector depends only on thecurrent state, (10) indicates that the observed vector ct at each frame depends on the entire state sequence. This implies that it is diﬃcult to apply theeﬃcient decoding and training algorithms used in
the HMM framework (such as the Viterbi algorithm)
and so some approximations and brute-force methods are necessary to perform training and decoding[3]. Thus, the decoding algorithm is not guaranteed
to ﬁnd the optimal state sequence and the training
algorithm is not guaranteed to converge to a local
optimal solution. Note that this also applies to theminimum generation error (MGE) training framework [4], which uses (10) in which Vs is replaced by
an identity matrix as the training criterion.

(7)

where p(c|s, λ) is deﬁned asp(c|s, λ) ∝ N (W c; μs , Us )Fig. 1 Illustrations that show the diﬀerence between HMM and trajectory HMM.

(13)

3.1 Model
While the HMM is only capable of describing
piecewise stationary sequences of data vectors, thetrajectory HMM is capable of describing continuously varying sequences of data vectors, governedby discrete hidden states (Fig. 1). This feature
is notable in that it can be used to model manykinds of time series data that are continuous in nature. However, the weakness of this model is thatit causes a diﬃculty as regards parameter inference.
To remedy this weakness, we propose introducing a
conceptually similar framework based on a diﬀerent
formulation, which is advantageous in that it alleviates the diﬃculty related to parameter inference.
Instead of treating o as a function of c, we treato as a latent variable that is related to c through asoft constraint o  W c. The relationship o  W c
can be expressed through the conditional distribution p(c|o). For example, we can deﬁne p(c|o) as1Tp(c|o) ∝ exp − (W c − o) Λ(W c − o) , (14)2where Λ is a constant positive deﬁnite matrix thatcan be set arbitrarily. Indeed, this probability density function becomes larger as o approaches W c.
By completing the square in the exponent of (14)with respect to c, we can write p(c|o) asp(c|o) = N (c; mc|o , Λ−1c|o ),where c is treated as the observed data.

日本音響学会講演論文集

Latent Trajectory HMM- 264 -(15)2015年9月

{λ̂, ŝ} = argmax log p(c|s, λ)p(s).

whereλ,smc|o = Ho,(16)T

H = (W ΛW )−1TW Λ,Λc|o = W T ΛW .

(17)
(18)

By using this and p(o|s, λ) deﬁned in (2), we canwrite p(c|s, λ) asp(c|s, λ) = p(c|o)p(o|s, λ)do,(19)
in a diﬀerent way from (10). Geometrically, this canbe viewed as a marginal distribution of the set of the
projected values of o onto the subspace o = W c.
From (2) and (15), the joint likelihood p(c, o|s, λ)can be written asp(c, o|s, λ) = p(c|o)p(o|s, λ)(20)T1c −mc −m,Λx o∝ exp −xx
2 owhere0mx = Λ−1x
U −1 μs

Λx =s,(21)−W T Λ.
H Λc|o H + Us−1

Λc|o−ΛWT(22)Thus, p(x|s, λ) = N (x; mx , Λ−1x ) where x =[cT , oT ]T . By using the blockwise matrix inversion
formula, Λ−1x is given asΣccΣx = Λ−1x = ΣocΣcoΣoo ,(23)whereTΣcc = Λ−1
c|o + HUs H ,(24)Σco = HUs ,(25)TΣoc = Us H ,Σoo = Us .

p(c|s, λ) = N (c; Hμs , Σcc ).

c(29)3.2 Decoding and training algorithms
As with the trajectory HMM framework, thepresent framework uses p(c|s, λ) for feature sequence generation, state decoding and parameter
training in a consistent manner. The problems ofstate decoding and parameter training can be formulated as the following optimization problems:s日本音響学会講演論文集

By leaving only the terms that depend on s and λ,Q(s, λ) can be written ass,λQ(s, λ) = Eo|c,s,λ [log p(o|s, λ)] + log p(s)1(33)= − log |Us | + Tr Us−1 R2
−1T −1− 2μTs Us ō + μs Us μs + log p(s),whereō = μs + Σoc Σ−1cc (c − Hμs ),(34)TR = Σoo − Σoc Σ−1cc Σco + ōō .

(35)Here, the prime mark indicates the values obtainedusing the model parameters updated at the previousiteration. Since Us is a block diagonal matrix, asgiven in (5), (33) can be decomposed into the sumof T individual terms:s,λQ(s, λ) = −T1{log |Ust | + Tr[Us−1Rt ]t
2 t=1−1T−1− 2μTst Ust ōt + μst Ust μst } +⎤⎡R1ō1ō = ⎣ ... ⎦ , R = ⎣ōT∗⎡(28)

Obviously, the solution to this is Hμs .

ŝ = argmax log p(c|s, λ)p(s)Q(s, λ) := Eo|c,s,λ [log p(c, o|s, λ)] + log p(s).
(32)

T


log πst−1 ,st ,t=1

where

We call this model the “latent trajectory HMM.”
With this framework, given a state sequence s anda parameter set λ, c is generated according toĉ = argmax p(c|s, λ).

Since the decoding problem (30) is a subproblem of
the training problem (31), here we only derive analgorithm for solving (31).
By regarding the set consisting of c and o as thecomplete data, this problem can be viewed as an incomplete data problem, which can be dealt with using the Expectation-Maximization (EM) algorithm.
The likelihood of s and λ given the complete data
is given by (20). By taking the conditional expectation of log p(c, o|s, λ) with respect to o given c,s = s and λ = λ , and then adding log p(q), we
obtain an auxiliary function

(26)
(27)

Hence, (19) can be written as

(31)

(30)

..

∗.

⎤⎦.

(36)

RTWith ﬁxed λ, Q(s, λ) can be maximized with respect to s by employing the Viterbi algorithm. Withﬁxed s, Q(s, λ) is maximized with respect to λ when
1[st = i]ōt,(37)μi = t1[st = i]t
1[st = i](ōt − μi )(ōt − μi )Tt,(38)
Ui =1[st = i]t
1[st−1 = i, st = j],(39)
πi,j = t 1[st−1 = i]

- 265 -t

2015年9月

if t 1[st = i] = 0 where i and j denote state indicesand 1[·] denotes an indicator function that takes thevalue 1 if its argument is true and 0 otherwise.
Overall, the parameter training algorithm can besummarized as follows:

800
600
400

Log-likelihood(E-step) Substitute s and λ into s and λ and recompute ō and R using (34) and (35).
(M-step) Update λ using (37)–(39) and ﬁnd s =argmaxs Q(s, λ) using the Viterbi algorithm.

-400-600
ProposedTrajectory HMM-1000
0600

4002008001000

Time (min)Fig. 2Evolutions of the log-likelihoods.

8Experiments

Conclusions

Inspired by the trajectory HMM framework proposed by Zen et al., this paper proposed a proba-

日本音響学会講演論文集

6Frequency (kHz)To conﬁrm the generalization ability of the present
model and the convergence speed of the present
training algorithm, we conducted parameter training experiments using mel-cepstrum sequences of
25 speech data excerpted from the ATR speech
database as experimental data. We chose the parameter training algorithm for the original trajectory HMM [3] as the baseline method, which uses
the method of steepest ascent for updating λ andthe “delayed decision Viterbi algorithm” for updating s. Since the degrees of freedom of the conventional and present models are exactly the same when
the numbers of hidden states are the same, the difference of the log-likelihood scores would reﬂect thediﬀerence in their generalization abilities. For boththe proposed and conventional models, the numbersof hidden states were set at 14.
Fig. 2 shows the evolution of the log-likelihoods
during the parameter training of the proposed andconventional models. As Fig. 2 shows, the presentalgorithm converged faster than the conventional algorithm. This reveals the eﬀectiveness of the combined use of eﬃcient statistical inference techniquessuch as the EM algorithm and the dynamic programming principle by the proposed algorithm. It
is also worth noting that the converged value of thelog-likelihood obtained with the proposed algorithm
was greater than that obtained with the conventional algorithm. This implies the possibility that,compared with the conventional model, the proposed model has a higher ability to ﬁt an arbitrary
set of feature sequences, given that the degrees-offreedom of the two models were the same. Fig. 3shows an example of parameter generation using theproposed model. After training λ using 25 speechdata, a feature sequence ĉ was generated according
to (29) given a state sequence s. The graph shows
the spectrogram constructed using ĉ obtained using
the trained λ and the state sequence s labeled fromthe speech sample. As Fig. 3 shows, the proposedmodel was able to represent the continuously timevarying nature of speech spectrograms reasonably
well, showing that it has a similar property to thetrajectory HMM.

50-200-800

Note that if λ is ﬁxed, the above algorithm reduces to a state decoding algorithm. One reasonable way to initialize s would be to search fors = argmaxs p(o|s, λ) using the Viterbi algorithm.

4200

4

2

00

0.2

0.4

0.6

0.8

1.0

1.2

Time (s)Fig. 3

Example of parameter generation.

bilistic generative model for describing continuously
time-varying sequences of data vectors governed bydiscrete hidden states. The proposed model is advantageous over the conventional trajectory HMM
in that it makes it possible to derive convergenceguaranteed and eﬃcient algorithms for parameter
training and state decoding. Interesting future
work involves incorporating the proposed model into
the factorial HMM formulation to develop a new
method for audio source separation.

Acknowledgement
We thank Dr. Tomoki Toda (NAIST) for fruitful discussions and Mr. Tomohiko Nakamura (University ofTokyo) for his help with conducting the experiments.
This work was supported by JSPS KAKENHI Grant
Numbers 26730100 and 26280060.

References[1] Yoshimura et al., Proc. EUROSPEECH, pp. 2347–2350, 1999.
[2] Tokuda et al., Proc. ICASSP, pp. 1315–1318, 2000.
[3] Zen et al., Computer Speech & Language, vol. 21,no. 1, pp. 153–173, 2007.
[4] Wu & Wang, Proc. ICASSP, pp. 89–92, 2006.
[5] Toda et al., IEEE Trans. ASLP, vol. 15, no. 8,pp. 2222–2235, 2007.
[6] Toda et al., Speech Communication, vol. 50, no. 3,pp. 215–227, 2007.
[7] Zhang & Renals, IEEE Signal Processing Letters,
vol. 15, pp. 245–248, 2008.
[8] Nakano et al., Proc. WASPAA, pp. 325–328, 2011.
[9] Higuchi et al., Proc. Interspeech, pp. 850–854, 2014.
[10] Higuchi & Kameoka, Proc. MLSP, 2014.
[11] Higuchi & Kameoka, Proc. GlobalSIP, 2014.

- 266 -2015年9月