Residual Deep Gaussian Processes on Manifolds
Authors: Kacper Wyrwal
Andreas Krause
Viacheslav Borovitskiy
Venue: ICLR 2025 Oral (2025)
Keywords: Gaussian processes
manifolds
deep Gaussian processes
probabilistic methods
variational inference
uncertainty quantification
geometric learning
Primary area: probabilistic methods (Bayesian methods, variational inference, sampling, UQ, etc.)

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.

Published as a conference paper at ICLR 2025
RESIDUAL DEEP GAUSSIAN PROCESSES ON MANIFOLDS
Kacper Wyrwal
Andreas Krause ETH Zürich
Viacheslav Borovitskiy ETH Zürich
ETH Zürich
University of Edinburgh
ABSTRACT
We propose practical deep Gaussian process models on Riemannian manifolds, similar in spirit to residual neural networks. With manifold-to-manitf hidden layers and an arbitrary last layer, they can model manifold- and scalar-valued fune- tions, as well as vector fields. We target data inherently supported on manifolds, which too complex for shallow Gaussian processes thereon. For example, while the latter perform well on high-altitude wind data, they struggle with the more in- tricate, nonstationary patterns low altitudes. Our models significantly improve performance in these settings, enhancing prediction quality and uncertainty cali- bration, and remain robust to overfitting, reverting to shallow models when addi- tional complexity is unneeded. We further showcase our models Bayesian op- timisation 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 proxy manifold well enough.
INTRODUCTION
Gaussian processes (GPs) are a widely adopted model class for learning functions within the Bayesian framework (Rasmussen and Williams, 2006). 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 (Snoek et al., 2012), active (Krause et al., 2008) and reinforcement (Kamthe and Deisenroth, 2018) learning.
In recent years, substantial work went into developing the analogues practical GP models on var- ious non-Fuclidean domains (Borovitskiy et al., 2021; 2023; 2020; Fichera et al., 2023). By virtue of being geometry-aware these analogues have demonstrated improved performance in variety of tasks on non- non-Fuclidean spaces. Their notable applications include Bayesian optimisation mani- folds for robotics (Jaquier et al., 2022), traffic flow interpolation on road networks (Borovitskiy et al., 2021), and wind velocity prediction on the globe (Hutchinson et al., 2021; Robert-Nicoud al., 2024). These models were also used to speed up inference for Euclidean GPs by transferring data to hypersphere and leveraging the attractive structure GPs thereon possess (Dutordoir et al., 2020).
Despite their advantages, GPs can sometimes fall short in modelling complex, irregular functions. To address this, deep Gaussian processes have been introduced as sequential composition of GPs (Damianou and Lawrence, 2013), providing improved flexibility through their layered structure (Dai al., 2016; Mattos al., 2016). Many techniques developed for shallow GPs, such as variational inference (Salimbeni and Deisenroth, 2017) and efficient sampling techniques (Wilson et al., 2020), 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 na- ture. often competing with Bayesian neural networks and deep ensembles. However, there has been limited work on expressive uncertainty-quant quanttffiiiiiiiinggggg 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, deep GP on 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 handling manifold-input GPs, outputs on manifolds conflict with the fundamental
Correspondence to wyrwal. cacceergmaa 1 com and v viaches lav. waw.ooovvvvoooooiooooooooooooooe t y@gma 1 com. Code available at https: / / thub. com kacprruy / resi dual deep-gge gps.

---

Published as a conference paper at ICLR 2025
Hidden Layer L
GVF
exp
GP
Schematic illustration of scalar- residual deep GP with L hidden layers. The last layer is scalar-valued GP on the manifold. If it is not present, the model is manifold-vale valued. If replaced with Gaussian vector field (GVF), the model a vector field on the manifold.
Figure 1:
concept GP, which dictates that outputs must be Gaussian and thus Euclidean. 1 Designing GPs with inputs or outputs on manifold is challenging; with both, is even more
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-te -o-mmdofooo GP layer directly, we represent it as a Gaussian vector field (GVF) combined with an exponential map. The former represents displacement vector, deviation from the identity map, translates the input by the given displacement vector. The mean of layer is always the output of the previous layer. We visualise this architecture in Figure 1, with sphere the manifold of interest. Notably, by changing the last layer only, one can get manifold-valuee or rector-valued deep GPs. As Section 3 will show, our residual deep Gaussian processes generalise the deep GP architecture of Salimbeni and Deisenroth (2017), perhaps the most successful architecture in the Euclidean case.
residual, while the latter
We examine residual deep GPs through synthetic and real-world experiments, demonstrating our models' superior performance over shallow geometry-aware GPss on tasks where complex data in- herently lies on 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 S,, due to their importance in key applications such as climate modelling and robotics, well as their particularly simple structure that allows for more powerful and spe- cialised GVFs (Robert-Nicoud et al., 2024). However, the applicability of our model extends to all Riemannian manifolds, including ones represented by meshes, indicating 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 there always exist a mean ja: X - I on the input domain X of f and kernel k: X x X + IR such that f() (((xx,, (xx, x)) for all S X. In this case, we write f GP k).
GPs define useful priors for learning functions from noisy observations y E R. at given input locations æ X within the Bayesian framework. fact, if the observation likelihood assumed Gaussian y f'xx) (m)())) E), then the posterior is a GP (Rasmussen and Williams, 2006) with
Prin (-) v(-) + k(, (k)()
Σ) (y m(x:)),
k = k( - kC, (k(2, + Σ) kxx ').
(2)
The input domain X can, principle, be any set. However, good kernel k: X X + R is nec- essary define practical GP models X. If X = R.d, the most widely used kernels from the Matérn family (Rasmussen and Williams, 2006), including. as a limiting case, the especially popular squared exponential kernel.2 These kernels attractive because they implement tive biases: (1) model behaviour should not change under any symmetries of the data-translationss in the case of Radand out that the Matérn family can be generalised to various non-Euclidean domains X in such way that it still implements (1) and (2). We now discuss such generalisation Riemannian manifolds
natural induc-
and (2) the unknown function possesses certain degree of smoothness. It turns
Nevertheless. some practical heuristics for this exist in the literature (Mallasto and Feragen, 2018).
"Thi kernel has many names; it also known as the RBF, Gaussian, heat, or diffusion kernel.
For details, see Azangulov et al. (2024), Borovitskiy al. (2020), and Hutchinson al. (2021).
2

---

Published as a conference paper at ICLR 2025
2.1 GAUSSIAN PROCESSES ON RIEMANNIAN MANIFOLDS
A principled way of generalising the family of Matérn kemels to Riemannian manifolds was pro- posed by Lindgren al. (2011) based the ideas dating back to Whittle (1963). Borovitskiy al. (2020) showed that the resulting kernels can be represented as the following infinite series
Ku,,, (x, Cv.r 1=0
o²
Σ t2 un (2)))))))))))))))
DC
PUEN (2)N(A)
2v
V 10000
(3)
==
where- - 0j are the eigempairs of the Laplace- Beltrami operator on X, is the dimension of X, and Ce,s a normalisation constant ensuring wix fx Kiu,, (2,,)) x)) dr a². The infinite sum must be truncated for computational tractability; nevertheless, the rapid decay of the coefficients 2 (Ay)
makes this a sensible approximation with convergence guarantees (Rosa et al., 2023).
Such models were used, e.g., in robotics (Jaquier al., 2022; 2024) and medical (Coveney al., 2020) applications. Their vector field counterparts, to which we return later in the paper, were used for modelling wind velocities on the globe (Hutchinson et al., 2021; Robert-Nicoud et al., 2024).
These models tend to perform well, especially when data is scarce and uncertainty quantification is crucial, but they can struggle capture complex irregular patterns. One potential way 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 composition F = f .o of multiple shallow GPs s (Damianou and Lawrence, 2013). To allow for richer structure, the hidden layers I', for 1 < L 1 1, are typically vector-valued GPs fl: E2 Fd i.c. vectors of scalar-valued GPs stacked together and potentially correlated with each other (Álvarez et al., 2012). The resulting random function F is itself not a GP. Thus, even for Gaussian likelihoods p(y F F())), inference for F is intractable.
To overcome this, various approximate inference techniques for deep GPs were proposed, perhaps the most popular being doubly stochastic variational inference (Salimbeni and Deisenroth, 2017). In it, the intractable posterior F' y approximated in terms of the KL divergence metric DKL by the elements of a certain variational family of tractable distributions. This family itself consists of deep GPs whose layers are sparse GPs (Hensman al., 2013; Titsias, 2009) which we now discuss.
Sparse GPs were originally proposed as variational family for approximating shallow GPs. For them, approximate inference helps scale to big datasets or accommodate non- non-Gaussian likelihoods. Take some GP(u, k). A sparse GP f=.,,5 is family of GPs parameterised by set of m inducing locations z S X, as well as a mean vector m E R.m and a covariance matrix S Rmxm which determine the corresponding inducing variable distribution q(u) N(m, Specifically,
p(f) = Eu-s(u) p(f(-) u, u,
q(u) (m,,)),
(4)
where p(f(.) u. =) is the prior f conditioned on f(z) u. Intuitively, are the pseudo-inputs, are random pseudo-ol soauudo-bservationn. and p(f2, m.)) is kind of pseudo-posterior. is a GP with (x) mx-) k(. 2) 2) (m m(2)),
2) '). These can be readily generalised the vector-valued setting, with m Rmd and S E Rmdxmmd
s.. = k(-,
(6)
For a deep GP, the respective variational family is the distribution of a composition of sparse GPs:
Fe fatm mm .... fz1 3m1mm St
θ = {2' m .'
(7)
The variational parameters O are found by minimising (p(( y) P(Fo)), which is equivalent (Salimbeni and Deisenroth, 2017) to maximising the following evidence lower bound (ELBO)
ELBO- E Erte) log logp(y)
-'_ DKLluuu p(u'),,
"
where q(u') ~ (mmmm, exactly using the formula for DKL between two Gaussian vectors. The first term is intractable, but can be efficiently approximated by drawing sample from p(Fe) this can be done in layerwise andisst-nddddd and subsampling the sum over a mini-batch of inputs æ. This way, optimisation proceeds by stochastic gradient descent.
S ( and p(u') (z')(=')) (2',, , z²)). The second term can be computed
fashion, as Salimbeni and Deisenroth (2017) suggest-
3

---

Published as a conference paper at ICLR 2025
3 RESIDUAL DEEP GAUSSIAN PROCESSES ON MANIFOLDS
In this section, we introduce the new model class of residual deep Gaussian processes on manifolds. generalises the notion of deep GPs to Riemannian manifolds, allowing for the modelling of scalar- and vector- 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 finding a practical notion of conifold-to-mmnine GPs f to serve as its hidden layers:
== lt fl ..... f
f': X for 1 1-1
where, to simplify exposition, we assume that the last layer f is real-valued, although 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 x X 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 Salimbeni and Deisenroth (2017). There, each layer f"; Rd + R2 is of the form
(((x) g'((),,
(10)
where g is zero-mean GP. That is, each layer f displaces its input x by residual vector (x) (x) the difference to the identity transform, much like a residual connection in neural networks (He et al., 2016), which modelled by GP. On a manifold X A Rd, when x E X and g (xr) E Rd, the addition operation in (x) = x + of (Ir) is undefined. However, there is a natural generalisation:
('(x) expx (z)).
(11)
Here, exp;: T,X -> X is the exponential map, the canonical mapping of the tangent space TxX x E X-i.e., the linear space of vectors tangent to X at point b--back to X itself. That is, a point is still displaced by the vector gz(x), but in geometrically sound manner.
The beauty of Equation (11) is that it reduces modelling f' to modelling g'. The latter is vector- valued, and thus compatible with the traditional notion of Gaussian, making the problem conceptu- ally much simpler. Still, major technical difficulty remains: for different inputs x. the value g' (x) must lie in the different spaces T.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 manifold X is a function that takes each x E X to an element of the tangent space T.X of X at the point æ. If X is* submanifold of 3 the difference between a vector field and a general vector-valued function on X is very intuitive: the latter, a vector attached to point x E X can point in any direction, while in the former it must always be tangential to the manifold X at z. This difference can be seen in Figure 2a, which features a vector-valued function on the left and an actual vector field on the right.
like the 2-sphere is a submanifold
A Gaussian vector field (GVF) can be thus thought of as vector- vector-valued GP whose outputs always happen to be tangential vectors. This notion is rigorously formalised in the appendices of Hutchinson et al. (2021). 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 Hutchinson et al. (2021) propose a simple idea, to build a GVF g from any given vector-valued GP h: X c IR RD by projecting its outputs onto the appropriate tangent spaces. Such projection P(.) DR + X exists because, if X a submanifold of R2, then any tangent space T X can be identified with a linear subspace of PR Thus, g(x) P.h((z) defines a random vector field (see Figure 2a), which turns out to be Gaussian because of the linearity of Po)-
Nash embedding theorem proves this always the case for large enough ambient space dimension De N.
4

---

Published as a conference paper at ICLR 2025
Pe
(a) Projected GVF
(b) Coordinate frame GVF
(c) Hodge GVF
Figure 2: Gaussian vector field constructions on the sphere. In (b), orange vectors depict the frame.
Coordinate-frame--ass GVFs Given any coordinate frame {ei 1 that is, set of functions such that {e,(-))(= is a linear basis of T and vector-valued GP h: X - R with compo- nents h;: + I, the equation = (x()));(z) defines GVE. This is shown in Figure 2b.
Hodge GVFs Most recently, Robert-Nicoud et al. (2024) extended the generalisation of Matérn GPs from Section 2.1 to the setting of vector fields on compact manifolds. They derive an analogue of Eqqation (3), representing the respective Hodge Matérn kernels, infinite series
k (x, CRR.. f=0
Dun (A))8)) (z()))
(12)
Here, 8y are the eigenfields of the Hodge Laplacian on X that correspond to the eigenvalues - -j-- is the tensor product, and & 2x is exactly in Equation (3). This family of kernels be made more expressive by using different hyperparameters o², and for different rypes of eigen- fields sy: the pure-divergence 8j. the pure-curl 8,. and the harmonic sy. The result is called Hodge- compositional Matérn kernels (Robert-Nicoud represent the inductive bias of divergence-frn vector fields, i.e. vector fields having no "sinks" and "sources", like the wind velocity field at certain altitudes. This can be
al., 2024; Yang et al., 2023). For example, it can
in Figure 2c.
The first two constructions are universal. This means that by choosing an appropriate h, one can obtain any possible GVF, Although this might seem advantageous, this also a major curse, is often unclear which particular h to take to get good inductive biases. What is more, simple solutions, such h with IID components, may lead undesirable artefacts (Robert-Nicoud al., 2024). On the other hand, the third construction is canonical, in the same way as the Matérn family is canonical 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, Robert-Nicoud et (2024) only provide practical expressions for the eigenfields 5n when X is the circle S,, the 2-sphere S2, any finite product of those. Thus, for the manifolds that go beyond this simple form, other GVF constructions have to 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, mak- ing exact Bayesian inference with them impossible. However, the doubly stochastic variational inference approach described in Section 2.2 is applicable to them after a few adjustments we detail below. What more, for compact manifolds, the approach be further modified to certain interdomain inducing variables, which, as Section 4 shows, tends to offer superior performance.
Doubly stochastic variational inference Consider the analogue of the variational family p(Fe) in Equation (7) with the sparse GP layers fz' s replaced by
fa = (x
(13)
where 92' ml S are sparse Gaussian vector fields. Again, intuitively and in practice, treating the manifold X as submanifold of RD, GVFs can be thought of, and worked with, special kind of
'Rohert-Nicoud et al. (2024) notice that sometimes, by analytically summing the terms corresponding to the eigenvalue >,, the series convergence be sped up by many orders of magnitude, improving efficiency.
5

---

Published as a conference paper at ICLR 2025
vector-valued GPs. For 9z', ml si however, 2' (2),..... .nn) with inducing locations 2f X, and m/ E T2,X, im(S') T x x... T X, (14)
m' = mm)
where im denotes the image of linear operator, and the last two constraints ensure that the random pseudo-observations u ~ (Mm, S) are tangent vectors which lie in the appropriate tangent spaces.
To satisfy the aforementioned constraints during optimisation, one can represent
m = Pamm
s = P
Ps = P... mm
(15)
with arbitrary m E RMD, arbitrary positive semi-definite S E RmDDxmD, and P.: RD - T,X denoting the projection of vectors in F. onto the tangent space T.X, as discussed in Section 3.2.
Instead of using such representations, one can fix a (locally) smooth frame and optimise the coeffi- cients of mn' and S represented in this frame. In any case, one needs make sure z always remain on the manifold X during optimisation, which can be done by using specialised líbraries, such as PYMANOPT (Townsend et al., 2016) or GEOOPT (Kochurov et al., 2020). For oiwdimensional man- ifolds, a fixed grid on X or set of cluster centroids can be effective alternative to optimising z/j.
Finally, to approximate ELBO, we need to sample F(x;)) ~ ((Fo(r))) As in Salimbeni and Deisen- roth (2017), this can be done sequentially. Specifically, if Fxx,) denotes the desired sample, then
f(x))= ft, f = f (Is.m's s (f( for 1 S S L,
f = Di.
(16)
and, given f fl-1 each individual 92 m² ss (ft- can be sampled in the usual manner.
Interdomain inducing variables 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 constructed by replacing the inducing locations = by inducing linear functionals (= (Ci,..... ... Sm).
Each &j takes in vector field and outputs real number. These Ç define sparse GVF through
pl.c.m.s(.) = Euvag(u) p(g(-) u, C),
g(u) = =N(m,)) S),
(17)
where p(g(.) u, <) is the prior g conditioned on ((g) = u, where ((9) = (gI()). Sm (g)). For example, linear functionals of the form (())= here, {cr is a coordinate frame- can be used to recover the usual doubly stochastic variational inference considered above.
The mean and covariance of m, s given by Equation (5), with replaced by Ç and k(C, 5') = Cov(S(g), ((g)), see, for example, Lázaro-Gredill and Figueiras Vidall (2009) and van der Wilk et al. (2020).
k( :) = (())(((()))))))))))))
(18)
Now, the kernel of g can be expressed as (z) @ (y²)) where {o,} 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 (-) = (s, (5, 1a/a; aa yield very simple covariance matrices
(19) In particular, k(C,C) () is diagonal, making it trivial invert. Dutordoir et al. (2020) 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 for further practical and thcoretical considerations regarding this variational family.
k), = (.),
k(G,, = 8,.//4-
Posterior mean, variance, and samples Expectations E Fo(x)) and variances Var (x( of the approximate posterior Fe cannot be computed exactly. Instead, they can be estimated by appropri- ate Monte Carlo averages, with Equation (16) providing a way to sample Fz(x.. However, since these estimates ignore the correlation between Fo Fo(x) and Fo(z), ), they are not continuous as func- tions of z. When continuity or differentiability of IE Fo()) and Var Fo(x) are desirable, another method can be used. The key idea in this case is to draw (approximate) samples from Fo(-) 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 Wilson et al. (2020) and Wilson et al. (2021) in a sequential manner, akin to Equation (16). This approach is useful for visualisation, performance metric estimation, and for working with downstream quanti- ties, such as acquisition functions in Bayesian optimisation, for which differentiability is key for efficiently finding their maxima.

---

Published as a conference paper at ICLR 2025
Residual (IL)
Baseline
Hodge (IV)
Residual (IV)
N 100
200
400
N 800
1600
1.5
-20
2.5
3
5
2
4
1
2 # Layers
3
4
5
1
2
3
4
5
2
3
4
5
# Layers
Layers
Layers
# Layers
Figure 3: NLPD of different residual deep GP variants and the baseline model, the regression problem for the synthetic benchmark function visualised in Figure 4a. Different subplots correspond different training set sizes N. The solid lines represent the mean, while the shaded areas represent the + standard deviation region around it. All statistics are computed over 5 randomised runs.
4 EXPERIMENTS
We begin this section by examining how various GVF and variational family choices impact the re- gression performance of residual deep GPs in synthetic experiments, discussed in Section 4.1. Throughout, we compare our models to a baseline with Euclidean hidden layers. Next, in the robotics-inspired experiments of Section 4.2, we demonstrate that residual deep GPs can signifi- cantly enhance Bayesian optimisation on a manifold when the optimised function is irregular. Fol- lowing this, in Section 4.3, we show state-of-the-un predictive and uncertainty calibration perfor- mance of residual deep GPs in wind velocity modelling on the globe, achieving interpretable patterns even low altitudes where data is more complex and irregular. Finally, in Section 4.4, we explore potential avenues for using residual deep GPs to accelerate inference on inherently Euclidean data.
4.1 SYNTIETIC EXAMPLES
Setup Deep GPs bear the promise of outperforming their shallow counterparts in modelling com- plex, irregular functions. To test this, construct benchmark function f on the 2-sphere Sn with multiple singularities, which is visualised in Figure 4a. We take N E (100, 200, 400, 800, 1600) training inputs æ Fibonacci lattice on S2 and put y (=) regress f" from x and y. On this problem, we compare different modifications of the residual deep GPs amongst themselves and 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 times with different random seeds for the observation noise E.
(0, 10 10-1I). Then,
The baseline is deep GP with Fuclidean (rather than manniold-to-mmifoodd) layers. It constructed by composing vector-valued Matérn GP whose signature X - R with the Euclidean deep GP of Salimbeni and Deisenroth (2017) on 3 For the residual deep GPs, 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.
Models
Results The NLPD values presented in Figure 3. The MSE values exhibit the same trends and can be found in Figure 8 in Appendix A. 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 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
Motivation GPs are a widely used model class for Bayesian optimisation, a technique for op- timising expensive-tus tovelllutte black-box functions that leverages uncertainty estimates to balance exploration and exploitation (Shahriari et al., 2016). In robotics, such problems arise, for exam- ple, when control policy needs to be fine-tuned specific real-world environment. This task
7

---

Published as a conference paper at ICLR 2025
Deep
Shallow
log Regret
-1
-1
S2 Irregular
83 Ackley
-2
-2
50
100 150
200
50
100
150
200
Step
Step
(a) Irregular benchmark function.
(b) Bayesian optimisation performance.
The irregular benchmark function, and Bayesian optimisation performance comparison. The target functions for Bayesian optimisation are: the aforementioned benchmark function, modi- fied have single global minimum (S2 Irregular), and the smooth Ackley function on the 3-sphere (Sz Ackley). In (b), the solid lines represent the median regret, while the shaded areas around them span +1 standard deviation. The statistics are computed over randomised runs.
Figure 4:
was shown to benefit from treating the optimisation space as manifold and using geometry- aware Gaussian processes to drive Bayesian optimisation (Jaquier et al., 2022). 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 Sec- tion 4.1, visualised in Figure 4a, modified have only one global minimum. The second is the much more regular Ackley function, projected onto S3, one of the benchmarks in Jaquier et al. (2022). In each Bayesian optimisation run, we perform the first 180 iterations using shallow geometry- GP, followed by 20 iterations using residual deep GP-both employing the expected improvement acquisition function (see, e.g., Frazier (2018)). In this experiment, we showease the coordinate- frame- based GVFs, as described in Appendix A. We do not use deep GPs in the initial iterations because, as intuition suggests and Section 4.1 affirms, deep GPs start outperforming their shallow counterparts only when data becomes more abundant. Although deep GPs show comparable per- formance even for small datasets, training them more computationally demanding, making their 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 GP models.
Results The optimisation performance, measured in terms of the logarithm of regret, is reported in Figure 4b. We find that residual deep GPs significantly improve performance in the Bayesian optimisation of the irregular function. Specifically, switching to 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 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 Jaquier al. (2022), with nearly identical median regret trajectories. This outcome aligns with our expectations, since the region around the minimum, explored during the initial 180 iterations, is smooth and thus modelled equally well by both deep and shallow models.
4.3 WIND INTERPOLATION ON THE GLOBE
Motivation Non-Euclidean geometry has 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 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 (Hutchinson et al., 2021; Robert-Nicoud et al., 2024). 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 interpolating the monthly average wind velocity from the ERA5 dataset (Hersbach et al., 2023), from set of locations the Acolus satellite track (Reitebuch,
8

---

Published as a conference paper at ICLR 2025
(a) The ground truth wind ve- locities as black arrows and the training locations along the Aco- lus satellite track as red points.
(b) Difference between the prediction and the ground truth wind veloci- ties, shown as black arrows, and the predictive uncertainty, shown using colour scale from purple (lowest) to yellow (highest), for 3-layer residual deep GP and wind velocities for July 2010, at 0.1 km altitude.
Figure 5: Using residual deep GPs for probabilistic wind velocity modelling on the surface of Earth.
2012), simulating a practical setting of weather-analysing satellite. We use ERA5 data from Jan- uary to December 2010, as in Robert-Nicoud et al. (2024), and sample the Acolus track locations every minute for 24-hour period from 9:00 am, January 1st, 2019. We choose 24-hour period instead of a 1-hour period, as in Hutchinson et al. (2021), because in that time frame, the satellite produces 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 5a. 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 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; MSE follows similar trends and can found Figure in Appendix A. We find that residual deep GPs improve upon the state-of-the-ar shallow Hodge GVFs (I-layer models in the plots, Robert-Nicoud et al. (2024)) both in prediction quality and uncertainty calibra- tion, evidenced by the significantly lower NLPD and MSE. Furthermore, Figure 5b shows that the un- certainty estimates of deep model at the lowest altitude are interpretable. Indeed, regions of high un-
5.5 km
2.0 km
0.1 km
certainty follow regions of irregular wind currents, such as boundaries where multiple currents meet or continental boundaries, as well as areas of season-
NLPD
ally high winds, such as India during the peak mon- soon season. At the same time, low uncertainty is as- signed to regions with constant-like currents. This is unachievable for shallow GVFs since their posterior covariance depends only the locations of the ob- servations, which are rather uniformly dense in our setup, and not the observations themselves. Ad- ditionally, Figure 5b shows predictive uncertainty to be well-calibrated, with areas of highest error corre- sponding to regions of high uncertainty.
# Layers
Figure 6: NL.D of residual deep GPs on the wind modelling task across three altitude lev- cls. Solid lines give the mean NLPD, while the shaded regions around it span +1 stan- dard deviation. Statistics are computed from 12 runs one for each month of 2010.
4.4 ACCELERATING INFERENCE FOR EUCLIDEAN DATA
Motivation Inspired by their connection infinitely wide neural networks, Dutordoir et al. (2020) 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

---

Published as a conference paper at ICLR 2025
Fuclidean (PI)
Residnal (IV)
Yacht
Energy B=691, d=8
Concrete B=927, d=8
Kin8mn
Power
B 277, d=6
B 10000, d=8
B 1000, d=4
2
5.0
Time (s)
1
2.5
2
0
0.0
0
2
3
2
4
2
3
5
1
3
4
# Layers
Layers
# Layers
# Layers
# Layers
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 UCI datasets, with dimension and batch size B. Solid lines show the mean, computed by averaging over 100 training steps, while the shaded areas span 1 standard deviation. However, they are often too narrow to be visible.
using the shallow analogue of the interdomain inducing variables applied to data mapped from Rd to the proxy manifold of S. can be significantly faster than dotucing-location based variational ference for 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 Dutordoir et al. (2020), and we use the same mapping from IR to the proxy manifold Sa. We use projected GVFs to accommo- date
hyperspheres. Overall, our experimental setup follows Dutordoir et al. (2020), except that working with deep models required us to optimise ELBB directly instead of marginalising out the variational mean and covariance S in Titsias (2009). Additionally, be- cause of high memory requirements of L-BFGS (Nocedal, 1980) arising from the lack of marginali- sation and depth, we switch Adam (Kingma and Ba, 2015). We use single Intel 17-13700H CPU.
Results We compare the variational inference speed, measured by wall-clock time for a single optimisation step, Figure The advantage for shallow I-layer GPs increases significantly with more layers, offering considerable edge in deep models. However, predictive performance com- parisons in Figure 17 in Appendix 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 the mapping from Rd to 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 more efficient alternative to Euclidean deep GPs. Achieving this, however, will require further work, which is beyond the scope of this paper.
optimisation. Also, our choice of
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 dis- cussed two variational inference techniques, including one tailored to the structure of Gaussian vector fields compact manifolds and based interdomain inducing variables. We evaluated our models in synthetic experiments, examining the impact of Gaussian vector field and variational fam- ily choices. These experiments supported favouring Hodge Gaussian vector fields and interdomain inducing variables. They also demonstrated that increasing the number of layers virtually never de- grades our models' performance, though can quickly saturate and plateau. We hypothesize that larger datasets will slow saturation, necessitating more layers' complexity. However, we leave this for future work to explore. In coootics-motivatee stylised experiment, our models significantly en- hanced Bayesian optimisation for irregular function on the sphere. For probabilistic interpolation of wind velocities, we achieved of performance, surpassing the recently proposed shal- low Hodge Gaussian vector fields. Finally, 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 Fuclidean 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.
Mapping d x + (21,-., Fz. b)/ (11, Fd. b)ll Sa, where E is bias term, J.)).
10

---

Published as a conference paper at ICLR 2025
ACKNOWLEDGMENTS
VB was supported by ELSA (European Lighthouse on Secure and Safe AI) funded by the European Union under grant agreement No. 101070617. KW thanks Edoardo Ponti for his mentorship. The Blender rendering scripts we used for plotting were adapted from Terenin (2022).
REFERENCES
M. A. Álvarez, L. Rosasco, N. D. Lawrence, et al. Kernels for Vector- Vector-Valued Functions: A Review. Found. Trends Mach. Learn., 4(3):195-266, 2012. Cited on page 3.
I. Azangulov, A. Smolensky, A. Terenin, and V. Borovitskiy. Stationary Kernels and Gaussian Pro- cesses on Lie Groups and their Homogeneous Spaces I: the compact case, Journal of Machine Learning Research, 25(280)11552 2024. Cited on pages 2, 14.
V. Borovitskiy, I. Azangulov, A. Terenin, P. Mostowsky, M. Deisenroth, and N. Durrande. Matérn Gaussian processes on graphs. International Conference on Artificial Intelligence and Statis- lics, 2021. Cited on page 1.
V. Borovitskiy, M. R. Karimi, V. R. Somnath, and A. Krause. Isotropic Gaussian Processes on Fi- nite Spaces of Graphs. In International Conference on Artificial Intelligence and Statistics, 2023. Cited on page 1.
V. Borovitskiy, A. Terenin, P. Mostowsky, and M. P. Deisenroth. Matérn Gaussian processes Riemannian manifolds. In Advances in Neural Information Processing Systems, 2020. Cited on pages 1-3.
S. Coveney, Corrado, C. H. Roney, D. O'Hare. E. Williams, M. D. O'Neill, S. A. Niederer, R. H. Clayton, J. E. Oakley, and R. D). Wilkinson. Gaussian process manifold interpolation for probabilistic atrial activation maps and uncertain conduction velocity. Philosophical Transactions of the Royal Society A, 378(213))22110833)
2020. Cited on page 3.
Z. Dai, A. C. Damianou, J. I. González, and N. D. Lawrence. Variational Auto-encoded Deep Gaus- sian Processes. In International Conference on Learning Representations, 2016. Cited on page 1.
A. Damianou and N. D. Lawrence. Deep Gaussian Processes. International Conference on Arti- ficial Intelligence and Statistics, 2013. Cited on pages 3.
E. De Vito, N. Mücke, and L. Rosasco. Reproducing kernel Hilbert spaces on manifolds: Sobolev and diffusion spaces. Analysis and Applications, (0(3))363-396) 2021. Cited on page 14.
V. Dutordoir, N. Durrande, and J. Hensman. Sparse Gaussian Processes with Spherical Harmonic Features. International Conference on Machine Learning, 2020. Cited on pages 1, 6, 10, 14.
B. Fichera, V. Borovitskiy, Krause, and A. Billard, Implicit Manifold Gaussian Process Regres- sion. In Advances in Neural Information Processing Systems, 2023. Cited page 1.
P. I. Frazier. A tutorial on Bayesian optimization. arXiv preprint arXiv: 2018. Cited on page 8.
K. He, X. Zhang, S. Ren, and J. Sun. Deep Residual Learning for Image Recognition. In Conference on Computer Vision and Pattern Recognition, 2016. Cited on page 4.
Hensman, N. Fusi, and N. D. Lawrence. Gaussian processes for big data, 2013. Cited on page 3.
H. Hersbach, B. Bell, Berrisford, Biavati, A. Horányi, Muñoz Sabater, Nicolas, C. Peubey, R. Radu, I. Rozum, D. Schepers, A. Simmons, C. Soci, D. Dee, and J.-N. Thépaut. ERA5 monthly averaged data on single levels from 1940 to present, 2023. Accessed on 13-Sep-2024. Cited on page 8.
M. Hutchinson, A. Terenin, V. Borovitskiy, S. Takao, Y. Teh, and M. Deisenroth. Vector-valued Gaussian Processes on Riemannian Manifolds via Gauge Independent Projected Kernels. Ad- vances in Neural Information Processing Systems, 2021. Cited pages 1-4, 8,

---

Published as a conference paper at ICLR 2025
N. Jaquier, V. Borovitskiy, A. Smolensky, A. Terenin, T. Asfour, and L. Rozo. Geometry- aware Bayesian Optimization in Robotics using Riemannian Matérn Kernels. In Conference on Robot Learning. 2022. Cited pages 1, 3, 8, 15-17.
N. Jaquier, L. Rozo, M. González-Duque, V. Borovitskiy, and Asfour. Bringing motion tax- onomies to continuous domains via GPLVM on hyperbolic manifolds, In International Confer- ence on Machine Learning. 2024. Cited on page 3.
Kamthe and Deisenroth. Data-Efficient Reinforcement Learning with Probabilistic Model Predictive Control. In International Conference on Artificial Intelligence and Statistics, 2018.
Cited on page
D. Kingma and J. Ba. Adam: A Method for Stochastic Optimization. In International Conference for Learning Representations, 2015. Cited pages 10, 16, 17, 19.
M. Kochurov, R. Karimov, and S. Kozlukov. Geoopt: Riemannian Optimization in PyTorch. In International Conference on Machine Learning, 2020. Cited on page 6.
A. Krause, A. Singh, and C. Guestrin. Near-Optimal Sensor Placements in Gaussian Processes: Theory, Efficient Algorithms and Empirical Studies. Journal of Machine Learning Research, 9(8):235-22) 284,,,4,, 2008. Cited on page 1.
M. Lázaro-Gredilla and A. R. Figueiras- Vidal. Inter-domain Gaussian Processes for Sparse Infer- ence using Inducing Features. Advances Neural Information Processing Systems, 2009. Cited page 6.
F. Lindgren, H. Rue, and J. Lindström. An explicit link between Gaussian fields and Gaussian Markov random fields: the stochastic partial differential equation approach. Journal of the Royal Statistical Society: Series B (Statistical Methodology), 73(4):423-495 2011. Cited page 3.
A. Mallasto and A. Feragen. Wrapped Gaussian process regression on Riemannian manifolds. In Conference on Computer Vision and Pattern Recognition, 2018. Cited on pages 2, 19, 22.
A. G. d. G. Matthews. Scalable Gaussian process inference using variational methods. In 2016. Cited on page 14.
C. L. C. Mattos, Z. Dai, A. C. Damianou, J. Forth, G. D. A. Barreto, and N. D. Lawrence. Recurrent Gaussian Processes. In International Conference Learning Representations, 2016. Cited on page 1.
J. Nocedal. Updating Quasi-Newton Matrices With Limited Storage. Mathematics of Computation, 35(151):773-78, 1980. Cited on page 10.
C. Rasmussen and C. K. Williams. Gaussian Processes Machine Learning. MIT Press, 2006.
Cited on pages 1, 2.
O. Reitebuch. The Spaceborne Wind Lidar Mission ADM-Acolus. In Atmospheric Physics, pages 815-827. 2012. Cited on page
D. Robert-Nicoud, A. Krause, and V. Borovitskiy. Intrinsic Gaussian Vector Fields on Manifolds. In International Conference Artificial Intelligence and Statistics, 2024. Cited pages 1-3, 5, 8, 9, 14.
P. Rosa, V. Borovitskiy, A. Terenin, and J. Rousseau. Posterior Contraction Rates for Matérn Gaus- sian Processes on Riemannian Manifolds. In Advances in Neural Information Processing Systems, 2023. Cited on page 3.
H. Salimbeni and M. P. Deisenroth. Doubly Stochastic Variational Inference for Deep Gaussian Processes. In Advances in Neural Information Processing Systems, 2017. Cited on pages 1-4, 6, 7, 22.
B. Shahriari, K. Swersky, Z. Wang, R. P. Adams, and N. de Freitas. Taking the Human Out of the Loop: A Review of Bayesian Optimization. Proceedings of the IEEE, 104(1):148-175, 2016. Cited on page?
12

---

Published as a conference paper at ICLR 2025
J. Snoek, H. Larochelle, and R. P. Adams. Practical Bayesian optimization of machine learning algorithms. In Advances in Neural Information Processing Systems, 2012. Cited on page 1.
A. Terenin. Gaussian processes and statistical decision-making in non-Euclidean spaces. arXiv preprint arXiv:2202.10611 2022. Cited page 11.
M. Titsias. Variational Learning of Inducing Variables in Sparse Gaussian Processes, Interna- tional Conference on Artificial Intelligence and Statistics, 2009. Cited on pages 3, 10.
J. Townsend, N. Koep, and S. Weichwald. Pymanopt: A Python Toolbox for Optimization on Man- ifolds using Automatic Differentiation. Journal of Machine Learning Research, 17(137):1-5, 2016. Cited on pages 6, 17.
M. van der Wilk, V. Dutordoir, T. A. Artemev, V. Adam, and Hensman. A Framework for Interdomain and Multioutput Gaussian Processes. arXiv preprint 2XXv:2003.01115, 2020. Cited on page
Whittle. Stochastic processes in several dimensions. Bulletin of the International Statistical Insti- lule, 40(2):974-994, 1963. Cited on page 3.
J. Wilson, V. Borovitskiy, A. Terenin, P. Mostowsky, and M. Deisenroth. Efficiently sampling func- tions from Gaussian process posteriors. In International Conference on Machine Learning, 2020.
Cited on pages 1, 6.
T. Wilson, V. Borovitskiy, A. Terenin, P. Mostowsky, and M. P. Deisenroth. Pathwise Conditioning of Gaussian Processes. Journal of Machine Learning Research, 22(105):1-47, 2021. Cited on page
M. Yang, V. Borovitskiy, and Isufi. ComposSoopoiiiio Edge Gaussian Processes. Interna- tional Conference on Artificial Intelligence and Statistics, 2023. Cited on page 5.
13

---

Published as a conference paper at ICLR 2025
A
ADDITIONAL EXPERIMENTAL DETAILS
A.1 IMPLEMENTATION
Efficient kernel evaluation with the addition theorem In our implementation of manifold Matérn kemels, we utilise the addition theorem for spherical harmonics (De Vito et al., 2021; Dutordoir et al., 2020) to accelerate kernel computation." 7 On Sd. the eigenfunctions of the Laplace- Beltrami operator are known to be certain special functions called spherical harmonics. The addi- tion theorem gives relation between all spherical harmonies corresponding to the (k: + 1) largest eigenvalue of the negative Laplace-Bees Beltrami operator -A, denoted ( and Gegenbauer
polynomials c other special functions:
Σ (1,(())0 & (o j=1
x²)
(20)
with the dot product computed after embedding S in Id+ as the unit sphere centred at the origin, d-1 and C,,d some known absolute constants. Thus, when computing the scalar Matérn kernel 2
on S.a we truncate the infinite sum in Equation (3) to include all spherical harmonics up to the
(K + 1)-th eigenvalue, and apply Equation (20). This gives the following formula
0²
K Dun(Au)en,dCC (x x').
Ku,no (2, (2)
Σ k=0
((o))
(21)
CVN
In case of the Hodge Matérn kernel on $2, we also apply addition theorem, only noting that here Φk.j are replaced with the normalised vector spherical harmonics:
Vdy., for the divergence-only kernel,
VAn
and Vddu. for the curl-only kernel, where is the Hodge star operator (Robert-Nicoud al., 2024).
VAx In this case, Ck kta Thus,
K
D (Pun(AR) ki t
Adly o (x.
(V, & Vor (c.
(22)
Cdiv
dk
A==
heurl k vunno (2,,02(2,2))
2')
Dun (Ak) ks t
(+V t )c (x x').
(23)
Ceurl
V.R
a
The Nodge-compositionaa kernel, which we typically use, is the sum pxx + curl
Accelerated training with whitened inducing variables In all the models we test, approximate inference requires variational mean vector m and a variational covariance matrix S that parame- terise q(u)= = NN(m,S) However, to accelerate convergence during training, instead working work with whitened inducing variables u' L L''u,
with the inducing variables u directly, where is the lower Cholesky factor of k(z, 2) or k((,C) (Matthews, 2016). Thus, in practice, denoting the whitened variational mean and covariance of q(u') m' and S', version of Equation (5):
modified
Hz,m' .s (-) =+(-) +
- m(z))
(24)
(u(-) + k(,,)L - (m' -- -L-))()))),,
(25)
k (. k(-,-) - ()2 2) LS'L ( (22, (1z,))
(26)
-k(--) - -k(, -T S') -k),-)'(z,,)).
(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 x 1. For kernels of output layers we
The analogues of this theorem hold for many other manifolds, see, c.g., Azangulov et al. (2024).
S is parametrised by its lower Cholesky factor ensure positive definiteness during optimisation.
14

---

Published as a conference paper at ICLR 2025
Residual (IL)
Baseline
Hodge (TV)
Residual (IV)
N 100
N 200
N 400
N 800
N 1600
0.0050
MSE
0.0025
0.0000
1
3
5
5
2
Layers
# Layers
# Layers
# Layers
Layers
Figure 8: MSE of different residual deep GP variants and the baseline model, on the regression problem for the synthetic benchmark function in Figure 4a. Different subplots correspond dif- ferent training set sizes N. The solid lines represent the mean, while the shaded areas represent the =o region around it, where the standard deviation, and all statistics were computed over randomised runs.
initialise the variance o² 1.0, while for kernels in hidden layers of an L-layer deep GPs
set o² = T.-T 10- at the start of training. In Section 4.1, Section 4.4, and Section 4.3 we initialise the smoothness parameter to = 2. while in Section 4.2 we set it to v = 5 2 replicate the setup in Jaquier et al. (2022). We optimise the smoothness of the manifold Matérn kernels during training. owing to their differentiability with respect to except in Section 4.2, where we fix to match the setup in Jaquier et al. (2022). Wherever we employ Hodge GVFs we use the Hodge compositional kernel, with separate V, K, a² for the curl-free and divergence-frr 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.
Vector-valued GPs We model vector-valued GPs as a set independent scalar-valued GPs stacked into vector. We utilise this construction in Euclidean deep GPs and in residual deep GPs with projected GVFs.
Inducing locations Following Salimbeni and Deisenroth (2017), for all models utilising the vari- ational family based on inducing locations z', we initialise z 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 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, use samples from the posterior. In evaluation, use 10 samples from the posterior to approximate the MSE and NLPD. In visualisation, we also use samples from the posterior approximate the predictive mean and standard deviation.
A.3 SYNTHETIC EXPERIMENTS
Data To examine the influence of data density model performance in controlled manner, generate the training sets as approximately uniform grid points on the 2-sphere S2. This is done using the Fibonacci lattice method, which, for grid of points, gives the colatitude and longitude of the i-th point as
colatitude = arccos (1
2i +
longitude = 2mi
1 + V5
Φ
(28)
2
Because this method gives consistent coverage of S2 also it generate the test set of 5000
points. We define the benchmark function f" S2 + R by
(w) (Y23 (x) (Y R))x,
(29)
15

---

Published as a conference paper at ICLR 2025
where
105
(0, =
sin (in())))
(30)
32m
(0,, = 15
sin sin
(31)
8T
R(x) (I1 -3, 22)
(32)
(xx) = (atan2 I1), arccos(z3)).
(33)
with x being an element of S2 embedded into R² as the unit sphere centred at the origin. Y2,3, V1.2 are spherical harmonics, smooth functions of their parameters. Singularities in f are caused by the composition with io, which converts x from Cartesian to spherical coordinates, but swaps the positions of the colatitude / and longitude σ. The function Y2.3 o yp has singularities at the poles, while the function Y1, oo 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 way that all models have almost the same total number of optimisable parameters. Specifically, we use the following formulae for the number of variational parameters in hidden layer, given that each GVF has m inducing variables
Chodge
(m² + m2 + m
(34)
covariance
mean
Oeuulideen
Oppojecctd
3 ((m² + m)/2 2 m)
(35)
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 spherical harmonics up to the 7-th negative eigenvalue. For Hodge GVFs, we use 70 interdomain inducing variables, corre- sponding vector spherical harmonics up the 6-th negative eigenvalue of the Hodge Laplacian. Despite this, residual deep GPs with Hodge GVFs are still at a disadvantage: single Hodge GVF has 2555 variational parameters and 2 kernel parameters, single projected GVF and manifold vector- vectue GP has 3822 variational parameters and 3 3 kernel parameters, while one Euclidean vector-valued GP has 3822 variational parameters and 3-4 kernel parameters (we use automatic rel- evance determination (ARD), we get 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, that each of them have additional 3 49 optimisable parameters.
Training and evaluation We optimise all models using the Adam optimiser (Kingma and Ba, 2015) for 1000 iterations with learning rate set 0.01.
Additional results extended comparison between the baseline model and Hodge residual deep GP with up to layers and 10 randomised runs.
The MSE comparison is presented in Figure 8. Figure 9 and Figure 10 show an
A.4 GROMETRY- BAYESIAN OPTIMISATION
Data To obtain an irregular function fr on S2, we modify the target function from Section 4.1 to have only one global minimum near singularity point. Specifically, f" was defined by
(*(x) = (Y2,3 o (x) (I3 + 1) (1 arccos(za)),
(36)
where Y2.3 and ip are as in Appendix A.3. The absence of Yi.2 removes the singularities around the equator, while the added scaling factors create 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, in Figure 4b we present the models that performed best GVFs
in this examination: in the left subplot, a 2-layer residual deep GP using coordinate-fram and inducing locations; in the right subplot, 3-layer residual deep GP using projected GVFs and inducing locations. We set the hyper-prion for the shallow model according to Jaquier et al. (2022). We do not hyper-priors for the deep models in this experiment.
16

---

Published as a conference paper at ICLR 2025
Hodge (TV)
Baseline
N = 400
=
800
N = 1600
-15
NL
2.0
-2.5
1
2 3
6 # Layers
4
5
7
9
10
1
2 3
4 # Layers
5
6
7
8
9
10
2 3 4 # Layers
5
6
7
8
9
10
Figure 9: NLPD of the baseline model and Hodge residual deep GP on the regression problem for the synthetic benchmark function in Figure 4a. Different subplots correspond to different training set sizes N. The solid lines represent the mean, while the shaded areas represent the to region around it, where a is the standard deviation; all statistics are computed over 10 randomised runs.
Hodge (IV)
Baseline
N = 400
N 800
N = 1600
0.0050
MSE
0.0025
0.0000
1
2
3
5
7 9 1
10
3
5 7 9 1
10
3
7 Layers
5
9
10
Layers
Layers
Figure 10: MSE of the baseline model and Hodge residual deep GP the regression problem for the synthetic benchmark function in Figure 4a. Different subplots correspond to different training set sizes N. The solid lines represent the mean, while the shaded areas represent the to region around it, where a is the standard deviation; all statistics are computed over 10 randomised runs.
Optimisation process Replicating the setup of Jaquier et al. (2022), we begin each Bayesian optimisation by sampling initial observations uniformly at random the hypersphere. At each step, we minimise the expected improvement acquisition function using the first-order geometry- aware gradient optimisation implemented in PYMANOPT (Townsend et al., 2016). We approximate the expected improvement acquisition function for deep models using Monte Carlo averages driven by pathwise sampling, described at the very end of Section 3.3. After each optimisation step, reinitialise the model and fit it to the data for 500 iterations, using the Adam optimiser (Kingma and Ba, 2015) with learning rate of 0.01 for the deep model, and the BFGS optimiser for the shallow model.
Additional results and analysis In Figure 4b, we see that the variance of the logarithm of regret increases after switching from the shallow GP the residual deep GP. This increase can be explained by variation in the initial 180 points acquired by the shallow GP between independent runs. The locations of these points determine the quality of fit of the deep GP. Indeed, if the points cluster near the true minimum, the deep model often achieves improvement in 1 or 2 iterations. With fewer points around the optimum the fit is poorer and more iterations are required make an improvement. This variance in the number of steps before a better observation is acquired increases the variance of the regret.
We also examined the performance of residual deep GPs on Bayesian optimisation of the irregular target function without the initial shallow GP stage. We equip the last layer of the deep with the same hyper-priors the shallow GP: however, instead of using the BFGS optimiser, BEGS due to memory constraints. Note that in the original experiment we did not use hyper-prion quasi-newton optimiser for the deep model. We use them here because we found that they are important for effective exploration with shallow GPs and they serve the same purpose for our deep GP.
use L-
17

---

Published as a conference paper at ICLR 2025
5.5 km
2.0 km
km
0.3
MSE
0.2
0.1
2
3
5
Layers
Figure 11: MSE of residual deep GPs the wind modelling task across three altitude levels. Solid lines give the mean MSE, while the shaded regions around it span =11 standard deviation. Statistics computed from 12 runs-one for each month of 2010.
log10 Regret
Deep only
-1
Deep only (outliers)
-2
Shallow Doep
1
50
100
150
200
Step
Figure 12: Comparison of Bayesian optimisation performed with shallow GP followed by resid- ual deep GP vs only with residual deep GP on the irregular target function. Solid lines show me- dian logarithm of regret while the shaded areas extend one standard deviation above and below. Blue dotted lines show three optimisation runs with the deep GP only which did not escape a lo- cal minimum -teeeseeeteeeeeeeeeee runs contribute strongly to the large variance of its regret. Grey dotted line indicates transition from the shallow GP to the deep GP at iteration 180.
We report the logarithm of regret achieved by using the residual deep GP model from the very first iteration in Figure We find that in 12 out of 15 the runs, our model improves upon the shallow GP, often even before the 100th iteration. This was expected, since we have seen that residual deep GPs recover shallow solutions when data is not abundant enough to capture target complexity. We also that the variance of the regret is considerably larger than for the shallow GP. This is largely caused by the outlier runs, indicated with dotted blue line, where the model gets stuck in local minimum. As this experiment is fairly sensitive to the setup setting, this could be due the fact that the baseline an exact GP, while our model recovers a sparse GP, and the fact that our model uses the L-BFGS optimiser, while the exact model uses the BFGS optimiser. Thus, using deep GPs exclusively can be more sample-efficient than employing the initial shallow GP stage; however, with the current setup appears that this poses an increased risk of being stuck in local minimum.
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 the FRA5 dataset. Our test set is grid 5000 points, same as in Appendix A.3. To each location in the test set, also assign the wind velocity from the closest location in the ERAS dataset.
Variational parameters For each GVF within the layers of the tested models we use 198 inter- domain inducing variables. They correspond to all vector spherical harmonics up and including, the 10th negative eigenvalue the Hodge Laplacian. This choice is arbitrary and simply serves to balance quality of fit with training time.
18

---

Published as a conference paper at ICLR 2025
Altitude km
Altitude 2.0 km
Altitude km
0.4
0.3
# Layers
MSB
0.2
2
0.1
4
5
101112
101112
101112
Month
Month
Month
Figure 13: MSE of residual deep GPs on the wind modelling task across the 12 months of 2010.
Altitude km
Altitude 2.0 km
Altitude 0.1 km
# Layers
NLPD
0
2
4
5
1 5 101112
1
6
9 101112
2
7 101112
Month
Month
Month
Figure 14: NLPD of residual deep GPs on the wind modelling task across the 12 months of 2010.
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 described in Equation (16), we visualise the predictive uncertainty at point which is computed as 1 10 16 Σn= ||S |]. where S is the posterior covariance matrix of the last layer given the nith sample from the penultimate layer, and is the Frobenius norm.
Additional results The MSE comparison is presented in Figure 11. The individual results for each month are shown in Figure 14 and Figure 13. A larger version of Figure 5b as well as its analogs for other altitudes are presented in Figure 18. Figures 19 to 21 present the ground truth vector field, predictive mean and uncertainty, and one posterior sample for July 2010 for the three different altitudes.
Additionally, we compare our model with the baseline model from Section 4.1. The final layer of the baseline is coordinate-framd GVF with independent Matérn 5/2 kernels, where the frame given by the gradients of the spherical coordinates (with singularities taken care of), this is motivated in part by the approach of Mallasto and Feragen (2018). The results are shown Figure 15 and Figure 16.
A.6 REGRESSION UCI DATASETS
Data In the mapping from Id to Sg, we set b to 1. This is done for all datasets and the bias is kept constant during training. Indeed, 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 (Kingma and Ba, 2015) with the learning rate set 0.01, such that the learning curves plateau
19

---

Published as a conference paper at ICLR 2025
Baseline
Hodge
5.5 km
2.0 km
0.1 km
NLPD
4 # Layers
5
5
# Layers
# Layers
Figure 15: Comparison of the NLPD of residual deep GPs and the baseline model for wind field modelling across the 0.1 km, 1.0 km, and 5.0 km altitudes.
Baseline
Hodge
5.5 km
2.0 km
0.1 km
0.2
MSE
0.1
5
# Layers
# Layers
# Layers
Figure 16: Comparison of the MSE of residual deep GPs and the baseline model for wind field modelling across the km, 1.0 km, and 5.0 km altitudes.
for both Euclidean deep GPs and residual deep GPs. Each iteration consists of gradient step using batch of data. When the size of the training set 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 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 17.
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 the true manifold Matérn kernel.
A potential reason for this phenomenon may be identified by examining the kernel of a posterior sparse Matémn GP f GP(u, more specifically, its whitened reparameterisation (see Ap- pendix A.1), as that makes equations cleaner. In this, will be helpful to define
(1 02 Cr. Cun
02 Cux Den (Snm()))))))))
Vnz())
V
(37)
which allows us to express the Matém keree approximated with K + eigenfunctions as kunn = (.):(.) (.)( Now, recalling Equation (19), and denoting the number of inter-
20

---

Published as a conference paper at ICLR 2025
Euclidean (IL)
Residual (IV)
Yacht
Energy R=69 691, d=8
Concrete 927, d=8
Kin8mn
Power
B=277 277, d=6
1000, d=8
B
1000, d=4
-2
NLPD
10
10
-3
0
1
2
1
3
5
Layers
# Layers
# Layers
Layers
Layers
Euclidean (IL)
Residual (IV)
Yacht
Energy
Concrete 927, d=8
Kin8mn
Power
B 277, d=6
33=69,15 8=8
1000, d=8
B
1000,
MSE
500
10
50
1000
0
3
5
Layers
# Layers
# Layers
# Layers
# Layers
Figure 17: 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 mapped from Rd to Sa via + (z. b)ll(a,, b)ll- Solid lines give the mean MSE and shaded regions around them span +1 standard deviation. All statistics were computed from 5 randomised runs.
domain inducing variables by M, we see
1/2
1
-(:
= =(>00-) ((o())...,,,M(()) diag
(38)
2 (m))
Puin(Arr
Tox
Co,
(
σ²
(00())) (m)()))) diag
(
Cose
Der (A)::
Cuk
Denu(das)
(39)
(.)
(40)
which we can substitute into Equation (27)
kg.,m f..m/s'// (t,'' kc.,.') k(,, - -
(41)
V (. P. (.:((') - T
(42)
= D TTTLK (-) (+): + L (-)S (.:
(43)
For = 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 m' and S'.
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 S' to approximate the true posterior covariance by introducing dependencies between basis coefficients. In practice, we observe that this lead to difficulty if the V V.. (-)S' T (+") is already at the desired variance but S' must still be adjusted approximate the true covariance, while V MM+KK( ( MM+iKK (+') is still too large. In this case, the first mechanism pushes ( M+IKK( downwards by lowering the prior variance o², which necessarily also reduces P (-)S' Po.M(() Consequently, S' must increase to compensate for this. This process results in oug-of-war between the variational parameters and kernel byperparameters, which seems to make optimisation difficult, and thus a possible reason for the drop in performance K > M.
T
21

---

Published as a conference paper at ICLR 2025
One remedy is to set K M, which is what actually do in 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 minimal cost. Our extension expands the S' matrix with a parametrised diagonal D' corresponding to the K . eigenfunctions previously used in the variational family but used in the kernel, giving
(:.m,s((..)) -V -MM+IK( M M + V V..) v W I MM+
(44)
(45)
(.) * D'V W N+IK(())
(46)
This eliminates the aforementioned conflict, allowing the variational parameters to affect both terms in similar way.
Ostensibly, when K M is large, it significantly impacts performance. However, as we have seen in Appendix A. this can be avoided by applying the addition theorem, which is possible if the parameters of D' 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 & M+1:K( (-)xx(()D' D')D's M (+) is simply substituted for v W (-) M+1x((-)).
C RELATION TO WRAPPED GAUSSIAN PROCESSES
Wrapped GPs of Mallasto and Feragen (2018) are used to model function of Euclidean data taking values on Riemannian manifold. They are constructed by choosing a base point function, which assigns point the manifold to each Euclidean input point, and (although the latter is implicit in the original paper). Inference is done by lifting the training labels from the manifold to the tangent space at the base points assigned to their corresponding inputs, performing inference in the tangent space, and projecting the posterior GVF into the manifold using the exponential map. The base point function is chosen either as constant mapping point minimising the squared distance from the training points (i.e. empirical Fréchet mean) or as an auxiliary regression function.
pardinatte-ramee GVF prior
We derived the manifold-to-manifoll layers of our model as a generalisation of the linear mean construction of Salimbeni and Deisenroth (2017); however, we may also build them based on the ideas of wrapped GPs. The key difference is that our construction is wrapped GPs are Fuclidean-input-aa and the input and output manifolds are identical, allowing lay- ers to be composed. Thus, the first non-trivial modification is to replace the Euclidean domain with the manifold domain. This requires adapting the GVF from Euclidean kernels to manifold kernels. Furthermore, instead of using auxiliary regression function Fréchet mean, the natural choice yielding a generalisation of the linear mean is the identity map for the base point function. This yields manifold-te to -ol GP: however, enable doubly stochastic variational inference, the next modification is to replace the exact GVFs with sparse GVFs using inducing points or interdo- main inducing variables. With these modifications, obtain our cerffooddtt-mmaiffon GPs, which can be composed sequentially to yield residual deep GPs.
-nhut-wheeass
22

---

Published as a conference paper at ICLR 2025
(a) 5.0 km altitude.
(b) 2.0 km altitude.
(c) 0.1 km altitude.
Figure 18: Difference between the prediction and the ground truth wind velocities, shown as black arrows, and the predictive uncertainty, shown using a colour scale from purple (lowest) to yellow (highest), for a 3-layer residual deep GP and wind velocities for July 2010, at three altitude levels.
23

---

Published as a conference paper at ICLR 2025
(a) Ground truth.
(b) Predictive mean and uncertainty.
(c) Posterior sample.
Figure 19: Ground truth wind velocity data altitude of 5.0 km from July 2010, and the corre- sponding posterior mean, uncertainty. and sample from a 3-layer residual deep GP. The mean and sample shown black arrows, while the predictive uncertainty is shown using colour scale from purple (lowest) to yellow (highest).
24

---

Published as a conference paper at ICLR 2025
(a) Ground truth.
(b) Predictive mean and uncertainty.
(c) Posterior sample.
Figure 20: Ground truth wind velocity data altitude of 2.0 km from July 2010, and the corre- sponding posterior mean, uncertainty. and sample from a 3-layer residual deep GP. The mean and sample shown black arrows, while the predictive uncertainty is shown using colour scale from purple (lowest) to yellow (highest).
25

---

Published as a conference paper at ICLR 2025
(a) Ground truth.
(b) Predictive mean and uncertainty.
(c) Posterior sample.
Figure 21: Ground truth wind velocity data altitude of 0.1 km from July 2010, and the corre- sponding posterior mean, uncertainty. and sample from a 3-layer residual deep GP. The mean and sample shown black arrows, while the predictive uncertainty is shown using colour scale from purple (lowest) to yellow (highest).
26


Reviews:
[{"review_id": "mU8mtX4fel", "reviewer": "Reviewer_a17f", "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.", "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.\nAlso, 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? \n\n2) 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.", "limitations": "", "rating": 8, "confidence": 3, "soundness": 3, "presentation": 4, "contribution": 3, "strengths": "--- \n\nExcellent presentation\n\nBetween 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.\n\n--- \n\nClear and consequential applications\n\nThe 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.", "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.\n\n--- \n\nNeed clearer comparison to prior manifold->manifold GPs.\n\nIt 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. \n\n--- \n\nRegarding BO on Ackley.\n\"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.\"\nBut 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?\n\n---\n\n Section 4.4 \n\nI'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.\nHowever, I have methodological concerns about this study, such that I fail to be convinced that a negative result is indeed appropriate. \ni) 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).\nii) 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. \niii) 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.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "boVvnrt8Gr", "reviewer": "Reviewer_VHwa", "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.\n\nDisclaimer: 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.", "questions": "See weaknesses.", "limitations": "", "rating": 8, "confidence": 2, "soundness": 4, "presentation": 4, "contribution": 3, "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.\n\n2. 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.\n\n3. 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.", "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)?", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "y1dWzevhxv", "reviewer": "Reviewer_GjkX", "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.\n\nThe paper concludes with a set of experiments, both synthetic experiments and a real experiment on a manifold.", "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.\n- 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?\n- 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.\n- 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. \n- 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.\n  \n[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). : .", "limitations": "", "rating": 6, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 2, "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!\n\nThe derivation of the model is well done keeping a good balance between details and intuition of the approach.", "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.\n\nIf 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.\n\nI 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.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "DEAP9NzDnk", "reviewer": "Reviewer_ezvQ", "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.", "questions": "See weaknesses", "limitations": "", "rating": 8, "confidence": 4, "soundness": 3, "presentation": 4, "contribution": 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.\n\n2. I think the paper is quite well-written and very clear in its presentation.", "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:\n1. 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.\n\n2. 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.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "Wr3MM3vZAq", "reviewer": "Reviewer_ymXL", "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.", "questions": "See the weakness", "limitations": "", "rating": 8, "confidence": 4, "soundness": 4, "presentation": 3, "contribution": 4, "strengths": "* The paper is well-written and easy to follow\n\n* The theoretical development is novel and solid\n\n* The experiments are well-designed and thorough, starting from synthetic manifold data, progressing to real-world climate examples, with additional analyses of Euclidean data\n\n* The appendix provides informative supplementary experiments that enhance understanding of the details", "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", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}]

Author rebuttal:

Decision: Accept (Oral)