# Fully Bayesian Autoencoders with Latent Sparse Gaussian Processes

Ba-Hien Tran  Babak Shahbaba  Stephan Mandt  Maurizio Filippone

###### Abstract

Autoencoders and their variants are among the most widely used models in
representation learning and generative modeling. However, autoencoder-based
models usually assume that the learned representations are i.i.d. and fail to
capture the correlations between the data samples. To address this issue, we
propose a novel Sparse Gaussian Process Bayesian Autoencoder (SGP-BAE) model
in which we impose fully Bayesian sparse Gaussian Process priors on the latent
space of a Bayesian Autoencoder. We perform posterior estimation for this
model via stochastic gradient Hamiltonian Monte Carlo. We evaluate our
approach qualitatively and quantitatively on a wide range of representation
learning and generative modeling tasks and show that our approach consistently
outperforms multiple alternatives relying on Variational Autoencoders.

Machine Learning, ICML

Table 1: A summary of related methods. Here,
${{\bm{\mathbf{\theta}}}},{{\bm{\mathbf{u}}}},{{\bm{\mathbf{S}}}}$ refer to gp
hyper-parameters, inducing variables and inducing inputs, respectively. $N$
and $B$ are the number of data points and the mini-batch size, whereas $M,H\ll
N$ are the number of inducing points and the low-rank matrix factor,
respectively. The colors denoting the methods shall be used consistently
throughout the paper.

Model |  | Scalable
---
(Minibatching)
| Non i.i.d.
---
data
| Free-form
---
posterior
| gp
---
complexity
| Arbitrary kernel
---
& data type
| Learnable
---
gp
| Inference
---
${{\bm{\mathbf{\theta}}}},{{\bm{\mathbf{u}}}},{{\bm{\mathbf{S}}}}$
Reference
( ) vae | ✓ | ✗ | ✗ | - | - | - | - | Kingma & Welling (2014)
( ) cvae | ✓ | ✗ | ✗ | - | - | - | - | Sohn et al. (2015)
( ) gppvae | ✓ | ✓ | ✗ | $\mathcal{O}(NH^{2})$ | ✗ | ✗ | ✗ | Casale et al. (2018)
( ) gpvae | ✗ | ✓ | ✗ | $\mathcal{O}(N^{3})$ | ✓ | ✗ | ✗ | Pearce (2019)
( ) gp-vae | ✓ | ✓ | ✗ | $\mathcal{O}(N)$ | ✗ | ✗ | ✗ | Fortuin et al. (2020)
( ) sgp-vae | ✓ | ✓ | ✗ | $\mathcal{O}(BM^{2}+M^{3})$ | ✓ | ✓ | ✗ | Ashman et al. (2020)
( ) svgp-vae | ✓ | ✓ | ✗ | $\mathcal{O}(BM^{2}+M^{3})$ | ✓ | ✓ | ✗ | Jazbec et al. (2021)
( ) bae | ✓ | ✗ | ✓ | - | - | - | - | This work
( ) gp-bae | ✗ | ✓ | ✓ | $\mathcal{O}(N^{3})$ | ✓ | ✗ | ✗ | This work
( ) sgp-bae | ✓ | ✓ | ✓ | $\mathcal{O}(BM^{2}+M^{3})$ | ✓ | ✓ | ✓ | This work

##  1 Introduction

The problem of learning representations of data that are useful for downstream
tasks is a crucial factor in the success of many machine learning applications
(Bengio et al., 2013). Among the numerous proposed solutions, modeling
approaches that evolved from autoencoders (aes) (Cottrell et al., 1989) are
particularly appealing, as they do not require annotated data and have proven
effective in unsupervised learning tasks, such as data compression and
generative modeling (Tomczak, 2022; Yang et al., 2022). aes are neural
networks consisting of an encoder that maps input data to a set of lower-
dimensional latent variables and a decoder that maps the latent variables back
to the observations.

In applications where data is scarce or uncertainty quantification and
learning interpretable representations are crucial, it is beneficial to treat
these models in a Bayesian manner (Mackay, 1992; Neal, 1996; Wilson &
Izmailov, 2020; Izmailov et al., 2021a) by imposing meaningful prior
distributions over both their latent spaces as well as their parameters, i.e.,
the parameters of the encoder/decoder. However, a fully Bayesian treatment of
these models, which we refer to as Bayesian autoencoders (baes), is
challenging due to the huge number of local (per-datapoint) latent variables
to be inferred. This difficulty has motivated the development of Variational
Autoencoders (vaes) (Kingma & Welling, 2014; Rezende et al., 2014b) that make
inference tractable by amortization. Note that vaes do not capture Bayesian
uncertainty in the decoder parameters.

A critical limitation of standard vaes is that latent encodings are assumed to
be factorized across latent dimensions. However, capturing latent correlations
is often necessary to model structured data. For example, in autonomous
driving or medical imaging applications, high-dimensional images are
correlated in time. Spatio-temporal dependencies between samples are also
common in environmental and life sciences datasets.

To address this limitation, several works (Pearce, 2019; Casale et al., 2018)
have attempted to introduce Gaussian process (gp) priors over the latent space
of vaes that capture correlations between pairs of latent variables through a
kernel function. While gp priors outperform conventional priors on many tasks,
they also introduce computational challenges such as $\mathcal{O}(N^{3})$
complexity for gp inference, where $N$ is the number of data instances.
Recently, Jazbec et al. (2021) proposed the svgp-vae model to tackle this
computational issue by relying on sparse approximations which summarize the
dataset into a set of points referred to as inducing points (Quiñonero-Candela
& Rasmussen, 2005).

Although the svgp-vae model (Jazbec et al., 2021) has achieved promising
results, it has several significant drawbacks. First, similarly to vaes, svgp-
vae is strongly tied to a variational inference (vi) formulation (Jordan et
al., 1999; Zhang et al., 2018), which can lead to poor approximations due to
vi making strong assumptions on both the factorization and functional form of
the posterior. Second, the svgp-vae model follows the common practice in the
sparse gp literature of optimizing the inducing points and kernel
hyperparameters based on the marginal likelihood. However, this approach does
not account for uncertainty in the inducing inputs and hyperparameters. It is
well known that this can result in biased estimates and underestimated
predictive uncertainties (Rossi et al., 2021; Lalchand et al., 2022). In this
work, we propose a novel Sparse Gaussian Process Bayesian Autoencoder (sgp-
bae) model that addresses these issues by providing scalable and flexible
inference through a fully Bayesian treatment.

#### Contributions.

Specifically, first, we develop a fully Bayesian autoencoder (bae) model (§
3), where we adopt a Bayesian treatment for both the local (per-datapoint)
latent variables and the global decoder parameters. This approach differs from
vaes in that it allows specifying any prior over the latent space while
decoupling the model from the inference. As a result, we can rely on powerful
alternatives to vi to carry out inference, and we adopt stochastic gradient
Hamiltonian Monte Carlo (sghmc) (Chen et al., 2014) as a scalable solution. To
achieve this, we propose an amortized Markov chain Monte Carlo (mcmc) approach
for our bae model by using an implicit stochastic network as the encoder to
learn to draw samples from the posterior of local latent variables. Our
approach addresses the prohibitively expensive inference cost induced by the
local latent variables and avoids making strong assumptions on the form of the
posterior. Second, when imposing a gp prior over the latent space, we propose
a novel scalable sgp-bae model (§ 4) in which the inducing points and kernel
hyperparameters of the sparse gp prior on the latent space are treated in a
fully Bayesian manner. This model offers attractive features such as high
scalability, richer modeling capability, and improved prediction quality.
Third, we extend sgp-bae model to allow deep gp prior (Damianou & Lawrence,
2013) and handling missing data. To the best of our knowledge, this is the
first work to consider deep gp prior to aes. Finally, we conduct a rigorous
evaluation of our sgp-bae model through a variety of experiments on dynamic
representation learning and generative modeling. The results demonstrate
excellent performance compared to existing methods of combining gps and aes (§
5).

##  2 Related Work

#### Autoencoders.

aes (Cottrell et al., 1989) are a powerful framework for representation
learning by projecting data into a lower-dimensional latent space using an
encoder-decoder architecture. vaes (Kingma & Welling, 2014; Rezende et al.,
2014b) elegantly combines ae with variational inference enabling the model to
generate new data and allowing for the specification of any prior on the
latent space. To improve the performance of vaes, recent works have attempted
to employ flexible priors such as mixture distributions (Tomczak & Welling,
2018; Bauer & Mnih, 2019), normalizing flows (Chen et al., 2017), non-
parametric model (Nalisnick & Smyth, 2017), or Gaussian processes (Casale et
al., 2018). In this work, we take an orthogonal view as we consider Bayesian
autoencoders, which recently have shown promising performance (Tran et al.,
2021), but still lack a mechanism to impose priors on the latent space.

#### Gaussian process priors for AE models.

The earliest attempts to combine ae models with gps are the gp prior vae
(gppvae) (Casale et al., 2018) and gp-vae (Pearce, 2019). Both these models
are not scalable for generic kernel choices and data types. gppvae is
restricted to a view-object gp product kernel and employs a Taylor
approximation for gps, while gp-vae relies on exact gp inference. Recently,
Fortuin et al. (2020) and Zhu et al. (2022) propose gp-vae and markovian-
gpvae, respectively, that are indeed scalable (even in linear time complexity)
by exploiting the Markov assumption but they only work exclusively for time-
series data. Most closely to our method is the approach of Jazbec et al.
(2021) (svgp-vae), Ashman et al. (2020) (sgp-vae) and Ramchandran et al.
(2021), where they utilize inducing point methods Titsias (2009); Hensman et
al. (2013) to make gps scalable. However, all these methods strongly relied on
vaes and a variational formulation for gps. In this work, we take a completely
different route, as we aim to treat sparse gps and aes in a fully Bayesian way
while enjoying scalability thanks to recent advances in stochastic-gradient
mcmc sampling. Table 1 compares our proposed models with the key related
works.

##  3 Imposing Distributions over the Latent Space of Bayesian Autoencoders

(a) bae.

(b) sgp-bae.

Figure 1:  The graphical models of vanilla bae (a), and the proposed sgp-bae
with a fully Bayesian sparse gp prior imposed on the latent space. The solid
lines denote the generative part, whereas the dashed lines denote the encoding
part. The cyclic thick line represents that the latent gp correlates with
every latent code.

We are interested in unsupervised learning problems with a high-dimensional
dataset consisting of $N$ data points
${{\bm{\mathbf{Y}}}}\stackrel{{\scriptstyle\text{def}}}{{=}}[{{\bm{\mathbf{y}}}}_{1},\cdots,{{\bm{\mathbf{y}}}}_{N}]^{\top}\in\mathbb{R}^{N\times
P}$. Each data point has a corresponding low-dimensional auxiliary data entry,
summarized as
${{\bm{\mathbf{X}}}}\stackrel{{\scriptstyle\text{def}}}{{=}}[{{\bm{\mathbf{x}}}}_{1},\cdots,{{\bm{\mathbf{x}}}}_{N}]^{\top}\in\mathbb{R}^{N\times
D}$. For instance, ${{\bm{\mathbf{y}}}}_{i}$ could be a video frame and
${{\bm{\mathbf{x}}}}_{i}$ the corresponding time stamp. As another example,
consider electronic health record (EHR) data, where the auxiliary data could
include patient covariate information, such as age, height, weight, sex, and
time since remission. Finally, we denote by
${{\bm{\mathbf{Z}}}}\stackrel{{\scriptstyle\text{def}}}{{=}}[{{\bm{\mathbf{z}}}}_{1},\cdots,{{\bm{\mathbf{z}}}}_{N}]^{\top}\in\mathbb{R}^{N\times
C}$ the low-dimensional latent representation of the data, meaning that each
latent variable ${{\bm{\mathbf{z}}}}_{i}$ lives in a $C$-dimensional latent
space. We aim to train a model that is able to (1) generate
${{\bm{\mathbf{Y}}}}$ conditioned on ${{\bm{\mathbf{X}}}}$, and (2) provide
useful and interpretable low-dimensional representations of
${{\bm{\mathbf{Y}}}}$.

#### Model setup.

In this work, we consider a model based on aes, and we aim to treat this in a
fully Bayesian manner. This treatment promises improved predictions, reliable
uncertainty estimates, and increased robustness under data sparsity (Mackay,
1992; Izmailov et al., 2021b). One difficulty in doing so is that the prior
distribution over the latent variables would be determined by the prior over
the weights of the encoder and not a distribution of interest. In many
applications, including the ones considered here, this is undesirable, and the
goal is to impose a certain prior distribution over the latent representation
in a similar vein as vaes and their variants. Therefore, we propose to treat
the entire ae in a fully Bayesian manner except for the encoder, and to design
the encoder in such a way that it maps data ${{\bm{\mathbf{y}}}}_{i}$ into
corresponding codes ${{\bm{\mathbf{z}}}}_{i}$ while allowing these mappings to
be compatible with posterior samples over the latent codes. For the encoder,
as we will elaborate on shortly, we will employ a so-called stochastic
inference network to learn to draw posterior samples of latent variables
${{\bm{\mathbf{z}}}}_{i}$ given high-dimensional inputs
${{\bm{\mathbf{y}}}}_{i}$, while for the decoder and the latent variables we
employ scalable mcmc techniques.

#### Bayesian treatment of the latent space and the decoder.

In order to retain a fully Bayesian treatment of the latent space and the
decoder, we impose a prior $p({{\bm{\mathbf{\varphi}}}})$ over the decoder’s
parameters, ${{\bm{\mathbf{\varphi}}}}$. In addition, another prior
$p({{\bm{\mathbf{Z}}}}\,|\,{{\bm{\mathbf{X}}}};{{\bm{\mathbf{\theta}}}})$ is
imposed on the latent variables. This prior is conditioned on the auxiliary
data ${{\bm{\mathbf{X}}}}$ and characterized by a set of hyper-parameters
${{\bm{\mathbf{\theta}}}}$. For example, one may employ an uninformative prior
such as an isotropic Gaussian commonly used in standard vaes, but in the next
section we will consider structured priors such as gps and their deep version
as well. We assume that the observed data ${{\bm{\mathbf{Y}}}}$ is fully
factorized and conditional on the latent variables ${{\bm{\mathbf{Z}}}}$ and a
decoder network with parameters ${{\bm{\mathbf{\varphi}}}}$, i.e.,
$p({{\bm{\mathbf{Y}}}}\,|\,{{\bm{\mathbf{Z}}}},{{\bm{\mathbf{\varphi}}}})=\prod_{i=1}^{N}p({{\bm{\mathbf{y}}}}_{i}\,|\,{{\bm{\mathbf{z}}}}_{i},{{\bm{\mathbf{\varphi}}}})$.
The full joint distribution of the model is as follows:

$\displaystyle\vspace{-1ex}p({{\bm{\mathbf{\varphi}}}},{{\bm{\mathbf{Z}}}},{{\bm{\mathbf{Y}}}}\,|\,{{\bm{\mathbf{X}}}})$
$\displaystyle=p({{\bm{\mathbf{\varphi}}}})p({{\bm{\mathbf{Z}}}}\,|\,{{\bm{\mathbf{X}}}};{{\bm{\mathbf{\theta}}}})p({{\bm{\mathbf{Y}}}}\,|\,{{\bm{\mathbf{Z}}}},{{\bm{\mathbf{\varphi}}}}).$
(1)

We wish to infer the posterior over the latent variables and decoder
parameters, which is given by Bayes’ rule:

$\displaystyle
p({{\bm{\mathbf{\varphi}}}},{{\bm{\mathbf{Z}}}}\,|\,{{\bm{\mathbf{Y}}}},{{\bm{\mathbf{X}}}})=\frac{p({{\bm{\mathbf{\varphi}}}},{{\bm{\mathbf{Z}}}},{{\bm{\mathbf{Y}}}}\,|\,{{\bm{\mathbf{X}}}})}{p({{\bm{\mathbf{Y}}}}\,|\,{{\bm{\mathbf{X}}}})},$
(2)

where $p({{\bm{\mathbf{Y}}}}\,|\,{{\bm{\mathbf{X}}}})=\int
p({{\bm{\mathbf{\varphi}}}},{{\bm{\mathbf{Z}}}},{{\bm{\mathbf{Y}}}}\,|\,{{\bm{\mathbf{X}}}})d{{\bm{\mathbf{\varphi}}}}d{{\bm{\mathbf{Z}}}}$
is the marginal likelihood. The generative process of data samples from this
bae model is illustrated in Fig. 1(a).

Characterizing the posterior distribution over
${{\bm{\mathbf{\varphi}}}},{{\bm{\mathbf{Z}}}}$ is analytically intractable
and requires approximations. Given the success of scalable mcmc techniques to
obtain samples from the posterior of model parameters in deep learning models
(Zhang et al., 2020; Tran et al., 2022), in this work, we propose to follow
this practice to obtain samples from
$p({{\bm{\mathbf{\varphi}}}},{{\bm{\mathbf{Z}}}}\,|\,{{\bm{\mathbf{Y}}}},{{\bm{\mathbf{X}}}})$.
In particular, we employ sghmc (Chen et al., 2014), which can scale up to
large datasets by relying on noisy but unbiased estimates of the energy
function (log-posterior)
$U({{\bm{\mathbf{\varphi}}}},{{\bm{\mathbf{Z}}}};{{\bm{\mathbf{X}}}},{{\bm{\mathbf{Y}}}})\stackrel{{\scriptstyle\text{def}}}{{=}}-\log
p({{\bm{\mathbf{\varphi}}}},{{\bm{\mathbf{Z}}}},{{\bm{\mathbf{Y}}}}\,|\,{{\bm{\mathbf{X}}}})$
and without the need to evaluate the energy function over the entire data set.
More precisely, when the prior over the latent codes is fully factorized, we
can approximate this energy function using mini-batches of size $B$ as
follows:

$\displaystyle{U}({{\bm{\mathbf{\varphi}}}},{{\bm{\mathbf{Z}}}};{{\bm{\mathbf{X}}}},{{\bm{\mathbf{Y}}}})$
$\displaystyle\approx\tilde{U}({{\bm{\mathbf{\varphi}}}},{{\bm{\mathbf{Z}}}};{{\bm{\mathbf{X}}}},{{\bm{\mathbf{Y}}}})=-\log
p({{\bm{\mathbf{\varphi}}}})-$ (3)
$\displaystyle-\frac{N}{B}\sum_{i\in\mathcal{I}_{B}}$
$\displaystyle\big{[}\log
p({{\bm{\mathbf{z}}}}_{i}\,|\,{{\bm{\mathbf{x}}}}_{i};{{\bm{\mathbf{\theta}}}})+\log
p({{\bm{\mathbf{y}}}}_{i}\,|\,{{\bm{\mathbf{z}}}}_{i},{{\bm{\mathbf{\varphi}}}})\big{]}$

where $\mathcal{I}_{B}$ is a set of $B$ random indices. The exact procedure
for generating samples from the posterior over
${{\bm{\mathbf{\varphi}}}},{{\bm{\mathbf{Z}}}}$ using sghmc can be found in
Fig. 6.

#### Encoder as a stochastic inference network.

sghmc can be challenging to implement on probabilistic models with many latent
variables due to the high computational burden of iteratively refining the
approximate posterior for each latent variable. Additionally, it can be
difficult to evolve the latent variables for each new test sample. To address
these challenges, we propose using a stochastic neural network as an inference
network to efficiently generate latent codes similar to those generated by the
posterior distribution, inspired by amortized inference techniques (Kingma &
Welling, 2014; Wang & Liu, 2016; Feng et al., 2017; Shi et al., 2019) and mcmc
distillation (Korattikara Balan et al., 2015; Wang et al., 2018; Li et al.,
2017).

More specifically, instead of storing every latent code, we use an inference
network
${{\bm{\mathbf{z}}}}_{i}=f_{{{\bm{\mathbf{\phi}}}}}({{\bm{\mathbf{y}}}}_{i};{\bm{\mathbf{\varepsilon}}})$
with parameters ${{\bm{\mathbf{\phi}}}}$ that generates a corresponding latent
code ${{\bm{\mathbf{z}}}}_{i}$ given an input ${{\bm{\mathbf{y}}}}_{i}$ and a
random seed ${\bm{\mathbf{\varepsilon}}}$. The random seed
${\bm{\mathbf{\varepsilon}}}$ is drawn from a distribution
$q({\bm{\mathbf{\varepsilon}}})$ that is easy to sample from, such as a
uniform or standard Gaussian distribution. The inference network $f_{\phi}$
serves as an encoder by generating posterior samples of the latent code
${{\bm{\mathbf{z}}}}_{i}$ given the observed input ${{\bm{\mathbf{y}}}}_{i}$.
Our approach using this encoder differs from that of vaes, as we do not make
any assumptions about the form of the posterior of the latent variables
${{\bm{\mathbf{z}}}}_{i}$.

We incrementally refine the encoder $f_{{{\bm{\mathbf{\phi}}}}}$ such that its
outputs mimic the sghmc dynamics. Specifically, after every $K$ iterations of
sampling the decoder parameters and the latent codes using sghmc, we adjust
the encoder parameters ${{\bm{\mathbf{\phi}}}}$ based on the following
objective:

$\displaystyle\mathcal{L}\left(\\{{{\bm{\mathbf{z}}}}_{i},{{\bm{\mathbf{y}}}}_{i}\\}_{i\in\mathcal{I}_{B}};{{\bm{\mathbf{\phi}}}}\right)\stackrel{{\scriptstyle\text{def}}}{{=}}\sum_{i\in\mathcal{I}_{B}}\Big{\|}f_{{{\bm{\mathbf{\phi}}}}}({{\bm{\mathbf{y}}}}_{i};\bm{\varepsilon}_{i})-{{\bm{\mathbf{z}}}}_{i}^{(k)}\Big{\|}_{2}^{2},$
(4)

where ${{\bm{\mathbf{z}}}}_{i}^{(k)}$ is the $k$-th posterior sample from
sghmc of the latent variable ${{\bm{\mathbf{z}}}}_{i}$, and it is used as
label to update ${{\bm{\mathbf{\phi}}}}$. As the analytic solution of Eq. 4 is
intractable, we perform $J$ steps of gradient descent to update
${{\bm{\mathbf{\phi}}}}$ using an optimizer such as Adam (Kingma & Ba, 2015).
The inference procedure for baes is described in Algorithm 1.

1

2

Input: Dataset $\\{{{\bm{\mathbf{X}}}},{{\bm{\mathbf{Y}}}}\\}$, mini-batch
size $B$, # sghmc iterations $K$, # encoder interations $J$

3 Initialize the autoencoder parameters ${{\bm{\mathbf{\phi}}}}$ and
${{\bm{\mathbf{\varphi}}}}$

4 while  _${{\bm{\mathbf{\varphi}}}}$ , ${{\bm{\mathbf{Z}}}}$ have not
converged_ do

5 Sample a mini-batch of $B$ random indices $\mathcal{I}_{B}$

6 Sample random seed $\\{\bm{\varepsilon}_{i}\\}_{i=1}^{B}\sim
q({\bm{\mathbf{\varepsilon}}})$

7 Initialize the latent codes from the encoder
$\\{{{\bm{\mathbf{z}}}}_{i}\\}_{i=1}^{B}=f_{{{\bm{\mathbf{\phi}}}}}(\\{{{\bm{\mathbf{y}}}}_{i}\\}_{i\in\mathcal{I}_{B}},\\{\bm{\varepsilon}_{i}\\}_{i=1}^{B})$

8 for  _K iterations_ do

9 Compute energy function $\tilde{U}$ using Eq. 3

10 Sample from posterior
$p({{\bm{\mathbf{\varphi}}}},{{\bm{\mathbf{Z}}}}\,|\,{{\bm{\mathbf{Y}}}},{{\bm{\mathbf{X}}}})$:
${{\bm{\mathbf{\varphi}}}},\\{{{\bm{\mathbf{z}}}}\\}_{i=1}^{B}\leftarrow$
SGHMC(${{\bm{\mathbf{\varphi}}}},\\{{{\bm{\mathbf{z}}}}\\}_{i=1}^{B};\nabla_{{{\bm{\mathbf{\varphi}}}},{{\bm{\mathbf{z}}}}}\tilde{U}$)

11 for  _J iterations_ do

12 Compute objective function $\mathcal{L}$ using Eq. 4

13 Update encoder:
${{\bm{\mathbf{\phi}}}}\leftarrow\texttt{Optimizer}({{\bm{\mathbf{\phi}}}};\nabla_{{{\bm{\mathbf{\phi}}}}}\mathcal{L})$

14

15

Algorithm 1 Inference for baes with sghmc

##  4 Scalable Gaussian Process Prior for Bayesian Autoencoders

In the previous section, we introduced our novel version of a bae where we
imposed a simple fully-factorized prior over the latent space, such as
isotropic Gaussians. However, in many applications, such priors are incapable
of appropriately modeling the correlation nature of the data. For example, it
is sensible to model structured data evolving over time with a bae with a
prior over the latent space in the form of a gp with the auxiliary data as
input. In this section, we consider these scenarios precisely by introducing
gp priors in the latent space, which allows us to model sample covariances as
a function of the auxiliary data. We then discuss the scalability issues
induced by the use of gp priors, and we propose sgp-bae where we recover
scalability thanks to sparse approximations to gps (Quiñonero-Candela &
Rasmussen, 2005; Rossi et al., 2021). In this model, we carry out fully
Bayesian inference of the decoder, as well as the inducing inputs and
covariance hyperparameters of the sparse gps, while we optimize the stochastic
inference network implementing the encoder.

###  4.1 Gaussian process prior

We assume $C$ independent latent functions $f^{[1]},...,f^{[C]}$, which
results in each ${{\bm{\mathbf{z}}}}_{i}$ being evaluated at the corresponding
${{\bm{\mathbf{x}}}}_{i}$, i.e.,
${{\bm{\mathbf{z}}}}_{i}=\left[f^{[1]}({{\bm{\mathbf{x}}}}_{i}),\cdots,f^{[C]}({{\bm{\mathbf{x}}}}_{i})\right]$.
We assume that each function is drawn from a zero-mean gp prior with a
covariance function
$\kappa({{\bm{\mathbf{x}}}},{{\bm{\mathbf{x}}}}^{\prime};{{\bm{\mathbf{\theta}}}})$:

$\displaystyle
p({{\bm{\mathbf{Z}}}}\,|\,{{\bm{\mathbf{X}}}};{{\bm{\mathbf{\theta}}}})=\prod_{c=1}^{C}\mathcal{N}\left({{\bm{\mathbf{z}}}}^{[c]}_{1:N}\,|\,{\bm{\mathbf{0}}},{{\bm{\mathbf{K}}}}_{{{\bm{\mathbf{x}}}}{{\bm{\mathbf{x}}}}\,|\,{{\bm{\mathbf{\theta}}}}}\right),$
(5)

where the $c$-th latent channel of all latent variables,
${{\bm{\mathbf{z}}}}^{[c]}_{1:N}\in\mathbb{R}^{N}$ (the $c$-th column of
${{\bm{\mathbf{Z}}}}$), has a correlated Gaussian prior with covariance
${{\bm{\mathbf{K}}}}_{{{\bm{\mathbf{x}}}}{{\bm{\mathbf{x}}}}\,|\,{{\bm{\mathbf{\theta}}}}}\in\mathbb{R}^{N\times
N}$ obtained by evaluating
$\kappa({{\bm{\mathbf{x}}}}_{i},{{\bm{\mathbf{x}}}}_{j};{{\bm{\mathbf{\theta}}}})$
over all input pairs of ${{\bm{\mathbf{X}}}}$. Here, the latent function
values are informed by all ${{\bm{\mathbf{y}}}}$ values according to the
covariance of the corresponding auxiliary input ${{\bm{\mathbf{x}}}}$. One can
recover to the fully factorized
$\mathcal{N}({\bm{\mathbf{0}}},{{\bm{\mathbf{I}}}})$ prior on the latent space
by simply setting
${{\bm{\mathbf{K}}}}_{{{\bm{\mathbf{x}}}}{{\bm{\mathbf{x}}}}\,|\,{{\bm{\mathbf{\theta}}}}}={{\bm{\mathbf{I}}}}$.

This gp prior over the latent space of baes introduces fundamental scalability
issues. First, we have to compute the inverse and log-determinant of the
kernel matrix
${{\bm{\mathbf{K}}}}_{{{\bm{\mathbf{x}}}}{{\bm{\mathbf{x}}}}\,|\,{{\bm{\mathbf{\theta}}}}}$,
which results in $\mathcal{O}(N^{3})$ time complexity. This is only possible
when $N$ is of moderate size. Second, it is impossible to employ a mini-
batching inference method like sghmc since the energy function
$U({{\bm{\mathbf{\varphi}}}},{{\bm{\mathbf{Z}}}};{{\bm{\mathbf{X}}}},{{\bm{\mathbf{Y}}}})$
does not decompose as a sum over all the observations.

###  4.2 Bayesian sparse Gaussian processes

In order to keep the notation uncluttered, we focus on a single channel and
suppress the superscript index $c$. Given a set of latent function evaluations
over the dataset, ${{\bm{\mathbf{f}}}}=[f_{1},\cdots,f_{N}]^{\top}$, we assume
that the latent codes are stochastic realizations based on
${{\bm{\mathbf{f}}}}$ and additive Gaussian noise i.e.,
$\mathcal{N}({{\bm{\mathbf{Z}}}}\,|\,{{\bm{\mathbf{f}}}},\sigma^{2}{{\bm{\mathbf{I}}}})$.
Sparse gps (Quiñonero-Candela & Rasmussen, 2005) are a family of approximate
models that address the scalability problem by introducing a set of $M\ll N$
inducing points ${{\bm{\mathbf{u}}}}=(u_{1},\cdots,u_{M})$ at corresponding
inducing inputs
${{\bm{\mathbf{S}}}}=\\{{{\bm{\mathbf{s}}}}_{1},\cdots,{{\bm{\mathbf{s}}}}_{M}\\}$
such that $u_{i}=f({{\bm{\mathbf{s}}}}_{i})$. We assume that these inducing
variables follow the same gp as the original process, resulting in the
following joint prior:

$\displaystyle
p({{\bm{\mathbf{f}}}},{{\bm{\mathbf{u}}}})=\mathcal{N}\bigg{(}{\bm{\mathbf{0}}},\begin{bmatrix}{{\bm{\mathbf{K}}}}_{{{\bm{\mathbf{X}}}}{{\bm{\mathbf{X}}}}\,|\,{{\bm{\mathbf{\theta}}}}}&{{\bm{\mathbf{K}}}}_{{{\bm{\mathbf{X}}}}{{\bm{\mathbf{S}}}}\,|\,{{\bm{\mathbf{\theta}}}}}\\\
{{\bm{\mathbf{K}}}}_{{{\bm{\mathbf{S}}}}{{\bm{\mathbf{X}}}}\,|\,{{\bm{\mathbf{\theta}}}}}&{{\bm{\mathbf{K}}}}_{{{\bm{\mathbf{S}}}}{{\bm{\mathbf{S}}}}\,|\,{{\bm{\mathbf{\theta}}}}}\end{bmatrix}\bigg{)},$
(6)

where the covariance matrices
${{\bm{\mathbf{K}}}}_{{{\bm{\mathbf{S}}}}{{\bm{\mathbf{S}}}}\,|\,{{\bm{\mathbf{\theta}}}}}$
and
${{\bm{\mathbf{K}}}}_{{{\bm{\mathbf{X}}}}{{\bm{\mathbf{S}}}}\,|\,{{\bm{\mathbf{\theta}}}}}$
are computed between the elements in ${{\bm{\mathbf{S}}}}$ and
$\\{{{\bm{\mathbf{X}}}},{{\bm{\mathbf{S}}}}\\}$, respectively.

#### Fully Bayesian sparse gps.

The fully Bayesian treatment of sparse gps requires priors
$p_{\tau}({{\bm{\mathbf{\theta}}}})$ and $p_{\xi}({{\bm{\mathbf{S}}}})$ over
covariance hyper-parameters and inducing inputs, respectively, with $\tau$ and
$\xi$ as prior hyper-parameters. With these assumptions, we term this model as
Bayesian sparse Gaussian process autoencoder (sgp-bae), and the corresponding
generative model is illustrated in Fig. 1(b).

By defining
${{\bm{\mathbf{\Psi}}}}\stackrel{{\scriptstyle\text{def}}}{{=}}\\{{{\bm{\mathbf{\varphi}}}},{{\bm{\mathbf{u}}}},{{\bm{\mathbf{S}}}},{{\bm{\mathbf{\theta}}}}\\}$,
we can rewrite the full joint distribution of parameters in sgp-bae:

$\displaystyle
p({{\bm{\mathbf{\Psi}}}},{{\bm{\mathbf{f}}}},{{\bm{\mathbf{Z}}}},{{\bm{\mathbf{Y}}}}\,|\,{{\bm{\mathbf{X}}}})=$
(7)
$\displaystyle=\underbrace{p({{\bm{\mathbf{\Psi}}}})}_{\begin{subarray}{c}\text{Priors
on inducing}\\\ \text{inputs \& variables,
decoder}\end{subarray}}\underbrace{p({{\bm{\mathbf{f}}}}\,|\,{{\bm{\mathbf{u}}}},{{\bm{\mathbf{X}}}},{{\bm{\mathbf{S}}}},{{\bm{\mathbf{\theta}}}})p({{\bm{\mathbf{Z}}}}\,|\,{{\bm{\mathbf{f}}}};\sigma^{2}{{\bm{\mathbf{I}}}})}_{\text{Sparse
GP prior on latent
space}}\underbrace{p({{\bm{\mathbf{Y}}}}\,|\,{{\bm{\mathbf{Z}}}},{{\bm{\mathbf{\varphi}}}})}_{\begin{subarray}{c}\text{Likelihood
of}\\\ \text{observed data}\end{subarray}},$

where
$p({{\bm{\mathbf{\Psi}}}})=p({{\bm{\mathbf{\varphi}}}})p_{\tau}({{\bm{\mathbf{\theta}}}})p_{\xi}({{\bm{\mathbf{S}}}})p({{\bm{\mathbf{u}}}}\,|\,{{\bm{\mathbf{S}}}},{{\bm{\mathbf{\theta}}}})$.
Here,
$p({{\bm{\mathbf{u}}}}\,|\,{{\bm{\mathbf{S}}}},{{\bm{\mathbf{\theta}}}})=\mathcal{N}(\bm{0,{{\bm{\mathbf{K}}}}_{{{\bm{\mathbf{S}}}}{{\bm{\mathbf{S}}}}\,|\,{{\bm{\mathbf{\theta}}}}}})$,
and
$p({{\bm{\mathbf{f}}}}\,|\,{{\bm{\mathbf{u}}}},{{\bm{\mathbf{X}}}},{{\bm{\mathbf{S}}}},{{\bm{\mathbf{\theta}}}})=\mathcal{N}({{\bm{\mathbf{K}}}}_{{{\bm{\mathbf{X}}}}{{\bm{\mathbf{S}}}}\,|\,{{\bm{\mathbf{\theta}}}}}{{\bm{\mathbf{K}}}}_{{{\bm{\mathbf{S}}}}{{\bm{\mathbf{S}}}}\,|\,{{\bm{\mathbf{\theta}}}}}^{-1}{{\bm{\mathbf{u}}}},{{\bm{\mathbf{K}}}}_{{{\bm{\mathbf{X}}}}{{\bm{\mathbf{X}}}}\,|\,{{\bm{\mathbf{\theta}}}}}-{{\bm{\mathbf{K}}}}_{{{\bm{\mathbf{X}}}}{{\bm{\mathbf{S}}}}\,|\,{{\bm{\mathbf{\theta}}}}}{{\bm{\mathbf{K}}}}_{{{\bm{\mathbf{S}}}}{{\bm{\mathbf{S}}}}\,|\,{{\bm{\mathbf{\theta}}}}}^{-1}{{\bm{\mathbf{K}}}}_{{{\bm{\mathbf{S}}}}{{\bm{\mathbf{X}}}}\,|\,{{\bm{\mathbf{\theta}}}}})$.
We assume a factorization
$p({{\bm{\mathbf{Z}}}}\,|\,{{\bm{\mathbf{f}}}};\sigma^{2}{{\bm{\mathbf{I}}}})=\prod_{i=1}^{N}p({{\bm{\mathbf{z}}}}_{i}\,|\,f_{i};\sigma^{2})$
and make no further assumtions about the other distribution.

#### Scalable inference objective.

We wish to infer the set of variables ${{\bm{\mathbf{\Psi}}}}$ and the latent
codes ${{\bm{\mathbf{Z}}}}$. To do so, we have to marginalize out the latent
process ${{\bm{\mathbf{f}}}}$ from the full joint distribution above. In
particular, we have:

$\displaystyle\log
p({{\bm{\mathbf{\Psi}}}},{{\bm{\mathbf{Z}}}},{{\bm{\mathbf{Y}}}}\,|\,{{\bm{\mathbf{X}}}})=\log
p({{\bm{\mathbf{\Psi}}}})+$ (8) $\displaystyle+\log\int
p({{\bm{\mathbf{f}}}}\,|\,{{\bm{\mathbf{\Psi}}}},{{\bm{\mathbf{X}}}})p({{\bm{\mathbf{Z}}}}\,|\,{{\bm{\mathbf{f}}}},\sigma^{2}{{\bm{\mathbf{I}}}})d{{\bm{\mathbf{f}}}}+\log
p({{\bm{\mathbf{Y}}}}\,|\,{{\bm{\mathbf{Z}}}},{{\bm{\mathbf{\varphi}}}}).$

This objective should be decomposed over observations to sample from the
posterior over all the latent variables using a scalable method such as sghmc.
As discussed by Rossi et al. (2021), this can be done effectively by imposing
independence in the conditional distribution (Snelson & Ghahramani, 2005),
i.e., by parameterizing
$p({{\bm{\mathbf{f}}}}\,|\,{{\bm{\mathbf{\Psi}}}},{{\bm{\mathbf{X}}}})=\mathcal{N}\big{(}{{\bm{\mathbf{K}}}}_{{{\bm{\mathbf{x}}}}{{\bm{\mathbf{s}}}}\,|\,{{\bm{\mathbf{\theta}}}}}{{\bm{\mathbf{K}}}}_{{{\bm{\mathbf{s}}}}{{\bm{\mathbf{s}}}}\,|\,{{\bm{\mathbf{\theta}}}}}^{-1}{{\bm{\mathbf{u}}}},\textrm{diag}\big{[}{{\bm{\mathbf{K}}}}_{{{\bm{\mathbf{x}}}}{{\bm{\mathbf{x}}}}\,|\,{{\bm{\mathbf{\theta}}}}}-{{\bm{\mathbf{K}}}}_{{{\bm{\mathbf{x}}}}{{\bm{\mathbf{s}}}}\,|\,{{\bm{\mathbf{\theta}}}}}{{\bm{\mathbf{K}}}}_{{{\bm{\mathbf{s}}}}{{\bm{\mathbf{s}}}}\,|\,{{\bm{\mathbf{\theta}}}}}^{-1}{{\bm{\mathbf{K}}}}_{{{\bm{\mathbf{s}}}}{{\bm{\mathbf{x}}}}\,|\,{{\bm{\mathbf{\theta}}}}}\big{]}\big{)}$.
With this approximation, the log joint marginal becomes as follows:

$\displaystyle\log$ $\displaystyle
p({{\bm{\mathbf{\Psi}}}},{{\bm{\mathbf{Z}}}},{{\bm{\mathbf{Y}}}}\,|\,{{\bm{\mathbf{X}}}})\approx\log
p({{\bm{\mathbf{\Psi}}}})+$ (9) $\displaystyle\vspace{-3ex}+\sum_{n=1}^{N}$
$\displaystyle\bigg{\\{}\log\mathbb{E}_{p(f_{n}\,|\,{{\bm{\mathbf{\Psi}}}},{{\bm{\mathbf{X}}}})}[p({{\bm{\mathbf{z}}}}_{n}\,|\,f_{n};\sigma^{2})]+\log
p({{\bm{\mathbf{y}}}}_{n}\,|\,{{\bm{\mathbf{z}}}}_{n},{{\bm{\mathbf{\varphi}}}})\bigg{\\}}$
$\displaystyle\stackrel{{\scriptstyle\text{def}}}{{=}}-U({{\bm{\mathbf{\Psi}}}},{{\bm{\mathbf{Z}}}};{{\bm{\mathbf{X}}}},{{\bm{\mathbf{Y}}}}).$

We can now carry out inference for this sgp-bae model by plugging a mini-
batching approximation of this energy function into Algorithm 1 and Algorithm
1 of Algorithm 1. By using this sparse approximation, we reduce the
computational complexity of evaluating the gp prior down to
$\mathcal{O}(BM^{2})$, and sgp-bae can be readily applied to a generic dataset
and arbitrary gp kernel function.

#### Extension to deep Gaussian processes.

We can easily extend sgp-bae to deep Gaussian process priors (Damianou &
Lawrence, 2013) to model much more complex functions in the latent space of
baes. To the best of our knowledge, the use of deep gps has not been
considered in previous work. We assume a deep Gaussian process prior
$f^{(L)}\circ f^{(L-1)}\circ\cdots\circ f^{(1)}$, where each $f^{(l)}$ is a
gp. Each layer is associated with a set of kernel hyperparameters
${{\bm{\mathbf{\theta}}}}^{(l)}$, inducing inputs ${{\bm{\mathbf{S}}}}^{(l)}$
and inducing variables ${{\bm{\mathbf{u}}}}^{(l)}$.The set of variables to be
inferred is
${{\bm{\mathbf{\Psi}}}}=\\{{{\bm{\mathbf{\varphi}}}}\\}\cup\\{{{\bm{\mathbf{u}}}}^{(l)},{{\bm{\mathbf{S}}}}^{(l)},{{\bm{\mathbf{\theta}}}}^{(l)}\\}_{l=1}^{L}$.
The joint distribution is as follows:

$\displaystyle p({{\bm{\mathbf{\Psi}}}},$
$\displaystyle\\{{{\bm{\mathbf{f}}}}^{(l)}\\}_{l=1}^{L},{{\bm{\mathbf{Z}}}},{{\bm{\mathbf{Y}}}}\,|\,{{\bm{\mathbf{X}}}})=$
(10) $\displaystyle=p({{\bm{\mathbf{\Psi}}}})$
$\displaystyle\underbrace{\prod_{l=1}^{L}p({{\bm{\mathbf{f}}}}^{(l)}\,|\,{{\bm{\mathbf{f}}}}^{(l-1)},{{\bm{\mathbf{\Psi}}}})p({{\bm{\mathbf{Z}}}}\,|\,{{\bm{\mathbf{f}}}}^{(L)};\sigma^{2}{{\bm{\mathbf{I}}}})}_{\text{Deep
GP prior
}}p({{\bm{\mathbf{Y}}}}\,|\,{{\bm{\mathbf{Z}}}},{{\bm{\mathbf{\varphi}}}}),$

where we omit the dependency on ${{\bm{\mathbf{X}}}}$ for notational brevity.
To perform inference, the hidden layers ${{\bm{\mathbf{f}}}}^{(l)}$ have to be
marginalized and propagated up to the final layer $L$ (Salimbeni & Deisenroth,
2017). The marginalization can be approximated by quadrature (Hensman et al.,
2015) or through Monte Carlo sampling (Bonilla et al., 2019). Detailed
derivations of this extension can be found in Appendix C.

Table 2:  Reconstructions of the latent trajectories of moving ball. In the
first column, frames of each test video are overlayed and shaded by time.
Ground truth trajectories are illustrated in orange, while predicted
trajectories are depicted in blue. We use $M=10$ inducing points for the
methods employed with sparse gps.

gt video |  ( ) vae |  ( ) gpvae |  ( ) svgp-vae |  ( ) bae |  ( ) gp-bae |  ( ) sgp-bae
---|---|---|---|---|---|---
|  |  |  |  |  |

&

| VAE |  | GPVAE |  | SVGP-VAE
---|---|---|---|---|---
| BAE |  | GP-BAE |  |  SGP-BAE

Figure 2: Performance of autoencoder models as a function of the number of
inducing points.

#### Extension for missing data.

In practice, real-world data may be sparse, with many missing and few
overlapping dimensions across the entire dataset. We can easily extend sgp-bae
to handle such datasets. We assume that any possible permutation of observed
features is potentially missing, such that each high-dimensional observation
${{\bm{\mathbf{y}}}}_{n}={{\bm{\mathbf{y}}}}^{o}_{n}\cup{{\bm{\mathbf{y}}}}^{u}_{n}$
contains a set of observed features ${{\bm{\mathbf{y}}}}^{o}_{n}$ and
unobserved features ${{\bm{\mathbf{y}}}}^{u}_{n}$. The likelihood term of the
inference objective (§ 4.2) factorizes across data points and dimensions, so
there is no major modification in this objective, as the summation of the
likelihood term should be done only over the non-missing dimensions, i.e.:

$\displaystyle
p({{\bm{\mathbf{Y}}}}\,|\,{{\bm{\mathbf{Z}}}},{{\bm{\mathbf{\varphi}}}})=\sum_{n=1}^{N}\log
p({{\bm{\mathbf{y}}}}_{n}^{o}\,|\,{{\bm{\mathbf{z}}}}_{n},{{\bm{\mathbf{\varphi}}}}).$
(11)

##  5 Experiments

In this section, we provide empirical evidence that our sgp-bae outperforms
alternatives of combination between gp priors and ae models on synthetic data
and real-world high-dimensional benchmark datasets. Throughout all
experiments, unless otherwise specified, we use the radial basis function
(rbf) kernel with automatic relevance determination (ard) with marginal
variance and independent lengthscales $\lambda_{i}$ per feature (MacKay,
1996). We place a lognormal prior with unit variance and means equal to and
$0.05$ for the lengthscales and variance, respectively. Since the auxiliary
data of most of the considered datasets are timestamps, we impose a non-
informative uniform prior on the inducing inputs. We observe that this prior
works well in our experiments. We set the hyperparameters of the number of
sghmc and optimizer steps $J=30$, and $K=50$, respectively. The details for
all experiments are available in Appendix D.

###  5.1 Synthetic moving ball data

We begin our empirical evaluation by considering the moving ball dataset
proposed by Pearce (2019). This dataset comprises grayscale videos showing the
movement of a ball. The two-dimensional trajectory of the ball is simulated
from a gp characterized by an rbf kernel. Our task is to reconstruct the
underlying trajectory in the 2D latent space from the frames in pixel space.
Unlike Jazbec et al. (2021), we generate a fixed number of $35$ videos for
training and another $35$ videos for testing, rather than generating new
training videos at each iteration, regardless of the fact that tens of
thousands of iterations are performed. This new setting is more realistic and
is designed to show the data efficiency of the considered methods. It is still
possible to perform full gp inference on such a small dataset. Therefore, we
consider full GP-based methods, such as Gaussian Process Bayesian Autoencoder
(gp-bae) and gp-vae (Pearce, 2019), as oracles for the sparse variants.
Because the dataset is quite small, we perform full-batch training/inference.

#### Are Bayesian AEs better than variational counterparts?

In this experiment, we show that, by relaxing strong assumptions on the
posterior of latent space and taking advantage of a powerful scalable mcmc
method, baes consistently outperform vaes. Fig. 2 illustrates the performance
of the considered methods in terms of root mean squared error (rmse). The
results show that our gp-bae model performs much better than gp-vae (Pearce,
2019) though both models use the same full gp priors. In addition, by treating
inducing inputs and kernel hyper-parameters of sparse gps in a Bayesian
fashion, sgp-bae offers a rich modeling capability. This is evident in the
improved performance of sgp-bae compared to sparse variational Gaussian
process vae (svgp-vae) (Jazbec et al., 2021) and is able to reach closely to
the gp-bae performance despite using a small number of inducing points. These
numerical results align with the qualitative evaluation of the reconstructed
trajectories shown in Table 2. As expected, the standard bae and vae endowed
with a $\mathcal{N}({\bm{\mathbf{0}}},{{\bm{\mathbf{I}}}})$ prior on the
latent space completely fail to model the trajectories faithfully. In
contrast, gp-bae and sgp-bae are able to match them closely.

#### Benefits of being fully Bayesian.

Our sgp-bae model has the same advantage as the svgp-vae (Jazbec et al., 2021)
in that it allows for an arbitrary gp kernel, and the kernel hyperparameters
and inducing inputs can be inferred jointly during training or inference. In
contrast, other methods either use fixed gp priors (Pearce, 2019) or employ a
two-stage approach (Casale et al., 2018), where the gp hyperparameters are
optimized separately from the aes. svgp-vae optimizes these hyperparameters
using the common practice of maximization of the marginal likelihood, ml-ii.
This results in a point estimate of the hyperparameters but may be prone to
overfitting, especially when the training data is small while there are many
hyperparameters. A distinct advantage of sgp-bae over svgp-vae is that it is
fully Bayesian for the gp hyperparameters and inducing points. This not only
improves the quality of the predictions but also offers sensible uncertainty
quantification. Fig. 3 illustrates the posterior of the lengthscale. By using
a sufficient number of inducing points, sgp-bae is able to recover the true
lengthscale. When using too few inducing points, the model tends to estimate a
larger lengthscale. This is expected as the effective lengthscale of the
observed process in the subspace spanned by these few inducing points is
larger. We also observe that the more inducing points are used, the more
confident the posterior is over the lengthscale. Our method also produces a
sensible posterior distribution on the inducing positions, as shown in Fig. 4.
The estimated inducing positions are evenly spaced over the time dimension,
which is reasonable since the latent trajectories are generated from
stationary gps.

Figure 3: The posterior of the lengthscale corresponding to using a different
number of inducing points, $M$. The red line denotes the true lengthscale.
Figure 4: The posterior of the inducing position.

###  5.2 Conditional generation of rotated MNIST

In the next experiment, we consider a large-scale benchmark of conditional
generation. We follow the experimental setup from (Casale et al., 2018; Jazbec
et al., 2021), where they use a rotated mnist dataset ($N=4050$). İn
particular, we are given a set of images of digit three that have been rotated
at different angles in $[0,2\pi)$. Our goal is to generate a new image rotated
at an unseen angle. As it is not trivial to apply full gp for such a large
dataset, we omit gp-vae (Pearce, 2019) and gp-bae baselines. We consider the
gppvae model (Casale et al., 2018), which employs a low-rank approximation for
the gp, and the svgp-vae model (Jazbec et al., 2021) as baselines. For both
svgp-vae and sgp-bae, we use a number of inducing points of $M=32$ and a mini-
batch size of $B=256$. We also compare to Conditional Variational Autoencoder
(cvae) (Sohn et al., 2015), which allows conditional generation tasks.
Following Jazbec et al. (2021), we consider an extension of the sparse gp
model (Hensman et al., 2013), named deep-svigp, that utilizes a deep
likelihood parameterized by a neural network. As shown in Table 4, our sgp-bae
model generates images that are more visually appealing and most faithful to
the ground truth when compared to other approaches.

Table 3: Results on the rotated mnist digit 3 dataset. Here, we report the
mean and standard deviation computed from 4 runs.

| model |  mse ($\downarrow$)
---|---|---
( ) |  cvae (Sohn et al., 2015) |  0.0819 $\pm$ 0.0027
( ) |  gppvae (Casale et al., 2018) |  0.0351 $\pm$ 0.0005
( ) |  svgp-vae (Jazbec et al., 2021) |  0.0257 $\pm$ 0.0004
( ) | sgp-bae (ours) | 0.0228 $\pm$ 0.0004
|  deep-svigp (Jazbec et al., 2021) |  0.0236 $\pm$ 0.0010

Table 4: Conditionally generated mnist images. The right most column depicts
the epistemic uncertainty obtained by our sgp-bae model.

gt image |  deep-svigp |  ( ) cvae |  ( ) gppvae |  ( ) svgp-vae |  ( ) sgp-bae |  ( ) var.
---|---|---|---|---|---|---
|  |  |  |  |  |
|  |  |  |  |  |
|  |  |  |  |  |

| CVAE |  | SVGP-VAE
---|---|---|---
| DEEP-SVIGP |  | SGP-BAE (Ours)

Figure 5: Comparison of test mse on the rotated mnist dataset as function of
training time.

Table 5: A comparison between methods of multi-output gp models and gp
autoencoders on the eeg and jura datasets.

dataset | metric | igp | gpar |  ( ) sgp-vae |  ( ) sgp-bae (ours) | dsgp-bae (ours)
---|---|---|---|---|---|---
eeg |  smse ($\downarrow$) |  1.70 $\pm$ 0.14 |  0.28 $\pm$ 0.00 |  0.52 $\pm$ 0.05 | 0.22 $\pm$ 0.01 | 0.21 $\pm$ 0.01
nll ($\downarrow$) |  2.59 $\pm$ 0.02 | 1.68 $\pm$ 0.01 |  1.98 $\pm$ 0.02 |  1.96 $\pm$ 0.08 |  2.25 $\pm$ 0.13
jura |  mae ($\downarrow$) |  0.57 $\pm$ 0.00 | 0.42 $\pm$ 0.01 |  0.54 $\pm$ 0.01 | 0.45 $\pm$ 0.03 | 0.44 $\pm$ 0.02
nll ($\downarrow$) |  1563.42 $\pm$ 166.55 |  17.81 $\pm$ 1.06 |  1.02 $\pm$ 0.01 | 0.91 $\pm$ 0.04 | 0.85 $\pm$ 0.04

#### Performance of conditional generation.

Table 3 presents the quantitative evaluation of the conditionally generated
images in terms of mse. Our sgp-bae model clearly outperforms the other
competing methods. It is worth noting that deep-svigp (Hensman et al., 2013)
does not use an amortization mechanism, and its performance is considered to
be an upper bound for that of svgp-vae. As discussed by Jazbec et al. (2021),
deep-svigp can be used for conditional generation tasks, where the goal is to
impose a single gp over the entire dataset, and therefore amortization is not
necessary. However, this model cannot be used in tasks where inference has to
be amortized across multiple gps, such as learning interpretable data
representations.

#### Computational efficiency.

Similarly to the competing methods that use sparse approximations such as
svgp-vae and deep-svigp, each iteration of gp-bae involves the computation of
the inverse covariance matrix, resulting in a time complexity of
$\mathcal{O}(M^{3})$. Fig. 5 shows the convergence in terms of test mse for
the competing methods and our sgp-bae, trained for a fixed training time
budget. We omit gppvae (Casale et al., 2018) from this comparison, as reported
by Jazbec et al. (2021), which is significantly slower than the sparse
methods. This result demonstrates that sgp-bae converges remarkably faster in
terms of wall-clock time while achieving better final predictive performance.

#### Epistemic uncertainty quantification.

Unlike the competing methods, our sgp-bae model can capture epistemic
uncertainty by treating the decoder’s parameters in a Bayesian manner. This
capability can be helpful in out-of-distribution settings or when uncertainty
quantification on reconstructed or generated images is desired. As shown in
the rightmost column of Table 4, sgp-bae can provide sensible epistemic
uncertainty quantification. Our model exhibits increased uncertainty for
semantically and visually challenging pixels, such as the boundaries of the
digits.

###  5.3 Missing data imputation

Next, we consider the task of imputing missing data on multi-output spatio-
temporal datasets. We compare our method against Sparse gp-vae (sgp-vae)
(Ashman et al., 2020) and multi-output gp methods such as Independent gps
(igp) and Gaussian Process Autoregressive Regression (gpar) (Requeima et al.,
2019). Additionally, we consider the dsgp-bae model, which is an extension of
our sgp-bae model endowed with $3$-layer gp priors. For a fair comparison, we
treat partially missing data as zeros to feed into the inference network
(encoder) for sgp-vae and our sgp-bae model. We leave the adoption of partial
inference networks (Ashman et al., 2020) to our model for future work. We
follow the experimental setup of Requeima et al. (2019) and Ashman et al.
(2020) and use two standard datasets for this comparison.

#### Electroencephalogram (EEG).

This dataset comprises $256$ measurements taken over one second. Each
measurement consists of seven electrodes, FZ and F1-F6, placed on the
patient’s scalp (${{\bm{\mathbf{x}}}}_{n}\in\mathbb{R}^{1}$,
${{\bm{\mathbf{y}}}}_{n}\in\mathbb{R}^{7}$). The goal is to predict the last
$100$ samples for electrodes FZ, F1, and F2, given that the first $156$
samples of FZ, F1, and F2 and the whole signals of F3-F6 are observed.
Performance is measured with the standardized mean squared error (smse) and
negative log-likelihood (nll).

#### Jura.

This is a geospatial dataset consisting of $359$ measurements of the topsoil
concentrations of three metals — Nickel, Zinc, and Cadmium — collected from a
region of Swiss Jura (${{\bm{\mathbf{x}}}}_{n}\in\mathbb{R}^{2}$,
${{\bm{\mathbf{y}}}}_{n}\in\mathbb{R}^{3}$). The dataset is split into a
training set comprised of Nickel and Zinc measurements for all $359$ locations
and Cadmium measurements for just $259$ locations. The task is to predict the
Cadmium measurements at the remaining $100$ locations conditioned on the
observed training data. Performance is evaluated with the mean absolute error
(mae) and negative log-likelihood.

Table 5 compares the performance of sgp-bae to the competing methods. As
mentioned in Ashman et al. (2020), gpar is the state-of-the-art method for
these datasets. We find that our sgp-bae and dsgp-bae models perform
comparably with gpar on the eeg dataset but better on the jura dataset. A
significant advantage of gp autoencoder methods is that they model $P$ outputs
using only $C$ latent gps, while gpar uses $P$ gps. This can be beneficial
when the dimensionality of the data, $P$, is very high. Similarly to the
previous experiments, sgp-bae consistently performs better than variational
approximation-based methods such as sgp-vae (Ashman et al., 2020). As
expected, igp is the worst-performing method due to its inability to model the
correlations between output variables.

##  6 Conclusions

We have introduced our novel sgp-bae that integrates fully Bayesian sparse gps
on the latent space of Bayesian autoencoders. Our proposed model is generic,
as it allows an arbitrary gp kernel and deep gps. The inference for this model
is carried out by a powerful and scalable sghmc sampler. Through a rigorous
experimental campaign, we have demonstrated the excellent performance of sgp-
bae on a wide range of representation learning and generative modeling tasks.

#### Limitations and future works.

While our model’s ability to learn disentangled representations has been
demonstrated through empirical evidence, there is a need to establish a
theoretical guarantee for the disentanglement of its latent space.
Furthermore, it would be useful to study the amortization gap (Cremer et al.,
2018; Krishnan et al., 2018; Marino et al., 2018) of our model. Additionally,
exploring more informative priors for the decoder’s parameters (Tran et al.,
2021), beyond the isotropic Gaussian prior used in this work, is also
worthwhile. As a further direction for future research, it would be
interesting to apply the model to downstream applications where modeling
correlations between data points and uncertainty quantification are required,
such as in bioinformatics and climate modeling.

## Acknowledgements

We thank Simone Rossi and Dimitrios Milios for helpful discussions. SM
acknowledges support by the National Science Foundation (NSF) under an NSF
CAREER Award, award numbers 2003237 and 2007719, by the Department of Energy
under grant DE-SC0022331, by the HPI Research Center in Machine Learning and
Data Science at UC Irvine, and by gifts from Qualcomm and Disney. MF
gratefully acknowledges support from the AXA Research Fund and the Agence
Nationale de la Recherche (grant ANR-18-CE46-0002 and ANR-19-P3IA-0002).

## References

  * Ashman et al. (2020) Ashman, M., So, J., Tebbutt, W., Fortuin, V., Pearce, M., and Turner, R. E.  Sparse Gaussian Process Variational Autoencoders.  _arXiv preprint arXiv:2010.10177_ , 2020.
  * Bauer & Mnih (2019) Bauer, M. and Mnih, A.  Resampled Priors for Variational Autoencoders.  In _The 22nd International Conference on Artificial Intelligence and Statistics, AISTATS 2019, 16-18 April 2019, Naha, Okinawa, Japan_ , volume 89 of _Proceedings of Machine Learning Research_ , pp. 66–75. PMLR, 2019.
  * Bengio et al. (2013) Bengio, Y., Courville, A., and Vincent, P.  Representation Learning: A Review and New Perspectives.  _IEEE Transactions on Pattern Analysis and Machine Intelligence_ , 35(8):1798–1828, 2013.
  * Bishop (2006) Bishop, C. M.  _Pattern Recognition and Machine learning_ , volume 4.  Springer, 2006.
  * Bonilla et al. (2019) Bonilla, E. V., Krauth, K., and Dezfouli, A.  Generic Inference in Latent Gaussian Process Models.  _Journal of Machine Learning Research_ , 20(117):1–63, 2019.
  * Casale et al. (2018) Casale, F. P., Dalca, A., Saglietti, L., Listgarten, J., and Fusi, N.  Gaussian Process Prior Variational Autoencoders.  In _Advances in Neural Information Processing Systems_ , volume 31. Curran Associates, Inc., 2018.
  * Chen et al. (2014) Chen, T., Fox, E., and Guestrin, C.  Stochastic Gradient Hamiltonian Monte Carlo.  In _Proceedings of the 31st International Conference on Machine Learning, ICML 2014_ , Proceedings of Machine Learning Research, pp. 1683–1691, Bejing, China, 22–24 Jun 2014. PMLR.
  * Chen et al. (2017) Chen, X., Kingma, D. P., Salimans, T., Duan, Y., Dhariwal, P., Schulman, J., Sutskever, I., and Abbeel, P.  Variational Lossy Autoencoder.  In _5th International Conference on Learning Representations, ICLR 2017, Toulon, France, April 24-26, 2017, Conference Track Proceedings_. OpenReview.net, 2017.
  * Cottrell et al. (1989) Cottrell, G. W., Munro, P., and Zipser, D.  Image Compression by Back Propagation: A Demonstration of Extensional Programming.  _Models of Cognition_ , pp. 208–240, 1989.
  * Cremer et al. (2018) Cremer, C., Li, X., and Duvenaud, D.  Inference Suboptimality in Variational Autoencoders.  In _6th International Conference on Learning Representations, ICLR 2018, Vancouver, BC, Canada, April 30 - May 3, 2018, Conference Track Proceedings_. OpenReview.net, 2018.
  * Damianou & Lawrence (2013) Damianou, A. and Lawrence, N. D.  Deep Gaussian Processes.  In _Proceedings of the Sixteenth International Conference on Artificial Intelligence and Statistics_ , volume 31 of _Proceedings of Machine Learning Research_ , pp. 207–215, Scottsdale, Arizona, USA, 29 Apr–01 May 2013. PMLR.
  * Feng et al. (2017) Feng, Y., Wang, D., and Liu, Q.  Learning to Draw Samples with Amortized Stein Variational Gradient Descent.  In _Proceedings of the Thirty-Third Conference on Uncertainty in Artificial Intelligence, UAI 2017, Sydney, Australia, August 11-15, 2017_. AUAI Press, 2017.
  * Fortuin et al. (2020) Fortuin, V., Baranchuk, D., Raetsch, G., and Mandt, S.  GP-VAE: Deep Probabilistic Time Series Imputation.  In _Proceedings of the Twenty Third International Conference on Artificial Intelligence and Statistics_ , volume 108 of _Proceedings of Machine Learning Research_ , pp. 1651–1661. PMLR, 26–28 Aug 2020.
  * Garnelo et al. (2018) Garnelo, M., Schwarz, J., Rosenbaum, D., Viola, F., Rezende, D. J., Eslami, S., and Teh, Y. W.  Neural Processes.  _arXiv preprint arXiv:1807.01622_ , 2018.
  * Gelman & Rubin (1992) Gelman, A. and Rubin, D. B.  Inference from Iterative Simulation using Multiple Sequences.  _Statistical Science_ , 7(4):457–472, 1992.
  * Hensman et al. (2013) Hensman, J., Fusi, N., and Lawrence, N. D.  Gaussian Processes for Big Data.  In _Proceedings of the Twenty-Ninth Conference on Uncertainty in Artificial Intelligence_ , UAI’13, pp. 282–290, Arlington, Virginia, USA, 2013\. AUAI Press.
  * Hensman et al. (2015) Hensman, J., de G. Matthews, A. G., and Ghahramani, Z.  Scalable variational gaussian process classification.  In _Proceedings of the Eighteenth International Conference on Artificial Intelligence and Statistics, AISTATS 2015, San Diego, California, USA, May 9-12, 2015_ , volume 38 of _JMLR Workshop and Conference Proceedings_. JMLR.org, 2015.
  * Izmailov et al. (2021a) Izmailov, P., Vikram, S., Hoffman, M. D., and Wilson, A. G.  What Are Bayesian Neural Network Posteriors Really Like?  In _Proceedings of the 38th International Conference on Machine Learning, ICML 2021, 18-24 June 2021, Virtual Event_ , 2021a.
  * Izmailov et al. (2021b) Izmailov, P., Vikram, S., Hoffman, M. D., and Wilson, A. G. G.  What Are Bayesian Neural Network Posteriors Really Like?  In _Proceedings of the 38th International Conference on Machine Learning_ , volume 139 of _Proceedings of Machine Learning Research_ , pp. 4629–4640. PMLR, 18–24 Jul 2021b.
  * Jazbec et al. (2021) Jazbec, M., Ashman, M., Fortuin, V., Pearce, M., Mandt, S., and Rätsch, G.  Scalable Gaussian Process Variational Autoencoders .  In _Proceedings of The 24th International Conference on Artificial Intelligence and Statistics_ , volume 130 of _Proceedings of Machine Learning Research_ , pp. 3511–3519. PMLR, 13–15 Apr 2021.
  * Jordan et al. (1999) Jordan, M. I., Ghahramani, Z., Jaakkola, T. S., and Saul, L. K.  An Introduction to Variational Methods for Graphical Models.  _Machine learning_ , 37(2):183–233, 1999.
  * Kingma & Ba (2015) Kingma, D. P. and Ba, J.  Adam: A Method for Stochastic Optimization.  In _International Conference on Learning Representations_ , 2015.
  * Kingma & Welling (2014) Kingma, D. P. and Welling, M.  Auto-Encoding Variational Bayes.  In _International Conference on Learning Representations_ , 2014.
  * Korattikara Balan et al. (2015) Korattikara Balan, A., Rathod, V., Murphy, K. P., and Welling, M.  Bayesian Dark Knowledge.  In _Advances in Neural Information Processing Systems_ , volume 28. Curran Associates, Inc., 2015.
  * Krishnan et al. (2018) Krishnan, R., Liang, D., and Hoffman, M.  On the Challenges of Learning with Inference Networks on Sparse, High-dimensional Data.  In Storkey, A. and Perez-Cruz, F. (eds.), _Proceedings of the Twenty-First International Conference on Artificial Intelligence and Statistics_ , volume 84 of _Proceedings of Machine Learning Research_ , pp. 143–151. PMLR, 09–11 Apr 2018.
  * Lalchand et al. (2022) Lalchand, V., Bruinsma, W., Burt, D., and Rasmussen, C. E.  Sparse Gaussian Process Hyperparameters: Optimize or Integrate?  In _Advances in Neural Information Processing Systems_ , volume 35 (to appear), 2022.
  * Lawrence (2005) Lawrence, N. D.  Probabilistic Non-linear Principal Component Analysis with Gaussian Process Latent Variable Models.  _Journal of Machine Learning Research_ , 6:1783–1816, 2005\.
  * Li et al. (2017) Li, Y., Turner, R. E., and Liu, Q.  Approximate Inference with Amortised MCMC.  _arXiv preprint arXiv:1702.08343_ , 2017.
  * MacKay (1996) MacKay, D. J.  Bayesian methods for backpropagation networks.  In _Models of neural networks III_ , pp. 211–254. Springer, 1996\.
  * Mackay (1992) Mackay, D. J. C.  _Bayesian Methods for Adaptive Models_.  PhD thesis, California Institute of Technology, USA, 1992.  UMI Order No. GAX92-32200.
  * Marino et al. (2018) Marino, J., Yue, Y., and Mandt, S.  Iterative Amortized Inference.  In _International Conference on Machine Learning_ , pp. 3403–3412. PMLR, 2018.
  * Nalisnick & Smyth (2017) Nalisnick, E. T. and Smyth, P.  Stick-Breaking Variational Autoencoders.  In _5th International Conference on Learning Representations, ICLR 2017, Toulon, France, April 24-26, 2017, Conference Track Proceedings_. OpenReview.net, 2017.
  * Neal (1996) Neal, R. M.  _Bayesian Learning for Neural Networks (Lecture Notes in Statistics)_.  Springer, 1 edition, August 1996.
  * Neal (2011) Neal, R. M.  _MCMC Using Hamiltonian Dynamics_ , chapter 5.  CRC Press, 2011.  doi: 10.1201/b10905-7.
  * Pearce (2019) Pearce, M.  The Gaussian Process Prior VAE for Interpretable Latent Dynamics from Pixels.  In _Symposium on Advances in Approximate Bayesian Inference, AABI 2019, Vancouver, BC, Canada, December 8, 2019_ , volume 118 of _Proceedings of Machine Learning Research_ , pp. 1–12. PMLR, 2019.
  * Quiñonero-Candela & Rasmussen (2005) Quiñonero-Candela, J. and Rasmussen, C. E.  A Unifying View of Sparse Approximate Gaussian Process Regression.  _Journal of Machine Learning Research_ , 6(65):1939–1959, 2005.
  * Ramchandran et al. (2021) Ramchandran, S., Tikhonov, G., Kujanpää, K., Koskinen, M., and Lähdesmäki, H.  Longitudinal Variational Autoencoder .  In _Proceedings of The 24th International Conference on Artificial Intelligence and Statistics_ , volume 130 of _Proceedings of Machine Learning Research_ , pp. 3898–3906. PMLR, 13–15 Apr 2021.
  * Rasmussen & Williams (2006) Rasmussen, C. E. and Williams, C.  _Gaussian Processes for Machine Learning_.  MIT Press, 2006.
  * Requeima et al. (2019) Requeima, J., Tebbutt, W., Bruinsma, W., and Turner, R. E.  The Gaussian Process Autoregressive Regression Model (GPAR).  In _Proceedings of the Twenty-Second International Conference on Artificial Intelligence and Statistics_ , volume 89 of _Proceedings of Machine Learning Research_ , pp. 1860–1869. PMLR, 16–18 Apr 2019.
  * Rezende et al. (2014a) Rezende, D. J., Mohamed, S., and Wierstra, D.  Stochastic backpropagation and approximate inference in deep generative models.  In _Proceedings of the 31th International Conference on Machine Learning, ICML 2014, Beijing, China, 21-26 June 2014_ , volume 32 of _JMLR Workshop and Conference Proceedings_ , pp. 1278–1286. JMLR.org, 2014a.
  * Rezende et al. (2014b) Rezende, D. J., Mohamed, S., and Wierstra, D.  Stochastic Backpropagation and Approximate Inference in Deep Generative Models.  In _Proceedings of the 31th International Conference on Machine Learning, ICML 2014_ , volume 32 of _Proceeding of Machine Learning Research_ , pp. 1278–1286, Beijing, China, 21-26 June 2014b. PMLR.
  * Rossi et al. (2021) Rossi, S., Heinonen, M., Bonilla, E., Shen, Z., and Filippone, M.  Sparse Gaussian Processes Revisited: Bayesian Approaches to Inducing-Variable Approximations .  In _Proceedings of The 24th International Conference on Artificial Intelligence and Statistics_ , volume 130 of _Proceedings of Machine Learning Research_ , pp. 1837–1845. PMLR, 13–15 Apr 2021.
  * Salimbeni & Deisenroth (2017) Salimbeni, H. and Deisenroth, M.  Doubly Stochastic Variational Inference for Deep Gaussian Processes.  In _Advances in Neural Information Processing Systems_ , volume 30. Curran Associates, Inc., 2017.
  * Shi et al. (2019) Shi, J., Khan, M. E., and Zhu, J.  Scalable Training of Inference Networks for Gaussian-Process Models.  In _Proceedings of the 36th International Conference on Machine Learning_ , volume 97 of _Proceedings of Machine Learning Research_ , pp. 5758–5768. PMLR, 09–15 Jun 2019.
  * Snelson & Ghahramani (2005) Snelson, E. and Ghahramani, Z.  Sparse Gaussian Processes using Pseudo-inputs.  In _Advances in Neural Information Processing Systems_ , volume 18. MIT Press, 2005.
  * Sohn et al. (2015) Sohn, K., Lee, H., and Yan, X.  Learning Structured Output Representation using Deep Conditional Generative Models.  In _Advances in Neural Information Processing Systems_ , volume 28. Curran Associates, Inc., 2015.
  * Springenberg et al. (2016) Springenberg, J. T., Klein, A., Falkner, S., and Hutter, F.  Bayesian Optimization with Robust Bayesian Neural Networks.  In _Advances in Neural Information Processing Systems_ , volume 29, pp. 4134–4142. Curran Associates, Inc., 2016.
  * Titsias (2009) Titsias, M.  Variational Learning of Inducing Variables in Sparse Gaussian Processes.  In _Proceedings of the Twelth International Conference on Artificial Intelligence and Statistics_ , volume 5 of _Proceedings of Machine Learning Research_ , pp. 567–574, Hilton Clearwater Beach Resort, Clearwater Beach, Florida USA, 16–18 Apr 2009. PMLR.
  * Tomczak (2022) Tomczak, J. M.  Deep Generative Modeling.  In _Deep Generative Modeling_ , pp. 1–12. Springer, 2022.
  * Tomczak & Welling (2018) Tomczak, J. M. and Welling, M.  VAE with a VampPrior.  In _International Conference on Artificial Intelligence and Statistics, AISTATS 2018, 9-11 April 2018, Playa Blanca, Lanzarote, Canary Islands, Spain_ , volume 84 of _Proceedings of Machine Learning Research_ , pp. 1214–1223. PMLR, 2018.
  * Tran et al. (2021) Tran, B.-H., Rossi, S., Milios, D., Michiardi, P., Bonilla, E. V., and Filippone, M.  Model Selection for Bayesian Autoencoders.  In _Advances in Neural Information Processing Systems_ , volume 34, pp. 19730–19742. Curran Associates, Inc., 2021.
  * Tran et al. (2022) Tran, B.-H., Rossi, S., Milios, D., and Filippone, M.  All You Need is a Good Functional Prior for Bayesian Deep Learning.  _Journal of Machine Learning Research_ , 23(74):1–56, 2022.
  * van der Maaten & Hinton (2008) van der Maaten, L. and Hinton, G.  Visualizing Data using t-SNE.  _Journal of Machine Learning Research_ , 9(86):2579–2605, 2008.
  * Wang & Liu (2016) Wang, D. and Liu, Q.  Learning to Draw Samples: With Application to Amortized MLE for Generative Adversarial Learning.  _arXiv preprint arXiv:1611.01722_ , 2016.
  * Wang et al. (2018) Wang, K.-C., Vicol, P., Lucas, J., Gu, L., Grosse, R., and Zemel, R.  Adversarial Distillation of Bayesian Neural Network Posteriors.  In _Proceedings of the 35th International Conference on Machine Learning_ , volume 80 of _Proceedings of Machine Learning Research_ , pp. 5190–5199. PMLR, 10–15 Jul 2018.
  * Wilson & Izmailov (2020) Wilson, A. G. and Izmailov, P.  Bayesian Deep Learning and a Probabilistic Perspective of Generalization.  In _Advances in Neural Information Processing Systems_ , volume 33, pp. 4697–4708. Curran Associates, Inc., 2020.
  * Yang et al. (2022) Yang, Y., Mandt, S., and Theis, L.  An Introduction to Neural Data Compression.  _arXiv preprint arXiv:2202.06533_ , 2022.
  * Zhang et al. (2018) Zhang, C., Bütepage, J., Kjellström, H., and Mandt, S.  Advances in variational inference.  _IEEE transactions on pattern analysis and machine intelligence_ , 41(8):2008–2026, 2018.
  * Zhang et al. (2020) Zhang, R., Li, C., Zhang, J., Chen, C., and Wilson, A. G.  Cyclical Stochastic Gradient MCMC for Bayesian Deep Learning.  In _8th International Conference on Learning Representations, ICLR 2020, Addis Ababa, Ethiopia, April 26-30, 2020_. OpenReview.net, 2020.
  * Zhu et al. (2022) Zhu, H., Rodas, C. B., and Li, Y.  Markovian Gaussian Process Variational Autoencoders.  _arXiv preprint arXiv:2207.05543_ , 2022.

##  Appendix A A taxonomy of latent variable models

By considering the characteristics of the prior distribution on latent
variables, the likelihood function, input dependencies, and Bayesian
treatments, we can draw connections between our proposed models and other
latent variable models. Fig. 6 summarizes these relationships. Here, we assume
an isotropic Gaussian likelihood with precision $\beta$ for the high-
dimensional observed data ${{\bm{\mathbf{y}}}}_{n}$ as used in our
experiments.

Probabilistic principal component analysis (pca) (Bishop, 2006) is a simple
latent variable model that imposes an isotropic Gaussian prior over the latent
space and linear projection from the latent variables to the observed data.
Variational Autoencoders (vaes) (Kingma & Welling, 2014; Rezende et al.,
2014a) build upon this by introducing a non-linear parametric mapping to the
observed data, while Gaussian process latent variable models (gplvms) utilize
a non-parametric Gaussian Process (gp) mapping. Conditional Variational
Autoencoder (cvae) (Sohn et al., 2015) is an extension of vaes that utilizes
an input-dependent prior in the latent space for conditional generation tasks.
However, this model does not account for correlations between latent
representations. Gaussian Process vaes (Casale et al., 2018; Pearce, 2019;
Jazbec et al., 2021) overcome this problem by imposing gp priors on the latent
space. Our model, Sparse Gaussian Process Bayesian Autoencoder (sgp-bae),
further enhances the modeling capabilities of these models by using a fully
Bayesian approach for the latent variables, decoder, and gp hyper-parameters
in a fully Bayesian manner, and decoupling the model from the variational
inference formulation.

Latent neural processes (Garnelo et al., 2018) can be seen as an extension of
cvae. However, this model follows a meta-learning approach and splits the data
into a context set,
$\left\\{{{\bm{\mathbf{x}}}}_{n}^{[C]},{{\bm{\mathbf{y}}}}_{n}^{[C]}\right\\}_{n=1}^{N_{C}}$,
and a target set,
$\left\\{{{\bm{\mathbf{x}}}}_{n}^{[T]},{{\bm{\mathbf{y}}}}_{n}^{[T]}\right\\}_{n=1}^{N_{T}}$.
This model uses a global latent variable ${{\bm{\mathbf{z}}}}$ to capture the
global properties of the context data. The likelihood is conditioned on both
new target input ${{\bm{\mathbf{x}}}}_{n}^{[T]}$ and the global latent
variable ${{\bm{\mathbf{z}}}}$.

Figure 6: Connections between our proposed models and other latent variables
models. References are [a] for Bishop (2006), [b] for Lawrence (2005), [c] for
Damianou & Lawrence (2013), [d] for Kingma & Welling (2014); Rezende et al.
(2014a), [e] for Sohn et al. (2015), [f] for Casale et al. (2018); Pearce
(2019); Jazbec et al. (2021), and [g] for Garnelo et al. (2018).

##  Appendix B Details on stochastic gradient Hamiltonian Monte Carlo

Hamiltonian Monte Carlo (hmc) (Neal, 2011) is a highly-efficient Markov Chain
Monte Carlo (MCMC) method used to generate samples ${{\bm{\mathbf{\theta}}}}$
from a potential energy function $U({{\bm{\mathbf{\theta}}}})$. hmc introduces
auxiliary momentum variables ${{\bm{\mathbf{r}}}}$ and samples are then
generated from the joint distribution
$p({{\bm{\mathbf{\theta}}}},{{\bm{\mathbf{r}}}})$ using the Hamiltonian
dynamics:

$\displaystyle\begin{cases}d{{\bm{\mathbf{\theta}}}}&={{\bm{\mathbf{M}}}}^{-1}{{\bm{\mathbf{r}}}}dt,\\\
d{{\bm{\mathbf{r}}}}&=-\nabla U({{\bm{\mathbf{\theta}}}})dt,\end{cases}$ (12)

where ${{\bm{\mathbf{M}}}}$ is an arbitrary mass matrix that serves as a
preconditioner. The continuous system described by hmc is approximated using
$\eta$-discretized numerical integration, and subsequent Metropolis steps are
applied to account for errors that may result from this approximation.

However, hmc becomes computationally inefficient when applied to large
datasets, as it requires the calculation of the gradient $\nabla
U({{\bm{\mathbf{\theta}}}})$ on the entire dataset. To address this issue,
(Chen et al., 2014) proposed an extension of hmc called stochastic gradient
Hamiltonian Monte Carlo (sghmc), which uses a noisy but unbiased estimate of
the gradient $\nabla\tilde{U}({{\bm{\mathbf{\theta}}}})$ computed from a mini-
batch of the data. The discretized Hamiltonian dynamics are then updated using
this estimate of the gradient as follows:

$\displaystyle\begin{cases}\Delta{{\bm{\mathbf{\theta}}}}&={{\bm{\mathbf{M}}}}^{-1}{{\bm{\mathbf{r}}}},\\\
\Delta{{\bm{\mathbf{r}}}}&=-\eta\nabla\tilde{U}({{\bm{\mathbf{\theta}}}})-\eta{{\bm{\mathbf{C}}}}{{\bm{\mathbf{M}}}}^{-1}{{\bm{\mathbf{r}}}}+\mathcal{N}(0,2\eta({{\bm{\mathbf{C}}}}-\tilde{{{\bm{\mathbf{B}}}}})),\end{cases}$
(13)

where $\eta$ is the step size, ${{\bm{\mathbf{C}}}}$ is a user-defined
friction matrix, $\tilde{{{\bm{\mathbf{B}}}}}$ is the estimate for the noise
of the gradient evaluation. To select suitable values for these
hyperparameters, we use a scale-adapted version of sghmc proposed by
Springenberg et al. (2016), in which the hyperparameters are automatically
adjusted during a burn-in phase. After this period, the values of the
hyperparameters are fixed.

#### Estimating ${{\bm{\mathbf{M}}}}$.

In our implementation of sghmc, we set the mass matrix
${{\bm{\mathbf{M}}}}^{-1}=\mathrm{diag}\left(\hat{V}_{{{\bm{\mathbf{\theta}}}}}^{-1/2}\right)$,
where $\hat{V}_{{{\bm{\mathbf{\theta}}}}}$ is an estimate of the uncentered
variance of the gradient,
$\hat{V}_{{{\bm{\mathbf{\theta}}}}}\approx\mathbb{E}[(\nabla\tilde{U}({{\bm{\mathbf{\theta}}}}))^{2}]$.
which can be estimated by using exponential moving average as follows:

$\Delta\hat{V}_{{{\bm{\mathbf{\theta}}}}}=-\tau^{-1}\hat{V}_{{{\bm{\mathbf{\theta}}}}}+\tau^{-1}\nabla(\tilde{U}({{\bm{\mathbf{\theta}}}}))^{2},$
(14)

where $\tau$ is a parameter vector that specifies the moving average windows.
This parameter can be chosen adaptively using the method proposed by
Springenberg et al. (2016) as follows:

$\Delta\tau=-g_{{{\bm{\mathbf{\theta}}}}}^{2}\hat{V}^{-1}_{{{\bm{\mathbf{\theta}}}}}\tau+1,\quad\text{and},\quad\Delta
g_{{{\bm{\mathbf{\theta}}}}}=-\tau^{-1}g_{{{\bm{\mathbf{\theta}}}}}+\tau^{-1}\nabla\tilde{U}({{\bm{\mathbf{\theta}}}}),$
(15)

where $g_{{{\bm{\mathbf{\theta}}}}}$ is a smoothed estimate of the gradient
$\nabla U({{\bm{\mathbf{\theta}}}})$.

#### Estimating $\tilde{{{\bm{\mathbf{B}}}}}$.

Ideally, the estimate of the noise of the gradient evaluation,
$\tilde{{{\bm{\mathbf{B}}}}}$, should be the empirical Fisher information
matrix of $U({{\bm{\mathbf{\theta}}}})$ However, this matrix is
computationally expensive to calculate, so we instead use a diagonal
approximation given by
$\tilde{{{\bm{\mathbf{B}}}}}=\frac{1}{2}\eta\hat{V}_{{{\bm{\mathbf{\theta}}}}}$.
This approximation is already obtained from the step of estimating
${{\bm{\mathbf{M}}}}$.

#### Choosing ${{\bm{\mathbf{C}}}}$.

In practice, it is common to set the friction matrix ${{\bm{\mathbf{C}}}}$ to
be equal to the identity matrix, i.e.
${{\bm{\mathbf{C}}}}=C{{\bm{\mathbf{I}}}}$, which means that the same
independent noise is applied to each element of ${{\bm{\mathbf{\theta}}}}$.

#### The discretized Hamiltonian dynamics.

By substituting
${{\bm{\mathbf{v}}}}:=\eta\hat{V}_{{{\bm{\mathbf{\theta}}}}}^{-1/2}{{\bm{\mathbf{r}}}}$,
the dynamics Eq. 13 become

$\displaystyle\begin{cases}\Delta{{\bm{\mathbf{\theta}}}}&={{\bm{\mathbf{v}}}},\\\
\Delta{{\bm{\mathbf{v}}}}&=-\eta^{2}\hat{V}_{{{\bm{\mathbf{\theta}}}}}^{-1/2}\nabla\tilde{U}({{\bm{\mathbf{\theta}}}})-\eta
C\hat{V}_{{{\bm{\mathbf{\theta}}}}}^{-1/2}{{\bm{\mathbf{v}}}}+\mathcal{N}(0,2\eta^{3}C\hat{V}_{{{\bm{\mathbf{\theta}}}}}^{-1}-\eta^{4}{{\bm{\mathbf{I}}}}).\end{cases}$
(16)

Following (Springenberg et al., 2016), we choose $C$ such that $\eta
C\hat{V}_{{{\bm{\mathbf{\theta}}}}}^{-1/2}=\alpha{{\bm{\mathbf{I}}}}$. This is
equivalent to using a constant momentum coefficient of $\alpha$. The final
discretized dynamics are then

$\displaystyle\begin{cases}\Delta{{\bm{\mathbf{\theta}}}}&={{\bm{\mathbf{v}}}},\\\
\Delta{{\bm{\mathbf{v}}}}&=-\eta^{2}\hat{V}_{{{\bm{\mathbf{\theta}}}}}^{-1/2}\nabla\tilde{U}({{\bm{\mathbf{\theta}}}})-\alpha{{\bm{\mathbf{v}}}}+\mathcal{N}(0,2\eta^{2}\alpha\hat{V}_{{{\bm{\mathbf{\theta}}}}}^{-1/2}-\eta^{4}{{\bm{\mathbf{I}}}}).\end{cases}$
(17)

##  Appendix C Details of the scalable sampling objective for sparse Gaussian
processes

gps (Rasmussen & Williams, 2006) are one of the main workhorses of Bayesian
nonparametric statistics and machine learning, as they provide a flexible and
powerful tool for imposing a prior distribution over functions:

$\displaystyle
f({{\bm{\mathbf{x}}}})\sim\mathcal{GP}(m({{\bm{\mathbf{x}}}}),\kappa({{\bm{\mathbf{x}}}},{{\bm{\mathbf{x}}}}^{\prime};{{\bm{\mathbf{\theta}}}})),$
(18)

where $f({{\bm{\mathbf{x}}}}):\mathbb{R}^{D}\rightarrow\mathbb{R}$ maps
$D$-dimensional inputs into one-dimensional outputs. gps are fully
characterized by their mean and their covariance:

$\displaystyle\mathbb{E}[f({{\bm{\mathbf{x}}}})]=m({{\bm{\mathbf{x}}}}),\quad\quad\text{cov}[f({{\bm{\mathbf{x}}}}),f({{\bm{\mathbf{x}}}}^{\prime})]=\kappa({{\bm{\mathbf{x}}}},{{\bm{\mathbf{x}}}}^{\prime};{{\bm{\mathbf{\theta}}}}),$
(19)

where $m:\mathbb{R}^{D}\rightarrow\mathbb{R}$ is the mean and
$k:\mathbb{R}^{D}\times\mathbb{R}^{D}\rightarrow\mathbb{R}$ is the kernel
function with hyperparameters ${{\bm{\mathbf{\theta}}}}$. gps indeed can be
viewed as an extension of a multivariate Gaussian distribution to infinitely
many dimensions, where any fixed set of inputs
${{\bm{\mathbf{X}}}}\in\mathbb{R}^{N\times D}$ follows the Gaussian
distribution:

$\displaystyle
p({{\bm{\mathbf{f}}}})=\mathcal{N}({\bm{\mathbf{0}}},{{\bm{\mathbf{K}}}}_{{{\bm{\mathbf{x}}}}{{\bm{\mathbf{x}}}}\,|\,{{\bm{\mathbf{\theta}}}}}).$
(20)

Here, we assume a zero-mean gp prior; ${{\bm{\mathbf{f}}}}\in\mathbb{R}^{N}$
are the realizations gp random variables at the inputs ${{\bm{\mathbf{X}}}}$;
and
${{\bm{\mathbf{K}}}}_{{{\bm{\mathbf{X}}}}{{\bm{\mathbf{X}}}}\,|\,{{\bm{\mathbf{\theta}}}}}$
is the covariance matrix obtained by evaluating
$\kappa({{\bm{\mathbf{x}}}},{{\bm{\mathbf{x}}}}^{\prime};{{\bm{\mathbf{\theta}}}})$
over all input pairs ${{\bm{\mathbf{x}}}}_{i}$, ${{\bm{\mathbf{x}}}}_{j}$ (we
will drop the explicit parameterization on ${{\bm{\mathbf{\theta}}}}$ for the
sake of notation). From now on, in order to keep the notation uncluttered, we
focus on a single channel of latent space of autoencoders (aes). We assume
that the latent codes ${{\bm{\mathbf{Z}}}}$ are stochastic realizations based
on ${{\bm{\mathbf{f}}}}$ and additive Gaussian noise i.e.
${{\bm{\mathbf{Z}}}}\,|\,{{\bm{\mathbf{f}}}}\sim\mathcal{N}({{\bm{\mathbf{f}}}},\sigma^{2}{{\bm{\mathbf{I}}}})$.
We further assume a full factorization
$p({{\bm{\mathbf{Z}}}}\,|\,{{\bm{\mathbf{f}}}};\sigma^{2}{{\bm{\mathbf{I}}}})=\prod_{n=1}^{N}p({{\bm{\mathbf{z}}}}_{n}\,|\,f_{n};\sigma^{2})$.

Though gps provide an elegant mechanism to handle uncertainties, their
computational complexity grows cubically with the number of inputs, i.e.
$\mathcal{O}(N^{3})$. This problem is commonly tackled by sparse gps, which
are a family of approximate models that address the scalability issue by
introducing a set of $M$ inducing variables
${{\bm{\mathbf{u}}}}=(u_{1},...,u_{M})^{\top}\in\mathbb{R}^{M\times 1}$ at the
corresponding inducing inputs
${{\bm{\mathbf{S}}}}=[{{\bm{\mathbf{s}}}}_{1}^{\top},...,{{\bm{\mathbf{s}}}}_{M}^{\top}]^{\top}\in\mathbb{R}^{M\times
D}$ with $u_{m}\equiv f({{\bm{\mathbf{s}}}}_{m})$. The inducing points
${{\bm{\mathbf{S}}}}$ can be interpreted as a compressed version of the
training data where the number of inducing points $M$ acts as a trade-off
parameter between the goodness of the approximation and scalability. The
inducing variables are assumed to be drawn from the same gp as the original
process, i.e.:

$\displaystyle p({{\bm{\mathbf{f}}}},{{\bm{\mathbf{u}}}})$
$\displaystyle=p({{\bm{\mathbf{u}}}})p({{\bm{\mathbf{f}}}}\,|\,{{\bm{\mathbf{u}}}}),\text{with}$
(21) $\displaystyle p({{\bm{\mathbf{u}}}})$
$\displaystyle=\mathcal{N}({\bm{\mathbf{0}}},{{\bm{\mathbf{K}}}}_{{{\bm{\mathbf{s}}}}{{\bm{\mathbf{s}}}}}),$
(22) $\displaystyle p({{\bm{\mathbf{f}}}}\,|\,{{\bm{\mathbf{u}}}})$
$\displaystyle=\mathcal{N}({{\bm{\mathbf{K}}}}_{{{\bm{\mathbf{x}}}}{{\bm{\mathbf{s}}}}}{{\bm{\mathbf{K}}}}_{{{\bm{\mathbf{s}}}}{{\bm{\mathbf{s}}}}}^{-1}{{\bm{\mathbf{u}}}},{{\bm{\mathbf{K}}}}_{{{\bm{\mathbf{x}}}}{{\bm{\mathbf{x}}}}}-{{\bm{\mathbf{K}}}}_{{{\bm{\mathbf{x}}}}{{\bm{\mathbf{s}}}}}{{\bm{\mathbf{K}}}}_{{{\bm{\mathbf{s}}}}{{\bm{\mathbf{s}}}}}^{-1}{{\bm{\mathbf{K}}}}_{{{\bm{\mathbf{s}}}}{{\bm{\mathbf{x}}}}}),$
(23)

where the covariance matrices
${{\bm{\mathbf{K}}}}_{{{\bm{\mathbf{s}}}}{{\bm{\mathbf{s}}}}},{{\bm{\mathbf{K}}}}_{{{\bm{\mathbf{x}}}}{{\bm{\mathbf{s}}}}}$
are computed between the elements in ${{\bm{\mathbf{S}}}}$ and
$\\{{{\bm{\mathbf{X}}}},{{\bm{\mathbf{S}}}}\\}$, respectively, and
${{\bm{\mathbf{K}}}}_{{{\bm{\mathbf{s}}}}{{\bm{\mathbf{x}}}}}={{\bm{\mathbf{K}}}}_{{{\bm{\mathbf{s}}}}{{\bm{\mathbf{x}}}}}^{\top}$.

There is a line of works using variational techniques for sparse gps priors
for vaes (Jazbec et al., 2021; Ashman et al., 2020). More specifically, they
rely on the variational framework proposed by Titsias (2009); Hensman et al.
(2013), enabling the use of gp priors on very large datasets. However, the
posterior of the inducing variables ${{\bm{\mathbf{u}}}}$ under this framework
involves constraining to have a parametric form (usually a Gaussian).

In this work, we take a different route as we treat the inducing variables
${{\bm{\mathbf{u}}}}$, inducing inputs ${{\bm{\mathbf{S}}}}$ as well as the
kernel hyperparameters ${{\bm{\mathbf{\theta}}}}$ in a fully Bayesian way. As
discussed in the main paper, we wish to infer these variables together with
the decoders’ parameters and the latent codes by using a powerful and scalable
sghmc (Chen et al., 2014) sampler. To do so, the sampling objective § 4.2
should be decomposed over all data points. The main obstacle is the joint
distribution
$p({{\bm{\mathbf{Z}}}},{{\bm{\mathbf{u}}}}\,|\,{{\bm{\mathbf{\theta}}}})=\mathbb{E}_{p({{\bm{\mathbf{f}}}}\,|\,{{\bm{\mathbf{u}}}})}[p({{\bm{\mathbf{Z}}}}\,|\,{{\bm{\mathbf{f}}}};\sigma^{2}{{\bm{\mathbf{I}}}})]$,
which has a complicated form due to the expectation under the conditional
$p({{\bm{\mathbf{f}}}}\,|\,{{\bm{\mathbf{u}}}})$. As discussed by Rossi et al.
(2021), this problem can be solved effectively by the fully independent
training conditionals (fitc; see Quiñonero-Candela & Rasmussen, 2005), i.e.,
parameterizing $p({{\bm{\mathbf{f}}}}\,|\,{{\bm{\mathbf{u}}}})$ as follows:

$\displaystyle p({{\bm{\mathbf{f}}}}\,|\,{{\bm{\mathbf{u}}}})$
$\displaystyle\approx\mathcal{N}({{\bm{\mathbf{f}}}};{{\bm{\mathbf{K}}}}_{{{\bm{\mathbf{x}}}}{{\bm{\mathbf{s}}}}}{{\bm{\mathbf{K}}}}_{{{\bm{\mathbf{s}}}}{{\bm{\mathbf{s}}}}}^{-1}{{\bm{\mathbf{u}}}},\text{diag}[{{\bm{\mathbf{K}}}}_{{{\bm{\mathbf{x}}}}{{\bm{\mathbf{x}}}}}-{{\bm{\mathbf{K}}}}_{{{\bm{\mathbf{x}}}}{{\bm{\mathbf{s}}}}}{{\bm{\mathbf{K}}}}_{{{\bm{\mathbf{s}}}}{{\bm{\mathbf{s}}}}}^{-1}{{\bm{\mathbf{K}}}}_{{{\bm{\mathbf{s}}}}{{\bm{\mathbf{x}}}}}])$
(24)
$\displaystyle=\prod_{n=1}^{N}p(f_{n}\,|\,{{\bm{\mathbf{u}}}})=\prod_{n=1}^{N}\mathcal{N}(f_{n};\tilde{\mu}_{n},\tilde{\sigma}^{2}_{n}),$
(25)

where

$\displaystyle\tilde{\mu}_{n}$
$\displaystyle=\kappa({{\bm{\mathbf{x}}}}_{n},{{\bm{\mathbf{S}}}}){{\bm{\mathbf{K}}}}_{{{\bm{\mathbf{s}}}}{{\bm{\mathbf{s}}}}}^{-1}{{\bm{\mathbf{u}}}},$
(26) $\displaystyle\tilde{\sigma}^{2}_{n}$
$\displaystyle=\kappa({{\bm{\mathbf{x}}}}_{n},{{\bm{\mathbf{x}}}}_{n})-\kappa({{\bm{\mathbf{x}}}}_{n},{{\bm{\mathbf{S}}}}){{\bm{\mathbf{K}}}}_{{{\bm{\mathbf{s}}}}{{\bm{\mathbf{s}}}}}^{-1}\kappa({{\bm{\mathbf{S}}}},{{\bm{\mathbf{x}}}}_{n}).$
(27)

We then can decompose the log of joint distribution over all data points as
follows:

$\displaystyle\log
p({{\bm{\mathbf{Z}}}},{{\bm{\mathbf{u}}}}\,|\,{{\bm{\mathbf{\theta}}}})$
$\displaystyle=\log\mathbb{E}_{p({{\bm{\mathbf{f}}}}\,|\,{{\bm{\mathbf{u}}}},{{\bm{\mathbf{\theta}}}})}[p({{\bm{\mathbf{z}}}}\,|\,{{\bm{\mathbf{f}}}};\sigma^{2}{{\bm{\mathbf{I}}}})]$
(28) $\displaystyle=\log\int
p({{\bm{\mathbf{f}}}}\,|\,{{\bm{\mathbf{u}}}},{{\bm{\mathbf{\theta}}}})p({{\bm{\mathbf{Z}}}}\,|\,{{\bm{\mathbf{f}}}},\sigma^{2}{{\bm{\mathbf{I}}}})d{{\bm{\mathbf{f}}}}$
(29)
$\displaystyle=\log\int\prod_{n=1}^{N}p({{\bm{\mathbf{z}}}}_{n}\,|\,f_{n};\sigma^{2})p(f_{n}\,|\,{{\bm{\mathbf{u}}}},{{\bm{\mathbf{\theta}}}})df_{1}...df_{n}$
(30)
$\displaystyle=\log\prod_{n=1}^{N}\int\mathcal{N}({{\bm{\mathbf{z}}}}_{n};f_{n},\sigma^{2})\mathcal{N}(f_{n};\tilde{\mu}_{n},\tilde{\sigma}_{n}^{2})df_{n}$
(31)
$\displaystyle=\log\prod_{n=1}^{N}\mathcal{N}({{\bm{\mathbf{z}}}}_{n};\tilde{\mu}_{n},\tilde{\sigma}_{n}^{2}+\sigma^{2})$
(32)
$\displaystyle=\sum_{n=1}^{N}\log\mathcal{N}({{\bm{\mathbf{z}}}}_{n};\tilde{\mu}_{n},\tilde{\sigma}_{n}^{2}+\sigma^{2}),$
(33)

where $\tilde{\mu}_{n},\tilde{\sigma}_{n}^{2}$ are given by Eq. 26 and Eq. 27,
respectively.

##  Appendix D Details of the extension to deep Gaussian processes

We assume a deep Gaussian process prior $f^{(L)}\circ
f^{(L-1)}\circ\cdots\circ f^{(1)}$ (Damianou & Lawrence, 2013), where each
$f^{(l)}$ is a gp. For the sake of notation, we use
${{\bm{\mathbf{\Psi}}}}^{(l)}$ to indicate the set of kernel hyper-parameters
and inducing inputs of the $l$-th layer and ${{\bm{\mathbf{f}}}}^{(0)}$ as the
input ${{\bm{\mathbf{X}}}}$. We additionally denote
${{\bm{\mathbf{\Psi}}}}=\\{{{\bm{\mathbf{\varphi}}}}\\}\cup\left\\{{{\bm{\mathbf{u}}}}^{(l)},{{\bm{\mathbf{\Psi}}}}^{(l)}\right\\}_{l=1}^{L}$,
where ${{\bm{\mathbf{\varphi}}}}$ is the decoder’s parameters.

The full joint distribution is as follows:

$\displaystyle
p\left({{\bm{\mathbf{\Psi}}}},\\{{{\bm{\mathbf{f}}}}^{(l)}\\}_{l=1}^{L},{{\bm{\mathbf{Z}}}},{{\bm{\mathbf{Y}}}}\,\middle|\,{{\bm{\mathbf{X}}}}\right)=p({{\bm{\mathbf{\Psi}}}})\underbrace{\prod_{l=1}^{L}p\left({{\bm{\mathbf{f}}}}^{(l)}\,|\,{{\bm{\mathbf{f}}}}^{(l-1)},{{\bm{\mathbf{\Psi}}}}\right)p\left({{\bm{\mathbf{Z}}}}\,\middle|\,{{\bm{\mathbf{f}}}}^{(L)};\sigma^{2}{{\bm{\mathbf{I}}}}\right)}_{\text{Deep
GP prior
}}p({{\bm{\mathbf{Y}}}}\,|\,{{\bm{\mathbf{Z}}}},{{\bm{\mathbf{\varphi}}}}).$
(34)

Here we omit the dependency on ${{\bm{\mathbf{X}}}}$ for notational brevity.

We wish to infer the posterior over ${{\bm{\mathbf{\Psi}}}}$ and the latent
variables ${{\bm{\mathbf{Z}}}}$. To do this, the hidden layers
${{\bm{\mathbf{f}}}}^{(l)}$ have to be marginalized and propagated up to the
final layer $L$ (Salimbeni & Deisenroth, 2017). In particular, the log
posterior is as follows:

$\displaystyle\log
p({{\bm{\mathbf{\Psi}}}},{{\bm{\mathbf{Z}}}}\,|\,{{\bm{\mathbf{Y}}}},{{\bm{\mathbf{X}}}})=\log
p({{\bm{\mathbf{\Psi}}}})+\log\mathbb{E}_{p\left(\left\\{{{\bm{\mathbf{f}}}}^{(l)}\right\\}_{l=1}^{L}\,\middle|\,\left\\{{{\bm{\mathbf{u}}}}^{(l)},{{\bm{\mathbf{\Psi}}}}^{(l)}\right\\}_{l=1}^{L}\right)}p\left({{\bm{\mathbf{Z}}}}\,\middle|\,{{\bm{\mathbf{f}}}}^{(L)};\sigma^{2}{{\bm{\mathbf{I}}}}\right)+\log
p({{\bm{\mathbf{Y}}}}\,|\,{{\bm{\mathbf{Z}}}},{{\bm{\mathbf{\varphi}}}})-\log
C,$ (35)

where $C$ is a normalizing constant.

The above posterior is intractable, but we have obtained the form of its (un-
normalized) log joint, from which we can sample using the hmc method.
Unfortunately, the expectation term is intractable, but we can obtain its
estimates using Monte Carlo sampling as follows:

$\displaystyle\log\mathbb{E}_{p\left(\left\\{{{\bm{\mathbf{f}}}}^{(l)}\right\\}_{l=1}^{L}\,\middle|\,\left\\{{{\bm{\mathbf{u}}}}^{(l)},{{\bm{\mathbf{\Psi}}}}^{(l)}\right\\}_{l=1}^{L}\right)}p\left({{\bm{\mathbf{Z}}}}\,\middle|\,{{\bm{\mathbf{f}}}}^{(L)};\sigma^{2}{{\bm{\mathbf{I}}}}\right)\approx$
$\displaystyle\qquad\approx\log\mathbb{E}_{p\left(\left\\{{{\bm{\mathbf{f}}}}^{(l)}\right\\}_{l=2}^{L}\,\middle|\,\widetilde{{{\bm{\mathbf{f}}}}}^{(1)},\left\\{{{\bm{\mathbf{u}}}}^{(l)},{{\bm{\mathbf{\Psi}}}}^{(l)}\right\\}_{l=2}^{L}\right)}p\left({{\bm{\mathbf{Z}}}}\,\middle|\,{{\bm{\mathbf{f}}}}^{(L)};\sigma^{2}{{\bm{\mathbf{I}}}}\right),\qquad\widetilde{{{\bm{\mathbf{f}}}}}^{(1)}\sim
p\left({{\bm{\mathbf{f}}}}^{(1)}\,\middle|\,{{\bm{\mathbf{u}}}}^{(1)},{{\bm{\mathbf{\Psi}}}}^{(1)},{{\bm{\mathbf{f}}}}^{(0)}\right),$
(36)
$\displaystyle\qquad\approx\log\mathbb{E}_{p\left(\left\\{{{\bm{\mathbf{f}}}}^{(l)}\right\\}_{l=3}^{L}\,\middle|\,\widetilde{{{\bm{\mathbf{f}}}}}^{(2)},\left\\{{{\bm{\mathbf{u}}}}^{(l)},{{\bm{\mathbf{\Psi}}}}^{(l)}\right\\}_{l=3}^{L}\right)}p\left({{\bm{\mathbf{Z}}}}\,\middle|\,{{\bm{\mathbf{f}}}}^{(L)};\sigma^{2}{{\bm{\mathbf{I}}}}\right),\quad\widetilde{{{\bm{\mathbf{f}}}}}^{(2)}\sim
p\left({{\bm{\mathbf{f}}}}^{(2)}\,\middle|\,{{\bm{\mathbf{u}}}}^{(2)},{{\bm{\mathbf{\Psi}}}}^{(2)},\widetilde{{{\bm{\mathbf{f}}}}}^{(1)}\right),$
(37) $\displaystyle\qquad\approx\ldots$
$\displaystyle\qquad\approx\log\mathbb{E}_{p\left({{\bm{\mathbf{f}}}}^{(L)}\,\middle|\,\widetilde{{{\bm{\mathbf{f}}}}}^{(L-1)},{{\bm{\mathbf{u}}}}^{(L)},{{\bm{\mathbf{\Psi}}}}^{(L)}\right)}p\left({{\bm{\mathbf{Z}}}}\,\middle|\,{{\bm{\mathbf{f}}}}^{(L)};\sigma^{2}{{\bm{\mathbf{I}}}}\right),\quad\widetilde{{{\bm{\mathbf{f}}}}}^{(L-1)}\sim
p\left({{\bm{\mathbf{f}}}}^{L-1}\left|{{\bm{\mathbf{u}}}}^{(L-1)},{{\bm{\mathbf{\Psi}}}}^{(L-1)},\widetilde{{{\bm{\mathbf{f}}}}}^{(L-2)}\right.\right),$
(38)
$\displaystyle\qquad=\sum_{n=1}^{N}\log\mathbb{E}_{p\left(f_{n}^{L}\,\middle|\,\widetilde{f}_{n}^{L-1},{{\bm{\mathbf{u}}}}^{(L)},{{\bm{\mathbf{\Psi}}}}^{(L)}\right)}p\left({{\bm{\mathbf{z}}}}_{n}\,\middle|\,f_{n}^{(L)};\sigma^{2}\right)$
(39)

Notice that each step of the approximation is unbiased due to the layer-wise
factorization of the joint distribution (Eq. 34). We can obtain a closed form
of the last-layer expectation as ${{\bm{\mathbf{z}}}}_{n}\,|\,f_{n}^{(L)}$ is
a Gaussian. In the case of using a different distribution, this expectation
can be approximated using quadrature (Hensman et al., 2015).

##  Appendix E Experimental details

In this section, we present details on implementation and hyperparameters used
in our experimental campaign. All experiments were conducted on a server
equipped with a Tesla T4 GPU having $16$ GB RAM. Throughout all experiments,
unless otherwise stated, we impose an isotropic Gaussian prior over the
decoder’s parameters. We use the radial basis function (rbf) kernel with
automatic relevance determination (ard) with marginal variance and independent
lengthscales $\lambda_{i}$ per feature (MacKay, 1996). We place a lognormal
prior with unit variance and means equal to and $0.05$ for the lengthscales
and variance, respectively. Since the auxiliary data of most of the considered
datasets are timestamps, we impose a non-informative uniform prior on the
inducing inputs. We observe that this prior works well in our experiments. The
lengthscales are initialized to $1$, whereas the inducing points are
initialized by a $k$-means algorithm as commonly used in practice (Hensman et
al., 2015). For inference, we use an adaptive version of sghmc (Springenberg
et al., 2016) in which the hyper-parameters are automatically tuned during a
burn-in phase. The hyper-parameters are tuned according to the performance on
the validation set.

The random seed for the stochastic inference network is drawn from an
isotropic Gaussian distribution, i.e.
$q({\bm{\mathbf{\varepsilon}}})=\mathcal{N}({\bm{\mathbf{0}}},{{\bm{\mathbf{I}}}})$.
In case the encoder is a multilayer perceptron (mlp), we concatenate the
random seeds and the input into a long vector. The dimension of the random
seeds is the same as that of the input. If the encoder is a convolutional
neural network, we spatially stack the random seeds and the input. We use an
Adam optimizer (Kingma & Ba, 2015) with a learning rate of $0.001$ for
optimizing the encoder network. We set the hyperparameters of the number of
sghmc and encoder optimization steps $J=30$, and $K=50$, respectively.

###  E.1 Moving ball experiment

We use the same network architectures as in Jazbec et al. (2021); Pearce
(2019). We follow the data generation procedure of Jazbec et al. (2021), in
which a squared-exponential gp kernel with a lengthscale $l=2$ was used.
Notice that, unlike Jazbec et al. (2021), we generate a fixed number of $35$
videos for training and another $35$ videos for testing rather than generating
new training videos at each iteration, regardless of the fact that tens of
thousands of iterations are performed. This new setting is more realistic and
is designed to show the data efficiency of the considered methods. The number
of frames in each video is $30$. The dimension of each frame is $32\times 32$.
Table 6 presents the default hyperparameters used for our sgp-bae and Gaussian
Process Bayesian Autoencoder (gp-bae) models. For the competing methods, we
follow the setups specified in Jazbec et al. (2021).

###  E.2 Rotated MNIST experiment

For the rotated MNIST experiment, we follow the same setup as in Jazbec et al.
(2021) and Casale et al. (2018). In particular, we use the gp kernel proposed
by Casale et al. (2018) and a neural network consisting of three convolutional
layers followed by a fully connected layer in the encoder and vice-versa in
the decoder. The details of hyperparameters used for our models are presented
in Table 7. For the competing methods, we refer to Jazbec et al. (2021).

###  E.3 Missing imputation experiment

We follow the experimental setting in Ashman et al. (2020), in which squared
exponential gp kernels are used. Notice that, to ensure a fair comparison, we
handle partially missing data by treating it as zero and feeding it into the
inference network (encoder) for sgp-vae (Ashman et al., 2020) and our sgp-bae
model. This is because it is not trivial to adapt partial inference networks
(Ashman et al., 2020) to our stochastic inference network, and we leave this
for future work. Table 8 and Table 9 show the default hyperparameters used for
our models. For multi-output gp baselines, we refer to Requeima et al. (2019).

Table 6: Parameter settings for the moving ball experiment.

parameter | value
---|---
num. of feedforward layers in encoder | $2$
num. of feedforward layers in decoder | $2$
width of a hidden feedforward layer | $500$
dimensionality of latent space | $2$
activation function | tanh
num. of inducing points | $10$
mini-batch size | full
step size | $0.005$
momentum | $0.05$
num. of burn-in steps | $1500$
num. of samples | $100$
thinning interval | $400$

Table 7: Parameter settings for the rotated mnist experiment.

parameter | value
---|---
num. of conv. layers in encoder | $3$
num. of conv. layers in decoder | $3$
num. of filters per conv. channel | $8$
filter size | $3\times 3$
num. of feedforward layers in encoder | $1$
num. of feedforward layers in decoder | $1$
activation function | elu
dimensionality of latent space | $16$
num. of inducing points | $32$
mini-batch size | $512$
step size | $0.01$
momentum | $0.01$
num. of burn-in steps | $1500$
num. of samples | $200$
thinning interval | $800$

Table 8: Parameter settings for the jura experiment.

parameter | value
---|---
num. of feedforward layers in encoder | $1$
num. of feedforward layers in decoder | $2$
width of a hidden encoder layer | $20$
width of a hidden decoder layer | $5$
dimensionality of latent space | $2$
activation function | relu
num. of inducing points | $128$
mini-batch size | $100$
step size | $0.003$
momentum | $0.05$
num. of burn-in steps | $1500$
num. of samples | $50$
thinning interval | $180$

Table 9: Parameter settings for the eeg experiment.

parameter | value
---|---
num. of feedforward layers in encoder | $1$
num. of feedforward layers in decoder | $2$
width of a hidden encoder layer | $20$
width of a hidden decoder layer | $5$
dimensionality of latent space | $3$
activation function | relu
num. of inducing points | $128$
mini-batch size | $100$
step size | $0.003$
momentum | $0.05$
num. of burn-in steps | $1500$
num. of samples | $50$
thinning interval | $180$

##  Appendix F Additional results

###  F.1 Latent space visualization

Fig. 7 illustrates two-dimensional t-SNE (van der Maaten & Hinton, 2008)
visualization of latent vectors ($C=16$) for the rotated mnist data obtained
by our sgp-bae model. It is evident that the clusters of embeddings are
organized in a clear and orderly manner according to the angles they
represent. Specifically, the embeddings of rotated images are arranged in a
continuous sequence from $0$ to $2\pi$ in a clockwise direction.

###  F.2 Convergence of SGHMC

We follow the common practice using $\hat{R}$ potential scale-reduction
diagnostic (Gelman & Rubin, 1992) to assess the convergence of Markov chain
Monte Carlo (mcmc) on Bayesian neural networks (bnns) (Izmailov et al., 2021b;
Tran et al., 2022). Given two or more chains, $\hat{R}$ estimates the ratio
between the chain variance and the average within-chain variance. For the
large-scale mnist experiment, we compute the $\hat{R}$ statistics on the
predictive posterior over four independent chains and obtain a median value of
less than $1.1$, which suggests good convergence (Izmailov et al., 2021b). In
addition, we present traceplots of the predictive posterior in Fig. 8, which
also support the conclusion of good mixing.

###  F.3 Visualization of EEG data

Fig. 9 depicts the predictive mean and uncertainty estimation for the missing
values of the eeg dataset.

Figure 7: Visualization of t-SNE embeddings (van der Maaten & Hinton, 2008) of
sgp-bae latent vectors on the training data of the rotated mnist. Each image
embedding is colored with respect to its associated angle. Here, we use a
perplexity parameter of $50$ for t-SNE.

Figure 8: Trace plots for four test points on the rotated mnist dataset. Here,
we simulate $4$ chains with $200$ samples for each.

(a) Independent Gaussian Processes (igp)

(b) Gaussian Process Autoregressive Model (gpar)

(c) sgp-bae

| Observed values |  | Missing values |  | Predicted mean |  |  $\pm$ 3 standard deviation
---|---|---|---|---|---|---|---

Figure 9: Visualization of predictions for missing data of the eeg dataset.
Each panel shows one of the three channels with missing data (orange crosses)
and observed data (black points).