Probabilistic Transformer for Time Series Analysis
Year: 2021

Abstract:
Generative modeling of multivariate time series has remained challenging partly due to the complex, non-deterministic dynamics across long-distance time steps. In this paper, we propose deep probabilistic methods that combine state-space models (SSMs) with transformer architectures. In contrast to previously proposed SSMs, our approaches use attention mechanism to model non-Markovian dynamics in the latent space and avoid recurrent neural networks entirely. We also extend our models to include several layers of stochastic variables organized in a hierarchy for further expressiveness. Compared to transformer models, ours are probabilistic, non-autoregressive, and capable of generating diverse long-term forecasts with accounted uncertainty. Extensive experiments show that our models consistently outperform competitive baselines on various tasks and datasets, including time series forecasting and human motion prediction.

heading: 1 Introduction
text: Generative modeling of multivariate time series is a challenging problem with wide-ranging applications in demand forecasting [15, 76], autonomous driving [2, 16], robotics [29, 67], and health care [20, 21, 59]. Despite remarkable progress in recent years, models that predict high-dimensional future observations from a few past examples have remained intractable, partly due to the complex, non-deterministic temporal dynamics across long-distance time steps. Given a sequence of human poses, for example, such models must internally figure out the involved dynamics of various body components across space and time while maintaining the inherent uncertainty of multiple plausible futures, even though only one such future is observed.
Among proposed probabilistic approaches, state space models (SSMs) provide a principled framework for learning and drawing inference from sequential inputs [27, 66]. While autoregressive models feed its predictions back into the dynamics model without any compressed representation of data, SSMs model stochastic transitions between abstract states using latent variables, allowing for efficient state-to-state sampling without the need to render high-dimensional observations. Gaussian linear dynamical systems (LDSs), one of the best known SSMs [92], for example, postulate linear state transitions and enjoy exact inference via the celebrated Kalman filter algorithm.
While early extensions of LDSs focus on linearization [46] and unscented transform [88], recent work that marry state space models with deep neural networks offers much more flexibility to model complex dependencies across different time steps. Some approaches retain the Markovian dynamics of LDSs and only replace their linear observation models with feed-forward networks [23, 31, 47, 71], whereas others favor nonlinear state transitions and parametrize such dependencies via recurrent neural networks (RNNs) [22, 23, 30, 39, 51, 75]. Despite differences, both Markovian transitions and RNNs are often not capable of capturing long-range dependencies in highly structured sequential inputs [36, 100], limiting the capacity of the corresponding SSMs.
35th Conference on Neural Information Processing Systems (NeurIPS 2021).
In this work, we propose to combine the complementary strengths of SSMs and transformer architectures [85], a powerful mechanism for modeling long-term interactions that enjoys success across a variety of sequence modeling tasks [26, 48, 99]. In contrast to most SSMs, our models make extensive use of attention mechanism [5, 85] between latent variables to model non-Markovian dynamics (see Figure 1). Compared to transformer-based methods, our models are probabilistic, nonautoregressive in a similar fashion to LDSs, and capable of generating diverse long-term forecasts with uncertainty estimates.
Our main contributions are threefold. First, we propose novel SSMs based on transformer architectures for multivariate time series, which include generative models and inference procedures based on variational inference [49, 74]. Second, we extend our models to include several layers of stochastic latent variables organized in a hierarchy for further expressiveness. Third, we conduct extensive experiments on time series forecasting and human motion prediction and demonstrate that our Probabilistic Transformer (ProTran) performs remarkably well compared to various state-of-theart baselines.

heading: 2 Preliminaries

heading: 2.1 Variational State Space Models
text: Let {x(i)1:Ti} N i=1 consist of N univariate time series where x (i) 1:Ti = (x (i) 1 ,x (i) 2 , · · ·x (i) Ti ) and x(i)t denotes the vaue of the i-th time series at time t. We consider the multivariate form x1:T = (x1,x2, . . . ,xT ) where xt = (x (1) t , . . .x (N) t ) ∈ RN . Conditioning on observed values up to time C, we aim to produce distributional forecasts into the future p(xC+1:T |x1:C). For clarity, we refer to x1:C and xC+1:T as contexts and targets, respectively.
We are interested in probabilistic models parametrized by θ of the form
pθ(x1:T |x1:C) = ∫ pθ(x1:T | z1:T )pθ(z1:T |x1:C)dz1:T (1)
where z1:T = (z1, z2, . . . , zT ) denotes the corresponding sequence of latent variables, sometimes referred to as states. In other words, we assume a generative model that can be decomposed into a transition model pθ(z1:T |x1:C) between the latent variables conditioned on the contexts, and an emission model pθ(x1:T | z1:T ) from the latent variables to observable outputs. In particular, we
further impose several assumptions on both models: 1
pθ(z1:T |x1:C) = T∏ t=1 pθ(zt | z1:t−1,x1:C), pθ(x1:T | z1:T ) = T∏ t=1 pθ(xt | zt). (2)
As demonstrated in Figure 1(b), the latent variable zt+1 depends not only on zt but also on all of its preceding latent variables, including zt−1, in contrast to linear dynamical systems (LDSs). In addition, the transition and emission models allow for non-linearity via neural network parametrizations. These assumptions aim to maximize model capacity for real-world applications with complex emissions or temporal dependencies.
However, neither x1:t−1 nor z1:t−1 are included in the emission model p(xt | z1:T ,x1:C). Such assumptions are important, as it has been argued previously that a leakage of information from the latent space in autoregressive models can hinder long-term predictions [23, 47]. While all ground truth observations are available during training, the entire sequence has to be generated sequentially at test time, making the dependencies on x1:t−1 prone to accumulated errors over multiple time steps. By letting the latent variable zt capture all information needed to render xt, we also avoid the computational costs associated with repeatedly decoding and encoding xt in multi-step predictions.
The inclusion of nonlinear state transitions and observation models necessarily requires approximate inference. We follow the stochastic variational inference framework [49, 74] and assume that the variational posterior parametrized by φ can be decomposed auto-regressively as qφ(z1:T |x1:T ) =∏ t qφ(zt | z1:t−1,x1:T ), which leads to a lower bound on the log likelihood:
log pθ(x1:T |x1:C) ≥ T∑ t=1 (Eq [log pθ(xt|zt)]− KL(qφ(zt|z1:t−1,x1:T ) ‖ pθ(zt|z1:t−1,x1:C))) ,
(3)
where KL is the Kullback-Leibler divergence.
For computational stability, we assume homoscedasticity and choose Laplace distribution with scale parameter β as a parametric form for pθ(xt | zt), i.e. we optimize for L1 reconstruction loss with a cross-validated factor β for the KL term, following similar variational autoencoder (VAE) work [24, 41, 86]. Such an assumption does not necessarily limit the capacity of our models, as powerful stochastic transitions and flexible emission models can theoretically characterize arbitrary noise covariance [66]. Incorporating structured probabilistic outputs such as Gaussian copulas [75] or normalizing flows [23] can potentially further improve our model performance.

heading: 2.2 Transformer Architectures
text: Central to our models and other transformer-based approaches [48, 85] is the notion of attention [5], which allows the models to focus on important parts within a context. Multi-head attention, for example, maps a sequence of queries Q ∈ R`q×d of length `q to a sequence of outputs O = [O1, . . . ,OH ] ∈ R`q×d of the same size by attending over given `k key-value pairs K ∈ R`k×d, V ∈ R`k×d:
Oh = Attention(Qh,Kh,Vh) = Softmax
( QhK
T h√
d
) Vh, (4)
where Qh = QW Q h , Kh = KW K h , Vh = VW V h are projected queries, keys, and values corresponding to head h ∈ [1,H] with learning parameters WQh ,WKh ,WVh , respectively. In case Q = K = V, we refer to such an attention mechanism as self-attention.
Given fully observed sequences of inputs, the mapping can be computed efficiently without any imposed sequential order often seen in recurrent neural networks [19, 42]. More importantly, the direct connections between long-distance time steps are baked into the mechanism as information from previous time steps is easily accessible without being compressed into a fixed representation, easing optimization and learning of long-term dependencies [5, 85].
1For notational simplicity, we assume x0 = z1:0 = ∅ and p(z |x0) = p(z).
Without recurrence, Transformer [85] encodes information about each time step t with pred efined sinusoidal positional embeddings Position(t) = [pt(1), . . . , pt(d)] ∈ Rd where the i-th embedding is given by pt(i) = sin(t · ci/d) for even i and pt(i) = cos(t · ci/d) for odd i and c is some large constant. Empirical results show that such positional embeddings are also important to our models.

heading: 3 Probabilistic Transformer
text: In this section, we first present our single-layered model and subsequently its multi-layered extension for a hierarchy of stochastic latent variables. As alluded earlier, our model consists of a generative model and an inference model that share information and parameters extensively.

heading: 3.1 Single-Layered Probabilistic Transformer
text: Generative Model. Given some contexts x1:C , we first apply a linear projection and combine it with a positional embedding to obtain h1:C ∈ Rd, i.e.
ht = LayerNorm(MLP(xt) + Position(t)), (5)
where LayerNorm and MLP denote layer normalizations [4] and multi-layer perceptrons, respectively. While a traditional transformer model often dedicates an entire encoder for the same purpose [55, 72], we find such a simple mapping works sufficiently well in conjunction with the contextattention module of the corresponding decoder.
As implied in Equation (2), our latent dynamics decomposes auto-regressively. At each time step, we parametri ze the distribution pθ(zt | z1:t−1,x1:C) by a Gaussian with parameters resulting from two sequential steps of attention: a self-attention over the previously inferred states z1:t−1 and another attention over the projected contexts h1:C . These two operations mirror those found in the decoder of Transformer [85], with the stochastic latent variables replacing its decoder inputs.
Unfortunately, using stochastic samples of zt as attention queries is problematic, as purely stochastic transitions make it difficult for the model to reliably retain information across multiple time steps [17, 30, 39]. We therefore encapsulate the latent variables in hidden representations wt that also has a deterministic component. Combined with the attention steps, such representations help model long-range temporal dependencies while accounting for the stochasticity of future observations.
Starting with a learnable, context-agnostic representation w0, we recursively update wt using a stochastic sample from pθ(zt | z1:t−1,x1:C) and th e positional embedding for the current time step t. The generating process for the time step t can be summarized by the following pseudocode:
w̄t = LayerNorm(wt−1 + Attention(wt−1,w1:t−1,w1:t−1)) (6) ŵt = LayerNorm(w̄t + Attention(w̄t,h1:C ,h1:C)) (7) zt = Sample(N (zt;MLP(ŵt), Softplus(MLP(ŵt)))) (8) wt = LayerNorm(ŵt + MLP(zt) + Position(t)), (9)
where Sample and Softplus are the Gaussian sampling and approximating rectifier operators.
Each stochastic sample of w1:T is then mapped to a sequence of x1:T via a multi-layer perceptron. We emphasize that our generation procedure in the latent space is more efficient than others in the observation space, which requires encoding and decoding high-dimensional inputs repeatedly.
Inference Model. We parametrize the approximate posterior qφ(zt | z1:t−1,x1:T ) at time step t in a simi lar fashion to the prior pθ(zt | z1:t−1,x1:C). Indeed, these parametrizations share most parameters and are done simultaneously in the same recursive loop, following the exact same steps in Equation (6) and Equation (7) (see Figure 1). We note that similar sharing techniques between the generative and inference processes have emerged as a common theme among recent successful VAE models [17, 62, 83].
While the prior only has access to the conditioning observations x1:C , the approximate posterior should take into account all observations during training, including the targets xC+1:T . Due to the inherent unidirectional aspect of RNNs, previous work that uses RNNs to parametrize the approximate posterior often disregards such a property [22, 30, 51] and often resorts to a filtering routine
p(zt | z1:t−1,x1:t). In contrast, our inference procedure resembles more of the smoothing process of LDSs, factoring in both past and future observations via another application of self-attention:
kt = Attention(h1:T ,h1:T ,h1:T )) (10) zt = Sample(N (zt;MLP([ŵt,kt]), Softplus(MLP([ŵt,kt]))). (11)
Here, we replace Equation (8) in the generative model with Equation (11), where the hidden representation kt summarizing all information relevant to the current tim estep t has been concatenate to the latent-and-context-aware representation ŵt preceding the Gaussian parametrization.
The generative model and the inference model are trained end-to-end with a single stochastic variational inference objective stated in Equation (3). Such a variational bound includes the reconstruction loss for x1:C and the KL term for z1:C . Alternatively, we can exclude these terms from the objective, which is equivalent to starting the inference process from t = C + 1 instead of t = 1.
Our models incur a time complexity of O(T 2d) and a memory cost of O(T 2d), where T is the total sequence length and d is the dimensionality of the latent space. The recursive latent dynamics also does not allow use the take full advantange of parallelizable attentions. However, we find that our models are still efficient in practice, especially for reasonably small values of T .

heading: 3.2 Multi-Layered Extension for Probabilistic Transformer
text: Inspired by recent work on hierarchical VAEs for non-sequential inputs [17, 80, 83, 101], we extend our proposed model to include several layers of latent variables, aiming to further increase its flexibility for modelling sequential data.
We represent each time step t with a Ma rkov chain of L latent variables z(1:L)t = (z (1) t , . . . , z (L) t ) for simplicity (see Figure 1). The generative and inference model also decompose auto-regressively across different time steps and may exhibit non-Markovian dynamics:
pθ ( x1:T , z (1:L) 1:T |x1:C ) = ( T∏ t=1 pθ ( xt | z(L)t ))( L∏ `=1 T∏ t=1 pθ ( z (`) t | z (`) 1:t−1, z (`−1) 1:T ,x1:C )) (12)
qφ ( z (1:L) 1:T |x1:T ) = L∏ `=1 T∏ t=1 qφ ( z (`) t | z (`) 1:t−1, z (`−1) 1:T ,x1:T ) . (13)
Intuitively, we generate samples x1:T conditioning on x1:C by following the latent dynamics from the bottom up and using the generative process described earlier within each layer. Analogously, inference proceeds in the same order, resulting in a variational bound similar to Equation (3):
log pθ(x1:T |x1:C) ≥ T∑ t=1 Eq [ log pθ(x (L) t | zt) ] (14)
− L∑ `=1 KL(qφ(z (`) t | z (`) 1:t−1, z (`) 1:T ,x1:T ) ‖ pθ(z (`) t | z (`) 1:t−1, z (`) 1:T ,x1:C)). (15)
As before, we parametrize the prior pθ(z (`) t | z (`) 1:t−1, z (`) 1:T ,x1:C) using self-attention over the inferred latent variables from previous time steps w(`)t−1 on the same layer and another attention over contexts h1:C . In this case, however, we include an additional self-attention over all latent variables from the layer immediately below it (see Equation (16)):
w̃ (`) t = LayerNorm(w (`) t−1 + Attention(w (`) t−1,w (`−1) 1:T ,w (`−1) 1:T )) (16)
w̄ (`) t = LayerNorm(w̃ (`) t + Attention(w̃ (`) t ,w (`) 1:t−1,w (`) 1:t−1)) (17)
ŵ (`) t = LayerNorm(w̄ (`) t + Attention(w̄ (`) t ,h1:C ,h1:C)) (18)
z (`) t = Sample(N (z (`) t ;MLP(ŵ (`) t ), Softplus(MLP(ŵ (`) t )))) (19)
w (`) t = LayerNorm(ŵ (`) t + MLP(z (`) t ) + Position(t)), (20)
Stacking multiple layers of latent variables increases model expressiveness, but it also result in a linear increase in running time and the number of parameters. The time complexity for the Llayers transformer is O(LT 2d), while the space complexity remains D(T 2d) due to the Markovian structure of the chain z(1:L)t at each time step t. In our experiments, we restrict the number of layers of our hierachical models to two or three.

heading: 4 Related Work
text: Deep State Space Models. Deep neural networks have been extensively combined with state space models, resulting in flexible, yet principledly motivated latent variable approaches. While some work keep the linear state transition intact to leverage the efficient Kalman filer algorithms [23, 31, 47, 71], more expressive, nonlinear latent dynamics parametrized by neural networks have been proposed [51, 52]. All such models are limited to the Markovian dynamics of LDSs, which hinders learning of long-range dependencies. The limitation is often alleviated by combining the stochastic transitions with a deterministic RNN that enables access to all past states [3, 8, 22, 30, 39, 77]. Our models are similarly non-Markovian, but the dependencies on the past states are done via attention, which allows for easy connections between long-distance time steps. In addition, while most existing deep SSMs represent each time step with a single latent variable, our models include several layers of hierarchical latent variables with tractable inference mechanism.
Attentive Recurrent Networks. Attention mechanism has also been widely adopted in recent time series work using sequence-to-sequence models [1, 28] or transformer architectures [14, 55, 57, 72, 81, 94]. While our models are equipped with latent variables, these transformer approaches [55, 72] lack inference mechanism and are susceptible to feeding back observation noise into the dynamics model at test time. Our work, however, can be considered as an extension of the attentive state space model proposed in [1], with discrete latent states replaced by their continuous analogs. Recent developments in natural language processing [58, 60, 90] also combine transformer and VAE; however, these approaches often use a time-agnostic latent variable, in contrast to our SSM formulation.
Time Series Forecasting. Traditional univariate time series models, such as Box-Jenkins methods [12] and exponential smoothing [43], often assume independence between any collection of time series [76]. While multivariate extensions of the classical approaches, including vector autoregression [82] and multivariate GARCH [7], do not require such a strong assumption, they come with many others such as stationarity and homocesdasticity, demand manual selection of covariates and models, and do not scale well to even a moderate number of time series [40, 69].
Deep learning methods for time series forecasting have recently emerged as an expressive, scalable framework for industrial applications [10, 68, 79, 91]. While early work focus on point forecasts [53, 70, 96], recent approaches often employ recurrent neural networks with probabilistic forecasts parametrized directly [76], using quantile functions [33], Gaussian copulas [75], normalizing flows [23], or diffusion models [73]. In contrast, our models are entirely devoid of such recurrent architectures and rely on latent variables to output distributional forecasts.
Human Motion Prediction. Despite being almost identical in formulation, human motion prediction has often been studied independently from time series forecasting. While some work deterministically generate future motions or video frames [13, 32, 34, 56], stochastic prediction has also been proposed with deep neural networks often outperforming traditional methods such as hidden Markov models [93] or Gaussian processes [89] on complex motion datasets [13, 32, 45, 54, 63]. In contrast to earlier work [95, 97] that employ a global latent variable across different time steps via conditional VAE [49], we leverage the principled framework of state space models for learning and inference of hierarchical, time-dependent latent variables.

heading: 5 Experiments
text: We present our experiment results on two tasks, namely, time series forecasting and human motion prediction. These tasks are often studied independently, despite being almost identical as conditional prediction problems.

heading: 5.1 Time-series Forecasting
text: Datasets & Covariates. Following the experiment setup in [72, 73, 75], we evaluate our models and multiple competitive baselines on five popular public datasets: SOLAR, ELECTRICITY, TRAFFIC, TAXI, and WIKIPEDIA. The data is recorded with hourly or daily frequency and shows seasonal patterns of different frequencies (see Appendix A for more dataset details). As in [72, 73], the covariates include lagged inputs, fixed time embeddings (e.g. day of week, hour of day), and learnable time-series embeddings. The inputs are scaled using the conditioning examples before being fed into the model, and the predictions are rescaled appropriately afterward.
Metrics. Following [23, 72, 75], we evaluate our model and all baselines using continuous ranked probability score (CRPS) [65] summed across time series, denoted by CRPSsum. Given a univariate distribution function F and an observation x, CRPS is defined as
CRPS(F ,x) = ∫ R (F (z)− 1{x≤z})2dz,
where 1{x≤z} is the indicator function. As argued in de Bézenac et al. [23], CRPSsum is a proper scoring rule [35] and can be computed without analytical forecast distributions. We compute the metrics in a rolling fashion and use 100 samples for the distributional forecasts, similar to the aforementioned work.
Baselines. We benchmark our models against various baselines, including (1) VES [43], an innovation state space model; (2) VAR-Lasso and VAR [61], two multivariate linear autoregressive models with and without Lasso regularization; (3) GARCH [84], a multivariate conditional heteroskedastic model; (4) DeepAR [76], an autoregressive recurrent neural network; LSTM-Copula and GP-Copula [75], two RNN-based models that use Gaussian copula to model nonlinearity; (5) KVAE [51], a variational approach based on linear dynamics; (6) NKF [23], a normalizing-flow model coupled with Kalman filters; (7) Transformer [72], a transformer-based model based on masked autoregressive flow; and (8) TimeGrad [73], a recent autoregressive approach that uses a diffusion model.
Implementations. We use 8-head attentions and 2-layers MLPs to parametrize the generative and inference models. The stochastic latent variables zt are 16-dimensional while the hidden representations wt are in R128. Our probabilistic transformers for SOLAR and ELECTRICITY have one stochastic layer while those for the other datasets of higher dimensional observations employ two layers. We re-
port the numbers of parameters of our models in Table 4 in Appendix C, which are all comparable to those of the state-of-the-art approaches. See Appendix D for more details about hyper-parameters and training processes.
Accuracy Comparison. Table 1 shows that our models perform competitively across all five highdimensional time series datasets, achieving CRPSsum comparable to the best methods on ELECTRICITY and WIKIPEDIA while outperforming all baselines, including a non-SSM transformer-based approach [72], by significant margins on SOLAR, TRAFFIC and TAXI. Further analyses with other metrics, including CRPS and MSE, in Appendix B also help confirm our findings.
Qualitative Results. Figure 2 shows that the distribution forecasts generated by our model follow closely the ground truths, which is consistent with our accuracy results. In addition, the model appears to capture the uncertainty of future forecasts to some extent; observations of large magnitudes and far into the future seem to correctly have higher variance estimates.
Ablation Study. We include a small scale ablation study on the TRAFFIC dataset to investigate which components of our models are essential. Table 2 suggests that removing the stochasticity from wt has most impacts on model performance, implying that incorprating latent variables into a transformer is indeed useful. Other aspects such as context attention or multiple layers of stochastic variables do not show dramatic effects in this study; however, they do contribute performance gains.

heading: 5.2 Human Motion Prediction
text: Datasets. Following the experiment setup in [97], we evaluation our models on two public motion capture datasets: Human3.6M[44] and HumanEva-I [78]. While Human3.6 is a large-scale dataset with 3.6 million video frames recorded at 50Hz, HumanEva-I is smaller with only 3 subjects and recorded at 60Hz. We follow the preprocessing steps of previous work [64, 97] and obtain a 17-joint skeleton for Human3.6 and a 15-joint skeleton for HumanEva-I. As in [97], we predict future motion for 2 seconds conditioning on observed motion of 0.5 seconds and 1 second conditioning on 0.25 seconds for Human3.6 and HumanEva-I, respectively.
Metrics. Following previous work on trajectory forecasting [2, 37], we adopt two popular metrics, namely, average displacement error (ADE) and final displacement error (FDE). ADE measures the average L2 distance over all time steps between the ground truth motion and the closest sample, while FDE only consider such distance for the final pose.
Baselines. We compare our models against 9 models, including ERD [32] and acLSTM [56], two deterministic RNN-based approaches; MT-VAE [95] and Pose-Knows [87], two conditional VAE models; HP-GAN [6], a conditional GAN; Best-Many [11], GMVAE [25], DeliGAN [38]. and DSP [98], four approaches optimizing for diversity objectives. The results for these baselines are reported as in [97].
Implementations. Similar to the previous experiments, we use 8-head attentions and 2-layers MLPs. Since Human3.6M is significantly more complex and multi-modal than the time series forecasting datasets, we make use of 3 stochastic layers, as opposed to 2 layers for HumanEva-I. For Human3.6M, the context and target observations are significantly longer and set up for long-term predictions, so we only infer latent variables for target observations. Appendix C also contains further details about our models and their number of parameters.
Quantitative Results. Table 3 shows that our models convincingly outperform all baselines based on both metrics ADE and FDE, with the gains significantly higher for the larger dataset Human3.6M. We emphasize that our favorable performance is evaluated using random samples, while the closest competitor, DLow [97], relies on a separate model for selecting samples to promote diversity, which can potentially be combined with our probabilistic transformer for further improvements.
Qualitative Results. We show in Figure 3 human pose predictions made by our model that are most similar to the corresponding ground truths among a collection of such stochastic predictions. The similarities between the body-part movements in both sequences suggest that our model has been able to capture the temporal dynamics quite well.

heading: 6 Conclusion & Discussion
text: In this work, we have introduced generative models for multivariate time series that combines strengths of state space models and transformer architectures. In contrast to previous work, our models do not rely on recurrent neural networks but make extensive use of attention mechanism. We also extend our models to include hierarchical latent variables, inspired by recent developments of VAEs for non-sequential data [17, 83]. Empirical experiments show that our models perform remarkably well on time series forecasting and human motion prediction.
Our models do not come without limitations, however. As in other transformer-based approaches, the reliance on attention incurs a quadratic time and memory complexity. While we do not find it problematic in our experiments, the limitation necessarily hinders applications of our models in tasks characterized by long-term dependencies such as language modelling or music generation [36]. Fortunately, recent work on sparse transformer [9, 18, 50, 55] can potentially address the issue, and we leave such an investigation for future work.
Probabilistic time series forecasting is a fundamental research problem with wide-ranging applications in society. Although we have not explored healthcare applications of our work, previously proposed methods with similar formulations have demonstrated potentials of forecasting techniques [1, 81] in diagnoses or disease control.

heading: 7 Acknowledgment
text: Financial support is gratefully acknowledged from a Xerox PARC Faculty Research Award, National Science Foundation Awards 1455172, 1934985, 1940124, and 1940276, USAID, and Cornell University Atkinson Center for a Sustainable Future.
,

Review summaries:
## Summary The authors propose a probabilistic generative non-autoregressive transformer based state space model (SSMs) and corresponding inference procedures for them. Starting with a detailed description of the related probabilities within their models. Followed by a description of the architecture split up into Generative Model and Inference Model also serving a Multi-Layered Extension. In related work several alternative forecasting method types are presented. They compare their methods with 11 of these methods.

The submission presents a probabilistic extension to the Transformer by augmenting the intermediate Transformer representations via stochastic latent variables. The attention mechanism allows for accessing the past and the future steps directly in contrast to the indirect hidden state of RNN-based state-space models (SSMs). The building blocks of the Transformer remain untouched. The attention operation is applied to the deterministic representations as in the Transformer. The random variables are then calculated by using these updated deterministic representations. The authors also propose a hierarchical variant of their model to enable stacking several layers as commonly observed in the Transformer architectures. 

The paper proposes a probabilistic forecasting model combining transformers with state space models in a non-trivial fashion. It uses variational inference for inference of the multiple, hierarchically stacked state variables (instead of other related works that assume linear state transitions and can then use kalman filtering). The empirical evaluations show the practical relevance of this approach.