Residual Deep Gaussian Processes on Manifolds

## Abstract

Abstract We propose practical deep Gaussian process models on Riemannian manifolds, similar in spirit to residual neural networks.
With manifold-to-manifold hidden layers and an arbitrary last layer, they can model manifold- and scalar-valued functions, as well as vector fields.
We target data inherently supported on manifolds, which is too complex for shallow Gaussian processes thereon.
For example, while the latter perform well on high-altitude wind data, they struggle with the more intricate, nonstationary patterns at low altitudes.
Our models significantly improve performance in these settings, enhancing prediction quality and uncertainty calibration, and remain robust to overfitting, reverting to shallow models when additional complexity is unneeded.
We further showcase our models on Bayesian optimisation problems on manifolds, using stylised examples motivated by robotics, and obtain substantial improvements in later stages of the optimisation process.
Finally, we show our models to have potential for speeding up inference for non-manifold data, when, and if, it can be mapped to a proxy manifold well enough.

## 1 Introduction

Gaussian processes (GPs) are a widely adopted model class for learning functions within the Bayesian framework rasmussen2006 ().
They offer accurate uncertainty estimates and perform well even when data is scarce.
Consequently, GPs have found success in decision-making tasks, where well-calibrated uncertainty is key, including Bayesian optimisation snoek2012 (), active krause2008 () and reinforcement kamthe2018 () learning.

In recent years, substantial work went into developing the analogs of practical GP models on various non-Euclidean domains borovitskiy2020 (; ; ; ).
By virtue of being geometry-aware, these analogs have demonstrated improved performance in a variety of tasks on non-Euclidean spaces.
Their notable applications include Bayesian optimisation on manifolds for robotics jaquier2022 (), traffic flow interpolation on road networks borovitskiy2021 (), and wind velocity prediction on the globe hutchinson2021 (; ).
These models were also used to speed up inference for Euclidean GPs by transferring data to a hypersphere and leveraging the attractive structure GPs thereon possess dutordoir2020 ().

Despite their advantages, GPs can sometimes fall short in modelling complex, irregular functions.
To address this, deep Gaussian processes have been introduced as a sequential composition of GPs damianou2013 (), providing improved flexibility through their layered structure dai2016 (; ).
Many techniques developed for shallow GPs, such as variational inference salimbeni2017 () and efficient sampling techniques wilson2020 (), can be adapted for deep GPs, enabling them to be efficiently trained and deployed on large datasets.
This scalability is vital when dealing with data with complex, irregular patterns.

Deep Gaussian processes have demonstrated success in handling complex data of Euclidean nature, often competing with Bayesian neural networks and deep ensembles.
However, there has been limited work on expressive uncertainty-quantifying models on manifolds beyond shallow GPs.
This gap leads to the natural question: how can we construct deep Gaussian processes on manifolds?

By analogy with the Euclidean case, a deep GP on a manifold should be a composition of GP layers which take inputs and produce outputs on the manifold of interest.
While significant advancements have been made in handling manifold-input GPs, outputs on manifolds conflict with the fundamental concept of a GP, which dictates that outputs must be Gaussian and thus Euclidean.(^†^†Nevertheless, some practical heuristics for this exist in the literature mallasto2018 ().)
Designing GPs with inputs or outputs on a manifold is challenging; with both, it is even more so.

Our solution to this problem is in part inspired by residual neural networks, thus we term our models *residual deep Gaussian processes*.
Instead of constructing a manifold-to-manifold GP layer directly, we represent it as a *Gaussian vector field* (GVF) combined with an exponential map.
The former represents a displacement vector, a deviation from the identity map, or a *residual*, while the latter translates the input by the given displacement vector.
The mean of a layer is always the output of the previous layer.
We visualise this architecture in [Figure 1](#S1.F1), with sphere as the manifold of interest.
Notably, by changing the last layer only, one can get manifold-valued or vector-valued deep GP.
As [Section 3](#S3) will show, our residual deep Gaussian processes generalise the deep GP architecture of salimbeni2017 (), perhaps the most successful architecture in the Euclidean case.

Figure: Figure 1: Schematic illustration of a scalar-valued residual deep GP with $L$ hidden layers. The last layer is a scalar-valued GP on the manifold. If it is not present, the model is manifold-valued. If it is replaced with a Gaussian vector field (GVF), the model is a vector field on the manifold.
Refer to caption: /html/2411.00161/assets/figures/residual_deep_gp_schematic/a.png

We examine residual deep GPs through synthetic and real-world experiments, demonstrating our models’ superior performance over shallow geometry-aware GPs on tasks where complex data inherently lies on a manifold.
Additionally, we show that our models offer prospective avenues for accelerating inference for inherently Euclidean data in the context of deep GPs.
Our main focus is on hypersphere manifolds $\mathbb{S}_{d}$, due to their importance in key applications such as climate modelling and robotics, as well as their particularly simple structure that allows for more powerful and specialised GVFs robertnicoud2024 ().
However, the applicability of our model extends to all Riemannian manifolds, including ones represented by meshes, indicating an even broader potential.

## 2 Background

Mathematically, a Gaussian process (GP) is a real-valued random function $f$ whose marginals are jointly Gaussian.
The same term is also used for the respective distribution over functions.
For such $f$ there always exist a mean $\mu:X\rightarrow\mathbb{R}$ on the input domain $X$ of $f$ and a kernel $k:X\times X\rightarrow\mathbb{R}$ such that $f(\boldsymbol{x})\sim\mathcal{N}(\mu(\boldsymbol{x}),k(\boldsymbol{x},\boldsymbol{x}))$ for all $\boldsymbol{x}\subseteq X$.
In this case, we write $f\sim\mathcal{GP}(\mu,k)$.

GPs define useful priors for learning functions from noisy observations $\boldsymbol{y}\in\mathbb{R}^{n}$ at given input locations ${\boldsymbol{x}\subseteq X}$ within the Bayesian framework.
In fact, if the observation likelihood is assumed Gaussian $\boldsymbol{y}\mid f(\boldsymbol{x})\sim\mathcal{N}(f(\boldsymbol{x}),\mathbf{\Sigma})$, then the posterior is a GP rasmussen2006 () with

$$ $\displaystyle\mu_{f\mid\boldsymbol{y}}(\cdot)$ $\displaystyle=\mu(\cdot)+k(\cdot,\boldsymbol{x})(k(\boldsymbol{x},\boldsymbol{x})+\mathbf{\Sigma})^{-1}\left(\boldsymbol{y}-\mu(\boldsymbol{x})\right),$ (1) $\displaystyle k_{f\mid\boldsymbol{y}}(\cdot,\cdot^{\prime})$ $\displaystyle=k(\cdot,\cdot^{\prime})-k(\cdot,\boldsymbol{x})(k(\boldsymbol{x},\boldsymbol{x})+\mathbf{\Sigma})^{-1}k(\boldsymbol{x},\cdot^{\prime}).$ (2) $$

The input domain $X$ can, in principle, be any set.
However, a good kernel $k:X\times X\to\mathbb{R}$ is necessary to define practical GP models on $X$.
If $X=\mathbb{R}^{d}$, the most widely used kernels are from the Matérn family rasmussen2006 (), including, as a limiting case, the especially popular squared exponential kernel.(^†^†This kernel has many names; it is also known as the RBF, Gaussian, heat, or diffusion kernel.)
These kernels are attractive because they implement two natural inductive biases: (1) model behaviour should not change under any symmetries of the data—translations, in the case of $\mathbb{R}^{d}$—and (2) the unknown function possesses a certain degree of smoothness.
It turns out that the Matérn family can be generalised to various non-Euclidean domains $X$ in such a way that it still implements (1) and (2).
We now discuss such a generalisation to Riemannian manifolds.

### 2.1 Gaussian Processes on Riemannian Manifolds

A principled way of generalising the family of Matérn kernels to Riemannian manifolds was proposed by lindgren2011 () based on the ideas dating back to whittle1963 ().
borovitskiy2020 () showed that the resulting kernels can be represented as the following infinite series

$$ $\displaystyle\!\!\!k_{\nu,\kappa,\sigma^{2}}(x,x^{\prime})=\frac{\sigma^{2}}{C_{\nu,\kappa}}\sum_{j=0}^{\infty}\Phi_{\nu,\kappa}(\lambda_{j})\phi_{j}(x)\phi_{j}(x^{\prime}),$ $\displaystyle\Phi_{\nu,\kappa}(\lambda)=\begin{cases}\left(\frac{2\nu}{\kappa^{2}}+\lambda\right)^{-\nu-\frac{d}{2}}&\nu<\infty\\ e^{-\frac{\kappa^{2}}{2}\lambda}&\nu=\infty\end{cases}$ (3) $$

where $-\lambda_{j},\phi_{j}$ are the eigenpairs of the Laplace–Beltrami operator on $X$, $d$ is the dimension of $X$, and $C_{\nu,\kappa}$ is a normalisation constant ensuring $\frac{1}{\mathrm{vol}X}\int_{X}k_{\nu,\kappa,\sigma^{2}}(x,x)\mathop{}\!\mathrm{d}x=\sigma^{2}$.
The infinite sum must be truncated for computational tractability; nevertheless, the rapid decay of the coefficients $\Phi_{\nu,\kappa}(\lambda_{j})$ makes this a sensible approximation with convergence guarantees rosa2023 ().

Such models were used, e.g., in robotics jaquier2022 (; ) and medical coveney2020 () applications.
Their vector field counterparts, to which we return later in the paper, were used for modelling wind velocities on the globe hutchinson2021 (; ).

These models tend to perform well, especially when data is scarce and uncertainty quantification is crucial, but they can struggle to capture complex irregular patterns.
One potential way to improve on this is to consider *deep Gaussian processes*, which—in the Euclidean case—we now review.

### 2.2 Deep Gaussian Processes and Approximate Inference

A deep Gaussian process $F$ is a composition $F=f^{L}\circ\cdots\circ f^{1}$ of multiple shallow GPs $f^{l}$ damianou2013 ().
To allow for richer structure, the hidden *layers* $f^{l}$, for $1\leq l\leq L-1$, are typically *vector-valued* GPs $f^{l}:\mathbb{R}^{d}\to\mathbb{R}^{d}$, i.e. vectors of scalar-valued GPs stacked together and potentially correlated with each other alvarez2012 ().
The resulting random function $F$ is itself not a GP.
Thus, even for Gaussian likelihoods $p(\boldsymbol{y}\mid F(\boldsymbol{x}))$, inference for $F\mid\boldsymbol{y}$ is intractable.

To overcome this, various approximate inference techniques for deep GPs were proposed, perhaps the most popular being *doubly stochastic variational inference* salimbeni2017 ().
In it, the intractable posterior $F\mid\boldsymbol{y}$ is approximated in terms of the KL divergence metric $D_{\mathrm{KL}}$ by the elements of a certain *variational family* of tractable distributions.
This family itself consists of deep GPs whose layers are *sparse GPs* titsias2009 (; ) which we now discuss.

Sparse GPs were originally proposed as a variational family for approximating shallow GPs.
For them, approximate inference helps scale to big datasets or accommodate non-Gaussian likelihoods.
Take some $f\sim\mathcal{GP}(\mu,k)$.
A sparse GP $f_{\boldsymbol{z},\boldsymbol{m},\boldsymbol{S}}$ is a family of GPs parameterised by a set of $m$ *inducing locations* $\boldsymbol{z}\subseteq X$, as well as a mean vector $\boldsymbol{m}\in\mathbb{R}^{m}$ and a covariance matrix $\mathbf{S}\in\mathbb{R}^{m\times m}$ which determine the corresponding *inducing variable* distribution $q(\boldsymbol{u})=\mathcal{N}(\boldsymbol{m},\mathbf{S})$.
Specifically,

$$ $\displaystyle p(f_{\boldsymbol{z},\boldsymbol{m},\boldsymbol{S}}(\cdot))={\operatorname*{\mathbb{E}}}_{\boldsymbol{u}\sim q(\boldsymbol{u})}\,p(f(\cdot)\mid\boldsymbol{u},\boldsymbol{z}),$ $\displaystyle q(\boldsymbol{u})=\mathcal{N}(\boldsymbol{m},\mathbf{S}),$ (4) $$

where $p(f(\cdot)\mid\boldsymbol{u},\boldsymbol{z})$ is the prior $f$ conditioned on $f(\boldsymbol{z})=\boldsymbol{u}$.
Intuitively, $\boldsymbol{z}$ are the pseudo-inputs, $\boldsymbol{u}$ are random pseudo-observations, and $p(f_{\boldsymbol{z},\boldsymbol{m},\boldsymbol{S}}(\cdot))$ is a kind of pseudo-posterior.
It is a GP with

$$ $\displaystyle\mu_{\boldsymbol{z},\boldsymbol{m},\boldsymbol{S}}(\cdot)$ $\displaystyle=\mu(\cdot)+k(\cdot,\boldsymbol{z})k(\boldsymbol{z},\boldsymbol{z})^{-1}(\boldsymbol{m}-\mu(\boldsymbol{z})),$ (5) $\displaystyle k_{\boldsymbol{z},\boldsymbol{m},\boldsymbol{S}}(\cdot,\cdot^{\prime})$ $\displaystyle=k(\cdot,\cdot^{\prime})-k(\cdot,\boldsymbol{z})k(\boldsymbol{z},\boldsymbol{z})^{-1}(k(\boldsymbol{z},\boldsymbol{z})-\mathbf{S})k(\boldsymbol{z},\boldsymbol{z})^{-1}k(\boldsymbol{z},\cdot^{\prime}).$ (6) $$

These can be readily generalised to the vector-valued setting, with $\boldsymbol{m}\in\mathbb{R}^{md}$ and $\mathbf{S}\in\mathbb{R}^{md\times md}$.

For a deep GP, the respective variational family is the distribution of a composition of sparse GPs:

$$ $\displaystyle F_{\boldsymbol{\theta}}=f_{\boldsymbol{z}^{L},\boldsymbol{m}^{L},\mathbf{S}^{L}}\circ\ldots\circ f_{\boldsymbol{z}^{1},\boldsymbol{m}^{1},\mathbf{S}^{1}},$ $\displaystyle\boldsymbol{\theta}=\left\{\boldsymbol{z}^{l},\boldsymbol{m}^{l},\mathbf{S}^{l}\right\}_{l=1}^{L}.$ (7) $$

The variational parameters $\boldsymbol{\theta}$ are found by minimising $D_{\mathrm{KL}}(p(F\mid\boldsymbol{y})\,\|\,p(F_{\boldsymbol{\theta}}))$, which is equivalent salimbeni2017 () to maximising the following evidence lower bound (ELBO)

$$ $\mathrm{ELBO}={\sum}_{i=1}^{n}{\operatorname*{\mathbb{E}}}_{F(x_{i})\sim p(F_{\boldsymbol{\theta}}(x_{i}))}\log p(y_{i}|F(x_{i}))-{\sum}_{l=1}^{L}D_{\mathrm{KL}}(q(\boldsymbol{u}^{l})\,\|\,p(\boldsymbol{u}^{l})),$ (8) $$

where $q(\boldsymbol{u}^{l})\sim\mathcal{N}(\boldsymbol{m}^{l},\mathbf{S}^{l})$ and $p(\boldsymbol{u}^{l})\sim\mathcal{N}(\mu^{l}(\boldsymbol{z}^{l}),k^{l}(\boldsymbol{z}^{l},\boldsymbol{z}^{l}))$.
The second term can be computed exactly using the formula for $D_{\mathrm{KL}}$ between two Gaussian vectors.
The first term is intractable, but can be efficiently approximated by drawing a sample from $p(F_{\boldsymbol{\theta}})$—this can be done in a layerwise fashion, as salimbeni2017 () suggest—and subsampling the sum over a mini-batch of inputs $\boldsymbol{x}$.
This way, optimisation proceeds by stochastic gradient descent.

## 3 Residual Deep Gaussian Processes on Manifolds

In this section, we introduce the new model class of *residual deep Gaussian processes on manifolds*.
It generalises the notion of deep GPs to Riemannian manifolds, allowing for the modelling of scalar- and vector-valued functions, vector fields, and functions taking values in the input manifold itself.

### 3.1 The Architecture

Let $X$ be a Riemannian manifold.
The key challenge in building a deep Gaussian process $F$ on $X$ is finding a practical notion of manifold-to-manifold GPs $f^{l}$ to serve as its hidden layers:

$$ $\displaystyle F=f^{L}\circ f^{L-1}\circ\cdots\circ f^{2}\circ f^{1},$ $\displaystyle f^{l}:X\to X\mathrel{\text{missing}}{~{}for~{}}1\leq l\leq L-1,$ (9) $$

where, to simplify exposition, we assume that the last layer $f^{L}$ is real-valued, although it can just as well be $X$-valued, vector-valued, or it can be a vector field, depending on the problem at hand.

While building a GP with inputs in $X$ amounts to finding an appropriate kernel ${k:X\times X\to\mathbb{R}}$, handling outputs in $X$ requires redefining the inherently Euclidean notion of a *Gaussian*.
We aim to circumvent this difficulty.
To explain how, we start by considering the popular Euclidean deep GP architecture of salimbeni2017 ().
There, each layer $f^{l}:\mathbb{R}^{d}\to\mathbb{R}^{d}$ is of the form

$$ $f^{l}(x)=x+g^{l}(x),$ (10) $$

where $g^{l}$ is a zero-mean GP.
That is, each layer
$f^{l}$ displaces its input $x$ by a *residual* vector $g^{l}(x)=f^{l}(x)-x$, the difference to the identity transform, much like a residual connection in neural networks he2016 (), which is modelled by a GP.
On a manifold $X\not=\mathbb{R}^{d}$, when $x\in X$ and $g^{l}(x)\in\mathbb{R}^{d}$, the addition operation in $f^{l}(x)=x+g^{l}(x)$ is undefined.
However, there is a natural generalisation:

$$ $f^{l}(x)=\exp_{x}\left(g^{l}(x)\right).$ (11) $$

Here, $\exp_{x}:T_{x}X\to X$ is the *exponential map*, the canonical mapping of the *tangent space* $T_{x}X$ at $x\in X$—i.e., the linear space of vectors tangent to $X$ at point $x$—back to $X$ itself.
That is, a point $x$ is still displaced by the vector $g^{l}(x)$, but in a geometrically sound manner.

The beauty of [Equation 11](#S3.E11) is that it reduces modelling $f^{l}$ to modelling $g^{l}$.
The latter is vector-valued, and thus compatible with the traditional notion of a Gaussian, making the problem conceptually much simpler.
Still, a major technical difficulty remains: for different inputs $x$, the value $g^{l}(x)$ must lie in the different spaces $T_{x}X$.
The mappings behaving like this are called *vector fields*, and their random Gaussian counterparts are called *Gaussian vector fields*, which we proceed to discuss.

### 3.2 Key Building Blocks: Gaussian Vector Fields

A vector field on a manifold $X$ is a function that takes each $x\in X$ to an element of the tangent space $T_{x}X$ of $X$ at the point $x$.
If $X$ is(^†^†*Nash embedding theorem* proves this is always the case for a large enough ambient space dimension $D\!\in\mathbb{N}$.) a submanifold of $\mathbb{R}^{D}$—like the $2$-sphere is a submanifold of $\mathbb{R}^{3}$—the difference between a vector field and a general vector-valued function on $X$ is very intuitive: in the latter, a vector attached to a point $x\in X$ can point in any direction, while in the former it must always be tangential to the manifold $X$ at $x$.
This difference can be seen in [Figure 2(a)](#S3.F2.sf1), which features a vector-valued function on the left and an actual vector field on the right.

A Gaussian vector field (GVF) can be thus thought of as a vector-valued GP whose outputs always happen to be tangential vectors.
This notion is rigorously formalised in the appendices of hutchinson2021 ().
However, for simplicity, we do not dwell on the formalism here.
Instead, we proceed to discuss three practicable GVF constructions that have been put forward in recent research.

##### Projected GVFs

hutchinson2021 () propose a simple idea, to build a GVF $g$ from any given vector-valued GP $\boldsymbol{h}:X\subset\mathbb{R}^{D}\to\mathbb{R}^{D}$ by *projecting* its outputs onto the appropriate tangent spaces.
Such a projection $P_{(\cdot)}:\mathbb{R}^{D}\rightarrow T_{(\cdot)}X$ exists because, if $X$ is a submanifold of $\mathbb{R}^{D}$, then any tangent space $T_{(\cdot)}X$ can be identified with a linear subspace of $\mathbb{R}^{D}$.
Thus, $g(x)=P_{x}\boldsymbol{h}(x)$ defines a random vector field (see [Figure 2(a)](#S3.F2.sf1)), which turns out to be Gaussian because of the linearity of $P_{(\cdot)}$.

##### Coordinate-frame-based GVFs

Given any *coordinate frame* $\{e_{i}\}_{i=1}^{d}$—that is, a set of functions such that $\{e_{i}(\cdot)\}_{i=1}^{d}$ is a linear basis of $T_{(\cdot)}X$—and a vector-valued GP $\boldsymbol{h}:X\to\mathbb{R}^{d}$ with components $h_{i}:X\to\mathbb{R}$, the equation $g(x)\!=\!\sum_{i=1}^{d}h_{i}(x)e_{i}(x)$ defines a GVF.
This is shown in [Figure 2(b)](#S3.F2.sf2).

##### Hodge GVFs

Most recently, robertnicoud2024 () extended the generalisation of Matérn GPs from [Section 2.1](#S2.SS1) to the setting of vector fields on compact manifolds.
They derive an analogue of [Equation 3](#S2.E3), representing the respective *Hodge Matérn kernels*, as an infinite series(^†^†robertnicoud2024 () notice that sometimes, by analytically summing the terms corresponding to the same eigenvalue $\lambda_{n}$, the series convergence can be sped up by many orders of magnitude, improving efficiency.)

$$ $\boldsymbol{k}_{\nu,\kappa,\sigma^{2}}(x,x^{\prime})=\frac{\sigma^{2}}{C_{\nu,\kappa}}\sum_{j=0}^{\infty}\Phi_{\nu,\kappa}(\lambda_{j})s_{j}(x)\otimes s_{j}(x^{\prime}).$ (12) $$

Here, $s_{j}$ are the eigenfields of the Hodge Laplacian on $X$ that correspond to the eigenvalues $-\lambda_{j}$, $\otimes$ is the tensor product, and $\Phi_{\nu,\kappa}$ is exactly as in [Equation 3](#S2.E3).
This family of kernels can be made more expressive by using different hyperparameters $\sigma^{2}$, $\kappa$, and $\nu$ for different *types* of eigenfields $s_{j}$: the *pure-divergence* $s_{j}$, the *pure-curl* $s_{j}$, and the harmonic $s_{j}$.
The result is called *Hodge-compositional* Matérn kernels robertnicoud2024 (; ).
For example, it can represent the inductive bias of divergence-free vector fields, i.e. vector fields having no “sinks” and “sources”, like the wind velocity field at certain altitudes.
This can be seen in  [Figure 2(c)](#S3.F2.sf3).

Figure: (a) Projected GVF
Refer to caption: /html/2411.00161/assets/figures/gvf_constructions/aa.png

The first two constructions are *universal*.
This means that by choosing an appropriate $\boldsymbol{h}$, one can obtain any possible GVF.
Although this might seem advantageous, this is also a major curse, as it is often unclear which particular $\boldsymbol{h}$ to take to get good inductive biases.
What is more, simple solutions, such as $\boldsymbol{h}$ with IID components, may lead to undesirable artefacts robertnicoud2024 ().
On the other hand, the third construction is *canonical*, in the same way as the Matérn family is canonical in the scalar Euclidean case, and it is based on the same simple and natural inductive biases.

Hodge GVFs seem to be the most attractive building blocks for deep GPs.
However, although they are generally applicable *in theory*, robertnicoud2024 () only provide practical expressions for the eigenfields $s_{n}$ when $X$ is the circle $\mathbb{S}_{1}$, the 2-sphere $\mathbb{S}_{2}$, or any finite product of those.
Thus, for the manifolds that go beyond this simple form, other GVF constructions have to be used.
With this, we finish introducing our models, and proceed to discuss Bayesian inference for them.

### 3.3 Inference

Like their Euclidean counterparts, residual deep GPs constitute complex non-Gaussian priors, making exact Bayesian inference with them impossible.
However, the doubly stochastic variational inference approach described in [Section 2.2](#S2.SS2) is applicable to them after a few adjustments we detail below.
What is more, for compact manifolds, the approach can be further modified to use certain interdomain inducing variables, which, as [Section 4](#S4) shows, tends to offer superior performance.

##### Doubly stochastic variational inference

Consider the analogue of the variational family $p(F_{\boldsymbol{\theta}})$ in [Equation 7](#S2.E7) with the sparse GP layers $f_{\boldsymbol{z}^{l},\boldsymbol{m}^{l},\mathbf{S}^{l}}$ replaced by

$$ $f_{\boldsymbol{z}^{l},\boldsymbol{m}^{l},\mathbf{S}^{l}}(\cdot)=\exp_{(\cdot)}\left(g_{\boldsymbol{z}^{l},\boldsymbol{m}^{l},\mathbf{S}^{l}}(\cdot)\right),$ (13) $$

where $g_{\boldsymbol{z}^{l},\boldsymbol{m}^{l},\mathbf{S}^{l}}$ are *sparse Gaussian vector fields*.
Again, intuitively and in practice, treating the manifold $X$ as a submanifold of $\mathbb{R}^{D}$, GVFs can be thought of, and worked with, as a special kind of vector-valued GPs.
For $g_{\boldsymbol{z}^{l},\boldsymbol{m}^{l},\mathbf{S}^{l}}$, however, $\boldsymbol{z}^{l}=(z^{l}_{1},\ldots,z^{l}_{m})$ with inducing locations $z^{l}_{j}\in X$, and

$$ $\displaystyle\boldsymbol{m}^{l}=(\boldsymbol{m}^{l}_{1},\ldots,\boldsymbol{m}^{l}_{m}),$ $\displaystyle\boldsymbol{m}^{l}_{j}\in T_{z^{l}_{j}}X,$ $\displaystyle\operatorname{im}(\mathbf{S}^{l})\subseteq T_{z^{l}_{1}}X\times\ldots\times T_{z^{l}_{m}}X,$ (14) $$

where $\operatorname{im}$ denotes the *image* of a linear operator, and the last two constraints ensure that the random pseudo-observations $\boldsymbol{u}\sim\mathcal{N}(\boldsymbol{m},\mathbf{S})$ are tangent vectors which lie in the appropriate tangent spaces.

To satisfy the aforementioned constraints during optimisation, one can represent

$$ $\displaystyle\boldsymbol{m}^{l}=P_{\boldsymbol{z}^{l}}\widetilde{\boldsymbol{m}}^{l},$ $\displaystyle\mathbf{S}^{l}=P_{\boldsymbol{z}^{l}}\widetilde{\mathbf{S}}^{l}P_{\boldsymbol{z}^{l}}^{\top},$ $\displaystyle P_{\boldsymbol{z}^{l}}=P_{z^{l}_{1}}\oplus\ldots\oplus P_{z^{l}_{m}},$ (15) $$

with arbitrary $\widetilde{\boldsymbol{m}}^{l}\in\mathbb{R}^{mD}$, arbitrary positive semi-definite $\widetilde{\mathbf{S}}^{l}\in\mathbb{R}^{mD\times mD}$, and $P_{z}:\mathbb{R}^{D}\to T_{z}X$ denoting the projection of vectors in $\mathbb{R}^{D}$ onto the tangent space $T_{z}X$, as discussed in [Section 3.2](#S3.SS2).

Instead of using such representations, one can fix a (locally) smooth frame and optimise the coefficients of $\boldsymbol{m}^{l}$ and $\mathbf{S}^{l}$ represented in this frame.
In any case, one needs to make sure $z^{l}_{j}$ always remain on the manifold $X$ during optimisation, which can be done by using specialised libraries, such as Pymanopt townsend2016 () or Geoopt kochurov2020 ().
For low-dimensional manifolds, a fixed grid on $X$ or a set of cluster centroids can be an effective alternative to optimising $z^{l}_{j}$.

Finally, to approximate ELBO, we need to sample $F(x_{i})\sim p(F_{\boldsymbol{\theta}}(x_{i}))$.
As in salimbeni2017 (), this can be done sequentially.
Specifically, if $\hat{F}(x_{i})$ denotes the desired sample, then

$$ $\displaystyle\hat{F}(x_{i})=\hat{f}^{L}_{i},$ $\displaystyle\hat{f}^{l}_{i}=\exp_{\hat{f}^{l-1}_{i}}(g_{\boldsymbol{z}^{l},\boldsymbol{m}^{l},\mathbf{S}^{l}}(\hat{f}^{l-1}_{i}))\mathrel{\text{missing}}{~{}for~{}}1\leq l\leq L,$ $\displaystyle\hat{f}^{0}_{i}=x_{i},$ (16) $$

and, given $\hat{f}^{l-1}_{i}$, each individual $g_{\boldsymbol{z}^{l},\boldsymbol{m}^{l},\mathbf{S}^{l}}(\hat{f}^{l-1}_{i})$ can be sampled in the usual manner.

##### Interdomain inducing variables on manifolds

On compact manifolds, an alternative variational family can be used that speeds up the inference and can often lead to better predictive performance.
It is constructed by replacing the inducing locations $\boldsymbol{z}$ by *inducing linear functionals* $\boldsymbol{\zeta}\!=\!(\zeta_{1},\dots,\zeta_{m})$.
Each $\zeta_{j}$ takes in a vector field $g$ and outputs a real number.
These $\boldsymbol{\zeta}$ define a sparse GVF through

$$ $\displaystyle p(g_{\boldsymbol{\zeta},\boldsymbol{m},\boldsymbol{S}}(\cdot))={\operatorname*{\mathbb{E}}}_{\boldsymbol{u}\sim q(\boldsymbol{u})}\,p(g(\cdot)\mid\boldsymbol{u},\boldsymbol{\zeta}),$ $\displaystyle q(\boldsymbol{u})=\mathcal{N}(\boldsymbol{m},\mathbf{S}),$ (17) $$

where $p(g(\cdot)\mid\boldsymbol{u},\boldsymbol{\zeta})$ is the prior $g$ conditioned on $\boldsymbol{\zeta}(g)=\boldsymbol{u}$, where $\boldsymbol{\zeta}(g)=(\zeta_{1}(g),\ldots,\zeta_{m}(g))$.
For example, linear functionals of the form $\zeta(g)\!\!\!=\!\!\!\left\langle g(z),e_{i}(z)\right\rangle_{T_{z}X}$—here, $\{e_{i}\}_{i=1}^{d}$ is a coordinate frame—can be used to recover the usual doubly stochastic variational inference considered above.

The mean and covariance of $g_{\boldsymbol{\zeta},\boldsymbol{m},\boldsymbol{S}}$ are given by [Equation 5](#S2.E5), with $\boldsymbol{z}$ replaced by $\boldsymbol{\zeta}$ and

$$ $\displaystyle k(\boldsymbol{\zeta},\cdot)=\operatorname{Cov}\left(\boldsymbol{\zeta}(g),g(\cdot)\right),$ $\displaystyle k(\boldsymbol{\zeta},\boldsymbol{\zeta}^{\prime})=\operatorname{Cov}\left(\boldsymbol{\zeta}(g),\boldsymbol{\zeta}(g)\right),$ (18) $$

see, for example, lazarogredilla2009 (; ).

Now, if the kernel of $g$ can be expressed as $\sum a_{j}\phi_{j}(x)\otimes\phi_{j}(x^{\prime})$ where $\left\{\phi_{j}\right\}$ is an orthonormal basis—this is obviously so for Hodge GVFs, but is also often the case for other GVFs on compact manifolds—the inducing functionals $\zeta_{j}(\cdot)=\left\langle\cdot,\phi_{j}\right\rangle_{L^{2}}/a_{j}$ yield very simple covariance matrices

$$ $\displaystyle k(\zeta_{j},\cdot)=\phi_{j}(\cdot),$ $\displaystyle k(\zeta_{i},\zeta_{j})=\delta_{i,j}/a_{i}.$ (19) $$

In particular, $k(\boldsymbol{\zeta},\boldsymbol{\zeta})$ is diagonal, making it trivial to invert.
dutordoir2020 () report that this can yield significant acceleration in practice.
For residual deep GPs, this phenomenon affects every individual layer, thus making the cumulative effect even more pronounced.
We refer the reader to [Appendix B](#A2) for further practical and theoretical considerations regarding this variational family.

##### Posterior mean, variance, and samples

Expectations $\operatorname*{\mathbb{E}}F_{\boldsymbol{\theta}}(x)$ and variances $\operatorname{Var}F_{\boldsymbol{\theta}}(x)$ of the approximate posterior $F_{\boldsymbol{\theta}}$ cannot be computed exactly.
Instead, they can be estimated by appropriate Monte Carlo averages, with [Equation 16](#S3.E16) providing a way to sample $F_{\boldsymbol{\theta}}(x)$.
However, since these estimates ignore the correlation between $F_{\boldsymbol{\theta}}(x)$ and $F_{\boldsymbol{\theta}}(x^{\prime})$, they are not continuous as functions of $x$.
When continuity or differentiability of $\operatorname*{\mathbb{E}}F_{\boldsymbol{\theta}}(x)$ and $\operatorname{Var}F_{\boldsymbol{\theta}}(x)$ are desirable, another method can be used.
The key idea in this case is to draw (approximate) samples from $F_{\boldsymbol{\theta}}(\cdot)$ which happen to be actual functions, for example linear combinations of some analytic basis functions or compositions of such.
This can be done by applying the *pathwise conditioning* of wilson2020 (; ) in a sequential manner, akin to [Equation 16](#S3.E16).
This approach is useful for visualisation, performance metric estimation, and for working with downstream quantities, such as acquisition functions in Bayesian optimisation, for which differentiability is key for efficiently finding their maxima.

## 4 Experiments

We begin this section by examining how various GVF and variational family choices impact the regression performance of residual deep GPs in synthetic experiments, as discussed in [Section 4.1](#S4.SS1).
Throughout, we compare our models to a baseline with Euclidean hidden layers.
Next, in the robotics-inspired experiments of [Section 4.2](#S4.SS2), we demonstrate that residual deep GPs can significantly enhance Bayesian optimisation on a manifold when the optimised function is irregular.
Following this, in [Section 4.3](#S4.SS3), we show state-of-the-art predictive and uncertainty calibration performance of residual deep GPs in wind velocity modelling on the globe, achieving interpretable patterns even at low altitudes where data is more complex and irregular.
Finally, in [Section 4.4](#S4.SS4), we explore potential avenues for using residual deep GPs to accelerate inference on inherently Euclidean data.

### 4.1 Synthetic Examples

Figure: Figure 3: NLPD of different residual deep GP variants and the baseline model, on the regression problem for the synthetic benchmark function visualised in [Figure 4(a)](#S4.F4.sf1). Different subplots correspond to different training set sizes $N$. The solid lines represent the mean, while the shaded areas represent the $\pm 1$ standard deviation region around it. All statistics are computed over $5$ randomised runs.
Refer to caption: /html/2411.00161/assets/x1.png

##### Setup

Deep GPs bear the promise of outperforming their shallow counterparts in modelling complex, irregular functions.
To test this, we construct a benchmark function $f^{*}$ on the $2$-sphere $\mathbb{S}_{2}$ with multiple singularities, which is visualised in [Figure 4(a)](#S4.F4.sf1).
We take $N\in\left\{100,200,400,800,1600\right\}$ training inputs $\boldsymbol{x}$ on a Fibonacci lattice on $\mathbb{S}_{2}$ and put $\boldsymbol{y}=f^{*}(\boldsymbol{x})+\boldsymbol{\varepsilon}$, $\boldsymbol{\varepsilon}\sim\mathcal{N}(\boldsymbol{0},10^{-4}\mathbf{I})$.
Then, we regress $f^{*}$ from $\boldsymbol{x}$ and $\boldsymbol{y}$.
On this problem, we compare different modifications of the residual deep GPs amongst themselves and to a baseline, in terms of the negative log predictive density (NLPD) and the mean squared error (MSE) metrics on the test set of $5000$ points, also on a Fibonacci lattice.
All runs are conducted $5$ times with different random seeds for the observation noise $\boldsymbol{\varepsilon}$.

##### Models

The baseline is a deep GP with Euclidean (rather than manifold-to-manifold) layers.
It is constructed by composing a vector-valued Matérn GP whose signature is $X\to\mathbb{R}^{3}$ with the Euclidean deep GP of salimbeni2017 () on $\mathbb{R}^{3}$.
For the residual deep GPs, we consider two different types of GVFs, projected and Hodge, and two types of variational families, the one based on inducing locations (IL) and the one based on interdomain variables (IV).
To ensure comparability, we match the number of optimised parameters between models as closely as possible.

##### Results

The NLPD values are presented in [Figure 3](#S4.F3). The MSE values exhibit the same trends and can be found in [Figure 8](#A1.F8) in [Appendix A](#A1).
We observe three key patterns.
First, residual deep GPs are never worse than their shallow counterparts, recovering the single-layer solution when data is sparse.
Second, as data becomes more abundant and thus captures more complexity of $f^{*}$, residual deep GPs outperform the shallow GPs.
Third, the IV variational family almost always improves over the IL one, and the best model—with considerable margin—is obtained by combining Hodge GVFs with the IV variational family.
The residual deep GP based on projected GVFs and using the IV variational family is the second best model, which still outperforms the baseline in most cases.

### 4.2 Geometry-aware Bayesian Optimisation

Figure: (a) Irregular benchmark function.
Refer to caption: /html/2411.00161/assets/figures/synthetic/irregular_function.png

##### Motivation

GPs are a widely used model class for Bayesian optimisation, a technique for optimising expensive-to-evaluate black-box functions that leverages uncertainty estimates to balance exploration and exploitation shahriari2016 ().
In robotics, such problems arise, for example, when a control policy needs to be fine-tuned to a specific real-world environment.
This task was shown to benefit from treating the optimisation space as a manifold and using geometry-aware Gaussian processes to drive Bayesian optimisation jaquier2022 ().
The functions upon which the technique was tested are rather regular, which is not always the case in reality, especially when dealing with increasingly complex systems.
Motivated by this challenge, we explore if residual deep GPs can offer improved performance in optimising complex irregular functions on manifolds.

##### Setup

We consider two target functions to optimise.
The first is the irregular function from [Section 4.1](#S4.SS1), visualised in [Figure 4(a)](#S4.F4.sf1), modified to have only one global minimum.
The second is the much more regular Ackley function, projected onto $\mathbb{S}_{3}$, one of the benchmarks in jaquier2022 ().
In each Bayesian optimisation run, we perform the first 180 iterations using a shallow geometry-aware GP, followed by 20 iterations using a residual deep GP—both employing the *expected improvement* acquisition function (see, e.g., frazier2018 ()).
In this experiment, we showcase the coordinate-frame-based GVFs, as described in [Appendix A](#A1).
We do not use deep GPs in the initial iterations because, as intuition suggests and [Section 4.1](#S4.SS1) affirms, deep GPs start outperforming their shallow counterparts only when data becomes more abundant.
Although deep GPs show comparable performance even for small datasets, training them is more computationally demanding, making their use less efficient in the early stages of optimisation.
We repeat each run 15 times to account for the stochasticity of initialisation, optimisation of the acquisition function, and training of GP models.

##### Results

The optimisation performance, measured in terms of the logarithm of regret, is reported in [Figure 4(b)](#S4.F4.sf2).
We find that residual deep GPs significantly improve performance in the Bayesian optimisation of the irregular function.
Specifically, switching to a residual deep GP in the latter stages of optimisation greatly, and often immediately, reduces the gap between the true optimum and the found optimum.
This trend is consistent across most runs, with only one outlier showing no improvement due to insufficient data collection near the singularity during the initial shallow GP phase.
In contrast, for the Ackley function, we observe no substantial difference in performance between the two methods: both approaches replicate results from jaquier2022 (), with nearly identical median regret trajectories.
This outcome aligns with our expectations, as the region around the minimum of the Ackley function is smooth and can be effectively modelled by a shallow GP.

### 4.3 Wind interpolation on the globe

Figure: (a) The ground truth wind velocities as black arrows and the training locations along the Aeolus satellite track as red points.
Refer to caption: /html/2411.00161/assets/figures/wind/aeolus_track_and_ground_truth_3d_1500pts.png

##### Motivation

Non-Euclidean geometry has a particularly pronounced effect on vector fields, such as wind velocity fields on the globe, more so than on scalar functions.
For instance, the famous hairy ball theorem states that a smooth vector field must always have a zero somewhere on the $2$-sphere, i.e. there must always be a location where wind does not blow.
Wind interpolation is thus an attractive use case for geometry-aware Gaussian vector fields, where they have been shown to perform well hutchinson2021 (; ).
Here, we show that residual deep GPs can improve the performance of probabilistic wind velocity models when the data contains complex irregular patterns, which naturally occur in wind fields at lower altitudes.

##### Setup

We consider the task of interpolating the monthly average wind velocity from the ERA5 dataset era5 (), from a set of locations on the Aeolius satellite track reitebuch2012 (), simulating a practical setting of a weather-analysing satellite.
We use ERA5 data from January to December 2010, as in robertnicoud2024 (), and sample the Aeolius track locations every minute for a 24-hour period from 9:00 am, January 1st, 2019.
We choose a 24-hour period instead of a 1-hour period, as in hutchinson2021 (), because in that time frame, the satellite produces a denser set of observations, crucial for capturing the complexity of wind behaviour at low altitudes.
The ground truth vector field and the input locations are visualised in [Figure 5(a)](#S4.F5.sf1).
To assess how decreasing regularity of data associated with decreasing altitude affects our model, we consider data at three altitudes: approximately 5.5 km, 2 km, and 0.1 km.
In our models, we use Hodge GVFs in hidden layers and as the last layer, and interdomain inducing variables for inference.

##### Results

We report regression performance, in terms of NLPD, in [Figure 6](#S4.F6); MSE follows similar trends and can be found in [Figure 9](#A1.F9) in [Appendix A](#A1).
We find that residual deep GPs improve upon the state-of-the-art shallow Hodge GVFs (1-layer models in the plots, robertnicoud2024 ())

Figure: Figure 6: NLPD of residual deep GPs on the wind modelling task across three altitude levels. Solid lines give the mean NLPD, while the shaded regions around it span $\pm 1$ standard deviation. Statistics were computed from 12 runs—one for each month of 2010.
Refer to caption: /html/2411.00161/assets/x3.png

both in prediction quality and uncertainty calibration, as evidenced by the significantly lower NLPD and MSE.
Furthermore, [Figure 5(b)](#S4.F5.sf2) shows that the uncertainty estimates of our deep model at the lowest altitude are interpretable.
Indeed, regions of high uncertainty follow regions of irregular wind currents, such as boundaries where multiple currents meet or continental boundaries, as well as areas of seasonally high winds, such as India during the peak monsoon season.
At the same time, low uncertainty is assigned to regions with constant-like currents.
This is unachievable for shallow GVFs since their posterior covariance depends only on the *locations* of the observations, which are rather uniformly dense in our setup, and not on the observations themselves.
Additionally, [Figure 5(b)](#S4.F5.sf2) shows predictive uncertainty to be well-calibrated, with areas of highest error corresponding to regions of high uncertainty.

### 4.4 Accelerating Inference for Euclidean Data

Figure: Figure 7: Wall clock time taken by one training step of Euclidean deep GPs with inducing locations, and residual deep GPs with interdomain variables. We consider 5 UCI datasets, with dimension $d$ and batch size $B$. Solid lines show the mean, computed by averaging over 100 training steps, while the shaded areas span $\pm 1$ standard deviation. However, they are often too narrow to be visible.
Refer to caption: /html/2411.00161/assets/x4.png

##### Motivation

Inspired by their connection to infinitely wide neural networks, dutordoir2020 () showed that geometry-aware GPs on hyperspheres can be applied to inherently Euclidean data to accelerate variational inference.
Specifically, they reported that approximate variational inference using the shallow analogue of the interdomain inducing variables applied to data mapped(^†^†Mapping $\mathbb{R}^{d}\ni\boldsymbol{x}\mapsto(x_{1},..,x_{d},b)/\left\lVert(x_{1},..,x_{d},b)\right\rVert\in\mathbb{S}_{d+1}$, where $b\in\mathbb{R}$ is a *bias* term, $\boldsymbol{x}=(x_{1},..,x_{d})$.) from $\mathbb{R}^{d}$ to the proxy manifold of $\mathbb{S}_{d+1}$ can be significantly faster than inducing-location-based variational inference for a Euclidean GP on the original data, while achieving competitive predictive performance.
We investigate whether this result can be extended to the case of deep Gaussian processes.

##### Setup

As Euclidean data, we use the same UCI datasets as dutordoir2020 (), and we use the same mapping from $\mathbb{R}^{d}$ to the proxy manifold $\mathbb{S}_{d+1}$.
We use projected GVFs to accommodate arbitrarily-dimensional hyperspheres.
Overall, our experimental setup follows dutordoir2020 (), except that working with deep models required us to optimise ELBO directly instead of marginalising out the variational mean $\boldsymbol{m}$ and covariance $\mathbf{S}$ as in titsias2009 ().
Additionally, because of high memory requirements of L-BFGS nocedal1980 () arising from the lack of marginalisation and depth, we switch to Adam kingma2015 ().
We use a single Intel i7-13700H CPU.

##### Results

We compare the variational inference speed, measured by wall-clock time for a single optimisation step, in [Figure 7](#S4.F7).
The advantage for shallow 1-layer GPs increases significantly with more layers, offering a considerable edge in deep models.
However, predictive performance comparisons in [Figure 12](#A1.F12) in [Appendix A](#A1) do not show such an optimistic picture: Euclidean deep GPs always outperform residual deep GPs with the same number of layers in terms of NLPD and MSE.
Possibly, this might be due to the aforementioned differences in optimisation.
Also, our choice of the mapping from $\mathbb{R}^{d}$ to a proxy manifold and the choice of the proxy manifold itself might be overly simplistic and thus hinder performance.
We hypothesise that better mappings or optimisation could potentially make the tested approach a more efficient alternative to Euclidean deep GPs.
Achieving this, however, will require further work, which is beyond the scope of this paper.

## 5 Conclusion

In this paper, we proposed a novel model class of *residual deep Gaussian processes* on manifolds.
We reviewed practical Gaussian vector field constructions for building their hidden layers and discussed two variational inference techniques, including one specifically tailored to the structure of Gaussian vector fields on compact manifolds and based on interdomain inducing variables.
We evaluated our models in synthetic experiments, examining the impact of Gaussian vector field and variational family choices.
These experiments supported favouring Hodge Gaussian vector fields and interdomain inducing variables.
In a robotics-motivated experiment, our models significantly enhanced Bayesian optimisation for an irregular function on the sphere.
For probabilistic interpolation of wind velocities, we achieved state-of-the-art performance, surpassing the recently proposed shallow Hodge Gaussian vector fields.
Finally, we showed interdomain inducing variables to be superior in terms of inference time, compared to doubly stochastic variational inference for Euclidean deep Gaussian processes.
This indicates potential future benefits for Euclidean data if suitable mappings from manifold data to proxy manifolds are found.
We believe residual deep Gaussian processes will provide a powerful toolset for applications in climate modelling, robotics, and beyond.

#### Acknowledgments

VB was supported by ELSA (European Lighthouse on Secure and Safe AI) funded by the European Union under grant agreement No. 101070617.
The Blender rendering scripts we used for plotting were adapted from terenin2022 ().

## Appendix A Additional Experimental Details

### A.1 Implementation

##### Efficient kernel evaluation with the addition theorem

In our implementation of manifold Matérn kernels, we utilise the *addition theorem for spherical harmonics* devito2021 (; ) to accelerate kernel computation.(^†^†The analogues of this theorem hold for many other manifolds, see, e.g., azangulov2022 ().)
On $\mathbb{S}_{d}$, the eigenfunctions of the Laplace–Beltrami operator are known to be certain special functions called *spherical harmonics*.
The addition theorem gives a relation between all spherical harmonics corresponding to the $(k+1)$-th largest eigenvalue of the negative Laplace–Beltrami operator $-\Delta$, denoted $\{\phi_{k,j}\}_{j=1}^{J}$, and Gegenbauer polynomials $C_{k}^{(\alpha)}$, other special functions:

$$ $\sum_{j=1}^{J}\phi_{k,j}(x)\phi_{k,j}(x^{\prime})=c_{k,d}C_{k}^{(\alpha)}(x\cdot x^{\prime})$ (20) $$

with the dot product computed after embedding $\mathbb{S}_{d}$ in $\mathbb{R}^{d+1}$ as the unit sphere centred at the origin, $\alpha=\frac{d-1}{2}$, and $c_{k,d}$ some known absolute constants.
Thus, when computing the scalar Matérn kernel on $\mathbb{S}_{d}$, we truncate the infinite sum in [Equation 3](#S2.E3) to include all spherical harmonics up to the $(K+1)$-th eigenvalue, and apply [Equation 20](#A1.E20).
This gives the following formula

$$ $k_{\nu,\kappa,\sigma^{2}}(x,x^{\prime})=\frac{\sigma^{2}}{C_{\nu,\kappa}}\sum_{k=0}^{K}\Phi_{\nu,\kappa}(\lambda_{k})c_{k,d}C_{k}^{(\alpha)}(x\cdot x^{\prime}).$ (21) $$

In case of the Hodge Matérn kernel on $\mathbb{S}_{2}$, we also apply addition theorem, only noting that here $\phi_{k,j}$ are replaced with the normalised vector spherical harmonics: $\frac{\nabla\phi_{k,j}}{\sqrt{\lambda_{k}}}$ for the divergence-only kernel, and $\frac{\star\nabla\phi_{k,j}}{\sqrt{\lambda_{k}}}$ for the curl-only kernel, where $\star$ is the Hodge star operator robertnicoud2024 ().
In this case, $c_{k,d}=\frac{k+\alpha}{\alpha}$.
Thus,

$$ $\displaystyle\boldsymbol{k}^{\mathrm{div}}_{\nu,\kappa,\sigma^{2}}(x,x^{\prime})$ $\displaystyle=\frac{\sigma^{2}}{C^{\mathrm{div}}_{\nu,\kappa}}\sum_{k=0}^{K}\frac{\Phi_{\nu,\kappa}(\lambda_{k})}{\lambda_{k}}\frac{k+\alpha}{\alpha}(\nabla_{x}\otimes\nabla_{x^{\prime}})C_{k}^{(\alpha)}(x\cdot x^{\prime})$ (22) $\displaystyle\boldsymbol{k}^{\mathrm{curl}}_{\nu,\kappa,\sigma^{2}}(x,x^{\prime})$ $\displaystyle=\frac{\sigma^{2}}{C^{\mathrm{curl}}_{\nu,\kappa}}\sum_{k=0}^{K}\frac{\Phi_{\nu,\kappa}(\lambda_{k})}{\lambda_{k}}\frac{k+\alpha}{\alpha}(\star\nabla_{x}\otimes\star\nabla_{x^{\prime}})C_{k}^{(\alpha)}(x\cdot x^{\prime}).$ (23) $$

The Hodge-compositional kernel, which we typically use, is the sum $\boldsymbol{k}^{\mathrm{div}}_{\nu,\kappa_{1},\sigma_{1}^{2}}+\boldsymbol{k}^{\mathrm{curl}}_{\nu,\kappa_{2},\sigma_{2}^{2}}$.

##### Accelerated training with whitened inducing variables

In all the models we test, approximate inference requires a variational mean vector $\boldsymbol{m}$ and a variational covariance matrix $\mathbf{S}$ that parameterise $q(\boldsymbol{u})=\mathcal{N}(\boldsymbol{m},\mathbf{S})$.(^†^†$\mathbf{S}$ is parametrised by its lower Cholesky factor to ensure positive definiteness during optimisation.)
However, to accelerate convergence during training, instead of working with the inducing variables $\boldsymbol{u}$ directly, we work with *whitened* inducing variables $\boldsymbol{u}^{\prime}=\mathbf{L}^{-1}\boldsymbol{u}$, where $\mathbf{L}$ is the lower Cholesky factor of $k(\boldsymbol{z},\boldsymbol{z})$ or $k(\boldsymbol{\zeta},\boldsymbol{\zeta})$ matthews2016 ().
Thus, in practice, denoting the whitened variational mean and covariance of $q(\boldsymbol{u}^{\prime})$ as $\boldsymbol{m}^{\prime}$ and $\mathbf{S}^{\prime}$, we use a modified version of [Equation 5](#S2.E5):

$$ $\displaystyle\mu_{\boldsymbol{z},\boldsymbol{m^{\prime}},\boldsymbol{S^{\prime}}}(\cdot)$ $\displaystyle=\mu(\cdot)+k(\cdot,\boldsymbol{z})k(\boldsymbol{z},\boldsymbol{z})^{-1}(\mathbf{L}\boldsymbol{m}^{\prime}-\mu(\boldsymbol{z}))$ (24) $\displaystyle=\mu(\cdot)+k(\cdot,\boldsymbol{z})\mathbf{L}^{-T}(\boldsymbol{m}^{\prime}-\mathbf{L}^{-1}\mu(\boldsymbol{z})),$ (25) $\displaystyle k_{\boldsymbol{z},\boldsymbol{m}^{\prime},\boldsymbol{S}^{\prime}}(\cdot,\cdot^{\prime})$ $\displaystyle=k(\cdot,\cdot^{\prime})-k(\cdot,\boldsymbol{z})k(\boldsymbol{z},\boldsymbol{z})^{-1}(k(\boldsymbol{z},\boldsymbol{z})-\mathbf{L}\mathbf{S}^{\prime}\mathbf{L}^{T})k(\boldsymbol{z},\boldsymbol{z})^{-1}k(\boldsymbol{z},\cdot^{\prime})$ (26) $\displaystyle=k(\cdot,\cdot^{\prime})-k(\cdot,\boldsymbol{z})\mathbf{L}^{-T}(\mathbf{I}-\mathbf{S^{\prime}})\mathbf{L}^{-1}k(\boldsymbol{z},\cdot^{\prime}).$ (27) $$

### A.2 Models

##### Mean and kernel

In all models, we equip constituent GPs with a zero mean and an appropriate variant of the Matérn kernel, initialised with length scale $\kappa=1$.
For kernels of output layers we initialise the variance to $\sigma^{2}=1.0$, while for kernels in hidden layers of an $L$-layer deep GPs we set $\sigma^{2}=\frac{10^{-4}}{L-1}$ at the start of training.
In [Section 4.1](#S4.SS1), [Section 4.4](#S4.SS4), and [Section 4.3](#S4.SS3) we initialise the smoothness parameter to $\nu=\frac{3}{2}$, while in section [Section 4.2](#S4.SS2) we set it to $\nu=\frac{5}{2}$ to replicate the setup in jaquier2022 ().
We optimise the smoothness of the manifold Matérn kernels during training, owing to their differentiability with respect to $\nu$, except in [Section 4.2](#S4.SS2), where we fix $\nu$ to match the setup in jaquier2022 ().
Wherever we employ Hodge GVFs we use the Hodge compositional kernel, with separate $\nu,\kappa,\sigma^{2}$ for the curl-free and divergence-free parts.
In models utilising interdomain variables, we use the same number of spherical harmonics for the kernel and inducing variables, as per our discussion in [Appendix B](#A2).

##### Vector-valued GPs

We model vector-valued GPs as a set of *independent* scalar-valued GPs stacked into a vector.
We utilise this construction in Euclidean deep GPs and in residual deep GPs with projected GVFs.

##### Inducing locations

Following salimbeni2017 (), for all models utilising the variational family based on inducing locations $\boldsymbol{z}^{l}$, we initialise $\boldsymbol{z}^{l}$ for every layer to be the centers of the clusters found via k-means clustering of training data.
In residual deep GPs, we further project these locations onto the sphere, and we do not optimise them during training.
In Euclidean deep GPs, we do not normalise the inducing locations and optimise them jointly with all other parameters.

##### Approximation in training and evaluation of deep models

In all experiments, to approximate the ELBO in deep models during training, we use $3$ samples from the posterior.
In evaluation, we use $10$ samples from the posterior to approximate the MSE and NLPD.
In visualisation, we also use ten samples from the posterior to approximate the predictive mean and standard deviation.

### A.3 Synthetic Experiments

Figure: Figure 8: MSE of different residual deep GP variants and the baseline model, on the regression problem for the synthetic benchmark function in [Figure 4(a)](#S4.F4.sf1). Different subplots correspond to different training set sizes $N$. The solid lines represent the mean, while the shaded areas represent the $\pm\sigma$ region around it, where $\sigma$ is the standard deviation, and all statistics were computed over $5$ randomised runs.
Refer to caption: /html/2411.00161/assets/x5.png

##### Data

To examine the influence on data density on model performance in a controlled manner, we generate the training sets as approximately uniform grid of points on the 2-sphere $\mathbb{S}_{2}$.
This is done using the Fibonacci lattice method, which, for a grid of $n$ points, gives the colatitude and longitude of the $i$-th point as

$$ $\displaystyle\mathrm{colatitude}=\arccos\left(1-\frac{2i+1}{n}\right),$ $\displaystyle\mathrm{longitude}=\frac{2\pi n}{\phi},$ $\displaystyle\phi=\frac{1+\sqrt{5}}{2}.$ (28) $$

Because this method gives consistent coverage of $\mathbb{S}_{2}$ we also use it to generate the test set of 5000 points.
We define the benchmark function $f^{\ast}:\mathbb{S}_{2}\rightarrow\mathbb{R}$ by

$$ $f^{\ast}(x)=(Y_{2,3}\circ\varphi)(x)+(Y_{1,2}\circ\varphi\circ R)(x),$ (29) $$

where

$$ $\displaystyle Y_{2,3}(\theta,\phi)$ $\displaystyle=\sqrt{\frac{105}{32\pi}}\sin^{3}\theta\sin(3\phi),$ (30) $\displaystyle Y_{1,2}(\theta,\phi)$ $\displaystyle=\sqrt{\frac{15}{8\pi}}\sin\theta\sin(2\phi)$ (31) $\displaystyle R(x)$ $\displaystyle=(x_{1},-x_{3},x_{2})$ (32) $\displaystyle\varphi(x)$ $\displaystyle=(\operatorname{atan2}(x_{2},x_{1}),\arccos(x_{3})),$ (33) $$

with $x$ being an element of $\mathbb{S}_{2}$ embedded into $\mathbb{R}^{3}$ as the unit sphere centred at the origin.
$Y_{2,3},Y_{1,2}$ are spherical harmonics, smooth functions of their parameters.
Singularities in $f^{\ast}$ are caused by the composition with $\varphi$, which converts $x$ from Cartesian to spherical coordinates, but swaps the positions of the colatitude $\theta$ and longitude $\phi$.
The function $Y_{2,3}\circ\varphi$ has singularities at the poles, while the function $Y_{1,2}\circ\varphi\circ R$ has singularities around the equator.

##### Variational parameters

To make the comparison between models as fair as possible, we set the number of inducing variables for each model in such a way that all models have almost the same total number of optimisable parameters.
In projected GVFs of the residual deep GPs and in vector-valued GPs of the baseline models, we equip each scalar GP with 49 inducing variables, which corresponds to spherical harmonics up to the 7-th negative eigenvalue.
For Hodge GVFs, we use 70 interdomain inducing variables, corresponding to vector spherical harmonics up to the 6-th negative eigenvalue of the Hodge Laplacian.
Despite this, residual deep GPs with Hodge GVFs are still at a disadvantage: a single Hodge GVF has $2555$ variational parameters and $2\cdot 3$ kernel parameters, a single projected GVF and manifold vector-valued GP has $3822$ variational parameters and $3\cdot 3$ kernel parameters, while one Euclidean vector-valued GP has $3822$ variational parameters and $3\cdot 4$ kernel parameters (we use automatic relevance determination (ARD), we get $4$ as the sum of $3$ length scale parameters and $1$ prior variance parameter).
Furthermore, we optimise the inducing locations in the Euclidean vector-valued GPs of the baseline model, so that each of them have an additional $3\cdot 49$ optimisable parameters.

##### Training and evaluation

We optimise all models using the Adam optimiser kingma2015 () for 1000 iterations with learning rate set to $0.01$.

##### Additional results

The MSE comparison is presented in [Figure 8](#A1.F8).

### A.4 Geometry-aware Bayesian optimisation

Figure: Figure 9: MSE of residual deep GPs on the wind modelling task across three altitude levels. Solid lines give the mean MSE, while the shaded regions around it span $\pm 1$ standard deviation. Statistics are computed from 12 runs—one for each month of 2010.
Refer to caption: /html/2411.00161/assets/x6.png

##### Data

To obtain an irregular function $f^{\ast}$ on $\mathbb{S}_{2}$, we modify the target function from [Section 4.1](#S4.SS1) to have only one global minimum near a singularity point.
Specifically, $f^{\ast}$ was defined by

$$ $f^{\ast}(x)=(Y_{2,3}\circ\varphi)(x)\cdot(x_{3}+1)\cdot(1-\arccos(x_{3})),$ (34) $$

where $Y_{2,3}$ and $\varphi$ are as in [Section A.3](#A1.SS3).
The absence of $Y_{1,2}$ removes the singularities around the equator, while the added scaling factors create a minimum at the north pole and ensure it is global.

##### Models

After a preliminary examination, we found that differences between GVFs variants on this task were not significant.
Nevertheless, [Figure 4(b)](#S4.F4.sf2) we present the models that performed best in this examination: in the left subplot, a 2-layer residual deep GP using coordinate-frame GVFs and inducing locations; in the right subplot, a 3-layer residual deep GP using projected GVFs and inducing locations.

##### Optimisation process

Replicating the setup of jaquier2022 (), we begin each Bayesian optimisation by sampling 5 initial observations uniformly at random on the hypersphere.
At each step, we minimise the expected improvement acquisition function using the first-order geometry-aware gradient optimisation implemented in Pymanopt townsend2016 ().
We approximate the expected improvement acquisition function for deep models using Monte Carlo averages driven by pathwise sampling, as described at the very end of [Section 3.3](#S3.SS3).
After each optimisation step we reinitialise the model and fit it to the data using the Adam optimiser kingma2015 () for 500 iterations with learning rate set to $0.01$.

### A.5 Wind interpolation on the globe

##### Data

To each location sampled along the track of the Aeolus satellite, we assign the wind velocity from the closest location present in the ERA5 dataset.
Our test set is a grid of 5000 points, same as in [Section A.3](#A1.SS3).
To each location in the test set, we also assign the wind velocity from the closest location in the ERA5 dataset.

##### Variational parameters

For each GVF within the layers of the tested models we use $198$ interdomain inducing variables.
They correspond to all vector spherical harmonics up to, and including, to the $10$th negative eigenvalue of the Hodge Laplacian.
This choice is arbitrary and simply serves to balance quality of fit with training time.

##### Training and evaluation

We fit the models to data using the Adam optimiser for 1000 iterations with the learning rate set to $0.01$.
To evaluate the models, we compute the MSE and NLPD via Monte Carlo sampling as described in [Equation 16](#S3.E16), we visualise the predictive uncertainty at a point $x_{i}$ which is computed as $\frac{1}{10}\sum_{n=1}^{10}\|\mathbf{S}^{n}_{i}\|$, where $\mathbf{S}^{n}_{i}$ is the posterior covariance matrix of the last layer given the $n$-th sample from the penultimate layer, and $\|\cdot\|$ is the Frobenius norm.

##### Additional results

The MSE comparison is presented in [Figure 9](#A1.F9).
The individual results for each month are shown in [Figure 11](#A1.F11) and [Figure 10](#A1.F10).
A larger version of [Figure 5(b)](#S4.F5.sf2) as well as its analogs for other altitudes are presented in [Figure 13](#A2.F13).
[Figures 14](#A2.F14), [15](#A2.F15) and [16](#A2.F16) present the ground truth vector field, predictive mean and uncertainty, and one posterior sample for July 2010 for the three different altitudes.

Figure: Figure 10: NLPD of residual deep GPs on the wind modelling task across the 12 months of 2010.
Refer to caption: /html/2411.00161/assets/x7.png

Figure: Figure 11: MSE of residual deep GPs on the wind modelling task across the 12 months of 2010.
Refer to caption: /html/2411.00161/assets/x8.png

### A.6 Regression on UCI datasets

Figure: Figure 12: NLPD and MSE of residual deep GPs with spherical harmonic features and Euclidean deep GPs with inducing points on five UCI datasets. Residual deep GPs had their inputs projected from $\mathbb{R}^{n}$ to $\mathbb{S}_{n+1}$ via $x\mapsto(x,b)/\|(x,b)\|$. Solid lines give the mean MSE and shaded regions around them span $\pm 1$ standard deviation. All statistics were computed from 5 randomised runs.
Refer to caption: /html/2411.00161/assets/x9.png

##### Data

In the projection from $\mathbb{R}^{d+1}$ to $\mathbb{S}^{d}$, we set the $b$ to $1$.
This is done for all datasets and the bias is kept constant during training.
Indeed, in our initial examinations, we found that learning the bias often seemed to result in overfitting and worse performance.

##### Training and evaluation

We train each model for 5000 iterations using the Adam optimiser kingma2015 () and the learning rate set to $0.01$.
Each iteration consists of a gradient step using a batch of data.
When the size of the training set is smaller than 1000 data points—that is, for the Yacht, Concrete, and Energy datasets—a batch is the entire dataset.
For the Kin8mn and Power datasets, whose training sets are considerably larger, a batch of 1000 data points is sampled with replacement from the training set.

##### Additional results

The test NLPD and MSE of all models can be seen in [Figure 12](#A1.F12).

## Appendix B More on Interdomain Inducing Variables on Manifolds

Empirically, we find that when the number of eigenfunctions $K$ used to approximate the manifold Matérn kernels exceeds the number of interdomain inducing variables, performance of residual deep GPs deteriorates.
This can be surprising, since a higher $K$ yields a better approximation of the true manifold Matérn kernel.

A potential reason for this phenomenon may be identified by examining the kernel of a posterior sparse Matérn GP $f\sim\mathcal{GP}(\mu,k)$—more specifically, its whitened reparameterisation (see [Section A.1](#A1.SS1)), as that makes equations cleaner.
In this, it will be helpful to define

$$ $\mathbf{\Psi}_{i:j}(\cdot)=\left(\sqrt{\frac{\sigma^{2}}{C_{\nu,\kappa}}\Phi_{\nu,\kappa}(\lambda_{i})}\phi_{i}(\cdot),\dots,\sqrt{\frac{\sigma^{2}}{C_{\nu,\kappa}}\Phi_{\nu,\kappa}(\lambda_{j})}\phi_{j}(\cdot)\right)^{\top},$ (35) $$

which allows us to express the Matérn kernel approximated with $K+1$ eigenfunctions as $k_{\nu,\kappa,\sigma}(\cdot,\cdot^{\prime})=\mathbf{\Psi}_{0,K}(\cdot)^{\top}\mathbf{\Psi}_{0,K}(\cdot^{\prime})$.
Now, recalling [Equation 19](#S3.E19), and denoting the number of interdomain inducing variables by $M$, we see

$$ $\displaystyle k(\cdot,\boldsymbol{\zeta})\mathbf{L}^{-\top}$ $\displaystyle=\left(\phi_{0}(\cdot),\dots,\phi_{M}(\cdot)\right)^{\top}\mathrm{diag}\left(\frac{1}{\frac{\sigma^{2}}{C_{\nu,\kappa}}\Phi_{\nu,\kappa}(\lambda_{0})},\dots,\frac{1}{\frac{\sigma^{2}}{C_{\nu,\kappa}}\Phi_{\nu,\kappa}(\lambda_{M})}\right)^{-1/2}$ (36) $\displaystyle=\left(\phi_{0}(\cdot),\dots,\phi_{M}(\cdot)\right)^{\top}\mathrm{diag}\left(\sqrt{\frac{\sigma^{2}}{C_{\nu,\kappa}}\Phi_{\nu,\kappa}(\lambda_{0})},\dots,\sqrt{\frac{\sigma^{2}}{C_{\nu,\kappa}}\Phi_{\nu,\kappa}(\lambda_{M})}\right)$ (37) $\displaystyle=\mathbf{\Psi}_{0,M}(\cdot)^{\top},$ (38) $$

which we can substitute into [Equation 27](#A1.E27)

$$ $\displaystyle k_{\boldsymbol{\zeta},\boldsymbol{m}^{\prime},\mathbf{S}^{\prime}}(\cdot,\cdot^{\prime})$ $\displaystyle=k(\cdot,\cdot^{\prime})-k(\cdot,\boldsymbol{\zeta})\mathbf{L}^{-\top}(\mathbf{I}-\mathbf{S^{\prime}})\mathbf{L}^{-1}k(\boldsymbol{\zeta},\cdot^{\prime}).$ (39) $\displaystyle=\mathbf{\Psi}_{0:K}^{\top}(\cdot)\mathbf{\Psi}_{0:K}(\cdot^{\prime})-\mathbf{\Psi}_{0:M}^{\top}(\cdot)(\mathbf{I}-\mathbf{S^{\prime}})\mathbf{\Psi}_{0:M}(\cdot^{\prime})$ (40) $\displaystyle=\mathbf{\Psi}_{M+1:K}^{\top}(\cdot)\mathbf{\Psi}_{M+1:K}(\cdot^{\prime})+\mathbf{\Psi}_{0:M}^{\top}(\cdot)\mathbf{S^{\prime}}\mathbf{\Psi}_{0:M}(\cdot^{\prime})$ (41) $$

For $K=M$, posterior covariance reduces to the second term only, which is determined by the kernel hyperparameters and the variational covariance matrix.
However, for $M>K$, the first term contributes additional variance that can only be reduced by changing the hyperparameters of the prior, like length scale and prior variance, rather than the variational parameters $\boldsymbol{m}^{\prime}$ and $\mathbf{S}^{\prime}$.

With this particular setup, there are two forces at play during optimisation: one which lowers the prior variance to match the posterior variance, and another which modifies $\mathbf{S}^{\prime}$ to approximate the true posterior covariance by introducing dependencies between basis coefficients.
In practice, we observe that this can lead to a difficulty if the $\mathbf{\Psi}_{0:M}^{T}(\cdot)\mathbf{S^{\prime}}\mathbf{\Psi}_{0:M}(\cdot^{\prime})$ is already at the desired variance but $\mathbf{S^{\prime}}$ must still be adjusted to approximate the true covariance, while $\mathbf{\Psi}_{M+1:K}^{T}(\cdot)\mathbf{\Psi}_{M+1:K}(\cdot^{\prime})$ is still too large.
In this case, the first mechanism pushes $\mathbf{\Psi}_{M+1:K}^{T}(\cdot)\mathbf{\Psi}_{M+1:K}(\cdot^{\prime})$ downwards by lowering the prior variance $\sigma^{2}$, which necessarily also reduces $\mathbf{\Psi}_{0:M}^{T}(\cdot)\mathbf{S^{\prime}}\mathbf{\Psi}_{0:M}(\cdot^{\prime})$.
Consequently, $\mathbf{S}^{\prime}$ must increase to compensate for this.
This process results in a tug-of-war between the variational parameters and kernel hyperparameters, which seems to make optimisation difficult, and is thus a possible reason for the drop in performance as $K>M$.

One remedy is to set $K=M$, which is what we actually do in our experiments.
However, since this comes at some cost to the kernel approximation, we propose that future work can consider an extended variational family which, in our preliminary tests, helped mitigate this issue at a minimal cost.
Our extension expands the $\mathbf{S}^{\prime}$ matrix with a parametrised diagonal $\mathbf{D}^{\prime}$ corresponding to the $K-M$ eigenfunctions not previously used in the variational family but used in the kernel, giving

$$ $\displaystyle k_{\boldsymbol{\zeta},\boldsymbol{m}^{\prime},\mathbf{S}^{\prime}}(\cdot,\cdot^{\prime})$ (42) $\displaystyle=$ $\displaystyle\mathbf{\Psi}_{M+1:K}^{T}(\cdot)\mathbf{\Psi}_{M+1:K}(\cdot^{\prime})+\mathbf{\Psi}_{0:M}^{\top}(\cdot)\mathbf{S^{\prime}}\mathbf{\Psi}_{0:M}(\cdot^{\prime})+\mathbf{\Psi}_{M+1:K}^{T}(\cdot)(\mathbf{D^{\prime}}-\mathbf{I})\mathbf{\Psi}_{M+1:K}(\cdot^{\prime})$ (43) $\displaystyle=$ $\displaystyle\mathbf{\Psi}_{0:M}^{\top}(\cdot)\mathbf{S^{\prime}}\mathbf{\Psi}_{0:M}(\cdot^{\prime})+\mathbf{\Psi}_{M+1:K}^{\top}(\cdot)\mathbf{D^{\prime}}\mathbf{\Psi}_{M+1:K}(\cdot^{\prime}).$ (44) $$

This eliminates the aforementioned conflict, allowing the variational parameters to affect both terms in a similar way.

Ostensibly, when $K-M$ is large, it significantly impacts performance.
However, as we have seen in [Section A.1](#A1.SS1), this can be avoided by applying the addition theorem, which is possible
if the parameters of $\mathbf{D}^{\prime}$ corresponding to eigenfunctions with the same eigenvalues are kept equal.
With this method, the number of additional variational parameters is minimal, and, comparing our extended variational family with its original variant, there is practically no increase in computation time, as $\mathbf{\Psi}_{M+1:K}^{T}(\cdot)\mathbf{D^{\prime}}\mathbf{\Psi}_{M+1:K}(\cdot^{\prime})$ is simply substituted for $\mathbf{\Psi}_{M+1:K}^{T}(\cdot)\mathbf{\Psi}_{M+1:K}(\cdot^{\prime})$.

Figure: (a) 5.0 km altitude
Refer to caption: /html/2411.00161/assets/figures/wind/uncertainty_and_difference-level=0.png

Figure: (a) Ground truth.
Refer to caption: /html/2411.00161/assets/figures/wind/ground_truth-level=0.png

Figure: (a) Ground truth.
Refer to caption: /html/2411.00161/assets/figures/wind/ground_truth-level=7.png

Figure: (a) Ground truth.
Refer to caption: /html/2411.00161/assets/figures/wind/ground_truth-level=15.png

## Peer Reviews

confidence: 3
contribution: 3
initial_rating: 6
presentation: 4
questions: 1) Can you please make more clear for readers without expertise in differential geometry how your approach to defining a manifold-valued stochastic process differs from that of Mallasto and Feragen, who also use an exponential map on the tangent space at a given point? Please emphasize the practical implications of the differences.
Also, is it the case that that method could not support a deep GP, or that there would be considerable work required to make a deep GP out of it? 

2) Does your work have implications for spaces that have much more structure than a manifold, like vector or function-valued mappings? Or do we recover already-used methods when deploying it in such context.
soundness: 3
strengths: --- 

Excellent presentation

Between the excellent illustrative figures such as figs 1/2 and the compelling narrative flow, this paper is a pleasure to read, which is especially outstanding given the proclivity of some authors to lean into the formalism when discussing ML on manifolds.

--- 

Clear and consequential applications

The authors make a good case that Deep GPs can be superior to single-layer GPs, and that manifold data are important, such that getting deep GPs for manifolds is an important problem to solve.
summary: The authors propose a technique for manifold->manifold GP regression, which they propose to use as a component for deep GP regression on manifolds. They illustrate their procedure via a battery of numerical experiments on real and synthetic problems, including regression and Bayesian optimization.
weaknesses: Since this paper uses numerical experiments rather than theory to introduce its methodology, the experiments are crucial to get right. I have some comments on them below.

--- 

Need clearer comparison to prior manifold->manifold GPs.

It would be help clarify the novelty and contribution of this article to make a more direct comparison to prior work with GPs on manifolds. Probably it is obvious to the authors, but it would be helpful to say exactly what the shortcomings of Mallasto and Feragen's 2018 work, say, which this work addresses. It would also be interesting to compare this methodology to theirs, or to explain why this is not feasible. See Questions for more. 

--- 

Regarding BO on Ackley.
"This outcome aligns with our expectations, as the region around the minimum of the Ackley function is smooth and can be effectively modelled by a shallow GP."
But the model is modeling the entire function, not just a region around the optimum, right? And the ackley function exhibits significant nonstationarity, so we would expect a deep GP to do a better job modeling it. How do the two methods compare in terms of out-of-sample MSE on this function? Does the deepGP do a better job modeling it, but just not optimizing it?

---

 Section 4.4 

I'd like to first commend the authors for sharing a negative result; the failure to investigate limitations of methods is doubtless a serious drag on machine learning research.
However, I have methodological concerns about this study, such that I fail to be convinced that a negative result is indeed appropriate. 
i) The per-iteration complexity is interesting, but it's not the ultimate issue of concern, which is wall-time to a given accuracy (as iterations required may differ between methods).
ii) The use of different optimization methods makes the comparison difficult; I understand the idea may have been to make a comparison to the existing baseline rather than a different one, but it's crucial to identify where any advantage or lake thereof lies by using the Euclidean GP with Adam as well. 
iii) It would be great if you could clarify on where the memory requirements for L-BFGS are coming from; a priori, one would expect that, since Adam keeps track of two vectors of the same size as the parameter vector, using L-BFGS with a history of size 2 would have similar memory requirements, and it seems like it would be possible to use a slightly larger one, like size 10.

---

confidence: 2
contribution: 3
initial_rating: 8
presentation: 4
questions: See weaknesses.
soundness: 4
strengths: 1. This paper proposes a simple and effective way to construct manifold-to-manifold deep GPs by modeling the residuals, which is both novel and interesting.

2. The paper is well-written with a good logic flow. It provides a clear introduction of the relevant technical background and a discussion of related works, which makes it easy for readers to understand the context of this work. The description and derivation of the model is also clear and easy to follow. The beautiful figures are a nice addition to the excellent content of the paper.

3. The proposed model is thoroughly evaluated on a wide range of synthetic and real-world experiments, including a careful ablation study of different choices of Gaussian vector fields and variational families. Empirical results show that the model consistently provides improved predictive accuracy and robust uncertainty estimation compared to shallow manifold GP baselines, which demonstrates its usefulness and wide applicability.
summary: This paper proposes a manifold-to-manifold deep GP inspired by residual neural networks. It only requires modifying the last layer with a Gaussian vector field that represents the residual from the identity map. The input to that layer is then translated by the residual using an exponential map. The authors demonstrate that this model is effective on a wide range of synthetic and real-world experiments.

Disclaimer: While I have expertise in GPs, I don’t know much about non-Euclidean spaces or manifolds except for a high level understanding. Therefore, I might underestimate or overestimate the technical contributions of this paper. This is reflected in my low confidence score.
weaknesses: It seems that the deep GPs saturate quite quickly as depth increases in practice as shown in Figures 3 and 6, especially when the number of data points is small (which is the setting where GPs are typically employed). More generally, deep GPs can be difficult to train and not very scalable in practice, due to, e.g., mode collapse in variational inference. I guess residual deep GPs also share these limitations or perhaps even worse due to the additional components for handling manifolds. In the conclusion section, could the authors comment more on the limitations of residual deep GP (e.g., balance between dataset size, effective depth, scalability, accuracy, uncertainty quality, robustness, approximation, etc) and provide some suggestions as to when residual deep GPs should be preferred over, say, manifold-to-manifold neural networks (and potentially their ensembles if uncertainty estimation is required)?

---

confidence: 3
contribution: 2
initial_rating: 6
presentation: 3
questions: - While the work is theoretically sound I would love to see more analysis of why it works well. My conclusions from the paper are two, the first, if my data lies on a manifold then of course I should constrain my prior to reflect this. However, it does also seem to be the case that this approach performs better in general. I would love to try and get an intuition into this. We know after 10 years of working with composite GPs that simply they do not work, there are too many symmetries in the structure and together with the lax approximations this makes them very hard to use. What of this does this approach address? I would love to get some insights into this from the authors.
- In general it would be great to have more explanation of the model set-up for the different experiment, what are the thoughts behind specific model choices etc. As far as I understand from the experiment you perform modelling on the hyper-sphere, why is this a motivated choice?
- I would love to see an experiment where you actually have known manifold structures that are combined in a composition. What would be interesting to see is if you can correctly factorise the uncertainty across the different layers in the composition.
- It would be great if the authors could make more comments related to paper [1] which is at the moment only referenced as a footnote comment. 
- Would be great if you could expand the caption to Figure 2 to include explanations to both a and b. Right now there is a lot of jumping between pages to get that caption.
  
[1] Mallasto, A., & Feragen, A. (2018). Wrapped gaussian process regression on riemannian manifolds. In , 2018 IEEE/CVF Conference on Computer Vision and Pattern Recognition (pp. 5580–5588). : .
soundness: 3
strengths: This is a very well written paper, there are quite a few concepts to explain to get the background and the authors do this in a splendid manner. A lot of these topics are quite involved and it would have been easy to stray of the necessary path but the authors do an excellent job of not doing this. Very concrete and to the point, well done!

The derivation of the model is well done keeping a good balance between details and intuition of the approach.
summary: This paper outlines the idea of using compositions of Gaussian process priors on manifold structures. The basic idea is to take a normal Gaussian process with identity mean function, here referred to as a residual GP, and at each layer use the exponential map to constrain this to the tangent space of the manifold at x. This Riemanninan formulation has quite a nice connection to the identity mean function as the map itself is now defined by the location as in Eq 10. The paper then continues to fit the proposed modelling approach in how it fits into the literature of the doubly stochastic inference mechanism typically used for composite GPs and with interdomain inducing points to formulate the variational bound.

The paper concludes with a set of experiments, both synthetic experiments and a real experiment on a manifold.
weaknesses: A substantive assessment of the weaknesses of the paper. Focus on constructive and actionable insights on how the work could improve towards its stated goals. Be specific, avoid generic remarks. For example, if you believe the contribution lacks novelty, provide references and an explanation as evidence; if you believe experiments are insufficient, explain why and exactly what is missing, etc.

If I should play devils advocate with this paper the authors do not do enough to highlight the novelty of their approach neither highlight the benefits except for experimentally. We know that composite GPs are horrific to train, while the doubly stochastic bound often allows us to train the model it often leads to very questionable uncertainties. The authors do not do enough to highlight how, or if, this approach addresses these issues. These might come as very generic comments as it more relates to composite GPs in general but this is also how the authors mainly "sell" the paper, that it is a benefit in general to have a manifold structure even if I do not know the manifold of the actual data. It would be great to hear the authors address this.

I also find the Bayes Opt experiment a little bit disappointing, the idea of switching model after a specific set of iterations makes sense but how do I know where to switch, is this the optimal choice that is shown and if so what happens in other locations, both earlier and later? It would be good to see a more in-depth comparison here, what is the result if I use the composite GP from iteration 1? What happens if I swap to another shallow model at the same time? How are the hyper-parameters dealt with, do you marginalize them at each iteration? What is also worrying is the variance in the composite GP doe seem to go up quite substantially even while having 175 points to start from. It would be great to see a few different functions not just Ackley to see if the behavior is consistent across experiments.

---

confidence: 4
contribution: 3
initial_rating: 6
presentation: 4
questions: See weaknesses
soundness: 3
strengths: 1. I think this work is quite a natural and well-motivated extension of recent works on GPs on Manifolds. I find the perspective of defining hidden layer outputs as a residual(deviation from identity maps) to be quite ingenious since it leads to natural generalization to manifolds. As far as I know, this is a novel contribution.

2. I think the paper is quite well-written and very clear in its presentation.
summary: This work generalizes Deep Gaussian Processes to Manifolds. They address the key challenge of modeling each hidden layer being a manifold to manifold map by interpreting the output of each hidden layer as a Gaussian Vector Field combined with an exponential map, allowing each layer to learn a residual displacement on the manifold. From pre-existing works, they outline three different ways for constructing Gaussian vector fields (GVFs) on manifolds: Projected GVFs, Coordinate-frame-based GVFs and Hodge GVFs. They perform inference with this model by optimizing a variational objective.
weaknesses: My main concern is that I found the experimental results to be a bit mixed, some of the choices to be a bit confusing and the main claim of the paper(in 015 "Our models significantly improve performance in these settings ...")  to be not fully validated:
1. The choice of GVF across the experiments seems a bit ad-hoc. It's understandable that the Hodge GVF can only be used when the eigenfunctions of the Laplacian are known but 4.2 uses coordinate-based GVF and 4.4 uses projected gvf. This choice seems arbitrary and it would have been nice to see both the gvfs across all experiments.

2. Euclidean GPs are not considered as baselines in experiments 4.2 and 4.3. Additionally, in experiment 4.1 I am a bit concerned about the fairness of the comparison. It is stated that "we match the number of optimized parameters between models as closely as possible". I am not sure what this means and it would have been nice to make this concrete by having plots of wall-clock times for optimization of the different GPs or something similar. In experiment 4.4 euclidean GPs outperform the manifold variants by a margin. Additionally in experiment 4.1 while Figure 3 paints a flattering picture of the method proposed by the authors, figure 8 seems to paint quite a different picture. The baseline either matches (or outperforms?, It is hard to tell just from the image) the Hodge GVF in terms of MSE. Also, the trend looks like with additional layers deep Euclidean GPs would outperform all the manifold GP variants. As it stands I am not convinced that Residual Deep Gaussian Processes consistently outperform Euclidean Deep Gaussian Processes.

---

confidence: 4
contribution: 4
initial_rating: 8
presentation: 3
questions: See the weakness
soundness: 4
strengths: * The paper is well-written and easy to follow

* The theoretical development is novel and solid

* The experiments are well-designed and thorough, starting from synthetic manifold data, progressing to real-world climate examples, with additional analyses of Euclidean data

* The appendix provides informative supplementary experiments that enhance understanding of the details
summary: The authors proposed practical deep Gaussian process models on Riemannian manifolds by leveraging Gaussian vector fields and manifold-specific operations. Specifically, the proposed interdomain variables-based residual deep Gaussian processes show considerable improvement over the baseline for manifold data. Additionally, this approach can significantly accelerate inference (~2x) for Euclidean data while maintaining comparable performance. The claims made by the authors are well-supported by both rigorous theoretical development and extensive experimental validation.
weaknesses: * The paper is technically dense and may be challenging for readers unfamiliar with Riemannian geometry. It would be helpful if the authors could provide a more detailed background on Gaussian processes in Riemannian geometry in the appendix. For practical use, I also encourage the authors to consider releasing the code for public use in the future

Meta-review: The paper 'Residual Deep Gaussian Processes on Manifolds' was reviewed by 5 reviewers who gave it an average score of 7.6 (final scores: 6+8+8+8+8). The reviewers appreciated the presentation, even if they also pointed out that it is rather dense and technical. The experiments were found a bit mixed, but still interesting. All reviewers argue for accepting the paper.

As recommended by reviewer VHwa, it would be good to include the requested discussion in the camera-ready paper.

Decision: Accept (Oral)