Tr

# Hierarchical Gaussian Processes with Wasserstein-2 Kernels

Sebastian G. Popescu
Imperial College London
<EMAIL_ADDRESS>
&David J. Sharp
Imperial College London
<EMAIL_ADDRESS>
&James H. Cole
University College London
<EMAIL_ADDRESS>
Ben Glocker
Imperial College London
<EMAIL_ADDRESS>

###### Abstract

Stacking Gaussian Processes severely diminishes the model’s ability to detect
outliers, which when combined with non-zero mean functions, further
extrapolates low non-parametric variance to low training data density regions.
We propose a hybrid kernel inspired from Varifold theory, operating in both
Euclidean and Wasserstein space. We posit that directly taking into account
the variance in the computation of Wasserstein-2 distances is of key
importance towards maintaining outlier status throughout the hierarchy. We
show improved performance on medium and large scale datasets and enhanced out-
of-distribution detection on both toy and real data.

_Keywords_ Gaussian Processes $\cdot$ Out-of-distribution Detection $\cdot$
Uncertainty Quantification

##  1 Introduction

Gaussian Processes (GP) are data-efficient Bayesian non-parametric models that
offer calibrated uncertainty quantification and are robust to overfitting,
recently finding applicability in data-scarce domains (Timonen et al., 2019;
Wang et al., 2020) or where uncertainty is of utmost importance (Chen et al.,
2014). Deep Gaussian Processes (DGP) (Damianou and Lawrence, 2013) are multi-
layered generalizations of GPs that share data-efficiency learning
capabilities and show enhanced representational capacity. However, the
reliability of uncertainties propagated throughout the hierarchy is up for
debate. Recent work (Ustyuzhaninov et al., 2019) has questioned the
uncertainties present in the hidden layers of DGP, showing that approximate
inference schemes using factorized variational Gaussian distributions result
in all but the last GP collapsing to deterministic transformations in the case
of noise-free data. Such pathological behaviour should be avoided as it
undermines the utility of multi-layered GP. We propose desideratas for useful
properties that a kernel must have in a DGP and based on this we introduce an
alternative kernel operating in both Euclidean and Wasserstein-2 space,
showcasing that the explicit distance between multivariate Gaussians as
present in hidden layers is of vital importance towards correctly propagating
non-parametric uncertainty forward. Prior work (Bachoc et al., 2017, 2018; Thi
Thien Trang et al., 2019) has only focused on shallow GPs operating on inputs
which do not have an explicit distribution, thus using Wasserstein-2 based
kernels to take advantage of the distributional nature of the input. We did
not find any other work that expanded the concept of Wasserstein-2 based
kernels to DGPs. In the remainder of this section we introduce two main
motivations behind this paper.

##### Distributional Uncertainty.

Kendall and Gal (2017) propose a dichotomy of classifying sources of
uncertainty into aleatoric or epistemic, pending on their reducibility status.
More specifically, aleatoric uncertainty is caused by inherent randomness in
the data, hence it cannot be reduced. The classical example of coin flipping
is the most obvious one to give, even though one can also argue that coin
flipping is deterministic given a complex system taking into consideration
strength of coin tosser, initial side, wind conditions and physics of coins
tossing. Conversely, epistemic uncertainty refers to the epistemic state of
the agent/model, or more colloquially the ignorance with respect to the
current state. This ignorance can be reduced by feeding the agent/model with
more information/data. Besides the dichotomy consisting of aleatoric and
epistemic uncertainty, reliably highlighting certain inputs which have
undergone a domain shift (Lakshminarayanan et al., 2017) or out-of-
distribution (OOD) samples (Hendrycks and Gimpel, 2016) has garnered a lot of
interest in the past years. Succinctly, the aim is to measure the degree to
which a model knows when it does not know, or more precisely if a network
trained on a specific dataset is evaluated at testing time on a completely
different dataset (potentially from a different modality or another
application domain), then the expectation is that the network should output
high predictive uncertainty on this set of data points that are very far from
the training data manifold. We remind ourselves that epistemic uncertainty can
be reduced by adding more data. By this logic, epistemic uncertainty cannot be
reduced outside the data manifold of our dataset since we don’t add data
points which do not stem from the same data generative pipeline (this is not
true in the case of OOD detection models which explicitly use OOD samples
during training/testing (Liang et al., 2017; Hafner et al., 2020)). Hence,
epistemic uncertainty can only be reduced inside the data manifold and should
be zero outside the data manifold (under the assumption that the model is
_distance-aware_ , a concept subsequently introduced). Conversely, our choosen
measure for OOD detection should grow outside the data manifold and be close
to or 0 inside the data manifold. Malinin and Gales (2018) introduced for the
first time the separation of total uncertainty into three components:
epistemic, aleatoric and distributional uncertainty. To make the distinction
clearer, the authors argue that aleatoric uncertainty is a "known-unknown",
whereby the model confidently states that an input data point is hard to
classify (class overlap). Contrary, distributional uncertainty is an "unknown-
unknown" due to the fact that the model is unfamiliar with the input space
region that the test data comes from, thereby not being able to make confident
predictions.

##### Do Hierarchical GPs reliably propagate distributional uncertainty?

Since their inception, research on deep GPs has predominantly focused on
enhancing their approximate inference framework. Little attention has been
paid so far in how deep GPs fare on OOD detection. Domingues et al. (2018)
have proposed an autoencoder architecture with both encoder and decoder
governed by deep GPs to detect outliers. However, this line of inquiry stems
from reconstruction-based methods for outlier detection (Baur et al., 2021),
whereby at training time an autoenconder is trained on normative data,
subsequently at testing time the model is presented with OOD samples. The
differences between reconstruction and OOD sample represent the regions in
input space which are anomalous. Conversely, in our line of research we are
interested in accurate predictive models which also have OOD detection
capabilities. Deep GPs have thus so far not been tested if they are reliable
OOD detection models. This will represent a major focus of this paper.

##### Contributions

To summarize, this paper makes the following advancements:

  * •

We provide empirical evidence that as the number of inducing point increases,
the non-parametric variance in hidden layers collapses in the case of zero-
mean function DGP.

  * •

We show that using Principal Component Analysis (PCA) mean functions solves
this problem, but this procedure induces low variance perpendicular to the
eigenvectors in the case of encoding down high-dimensional data.

  * •

We derive theoretical requirements for DGPs to propagate outliers with high
total uncertainty, respectively high non-parametric variance throughout the
hierarchy.

  * •

We propose a variant of DGPs that operates on both Euclidean and Wasserstein-2
space. We demonstrate that this formulation guarantees higher non-parametric
variance for OOD points compared to in-distribution points, alongside
increased chances of not suffering from _feature-collapse_.

  * •

We empirically show enhanced out-of-distribution detection capabilities,
alongside improved test time results on a wide range of datasets for our
proposed model.

##  2 Background

In this section we introduce the necessary theoretical background needed for
this paper.

###  2.1 Primer on Gaussian Processes

A stochastic process assigns a random variable $f(x)$ to every element of an
index set $\mathbb{X}$, more precisely $f=\\{f(x)\\}_{x\in\mathbb{X}}$.
Stochastic processes can be defined via their finite dimensional marginal
distributions, more concretely the stochastic process can be defined as a
random function $f:\mathbb{X}\to\mathbb{Y}$, with any finite sequence
$x_{1:n}=\\{x_{1},\cdots,x_{n}\\}$ having a joint marginal distribution over
their values
$y_{1:n}=\\{y_{1},\cdots,y_{n}\\}=\\{f(x_{1}),\cdots,f(x_{n})\\}$. Stochastic
processes must respect two key conditions from the Kolmogorov Extension
Theorem, which we now outline. Consider finite sequence $x_{1:n}\in\mathbb{X}$
with corresponding probability measure $p_{x_{1:n}}$ on
$\left(\mathbb{R}^{D}\right)^{n}$.

##### Exchangeability.

For all permutations $\pi$ of $x_{1:n}=\\{x_{1},\cdots,x_{n}\\}$ we impose
$p_{x_{1:n}}(y_{1:n})=p_{\pi\left(x_{1:n}\right)}(\pi\left(y_{1:n}\right))$
which translates into invariance of the joint distribution with respect to
permutations of elements in any given finite sequence.

##### Consistency.

$p_{x_{1:m}}(y_{1:m})=\int p_{x_{1:n}}(y_{1:n})dy_{m+1:n}$ which ensures that
if we marginalize a subset of elements, the remainder will remain unchanged.

Random functions are intrinsically complex mathematical objects with a
Gaussian process representing the simplest random function due to its complete
characterization by its mean and covariance functions (Pugachev, 2013). A
Gaussian Process can be seen as a generalization of multivariate Gaussian
random variables to infinite sets. We define this statement in more detail
now. We consider $f(x)$ to be a stochastic field, with $x\in\mathbb{R}^{D}$
and we define $m(x)=\mathbb{E}\left[f(x)\right]$ and
$C(x_{i},x_{j})=Cov\left[f(x_{i}),f(x_{j})\right]$. We denote a Gaussian
Process (GP) $f(x)$ as:

$f(x)\sim\mathcal{GP}\left(m(x),C\left(x_{i},x_{j}\right)\right)$ (1)

Throughout this paper, we will generally consider zero-mean processes, hence
our GPs are strictly characterized by their covariance functions. The latter
have the condition to generate non-negative-definite covariance matrices, more
specifically they have to satisfy:

$\sum_{i,j}a_{i}a_{j}C\left(x_{i},x_{j}\right)\geq 0$ (2)

for any finite set $\\{x_{1},\cdots,x_{n}\\}$ and any real valued coefficients
$\\{a_{1},\cdots,a_{n}\\}$. Throughout this paper we will only consider
second-order stationary processes which have constant means and
$Cov\left[f(x_{i}),f(x_{j})\right]=C\left(\|x_{i}-x_{j}\|\right)$. We can see
that such covariance functions are invariant to translations.

The squared exponential/radial basis function kernel is a commonly used
stationary kernel, which has the following form:

$k^{SE}(x_{i},x_{j})=\sigma^{2}\exp{\left[\sum_{d=1}^{D}-\frac{\left(x_{i,d}-x_{j,d}\right)^{2}}{l^{2}}\right]}$
(3)

, where we have written its definition in the anisotropic case. Intuitively,
the lengthscale values $\\{l_{1}^{2},\cdots,l_{D}^{2}\\}$ represent the
strength along a particular dimension of input space by which successive
values are strongly correlated with correlation invariably decreasing as the
distance between points increases. Such a covariance function has the property
of Automatic Relevance Determination (ARD) (Neal, 2012). However, the usage of
ARD to determine the relative importance of a given feature in comparison to
others has been recently challenged (Paananen et al., 2019). Lastly, the
kernel variance $\sigma^{2}$ controls the variance of the process, more
specifically the amplitude of function samples.

A GP has the following joint distribution over finite subsets
$\mathbb{X}_{1}\in\mathbb{X}$ with function values $f(X_{1})$. Analogously for
$\mathbb{X}_{2}$, with their union being denoted as
$x=\\{x_{1},\cdots,x_{n}\\}$.

$p\begin{pmatrix}f(x_{1})\\\
f(x_{2})\end{pmatrix}=\mathcal{N}\left[\begin{pmatrix}m(x_{1})\\\
m(x_{2})\end{pmatrix},\begin{pmatrix}k(x_{1},x_{1}),k(x_{1},x_{2})\\\
k(x_{2},x_{1}),k(x_{2},x_{2})\end{pmatrix}\right]$ (4)

The following observation model is used:

$p(y|f,x)=\prod_{i=1}^{N}p(y_{i}\mid f(x_{i}))$ (5)

, where given a supervised learning scenario, the dataset
$D=\\{x_{i},y_{i}\\}_{i=1,\cdots,n}$ can be shorthand denoted as
$D=\\{x,y\\}$. In the case of probabilistic regression, we make the assumption
that the noise is additive, independent and Gaussian, such that the latent
function $f(x)$ and the observed noisy outputs $y$ are defined by the
following equation:

$y_{i}=f(x_{i})+\epsilon_{i}\text{, where
}~{}\epsilon_{i}\sim\mathcal{N}\left(0,\sigma^{2}_{noise}\right)$ (6)

To train a GP for regression tasks, one performs Marginal Likelihood
Maximization of Type 2 over the following equation:

$p(y)=\mathcal{N}\left(y\mid m,K+\sigma^{2}_{noise}\mathbb{I}_{n}\right)$ (7)

by treating the kernel hyperparameters as point-mass. In the case of priors
being placed on the hyperparameters, Markov Chain Monte Carlo (MCMC) methods
can be used (Lalchand and Rasmussen, 2020).

We are interested in finding the posterior $p\left(f(x^{*})\mid y\right)$
since the goal is to predict for unseen data points $x^{*}$ which are
different than the training set. We know that the joint prior over training
and testing set latent functions is given by:

$p\begin{pmatrix}f(x)\\\
f(x^{*})\end{pmatrix}=\mathcal{N}\left[\begin{pmatrix}m(x)\\\
m(x^{*})\end{pmatrix},\begin{pmatrix}k(x,x)&k(x,x^{*})\\\
k(x^{*},x)&k(x^{*},x^{*})\end{pmatrix}\right]$ (8)

and using marginalization properties of multivariate normals we can then write
the following : $f(x^{*})\sim\mathcal{N}\left(m(x^{*}),k(x^{*},x^{*})\right)$.
With this in mind, and also reiterating that the marginal likelihood is also a
multivariate normal given by
$y\sim\mathcal{N}\left(m(x),k(x,x,)+\sigma^{2}_{noise}\mathbb{I}_{n}\right)$,
we can then write their joint prior since both of them are Gaussian:

$p\begin{pmatrix}y\\\
f(x^{*})\end{pmatrix}=\mathcal{N}\left[\begin{pmatrix}m(x)\\\
m(x^{*})\end{pmatrix},\begin{pmatrix}k(x,x)+\sigma^{2}_{noise}\mathbb{I}_{n}&k(x,x^{*})\\\
k(x^{*},x)&k(x^{*},x^{*})\end{pmatrix}\right]$ (9)

Now we can simply apply the conditional rule for multivariate Gaussians to
obtain:

$\displaystyle p\left(f(x^{*})\mid y\right)$
$\displaystyle=\mathcal{N}(f(x^{*})\mid
m(x^{*})+K_{f^{*}f}\left[K_{ff}+\sigma^{2}_{noise}\mathbb{I}_{n}\right]^{-1}\left[y-m(x)\right],$
(10) $\displaystyle\hskip
28.45274ptK_{f^{*}f^{*}}-K_{f^{*}f}\left[K_{ff}+\sigma^{2}_{noise}\mathbb{I}_{n}\right]^{-1}K_{ff^{*}})$

The usage of GP in real-life datasets is hindered by the $\mathbb{O}(n^{3})$
time, $\mathbb{O}(n^{2})$ memory for training, where $n$ is the number of data
points in the training set. In the next subsection we will see how to avert
having to incur these expensive computational budgets.

###  2.2 Variational Free Energy approximation

Titsias (2009) introduced the first variational lower bound comprising a
probabilistic regression model over inducing points. More specifically, the
authors applied variational inference (VI) in an augmented probability space
that comprises training set latent function values $F$ alongside inducing
point latent function values $U$:

$p(y,F,U)=p(y\mid F)p(F\mid U)p(U)$ (11)

, where the conditional of the training set latent functions based on inducing
points is taken to be exact $p(F\mid U)=\mathcal{N}\left(F\mid
K_{fu}K_{uu}^{-1}U,K_{ff}-K_{fu}K_{uu}^{-1}K_{uf}\right)$.

In terms of doing exact inference in this new model, respectively computing
the posterior $p(f|y)$ and the marginal likelihood $p(y)$, it remains
unchanged even with the augmentation of the probability space by $U$ as we can
marginalize $p(F)=\int p(F,U)~{}dU$ due to the marginalization properties of
Gaussian processes. Succinctly, $p(F)$ is not changed by modifying the values
of $U$, even though $p(F\mid U)$ and $p(U)$ do indeed change. This translates
into the fundamental difference between variational parameters $U$ and
hyperparameters of the model
$\\{\sigma^{2}_{noise},\sigma^{2},l_{1}^{2},\cdots,l_{D}^{2}\\}$, whereby the
introduction of more variational parameters does not change the fundamental
definition of the model before probability space augmentation.

In Titsias (2009), the authors have derived a collapsed variational lower
bound which enties using the entire training set to derive the lower bound to
the log marginal likelihood, which is used for training the hyperparameters
$\\{\sigma^{2}_{noise},\sigma^{2},l_{1}^{2},\cdots,l_{D}^{2}\\}$ of the GP,
respectively for deriving the optimal posterior for $U$. Doing exact inference
in the VFE formulation of sparse GP (SGP) is of sound theoretical interest,
however it becomes infeasible for even medium sized datasets. In this
subsection we will see how we can train VFE formulated SGP using minibatches
by keeping an explicit representation over $q(U)$.

Stochastic Variational Inference (SVI) (Hoffman et al., 2013) enables the
application of VI for extremely large datasets, by virtue of performing
inference over a set of global variables, which induce a factorisation in the
observations and latent variables, such as in the Bayesian formulation of
Neural Networks with distributions (implicit or explicit) over matrix weights.
GP do no exhibit these properties, but by virtue of the approximate prior over
testing and training latent functions for SGP approximations with inducing
points $U$, which we remind here for the VFE scenario:

$p(f,f^{*})\approx q(f,f^{*})=\int p(f\mid U)p(f^{*}\mid U)p(U)~{}dU$ (12)

This translates into a fully factorized model with respect to observations at
training and testing time, conditioned on the global variables $U$. One can
easily see that by integrating out $U$ in the Titsias bound, this key property
for SVI dissipates. In the remainder of this subsection, we will derive a
lower bound to the Titsias bound that satisfies the SVI requirements.

Our goal is to approximate the true posterior distribution $p(F,U\mid
y)=p(F\mid U,Y)p(U\mid Y)$ by introducing a variational distribution $q(F,U)$
and minimizing the Kullback-Leibler divergence:

$KL\left[q(F,U)\|p(F,U\mid y)\right]=\int q(F,U)\log\frac{q(F,U)}{p(F,U\mid
y)}~{}dF~{}dU$ (13)

where the approximate posterior is factorized as $q(F,U)=p(F\mid U)q(U)$ and
$q(U)$ is an unconstrained variational distribution over $U$. Following the
standard VI framework we need to maximize the following variational lower
bound on the log marginal likelihood:

$\displaystyle\log p(y)$ $\displaystyle\geq\int p(F\mid
U)q(U)\log\frac{p(y\mid F)p(F\mid U)p(U)}{p(F\mid U)p(U)}~{}dF~{}dU$ (14)
$\displaystyle\geq\int q(U)\left[\int\log p(Y\mid F)p(F\mid
U)~{}dF+\log\frac{p(U)}{q(U)}\right]~{}dU$ (15)

We can now solve for the integral over $F$:

$\displaystyle\int\log p(y|F)p(F|U)~{}dF$
$\displaystyle=\mathbb{E}_{p(F|U)}\left[-\frac{n}{2}\log(2\pi\sigma^{2}_{noise})-\frac{1}{2\sigma^{2}_{noise}}Tr\left[yy^{\top}-2yF^{\top}+FF^{\top}\right]\right]$
(16)
$\displaystyle=-\frac{n}{2}\log(2\pi\sigma^{2}_{noise})-\frac{1}{2\sigma^{2}_{noise}}Tr[yy^{\top}-2y\left(K_{fu}K_{uu}^{-1}U\right)^{\top}+$
(17)
$\displaystyle\left(K_{fu}K_{uu}^{-1}U\right)\left(K_{fu}K_{uu}^{-1}U\right)^{\top}+K_{ff}-Q_{ff}]$
$\displaystyle=\log\mathcal{N}\left(y|K_{fu}K_{uu}^{-1}U,\sigma^{2}_{noise}\mathbb{I}_{n}\right)-\frac{1}{2\sigma^{2}_{noise}}Tr\left[K_{ff}-Q_{ff}\right]$
(18)

We can now rewrite our variational lower bound as follows:

$\log p(y)\geq\int q(U)\log\frac{\mathcal{N}\left(y\mid
K_{fu}K_{uu}^{-1}U,\sigma^{2}_{noise}\mathbb{I}_{n}\right)p(U)}{q(U)}~{}dU-\frac{1}{2\sigma^{2}_{noise}}Tr\left[K_{ff}-Q_{ff}\right]$
(19)

The variational posterior is explicit in this case, respectively
$q(F,U)=p(F\mid U;X,Z)q(U)$, where $q(U)=\mathcal{N}(U\mid m_{U},S_{U})$.
Here, $m_{U}$ and $S_{U}$ are free variational parameters. Due to the Gaussian
nature of both terms we can marginalize $U$ to arrive at $q(F)=\int p(F\mid
U)q(U)=\mathcal{N}(F\mid\tilde{U}(x),\tilde{\Sigma}(x))$, where:

$\displaystyle\tilde{U}(x)$ $\displaystyle=K_{fu}K_{uu}^{-1}m_{U}$ (20)
$\displaystyle\tilde{\Sigma}(x)$
$\displaystyle=K_{ff}-K_{fu}K_{uu}^{-1}\left[K_{uu}-S_{U}\right]K_{uu}^{-1}K_{uf}$
(21)

The lower bound can be re-expressed as follows:

$\displaystyle\log p(y)$ $\displaystyle\geq\int
q(U)\log\mathcal{N}\left(y|K_{fu}K_{uu}^{-1}U,\sigma^{2}\mathbb{I}\right)~{}dU-
KL\left[q(U)\|p(U)\right]-\frac{1}{2\sigma^{2}_{noise}}Tr\left[K_{ff}-Q_{ff}\right]$
(22)
$\displaystyle\geq\mathcal{N}\left(y|K_{fu}K_{uu}^{-1}m_{U},\sigma^{2}_{noise}\mathbb{I}_{n}\right)-\frac{1}{2\sigma^{2}_{noise}}Tr\left[K_{fu}K_{uu}^{-1}S_{U}K_{uu}^{-1}K_{uf}\right]$
(23) $\displaystyle\hskip
28.45274pt-\frac{1}{2\sigma^{2}_{noise}}Tr\left[K_{ff}-Q_{ff}\right]-KL\left[q(U)\|p(U)\right]=\mathcal{L}_{SVGP}$

, where we can easily see that the last equation factorized with respect to
individual observations. This lower variational bound will be denoted as the
Sparse Variational GP (SVGP) or uncollapsed SVGP bound and will constitute a
stepping-stone in training our proposed model.

###  2.3 Approximation of Gaussian Processes with noisy inputs

In subsection 2.1 we have seen that the predictive distribution of the output
at a unseen noise-free input is Gaussian. In Girard (2004), the authors
propose several approximation techniques devising approximations of the non-
Gaussian predictive distribution when the input is noisy (or a random variable
in a more general scenario as we shall later see). In this subsection we will
attempt to make a summary of aforementioned work, keeping in mind that this
approximation techniques will play a key role in analyzing properties of deep
GP and how they propagate uncertainty through the hierarchy.

We now consider the case where the input is corrupted by noise
$\epsilon_{x^{*}}\sim\mathcal{N}\left(0,\Sigma_{x^{*}}\right)$ such that
$x^{*}=\mu_{x^{*}}+\epsilon_{x^{*}}$. Hence, our aim is to make predictions at
$x^{*}\sim\mathcal{N}\left(\mu_{x^{*}},\Sigma_{x^{*}}\right)$, which involves
the following integral:

$p\left(y^{*}\mid D,\mu_{x^{*}},\Sigma_{x^{*}}\right)=\int p\left(y^{*}\mid
D,x^{*}\right)p\left(x^{*}\mid\mu_{x^{*}},\Sigma_{x^{*}}\right)~{}dx^{*}$ (24)

where $p\left(y^{*}\mid
D,x^{*}\right)=\frac{1}{\sqrt{2\pi\sigma^{2}(x^{*})}}\exp{\left[-\frac{1}{2}\frac{\left(y^{*}-\mu(x^{*})\right)^{2}}{\sigma^{2}(x^{*})}\right]}$,
respectively $\mu(x^{*})=K_{f^{*}f}K_{ff}^{-1}y$ and
$\sigma^{2}(x^{*})=K_{f^{*}f^{*}}-Q_{f^{*}f^{*}}$. In this subsection, for
notation purposes, we will consider
$K_{x_{i},x_{j}}=K^{SE}(x_{i},x_{j})+\delta_{i,j}\sigma^{2}_{noise}$. One can
notice that $p\left(y^{*}\mid D,x^{*}\right)$ is a non-linear function of
$x^{*}$, so we cannot integrate out $x^{*}$. We will now consider
approximations to the solution of this integral.

##### I Monte Carlo Integration

uses random sampling of a function to numerically compute an estimate of its
integral.

$p^{MC}\left(y^{*}\mid
D,\mu_{x^{*}},\Sigma_{x^{*}}\right)\approx\frac{1}{T}\sum\limits_{t=1}^{T}p\left(y^{*}\mid
D,x_{t}^{*}\right)=\frac{1}{T}\sum\limits_{t=1}^{T}\mathcal{N}\left(y^{*}\mid\mu(x^{*}),\sigma^{2}(x^{*})\right)$
(25)

It will result in a mixture of T Gaussians with equal mixing weights and as
$T\to\infty$ we have $p^{MC}\left(y^{*}\mid
D,\mu_{x^{*}},\Sigma_{x^{*}}\right)\approx p\left(y^{*}\mid
D,\mu_{x^{*}},\Sigma_{x^{*}}\right)$. We will later see that the case $T=1$
will be default method to computing this type of integrals within the
framework of Deep Gaussian Processes.

##### II Gaussian Approximation

involves approximating the non-Gaussian predictive equation with a Gaussian:

$p\left(y^{*}\mid
D,\mu_{x^{*}},\Sigma_{x^{*}}\right)=\mathcal{N}\left(y^{*}\mid
m(u_{x^{*}},\Sigma_{x^{*}}),v(u_{x^{*}},\Sigma_{x^{*}})\right)$ (26)

We can now explicitly write the first two moments:

$\displaystyle m(u_{x^{*}},\Sigma_{x^{*}})$ $\displaystyle=\int
y^{*}\left[\int p\left(y^{*}\mid D,x^{*}\right)p\left(x^{*}\mid
u_{x^{*}},\Sigma_{x^{*}}\right)dx^{*}\right]dy^{*}$ (27)
$\displaystyle=\int\left[\int y^{*}p\left(y^{*}\mid
D,x^{*}\right)dy^{*}\right]p\left(x^{*}\mid
u_{x^{*}},\Sigma_{x^{*}}\right)dx^{*}$ (28)
$\displaystyle=\int\mu(x^{*})p\left(x^{*}\mid
u_{x^{*}},\Sigma_{x^{*}}\right)dx^{*}$ (29)
$\displaystyle=\mathbb{E}_{x^{*}}\left[\mu(x^{*})\right]$ (30) $\displaystyle
v(u_{x^{*}},\Sigma_{x^{*}})$
$\displaystyle=\int\left(y^{*}\right)^{2}\left[\int p\left(y^{*}\mid
D,x^{*}\right)p\left(x^{*}\mid
u_{x^{*}},\Sigma_{x^{*}}\right)dx^{*}\right]dy^{*}-m(u_{x^{*}},\Sigma_{x^{*}})^{2}$
(31) $\displaystyle=\int\left[\int\left(y^{*}\right)^{2}p\left(y^{*}\mid
D,x^{*}\right)dy^{*}\right]p\left(x^{*}\mid
u_{x^{*}},\Sigma_{x^{*}}\right)dx^{*}-m(u_{x^{*}},\Sigma_{x^{*}})^{2}$ (32)
$\displaystyle=\int\left[\sigma^{2}(x^{*})+\mu(x^{*})^{2}\right]p\left(x^{*}\mid
u_{x^{*}},\Sigma_{x^{*}}\right)dx^{*}-m(u_{x^{*}},\Sigma_{x^{*}})^{2}$ (33)
$\displaystyle=\mathbb{E}_{x^{*}}\left[\sigma^{2}(x^{*})\right]+\mathbb{E}_{x^{*}}\left[\mu(x^{*})^{2}\right]-m(u_{x^{*}},\Sigma_{x^{*}})^{2}$
(34)
$\displaystyle=\mathbb{E}_{x^{*}}\left[\sigma^{2}(x^{*})\right]+V_{x^{*}}\left[\mu(x^{*})\right]$
(35)

, where we now consider just predictive distributions over a single data point
$x^{*}$, hence we update predictive equations to correspond to the univariate
prediction case:

$\displaystyle\mu(x^{*})$ $\displaystyle=K_{x^{*}f}K_{ff}^{-1}y$ (36)
$\displaystyle\sigma^{2}(x^{*})$
$\displaystyle=K_{x^{*}x^{*}}-\sum\limits_{i,j=1}^{N}K_{ij}^{-1}K(x^{*},x_{i})K(x^{*},x_{j})$
(37)

To recap, we are interested in finding the first two moments of the Gaussian
approximation in equation (26), which we now write explicitly:

$\displaystyle m(u_{x^{*}},\Sigma_{x^{*}})$
$\displaystyle=\mathbb{E}_{x^{*}}\left[K_{x^{*}f}K_{ff}^{-1}y\right]$ (38)
$\displaystyle v(u_{x^{*}},\Sigma_{x^{*}})$
$\displaystyle=\mathbb{E}_{x^{*}}\left[K_{x^{*}x^{*}}-\sum\limits_{i,j=1}^{N}K_{ij}^{-1}K(x^{*},x_{i})K(x^{*},x_{j})\right]$
(39)
$\displaystyle+\mathbb{E}_{x^{*}}\left[\left[K_{ff}^{-1}y\right]_{i,:}\left[K_{ff}^{-1}y\right]_{j,:}K_{x^{*},x_{i}}K_{x^{*},x_{j}}\right]$
$\displaystyle-m(u_{x^{*}},\Sigma_{x^{*}})^{2}$

, where $\left[K_{ff}^{-1}y\right]_{i,:}$ takes the i-th row. Therefore, we
need to compute the following integrals:

$\displaystyle\Psi=\mathbb{E}_{x^{*}}\left[K_{x^{*},x^{*}}\right]$
$\displaystyle=\int K(x^{*},x^{*})p(x^{*})~{}dx^{*}$ (40)
$\displaystyle\Psi_{i}=\mathbb{E}_{x^{*}}\left[K_{x^{*},x_{i}}\right]$
$\displaystyle=\int K(x^{*},x_{i})p(x^{*})~{}dx^{*}$ (41)
$\displaystyle\Psi_{i,j}=\mathbb{E}_{x^{*}}\left[K_{x^{*},x_{i}}K_{x^{*},x_{j}}\right]$
$\displaystyle=\int K(x^{*},x_{i})K(x^{*},x_{j})p(x^{*})~{}dx^{*}$ (42)

These integrals will have an analytic solution in the case of specific
kernels. We will first focus on the general case, whereby we devise
approximations for any arbitrary kernel, subsequently calculating the exact
moments for squared exponential kernels.

####  2.3.1 Approximate Moments

Moment approximations involves approximating the integrand by a Taylor
expansion.

##### General case.

To illustrate this method, we take a simple scenario where $x$ is a random
variable such that $\mathbb{E}_{x}\left[x\right]=\mu_{x}$ and
$V_{x}\left[x\right]=\Sigma_{x}$. For a well-behaved function $f$ with
$y=f(x)$ and sufficiently small $\sigma_{x}=\sqrt{\Sigma_{x}}$ we can consider
the Taylor expansion of $f(x)$ at $\mu_{x}$ up to second order to be an
accurate representation of $f(x)$:

$y=f(x)\approx
f(\mu_{x})+\left(x-\mu_{x}\right)f^{{}^{\prime}}(\mu_{x})+\frac{1}{2}\left(x-\mu_{x}\right)^{2}f^{{}^{\prime\prime}}(\mu_{x})$
(43)

Using the Taylor expansion, we can calculate the first two moments as follows:

$\displaystyle\mathbb{E}_{x}\left[y\right]$
$\displaystyle\approx\mathbb{E}_{x}\left[f(\mu_{x})\right]+\mathbb{E}_{x}\left[\left(x-\mu_{x}\right)f^{{}^{\prime}}(\mu_{x})\right]+\mathbb{E}_{x}\left[\frac{1}{2}\left(x-\mu_{x}\right)^{2}f^{{}^{\prime\prime}}(\mu_{x})\right]$
(44)
$\displaystyle\approx\mathbb{E}_{x}\left[f(\mu_{x})\right]+\frac{1}{2}\Sigma_{x}f^{{}^{\prime\prime}}(\mu_{x})$
(45) $\displaystyle V_{x}\left[y\right]$
$\displaystyle\approx\mathbb{E}_{x}\left[y^{2}\right]-\mathbb{E}_{x}\left[y\right]^{2}$
(46) $\displaystyle\approx
f(\mu_{x})^{2}+\Sigma_{x}f^{{}^{\prime}}(\mu_{x})^{2}+f(\mu_{x})f^{{}^{\prime\prime}}(\mu_{x})\Sigma_{x}-f(\mu_{x})^{2}-f(\mu_{x})f^{{}^{\prime\prime}}(\mu_{x})\Sigma_{x}$
(47) $\displaystyle\approx\Sigma_{x}f^{{}^{\prime}}(\mu_{x})^{2}$ (48)

##### GP with noisy inputs.

We now apply the above results to equations (30) and (35), obtaining:

$\displaystyle m(u_{x^{*}},\Sigma_{x^{*}})$
$\displaystyle=\mathbb{E}_{x^{*}}\left[\mu(x^{*})\right]$ (49)
$\displaystyle=\mu(\mathbb{E}_{x^{*}}(x^{*}))+\frac{1}{2}V_{x^{*}}(x^{*})\mu^{{}^{\prime\prime}}(\mathbb{E}_{x^{*}}\left[x^{*}\right])$
(50) $\displaystyle v(u_{x^{*}},\Sigma_{x^{*}})$
$\displaystyle=\mathbb{E}_{x^{*}}\left[\sigma^{2}(x^{*})\right]+\mathbb{E}_{x^{*}}\left[\mu(x^{*})^{2}\right]-m(u_{x^{*}},\Sigma_{x^{*}})^{2}$
(51)
$\displaystyle=\sigma^{2}(\mathbb{E}_{x^{*}}\left[x^{*}\right])+V_{x^{*}}\left[\frac{1}{2}\sigma^{2}(\mathbb{E}_{x^{*}}\left[x^{*}\right])^{{}^{\prime\prime}}+\mu^{{}^{\prime}}(\mathbb{E}_{x^{*}}\left[x^{*}\right])^{2}\right]$
(52)

####  2.3.2 Exact Moments

We can calculate the exact moments in the case of squared exponential kernels
of the following form:

$K\left(x_{i},x_{j}\right)=\sigma^{2}\exp{\left[-\frac{1}{2}\left[x_{i}-x_{j}\right]^{T}W^{-1}\left[x_{i}-x_{j}\right]\right]}$
(53)

, where in most cases $W=diag\left[l_{1}^{2},\cdots,l_{D}^{2}\right]$.

We remind ourselves that we need to compute the following equations:

$\displaystyle m\left(\mu_{x^{*}},\Sigma_{x^{*}}\right)$
$\displaystyle=\sum\limits_{i=1}^{N}\left[yK_{ff}^{-1}\right]_{i,:}\Psi_{i}$
(54) $\displaystyle v\left(\mu_{x^{*}},\Sigma_{x^{*}}\right)$
$\displaystyle=\Psi-\sum\limits_{i,j=1}^{N}\left[K_{ij}^{-1}-\left[yK_{ff}^{-1}\right]_{i,:}\left[yK_{ff}^{-1}\right]_{j:}\right]\Psi_{ij}-m\left(\mu_{x^{*}},\Sigma_{x^{*}}\right)^{2}$
(55)

We now need to derive the analytic formulas for $\Psi$, $\Psi_{i}$ and
$\Psi_{ij}$. We will make use of the following identity.

######  Definition 2.1 (Product of Multivariate Gaussians).

For two multivariate Gaussian distributions, $\mathcal{N}\left(x\mid
a,A\right)$ and $\mathcal{N}\left(x\mid b,B\right)$, their product can be
expressed as follows:

$zN\left(x\mid c,C\right)=\mathcal{N}\left(x\mid
a,A\right)\mathcal{N}\left(x\mid b,B\right)$ (56)

where $c=C\left[A^{-1}a+B^{-1}b\right]$ and
$C=\left[A^{-1}+B^{-1}\right]^{-1}$

Armed with this toolkit we can proceed:

$\displaystyle\Psi$
$\displaystyle=\mathbb{E}_{x^{*}}\left[K\left(x^{*},x^{*}\right)\right]=\sigma^{2}$
(57) $\displaystyle\Psi_{i}$
$\displaystyle=\mathbb{E}_{x^{*}}\left[K\left(x^{*},x_{i}\right)\right]$ (58)
$\displaystyle=\int
K\left(x^{*},x_{i}\right)\mathcal{N}\left(x^{*}\mid\mu_{x^{*}},\Sigma_{x^{*}}\right)dx^{*}$
(59) $\displaystyle=c\int\mathcal{N}\left(x^{*}\mid
x^{i},W\right)\mathcal{N}\left(x^{*}\mid\mu_{x^{*}},\Sigma_{x^{*}}\right)dx^{*}$
(60) $\displaystyle=c\mathcal{N}\left(\mu_{x^{*}}\mid
x_{i},W+\Sigma_{x^{*}}\right)$ (61) $\displaystyle\Psi_{ij}$
$\displaystyle=\mathbb{E}_{x^{*}}\left[K\left(x^{*},x_{i}\right)K\left(x^{*},x_{j}\right)\right]$
(62) $\displaystyle=c^{2}\int\mathcal{N}\left(x^{*}\mid
x^{i},W\right)\mathcal{N}\left(x^{*}\mid
x^{j},W\right)\mathcal{N}\left(x^{*}\mid\mu_{x^{*}},\Sigma_{x^{*}}\right)dx^{*}$
(63) $\displaystyle=c^{2}\mathcal{N}\left(x_{i}\mid
x_{j},2W\right)\int\mathcal{N}\left(x^{*}\mid\frac{x_{i}+x_{j}}{2},\frac{W}{2}\right)\mathcal{N}\left(x^{*}\mid\mu_{x^{*}},\Sigma_{x^{*}}\right)dx^{*}$
(64) $\displaystyle=c^{2}\mathcal{N}\left(x_{i}\mid
x_{j},2W\right)\mathcal{N}\left(\mu_{x^{*}}\mid\frac{x_{i}+x_{j}}{2},\Sigma_{x^{*}}+\frac{W}{2}\right)$
(65)

where $c=\left(2\pi\right)^{\frac{D}{2}}\mid W\mid^{\frac{1}{2}}\sigma^{2}$.
We can now write the full explicit form for the two exact moments:

$\displaystyle m\left(\mu_{x^{*}},\Sigma_{x^{*}}\right)$
$\displaystyle=\sum\limits_{i=1}^{N}\left[yK_{ff}^{-1}\right]_{i,:}c\mathcal{N}\left(\mu(x^{*})\mid
x_{i},W+\Sigma_{x^{*}}\right)$ (66) $\displaystyle
v\left(\mu_{x^{*}},\Sigma_{x^{*}}\right)$
$\displaystyle=\sigma^{2}-c^{2}\sum\limits_{i,j=1}^{N}\left[yK_{ff}^{-1}\right]_{i,:}\left[yK_{ff}^{-1}\right]_{j,:}\mathcal{N}\left(x_{i}\mid
x_{j},2W\right)$ (67) $\displaystyle\hskip
28.45274pt\mathcal{N}\left(\mu_{x^{*}}\mid\frac{x_{i}+x_{j}}{2},\Sigma_{x^{*}}+\frac{W}{2}\right)-m\left(\mu_{x^{*}},\Sigma_{x^{*}}\right)^{2}$

We follow Quinonero-Candela et al. (2003) and re-write the explicit first two
moments as:

$\displaystyle m\left(\mu_{x^{*}},\Sigma_{x^{*}}\right)$
$\displaystyle=q^{T}K_{ff}^{-1}y$ (68) $\displaystyle
v\left(\mu_{x^{*}},\Sigma_{x^{*}}\right)$
$\displaystyle=\sigma^{2}+Tr\left[\left(K_{ff}^{-1}yy^{T}K_{ff}^{-T}-K_{ff}^{-1}\right)Q\right]-Tr\left[y^{T}K_{ff}^{-1}y\right]$
(69)

, where we have:

$\displaystyle q_{i}$ $\displaystyle=\mid
W^{-1}\Sigma_{x^{*}}+\mathbb{I}\mid^{-\frac{1}{2}}\exp{\left[-\frac{1}{2}\left(\mu_{x^{*}}-x_{i}\right)^{T}\left(\Sigma_{x^{*}}+W\right)^{-1}\left(\mu_{x^{*}}-x_{i}\right)\right]}$
(70) $\displaystyle Q_{ij}$ $\displaystyle=\mid
W^{-1}\Sigma_{x^{*}}+\mathbb{I}\mid^{-\frac{1}{2}}\exp\left[-\frac{1}{2}\left(\frac{x_{i}+x_{j}}{2}-\mu_{x^{*}}\right)^{T}\left(\Sigma_{x^{*}}+\frac{1}{2}W\right)^{-1}\left(\frac{x_{i}+x_{j}}{2}-\mu_{x^{*}}\right)\right]$
(71) $\displaystyle\hskip
28.45274pt\exp\left[-\frac{1}{2}\left(x_{i}-x_{j}\right)^{T}\left(2W\right)^{-1}\left(x_{i}-x_{j}\right)\right]$

###  2.4 Deep Gaussian Processes

We can view DGPs as a composition of functions $f_{L}(x)=f_{L}\circ...\circ
f_{1}(x)$ with
$f_{l}=\mathcal{GP}\left(m_{l},k_{l}\left(\cdot,\cdot\right)\right)$. Assuming
a general likelihood function, we can write the joint prior as:

$p\left(y,\\{f_{l}\\}_{l=1}^{L};X\right)=\underbrace{p(y|f_{L})}_{\text{likelihood}}\underbrace{\prod_{l=1}^{L}p(f_{l}|f_{l-1})}_{\text{prior}}$
(72)

with $p\left(f_{l}\mid
f_{l-1}\right)\sim\mathcal{GP}\left(m_{l}(f_{l-1}),k_{l}\left(f_{l-1},f_{l-1}\right)\right)$,
where we use squared exponential kernels
$k^{SE}_{l}(f_{l,i},f_{l,j})=\sigma^{2}_{l}\exp{\left[\sum_{d=1}^{D_{l}}-\frac{\left(f_{l,i,d}-f_{l,j,d}\right)^{2}}{l^{2}_{l,d}}\right]}$,
where $D_{l}$ represents the number of dimensions of $F_{l}$ and we introduce
layer specific kernel hyperparameters
$\theta_{l}=\\{\sigma^{2}_{l},l^{2}_{l,1},\cdots,l^{2}_{l,D_{l}}\\}$.

Analytically integrating this Bayesian hierarchical model is intractable as it
requires integrating Gaussians present in a non-linear way. To enable faster
inference over our model we can augment each layer $l$ with $M_{l}$ inducing
locations $Z_{l-1}$, respectively inducing values $U_{l}$ resulting in the
following augmented joint prior
$p\left(y,\\{f_{l}\\}_{l=1}^{L},\\{U_{l}\\}_{l=1}^{L};X,\\{Z_{l}\\}_{l=0}^{L-1}\right)$:

$\underbrace{p(y|f_{L})}_{\text{likelihood}}\underbrace{\prod_{l=1}^{L}p(f_{l}|f_{l-1},U_{l};Z_{l-1})p(U_{l})}_{\text{prior}}$
(73)

To perform stochastic variational inference (Hoffman et al., 2013) we
introduce an approximate posterior:
$q\left(\\{U_{l}\\}_{l=1}^{L}\right)=\prod\limits_{l=1}^{L}\mathcal{N}\left(U_{l}\mid
m_{U_{l}},S_{U_{l}}\right)$. Using a similar derivation as for sparse
variational GPs (SVGP), we can arrive at the evidence lower bound:

$\displaystyle\log{p(y)}$
$\displaystyle\geq\mathbb{E}_{q(\prod\limits_{l=1}^{L}q\left(f_{l}\mid
f_{l-1}\right))}\left[\log{p\left(y\mid f_{L}\right)}\right]-$ (74)
$\displaystyle\hskip
28.45274pt\sum\limits_{l=1}^{L}KL\left[q(U_{l})\|p(U_{l})\right]=\mathcal{L}_{DGP}$

, where $q\left(f_{l}\mid
f_{l-1}\right)=\mathcal{N}\left(f_{l}\mid\tilde{U_{l}}(f_{l-1}),\tilde{\Sigma_{l}}(f_{l-1})\right)$,
respectively:

$\displaystyle\tilde{U_{l}}(f_{l-1})$
$\displaystyle=m_{l}(f_{l-1})+K_{fu}K_{uu}^{-1}\left[m_{U_{l}}-m_{l}(Z_{l-1})\right]$
(75) $\displaystyle\tilde{\Sigma_{l}}(f_{l-1})$
$\displaystyle=K_{ff}-K_{fu}K_{uu}^{-1}\left[K_{uu}-S_{U_{l}}\right]K_{uu}^{-1}K_{uf}$
(76)

This composition of functions is approximated via Monte Carlo integration as
in Salimbeni and Deisenroth (2017).

###  2.5 Wasserstein-2 distance kernels on Gaussian measures

The Wasserstein space on $\mathbb{R}$ can be defined as the set
$W_{2}(\mathbb{R})$ of probability measures on $\mathbb{R}$ with a finite
moment of order two. We denote by $\Pi(\mu,\nu)$ the set of all probability
measures $\Pi$ over the product set $\mathbb{R}\times\mathbb{R}$ with
marginals $\mu$ and $\nu$, which are probability distributions in
$W_{2}(\mathbb{R})$. The transportation cost between two measures $\mu$ and
$\nu$ is defined as
$T_{2}(\mu,\nu)=\inf_{\pi\in\Pi(\mu,\nu)}\int[x-y]^{2}d\pi(x,y)$. This
transportation cost allows us to endow the set $W_{2}(\mathbb{R})$ with a
metric by defining the quadratic Wasserstein distance between $\mu$ and $\nu$
as $W_{2}(\mu,\nu)=T_{2}(\mu,\nu)^{1/2}$. More details on Wasserstein spaces
and their links with optimal transport problems can be found in Villani
(2008).

######  Theorem 2.1 (Theorem IV.1. in Bachoc et al. (2017)).

Let $k_{W}:W_{2}(\mathbb{R})\times W_{2}(\mathbb{R})\rightarrow\mathbb{R}$ be
the Wasserstein-2 RBF kernel defined as following:

$k^{W_{2}}(\mu,\nu)=\sigma^{2}\exp\frac{-W_{2}^{2}(\mu,\nu)}{l^{2}}$ (77)

then $k^{W_{2}}(\mu,\nu)$ is a positive definite kernel for any $\mu,\nu\in
W_{2}(\mathbb{R})$, respectively $\sigma^{2}$ is the kernel variance, $l^{2}$
being the lengthscale.

For completeness a proof following the layout in Thi Thien Trang et al. (2019)
is provided in Appendix LABEL:apd:pos_def_kernels.

Multiplication of positive definite kernels results again in a positive
definite kernel, hence we arrive at the automatic relevance determination
kernel based on Wasserstein-2 distances:

$k^{W_{2}}([\mu_{d}]_{d=1}^{D},[\nu_{d}]_{d=1}^{D})=\sigma^{2}\exp\sum_{d=1}^{D}\frac{-W_{2}^{2}(\mu_{d},\nu_{d})}{l_{d}^{2}}$
(78)

##### Wasserstein-2 Distance between Gaussian distributions

Gaussian measures fulfill the condition of finite second order moment, thereby
being a clear example of probability measures for which we can compute
Wasserstein metrics. The Wasserstein-2 distance between two multivariate
Gaussian distributions $\mathcal{N}(m_{1},\Sigma_{1})$ and
$\mathcal{N}(m_{2},\Sigma_{2})$, which have associated Gaussian measures and
implicitly the Wasserstein metric is well defined for them, has been shown to
have the following form $\parallel
m_{1}-m_{2}\parallel_{2}^{2}+\Tr\Big{[}\Sigma_{1}+\Sigma_{2}-2\Big{(}\Sigma_{1}^{1/2}\Sigma_{2}\Sigma_{1}^{1/2}\Big{)}^{1/2}\Big{]}$
(Dowson and Landau, 1982), which in the case of univariate Gaussians
simplifies to $|m_{1}-m_{2}|^{2}+|\sqrt{\Sigma_{1}}-\sqrt{\Sigma_{2}}|^{2}$.
This last formulation will be used throughout this paper.

###  2.6 Varifold theory or kernels between generalized measure spaces

Varifolds were first introduced in Geometric Measure Theory (GMT) (Almgren,
1966), in which shapes like submanifolds are represented in certain spaces of
generalized measures such as currents, varifolds or normal cycles. This
general framework was adapted to suit requirements for computational anatomy
matching in Charon and Trouvé (2013) and without loss of generality,
throughout the remainder of this subsection we will focus solely on surfaces
embedded in $\mathbb{R}^{3}$ and how to match them.

We commence by introducing elementary definitions from GMT. Subsequently
introduced methods are capable of working with smooth surfaces alongside
piecewise smooth or discrete objects. The umbrella term of rectifiable subsets
in GMT encompasses the aforementioned scenarios. We denote
$\mathbb{M}\in\mathbb{R}^{3}$ a rectifiable surface if for almost every
$x\in\mathbb{M}$, there exists a tangent space
$\mathbb{T}_{x}\mathbb{M}\in\mathbb{R}^{2}$.

To make the abstract mathematical framework more intuitive, we consider the
problem of partial matching between shapes (Antonsanti et al., 2021). Hence,
the aim is to find the optimal deformation to register a source $S$ to a
target shape $T$, where $S,T\in\mathbb{R}^{m}$, hence they are a finite union
of $m$-dimensional submanifolds of the ambient space $E=\mathbb{R}^{d}$. We
consider them to be m-rectifiable, which in the case of hyperspheres implies
that any point $x\in S$ ($x\in T$) has a normal vector $\tau_{x}S\in S^{d-1}$
($\tau_{x}T\in\mathbb{S}^{d-1}$).

######  Definition 2.2.

A varifold is a distribution on the product
$\mathbb{R}^{k}\times\mathbb{S}^{k-1}$, namely a continuous linear form on a
given space $W$ of smooth functions
$\mathbb{R}^{k}\times\mathbb{S}^{k-1}\to\mathbb{R}$. If $M$ is a surface, we
associate the varifold $\\{M\\}\in W^{{}^{\prime}}$, for any function $w\in
W$, by the surface integral:
$\\{M\\}(w)=\int_{M}w\left(x,\tau_{x}T\right)dS(x)$

######  Proposition 2.1 (Proposition 4.1. in Charon and Trouvé (2013)).

Given positive-definite real kernel $k_{e}$ on $\mathbb{R}^{d}$ s.t. $k_{e}$
is continuous and for all $x\in\mathbb{R}^{d}$, $k_{e}(x,\cdot)$ vanishes at
infinity, respectively a positive-definite real kernel $k_{t}$ defined on the
manifold $\mathbb{S}^{d-1}$ that is also continuous. Then the Reproducing
Kernel Hilbert Space (RKHS) $W$ associated to the positive-definite kernel
$k_{e}\otimes k_{t}$ is continuously embedded into the space
$C_{0}\left(\mathbb{R}^{d}\times\mathbb{S}^{d-1}\right)$, the space of
continuous functions on $\mathbb{R}^{d}\times\mathbb{S}^{d-1}$ that decay to 0
at infinity.

Its dual space $W^{{}^{\prime}}$ is the space of varifolds. The reproducing
kernel of $W$ is given by $k_{e}\otimes k_{t}$, with
$k_{e}(x,y)=\exp{\left(-\mid x-y\mid^{2}\right)}/\sigma_{w}^{2}$ and
$k_{t}(u,v)=\exp{\langle u,v\rangle_{R}^{d}}$ and we associate the following
canonical function $w_{S}\in W$:

$w_{S}(y,\tau)=\int_{S}k_{e}(y,x)k_{t}(\tau,\tau_{x}S)dx$ (79)

We can finally define the scalar product between canonical functions $w_{S}$
and $w_{T}$ as follows:

$\langle
w_{S},w_{T}\rangle=\int_{S}\int_{T}k_{e}(x,y)k_{t}(\tau_{x}S,\tau_{y}T)dx~{}dy$
(80)

######  Remark 2.1.

Intuitively, the inner product associated to RKHS $W$ keeps track of both
Euclidean differences between 3D points on the surface mesh, but also takes
into account geometric information surrounding those points, as given by the
points’ surface normal vectors (Figure 1).

Figure 1: Similarities between Varifold kernel and Hybrid kernel. Left:
Comparing two manifolds involves taking 3D points in Euclidean space (denoted
by circles) and their associated surface normal vectors (denoted by directed
arrows), subsequently computing their distance in their respective domains. We
only show the comparison between two points for aesthetical reasons, however
it should be noted that in practice a multitude of points and associated
surface normal vectors are used; Right: Considering two data points $S$ and
$T$, their predictive distribution at a certain hidden layer of a DDGP would
consist of mixture of Gaussians (p.d.f. of mixture components are depicted
either in pink or green). The hybrid kernel takes samples from mixture
components and their underlying Gaussian distributions.

##  3 Propagation of variance in DGP

Ustyuzhaninov et al. (2019) have proved for noiseless data that the various
stochastic layers in the hierarchy of a DGP collapse to deterministic
transformations. In this section, we analyze what are the necessary
theoretical requirements for a zero-mean DGP to maintain comparatively higher
total uncertainty for an OOD data point throughout the hierarchy and draw some
insights into what type of uncertainty (e.g., parametric/epistemic or non-
parametric/distributional) is propagate forward in each case. We analyze two
different case studies that might occur in practice. The first one involves
inducing points $\\{Z_{0},\cdots,Z_{L}\\}$ that are evenly spread around 0
with a maximum range of $\left[-3\sigma_{l},3\sigma_{l}\right]$ at each layer
such that samples stemming from the previous layers are within the
neighbourhood of said inducing points. Secondly, we will consider the scenario
where we have two separate groups of inducing points
$Z_{l,1}\subsetneq\left[-\infty,-c\right]$ and
$Z_{l,2}\subsetneq\left[c,\infty\right]$ for $c\in\mathbb{R^{+}}$. This
situation can occur in practice for a binary classification task, where at the
respective layer, the inducing points can reliably separate the two classes.
We now detail what are the necessary conditions in the first scenario to
predict higher total uncertainty for OOD data points in comparison to in-
distribution points.

Figure 2: Collapse of non-parametric variance in DGP. Layer-wise decomposition
of uncertainty into parametric and non-parametric for a zero-mean DGP,
alongside first order derivatives.

######  Proposition 3.1.

We consider the approximate posterior DGP as a composition of functions
$q_{L}=q_{L}\circ...\circ q_{1}$ with $q_{l}$ being given by
$q(F_{l}|F_{l-1})=\int
p(F_{l}|U_{l},F_{l-1})q(U_{l})=\mathcal{N}(F_{l}\mid\tilde{U_{l}}(F_{l-1}),\tilde{\Sigma_{l}}(F_{l-1}))$.
We consider an OOD data point $x_{out}$ in input space such that
$\tilde{\Sigma_{1}}(x_{out})=\sigma^{2}_{1}$ and $\tilde{U_{1}}(x_{out})=0$,
respectively $x_{in-d}$ to be in-distribution, with
$\tilde{\Sigma}_{1}(x_{in-d})=V_{in-d}\leq\sigma^{2}_{1}$ and
$\tilde{U}_{1}(x_{in-d})=M_{in-d}$. We assume that $Z_{2}$ are equidistantly
placed between $[-3\sigma_{1},3\sigma_{1}]$. The variance in the second layer
of $x_{out}$ will be higher than $x_{in-d}$ if the following holds:

$\frac{\left(\frac{\partial\tilde{U}_{2}(F_{1})}{\partial
F_{1}}\right)^{2}\Bigr{|}_{\begin{subarray}{c}F_{1}=M_{in-d}\end{subarray}}}{\left(\frac{\partial\tilde{U}_{2}(F_{1})}{\partial
F_{1}}\right)^{2}\Bigr{|}_{\begin{subarray}{c}F_{1}=0\end{subarray}}}\leq\frac{\sigma^{2}_{1}}{V_{in-d}}$
(81)

###### Proof.

We consider the approximate posterior DGP as a composition of functions
$q_{L}=q_{L}\circ...\circ q_{1}$ with a certain $q_{l}$ being given by:

$q(F_{l}|F_{l-1})=\int
p(F_{l}|U_{l},F_{l-1})q(U_{l})=\mathcal{N}(F|\tilde{U_{l}}(f_{l-1}),\tilde{\Sigma_{l}}(f_{l-1}))$
(82)

Following the layout introduced in Girard et al. (2002) for obtaining Gaussian
approximations of GP with uncertain inputs, we have the following adapted case
for a two layer DGP, with the forthcoming analysis being easily extended to
more layers: U

$q(F_{2})(x)=\int p(F_{2}|F_{1})q(F_{1}(x))~{}dF_{1}$ (83)

Using the framework described in subsection 2.3.1 (equations (50) and (52))
and adapting to our case at hand, we have the following approximate moments
for $q(F_{2})(x)$:

$\displaystyle m(F_{2}(x))$ $\displaystyle=\tilde{U}_{2}(\tilde{U}_{1}(x))$
(84) $\displaystyle v(F_{2}(x))$
$\displaystyle=\tilde{\Sigma}_{2}(\tilde{U}_{1}(x))+\tilde{\Sigma}_{1}(x)\Bigg{[}\frac{1}{2}\frac{\partial^{2}\tilde{\Sigma}_{2}(F_{1})}{\partial^{2}F_{1}}\Bigr{|}_{\begin{subarray}{c}F_{1}=\tilde{U}_{1}(x)\end{subarray}}+\left(\frac{\partial\tilde{U}_{2}(F_{1})}{\partial
F_{1}}\right)^{2}\Bigr{|}_{\begin{subarray}{c}F_{1}=\tilde{U}_{1}(x)\end{subarray}}\Bigg{]}$
(85)

, where for the purposes of our derivation here we only considered a first
order Taylor expansion for the approximate mean. In Ustyuzhaninov et al.
(2019), the authors have empirically shown that
$\lim_{M\to\infty}\frac{\partial^{2}\tilde{\Sigma_{2}}(F_{1})}{\partial^{2}F_{1}}=0$,
hence that term is dropped. Intuitively, it follows a sinusoidal path with the
amplitude converging to 0 as the number of inducing points is increased.

We remind ourselves that in the proposition statement we made the following
set of assumptions: inducing point locations for the second layer $Z_{1}$ are
equidistantly placed in the interval $[-3\sigma_{1},3\sigma_{1}]$; the
variance of the inducing point values $S_{U_{2}}$ are taken to be equal to
simplify the problem.

For $x_{out}$ and $x_{in-d}$ using equation 85 we obtain:

$\displaystyle v(F_{2}(x_{out}))$
$\displaystyle=\tilde{\Sigma}_{2}(0)+\sigma^{2}_{1}\left(\frac{\partial\tilde{U}_{2}(F_{1})}{\partial
F_{1}}\right)^{2}\Bigr{|}_{\begin{subarray}{c}F_{1}=0\end{subarray}}$ (86)
$\displaystyle v(F_{2}(x_{in-d}))$
$\displaystyle=\tilde{\Sigma}_{2}(M_{in-d})+V_{in-d}\left(\frac{\partial\tilde{U}_{2}(F_{1})}{\partial
F_{1}}\right)^{2}\Bigr{|}_{\begin{subarray}{c}F_{1}=M_{in-d}\end{subarray}}$
(87)

The desired behaviour for the second layer is:

$\tilde{\Sigma}_{2}(M_{in-d})+V_{in-d}\left(\frac{\partial\tilde{U}_{2}(x)}{\partial
x}\right)^{2}\Bigr{|}_{\begin{subarray}{c}F_{1}=M_{in-d}\end{subarray}}\leq\tilde{\Sigma}_{2}(0)+\sigma^{2}_{1}\left(\frac{\partial\tilde{U}_{2}(x)}{\partial
x}\right)^{2}\Bigr{|}_{\begin{subarray}{c}F_{1}=0\end{subarray}}$ (88)

Since we have assumed that the inducing point locations are equidistantly
placed in a grid with centre zero, which is a sufficiently similar situation
to what one might observe in practice, then
$\tilde{\Sigma}_{2}(0)\approx\tilde{\Sigma}_{2}(M_{in-d})$ as we have also
assumed that the posterior variances of the inducing point values are equal.
Hence, the inequality $v(F_{2}(x_{in-d}))\leq v(F_{2}(x_{out}))$ holds if and
only if:

$\frac{\left(\frac{\partial\tilde{U}_{2}(F_{1})}{\partial
F_{1}}\right)^{2}\Bigr{|}_{\begin{subarray}{c}F_{1}=M_{in-d}\end{subarray}}}{\left(\frac{\partial\tilde{U}_{2}(F_{1})}{\partial
F_{1}}\right)^{2}\Bigr{|}_{\begin{subarray}{c}F_{1}=0\end{subarray}}}\leq\frac{\sigma^{2}_{1}}{V_{in-d}}$
(89)

∎

######  Remark 3.1.

Intuitively, this inequality can be easily satisfied if
$\left(\frac{\partial\tilde{U}_{2}(F_{1})}{\partial
F_{1}}\right)^{2}\Bigr{|}_{\begin{subarray}{c}F_{1}=0\end{subarray}}$ reaches
its absolute peak at $F_{1}=0$ since $1\leq\frac{\sigma^{2}_{1}}{V_{in-d}}$
holds in general. However, this requirement is not guaranteed to be satisfied
in practice. On a toy regression task we can observe that the overall
uncertainty is correctly predicted as higher for OOD points due to high
absolute values of the derivative at 0 (i.e., where OOD points get mapped in
hidden layers on average) in hidden layers’ space (Figure 2).

In the first investigated case, we have thus established the necessary
condition for the reliable propagation of total variance so as to guarantee
that OOD data points have higher total uncertainty compared to in-distribution
data points at each layer in the hierarchy. However, a reasonable question is
whether DGPs in this case are _distance-aware_ (see Definition 4.1) at each
layer? From Figure 2 we can see that the non-parametric uncertainty (which
stands as a proxy for distributional uncertainty; the two terms will be used
interchangeably throughout the remainder of the paper) collapses to zero
across the whole spectrum of data points in the hidden layers. Therefore, DGPs
in this scenario are not _distance-aware_ , effectively rendering the non-
parametric part of layer-specific SVGP obsolete, thus resulting in a fully
parametric DGP from the first hidden layer up to the output layer (for in-
depth insight into why this pathology occurs we refer the reader to Figure 6
and accompanying text). The second scenario for DGPs will investigate the case
where the non-parametric uncertainty for an OOD data point does get correctly
propagated throughout the layers, thus making use of the non-parametric part
of SVGP in each layer of the hierarchy. We lay out the necessary conditions in
the following proposition. Before we proceed, we remind ourselves that a SVGP
can be decomposed as
$F(\cdot)\sim\underbrace{K_{fu}K_{uu}^{-1}U}_{G(\cdot)}+\underbrace{\mathcal{N}\left(0,K_{ff\cdot
u}\right)}_{H(\cdot)}$.

######  Proposition 3.2.

Under the assumption that for any OOD point in input space $x_{out}$ we have
$G_{l}(x_{out})=K_{fu}K_{uu}^{-1}U_{l}\approx 0$ at each layer in a DGP, we
can instead focus solely on the non-parametric components of DGP as a
composition of $p(H_{L})=p(H_{L}\mid H_{L-1})\circ...\circ p(H_{1}\mid H_{0})$
with $p\left(H_{l}\mid H_{l-1}\right)=\mathcal{N}\left(H_{l}\mid 0,K_{ff\cdot
u}\right)$ and $H_{0}=X$. We consider two inducing points in the first hidden
layer $\\{Z_{1,1}=-c,Z_{1,2}=c\\}$ with $c\geq 0$ and an OOD data point
$x_{out}$ in input space such that
$p(H_{1}(x_{out}))\sim\mathcal{N}\left(0,\sigma^{2}_{1}\right)$. The
approximate variance of the non-parametric component pertaining to the second
layer will satisfy $v\left(H_{2}(x_{out})\right)\geq K_{ff\cdot u}$ if the
following holds $l^{2}_{2}\geq 2c^{2}$, where the Schur complement is
evaluated at $F_{1}=0$.

###### Proof.

We commence by writing our Bayesian hierarchical formulation so that the link
with approximation of Gaussian Processes’ moments with noisy inputs (see
section 2.3.1) is made clearer. Under our set of assumptions we have:

$\displaystyle p(H_{1}\mid H_{0})$ $\displaystyle=\mathcal{N}\left(H_{1}\mid
0,K_{ff}-K_{fu}K_{uu}^{-1}K_{uf}\right)$ (90)
$\displaystyle=\mathcal{N}\left(H_{1}\mid 0,\sigma^{2}_{1}\right)$ (91)
$\displaystyle p(H_{2}\mid H_{1})$
$\displaystyle=\int\mathcal{N}\left(H_{2}\mid
0,K_{ff}-K_{fu}K_{uu}^{-1}K_{uf}\right)\mathcal{N}\left(H_{1}\mid
0,\sigma^{2}_{1}\right)~{}dH_{1}$ (92)

We can approximate the non-Gaussian predictive equation using the framework of
exact moments introduced in subsection 2.3.2 since we are only considering
squared exponential kernels. We adapt equation (69) to suit our scenario:

$\displaystyle
v\left(0,\sigma^{2}_{1}\right)=\sigma^{2}_{2}-Tr\left[K_{uu}^{-1}Q\right]$
(93)

where we have:

$\displaystyle K$
$\displaystyle=\begin{pmatrix}\sigma^{2}_{2}&\sigma^{2}_{2}\exp{\left[-\frac{2c^{2}}{l^{2}}\right]}\\\
\sigma^{2}_{2}\exp{\left[-\frac{2c^{2}}{l^{2}}\right]}&\sigma^{2}_{2}\end{pmatrix}$
(94) $\displaystyle Q$
$\displaystyle=\frac{1}{\sqrt{\frac{2\sigma^{2}_{1}}{l^{2}}+1}}\begin{pmatrix}\exp{\left[-\frac{c^{2}}{l^{2}+2\sigma^{2}_{1}}\right]}&\exp{\left[-\frac{c^{2}}{l_{2}}\right]}\\\
\exp{\left[-\frac{c^{2}}{l_{2}}\right]}&\exp{\left[-\frac{c^{2}}{l^{2}+2\sigma^{2}_{1}}\right]}\end{pmatrix}$
(95)

We can now write the exact variance of the Gaussian approximation of
$p(H_{2})$ as:

$v\left(0,\sigma^{2}_{1}\right)=\sigma^{2}_{2}-\frac{2\left[\exp{\left[-\frac{c^{2}}{l^{2}+2\sigma^{2}_{1}}\right]}-\exp{\left[-\frac{3c^{2}}{l^{2}}\right]}\right]}{\sigma^{2}_{2}\left[1-\exp{\left(-\frac{4c^{2}}{l^{2}}\right)}\right]\sqrt{\frac{2\sigma^{2}_{1}}{l_{2}}+1}}$
(96)

which has the following derivative with respect to $\sigma^{2}_{1}$:

$\frac{\partial
v\left(0,\sigma^{2}_{1}\right)}{\partial\sigma^{2}_{1}}=\sigma^{2}_{2}\left[1-\exp{\left(-\frac{4c^{2}}{l^{2}}\right)}\right]2\exp{\left[-\frac{c^{2}}{l^{2}+2\sigma^{2}_{1}}\right]}\left[\frac{2\sigma^{2}_{1}+l^{2}}{l^{2}}\right]^{-\frac{1}{2}}\left[2c^{2}\left[l^{2}+2\sigma^{2}_{l-1}\right]-\frac{1}{2\sigma^{2}_{1}+l^{2}}\right]$
(97)

We first consider the derivative evaluated at $\sigma^{2}_{1}=0$:

$\frac{\partial
v\left(0,\sigma^{2}_{1}\right)}{\partial\sigma^{2}_{1}}\Bigr{|}_{\begin{subarray}{c}\sigma^{2}_{1}=0\end{subarray}}=\sigma^{2}_{2}\left[1-\exp{\left(-\frac{4c^{2}}{l^{2}}\right)}\right]\left[-\frac{c^{2}}{l^{2}}\right]\left[\frac{2c^{2}}{l^{4}}-\frac{1}{l^{2}}\right]$
(98)

The first two terms are always positive in the case that $K_{uu}$ is
invertible, leaving us with just the last term. Hence, we have that:

$\frac{\partial
v\left(0,\sigma^{2}_{1}\right)}{\partial\sigma^{2}_{1}}\Bigr{|}_{\begin{subarray}{c}\sigma^{2}_{1}=0\end{subarray}}\geq
0\iff l^{2}\geq 2c^{2}$ (99)

Similarly, one can investigate whether the derivative is positive at any other
value except 0:

$\frac{\partial
v\left(0,\sigma^{2}_{1}\right)}{\partial\sigma^{2}_{1}}\Bigr{|}_{\begin{subarray}{c}\sigma^{2}_{1}=\tilde{\sigma^{2}_{1}}\end{subarray}}\geq
0\iff\tilde{\sigma^{2}_{1}}+l^{2}\geq 2c^{2}$ (100)

One can easily see that at any value of $\sigma^{2}_{1}$ the derivative will
be positive only if $l^{2}\geq 2c^{2}$. Hence, the minimum values will be at
$v(0,0)=K_{ff}-Q_{ff}$, where the Schur complement is evaluated at point
$F_{1}=0$. ∎

######  Remark 3.2.

Reliable propagation of non-parametric uncertainty for OOD data points in DGP
can only occur if the parametric posterior at each layer
$\mathcal{N}\left(G_{l}\mid
K_{fu}K_{uu}^{-1}m_{U_{l}},K_{fu}K_{uu}^{-1}S_{U_{l}}K_{uu}^{-1}K_{uf}\right)$
manages to reliably map in-distribution points stemming from the previous
layer in the intervals $I_{in-d}^{1}=[-\infty,-3\sigma_{l-1}]$, respectively
$I_{in-d}^{2}=[3\sigma_{l-1},\infty]$. This scenario can occur for example in
simple binary classification tasks, where the parametric posteriors at each
layer have learned to reliably separate the two classes in hidden layer
features. Considering that OOD data points from the previous layer have a
predictive distribution $\mathcal{N}\left(0,\sigma^{2}_{l-1}\right)$, upon
sampling $95\%$ of samples from this distribution will be in the interval
$I_{OOD}=[-3\sigma_{l-1},3\sigma_{l-1}]$. Under an assumption of well-
optimized inducing points’ locations, $Z_{l-1}$ will invariably find
themselves in either $I_{in-d}^{1}$ or $I_{in-d}^{2}$, with no locations in
$I_{OOD}$. Using notation from Proposition 3.2, $c=3\sigma_{l-1}$ and the
approximate variance $v\left(0,\sigma^{2}_{l-1}\right)\geq K_{ff\cdot u}$
where the Schur complement is evaluated at $F_{l-1}=0$ only if $l_{l}^{2}\geq
2c^{2}$.

To summarize our results so far, in this section we have shown the necessary
conditions for DGPs to maintain a comparatively higher total uncertainty for
OOD points in comparison to in-distribution points. Through empirical
examples, we have seen that the non-parametric variance collapses to zero even
in toy regression datasets, thus making the DGP rely solely on its parametric
components. Next, we have proposed a scenario in which the posterior DGP can
reliably propagate non-parametric variance for OOD points. However, this
latter scenario is not guaranteed to occur in practice. In the next section,
we propose a variant of DGPs which are guaranteed to propagate non-parametric
variance for OOD points.

##  4 Distributional Deep Gaussian Processes

In this section we commence with a motivation behind the hybrid kernel
(Euclidean $\&$ Wasserstein space), subsequently introducing the generative
process and inference framework for our newly introduced model. Moreover, we
derive the theoretical requirements of correctly propagating outliers through
the hierarchy of this probabilistic formulation and show that it’s easier
satisfied in comparison to DGPs.

###  4.1 Desiderata for kernels in DGP

##### I Maintenance of outlier status.

While there are many definitions of what an outlier constitutes (Ruff et al.,
2021), for the purposes of this work we follow the most basic one,
respectively “An anomaly is an observation that deviates considerably from
some concept of normality.” More concretely, it can be formalised as follows:
our data resides in $X\in\mathbb{R}^{D}$, an anomaly/outlier is a data point
$x\in X$ that lies in a low probability region under $\mathcal{P}$ such that
the set of anomalies/outliers is defined as $A=\\{x\in X\mid
p(x)\leq\xi\\},~{}\xi\geq 0$, with $\xi$ is a threshold under which we
consider data points to deviate sufficiently from what normality constitutes.
GP constitute an easy plug-in method to define what normality constitutes, or
any other kernel-based algorithm for that matter, due to the mathematical
formalism embedded in the kernel.

######  Definition 4.1 (Definition 1 in Liu et al. (2020)).

We consider the predictive distribution for unseen point $p\left(y^{*}\mid
x^{*}\right)$ at testing time, for model trained on
$\mathbb{X}_{in-d}\in\mathbb{X}$, with the data manifold being equipped with
metric $\|\cdot\|_{X}$. Then, we can affirm that $p\left(y^{*}\mid
x^{*}\right)$ is _distance-aware_ if there exists a summary statistic
$u\left(x^{*}\right)$ of $p\left(y^{*}\mid x^{*}\right)$ that embeds the
distance between $\mathbb{X}_{in-d}$ and $x^{*}$:

$u\left(x^{*}\right)=v\left[\mathbb{E}_{x\sim\mathbb{X}_{in-d}}\left[\|x^{*}-x\|_{X}^{2}\right]\right]$
(101)

where $v$ is a monotonic function that increases with distance.

##### GP predictive variance as distributional uncertainty.

GPs are clearly _distance-aware_ provided we use a translation-invariant
kernel. The summary statistics for an unseen point is given by
$u(x^{*})=K_{f^{*}f^{*}}-K_{f^{*}f}K_{ff}^{-1}K_{ff^{*}}$, which is
monotonically increasing as a function of distance. Throughout this paper, we
will use the non-parametric variance of sparse GPs as a proxy for
distributional uncertainty, which will be used to assess if inputs are in or
outside the distribution.

Definition 4.1 does not make any assumptions related to the architecture of
the model from which the predictive distribution stems. In practice we would
have the following composition to arrive at the logits $logit(x^{*})=f\circ
enc\left(x^{*}\right)$, where $enc\left(\cdot\right)$ represents a network
that outputs the representation learning layer and $f\left(\cdot\right)$ is
the output layer. In Liu et al. (2020) the authors propose the following two
conditions to ensure that the composition is _distance-aware_ :

  * •

$f(\cdot)$ is _distance-aware_

  * •

$\mathbb{E}_{x\sim\mathbb{X}_{in-d}}\left[\|x^{*}-X\|_{X}^{2}\right]\approx\mathbb{E}_{x\sim\mathbb{X}_{in-d}}\left[\|enc(x^{*})-enc(X)\|_{enc(X)}^{2}\right]$

The last condition means that distances between data points in input space
should be correlated with distances in the representation learning layer,
which is equipped with a $\|\cdot\|_{enc(X)}$ metric. Whereas we have seen
that for example a Gaussian Process satisfies the _distance-aware_ condition
for the last layer predictor, we are still left with the question on how to
maintain distances in the representation learning layer correlated to
distances in the input layer.

As we have previously stated, GPs are _distance-aware_. Thus, they can
reliably notice departures from the training set manifold. For SVGP we
decompose the model uncertainty into two components:

$\displaystyle h(\cdot)$
$\displaystyle=\mathcal{N}(h|0,K_{ff}-K_{fu}K_{uu}^{-1}K_{uf})$ (102)
$\displaystyle g(\cdot)$
$\displaystyle=\mathcal{N}(g|K_{fu}K_{uu}^{-1}m_{U},K_{fu}K_{uu}^{-1}S_{U}K_{uu}^{-1}K_{uf})$
(103)

The $h(\cdot)$ variance captures the shift from within to outside the data
manifold and will be denoted as _distributional uncertainty_ (interchangeably
also denoted as non-parametric variance). The variance $g(\cdot)$ is coined as
_epistemic uncertainty_ (interchangeably also denoted as parametric variance
or within-data uncertainty) encapsulating uncertainty present inside the data
manifold. To capture the overall uncertainty in $h(\cdot)$, thereby also
capturing the spread of samples from it, we can calculate its differential
entropy as:

$h(h)=\frac{n}{2}\log{2\pi}+\frac{1}{2}\log{\mid
K_{ff}-K_{fu}K_{uu}^{-1}K_{uf}\mid}+\frac{1}{2}n$ (104)

In practice we only use the diagonal terms of the Schur complement, hence the
log determinant term will considerably simplify. Intuitively, if terms on the
diagonal of the Schur complement have higher values, so will the
distributional differential entropy. This OOD measure in logit space will be
used throughout the rest of the paper for measuring departures from training
set manifold. Distributional uncertainty as defined in equation (102) can be
used as proxy for OOD detection due to its _distance-awareness_ properties.
Therefore we desire DGP with similar OOD detection properties as it shallow
counterpart at every layer in its hierarchy.

##### II Smoothness in hidden layers.

Imposing smoothness in neural networks’ hidden layers has been linked to
increased generalization, robustness to input perturbations or reliable
uncertainty estimates (Rosca et al., 2020). With the latter property in mind,
recent work (van Amersfoort et al., 2021; Liu et al., 2020) have proposed bi-
Lipschitz constraints for extracting features to be fed into a GP. This
paradigm is based on the notion that objects which were close in previous
layers should also be close going forward in the hierarchy, thus avoiding
_feature-collapse_ (van Amersfoort et al., 2021). Consequently, kernels in DGP
should provide a certain guarantee of smoothness, respectively it should guard
against _feature-collapse_.

###  4.2 Hybrid kernel

The term _Distributional Gaussian Process_ was first introduced in (Bachoc et
al., 2017) to describe a shallow GP that operates on measures using a
Wasserstein-2 based kernel ( equation (78)) and for purposes of
differentiating it from subsequent models it will be coined as _DistGP-Bahoc_
throughout the remainder of the paper. We first investigate prior function
space properties if the latter formulation is extended to the multi-layered
case and then introduce our proposed kernel and generative process.

Considering the case of Euclidean input data, the fist layer is governed by a
GP $p(F_{1})=\mathcal{N}(0,K_{F_{0}F_{0}})$. We immediately see that the prior
has the capability to capture correlations between data points. Adding a new
layer given by _DistGP-Bahoc_ operating directly on $p(F_{1})$, where we only
take diagonal terms of its covariance matrix to compute
$K^{W_{2}}_{F_{1,i},F_{1,j}}=\sigma^{2}_{2}\exp\frac{-W_{2}^{2}(\mathcal{N}(0,\sigma^{2}_{1}),\mathcal{N}(0,\sigma^{2}_{1}))}{l^{2}_{2}}=\sigma^{2}_{2}$.
Therefore, correlations between different data points are ignored as we arrive
at the prior on $p\left(F_{2}\right)$:

$p(F_{2})=\mathcal{N}\left[\begin{pmatrix}0\\\ \vdots\\\
0\end{pmatrix},\begin{pmatrix}\sigma^{2}_{2}&\dots&\sigma^{2}_{2}\\\
\vdots&\ddots&\vdots\\\
\sigma^{2}_{2}&\dots&\sigma^{2}_{2}\end{pmatrix}\right]$ (105)

Consequently, this type of construction will exhibit an overly-correlated
prior which will cause samples to collapse to set of finite values, thereby
not satisfying Desiderata II. However, we desire a generative process that is
without pathologies in the zero-mean function case.

Sampling hidden layers’ features in DGP is of utmost importance towards
introducing correlations between data points. On the other hand, the
Wasserstein-2 kernel (equation (78)) directly applied on hidden layers’
resulting multivariate Gaussians is attractive as it takes into account
differences between variance terms. The hindsight is that differences between
low and high variance terms in the multivariate normals would enable a finer
separation of data points which were OOD compared to in-distribution as
detected by the previous layer. With these observations in mind, we argue that
both stochastic (sampling hidden layer features in DSVI framework for DGPs)
and deterministic versions (Wasserstein-2 kernels directly applied on hidden
layer GP) are important from entirely complementary viewpoints. We seek to
combine these two formulations in a single one that preservers both key
properties, we can now write the generative process of this new probabilistic
framework coined Distributional Deep Gaussian Processes (DDGP) for 2 layers:

$\displaystyle F_{1}$ $\displaystyle\sim\mathcal{N}\left(0,K_{ff}\right)$
(106) $\displaystyle F_{1}^{sth}$
$\displaystyle=m\left(F_{1}\right)+\sqrt{v\left(F_{1}\right)}\epsilon,~{}\epsilon\sim\mathcal{N}\left(0,\mathbb{I}_{n}\right)$
(107) $\displaystyle F_{1}^{det}$
$\displaystyle=\mathcal{N}\left(m(F_{1}),diag\left[v(F_{1})\right]\right)$
(108) $\displaystyle F_{2}$
$\displaystyle\sim\mathcal{N}\left(0,k_{hybrid}\left(\\{F_{1}^{sth},F_{1}^{det}\\},\\{F_{1}^{sth},F_{1}^{det}\\}\right)\right)$
(109)

More intuitively, this generative process implies keeping track of a
_stochastic_ , respectively _deterministic_ component of the same SVGP at any
given hidden layer. In Figure 3 we provide a schematic of the differences in
the generative process of DGPs and DDGPs.

Figure 3: Conceptual differences between DGPs and DDGPs. DGPs operate in
Euclidean space of sampled hidden layer features. DDGPs also operate on
Wasserstein space of underlying hidden layer GPs, with inducing locations
following a Gaussian distribution.

We took inspiration from Varifold theory (see subsection 2.6 for details) and
devised a kernel operating on both Euclidean and Wasserstein-2 space, thereby
incorporating both stochastic and deterministic pathways. Considering again
two data points in the prior space of $F_{1}$, respectively
$p\begin{pmatrix}F_{1,i}\\\
F_{1,j}\end{pmatrix}\sim\mathcal{N}\left[\begin{pmatrix}0\\\
0\end{pmatrix},\begin{pmatrix}\sigma_{1}^{2}&K^{E}_{i,j}\\\
K^{E}_{j,i}&\sigma_{1}^{2}\end{pmatrix}\right]$ and their associated samples
$\\{f_{1,i},f_{1,j}\\}$. To compute Wasserstein-2 distances using the previous
layer’s MVN, we take the diagonal of the covariance matrix and we denote
$\mu_{i}=\mathcal{N}(0,\sigma_{1}^{2})$;
$\mu_{j}=\mathcal{N}(0,\sigma_{1}^{2})$. We can now introduce the hybrid
kernel:

$k^{H}\left(\mu_{i},\mu_{i}\right)=k^{SE}(x_{i},x_{j})\exp\sum_{d=1}^{D}\frac{-W_{2}^{2}(\mu_{i,d},\mu_{j,d})}{l_{d}^{2}}$
(110)

Prior samples of $F_{2}$ will be correlated due to the Euclidean part of this
kernel, whereas the Wasserstein-2 component will always be equal to 1, hence
neither introducing or removing extra correlations between data points. To
make this clear, we can explicitly calculate it:

$\displaystyle p\begin{pmatrix}F_{2,i}\\\ F_{2,j}\end{pmatrix}$
$\displaystyle\sim\mathcal{N}\left[\begin{pmatrix}0\\\
0\end{pmatrix},\begin{pmatrix}\sigma_{2}^{2}\exp-\frac{-W_{2}^{2}(\mu_{i},\mu_{i})}{l^{2}}&K^{E}_{i,j}\exp-\frac{-W_{2}^{2}(\mu_{i},\mu_{j})}{l^{2}}\\\
K^{E}_{j,i}\exp-\frac{-W_{2}^{2}(\mu_{j},\mu_{i})}{l^{2}}&\sigma_{2}^{2}\exp-\frac{-W_{2}^{2}(\mu_{j},\mu_{j})}{l^{2}}\end{pmatrix}\right]$
(111) $\displaystyle\sim\mathcal{N}\left[\begin{pmatrix}0\\\
0\end{pmatrix},\begin{pmatrix}\sigma_{2}^{2}&K^{E}_{i,j}\\\
K^{E}_{j,i}&\sigma_{2}^{2}\end{pmatrix}\right]$ (112)

since the Wasserstein-2 distance between two equal distributions ($\mu_{i}$
and $\mu_{j}$) is zero.

######  Remark 4.1.

In prior space, a hierarchical construction based on hybrid kernels (equation
(110)) for hidden layers is equivalent to a DGP. In posterior space, this no
longer holds true with several different properties which we will outline
subsequently.

##### Connections between Varifold kernel and Hybrid Kernel.

Whereas usually during training we only sample once from $q(F_{l})$ at each
layer in order to obtain the samples that get used in the Euclidean part of
the kernel, one could sample arbitrary times, thereby obtaining a collection
of Gaussian distributions. In this case, we can re-express our hybrid kernel
in the multi-sample scenario:

$\sum_{i=1}^{S}\sum_{j=1}^{S}k^{SE}(x_{i},y_{j})\exp\sum_{d=1}^{D}\frac{-W_{2}^{2}(\mu_{i,d},\mu_{j,d})}{l_{d}^{2}}$
(113)

We can observe similarities with the varifold kernel in equation (80), where
we compute a distance between submanifolds, taking into consideration
Euclidean distances between discrete points on the manifolds’ surface,
combined with information pertaining to the points’ geometry. For hybrid
kernels as applied to DDGPs, we compare two sets of multivariate Gaussians,
taking into account Euclidean differences between points sampled from the
underlying distributions, but also considering geometric differences between
these distributions as captures by the Wasserstein-2 distance. A visual
depiction of these different but similar objectives is provided in figure 1.

Figure 4: Conceptual difference between euclidean and hybrid kernel. Explicit
difference between moments in Wasserstein-2 distance ensures _distance-
awareness_.

##### Hybrid kernels and OOD detection.

The key motivation behind the introduction of hybrid kernels can be succinctly
described by a simple toy example. Considering two-layered DGPs and DDGPs, we
assume an in-distribution point $x_{in-d}$ to have low total variance in
hidden layer $F_{1}$, respectively an OOD point $x_{ood}$ to have high total
variance. In the DGP case, upon sampling from $q(F_{1}(x_{in-d}))$ and
$q(F_{1}(x_{ood}))$ we can end up with samples which are equally distant with
respect to inducing points’ location $Z_{1}$. If this occurs, then non-
parametric variance (proxy for distributional uncertainty) will be equal for
$x_{in-d}$ and $x_{ood}$ in $F_{2}$. Hence, what was initially flagged as OOD
in the first hidden layer will be considered as in-distribution by the second
hidden layer, thus not satisfying Desiderata I. Conversely, in the DDGP case
and under the assumption that the variance of distributional inducing points’
locations $Z_{2}$ is almost equal in distribution to the total variance of in-
distribution points in $F_{1}$, the Wasserstein-2 component of the hybrid
kernel will notice that there is a higher distance between the now
distributional inducing point location and $x_{ood}$, as opposed to the
distance of the former with $x_{in-d}$. Then, the non-parametric variance of
$x_{ood}$ will be higher than that of $x_{in-d}$, thus satisfying Desiderata
I. A visual depiction of this case study is illustrated in Figure 4.

###  4.3 Evidence lower bound

In this subsection we outline the generative process of constructing multi-
layered architectures using the hybrid kernel, which we will coin as
Distributional Deep Gaussian Processes (DDGP). First, we remind ourselves that
a DGP is defined as a stack of shallow GP acting as the prior:

$\underbrace{p(Y|F_{L})}_{\text{likelihood}}\underbrace{\prod_{l=1}^{L}p(F_{l}|F_{l-1},U_{l};Z_{l-1})p(U_{l})}_{\text{prior}}$
(114)

, where for brevity of notation we denote $F_{0}=X$. The SVGPs between hidden
layers are treated as being noiseless. As the prior is analytically
intractable to integrate, Salimbeni and Deisenroth (2017) have suggested to
sample from each hidden layer of the DGP in order to obtain unbiased
stochastic gradients. The difference between DGP and DDGP resides in the fact
that for $l\geq 2$ we not only sample from the marginal of the data points
$q(F_{l};F_{l-1})$, but we also explicitly use this marginal in the
computation of the next layer’s kernel. More concretely, we have the following
joint density prior $p(Y,\\{F_{l},U_{l}\\}_{l=1}^{L})$:

$\underbrace{p(Y|F_{L})}_{likelihood}\underbrace{\prod_{l=2}^{L}p(F_{l}|F_{l-1},U_{l};Z_{l-1})p(U_{l})}_{\text{Wasserstein
\& Euclidean space}}\underbrace{p(F_{1}|U_{1};X)p(U_{1})}_{\text{Euclidean
space}}$ (115)

We introduce a factorized posterior between layers and dimensions of the
following form : $q(\\{F_{l}\\}_{l=1}^{L},\\\
\\{U_{l}\\}_{l=1}^{L})=p(F_{L}|U_{L};Z_{L-1})\prod_{l=1}^{L}q(U_{l})$, where
$q(U_{l})$ is taken to be a multivariate Gaussian with mean $m_{U_{l}}$ and
variance $S_{U_{l}}$. The posterior over $q(F_{1})$ is similar to the DGP
case, with the difference residing in $l\geq 2$, where we have the following
conditional DistGP:

$p(F_{l}|U_{l};Z_{l-1},F_{l-1})=\mathcal{N}(F_{l}|K_{fu}^{H}K^{H}_{uu}\raisebox{4.95134pt}{$\scriptscriptstyle-\\!1$}U_{l},K_{ff}^{H}-Q_{ff}^{H})$
(116)

After sampling from $q(F_{l-1})$ and
$Z_{l-1}\sim\mathbf{N}(M_{Z_{l-1}},\Sigma_{Z_{l-1}})$, we proceed to integrate
out $U_{l}$, arriving at $q(F_{l})$:

$\mathcal{N}(F_{l}|K_{fu}^{H}K^{H}_{uu}\raisebox{4.95134pt}{$\scriptscriptstyle-\\!1$}m_{U_{l}},K_{ff}^{H}-Q_{ff}^{H}+K_{fu}^{H}K^{H}_{uu}\raisebox{4.95134pt}{$\scriptscriptstyle-\\!1$}S_{U_{l}}K^{H}_{uu}\raisebox{4.95134pt}{$\scriptscriptstyle-\\!1$}K_{uf}^{H})$
(117)

where in the computation of the hybrid kernels we utilise $q(F_{l-1})$ for
computing the Wasserstein-2 distance, respectively samples from $q(F_{l-1})$
in the Euclidean part. It is worth highlighting the fact that for $l\geq 2$,
$Z_{l}$ are MVN. Further derivations will lead us to the DDGP ELBO:

$\mathcal{L}_{DDGP}=\mathbb{E}_{\\{q(F_{l},U_{l}\\}_{l=1}^{L}}p(Y|F_{L})-\sum_{l=2}^{L}\mathbb{E}_{Z_{l}}\left[KL\left[q(U_{l})\|p(U_{l})\right]\right]-KL\left[q(U_{1})\|p(U_{1})\right]$
(118)

###  4.4 Propagation of variance in DDGP

In this subsection, we investigate what are the necessary requirements for
reliable propagation of non-parametric variance through the hierarchy for OOD
data points.

######  Proposition 4.1.

We consider the approximate posterior DDGP as a composition of functions
$q_{L}=q_{L}\circ...\circ q_{1}$ with $q_{l}$ being given by
$q(F_{l}|F_{l-1})=\int
p(F_{l}|U_{l},F_{l-1})q(U_{l})=\mathcal{N}(F_{l}\mid\tilde{U_{l}}(F_{l-1}),\tilde{\Sigma_{l}}(F_{l-1}))$.
We consider an OOD data point $x_{out}$ in input space such that
$\tilde{\Sigma_{1}}(x_{out})=\sigma^{2}_{1}$ and $\tilde{U_{1}}(x_{out})=0$,
respectively $x_{in}$ to be in-distribution, with
$\tilde{\Sigma}_{1}(x_{in-d})=V_{in-d}\leq\sigma^{2}_{1}$ and
$\tilde{U}_{1}(x_{in-d})=M_{in-d}$. We assume that $M_{Z_{1}}$ are
equidistantly placed between $[-3\sigma_{1},3\sigma_{1}]$ and
$\Sigma_{Z_{1}}\leq\sigma^{2}_{1}$. The approximate variance in the second
layer of $x_{out}$ will be higher than $x_{in-d}$ if the following hold:

$\sigma^{2}_{1}>>\Sigma_{Z_{2}}~{}\textit{and}~{}V_{in-d}\approx\Sigma_{Z_{2}}$

###### Proof.

From the proof of Proposition 3.1 we remind the desired behaviour for the
second layer:

$\tilde{\Sigma}_{2}(M_{in-d})+V_{in-d}*\left(\frac{\partial\tilde{U}_{2}(F_{1})}{\partial
F_{1}}\right)^{2}\Bigr{|}_{\begin{subarray}{c}F_{1}=M_{in-d}\end{subarray}}\leq\tilde{\Sigma}_{2}(0)+\sigma^{2}*\left(\frac{\partial\tilde{U}_{2}(F_{1})}{\partial
F_{1}}\right)^{2}\Bigr{|}_{\begin{subarray}{c}F_{1}=0\end{subarray}}$ (119)

In contrast to zero mean DGP, where under current assumptions we obtained
$\tilde{\Sigma}_{2}(M_{in-d})\approx\tilde{\Sigma}_{2}(0)$, this is not the
case for DDGP. Expanding on the variance terms $\tilde{\Sigma}_{2}(M_{in-d})$
and $\tilde{\Sigma}_{2}(0)$ looking at their non-parametric variance:

$\sigma^{2}_{2}-K_{M_{in-d}u}^{H}K_{uu}^{H}\raisebox{4.95134pt}{$\scriptscriptstyle-\\!1$}K_{uM_{in-d}}^{H}\leq\sigma^{2}_{2}-K_{0u}^{H}K_{uu}^{H}\raisebox{4.95134pt}{$\scriptscriptstyle-\\!1$}K_{u0}^{H}$
(120)

To simplify the notation we can consider just one inducing point, lying
halfway between 0 and $M_{in-d}$ at $Z_{m}$ and with
$Z_{var}<<\sigma^{2}_{1}$. Using the hybrid kernel formulation and re-
arranging the terms in the above equation considering the matrices are scalars
we arrive at:

$\displaystyle\sigma^{2}_{2}$
$\displaystyle\exp{\frac{-W_{2}^{2}(\mathcal{N}(0,\sigma^{2}_{1}),\mathcal{N}(Z_{m},Z_{var}))}{l^{2}_{2}}}\exp{\frac{-(0-Z_{m})^{2}}{l^{2}_{2}}}\leq$
(121)
$\displaystyle\sigma^{2}_{2}\exp\frac{-W_{2}^{2}(\mathcal{N}(M_{in-d},V_{in-d}),\mathcal{N}(Z_{m},Z_{var}))}{l^{2}_{2}}\exp\frac{-(M_{in-d}-Z_{m})^{2}}{l^{2}_{2}}$

which after eliminating redundant terms and taking the logarithm we obtain:

$W_{2}^{2}(\mathcal{N}(0,\sigma^{2}_{1}),\mathcal{N}(Z_{m},Z_{var}))+(0-Z_{m})^{2}\geq
W_{2}^{2}(\mathcal{N}(M_{in-d},V_{in-d}),\mathcal{N}(Z_{m},Z_{var}))+(M_{in-d}-Z_{m})^{2}$
(122)

The above inequality holds if and only if $\sigma^{2}_{1}>>Z_{var}$ and
$V_{in-d}\approx Z_{var}$, thus proving our proposition.

∎

######  Remark 4.2.

In opposition to DGP which require a specific pattern of first order
derivatives surrounding 0, DDGP rely less on this constraint as the above
inequality is already providing a significant difference in non-parametric
variance between outlier and inlier points, thus satisfying Desiderata I. For
any layer $l$, $\Sigma_{Z_{l}}$ will converge towards values of the parametric
posterior variance of in-distribution points, since the inducing points are
encouraged to gravitate towards the in-distribution manifold at layer $l$. In
well optimized scenarios, the parametric posterior variance will be reduced
compared to the prior variance ($\sigma_{l-1}^{2}$). Considering this
observations, the conditions laid out in Proposition 4.1 will be satisfied in
practice.

###  4.5 _Feature-collapse_ in (D)DGP

Dunlop et al. (2018) provide an in-depth analysis of function space properties
DGP have, which DDGPs are equivalent to in prior space (see subsection 4.2 ).
We are interested to find what are the necessary requirements, given a few
constraints, of the function space values to collapse in the posterior for in-
distribution data points.

We now introduce some notation conventions. For $D\in\mathcal{R}^{d}$ and
$\\{D_{l}\in\mathcal{R}^{d_{l}}\\}_{l=1}^{L}$, where $d_{l}$ is the number of
dimensions in the l-th layer of the hierarchy. We consider the functions
$f_{l}:D_{l-1}\to\mathcal{R}^{d_{l}}$, with a DGP being defined as composition
of functions: $f_{n}\left(\cdots
f_{2}\left(f_{1}\left(x\right)\right)\right)$.

For a DGP at layer $n$ we assume the following posterior conditional
formulation for two data points $x,~{}x^{*}\in D_{in-d}$, respectively:

$p\begin{pmatrix}f_{n}(x)\\\
f_{n}(x^{*})\end{pmatrix}=\mathcal{N}\left[\begin{pmatrix}K_{xu}K_{uu}^{-1}U_{n}\\\
K_{x^{*}u}K_{uu}^{-1}U_{n}\end{pmatrix},\begin{pmatrix}K_{xx}-Q_{xx}&K_{xx^{*}}-Q_{xx^{*}}\\\
K_{x^{*}x}-Q_{x^{*}x}&K_{x^{*}x^{*}}-Q_{x^{*}x^{*}}\end{pmatrix}\right]$ (123)

For subsequent derivations, we denote
$\mu(x)=\mathcal{N}\left(K_{xu}K_{uu}^{-1}U_{n},K_{xx}-Q_{xx}\right)$.

######  Proposition 4.2.

We assume $f_{0}$ to be bounded on bounded sets almost-surely. For a DGP, if
at each layer $l$ we have satisfied the following inequality
$\frac{\sigma^{4}_{l}}{4l_{l}^{4}}\|K_{uu}^{-1}U_{l}\|_{2}^{2}\left(f_{l}(x)-f_{l}(x^{*})\right)^{2}<1$,
alongside $Q_{ff}\approx K_{ff}$, given two in-distributions points $x$ and
$x^{*}$ we have the following result:

$P\left(\|f_{n}(x)-f_{n}(x^{*})\|_{2}\to 0\mid\\{U_{l}\\}_{l=1}^{n}\right)=1$

###### Proof.

In subsequent derivations, we assume that the set of inducing points
$\\{Z_{l}\\}_{l=1}^{L}$ are sufficiently large and well-optimized so that for
points $x\in D$ we have $K_{ff}\approx K_{fu}K_{uu}^{-1}K_{uf}$.

We commence the derivation:

$\displaystyle\mathbb{E}\left[\|f_{n}(x)-f_{n}(x^{*})\|_{2}^{2}\mid
f_{n-1},U_{n}\right]$ $\displaystyle=\sum_{d=1}^{D_{l}}\mathbb{E}\left[\mid
f_{n}^{d}(x)-f_{n}^{d}(x^{*})\mid^{2}\mid f_{n-1},U_{n}\right]$ (124)
$\displaystyle=\sum_{d=1}^{D_{l}}\mid\left(K_{fu}-K_{f^{*}u}\right)K_{uu}^{-1}U_{n}^{d}\mid^{2}$
(125)

We can apply Cauchy-Schwarz inside the brackets to obtain the following
inequality:

$\displaystyle\leq\sum_{d=1}^{D_{l}}\|K_{fu}-K_{f^{*}u}\|_{2}^{2}\|K_{uu}^{-1}U_{n}^{d}\|_{2}^{2}$
(126)
$\displaystyle\leq\sum_{d=1}^{D_{l}}\|\sigma^{2}_{n}\exp{-\frac{\left(x-Z_{n}\right)^{2}}{2l_{n}^{2}}}-\sigma^{2}_{n}\exp{-\frac{\left(x^{*}-Z_{n}\right)^{2}}{2l_{n}^{2}}}\|_{2}^{2}\|K_{uu}^{-1}U_{n}^{d}\|_{2}^{2}$
(127)

We assume that $x=x^{*}+h$ and $x^{*}$ is a high-density data point with
regards to $Z_{n}$, respectively $\mid\mid x-Z_{n}\mid\mid_{2}^{2}\approx 0$.
Now we can re-arrange the last equation as:

$\displaystyle\leq\sum_{d=1}^{D_{l}}\|\sigma^{2}_{n}\exp{-\frac{\left(h+x^{*}-Z_{n}\right)^{2}}{2l_{n}^{2}}}-\sigma^{2}_{n}\|_{2}^{2}\|K_{uu}^{-1}U_{n}^{d}\|_{2}^{2}$
(128)
$\displaystyle\leq\sum_{d=1}^{D_{l}}\sigma^{4}_{n}\mid\exp{-\frac{\left(h\right)^{2}}{2l_{n}^{2}}}-1\mid^{2}\|K_{uu}^{-1}U_{n}^{d}\|_{2}^{2}$
(129)

We have the following general result $1-\exp{-x}<x,~{}\text{for}~{}x>-1$ which
we apply in our case since $\frac{\left(h\right)^{2}}{2l^{2}}>-1$, we obtain:

$\mathbb{E}\left[\|f_{n}(x)-f_{n}(x^{*})\|_{2}^{2}\mid
f_{n-1},U_{n}\right]\leq\sum_{d=1}^{D_{l}}\frac{\sigma^{4}_{n}}{4l_{n}^{4}}\|K_{uu}^{-1}U_{n}^{d}\|_{2}^{2}\left(h\right)^{2}\left(h\right)^{2}$
(130)

which we can re-express with regards to previous layer features to make
subsequent derivations easier to follow:

$\mathbb{E}\left[\|f_{n}(x)-f_{n}(x^{*})\|_{2}^{2}\mid
f_{n-1},U_{n}\right]\leq\sum_{d=1}^{D_{l}}\frac{\sigma_{n}^{4}}{4l_{n}^{4}}\|K_{uu}^{-1}U_{n}^{d}\|_{2}^{2}\left(h\right)^{2}\left(f_{n-1}(x)-f_{n-1}(x^{*})\right)^{2}$
(131)

We can now apply induction and the tower property of conditional expectation
to arrive at:

$\mathbb{E}\left[\|f_{n}(x)-f_{n}(x^{*})\|_{2}^{2}\mid\\{U_{l}\\}_{l=1}^{n}\right]\leq\prod_{l=1}^{n}\frac{\sigma_{l}^{4}}{4l_{l}^{4}}\|K_{uu}^{-1}U_{l}\|_{2}^{2}\left(h\right)^{2}\left(f_{0}(x)-f_{0}(x^{*})\right)^{2}$
(132)

By Markov’s inequality, for any $\epsilon>0$ we have that:

$P\left(\|f_{l+1}(x)-f_{l+1}(x^{*})\|_{2}\geq\epsilon\mid\\{U_{l}\\}_{l=1}^{n}\right)\leq\frac{1}{\epsilon^{2}}\prod_{l=1}^{n}\frac{\sigma_{l}^{4}}{4l_{l}^{4}}\|K_{uu}^{-1}U_{l}\|_{2}^{2}\left(h\right)^{2}\left(f_{0}(x)-f_{0}(x^{*})\right)^{2}$
(133)

We can apply the first Borel-Cantelli lemma to obtain:

$P\left(\cap_{l=1}^{\infty}\cup_{m=l}^{\infty}\|f_{m}(x)-f_{m}(x^{*})\|_{2}\geq\epsilon\mid\\{U_{l}\\}_{l=1}^{n}\right)=0$
(134)

Lastly, we can express the following:

$\displaystyle P\left(\|f_{n}(x)-f_{n}(x^{*})\|_{2}\to
0\mid\\{U_{l}\\}_{l=1}^{n}\right)$
$\displaystyle=P(\cap_{k=1}^{\infty}\cup_{l=1}^{\infty}\cap_{m=l}^{\infty}\|f_{m}(x)-f_{m}(x^{*})\|_{2}\leq\frac{1}{k}\mid\\{U_{l}\\}_{l=1}^{n})$
(135)
$\displaystyle=1-P(\cup_{k=1}^{\infty}\cap_{l=1}^{\infty}\cup_{m=l}^{\infty}\|f_{m}(x)-f_{m}(x^{*})\|_{2}\geq\frac{1}{k}\mid\\{U_{l}\\}_{l=1}^{n})$
(136) $\displaystyle~{}\geq
1-\sum\limits_{k=1}^{\infty}P(\cap_{l=1}^{\infty}\cup_{m=l}^{\infty}\|f_{m}(x)-f_{m}(x^{*})\|_{2}\geq\frac{1}{k}\mid\\{U_{l}\\}_{l=1}^{n})$
(137) $\displaystyle=1$ (138)

From which we obtain the proof of our proposition, respectively
$P\left(\|f_{n}(x)-f_{n}(x^{*})\|_{2}\to 0\mid\\{U_{l}\\}_{l=1}^{n}\right)=1$
∎

######  Remark 4.3.

From the constraint
$\frac{\sigma_{l}^{4}}{4l_{l}^{4}}\|K_{uu}^{-1}U_{l}\|_{2}^{2}\left(f_{l}(x)-f_{l}(x^{*})\right)^{2}<1$
to ensure _feature-collapse_ , we can see that in the converse case, of
avoiding _feature-collapse_ , this can be assured by either increasing the
kernel variance $\sigma^{2}_{l}$ or $\|K_{uu}^{-1}U_{l}\|_{2}^{2}$, which are
both responsible for the amplitude of sampled functions. Moreover, another
pathway to avoid _feature-collapse_ is to increase the lengthscale
$l^{2}_{l}$. Intuitively, this will ensure that unseen test points close to an
inducing point will reliably have their predicted mean almost equal to that
respective inducing point’s value. Conversely, for lower lengthscale values,
the predicted mean for an unseen points will be a weighted version of inducing
points’ output values across a larger distance, hence favouring similar values
for unseen points, which will eventually lead to _feature-collapse_.

Now we are interested in deriving similar constraints for DDGPs to suffer from
_feature-collapse_ and possibly seeing if there are any differences.

######  Proposition 4.3.

We assume $f_{0}$ to be bounded on bounded sets almost-surely. For a DDGP, if
at each layer $l$ except the first one we have satisfied the following
inequality
$\frac{\sigma^{4}_{l}}{4l^{4}_{l}}\|K_{uu}^{-1}U_{l}\|_{2}^{2}[f_{l}(x)-f_{l}(x*)+\frac{\|m(f_{l}(x))-m(f_{l}(x^{*}))\|_{2}^{2}+\|v(f_{l}(x))-v(f_{l}(x^{*}))\|_{2}^{2}}{f_{l}(x)-f_{l}(x*)}]^{2}\leq
1$, given two in-distributions points $x$ and $x^{*}$ s.t. $Q_{ff}\approx
K_{ff}$ we have the following result:

$P\left(\|f_{n}(x)-f_{n}(x^{*})\|_{2}\to 0\mid\\{U_{l}\\}_{l=1}^{n}\right)=1$

###### Proof.

We adapt our previous proof for Proposition 4.2 to suit the hybrid kernel
specific to DDGPs. We have the following upper bound:

$\displaystyle\mathbb{E}\left[\|f_{n}(x)-f_{n}(x^{*})\|_{2}^{2}\mid
f_{n-1},U_{n}\right]\leq$
$\displaystyle\sum_{d=1}^{D_{n}}\|\sigma^{2}_{n}\exp{-\frac{\left(x-Z_{n}\right)^{2}+W_{2}^{2}(\mu(x),\mu(Z_{n}))}{2l_{n}^{2}}}$
(139)
$\displaystyle-\sigma^{2}_{n}\exp{-\frac{\left(x^{*}-Z_{n}\right)^{2}+W_{2}^{2}(\mu(x^{*}),\mu(Z_{n}))}{2l_{n}^{2}}}\|_{2}^{2}\|K_{uu}^{-1}U_{n}^{d}\|_{2}^{2}$

We assume that $x=x^{*}+h$ and $x^{*}$ is a high-density data point with
regards to $Z_{n}$, respectively $\mid\mid x-Z_{n}\mid\mid_{2}^{2}\approx 0$.
For their first two moments we have equivalent assumptions
$m(x)=m(x^{*})+m(h)$ and $v(x)=v(x^{*})+v(h)$, with
$\mu(x^{*})\stackrel{{\scriptstyle d}}{{\approx}}\mu(Z_{n})$ We can re-arrange
the last equation as:

$\leq\sum_{d=1}^{D_{n}}\|\sigma^{2}_{n}\exp{-\frac{\left(h+x^{*}-Z_{n}\right)^{2}+\|m(h)\|_{2}^{2}+\|v(h)\|_{2}^{2}}{2l_{n}^{2}}}-\sigma^{2}\|_{2}^{2}\|K_{uu}^{-1}U_{n}^{d}\|_{2}^{2}$
(140)

We have the following general result $1-\exp{-x}<x,~{}\text{for}~{}x>-1$ which
we apply in our case since $\frac{\left(h\right)^{2}+\mid m(h)\mid^{2}+\mid
v(h)\mid^{2}}{2l^{2}}>-1$:

$\mathbb{E}\left[\|f_{n}(x)-f_{n}(x^{*})\|_{2}^{2}\mid
f_{n-1},U_{n}\right]\leq\sum_{d=1}^{D_{n}}\frac{\sigma_{n}^{4}}{4l_{n}^{4}}\|K_{uu}^{-1}U_{n}^{d}\|_{2}^{2}\left(h^{2}+\|m(h)\|_{2}^{2}+\|v(h)\|_{2}^{2}\right)^{2}$
(141)

We can see that we can re-express the final term of above equation as follows:

$\displaystyle\left(h^{2}+\underbrace{\|m(h)\|_{2}^{2}+\|v(h)\|_{2}^{2}}_{c}\right)^{2}$
$\displaystyle=\left(h^{2}+c\right)^{2}$ (142)
$\displaystyle=h^{4}+2h^{2}c+c^{2}$ (143)
$\displaystyle=h^{2}\left[h^{2}+2c+\frac{c^{2}}{h^{2}}\right]$ (144)
$\displaystyle=h^{2}\left[h+\frac{c}{h}\right]^{2}$ (145)

We can apply this side derivation to equation (141):

$\mathbb{E}\left[\|f_{n}(x)-f_{n}(x^{*})\|_{2}^{2}\mid
f_{n-1},U_{n}\right]\leq\sum_{d=1}^{D_{n}}\frac{\sigma_{n}^{4}}{4l_{n}^{4}}\|K_{uu}^{-1}U_{n}^{d}\|_{2}^{2}h^{2}\left[h+\frac{\|m(h)\|_{2}^{2}+\|v(h)\|_{2}^{2}}{h}\right]^{2}h^{2}$
(146)

We can now apply induction and the tower property of conditional expectation
to arrive at:

$\mathbb{E}\left[\|f_{n}(x)-f_{n}(x^{*})\|_{2}^{2}\mid\\{U_{l}\\}_{l=2}^{n}\right]\leq\prod_{l=2}^{L}\frac{\sigma_{l}^{4}}{4l_{l}^{4}}\|K_{uu}^{-1}U_{l}^{d}\|_{2}^{2}\left[h+\frac{\|m(h)\|_{2}^{2}+\|v(h)\|_{2}^{2}}{h}\right]^{2}\left(f_{1}(x)-f_{1}(x^{*})\right)^{2}$
(147)

Using the upper bound on $\left(f_{1}(x)-f_{1}(x^{*})\right)^{2}$ from the
proof of Proposition 4.2, we obtain:

$\displaystyle\mathbb{E}\left[\|f_{n}(x)-f_{n}(x^{*})\|_{2}^{2}\mid\\{U_{l}\\}_{l=1}^{n}\right]$
$\displaystyle\leq\frac{\sigma^{4}_{1}}{4l_{1}^{4}}\|K_{uu}^{-1}U_{1}^{d}\|_{2}^{2}\prod_{l=2}^{L}\frac{\sigma_{l}^{4}}{4l_{l}^{4}}\|K_{uu}^{-1}U_{l}^{d}\|_{2}^{2}\left[h+\frac{\|m(h)\|_{2}^{2}+\|v(h)\|_{2}^{2}}{h}\right]^{2}\left(f_{0}(x)-f_{0}(x^{*})\right)^{2}$
(148)

Similar derivations to the proof of Proposition 4.2 can be applied to arrive
at the final result. ∎

######  Remark 4.4.

The constraint for DDGPs to suffer from _feature-collapse_
$\frac{\sigma^{4}_{l}}{4l_{l}^{4}}\|K_{uu}^{-1}U_{l}\|_{2}^{2}[f_{l}(x)-f_{l}(x*)+\frac{\|m(f_{l}(x))-m(f_{l}(x^{*}))\|_{2}^{2}+\|v(f_{l}(x))-v(f_{l}(x^{*}))\|_{2}^{2}}{f_{l}(x)-f_{l}(x*)}]^{2}\leq
1$, given two in-distributions points $x$ and $x^{*}$ s.t. $Q_{ff}\approx
K_{ff}$ is similar to the one for DGPs expect the final term in the brackets.
We can see that for any given value of $f_{l}(x)-f_{l}(x^{*})$, since
$\|m(f_{l}(x))-m(f_{l}(x^{*}))\|_{2}^{2}+\|v(q(f_{l}(x))-v(f_{l}(x^{*}))\|_{2}^{2}$
is always positive, it will result in
$[f_{l}(x)-f_{l}(x*)+\frac{\|m(f_{l}(x))-m(f_{l}(x^{*}))\|_{2}^{2}+\|v(f_{l}(x))-v(f_{l}(x^{*}))\|_{2}^{2}}{f_{l}(x)-f_{l}(x*)}]^{2}\geq\left[f_{l}(x)-f_{l}(x*)\right]^{2}$.
This means that DDGPs are more suitable to guard against _feature-collapse_
than DGPs ceteris paribus. Moreover, we posit that the dependence on moment
differences between two data points will embed DDGPs with _smoother_
properties in comparison to DGPs, since the moments are computed using squared
exponential kernels, which will produce smoother functions. With these two
properties in mind, we can affirm that DDGPs satisfy Desiderata II.

Figure 5: Progressive collapse of distributional uncertainty in DGPs.
Distributional uncertainty of DGPs trained on a toy dataset with data between
-5.0 and 5.0 with a varying number of inducing points. X axis denotes input
space.

##  5 Experiments

##### Collapse of non-parametric variance in DGP.

Duvenaud et al. (2014) highlighted a certain pathology when stacking GP,
respectively the increasingly more non-injective mappings as the network
increases in number of layers. In this section we highlight another pathology
of DGP with zero-mean functions, namely the progressive collapse of non-
parametric uncertainty for larger number of inducing points. This erroneous
behaviour is also present in relatively small networks with just two hidden
layers (Figure 5).

To understand what is causing this pathology, we take a simple case study of a
DGP with two hidden layers trained on a toy regression dataset (Figure 6).
Taking a clear outlier in input space, say the data point situated at -7.5, it
gets correctly identified as an outlier in the mapping from input space to
hidden layer space as given by its distributional variance. However, its
outlier property dissipates in the next layer after sampling, as it gets
mapped to regions where the next GP assigns inducing point locations. This is
due to points inside the data manifold getting confidently mapped between -2.5
and 2.5 in hidden layer space. Consequently, what was initially correctly
identified as an outlier will now have its final distributional uncertainty
close to zero. Adding further layers, will only compound this pathology.

Figure 6: Rationale behind collapse of non-parametric uncertainty. Layer-wise
decomposition of uncertainty into parametric and non-parametric for a zero-
mean DGP. OOD points in input space (red dots) get sampled in hidden layers
close to in-distribution points (green dots), hence close to inducing point
locations.

##### Interplay between PCA mean functions and non-parametric uncertainty.

In this subsection we investigate the effect that PCA mean function has on
distributional uncertainty. In the case that the hidden layers’ size is equal
or larger than input size, then DGPs will retain their OOD capabilities.
Conversely, if the hidden layers’ size is considerably lower than input size
(such as image classification) we posit that it will negatively affect OOD
capabilities (Figure 5). We use the “banana” dataset and train our proposed
models with 2 hidden layers and 10 inducing points. The DGP distributional
uncertainty is guided by the direction of the PCA projection of the data,
thereby widely extrapolating low uncertainty regions to spatial locations
where there is no data. For DDGP the PCA projection only slightly influences
distributional variance in the last layer (Figure 7).

(a) DGP (PCA mean function)

(b) DDGP (PCA mean function)

Figure 7: PCA mean functions influence distributional uncertainty in DGPs.
Layer-wise moments of specified models trained on the “banana” dataset.
Distributional uncertainty in hidden layers is lower on the data manifold
compared to outside the manifold for DDGPs. Conversely, DGPs widely
extrapolate low distributional uncertainty outside the data manifold.

##### Smoothness properties.

From Proposition 4.2 and 4.3 we would expect that DDGPs exhibit smoother
hidden layers features. Using our previous experiment on the “banana” dataset,
we keep track of correlations between two focus points with respect to full
input space (Figure 8). The increased “locality” of features in layer 3 is
striking in the case of DDGPs. Moreover, we notice even in layer 2 a slight
difference in smoothness between the two models.

Figure 8: Smoothness comparison between DGPs and DDGPs. Correlation of two
focus data points (purple dots) with respect to full input space spectrum at
hidden layers for DGPs (PCA mean function) and DDGPs (PCA mean function)
trained on “banana” dataset. DDGPs (right panel) exhibit smoother and more
“local” correlations.

##### Performance on real-life data.

Besides testing whether DDGPs surpass DGPs in terms of OOD detection, it is
also of vital importance to see whether DDGPs manage to obtain calibrated
predictive distributions on medium and large scale datasets. For this we
evaluate our models on datasets from the UCI machine learning dataset
repository. All of the experiments were run with the same initializations of
variational parameters and kernel hyperparameters, with the goal of comparing
DDGPs with DGPs for 2,3, respectively 4 layers. We used 100 inducing points at
each layer, with 5 dimensions per hidden layer and PCA mean functions. Further
details regarding training and initialization are given in Appendix
LABEL:apd:uci_details.

(a) Test log-likelihood regression datasets

(b) Test log-likelihood classification datasets

Figure 9: Results on UCI datasets. Top: Test log-likelihood values on
regression datasets. Bottom: Test log-likelihood values on classification
datasets. All subplots: 10 different runs of each model with different
initialization seeds are taken into the composition of each boxplot. Higher
values (to the right) indicate better model fit.

In terms of predictive log-likelihood values on test data, DDGPs generally
surpass their Euclidean only counterpart on regression datasets, with
diminished performance seen on “WineRed” and “WineWhite”. In terms of
predictive testing set results on classification datasets, the behaviour is
roughly similar, with noticeable improved performance of DDGPs on “Avila”,
“Letters” and “MNIST”. Conversely, the same can be said for DGPs on
“Wisconsin”. On remaining classification datasets the behaviour is similar
(Figure 9).

##### Out-of-distribution detection.

In this subsection we adapt DDGP to the architecture introduced in Deep
Convolutional GP (DCGP) (Blomqvist et al., 2018), obtaining the Distributional
Deep Convolutional GP (DistDCGP). The goal is to investigate their OOD
detection capabilities on real data. To quantify the departure from normative
data, we need a scalar measure encompassing the uncertainty present across all
classes. For this we use the distributional differential entropy introduced in
equation (104). This has the effect to measure the overall distributional
uncertainty in logit space, with higher values representing more uncertainty.

In these experiments we assess the capacity of our model to detect dataset
shift by training it on MNIST and looking at the uncertainty measures computed
on the testing set of MNIST and NotMNIST and SVHN. The hypothesis is that we
ought to see higher differential entropy for distributional uncertainty for
the digits stemming from a wide array of fonts present in NotMNIST and SVHN as
none of the fonts are handwritten. We additionally compare with FashionMNIST,
to discern if models are capable of flagging anomalous inputs when deployed in
a user-facing system.

Model |  Hidden Layers |  MNIST vs. NotMNIST AUC |  MNIST vs. SVHN AUC |  MNIST vs. FashionMNIST AUC
---|---|---|---|---
Inducing Points | 50 | 100 | 50 | 100 | 50 | 100
DCGP | 1 | $0.71$ | $0.69$ | 1.0 | 1.0 | 0.98 | $0.98$
DistDCGP | 1 | 0.87 | 0.95 | 1.0 | 1.0 | 0.98 | 1.0
DCGP | 2 | $0.64$ | $0.59$ | $0.97$ | $0.99$ | $0.94$ | $0.92$
DistDCGP | 2 | 0.97 | 0.84 | 1.0 | 1.0 | 1.0 | 0.98
DCGP | 3 | $0.64$ | $0.70$ | $0.97$ | $0.89$ | $0.94$ | $0.84$
DistDCGP | 3 | 0.95 | 0.86 | 1.0 | 1.0 | 0.96 | 0.92
Table 1: OOD detection results. Performance of OOD detection based on
distributional uncertainty differential entropy. Models with different
inducing point numbers in hidden layers are trained on MNIST (normative data).

We notice an abrupt loss in AUC scores on SVHN and FashionMNIST for 3 layers
DCGP when transitioning from 50 to 100 inducing points. This reinforces the
observed pathology that increasing inducing points results in loss of OOD
capabilities. On NotMNIST we notice a substantial discrepancy between the two
in favor of DistDCGP (Table 1). Overall, DistDCGP manages to surpass DCGP in
any given comparison regardless of network architecture. Further figures
showcasing differences in their distributional differential entropy histograms
are provided in Appendix LABEL:apd:additional_ood_detection.

##  6 Discussion

We have shown under a set of assumptions that total uncertainty can be
correctly propagated under the DSVI framework introduced in Salimbeni and
Deisenroth (2017) only in the case that the first order derivative of the
predictive mean function attains maximum values at the point where previously
identified OOD points get mapped to. Intuitively, total uncertainty for OOD
points will be higher compared to in-distribution points, if the parametric
component of the current layer’s SVGP provides a diverse set of function
samples when evaluated at points where the OOD point got mapped in the
previous layer.

The collapse of distributional uncertainty as we increase the number of
inducing points in the hidden layers is caused by the progressive mapping of
data points which in lower layers were flagged as OOD to regions close to
inducing points in the next layer in the case that in-distribution data points
are mapped in intervals centred around zero in hidden layer space. This leads
to the deactivation of $h(\cdot)$ (102) (non-parametric component in SVGP) in
all layers, except the first. Hence, in these situations DGPs rely just on the
variance stemming from $g(\cdot)$ (103) (parametric component of SVGP).

The PCA based projections stemming from the training data as introduced in
Salimbeni and Deisenroth (2017) help to solve these aforementioned issues.
However, these also come with a hidden cost, as the PCA projections will
wildly extrapolate the regions of low distributional uncertainty in directions
perpendicular to the eigenvectors in the case of high-dimensional datasets
that are projected down in hidden layer space (see experiments on “banana"
dataset). Our proposed framework solves this problem as the Wasserstein-2 part
of our hybrid kernel explicitly takes into account differences between the
variance of hidden layer SVGPs. Intuitively, even if the mean component of the
SVGP is dominated by the PCA projection, which could coincide with the
projection values of the training set, if the variance component is higher
then the Wasserstein-2 distance will highlight this domain shift with higher
distributional uncertainty. Besides this intuition, we have also provided
theoretical guarantees that total uncertainty is higher for OOD points
compared to in-distribution points for DDGPs in the case that the inducing
locations’ variance is close in distribution to total uncertainty of in-
distribution points at that hidden layer. This is a scenario which happens in
practice as the variance of inducing locations will be optimized to match that
of the in-distribution points. Future work should explore methods to
explicitly enforce this match in distribution between inducing locations and
in-distribution predictive distributions. Recent work (Ustyuzhaninov et al.,
2019; Ober and Aitchison, 2021) has proposed to propagate the initial inducing
locations in input space through the hierarchy, thus obtaining predictive
distributions for these “global" inducing points at each layer. Under this
construction, the distribution match between inducing locations and in-
distribution points will be guaranteed at each layer in DDGPs.

Besides theoretical guarantees on correct propagation of distributional
uncertainty in DDGPs, we have also shown that DDGPs have theoretical
guarantees to better guard against _feature-collapse_. Moreover, the
Wasserstein-2 distance based on predictive moments will cause DDGPs to be
smoother compared to DGPs (see Figure 8). Due to all these attractive
properties of DDGPs, it is perhaps unsurprising that on OOD detection tasks
our proposed model has managed to surpass baseline models. Additionally, we
analyze their sensitivity to input perturbations (see Appendix
LABEL:apd:input_perturbation).

Moreover, we discover another pathology pertaining to DGPs, respectively that
when deploying very wide (i.e., high-dimensional hidden layers) DGPs, samples
can collapse to a thin hyper-sphere, thus requiring increasing inducing
locations to cover it. In Appendix LABEL:apd:fail_high_dimensions we offer
more details on this pathology, also showing that hybrid kernels can guard
against it.

In conclusion, we have shown that DDGPs are better suited for outlier
detection on both toy and real datasets, while also showing slight
improvements at testing time in medium and large scale datasets. Future work
should consider adapting DDGPs to safety-critical domains such as biomedical
image segmentation (Czolbe et al., 2021) or semantic segmentation (Franchi et
al., 2020).

## References

  * Timonen et al. (2019) Juho Timonen, Henrik Mannerström, Aki Vehtari, and Harri Lähdesmäki.  lgpr: An interpretable nonparametric method for inferring covariate effects from longitudinal data.  _arXiv preprint arXiv:1912.03549_ , 2019.
  * Wang et al. (2020) Dennis Wang, James Hensman, Ginte Kutkaite, Tzen S Toh, Ana Galhoz, Jonathan R Dry, Julio Saez-Rodriguez, Mathew J Garnett, Michael P Menden, Frank Dondelinger, et al.  A statistical framework for assessing pharmacological responses and biomarkers using uncertainty estimates.  _Elife_ , 9:e60352, 2020.
  * Chen et al. (2014) Tongtong Chen, Bin Dai, Ruili Wang, and Daxue Liu.  Gaussian-process-based real-time ground segmentation for autonomous land vehicles.  _Journal of Intelligent & Robotic Systems_, 76(3):563–582, 2014.
  * Damianou and Lawrence (2013) Andreas Damianou and Neil Lawrence.  Deep gaussian processes.  In _Artificial Intelligence and Statistics_ , pages 207–215, 2013\.
  * Ustyuzhaninov et al. (2019) Ivan Ustyuzhaninov, Ieva Kazlauskaite, Markus Kaiser, Erik Bodin, Neill DF Campbell, and Carl Henrik Ek.  Compositional uncertainty in deep gaussian processes.  _arXiv preprint arXiv:1909.07698_ , 2019.
  * Bachoc et al. (2017) François Bachoc, Fabrice Gamboa, Jean-Michel Loubes, and Nil Venet.  A gaussian process regression model for distribution inputs.  _IEEE Transactions on Information Theory_ , 64(10):6620–6637, 2017.
  * Bachoc et al. (2018) Francois Bachoc, Alexandra Suvorikova, David Ginsbourger, Jean-Michel Loubes, and Vladimir Spokoiny.  Gaussian processes with multidimensional distribution inputs via optimal transport and hilbertian embedding.  _arXiv preprint arXiv:1805.00753_ , 2018.
  * Thi Thien Trang et al. (2019) Bui Thi Thien Trang, Jean-Michel Loubes, Laurent Risser, and Patricia Balaresque.  Distribution regression model with a reproducing kernel hilbert space approach.  _Communications in Statistics-Theory and Methods_ , pages 1–23, 2019\.
  * Kendall and Gal (2017) Alex Kendall and Yarin Gal.  What uncertainties do we need in bayesian deep learning for computer vision?  _arXiv preprint arXiv:1703.04977_ , 2017.
  * Lakshminarayanan et al. (2017) Balaji Lakshminarayanan, Alexander Pritzel, and Charles Blundell.  Simple and scalable predictive uncertainty estimation using deep ensembles.  _Advances in neural information processing systems_ , 30, 2017.
  * Hendrycks and Gimpel (2016) Dan Hendrycks and Kevin Gimpel.  A baseline for detecting misclassified and out-of-distribution examples in neural networks.  _arXiv preprint arXiv:1610.02136_ , 2016.
  * Liang et al. (2017) Shiyu Liang, Yixuan Li, and Rayadurgam Srikant.  Enhancing the reliability of out-of-distribution image detection in neural networks.  _arXiv preprint arXiv:1706.02690_ , 2017.
  * Hafner et al. (2020) Danijar Hafner, Dustin Tran, Timothy Lillicrap, Alex Irpan, and James Davidson.  Noise contrastive priors for functional uncertainty.  In _Uncertainty in Artificial Intelligence_ , pages 905–914. PMLR, 2020.
  * Malinin and Gales (2018) Andrey Malinin and Mark Gales.  Predictive uncertainty estimation via prior networks.  In _Advances in Neural Information Processing Systems_ , pages 7047–7058, 2018.
  * Domingues et al. (2018) Rémi Domingues, Pietro Michiardi, Jihane Zouaoui, and Maurizio Filippone.  Deep gaussian process autoencoders for novelty detection.  _Machine Learning_ , 107(8):1363–1383, 2018.
  * Baur et al. (2021) Christoph Baur, Stefan Denner, Benedikt Wiestler, Nassir Navab, and Shadi Albarqouni.  Autoencoders for unsupervised anomaly segmentation in brain mr images: a comparative study.  _Medical Image Analysis_ , page 101952, 2021.
  * Pugachev (2013) Vladimir Semenovich Pugachev.  _Theory of Random Functions: And Its Application to Control Problems_ , volume 5.  Elsevier, 2013.
  * Neal (2012) Radford M Neal.  _Bayesian learning for neural networks_ , volume 118.  Springer Science & Business Media, 2012.
  * Paananen et al. (2019) Topi Paananen, Juho Piironen, Michael Riis Andersen, and Aki Vehtari.  Variable selection for gaussian processes via sensitivity analysis of the posterior predictive distribution.  In _The 22nd International Conference on Artificial Intelligence and Statistics_ , pages 1743–1752. PMLR, 2019.
  * Lalchand and Rasmussen (2020) Vidhi Lalchand and Carl Edward Rasmussen.  Approximate inference for fully bayesian gaussian process regression.  In _Symposium on Advances in Approximate Bayesian Inference_ , pages 1–12. PMLR, 2020.
  * Titsias (2009) Michalis Titsias.  Variational learning of inducing variables in sparse gaussian processes.  In _Artificial Intelligence and Statistics_ , pages 567–574, 2009\.
  * Hoffman et al. (2013) Matthew D Hoffman, David M Blei, Chong Wang, and John Paisley.  Stochastic variational inference.  _The Journal of Machine Learning Research_ , 14(1):1303–1347, 2013.
  * Girard (2004) Agathe Girard.  _Approximate methods for propagation of uncertainty with Gaussian process models_.  University of Glasgow (United Kingdom), 2004.
  * Quinonero-Candela et al. (2003) Joaquin Quinonero-Candela, Agathe Girard, and Carl Edward Rasmussen.  Prediction at an uncertain input for gaussian processes and relevance vector machines-application to multiple-step ahead time-series forecasting.  2003\.
  * Salimbeni and Deisenroth (2017) Hugh Salimbeni and Marc Deisenroth.  Doubly stochastic variational inference for deep gaussian processes.  In _Advances in Neural Information Processing Systems_ , pages 4588–4599, 2017.
  * Villani (2008) Cédric Villani.  _Optimal transport: old and new_ , volume 338.  Springer Science & Business Media, 2008.
  * Dowson and Landau (1982) D. Dowson and B. Landau.  The fréchet distance between multivariate normal distributions.  _Journal of Multivariate Analysis_ , 12:450–455, 1982.
  * Almgren (1966) Frederick J Almgren.  _Plateau’s problem: an invitation to varifold geometry_ , volume 13.  American Mathematical Soc., 1966.
  * Charon and Trouvé (2013) Nicolas Charon and Alain Trouvé.  The varifold representation of nonoriented shapes for diffeomorphic registration.  _SIAM journal on Imaging Sciences_ , 6(4):2547–2580, 2013.
  * Antonsanti et al. (2021) Pierre-Louis Antonsanti, Joan Glaunès, Thomas Benseghir, Vincent Jugnon, and Irène Kaltenmark.  Partial matching in the space of varifolds.  In _International Conference on Information Processing in Medical Imaging_ , pages 123–135. Springer, 2021.
  * Girard et al. (2002) A Girard, CE Rasmussen, and R Murray-Smith.  Gaussian process priors with uncertain inputs: Multiple-step ahead prediction. univ. glasgow.  Technical report, Glasgow, Technical Report TR-2002-119, 2002.
  * Ruff et al. (2021) Lukas Ruff, Jacob R Kauffmann, Robert A Vandermeulen, Grégoire Montavon, Wojciech Samek, Marius Kloft, Thomas G Dietterich, and Klaus-Robert Müller.  A unifying review of deep and shallow anomaly detection.  _Proceedings of the IEEE_ , 2021.
  * Liu et al. (2020) Jeremiah Zhe Liu, Zi Lin, Shreyas Padhy, Dustin Tran, Tania Bedrax-Weiss, and Balaji Lakshminarayanan.  Simple and principled uncertainty estimation with deterministic deep learning via distance awareness.  _arXiv preprint arXiv:2006.10108_ , 2020.
  * Rosca et al. (2020) Mihaela Rosca, Theophane Weber, Arthur Gretton, and Shakir Mohamed.  A case for new neural network smoothness constraints.  2020\.
  * van Amersfoort et al. (2021) Joost van Amersfoort, Lewis Smith, Andrew Jesson, Oscar Key, and Yarin Gal.  On feature collapse and deep kernel learning for single forward pass uncertainty.  _arXiv preprint arXiv:2102.11409_ , 2021.
  * Dunlop et al. (2018) Matthew M Dunlop, Mark A Girolami, Andrew M Stuart, and Aretha L Teckentrup.  How deep are deep gaussian processes?  _Journal of Machine Learning Research_ , 19(54):1–46, 2018.
  * Duvenaud et al. (2014) David Duvenaud, Oren Rippel, Ryan Adams, and Zoubin Ghahramani.  Avoiding pathologies in very deep networks.  In _Artificial Intelligence and Statistics_ , pages 202–210, 2014\.
  * Blomqvist et al. (2018) Kenneth Blomqvist, Samuel Kaski, and Markus Heinonen.  Deep convolutional gaussian processes.  _arXiv preprint arXiv:1810.03052_ , 2018.
  * Ober and Aitchison (2021) Sebastian W Ober and Laurence Aitchison.  Global inducing point variational posteriors for bayesian neural networks and deep gaussian processes.  In _International Conference on Machine Learning_ , pages 8248–8259. PMLR, 2021.
  * Czolbe et al. (2021) Steffen Czolbe, Kasra Arnavaz, Oswin Krause, and Aasa Feragen.  Is segmentation uncertainty useful?  In _International Conference on Information Processing in Medical Imaging_ , pages 715–726. Springer, 2021.
  * Franchi et al. (2020) Gianni Franchi, Andrei Bursuc, Emanuel Aldea, Severine Dubuisson, and Isabelle Bloch.  One versus all for deep neural network incertitude (ovnni) quantification.  _arXiv preprint arXiv:2006.00954_ , 2020.