Title: RECURRENT HIDDEN SEMI-MARKOV MODEL

Abstract: Segmentation and labeling of high dimensional time series data has wide applications in behavior understanding and medical diagnosis. Due to the difficulty of obtaining a large amount the label information, realizing this objective in an unsupervised way is highly desirable. Hidden Semi-Markov Model (HSMM) is a classical tool for this problem. However, existing HSMM and its variants typically make strong generative assumptions on the observations within each segment, thus their ability to capture the nonlinear and complex dynamics within each segment is limited. To address this limitation, we propose to incorporate the Recurrent Neural Network (RNN) as the generative process of each segment, resulting the Recurrent HSMM (R-HSMM). To accelerate the inference while preserving accuracy, we designed a structure encoding function to mimic the exact inference. By generalizing the penalty method to distribution space, we are able to train the model and the encoding function simultaneously. We also demonstrate that the R-HSMM significantly outperforms the previous state-of-the-art on both the synthetic and real-world datasets.

1. INTRODUCTION:

Segmentation and labeling of time series data is an important problem in machine learning and signal processing. Given a sequence of observations {x 1 , x 2 , . . . , x T }, we want to divide the T observations into several segments and label each segment simultaneously, where each segment consists of consecutive observations. The supervised sequence segmentation or labeling techniques have been well studied in recent decades (Sutskever et al., 2014; Kong et al., 2015; Chen et al., 2015) . However, for complicated signals, like human activity sensor data, accurately annotating the segmentation boundary or the activity type would be prohibitive. Therefore, it is urgent to develop unsupervised algorithms that can jointly learn segmentation and labeling information directly from the data without supervisions. Figure 1 provides an illustration which we are focus on.

The Hidden Semi-Markov Model (HSMM) (Murphy, 2002 ) is a powerful model for such task. It eliminates the implicit geometric duration distribution assumptions in HMM (Yu, 2010) , thus allows the state to transit in a non-Markovian way. Most of the HSMM variants make strong parametric assumptions on the observation model (Rabiner, 1989; Johnson & Willsky, 2013; Yu, 2010) . This makes the learning and inference simple, but ignores the nonlinear and long-range dependency within a segment. Take the human activity signals as an example. The movements a person performs at a certain time step would rely heavily on the previous movements, like the interleaving actions of left hand and right hand in swimming, or more complicated dependency like shooting after jumping in playing basketball. Some models have been proposed to tackle this problem (Ghahramani & Hinton, 2000; Fox et al., 2009; Linderman et al., 2016) , but are limited in linear case.

Since people have justified RNN's ability in modeling nonlinear and complicated dependencies (Sutskever et al., 2014; Du et al., 2016) , we introduce the recurrent neural emission model into HSMM for capturing various dependencies within each segment to address such issue. However, the flexibility of recurrent neural model comes with prices: it makes the exact Expectation-Maximization (EM) algorithm computationally too expensive.

To speed up the learning and inference, we exploit the variational encoder (VAE) framework (Kingma & Welling, 2013) . Specifically, we propose to use bidirectional RNN (bi-RNN) encoder. Such The first two rows are our algorithms which almost exact locate every segment.

architecture will mimic the forward-backward algorithm, and hence is expected to capture similar information as in exact posterior calculation.

It should be emphasized that due to the discrete nature of the latent variables in our model, the algorithm proposed in Kingma & Welling (2013) and its extension on time-series models (Gao et al., 2016; Krishnan et al., 2015) are not directly applicable. There are plenty of work proposed based on stochastic neuron (Tang & Salakhutdinov, 2013; Bengio et al., 2013; Mnih & Gregor, 2014; Raiko et al., 2014; Gu et al., 2015; Chung et al., 2016) to remedy such issue. However, none of these off-the-shelf methods are easy to achieve good performance according to our experiment: the hundreds or thousands layers of stochastic neuron (which is equal to the length of sequence), together with the switching generative RNN, make the encoding function very sensitive, and thus, extremely difficult to train fully on unsupervised setting. We propose a solution, stochastic distributional penalty method, which introduces auxiliary distributions to separate the decoding R-HSMM and encoding bi-RNN in training procedure, and thus, reduces the learning difficulty for each component. This novel algorithm is general enough and can be applied to other VAE with discrete latent variables, which can be of independent interest. We emphasize that the proposed algorithm is maximizing exact the nagative Helmholtz variational free energy. It is different from Johnson et al. (2016) in which a lower bound of the variational free energy is proposed as the surrogate to be maximized for convenience.

We experimentally justified our algorithm on the synthetic datasets and three real-world datasets, namely the segmentation tasks for human activity, fruit fly behavior and heart sound records. The R-HSMM with Viterbi exact inference significantly outperforms basic HSMM and its variants, demonstrating the generative model is indeed flexible. Moreover, the trained bi-RNN encoder also achieve similar state-of-the-art performances to the exact inference, but with 400 times faster inference speed, showing the proposed structured encoding function is able to mimic the exact inference efficiently.

2. MODEL ARCHITECTURE:

Given a sequence x = [x 1 , x 2 , . . . , x |x| ], where x t ∈ R m is an m dimensional observation at time t, our goal is to divide the sequence into meaningful segments. Thus, each observation x t will have the corresponding label z t ∈ Z, where Z = {1, 2, . . . , K} is a finite discrete label set and K is predefined. The label sequence z = [z 1 , z 2 , . . . , z |x| ] should have the same length of x.

Besides labels, HSMM will associate each position t with additional variable d t ∈ D = {1, 2, . . . , D}, where d t is known as duration variable and D is the maximum possible duration. The duration variable can control the number of steps the current hidden state will remain. We use d to denote the duration sequence. We also use notation x t1:t2 to denote the substring [x t1 , x t1+1 , . . . , x t2 ] of x. Without ambiguity, we use z as a segment label, and d as the duration.

B A π d 1 d 2 d 3 d 4 d 5 z 1 z 2 z 3 z 4 z 5 x 1 x 2 x 3 x 4 x 5 θ (a) Classical Hidden Semi-Markov Model. B d 1 d 2 d 3 d 4 d 5 z 1 z 2 z 3 z 4 z 5 A π θ rnn x 1 x 2 x 3 x 4 x 5 h 1 h 2 h 3 h 4 h 5 (b) Recurrent Hidden Semi-Markov Model.

Figure 2 : Graphical models of HSMM and R-HSMM. Different from classical HSMM, the R-HSMM has two-level emission structure with recurrent dependency.

In this paper, we focus on one of the variants of HSMM, namely the explicit duration HMM (EDHMM) (Rabiner, 1989) , and use Decreasing Count Variables (Chiappa, 2014) for the notation.

Explicit Duration Hidden Markov Model. Similar to HMM, this model treats the pair of (z, d) as 'macro hidden state'. The probability of initial macro state is defined as P (z, d) = P (z)P (d|z). We use the notation π z P (z) and P (d|z) B z,d to parametrize the initial probability and duration probability, respectively. A i,j P (z t = i|z t-1 = j, d t-1 = 1) is the state transition probability on the segment boundary. Here π ∈ R K is in K-dimensional simplex. For each hidden state z, the corresponding rows B z,: and A z,: are also in probability simplex. Here we assume the multinomial distribution for P (d|z).

In EDHMM, the transition probability of macro hidden state P (z t , d t |z t-1 , d t-1 ) is decomposed by P (z t |z t-1 , d t-1 )P (d t |z t , d t-1 ) and thus can be defined as:

P (z t |z t-1 , d t-1 ) = A zt-1,zt if d t-1 = 1 I (zt=zt-1) if d t-1 > 1 ; P (d t |z t , d t-1 ) = B zt,dt if d t-1 = 1 I (dt=dt-1-1) if d t-1 > 1

.

(1) The graphical model is shown in Figure 2a .

Recurrent Hidden Semi-Markov Model. For the simplicity of explanation, we focus our algorithm on the single sequence first. It is straightforward to apply the algorithm for dataset that has multiple sequences. Given the parameters {π, A, B}, the log-likelihood of a single observation sequence x can be written as below,

EQUATION

where P (x|z, d) is the emission probability. To define P (x|z, d), we further denote the sequence variable s = [s 1 , s 2 , . . . , s |s| ] to be the switching position (or the beginning) of segments. Thus s 1 = 1, and s i = s i-1 + d si-1 and |s| is the number of segments. Traditional HSMM assumes P (x|z, d) = |x| t=1 P (x t |z t ), which ignores the dependency and some degree of dynamics exhibited in each segment. While in this paper, we use RNN as the generative model to capture the dependent emission probability. Specifically, for the i-th segment starting from position s i , the corresponding generative process is

EQUATION

where we assume that the dependency of history before time step j can be captured by a vector h j ∈ R h . As in RNN, we use a recurrent equation to formulate the history vector,

EQUATION

Finally, in this model, P (x|z, d) = |s| i=1 P (x si:si+ds i -1 |z si , d si ) is computed by the product of generative probabilities for each segment. In Eq. 4, W ∈ R m×h is a weight matrix capturing the last observation x t-1 , and V ∈ R h×h is for the propagation of history h t-1 . The b is a bias term. The superscript z si indexes the RNN we used for the corresponding segment. The segments with different labels are generated using different RNNs. So we should maintain K RNNs. σ(•) is a nonlinear activation function. We use tanh in our experiments.

At the time step t, we assume a diagonal multivariate gaussian distribution over the conditional likelihood, where the mean and covariance matrix are the output of RNN, i.e.,

EQUATION

The matrices U µ , U Σ ∈ R m×h are used for parametrizing the mean and covariance at each time step j, given the history information. b µ , b Σ ∈ R m are bias terms. For simplicity, let's use θ rnn = {θ

(1)

EQUATION

rnn , . . . , θ

(K)

rnn } to denote the collection of parameters in each RNN. On the boundary case, i.e., the starting point of each segment, we simply set h t = 0, which can be viewed as the setting according to the prior knowledge (bias terms) of RNN.

The above formulation indicates that the generative model P (x t |h t , z si ) depends not only on the last step observation x t-1 , but also the last hidden state h t-1 , which is together captured in Eq. 4. In summary, we denote all the parameters in the proposed R-HSMM as θ = {π, A, B, θ rnn }. The corresponding graphical model is shown in Figure 2b .

3. SEQUENTIAL VARIATIONAL AUTOENCODER:

To obtain the posterior or MAP in the proposed R-HSMM, the classical forward-backward algorithm or Viterbi algorithm needs to solve one dynamic programming per sample, which makes the inference costly, especially for the long sequence with thousands of timestamps. So instead, we treat the Bayesian inference from optimization perspective, and obtain the posterior by maximizing the negative Helmholtz variational free energy (Williams, 1980; Zellner, 1988; Dai et al., 2016) , max

EQUATION

over the space of all valid densities P. To make the optimization (6) tractable, the variational autoencoder restricts the feasible sets to be some parametrized density Q ψ , which can be executed efficiently comparing to the forward-backward algorithm or Viterbi algorithm. However, such restriction will introduce extra approximation error. To reduce the approximation error, we use a structured model, i.e., bidirectional RNN, to mimic the dynamic programming in forward-backward algorithm. Specifically, in the forward-backward algorithm, the forward message α t (z t , d t ) and backward message β t (z t , d t ) can be computed recursively, and marginal posterior at position t depends on both α t (z t , d t ) and β t (z t , d t ). Similarly, in bi-RNN we embed the posterior message with RNN's latent vector, and marginal posterior is obtained from the latent vectors of two RNNs at the same time step t.

Let ψ = {ψ---→ RNN1 , ψ← --- RNN2 , W z ∈ R h×K , W d ∈ R h×D } be the parameters of bi-RNN encoder, the Q ψ is decomposed as: Q ψ (z, d|x) = Q(z 1 |h 1 ; ψ)Q(d 1 |z 1 , h 1 ; ψ) |x| t=2 Q(z t |d t-1 , h t ; ψ)Q(d t |z t , d t-1 , h t ; ψ), (7) where h t = [ ---→ RNN 1 (x 1:t ), ←--- RNN 2 (x t:|x| )] is computed by bi-RNN. We use multinomial distributions Q(z t |h t ; ψ) = M(softmax(W z h t )) and Q(d t |z t , h t ; ψ) = M(softmax(W d h t )).

The dependency over d t-1 ensures that the generated segmentation (z, d) is valid according to Eq. 1. For example, if we sampled duration d t-1 > 1 from Q ψ at time t -1, then d t and z t should be deterministic. In our experiment, we use LSTM (Hochreiter & Schmidhuber, 1997) as the recursive units in bi-RNN.

Since with any fixed θ, the negative Helmholtz variational free energy is indeed the lower bound of the marginal likelihood, i.e., log

EQUATION

) therefore, we can treat it as a surrogate of the marginal log-likelihood and learn θ jointly with approximate inference, i.e.,

EQUATION

It should be emphasized that due to the discrete nature of latent variables in our model, the algorithm proposed in Kingma & Welling (2013) is not directly applicable, and its extension with stochastic neuron reparametrization (Bengio et al., 2013; Raiko et al., 2014; Gu et al., 2015; Chung et al., 2016) cannot provide satisfied results for our model according to our experiments. Therefore, we extend the penalty method to distribution space to solve optimization (9).

Algorithm 1 Learning sequential VAE with stochastic distributional penalty method

1: Input: sequences {x (n) } N n=1

2: Randomly initialize ψ (0) and θ 0 = {π 0 , A 0 , B 0 , θ 0 rnn } 3: for λ = 0, . . . , ∞ do 4:

for t = 0 to T do 5: Sample {x (n) } M

n=1 uniformly from dataset with mini-batch size M .

6::

Get {z (n) , d (n) } M n=1 with θ t by dynamic programming in (13).

7:

Update π t+1 , A t+1 , B t+1 using rule ( 16).

8:

θ t+1 rnn = θ t rnn , -γ t 1 M M n=1 ∇ θ t rnn Lλ (θ, ψ|x (n) )

9:

ψ t+1 = ψ t -η t 1 M M n=1 ∇ ψ t Lλ (θ, ψ|x (n) )

bi-rnn sequence to sequence learning 10:

end for 11: end for

4. LEARNING VIA STOCHASTIC DISTRIBUTIONAL PENALTY METHOD:

As we discussed, learning the sequential VAE with stochastic neuron reparametrization in unsupervised setting is extremely difficult, and none the off-the-shelf techniques can provide satisfied results. In this section, we introduce auxiliary distribution into (9) and generalize the penalty method Bertsekas (1999) to distribution space. Specifically, we first introduce an auxiliary distribution Q(z, d|x) for each x and reformulate the optimization (9) as

EQUATION

.t. KL Q(z, d|x (n) )||Q ψ (z, d|x (n) ) = 0, ∀x (n) , n = 1, . . . , N.

We enforce the introduced Q(z, d|x) equals to Q ψ (z, d|x) in term of KL-divergence, so that the optimization problems ( 9) and (10) are equivalent. Because of the non-negativity of KL-divergence, itself can be viewed as the penalty function, we arrive the alternative formulation of (10) as

EQUATION

where

Lλ (θ, ψ|x) = E Q(z,d|x) log P θ (x, z, d) -log Q(z, d|x i ) -λKL Q(z, d|x)||Q ψ (z, d|x)

and λ ≥ 0. Obviously, as λ → ∞, KL( Q(z, d|x)||Q ψ (z, d|x)) must be 0, otherwise the L∞ (θ, ψ|x) will be -∞ for arbitrary θ, ψ. Therefore, the optimization (11) will be equivalent to problem (10) . Following the penalty method, we can learn the model with λ increasing from 0 to ∞ gradually. The entire algorithm is described in Algorithm 1. Practically, we can set λ = {0, ∞} and do no need the gradually incremental, while still achieve satisfied performance. For each fixed λ, we optimize Q and the parameters θ, ψ alternatively. To handle the expectation in the optimization, we will exploit the stochastic gradient descent. The update rules for θ, ψ and Q derived below.

4.1. UPDATING Q:

In fact, fix λ, Q ψ and P θ in optimization (11), the optimal solution Q * (z, d|x) for each x has closed-form.

Theorem 1 Given fixed λ, Q ψ and

P θ , Q * (z, d|x) ∝ Q ψ (z, d|x) λ 1+λ P θ (x, z, d) 1 1+λ

achieves the optimum in (11) . Proof The proof is straightforward. Take the functional derivative of L w.r.t. Q and set it to zeros, d|x) . Take exponential on both sides, we achieve the conclusion. In fact, because we are using the stochastic gradient for updating θ and ψ later, Q * (z, d|x) is never explicitly computed and only samples from it are required. Recall the fact that Q ψ (z, d|x) has a nice decomposition 7, we can multiply its factors into each recursion step and still get the same complexity as original Viterbi algorithm for MAP or sampling. Specifically, let's define α t (j, r) to be the best joint log probability of prefix x 1:t and its corresponding segmentation which has the last segment with label j and duration r, i.e., α t (j, r) max

∇ Q L = log P θ (x, z, d) + λ log Q(z, d|x) -(1 + λ) log Q(z, d|x) = 0, we obtain 1 1+λ log P θ (x, z, d) + λ 1+λ log Q(z, d|x) = log Q(z,

z1:t,d1:t log Q(z 1:t , d 1:t |x 1:t ), s.t. z t = j, d t = d t-r = 1, d t-r+1 = r (12)

here t ∈ {1, 2, . . . , |x|}, j ∈ Z, r ∈ D. Then we can recursively compute the entries in α as below:

α t (j, r) =                    α t-1 (j, r -1) + 1 1+λ log( Bj,r Bj,r-1 P (x t |x t-r+1:t-1 , z = j)) r > 1, t > 1 + λ 1+λ log Q ψ (dt-r+1=r|z=j,x) Q ψ (dt-r+1=r-1|z=j,x) ; max i∈Z\j max r ∈D α t-1 (i, r ) + 1 1+λ log(A i,j B j,1 P (x t |z = j)) r = 1, t > 1 + λ 1+λ log Q ψ (z t-r+1 = j, d t-r+1 = r|x); λ 1+λ log Q ψ (z 1 = j, d 1 = r|x) + 1 1+λ log(π j B j,1 P (x 1 |z = j)); r = 1, t = 1 0.

otherwise (13) To construct the MAP solution, we also need to keep a back-tracing array β t (j, r) that records the transition path from α t-1 (i, r ) to α t (j, r). The sampling from Q(z, d|x) also can be completed with almost the same style forwards filtering backwards sampling algorithm, except replacing the max-operator by sum-operator in α propagation Murphy (2012) .

Without considering the complexity of computing emission probabilities, the dynamic programming needs time complexity O |x|K 2 + |x|KD (Yu & Kobayashi, 2003) and O(|x|K) memory. We explain the details of optimizing the time and memory requirements in Appendix A.

Remark: When λ = ∞, the Q(z, d|x) will be exactly Q ψ (z, d|x) and the algorithm will reduce to directly working on Q ψ (z, d|x) without the effect from P θ (x, z, d). Therefore, it is equivalent to obtaining MAP or sampling of the latent variables z, d from Q ψ (z, d|x), whose cost is O(|x|K). In practical, to further accelerate the computation, we can follow such strategy to generate samples when λ is already large enough, and thus, the effect of P θ (x, z, d) is negligible.

4.2. UPDATING θ AND ψ:

With the fixed Q(z, d|x), we can update the θ and ψ by exploiting stochastic gradient descent algorithm to avoid scanning the whole training set. Sample a mini-batch of sequences {x n } M n=1 with size M N , we proceed to update {θ, ψ} by optimizing the Monte Carlo approximation of (11),

EQUATION

where {z (n) , d (n) } is the MAP or a sample of Q(z, d|x (n) ). Note that the two parts related to θ and ψ are separated now, we can optimize them easily.

Update θ: Finding parameters to maximize the likelihood needs to solve the constrained optimization shown below

EQUATION

where {π, A, B} are constrained to be valid probability distribution. We use stochastic gradient descent to update θ rnn in totally K RNNs. For parameters π, A, B which are restricted to simplex, the stochastic gradient update will involve extra projection step. To avoid such operation which may be costly, we propose the closed-form update rule derived by Lagrangian, 

EQUATION

ψ: Given fixed λ, log Q ψ (z (n) , d (n) |x (n)

) is essentially the sequence to sequence likelihood, where the input sequence is x and output sequence is {z, d}. Using the form of Q ψ in Eq 7, this likelihood can be decomposed by positions. Thus we can conveniently train a bi-RNN which maximize the condition likelihood of latent variables by stochastic gradient descent.

Remark: We can get multiple samples {z, d} for each x from Q(z, d|x) to reduce the variance in stochastic gradient. In our algorithm, the samples of latent variable come naturally from the auxiliary distributions (which are integrated with penalty method), rather than the derivation from lower bound of objective (Tang & Salakhutdinov, 2013; Raiko et al., 2014; Mnih & Rezende, 2016) .

5. EXPERIMENTS:

Baselines We compare with classical HSMM and two popular HSMM variants. The first one is Hierarchical Dirichlet-Process HSMM (HDP-HSMM) (Johnson & Willsky, 2013) , which is the nonparametric Bayesian extension to the traditional HSMM that allows infinite number of hidden states; the second one is called subHSMM (Johnson & Willsky, 2014) , which uses infinite HMM as the emission model for each segment. This model also has two-level of latent structure. It considers the dependency within each segment, which is a stronger algorithm than HDP-HSMM. We also compare with the CRF autoencoder (CRF-AE) (Ammar et al., 2014) , which uses markovian CRF as recognition model and conditional i.i.d.model for reconstruction. Comparing to HSMM, this model ignores the segmentation structures in modeling and is more similar to HMM. Evaluation Metric We evaluate the performance of each method via the labeling accuracy. Specifically, we compare the labels of each single observations in each testing sequence. Since the labels are unknown during training, we use KM algorithm (Munkres, 1957) to find the best mapping between predicted labels and ground-truth labels. Settings Without explicitly mentioned, we use leave-one-sequence-out protocol to evaluate the methods. Each time we test on one held-out sequence, and train on other sequences. We report the mean accuracy in Table 1 . We set the truncation of max possible duration D to be 400 for all tasks. We also set the number of hidden states K to be the same as ground truth.

For the HDP-HSMM and subHSMM, the observation distributions are initialized as standard Multivariate Gaussian distributions. The duration is modeled by the Poisson distribution. We tune the concentration parameters α, γ ∈ {0.1, 1, 3, 6, 10}. The hyperparameters are learned automatically. For subHSMM, we tune the truncation threshold of the second level infinite HMM from {2 . . . 15}.

For CRF-AE, we extend the origin model for the continuous observations, and learn all parameters similar to M. Schmidt (2008) . We use mixture of Gaussians to model the emission, where the number of mixtures is tuned in {1, . . . , 10}.

For the proposed R-HSMM, we use Adam (Kingma & Ba, 2014) to train the K generative RNN and bi-RNN encoder. To make the learning tractable for long sequences, we use back propagation through time (BPTT) with limited budget. We also tune the dimension of hidden vector in RNN, the L 2 -regularization weights and the stepsize. We implemented with CUDA that parallelized for different RNNs, and conduct experiments on K-20 enabled cluster. We include both the R-HSMM with the exact MAP via dynamic programming (rHSMM-dp) and sequential VAE with forward pass (rHSMM-fw) in experiments. In all tasks, the rHSMM-fw achieves almost the same performance to rHSMM-dp, but 400 times faster, showing the bi-RNN is able to mimic the forward-backward algorithm very well with efficient computation.

5.1. SEGMENTATION ACCURACY Synthetic Experiments:

We first evaluate the proposed method on two 1D synthetic sequential data sets. The first data set is generated by a HSMM with 3 hidden states, where π, A, B are designed beforehand. A segment with hidden state z is a sine function λ z sin(ω z x + 1 ) + 2 , where 1 and 2 are Gaussian random noises. Different hidden states use different scale parameters λ z and frequency parameters ω z . The second data set also has 3 hidden states, where the segment with hidden state z is sampled from a Gaussian process (GP) with kernel function k z (x, y). Different hidden states employ different kernel functions. The specific kernel functions used here are k

1 (x, y) = exp{-min(| x-y | , |x + y|) 2 /10}, k 2 (x, y) = exp{-(x -y) 2 /10} and k 3 (x, y) = (5 -|x -y|)I{(5 -|x -y|) < 5}.

For both of the Sine and GP data sets, the duration of a segment is randomly sampled from a distribution defined on {1, ..., 100}, which depends on the hidden states. Thus, the segmentation task corresponds to finding out different functions embedded in the sequences. We visualize the segmentation results of ground truth and three competitors on Sine and GP data sets in Figure 1a and Figure 1b respectively, and report the numerical results in Table 1 . As we can see, R-HSMM provides much better results on even small segments, dramatically outperforms HSMM variants and CRF-AE. Also note that, the sine function depicts short term dependencies, while Gaussian process has long dependency that determined by the kernel bandwidth. This demonstrates the ability of R-HSMM in capturing the long or short term dependencies.

Human activity This dataset which is collected by Reyes-Ortiz et al. ( 2016) consists of signals collected from waist-mounted smartphone with accelerometers and gyroscopes. Each of the volunteers is asked to perform a protocol of activities composed of 12 activities (see Figure 3a for the details). Since the signals within an activity type exhibit high correlation, it is natural for RNN to model this dependency. We use these 61 sequences, where each sequence has length around 3000.

Each observation is a 6 dimensional vector, consists of triaxial measures from accelerometers and gyroscopes.

Figure 3a shows the ground truth and the segmentation results of all methods. Both rHSMM-dp and rHSMM-fw almost perfectly recover the true segmentation. They can also capture the transition activity types, e.g., stand to lie or sit to lie. The HSMM, HDP-HSMM and CRF-AE makes some fragmental but periodical segmentations for walking, caused by lacking the dependency modeling within a segment. The subHSMM also has similar problem, possibly due to the limited ability of HMM generative model.

Drosophila Here we study the behavior patterns of drosophilas. The data was collected by Kain et al. (2013) with two dyes, two cameras and some optics to track each leg of a spontaneously behaving fruit fly. The dimension of observation in each timestamp is 45, which consists of the raw features and some higher order features. See Figure 3b for the detail of the 11 behavior types. We perform leave-one-sequence-out experiment on 10 sequences of length 10000 each. Figure 3b shows the segmentation results on the prefix of one sequence, while Table 1 gives the mean accuracy on all sequences. Different from the previous experiment, where the human activity signals are relatively smooth, here the signals depict high variance. Different activities exhibit quite different duration and patterns. Also, the activity types changes frequently. The R-HSMM almost captured each changing point of activities with both long and short durations. The corresponding mean accuracy also outperforms the baselines. However, we observed there are some correct segmentations with wrong labels. This happens mostly to the short segments, in which the RNN doesn't have enough history established for distinguishing similar activity types.

Physionet The heart sound records, usually represented graphically by phonocardiogram (PCG), are key resources for pathology classification of patients. We collect data from PhysioNet Challenge 2016 (Springer et al., 2015) , where each observation has been labeled with one of the four states, namely Diastole, S1, Systole and S2. We experiment with both the raw signals and the signals after feature extraction. Regarding the raw signals (Heart dataset), we collect 7 1-dimensional sequences of length around 40000. The feature-rich dataset (PN-Full) contains 2750 sequences, where each of them consists of 1500 4-dimensional observations. We do 5-fold cross validation for PN-Full. The visualization of segmentation results are shown in Appendix B.4. As the results shown in Table 1 , our algorithm still outperforms the baselines significantly. Also for such long raw signal sequences, the speed advantage of bi-RNN encoder over Viterbi is more significant. Viterbi takes 8min to do one inference, while bi-RNN only takes several seconds. Our framework is also flexible to incorporate prior knowledge, like the regularity of heart state transition into HSMM.

5.2. RECONSTRUCTION:

In this section, we examine the ability of learned generative model by visualizing the reconstructed signals. Given a sequence x, we use recognition model to get the latent variables z and d, then use learned K generative RNNs to generate signals within each segment. For the ease of visualization, we show the results on 1D signal dataset in Fig. 4a and Fig. 4b .

From Fig. 4 we can see the generative RNN correctly captures different characteristics from signals of different segment labels, such as different frequencies and scales in Sine dataset, or the different variance patterns in GP dataset. This is essential to distinguish between different segments.

6. CONCLUSION:

We presented the R-HSMM, a generalization of HSMM by incorporating recurrent neural generative model as the emission probability. To eliminate the difficulty caused by such flexible and powerful model in inference, we introduced the bi-RNN as the encoding distribution via the variational autoencoder framework to mimic the forward-backward algorithm. To deal with the difficulty of training VAE containing discrete latent variables, we proposed a novel stochastic distributional penalty method. We justified the modeling power of the proposed R-HSMM via segmentation accuracy and reconstruction visualization. From the comprehensive comparison, the proposed model significantly outperforms the existing models. It should be emphasized that the structured bi-RNN encoder yields similar performance as the exact MAP inference, while being 400 times faster. Future work includes further speeding up of our algorithm, as well as generalizing our learning algorithm to other discrete variational autoencoder.

Published as a conference paper at ICLR 2017 

B MORE EXPERIMENT RESULTS:



B.1 SYNTHETIC DATASETS:

The reconstructed signals from the original signals are shown in Fig. 5 and Fig. 6 for sine dataset gaussian Process dataset respectively. We can see the reconstructed signal almost recovered the original signal. The RNN captured the key differences of states, such as the frequency and scale; while in gaussian process dataset, it also recovered the complicated pattern involving long term dependencies.

We show the confusion matrix of all methods on synthetic sine and gaussian process dataset in Figure 7 and Figure 8 respectively.

B.2 HUMAN ACTIVITY:

The confusion matrices of our method and two baseline algorithms are shown in Figure 9 .

In Figure 10 , we also show several other segmentation results on different testing sequences.

B.3 DROSOPHILA:

The confusion matrices of our method and two baseline algorithms are shown in Figure 11 .

Since each sequence is too long to be clearly shown in one figure, we split the segmentation results of one sequence into four parts, and show them in Figure 12 .

B.4 HEART SOUND RECORDS:

The confusion matrices of our method and two baseline algorithms are shown in Figure 13 . Also, we split the segmentation results of one sequence into four parts, and show them in Figure 14 . 

:



Figures and Tables:

Figure 1: Synthetic experiment results. Different background colors represent the segmentations with different labels. In the top row, the black curve shows the raw signal. (a) The Sine data set is generated by a HSMM with 3 hidden states, where each one has a corresponding sine function; (b) Similar to 1a, but the segments are generated from Gaussian processes with different kernel functions. The first two rows are our algorithms which almost exact locate every segment.

st = r and z (n) st = j) M n=1 |s (n) | Since we already have the segmentation solution, the total number of samples used for training is equal to the number of observations in dataset. The different RNNs use different parameters, and train on different parts of observations. This makes it easy for parallelized training.



Figure 3: Segmentation results on Human activity and Drosophila datasets. Different background colors represent the segmentations with different labels. In the top row, the black cure shows the signal sequence projected to the first principle component. The following two rows are our algorithms which almost exact locate every segment. (a) The Human activity data set contains 12 hidden states, each of which corresponds to a human action; (b) The Drosophila data set contains 11 hidden states, each of which corresponds to a drosophila action.

Reconstruction illustration on GP dataset.

Figure 4: Reconstruction illustration. The generative RNNs (decoders) are asked to reconstruct the signals from only the discrete labels and durations (which are generated from encoder).

Figure 5: More reconstruction illustration on Sine dataset.

Figure 7: Confusion matrix on Synthetic Sine dataset.

Figure 9: Confusion matrix on Human Activity dataset.

Figure 10: More segmentation results on Human Activity dataset.

Error rate of segmentation. We report the mean and standard deviation of error rate.

Table content: <table><tr><td>Methods</td><td>SINE</td><td>GP</td><td>HAPT</td><td>Drosophila</td><td>Heart</td><td>PN-Full</td></tr><tr><td>rHSMM-dp</td><td colspan="6">2.67 ± 1.13% 12.46 ± 2.79% 16.38 ± 5.03% 36.21 ± 1.37% 33.14 ± 7.87% 31.95 ± 4.32%</td></tr><tr><td>rHSMM-fw</td><td colspan="6">4.02 ± 1.37% 13.13 ± 2.89% 17.74 ± 7.64% 35.79 ± 0.51% 33.36 ± 8.10% 32.34 ± 3.97%</td></tr><tr><td>HSMM</td><td colspan="2">41.85 ± 2.38% 41.15 ± 1.99%</td><td>41.59 ± 8.58%</td><td colspan="3">47.37 ± 0.27% 50.62 ± 4.20 % 45.04 ± 1.87%</td></tr><tr><td>subHSMM</td><td colspan="2">18.14 ± 2.63% 24.81 ± 4.63%</td><td>22.18 ± 4.45%</td><td colspan="3">39.70 ± 2.21% 46.67 ± 4.22% 43.01 ± 2.35%</td></tr><tr><td colspan="3">HDP-HSMM 42.74 ± 2.73% 41.90 ± 1.58%</td><td>35.46 ± 6.19%</td><td colspan="3">43.59 ± 1.58% 47.56 ± 4.31% 42.58 ± 1.54%</td></tr><tr><td>CRF-AE</td><td colspan="6">44.87 ± 1.63% 51.43 ± 2.14% 49.26 ± 10.63% 57.62 ± 0.22% 53.16 ± 4.78% 45.73 ± 0.66%</td></tr></table>

Peer review (iclr2017):
This paper presents a novel model for unsupervised segmentation and classification of time series data. A recurrent hidden semi-markov model is proposed. This extends regular hidden semi-markov models to include a recurrent neural network (RNN) for observations. Each latent class has its own RNN for modeling observations for that category. Further, an efficient training procedure based on a variational approximation. Experiments demonstrate the effectiveness of the approach for modeling synthetic and real time series data. This is an interesting and novel paper. The proposed method is a well-motivated combination of duration modeling HMMs with state of the art observation models based on RNNs. The combination alleviates shortcomings of standard HSMM variants in terms of the simplicity of the emission probability. The method is technically sound and demonstrated to be effective. It would be interesting to see how this method compares quantitatively against CRF-based methods (e.g. Ammar, Dyer, and Smith NIPS 2014). CRFs can model more complex data likelihoods, though as noted in the response phase there are still limitations. Regardless, I think the merits of using RNNs for the class-specific generative models are clear.