Implicit Manifold Gaussian Process Regression
Authors: Bernardo Fichera
Viacheslav Borovitskiy
Andreas Krause
Aude Billard
Venue: NeurIPS 2023 poster (2023)
Keywords: Gaussian process
manifolds
manifold learning
uncertainty
regression
graph Laplacian

Abstract: Gaussian process regression is widely used because of its ability to provide well-calibrated uncertainty estimates and handle small or sparse datasets. However, it struggles with high-dimensional data. One possible way to scale this technique to higher dimensions is to leverage the implicit low-dimensional manifold upon which the data actually lies, as postulated by the manifold hypothesis. Prior work ordinarily requires the manifold structure to be explicitly provided though, i.e. given by a mesh or be known to be one of the well-known manifolds like the sphere. In contrast, in this paper we propose a Gaussian process regression technique capable of inferring implicit structure directly from data (labeled and unlabeled) in a fully differentiable way. For the resulting model, we discuss its convergence to the Matérn Gaussian process on the assumed manifold. Our technique scales up to hundreds of thousands of data points, and improves the predictive performance and calibration of the standard Gaussian process regression in some high-dimensional settings.

Implicit Manifold Gaussian Process Regression
Bernardo Fichera'
Viacheslav Borovitskiy
Andreas Krause
Aude Billard
EPFL
ETH Zürich
Abstract
Gaussian process regression widely used because of its ability to provide well-calibrated uncertainty estimates and handle small or sparse datasets. However, it struggles with nigh-dimensional data. One possible way to scale this technique higher dimensions is to leverage the implicit fow-dimensional manifold upon which the data actually lies, as postulated by the manifold hypothesis. Prior work ordinar- ily requires the manifold structure to be explicitly provided though, i.e. given by a mesh or known be one of the well-known manifolds like the sphere. con- trast, in this paper propose Gaussian process regression technique capable of inferring implicit structure directly from data (labeled and unlabeled) in fully dif- ferentiable way. For the resulting model, discuss its convergence to the Matérn Gaussian process on the assumed manifold. Our technique scales up to hundreds thousands of data points, and improves the predictive performance and calibration the standard Gaussian process regression in some high- nigh-dimensional settings.
Introduction
Gaussian processes are among the most adopted models for learning unknown functions within the Bayesian framework, Their data efficiency and aptitude for uncertainty quantification make them appealing for modeling and decision-making applications in the fields like robotics (Deisenroth and Rasmussen, 2011) 1), geostatistics (Chilès and Delfiner, 2012), numerics (Hennig et al., 2015), etc.
The most widely used Gaussian process models, like squared exponential and Matérn (Rasmussen and Williams, 2006), impose the simple assumption of differentiability of the unknown function while also respecting the geometry of IRd by virtue of being stationary or isotropic. Such simple assumptions make uncertainty estimates reliable, albeit too conservative times. The same simplicity makes these models struggle from the curse of dimensionality. We hypothesize that it is still possible leverage these simple priors for real world high dimensional problems granted that they are adapted to the implicit low-dimensional submanifolds where the data actually lies, as illustrated by Figure 1.
Prediction
Unvertuinty
(a) Euclidean, prediction
(b) Ours, prediction Ec)cucci,,,n, uncertainty
(d) Ours, uncertainty
Figure Fuclidean (standard Matérn-5/2 kernel) vs ours (implicit manifold) Gaussian process regression for data that lies dumbbell-shaped curve (1-dimensional manifold) assumed unknown. The data contains a small set of labeled points and large set of unlabeled points. Our technique recognizes that the lines in the middle intrinsically far away from each other, giving much better model on and near the manifold. Far away from the manifold reverts to the Euclidean model.
Code available https:://giiiii http:/ghhwb..ccc////aaa////aaaaa gp.
37th Conference on Neural Information Processing Systems (NcurIPS 2023).

---

Recent works in machine learning generalized Matérn Gaussian processes for modeling functions on -uclluei domains such as manifolds or graphs (Azangulov al., 2022: 2023; Borovitskiy et al., 2021; 2023; 2020). Crucially, this line of work assumes known geometry (e.g. a manifold or a graph) beforehand. In this work we aim to widen the applicability of Gaussian processes for higher dimensional problems by automatically learning the implicit low-dimensional manifold upon which the data lies, the existence of which is suggested by the manifold hypothesis. We propose a new model which learns this structure and approximates the Matérn kernel on the implicit manifold.
Our approach can operate in both supervised and cumi-suprrviee settings, with the emphasis on the latter: uncovering the implicit manifold may require lot of samples from it, however these samples need not be labeled, and unlabeled data usually more abundant. Taking inspiration in the manifold learning results of Coifman and Lafon (2006), Dunson et al. (2021) and others we approximate the unknown manifold by an appropriately weighted nearest neighbor graph. Then we use graph Matérn kernels thereon as approximations to the manifold Matérn kernels of Borovitskiy et al. (2020), extending them to the vicinity of the manifold in the ambient Rd an appropriate way.
1.1 Related Work and Contribution
High-dimensional Gaussian process regression is an area of active research, primarily motivated by decision making applications like Bayesian optimization. There are three main directions in this area: (1) selecting a small subset of input dimensions, (2) learning a small number of new features by linearly projecting the inputs and (3) learning non-linee non-linear features. Our technique belongs to the third direction. Further details on the area can be found in the recent review by Binois and Wycoff (2022).
A close relative of our technique in the literature is described in Dunson et al. (2022). It targets the low-dimensional setting where the inputs are densely sampled on the underlying surface. is based on the heat (diffusion) kernels on graphs as in Kondor and Lafferty (2002) and uses the Nyström method to extend kernels to Rd, both of which may incur a high computational cost. Another close relative is the concurrent work by Peach et al. (2023) modeling vector fields on implicit manifolds.
We are targeting the high-dimensional setting. Here, larger datasets of partly labeled points are often needed to infer geometry. Because of this, we emphasize computational efficiency by leveraging sparse precision matrix structure of Matérn kernels (as opposed to the heat kernels) and use KNN for sparsifying the graph and accelerating the Nyström method. This results in linear computational complexity with respect to the number of data points. Furthermore, the model we propose fully differentiable, which may be used to find both kernel and geometry hyperparameters by maximizing the marginal likelihood. Finally, to get reasonable predictions on the whole ambient space Rd, we combine the prediction of the geometric model with the prediction of a classical Euclidean Gaussian process, weighting these by the relative distance to the manifold.
The geometric model is differentiable with respect to its kernel., likelihood- and geometry-related hyperparameters, with gradient evaluation cost being linear with respect to the number of data points. After training, we can efficiently compute the predictive mean and kernel as well as sample the predictive model, providing the basic computational primitives needed for the downstream applications like Bayesian optimization. We evaluate our technique on a synthetic low-dimensional example and test in a high-dimensional large dataset setting of predicting rotation angles of rotated MNIST images, improving over the standard Gaussian process regression.
2 Gaussian Processes
A Gaussian process f ~ Gm,m, k) is distribution over functions on set X. It determined by the mean function m(x) = E f(x) and the covariance function (kernel) k(xx, x') = (fv(f()),
Given data X, y. where X = (æ) 2n.) and y = (Wi, ..: Vnn) with 2;; E X yi E IR, one usually assumes yi = =xx +6; where 6; ~ (O,oo)) is IID noise and f ~ GP k)) is some prior Gaussian process, whose mean is assumed to be zero in order simplify notation. The posterior distribution fI y is then another Gaussian process f y ~ GP k) with (Rasmussen and Williams, 2006)
mm(-) = {() (Kxx + o y, k(,,-) = K(,,) (Xxx(Xxx +o21) Kx('), (1) where the matrix Kxx has entries k( 2.)), the vector Kx() = KTix has components k( :). If needed, can efficiently sample f using pathwise conditioning (Wilson et al., 2020; 2021)
2

---

(a) Kernel, sphere $
(b) Kernel, graph P 2
(c) Sample, sphere 82
(d) Sample, graph P82
Figure 2: Kernel values )) and samples for the Matérn Gaussian processes the sphere manifold 82 and for the approximating Matémn5/2 process on a geodesic polyhedron graph P S2.
Matérn Gaussian processes-inchuding the limiting - 0o case, squared exponential Gaussian processes- are the most popular family of models for X' = Rd. These have zero mean and kernels
k (x, x') = 0² 222-0 T(v)
T x'1l
x x'll
22
Kr
2v
(2)
where K2 the modified Bessel function the second kind (Gradshteyn and Ryzhik, 2014) and v, K. o² are the hyperparameters responsible for smoothness, length scale and variance, respectively. We proceed to describe how Matérn processes can be generalized to inputs æ lying on explicitly given manifolds graphs instead of the Euclidean space IR
2.1 Matérn Gaussian Processes on Explicit Manifolds and Graphs
For domain which is Riemannian manifold, an obvious and natural idea for generalizing Matérn Gaussian processes could be to substitute the Euclidean distances ||a x' Equation (2) with the geodesic distance. However, this approach results in ill-defined kernels that fail to be positive semi-definite (Feragen et al., 2015; Gneiting, 2013).
Another direction for generalization is based on the stochastic partial differential equation (SPDE) characterization of Matérn processes first described by Whittle (1963): GP xa) solves
2v A2
+
Agd
f=W,
(3)
where Ag is the standard Laplacian operator and W is the Gaussian white noise with variance proportional to o². If taken to be the definition, this characterization can be easily extended to general Riemannian manifolds X' = M by substituting Aga with the Laplace-Beltrami Beltrami operator AM. taking d = dim M and substituting W with the appropriate generalization of the Gaussian white noise (Lindgren al., 2011). Based on this idea, Borovitskiy al. (2020) showed that on compact Riemannian manifolds, Matérn Gaussian processes are the zero-mean processes with kernels
4/2
ky (I, x') = Cu.
n² +
f((x)(((').
(4)
1=0
where fi fi eigenvalues and eigenfunctions of the Laplace- Beltrrai operator and Cv.. is the normalizing constant ensuring that & .x kvvn. (æ, 2) dz = =² This, alongside with considerations from Azangulov et al. (2022) allows one to practically compute kv,r, 2 for many compact manifolds.
If the domain X' is weighted undirected graph G, we can also use Equation (3) to define Matérn Gaussian processes on & (Borovitskiy et al., 2021). this case, Ag is substituted with the minus graph Laplacian -g and W ~ N(0, o I) is the vector of IID Gaussians. Here, SPDE transforms into a stochastic linear system, whose solution of the same form as Equation (4) but with a finite sum instead of the infinite series, with d 0 because there is no canonical notion of dimension for graphs and with At, fi being the eigenvalues and eigenvectors- functions on the node set-of the matrix Ag. These processes illustrated on Figure 2.
3

---

3 Implicit Manifolds and Gaussian Processes on Them
Consider a dataset X (N) 2; R partially labeled with labels y1, Un E n N. Assume that x are IID randomly sampled from compact Riemannian submanifold M S M As by Section 2.1, the manifold M is associated to family of Matérn Gaussian processes tailored to its geometry. We do not assume know M, only the fact that it exists, hence the question is: how can we recover the kernels of the aforementioned geometry-aware processes from the observed dataset?
It is clear from Equation (4) that to recover ku,x,0" need to get the eigenpairs - ft of the Laplace-Bell Beltramii operator on M. Naturally, for a finite dataset this can only be done approximately. We proceed to discuss the relevant theory of Saplace-Beltramt eigenpair approximation.
3.1 Background on Approximating the Eigenpairs of the Laplace-Beltrami Operator
There exists number of theoretical and empirical results on eigenpair approximation. Virtually all of them study approximating the implicit manifold by some kind of weighted undirected graph' with node set [., ØN) and weights that are somehow determined by the Fuclidean distances ||z z. The eigenvalues of the graph Laplacian on this graph are supposed to approximate the as functions on the
operator, while the eigenvectors- node set-approiimate the values of the eigenfunctions of the Laplace-Beltrami operator at æ; M. To approximate eigenfunctions elsewhere, any sort of continuous (smooth) interpolation suffices.
cigenvalues of the Haplace-Beltrane
There are three popular notions of graph Laplacian. Let us denote the adjacency matrix of the weighted graph by A and define D to be the diagonal degree matrix with D;; Σ, A.j. Then
A Aun A
sym - 1/2 AD 1/2
rw - 'A.
(5)
unnommalized
symmetric normalized
random walk normalized
The first two of these are symmetric positive semi-definite matrices, the third is, generally speaking, non-symmetric. However, from the point of view of linear operators, all of them can be considered symmetric (self-adjoint) positive semi-definite: the first two with respect to the standard Fuclidean inner product (., +), and the third one with respect to the modified inner product (v, u)D = (Dv, u). Thus for each there exists orthonormal basis of eigenvectors and eigenvalues are non-negative."
The most common way define the graph by setting Auj exp( -|,,,-1;11///0 for an a 0. If x; are IID samples from the uniform distribution on the manifold M, then all of the graph Laplacians, each multiplied by an appropriate power of a, converge the Laplace-Beltraml operator, both pointwise (Hein et al., 2007) and spectrally (García Trillos et al., 2020). i.c. in the sense of eigenpair convergence, at least at the node set of the graph. However, the inputs æ; are sampled non-un non-uniformly, graph Laplacians, at best, converge to different continuous limits, none of which Saplace-Beltrami coincides with the Laplace- Beltrami operator (Hein al., 2007).
Coifman and Lafon (2006) proposed a clever trick to handle non-uniformly sampled data I N... Starting with Á and D defined in the same way as A and D before, they define A = 1D-1AD Intuitively, this corresponds normalizing by the kernel density estimator cancel out the unknown density. The corresponding A then converges pointwise to the Laplace-Beitrami operator (Hein al., 2007), though A un and Aaym do not: they converge different continuous limits. Dunson et al. (2021, Theorem 2) show that, under technical regularity assumptions, eigenvalues d and renormalized eigenvectors of Arw converge the respective eigenvalues and eigenfunctions of the Laplace BentreamBBlraiii operator, regardless of the sampling density of z;,
Both in the simple case and in the sampling density independent case, the graphs and their respective Laplacians turn out to be dense, requiring a lot of memory to store and being inefficient to operate with. To make computations efficient, sparse graphs such as KNN graphs much more preferable over the dense graphs. Spectral convergence for KNN graphs is studied, for example, in Calder and Trillos (2022), for Aij h(1l= with a compactly supported regular function and with
'Otee possibilities include linear (classical PCA) or quadratic (Pavutnitskly et al., 2022) approximations. See the books by Lee and Verleysen (2007) and Ma and Fu (2011) for additional context.
NNaturally, for the random walk normalized Laplacian Arw the orthonormality with respect to (:, :) *Garefa Trillos et al. (2020) do not explicitly study Asym Since 4.ym and A.w are similar matrices, they share eigenvalues, eigenvalue convergence for Asym is trivial, the eigenvector convergence, however, is not.

---

limit depending on the sampling density. Unfortunately, we are unaware of any spectral convergence results in the literature that hold for KNN graphs and are independent of the data sampling density.
3.2 Approximating Matérn Kernels on Manifolds
Here we incorporate various convergence results, including but not limited to the ones described in Section 3.1, proving that all spectral convergence results imply the convergence of graph Matérn kernels to the respective manifold Matérn kernels.
Proposition 1. Denote the eigenpairs by At fi for a graph Laplacian and by X) Si for the Laplace- Beltrami operator. Fix O. Assume that, with probability least 1 5. for all E fo a small enough and for N large enough we have { xX1 < E and fi(æ.) - f/ (æ,)) < E. Then, with probability at least 1 5, we have k, (2); O (2:; T.) as a 0, N, L oc, where
o²
2
- -vddi((M)/2
k p.n,o (x),,
t
fi(z))fiaa;)
C2,N
Proof. First prove that the tail of the series in Equation (4) converges uniformly zero, then combine this with eigenpair bounds. See details in Appendix A.
Remark, The convergence in æ, E M can be lifted to pointwise convergence for all æ E M eigenvectors are interpolated Lipschitz-continuuusly. simply because the eigenfunctions are smooth.
Inspired by this theory, we proceed to present the implicit manifold Gaussian process model.
3.3 Implicit Manifold Gaussian Process
Guided by the theory we described the previous section we are now ready to formulate the implicit manifold Gaussian process model. Given the dataset I1; EN IR and f/..., Ur. IR put
||=;;
xz 112 1²).
E $000m 00
Aim
i i= j.
A A - À.j- Sx Sx(x,, x.) exp
D.
(7)
4a²
i i
Here (2,; ..y) if n; is of the K nearest neighbors of a or vice versa and S (a:; 0.)) otherwise; all matrices are of size N x N and depend on a and K as hyperparameters. Thanks to the coefficient (x 2.y) that performs KNN sparsification, the matrix A is sparse when K < N.
Then we consider the operator A.w = I-D- -1 A defined by Equation (5), whose matrix is also sparse.
Denoting its sigenvalues-orden from the smallest to the largest- by 0 = Åo S A < ... S A and its eigenvectors- gngervecors--orhhooonmee orthonormal under the modified inner product (- :) p and regarded as functions on the node set of the graph- Jo.... _N-1., we define Matérn kernel on graph nodes n; by
X M,K 0 23))
0²
Σ 1=0
1--1
Du.n (M))))))
P (A)
(2v Ar
+
(8)
CU.
where L does not need to be equal to the actual number N of eigenpairs. Doing so means truncating the high frequency eigenvectors (fr for / large), which always contribute less the sum because they correspond to smaller values of & (At). This can massively reduce the computational costs.
By Proposition 1, Equation (8) approximates the manifold Matérn kernel with smoothness '' = - dim(M) /2. We adopt such a reparametrization because does not require estimating the a priori unknown dim(MM). This, however, makes the typical assumption of (1/2, 8/2, 5/2) in- adequate. We chose a particular graph Laplacian normalization, namely the random walk normalized graph Laplacian Arw. to approximate the true Laplace- Beltrami operator regardless of the potential non-uniform sampling of æ NN, based on the theoretical insights described in Section 3.1.
The kernel in Equation (8) is only defined the set of nodes I,, next step is to extend it to the whole space Rdd Extending kernels is usually a difficult problem because one has to worry about positive
'Note: An 1, if the graph has loops. Assuming As 0 would lead to discontinuities later stages.
5

---

اممامم
امن
(a) Extended eigenvector f.(+)
(b) Extended kemel 1.x
n, (,,
(c) Function (dist((;.))
Figure 3: Different quantities connected to kernel extension. Notice that the values on subfigures (a)
and (b) artificially restricted to the set dist((,,M)) 3a maintain numerical stability.
semi demi-definitnness. To work around it, we extend the features ft. For this, we use Nyström method: we allow the first argument of SK be an arbitrary vector from R4, defining Sx (ar, myy) if æ is one of the nearest neighbors of æ among æ1..., ZN. This allows us to extend A, D, A and D as
| x,
N
A(x,0,) = Sx(xx, æt,) xx) exp
D(Ø) = Σ
2() = E A(x,0;), (9) 2;EKNN(#)
4a
==1
A(x,2,) = Ålæ, æzi)
N
D(x) = Σ A(x,
Σ A(x, (10)
D(x) D(x,)
,=1
2,EKNN(m)
where KNN(②) is the set of the K nearest neighbors from æ among 21,
With this, we define
1 A(,,0)) 1 D(x) f((ar))
A(x, D(x) f((2)).
f1(x))
Σ EKNN(())
(11)
1
It is easy check that this extension perfectly reproduces the values fi(x)). simply because A(x, æz)) coincides with Aij when 2; and because (i,
is the eigenvector of A corre- sponding to the eigenvalue - A₁. It also allows us to extend the kernel 1-X "Mn,, as well, such that k- k "xxKo= (x, æ') is defined for arbitrary æ, x' E Idd and Equation (8) still holds for the nodes Ui...,aN
We visualize an extended eigenvector and an extended kernel in Figures 3a and 3b.
For æ far away from the nodes æ; the values of Dxx) become very small, making the extension procedure numerically unstable. Furthermore, the geometric model is generally not so relevant far away from the manifold. Because of this, the final predictive model ((x) ~ (P(m(r), k(P)) combines
the geometric model f(m) ~ (m(m(()) (n)_))))))()))))))))))))) ) posterior under the kernel y.X "2,x,0² we defined
just above-wirt the standard Euclidean model f(e) GP(m( ((()) the posterior under the standard Euclidean Gaussian process in R.d, for instance with the squared exponential kernel:
(3a))
(( (=)) (m) + (o)
r(x) exp
(3x)) dist M)) (12)
where 7 is bump function that is zero outside the 3a neighborhood of the manifold, illustrated in Figure 3c. Here dist(x, Je.M .M) can be computed as the distance from æ to its nearest neighbor node zy..
Efficient Training of the Implicit Manifold Gaussian Processes
Here describe how perform the implicit manifold Gaussian process regression efficiently, being able to handle hundreds of thousands points, in both supervised and semi-supervised regimes.
In all cases need efficient (approximate) KNN build graph, extend the kernel beyond the nodes x; and combine the geometric model with the standard Euclidean one. For this we use FAISS (Johnson al., 2019). The resulting sparse matrices, such the Laplacian Arw. represent as black box functions capable of performing matrix-vectis multiplications for any given input vector.
'In practice it makes sense to compute dist(æ, as the average of distances between æ and its K nearest neighbors to smoothen the resulting y(x) this computationally cheap given efficient KNN implementation.

---

After hyperparameters are found-we we will return to their search later-we need to compute the eigenpairs Ar, fi of Arw- For this we run Lanczos algorithm (Meurant, 2006) to evaluate the eigenpairs X ft of the symmetric matrix syms putting À = A and fi D D-1 flym because the matrices Aw and Asy are similar, i.e. Arw D D-1/22 A sym D1/2, Importantly, Lanczos algorithm only relies on matrix-vector products with Asym- We only compute few hundred of cigenpairs, asking Lanczos to provide twice or trice as many and disregarding the rest.
When there is lot of labeled data, we approximate the classical Euclidean (e.g. squared exponential) kernel using random Fourier features approximation (Rahimi and Recht, 2007), this allows linear computational complexity scaling with respect to the number of data points. The number of Fourier features taken be equal L, with the same L as Equation (8).
As it was already mentioned, we need to find hyperparameters θ (,, that determine the graph, Gaussian process prior and the noise variance that fit the observations best. The parameter we assume manually fixed. To avoid nonsensical parameter values-a common difficulty often occurring when the data is scarce one might want to assume prior p(0) on O. Some specific choices of p(0) are discussed in Appendix C. Then 0 is maximum posterioni (MAP) estimate:
i arg max, log p(y 0, X) + log p(0).
(13)
To simplify hyperparameter initialization and align with zero prior mean assumption it makes sense to preprocess y, to be centered and normalized.
To solve the optimization problem Equation (13) we use restarted gradient descent. Repeatedly evaluating the gradient of log pty 0, X) is the main computational bottleneck. 'he key idea for doing this efficiently- viable for integer values of v-is reduce matrix-vector products with Matérn kernels' precision to iterated matrix-vector products with the Laplacian, which is sparse. First, describe this in detail in the noiseless supervised setting, where the idea most directly applicable.
4.1 Noiseless Supervised Learning
Here we assume that all inputs are labeled, i.c. N = and all observations are noiseless, i.c. of = 0.
Denoting by Pxx Kxx the precision matrix, the log-likelihood log p(y 0,X), up to multiplicative constant and additive constant irrelevant for optimization, is given by
(()) = - log det(Kxx)
TK = log det det(Pxx) y Pxxxy.
(14)
Its gradient may be given and then subsequently approximated (Gardner al., 2018) by
aC( = tr P XX 00 00
pxxx
0Pxx T
0 2
xx
(15)
00
d0
00
where z is a random vector consisting of IID variables that are either or - with probability 1/2. The first term on the right-hand side is the stochastic estimate of the trace of Hutchinson (1989). Since the kernels from Section 3.3 coincide with graph Matérn kernels on the nodes Tú, we have
0² Kxx Cun
Σ (n)(fi
A Afr.
D. Spm
(16)
However, as the graph bandwidth is one of the hyperparameters optimize over, using Equa- tion (16) would entail repeated eigenpair computations and differentiating through this procedure. Because of this, we use an alternative way to compute matrix-vector products Pxx* detailed below.
Proposition 2. Assuming v E N, the precision matrix P Xxx of F., (x:; 2z)) can be given by
-2
2v
2v
Pxx
D
I A
(17)
rw
C
Proof. See Appendix A.
'Though automatic differentiability could in principle work for iterative methods like the Lanczos algorithm. the amount of memory required for storing the gradients of the intermediate steps quickly becomes prohibitive.
7

---

Using Proposition 2 to evaluate matrix- products Pxx u and conjugate gradients (Meurant, 2006) to solve z P using only the matrix-vector products, we can efficiently evaluate the right- hand side of Equation (15), with linear costs with respect to N, assuming that the graph is sparse. Preconditioning (Wenger et al., 2022) can be used to further improve the efficiency of the solve.
4.2 Noiseless Semi-Supervised Learning
Here we assume that inputs are partly unlabeled, i.e. N * n. while observations are still noiseless, ie. o² = Denote Z be the labeled part of X. Then the matrix Kxx in the log- ligellilooood given by Equation (14) should be substituted with Kzz. However, while Pxx can be represented using Equation (17), the precision Pzz = Kzz cannot. We thus compute it as the Schur complement:
A B Pxx = C D
Pzz = A BD,
(18)
where partitioning of Pxx corresponds to partitioning X into Z and the rest. Evaluating matrix- vector product P zzzu requires a solve of D D- DCC).. This solve can also be performed using conjugate gradients, keeping the computational complexity linear in N but increasing the constants.
4.3 Handling Noisy Observations
Finally, assume noisy observations, i.e. of O. The inputs be partially unlabeled, i.e. N of
In this case, matrix Kxx in the log-likelihood given by Equation (14) should be substituted with Kzz + L To reduce this to the previously considered cases, we use the Taylor expansion
(Kzz + ~ KZZ KzZ a K zz
Pzz a² P + o
(19)
In practice, we only use the first two terms on the right-hand side an approximation. This allows to retain linear computational complexity scaling with respect to N but increases the constants.
4.4 Resulting Algorithm
Here we provide concise summary of the implicit manifold Gaussian process regression algorithm.
Step 1: KNN-index. Construct the KNN index on the points 31,.., N. This allows linear time evaluation of any matrix-vector product with A, and thus also with A, A yucr Pxx for E N. etc.
Step 2: hyperparameter optimization. Find the hyperparameters i that solve Equation (13). Assuming V = & E N is manually fixed, this relies only on matrix-vector products with Arw-
Step 3: computing the eigenpairs. Fixing the graph bandwidth & found on Step 2, compute the eigenpairs At, fi corresponding to the L smallest eigenvalues \₁. For large N, use Lanczos algorithm.
After the steps above are finished, Equations (8) and (11) define the geometric kernel A,A,aa (2x, x') for arbitrary x, x IRd. Then the respective prior GP X-X (D, can be conditioned by the labeled data in the standard way, yielding the posterior f(m) ~ GP(m k((m).. To get sensible predictions far away from the data, the geometric model (m)) is convexly combined with an independently trained classical Gaussian process model, as given by Equation (12). The resulting predictive model is still a Gaussian process, sum of two appropriately weighted independent Gaussian processes.
The number of neighbors K, the number of eigenpairs L and the smoothness are assumed to be manually fixed parameters. Higher values of K and L improve the quality of approximation of the manifold kernel, which is often linked to better predictive performance, but requires more computational resources. The parameter v can be picked using cross validation or prior knowledge. Small integer values of v reduce computational costs, but may be inadequate for higher dimensions of the assumed manifold due to the v/ = v dim(M)/2 /2 link with the manifold kernel smoothness v/.
Remark.
5 Experiments
We start in Section 5.1 by examining simple synthetic example to gain intuition on how noise- sensitive the technique is. Then in Section 5.2 we consider real datasets, showing improvements in higher dimensions. More experiments, results, and additional discussion can be found in Appendix B
8

---

(a) Ground truth
(b) Noiseless, prediction
(c) Low noise, prediction (d) High noise, prediction
Figure 4: The ground truth function the dumbbell manifold and the predictions of the implicit
manifold Gaussian process regression (IMGP) under different levels noise.
5.1 Synthetic Examples
We consider one dimensional manifold resembling the shape of a dumbbell which already appeared in Figures and 3. The unknown function f. is defined by fixing a point x² in the top left part of the dumbbell, and computing sin (n(x"",. .)) where d(-,-) denotes the geodesic (intrinsic) distance between pair of points on the manifold. This function illustrated in Figure 4a.
To measure performance we primarily rely on measuring negative log-likelihood (NLL) on the dense mesh of test locations. We do this because such metric is able combine accuracy and calibration simultaneously. Additionally, we present the root mean square error (RMSE).
We investigate semi- supervisdd setting where the number of unlabeled points is large (N - 15466 and the number of labeled points is small (n = 10). We contaminate the inputs with noise, putting X noiseless + NO, x I) and do the same with the outputs, putting y f. (X) + N(0, og.1 for various values of xx By > O. Specifically, we refer to as the noiseless setting, the low noise setting and the high noise setting, respectively.
consider ox dy B E {0. 0.01, 0.05) to which
The results for these are visualized in Figures 4b to 4d with performance metrics reported in Table 1. The implicit manifold Gaussian process regression referred to as IMGP (we use 1) and is compared with the standard Euclidean Matérn-5/2 Gaussian process. IMGP performs much better in the noiseless and the low noise settings. The high noise is enough damage the calibration of IMGP, as it ties with the baseline model: NLL is slightly worse and RMSE is slightly better.
In Appendix B.1 show how performance depends on the fraction n. n/N of labeled data points, the truncation level L and we discuss the choice of the K parameter in KNN, Additionally, in Ap- pendix B.2 we consider sensiiveny for 2D manifold.
5.2 High Dimensional Datasets
For the nigh-dimensional setting, we considered predicting rotation angles for MNIST-based datasets. Additionally, we examined high-dimensional dataset from the UCI MI. Repository, CT slices.
5.2.1
Setup
Datasets. We consider two MNIST-based datasets. The first one is created by extracting a single image per digit from the complete MNIST dataset. By randomly rotating these 10 images we obtained N 10000 training samples and 1000 testing samples. We call it Single Rotated MNIST (SR-MNIST). For the second dataset, we select 100 random samples from MNIST. By randomly rotating these, we generate N 100 000 training samples, most will be unlabeled, and 000 testing samples. We call it Multiple Rotated MNIST (MR-MNIST). The last dataset, CT slices, has dimensionality of d = 385, we split have N 24075 training samples and 24075 testing samples. Dataset names can be complemented by the fraction of labeled samples, e.g. MR-MNIST--105 refers to n = 10%N.
RMSE
NLL
0.01
8 0.05
B A=
8 0.01
0.05
Ruclidean Matérn-F/2
0.98
0.99
1.02 + 0.03
-2.17
-2.09
0.03
0.10
IMGP
0.33
0.34 + 0.02
1.00 = 0.02
-5.02
-4.19 = 0.1
-1.91 1.88
Table 1: Performance metrics for the dumbbell manifold with varying magnitude of noise B

---

MNIST
CT slices
Method
SR 10%
MR 1%
MR 10%
5%
10%
25%
EGP
0.01
0.00
-0.33 0.01
-0,80 0.02
-0.96 0.00
-1,20
S-IMGP
-1.42 0.01
2:24 ::
-0.88 +: 0.08
0.47 +: 0.06
-0.39 4: 0.08
-0.08 0.01
SS-IMGP
-1.52
-0.59
-0.79
26.1
12.7
1.03
0.68
S-IMGP (full)
0.64 -: 0.83
0.88 = 0.29
-0.02
SS-IMGP (full)
-2.48 + 0.08
-.33 1: 0.04
-1.99 0.04
Table 2: Negative log likelihood on test samples for real datasets. For RMSE see Tables 4 and 5.
Methods. We consider implicit manifold Gaussian processes in the supervised regime (S-IMGP) and in the semi-supervised regime (SS-IMGP). We compare them to the GPyTorch implementation of the Euclidean Matérn-5/2 Gaussian Process. We refer to it as the Euclidean Gaussian Process (EGP).
Additional details. We run 100 iterations of hyperparameter optimization using Adam with fixed learning rate of 0.01. For MNIST, with use IMGP with v = 2; for CT slices with v =
5.2.2 Results
Table 2 shows the negative log- legllikllhoodd metric for different datasets and methods on the test set. The respective RMSEs are presented in Appendix B.3. On SR-MNIST, IMGP outperforms EGP in both supervised and semi-super scenarios. MR-MNIST is more challenging. In the supervised setting for n = 1%N, S-IMGP is incapable of inferring the underlying manifold structure, performing worse than EGP. However, SS-IMGP, with more data to infer manifold from, performs best. For n 10%N, IMGP gets better grip of the dataset's geometry, outperforming EGP in both regimes.
For CT slices, regardless of both S-IMGP and SS-IMGP performed poorly. Looking for an explanation, we considered two modifications. First, we fixed the graph bandwidth á found in the algorithm's Step 2 (ef. Section 4.4), and re-optimized the other hyperparameters K, o² oo by maximizing the likelihood the eigenpair- based model (truncated to L = 2000 eigenpairs) computed in the algorithm's 's Step 3. This resulted in limited improvement but did not change the big picture in fact, values for S-IMGP and SS-IMGP in Table 2 for CT slices already include this modification.
Second, on top of this hyperparameter re-optimization, we tried computing the eigenpairs using torch. brich.linalg.eff instead of the Lanczos implementation in GPyTorch, taking the same number L 2000 of eigenpairs. The resulting methods S-IMGP (full) and SS-IMGP (full) showed considerable improvement over the baseline, as shown Table 2. This indicated an issue with the quality of eigenpairs derived from the Lanczos method which requires further investigation. We discuss this in Appendix B.3, together with the aforementioned hyperparameter re-optimization procedure.
Conclusion
In this work, we propose the implicit manifold Gaussian process regression technique. It is able unlabeled data to improve predictions and uncertainty calibration by learning the implicit manifold upon which the data lies, being inspired by the convergence of graph Matérn Gaussian processes to their manifold counterparts. This helps building better probabilistic models in higher dimensional settings where the standard Euclidean Gaussian processes usually struggle. This is supported by our experiments in a synthetic fow-dimensional setting and for high-dimensional datasets. Leveraging sparse structure of graph Matérn precision matrices and efficient approximate KNN, the technique able to scale to large datasets of hundreds of thousands points, which is especially important in high dimension. where a large number of unlabeled points is often needed learn the implicit manifold. The model is fully differentiable, making it possible to infer hyperparameters in the usual way.
Limitations. The quality of the constructed graph significantly influences the technique's 's perfor- mance. When dealing with data from complex manifolds or exhibiting highly non-uniform density, simplistic KNN strategies might fail to capture the manifold structure due their reliance on single graph bandwidth. In such scenarios, larger values of parameters K and L, or in high dimensions, of parameter v, may be beneficial but could substantially increase computational costs. Furthermore, larger datasets coupled with high parameter values can lead to numerical stability issues, for instance, in the Lanezos algorithm, calling for further improvements and research. Despite these challenges, method shows promise for advancing probabilistic modeling in higher dimensions.
10

---

Acknowledgements
We express gratitude Alexander Shulzhenko (St. Petersburg University, mentored by VB), whose initial work on a similar problem and accessible source material at https://githm: com/ AlexanderShulzine wtanko/ wwweeeo/eemeeiciiieeeeeeeee implicit--e paper, sparked inspiration for this project. We thank Dr. Alexander Terenin (Comell University) for generously making his Blender rendering scripts accessible to the public. These scripts, which com/
while not included in the
aided us in creating Figure 2, can be found in his Ph.D. thesis repository at attps://githi aterenin/ iterenin/phdthesee hwwppphhs......eeee BF and AB acknowledge support by the European Research Council (ERC), Advanced Grant agreement No 741945, Skill Acquisition in Humans and Robots. VB acknowledges support by ETH Zürich Postdoctoral Fellowship.
References
I. Azangulov, A. Smolensky, A. Terenin, and V. Borovitskiy. Stationary Kernels and Gaussian Processes on Lie Groups and their Homogeneous Spaces I: the compact case. arXiv preprint arXiv:2208. 14960, 2022. Cited on pages 2, 3.
I. Azangulov, A. Smolensky, A. Terenin, and V. Borovitskiy. Stationary Kernels and Gaussian Processes on Lie Groups and their Homogeneous Spaces II: non-compact symmetric spaces, arXiv preprint arXiv:2301. 13088, 2023. Cited on page 2.
M. Binois and N. Wycoff. A survey on high digeesinnniional Gaussian process modeling with application to Bayesian optimization, ACM Transactions Evolutionary Learning and Optimization, 2(2):1- 26, 2022. Cited on page 2.
V. Borovitskiy, I. Azangulov, A. Terenin, P. Mostowsky, M. Deisenroth, and N. Durrande. Matém Gaussian Processes on Graphs. In International Conference on Artificial Intelligence and Statistics, 2021. Cited on pages 3.
V. Borovitskiy, M. R. Karimi, R. Somnath, and A. Krause. Isotropic Gaussian Processes Finite Spaces Graphs. In International Conference on Artificial Intelligence and Statistics, 2023. Cited on page 2.
V. Borovitskiy, A. Terenin, P. Mostowsky, and M. P. Deisenroth. Matérn Gaussian Processes on Riemannian Manifolds. In Advances in Neural Information Processing Systems, 2020. Cited on pages 2, 3, 14.
J. Calder and N. G. Trillos. Improved spectral convergence rates for graph Laplacians on e-graphs and k-NN graphs. Applied and Computational Harmonic Analysis, 60:123-177 2022. Cited on page 4.
J.-P. Chilès and P. Delfiner. Geostatistics: Modeling Spatial Uncertainty. John Wiley & Sons, 2012. Cited on page 1.
R. R. Coifman and S. Lafon. Diffusion Maps. Applied and Computational Harmonic Analysis Special Issue: Diffusion Maps and Wavelets, 21(1):5-30, 2006. Cited on pages 2, 4.
E. De Vito, N. Mücke, and L. Rosasco. Reproducing kernel Hilbert spaces on manifolds: Sobolev and diffusion spaces. Analysis and Applications, 19(03):363-396, 2021. Cited on page 14.
M. Deisenroth and C. E. Rasmussen. PILCO: model- based and data-efficient approach to policy search. In International Conference on Machine Learning, 2011. Cited on page 1.
H. Donnelly. Eigenfunctions of the Laplacian on compact Riemannian manifolds. Asian Journal of Mathematics, 10(1): 2006. Cited on page 14.
D. B. Dunson, H.-T. Wu, and N. Wu. Graph based Gaussian processes restricted domains. Journal of the Royal Statistical Society Series B: Statistical Methodology, 84(2):414 439, 2022, Cited on page 2.
D. B. Dunson, H.-T. Wu, and N. Wu. Spectral convergence of graph Laplacian and heat kernel reconstruction in 7,00 from random samples. Applied and Computational Harmonic Analysis, 55:282-336, 2021. Cited on pages 2, 4.
11

---

A. Feragen, F. Lauze, and S. Hauberg. Geodesic exponential kernels: When curvature and linearity conflict. In Conference on Computer Vision and Pattern Recognition, 2015. Cited on page 3.
N. García Trillos, M. Gerlach, M. Hein, and D. Slepëev. Error estimates for spectral convergence of the graph Laplacian on random geometric graphs toward the Laplace Beltrami operator. Foundations of Computational Mathematics, 20(4):827-88) 887, 2020. Cited on page 4.
Gardner, G. Pleiss, K. Q. Weinberger, D. Bindel, and A. G. Wilson. GPyTorch: Blackbox Matrix- Matrix Gaussian Process Inference with GPU Accelcration. In Advances in Neural Information Processing Systems, 2018. Cited on pages 18.
T. Gneiting, Strictly and non-strictly positive definite functions on spheres. Bernoulli, 19(4): 1327- 1349, 2013. Cited on page 3.
I. S. Gradshteyn and I. M. Ryzhik. Table of Integrals, Series, and Products. Academic Press, 7th edi- tion, 2014. Cited on page 3.
M. Hein, J.-Y. Audibert, and U. von Luxburg. Graph Laplacians and their Convergence on Random Neighborhood Graphs. Journal of Machine Learning Research, 8(48):1335-1336, 2007. Cited on page 4.
P. Hennig, M. A. Osborne, and M. Girolami. Probabilistic numeries and uncertainty in computa- tions. Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences, 471 47102179)20150142,, 2015. Cited on page 1.
M. F. Hutchinson. A stochastic estimator of the trace of the influence matrix for Laplacian smoothing splines. Communications in on page 7.
and Computation, 1853):005- 1989. Cited
J. Johnson, M. Douze, and H. Jégou. Billion-scale similarity search with GPUs. JEEE Transactions on Big Data, ((3):535-547) 2019. Cited on page 6.
R. I. Kondor and J. Lafferty. Diffusion kernels on graphs and other discrete structures. In International Conference Machine Learning, 2002. Cited on page 2.
J. A. Lee and M. Verleysen. Nonlinear Dimensionality Reduction. Springer Science & Business Media, 2007. Cited on page 4.
Lindgren, H. Rue, and 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 498, 2011. Cited on page 3.
Y. Ma and Y. Fu. Manifold Learning Theory and Applications. CRC press, 2011. Cited on page 4.
G. Meurant. The Lanezos and conjugate gradient algorithms: from theory finite precision computa- fions. SIAM, 2006. Cited on pages 7, 8.
F. Pavutnitskiy, S. O. Ivanov, E. Abramov, V. Borovitskiy, A. Klochkov, V. Vyalov, A. Zaikovskii, and A. Petiushko. Quadric hypersurface intersection for manifold learning in feature space. In International Conference on Artificial Intelligence and Statisties, 2022. Cited on page 4.
R. L. Peach, M. Vinao-Carl, N. Grossman, M. David, E. Mallas, D. Sharp. P. A. Malhotra, P. Vandergheynst, and A. Gosztolai. Implicit Gaussian process representation of vector fields over arbitrary latent manifolds. arXiv preprint arXiv:2309. 2023, Cited on page 2.
A. Rahimi and B. Recht. Random features for large-scale kernel machines. In Advances in Neural Information Processing Systems, 2007. Cited on page
C. E. Rasmussen and C. K. Williams. Gaussian Processes for Machine Learning. Adaptive Computation and Machine Learning. MIT Press, 2006. Cited on pages 1, 2.
J. Wenger, G. Pleiss, P. Hennig, J. Cunningham, and J. Gardner. Preconditioning for scalable Gaus- sian process hyperparameter optimization. In International Conference on Machine Learning, pages 23751-23788 PMLR, 2022. Cited on page 8
12

---

P. Whittle, Stochastic Processes in Several Dimensions. Bulletin of the International Statistical Institute, 40:974-999 1963. Cited on page 3.
J. T. Wilson, V. Borovitskiy, A. Terenin, P. Mostowsky, and M. P. Deisenroth. Efficiently Sampling Functions from Gaussian Process Posteriors. In International Conference on Machine Learning. 2020. Cited on page
T. Wilson, V. Borovitskiy, 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 2.
13

---

A Theory
Proposition 1. Denote the eigenpairs by At fi for a graph Laplacian and by X)4 Si for the Laplace- Laplace-Beltrami operator. Fix > O. Assume that, with probability at least 1 1- 5, for all E > for small enough and for N large enough we have 1 XM < E and fi(2;)) - (
< e. Then, with
probability at least 1 8, we have k2 (Di 2: - 2 (xx, R.)) as - 0, N, I. - co, where
0² C2,e
2v
Kynnoa (2),
thi
fi(z))ff(z3))
(6)
1=0
Proof. Fix small & > O. We will prove that for a small enough and N, L large enough we have Iku,a.o7 (r,, xz)) - N. (x; xz)l < Ce for some C > 0 with probability at least 1 1 8. Since the assumption holds on the same event of probability convergence on the same event. In fact, probabilistic narrative is nonessential for what actually prove, and we do not use it below. To simplify notation, we replace El- by Σi=
5 for all e, this directly translates to the
het First, for any L E Σzo define the truncated version Ne,.,03 of the manifold kernel ku,n,o3 by
L
di/(MM/2
k ke x')
Σ 1=0
2v R2
(M)
(20)
Cunu
The manifold Matérn kernels are the reproducing kernels of Sobolev spaces, if the latter are defined appropriately (Borovitskiy et al., 2020). These are Mercer kernels (De Vito et al., 2021), hence, by Mercer's theorem, het k2,n.o kymnoe uniformly on .M, i.e. for L large enough we have
| 2,) (um.3(101,
sup x.
Rh K (xx, (x, x') - Kum.o (2) ( <E
(21)
Now suppose that is small enough and N is large enough so that
A XM| < 6', (i(())) JiM (ax)) < <'', e' = min (1, (AL)) (2 2
(22)
T.
for all E {1,.., L} and for all i E {1,.., N} with probability at least 5.
Assuming the manifold is connected, by Donnelly (2006) we have ||/M| S Cx (A d-1 for >
where CA > 0 is constant that depends on the geometry of the manifold. The case = 0 is special
because X 0. Since M is a constant function, we have
IMM S C' (
(23)
for all 2 O, where Cs here is potentially different from the C3 before. Assuming, without loss of generality, 1, we have
[f(():))( (Ma)))
s [f(=a))](((((>>>M))))))))
(24)
+ + - -(
(25)
< +
(26)
S ( -
+ [)( + (())M(a;)))
(27)
(1 2CxJe
< (e' 20' max ((MM) $ 1)) <
(28)
L
The function ((A)= = (2x + +2)
n (M) is Lipschitz: (a( - +())) S Cola - XI where
-0-din(M//2-1
-v-din(Mi/
2v Co sup 1 ((()): sup (v + dim(M)/22
(2V) n3
> +2)
(v dim(MM/))
(29)
X
X
x2
14

---

Define an auxiliary kernel with manifold eigenvalues and graph eigenfunctions by
W Mon,o (2,, x²)
of Cv.
2v
+ xp1)
-u-dim(M)/2
fi(x))
(30)
x22
I=0
Then
Cux of
k ounnog(0: (z)) - k kem (0; m;}
5
Σ>
22
+ 1M))
-u-dim(M)/2
(1 2Cx)e
(31)
RZ
L
1=0
$ {0)(1 + 20x)c.
(32)
Also, noting that [f(((e;)| S (x;) - f((20)1 + 11 S + +sx max - ,1, write
Cu. Tie (ur.o k.no (2ts; mz){
5
Σ
((1))) +
(33)
o3
1=0
C A -[1[(())]](({(
(34)
1m0
2C ((')² C max ( 1))
(35)
1=0
$ 2C (1 C%)e.
(36)
Finally,
ku.n.oz (Di, (2), 20,) - Mono (2; 2;} 1
k M.m (2;; (2) 20,) - kL un (20 xz))
(37)
Rb ummoo (2 125) - Mum (ar (23)
(38)
K kimm.o3 (0) vz,) - kino (20 xz;)
(39)
of (+(0)(+20s>) + 2C (1 + C3))..
S E +
(40)
Onn
This proves the claim.
Proposition 2. Assuming E N. the precision matrix Pxx of x k. (I xz) can be given by
Pxx
D
2
A
kw
(
2y
+ A Iw)
(17)
C
times
Proof. The covariance matrix Kxx is given by
Kxx = O² L-1
Σ
+())fff
( =
2v k2 +1)
(41)
C 1=
where S: D 1/2 s and faym the orthonormal eigenvectors of the symmetric normalized
Laplacian Aaym Denote
T²
sym K
(
(42)
Xx
C
1=0
then (K*)) xx
L-1 2
X1)
f
paym ( T
2 A sym)
We have
C
|==
C
0² 1--1
Kxx Cv, 1=0
Σ ((x))
D 1/2 paym
D-1/22 1/2
D
sym Xx
D
1/2
(43)
15

---

900
1.0
IMGP
-2.0
IMGP
0.
100
0.5
EGP
0.6
100
200
0.
2.5
1000
0.4
044
-3.0
200
0.2
0.0
-3.5
510
10
140
(a) RMSE depending on L (b) NLL depending L (c) RMSE depending on ʃ (d) NLL depending on
Figure 5: Root Mean Square Error (RMSE) and Negative Log-Likelihood (NLL) for increasing number of eigenpairs L (left panels) and increasing fraction f n%,N of labeled points (right
panels). The legend in (a) and (b) refers to the number of hyperparameter optimization iterations.
On the other hand,
2v 2
A
D
2v
D
1/2
D1/2
(44)
1w
C
k2
sym
CN
D1/2
22
I
+
A
D1/2 = =
(
sym xx
D
(45)
Cu.
62
sym
KxX
Pxx.
(46)
B Additional Experimental Results and Details
Here we provide additional results and details for synthetic examples and high-dimensional datasets.
B.1 ID Dumbbell Manifold
In this section, analyze our method's 's sensitivity to the number of labeled points, spectrum truncation, and neighbor count in graph construction, offering deeper insights into its behavior.
Eigenpairs Truncation. From a theoretical standpoint, utilizing the complete set of eigenpairs to construct the kernel should be beneficial. However, this assumption may not hold true in cases where the optimization problem is not fully converged. The trends of RMSE and NLL, as depicted in Figures 5a and 5b, illustrate the impact of increasing the number of eigenpairs for 100, 200, and 1000 iterations. In situations where hyperparameter optimization does not converge fully, the length scale parameter x fails to reach sufficiently high values necessary to generate appropriate spectral density decay, which in turn would properly weigh higher frequency eigenpairs. In such scenarios, truncating the spectrum similar to increasing the length scale, which might improve results in certain cases. In the ID scenario, due to its simplicity, we opted for modest number of eigenpairs, namely L 50. Exceeding this count did not yield any discernible improvements.
Dataset Size. We analyze the sensitivity to dataset size of our method (IMGP) against the Euclidean case (EGP) in the semi- -upervueee learning scenario. For fixed number of eigenpairs, Figures 5c and 5d show performance metrics (RMSE and NLL. depending on the percentage n% of labeled points. In a typical scenario where we have at our disposal fewer labeled points compared to the number of unlabeled points, our method outperforms EGP in both accuracy (RMSE) and uncertainty quantification (NLL). As n increases EGP starts to match the performance of IMGP.
50,
Number of neighbors. For the one dimensional case considered in this section, only three neighbors are necessary to capture the essential features of the manifold's geometry. Choosing number less than three would hinder the algorithm's 's ability to capture the underlying manifold structure. On the other hand, increasing the number of neighbors beyond this threshold does not affect the solution, provided that sufficient time is allocated for hyperparameter optimization to converge and the graph bandwidth becomes small enough to "correct" for all undesired edges the graph structure (for instance in the central region of the dumbbell). In high-dimensional problems where the dimension of the underlying manifold is unknown, this suggests to incrementally increase the number of neighhors until the loss function stops improving, if it is computationally feasible try multiple values of K.
16

---

400
(a) Ground truth
(b) IMGP pred. B = (c)IMG pred. 8
(dJMGP pred. 8 0.05
(c) Euclidean, prediction
(f) Euclidean, uncertainty
(g) IMGP, uncertainty
Figure 6: (a) Ground truth function on the complex 2D manifold and (b)-(d) predictions of the implicit manifold Gaussian process regression (IMGP) increasing level of sampling noise B. (e)-(d) Euclidean GP prediction and uncertainty and (g) IMGP uncertainty in noiseless scenario.
B.2 2D Dragon Manifold
We consider another synthetic setting, a complex 2D manifold, depicted in Figure 6. The ground truth function is visualized Figure it is the same as the one-dimensional case, the sine the geodesic distance to a point, which is located in the green area.
In the sumi-super learning scenario, Figures 6e and 6b offer comparison of the posterior mean between the Euclidean GP and IMGP, while Figures 6f and 6g illustrate the posterior standard deviation for both models-quee quite similar each other, in this regime. Similar what did for the 1D compact manifold in Section 5.1, we evaluated the performance of IMGP under varying levels of sampling noise. Figures 6b 6d display the IMGP predictions in the semi-supervised learning scenario as sampling noise increases. Similar to the ID case, our approach consistently outperforms the standard Euclidean GP, as evidenced in Table
Remark. We observed that for higher levels of sampling noise, the linear combination of posteriors, as described by Equation (12), significantly outperform the single geometric model.
B.3 High Dimensional Datasets
Rotated MNIST. For IMGP, we use 2. As discussed in Appendix B. the optimal number of eigenpairs varies considerably depending on the convergence the optimization problem. Given the limited number of iterations per run we opted to fix the number of eigenpairs at 20%NN and 2%.N, for SR-MNIST and MR-MNIST, respectively. Table 4 reports the complete results obtained for the rotated MNIST dataset, including both RMSE and NL... Notably, these emphasize the importance of unlabeled points, S-IMGP performs worse than both SS-IMGP and EGP on MR-MNIST-1 1%.
RMSE
NLL
A=
0.01
= 0.05
9 = 0.01
B = 0.05
Fuclidean Matérn-6/2
0.24 0.02
0.24 0.01
0.26 = 0.00
-0.85 + 0.46
0.16 1.58
1.26 1.80
IMGP
0.12 111::00 0.01
0.21 0.01
0.01
0.14
0.00
-130 0.09
Table 3: Results for a complex 2D manifold with varying magnitude of sampling noise.
17

---

CT slices. For IMGP, we use v = 3. In Table 5 we report results for IMGP-S (full) and IMGP-SS (full). These are based on the "exact" eigenpairs, computed by torch.linalg.eige as opposed the standard Lanczos implementation we use by default. Additionally, these include the hyperparameter re-optimization step, as described in Section 5.2 and discussed below. Regarding RMSE, all three compared methods-IMGP (full) and IMGP-SS (full) and EGP- exG-bxhibit similar performance, with a slight advantage for SS-IMGP (full) as the number of training points increases. When considering NLL, as previously observed with MNIST, SS-IMGP performs best in all settings. However, NLL decreases for larger values of n/ n/N, indicating a probable overfit in this regime.
Hyperparameter re-optimization. We observed this step to serve two important purposes: (I) fixing overly small values of the signal variance o², potentially caused by the absence covariance normalization in the optimization process and poor convergence; (2) adjusting the length scale parameter to take into account the loss of high-frequency components due to truncation.
Implementation. Currently, the implementation faces two significant limitations that might restrict its usability to high-memory hardware setups. Firstly, due to the absence of rich enough sparse matrix routines in PyTorch, we had to develop own custom implementation of differentiable sparse operators. We kept to high-level routines, which forced us to strike a balance between performance and memory efficiency. In particular, for matrix-vector sparse product operations, our approach relies on highly optimized vectorized code, delivering high performance on GPU at the expense of increased memory allocation. Secondly, in the PyTorch ecosystem. Our attempts of using the PyTorch implementation of the LOBPCG algorithm, torch.lobpcg. yielded relatively poor results. We had similar experience with the Lanezos implementation from GPyTorch (Gardner et al., 2018). In light of this, when extracting higher- frequency components was needed, we resorted two alternative solutions: PyTorch dense matrix cigen-decomposition torch. winal..inagg.aao eiih and the SciPy Arpack wrapper scipy.linll..eige approach, while benefiting from GPU acceleration, can be infeasible because of limited GPU memory. The second approach, known for its efficiency in memory usage due to its Krylov subspace-based nature, is constrained to CPU utilization, significantly impacting the algorithm's 's overall performance.
faced challenges with sparse eigen-solvers
The first
C Hyperparameter Priors and Initialization
Here we describe hyperpameter priors which might be of help when using implicit manifold Gaussian process regression.
Graph Bandwidth a. Graph Laplacian converges to the Laplace- Beltrami operator when a tends to zero, motivating smaller However, in a non-asymptotic setting it is impractical to have overly small as it will render the graph effectively disconnected and cause numerical instabilities. One appropriate prior could thus be gamma distribution, whose right tail discourages high values of and which, given an appropriate choice of the parameters, encourages a to align with the scale of pairwise distances between graph nodes. Specifically, we choose the parameters of the gamma distribution so as to (1) match its mode with the median (Q2) of pairwise average distance between K-th nearest neighbors because such would give reasonable weights in the KNN graph and (2) so to place its standard 0.95 confidence interval to the right of a certain lower bound ā we define further that ensures the graph is numerically not disconnected. Let
D = max EKNNN,, - æ;i1 where j = L,..,N}
(47)
RMSE
NLL
Dataset
n//N
S-IMGP
SS-IMGP
EGP
S-IMGP
SS-IMGP
EGP
SR- MNIST
10%
0.04 - 0.01
0.01 + 0.00
0.12 + 0.01
-1.42 + 0.01
-1.52 + 0.01
-0.54 + 0.01
-NNIT
15
0.74 = 0.02
0.43 1 0.02
0.74 4 0.02
2.24 1 0.20
0.59 1 0.01
-0.20 1: 0.01
MR- -NNIST
10%
0.14
1 0.05
0.03 1 0.00
0.13 1 0.00
0.68 0.08
0.79 1 0.00
0.43 0.01
Table 4: Results for the rotated MNIST dataset.
18

---

RMSE
NLL
R/N
S-IMGP (full)
SS-IMGP (full)
EGP
S-IMGP (full)
SS-IMGP (full)
EGP
5%
0.27 + 0.02
0.39 = 0.04
0.24
0.64 0.83
0.08
0.02
10%
0.22 +
0.19 4: 0.01
0.16
0.88 0.29
-2.35
-0.96
25%
0.14 + 0.01
0.08 + 0.02
0.09 + 0.01
-0.42 + 0.10
-1.99
-1.20 0.09
50%
0.08 + 0.01
0.07 + 0.01
0.08 + 0.04
0.96 + 0.11
-2.04=0.02
-1.02 = 0.04
Table 5: Results for Relative location of CT slices on axial axis (d 385, N 48150) from UCI
Machine Learning Repository.
1000
200-
Hiureram percentle
Froquency Prcbbbalty
perhetee Lower Brand
Pooquenccy Peduulittiyy
750
130
at Lover Bouad Produivillt Carilty Carve
Prchobilly Teeeity Orrse
500
100
250
0.3
1.5
2.0
3.0
3.5
0-
Avaogge Distance
Avorage Distaras Valuce
(a) Prior for the -uniisuperviee case
(b) Prior for the supervised case
Figure 7: The histogram of 'D and the prior for the bandwidth hyperparameter a.
The bandwidth's 's lower bound we compute as
d2 4og T
min a
V
(48)
dED
where 7 user-defined parameter.
Let 7 and 33 the shape and rate parameters of the gamma distribution. order to achieve (1), we define
1 pQ2 +
and
B == p
(49)
where p is used achieve (2) and is given by
4Q2
(50)
(Q2 a)2
Considering the S-MNIST dataset, Figure 7 shows the bandwidth prior distribution for the semi- supervised (labeled and unlabeled points) and the supervised (labeled points) scenarios.
Signal Variance o3. Assuming normalized ye. the signal variance of should be close 1. Because
of this a natural prior for o² is the truncated normal (onto the set of positive reals a² > 0), with mode at 1. The pre-truncation variance away from the mode, i.e. it can be chosen to be 1/9. Note that setting a prior over the signal variance parameter requires evaluating the normalization constant C.,"" for the kernel at each hyperparameter optimization step, something that can be avoided otherwise.
be chosen, for example, have 0 at standard deviations
Noise Variance o². Choosing a prior for the noise variance o² is heavily problem-dependent. Truncated normal (onto the set o² > 0) with mode 0 could be reasonable option. The pre- truncation variance can be chosen, for example, to be 1, 1/4 or 1/9, placing the value 1, which is the variance of the normalized observations ys. at 1, 2 or 3 standard deviations away from the mode.
M "Covariance matrix normalization constant Cs. can be approximated as Co.e. TI Eis e. P Pxxe: ei. where M is relatively small, e; random standard basis vectors and the solve is performed by running the conjugate gradients. Differentiability the model preserved. Note however that we did not use this normalization in our tests since the performance improvement we observed was not enough to justify the additional computation overhead. This trade-off can vary considerably from case to case, which is why in our inplementation, the covariance normalization is optional.
19

---

Length scale The interpretation and the scale of the length scale parameter is manifold-specific. This makes very difficult to come up with any reasonable prior. Because of this, we suggest actually leaving the length scale parameter free.
Parameter initialization When doing MAP estimation, one can initialize parameters randomly, sampling them from respective priors.
20


Reviews:
[{"review_id": "ikPxgHK7Yz", "reviewer": "Reviewer_ejKa", "summary": "This work extends the reach of Matern Gaussian processes to additionally learning the implicit low-dimensional (unknown) manifold the data lives on - the existence of such a manifold is suggested by the manifold hypothesis. The theory draws from existing Laplacian Matern Gaussian processes (Borovitskiy et al., 2020). The model they propose is differentiable wrt all model and geometry hypers. It is able to scale to thousands of data points by leveraging the sparse structure of Mather precision matrices. They provide a way to extend predictions to the whole R^{d} space by reverting to a Euclidean GP away from the manifold.  The experimental evaluation shows support for the implicit manifold Gaussian process.", "questions": "- The final predictive distribution is a mixture of Gaussians - one is the canonical Euclidean posterior predictive with SE-ARD kernel and the other is the geometry aware kernel underlying the same canonical posterior predictive equations? \n \n- What is the order of the steps in training - the log marginal likelihood depends entails the precision matrix which depends on the eigenpairs, but you say in line 206 that they are computed after the hyperparameters are found, these are mentioned in line 216 - is the \\hat{\\kappa} the number of neighbours of the KNN graph?", "limitations": "I think authors should devote a paragraph to this. I don't see limitations discussed anywhere.", "rating": 6, "confidence": 4, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "- The end-to-end differentiability where the kernel hyperparameters along with those that parameterise the underlying geometry can be learnt simultaneously using a single objective. \n- Scalability is achieved as the precision matrix corresponding to the KNN implied graph is sparse.  \n- The semi-supervsied case is interesting where a large amount of unlabelled data is leveraged to infer the underlying geometry of the manifold through a weighted nearest neighbour  graph.\nAll these points contribute to the quality and significance of the work.", "weaknesses": "- The exposition could benefit from an algorithm style pseudocode for the training and prediction steps, it should start with the interpretation of the data as a graph, computation of the elements needed for the Matern kernel etc - as I am a bit confused about the order of the steps. \n\n- While the details are all there in section 3.2, I think there need to be separate sections for construction of the matern kernel on graph nodes, computing the dependencies ie. the eigenpairs of the graph Laplacian, computing the kernel on arbitrary vectors in R^{d}, computing the predictive posterior. \n\n- I dont understand this line - pls explain or rewrite. Basically, why dont you just say how do you compute dim(M)?\n\nEssentially, from the theoretical results of Section 3.1 we borrow a particular weighting scheme, aiming to cancel out the possibly non-uniform density, and a specific choice of the graph Laplacian (the random walk normalized one).", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "*\"The exposition could benefit from an algorithm style pseudocode ..\", \".. I am a bit confused about the order of the steps.\", \"What is the order of the steps in training\"*\n\n* Thank you for mentioning this.\nWe will emphasize the general flow of the algorithm.\nTo clarify, the algorithm can be summarized as follows.\n1. We compute KNN index using FAISS. This is enough to define matrix-vector products (MVPs) with the weighted adjacency matrix of the graph.\n2. We find the optimal hyperparameters $\\hat{\\mathbf{\\theta}}$ by maximizing the likelihood (14). Here we assume that $\\nu$ is a small integer and use Proposition 2 to efficiently compute MVPs with the precision matrix using MVPs with the adjacency.\n3. We compute a set of eigenpairs corresponding to the smallest eigenvalues of the Laplaican (with fixed hyperparameters $\\hat{\\mathbf{\\theta}}$) using the Lanczos algorithm.\n4. We define the kernel as\n    $$\n    k(\\mathbf{x}, \\mathbf{x}') = \\frac{\\hat{\\sigma_f}^2}{C_{\\hat{\\nu}, \\hat{\\kappa}}} \\sum_{l=1}^L \\Phi_{\\hat{\\nu}, \\hat{\\kappa}}(\\lambda_l) f_l(\\mathbf{x}) f_l(\\mathbf{x}'), \\qquad \\Phi_{\\nu, \\kappa}(\\lambda) = \\left(\\frac{2 \\nu}{\\kappa^2} + \\lambda\\right)^{-\\nu}\n    $$\n    where $f_l$ are the Laplacian eigenvectors extended to the whole $\\mathbb{R}^d$ via (11).\n5. We compute the posterior $f^{(m)}$ corresponding to the kernel $k$. We also perform Gaussian process regression (including hyperparameter tuning) with the Euclidean squared exponential kernel, we call the result $f^{(e)}$. The final predictive model is then given by their weighted average (12).\n\n\n*\"While the details are all there in section 3.2, I think there need to be separate sections for construction of the matern kernel on graph nodes, computing the dependencies ie. the eigenpairs of the graph Laplacian, computing the kernel on arbitrary vectors in $R^{d}$, computing the predictive posterior.\"*\n\n* Thank you for suggesting how we can make the presentation more clear.\nThis suggestion is clearly connected to the one above and to the questions some other referees have.\nWe will revisit the presentation in Sections 3 and 4 to make them better structured, emphasize specific algorithm steps and their order, as well as to distinguish our own contributions there from the previous work.\nWe believe it will be fairly easy to do this and this will not require major changes.\n\n*\"I dont understand this line - pls explain or rewrite.\"*\n\n* Thank you for spotting this.\nWe realize now that this statement is unclear, the term \"weighting scheme\" is nor very intuitive nor is it introduced in any way.\nWhat is meant here is that, following Coifman et al. (2006), we (1) use the random-walk normalized Laplacian, as opposed to the unnormalized one or the symmetric normalized one and (2) use the normalized adjacency matrix $\\mathbf{A} = \\tilde{\\mathbf{D}}^{-1} \\tilde{\\mathbf{A}} \\tilde{\\mathbf{D}}^{-1}$ to counteract the presence of non-uniform density---using this adjacency is what we meant by the \"weighting scheme\" because it defines graph weights.\n\n*\"how do you compute dim(M)?\"*\n\n* There is plenty of ways to estimate $\\operatorname{dim}(M)$, from classical (Levina and Bickel, 2004) to modern (Denti et al. 2022).\nWe do not estimate $\\operatorname{dim}(M)$ though, for multiple reasons.\nFirst and foremost, because this merely results in  a reparameterization of the Matérn family: the meaning of parameter $\\nu$ becomes different, but in any case it is an unknown parameter to be chosen somehow.\nSecond, the possible values of $\\nu$ will be usually restricted by the problem size and computational resources at hand because larger $\\nu$ implies higher costs.\nA reasonable approach is to choose $\\nu$ by grid search where the grid consists of small integer values $1, 2, 3, \\ldots$.\n\n*The final predictive distribution is a mixture of Gaussians - one is the canonical Euclidean posterior predictive with SE-ARD kernel and the other is the geometry aware kernel underlying the same canonical posterior predictive equations?*\n\n* Almost, with one clarification, the final predictive distribution is a *Gaussian distribution*, not a mixture of Gaussians.\nIts mean and variance are convex combinations of the mean and variance of the canonical Euclidean posterior and the posterior under the geometry aware kernel similar to what you suggested.\nAnother point to note is that in our experiments we did not use the ARD kernel, just the plain kernel with one length scale although nothing prevents you from using it.\n\n*\"the log marginal likelihood depends entails the precision matrix which depends on the eigenpairs, but you say in line 206 that they are computed after the hyperparameters are found, these are mentioned in line 216\"*\n\n* You are right that the precision matrix can be (approximately!) expressed in terms of the eigenpairs corresponding to the smallest eigenvalues of the Laplacian.\nHowever, thanks to Proposition 2, for integer values of $\\nu$ the precision matrix can also be represented (exactly!) as a polynomial of the Laplacian. This allows us to efficiently compute exact matrix-vector products with the precision without evaluating the eigenpairs on each optimization step, and to find the hyperparameters without inefficient differentiation through eigenpair-computing routines.\n\n*\"is the $\\hat{\\kappa}$ the number of neighbors of the KNN graph?\"*\n\n* No, $\\hat{\\kappa}$ is the length scale of the geometric kernel, after optimization.\n\n*\"I think authors should devote a paragraph to this. I don't see limitations discussed anywhere.\"*\n\n* Thank you for mentioning this. Please see the general response."}, {"review_id": "ebnCoF1xO9", "reviewer": "Reviewer_2ieq", "summary": "The authors propose a methodology to extend Matern processes on implicit manifolds, which are modeled by $K$-NN graphs, and the kernel relies on the set of eigenvalues/vectors of the associated graph-Laplacian. An approximation of the eigenfunctions based on the eigenvectors is provided, which together with an additional trick allows to extend the process in the vicinity of the implicit manifold. The hyperparameters are learned in supervised and semi-supervised scenarios via a differentiable objective, and the performance is demonstrated in the experiments.", "questions": "Q1. Lines 98-100: Is there any reference that plugging in Eq. 3 the graph Laplacian and the Gaussian noise gives a Matern process on the graph?\n\nQ2. Why the dimension of the implicit manifold (graph) is 0? I suppose that the underlying manifold has some dimensionality. I think that there are approaches relying on graphs that try to estimate this dimensionality.\n\nQ3. How sensitive is the method with respect to the graph construction? In particular, how the number of neighbors $K$ influences the result? Can this be estimated somehow e.g. by cross-validation?\n\nQ4. How the method behaves if there are more than one connected components in the dataset? \n\nQ5. Lines 144-150: The graph Laplacian converges to the Laplace-Beltrami (in terms of eigenvalues/vectors/functions) for $K$-NN graphs as long as the data distribution is uniform on the manifold? While for non-uniform distributions such a convergence result does not exist?\n\nQ6. Lines 177-179: I think that here the gist is unclear.\n\nQ7. Which Nyström method is used?\n\nQ8. Fig 3. I think it is interesting to see the actual eigenvector on the data (e.g. color the points) and how this correlates to the extended version.\n\nQ9. Eq. 11 the eigenvalue should be probably $\\lambda_l$ and not $\\lambda_j$?", "limitations": "The authors do not discuss the limitations of their work, and in particular, the construction of the graph which I think is critical for the performance of the model. The methodology that is propose does not have a direct negative societal impact.", "rating": 5, "confidence": 4, "soundness": 2, "presentation": 3, "contribution": 3, "strengths": "- The technical part of the paper seems correct and sensible, but I have not checked the derivations in details.\n- The idea is a direct and relatively simple extension of previous work.\n- The results mainly due to the synthetic experiment seem convincing.\n- The paper in general is well-written.", "weaknesses": "- The main drawback of graph-based methods is typically the construction of the graph and how well this recovers the actual structure of the underlying manifold. I suppose that the current method has the same limitation that is not clearly discussed/analyzed in the paper.\n- As the methodology applies to general manifolds, I would expect to see more empirical results. For example, some synthetic 3d surfaces and more complicated high-dimensional real-world datasets as the rotation of MNIST is closer to a synthetic experiment. \n", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "*\"the construction of the graph and how well this recovers the actual structure of the underlying manifold\"*\n\n* This is indeed a very appropriate valid concern.\nWe apologize for not explicitly discussing the limitations in the paper.\nWe sketched a paragraph on this (see general response) and aim to include it into the new version of the paper.\nKNN graph can fail to capture the geometry.\nBecause of the asymptotic convergence, it should not be a problem when there is a lot of (potentially unlabeled) data.\nThis is why we emphasize the semi-supervised case where large number of unlabeled data points can facilitate small data learning from the labeled ones.\nFurthermore, our framework, being fully differentiable, sets up the groundwork for more sophisticated techniques that might include *graph learning* as in, for example, Kazi et al (2022).\nIn this particular paper, however, we wanted to describe the simple model and lay groundwork for the possible extensions.\n\n*\"The paper would greatly benefit from additional empirical findings.\"*\n\n* In addressing the high-dimensional context, we employed the rotated MNIST dataset where we have good reasons to believe that the manifold hypothesis holds.\nOther high-dimensional datasets may lack an explicit enough underlying manifold structure.\nThe attached PDF (see general response) contains a table presenting preliminary results in the supervised scenario for a random dataset from the UCI ML Repository.\nRMSE performance with IMGP surpassed that of EGP for larger sample sizes.\nHowever, NLL exhibited a less favorable trend compared to our observations for rotated MNIST.\nThis warrants further investigation.\nIt is possible that the chosen dataset, despite the high-dimensional feature space, does not inherently embody a manifold structure or such a structure could be highly irregular.\nThis is supported by the need for a substantial number of eigenpairs in achieving satisfactory outcomes in contrast to rotated MNIST i.e. slow spectrum decay.\nWe will discuss this example or a similar additional example in the paper as well as an example for a manifold in 3D (results in this case are also present in the attached pdf) that we excluded before to save space.\n\n*\"Q1\"*\n\n* Yes, this is discussed in detail in Borovitskiy et al. (2021).\nThis reference was accidentally removed from lines 98-100 in one of the final versions of the draft.\nThank you very much for spotting this!\n\n*\"Q2\"*\n\n* There is plenty of ways to estimate $\\operatorname{dim}(M)$, from classical(Levina and Bickel, 2004) to modern (Denti et al. 2022).\nWe do not estimate $\\operatorname{dim}(M)$ though, for multiple reasons.\nFirst and foremost, because this merely results in  a reparameterization of the Matérn family: the meaning of parameter $\\nu$ becomes different, but in any case it is an unknown parameter to be chosen somehow.\nSecond, the possible values of $\\nu$ will be usually restricted by the problem size and computational resources at hand because larger $\\nu$ implies higher costs.\nA reasonable approach is to choose $\\nu$ by grid search where the grid consists of small integer values $1, 2, 3, \\ldots$.\n\n*\"Q3\"*\n\n* Larger values of $K$ should in principle---and usually in practice, as we observed---lead to better results: the weighting of the graph does the heavy-lifting of capturing the geometry.\nAs a rule of thumb you should use $K$ as large as you can afford.\n\n*\"Q4\"*\n\n* The KNN graph we are using is always connected. \nWe further use such bandwidth prior that prevents weights from being overly small.\nTo consider datasets with well separated connected components, one can consider clustering them and then applying our technique to each of the clusters separately.\nThis would correspond to the geometric assumption of having a disconnected manifold.\nThis would, however, often make little sense because the aforementioned implies zero correlation between points on different connected components, i.e. this would correspond to modeling the connected components independently.\nSome datasets may posses a hierarchical structure: have many (small) disconnected components which, as whole components, covary with one another.\nThis might require some modified geometric assumptions (and model) to handle.\n\n*\"Q5\"*\n\n* We expect that convergence still holds for KNN graphs even in the case of a non-uniform density (when using Coifman's normalized adjacency matrix $\\mathbf{A} = \\tilde{\\mathbf{D}}^{-1} \\tilde{\\mathbf{A}} \\tilde{\\mathbf{D}}^{-1}$), but such a result cannot be found in the literature, to the best of our knowledge.\nWhen the distribution is uniform and using the usual adjacency matrix ($\\mathbf{A} = \\tilde{\\mathbf{A}}$) the convergence for KNN graphs is studied, for example, in Calder and Trillos (2022).\n\n*\"Q6\"*\n\n* Thank you for spotting this.\nWe realize now that this statement is unclear, the term \"weighting scheme\" is nor very intuitive nor is it introduced in any way.\nWhat is meant here is that, following Coifman et al. (2006), we (1) use the random-walk normalized Laplacian, as opposed to the unnormalized one or the symmetric normalized one and (2) use the normalized adjacency matrix $\\mathbf{A} = \\tilde{\\mathbf{D}}^{-1} \\tilde{\\mathbf{A}} \\tilde{\\mathbf{D}}^{-1}$ to counteract the presence of non-uniform density---using this adjacency is what we meant by the \"weighting scheme\" because it defines graph weights.\n\n*\"Q7\"*\n\n* Like, for example, in Section 3 of https://proceedings.neurips.cc/paper_files/paper/2003/file/cf05968255451bdefe3c5bc64d550517-Paper.pdf.\n\n*\"Q8\"*\n\n* Thank you for this suggestion. We did not include such a picture because of continuity: the values of the eigenvector match the values of the extension close to the manifold. We are happy to add the picture you requested to the appendix.\n\n*\"Q9\"*\n\n* This is a typo indeed, thank you very much for spotting this!\n\n*\"The authors do not discuss the limitations of their work\"*\n\n* Thank you for mentioning this. Please see the general response."}, {"review_id": "4dKZ61M9by", "reviewer": "Reviewer_g5fB", "summary": "Authors propose a novel methodology for doing GP regression which is able to learn the implicit structure from the data. This is particularly useful in high-dimensional problems, where the data  lies on low-dimensional manifold. The proposed methods allows to learn the implicit manifold in a fully differentiable way, using nearest neighbour graph.\n\nThe method is based on Matern GPs on manifolds and graphs (Borovitskiy, 2020, 2021). To approximate the eigenvalues of the Laplace-Beltrami operator on the implicit data manifold, authors use random walk normalized graph Laplacian (on KNN graph) weighted to overcome possible non-uniform sampling density.\n\nThis work leverages efficient approximation to the KNN, sparse structure of the precision matrices, and RFF kernel approximation to be able to scale to large datasets.\n\nThe model is tested on a synthetic dumbbell-shape manifold and on rotated MNIST data.", "questions": "1. Could you please clarify on the following: There are no results on spectral convergence for KNN graphs independent of data sampling density. You use a weighing scheme to deal with non-uniform data sampling. You do not claim that it implies spectral convergence. How much of a problem is it, both theoretically and practically?\n\nSuggestions:\n1. In figure 2, perhaps use different colours foe kernel and sample plots.\n2. eq(11) I believe $\\lambda_i$ should be $\\lambda_l$?", "limitations": "The limitations are adequately addressed.", "rating": 7, "confidence": 4, "soundness": 4, "presentation": 3, "contribution": 3, "strengths": "This is, to the best of my knowledge, a novel methodology for GP regression. This work is build on solid foundation and the method is derived naturally from Matern GPs on manifolds and graphs. Authors cleverly use a variety of techniques to keep good scalability. \n\nIn my opinion, the main contribution of this work is in the elegant assembly of the right methods, which gives a practically usable model for a well-known problem.\n\nThe manuscript is written clearly and is easy to follow.", "weaknesses": "Due to the incremental nature of the paper, it is not always clear what is a novel contribution of this work and what is a background from previous works (example, last paragraph of 2.1). \n\nThe experiments are somewhat limited. The method is tested on one synthetic dataset and a few versions of rotated MNIST. In the MNIST dataset, however, the added rotation must be inducing a lot of structure. It would be beneficial to showcase the method on a real dateset where an implicit manifold is already existing in the data.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "*\"it is not always clear what is a novel contribution of this work and what is a background from previous works (example, last paragraph of 2.1).\"*\n\n* Thank you for mentioning this point.\nIn light of this and some other comments we intend to revisit the presentation in Sections 3 and 4 to make them better structured, emphasize specific algorithm steps and their order, as well as to distinguish our own contributions there from the previous work.\nWe believe it will be fairly easy to do this and this will not require major changes.\n\n\"*The experiments are somewhat limited. The method is tested on one synthetic dataset and a few versions of rotated MNIST. In the MNIST dataset, however, the added rotation must be inducing a lot of structure. It would be beneficial to showcase the method on a real dateset where an implicit manifold is already existing in the data.*\"\n\n* In addressing the high-dimensional context, we employed the rotated MNIST dataset where we have good reasons to believe that the manifold hypothesis holds. Other high-dimensional datasets may lack an explicit enough underlying manifold structure. The attached PDF (see general response) contains a table presenting preliminary results in the supervised scenario for a random dataset from the UCI Machine Learning Repository. RMSE performance with IMGP surpassed that of EGP for larger sample sizes. However, NLL exhibited a less favorable trend compared to our observations for rotated MNIST. This warrants further investigation. It is possible that the chosen dataset, despite the high-dimensional feature space, does not inherently embody a manifold structure or such a structure could be highly irregular. This is supported by the need for a substantial number of eigenpairs in achieving satisfactory outcomes in contrast to rotated MNIST i.e. slow spectrum decay. We will discuss this example or a similar additional example in the paper as well as an example for a manifold in 3D---the results for the latter were actually excluded from the draft to save space. We think it might make sense to apply the technique in the contexts like modeling 3d structures of molecules where there is inherent symmetry (invariance or equivariance to rotations and translations), as an interesting direction for further work.\n\n*\"Could you please clarify on the following: There are no results on spectral convergence for KNN graphs independent of data sampling density. You use a weighing scheme to deal with non-uniform data sampling. You do not claim that it implies spectral convergence. How much of a problem is it, both theoretically and practically?\"*\n\n* We hypothesize that convergence should still hold true.\nThere is apparently no proof of this fact in the literature though.\nIt does not seem to be out of reach for the proof techniques used for similar problems, but such a proof would likely turn out to be long and heavy due to multiple \"moving parts\".\n\n* Practically, we do not expect this to be a problem either.\nAlthough we use convergence results to motivate the technique, the graph construction itself reflects the geometry of a point cloud in a rather intuitive way.\nWhat is more, it easy to imagine that the graph is capable of representing structures beyond manifolds in the rigorous mathematical sense, e.g. manifolds of different dimensions glued together (like presented on Figures 2--4 in Dunson et al. (2022).\n\n*\"In figure 2, perhaps use different colours foe kernel and sample plots.\"*\n\n* Thank you, we agree this should improve clarity.\nWe will make this change for the camera ready version.\n\n*\"eq(11) I believe  should be ?\"*\n\n* This is an actual typo.\nThank you very much for reporting this!"}, {"review_id": "qrOyi3MdXL", "reviewer": "Reviewer_DAU4", "summary": "The curse of dimensionality demonstrates in Gaussian process models in that they’re default kernel choices that often depend on the Euclidean distance between 2 points. Euclidean distance is a poor metric especially in high-dimensional settings, where we expect data points to lie on an, often unknown, manifold. This paper combines findings from multiple precursor works in assembling a tool set to for Gaussian process regression with a Metern kernel, while implicitly inferring the underlying manifold structure. ", "questions": "None", "limitations": "- I believe the presentation of the paper could be improved by a clear separation of predecessor works and original contribution: for example, Section 3.1 is almost solely composed of existing work except for the theoretical contribution of the Matern kernel convergence — an overall marginal contribution in the paper. \n- Table 3: “RMSE” in the table seems confusing for prediction of rotating angles. Is it so that the author is referring to the shortest arc as a metric for distance in rotations? \n- I believe that the paper could benefit from a illustration of what implicit manifold is learned from the rotated MNIST data: for example, a display of what the kernel values matrix looks like along the trajectory of a rotation compared to the Euclidean kernel matrix, or something else showing the geodesics of the points on the rotation trajectory. ", "rating": 7, "confidence": 4, "soundness": 4, "presentation": 3, "contribution": 3, "strengths": "I believe that while the paper mostly combines existing theoretical results in manifold learning and numerical analysis, it presents a compelling, novel framework of tackling manifold learning with Gaussian process regression. \n\nThe paper is largely self-sufficient: it tackles the complex problem of learning manifolds in the context of Gaussian process by proposing solutions to all facets of the problem: characterizing manifolds with graph Laplacian, gradient-based learning and scalability. ", "weaknesses": "I think the paper lacks a presentation of _how_ learning the manifold structure aids in the predictive performance of the rotated MNIST dataset. While one can clearly see _why_ the learning of an implicit manifold helps in predicting rotation angles for handwritten digit, as rotation expressly traverses around a nontrivial manifold, the benefits of this work go beyond simple improvements in prediction, but interpretability. While the authors did a very good job in demonstrating how the model correctly inferred the dumbbell-shaped manifold, I believe similar illustrations can also be done in the MNIST rotations. \n\n", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "*\"I believe the presentation of the paper could be improved by a clear separation of predecessor works and original contribution: for example, Section 3.1 is almost solely composed of existing work except for the theoretical contribution of the Matern kernel convergence — an overall marginal contribution in the paper.\"*\n\n* Thank you for mentioning this point.\nIn light of this and some other comments we intend to revisit the presentation in Sections 3 and 4 to make them better structured, emphasize specific algorithm steps and their order, as well as to distinguish our own contributions there from the previous work.\nWe believe it will be fairly easy to do this and this will not require major changes.\n\n\"*Table 3: “RMSE” in the table seems confusing for prediction of rotating angles. Is it so that the author is referring to the shortest arc as a metric for distance in rotations.*\"\n\n* This is a good catch! We actually do use RMSE here although arc length is more appropriate.\nThis, however, does not affect the results because the angle of rotation in the dataset was limited to $\\pm 45$ degrees.\n\n\"*I believe that the paper could benefit from a illustration of what implicit manifold is learned from the rotated MNIST data: for example, a display of what the kernel values matrix looks like along the trajectory of a rotation compared to the Euclidean kernel matrix, or something else showing the geodesics of the points on the rotation trajectory.*\"\n\n* Thank you for suggesting the ideas for visualizing the learned manifold.\nSince our method uses the eigenfunctions of the Laplacian as features, it is closely related to the *Laplacian eigenmaps* and the *diffusion maps* dimensionality reduction techniques.\nUsing one of these techniques is most likely the best way of visualizing the manifold we are learning.\nActually, a manifold of rotated images is visualized right in the original paper on diffusion maps: see Figure 2 in Coifman et al. (2006).\nNote though that such a visualization only relies on the first two eigenvectors of the Laplacian---the number of vectors determines the dimension of the \"picture\" and thus should be no larger than 3---and our kernel uses much more of them.\nWe will mention this in the paper."}]

Author rebuttal:
We thank the referees for their valuable summaries and insights which will help us improve the paper.
We have replied to every review with detailed comments.

An essential suggestion made by many referees was to add an explicit discussion of limitations.
One major limitation is that there is a trade-off between the complexity of geometry and non-uniformity of sampling data on one side and the number of unlabeled data points needed to capture meaningful structure on the other side.
The model works well when the former is not too pronounced and the latter is large.
Further, the model could struggle for clustered data.
Possible avenues for improving on the mentioned drawbacks include learning a graph in a more sophisticated way (e.g. adaptive bandwidth, differentiable graph learning) or clustering and combining models.
Another limitation is that using larger values of $\nu$ and $K$ result that may be desirable sometimes imply higher costs and are not always feasible.
We fully agree that including this explicitly would improve paper quality and aim to include it in the camera-ready version.

Meta-review: Four knowledgeable referees support accept, and I accept.
Decision: Accept (poster)