<figure><img src="image_1.png" /></figure>

HAL
open science

# Dynamical Variational Autoencoders: A Comprehensive Review

Laurent Girin, Simon Leglaive, Xiaoyu Bie, Julien Diard, Thomas Hueber, Xavier Alameda-Pineda

► To cite this version:

Laurent Girin, Simon Leglaive, Xiaoyu Bie, Julien Diard, Thomas Hueber, et al. Dynamical Variational Autoencoders: A Comprehensive Review. Foundations and Trends in Machine Learning, 2021, 15 (1-2), pp.1-175. 10.1561/2200000089. hal-02926215v2

HAL Id: hal-02926215

https://inria.hal.science/hal-02926215v2

Submitted on 5 Jul 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.
The version of record is available at:
http://dx.doi.org/10.1561/2200000089
# Dynamical Variational Autoencoders: A Comprehensive Review

Laurent Girin<fnref n="1" />, Simon Leglaive<fnref n="2" />, Xiaoyu Bie<fnref n="3" />, Julien Diard<fnref n="4" />,
Thomas Hueber<fnref n="1" /> and Xavier Alameda-Pineda<fnref n="3" />

## ABSTRACT

Variational autoencoders (VAEs) are powerful deep generative models widely used to represent high-dimensional complex data through a low-dimensional latent space learned in an unsupervised manner. In the original VAE model, the input data vectors are processed independently. Recently, a series of papers have presented different extensions of the VAE to process sequential data, which model not only the latent space but also the temporal dependencies within a sequence of data vectors and corresponding latent vectors, relying on recurrent neural networks or state-space models. In this paper, we perform a literature review of these models. We introduce and discuss a general class of models, called dynamical variational autoencoders (DVAEs), which

<footnote n="1">¹Univ. Grenoble Alpes, CNRS, Grenoble-INP, GIPSA-lab;
laurent.girin@grenoble-inp.fr, thomas.hueber@grenoble-inp.fr</footnote>

<footnote n="2">²CentraleSupélec, IETR; simon.leglaive@centralesupelec.fr</footnote>

<footnote n="3">³Inria, Univ. Grenoble Alpes, CNRS, LJK; xiaoyu.bie@inria.fr,
xavier.alameda-pineda@inria.fr</footnote>

<footnote n="4">⁴Univ. Grenoble Alpes, CNRS, LPNC;
julien.diard@univ-grenoble-alpes.fr</footnote>
encompasses a large subset of these temporal VAE extensions. Then, we present in detail seven recently proposed DVAE models, with an aim to homogenize the notations and presentation lines, as well as to relate these models with existing classical temporal models. We have reimplemented those seven DVAE models and present the results of an experimental benchmark conducted on the speech analysis-resynthesis task (the PyTorch code is made publicly available). The paper concludes with a discussion on important issues concerning the DVAE class of models and future research guidelines.
# Contents

<table><tr><td>1</td><td>Introduction</td><td>7</td></tr><tr><td>1.1</td><td>Deep Dynamical Bayesian Networks</td><td>8</td></tr><tr><td>1.2</td><td>Variational inference and VAEs</td><td>10</td></tr><tr><td>1.3</td><td>Dynamical VAEs</td><td>12</td></tr><tr><td>1.4</td><td>Aim, contributions, and outline of the paper</td><td>14</td></tr><tr><td>2</td><td>Variational Autoencoders</td><td>19</td></tr><tr><td>2.1</td><td>Principle</td><td>19</td></tr><tr><td>2.2</td><td>VAE generative model</td><td>20</td></tr><tr><td>2.3</td><td>Learning with variational inference</td><td>23</td></tr><tr><td>2.4</td><td>VAE inference model</td><td>28</td></tr><tr><td>2.5</td><td>VAE training</td><td>28</td></tr><tr><td>3</td><td>Recurrent Neural Networks and State Space Models</td><td>31</td></tr><tr><td>3.1</td><td>Recurrent Neural Networks</td><td>31</td></tr><tr><td>3.2</td><td>State Space Models</td><td>34</td></tr><tr><td>4</td><td>Definition of Dynamical VAEs</td><td>39</td></tr><tr><td>4.1</td><td>Generative model</td><td>39</td></tr><tr><td>4.2</td><td>Inference model</td><td>46</td></tr><tr><td>4.3</td><td>VLB and training of DVAEs</td><td>50</td></tr><tr><td>4.4</td><td>Additional dichotomy for autoregressive DVAE models</td><td>52</td></tr></table>
<table>
  <tbody>
    <tr>
      <td>4.5</td>
      <td>DVAE summary</td>
      <td>54</td>
    </tr>
    <tr>
      <td><b>5</b></td>
      <td><b>Deep Kalman Filters</b></td>
      <td><b>57</b></td>
    </tr>
    <tr>
      <td>5.1</td>
      <td>Generative model</td>
      <td>58</td>
    </tr>
    <tr>
      <td>5.2</td>
      <td>Inference model</td>
      <td>59</td>
    </tr>
    <tr>
      <td>5.3</td>
      <td>Training</td>
      <td>62</td>
    </tr>
    <tr>
      <td><b>6</b></td>
      <td><b>Kalman Variational Autoencoders</b></td>
      <td><b>65</b></td>
    </tr>
    <tr>
      <td>6.1</td>
      <td>Generative model</td>
      <td>66</td>
    </tr>
    <tr>
      <td>6.2</td>
      <td>Inference model</td>
      <td>68</td>
    </tr>
    <tr>
      <td>6.3</td>
      <td>Training</td>
      <td>70</td>
    </tr>
    <tr>
      <td><b>7</b></td>
      <td><b>STOchastic Recurrent Networks</b></td>
      <td><b>71</b></td>
    </tr>
    <tr>
      <td>7.1</td>
      <td>Generative model</td>
      <td>71</td>
    </tr>
    <tr>
      <td>7.2</td>
      <td>Inference model</td>
      <td>74</td>
    </tr>
    <tr>
      <td>7.3</td>
      <td>Training</td>
      <td>76</td>
    </tr>
    <tr>
      <td><b>8</b></td>
      <td><b>Variational Recurrent Neural Networks</b></td>
      <td><b>79</b></td>
    </tr>
    <tr>
      <td>8.1</td>
      <td>Generative model</td>
      <td>79</td>
    </tr>
    <tr>
      <td>8.2</td>
      <td>Inference model</td>
      <td>81</td>
    </tr>
    <tr>
      <td>8.3</td>
      <td>Training</td>
      <td>82</td>
    </tr>
    <tr>
      <td>8.4</td>
      <td>Improved VRNN and VRNN applications</td>
      <td>83</td>
    </tr>
    <tr>
      <td><b>9</b></td>
      <td><b>Stochastic Recurrent Neural Networks</b></td>
      <td><b>85</b></td>
    </tr>
    <tr>
      <td>9.1</td>
      <td>Generative model</td>
      <td>85</td>
    </tr>
    <tr>
      <td>9.2</td>
      <td>Inference model</td>
      <td>87</td>
    </tr>
    <tr>
      <td>9.3</td>
      <td>Training</td>
      <td>89</td>
    </tr>
    <tr>
      <td><b>10</b></td>
      <td><b>Recurrent Variational Autoencoders</b></td>
      <td><b>91</b></td>
    </tr>
    <tr>
      <td>10.1</td>
      <td>Generative model</td>
      <td>91</td>
    </tr>
    <tr>
      <td>10.2</td>
      <td>Inference model</td>
      <td>95</td>
    </tr>
    <tr>
      <td>10.3</td>
      <td>Training</td>
      <td>98</td>
    </tr>
    <tr>
      <td><b>11</b></td>
      <td><b>Disentangled Sequential Autoencoders</b></td>
      <td><b>99</b></td>
    </tr>
    <tr>
      <td>11.1</td>
      <td>Generative model</td>
      <td>99</td>
    </tr>
    <tr>
      <td>11.2</td>
      <td>Inference model</td>
      <td>100</td>
    </tr>
    <tr>
      <td>11.3</td>
      <td>Training</td>
      <td>102</td>
    </tr>
  </tbody>
</table>
<table>
  <tbody>
    <tr>
      <td><strong>12 Brief tour of other models</strong></td>
      <td><strong>105</strong></td>
    </tr>
    <tr>
      <td>12.1 Models related to DKF</td>
      <td>105</td>
    </tr>
    <tr>
      <td>12.2 Models related to STORN, VRNN, and SRNN</td>
      <td>107</td>
    </tr>
    <tr>
      <td>12.3 Other models</td>
      <td>112</td>
    </tr>
    <tr>
      <td><strong>13 Experiments</strong></td>
      <td><strong>115</strong></td>
    </tr>
    <tr>
      <td>13.1 DVAE architectures</td>
      <td>115</td>
    </tr>
    <tr>
      <td>13.2 Experimental protocol</td>
      <td>118</td>
    </tr>
    <tr>
      <td>13.3 Results on speech data</td>
      <td>121</td>
    </tr>
    <tr>
      <td>13.4 Results on 3D human motion data</td>
      <td>133</td>
    </tr>
    <tr>
      <td>13.5 Conclusion</td>
      <td>136</td>
    </tr>
    <tr>
      <td><strong>14 Discussion</strong></td>
      <td><strong>139</strong></td>
    </tr>
    <tr>
      <td>14.1 Fundamental motivation for DVAEs</td>
      <td>139</td>
    </tr>
    <tr>
      <td>14.2 DVAE outcome: A story of flexibility</td>
      <td>140</td>
    </tr>
    <tr>
      <td>14.3 VAE improvements and extensions applicable to DVAEs</td>
      <td>144</td>
    </tr>
    <tr>
      <td>14.4 Perspectives on source coding</td>
      <td>155</td>
    </tr>
    <tr>
      <td><strong>Acknowledgements</strong></td>
      <td><strong>159</strong></td>
    </tr>
    <tr>
      <td><strong>Appendices</strong></td>
      <td><strong>161</strong></td>
    </tr>
    <tr>
      <td><strong>A Marginalization of h<sub>1:T</sub> in STORN</strong></td>
      <td><strong>163</strong></td>
    </tr>
    <tr>
      <td><strong>B DVAE implementation with speech data</strong></td>
      <td><strong>165</strong></td>
    </tr>
    <tr>
      <td>B.1 DKF</td>
      <td>166</td>
    </tr>
    <tr>
      <td>B.2 STORN</td>
      <td>166</td>
    </tr>
    <tr>
      <td>B.3 VRNN</td>
      <td>167</td>
    </tr>
    <tr>
      <td>B.4 SRNN</td>
      <td>167</td>
    </tr>
    <tr>
      <td>B.5 RVAE</td>
      <td>168</td>
    </tr>
    <tr>
      <td>B.6 DSAE</td>
      <td>168</td>
    </tr>
    <tr>
      <td><strong>C DVAE implementation with 3D human motion data</strong></td>
      <td><strong>171</strong></td>
    </tr>
    <tr>
      <td>C.1 DKF</td>
      <td>171</td>
    </tr>
    <tr>
      <td>C.2 STORN</td>
      <td>172</td>
    </tr>
    <tr>
      <td>C.3 VRNN</td>
      <td>172</td>
    </tr>
  </tbody>
</table>
<table><tr><td>C.4 SRNN</td><td>172</td></tr><tr><td>C.5 RVAE</td><td>173</td></tr><tr><td>C.6 DSAE</td><td>173</td></tr><tr><td><strong>References</strong></td><td><strong>175</strong></td></tr></table>
1

Introduction

Deep Generative Models (DGMs) constitute a large family of probabilis-
tic models that are currently of high interest in the machine learning and
signal processing. They result from the combination of conventional (i.e.,
non-deep) generative probabilistic models and Deep Neural Networks
(DNNs). For both conventional models and DGMs, different noncon-
flicting taxonomies can be established due to the domain richness and
percolation across the different approaches. Nevertheless, these models
can be grossly classified into the following two categories. Using the
terminology of Diggle and Gratton (1984), the first category corresponds
to *prescribed* models for which the probability density function (pdf)
of the generative model is defined explicitly, generally through a para-
metric form. The second category corresponds to *implicit* models that
can generate data “directly,” without using an explicit formulation and
manipulation of a pdf model. Generative adversarial networks (GANs)
are a popular example of this second category (Goodfellow *et al.*, 2014;
Goodfellow *et al.*, 2016; Goodfellow, 2016).

In the present review, we focus on the first category, in which a
parametric pdf model is used. A suitable feature of generative models
based on an explicit formulation of the pdf is that they can be easily
plugged into a more general Bayesian framework, not only for gener-
ating data but also for modeling the data structure (without actually
generating them) in various applications (e.g., data denoising or data
transformation). In any case, the pdf model must be as close as possible
to the true pdf of the data, which is generally unknown. To achieve this
aim, the model must be trained from data, and model parameters are
generally estimated by following the maximum likelihood methodology
(Goodfellow et al., 2016; Bishop, 2006; Koller and Friedman, 2009).
These principles are valid for both conventional generative models and
DGMs; however, in the case of DGMs, the pdf parameters are generally
the output of DNNs, which makes model training potentially difficult.

## 1.1 Deep Dynamical Bayesian Networks

In the present review, we focus on an important subfamily of DGMs,
namely the deep dynamical Bayesian networks (DDBNs), which are
built on the following models:

* Bayesian networks (BNs) are a popular class of probabilistic models for which i) the dependencies among all involved random variables are explicitly represented by conditional pdfs (i.e. BNs are prescribed models), and ii) these dependencies can be schematically represented using a directed acyclic graph (Bishop, 2006; Koller and Friedman, 2009). The structure of these dependencies often reflects (or originates from) an underlying hierarchical generative process.

* Dynamical Bayesian networks are BNs that include temporal dependencies and are widely used to model dynamical systems and/or data sequences. Dynamical BNs are BNs “repeated over time”; that is, they exhibit a repeating dependency structure (a time-slice at discrete time t) and some dependencies across these time-slices (the dynamical models). Recurrent neural networks (RNNs) and state-space models (SSMs) can be considered special cases of dynamical BNs. In fact, a temporal dependency in a dynamical BN is often implemented either as a deterministic recursive process, as in RNNs, or as a first-order Markovian process,
<figure><img src="image_2.png" /><figcaption>**Figure 1.1:** A graphical taxonomy of generative probabilistic models. Only the branch corresponding to the models covered in this review is detailed. Nodes represent classes of models, and arrow labels specify some of the relationships between the classes of models. Please refer to the text for details and acronym definitions. The numbered gray circles indicate the section number in which the corresponding class of models is detailed.</figcaption></figure>
as in usual SSMs.

* Deep Bayesian networks combine BNs with DNNs. DNNs are used to generate the parameters of the modeled distributions. This enables them to be high-dimensional and highly multi-modal while having a reasonable number of parameters. In short, deep BNs have can appropriately combine the “explainability” of Bayesian models with the modeling power of DNNs.

DDBNs are thus a combination of all these aspects, as illustrated in Figure 1.1. They can be equally seen as dynamical versions of deep BNs (i.e., deep BNs including temporal dependencies) or deep versions of dynamical BNs (i.e., dynamical BNs mixed with DNNs). As an extension of dynamical BNs, DDBNs are expected to be powerful tools for modeling dynamical systems and/or data sequences. However, as mentioned above, the combination of probabilistic modeling with DNNs in deep BNs can result in a complex and costly model training. This is an even more serious issue for DDBNs, in which the repeating structure due to temporal modeling adds a level of complexity.

## 1.2 Variational inference and VAEs

Recently, the application of the variational inference methodology (Jordan et al., 1999; Bishop, 2006, Chapter 10; Šmídl and Quinn, 2006; Murphy, 2012, Chapter 21) to a fundamental deep BN architecture –a low-dimensional to a high-dimensional generative feed-forward DNN– has led to efficient inference and training of the resulting model, called a variational autoencoder (VAE) (Kingma and Welling, 2014). A similar approach was proposed the same year (Rezende et al., 2014).<fnref n="1" /> The VAE is directly connected to the concepts of a *latent variable* and *unsupervised representation learning*: the observed random variable representing the data of interest is assumed to be generated from an unobserved latent variable through a probabilistic process. Often, this latent variable is of lower dimension than the observed data (which can be high-dimensional) and is assumed to “encode” the observed data in

<footnote n="1">¹Kingma and Welling (2014) and Rezende et al. (2014) were both pre-published in 2013 as ArXiv papers. Connections also exist with Mnih and Gregor (2014).</footnote>
a compact manner so that new data can be generated from new values of the latent variable. Moreover, one wishes to extract a latent representation that is *disentangled* (i.e., different latent coefficients encode different properties or different factors of variation in the data). When successful, this provides good interpretability and control of the data generation/transformation process.

The automatic discovery of a latent space structure is part of the model training process. The inference process, which is defined in the present context as the estimation of latent variables from the observed data, also plays a major role. As presented in detail later, in a deep BN, the exact posterior distribution (i.e., the posterior distribution of the latent variable given the observed variable corresponding to the generative model) is generally not tractable. It is thus replaced with a parametric approximate posterior distribution (i.e., an inference model) that is implemented with a DNN. As the observed data likelihood function is also not tractable, the model parameters are estimated by chaining the inference model (also known as the *encoder* in the VAE framework) and the generative model (the *decoder*) and maximizing a lower bound of the log-likelihood function, called the variational lower bound (VLB), over a training dataset.<fnref n="2" /> Hereinafter, we refer to this general variational inference and training methodology as *the VAE methodology*.

In summary, the VAE methodology enables deep unsupervised representation learning while providing efficient inference and parameter estimation in a Bayesian framework. As a result, the seminal papers by Kingma and Welling (2014) and Rezende et al. (2014) have had and continue to have a strong impact on the machine learning community. VAEs have been applied to many signal processing problems, such as the generation and transformation of images and speech signals (we provide a few references in Chapter 2).

<footnote n="2">²The idea of using an artificial neural network to approximate an inference model and chaining the encoder and decoder dates back to the early studies of Hinton *et al.* (1995) and Dayan *et al.* (1995). However, the algorithms presented in these papers for model training are different from the one used to optimize the VAE.</footnote>
## 1.3 Dynamical VAEs

As a deep BN, the original VAE proposed by Kingma and Welling (2014) did not include temporal modeling. This means that each data vector was processed independently of the other data vectors (and the corresponding latent vector was also processed independently of the other latent vectors). This is clearly suboptimal for the modeling of correlated (temporal) vector sequences.

In the years following the publication of Kingma and Welling (2014) and Rezende et al. (2014), the VAE methodology was extended and successfully applied to several more complex deep BNs. In particular, it was applied to deep BNs with a temporal model (i.e., DDBNs) dedicated to the modeling of sequential data exhibiting temporal correlation. In the present review, we are particularly interested in the models presented in the following papers: (Bayer and Osendorfer, 2014; Krishnan et al., 2015; Chung et al., 2015; Gu et al., 2015; Fraccaro et al., 2016; Krishnan et al., 2017; Fraccaro et al., 2017; Goyal et al., 2017; Hsu et al., 2017b; Li and Mandt, 2018; Leglaive et al., 2020). In addition to including temporal dependencies, the unsupervised representation learning essence of the VAE is preserved and cherished in these studies. These DDBNs combine the observed and latent variables and aim at modeling not only data dynamics but also discovering the latent factors governing them.

To achieve this aim, these models are trained using the VAE methodology (i.e., design of an inference model and maximization of the corresponding VLB). We can thus encompass these models under the common class and terminology of *variational DDBNs* (i.e., DDBNs immersed in the VAE framework). In the following of the paper, as well as the title, we prefer to refer to them as *dynamical VAEs* (DVAEs) (i.e., VAEs including a temporal model for modeling sequential data). This is simply because we assume that the term “VAE” is currently more popular than the term “DBN,” and “dynamical VAEs” gives a more speaking-first evocation of these models, compared to “variational DDBNs.” This convergence of DDBNs and VAEs into DVAEs is illustrated in Figure 1.1.

In practice, these different DVAE models vary in how they define the dependencies between the observed and latent variables, how they
define and parameterize the corresponding generative pdfs, and how they
define and parameterize the inference model. They also differ in how
they combine the variables with RNNs to model temporal dependencies,
at both generation and inference. In contrast, they are all characterized
by the following common set of features.

First, as stated above, they are all trained using the VAE methodol-
ogy, possibly with a few adaptations and refinements. In this paper, we
do not review models based on GANs and, more generally, on adversarial
training. Examples of extensions of “static” GANs to sequence modeling
and generation can be found in the literature (Mathieu et al., 2016;
Villegas et al., 2017; Denton and Birodkar, 2017; Tulyakov et al., 2018;
Lee et al., 2018). This approach is particularly popular for separating
content and motion in videos.

Second, even if the observed random vectors can be continuous or
discrete, as in the original VAE formulation, they all feature continuous
latent random variables. In the present review, we do not consider the
case of discrete latent random variables. The latter can be incorporated
in DVAE models, in the line with the case of, for example, conditional
VAEs (Sohn et al., 2015; Zhao et al., 2017). Temporal models with binary
observed and latent random variables have been proposed (Boulanger-
Lewandowski et al., 2012; Gan et al., 2015). These models are based
on restricted Boltzman machines (RBMs) or sigmoid belief networks
(SBNs) combined with RNNs. A detailed analysis of such models is
beyond the scope of the present review.

Third, all DVAE models we consider feature a *discrete-time sequence of (continuous or discrete) observed random vectors associated with a corresponding discrete-time sequence of (continuous) latent random vectors*. In other words, these models function in a *sequence-to-sequence mode for both encoding and decoding*. Thus, we do not focus on VAE-based models specifically designed for text and dialogue generation (Bowman *et al.*, 2016; Miao *et al.*, 2016; Serban *et al.*, 2016; Serban *et al.*, 2017; Yang *et al.*, 2017; Semeniuta *et al.*, 2017; Hu *et al.*, 2017; Zhao *et al.*, 2018; Jang *et al.*, 2019) or (2D) image modeling (Gulrajani *et al.*, 2016; Chen *et al.*, 2017; Lucas and Verbeek, 2018; Shang *et al.*, 2018). These models generally have a many-to-one encoder and a one-to-many decoder; that is, a long sequence of data (e.g., words or pixels) is
encoded into a single latent vector, which is in turn decoded into a whole
data sequence (see also Roberts et al. (2018) and Pereira and Silveira
(2018) for examples on music score modeling and anomaly detection
in energy time series, respectively). Even if those models can include
a hierarchical structure at encoding and/or at decoding, they do not
consider a temporal sequence of latent vectors.

All these latent-variable deep temporal models, the ones we detail and unify in the DVAE class, and the ones we do not detail, remain strongly connected, with a similar overall encoding-decoding architecture and possibly a similar inference and training VAE methodology. Therefore, we must keep in mind that some of the propositions made in the literature for one type of model can be adapted and be beneficial to the other.

1.4 Aim, contributions, and outline of the paper

This paper aims to provide a comprehensive overview of DVAE models.
The contributions of this paper are detailed as follows.

**We provide a formal definition of the general class of DVAEs.**
We describe its main properties and characteristics and how this class is
related to previous classical models, such as VAEs, RNNs, and SSMs. We
discuss the structure of dependencies between the observed and latent
random variables in DVAE pdfs, as well as how these dependencies are
implemented with neural networks. We discuss the design of inference
models considering the general methodology used to identify the actual
dependencies of the latent variables at inference time. We also discuss
the VLB computation for training DVAEs. All these points are presented
in Chapter 4. To the best of our knowledge, this is the first time this
class of models has been presented in such a general and unified manner.

**We provide a detailed and complete technical description of seven DVAE models selected from the literature.** In Chapter 5, we start with the deep Kalman filter (DKF) (Krishnan et al., 2015; Krishnan et al., 2017), which is a basic combination of an SSM with DNNs. Then, we examine the Kalman variational autoencoder (KVAE) (Fraccaro et al., 2017) in Chapter 6, the stochastic recurrent neural
network (STORN) (Bayer and Osendorfer, 2014) in Chapter 7, the variational recurrent neural network (VRNN) (Chung et al., 2015; Goyal et al., 2017) in Chapter 8, another type of stochastic recurrent neural network (SRNN) (Fraccaro et al., 2016) in Chapter 9, the recurrent variational autoencoder (RVAE) (Leglaive et al., 2020) in Chapter 10, and finally the disentangled sequential autoencoder (DSAE) (Li and Mandt, 2018) in Chapter 11.

We have spent effort on consistency of presentation. For all seven models that we detail, we first present the generative equations in time-step form and then for an entire data sequence. Then, we present the structure of the exact posterior distribution of the latent variables given the observed data and present the inference model as proposed in the original papers. Finally, we present the corresponding VLB. In the original papers, some parts of this complete picture are often overviewed or even missing (not always the same parts), independently of the authors’ goodwill, because of lack of space.

We discuss the links, similarities, and differences of the selected DVAE models. We comment on the choices of the authors of the reviewed papers regarding the inference model, its relation to the exact posterior distribution, and implementation issues. In the present review, we discuss only high-level implementation issues related to the general structure of the neural network that implements a given DVAE at generation or inference (e.g., the type of RNN). We do not discuss practical implementation issues (e.g., the number of layers), which are too low-level in the present technical review context.

We have also spent some effort making the notations homogeneous across all models. This is valid for both the review of the seven detailed models and the other sections of the paper, including the general presentation of the DVAE class of models in Chapter 4. For some models, we have changed the time indexation notation, and in some instances, the names of some variables compared to the original papers. We have taken great care to do that consistently in the generative model, the inference part, and the VLB so that these notation changes do not affect the essence and functioning of the model. Together with consistency of presentation, this enables us to better put in evidence the commonalities and differences across models and make their comparison easier.
Notation remarks are specified in independent dedicated paragraphs throughout the paper to facilitate connections with the original papers.

In complement to the detailed review of the seven selected models, **we provide a more rapid overview of other DVAE models presented in the recent literature in Chapter 12.**

**We relate the recent developments in DVAEs to the history and technical background of the classical models DVAEs are built on, namely VAEs, RNNs, and SSMs.** Although there are already many papers on VAEs, including tutorials, we present them in Chapter 2 because all subsequent DVAE models rely on the VAE methodology. Then, as the introduction of temporal models in the VAE framework is closely linked to RNNs and SSMs, we briefly present these two classes of models in Chapter 3. The unified notation that we use will help readers from different communities (e.g., machine learning, signal processing, and control theory) who are not familiar with the relations among VAEs, RNNs, and SSMs to discover them comfortably.

**We provide a quantitative benchmark of the selected DVAE models in an analysis-resynthesis task, as well as qualitative examples of data generation.** We have reimplemented the seven DVAE models detailed in this review and evaluated them on two different datasets (speech signals and 3D human motion data). This benchmark is presented in Chapter 13.

The performance comparison of the different models from the literature review is a difficult task for many reasons. First, all models are not evaluated on the same data. Then, a newly proposed model generally performs better than some previously proposed model(s), at least on some aspect(s), but this can depend on model tuning, task, data, and experimental setup. Moreover, the comparison performed with a subset of previous models is incomplete in essence. In short, an extended benchmark of DVAE models is not yet available in the literature. Conducting an extended benchmark is a huge endeavor, as there are many possible configurations for the models and many tasks for evaluating them. In particular, it is not yet clear how to evaluate the degree of “disentanglement” of the extracted latent space. The presented
experiments are a first step in that direction. We plan to exploit and
compare the models more extensively and on more complex tasks in
future studies. For example, we compared three models in the recently
proposed DVAE-based unsupervised speech enhancement method (Bie
et al., 2021).

The code reimplementing the seven DVAE models and
used on the benchmark task is made available to the com-
munity. A link to the open-source code and the best-trained models
can be found at https://team.inria.fr/robotlearn/dvae/. We have also
taken care, in the code, to follow the unified presentation and nota-
tion used in the paper, making it, hopefully, a useful and pedagogical
resource.

**We provide a discussion to put the DVAE class of models into perspective.** We summarize the outcome of this review and discuss the future challenges and possible improvements of VAEs and DVAEs. This is presented in Chapter 14.

In summary, we believe that comparison of models across papers is a
difficult task in essence, regardless of the efforts spent by the authors of
the original papers, because of the use of different notations, presentation
lines, missing information, etc. We hope that the present review paper
and accompanying code will enable the readers to access the technical
substance of the different DVAE models, their connections with classical
models, their cross-connections, and their unification in the DVAE class
more rapidly and “comfortably” than by analyzing and comparing the
original papers by themselves.

We wish the reader to enjoy this DVAE tour.
<figure><img src="image_3.png" /></figure>
In this section, we present the VAE and the associated methodology for model training and approximate posterior distribution estimation (i.e., inference) with variational methods (Kingma and Welling, 2014; Rezende et al., 2014). An extended tutorial on VAEs can be found in Kingma and Welling's (2019) paper.

## 2.1 Principle

For clarity of presentation, let us start with an autoencoder (AE). As illustrated in Figure 2.1, an AE is a DNN that is trained to replicate an input vector $x \in \mathbb{R}^F$ at the output (Hinton and Salakhutdinov, 2006; Vincent et al., 2010). At training time, the target output is thus set equal to $x$, and at test time, the output $\hat{x}$ is an estimated value of $x$ (i.e., we have $\hat{x} \approx x$). An AE usually has a diabolo shape. The left part of the AE, the *encoder*, provides a low-dimensional latent representation $z \in \mathbb{R}^L$ of the data vector $x$, with $L \ll F$, at the so-called bottleneck layer. The right part of the AE, the *decoder*, tries to reconstruct $x$ from $z$. So far, everything is deterministic: At test time, each time the AE is fed with a specific input vector $x_0$, it will provide the same corresponding output $\hat{x}_0$.
<figure><img src="image_4.png" /><figcaption>**Figure 2.1:** Schematic representation of an AE. The left trapezoid represents a high-to-low-dimensional encoder DNN (denoted $e_z$), and the right trapezoid represents a low-to-high-dimensional decoder DNN (denoted $d_x$). Calculation of the latent variable **z** and output $\hat{\mathbf{x}}$ from input **x** is deterministic. In line with probabilistic graphical models, deterministic variables are represented within diamonds.</figcaption></figure>

The VAE was initially proposed by Kingma and Welling (2014) and Rezende et al. (2014). It can be seen as a probabilistic version of an AE, where the output of the decoder is not directly a value of **x** but the parameters of a probability distribution of **x**. As shown below, the same probabilistic formulation applies to the encoding of **z**. The resulting probabilistic model can be used to generate new data from new values of **z**. It can also be used to transform existing data within an encoding-modification-decoding scheme. For instance, the seminal papers on VAEs and many subsequent ones have considered image generation and transformation. Examples of speech/music signals transformation based on a VAE can be found in the literature (Blaauw and Bonada, 2016; Hsu et al., 2017a; Esling et al., 2018; Roche et al., 2019; Bitton et al., 2020). Finally, it can be employed as a prior distribution of **x** in more complex Bayesian models for, for example, speech enhancement (Bando et al., 2018; Leglaive et al., 2018; Pariente et al., 2019; Leglaive et al., 2019) or source separation (Kameoka et al., 2018).

For clarity of presentation, at this point, it is convenient to separate the presentation of the VAE decoder (i.e., the generative model) and that of the VAE encoder (i.e., the inference model).

## 2.2 VAE generative model

In the following, $N(\cdot; \mu, \Sigma)$ denotes a multivariate Gaussian distribution with mean vector $\mu$ and covariance matrix $\Sigma$, $\text{diag}\{\cdot\}$ is the operator that forms a diagonal matrix from a vector by putting the vector entries
<figure><img src="image_5.png" /><figcaption>**Figure 2.2:** Schematic representation of the VAE: Encoder (left) and decoder (right). Dashed lines represent a sampling process. In line with probabilistic graphical models, latent variables are represented within empty circles and observed variables are represented within shaded circles. When the encoder and decoder are cascaded, using the same variable name x at both the input and output is an abuse of notation, but this is done to be more consistent with the separate encoder and decoder equations.</figcaption></figure>

on the diagonal, $\mathbf{0}_L$ is the zero-vector of size $L$, and $\mathbf{I}_L$ is the identity matrix of size $L$. $p_{\theta_x}(\mathbf{x})$ is a generic notation for a parametric pdf of the random variable $\mathbf{x}$, where $\theta_x$ is the set of parameters. It is equivalent to $p(\mathbf{x}; \theta_x)$.

Formally, the VAE decoder is defined by

$$p_{\theta}(\mathbf{x}, \mathbf{z}) = p_{\theta_x}(\mathbf{x}|\mathbf{z})p_{\theta_z}(\mathbf{z}), \qquad (2.1)$$

with

$$p_{\theta_z}(\mathbf{z}) = \mathcal{N}(\mathbf{z}; \mathbf{0}_L, \mathbf{I}_L), \qquad (2.2)$$

and $p_{\theta_x}(\mathbf{x}|\mathbf{z})$ is a parametric conditional distribution, the parameters of which are a nonlinear function of $\mathbf{z}$ modeled by a DNN. This DNN is called the *decoder network*, or the *generation network*, and is parametrized by a set of weights and biases denoted $\theta_x$. In the standard VAE, the set of parameters $\theta_z$ is empty, but we write it explicitly to be coherent with the rest of the paper, and we have here $\theta = \theta_x \cup \theta_z = \theta_x$. The decoder network is illustrated in Figure 2.2 (right).

The VAE model and associated variational methodology was introduced by Kingma and Welling (2014) in the general framework of parametric distributions, independently of the practical choice of the pdf $p_{\theta_x}(\mathbf{x}|\mathbf{z})$ (and to a lesser extend of $p_{\theta_z}(\mathbf{z})$). The observed variable $\mathbf{x}$ can be a continuous or discrete random variable with any arbitrary conditional distribution. The Gaussian case was then presented by Kingma and Welling (2014) as a major example. Of course, other pdfs (rather
than Gaussian) can be used depending on the nature of the data vector **x**. For example, Gamma distributions better fit the natural statistics of speech/audio power spectra (Girin et al., 2019). For simplicity of presentation and consistency across models, in the present review, $p_{\theta_x}(\mathbf{x}|\mathbf{z})$ is assumed to be a Gaussian distribution with diagonal covariance matrix for all models; that is,

$$p_{\theta_x}(\mathbf{x}|\mathbf{z}) = \mathcal{N}(\mathbf{x}; \boldsymbol{\mu}_{\theta_x}(\mathbf{z}), \mathrm{diag}\{\boldsymbol{\sigma}_{\theta_x}^2(\mathbf{z})\}) \quad (2.3)$$

$$= \prod_{f=1}^{F} p_{\theta_x}(x_f|\mathbf{z}) = \prod_{f=1}^{F} \mathcal{N}(x_f; \boldsymbol{\mu}_{\theta_x}, f(\mathbf{z}), \boldsymbol{\sigma}_{\theta_x}^2, f(\mathbf{z})), \quad (2.4)$$

where the subscript *f* denotes the *f*-th entry of a vector, and $\boldsymbol{\mu}_{\theta_x}$ : $\mathbb{R}^L \mapsto \mathbb{R}^F$ and $\boldsymbol{\sigma}_{\theta_x}$ : $\mathbb{R}^L \mapsto \mathbb{R}_+^F$ are nonlinear functions of **z** modeled by the decoder DNN. Although from a mathematical perspective, we could choose to work with full covariance matrices, assuming diagonal covariance matrices is preferable for computational reasons, since the number of free parameters of a covariance matrix grows quadratically with the variable dimension. This is problematic not only because we need to learn the neural network that computes all these parameters but also because covariance matrices often need to be inverted. The use of full covariance matrices also requires choosing an appropriate representation, for instance, based on the Cholesky decomposition. Please refer to Section 2.5.1 of Kingma and Welling (2019) for an extended discussion on this topic.

For maintaining consistency with the presentation of the other models in the next sections, we gather into $d_x$ the functions implemented by the decoder DNN; that is,

$$[\boldsymbol{\mu}_{\theta_x}(\mathbf{z}), \boldsymbol{\sigma}_{\theta_x}(\mathbf{z})] = d_x(\mathbf{z}). \quad (2.5)$$

A VAE decoder can be considered a generalization of the probabilistic principal component analysis (PPCA) (Tipping and Bishop, 1999) with a nonlinear (instead of linear) relationship between **z** and the parameters $\theta_x$. It can also be considered the generalization of a generative mixture models, with a continuous conditional latent variable instead of a discrete one (Kingma and Welling, 2019). Indeed, the marginal distribution of
$\mathbf{x}$, $p_{\theta}(\mathbf{x})$, is given by

$$p_{\theta}(\mathbf{x}) = \int p_{\theta_x}(\mathbf{x} | \mathbf{z}) p_{\theta_z}(\mathbf{z}) d\mathbf{z}. \qquad (2.6)$$

As any conditional distribution $p_{\theta_x}(\mathbf{x} | \mathbf{z})$ can provide a mode, $p_{\theta}(\mathbf{x})$ can be highly multimodal (in addition to being potentially high-dimensional). Unlike PPCA, in VAEs, the posterior distribution cannot be written analytically and has to be approximated, as discussed in the next section.

## 2.3 Learning with variational inference

Training the generative model defined in (2.1)–(2.3) amounts to estimating the parameters $\theta$ so as to minimize the Kullback-Leibler (KL) divergence between the true data distribution $p^*(\mathbf{x})$ and the model distribution (i.e., the marginal likelihood) $p_{\theta}(\mathbf{x})$:

$$\begin{align}
\min_{\theta} & \left\{ D_{\mathrm{KL}}(p^*(\mathbf{x}) \| p_{\theta}(\mathbf{x})) = \mathbb{E}_{p^*(\mathbf{x})} [\log p^*(\mathbf{x}) - \log p_{\theta}(\mathbf{x})] \right\} \nonumber \\
\Leftrightarrow & \max_{\theta} \mathbb{E}_{p^*(\mathbf{x})} [\log p_{\theta}(\mathbf{x})]. \tag{2.7}
\end{align}$$

This equivalence relation shows that this definition of model training actually corresponds to the maximum (marginal) likelihood parameter estimation. In practice, the true data distribution $p^*(\mathbf{x})$ is unknown, but we assume the availability of a training dataset $\mathbf{X} = \{\mathbf{x}_n \in \mathbb{R}^F\}_{n=1}^N$, where the training examples $\mathbf{x}_n$ are independent and identically distributed (i.i.d.) according to $p^*(\mathbf{x})$. Following the principle of empirical risk minimization, where the risk is here defined as the negative log-marginal likelihood, the intractable expectation in (2.7) is replaced by a Monte Carlo estimate:

$$\max_{\theta} \frac{1}{N} \sum_{n=1}^{N} \log p_{\theta}(\mathbf{x}_n), \quad \mathbf{x}_n \stackrel{\text{i.i.d.}}{\sim} p^*(\mathbf{x}). \qquad (2.8)$$

The estimated model parameters can then be used, for example, to generate new data from (2.1).

For many generative models with latent variables, directly solving this optimization problem is difficult, if not impossible when the marginal likelihood is analytically intractable, because of the integral
in (2.6), which cannot be computed in closed form. For the VAE, this intractability arises from the nonlinear relationship between the latent and observed variables, the latter being generated from the former through a DNN, which makes $p_{\theta_x}(x|z)$ in (2.6) a nonlinear function of $z$. One standard approach then involves leveraging the latent variable nature of the model in to maximize a lower bound of the intractable log-marginal likelihood (Neal and Hinton, 1998), which precisely depends on the *posterior distribution* of the latent variables or its approximation. This strategy leads to the expectation-maximization (EM) algorithm (Dempster et al., 1977) and its variants when the posterior distribution is intractable, such as Monte Carlo EM (Wei and Tanner, 1990) and variational EM (Jordan et al., 1999) algorithms.

As the name suggests, a VAE builds upon variational inference techniques, the general principles of which will now be briefly reviewed. Let $\mathcal{F}$ denote a *variational family* defined as a set of pdfs over the latent variables $z$. For any variational distribution of pdf $q(z) \in \mathcal{F}$, the following decomposition of the log-marginal likelihood holds (Neal and Hinton, 1998):

$$ \log p_{\theta}(\mathbf{x}) = \mathcal{L}(\theta, q(\mathbf{z}); \mathbf{x}) + D_{\text{KL}}(q(\mathbf{z}) \parallel p_{\theta}(\mathbf{z}|\mathbf{x})), \quad (2.9) $$

where $\mathcal{L}(\theta, q(\mathbf{z}); \mathbf{x})$ is referred to in the literature as the evidence lower bound (ELBO), the negative variational free energy, or the VLB, and is defined as

$$ \mathcal{L}(\theta, q(\mathbf{z}); \mathbf{x}) = \mathbb{E}_{q(\mathbf{z})} \left[ \log p_{\theta}(\mathbf{x}, \mathbf{z}) - \log q(\mathbf{z}) \right] \leq \log p_{\theta}(\mathbf{x}). \quad (2.10) $$

The inequality in (2.10) is obtained from (2.9) based on the fact that $D_{\text{KL}}(\cdot \parallel \cdot) \geq 0$. Equality holds (i.e., the VLB is tight to the log-marginal likelihood) if and only if the variational distribution $q(\mathbf{z})$ is equal to the exact posterior distribution $p_{\theta}(\mathbf{z}|\mathbf{x})$.

The EM algorithm (Dempster et al., 1977) is an iterative algorithm that consists in alternatively maximizing the VLB with respect to $q(\mathbf{z}) \in \mathcal{F}$ in the E-step and with respect to $\theta$ in the M-step (Neal and Hinton, 1998). From (2.9), we see that the E-step involves finding the variational distribution $q(\mathbf{z})$ in the variational family $\mathcal{F}$ that best approximates the true posterior $p_{\theta}(\mathbf{z}|\mathbf{x})$ according to the KL divergence
measure of fit:

$$q^*(\mathbf{z}) = \arg\max_{q \in \mathcal{F}} \mathcal{L}(\theta, q(\mathbf{z}); \mathbf{x}) = \arg\min_{q \in \mathcal{F}} D_{\text{KL}}(q(\mathbf{z}) \parallel p_{\theta}(\mathbf{z}|\mathbf{x})). \quad (2.11)$$

In the exact EM algorithm, the variational family $\mathcal{F}$ is unconstrained, so the solution to the E-step is given by the exact posterior distribution: $q^*(\mathbf{z}) = p_{\theta}(\mathbf{z}|\mathbf{x})$. As this optimal variational distribution over $\mathbf{z}$ is actually conditioned on $\mathbf{x}$, we will now use the notation $q(\mathbf{z}|\mathbf{x})$ instead of $q(\mathbf{z})$. The difficulty arises when the posterior distribution $p_{\theta}(\mathbf{z}|\mathbf{x})$ is intractable, which prevents us from solving the E-step analytically. Variational inference then consists in constraining the variational family $\mathcal{F}$ and resorting to optimization methods for solving the E-step (Jordan et al., 1999).

Seminal works on variational inference relied on the so-called mean-field approximation, which constrains the variational family $\mathcal{F}$ to be a set of completely factorized distributions (i.e., multivariate distributions over $\mathbf{z}$ that are written as a product of univariate marginal distributions over the entries of $\mathbf{z}$). All marginal posterior dependencies between different entries of $\mathbf{z}$ are ignored here, while the “structured” mean-field approximation (Saul and Jordan, 1996) partially restores some of them. Solving the E-step under the mean-field approximation leads to a set of closed-form coupled solutions for each univariate distribution in the factorization. This approach is also referred to as coordinate-ascent variational inference in the literature (Bishop, 2006; Blei et al., 2017). However, closed-form updates are usually only available for conjugate-exponential models (Winn and Bishop, 2005) when the distribution of each scalar latent variable, conditionally on its parents, belongs to the exponential family and is conjugate with respect to the distribution of these parent variables. Moreover, this coordinate-ascent approach does not scale well for high-dimensional and large-scale inference problems (Hoffman et al., 2013).

An alternative to the mean-field approximation is then to define the variational family as a set of distributions with a certain parametric form $q_{\lambda}(\mathbf{z}|\mathbf{x})$, where the parameters $\lambda$ govern the shape of the distribution. For example, we can define the Gaussian variational family where the
parameters $\lambda$ correspond to the mean vector and covariance matrix:

$$ \mathcal{F} = \{ q_{\lambda}(\mathbf{z}|\mathbf{x}) = \mathcal{N}(\mathbf{z}; \boldsymbol{\mu}, \boldsymbol{\Sigma}), \lambda = \{\boldsymbol{\mu}, \boldsymbol{\Sigma}\} \}. \quad (2.12) $$

As shown below, the optimal parameters $\lambda$ that maximize the VLB depend on $\mathbf{x}$ and $\theta$. This approach is called fixed-form or structured variational inference (Honkela et al., 2010; Salimans and Knowles, 2013). The VLB in (2.10) then becomes a function of both the generative model parameters $\theta$ and variational parameters $\lambda$:

$$ \mathcal{L}(\theta, \lambda; \mathbf{x}) = \mathbb{E}_{q_{\lambda}(\mathbf{z}|\mathbf{x})} \left[ \log p_{\theta}(\mathbf{x}, \mathbf{z}) - \log q_{\lambda}(\mathbf{z}|\mathbf{x}) \right]. \quad (2.13) $$

The E-step in (2.11) consequently reduces to a parametric optimization problem:

$$ \lambda^* = \arg\max_{\lambda} \mathcal{L}(\theta, \lambda; \mathbf{x}) = \arg\min_{\lambda} D_{\text{KL}}(q_{\lambda}(\mathbf{z}|\mathbf{x}) \parallel p_{\theta}(\mathbf{z}|\mathbf{x})). \quad (2.14) $$

As the objective function depends on the observed data vector $\mathbf{x}$ and the generative model parameters $\theta$, so does the solution $\lambda^*$. In fact, the optimal variational distribution depends on $\mathbf{x}$ through the parameters $\lambda^*$. The M-step remains unchanged in fixed-form variational inference; that is, it consists in updating the generative model parameters by maximizing $\mathcal{L}(\theta, \lambda; \mathbf{x})$ w.r.t. $\theta$, using the current estimate of the variational parameters. If the expectation in (2.13) and its gradient w.r.t. $\lambda$ can be computed analytically, the optimization problem of the E-step can be solved using gradient-based optimization methods.

In general, given a dataset of i.i.d. data vectors $\mathbf{X} = \{\mathbf{x}_1, ..., \mathbf{x}_N\}$, one needs to find the parameters $\Lambda = \{\lambda_1, ..., \lambda_N\}$ of the variational distributions $q_{\lambda_n}(\mathbf{z}_n|\mathbf{x}_n)$, $n = 1, ..., N$. Taking the same example as before, with $q_{\lambda_n}(\mathbf{z}_n|\mathbf{x}_n) = \mathcal{N}(\mathbf{z}_n; \boldsymbol{\mu}_n, \boldsymbol{\Sigma}_n)$, we have here $\lambda_n = \{\boldsymbol{\mu}_n, \boldsymbol{\Sigma}_n\}$. This problem is solved by maximizing the following total VLB, which is the sum (or equivalently, the mean) of the local VLB defined in (2.13) over each vector in the training dataset:

$$ \mathcal{L}(\theta, \Lambda; \mathbf{X}) = \sum_{n=1}^{N} \mathcal{L}(\theta, \lambda_n; \mathbf{x}_n). \quad (2.15) $$

To scale to large amounts of data, stochastic variational inference (Hoffman et al., 2013) relies on gradient-based stochastic optimization
(Robbins and Monro, 1951; Bottou, 2004) for maximizing the total VLB in (2.15) w.r.t. the generative model parameters $\theta$. The gradient of the total VLB, $\mathcal{L}(\theta, \Lambda; \mathbf{X})$, is the sum of the gradients of the local VLBs, $\mathcal{L}(\theta, \lambda_n; \mathbf{x}_n)$, defined for each sample $\mathbf{x}_n$ in the dataset. For large datasets, computing this sum to perform a single update of $\theta$ with a step of gradient ascent can be inefficient. Therefore, stochastic variational inference exploits a noisy stochastic estimate of the gradient, computed from a single example $\mathbf{x}_n$ or from a mini-batch of examples in the dataset. This is the same principle as that used in stochastic and mini-batch gradient descent (Bottou, 2004), such that stochastic variational inference inherits from the same convergence properties (Robbins and Monro, 1951).

However, the estimation of the complete set of variational parameters can remain expensive for large datasets. Thus, amortized variational inference makes a stronger assumption for defining the variational family, by introducing an inference model $f_{\phi}$ such that

$$ \lambda_n = f_{\phi}(\mathbf{x}_n), \qquad (2.16) $$

where $\phi$ is a set of parameters that is shared among all variational distributions $q_{\lambda_n}(\mathbf{z}_n|\mathbf{x}_n)$. This inference model is used to map the observation $\mathbf{x}_n$ to the local variational parameter $\lambda_n$. The variational family $\mathcal{F}$ then corresponds to the set of variational distributions parametrized by $\phi$, which are denoted by $q_{\phi}(\mathbf{z}_n|\mathbf{x}_n)$. For instance, for $q_{\phi}(\mathbf{z}_n|\mathbf{x}_n) = \mathcal{N}(\mathbf{z}_n; \boldsymbol{\mu}_n, \boldsymbol{\Sigma}_n)$, we have $\lambda_n = [\boldsymbol{\mu}_n, \boldsymbol{\Sigma}_n] = f_{\phi}(\mathbf{x}_n)$. This amortization principle corresponds to a stronger assumption for the variational family compared to nonamortized fixed-form and mean-field approximations. Therefore, the KL divergence between the exact posterior and its approximation is likely to be larger in the amortized case than in the previous cases. The total VLB for the complete training dataset then becomes a function of $\phi$:

$$ \mathcal{L}(\theta, \phi; \mathbf{X}) = \sum_{n=1}^{N} \mathbb{E}_{q_{\phi}(\mathbf{z}_n|\mathbf{x}_n)}[\log p_{\theta}(\mathbf{x}_n, \mathbf{z}_n) - \log q_{\phi}(\mathbf{z}_n|\mathbf{x}_n)]. \quad (2.17) $$

This means that the optimization of the set of local variational parameters $\Lambda = \{\lambda_1, ..., \lambda_N\}$ is replaced by the optimization of the shared
set of inference model parameters $\phi$. Hereinafter, we will use the term inference model to directly denote the variational distribution $q_{\phi}(\mathbf{z}_n|\mathbf{x}_n)$.

## 2.4 VAE inference model

VAEs belong to the family of amortized variational inference techniques, where the VLB in (2.17) is optimized using stochastic gradient-based optimization techniques. The VAE generative model $p_{\theta}(\mathbf{x}, \mathbf{z})$ has already been defined in (2.1)–(2.3). It involves a decoder neural network through $p_{\theta}(\mathbf{x}|\mathbf{z})$. To fully specify the VLB, which is required to learn the generative model parameters $\theta$, it is also necessary to define the inference model $q_{\phi}(\mathbf{z}|\mathbf{x})$, which approximates the intractable exact posterior $p_{\theta}(\mathbf{z}|\mathbf{x})$.

Similar to the generative model, the inference model for $q_{\phi}(\mathbf{z}|\mathbf{x})$ is defined by an encoder neural network. A common choice for the approximate posterior distribution $q_{\phi}(\mathbf{z}|\mathbf{x})$ is to use a Gaussian distribution:

$$q_{\phi}(\mathbf{z}|\mathbf{x}) = \mathcal{N}(\mathbf{z}; \boldsymbol{\mu}_{\phi}(\mathbf{x}), \mathrm{diag}\{\boldsymbol{\sigma}_{\phi}^2(\mathbf{x})\}) \quad (2.18)$$

$$= \prod_{l=1}^{L} q_{\phi}(z_l|\mathbf{x}) = \prod_{l=1}^{L} \mathcal{N}(z_l; \mu_{\phi,l}(\mathbf{x}), \sigma_{\phi,l}^2(\mathbf{x})), \quad (2.19)$$

where index $l \in \{1, ..., L\}$ is used to denote the $l$-th entry of the corresponding vectors, and $\boldsymbol{\mu}_{\phi}: \mathbb{R}^F \mapsto \mathbb{R}^L$ and $\boldsymbol{\sigma}_{\phi}: \mathbb{R}^F \mapsto \mathbb{R}_+^L$ are nonlinear functions of $\mathbf{x}$, modeled by a DNN called the encoder or recognition network, which is parametrized by a set of weights and biases denoted by $\phi$. The encoder network is illustrated in Figure 2.2 (left). As for the VAE generative model, for the sake of consistency with the presentation of the other models, we denote

$$[\boldsymbol{\mu}_{\phi}(\mathbf{x}), \boldsymbol{\sigma}_{\phi}(\mathbf{x})] = e_{\mathbf{z}}(\mathbf{x}), \quad (2.20)$$

where $e_{\mathbf{z}}$ is the nonlinear function implemented by the encoder DNN.

## 2.5 VAE training

In the VAE methodology (Kingma and Welling, 2014; Rezende et al., 2014), the VLB in (2.17) is optimized using stochastic gradient-based
optimization techniques to learn the generative and inference model parameters. For training the VAE, the encoder and decoder networks are cascaded, as illustrated in Figure 2.2, and the sets of parameter $\theta$ and $\phi$ are jointly estimated from the training data $\mathbf{X}$. This is different from an EM algorithm strategy, which would alternatively optimize the VLB w.r.t. $\phi$ and $\theta$ in the E- and M-steps, respectively. He et al. (2018) showed that this joint encoder-decoder training of the VAE can, however, be suboptimal.

The VLB in (2.17) can be reshaped as (Kingma and Welling, 2014)

$$ \mathcal{L}(\theta, \phi; \mathbf{X}) = \underbrace{\sum_{n=1}^{N} \mathbb{E}_{q_{\phi}(\mathbf{z}_n|\mathbf{x}_n)} [\log p_{\theta_x}(\mathbf{x}_n|\mathbf{z}_n)]}_{\text{Reconstruction accuracy}} - \underbrace{\sum_{n=1}^{N} D_{KL}(q_{\phi}(\mathbf{z}_n|\mathbf{x}_n) \| p_{\theta_z}(\mathbf{z}_n))}_{\text{Regularization}}. \quad (2.21) $$

The first term on the right-hand side of (2.21) is a reconstruction term that represents the average accuracy of the chained encoding-decoding process. For instance, if the generative model $p_{\theta_x}(\mathbf{x}|\mathbf{z})$ is chosen to be Gaussian with an identity covariance matrix, the reconstruction term is equal to the opposite of the mean-squared error (MSE) between the original data and decoder output, up to additive constants. The second term is a regularization one, which enforces the approximate posterior distribution $q_{\phi}(\mathbf{z}|\mathbf{x})$ to be close to the prior distribution $p_{\theta_z}(\mathbf{z})$. Provided that an independent Gaussian prior is used, this term forces $\mathbf{z}$ to be a disentangled data representation; that is, the $\mathbf{z}$ entries tend to be independent and encode a different characteristic (or factor of variation) of the data.

For usual distributions, the regularization term has an analytical expression as a function of $\theta$ and $\phi$. However, the expectation taken with respect to $q_{\phi}(\mathbf{z}_n|\mathbf{x}_n)$ in the reconstruction accuracy term is analytically intractable. Therefore, in practice, it is approximated using a Monte Carlo estimate with $R$ samples $\mathbf{z}_n^{(r)}$ independently and identically drawn
from $q_{\phi}(\mathbf{z}_n|\mathbf{x}_n)$ (for each index $n$):

$$ \mathbb{E}_{q_{\phi}(\mathbf{z}_n|\mathbf{x}_n)}[\log p_{\theta}(\mathbf{x}_n|\mathbf{z}_n)] \approx \frac{1}{R} \sum_{r=1}^{R} \log p_{\theta}(\mathbf{x}_n|\mathbf{z}_n^{(r)}). \quad (2.22) $$

The resulting Monte Carlo estimate of the VLB is given by

$$ \hat{\mathcal{L}}(\theta, \phi; \mathbf{X}) = \sum_{n=1}^{N} \frac{1}{R} \sum_{r=1}^{R} \log p_{\theta}(\mathbf{x}_n|\mathbf{z}_n^{(r)}) - \sum_{n=1}^{N} D_{KL}(q_{\phi}(\mathbf{z}_n|\mathbf{x}_n) \| p(\mathbf{z}_n)). \quad (2.23) $$

To optimize this objective function, we can typically resort to the (variants of) stochastic or mini-batch gradient descent (on the negative VLB) (Bottou, 2004). While the gradient of $\hat{\mathcal{L}}(\theta, \phi; \mathbf{X})$ w.r.t. $\theta$ can be easily computed using the standard backpropagation algorithm, that w.r.t. $\phi$ is problematic because the sampling operation from $q_{\phi}(\mathbf{z}_n|\mathbf{x}_n)$ is not differentiable w.r.t. $\phi$. The solution to this problem, proposed by Kingma and Welling (2014) and referred to as the *reparameterization trick*, consists in reparametrizing the sample $\mathbf{z}_n^{(r)}$ using a differentiable transformation of a sample $\boldsymbol{\epsilon}^{(r)}$ drawn from a standard Gaussian distribution, which does not depend on $\phi$:

$$ \mathbf{z}_n^{(r)} = \boldsymbol{\mu}_{\phi}(\mathbf{x}_n) + \operatorname{diag}\{\boldsymbol{\sigma}_{\phi}^2(\mathbf{x}_n)\}^{\frac{1}{2}} \boldsymbol{\epsilon}^{(r)}, \quad \boldsymbol{\epsilon}^{(r)} \sim \mathcal{N}(\mathbf{0}_L, \mathbf{I}_L). \quad (2.24) $$

Using this reparameterization trick, $\hat{\mathcal{L}}(\theta, \phi; \mathbf{X})$ is now differentiable w.r.t. $\phi$. This differentiable Monte Carlo approximation of the VLB is referred to as the stochastic gradient variational Bayes (SGVB) estimator (Kingma and Welling, 2014). The gradient of $\hat{\mathcal{L}}(\theta, \phi; \mathbf{X})$ w.r.t. $\phi$ is an unbiased estimate of the gradient of the exact VLB $\mathcal{L}(\theta, \phi; \mathbf{X})$ (Kingma and Welling, 2019). This property allows using very few samples to compute the SGVB estimator, which however impacts the variance of the estimator. Kingma and Welling (2014) suggested setting $R = 1$ provided that sufficiently large mini-batches are used for the gradient descent. This training procedure of a VAE model is now considered routine within deep learning toolkits, such as TensorFlow (Abadi et al., 2016) and PyTorch (Paszke et al., 2019).
# Recurrent Neural Networks and State Space Models

As mentioned earlier, DVAEs are formed of combinations of a VAE and temporal models. Most of these temporal models rely on RNNs and/or SSMs. We thus briefly present the basics of RNNs and SSMs in this chapter before moving on to DVAEs in the next chapters. An extended technical overview of RNNs and SSMs, as well as their applications, is beyond the scope of the present paper.

## 3.1 Recurrent Neural Networks

### 3.1.1 Principle and definition

RNNs have been and are still widely used for data sequence modeling and generation and sequence-to-sequence mapping. An RNN is a neural network that processes ordered vector sequences and uses a memory of past input/output data to condition the current output (Sutskever, 2013; Graves et al., 2013). This is achieved using an additional vector that recursively encodes the internal state of the network.

We denote by $x_{t1:t2} = \{x_t\}_{t=t_1}^{t=t_2}$ a sequence of vectors $x_t$ indexed from $t_1$ to $t_2$, where $t_1 \le t_2$. When $t_1 > t_2$, we assume $x_{t1:t2} = \emptyset$. We present RNNs in the general framework of nonlinear systems, which
transform an input vector sequence $\mathbf{u}_{1:T}$ into an output vector sequence $\mathbf{x}_{1:T}$, possibly through an internal state vector sequence $\mathbf{h}_{1:T}$. The input, output, and internal state vectors can have arbitrary (different) dimensions. If $\mathbf{u}_{1:T}$ is an “external” input sequence, the network can be considered as a “system,” as is usual in control theory ($\mathbf{u}_{1:T}$ being considered as a command to the system). If $\mathbf{u}_t = \emptyset$, the RNN is in the *undriven mode*. In contrast, if $\mathbf{u}_t = \mathbf{x}_{t-1}$, the RNN is in the *predictive mode*, or *sequence generation mode*, which is a usual mode when we are interested in modeling the evolution of a data sequence $\mathbf{x}_{1:T}$ “alone” (i.e., independently of any external input; in this case, $\mathbf{x}_{1:T}$ can be seen both as an input and an output sequence).

A basic single-layer RNN model is defined by

$$ \mathbf{h}_t = d_{\text{hid}}(\mathbf{W}_{\text{in}}\mathbf{u}_t + \mathbf{W}_{\text{rec}}\mathbf{h}_{t-1} + \mathbf{b}_{\text{hid}}), \quad (3.1) $$

$$ \mathbf{x}_t = d_{\text{out}}(\mathbf{W}_{\text{out}}\mathbf{h}_t + \mathbf{b}_{\text{out}}), \quad (3.2) $$

where $\mathbf{W}_{\text{in}}$, $\mathbf{W}_{\text{rec}}$ and $\mathbf{W}_{\text{out}}$ are weight matrices of appropriate dimensions; $\mathbf{b}_{\text{hid}}$ and $\mathbf{b}_{\text{out}}$ are bias vectors; and $d_{\text{hid}}$ and $d_{\text{out}}$ are nonlinear activation functions. We also define the initial internal state vector $\mathbf{h}_0$. This model is extendable to more complex recurrent architectures:

$$ \mathbf{h}_t = d_{\mathbf{h}}(\mathbf{u}_t, \mathbf{h}_{t-1}), \quad (3.3) $$

$$ \mathbf{x}_t = d_{\mathbf{x}}(\mathbf{h}_t), \quad (3.4) $$

where $d_{\mathbf{h}}$ and $d_{\mathbf{x}}$ denote any arbitrary complex nonlinear functions implemented with a DNN. We assume that this representation includes long short-term memory (LSTM) networks (Hochreiter and Schmidhuber, 1997) and gated recurrent unit (GRU) networks (Cho et al., 2014), which comprise additional internal variables called gates. For simplicity of presentation, these additional internal gates are not formalized in (3.3) and (3.4). The same is true for multi-layer RNNs, where several recursive layers are stacked on top of each other (Graves et al., 2013) (in this case, for the same reason, we do not report layer indexes in (3.3) and (3.4)). This is also true for combinations of multi-layer RNNs and LSTMs (i.e., multi-layer LSTM networks). In summary, we assume that (3.3) and (3.4) are a “generic” or “high-level” representation of an RNN of arbitrary complexity.
**Notation remark:** To clarify the presentation and links between the different models, we use the same generic notation $d_x$ for the generating function in (2.5) and (3.4), and will do that throughout the paper (and the same for $d_h$ and for $d_z$ later in the paper).

So far, the above RNNs are deterministic: given $u_{1:T}$ and $h_0$, $x_{1:T}$ is completely determined. Such networks are trained by optimizing a deterministic criterion, e.g. the MSE between the target output sequences from a training dataset and the corresponding actual output sequences obtained by the network. The training set of i.i.d. vectors used for VAE training is replaced with that of vector sequences, and consecutive vectors within a training sequence are generally correlated, which is the point of using a dynamical model.

### 3.1.2 Generative recurrent neural networks

Deterministic RNNs can easily be transformed into *generative* RNNs (GRNNs) by adding stochasticity at the output level. We just have to define a probabilistic observation model and replace the output data sequence with an output sequence of distribution parameters, similar to the VAE decoder:

$$h_t = d_h(u_t, h_{t-1}), \qquad (3.5)$$

$$[\mu_{\theta_x}(h_t), \sigma_{\theta_x}(h_t)] = d_x(h_t), \qquad (3.6)$$

$$p_{\theta_x}(x_t|h_t) = \mathcal{N}(x_t; \mu_{\theta_x}(h_t), \text{diag}\{\sigma_{\theta_x}^2(h_t)\}). \qquad (3.7)$$

Eq. (3.5) is the same recursive internal state model as (3.3). Eqs. (3.6) and (3.7) constitute the observation model. In (3.7), we use the Gaussian distribution for its generality and for the convenience of illustration, although any distribution can be used, just as for the VAE decoder. Again, one may choose a distribution that is more appropriate for the nature of the data. For example, Graves (2013) proposed using mixture distributions. The complete set of model parameters $\theta$ here includes $\theta_h$ and $\theta_x$, the parameters of the networks implementing $d_h$ and $d_x$, respectively. Because the output of $d_x$ in (3.6) is now two vectors of pdf parameters instead of a data vector in (3.4), its size is twice that of the deterministic RNN. When the internal state vector $h_t$ is of (much)
lower dimension than the output vector $\mathbf{x}_t$, the GRNN observation model becomes similar to the VAE decoder, except that, again, $\mathbf{h}_t$ has a deterministic evolution through time, whereas the latent state $\mathbf{z}$ of the VAE is stochastic and i.i.d., which is a fundamental difference.

Even if the generation of $\mathbf{x}_t$ is now stochastic, the evolution of the internal state is still deterministic. Let us denote $\mathbf{h}_t$ as a function $\mathbf{h}_t = \mathbf{h}_t(\mathbf{u}_{1:t})$ to make the deterministic relation between $\mathbf{u}_{1:t}$ and $\mathbf{h}_t$ explicit (for each time index $t$).<fnref n="1" /> We thus have $p_{\theta_x}(\mathbf{x}_t|\mathbf{h}_t) = p_{\theta_x}(\mathbf{x}_t|\mathbf{h}_t(\mathbf{u}_{1:t}))$. In the predictive mode, we have $p_{\theta_x}(\mathbf{x}_t|\mathbf{h}_t) = p_{\theta_x}(\mathbf{x}_t|\mathbf{h}_t(\mathbf{x}_{0:t-1}))$.<fnref n="2" /> Such stochastic version of the RNN can be trained with a statistical criterion (e.g., maximum likelihood). As for the VAE training, we search for the maximization of the observed data log-likelihood w.r.t. $\theta$ over a set of training sequences. For one sequence, with the conditional independence of successive data vectors, the data log-likelihood is given by

$$ \log p_{\theta_x}(\mathbf{x}_{1:T}|\mathbf{u}_{1:T}) = \sum_{t=1}^{T} \log p_{\theta_x}(\mathbf{x}_t|\mathbf{h}_t(\mathbf{u}_{1:t})). \quad (3.8) $$

## 3.2 State Space Models

### 3.2.1 Principle and definition

SSMs are a rich family of models that are widely used to model dynamical systems (e.g., in statistical signal processing, time-series analysis, and control theory) (Durbin and Koopman, 2012). Here, we focus on discrete-time, continuous-valued SSMs of the form

$$ [\mu_{\theta_z}(\mathbf{z}_{t-1}, \mathbf{u}_t), \sigma_{\theta_z}(\mathbf{z}_{t-1}, \mathbf{u}_t)] = d_z(\mathbf{z}_{t-1}, \mathbf{u}_t), \quad (3.9) $$

$$ p_{\theta_z}(\mathbf{z}_t|\mathbf{z}_{t-1}, \mathbf{u}_t) = \mathcal{N}(\mathbf{z}_t; \mu_{\theta_z}(\mathbf{z}_{t-1}, \mathbf{u}_t), \text{diag}\{\sigma_{\theta_z}^2(\mathbf{z}_{t-1}, \mathbf{u}_t)\}), \quad (3.10) $$

$$ [\mu_{\theta_x}(\mathbf{z}_t), \sigma_{\theta_x}(\mathbf{z}_t)] = d_x(\mathbf{z}_t), \quad (3.11) $$

$$ p_{\theta_x}(\mathbf{x}_t|\mathbf{z}_t) = \mathcal{N}(\mathbf{x}_t; \mu_{\theta_x}(\mathbf{z}_t), \text{diag}\{\sigma_{\theta_x}^2(\mathbf{z}_t)\}), \quad (3.12) $$

<footnote n="1">¹$\mathbf{h}_t(\mathbf{u}_{1:t})$ also depends on the initial internal state vector $\mathbf{h}_0$, but we omit this term as an argument of the function for conciseness.</footnote>

<footnote n="2">²Here, the first “input” $\mathbf{x}_0$ has to be set arbitrarily, just like $\mathbf{h}_0$. Alternately, one can directly start the generation process from an arbitrary internal state vector $\mathbf{h}_1$.</footnote>
where $d_z$ and $d_x$ are functions of arbitrary complexity, each being parameterized by a set of parameters denoted $\theta_z$ and $\theta_x$, respectively. As for the complete generative model, we have $\theta = \theta_x \cup \theta_z$, and we retain this notation hereinafter. At this point, $d_z$ and $d_x$ can be linear or nonlinear functions, and we will differentiate the two cases later. The observation model (3.11)–(3.12) is very similar to the GRNN observation model (3.6)–(3.7). However, $z_t$ is here a stochastic internal state vector in contrast to the deterministic internal state $h_t$ of the (G)RNN. The distribution of $z_t$, known as the state model or the dynamical model, is given by (3.9)–(3.10). It follows a first-order Markov model; that is, a temporal dependency is introduced where $z_t$ depends on the previous state $z_{t-1}$ and the corresponding input $u_t$ through the function $d_z$. In short, the above SSM can be considered a GRNN in which the deterministic internal state $h_t$ is replaced with a stochastic internal state $z_t$, as illustrated in Figure 3.1.

**Notation remark:** In the control theory literature, the input corresponding to the generation of $z_t$ is often denoted as $u_{t-1}$, or equivalently, $u_t$ is used to generate the next state $z_{t+1}$. This notation is arbitrary. In the present paper, we prefer to realign the temporal indices, so that the input $u_t$ is used to generate $z_t$, which in turn is used to generate $x_t$, to maintain better consistency through all presented models.

As for the complete sequence, given the dependencies represented in Figure 3.1, the joint distribution of all variables can be expressed as

$$p_{\theta}(\mathbf{x}_{1:T}, \mathbf{z}_{1:T}, \mathbf{u}_{1:T}) = \prod_{t=1}^{T} p_{\theta_x}(\mathbf{x}_t|\mathbf{z}_t) p_{\theta_z}(\mathbf{z}_t|\mathbf{z}_{t-1}, \mathbf{u}_t) p(\mathbf{u}_t), \quad (3.13)$$

from which we can deduce

$$p_{\theta_x}(\mathbf{x}_{1:T}|\mathbf{z}_{1:T}) = \prod_{t=1}^{T} p_{\theta_x}(\mathbf{x}_t|\mathbf{z}_t), \quad (3.14)$$

and

$$p_{\theta_z}(\mathbf{z}_{1:T}|\mathbf{u}_{1:T}) = \prod_{t=1}^{T} p_{\theta_z}(\mathbf{z}_t|\mathbf{z}_{t-1}, \mathbf{u}_t). \quad (3.15)$$
<figure><img src="image_6.png" /><figcaption>**Figure 3.1:** GRNN (left) against SSM (right): The two models have an identical structure, though the internal state of the GRNN is deterministic (represented with a diamond), whereas that of the SSM is stochastic (represented with a circle).</figcaption></figure>

Given the state sequence $z_{1:T}$, the observation vectors at different time frames are mutually independent. The prior distribution of $u_t$ also factorizes across time frames, but this is of limited interest here. To be complete, we should specify the model “initialization”: At $t = 1$, we need to define $z_0$, which can be set to an arbitrary deterministic value, or defined through a prior distribution $p_{\theta_z}(z_0)$ (which then must be added to the right-hand side of (3.13) and (3.15)), or we can set $z_0 = \emptyset$, in which case the first term of the state model in these equations is $p_{\theta_z}(z_1|u_1)$.

Solving the above SSM means that we run the inference process; that is, we estimate the state vector sequence $z_{1:T}$ from an observed data vector sequence $x_{1:T}$. The use of Gaussian distribution in (3.10) and (3.12) is a convenient choice that generally facilitates inference. More generally, these distributions are within the exponential family, so either exact or approximate inference algorithms can be applied, depending on the nature of $d_z$ and $d_x$. In the next subsection, we provide an example of a closed-form inference solution when $d_z$ and $d_x$ are linear functions.

### 3.2.2 Kalman filters

Some classical SSMs have been successfully used for decades for a wide set of applications. For example, when $d_x$ and $d_z$ are linear functions
of the form

$$ \mu_{\theta_z}(\mathbf{z}_{t-1}, \mathbf{u}_t) = \mathbf{A}_t \mathbf{z}_{t-1} + \mathbf{B}_t \mathbf{u}_t + \mathbf{m}_t, \quad \sigma_{\theta_z}^2(\mathbf{z}_{t-1}, \mathbf{u}_t) = \mathbf{\Lambda}_t, \quad (3.16) $$

$$ \mu_{\theta_x}(\mathbf{z}_t) = \mathbf{C}_t \mathbf{z}_t + \mathbf{n}_t, \quad \sigma_{\theta_x}^2(\mathbf{z}_t) = \mathbf{\Sigma}_t, \quad (3.17) $$

where $\mathbf{A}_t$, $\mathbf{B}_t$, $\mathbf{m}_t$, $\mathbf{\Lambda}_t$, $\mathbf{C}_t$, $\mathbf{n}_t$, and $\mathbf{\Sigma}_t$ are matrices and vectors of appropriate size, the SMM transforms into a linear-Gaussian linear dynamical system (LG-LDS). In this case, the inference has a very popular closed-form solution, known as a Kalman filter (Moreno and Pigazo, 2009). More precisely, a Kalman filter is the solution obtained when the past and present observations (outputs and inputs) are used at each time $t$ (i.e., causal inference). When a complete sequence of observations is used at each time $t$ (i.e., noncausal inference), the solution is referred to as a Kalman smoother, also obtainable in closed form. In practical problems, $\mathbf{x}_{1:T}$ is generally noisy, and the terms “filter” and “smoother” refer to the estimation of a “clean” state vector trajectory $\mathbf{z}_{1:T}$ from noisy observed data.

The Kalman filter is an iterative solution that alternates between a prediction step and an update step. The prediction step involves computing the predictive distribution, which is the posterior distribution of $\mathbf{z}_t$ given the observations up to time $t-1$. Starting from the joint distribution of all variables and exploiting the dependencies in the generative model, the predictive distribution can be expressed as (we omit the input $\mathbf{u}$ for simplicity of presentation)

$$ p(\mathbf{z}_t|\mathbf{x}_{1:t-1}) = \int p(\mathbf{z}_t|\mathbf{z}_{t-1})p(\mathbf{z}_{t-1}|\mathbf{x}_{1:t-1})d\mathbf{z}_{t-1}. \quad (3.18) $$

The update step involves integrating the new (current) observation $\mathbf{x}_t$ using Bayes' rule to obtain the so-called filtering distribution (up to some normalizing factor that does not depend on $\mathbf{z}_t$):

$$ p(\mathbf{z}_t|\mathbf{x}_{1:t}) \propto p(\mathbf{x}_t|\mathbf{z}_t) \int p(\mathbf{z}_t|\mathbf{z}_{t-1})p(\mathbf{z}_{t-1}|\mathbf{x}_{1:t-1})d\mathbf{z}_{t-1}. \quad (3.19) $$

The filtering distribution at time $t$ can be computed recursively from the filtering distribution at time $t-1$ (inside the integral). In the case of linear-Gaussian generative distributions, the filtering distribution is Gaussian, with parameters that can be computed recursively from the
parameters at time $t-1$ and the generative model parameters with basic matrix/vector operations. In practice, these parameters are computed in the following two steps: prediction step and update step. Finally, the mean vector of the filtering distribution, which is often used as the state estimate, is a linear form of the observation vector.

In the noncausal case, a similar two-step predictive/update recursive process can be computed, except that the recursion is processed in both forward (causal) and backward (anticausal) directions, leading to the *smoothing distribution*. A more detailed presentation of the Kalman filter and Kalman smoother is beyond the scope of the present paper.

### 3.2.3 Nonlinear Kalman filters

Nonlinear dynamical systems (NDS), sometimes abusively referred to as nonlinear Kalman filters, have also been extensively studied, well before the deep learning era. Principled extensions to the Kalman Filter have been proposed to deal with the nonlinearities (e.g., the extended Kalman filter and the unscented Kalman filter) (Wan and Van Der Merwe, 2000; Daum, 2005). The review of nonlinear Kalman filters is beyond the scope of the present paper, to retain the focus on DVAEs.
# Definition of Dynamical VAEs

In this section, we describe a general methodology for defining and training dynamical VAEs. Our goal is to encompass different models proposed in the literature, which we will describe in detail later. These models can be considered particular instances of this general definition, given simplifying assumptions. This section will prepare the readers to understand well the commonalities and differences among all models that we will review and may motivate future developments. We first define a DVAE in terms of a generative model and then present the general lines of inference and training in the DVAE framework.

## 4.1 Generative model

As already mentioned, DVAEs consider a sequence of observed random vectors $x_{1:T} = \{x_t \in \mathbb{R}^F\}_{t=1}^T$ and that of latent random vectors $z_{1:T} = \{z_t \in \mathbb{R}^L\}_{t=1}^T$. As opposed to the a “static” VAE and similarly to SSMs, these two data sequences are assumed to be temporally correlated and can have somewhat complex (cross-)dependencies across time. Defining a DVAE generative model involves specifying the joint distribution of the observed and latent sequential data, $p_\theta(x_{1:T}, z_{1:T})$, the parameters of which are provided by DNNs, which themselves depend on a set of
parameters $\theta$.

When the model works in the so-called *driven mode*, one additionally considers an input sequence of observed random vectors $\mathbf{u}_{1:T} = \{\mathbf{u}_t \in \mathbb{R}^U\}_{t=1}^T$, and in that case, $\mathbf{x}_{1:T}$ is considered the output sequence. In this case, to define the full generative model, we need to specify the joint distribution $p_{\theta}(\mathbf{x}_{1:T}, \mathbf{z}_{1:T}, \mathbf{u}_{1:T})$. However, in practice, we are usually only interested in modeling the generative process of $\mathbf{x}_{1:T}$ and $\mathbf{z}_{1:T}$ given the input sequence $\mathbf{u}_{1:T}$. Loosely speaking, the input sequence is assumed deterministic, while $\mathbf{x}_{1:T}$ and $\mathbf{z}_{1:T}$ are stochastic. Therefore, as is commonly observed in the DVAE literature (Krishnan et al., 2015; Fraccaro et al., 2016; Fraccaro et al., 2017), we will only focus on modeling the distribution $p_{\theta}(\mathbf{x}_{1:T}, \mathbf{z}_{1:T}|\mathbf{u}_{1:T})$.

In the following section, we will first omit $\theta$ when defining the general structure of dependencies in the generative model. We will specify the parameter notation later when introducing how RNNs are used to parametrize the model. In addition, we will consider the model in the driven mode (i.e., with $\mathbf{u}_{1:T}$ as input) as it is more general than that in the undriven mode (i.e., with no “external” input). The undriven mode equations can be obtained from the driven mode equations by simply removing $\mathbf{u}_{1:T}$.

### 4.1.1 Structure of dependencies in the generative model

As we will discuss in detail in Section 14.3, a DVAE can be considered a structured or hierarchical VAE in which both observed and latent variables are a set of ordered vectors, and the ordering is imposed by time. However, the natural order present in the data does not imply a unique possible structure of variable dependencies for a DVAE generative (or inference) model. In fact, in DVAEs, the joint distribution of the observed and latent vector sequences is usually defined using the chain rule; that is, it is written as a product of conditional distributions over the vectors at different time indices. When writing the chain rule, different orderings of the random vectors can be arbitrarily chosen. This is an important point because the choice of ordering when applying the chain rule yields different practical implementations, which result in different sampling processes.
A natural choice for ordering dependencies at generation is to use a *causal* model. In the present context, a generation (or inference) model is said to be causal if the distribution of the generated (or inferred) variable at time *t* depends only on its values at previous time indices and/or on the values of the other variables at time *t* and at previous time indices. If the dependency is only over future time indices, the model is said to be *anticausal*, and if the dependency combines the past, present, and future of the conditioning variables, the model is said to be *noncausal*.

Let us consider the following simple example:

$$p(\mathbf{x}_1, \mathbf{x}_2, \mathbf{z}_1, \mathbf{z}_2) = p(\mathbf{x}_2|\mathbf{x}_1, \mathbf{z}_1)p(\mathbf{z}_2|\mathbf{x}_1, \mathbf{z}_1)p(\mathbf{x}_1|\mathbf{z}_1)p(\mathbf{z}_1) \quad (4.1)$$

$$= p(\mathbf{x}_2|\mathbf{x}_1, \mathbf{z}_1, \mathbf{z}_2)p(\mathbf{x}_1|\mathbf{z}_1, \mathbf{z}_2)p(\mathbf{z}_2|\mathbf{z}_1)p(\mathbf{z}_1). \quad (4.2)$$

In (4.1), the sampling is causal because we alternate between sampling $\mathbf{z}_t$ and $\mathbf{x}_t$ from their past value or their past and present values, from $t=1$ to 2. In contrast, in (4.2), the sampling is not causal because we first have to sample the complete sequence of latent vectors $\mathbf{z}_{1:2}$ before sampling $\mathbf{x}_1$, and then $\mathbf{x}_2$. This principle generalizes to much longer sequences.

In the DVAE literature, causal modeling is the most popular approach. In what follows, we will therefore focus on causal modeling, but the general methodology is similar for noncausal modeling. To the best of our knowledge, only one noncausal model has been proposed in the literature: the RVAE model (Leglaive et al., 2020). In fact, both causal and noncausal versions of RVAE were proposed in this paper, and both versions will be presented in Section 10.

In (causal) DVAEs, the joint distribution of the latent and observed sequences is first factorized according to the time indices using the chain rule:

$$p(\mathbf{x}_{1:T}, \mathbf{z}_{1:T}|\mathbf{u}_{1:T}) = \prod_{t=1}^{T} p(\mathbf{x}_t, \mathbf{z}_t|\mathbf{x}_{1:t-1}, \mathbf{z}_{1:t-1}, \mathbf{u}_{1:t}). \quad (4.3)$$

The only assumption made in (4.3) is the causal dependence of $\mathbf{x}_t$ and $\mathbf{z}_t$ on the input sequence $\mathbf{u}_{1:T}$. Then, at each time index $p(\mathbf{x}_t, \mathbf{z}_t|\mathbf{x}_{1:t-1}, \mathbf{z}_{1:t-1}, \mathbf{u}_{1:t})$
is again factorized using the chain rule, so that

$$p(\mathbf{x}_{1:T}, \mathbf{z}_{1:T} | \mathbf{u}_{1:T}) = \prod_{t=1}^{T} p(\mathbf{x}_t | \mathbf{x}_{1:t-1}, \mathbf{z}_{1:t}, \mathbf{u}_{1:t}) p(\mathbf{z}_t | \mathbf{x}_{1:t-1}, \mathbf{z}_{1:t-1}, \mathbf{u}_{1:t}). \quad (4.4)$$

This equation is a generalization of (4.1), and again, it exhibits the alternate sampling of $\mathbf{z}_t$ and $\mathbf{x}_t$. Similarly to our remark in Section 3.2.1, for $t = 1$, the first terms of the products in (4.3) and (4.4) are $p(\mathbf{x}_1, \mathbf{z}_1 | \mathbf{u}_1)$ and $p(\mathbf{x}_1 | \mathbf{z}_1, \mathbf{u}_1) p(\mathbf{z}_1 | \mathbf{u}_1)$, respectively. This is consistent with our notation choice of $\mathbf{x}_{1:0} = \mathbf{z}_{1:0} = \emptyset$. Alternatively, we can define $\mathbf{z}_0$ as the initial state vector and consider $p(\mathbf{z}_0)$, $p(\mathbf{z}_1 | \mathbf{z}_0, \mathbf{u}_1)$, and so on, in these equations. Hereinafter, for each detailed model, we will present the joint distribution in the general form of a product over frames from $t = 1$ to $T$, and for conciseness, will not detail the model “initialization.”

As will be detailed later, the different models proposed in the literature make different conditional independence assumptions to simplify the dependencies in the conditional distributions of (4.4). For instance, the SSM family presented in Section 3.2 is based on the following conditional independence assumptions:

$$p(\mathbf{x}_t | \mathbf{x}_{1:t-1}, \mathbf{z}_{1:t}, \mathbf{u}_{1:t}) = p(\mathbf{x}_t | \mathbf{z}_t), \quad (4.5)$$

$$p(\mathbf{z}_t | \mathbf{x}_{1:t-1}, \mathbf{z}_{1:t-1}, \mathbf{u}_{1:t}) = p(\mathbf{z}_t | \mathbf{z}_{t-1}, \mathbf{u}_t). \quad (4.6)$$

We have already introduced the concept of the driven mode. In the causal context, we say that a DVAE is in the driven mode if $\mathbf{u}_{1:t}$ is used to generate either $\mathbf{x}_{1:t}$, $\mathbf{z}_{1:t}$, or both. A DVAE is in *predictive mode* if $\mathbf{x}_{1:t-1}$, or part of this sequence, typically $\mathbf{x}_{t-1}$, is used to generate either $\mathbf{x}_t$ or $\mathbf{z}_t$, or both. This corresponds to feedback or closed-loop control in control theory. This is also strongly related to the concept of *autoregressive process*, jointly found in the control theory, machine learning, signal processing, or time-series analysis literature (Papoulis, 1977; Frey, 1998; Durbin and Koopman, 2012; Hamilton, 2020). Therefore, in what follows, we indifferently use the terms *predictive DVAE* or *autoregressive DVAE* to qualify a DVAE in the predictive mode.

In its most general form (4.4), a DVAE is both in the driven and predictive modes; however, it can also be in only one of the two modes
(e.g., the above SSM is in the driven mode but not in the predictive mode), or even in none of them. In the literature, we did not encounter any DVAE in both modes at the same time. Moreover, there are models in the driven and nonpredictive modes that are converted to the undriven and predictive modes by replacing the control input $u_t$ with the previously generated output $x_{t-1}$, see (Fraccaro et al., 2016). Note that a model's behavior can be quite different under the various modes. This is consistent with the concept of using a model in an open loop or in a closed loop in control theory. The principle of these different modes has been poorly discussed in the DVAE literature, and it is interesting to clarify it at an early stage of the DVAE presentation.

### 4.1.2 Parameterization with (R)NNs

The factorization in (4.4) is a general umbrella for all (causal) DVAEs. As discussed above, each DVAE model will make different conditional independence assumptions, which will simplify the general factorization in various ways. Once the conditional assumptions are made, one can easily determine if there is a need to accumulate the past information (e.g., $z_t$ or $x_t$ depends on past observations $x_{1:t}$) or if a first-order Markovian relationship holds (e.g., $z_t$ and $x_t$ depend at most on $z_{t-1}$ and $x_{t-1}$). Usually, the former is implemented using RNNs, whereas feed-forward DNNs can be used to implement first-order Markovian dependency. Moreover, once the conditional assumptions are made, the remaining dependencies can be implemented in different ways. Therefore, the final family of distributions depends not only on the conditional independence assumptions but also on the networks that are used to implement the remaining dependencies.

Let us showcase this with a concrete example in which we have the following conditional independence assumptions:

$$p(z_t|\mathbf{x}_{1:t-1}, \mathbf{z}_{1:t-1}, \mathbf{u}_{1:t}) = p(z_t|\mathbf{x}_{1:t-1}, \mathbf{u}_t), \quad (4.7)$$

$$p(\mathbf{x}_t|\mathbf{x}_{1:t-1}, \mathbf{z}_{1:t}, \mathbf{u}_{1:t}) = p(\mathbf{x}_t|\mathbf{x}_{1:t-1}, \mathbf{z}_t). \quad (4.8)$$

Here, we assume that the generation of both $x_t$ and $z_t$ depends on $x_{1:t-1}$. In addition, the generation of $x_t$ also depends on $z_t$ and that of $z_t$ also depends on $u_t$. To accumulate the information of all past outputs
<figure><img src="image_7.png" /><figcaption>**Figure 4.1:** Two different implementations of a given factorization. The probabilistic graphical model (left) shows the dependencies between random variables and corresponds to the factorization in (4.7) and (4.8). Two possible implementations based on RNNs are shown: sharing the internal state variables (middle) or with two different internal state variables (right). We refer to the *compact* representation (left) and to the *developed* representations (middle and right). This terminology holds true for both the graphical representations and model formulations.</figcaption></figure>

$x_{1:t-1}$, one can use an RNN. In practice, the past information is accumulated in the internal state variable of the RNN, namely $h_t$, computed recurrently at each frame $t$. Among the many possible implementations, we consider two in this example: in the first implementation, illustrated in Figure 4.1 (middle), a single RNN internal state variable $h_t$ is used to generate both $x_t$ and $z_t$, while in the second implementation, illustrated in Figure 4.1 (right), two different internal state variables, $h_t$ and $k_t$, are used to generate $x_t$ and $z_t$ separately.

Assuming that all probability distributions are Gaussian, the first implementation can be expressed as

$$h_t = d_h(x_{t-1}, h_{t-1}; \theta_h), \quad (4.9)$$

$$[\mu_{\theta_z}(x_{1:t-1}, u_t), \sigma_{\theta_z}(x_{1:t-1}, u_t)] = d_z(h_t, u_t; \theta_{hz}), \quad (4.10)$$

$$p_{\theta_z}(z_t|\mathbf{x}_{1:t-1}, \mathbf{u}_t) = \mathcal{N}(z_t; \mu_{\theta_z}(\mathbf{x}_{1:t-1}, \mathbf{u}_t), \mathrm{diag}\{\sigma_{\theta_z}^2(\mathbf{x}_{1:t-1}, \mathbf{u}_t)\}), \quad (4.11)$$

$$[\mu_{\theta_x}(\mathbf{x}_{1:t-1}, \mathbf{z}_t), \sigma_{\theta_x}(\mathbf{x}_{1:t-1}, \mathbf{z}_t)] = d_x(h_t, z_t; \theta_{hx}), \quad (4.12)$$

$$p_{\theta_x}(\mathbf{x}_t|\mathbf{x}_{1:t-1}, \mathbf{z}_t) = \mathcal{N}(\mathbf{x}_t; \mu_{\theta_x}(\mathbf{x}_{1:t-1}, \mathbf{z}_t), \mathrm{diag}\{\sigma_{\theta_x}^2(\mathbf{x}_{1:t-1}, \mathbf{z}_t)\}), \quad (4.13)$$

where $d_h$, $d_z$, and $d_x$ are nonlinear functions implemented with DNNs. It is now clear that the parameters of the conditional distribution of
$\mathbf{z}_t$ are $\theta_z = \theta_h \cup \theta_{hz}$, whereas those of the conditional distribution of $\mathbf{x}_t$ are $\theta_x = \theta_h \cup \theta_{hx}$. Thus, the two conditional distributions share the recurrent parameters $\theta_h$. Regarding the second implementation, the generative process can be expressed as

$$ \mathbf{h}_t = d_h(\mathbf{x}_{t-1}, \mathbf{h}_{t-1}; \theta_h), \quad (4.14) $$

$$ [\mu_{\theta_z}(\mathbf{x}_{1:t-1}, \mathbf{u}_t), \sigma_{\theta_z}(\mathbf{x}_{1:t-1}, \mathbf{u}_t)] = d_z(\mathbf{h}_t, \mathbf{u}_t; \theta_{hz}), \quad (4.15) $$

$$ p_{\theta_z}(\mathbf{z}_t|\mathbf{x}_{1:t-1}, \mathbf{u}_t) = \mathcal{N}(\mathbf{z}_t; \mu_{\theta_z}(\mathbf{x}_{1:t-1}, \mathbf{u}_t), \text{diag}\{\sigma_{\theta_z}^2(\mathbf{x}_{1:t-1}, \mathbf{u}_t)\}), \quad (4.16) $$

$$ \mathbf{k}_t = d_k(\mathbf{x}_{t-1}, \mathbf{k}_{t-1}; \theta_k), \quad (4.17) $$

$$ [\mu_{\theta_x}(\mathbf{x}_{1:t-1}, \mathbf{z}_t), \sigma_{\theta_x}(\mathbf{x}_{1:t-1}, \mathbf{z}_t)] = d_x(\mathbf{k}_t, \mathbf{z}_t; \theta_{kx}), \quad (4.18) $$

$$ p_{\theta_x}(\mathbf{x}_t|\mathbf{x}_{1:t-1}, \mathbf{z}_t) = \mathcal{N}(\mathbf{x}_t; \mu_{\theta_x}(\mathbf{x}_{1:t-1}, \mathbf{z}_t), \text{diag}\{\sigma_{\theta_x}^2(\mathbf{x}_{1:t-1}, \mathbf{z}_t)\}). \quad (4.19) $$

We have an additional DNN-based nonlinear function $d_k$, and analogously, it is clear that the parameters of the conditional distribution of $\mathbf{z}_t$ are $\theta_z = \theta_h \cup \theta_{hz}$, whereas those of the conditional distribution of $\mathbf{x}_t$ are $\theta_x = \theta_h \cup \theta_{hx}$. In this case, the two conditional distributions do not share any parameter. To ease the notation, hereinafter, we will denote the parameters of $p(\mathbf{z}_t|\mathbf{x}_{1:t-1}, \mathbf{z}_{1:t-1}, \mathbf{u}_{1:t})$ and $p(\mathbf{x}_t|\mathbf{x}_{1:t-1}, \mathbf{z}_{1:t}, \mathbf{u}_{1:t})$ as $\theta_z$ and $\theta_x$, respectively, in (4.4), irrespective of whether or not they share some parameters. We will also use $\theta$ to denote $\theta_z \cup \theta_x$.

In the equations above, the operators $d_h$, $d_k$, $d_x$ and $d_z$ are nonlinear mappings parametrized by DNNs of arbitrary architecture. How to choose and design these architectures is beyond the scope of this paper, as it largely depends on the target application. To fix ideas, in the present example, $d_h$ and $d_k$ are RNNs, and $d_x$ and $d_z$ are feed-forward DNNs. In this paper, we will not discuss how to select the hyper-parameters of these networks, such as the number of layers, or the number of units per layer.

Note that (4.11) and (4.16) are exactly the same, meaning that the conditional distributions of $\mathbf{z}_t$ are the same for both models. The same remark holds for (4.13) and (4.19), defining the conditional distribution of $\mathbf{x}_t$. However, the computations performed to obtain the parameters $\theta_z$ and $\theta_x$ differ depending on the model. Clearly, we need to make a distinction between the two forms. We propose to call the form of a DVAE model or its graphical representation *compact* when only random variables appear (e.g., (4.11) and (4.13), and Figure 4.1 (left)).
In addition, we propose to call the form of a DVAE or its graphical representation *developed* when both random and deterministic variables appear (e.g., (4.9)–(4.13), (4.14)–(4.19) and Figure 4.1 (middle) and (right)). Each compact form can have different developed forms corresponding to different implementations. The distinction between the compact and developed forms is important as the optimization occurs on the parameters of the developed form, which is only a subgroup of all possible models satisfying the compact form. It is thus important to present the developed form of a model. However, the temporal dependencies of order higher than one are not directly visible in the developed graphical form, as they might be implicitly encoded in the internal state variables. Therefore, when reviewing DVAE models in the following chapters, we will always present both the compact and developed graphical representations.

## 4.2 Inference model

In the present DVAE context, the posterior distribution of the state sequence $z_{1:T}$ is $p_\theta(z_{1:T}|x_{1:T}, u_{1:T})$ in the driven mode or $p_\theta(z_{1:T}|x_{1:T})$ in the undriven mode. As for the standard VAE, this posterior distribution is intractable because of the nonlinearities in the generative model. In fact, having temporal dependencies only makes things even more complicated. Therefore, we also need to define an inference model $q_\phi(z_{1:T}|x_{1:T}, u_{1:T})$, which is an approximation of the intractable posterior distribution $p_\theta(z_{1:T}|x_{1:T}, u_{1:T})$. As for the standard VAE, this model is required not only for performing inference of the latent sequence $z_{1:T}$ from the observed sequences $x_{1:T}$ and $u_{1:T}$ but also for estimating the parameters of the generative model, as will be seen below. As for the standard VAE again, the inference model also uses DNNs to generate its parameters.

### 4.2.1 Exploiting D-separation

In a Bayesian network, and in a DVAE in particular, even though the computation of the posterior distribution is often intractable, there exists a general methodology to express its general form (i.e., to specify the
dependencies between the variables of a generative model at inference time). This methodology is based on the so-called *D-separation* property of Bayesian networks (Geiger et al., 1990; Bishop, 2006, Chapter 8). The general principle is that some of the conditioning variables in the expression of the posterior distribution of a given variable can vanish depending on whether the nodes between these conditioning variables and the given variable represent variables that are observed or unobserved and depending on the direction of the dependencies (i.e., the direction of the arrows of the graphical representation).

In detail, D-separation is based on the three principles derived for a Bayesian network with three random variables *a*, *b*, and *c*:

* A *tail-to-tail* (or common parent) node *c* corresponding to the structure *a* ← *c* → *b* makes the two other nodes *a* and *b* conditionally independent when it is observed. In short, we have $p(a, b|c) = p(a|c)p(b|c)$.

* A *head-to-tail* (or cascade) node *c* corresponding to the structure *a* → *c* → *b* or *a* ← *c* ← *b* makes the two other nodes *a* and *b* conditionally independent when it is observed. In short, we have $p(a, b|c) = p(a|c)p(b|c)$.

* A *head-to-head* (or V-structure) node *c* corresponding to the structure *a* → *c* ← *b* makes *a* and *b* conditionally dependent when it is observed, hence $p(a, b|c) \neq p(a|c)p(b|c)$.

D-separation consists in applying these three principles recursively to analyze larger Bayesian networks with any arbitrary structure. Let us consider a Bayesian network in which *A*, *B*, and *C* are arbitrary nonintersecting node sets. *A* and *B* are *D*-separated given *C* if all possible paths that connect any node in *A* to any node in *B* are blocked given *C*. A path is said to be *blocked* given a set of observed nodes *O* if it includes a node *c* such that either

* *c* is a tail-to-tail node and $c \in O$ (i.e., it is observed) or

* *c* is a head-to-tail node and $c \in O$ (i.e., it is observed) or

* *c* is a head-to-head node and $c \notin O$ (i.e., it is not observed).
Equivalently, A and B are D-separated given C if they are not connected by any active path (i.e., a path that is not blocked). Finally, if A and B are D-separated given C, we have $p(A, B|C) = p(A|C)p(B|C)$.

D-separation is helpful even for more conventional (i.e., nondeep) models because the algebraic derivation of a posterior distribution from a joint distribution is not always easy. In the present variational framework, we can exploit the above methodology to design the approximate posterior distribution $q_\phi$. It is reasonable to assume that a good candidate for $q_\phi$ will have the same structure as the exact posterior distribution in terms of variable dependency. In other words, if we cannot derive the exact posterior distribution, let us at least use an approximation that exhibits the same dependencies between variables so that it is fed with the same information. Yet, it is quite surprising to see that a significant proportion of the DVAE papers we have reviewed, especially the early papers, neither refer to this methodology nor consider looking at the form of the exact posterior distribution when designing an approximate distribution. In the early studies in particular, the formulation of $q_\phi$ is chosen quite arbitrarily and with no reference to the structure of the exact posterior distribution. In more recent papers however, the structure of $q_\phi$ generally follows that of the exact posterior distribution. We will come back on this point on a case-by-case basis when presenting the DVAE models of the literature in the next chapters.

### 4.2.2 Noncausal and causal inference

Being aware of this problem, we can now go back to the general form of the exact posterior distribution and factorize it as follows, applying again the chain rule the same way as we did for the generative model:

$$p_{\theta}(\mathbf{z}_{1:T}|\mathbf{x}_{1:T}, \mathbf{u}_{1:T}) = \prod_{t=1}^{T} p_{\theta_z}(\mathbf{z}_t|\mathbf{z}_{1:t-1}, \mathbf{x}_{1:T}, \mathbf{u}_{1:T}). \quad (4.20)$$

For the most general generative model defined in (4.4), the dependencies in each conditional distribution $p_{\theta_z}(\mathbf{z}_t|\mathbf{z}_{1:t-1}, \mathbf{x}_{1:T}, \mathbf{u}_{1:T})$ cannot be simplified. In other words, $\mathbf{z}_t$ depends on the past latent vectors $\mathbf{z}_{1:t-1}$ and on the complete sequences of observed vectors $\mathbf{x}_{1:T}$ and $\mathbf{u}_{1:T}$ (past,
current, and future time steps). The exact inference is thus a noncausal process, even if the generation is causal. This is reminiscent of the Kalman smoother (i.e., the noncausal solution to inference in LG-LDS, see Section 3.2.2). As discussed in the previous subsection, the inference model $q_{\phi}$ should here have the same most general structure as the exact posterior distribution of (4.20):

$$ q_{\phi}(\mathbf{z}_{1:T}|\mathbf{x}_{1:T}, \mathbf{u}_{1:T}) = \prod_{t=1}^{T} q_{\phi}(\mathbf{z}_t|\mathbf{z}_{1:t-1}, \mathbf{x}_{1:T}, \mathbf{u}_{1:T}). \quad (4.21) $$

Similar to the generative model, each conditional posterior distribution $q_{\phi}(\mathbf{z}_t|\mathbf{z}_{1:t-1}, \mathbf{x}_{1:T}, \mathbf{u}_{1:T})$ should accumulate information from past latent variables and past observations, but in contrast to the generative model, it should also accumulate information from present and future observations. Typically, this process is implemented with a bidirectional recurrent network.

Depending on the conditional independence assumptions made when defining the generative model, the posterior dependencies in $p_{\theta}(\mathbf{z}_t|\mathbf{z}_{1:t-1}, \mathbf{x}_{1:T}, \mathbf{u}_{1:T})$ can be simplified using the D-separation property of Bayesian networks described in the previous subsection. Thus, the posterior dependencies in $q_{\phi}(\mathbf{z}_t|\mathbf{z}_{1:t-1}, \mathbf{x}_{1:T}, \mathbf{u}_{1:T})$ can be simplified similarly. Of course, it is always possible to use an approximate posterior $q_{\phi}$ that does not follow the structure of the exact posterior distribution. In fact, it makes sense to use a simplified version if one wants to decrease the computational cost or satisfy other constraints. In particular, for online or incremental data processing, the inference can be forced to be a causal process by removing the dependencies of $q_{\phi}$ on the future observations (and future inputs). This is similar to the Kalman filter for an LG-LDS, see again Section 3.2.2. This will generally be at risk of degrading the inference performance. Again, we will return to these points when reviewing the DVAE models proposed in the literature.

### 4.2.3 Sharing variables and parameters at generation and inference

We can note a similarity between the (most general causal) generative distribution $p_{\theta_z}(\mathbf{z}_t|\mathbf{x}_{1:t-1}, \mathbf{z}_{1:t-1}, \mathbf{u}_{1:t})$ and the corresponding inference model $q_{\phi}(\mathbf{z}_t|\mathbf{z}_{1:t-1}, \mathbf{x}_{1:T}, \mathbf{u}_{1:T})$ in terms of random variable dependencies.
For instance, the general form of the dependency of $z_t$ on past latent vectors is the same at inference and generation: in both cases, $z_t$ depends on the complete past sequence $z_{1:t-1}$. Implementing this recurrence at inference and at generation can be made either with a single unique RNN or with two different RNNs, in line with what we discussed in Section 4.1.2. The same principle applies to $u_{1:t}$ and $x_{1:t}$, which are both used at generation and inference. Depending on which variables we consider, it can make sense to use the same RNN at generation and inference, meaning that the deterministic link between the realizations of random variables is the same at generation and at inference. If this is the case, the decoder and encoder share some network modules and thus $\theta$ and $\phi$ share some parameters. Note that this is not the case in standard VAEs.

Hereinafter, we will use $h_t$ to denote the internal state of the decoder and $g_t$ to denote that of the encoder if it is different from the internal state of the decoder. Otherwise, we will use $h_t$ for the encoder as well.

## 4.3 VLB and training of DVAEs

As for the standard VAE, training a DVAE is based on the maximization of the VLB. In the case of DVAEs, the VLB initially defined in (2.13) is extended to data sequences as follows:

$$ \mathcal{L}(\theta, \phi; \mathbf{x}_{1:T}, \mathbf{u}_{1:T}) = \mathbb{E}_{q_{\phi}(\mathbf{z}_{1:T}|\mathbf{x}_{1:T}, \mathbf{u}_{1:T})} [\log p_{\theta}(\mathbf{x}_{1:T}, \mathbf{z}_{1:T}|\mathbf{u}_{1:T})] \\ - \mathbb{E}_{q_{\phi}(\mathbf{z}_{1:T}|\mathbf{x}_{1:T}, \mathbf{u}_{1:T})} [\log q_{\phi}(\mathbf{z}_{1:T}|\mathbf{x}_{1:T}, \mathbf{u}_{1:T})]. \quad (4.22) $$

With the factorization in (4.21), the expectation in (4.22) can be expressed as a cascade of expectations taken with respect to conditional distributions over individual latent vectors at different time indices:

$$ \mathbb{E}_{q_{\phi}(\mathbf{z}_{1:T}|\mathbf{x}_{1:T}, \mathbf{u}_{1:T})} [\psi(\mathbf{z}_{1:T})] = \mathbb{E}_{q_{\phi}(\mathbf{z}_1|\mathbf{x}_{1:T}, \mathbf{u}_{1:T})} \left[ \mathbb{E}_{q_{\phi}(\mathbf{z}_2|\mathbf{z}_1, \mathbf{x}_{1:T}, \mathbf{u}_{1:T})} \cdots \right. \\ \left. \mathbb{E}_{q_{\phi}(\mathbf{z}_T|\mathbf{z}_{1:T-1}, \mathbf{x}_{1:T}, \mathbf{u}_{1:T})} [\psi(\mathbf{z}_{1:T})] \cdots \right] \right], \quad (4.23) $$

where $\psi(\mathbf{z}_{1:T})$ denotes any function of $\mathbf{z}_{1:T}$. Then, by injecting (4.4) and (4.21) into (4.22), and using the above cascade, we can develop the
VLB as follows:

$$
\begin{align}
\mathcal{L}(\theta, \phi; \mathbf{x}_{1:T}, \mathbf{u}_{1:T}) &= \mathbb{E}_{q_{\phi}(\mathbf{z}_{1:T}|\mathbf{x}_{1:T}, \mathbf{u}_{1:T})} \left[ \log p_{\theta}(\mathbf{x}_{1:T}, \mathbf{z}_{1:T}|\mathbf{u}_{1:T}) \right. \nonumber \\
&\qquad \left. - \log q_{\phi}(\mathbf{z}_{1:T}|\mathbf{x}_{1:T}, \mathbf{u}_{1:T}) \right] \nonumber \\
&= \sum_{t=1}^{T} \mathbb{E}_{q_{\phi}(\mathbf{z}_{1:t}|\mathbf{x}_{1:T}, \mathbf{u}_{1:T})} \left[ \log p_{\theta_x}(\mathbf{x}_t|\mathbf{x}_{1:t-1}, \mathbf{z}_{1:t}, \mathbf{u}_{1:t}) \right] \nonumber \\
&\quad - \sum_{t=1}^{T} \mathbb{E}_{q_{\phi}(\mathbf{z}_{1:t-1}|\mathbf{x}_{1:T}, \mathbf{u}_{1:T})} \left[ D_{\text{KL}}(q_{\phi}(\mathbf{z}_t|\mathbf{z}_{1:t-1}, \mathbf{x}_{1:T}, \mathbf{u}_{1:T}) \parallel p_{\theta_z}(\mathbf{z}_t|\mathbf{x}_{1:t-1}, \mathbf{z}_{1:t-1}, \mathbf{u}_{1:t})) \right]. \tag{4.24}
\end{align}
$$

To the best of our knowledge, this is the first time that the VLB is
presented in this most general form, which is valid for the entire class
of (causal) DVAE models.

As for the standard VAE, the VLB contains a reconstruction ac-
curacy term and a regularization term. However, in contrast to the
standard VAE, where the regularization term has an analytical form
for usual distributions, here, both the reconstruction accuracy and
regularization term require the computation of Monte Carlo estimates
(i.e., empirical averages) using samples drawn from $q_{\phi}(\mathbf{z}_{1:\tau}|\mathbf{x}_{1:T}, \mathbf{u}_{1:T})$,
where $\tau \in \{1, ..., T\}$ is an arbitrary index. Using the chain rule in (4.21),
we sample from the joint distribution $q_{\phi}(\mathbf{z}_{1:\tau}|\mathbf{x}_{1:T}, \mathbf{u}_{1:T})$ by sampling
recursively from $q_{\phi}(\mathbf{z}_t|\mathbf{z}_{1:t-1}, \mathbf{x}_{1:T}, \mathbf{u}_{1:T})$, going from $t=1$ to $t=\tau$.
Sampling each random vector $\mathbf{z}_t$ at a given time instant is straightfor-
ward, as $q_{\phi}(\mathbf{z}_t|\mathbf{z}_{1:t-1}, \mathbf{x}_{1:T}, \mathbf{u}_{1:T})$ is analytically specified by the chosen
inference model (e.g., Gaussian with mean and variance provided by an
RNN). We have to use a similar reparameterization trick as for standard
VAEs, so the sampling-based VLB estimator remains differentiable with
respect to $\phi$. The VLB can then be maximized with respect to both $\phi$
and $\theta = \theta_z \cup \theta_x$ using gradient-ascent-based algorithms. We recall that
for DVAEs, $\phi$ and $\theta$ can share parameters, which is different from the
"static" VAE, but perfectly alright for the optimization. Finally, the VLB
is here defined here for a single data sequence, but a common practice is
to average the VLB over a mini-batch of training data sequences before
updating the model parameters with gradient ascent.
## 4.4 Additional dichotomy for autoregressive DVAE models

A DVAE can be used to generate new data, for analysis-synthesis (by chaining the encoder and decoder), or for data transformation, by modifying the latent vector sequence in between analysis and synthesis. In the case of DVAE models functioning in the predictive mode (i.e., autoregressive DVAEs, see Section 4.1.1), these tasks can be processed in different manners, leading to an additional dichotomy of functioning modes. We describe these functioning modes in the next subsection before we see the implications for model training in the following subsection. Because these additional different modes concern the recursive part of the models, nonpredictive DVAEs are not concerned here.

### 4.4.1 Teacher forcing against generation mode

In practice, for autoregressive DVAE models, we have two generation modes, for the generation of both $x_t$ and $z_t$. A mode in which we assume that the ground-truth past observed vectors $x_{1:t-1}$ are used for generating the current vector ($x_t$ or $z_t$), and a mode in which the generated past observed vectors are used for generating the current vector. At this point, it is important to distinguish between the notation for the ground-truth value of the observed data vector $x_t$ and that for its modeled version produced by a DVAE, which we denote by $\hat{x}_t$. In practice, $x_{1:T}$ is a given data sequence that we want to model with a DVAE (or that we use for model training, as shown below), and $\hat{x}_{1:T}$ is the actual output of the DVAE.

This issue of either using the ground-truth past observed data vectors $x_{1:t-1}$ or reinjecting the previously generated vectors $\hat{x}_{1:t-1}$ at the input of a generative model is a classical problem of recursive models and, in particular, of RNNs. Yet it is poorly discussed in the DVAE literature. In the RNN literature, the first configuration is sometimes referred to as *teacher-forcing* (Williams and Zipser, 1989), as it is assumed that a teacher (or oracle) can provide the model with ground-truth values. Hereinafter, we will use this terminology. This is a classical configuration at training time, when the whole sequence $x_{1:T}$ is available and the model is tuned so that $\hat{x}_{1:T}$ fits $x_{1:T}$. However, this is unrealistic at generation
time, when the model produces a new sequence $\hat{x}_{1:T}$. Here, the second configuration must be used. We refer to this second configuration as the *generation mode*. Regarding the generation of $z_{1:T}$, the teacher-forcing concept is irrelevant as the concept of ground-truth values for latent vectors is questionable in essence. In practice, $z_{1:T}$ is either the output of the inference model (this is the case during DVAE training or in analysis-synthesis) or any arbitrary latent vector sequence (generated with $p_{\theta_z}$ or predefined). Note also that because of the recursivity of the generative process, data generation with a DVAE (strongly) depends on the initialization of the generative process. We do not detail this aspect in the present review.

If we now focus on the analysis-synthesis task, we have to chain the encoder and decoder. The encoder takes $x_{1:T}$ as the input and produces a sequence of latent vectors $z_{1:T}$. Then, the decoder uses $z_{1:T}$ to generate $\hat{x}_{1:T}$. Decoding can be performed with either teacher-forcing or generation mode. The former case is expected to produce a sequence $\hat{x}_{1:T}$ that is closer to $x_{1:T}$ than in the latter case, as it uses ground-truth values, whereas the generation mode uses approximate values. However, it suffers from the same “unrealistic” aspect as that used for data generation. This configuration can be used to evaluate the prediction power of DVAE models in an ideal (oracle) setting. In contrast, analysis-synthesis with the generation mode is expected to yield lower performance but is the natural configuration from an information-theoretic viewpoint. Here, we test the capability of the model to encode the information of a (generally high-dimensional) data sequence $x_{1:T}$ into a (generally low-dimensional) latent sequence $z_{1:T}$. From an application point of view, this corresponds to telecommunication or storage applications, where $x_{1:T}$ would be encoded into $z_{1:T}$, $z_{1:T}$ would be transmitted or stored, and then $\hat{x}_{1:T}$ would be decoded from $z_{1:T}$. In short, we apply DVAEs to source coding, and the DVAE turns into a codec, apart from quantization issues. Such coding/decoding scheme can be applied offline by using a noncausal inference model (with an optimal structure following that of the exact posterior distribution) or online with a (sub-optimal) causal inference model. If some amount of latency is tolerated, one can also use a noncausal inference model with a suitable lookahead. The interest of DVAE models for source coding is further discussed in
Section 14.4.

### 4.4.2 Train/test matching

As mentioned in the previous subsection, teacher-forcing is a conventional strategy used for training recursive models. However, in practice, when using autoregressive DVAEs for data generation or compression, the generation mode must be used. This leads to a mismatch between the training and testing conditions, a general problem in machine learning that leads to performance degradation compared to the case in which the same configuration is used for training and testing. Therefore, if the generation mode is used in a practical DVAE use-case, it can be beneficial to use the generation mode during model training as well, so that the training configuration matches the practical use-case configuration. In practice, this implies replacing $x_{1:t-1}$ with $\hat{x}_{1:t-1}$ in the conditioning variables in the VLB equations of Section 4.3. In Chapter 13, we illustrate this strategy in our experimental benchmark. We observe in our experiments that using the generation mode at both model training and testing leads to a significant gain in performance compared to the mismatched configuration, though the performance remains slightly lower than that in the case where teacher-forcing is used at both training and testing. More details are given in Section 13.3.3.

## 4.5 DVAE summary

Dynamical VAEs are constructed with various stochastic relationships among the control variables $u_{1:T}$, latent variables $z_{1:T}$, and observed variables $x_{1:T}$. We recall that a random variable **a** is called a parent of another random variable **b** when the realization of **a** is used to compute the parameters of the distribution of **b**. These parameters can be obtained with a linear or a nonlinear mapping of the realization of **a** (and possibly of other random variables). A DVAE model must contain two types of relationships:

* **Decoding link:** $z_t$ is always a parent of $x_t$. Graphically, there is always an arrow from $z_t$ to $x_t$ in the compact graphical representation. This is a fundamental characteristic inherited from the
standard VAE.

* **Temporal link:** At least one element in $z_{1:t-1}$ or in $x_{1:t-1}$ is parent to either $z_t$ or $x_t$. One of the simplest forms of a temporal link, namely $z_{t-1}$ is a parent of $z_t$, is a fundamental characteristic of first-order SSMs.

In a way, the “minimal DVAE” is the straightforward combination of a first-order SSM and a VAE, which is the DKF model that we will detail in Section 5. Other DVAEs include additional temporal links. Moreover, temporal links such as “$z_{t-1}$ is a parent of $x_t$” can be considered additional decoding links, in that $x_t$ is generated from $z_t$ and $z_{t-1}$. As for temporal links, in the papers that we detail in this overview, $z_t$ and/or $x_t$ depend either on $z_{t-1}$ and/or $x_{t-1}$, or on $z_{1:t-1}$ and/or $x_{1:t-1}$. In other words, the order of temporal dependencies is either 1 (implemented with a basic feed-forward neural network, such as a Multi-Layer Perceptron (MLP)) or infinity (implemented with an RNN). However, one can, in principle, use $N$-order temporal dependencies with $1 < N < \infty$, relying, for instance, on convolutional neural networks (CNNs) with finite-length receptive fields. In particular, temporal convolutional networks (TCNs) (Lea et al., 2016), which are based on dilated convolutions, are competitive with RNNs on several sequence modeling tasks, including generative modeling (Aksan and Hilliges, 2019).

Finally, as discussed before, a DVAE can be in the *driven mode*, in the *predictive mode*, in both, or in none of these modes. This is also modeled by parenthood relationships that may or may not exist.

* **Driving link:** A DVAE is said to be in the driven mode if $u_t$ is a parent of either $z_t$ or $x_t$, or both.

* **Predictive link:** A DVAE is said to be in the predictive mode if $x_{1:t-1}$, or part of this sequence, is a parent of either $z_t$ or $x_t$, or both. In practice, a predictive link is implemented either in the teacher-forcing mode (using the ground-truth past vector sequence) or in the generation mode (using the previously generated vector sequence). Generally, better performance is achieved if the same mode is used at training time and test time.
<figure><img src="image_8.png" /></figure>
Continuing from the previous section, we start our DVAE tour by combining SSMs with neural networks. Such a combination is not recent, see, e.g., (Haykin, 2004; Raiko and Tornio, 2009), but has been recently investigated under the VAE angle in two papers by the same authors (Krishnan et al., 2015; Krishnan et al., 2017). The resulting deep SSM is referred to as a DKF (Krishnan et al., 2015) or a deep Markov model (DMM) (Krishnan et al., 2017).<fnref n="1" /> Therefore, these papers do not provide a new concept in terms of models, but they provide a solution to the joint problem of inference and model parameter estimation in the VAE methodological framework applied to the SSM model architecture. In other words, this is, to the best of our knowledge, the first example of unsupervised training of a deep SSM by chaining an approximate inference model with a generative model and using the VLB maximization methodology. This training leads to an unsupervised discovery of the latent space that encodes the temporal

<footnote n="1">¹The same generative model is considered in both papers, but as we will detail later, the second paper proposes notable improvements regarding the inference model. The authors change the model name from DKF to DMM, maybe because the second denomination appears more general. In the present review, we retain the denomination DKF.</footnote>
dynamics of the data. The VAE methodology circumvents the difficulties encountered in previous approaches (Haykin, 2004; Raiko and Tornio, 2009) concerning the computational complexity and practicability of model parameter estimation, particularly allowing to move directly from single-layer neural networks to DNNs.

Hereinafter, for all detailed DVAE models, we first present the generative equations, then the inference model (and we discuss its choice by referring to the exact posterior distribution structure as deduced from the generative model), and finally the detailed form of the VLB used for model training together with clues about the optimization algorithm.

## 5.1 Generative model

We have already seen the generative equations of the DKF, as they are the same as (3.9)–(3.12), with the specificity that $d_z$ and $d_x$ are DNNs here.<fnref n="2" /> Krishnan et al. (2015) did not specify these DNNs; we can assume that basic feed-forward neural networks (i.e., MLPs) were used. In their second paper, Krishnan et al. (2017) implemented $d_x$ with a two-layer MLP and used a slightly more refined model for $d_z$: a gated linear combination of a linear model and an MLP for the mean parameter, where the gate is itself provided by an MLP, and the chaining of MLP, rectified linear unit (ReLU) activation, and Softmax activation layers for the variance parameter (see Section B.1). According to the authors, this allows “the model have the flexibility to choose a linear transition for some dimensions while having a nonlinear transition for the others.”

Even if we are still at an early point in our presentation of the different DVAE models, we can make a first series of remarks to clarify the links between the models we have discussed so far.

* The stochastic state $z_t$ of an SSM is similar in essence to the latent state of the VAE. In the present DKF case, where $d_x$ is implemented

<footnote n="2">²In fact, a Bernoulli distribution was considered for $x_t$ by Krishnan et al. (2015) and Krishnan et al. (2017), but we have already mentioned that different pdfs can be considered for $p_{\theta_x}(x_t|z_t)$ depending on the data, without affecting the fundamental issues of this review. Hence, we consider here a Gaussian distribution for a better comparison with the other models.</footnote>
with a DNN, if $\mathbf{z}_t$ is of reduced dimension compared to $\mathbf{x}_t$, the DKF observation model is identical to the VAE decoder.

* Consequently, a DKF can be viewed as a VAE decoder with a temporal (Markovian) model of the latent variable $\mathbf{z}$.

* A (deep) SSM can also be viewed as a “fully stochastic” version of a (deep) RNN, where stochasticity is introduced at both the observation model level (like a GRNN) and internal state level. As mentioned before, in an SSM, the deterministic internal state $\mathbf{h}_t$ of the (G)RNN is simply replaced with a stochastic state $\mathbf{z}_t$.

* In summary, DKF = deep SSM = “Markovian” VAE decoder = “fully stochastic” RNN. The graphical model of the DKF is given by the right-hand schema in Figure 3.1 (which does not make the DNNs apparent).

## 5.2 Inference model

Following the general line of Section 4.2, we first identify the structure of the SSM/DKF posterior distribution $p_{\theta}(\mathbf{z}_{1:T}|\mathbf{x}_{1:T}, \mathbf{u}_{1:T})$. Let us first recall that applying the chain rule enables us to rewrite this distribution as follows:

$$p_{\theta}(\mathbf{z}_{1:T}|\mathbf{x}_{1:T}, \mathbf{u}_{1:T}) = \prod_{t=1}^{T} p_{\theta}(\mathbf{z}_{t}|\mathbf{z}_{t-1}, \mathbf{x}_{1:T}, \mathbf{u}_{1:T}). \quad (5.1)$$

Then, D-separation can be used to simplify each term of the product. The structure presented in Figure 3.1 (right) shows that the $\mathbf{z}_{t-1}$ node “blocks” all information coming from the past and flowing to $\mathbf{z}_t$ (i.e., $\mathbf{z}_{1:t-2}$, $\mathbf{x}_{1:t-1}$, and $\mathbf{u}_{1:t-1}$). In other words, $\mathbf{z}_{t-1}$ has accumulated this past information or is a summary of this information. We thus have $p_{\theta}(\mathbf{z}_t|\mathbf{z}_{t-1}, \mathbf{x}_{1:T}, \mathbf{u}_{1:T}) = p_{\theta}(\mathbf{z}_t|\mathbf{z}_{t-1}, \mathbf{x}_{t:T}, \mathbf{u}_{t:T})$, and therefore (with $\mathbf{z}_0$ being arbitrarily set)

$$p_{\theta}(\mathbf{z}_{1:T}|\mathbf{x}_{1:T}, \mathbf{u}_{1:T}) = \prod_{t=1}^{T} p_{\theta}(\mathbf{z}_{t}|\mathbf{z}_{t-1}, \mathbf{x}_{t:T}, \mathbf{u}_{t:T}). \quad (5.2)$$
At each time $t$, the posterior distribution of $z_t$ depends on the previous latent state $z_{t-1}$ and on the present and future observations $x_{t:T}$ and inputs $u_{t:T}$ (it is thus a first-order Markovian causal process on $z_t$ combined with an anticausal process on $x_t$ and $u_t$).

Krishnan et al. (2015) indicated this structure and the fact that we should inspire from it to design the approximate posterior $q_\phi$. However, somewhat surprisingly, they proposed the following four different models:

* an instantaneous model: $q_\phi(z_t|x_t, u_t)$ parameterized by an MLP;

* a model with local past and future context: $q_\phi(z_t|x_{t-1:t+1}, u_{t-1:t+1})$ parameterized by an MLP;

* a model with the complete past context: $q_\phi(z_t|x_{1:t}, u_{1:t})$ parameterized by an RNN;

* a model with the whole sequence: $q_\phi(z_t|x_{1:T}, u_{1:T})$ parameterized by a bidirectional RNN.

We do not detail these implementations of $q_\phi$ here, as we will provide other detailed examples in the next chapters. One might wonder why $z_{t-1}$ is not present in the conditional variables of the approximate posterior, but this might just be an oversight from the authors. This is difficult to be determined from the paper, as the implementation is not detailed. The authors mention an RNN to model the dependencies on $x$ and $u$, but the implementation of the dependency on $z_{t-1}$ is not specified. The point is that, although the authors pointed out the dependency of the exact posterior on the present and future observations and inputs, they did not propose a corresponding approximate model.

**Notation remark:** Krishnan et al. (2015) denoted by $u_{t-1}$ the input at time $t$ in the generative model, as in many control theory papers published on SSMs, and not $u_t$, as we do (as pointed out in a previous footnote). However, when defining the four approximate posterior models, they did it exactly as we report here (i.e., with $u_t$ being synchronous to $z_t$ and $x_t$). We conjecture that this problem is just a notation mistake made by Krishnan et al. (2015), which we have implicitly corrected by using $u_t$ instead of $u_{t-1}$ as the input at time $t$ in the generative model.
Krishnan et al. (2017) proposed the same generative model, renamed it DMM and presented it in the undriven mode, in which $\mathbf{u}_{1:T}$ was simply removed. However, they largely clarified and improved on their previous paper regarding the inference model. They proposed a new series of inference models that clearly do or do not depend on $\mathbf{z}_{t-1}$ and varied the dependency on the observed data. They again considered the case of dependency on the past and present data sequence $\mathbf{x}_{1:t}$ and on the complete data sequence $\mathbf{x}_{1:T}$. More importantly, they now also consider the case of an inference model with a functional form that corresponds exactly to the form of the exact posterior distribution, namely $q_{\phi}(\mathbf{z}_t|\mathbf{z}_{t-1}, \mathbf{x}_{t:T})$. In this case, for a complete data sequence, we have

$$q_{\phi}(\mathbf{z}_{1:T}|\mathbf{x}_{1:T}) = \prod_{t=1}^{T} q_{\phi}(\mathbf{z}_t|\mathbf{z}_{t-1}, \mathbf{x}_{t:T}). \qquad (5.3)$$

This model is referred to as the deep Kalman smoother (DKS), as it combines information from the past, through $\mathbf{z}_{t-1}$, and information from the present and future observations.

For conciseness, we report the detailed inference equations only for the DKS and will comment on their extension to the other proposed inference models. The DKS is implemented with a backward RNN on $\mathbf{x}_t$, followed by an additional layer for combining the RNN output with $\mathbf{z}_{t-1}$:

$$\overleftarrow{\mathbf{g}}_t = e_{\overleftarrow{\mathbf{g}}} (\overleftarrow{\mathbf{g}}_{t+1}, \mathbf{x}_t), \qquad (5.4)$$

$$\mathbf{g}_t = \frac{1}{2} (\tanh(\mathbf{W}\mathbf{z}_{t-1} + \mathbf{b}) + \overleftarrow{\mathbf{g}}_t), \qquad (5.5)$$

$$[\boldsymbol{\mu}_{\phi}(\mathbf{g}_t), \boldsymbol{\sigma}_{\phi}(\mathbf{g}_t)] = e_{\mathbf{z}}(\mathbf{g}_t), \qquad (5.6)$$

$$q_{\phi}(\mathbf{z}_t|\mathbf{g}_t) = \mathcal{N}(\mathbf{z}_t; \boldsymbol{\mu}_{\phi}(\mathbf{g}_t), \mathrm{diag}\{\boldsymbol{\sigma}_{\phi}^2(\mathbf{g}_t)\}). \qquad (5.7)$$

In the above equations, $e_{\mathbf{z}}$ is a basic combining network, parameterized by $\phi_{\mathbf{z}}$, $\boldsymbol{\mu}_{\phi}(\mathbf{g}_t)$ is an affine function of $\mathbf{g}_t$, and $\boldsymbol{\sigma}_{\phi}^2(\mathbf{g}_t)$ is a Softplus of an affine function of $\mathbf{g}_t$. We thus have $\phi = \phi_{\overleftarrow{\mathbf{g}}} \cup \phi_{\mathbf{z}}$, assuming $\{\mathbf{W}, \mathbf{b}\} \in \phi_{\mathbf{z}}$ for simplicity.

Because of the recursivity in (5.4), we can see $\mathbf{g}_t$ as an unfolded deterministic function of $\mathbf{z}_{t-1}$ and $\mathbf{x}_{t:T}$, which we can rewrite as $\mathbf{g}_t =$
$g_t(\mathbf{z}_{t-1}, \mathbf{x}_{t:T})$,<fnref n="3" /> and we have $q_{\phi}(\mathbf{z}_t|\mathbf{g}_t) = q_{\phi}(\mathbf{z}_t|\mathbf{g}_t(\mathbf{z}_{t-1}, \mathbf{x}_{t:T}))$. For a complete data sequence, we have

$$q_{\phi}(\mathbf{z}_{1:T}|\mathbf{x}_{1:T}) = \prod_{t=1}^{T} q_{\phi}(\mathbf{z}_t|\mathbf{g}_t) = \prod_{t=1}^{T} q_{\phi}(\mathbf{z}_t|\mathbf{g}_t(\mathbf{z}_{t-1}, \mathbf{x}_{t:T})), \quad (5.8)$$

which is just a rewriting of (5.3). In short, the inference of $\mathbf{z}_t$ with a DKS requires a first backward pass from $\mathbf{x}_T$ up to $\mathbf{x}_t$ to compute $\overleftarrow{\mathbf{g}}_t$, which is then combined with $\mathbf{z}_{t-1}$.

As mentioned above, this inference model is extended to a noncausal (bidirectional) model regarding the observations, $q_{\phi}(\mathbf{z}_t|\mathbf{z}_{t-1}, \mathbf{x}_{1:T})$. This is done by adding a forward RNN on $\mathbf{x}_t$ and sending its output $\overleftarrow{\mathbf{g}}_t$ to the combining network, in addition to $\mathbf{z}_{t-1}$ and $\overleftarrow{\mathbf{g}}_t$. For conciseness, we do not report the corresponding detailed equations, but this more general model is represented in Figure 5.1. The other models proposed by Krishnan et al. (2017) can be deduced from this general model by removing some elements. In particular, DKS is obtained by simply removing the forward RNN. Moreover, models that do not depend on $\mathbf{z}_{t-1}$ are obtained by removing the arrows between $\mathbf{z}_{t-1}$ and $\mathbf{g}_t$ for all $t$.

## 5.3 Training

A comparison of the compact form of the DKF model in (3.13) with the general compact form of a DVAE in (4.4) shows that the DKF model makes the following conditional independence assumptions:

$$\begin{aligned} p_{\theta_x}(\mathbf{x}_t|\mathbf{x}_{1:t-1}, \mathbf{z}_{1:t}, \mathbf{u}_{1:t}) &= p_{\theta_x}(\mathbf{x}_t|\mathbf{z}_t); \\ p_{\theta_z}(\mathbf{z}_t|\mathbf{x}_{1:t-1}, \mathbf{z}_{1:t-1}, \mathbf{u}_{1:t}) &= p_{\theta_z}(\mathbf{z}_t|\mathbf{z}_{t-1}, \mathbf{u}_t). \end{aligned} \quad (5.9)$$

Using these two simplifications along with the inference model (5.3) (extended to be in the driven mode for the sake of generality), the VLB

<footnote n="3">³This function is also a function of the initial state $\overleftarrow{\mathbf{g}}_T$ of the backward RNN.</footnote>
<figure><img src="image_9.png" /><figcaption>**Figure 5.1:** Graphical model of DKF at inference time corresponding to the inference model $q_{\phi}(z_t|z_{t-1}, x_{1:T})$, in developed form (left) and compact form (right). The specific DKS model, which functional form $q_{\phi}(z_t|z_{t-1}, x_{t:T})$ perfectly corresponds to the form of the exact posterior distribution, is obtained by removing the forward RNN (and removing the blue arrows on the right-hand schema).</figcaption></figure>
in its most general form (4.24) can be simplified as follows:

$$ \mathcal{L}(\theta, \phi; \mathbf{x}_{1:T}, \mathbf{u}_{1:T}) = \sum_{t=1}^{T} \mathbb{E}_{q_{\phi}(\mathbf{z}_t|\mathbf{x}_{1:T}, \mathbf{u}_{1:T})} \left[ \log p_{\theta_x}(\mathbf{x}_t|\mathbf{z}_t) \right] \\ - \sum_{t=1}^{T} \mathbb{E}_{q_{\phi}(\mathbf{z}_{t-1}|\mathbf{x}_{1:T}, \mathbf{u}_{1:T})} \left[ D_{\text{KL}} \left( q_{\phi}(\mathbf{z}_t|\mathbf{z}_{t-1}, \mathbf{x}_{t:T}, \mathbf{u}_{t:T}) \parallel p_{\theta_z}(\mathbf{z}_t|\mathbf{z}_{t-1}, \mathbf{u}_t) \right) \right]. \tag{5.10} $$

The KL divergence in (5.10) can be computed analytically, while the two expectations are intractable. (Krishnan et al., 2015; Krishnan et al., 2017) provided no detailed information regarding how to approximate these expectations; they only mentioned that “stochastic backpropagation” is used, referring the reader to the papers of Kingma and Welling (2014) and Rezende et al. (2014), who introduced the reparameterization trick for standard “static” VAEs. However, due to the dynamical nature of the model, the sampling procedure required for stochastic backpropagation in DKF is more complicated than in standard VAEs. In particular, we do not have an analytical form for $q_{\phi}(\mathbf{z}_t|\mathbf{x}_{1:T}, \mathbf{u}_{1:T})$, and only have one for $q_{\phi}(\mathbf{z}_{\tau}|\mathbf{z}_{\tau-1}, \mathbf{x}_{\tau:T}, \mathbf{u}_{\tau:T})$. Therefore, we need to exploit the chain rule and the “cascade trick” to develop and then approximate the intractable expectations in (5.10). The first expectation in this VLB expression can be developed as follows:

$$ \begin{aligned} \mathbb{E}_{q_{\phi}(\mathbf{z}_t|\mathbf{x}_{1:T}, \mathbf{u}_{1:T})}[f(\mathbf{z}_t)] &= \mathbb{E}_{q_{\phi}(\mathbf{z}_{1:t}|\mathbf{x}_{1:T}, \mathbf{u}_{1:T})}[f(\mathbf{z}_t)] \\ &= \mathbb{E}_{q_{\phi}(\mathbf{z}_1|\mathbf{x}_{1:T}, \mathbf{u}_{1:T})}[ \\ & \quad \mathbb{E}_{q_{\phi}(\mathbf{z}_2|\mathbf{z}_1, \mathbf{x}_{2:T}, \mathbf{u}_{2:T})}[ \dots \\ & \quad \mathbb{E}_{q_{\phi}(\mathbf{z}_t|\mathbf{z}_{t-1}, \mathbf{x}_{t:T}, \mathbf{u}_{t:T})}[f(\mathbf{z}_t)] \dots ], \end{aligned} \tag{5.11} $$

where $f(\mathbf{z}_t)$ denotes an arbitrary function of $\mathbf{z}_t$. A similar procedure can be used to develop the second expectation in (5.10). Each intractable expectation in this cascade of expectations can then be approximated with a Monte Carlo estimate. It requires to sample $q_{\phi}(\mathbf{z}_{\tau}|\mathbf{z}_{\tau-1}, \mathbf{x}_{\tau:T}, \mathbf{u}_{\tau:T})$ iteratively from $\tau = 1$ to $t$, using the same reparametrization trick as in standard VAEs. Doing so, the VLB becomes differentiable w.r.t. $\theta = \theta_x \cup \theta_z$ and $\phi = \phi_{\bar{g}} \cup \phi_z$, and it can be optimized with gradient-ascent-based techniques.
The KVAE model was presented by Fraccaro *et al.* (2017). This model can be considered a variant of the DKF model, and hence as another deep SSM, where an additional random variable, denoted $a_t$, is inserted between the latent vector $z_t$ and the observed vector $x_t$, as illustrated in Figure 6.1. This enables us to separate the model into two parts: A deep feature extractor linking $a_t$ and $x_t$ and the dynamical model on $z_t$ with “new observations” $a_t$. As we will see below, this provides the model with interesting properties for inference and training.
<figure><img src="image_10.png" /><figcaption>Figure 6.1: KVAE’s graphical model.</figcaption></figure>

## 6.1 Generative model

The general formulation of the KVAE model is given by

$$[\mu_{\theta_z}(z_{t-1}, u_t), \sigma_{\theta_z}(z_{t-1}, u_t)] = d_z(z_{t-1}, u_t), \qquad (6.1)$$
$$p_{\theta_z}(z_t|z_{t-1}, u_t) = \mathcal{N}(z_t; \mu_{\theta_z}(z_{t-1}, u_t), \text{diag}\{\sigma_{\theta_z}^2(z_{t-1}, u_t)\}), \qquad (6.2)$$

$$[\mu_{\theta_a}(z_t), \sigma_{\theta_a}(z_t)] = d_a(z_t), \qquad (6.3)$$

$$p_{\theta_a}(a_t|z_t) = \mathcal{N}(a_t; \mu_{\theta_a}(z_t), \text{diag}\{\sigma_{\theta_a}^2(z_t)\}), \qquad (6.4)$$

$$[\mu_{\theta_x}(a_t), \sigma_{\theta_x}(a_t)] = d_x(a_t), \qquad (6.5)$$

$$p_{\theta_x}(x_t|a_t) = \mathcal{N}(x_t; \mu_{\theta_x}(a_t), \text{diag}\{\sigma_{\theta_x}^2(a_t)\}). \qquad (6.6)$$

In Fraccaro et al.’s (2017) paper, (6.1) and (6.3) are linear equations; that is, they are given as (3.16) and (3.17), respectively (with **a** in place of **x**, null bias vectors **n**<sub>t</sub> and **m**<sub>t</sub>, and time-invariant covariance matrices **Λ** and **Σ**). Therefore, we have $\theta_z = \{A_t, B_t\}_{t=1}^T \cup \{\Lambda\}$ and $\theta_a = \{C_t\}_{t=1}^T \cup \{\Sigma\}$, and the submodel on $\{u_t, z_t, a_t\}$ is a classical (non-deep) LG-LDS. In contrast, $d_x$ in (6.5) is implemented with a DNN, with parameter set $\theta_x$ (e.g., a basic MLP or a CNN for video sequence modeling). This network plays the role of a deep feature extractor, with the dimension of $a_t$ being possibly much smaller than that of $x_t$. The feature vector $a_t$ is expected to encode the properties of the “object(s)”
present in the observation $\mathbf{x}_t$, whereas $\mathbf{z}_t$ is expected to encode the dynamics of these objects, which is an important application of the previously described disentanglement concept. As we will see, in the KVAE case, this can be a great advantage for solving the dynamical part of the model.

The joint distribution of all random variables can be factorized as follows:

$$p_{\theta}(\mathbf{x}_{1:T}, \mathbf{a}_{1:T}, \mathbf{z}_{1:T}, \mathbf{u}_{1:T}) = \prod_{t=1}^{T} p_{\theta_x}(\mathbf{x}_t|\mathbf{a}_t)p_{\theta_a}(\mathbf{a}_t|\mathbf{z}_t)p_{\theta_z}(\mathbf{z}_t|\mathbf{z}_{t-1}, \mathbf{u}_t)p(\mathbf{u}_t), \qquad (6.7)$$

and we also have

$$p_{\theta_x}(\mathbf{x}_{1:T}|\mathbf{a}_{1:T}) = \prod_{t=1}^{T} p_{\theta_x}(\mathbf{x}_t|\mathbf{a}_t), \qquad (6.8)$$

$$p_{\theta_a}(\mathbf{a}_{1:T}|\mathbf{z}_{1:T}) = \prod_{t=1}^{T} p_{\theta_a}(\mathbf{a}_t|\mathbf{z}_t), \qquad (6.9)$$

$$p_{\theta_z}(\mathbf{z}_{1:T}|\mathbf{u}_{1:T}) = \prod_{t=1}^{T} p_{\theta_z}(\mathbf{z}_t|\mathbf{z}_{t-1}, \mathbf{u}_t). \qquad (6.10)$$

Given the state sequence $\mathbf{z}_{1:T}$, the features $\mathbf{a}_{1:T}$ are independent, and given the sequence of features, the observations $\mathbf{x}_{1:T}$ are independent.

Fraccaro et al. (2017) mentioned the classical limitation of LDS for modeling abrupt changes in trajectories. A classical solution to this problem is to include in the model a “switching strategy” between different models or different parameterizations of the model, see, e.g., the switching Kalman filter (Murphy, 1998; Fox et al., 2011). Fraccaro et al. (2017) proposed to define each LDS parameter at each time $t$ (e.g., matrix $\mathbf{A}_t$) as a linear combination of predefined matrices/vectors from a parameter bank, and the coefficients of the linear combination were estimated at each time $t$ from the past features $\mathbf{a}_{1:t-1}$ using an LSTM network. Although it is an interesting contribution on its own, we do not further consider this part of the KVAE model here, as it is loosely relevant to our model review. A similar transition model was proposed independently by Karl et al. (2017) as an instance of a deep variational Bayesian filter (DVBF), an extended class of SSM-based DVAE models
enriched with stochastic transition parameters (see also Watter et al. (2015)).

## 6.2 Inference model

For the KVAE model, the posterior distribution of all latent variables given the observed variables, $p_{\theta}(\mathbf{a}_{1:T}, \mathbf{z}_{1:T}|\mathbf{x}_{1:T}, \mathbf{u}_{1:T})$, can be factorized as follows:

$$ 
\begin{aligned}
p_{\theta}(\mathbf{a}_{1:T}, \mathbf{z}_{1:T}|\mathbf{x}_{1:T}, \mathbf{u}_{1:T}) &= p_{\theta}(\mathbf{z}_{1:T}|\mathbf{a}_{1:T}, \mathbf{x}_{1:T}, \mathbf{u}_{1:T})p_{\theta}(\mathbf{a}_{1:T}|\mathbf{x}_{1:T}, \mathbf{u}_{1:T}), \\
&= p_{\theta}(\mathbf{z}_{1:T}|\mathbf{a}_{1:T}, \mathbf{u}_{1:T})p_{\theta}(\mathbf{a}_{1:T}|\mathbf{x}_{1:T}, \mathbf{u}_{1:T}),
\end{aligned}
\quad (6.11) $$

where the simplification of the first term on the right-hand side results from D-separation. A keypoint that appears here is that, if the sequence of features $\mathbf{a}_{1:T}$ is known, then $p_{\theta}(\mathbf{z}_{1:T}|\mathbf{a}_{1:T}, \mathbf{u}_{1:T})$ has a closed-form solution, which is a Kalman filter or a Kalman smoother (see Section 3.2.2). This Kalman solution is classic, and it is not detailed here for conciseness. We can simply mention that it depends only on $\theta_{\mathbf{a}} \cup \theta_{\mathbf{z}}$ but not on $\theta_{\mathbf{x}}$. The other factor, $p_{\theta}(\mathbf{a}_{1:T}|\mathbf{x}_{1:T}, \mathbf{u}_{1:T})$, is more problematic.

Fraccaro et al. (2017) did not discuss the form of the exact posterior distribution, yet they proposed the following factorized inference model, which exploits the Kalman solution:

$$ q_{\phi}(\mathbf{a}_{1:T}, \mathbf{z}_{1:T}|\mathbf{x}_{1:T}, \mathbf{u}_{1:T}) = p_{\theta}(\mathbf{z}_{1:T}|\mathbf{a}_{1:T}, \mathbf{u}_{1:T})q_{\phi}(\mathbf{a}_{1:T}|\mathbf{x}_{1:T}) \quad (6.12) $$

$$ = p_{\theta}(\mathbf{z}_{1:T}|\mathbf{a}_{1:T}, \mathbf{u}_{1:T}) \prod_{t=1}^{T} q_{\phi}(\mathbf{a}_{t}|\mathbf{x}_{t}), \quad (6.13) $$

where

$$ q_{\phi}(\mathbf{a}_{t}|\mathbf{x}_{t}) = \mathcal{N}(\mathbf{a}_{t}; \boldsymbol{\mu}_{\phi}(\mathbf{x}_{t}), \mathrm{diag}\{\boldsymbol{\sigma}_{\phi}^{2}(\mathbf{x}_{t})\}) \quad (6.14) $$

is implemented with a fully-connected DNN:

$$ [\boldsymbol{\mu}_{\phi}(\mathbf{x}_{t}), \boldsymbol{\sigma}_{\phi}(\mathbf{x})] = e_{\mathbf{a}}(\mathbf{x}_{t}). \quad (6.15) $$

Eqs. (6.5), (6.6), (6.14), and (6.15) are identical to (2.5), (2.3), (2.18), and (2.20), respectively, with $\mathbf{a}_t$ substituting $\mathbf{x}$ and $\mathbf{z}_t$ substituting $\mathbf{z}$. This means that, with the proposed inference model, a KVAE is
<figure><img src="image_11.png" /><figcaption>**Figure 6.2:** KVAE's graphical model at inference time. The black arrows represent the Kalman filter solution (causal solution) of LG-LDS on {$u_{1:T}, z_{1:T}, a_{1:T}$}. The blue arrows represent the VAE encoder from $x_t$ to $a_t$. The inference solution for the complete KVAE model is a combination of these two.</figcaption></figure>

composed of a VAE modeling the relationship between $\mathbf{a}_t$ and $\mathbf{x}_t$, placed on top of an LG-LDS on $\mathbf{u}_t, \mathbf{z}_t$ and $\mathbf{a}_t$. Therefore, the inference solution of the VAE and that of LG-LDS can be combined for the solution of this combined model. This is illustrated in Figure 6.2.

This inference model is thus designed to benefit from both the VAE methodology and the well-known efficiency of the “simple” LG-LDS model for tracking data dynamics. The Kalman solution, which requires inverse matrix calculation, greatly benefits from the notable $\mathbf{x}_t$-to-$\mathbf{a}_t$ dimension reduction. In addition, although this was not exactly discussed in these terms by Fraccaro et al. (2017), one possible motivation for designing the KVAE model is that the joint learning of all parameters (see the next subsection) can encourage the feature extractor to provide $\mathbf{a}_t$ features that are well-suited to a linear dynamical model; that is, the nonlinear relations between observations $\mathbf{x}_t$ and dynamics $\mathbf{z}_t$ are (at least largely) captured by the VAE.
## 6.3 Training

The VLB for the KVAE model writes as follows:

$$ \mathcal{L}(\phi, \theta; \mathbf{x}_{1:T}, \mathbf{u}_{1:T}) = \mathbb{E}_{q(\mathbf{z}_{1:T}, \mathbf{a}_{1:T}|\mathbf{x}_{1:T}, \mathbf{u}_{1:T})} \left[ \log p_{\theta}(\mathbf{x}_{1:T}|\mathbf{z}_{1:T}, \mathbf{a}_{1:T}, \mathbf{u}_{1:T}) \right] \\ - D_{KL}(q_{\phi}(\mathbf{z}_{1:T}, \mathbf{a}_{1:T}|\mathbf{x}_{1:T}, \mathbf{u}_{1:T}) \| p_{\theta}(\mathbf{z}_{1:T}, \mathbf{a}_{1:T}|\mathbf{u}_{1:T})) \quad (6.16) $$

$$ = \mathbb{E}_{q_{\phi}(\mathbf{a}_{1:T}|\mathbf{x}_{1:T})p_{\theta}(\mathbf{z}_{1:T}|\mathbf{a}_{1:T}, \mathbf{u}_{1:T})} \left[ \log p_{\theta}(\mathbf{x}_{1:T}|\mathbf{a}_{1:T}) \right] \\ - D_{KL}(q_{\phi}(\mathbf{a}_{1:T}|\mathbf{x}_{1:T})p_{\theta}(\mathbf{z}_{1:T}|\mathbf{a}_{1:T}, \mathbf{u}_{1:T}) \| p_{\theta}(\mathbf{z}_{1:T}, \mathbf{a}_{1:T}|\mathbf{u}_{1:T})) \quad (6.17) $$

$$ = \mathbb{E}_{q_{\phi}(\mathbf{a}_{1:T}|\mathbf{x}_{1:T})} \left[ \log p_{\theta}(\mathbf{x}_{1:T}|\mathbf{a}_{1:T})/q_{\phi}(\mathbf{a}_{1:T}|\mathbf{x}_{1:T}) \right. \\ \left. - D_{KL}(p_{\theta}(\mathbf{z}_{1:T}|\mathbf{a}_{1:T}, \mathbf{u}_{1:T}) \| p_{\theta}(\mathbf{a}_{1:T}|\mathbf{z}_{1:T})p_{\theta}(\mathbf{z}_{1:T}|\mathbf{u}_{1:T})) \right], \quad (6.18) $$

where we use the proposed decomposition of both the generative and inference models.

In practice, one first samples from $q_{\phi}(\mathbf{a}_t|\mathbf{x}_t)$ for each $t$. These samples are fed to a standard Kalman smoother that computes $p_{\theta}(\mathbf{z}_{1:T}|\mathbf{a}_{1:T}, \mathbf{u}_{1:T})$. We can then easily sample from this distribution. This allows for jointly learning the parameters of the VAE (both the encoder and decoder) and those of the LG-LDS.
To the best of our knowledge, STORN (Bayer and Osendorfer, 2014) is the first DVAE model to combine an internal deterministic state $h_t$ and an internal stochastic state $z_t$. Bayer and Osendorfer (2014) presented STORN in the undriven and predictive modes (i.e., with $u_t = x_{t-1}$). Hereinafter, we retain this mode for STORN, VRNN, and SRNN, for an easier comparison, but these models can also be easily set up in the driven mode with an external input $u_t$.

## 7.1 Generative model

The STORN observation model is given by

$$h_t = d_{\text{hid}}(\mathbf{W}_{\text{in}}\mathbf{x}_{t-1} + \mathbf{W}_{\text{lat}}\mathbf{z}_t + \mathbf{W}_{\text{rec}}h_{t-1} + \mathbf{b}_{\text{hid}}), \quad (7.1)$$

$$[\mu_{\theta_x}(h_t), \sigma_{\theta_x}(h_t)] = d_{\text{out}}(\mathbf{W}_{\text{out}}h_t + \mathbf{b}_{\text{out}}), \quad (7.2)$$

$$p_{\theta_x}(\mathbf{x}_t|\mathbf{h}_t) = \mathcal{N}(\mathbf{x}_t; \mu_{\theta_x}(\mathbf{h}_t), \text{diag}\{\sigma_{\theta_x}^2(\mathbf{h}_t)\}). \quad (7.3)$$

Therefore, (7.2) and (7.3) are the same as those for a basic single-layer GRNN, but in STORN, $z_t$ forms an input additional to the internal state $h_t$. Moreover, STORN was originally presented in the above framework
of a single-layer RNN, but it can be easily generalized to a deep RNN, defined by (3.5)–(3.7), by inserting **z**<sub>t</sub> as an additional input to the network (and setting **u**<sub>t</sub> = **x**<sub>t-1</sub>):

$$ \mathbf{h}_t = d_h(\mathbf{x}_{t-1}, \mathbf{z}_t, \mathbf{h}_{t-1}), \quad (7.4) $$

$$ [\boldsymbol{\mu}_{\theta_x}(\mathbf{h}_t), \boldsymbol{\sigma}_{\theta_x}(\mathbf{h}_t)] = d_x(\mathbf{h}_t), \quad (7.5) $$

$$ p_{\theta_x}(\mathbf{x}_t|\mathbf{h}_t) = \mathcal{N}(\mathbf{x}_t; \boldsymbol{\mu}_{\theta_x}(\mathbf{h}_t), \text{diag}\{\boldsymbol{\sigma}_{\theta_x}^2(\mathbf{h}_t)\}). \quad (7.6) $$

In the following, we retain this latter more general formulation for easier comparison with the other models. We denote by $\theta_h$ and $\theta_{hx}$ the set of parameters of the networks implementing $d_h$ and $d_x$, respectively, and have $\theta_x = \theta_h \cup \theta_{hx}$.

In STORN, $\mathbf{z}_t$ is assumed i.i.d. with a standard Gaussian distribution:

$$ p_{\theta_z}(\mathbf{z}_{1:T}) = \prod_{t=1}^{T} p_{\theta_z}(\mathbf{z}_t) \quad \text{with} \quad p_{\theta_z}(\mathbf{z}_t) = \mathcal{N}(\mathbf{z}_t; \mathbf{0}, \mathbf{I}_L). \quad (7.7) $$

In short, there is no temporal model on the prior distribution of $\mathbf{z}_t$ and $\theta_z = \emptyset$ (and therefore $\theta = \theta_x$). Here, it is the temporal recursion on $\mathbf{h}_t$ and the use of $\mathbf{h}_t$ to generate $\mathbf{x}_t$ that makes STORN a member of the DVAE family. The graphical model of STORN is shown in Figure 7.1 (left).

**Notation remark:** To ensure homogeneous notations across models, we slightly change the notation used by Bayer and Osendorfer (2014) by “synchronizing” $\mathbf{x}_t$ and $\mathbf{h}_t$; that is, in our presentation of STORN, $\mathbf{x}_t$ is generated from $\mathbf{h}_t$ (which is generated from $\mathbf{x}_{t-1}$, $\mathbf{h}_{t-1}$, and $\mathbf{z}_t$). In contrast, in Eq. (4) in Bayer and Osendorfer’s (2014) paper, $\mathbf{x}_{t+1}$ is generated from $\mathbf{h}_t$ (which is generated from $\mathbf{x}_t$, $\mathbf{h}_{t-1}$, and $\mathbf{z}_t$). In other words, we replace $\mathbf{x}_t$ with $\mathbf{x}_{t-1}$. This change of notation does not change the model in essence but makes the comparison with other models easier.

Eq. (7.4) shows that the two states, $\mathbf{h}_t$ and $\mathbf{z}_t$, are intricate, and the interpretation of $\mathbf{h}_t$ as a deterministic state is now an issue. In fact, $\mathbf{h}_t$ is now a random variable, but it is not a “free” one, as it is a deterministic function of the latent random variables $\mathbf{z}_t$ and $\mathbf{x}_{t-1}$ and of its previous
<figure><img src="image_12.png" /><figcaption>**Figure 7.1:** STORN’s graphical model in developed form (left) and compact form (right).</figcaption></figure>

value $\mathbf{h}_{t-1}$. We present a proper treatment of this issue in Appendix A. The recurrence on $\mathbf{h}_t$ is unfolded to consider $\mathbf{h}_t$ as a deterministic function of $\mathbf{z}_{1:t}$ and $\mathbf{x}_{1:t-1}$, which we denote $\mathbf{h}_t = \mathbf{h}_t(\mathbf{x}_{1:t-1}, \mathbf{z}_{1:t})$,<fnref n="1" /> and we consider a Dirac probability distribution over $\mathbf{h}_t$, positioned at that function’s output value. These points were briefly discussed by Bayer and Osendorfer (2014). In Appendix A, it is shown that marginalizing the joint density $p_{\theta}(\mathbf{x}_{1:T}, \mathbf{z}_{1:T}, \mathbf{h}_{1:T})$ w.r.t. $\mathbf{h}_{1:T}$ leads to

$$p_{\theta}(\mathbf{x}_{1:T}, \mathbf{z}_{1:T}) = \prod_{t=1}^{T} p_{\theta_x}(\mathbf{x}_t|\mathbf{h}_t(\mathbf{x}_{1:t-1}, \mathbf{z}_{1:t}))p(\mathbf{z}_t). \quad (7.8)$$

From the above equation and (7.7), we deduce the conditional distribution

$$p_{\theta_x}(\mathbf{x}_{1:T}|\mathbf{z}_{1:T}) = \prod_{t=1}^{T} p_{\theta_x}(\mathbf{x}_t|\mathbf{h}_t(\mathbf{x}_{1:t-1}, \mathbf{z}_{1:t})). \quad (7.9)$$

Note that these data sequence densities factorize across time frames, but the whole history of the present and past latent variables and past outputs is necessary to generate the present output. This history is summarized in $\mathbf{h}_t(\mathbf{x}_{1:t-1}, \mathbf{z}_{1:t})$.

In line with the discussion in Section 4.1.2, an alternate description of STORN can be written, where we remove the internal deterministic

<footnote n="1">¹This function also depends on the initial vectors $\mathbf{x}_0$ and $\mathbf{h}_0$ (and we can set $\mathbf{x}_0 = \emptyset$), but we omit them for clarity of presentation.</footnote>
state **h** and express the model only in terms of the free random variables **x**<sub>1:T</sub> and **z**<sub>1:T</sub>:

$$p_{\theta}(\mathbf{x}_{1:T}, \mathbf{z}_{1:T}) = \prod_{t=1}^{T} p_{\theta_x}(\mathbf{x}_t | \mathbf{x}_{1:t-1}, \mathbf{z}_{1:t}) p(\mathbf{z}_t), \quad (7.10)$$

and

$$p_{\theta_x}(\mathbf{x}_{1:T} | \mathbf{z}_{1:T}) = \prod_{t=1}^{T} p_{\theta_x}(\mathbf{x}_t | \mathbf{x}_{1:t-1}, \mathbf{z}_{1:t}). \quad (7.11)$$

The above two equations are more general than (7.8) and (7.9), but they lose some information on the deterministic link between **x**<sub>1:t-1</sub> and **z**<sub>1:t</sub> in the process of generating **x**<sub>t</sub>. The compact graphical representation corresponding to this alternate formulation is given in Figure 7.1 (right).

## 7.2 Inference model

Following Section 4.2, it is easy to show that the exact posterior distribution of STORN takes the following form:

$$p_{\theta}(\mathbf{z}_{1:T} | \mathbf{x}_{1:T}) = \prod_{t=1}^{T} p_{\theta}(\mathbf{z}_t | \mathbf{z}_{1:t-1}, \mathbf{x}_{1:T}). \quad (7.12)$$

In fact, this expression is obtained by the chain rule, and it cannot be simplified by applying D-separation. This is because any vector in **z**<sub>1:t-1</sub> and **x**<sub>1:T</sub> is either a child, a parent, or a co-parent of **z**<sub>t</sub> in the graphical representation of STORN. In other words, each product term at time t in (7.12) depends on the past observed and latent state vectors that propagate through the internal state, and it depends on the present and future observed vectors, as **z**<sub>t</sub> propagates to them through the internal hidden states.

To complement the discussion on the form of the exact posterior distribution, we note the following:

$$p_{\theta}(\mathbf{z}_{1:T} | \mathbf{x}_{1:T}) \propto p_{\theta}(\mathbf{x}_{1:T}, \mathbf{z}_{1:T}) = p_{\theta_x}(\mathbf{x}_{1:T} | \mathbf{z}_{1:T}) p(\mathbf{z}_{1:T}). \quad (7.13)$$

Thus, combining the above equation with (7.8), we get

$$p_{\theta}(\mathbf{z}_{1:T} | \mathbf{x}_{1:T}) \propto \prod_{t=1}^{T} p_{\theta_x}(\mathbf{x}_t | \mathbf{h}_t(\mathbf{x}_{1:t-1}, \mathbf{z}_{1:t})) p(\mathbf{z}_t). \quad (7.14)$$
As $\mathbf{z}_t$ is present in all terms of $\mathbf{h}_t(\mathbf{x}_{1:t-1}, \mathbf{z}_{1:t})$ from $t$ to $T$, it is confirmed that (7.14), considered as a function of $\mathbf{z}_t$, depends not only on $\mathbf{x}_{1:t-1}$ but also on $\mathbf{x}_{t:T}$.

As for the practical inference in STORN, the approximate posterior distribution $q_{\phi}$ was chosen by Bayer and Osendorfer (2014) as generated by an additional forward RNN. Little information is available on the implementation. The parameters of $q_{\phi}$ are said to be generated from $\mathbf{x}_{1:t}$, which is slightly odd as $\mathbf{x}_{t+1}$ was generated from $\mathbf{z}_t$ (with their notations; see our remark in the previous subsection). There is thus a one-step lag between generation and inference, which is difficult to justify (in practice, we have found that this leads to significantly inferior inference performance). With our change in notation at generation, we somehow automatically compensate for this gap and assume that the “correct” detailed inference equations are given as

$$ \mathbf{g}_t = e_g(\mathbf{W}_{\text{in}}^{\text{enc}}\mathbf{x}_t + \mathbf{W}_{\text{rec}}^{\text{enc}}\mathbf{g}_{t-1} + \mathbf{b}_{\text{hid}}^{\text{enc}}), \quad (7.15) $$

$$ [\mu_{\phi}(\mathbf{g}_t), \sigma_{\phi}(\mathbf{g}_t)] = e_z(\mathbf{W}_{\text{out}}^{\text{enc}}\mathbf{g}_t + \mathbf{b}_{\text{out}}^{\text{enc}}), \quad (7.16) $$

$$ q_{\phi}(\mathbf{z}_t|\mathbf{g}_t) = \mathcal{N}(\mathbf{z}_t; \mu_{\phi}(\mathbf{g}_t), \text{diag}\{\sigma_{\phi}^2(\mathbf{g}_t)\}), \quad (7.17) $$

where $\mathbf{g}_t$ denotes the inference RNN internal state,<fnref n="2" /> and $e_g$ and $e_z$ are nonlinear activation functions. Similarly to the generative model, because of the recursivity in (7.15), $\mathbf{g}_t$ can be considered as an unfolded deterministic function of $\mathbf{x}_{1:t}$, which we note $\mathbf{g}_t = \mathbf{g}_t(\mathbf{x}_{1:t}).$<fnref n="3" /> For a complete data sequence, we have

$$ q_{\phi}(\mathbf{z}_{1:T}|\mathbf{x}_{1:T}) = \prod_{t=1}^{T} q_{\phi}(\mathbf{z}_t|\mathbf{g}_t) = \prod_{t=1}^{T} q_{\phi}(\mathbf{z}_t|\mathbf{g}_t(\mathbf{x}_{1:t})). \quad (7.18) $$

This inference model can be rewritten in compact form as

$$ q_{\phi}(\mathbf{z}_{1:T}|\mathbf{x}_{1:T}) = \prod_{t=1}^{T} q_{\phi}(\mathbf{z}_t|\mathbf{x}_{1:t}). \quad (7.19) $$

The corresponding graphical model is shown in Figure 7.2. Note that this inference model is inconsistent with the exact posterior distribution $p_{\theta}(\mathbf{z}_{1:T}|\mathbf{x}_{1:T})$ at several points: At each time $t$, it neither considers

<footnote n="2">²In Bayer and Osendorfer’s (2014) paper, it is denoted as $\mathbf{h}_t^r$.</footnote>

<footnote n="3">³This function is also a function of $\mathbf{g}_0$.</footnote>
<figure><img src="image_13.png" /></figure>

<figure><img src="image_14.png" /><figcaption>**Figure 7.2:** STORN's graphical model at inference time in developed form (left) and compact form (right). Golden arrows correspond to missing links on the proposed probabilistic dependencies (compared to the exact inference dependencies).</figcaption></figure>

future observations $\mathbf{x}_{t+1:T}$ nor past latent states $\mathbf{z}_{1:t-1}$. As discussed in Section 4.2.3, the internal states of the encoder and of decoder can be identical, or they can be different. In STORN, given the choice of the inference model, these internal states depend on different variables and, therefore, are different.

## 7.3 Training

A comparison of the compact form of STORN in (7.10) with the general compact form of a DVAE in (4.4) (simplified without $\mathbf{u}_{1:T}$) shows that STORN makes the following conditional independence assumption:

$$p_{\theta_z}(\mathbf{z}_t | \mathbf{x}_{1:t-1}, \mathbf{z}_{1:t-1}) = p(\mathbf{z}_t). \qquad (7.20)$$

Using this single simplification, the VLB given in its most general form in (4.24) becomes

$$\begin{aligned} \mathcal{L}(\theta, \phi; \mathbf{x}_{1:T}) = & \sum_{t=1}^{T} \mathbb{E}_{q_{\phi}(\mathbf{z}_{1:t} | \mathbf{x}_{1:T})} [\log p_{\theta_x}(\mathbf{x}_t | \mathbf{x}_{1:t-1}, \mathbf{z}_{1:t})] \\ & - \sum_{t=1}^{T} \mathbb{E}_{q_{\phi}(\mathbf{z}_{1:t-1} | \mathbf{x}_{1:T})} [D_{\text{KL}}(q_{\phi}(\mathbf{z}_t | \mathbf{z}_{1:t-1}, \mathbf{x}_{1:T}) \| p(\mathbf{z}_t))] . \end{aligned} \qquad (7.21)$$

This expression of the VLB relies on an inference model that is consistent with the exact posterior distribution (7.12). However, as discussed above,
STORN assumes an inference model of the form: $q_{\phi}(\mathbf{z}_t|\mathbf{z}_{1:t-1}, \mathbf{x}_{1:T}) = q_{\phi}(\mathbf{z}_t|\mathbf{x}_{1:t})$. Consequently, the VLB in (7.21) can be simplified as follows:

$$ \mathcal{L}(\theta, \phi; \mathbf{x}_{1:T}) = \sum_{t=1}^{T} \mathbb{E}_{q_{\phi}(\mathbf{z}_{1:t}|\mathbf{x}_{1:T})} \left[ \log p_{\theta_{\star}}(\mathbf{x}_{t}|\mathbf{x}_{1:t-1}, \mathbf{z}_{1:t}) \right] \\ - \sum_{t=1}^{T} D_{\text{KL}}(q_{\phi}(\mathbf{z}_t|\mathbf{x}_{1:t}) \parallel p(\mathbf{z}_t)). \qquad (7.22) $$

The KL divergence in this expression can be computed analytically, while the expectation is intractable and should be approximated by a Monte Carlo estimate, using samples drawn recursively from $q_{\phi}(\mathbf{z}_{1:t}|\mathbf{x}_{1:T})$ based on the inference model (7.18). As for DKF, using the reparameterization trick for this recursive sampling leads to an objective function, which is differentiable w.r.t. $\theta$ and $\phi$.
<figure><img src="image_15.png" /></figure>
The VRNN model was proposed by Chung et al. (2015) as a combination of a VAE and an RNN.

## 8.1 Generative model

The VRNN observation model is given by

$$ \mathbf{h}_t = d_{\mathbf{h}}(\varphi_{\mathbf{x}}(\mathbf{x}_{t-1}), \varphi_{\mathbf{z}}(\mathbf{z}_{t-1}), \mathbf{h}_{t-1}), \quad (8.1) $$

$$ [\mu_{\theta_x}(\mathbf{z}_t, \mathbf{h}_t), \sigma_{\theta_x}(\mathbf{z}_t, \mathbf{h}_t)] = d_x(\varphi_z(\mathbf{z}_t), \mathbf{h}_t), \quad (8.2) $$

$$ p_{\theta_x}(\mathbf{x}_t|\mathbf{z}_t, \mathbf{h}_t) = \mathcal{N}(\mathbf{x}_t; \mu_{\theta_x}(\mathbf{z}_t, \mathbf{h}_t), \text{diag}\{\sigma_{\theta_x}^2(\mathbf{z}_t, \mathbf{h}_t)\}), \quad (8.3) $$

where $\varphi_z$ and $\varphi_x$ are feature extractors, which were mentioned by Chung et al. (2015) to be important in practice. These feature extractors are DNNs parameterized by a set of weights and biases, denoted as $\tau$.

The generative distribution of $\mathbf{z}_t$ is given by

$$ [\mu_{\theta_z}(\mathbf{h}_t), \sigma_{\theta_z}(\mathbf{h}_t)] = d_z(\mathbf{h}_t), \quad (8.4) $$

$$ p_{\theta_z}(\mathbf{z}_t|\mathbf{h}_t) = \mathcal{N}(\mathbf{z}_t; \mu_{\theta_z}(\mathbf{h}_t), \text{diag}\{\sigma_{\theta_z}^2(\mathbf{h}_t)\}). \quad (8.5) $$
<figure><img src="image_16.png" /><figcaption>**Figure 8.1:** VRNN's graphical model in developed (left) and compact (right) forms.</figcaption></figure>

**Notation remark:** In Chung et al.'s (2015) paper, $\tau$ denotes the set of parameters for both feature extractors, which are denoted as $\varphi_{\tau}^{x}$ and $\varphi_{\tau}^{z}$, respectively, as well as for $d_{x}$ and $d_{z}$, which are denoted $\varphi_{\tau}^{\text{dec}}$ and $\varphi_{\tau}^{\text{prior}}$, respectively. Moreover, $d_{h}$ is denoted $d_{\theta}$. We find this a bit confusing and prefer to distinguish among $\tau$, $\theta_{x} = \theta_{h} \cup \theta_{hx}$, $\theta_{z} = \theta_{h} \cup \theta_{hz}$, and $\theta = \tau \cup \theta_{x} \cup \theta_{z}$. Moreover, one may also want to distinguish between $\tau_{x}$ and $\tau_{z}$ to clarify that the two feature extractors are different. We retain $\tau$ for simplicity. Besides, we replace $h_{t-1}$ in the paper by Chung et al. (2015) with $h_{t}$ to synchronize $h_{t}$, $z_{t}$, and $x_{t}$; that is, $x_{t}$ is generated from $h_{t}$ and $z_{t}$. This arbitrary reindexing of $h_{t}$ does not change the model conceptually.

In VRNN, we thus have multiple intrications of $h_t$ and $z_t$ in both the observation model and the distribution of $z_t$. The graphical model of VRNN is given in Figure 8.1 (left). The generative process starts with an initial internal state $h_1$, from which we generate $z_1$. From $h_1$ and $z_1$, we generate $x_1$. Then, $h_2$ is deterministically calculated from $h_1$, $z_1$, and $x_1$, and so on, except that hereinafter, $z_t$ is generated from $z_{t-1}$, $h_{t-1}$, and $x_{t-1}$. Using the “unfolding the recurrence” trick mentioned in the previous sections, we can here denote $h_t = h_t(x_{1:t-1}, z_{1:t-1})^1$ and have $p_{\theta_z}(z_t|h_t) = p_{\theta_z}(z_t|h_t(x_{1:t-1}, z_{1:t-1}))$ and $p_{\theta_x}(x_t|z_t, h_t) = p_{\theta_x}(x_t|z_t, h_t(x_{1:t-1}, z_{1:t-1}))$. This provides both $z_t$ and $x_t$ with an implicit temporal model.

<footnote n="1">¹This function also depends on $h_1$, which is omitted for clarity of presentation.</footnote>
As for a data sequence, when marginalizing the joint distribution of all variables w.r.t. $\mathbf{h}_{1:T}$ following the line of Appendix A, we get

$$p_{\theta}(\mathbf{x}_{1:T}, \mathbf{z}_{1:T}) = \prod_{t=1}^{T} p_{\theta_x}(\mathbf{x}_t|\mathbf{z}_t, \mathbf{h}_t(\mathbf{x}_{1:t-1}, \mathbf{z}_{1:t-1}))p_{\theta_z}(\mathbf{z}_t|\mathbf{h}_t(\mathbf{x}_{1:t-1}, \mathbf{z}_{1:t-1})) \quad (8.6)$$

$$= \prod_{t=1}^{T} p_{\theta}(\mathbf{x}_t, \mathbf{z}_t|\mathbf{h}_t(\mathbf{x}_{1:t-1}, \mathbf{z}_{1:t-1})). \quad (8.7)$$

Again, we have a factorization of the conditional densities over time frames. However, we do not have conditional independence of $\mathbf{x}_t$ and $\mathbf{z}_t$ conditionally to the state $\mathbf{h}_t(\mathbf{x}_{1:t-1}, \mathbf{z}_{1:t-1})$ due to the direct link from $\mathbf{z}_t$ to $\mathbf{x}_t$.

As for STORN, we can provide a more general alternate expression for $p_{\theta}(\mathbf{x}_{1:T}, \mathbf{z}_{1:T})$ that does not make the internal state explicit but only represents the general dependencies among the “free” random variables:

$$p_{\theta}(\mathbf{x}_{1:T}, \mathbf{z}_{1:T}) = \prod_{t=1}^{T} p_{\theta_x}(\mathbf{x}_t|\mathbf{x}_{1:t-1}, \mathbf{z}_{1:t})p_{\theta_z}(\mathbf{z}_t|\mathbf{x}_{1:t-1}, \mathbf{z}_{1:t-1}) \quad (8.8)$$

$$= \prod_{t=1}^{T} p_{\theta}(\mathbf{x}_t, \mathbf{z}_t|\mathbf{x}_{1:t-1}, \mathbf{z}_{1:t-1}). \quad (8.9)$$

The corresponding compact graphical model is shown in Figure 8.1 (right).

## 8.2 Inference model

The general form of the exact posterior distribution of VRNN is identical to the one of STORN; that is, it factorizes into

$$p_{\theta}(\mathbf{z}_{1:T}|\mathbf{x}_{1:T}) = \prod_{t=1}^{T} p_{\theta}(\mathbf{z}_t|\mathbf{z}_{1:t-1}, \mathbf{x}_{1:T}), \quad (8.10)$$

and, here also, no further simplification can be obtained from D-separation.

The approximate posterior distribution $q_{\phi}$ was chosen by Chung
et al. (2015) as

$$[\mu_{\phi}(\mathbf{x}_{t}, \mathbf{h}_{t}), \sigma_{\phi}(\mathbf{x}_{t}, \mathbf{h}_{t})] = e_{z}(\varphi_{x}(\mathbf{x}_{t}), \mathbf{h}_{t}), \qquad (8.11)$$

$$q_{\phi}(\mathbf{z}_{t}|\mathbf{x}_{t}, \mathbf{h}_{t}) = \mathcal{N}(\mathbf{z}_{t}; \mu_{\phi}(\mathbf{x}_{t}, \mathbf{h}_{t}), \text{diag}\{\boldsymbol{\sigma}_{\phi}^{2}(\mathbf{x}_{t}, \mathbf{h}_{t})\}), \quad (8.12)$$

where $e_z$ is the encoder DNN, parameterized by $\phi_z$. As for data sequence inference, we have

$$q_{\phi}(\mathbf{z}_{1:T}|\mathbf{x}_{1:T}) = \prod_{t=1}^{T} q_{\phi}(\mathbf{z}_{t}|\mathbf{x}_{t}, \mathbf{h}_{t}(\mathbf{x}_{1:t-1}, \mathbf{z}_{1:t-1})). \qquad (8.13)$$

In contrast to STORN, the same internal state $\mathbf{h}_t$ is here shared by the VRNN encoder and decoder, which, in our opinion, makes the approximate model more consistent with the exact posterior distribution. This makes the set of parameters $\theta_{\mathbf{h}}$ common to the encoder and decoder. Because the feature extractor $\varphi_{\mathbf{x}}(\mathbf{x}_t)$ is also used at the encoder, the same remark applies to its parameter set $\tau_{\mathbf{x}}$. In addition, the inference at time $t$ depends on past outputs and past latent states, which also makes the inference model closer to the exact posterior distribution. However, compared to the exact posterior, the future observations (from $t+1$ to $T$) are missing again. In short, here also, the approximate inference is causal, whereas the exact posterior distribution is noncausal. The graphical model corresponding to the VRNN approximate inference process is shown in Figure 8.2. The inference model can be rewritten in the following general form:

$$q_{\phi}(\mathbf{z}_{1:T}|\mathbf{x}_{1:T}) = \prod_{t=1}^{T} q_{\phi}(\mathbf{z}_{t}|\mathbf{z}_{1:t-1}, \mathbf{x}_{1:t}). \qquad (8.14)$$

## 8.3 Training

A comparison of the compact form of VRNN in (8.8) with the general compact form of a DVAE in (4.4) (simplified without $\mathbf{u}_{1:T}$) shows that VRNN does not make any conditional independence assumption in the generative model. In this sense, VRNN is the most general DVAE model we have seen so far. The expression of the VLB for VRNN should therefore be the one given in (4.24). However, as discussed above,
<figure><img src="image_17.png" /><figcaption>**Figure 8.2:** VRNN's graphical model at inference time in developed form (left) and compact form (right). Golden arrows correspond to missing links on the proposed probabilistic dependencies (compared to the exact inference dependencies).</figcaption></figure>

the inference model in VRNN is inconsistent with the exact posterior distribution, as the following conditional independence assumption is made: $q_{\phi}(\mathbf{z}_t|\mathbf{z}_{1:t-1}, \mathbf{x}_{1:T}) = q_{\phi}(\mathbf{z}_t|\mathbf{z}_{1:t-1}, \mathbf{x}_{1:t})$. Consequently, the VLB in (4.24) can be simplified as follows:

$$ \mathcal{L}(\theta, \phi; \mathbf{x}_{1:T}) = \sum_{t=1}^{T} \mathbb{E}_{q_{\phi}(\mathbf{z}_{1:t}|\mathbf{x}_{1:T})} [\log p_{\theta_x}(\mathbf{x}_t|\mathbf{x}_{1:t-1}, \mathbf{z}_{1:t})] \\ - \sum_{t=1}^{T} \mathbb{E}_{q_{\phi}(\mathbf{z}_{1:t-1}|\mathbf{x}_{1:T})} [D_{\text{KL}}(q_{\phi}(\mathbf{z}_t|\mathbf{z}_{1:t-1}, \mathbf{x}_{1:t}) \parallel p_{\theta_z}(\mathbf{z}_t|\mathbf{x}_{1:t-1}, \mathbf{z}_{1:t-1}))]. \quad (8.15) $$

As for the previously presented DVAE models, the KL divergence can be computed analytically, and intractable expectations are approximated by Monte Carlo estimates.

## 8.4 Improved VRNN and VRNN applications

To complement this VRNN section, we report that an improved version of VRNN was presented by Goyal et al. (2017). The authors pointed out the difficulty in learning meaningful latent variables when coupled with a strong autoregressive decoder. We further discuss this point and we provide a series of references in Chapter 14. Goyal et al. (2017) proposed to improve the inference and training of VRNN with the following three features.
First, possibly inspired by Krishnan et al. (2017), as well as (Fraccaro et al., 2016) (see Section 9.2), they introduced a backward RNN on **x**_t to feed the approximate posterior distribution q_φ(**z**<sub>1:T</sub>|**x**<sub>1:T</sub>), in line with DKS (see Section 5.2). Therefore, they accounted for the future observations in the inference process, as opposed to the original VRNN, moving toward a better compliance with the structure of the exact posterior distribution p_θ(**z**<sub>1:T</sub>|**x**<sub>1:T</sub>).

Second, they forced the latent variable **z**_t to contain relevant in-
formation about the future of the sequence by connecting **z**_t with the
internal state of the inference backward network (denoted **b**_t by Goyal
et al. (2017)). This was achieved by introducing an additional condi-
tional model p_ε(**b**_t|**z**_t) and adding it in the VLB. Similarly, they also
considered an additional conditional model p_ε(**x**_t|**b**_t).

Finally, they slightly modified the VRNN model itself by removing
the direct link between **z**_t and **x**_t; that is, they replaced (8.2)–(8.3)
with (7.5)–(7.6), while all other equations remained identical to the
VRNN equations. The authors reported that “[they] observed better
performance by avoiding the latent variables from directly producing
the next output.”

An adaptation of VRNN to automatic language translation was pro-
posed by Su et al. (2018). This is doubly interesting because this paper
considers a sequence of discrete inputs and outputs, which contrasts
with the “all continuous” models we focus on. This paper also contrasts
with the previously proposed VAE-based models for text/language pro-
cessing, which, as mentioned in the Introduction, usually consider a
single latent vector to encode the whole input sequence (a full sentence).
In Su et al.’s (2018) paper, it is the sequence of latent vectors **z**<sub>1:T</sub> that
encodes the semantic content of the sequence to translate “over time.”

Finally, we can also mention the study by Lee et al. (2018), which
uses VRNN for speech synthesis and adopts adversarial training. Besides
VRNN, all these papers illustrate the flexibility of the DVAE class of
models.
The SRNN model was proposed by Fraccaro et al. (2016), with an objective to “glue (or stack) a deterministic recurrent neural network and a state space model together to form a stochastic and sequential neural generative model.”

**Notation remark:** In Fraccaro et al.’s (2016) paper, $\mathbf{h}_t$ is denoted as $\mathbf{d}_t$, and the model is presented in the driven mode. We replace the external input $\mathbf{u}_t$ with $\mathbf{x}_{t-1}$ (i.e., predictive mode) for a better comparison with VRNN and STORN.

## 9.1 Generative model

The SRNN observation model is given by

$$ \mathbf{h}_t = d_{\mathbf{h}}(\mathbf{x}_{t-1}, \mathbf{h}_{t-1}), \qquad (9.1) $$

$$ [\mu_{\theta_x}(\mathbf{z}_t, \mathbf{h}_t), \sigma_{\theta_x}(\mathbf{z}_t, \mathbf{h}_t)] = d_x(\mathbf{z}_t, \mathbf{h}_t), \qquad (9.2) $$

$$ p_{\theta_x}(\mathbf{x}_t|\mathbf{z}_t, \mathbf{h}_t) = \mathcal{N}(\mathbf{x}_t; \mu_{\theta_x}(\mathbf{z}_t, \mathbf{h}_t), \text{diag}\{\sigma_{\theta_x}^2(\mathbf{z}_t, \mathbf{h}_t)\}). \qquad (9.3) $$

Eq. (9.1) is identical to (3.3) (with $\mathbf{x}_{t-1}$ instead of $\mathbf{u}_t$) and thus refers to the usual deterministic RNN. Eqs. (9.2) and (9.3) are quite similar
<figure><img src="image_18.png" /><figcaption>Figure 9.1: SRNN’s graphical model in developed (left) and compact (right) forms.</figcaption></figure>

to (8.2) and (8.3), respectively. Thus, the internal state $h_t$ remains deterministic here, and the latent state $z_t$ is integrated at the $d_x$ level. This justifies the “clear(er) separation of deterministic and stochastic layers” claimed by Fraccaro et al. (2016), compared to VRNN.

Fraccaro et al. (2016) also introduced an explicit temporal model on the distribution of $z_t$ (as opposed to implicit temporal dependency through $h_t$ in VRNN), in addition to the dependency on the internal state $h_t$:

$$[\mu_{\theta_z}(z_{t-1}, h_t), \sigma_{\theta_z}(z_{t-1}, h_t)] = d_z(z_{t-1}, h_t), \quad (9.4)$$

$$p_{\theta_z}(z_t|z_{t-1}, h_t) = \mathcal{N}(z_t; \mu_{\theta_z}(z_{t-1}, h_t), \text{diag}\{\sigma_{\theta_z}^2(z_{t-1}, h_t)\}). \quad (9.5)$$

Compared to VRNN, the arrow from $z_{t-1}$ to $h_t$ is replaced with an arrow from $z_{t-1}$ to $z_t$, leading to an explicit first-order Markovian dependency for $z_t$ (which is combined with the dependency on $h_t$). In addition, compared to VRNN, no feature extractor is mentioned in SRNN, so we have here $\theta = \theta_x \cup \theta_z$ (and again $\theta_x = \theta_h \cup \theta_{hx}$ and $\theta_z = \theta_h \cup \theta_{hz}$). The graphical model of SRNN is shown in Figure 9.1 (left). In Fraccaro et al.’s (2016) paper, both $d_x$ and $d_z$ are two-layer feed-forward networks. The function $d_h$ is a GRU RNN, so that, according to the authors, “the SSM can therefore utilize long-term information captured by the RNN.”

Using the same “unfolding the recurrence” trick as in the previous sections, we here denote $h_t$ as $h_t(x_{1:t-1})^1$ and have $p_{\theta_z}(z_t|z_{t-1}, h_t) =$

<footnote n="1">¹Again, we omit the initial term $h_1$ for clarity of presentation.</footnote>
$p_{\theta_z}(\mathbf{z}_t|\mathbf{z}_{t-1}, \mathbf{h}_t(\mathbf{x}_{1:t-1}))$ and $p_{\theta_x}(\mathbf{x}_t|\mathbf{z}_t, \mathbf{h}_t) = p_{\theta_x}(\mathbf{x}_t|\mathbf{z}_t, \mathbf{h}_t(\mathbf{x}_{1:t-1}))$. Again, if we follow the line of Appendix A, marginalizing the joint distribution of all variables w.r.t. $\mathbf{h}_{1:T}$ leads to

$$p_{\theta}(\mathbf{x}_{1:T}, \mathbf{z}_{1:T}) = \prod_{t=1}^{T} p_{\theta_x}(\mathbf{x}_t|\mathbf{z}_t, \mathbf{h}_t(\mathbf{x}_{1:t-1}))p_{\theta_z}(\mathbf{z}_t|\mathbf{z}_{t-1}, \mathbf{h}_t(\mathbf{x}_{1:t-1})). \quad (9.6)$$

As for VRNN, we have a factorization over time frames, but no in-dependence of $\mathbf{x}_t$ and $\mathbf{z}_t$ conditionally to the state $\mathbf{h}_t(\mathbf{x}_{1:t-1})$. As for STORN and VRNN, (9.6) can be reshaped into the following more general expression:

$$p_{\theta}(\mathbf{x}_{1:T}, \mathbf{z}_{1:T}) = \prod_{t=1}^{T} p_{\theta_x}(\mathbf{x}_t|\mathbf{x}_{1:t-1}, \mathbf{z}_t)p_{\theta_z}(\mathbf{z}_t|\mathbf{z}_{t-1}, \mathbf{x}_{1:t-1}). \quad (9.7)$$

The corresponding compact graphical model is shown in Figure 9.1 (right).

## 9.2 Inference model

For SRNN, because of the dependencies in the generative model, the general form of the exact posterior distribution is given by

$$p_{\theta}(\mathbf{z}_{1:T}|\mathbf{x}_{1:T}) = \prod_{t=1}^{T} p_{\theta}(\mathbf{z}_t|\mathbf{z}_{t-1}, \mathbf{x}_{1:T}). \quad (9.8)$$

At each time $t$, the posterior distribution of $\mathbf{z}_t$ depends on the previous latent state and the whole observation sequence.

Fraccaro et al. (2016) indicated this structure and proposed an approximate posterior distribution $q_{\phi}$ with the same structure. This is the second time this proper methodology is considered in the present review after the DKS in Section 5.2, but in the publication chronology, to the best of our knowledge, this was the first time. The dependency of $q_{\phi}$ on future observations, as well as on past observations through the future internal states, is implemented with a gated backward RNN. This network is denoted by $e_{\leftarrow\mathbf{g}}$ in the equations below, it is parameterized by $\phi_{\leftarrow\mathbf{g}}$, and $\leftarrow\mathbf{g}_t$ denotes its internal state, where the right-to-left arrow highlights the backward nature of the process. This backward RNN is
followed by a basic feed-forward network $e_z$, which is parameterized by $\phi_z$. Formally, $q_\phi$ can be written as

$$q_{\phi}(\mathbf{z}_{1:T}|\mathbf{x}_{1:T}) = \prod_{t=1}^{T} q_{\phi}(\mathbf{z}_t|\mathbf{z}_{t-1}, \overleftarrow{\mathbf{g}}_t), \qquad (9.9)$$

with

$$\overleftarrow{\mathbf{g}}_t = e_{\overleftarrow{\mathbf{g}}} ([\mathbf{h}_t, \mathbf{x}_t], \overleftarrow{\mathbf{g}}_{t+1}), \qquad (9.10)$$

$$[\boldsymbol{\mu}_{\phi}(\mathbf{z}_{t-1}, \overleftarrow{\mathbf{g}}_t), \boldsymbol{\sigma}_{\phi}(\mathbf{z}_{t-1}, \overleftarrow{\mathbf{g}}_t)] = e_z(\mathbf{z}_{t-1}, \overleftarrow{\mathbf{g}}_t), \qquad (9.11)$$

$$q_{\phi}(\mathbf{z}_t|\mathbf{z}_{t-1}, \overleftarrow{\mathbf{g}}_t) = \mathcal{N}(\mathbf{z}_t; \boldsymbol{\mu}_{\phi}(\mathbf{z}_{t-1}, \overleftarrow{\mathbf{g}}_t), \mathrm{diag}(\boldsymbol{\sigma}_{\phi}^2(\mathbf{z}_{t-1}, \overleftarrow{\mathbf{g}}_t))). \qquad (9.12)$$

We thus have here $\phi = \phi_{\overleftarrow{\mathbf{g}}} \cup \phi_z$.

**Notation remark:** In Fraccaro et al.'s (2016) paper, $\mathbf{h}_t$ is denoted by $\mathbf{d}_t$, $\overleftarrow{\mathbf{g}}_t$ is denoted by $\mathbf{a}_t$, and $q_\phi(\mathbf{z}_{1:T}|\mathbf{x}_{1:T})$ is denoted by $q_\phi(\mathbf{z}_{1:T}|\mathbf{d}_{1:T}, \mathbf{x}_{1:T})$. Because we have $\mathbf{h}_t = \mathbf{h}_t(\mathbf{x}_{1:t-1})$, we can stick to $q_\phi(\mathbf{z}_{1:T}|\mathbf{x}_{1:T})$.

The above equations show that inference requires a forward pass on the internal state $\mathbf{h}_t$ (which is shared by the encoder and decoder), its combination with $\mathbf{x}_t$, and a backward pass on the inference RNN, which makes $\overleftarrow{\mathbf{g}}_t$ a deterministic function of the whole data sequence $\mathbf{x}_{1:T}$. Similarly to $\mathbf{h}_t$, we can denote this function by $\overleftarrow{\mathbf{g}}_t(\mathbf{x}_{1:T})$ to make this latter point explicit; however, this would be poorly informative about the way $\overleftarrow{\mathbf{g}}_t$ depends on $\mathbf{x}_{1:T}$. The graphical model corresponding to the inference process in SRNN is shown in Figure 9.2. The inference model can be rewritten in the following general form:

$$q_{\phi}(\mathbf{z}_{1:T}|\mathbf{x}_{1:T}) = \prod_{t=1}^{T} q_{\phi}(\mathbf{z}_t|\mathbf{z}_{t-1}, \mathbf{x}_{1:T}). \qquad (9.13)$$

Fraccaro et al. (2016) stated that this smoothing process (combination of forward and backward RNNs on $\mathbf{x}_t$) can be replaced with a filtering process, by replacing (9.10)–(9.11) with an “instantaneous” DNN $e_z(\mathbf{z}_{t-1}, \mathbf{h}_t, \mathbf{x}_t)$.
<figure><img src="image_19.png" /><figcaption>**Figure 9.2:** SRNN's graphical model at inference time in developed form (left) and compact form (right). In this case, there are no missing links on the proposed probabilistic dependencies (compared to the exact inference dependencies). In Fraccaro et al.'s (2016) paper, the dependencies of $h_t$ were omitted in the inference graphical model for clarity. We make them explicit here to recall that $h_t$ follows the deterministic update (9.1).</figcaption></figure>

## 9.3 Training

A comparison of the compact form of SRNN in (9.7) with the general compact form of a DVAE in (4.4) (simplified without $u_{1:T}$) shows that the SRNN model makes the following conditional independence assumptions:

$$ 
\begin{aligned}
p_{\theta_x}(\mathbf{x}_t|\mathbf{x}_{1:t-1}, \mathbf{z}_{1:t}) &= p_{\theta_x}(\mathbf{x}_t|\mathbf{x}_{1:t-1}, \mathbf{z}_t); \\
p_{\theta_z}(\mathbf{z}_t|\mathbf{x}_{1:t-1}, \mathbf{z}_{1:t-1}) &= p_{\theta_z}(\mathbf{z}_t|\mathbf{x}_{1:t-1}, \mathbf{z}_{t-1}).
\end{aligned} 
\quad (9.14) $$

Using these two simplifications, along with the inference model (9.13) (which we recall is consistent with the exact posterior distribution), the
VLB in its most general form (4.24) can be simplified as follows:

$$ \mathcal{L}(\theta, \phi; \mathbf{x}_{1:T}) = \sum_{t=1}^{T} \mathbb{E}_{q_{\phi}(\mathbf{z}_t|\mathbf{x}_{1:T})} \left[ \log p_{\theta_x}(\mathbf{x}_t|\mathbf{x}_{1:t-1}, \mathbf{z}_t) \right] \\ - \sum_{t=1}^{T} \mathbb{E}_{q_{\phi}(\mathbf{z}_{t-1}|\mathbf{x}_{1:T})} \left[ D_{\text{KL}}(q_{\phi}(\mathbf{z}_t|\mathbf{z}_{t-1}, \mathbf{x}_{1:T}) \parallel p_{\theta_z}(\mathbf{z}_t|\mathbf{x}_{1:t-1}, \mathbf{z}_{t-1})) \right]. \tag{9.15} $$

Again, the KL divergence can be computed analytically, and intractable expectations are approximated by Monte Carlo estimates. The procedure to sample from $q_{\phi}(\mathbf{z}_t|\mathbf{x}_{1:T})$ and $q_{\phi}(\mathbf{z}_{t-1}|\mathbf{x}_{1:T})$ relies on the “cascade trick,” as for DKF (see Section 5.3).
10

Recurrent Variational Autoencoders

The RVAE model was introduced by Leglaive et al. (2020) to repre-
sent clean speech signals in a speech enhancement application. It was
combined with a Gaussian noise model with nonnegative Matrix fac-
torization of the variance within a Bayesian framework. The RVAE
parameters were estimated offline on a large dataset of clean speech
signals using the VAE methodology (maximization of the VLB). A
variational expectation-maximization (VEM) algorithm was used for
estimating the remaining parameters from a noisy speech signal, and
probabilistic Wiener filters were then derived for speech enhancement.
Here, we present only the RVAE model.

10.1 Generative model

The RVAE model was designed to model speech signals in the short-term
Fourier transform (STFT) domain. This implies that the model applies
to a sequence of *complex-valued* vectors. Therefore, the observation
model uses a multivariate zero-mean circular complex Gaussian distri-
bution (Neeser and Massey, 1993), denoted $N_c$, instead of the usual
multivariate real-valued Gaussian distribution. This observation model
has the following generic form:

$$ \sigma_{\theta_x}(z_T) = d_x(z_T), \qquad (10.1) $$

$$ p_{\theta_x}(\mathbf{x}_t|z_T) = N_c(\mathbf{x}_t; \mathbf{0}, \text{diag}\{\sigma_{\theta_x}^2(z_T)\}), \qquad (10.2) $$

where $T$ denotes a set of time frames, and the following three cases are considered: i) an instantaneous model: $T = \{t\}$, which only considers the current latent state vector to model the observation at time $t$; ii) a causal model: $T = \{1:t\}$, which considers the sequence of past and present latent state vectors; and iii) a noncausal model: $T = \{1:T\}$, which considers the complete sequence of latent state vectors.

This model can be adapted to real-valued observations with a usual Gaussian distribution:<fnref n="1" /> we just have to replace $N_c$ with $N$, replace $\mathbf{0}$ with a mean parameter $\mu_{\theta}(z_T)$ in (10.2), and add this mean parameter to the left-hand side of (10.1), as usual. This is what we have done hereinafter for easier comparison with the other models.

As in STORN, $z_t$ is assumed i.i.d. with a standard Gaussian distribution:

$$ p(\mathbf{z}_{1:T}) = \prod_{t=1}^{T} p(\mathbf{z}_t) \quad \text{with} \quad p(\mathbf{z}_t) = N(\mathbf{z}_t; \mathbf{0}, \mathbf{I}_L). \qquad (10.3) $$

Therefore, there is no explicit temporal model on $z_t$, and $x_t$ possibly depends on the past and future values of the latent state through (10.2). We have here $\theta_z = \emptyset$ and $\theta = \theta_x$. As case i) is strictly equivalent to the original VAE of Section 2, with no temporal model at all, we will now focus on cases ii) and iii).

Leglaive et al. (2020) only mentioned that cases ii) and iii) are implemented using a forward RNN and a bidirectional RNN, respectively, which take as input the sequence $z_{1:t}$ or $z_{1:T}$, respectively. The authors did not provide detailed implementation equations (though they provided a link to some supplementary material, including informative schemas). Let us write them now for easier comparison with the other models (for the same reason, we consider real-valued observations).

<footnote n="1">¹Or any other distribution for real-valued vectors, as already mentioned. In fact, under some conditions, the complex proper Gaussian distribution applied on STFT coefficients corresponds to a Gamma distribution on the squared magnitude of those coefficients (Girin et al., 2019).</footnote>
<figure><img src="image_20.png" /><figcaption>Figure 10.1: Causal RVAE’s graphical model in developed form (left) and compact form (right).</figcaption></figure>

**Causal case:** Let us start with the causal case, for which we have

$$h_t = d_h(z_t, h_{t-1}), \qquad (10.4)$$

$$[\mu_{\theta}(h_t), \sigma_{\theta}(h_t)] = d_{\chi}(h_t), \qquad (10.5)$$

$$p_{\theta}(\mathbf{x}_t|\mathbf{h}_t) = \mathcal{N}(\mathbf{x}_t; \mu_{\theta}(\mathbf{h}_t), \mathrm{diag}\{\sigma_{\theta}^2(\mathbf{h}_t)\}). \qquad (10.6)$$

Eq. (10.4) is similar to the RNN internal state update (3.3) or (9.1), with the major difference being that the latent state $z_t$ is used as an input instead of an external input $u_t$ or previous observation vector $x_{t-1}$. Alternatively, (10.4) can be viewed as a simplified version of the STORN or VRNN internal state updates (7.1) or (8.1), where only $z_t$ and $h_{t-1}$ (and not $x_{t-1}$) are used as inputs. Considering both the observation model and the prior latent-state model, the causal RVAE model is quite close to STORN. The two differences with STORN are that here $x_{t-1}$ is not reinjected as input to the internal state $h_t$ and an LSTM network is used instead of a single-layer RNN in the original STORN formulation.

The graphical model of RVAE (causal case) is shown in Figure 10.1. As is now usual in our developments, we rewrite $h_t = h_t(z_{1:t})^2$ and have $p_{\theta}(\mathbf{x}_t|\mathbf{h}_t) = p_{\theta}(\mathbf{x}_t|\mathbf{h}_t(z_{1:t}))$. For a complete data sequence, we have

$$p_{\theta}(\mathbf{x}_{1:T}, \mathbf{z}_{1:T}) = \prod_{t=1}^{T} p_{\theta}(\mathbf{x}_{t}|\mathbf{h}_{t}(\mathbf{z}_{1:t}))p(\mathbf{z}_{t}), \qquad (10.7)$$

<footnote n="2">²This is also a function of $h_0$, which we omit for clarity.</footnote>
<figure><img src="image_21.png" /><figcaption>**Figure 10.2:** Noncausal RVAE's graphical model in developed form (left) and compact form (right).</figcaption></figure>

which, as for the other models, can be reshaped into the following more general expression:

$$p_{\theta}(\mathbf{x}_{1:T}, \mathbf{z}_{1:T}) = \prod_{t=1}^{T} p_{\theta_x}(\mathbf{x}_t | \mathbf{z}_{1:t})p(\mathbf{z}_t). \qquad (10.8)$$

**Noncausal case:** All DVAE models we have seen so far are causal (at generation). The noncausal case presented by Leglaive et al. (2020) is the first noncausal DVAE model found in the literature. It is implemented with a combination of a forward RNN and a backward RNN on $\mathbf{z}_t$:

$$\vec{\mathbf{h}}_t = d_{\vec{\mathbf{h}}}(\mathbf{z}_t, \vec{\mathbf{h}}_{t-1}), \qquad (10.9)$$

$$\leftarrow_{\mathbf{h}} h_t = d_{\leftarrow_{\mathbf{h}}}(\mathbf{z}_t, \leftarrow_{\mathbf{h}} h_{t+1}), \qquad (10.10)$$

$$\mathbf{h}_t = [\vec{\mathbf{h}}_t, \leftarrow_{\mathbf{h}} h_t], \qquad (10.11)$$

$$[\mu_{\theta_x}(\mathbf{h}_t), \sigma_{\theta_x}(\mathbf{h}_t)] = d_x(\mathbf{h}_t), \qquad (10.12)$$

$$p_{\theta_x}(\mathbf{x}_t|\mathbf{h}_t) = \mathcal{N}(\mathbf{x}_t; \mu_{\theta_x}(\mathbf{h}_t), \mathrm{diag}\{\sigma_{\theta_x}^2(\mathbf{h}_t)\}). \qquad (10.13)$$

We thus have $\mathbf{h}_t = \mathbf{h}_t(\mathbf{z}_{1:T}).^3$ The graphical representation of the

<footnote n="3">³This function is also a function of the initial internal states $\vec{\mathbf{h}}_0$ and $\leftarrow_{\mathbf{h}} h_{T+1}$,</footnote>
noncausal RVAE model is shown in Figure 10.2. For a complete data sequence, we have

$$p_{\theta}(\mathbf{x}_{1:T}, \mathbf{z}_{1:T}) = \prod_{t=1}^{T} p_{\theta}(\mathbf{x}_t | \mathbf{h}_t(\mathbf{z}_{1:T})) p(\mathbf{z}_t), \quad (10.14)$$

which can be reshaped into

$$p_{\theta}(\mathbf{x}_{1:T}, \mathbf{z}_{1:T}) = \prod_{t=1}^{T} p_{\theta_x}(\mathbf{x}_t | \mathbf{z}_{1:T}) p(\mathbf{z}_t). \quad (10.15)$$

## 10.2 Inference model

As for the inference model, Leglaive et al. (2020) first remarked that, using the chain rule and D-separation, the posterior distribution of the latent vectors can be expressed as follows:

$$p_{\theta}(\mathbf{z}_{1:T} | \mathbf{x}_{1:T}) = \prod_{t=1}^{T} p_{\theta}(\mathbf{z}_t | \mathbf{z}_{1:t-1}, \mathbf{x}_{T'}), \quad (10.16)$$

where in the causal case, $T' = \{t: T\}$, and in the noncausal case, $T' = \{1: T\}$. For the causal generative model, the latent vector at a given time step depends (a posteriori) on past latent vectors and on present and future observations, whereas for the noncausal generative model, it also depends on the past observations. Therefore, the authors chose to define the variational distribution $q_{\phi}$ with the same form:

$$q_{\phi}(\mathbf{z}_{1:T} | \mathbf{x}_{1:T}) = \prod_{t=1}^{T} q_{\phi}(\mathbf{z}_t | \mathbf{z}_{1:t-1}, \mathbf{x}_{T'}). \quad (10.17)$$

As for the generative model, we now detail the implementation of the inference model.

**Causal case:** The inference corresponding to the causal generative model is implemented by combining a forward RNN on the latent

<footnote>which we omit for clarity.</footnote>
<figure><img src="image_22.png" /><figcaption>**Figure 10.3:** Graphical model of causal RVAE at inference time in developed form (left) and compact form (right).</figcaption></figure>

vectors and a backward RNN on the observations:

$$ \vec{\mathbf{g}}_t = e_{\vec{\mathbf{g}}} (\mathbf{z}_{t-1}, \vec{\mathbf{g}}_{t-1}), \qquad (10.18) $$

$$ \overleftarrow{\mathbf{g}}_t = e_{\overleftarrow{\mathbf{g}}} (\mathbf{x}_t, \overleftarrow{\mathbf{g}}_{t+1}), \qquad (10.19) $$

$$ \mathbf{g}_t = [\vec{\mathbf{g}}_t, \overleftarrow{\mathbf{g}}_t], \qquad (10.20) $$

$$ [\mu_{\phi}(\mathbf{g}_t), \sigma_{\phi}(\mathbf{g}_t)] = e_{\mathbf{z}}(\mathbf{g}_t), \qquad (10.21) $$

$$ q_{\phi}(\mathbf{z}_t|\mathbf{z}_{1:t-1}, \mathbf{x}_{t:T}) = \mathcal{N}(\mathbf{z}_t; \mu_{\phi}(\mathbf{g}_t), \mathrm{diag}\{\boldsymbol{\sigma}_{\phi}^2(\mathbf{g}_t)\}), \qquad (10.22) $$

$$ q_{\phi}(\mathbf{z}_{1:T}|\mathbf{x}_{1:T}) = \prod_{t=1}^{T} q_{\phi}(\mathbf{z}_{t}|\mathbf{z}_{1:t-1}, \mathbf{x}_{t:T}). \qquad (10.23) $$

This inference model is shown in Figure 10.3.

**Noncausal case:** The inference corresponding to the noncausal generative model is similar to the causal case, except that the RNN on the
<figure><img src="image_23.png" /><figcaption>**Figure 10.4:** Graphical model of noncausal RVAE at inference time in developed form (left) and compact form (right).</figcaption></figure>

observations is bidirectional:

$$ \vec{g}_t^z = e_{\vec{g}^z}(z_{t-1}, \vec{g}_{t-1}^z), \qquad (10.24) $$

$$ \vec{g}_t^x = e_{\vec{g}^x}(x_t, \vec{g}_{t-1}^x), \qquad (10.25) $$

$$ \vec{g}_t^x = e_{\vec{g}^x}(x_t, \vec{g}_{t+1}^x), \qquad (10.26) $$

$$ g_t = [\vec{g}_t^z, \vec{g}_t^x, \vec{g}_t^x], \qquad (10.27) $$

$$ [\mu_{\phi}(g_t), \sigma_{\phi}(g_t)] = e_z(g_t), \qquad (10.28) $$

$$ q_{\phi}(\mathbf{z}_t|\mathbf{z}_{1:t-1}, \mathbf{x}_{1:T}) = \mathcal{N}(\mathbf{z}_t; \mu_{\phi}(\mathbf{g}_t), \mathrm{diag}\{\boldsymbol{\sigma}_{\phi}^2(\mathbf{g}_t)\}), \qquad (10.29) $$

$$ q_{\phi}(\mathbf{z}_{1:T}|\mathbf{x}_{1:T}) = \prod_{t=1}^{T} q_{\phi}(\mathbf{z}_t|\mathbf{z}_{1:t-1}, \mathbf{x}_{1:T}). \qquad (10.30) $$

This inference model is shown in Figure 10.4.
## 10.3 Training

For conciseness, and because we focus on reviewing causal DVAEs, we only describe in this section the VLB for the causal RVAE model. The methodology to derive the VLB in the noncausal case is similar.

A comparison of the compact form of causal RVAE in (10.8) with the general compact form of a DVAE in (4.4) (simplified without $u_{1:T}$) shows that the causal RVAE model makes the following conditional independence assumptions:

$$ 
\begin{aligned}
p_{\theta_x}(\mathbf{x}_t|\mathbf{x}_{1:t-1}, \mathbf{z}_{1:t}) &= p_{\theta_x}(\mathbf{x}_t|\mathbf{z}_{1:t}); \\
p_{\theta_z}(\mathbf{z}_t|\mathbf{x}_{1:t-1}, \mathbf{z}_{1:t-1}) &= p(\mathbf{z}_t).
\end{aligned} 
\quad (10.31) 
$$

Using these two simplifications, along with the inference model (10.23) (which we recall is consistent with the exact posterior distribution), the VLB can be simplified as follows:

$$ 
\begin{aligned}
\mathcal{L}(\theta, \phi; \mathbf{x}_{1:T}) = & \sum_{t=1}^{T} \mathbb{E}_{q_{\phi}(\mathbf{z}_{1:t}|\mathbf{x}_{1:T})} [\log p_{\theta_x}(\mathbf{x}_t|\mathbf{z}_{1:t})] \\
& - \sum_{t=1}^{T} \mathbb{E}_{q_{\phi}(\mathbf{z}_{1:t-1}|\mathbf{x}_{1:T})} [D_{\text{KL}}(q_{\phi}(\mathbf{z}_t|\mathbf{z}_{1:t-1}, \mathbf{x}_{1:T}) \| p(\mathbf{z}_t))] .
\end{aligned} 
\quad (10.32) 
$$

As for the previous models, the KL divergence can be computed analytically, while the two intractable expectations are approximated by Monte Carlo estimates using samples drawn from $q_{\phi}(\mathbf{z}_{1:t}|\mathbf{x}_{1:T})$ and $q_{\phi}(\mathbf{z}_{1:t-1}|\mathbf{x}_{1:T})$ in a recursive manner.
# 11

## Disentangled Sequential Autoencoders

Li and Mandt (2018) proposed a hierarchical model called DSAE. This model introduces the idea of adding to the usual sequence of latent variables $z_{1:T}$ a sequence-level latent vector $v$ (denoted $f$ by Li and Mandt (2018)), which is assumed to encode the sequence-level characteristics of the data. Therefore, $z_t$ is assumed to encode time-dependent data features (e.g., the dynamics of an object in a video clip), and $v$ is assumed to encode “everything else” (e.g., object characteristics in a video clip).

### 11.1 Generative model

Li and Mandt (2018) only provided the following general form of the generative DSAE model for a complete data sequence:

$$p_{\theta}(\mathbf{x}_{1:T}, \mathbf{z}_{1:T}, \mathbf{v}) = p_{\theta_v}(\mathbf{v}) \prod_{t=1}^{T} p_{\theta_x}(\mathbf{x}_t|\mathbf{z}_t, \mathbf{v}) p_{\theta_z}(\mathbf{z}_t|\mathbf{z}_{1:t-1}). \quad (11.1)$$

More detailed information about the pdfs and implementation issues are found in annexes from the ArXiv version of the paper. The authors used different variants for different datasets according to the nature of the data (e.g., video clips or speech signals). Here, we only report the model
implemented for speech signals. The dynamical model $p_{\theta_z}(\mathbf{z}_t|\mathbf{z}_{1:t-1})$ is a Gaussian distribution whose parameters are provided by an LSTM network. The observation model $p_{\theta_x}(\mathbf{x}_t|\mathbf{z}_t, \mathbf{v})$ is a Gaussian distribution whose parameters are provided by a feed-forward DNN. With the simplified generic RNN formalism used for LSTM (see Section 3.1.1), we can thus write

$$ \mathbf{h}_t = d_h(\mathbf{z}_{t-1}, \mathbf{h}_{t-1}), \qquad (11.2) $$

$$ [\mu_{\theta_z}(\mathbf{h}_t), \sigma_{\theta_z}(\mathbf{h}_t)] = d_z(\mathbf{h}_t), \qquad (11.3) $$

$$ p_{\theta_z}(\mathbf{z}_t|\mathbf{h}_t) = \mathcal{N}(\mathbf{z}_t; \mu_{\theta_z}(\mathbf{h}_t), \mathrm{diag}\{\sigma_{\theta_z}^2(\mathbf{h}_t)\}), \qquad (11.4) $$

$$ [\mu_{\theta_x}(\mathbf{z}_t, \mathbf{v}), \sigma_{\theta_x}(\mathbf{z}_t, \mathbf{v})] = d_x(\mathbf{z}_t, \mathbf{v}), \qquad (11.5) $$

$$ p_{\theta_x}(\mathbf{x}_t|\mathbf{z}_t, \mathbf{v}) = \mathcal{N}(\mathbf{x}_t; \mu_{\theta_x}(\mathbf{z}_t, \mathbf{v}), \mathrm{diag}\{\sigma_{\theta_x}^2(\mathbf{z}_t, \mathbf{v})\}). \qquad (11.6) $$

The graphical representation of DSAE is shown in Figure 11.1. This model is similar to a DKF in the undriven mode conditioned on variable **v**, except that, in addition to this conditioning, the first-order Markov temporal model of DKF is replaced with a virtually infinite-order model owing to the LSTM. Although this is poorly discussed in the DVAE papers in general, it is an example of interesting model extensions that are easy to implement in the deep learning and VAE frameworks. As stated by Krishnan et al. (2015), “using deep neural networks, we can enhance Kalman filters with arbitrarily complex transition dynamics and emission distributions. [...] we can tractably learn such models by optimizing a bound on the likelihood of the data.”

## 11.2 Inference model

For DSAE, the posterior distribution of latent variables is given by

$$ p_{\theta}(\mathbf{z}_{1:T}, \mathbf{v}|\mathbf{x}_{1:T}) = p_{\theta}(\mathbf{v}|\mathbf{x}_{1:T})p_{\theta}(\mathbf{z}_{1:T}|\mathbf{v}, \mathbf{x}_{1:T}) \qquad (11.7) $$

$$ = p_{\theta}(\mathbf{v}|\mathbf{x}_{1:T}) \prod_{t=1}^{T} p_{\theta}(\mathbf{z}_t|\mathbf{z}_{1:t-1}, \mathbf{v}, \mathbf{x}_{1:t}) \qquad (11.8) $$

$$ = p_{\theta}(\mathbf{v}|\mathbf{x}_{1:T}) \prod_{t=1}^{T} p_{\theta}(\mathbf{z}_t|\mathbf{z}_{1:t-1}, \mathbf{v}, \mathbf{x}_{t:T}). \qquad (11.9) $$
<figure><img src="image_24.png" /><figcaption>**Figure 11.1:** DSAE's graphical model in developed form (left) and compact form (right).</figcaption></figure>

The simplification in the last line results from D-separation. This decomposition can be interpreted as follows: The whole sequence of observations $\mathbf{x}_{1:T}$ is used to estimate the “object” representation $\mathbf{v}$, and then $\mathbf{v}$, the present and future observations $\mathbf{x}_{t:T}$, and previous latent state vectors $\mathbf{z}_{1:t-1}$ are used to update the object dynamics.

As for the approximate posterior distribution $q_{\phi}$, Li and Mandt (2018) proposed two models. The first one, referred to as “factorized,” is expressed as

$$q_{\phi}(\mathbf{z}_{1:T}, \mathbf{v}|\mathbf{x}_{1:T}) = q_{\phi_{\mathbf{v}}}(\mathbf{v}|\mathbf{x}_{1:T}) \prod_{t=1}^{T} q_{\phi_{\mathbf{z}}}(\mathbf{z}_{t}|\mathbf{x}_{t}). \quad (11.10)$$

This model thus relies on an instantaneous frame-wise inference model $q_{\phi}(\mathbf{z}_t|\mathbf{x}_t)$ for encoding the latent vector dynamics. This approach is oversimplistic compared to the exact posterior distribution and yields inferior performance to that of the second inference model. We thus focus on the latter, which is referred to as “full,” and is given by

$$q_{\phi}(\mathbf{z}_{1:T}, \mathbf{v}|\mathbf{x}_{1:T}) = q_{\phi_{\mathbf{v}}}(\mathbf{v}|\mathbf{x}_{1:T})q_{\phi_{\mathbf{z}}}(\mathbf{z}_{1:T}|\mathbf{v}, \mathbf{x}_{1:T}). \quad (11.11)$$

As the authors mentioned, “The idea behind [this] structured approximation is that content may affect dynamics.” So far, this model has been compliant with the exact posterior distribution, as expressed by (11.7). From the information given in the annexes of the ArXiv version of the paper, we can write the detailed equations of the full inference
model as follows:

$$ \vec{\mathbf{g}}_t^\mathbf{v} = e_{\vec{\mathbf{g}}^\mathbf{v}}(\mathbf{x}_t, \vec{\mathbf{g}}_{t-1}^\mathbf{v}), \qquad (11.12) $$

$$ \vec{\mathbf{g}}_t^\mathbf{v} = e_{\vec{\mathbf{g}}^\mathbf{v}}(\mathbf{x}_t, \vec{\mathbf{g}}_{t+1}^\mathbf{v}), \qquad (11.13) $$

$$ \mathbf{g}^\mathbf{v} = [\vec{\mathbf{g}}_T^\mathbf{v}, \vec{\mathbf{g}}_1^\mathbf{v}], \qquad (11.14) $$

$$ [\mu_{\phi_v}(\mathbf{g}^\mathbf{v}), \sigma_{\phi_v}(\mathbf{g}^\mathbf{v})] = e_v(\mathbf{g}^\mathbf{v}), \qquad (11.15) $$

$$ q_{\phi_v}(\mathbf{v}|\mathbf{x}_{1:T}) = \mathcal{N}(\mathbf{v}; \mu_{\phi_v}(\mathbf{g}^\mathbf{v}), \text{diag}\{\boldsymbol{\sigma}_{\phi_v}^2(\mathbf{g}^\mathbf{v})\}), \qquad (11.16) $$

$$ \vec{\mathbf{g}}_t^\mathbf{z} = e_{\vec{\mathbf{g}}^\mathbf{z}}([\mathbf{x}_t, \mathbf{v}], \vec{\mathbf{g}}_{t-1}^\mathbf{z}), \qquad (11.17) $$

$$ \vec{\mathbf{g}}_t^\mathbf{z} = e_{\vec{\mathbf{g}}^\mathbf{z}}([\mathbf{x}_t, \mathbf{v}], \vec{\mathbf{g}}_{t+1}^\mathbf{z}), \qquad (11.18) $$

$$ \mathbf{g}_t^\mathbf{z} = [\vec{\mathbf{g}}_t^\mathbf{z}, \vec{\mathbf{g}}_t^\mathbf{z}], \qquad (11.19) $$

$$ [\mu_{\phi_z}(\mathbf{g}_t^\mathbf{z}), \sigma_{\phi_z}(\mathbf{g}_t^\mathbf{z})] = e_z(\mathbf{g}_t^\mathbf{z}), \qquad (11.20) $$

$$ q_{\phi_z}(\mathbf{z}_t|\mathbf{v}, \mathbf{x}_{1:T}) = \mathcal{N}(\mathbf{z}_t; \mu_{\phi_z}(\mathbf{g}_t^\mathbf{z}), \text{diag}\{\boldsymbol{\sigma}_{\phi_z}^2(\mathbf{g}_t^\mathbf{z})\}), \qquad (11.21) $$

and for the full sequence $\mathbf{z}_{1:T}$ we have

$$ q_{\phi_z}(\mathbf{z}_{1:T}|\mathbf{v}, \mathbf{x}_{1:T}) = \prod_{t=1}^{T} q_{\phi_z}(\mathbf{z}_t|\mathbf{v}, \mathbf{x}_{1:T}). \qquad (11.22) $$

Note that none of the two approximations proposed by the authors (fقطorized and full) actually follow the dependencies of the exact posterior distribution shown in (11.9). The graphical representation of the full inference model is shown in Figure 11.2.<fnref n="1" />

## 11.3 Training

To derive the VLB for the DSAE model, we apply the same strategy as that applied for the other models (i.e., inject the generative model and the approximate posterior in the VLB general formulation). We do so

<footnote n="1">¹In Li and Mandt’s (2018) paper, Appendix A, the schematic representation of the inference model given in Figure 9(b) is not consistent with the following sentence (reported with our notations): “Finally the parameters of $q_\phi(\mathbf{z}_{1:T}|\mathbf{v}, \mathbf{x}_{1:T})$ are computed by a simple RNN with input $[\vec{\mathbf{g}}_t^\mathbf{z}, \vec{\mathbf{g}}_t^\mathbf{z}]$ at time $t$. It is indeed inconsistent and a bit odd that $\mathbf{z}_{t-1}$ is not mentioned as an input of the $\mathbf{z}_t$ inference process, as is well apparent in Figure 9(b). We base the writing of (11.19)–(11.21) on their text and not on their figure.</footnote>