Thin and deep Gaussian processes
Authors: Daniel Augusto de Souza
Alexander V Nikitin
S. T. John
Magnus Ross
Mauricio A Álvarez
Marc Peter Deisenroth
João Paulo Pordeus Gomes
Diego Mesquita
César Lincoln Mattos
Venue: NeurIPS 2023 poster (2023)
Keywords: Gaussian Processes
Deep Gaussian Processes
non-stationary kernels
TL;DR: We propose an extension to deep Gaussian processes that allow for interpretability in terms of lengthscales fields

Abstract: Gaussian processes (GPs) can provide a principled approach to uncertainty quantification with easy-to-interpret kernel hyperparameters, such as the lengthscale, which controls the correlation distance of function values.However, selecting an appropriate kernel can be challenging.

Deep GPs avoid manual kernel engineering by successively parameterizing kernels with GP layers, allowing them to learn low-dimensional embeddings of the inputs that explain the output data.
Following the architecture of deep neural networks, the most common deep GPs warp the input space layer-by-layer but lose all the interpretability of shallow GPs. An alternative construction is to successively parameterize the lengthscale of a kernel, improving the interpretability but ultimately giving away the notion of learning lower-dimensional embeddings. Unfortunately, both methods are susceptible to particular pathologies which may hinder fitting and limit their interpretability.

This work proposes a novel synthesis of both previous approaches: {Thin and Deep GP} (TDGP). Each TDGP layer defines locally linear transformations of the original input data maintaining the concept of latent embeddings while also retaining the interpretation of lengthscales of a kernel. Moreover, unlike the prior solutions, TDGP induces non-pathological manifolds that admit learning lower-dimensional representations.

We show with theoretical and experimental results that i) TDGP is, unlike previous models, tailored to specifically discover lower-dimensional manifolds in the input data, ii) TDGP behaves well when increasing the number of layers, and iii) TDGP performs well in standard benchmark datasets.

Thin and Deep Gaussian Processes
Daniel Augusto de Souza" University College London
Alexander Nikitin
ST John Aalto University
Aalto University
Magnus Ross University of Manchester
Mauricio A. Álvarez University of Manchester
Mare Peter Deisenroth University College London
João P. P. Gomes Federal University of Ceará
Diego Mesquita
César Lincoln C. Mattos Federal University of Ceará
Getulio Vargas Foundation
Abstract
Gaussian processes (GPs) provide principled approach uncertainty quan- tification with easy-to-interpret kernel hyperparameters, such as the lengthscale, which controls the correlation distance of function values. However, selecting appropriate kernel can be challenging. Deep GPs avoid manual kernel engineer- ing by successively parameterizing kernels with GP layers, allowing them learn low-dimensional embeddings of the inputs that explain the output data. Follow- ing the architecture of deep neural networks, the most common deep GPs warp the imput space layer-by-layee but lose all the interpretability of shallow GPs. An alternative construction to successively parameterize the lengthscale of a ker- nel, improving the interpretability but ultimately giving away the notion of learn- ing lower-dimensional embeddings. Unfortunately, both methods are susceptible to particular pathologies which may hinder fitting and limit their interpretability. This work proposes novel synthesis of both previous approaches: Thin and Deep GP (TDGP). Each TDGP layer defines locally linear transformations of the orig- inal input data maintaining the concept of latent embeddings while also retaining the interpretation of lengthscales of kernel. Moreover, unlike the prior methods, TDGP induces non-pathological manifolds that admit learning lower-dimensional representations. We show with thcoretical and experimental results that TDGP is, unlike previous models, tailored to specifically discover dover-dimensional manifolds in the input data, ii) TDGP behaves well when increasing the number of layers, and iii) TDGP performs well in standard benchmark datasets.
Introduction
Gaussian processes (GPs) are probabilistic models whose nonparametric nature and interpretable hyperparameters make them appealing in many applications where uncertainty quantification and data efficiency matter, such as Bayesian optimization [II, spatiotemporal modeling (6), robotics and control [5). The key modeling choice for GP prior is its covariance or kernel function, which determines the class of functions it represents. However, due to finite data and a fixed-form kernel function, the GP may not inter- or extrapolate as desired. Stationary kernels, such as the commonly used squared exponential kernel and the Matérn family, assume the existence of constant charac- teristic lengthscale, which makes them unsuitable for modeling non-stationary data.
To construct more expressive kernels, can consider hierarchical GP models (deep GPs or DGP). The most common deep GP construction is a functional composition of GP layers with standard
Contact: daniel. daniel.souza souza. 210ucl. ac. uk.
37th Conference on Neural Information Processing Systems (NcurIPS 2023).

---

stationary kernels that results in non con-stationary Gaussaan process B 23 for clarity, we will refer this model type as compositional DGPs (CDGPs). However. CDGPs can show pathological behavior, where adding layers leads to a loss of representational ability (9 7l. Alternatively, it is also possible to extend "shallow" GPs by making the kernel lengthscales a function of the input (19), resulting in the deeply non-stationary GP [DNSGP, 22). Although this covariance function approach to DGPs does not degenerate with more layers, care must be taken to guarantee a positive semi-definite kernel matrix. Moreover, the induced space is not proper inner-product space (18). which hinders the learning of useful manifolds.
In this paper, we address the shortcomings of previous DGP paradigms by retaining the flexibility provided by learning hierarchical lengthscale fields while also enabling manifold learning. Instead of pursuing neural networks to enhance standard GPs 26], which may lead to overfitting [17], we re- tain the nonparametric formulation by using additional GPs to model a linear projection of each input onto a latent manifold. The key insight of our approach is that such projections are input-dependent and tailored towards more interpretable lower-dimensional manifolds with corresponding length- scale fields. The resulting Thin and Deep (TDG)) avoids the pitfalls of other DGP constructions while maintaining its hierarchical composition and modeling capacity beyond shallow GPs.
Our contributions are three-fold:
1.
We propose TDGP, new hierarchical architecture for DGPs that is highly interpretable and does not degenerate as the number of layers increases. Notably, TDGP is the only deep architecture that induces both a lengthscale field and data embeddings.
We prove that TDGPs and compositional DGPs are the limits of more general DGP construc- tion. Thus, we establish a new perspective on standard CDGPs while reaping the benefits of inducing a lengthscale field.
3. We demonstrate that TDGPs perform as well as or better than previous approaches. Our experi- ments also show that TDGP leans towards inducing -iwddimnniionll embeddings.
2 Background
Gaussian processes (GPs) are distributions over functions and fully characterized by a mean function m and kernel (covariance function) A (20). If function is GP distributed, we write f ~ GP(m, k). If not stated otherwise, we assume that the prior mean function is 0 everywhere, i.c., m() = 0. Typically, the kernel possesses few interpretable hyperparameters, such lengthscales or signal variances, estimated by the standard maximization of the marginal likelihood (20).
The squared exponential (SE) kernel is arguably the most commonly used covariance function in the GP literature and, in its most general form 255 can be written as
kse a, b) = 0² exp (a b) A a b)].
(1)
where a, b, E R. the constant o² > 0 defines the signal variance, and A E DxD is lengthscale matrix. Notably, the SE kernel is stationary, i.e., kgg + c, b + e) = ksg b) for any e E R.". Furthermore, when the lengthscale matrix A I the kernel is isotoopic, meaning that it can be written as a function usy (d²) of the squared distance d² |aa bl.3. For more details on the significance of the lengthscale parameter, we refer the reader to Appendix A.
Stationary kernels enforce invariances, which may always be desirable. However, stationary kernels can be used as building blocks to derive broader families of kernels (including non-stationary kernels), either by composing them with deformation functions or through mixtures of lengthscales.
Deformation kernels result from applying a deformation function T : RR + Eº to a and before feeding them to stationary kernel k in R Thus, a deformation kernel k, follows
(a, b) = ((T(a),
For a linear transformation T(x) Wx, we can interpret k,p as a stationary kernel with lengthscale matrix A= [WTW] However, for more intricate T(-). interpreting analyzing k can be very challenging. For instance, this is the case for deep kernel learning [DKL, 26) models, in which T (-)) is an arbitrary neural network.
2"Thin" refers to the graph-theoretic girth of the graphical model of proposed DGP construction.

---

Compositional deep GPs JCDGPs, 3) also rely on deformation kernels, but they use a GP prior to model r(:). The kernel of the T-)) process can also (recursively) be considered be a deformation kernel, thereby extending DGPs to arbitrary depths. However, stacking GP layers reduces the interpretability of DGPs, and their non-injective nature makes them susceptible pathologies (9).
Lengthscale mixture kernels [19] are a generalization of the process of convolving stationary ker- nels with different lengthscales from Higdon et al. [14]. For arbitrary isotropic kernels k, i.e. k(a, b) T& (I b)12), Paciorek and Schervish [19] construct non-stationary Almx as a func-
tion of field of lengthscale matrices
R
such that
kmmx (a, b) a ( + A(b)) (6),
(3)
where 5 = 2(a (A(() + A (b)) b). The explicit lengthscale field A makes this model more interpretable than general deformation kernels. However, Paciorek [8] notes that 8(a, b) may violate the triangle inequality and therefore does not induce a manifold over the input space. This departure from the properties of stationary kernels is due the matrix inside the quadratic form being function of both a and b. Another caveat is that the scale of kimx is also controlled by pre-factor term that depends on A, allowing for unintended effects (such unwanted long-range correlations), especially in rapidly varying lengthscale fields.
Deeply non-stationary GPs IDNSGPs, 7 22 use the lengthscale mixture kernel and parameterize the function A (x) with a warped GP prior to obtain deep GP model. Similar to CDGPs, this model can be extended in depth [22 by considering the kernel of A (xx) be non-stationary with its lengthscales stemming from another GP. A practical issue in these models is guaranteeing that A(x) is positive semi-definite. Therefore, in practice, A(x) is usually restricted to be diagonal.
3 Thin and Deep GPs (TDGPs)
As discussed in the previous section, deep GP constructions follow naturally from hierarchical ex- tensions of base kernel. Therefore, we arrange the presentation of TDGPs in three parts. First, we propose kernel that admits interpretations both in terms of its lengthscale and of its induced manifold. Second, we use this kernel to build novel type of deep GP model (TDGP). Third, we describe how to carry out inference for TDGPs. Finally, we discuss limitations of our approach.
Kernel. We address the drawbacks of the approaches based on deformation and lengthscale mix- ture kernels - i.e., the lack of interpretability and failure to induce a manifold, respectively - by proposing a synthesis of both methods, retaining their positives whilst mitigating some of their is- sues. Starting from the discussion in Seectionn 2| a squared exponential kernel with lengthscale matrix [WTW]- corresponds to deformation of an isotropic kernel:
kisT (0, b)
STF ((a b) WT b))
= T (Wa WWbT(WW Wb))
TSE Waa bbll)
(4)
Therefore, this is the same as applying linear deformation r(x) W..
We propose extending this non-linear transformation that locally linear by letting W vary as a function of that input, so (xx) = W W(w)æ This results in the TDGP kernel:
KTDGP (a, b) (T W (a)a ((b)bll
=Ws=((((a)a (W
(5)
Equation (5)
cannot be written as Mahalanobis distance like in Eq. (4) but in the neighborhood of
(x)) is a lengthscale matrix, thus allowing KYDGP to be implicitly parametrized by
a lengthscale field just like the lengthscale mixture kernels Paciorek and Schervish [19]. Hence,
it allows for better interpretability than the deformation kernels considered in compositional DGPs.
However, unlike the lengthscale mixture approach of BB. (3) our kernel in in (5) not introduce an input-dependent fre-factor, thereby avoiding pathologies when the lengthscale varies rapidly. Moreover, since the distance induced by the quadratic form obeys the triangle inequality, it induces
manifold in the input space. Hence, it addresses the two issues of lengthscale mixture kernels.

---

CDCP wii::
TDGP mil::
GP
GP
W
f
W
Figure 1: Graphical models for a CDGP (left) and TDGP (right) with layers. The colored boxes represent the architecture of a single layer; dashed nodes represent where layers connect. CDGP directly composes functions, whereas TDGP builds hierarchy of input-dependent lengthscale fields.
Model. an analogous manner to the compositional DGP and DNSGP, we present how to use this kernel to build hicrarchical GP model. Our L-layers deep model is described as follows:
p()(() h (.)) GP(( 0, TT& (
(a) h--1 (b)
where
h²(x) (n' D Q II GP Wey pen kw) d=1q=1
h hx(x)= I.
(W)())
(8)
Therefore, this model can be extended arbitrarily by deforming the kernel of the entries of W with another locally linear transformation with its own matrix W
compares the graphical models of TDGP and compositional DGP. We notice that TDGP learns hierarchical input-dependent lengthscale fields, instead of a straight composition of processes. We call this model Thin and Deep GP (TDGP) by the fact that our graphical model always has cycles of bounded length due to the connection of every hidden layer with the inputs X. Therefore, it has finite girth, in contrast to the unbounded girth of the compositional DGP graphical model. Importantly, however, TDGP are related to CDGPs as both can be seen as locally affine deformations of the input space (proof in Appenddix B).
Figure I
Theorem 3.1 (Relationship between TDGP and CDGP), Any L-layer CDGP prior over a func- tion f(x) (h (.). (x)...) is special case of generalized TDGP prior with equal depth defined over the augmented imput space ã = [ax, 1]T. Since linear deformations W ã in the augmented space correspond to affine transformations Wx + d in the original space, the special case of the CDGP model corresponds to a TDGP where the prior variance of W approaches zero.
Unlike DNSGP, which makes the positive semi-definite matrix A function and requires the use of warped GP priors, our hidden layers W are arbitrary matrices that admit a regular Gaussian prior. By placing zero-mean priors the entries of W, encourage the MLE estimate to maximally reduce the latent dimensionality of h. This is because the latent dimensionality becomes linked to the number of rows with non-zero variance in the prior, c.g., if the prior kernel variance i-th row of 4 tends to zero, the posterior over that row is highly concentrated at zero, eliminating the i-th latent dimension. In compositional DGP, this inductive bias also corresponds in the hidden layers. However, as discussed by Duvenaud et al. [9) and reproduced in Fig. 2 this choice introduces pathology that makes the derivatives of the process almost zero.
zero-mean prior
Inference. To estimate the posterior distributions and the hyperparameters of our model, we perform variational inference (VT). We introduce inducing points u for the last-layer GP f(() and inducing points V for the W() processes. Similar to how VI methods for DGPs were based on the non- deep GPLVM model [24), our method builds upon VI for the hyperparameters of square exponential kernels in shallow GPs as discussed by Damianou et al. [4. However, we replace the Gaussian prior p( (W) with GP prior p-W(.)) For instance, the variational distribution for a two-layer TDGP is:
Q
qf, W, u, V) = =ff u) N(u jun ESw)
II d=1
pl Wgd Uga) N N(vdd Mogar Evaa), (9)

---

where for the final layer the parameters Pu.. Σn are not estimated but replaced with their closed-form optimal solutions, [Apendix Coontains derivation of the ELBO and more details.
Limitations. For large numbers of layers, our VI scheme for TDGP may become computa- tionally intensive. By stacking additional hidden layers W(x), add D x Q GPs into the model; consequently, the number of variational parameters increases, which can slow down op- timization. More specifically, inference uses OLL x D Q x m²) parameters and takes time (L x D Q x (m² + m²))) to compute the ELBO. Additional details on runtime as a func- tion of data size and width can be found in the Appendix Ej Furthermore, without the addition of bias to the input data æ. TDGP performs locally linear transformations W W(x)aa ar meaning that the neighborhood of =
0 is kept unchanged for each layer.
Layer
Layer
Layer 10
Layer
Layer 50
TDGP DGP
My
Figure 2: Samples from a CDGP prior (top row) and TDGP prior (hottom row), hoth with zero mean, Fach column represents the number of layers. where one layer corresponds to regular shallow GP. When the number of layers increases, CDGP samples quickly become flat, i.e. the prior distribution no longer supports 'interest- ing' functions. Notably, this pathology does not occur with TDGP.
4 Experiments
To assess whether TDGP accurately captures non-smooth behavior and has competitive results with other GP models, we performed regression experiment on synthetic dataset with input-dependent linear projections and another experiment in varied empirical datasets. In this section, we examine the two-layer TDGP against a two-layer CDGP, two-layer DNSGP, DKL, and the shallow sparse GP (SGP). These models are evaluated on average negative log-predietive density (NLPD) and mean relative absolute error (MRAE). All metries are "the lower, the better". Importantly, throughout the experiments, we also emphasize (i) the interpretability of our model compared to the prior art and (ii) TDGP's inductive bias towards learning low-dimensional embeddings. all experiments, inputs and targets are normalized so that the training set has zero mean and unit variance. [Appeddxx D) contains more details of architecture, training, and initialization. We implemented the experiments in Python using GPflow [10), GPflux (8), and Keras (2). Code is available as supplementary material at whtttt:::0 / //ggth/gtttttttt com httmm//eectraaaccmo http:/ppcccrllaaaa////eeeeegggg
4.1 Synthetic experiment
Data. To assess intuition that TDGP leans towards inducing low-dimensional manifolds, show how well TDGP and competitors can fit a composite function f f= g o h, where g : R - R and h: R IR are non-linear functions. In this context, h acts as "funnel" inducing 1D manifold. For more details on f and how we sample from the input space, we refer the reader to Appendix D
Results. Figure 3 shows the posterior mean of TDGP and compet- ing methods for our synthetic dataset. Note TDGP accurately fits the target function (leftmost panel), while the other methods fail to capture the shape of several of the local maxima. Consequently, Table 1 shows that TDGP performs significantly better in terms of both NLPD and MRAE in the test data.
Table Results for the synthetic data. TDGP significantly outper- forms baselines.
NLPD
MRAE
SGP
-1.49
0.11
DKT.
36.82
0.21
Examining the latent space. TDGP outperformed competitors in predictive performance, but how do their latent spaces compare? To further understand the nuances that distinguish each model, Fig. 4
CDGP
-1.28
0.12
DNSGP
-1.46
0.12
TDGP
-3.57
0.00
plots their respective latent spaces and lengthscale fields - along with the true latent space induced by h as baseline. While TDGP induces both a latent space and a lengthscale field, it is important to highlight that the same does not hold for CDGP and DNSGP. Thus. Fig. 4foee not show latent space for DNSGP or lengthscale field for CDGP. Notably, TDGP's s latent space perfectly captures
5

---

Function outpot
1.5
NLPD
2 -2
-2
-2
Figure 3: Synthetic 2D dataset with one latent dimension. The leftmost plot shows the true function and the location of the training data red dots. The remaining plots show the mean predictions and NLPD of each method. TDGP presents the best fit for the true function, as observed by mean and NLPD plots.
space
DKL
CDGP
DNSGP
TDGP
Latent space
X
-1 0
Identity
Inverso Nengshscclles
X
X
0.5
0
11
Figure 4: Synthetic dataset: true latent space (left), latent spaces learned by cach model (top-right), and their inverse lengthscale fields (bottom-right), Models that do not induce latent space a lengthscale field marked in red Note that TDGP is the only model that allows both interpretations. Furthermore, TDGP's latent space perfectly captures the shape the ground truth.
shape of while DKL CDGP fail to so. Analyzing the lengthscale fields, we conclude that TDGP successfully learns a non-stationary kernel, unlike DNSGP.
Finally, Fig. shows that TDGP learns to give much higher importance to one of its latent di- mensions, supporting our intuition that TDGPs lean towards learning dow-dinensional mani- folds. While the shallow GP and CDGP also weigh one dimension higher than the other, this discrepancy less accentuated. tioning that measuring the relevance of latent di- mensions depends on the architecture we evalu- the inverses
1.0
largest
0.8
0.6
0.2
worth men-
0.0
Shallow SGP DKL
CDGP TDGP
(true)
Figure 5: Synthetic dataset: what is the effective di- mensionality of the inferred latent space? A dimen- sion's relevance given by its inverse lengthscale; we plot the relative relevance of the least relevant com- pared to the largest. Only TDGP matches the one- dimensionality of the true latent space.
ate. For DGPs and DKL, of the output layer's 's kernel lengthscales. Natu- rally, the same rationale applies to shallow GPs. For TDGP, the analogous variable is the kernel variance of each hidden layer's row, since the larger this value is, the more distant from zero the values that row are. Overall, these results suggest that TDGP leans towards learning sparser representations than its competitors.
4.2 Bathymetry case study
Data. As case-study, we also apply TDGP to the bathymetry dataset GEBCO, This dataset contains a global terrain model (elevation data) for ocean and land. As an example we selected an especially challenging subset of the data covering the Andes mountain range, ocean, for
-tatiioonn task,
and land (see Appendix D.2 for details). We subsample 1,000 points from this region and compare the methods via fold cross-validation.

---

Punetion output
Figure 6: Results for GEBCO bathymetry dataset. The leftmost plot shows the dataset, while the remaining plots show the mean predictions and NLPD of each method. Note how, despite having best mean posteriors, DKL has an overfitting problem, as shown in the NLPD. TDGP is the most balanced, especially when capturing the tramsitions hetween occan/mountains/panne
DKL
CDGP
DNSGP
TDGP
Latent space
1
-2
0
-R
1.0
lengthscaler
Inverse
1
X
0.5
-2
0.0
-2
-2 -1
Figure 7: GEBCO dataset: visualization of the latent spaces learned by each model (top) and their inverse lengthscale fields (bottom). Models that do not induce a latent space or a lengthscale field are marked in red crosses. TDGP is the only model that allows both interpretations. For expert users, the lengthscale field is the only which informative priors can he put.
Results. The NLPD and MRAE results are listed in Table 2] where we observed our method to be slightly better than oth- ers. However, more importantly, this dataset can be thoroughly examined and interpreted, Figure ssowws correlation plots for a point located on the lower slope of the Andes. We observe more sensible correlations for TDGP and CDGP compared to other methods the correlation is high in the locations across the slope. This plot also highlights recurring problem with DNSGP: Despite a high inverse lengthscale barrier the An- des, there is still correlation between the slope of the moun- tain and the sea level. Additionally, Fig. 7 shows the domain coloring of learned latent spaces and the sum of the eigenvalues of the lengthscale fields. Once again, note that only TDGP can be analyzed in both ways. This is an advantage in settings where expert priors on the smoothness exist, the lengthscale field instead of the less accessible latent space
Table 2: Performance of each model on the GEBCO dataset (avg Lower is better.
NLPD
MRAE
SGP
0..1
1.19
DKL
3.85
0.39
CDGP -0.44+0.12
0.83 183=0.3
DNSGP -0.31
1.77
TDGP
-0153
0.43
these be placed directly mapping, As expected, methods that learn the lengthscale field place high inverse lengthscale values around the mountain range and low values in the smooth oceans and plains.
4.3 Benchmark datasets
Data. We also compare the methods on four well-known regression datasets from the UCI repos- itory. To assess each model fairly, we adopt ten-fold separation of the datasets into training and testing.
Results. Figure 9 shows the average NLPD for each method in the UCI datasets along with their respective standard deviations. TDGP performs either on par with other methods or outperforms them. Figure 10 also shows the relative relevance of each method's latent dimensions. Similarly
7

---

Shalloo DKL
CDGP
DNSGP
TDGP
0.5
-2
0.0
-2
-2
-2
2 -1
Figure 8: Correlation plots GEBCO dataset for datapoint marked * (1-7, 0.9), which located close mountains. observe better interpretability with TDGP (the zone high correlation coincides with moun- tains) compared to DNSGP (high correlation extends into the sea and plains).
housing
concrete
energy
wine red
SGP
SGP
CDGP
X
CDGP
DNSGP
X
DNSGP
TDGP
--00
TDGP
0.02
0.29
0.56
0.84
0.01
0.18
0.55 NLPD
1.77 1.53
1.41 NLPD
1.09
1.96
2.39
NLPD
NLPD
Figure 9: Benchmark datasets: test NLPD (lower is better) each model for 10 folds. Each dot represents result of fold. The cross and bar represent the mean and standard deviation across all folds. TDGP is as good or better than the compared models.
to what observed in the synthetic experiments, TDGP's inductive bias leads to sharp split between relevant and irrelevant dimensions compared to the other approaches. Even for the cases with similar generalization error, TDGP shows better interpretability and bias towards learning lower-dimensional representations.
TDGP
concrete housing
1
100%
0%
Most
Least
Most
Least
Most
Least
Most
Least
100%
0%
Most
Least relv.
Most relv.
Least relv.
Most relv.
Least relv.
Most relv.
Least relv.
relv.
100%
0%
Most
Least reiv.
Most relv.
Loast relv.
Most reiv.
Least relv.
Most reiv.
Least reiv.
reiv.
energy
1
100%
Most
Least relv.
Input dimennsoon
Latent dimension
Latent dimension
Latent dimension
Figure 10: Benchmark datasets: comparison of the relevances of the latent dimensions identified by each model (mean and la interval). TDGP consistently identifies low-dimensional latent spaces (most dimensions are irrelevant).
Related work
Expressive non-stationary GPs can be constructed in variety of ways. We give brief overview of some contributions available in the literature below.
Convolution of SE and scale mixture kernels. Non-stationary covariance learning can be done by convolving two SE kernels with different lengthscales. This well-studied approach was later ex- tended to arbitrary kernels by Paciorek and Schervish [9] Gibbs [12] introduced this formulation,

---

representing the lengthscale as the sum of squares of an arbitrary function applied to each covari- ance input. Later, Higdon al. [14] introduced this approach in the spatial context, varying the focus points of the ellipse representing the 2 x 2 lengthscale matrix of a SE covariance with inde- pendent GPs, which constitutes a two-layer DGP 2D. Heinonen et al. (13) introduced a method to jointly vary the covariance amplitude and lengthscale, and the observation noise over the domain, using Hamiltonian Monte Carlo for inference. As discussed in Seccionn 2] all these methods inherit the limitations of their kernel, i.e. issues with the pre-factor and violation of the triangle inequality, which hinders manifold learning.
Compositional deep GPs. Compositional DGPs 15 construct priors over complex functions by repeatedly applying a sequence of non non-paraettric transformations the input, with each trans- formation being represented by a GP with a stationary covariance. One of the the most significant issues with CDGPs is that, unlike neural networks which compose linear functions through non- linear scalar activations, every latent layer is non-linear function. Therefore, several authors have suggested improvements to DGPs by adding skip connections that directly concatenate the hidden layer with the input space [9] the use of linear mean functions [16] 23 for improved perfor- mance and capability. CDGP research also explores improvements in inference, such as the use of mini-batching and more efficient approximations to intractable integrals [23). and the inclusion of auxiliary latent variables to allow for more flexible posteriors [2].. Since TDGP is based on de- formation kernels and is related to CDGPs, we believe it can leverage many of these advances in inference, which is left to future investigations.
Deep kernel learning. An alternative to CDGPs is to use layers of parametric transformations to learn arbitrarily complex kernels. Deep kernel learning [T (66) parameterizes such transformations with a deterministic neural network. Since these deep kernels have many parameters and are not constructed in a Bayesian way, they are prone to overfitting and miscalibrated uncertainties, and can perform poorly in practice [77]. Nevertheless, by using Bayesian injective warping functions
that are well-suited for spatial applications, Zammit-Mangion et al. [27 were able to achieve good
results in low-dimensional datasets.
6 Conclusion
This work presented Thin and Deep GP (TDGP), a new hierarchical architecture for DGPs. TDGP can recover non-stationary functions through a locally linear deformation of stationary kernels. Im- portantly, this construction allows us to interpret TDGP in terms of the latent embeddings and the lengthscale field it induces. Additionally, while TDGP shares connection with CDGP, our prior does not concentrate on "flat" samples when the number of layers increases - as is the case with CDGP. Furthermore, our experiments show that TDGP performs as well as better than previous approaches. Moreover, TDGP has a robust inductive bias toward learning low-dimensional embed- dings, which is instrumental for better interpretability.
That being said, we expect TDGP will be especially useful for geospatial modeling in cases where are modeling non-stationary functions and when have expert knowledge on how this function should vary over space/time, which could be inserted as priors over the locally linear transforma- tions. We also believe that recent improvements in inference available in the literature could greatly enhance the expressivity and ease of training of TDGP.
Acknowledgments
This work was supported in part by the CONFAP-CNPQ UK Academies program (grant UKA- 001 600000003000 100003.0.000 01.00/19). Diego Mesquita was supported by the Silicon Valley Community Founda- tion (SVCF) through the Ripple impact fund, the Fundação de Amparo à Pesquisa do Estado do Rio de Janeiro (FAPERJ) through the Jovem Cientista do Nosso Estado program, and the Fundação de Amparo Pesquisa do Estado de São Paulo (FAPESP) through the grant 2023/00815-6

---

References
[1] Roberto Calandra, Jan Peters, Carl Edward Rasmussen, and Mare Peter Deisenroth. "Man- ifold Gaussian processes for regression". In: International Joint Conference Neural Net- works (IJCNN). 2016.
[2] François Chollet. Keras. httts::/kkers 10. 2015,
[3] Andreas C. Damianou and Neil D. Lawrence. "Deep Gaussian processes". In: Artificial Intel-
ligence and Statistics (AISTATS). 2013.
[4] Andreas C. Damianou, Michalis K. Titsias, and Neil D. Lawrence. "Variational inference for latent variables and uncertain inputs in Gaussian processes". In: Journal of Machine Learning Research (JMLR) 17 (2016).
[5] Mare Peter Deisenroth and Carl Edward Rasmussen. "PILCO: A Model-Based and Data- Efficient Approach Policy Search". In: International Conference on Machine Learning (ICML). 2011.
[6] Peter J. Diggle and Paulo J. Ribeiro. "Gaussian models for geostatistical data". In: Model- based Geostatistics. Springer New York, 2007, pp. 46-78.
[7] Matthew M. Dunlop, Mark A. Girolami, Andrew M. Stuart, and Aretha L. Teckentrup. "How deep are deep Gaussian processes?" In: Journal of Machine Learning Research (JMLR) 19 (2018).
[8] Vincent Dutordoir, Hugh Salimbeni, Eric Hambro, John McLeod, Felix Leibfried, Artem Artemev, Mark van der Wilk, Mare Peter Deisenroth, James Hensman, and T. John. "GPflux: A library for Deep Gaussian Processes". In: (2021).
[9] David Duvenaud, Oren Rippel, Ryan Adams, and Zoubin Ghahramani. "Avoiding pathologies in very deep networks". In: Artificial Intelligence and Statistics (AISTATS). 2014.
[10] Alexander G. de G. Matthews, Mark van der Wilk, Tom Nickson, Keisuke Fujii, Alexis Bouk- ouvalas, Pablo Léon- Vool-Villagrá, Zoubin Ghahramani, and James Hensman. "GPflow: Gaus- sian Process Library using TensorFlow". In: Journal of Machine Learning Research (JLMR) 18 (2017).
[11]
Roman Garnett. Bayesian Optimization. Cambridge University Press, 2023.
[12] Mark N. Gibbs. "Bayesian Gaussian processes for regression and classification". PhD thesis. University of Cambridge, 1997.
[13] Markus Heinonen, Henrik Mannerström, Juho Rousu, Samuel Kaski, and Harri Lähdesmäki. "Non -tttintttionary Gaussian Process Regression with Hamiltonian Monte Carlo". In: Artificial Intelligence and Statistics (AISTATS). 2016.
[14] Dave Higdon, Jenise Swall, and John Kern. "Non-Stationary Spatial Modeling". In: 6. 1999.
[15] Neil D. Lawrence and Andrew J. Moore. "Hierarchical Gaussian process latent variable mod- els". In: International Conference on Machine learning (ICML). 2007.
[16] Miguel Lázaro-Gredilla. "Bayesian Warped Gaussian Processes". In: Advances in Neural In- formation Processing Systems (NeurIPS). 2012.
[17] Sebastian W. Ober, Carl Edward Rasmussen, and Mark van der Wilk. "The promises and pitfalls of deep kernel learning". In: Uncertainty in Artificial Intelligence (UAI). 2021.
[18] Christopher J. Paciorek. "Nonstationaa Gaussian Processes for Regression and Spatial Mod- elling". PhD thesis. Carnegie Mellon University, 2003.
[19] Christopher Paciorek and Mark Schervish. "Nonstationary Covariance Functions for Gaussian Process Regression". In: Advances in Neural Information Processing Systems (NeurIPS). 2003.
[20] Carl Edward Rasmussen and Christopher K. I. Williams. Gaussian Processes for Machine Learning. MIT Press, 2006.
[21] Hugh Salimbeni. "Deep Gaussian Processes: Advances in Models and Inference". PhD thesis. Imperial College London, 2019.
[22] Hugh Salimbeni and Marc Peter Deisenroth. "Deeply non-stationary Gaussian processes". In: 2nd Workshop on Bayesian Deep Learning (NeurIPS). 2017.
[23] Hugh Salimbeni and Mare Peter Deisenroth. "Doubly Stochastic Variational Inference
for Deep Gaussian Processes". In: Advances Neural Information Processing Systems
(NeurIPS), 2017.
10

---

[24] Michalis K. Titsias and Neil D. Lawrence. "Bayesian Gaussian Process Latent Variable
Model". In: Artificial Intelligence and Statistics (AISTATS). 2010.
[25] Michalis K. Titsias and Miguel Lázaro-Gredilla, "Variational Inference for Mahalanobis Dis- tance Metries in Gaussian Process Regression" In: Advances in Neural Information Process- ing Systems (NeurIPS). 2013.
[26] Andrew Gordon Wilson, Zhiting Hu, Ruslan Salakhutdinov, and Eric P. Xing. "Deep Kernel Learning.. In: Artificial Intelligence and Statistics (AISTATS). 2016.
[27] Andrew Zammit-Mangion Tin Lok James Ng. Quan Vu, and Maurizio Filipponee "Deep Compositional Spatial Models", In: Journal of the American Statistical Association 117.540
(2022).


Reviews:
[{"review_id": "CkC8eSJVZ9", "reviewer": "Reviewer_6DJk", "summary": "The paper presents a new deep Gaussian process (DGP) model, the thin and deep GP (TDGP), which does not suffer from a diminishing signal as the number of layers increases. The crux of the TDGP model is the covariance function that, for each layer $\\ell$, acts on a linear combination of a (non-linear) transformation of the outputs of the previous layer and the inputs $\\mathbf{x}$, $\\mathbf{h}^\\ell(\\mathbf{x}) = \\mathbf{W}^\\ell(\\mathbf{h}^{\\ell-1}(\\mathbf{x})) \\mathbf{x}$. This hierarchical construction is shown to cover the usual DGP models while not suffering from diminishing signals as the depth increases. To perform inference in this model, the authors adopt a mean-field VI scheme, which is demonstrated to be effective through a number of experiments on both synthetic and real-world data.", "questions": "1. Which specific variation of the DGP model is referred to as \"CDGP\"? \n2. In lines 138-139, you say that \"By placing zero-mean priors on the entries of $\\mathbf{W}$, we encode an inductive bias to maximally reduce the latent dimensionality of $\\mathbf{h}$.\" Why is this? Intuitively, since the prior is Gaussian, it shouldn't encourage sparsity.\n3. Your contribution seems to be an input transformation, which isn't tied to a specific covariance function, so is there a reason that you chose the squared exponential (or exponentiated quadratic) covariance function? The Matérn class is more often used in practice, I believe.\n4. Is a model with more than two layers still interpretable? If so, how would you interpret the latent space?\n5. In Figures 4 and 7, how is the inverse length scale computed for TDGP? Also, why can we not compute the inverse length scales for CDGP or DKL? And why can we not visualise the latent space for DNSGP?\n6. In lines 51-52, you mention that TDGP is the only deep architecture that induces data embeddings. What do you mean by \"data embedding\" here, and why do other models not have this property? For instance, even for the original DGP model, we can easily visualise the distribution of a particular sample at a particular layer.", "limitations": "The authors have adequately addressed limitations in a dedicated paragraph.", "rating": 7, "confidence": 3, "soundness": 3, "presentation": 4, "contribution": 3, "strengths": "1. The paper presents an interesting and original approach to DGPs, an important field that should interest the NeurIPS community.\n2. The paper is of high overall quality and very well-written.\n3. The presented model has the potential to be of considerable significance, given that it generalises previous DGP models while not suffering from the same pathologies.", "weaknesses": "1. While I think the paper presents a great idea with lots of potential, my main concern is the lack of a deeper analysis of the model. Given that the proposed model is essentially just a particular transformation of the inputs at each layer and a rather straightforward mean-field VI scheme, I would expect some insights into the behaviour of the model, even if these were just empirical. The current experiments do present some interesting insights, but, say, what happens to the model's performance as we increase the depth of the model? What if we change the width of the layers? What if we use different covariance functions than the squared exponential? How does the induced covariance matrix look at different layers of the model compared to, say, DNSGP? Why does the model encourage low-dimensional latent spaces, and does this happen too for models deeper than just two layers? A Gaussian prior shouldn't encourage sparsity, so this aspect of the model is particularly puzzling.\n2. As mentioned, a standard mean-field VI scheme seems very rough given that this is known to perform poorly for many models, including DGPs. To my knowledge, the current state-of-the-art is still the doubly stochastic VI scheme by Salimbeni & Deisenroth (2017), so it would make sense to try something similar (and unless I misunderstand something, it seems fairly straightforward to do this for the proposed model). \n3. One clarity issue (perhaps the only one) I find with the paper is that it is unclear to me which DGP model the authors refer to as \"CDGP\". Deep GPs have evolved dramatically from the original formulation of Damianou & Lawrence (2013), and the mentioned pathologies are not as pronounced anymore (to my knowledge). It is great that the authors compare against DNSGP, but comparing against the doubly stochastic DGP rather than the original formulation makes much more sense. As it is unclear if the authors actually compare to the current state-of-the-art, it is also unclear if the proposed model addresses actual problems.\n4. Experimentally, the model seems to work well, but the authors use only two layers for all models. This is particularly strange as a key selling point of the proposed model is that it doesn't degenerate as the number of layers increases. But it is also an odd choice since deeper models should work better (and DNSGP demonstrably does for two of the UCI datasets). It is also strange that MRAE statistics are only reported for the GEBCO dataset; these should also be reported for the UCI experiments.\n5. The model is claimed to be \"interpretable\", but this only seems to be the case for a two-layer model (i.e., one with a single latent space). For deeper models, the distance matrix is still a highly non-linear function of the inputs, which I don't think will be particularly interpretable.\n\n\nReferences:\n\n- Andreas Damianou & Neil Lawrence, \"Deep Gaussian processes\", AISTATS 2013.\n- Hugh Salimbeni & Marc Peter Deisenroth. \"Doubly Stochastic Variational Inference for Deep Gaussian Processes\", NeurIPS 2017.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for the review. We are glad you found our paper to be interesting, significant, and of high quality. We appreciate your concerns and questions and hope you will find them addressed in the following:\n## Approximate inference and kernel choice (W2, Q3)\nWe derived our ELBO in Appendix B.3 for the two-layer model. This gives a collapsed bound which is superior to the doubly stochastic (DS-VI) scheme of Salimbeni & Deisenroth (2017) in small to medium scale regression settings, but relies on Gaussian likelihood and squared exponential kernels. However, we can train our model using the DS-VI scheme for arbitrary likelihoods, kernel functions, and depths.\n\nAs seen in Figure 2 of the rebuttal, we ran experiments with DS-VI for the TDGP using both squared exponential and Matérn-3/2 kernels; the results are worse than our original results due to the stochasticity in the DS-VI ELBO.\n## Induced covariance matrix (W1) and deeper models (W4)\nWe include a comparison of the induced covariances at different depths in the rebuttal's Figure 5. CDGP's covariances eventually \"flatten\" out due to contraction of the image, and the prior covariance gets stuck. In contrast, DNSGP and TDGP covariance can never get stuck: even if a previous layer becomes constant, this is not an issue, as each layer is always a function of the input, as the range of one layer does not depend on the previous.\n\nThus, TDGP with increasing depth does not degenerate, as opposed to CDGP. Nevertheless, in the current work we focused on deriving and experimenting with the collapsed ELBO for a 2-layer TDGP. We emphasize that a similar ELBO can still be obtained for more layers if a SE kernel and a Gaussian likelihood is considered. Further experiments with deeper models are left for future investigations due to their complexity.\n## Clarification on CDGP (W3, Q1)\nBy CDGP we refer to the general prior model with stacked GP layers; this describes the bulk of DGP works since Damianou & Lawrence (2013). Most of the advances in DGPs have been in improving approximate posterior and inference for this model class. In our experimental comparison, we indeed use the doubly stochastic DGP with a linear mean function.\n\nWe will clarify the methods used for approximate inference in each model.\n## Other metrics (W4)\nAs for the metrics, we did not report MRAE statistics for the UCI experiments due to layout constraints and we believed NLPD would be the most relevant metric to show. However, we will add a table in the appendix with the MRAE scores.\n## Why TDGP induces sparsity (W1, Q2)\nThanks for bringing this to our attention, we acknowledge that lines 136-141 are not as clear as we would like. Indeed, a Gaussian prior does not encourage sparsity. However, every entry of the $W(x)$ matrix is an independent Gaussian process with its own kernel function (Eq. 8). Each kernel has a hyperparameter $\\sigma_{w_ij}^2$ which controls the prior variance of $w_{ij}(x)$. As usually in GP models, we optimize these hyperparameters with respect to the ELBO.\n\nIf $\\sigma_{w_ij}^2$ tends to zero, then the posterior of $w_{ij}(x)$ concentrates on the prior mean. If, as discussed in lines 136-141, we set the prior mean to zero, then $\\sigma_{w_ij}^2\\to 0$ corresponds to those entries of the $W(x)\\to 0$, so a whole row tending to zero can effectively remove a dimension of the latent space. This effect is explained in lines 199-202 and it can be visualized in practice as shown in Figures 5 and 10.\n\nWe will clarify this in the final version.\n## Interpretability in deeper models (W5, Q4)\nAs argued in lines 104-119, $h^\\ell$ being a locally linear transformation means we can analyze each $W^\\ell$ as a function of the original input space $x$ with lengthscale field $W^{\\ell-1}(x)$. Take a TDGP model with three layers, $f(W^2(W^1(x) \\cdot x) \\cdot x)$, as an example. When analyzing the last layer $f$, we can reproduce the same analysis of Figures 4, 7 and 8 and plot the output of the hidden layer $W^2(x)$ to analyze the lengthscale field of $f$: areas of high inverse lengthscale correspond to areas with high variance for the derivative of $f$ and vice versa. The same analysis can be repeated for $W^2$ and _its_ lengthscale field defined through $W^1$. This analysis holds since for TDGP each hidden layer is a locally stationary function with respect to the original input space.\n\nWe will add an extended description of how to interpret the lengthscale field to our appendix.\n## Computation of lengthscale field (Q5)\nAs mentioned in line 113, $[W(x)^\\intercal W(x)]^{-1}$ is the lengthscale field, and correspondingly the local inverse lengthscale is computed as $\\Phi(x) = W(x)^\\intercal W(x)$. In Figures 4 and 7, we plot the trace of this $\\Phi(x)$ for the last layer, which reflects the smoothness of the output around each point. We have clarified this in the text.\n\nThis cannot be computed for CDGP or DKL, as both use deformation kernels of the form $k(\\tau(a), \\tau(b))$, which cannot in general be represented as a function of $(a - b)\\Phi(a,b)(a - b)$ nor of $\\| \\Phi(a)^½ a - \\Phi(b)^½ b \\|$ for some appropriate inverse lengthscale function $\\Phi$.\n\n## Latent space (Q5, Q6)\n\n> why can we not visualise the latent space for DNSGP?\n\nIf we can write the kernel as $k(a,b) = k’ (\\tau(a), \\tau(b))$, then the deformation $\\tau(x)$ defines the latent space. However, the $\\delta$ term at the core of DNSGP (see our Eq. 3 and line 86) makes this impossible.\n\n> In lines 51-52, you mention that TDGP is the only deep architecture that induces data embeddings.\n\nWe meant that TDGP is the only architecture that induces data embeddings **and** lengthscale field. We have clarified in the text that this means _at the same time_. By “data embedding” we mean that the models with successive transformations of the original data through the layers of the model, as in the DKL and CDGP models (but not in DNSGP)."}, {"review_id": "dWv3Cyfb9r", "reviewer": "Reviewer_CrCc", "summary": "This paper considers the problem of non-stationary kernel design for Gaussian processes. Two main approaches for this problem are deformation kernels and length-scale mixture kernels. On one hand, deformation kernels, eg deep GPs, have found great success while trading off expressivity and the ability to learn latent manifolds in intermediate layers, for lack of interpretability and being prone to pathologies. On the other, length-scale mixture kernels retain interpretability but typically break the triangle inequality in the induced latent space, precluding a manifold structure. The authors propose thin and deep GPs, which aim to combine the benefits of both these approaches, while mitigating their drawbacks. The approach learns a position-dependent matrix, parameterized by a GP, which learns a pointwise linear transformation of the data. The output of such a layer can then be composed to learn a new matrix acting on the original data, enabling a deep learning structure. Numerical examples are provided which appear competitive with existing techniques.", "questions": "Have you tested the method with and without the data being augmented with a bias term? Has this been tested in prior methods as well?\n\nSince the matrices are position dependent, the dimension of their range can vary across the domain as a different number of their singular values. How does this influence the latent manifold interpretation that motivated your work?", "limitations": "Besides those listed above, the computational expense is listed as a major limitation. This is a common drawback of related work, and it may be helpful to list the memory and computational complexity of this method with L layers, as is commonly done in the support literature.", "rating": 4, "confidence": 4, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "This paper is clearly presented and combines existing ideas in a novel way. The numerical examples seem well selected, with a mixture of synthetic/analytic examples that demonstrate the utility of the method, with standard benchmark datasets allowing for meaningful comparisons on real data. A case study is also provided that is geared towards the strengths of the technique.", "weaknesses": "The main weakness of the paper in my view is Theorem 3.1. It appears that the proof that the CDGP is a generalization of deep GPs hinges upon the introduction of a new variable that plays the role of the deep GP, and then turning off their method's parameters by concentrating their priors around zero. This seems highly misleading unless the method is introduced in a way that uses the augmented space as default, or some numerical evidence is given that it should be turned off if appropriate. In this case, presumably the deep GP should also be tested with and without an augmented space for fair comparison. It seems that the augmented space is not used in any of the experiments in any of the models presented, since the local linearity leaving the data unchanged around zero is listed as a limitation.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for the comments. \n\nWe would like to highlight that the bias term is a simple design choice (used profusely throughout ML) that suffices to show TDGP is at least as expressive as CDGP. Both CDGP and TDGP with no bias are limiting cases of TDGP with bias.\n\n> Have you tested the method with and without the data being augmented with a bias term? Has this been tested in prior methods as well?\n\nEmpirical performance may often rely on factors that are alien to theoretical analysis — e.g., inference quality, initialization, learning rate scheduling. In particular, adding bias terms increases the number of parameters of each layer $\\ell$ in $Q_\\ell$, which makes inference more challenging and is the main reason we did not test it before the submission. We have run TDGP with bias on the synthetic experiment using doubly-stochastic inference. As expected, results were significantly worse when we included the bias (NLPD: 0.7 vs -0.3). However, if we reduce the number of parameters  restricting $W(\\cdot)$ to be diagonal, TDGP w/ bias outperforms TDGP w/o bias (NLPD: -2.7 vs 1.8). This corroborates our hypothesis.\n\n> Since the matrices are position dependent, the dimension of their range can vary across the domain as a different number of their singular values. How does this influence the latent manifold interpretation that motivated your work?\n\nFormally, to characterize a manifold, we need the number of singular values of $W(x)$ to be a constant for all $x$ in the domain. The restriction that the local dimensionality must be constant is usually not  enforced in most deep learning models. In particular, this also applies to CDGP. However, even if this quantity fluctuates, it is still reasonable and commonplace to use the latent embedding plots to interpret the model. Nonetheless, unlike CDGP, our model provides a straightforward way to quantify the local dimensionality. Therefore, we also include Figure 3 with the distribution of non-zero singular values in the data for different folds and datasets and also Figure 4 with the average singular values in the data, showing their correlation with Figure 10 in the paper.\n\n>This is a common drawback of related work, and it may be helpful to list the memory and computational complexity of this method with L layers, as is commonly done in the support literature\n\nThanks for the suggestion. We will include the spatial and temporal complexity analysis to the final manuscript, as they are only implied by the description of the model. Just as shown in Figure 1 of the rebuttal, we expect increases in depth to linearly increase the time spent in the optimization of the ELBO.\n"}, {"review_id": "2k2ES3eNqC", "reviewer": "Reviewer_HPwG", "summary": "This paper extends the CDGP to address two weaknesses in deep Gaussian processes. The solution is similar intuitively to a residual connection. Instead of letting each layer depend only on the last layer, they let each layer depend on the last layer and the input. This allows the model to induce a manifold and use a lengthscale field.", "questions": "Most of my questions are in the weaknesses.", "limitations": "They have a good discussion of limitations: mainly that it is computationally expensive and leads to requiring more parameters.", "rating": 4, "confidence": 4, "soundness": 3, "presentation": 2, "contribution": 3, "strengths": "This seems like an important problem, and the solution seems principled.", "weaknesses": "There should be more discussion of what it means to not induce a manifold, why violating the triangle inequality leads to that, and why not inducing a manifold is a drawback. I think inducing a manifold implies that one can make the latent space plots: why are these valuable? Why do they require the triangle inequality?\n\nThere should be more discussion of how a lengthscale field leads to interpretability.\n\nIn both the latent space and the inverse lengthscale plots, there should be some discussion about why we want to see these plots: for instance, in figure 7, we have the two sets of plots: what can we glean from them about the data/problem/etc? Figure 8 has a little bit (the zone of high correlation extends into the mountains), but this is not obvious for a non-expert application reader. You could have a short tutorial perhaps of how to read these plots in the supplement.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for your suggestions to improve our manuscript. We agree that elaborating further on the interpretation of the experiments can make it more accessible to the general public. In the next paragraphs, we elaborate further on the importance of both types of representation (manifold/latent space and lengthscale fields). While we cannot update the manuscript at this point, we will include this discussion in the camera-ready version along with a guide on how to read the plots.  We hope our answers have sufficiently addressed your concerns.\n\n>Manifolds, latent space plots and the triangle inequality\n\nIndeed, inducing a manifold implies there is some function $\\tau$ mapping inputs to a finite-dimensional latent space. We can use the plots to analyze how $\\tau$ changes the notion of neighborhood (or similarity) between points in the input space. For instance, Figure 7 reveals that TDGP learns to differentiate the mountains from both the landmass above and the water below, mapping the mountains to a different range of values. On the other hand, Figure 7 also shows that DKL and CDGP do not discriminate between the mountain range and the landmass. In general, we can use these plots to analyze $\\tau$ and help human experts validate/understand/debug our models. It is not clear how this could be done for models without an explicit $\\tau$, like DNSGP. In the case of DNSGP, the quadratic term in our Equation 3 ($\\delta$) is what makes it hard to write the kernel as $k(\\tau(a), \\tau(b))$. To rewrite $\\delta$ as a quadratic form in $\\tau(a)$ and $\\tau(b)$, we would need to redefine $\\tau$ to take $a$ and $b$ at the as inputs at the same time, as discussed by Paciorek [19, page 32]. If $\\delta$ induced a metric, writing DNSGP in terms of some $\\tau$ would be simpler. Importantly, while it might be possible to extract some $\\tau$ from DNSGP, how it is possible or how to do so is an open problem, as far as we are aware.\n\n>Importance of lengthscale fields for interpretability\n\nThe inverse lengthscales allow interpreting the prior smoothness of a function at each input point. For instance, areas in which $W(x) W(x)^\\intercal$ is “high” (sum of eigenvalues)  correspond to areas where the function should be less smooth, i.e. have increased variance in their derivatives. More concretely, the lengthscale plots in Figure 7 shows that TDGP captures the non-smoothness in the transitions between the landmass, the mountain and the water. A similar phenomenon happens for DNSGP. This is in contrast with the latent space plots which don't directly relate to properties of the learned function.\n\nWe will add an extended description of how to interpret the lengthscale field to our appendix, as this might not be widespread knowledge."}, {"review_id": "rKILPdyMxC", "reviewer": "Reviewer_DvNb", "summary": "Techniques such as Deep Kernel Learning and extensions such as Deep GPs have been proposed in order to overcome the flexibility constraints associated with standard shallow GPs. However, although these techniques are better suited to model non-stationary data, they are still susceptible to issues such as pathologies when extended to more than a few layers, as well as overfitting. These techniques also lack the interpretability of the latent space and associated length-scales that make standard GPs so appealing. In order to mitigate these limitations, the authors propose a novel approach that generalises over the more widely-used compositional DGPs via a hierarchical model that incorporates locally linear deformations of stationary kernels. A synthetic example is crafted in order to showcase how TDGP improves over competing methods for problems exhibiting heavy non-stationarity. This is further complemented by an analysis on a real-world dataset, as well as 4 benchmark datasets from the UCI repository.", "questions": "I do not have additional questions beyond the points raised in the *Weaknesses* section. I would appreciate if the authors could read through this and indicate whether they can improve upon these aspects in a future revision.\n\nThe paper is a solid addition to the literature on hierarchical GP models – besides proposing a novel method, it also gather and illustrates potential deficiencies in existing methods that are widely-used in the community. I also believe that the paper is already in a state that is fit for publication, and can only be improved further with additional minor iterations.", "limitations": "N/A", "rating": 7, "confidence": 4, "soundness": 3, "presentation": 4, "contribution": 3, "strengths": "- The paper is very well-written and a pleasure to read. I commend the authors for structuring the paper in a manner that clearly shows connections to related models, while also showcasing the key contributions of this work.\n- I particularly appreciated Sections 4.1 and 4.2, which capture the various dimensions along which TDGP can lead to improvements over other methods (predictive performance, uncertainty quantification, and interpretability).\n- The problem statement tackled in this paper is well-motivated, and the limitations of other techniques are clearly pointed out both verbally and visually in the experiments section. One possible improvement to this could be highlighting application domains where TDGP is expected to be most impactful – i.e., which application domains tend to non-stationary data where well-calibrated uncertainty is especially useful. While fairly minor, this could help clarify the significance of this work earlier on the paper.  ", "weaknesses": "- The authors comment in the *Limitations* section that the increase in the number of variational parameters could slow down optimisation – it could be worth supporting this statement by specifying the computational complexity, and maybe an analysis reporting wall clock time for one of the experiments in Section 4. \n- I would also be interested in whether the increased parametrisation of the model results in models that are less stable or difficult to consistently converge. This has implications on the practical utility of the model.\n- The paper focuses heavily on using the SE kernel throughout – can this be extended to other kernels or is the value-add of this extension not considered to be worthwhile for the paper?", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for your detailed feedback. We have run additional experiments to answer your questions (below). Besides including the additional results and complexity analysis to the final manuscript, we will also add a discussion on which domains we expect TDGP to be most useful, as suggested. \n\n>The authors comment in the Limitations section that the increase in the number of variational parameters could slow down optimisation – it could be worth supporting this statement by specifying the computational complexity, and maybe an analysis reporting wall clock time for one of the experiments in Section 4.\n\nThanks for the suggestion. For simplicity, assume the number of latent dimensions is constant across layers. Regarding memory complexity, TDGP uses $L\\times D\\times Q\\times m$ variational parameters for the mean and $L\\times D\\times Q\\times m^2$ for the covariance. As for the computational complexity one ELBO evaluation costs $O(L\\times D\\times Q\\times m^3)$ for the KL terms and $O(L\\times D\\times Q\\times m^2\\times n)$ for the expectation term. Although we did not link it properly in the main text, we have a wall-clock time analysis for all experiments in Appendix C, Tables 3 and 4. We will include the complexity analysis to the final manuscript. We will also add a reference to the wall-clock tables in the experiments section.\n\n>I would also be interested in whether the increased parametrisation of the model results in models that are less stable or difficult to consistently converge. This has implications on the practical utility of the model.\n\nAs discussed in our limitations, we expect models with a higher number of parameters to be more difficult to optimize. However, Figure 1 in the rebuttal shows that TDGP is reasonably stable for moderately-sized values of $Q$. Since you have also asked for wall-clock times (above), we also included a complimentary plot with time as a function of $Q$ in Figure 1. As TDGP has an architecture that is not comparable to previous DGP models, best practices for learning rate scheduling and pre-training for increased robustness are still unknown and we intend to explore this topic in the future.\n\n>The paper focuses heavily on using the SE kernel throughout – can this be extended to other kernels or is the value-add of this extension not considered to be worthwhile for the paper?\n\nThe inference scheme that we used in our experiments and as derived in Appendix B.3 is only valid for SE kernels due to the necessity of the computation of the $\\Psi$-statistics, which are convolutions of the kernel with Gaussian densities. However, using solely the SE kernel is commonplace in the Deep GP literature, as the presence of the hidden layer alleviates the drawbacks of such a simple kernel. For example, the experiments done by Salimbeni & Deisenroth (2017) only consider SE kernels.\n\nNonetheless, we can train our model using the doubly stochastic VI (DS-VI) scheme of Salimbeni & Deisenroth (2017), which allows us to benchmark different kernel choices. Note, however, that DS-VI is stochastic and black-box. Therefore, compared to our ELBO, DS-VI results in a loose lower-bound on the evidence. For completeness, we have run TDGP using DS-VI on the benchmark datasets for the SE kernel and the Matérn 3/2 kernel (Rebuttal, Figure 2). As expected, results were worse than with our tailored inference scheme."}]

Author rebuttal:
We thank the reviewers for their thoughtful feedback, as we believe addressing them will broaden the audience for this paper. We are also glad that most reviewers feel our ideas are clearly presented and constitute a solid and principled contribution to this problem space.

In particular, the considerations shared across reviewers which will be revised in the final manuscript are:
* further explanation on the interpretability of our model and lengthscale fields in general,
* inclusion of big-O notation complexity analysis,
* further explanation on the possible use of alternative kernels and inference methods.

As for additional plots and experiments, our attached PDF includes experimental analysis relating to: 
* (Figure 1) convergence and wall clock timing as a function of TDGP’s width,
* (Figure 2) additional experiments with alternative inference schemes and non-SE kernels,
* (Figures 3 and 4) an analysis of the local dimensionality of TDGP’s latent space,
* (Figure 5) an illustration of the behavior of the prior covariance for different numbers of hidden layers.

In each review’s personal rebuttal, we address in detail each question and concern.


Meta-review: New approaches to non-stationary kernel design and (variational) inference are presented. Contributions and limitations are described and it is well written. However, the rebuttal process has lead to better explanations/presentation suggestions; I strongly suggest updating the final version that incorporates it carefully.
Decision: Accept (poster)