Switching Poisson Gamma Dynamical Systems

1 Introduction

Temporal sequences are abundant in real world and analyzing them is always an utmost important task in machine learning. Among them, time-series count data has attracted wide attention to deal with a variety of real-world applications, such as text analysis, social network modeling and natural language processing. The widely-used dynamic models, such as hidden Markov models (HMMs) [Rabiner and Juang, 1986] and linear dynamic systems (LDSs) [Ghahramani and Roweis, 1999], have difficulty in modeling such data, which are often high dimensional, sparse, and overdispersed. To address this issue, several dynamic methods have been proposed, especially based on the Poisson-gamma structure [Zhou , 2012; Zhou et al., 2016]. Specifically, Poisson-gamma dynamical systems (PGDS) [Schein et al., 2016] is a typical dynamic model for count sequence analysis and performs well in capturing cross-factor temporal dependence, which has been wildly used in text analysis, international relation study and so on. Generally, current developments of dynamic models are always focusing on modeling more complex sequences.

∗ : Corresponding author

In this paper, we propose switching Poisson-gamma dynamical systems (SPGDS), which is a powerful dynamical model that can capture different dynamics among time-steps, so as to model complex sequential relationships efficiently. Specifically, we build a dynamic system with the assumption that the latent variables of each time-step are drawn from the gamma mixture distributions, whose shape parameter for each mixture component is factorized into a linear transformation of the latent unit of previous time-step. Combining the temporal structure and mixture model, SPGDS can not only benefit the transmitting of nonlinear and diverse temporal variation for ample representational capacity by gamma mixture distributions, but also enable our model to cluster diverse dynamics among time-steps into different patterns. We introduce a discrete indicator variable zt, called switching variable, to guide how the latent state θt varies from time t - 1 to time t, as illustrated in Fig. 1 (a). The switching mechanism in our model serves two benefits: (1) it results in better fitting and prediction performance on complex count sequences; (2) it provides an insight into which dynamical patterns contains within sequences, which can be of value in application and analysis [Fraccaro et al., 2017; Becker Ehmck et al., 2019]. With the assumption that the real-world data changes in dynamics at time t are causality, we assign zt to non-linearly depend on history input x1:t-1 via a Gumble-softmax based recurrent variational inference network. Based on this, we further develop a Weibulldistribution-based switching recurrent variational inference network for structured inference [Krishnan et al., 2017] of latent variable θt. This structured inference network enables SPGDS to learn rich latent representations and fast in outof-sample prediction. The detail structure of our inference network is shown in Fig. 1 (b). Moreover, we develop a minibatch based stochastic inference algorithm that combines stochastic gradient MCMC (SG-MCMC) [Patterson and W, 2013; Ma et al., 2015] and autoencoding variational inference, which accelerates our model in both training and testing phase for large scale sequences. Furthermore, to prove the flexibility and compatibility of our model with prevalent deep learning structures, we extend SPGDS into supervised SPGDS (sSPGDS), which incorporates the label information into the model and extract discriminative features to achieve enhanced performance on both data representation and classification.

2 RelatedWork

To model count sequence data, several dynamical model based on the Poisson-gamma construction has been proposed. The gamma process dynamic Poisson factor analysis (GPDPFA) [Acharya et al., ] assumes that the data comes from the Poisson distribution and models the count vector at each time step under Poisson factor analysis (PFA) [Zhou et al., 2012] as xt ∼ Pois (Φθt). It further smoothes the translation through time by assigning θt ∼ Gam(θt-1, βt). Poissongamma dynamical systems (PGDS) [Schein et al., 2016] further improves the ability to capture cross-factor temporal dependence by introducing a translation matrix as θt ∼ Gam(Πθt-1, β). Moreover, to capture long-range temporal dependencies and model long sequences, several multilayer probabilistic dynamic models are proposed. For instance, deep dynamic Poisson factor analysis (DDPFA) [Gong and Huang, 2017] combines recurrent neural network (RNN) [Martens and Sutskever, 2011] with PFA, to capture longrange temporal dependencies of the latent factor via RNN. Deep temporal sigmoid belief network (DTSBN) [Zhe et al., 2015] is an extension of deep sigmoid belief network (DSBN) [Gan et al., 2015] with sequential feedback loops on each layer. However, DTSBN restricts its hidden units to be binary, which limits its representational power. [Guo et al., 2018] extends PGDS into deep PGDS (DPGDS), a model with deep hierarchical latent structure and captures the correlations between the features across layers and over times using the gamma belief network [Zhou et al., 2016].

Although these methods exhibit attractive performance in describing temporal dependencies, due to the weights sharing across different time steps, there are still some difficulties for them in modeling sequences with very complex and highly nonlinear dynamics. However, switching linear dynamical systems (SLDS) [Linderman et al., 2016; Fraccaro et al., 2017; Becker Ehmck et al., 2019] is a widely used method to model those complex sequences by breaking down them into multiple simpler units and modeling them separately. But, it is difficult for SLDS to model count sequences for it is under the Gaussian assumption. In this paper, the proposed SPGDS could be seen as a novel switching extension of Poisson-gamma structure models, which not only fit the count sequence, but also inherits various virtues of switching dynamic models. Moreover, our model is also fast in out-of-sample prediction with the help of the structure variational inference network.

3 Switching Dynamical Systems for Count Sequences Modeling

In this section, we first propose a switching PGDS model for analysing count sequences. Then, we propose switching recurrent inference network to map the observation directly to the switching variables and latent representation of SPGDS, by which our model has benefited fast inference for out-ofsampling prediction.

3.1 Switching Poisson Gamma Dynamical Systems

We demonstrate the graphical representation of SPGDS in Fig. 1 (a). Assuming that dataset of V -dimensional sequentially observed multivariate count data x1, ..., xt are represented as a V × T count matrix X. The generation process of SPGDS can be expressed as

where the latent factors Φ, θ, {Πc}Cc=1 , {τc}Cc=1 and z are all positive. The input count data xt of time t are factored into the loading matrix Φ ∈ RV×K + and the corresponding hidden states θt ∈ RK+ . And δ ∈ R+ is the scaling parameter. We characterize the relationship between the hidden units of adjacent time-steps as multiple gamma distribution, so as to characterize the nonlinear sequential relationship between time-steps. zt = (zt1, ..., ztC) is a C-dimensional categorical variable which choosing parameters Πc and τc at time t from different gamma distribution, we call it switching variable. Marginalizing zt in (1), we get

which is a mixture of gamma distribution to characterize the complex and diverse sequential relationship between timesteps more accurately than other Poisson-gamma structure models. {Πc ∈ RK×K + }Cc=1 are the latent transition matrices that capture the various temporal dependencies between components, and C represents the number of components in mixture distribution. We denote {τc ∈ R+}Cc=1 as the scaling parameters that control the variousness of temporal amplitude variation of the hidden states. rt = (r1 t , ..., rC t )T is the parameter for categorical distribution. In this way, we can also label sequences into segmentations that exhibit different dynamics with zt. The vector θt has an expected value

that {Πc}Cc=1 play the role of transiting the latent presentation across time, and controlled by rt. Hence, our model is capable to capture the primitive patterns of temporal sequences concisely and accurately, when the data contains multiple dynamics.

To complete the diverse dynamic model, we introduce K factor weights in each components of mixture of gamma distribution:

where πck = (πc1k , . . . , πcKk ) is the kth column of Πc and {πck1k }Cc=1 can be interpreted as the C different probabilities for transitioning from component k to component k1 for different dynamics. For the latent state at the first timestep, we define its prior as θ1 ∼ Gam (τ0v, 1 c0 ). Moreover, we place Dirichlet priors over the feature factors and draw the other parameters from a noninformative gamma prior: φk = (φ1k , . . . , φKk ) ∼ Dir (η, . . . , η) and ξ (c) , β(c) ∼ Gam(ε0, ε0), τ (c) ∼ Gam(α0, 1/β0) .

In particular, when C = 1, SPGDS reduces to PGDS. For non-count observations, we use Bernoulli-Poisson distribution [Zhou, 2015] and Poisson randomized gamma distribution [Aaron et al., 2019] to link binary observation and Nonnegative-real-valued observation to latent poisson count, respectively. More details can be found in [Schein et al., 2016].

3.2 Switching Recurrent Inference Network for SPGDS

For efficient out-of-sample prediction, we develop a switching recurrent inference network, which will be used in hybrid SG-MCMC and variational inference method described in section 4, to map the observation directly to the latent variables. Specifically, we use Concrete distribution [Maddison et al., 2016] to approximate the categorical distributed switching variables, Weibull distribution [Zhang et al., 2018] to approximate gamma distributed conditioned latent representations.

Gumble-softmax based recurrent variational inference network for ztn: Assuming there are N count sequential data {x1n, ..., xTn}Nn=1 . It is clear that categorical variable ztn in control of the changes of latent states θtn from time t - 1 to t. In more general setting, the changes in dynamics at time t depend on the history of system and determined by the input ranging from 0 to t - 1 [Fraccaro et al., 2017; Becker Ehmck et al., 2019]. Here, we let ztn determined by a learnable function of x1:t-1,n and modulate it by a recurrent variational parameter inference network. We construct the autoencoding variational distribution as q(ztn) = Categorical(rtn) and transform the observation to rtn using recurrent structure as:

rtn = softmax (Wxcxt-1,n +Wccrt-1,n + bxπ) (2) To obtain samples from categorical distribution, and to backpropagate through the categorical latent variables, we use q(˜ztn) = Gumble - softmax(rtn) to approximate q(ztn) [Maddison et al., 2016; Jang et al., 2016]. It draws samples via˜

where λ denotes the softmax temperature and ϵctn is the standard uniform variable. As λ approaches 0, samples from the Gumbel-softmax distribution become one-hot and the Gumbel-softmax distribution q(ztn) becomes identical to the categorical distribution q(ztn).

Weibull distribution based switching recurrent variational inference network for θtn: Following [Zhang et al., 2018], we approximate the gamma distributed conditional posterior of θtn with Weibull distribution by assigning q (θtn|ztn) ∼ Weibull (ktn, λtn), the random sample ϵtn can be obtained by transforming standard uniform variables ϵtn as: θtn = λtn(- ln (1 - ϵtn))1/ktn , where ktn and λtn are the parameters of θtn and they are nonlinearly transformed from the hidden units ht as

where Whk ∈ RK×K ,Whλ ∈ RK×K , b1 ∈ RK , b2 ∈ RK . A nonlinear transformation deterministically covert htn from xtn and ht-1,n. To exploit the various temporal information, we propose a switching recurrent inference network considering diverse temporal dependence across time-steps, as illustrated in Fig. 1 (b). Therefore, htn can be expressed as

where f is a nonlinear transformation function. Wxh ∈ RV×K is a input-hidden weight matrix, and {W chh ∈ RK×K }Cc=1 denote hidden states connected matrices. The detail structure of our inference network is shown in Fig. 1 (b).

4 Hybrid SG-MCMC and Variational Inference

In this section, we provide a hybrid stochastic gradient MCMC and autoencoding variational inference for SPGDS, which is scalable in training phase and fast in testing phase. Aforementioned in Section 3, we develop a switching recurrent autoencoding variational inference network for switching variable zt and latent representation θt, which enables our model be fast in testing phase. For inferencing global parameters in SPGDS, including Φ, {Πc}Cc=1 and {τc}Cc=1 , the topic-layeradaptive stochastic gradient Riemannian (TLASGR) MCMC algorithm [Cong et al., 2017; Zhang et al., 2018], which is proposed to sample simplex-constrained global parameters in a mini-batch learning setting and improve its sampling efficiency by using the Fisher information matrix (FIM), can be extended to our model. More specifically, after sampling auxiliary latent counts using augmentable techniques as in [Guo et al., 2018], π(c) k , kth column of the transition matrix Πc of component c, can be sampled as

where ˜Z comes from the augmented latent counts and the definition of ρ, εn, [·]∠ andM c k are analogous to [Cong et al., 2017]'s setting. The update of Φ is the same as [Guo et al., 2018]. For {τc}Cc=1 , which are one-dimensional non-negative data, it is efficient to sample them with stochastic gradient Langevin dynamics [Welling and Teh, 2011] as

Given the global parameters Φ, {Πc}Cc=1 , {τc}Cc=1 , the task is to optimize the parameters of the switching recurrent inference network. As the usual strategy of autoencoding variational inference, this optimization can be achieved by minimizing the negative Evidence Lower Bound (ELBO). We can express the ELBO of our inference network as:

Denoting Ω as the parameters of the inference network: Ω = {Wxc,Wcc,Wxh, {W chh }Cc=1 ,Whk,Whλ}, the corresponding hybrid SG-MCMC and switching recurrent autoencoding variational inference method for SPGDS is described in Algorithm.1.

Furthermore, our model can be extended to a supervised version called supervised SPGDS (sSPGDS) to handle the categorization task of sequential data. We achieve this by Algorithm 1 Hybrid stochastic-gradient MCMC and autoencoding variational inference for SPGDS

Set mini-batch size as M, the width of layer K and hyperparameters.

Initialize inference model parameters Ω and generative model parameters D = {Φ, {Πc, τc}Cc=1}.

for iter = 1, 2, ... do

Randomly select a mini-batch of time sequential data to form a subset {xt,m}T,M t=1,m=1 ;

Draw random noise {ϵt,m, ˜ϵt,m} T,M t = 1,m= 1 from uniform distribution for sampling latent states θt,m and zt,m;

Calculate subgradient ∇ΩL (Ω, D; xt,m, ϵt,m, ˜ϵt,m) according to (7), and update Ω using the subgradient;

Update M c k according to eqn. (18) and eqn. (19) in [Cong et al., 2017]; then πck with (5);

Update τc with (6);

Update φk according to eqn. (15) in [Cong et al., 2017];

end for end for

concatenate the latent states across all time-steps to construct a T ×K-dimensional latent feature and add the softmax classifier on it [Wang et al., 2019]. Then, the loss function of the entire framework should be modified as

where Lg refers to ELBO of generative model shown in (7), Ls denotes the classification criterion, and ξ is a tradeoff parameter to balance aspects of generation and classification.

5 Experiments

In this section, We examine the performance of the proposed model on both unsupervised and supervised tasks.

5.1 UnsupervisedModels

We first examine the performance on fitting and prediction tasks of our proposed model, both synthetic datasets and realworld datasets are exploited here. Besides, we compare our model with some existing dynamic methods introduced in section 2, including HMM [Rabiner and Juang, 1986], LDS [Ghahramani and Roweis, 1999], GP-DPFA [Acharya et al., ], TSBN [Zhe et al., 2015] and PGDS [Schein et al., 2016]. We set the hyperparameters of GP-PFA, TSBN and PGDS the same as their original settings.

Synthetic Dataset

Inspired by [Gong and Huang, 2017], we consider three different multi-dimensional synthetic datasets:

Toydata 2: f1 (t) = t, f2 (t) = 2 mod (t, 3), f3 (t) = 200 exp (-t/3) for t = 1, . . . , 50 and f1 (t) = 2t + 30, f2 (t) = 3 mod (t, 2) + 50 and f3 (t) = 30t exp (-t) + 100 for t = 51, . . . , 100, where mod(t, n) denotes the modulo operation which returns the remainder after division of t by n.

Toydata 3: f1 = 5t, f2 = 10t, f3 = 10t + 2 for t = 1, . . . , 50 and f1 (t) = f1 (50) +mod(t, 2), f2 (t) = f2 (50) + 2 mod (t, 2) + 2, f3 (t) = f3 (50) + mod (t, 10) for t = 51, . . . , 100 and f1 (t) = mod(t, 3), f2 (t) = 2 mod (t, 2) + 2, f3 (t) = mod (t, 5) for t = 101, . . . , 150.

We set the number of latent states as K = 2, and compare the proposed SPGDS with PGDS, HMM and LDS on Mean Square Error (MSE) between the ground truth and the estimated value and Prediction Mean Square Error (PMSE), which is the MSE between the ground truth and the prediction in the next time-step. The best performance of different methods are listed in Tab. 1. Clearly, SPGDS has the best performance in fitting and prediction tasks on all datasets. We attribute these to our model's ability of capturing diverse temporal dependencies and modeling them separately in the latent space. Taking Toydata 2 as an example, SPGDS captures three dynamic patterns at different time-steps, which are t = 1 : 7, t = 50 : 51 and t = 7 : 50, 51 : 99, as shown in Fig. 2. Corresponding transition matrices and weights to these three dynamic patterns are shown in Tab. 2. For stable sequential variation as dynamic 1, the transition matrix more closely approaches a diagonal matrix and transition weight is approximately equal to one. Relatively, as dynamic 2 and 3, transition matrix and weight changes with various temporal dependencies.

Real-world Datasets

Following [Gong and Huang, 2017; Schein et al., 2016], five real-world datasets are used:

• Global Database of Events, Language, and Tone (GDELT): GDELT is an international relationship dataset, which is extracted from news corpora.
• Integrated Crisis Early Warning System (ICEWS): ICEWS is another international relationship dataset extracted from news corpora.
• State-of-the-Union transcript (SOTU): The SOTU dataset contains the text of the annual SOTU speech transcripts from 1790 to 2014.
• DBLP conference abstract (DBLP): DBLP corpus is a database of computer research papers.
• NIPS corpus (NIPS): The NIPS contains the text of every NIPS conference paper from 1987 to 2003.

For each of these datasets, we summarize it as a V × T matrix, as shown in Tab.3. Specifically, we set V = 1000 for all data by choosing the top 1000 most frequently used features. Similar as previous methods [Zhe et al., 2015], we evaluate the prediction performance of our model by calculating the precision and recall at top-M as in [Han et al., 2014], which is given by the fraction of the predicted top-M words, that matches the true ranking of the words. M is set as 50 here. We use three criterion MP, MR and PP. MP and MR are mean precision and mean recall over all years that appears in the training set, PP is the prediction precision for the final year. Moreover, we employ the setup in [Zhe et al., 2015] that the entire data of the last year is held-out, while the words of the each document for the documents in the previous years are randomly partitioned into 80%/20% split. The 80% portion is used to train the model, and the prediction at the next year is tested on the rest of 20% held-out words. We compare the proposed model with several related works, including GPDPFA, TSBN and PGDS, and the results are summarized in Tab. 3. Fig. 3 presents the different number of mixture components influence to the results in SPGDS on NIPS dataset. They have similar trends with different number of mixture components. Thus, cross validation where we set MSE as metric is used here to determine the number of mixture components in SPGDS for each dataset. We select Cg = 5, Ci = 5, Cs = 3, Cd = 2, Cn = 3 for datasets from left to right in Tab. 3 Clearly, SPGDS outperforms other methods on most of the evaluation criteria which we contribute to the superiority of SPGDS in modeling nonlinear sequential data. In Fig. 4, we present the visualization of top four latent factors inferred by SPGDS with three mixture components from the NIPS. The three dynamical patterns captured by our model, including the decline of research on neuron from 1987 to 1991, the relatively stable phase from 1991 to 1994, the decline of research on neural network and the raised of research on bayesian learning from 1994 to 2002. The sustainable growth of topics on "learning, network" indicates the increase in the number of papers accepted by NIPS.

5.2 SupervisedModels

To evaluate how well sSPGDS leverages the label information for feature learning, we compare its classification performance with a variety of algorithms on sequential MNIST dataset and permuted sequential MNIST dataset. For sequential MNIST, the pixels of MNIST digits [LeCun et al., 1998] are presented sequentially to the network and classification is performed at the end. By permuting the image sequences by a fixed random order, we can get permuted sequential MNIST dataset. Since the MNIST image are 28 × 28 pixels, they are reshaped into 784 × 1 sequences in sequential MNIST.

We compare our model with (a)RNN(relu), that is an RNN with RELU activations, (b) IRNN [Le et al., 2015], that is an RNN with RELU activations and with the recurrent weight matrix initialized to the identity, (c) Unitary evolution RNN (uRNN) [Arjovsky et al., 2016], that uses orthogonal and unitary matrices in RNN, (d) Full-Capacity Unitary RNN [Wisdom et al., 2016], which is a modified version of uRNN, (e) Skip RNN [Campos et al., 2018], an extending of existing RNNmodels, that learns to skip state updates and shortens the effective size of the computational graph, (f) long short-term

0.8

0.6

0.4

0.2

Table 4: Supervised Results on MNIST

memory (LSTM) [Greff et al., 2016] and (g) the gated recurrent unit (GRU) [Cho et al., 2014], which are the variants of RNN to address the gradient problem, (h) supervised Poisson gamma dynamical system, a supervised extension of PGDS. The classification accuracy of different methods are shown in Tab. 4, where the results of compared methods are provided in their corresponding papers. The latent dimension of models are 100 and Tensor Recurrent autoencoding network is developed based on RNN (relu). Clearly, our proposed sSPGDS can achieve a comparable performance among these methods and we contribute to the two characteristics of our model: (1) uncertainty is included in our hidden states [Chung et al., 2015] and (2) mixture distribution is assigned to our latent variables. These two characters enable sSPGDS to be robustness when track with complex and nonlinear sequential data. In addition, we show various dynamics captured by sSPGDS in Fig. 5. The number of mixture components is set as three. We reshape I(zt), which is the index for zt in (1), from t = 2 to t = 785 into 28 × 28 pixel and each pixel represents the index for dynamical patterns from the corresponding pixel in data to its next pixel. We visualize it by assign different colors to I(zt) with different values: I (zt) = 1 → blue, I (zt) = 2 → green, I (zt) = 3 → yellow, and show it in the second row of Fig. 5. As we can see, our model captures three different dynamical patterns including: changing within noise or target, changing from noise to target and changing from target to noise. Moreover, we also visualize the probability of p(I(zt) = 1) , p(I(zt) = 2) and p(I(zt) = 3) in third to fifth row of Fig. 5.

6 Conclusion

In this paper, we propose switching Poisson gamma dynamical systems (SPGDS) that takes advantage of gamma mixture distributions to model complex and nonlinea temporal dependencies, while capturing various dynamic patterns. To model the dependency of dynamic patterns among different timesteps and achieve fast out-of-sample prediction, a switching recurrent variational inference network is developed to infer the switching variable and latent representation. A minibatch based stochastic inference method that combines both stochastic-gradient MCMC and variational inference algorithm is developed to accelerate both training and testing for large scale sequences. In addition, we provide a supervised extension of SPGDS. Experiment results show that our model not only has an excellent fitting and prediction performance on unsupervised feature extraction tasks, but also achieves comparable classification performance on supervised tasks.

Acknowledgements

B. Chen acknowledges the support of the Program for Young Thousand Talent by Chinese Central Government, the 111 Project (No. B18039), NSFC (61771361), NSFC for Distinguished Young Scholars (61525105), Shaanxi Innovation Team Project, and the Innovation Fund of Xidian University. M. Zhou acknowledges the support of the U.S. National Science Foundation under Grant IIS-1812699.