Graph Random Features for Scalable Gaussian Processes
Authors: Matthew Zhang
Jihao Andreas Lin
Krzysztof Marcin Choromanski
Adrian Weller
Richard E. Turner
Isaac Reid
Venue: ICLR 2026 Poster (2026)
Keywords: kernels
graphs
Gaussian processes
Monte Carlo
inference
Primary area: probabilistic methods (Bayesian methods, variational inference, sampling, UQ, etc.)
TL;DR: Approximating graph node kernels using random walks unlocks efficient GPs on graphs.

Abstract: We study the application of graph random features (GRFs) – a recently-introduced stochastic estimator of graph node kernels – to scalable Gaussian processes on discrete input spaces. We prove that (under mild assumptions) Bayesian inference with GRFs enjoys $\mathcal{O}(N^{3/2})$ time complexity with respect to the number of nodes $N$, with probabilistic accuracy guarantees. In contrast, exact kernels generally incur $\mathcal{O}(N^{3})$. Wall-clock speedups and memory savings unlock Bayesian optimisation with over 1M graph nodes on a single computer chip, whilst preserving competitive performance.

Published as a conference paper at ICLR 2026
GRAPH RANDOM FEATURES FOR SCALABLE GAUSSIAN PROCESSES
Matthew Zhang Jihao Andreas Lin 12 Krzysztof Choromanski
Adrian Weller Richard E. Turner' Isaac Reid 13
'Uiiversity of Cambridge "Max Planck Institute for Intelligent Systems 'Google DeepMind 'Alan Turing Institute
matthew 2hang4738gmail com ir337@cam.ac..
ABSTRACT
We study the application of graph random features (GRFs) rocently-in- troduced stochastic estimator of graph node kernels to scalable Gaussian processes on discrete input spaces. We prove that (under mild assumptions) Bayesian inference with GRBs enjoys O (N3/22 time complexity with respect to the number of nodes N, with probabilistic accuracy guarantees. In contrast, exact kernels generally incur ((N3). Wall-clock speedups and memory savings unlock Bayesian optimisation with over 1M graph nodes on a single computer chip, whilst preserving competitive performance.
1 INTRODUCTION AND RELATED WORK
Gaussian processes (GPs) provide powerful framework for learning unknown functions in the presence of uncertainty (Rasmussen, 2003). In certain applications, kernels based on Euclidean distance may be unsuitable: for example, when modelling traffic congestion, since pairs of locations that are spatially close may not be connected by roads. In this case, kernels defined the nodes of graph may be more appropriate (Borovitskiy et al., 2021; Smola and Kondor, 2003). One can then perform inference and make principled predictions, including during Bayesian optimisation, using GPs on graphs.'
Scalability of GPs on graphs. Like their Euclidean cousins, exact GPs on graphs incur (NN)) time complexity with respect to the number of nodes N. This makes them impractical when working with very large graphs. To mitigate this, practitioners use techniques such as 'graph Fourier features', which approximate the kernel matrix with a truncated eigenvalue expansion, or specifie sparse kernel families (Borovitskly et al., 2021). The former loses high-frequency kernel information and the latter limits flexibility. Alternatively, one can use kernels for small, local subgraphs, at the cost of no longer performing inference on the whole of the graph 9 (Wan et al., 2023).
Graph random features. this paper, we propose to instead use the recently-intro duced class of gruph random features (GRFs) sparse, unbiased estimates of graph node kernels computed using random walks (Choromanski, 2023; Reid et al., 2023). GRFs are Monte Carlo estimators of power series of weighted adjacency matrices, analogous to Von Neumann's 's celebrated Russian Roulette estimator (Carter and Cashwell, 1975: Hendricks and Booth, 2006). GRFss enjoy strong concentration properties (Reid et al., 2024b). They are able to estimate ficxible class of graph node kernels including the popular diffusion and Matérn kernels by varying the so-called 'modulation function..
GRFs for scalable GPs. Reid et al. (2024a) previously suggested using GRFs for GPs, as part of a broader study of variance reduction techniques. However, their experiments focused exclusively on the diffusion kernel with small graphs, failing exploit the estimator's sparsity to accelcrate inference. Morcover, they limited their (chiefly theoretical) study to computing the posterior, omitting exploration of applications such as Bayesian optimisation.
We consider GPs defined on the nodes a fined graph. The input space is finite and we perform inference for finite set of random variables, one per node. The relationships between these variables determined by the structure of via graph node kernel. Whilst some might prefer to call this a Gaussiaa random field' or simply a 'multivariate Gaussian', in this paper we use 'GP on graph' for consistency with recent literature (Borovitskly et al., 2021; Reid et al., 2024a; Wan et al,, 2023).

---

Published as a conference paper at ICLR 2026
Random Walks
Graph Random Features
Graph GPs
((1)
(
Kg
p ar W
K
E
g
g
Figure 1: GRFs for scalable GPs graphs. The GRF algorithm constructs random feature (i) for node (I,..., using random walks. K = [ ( +=== is sparse approximation of the kernel matrix K, enabling efficient posterior inference in (N
Key contributions. We investigate graph random features (GRFs) for Gaussian processes (GPs), unlocking scalable Bayesian inference on graphs with > 1M nodes. See Figure 1.
We GRFs to construct sparse estimates of learnable graph node kernels, and these as covariance functions for GPs.
2. We prove that Bayesian inference with GRFs enjoys O(N3/2) time complexity with prob abilistic guarantees on approximation quality, compared to ((N3) for exact alternatives. In experiments, this translates to 50 wall-clock speedups on graphs with fewer than 10K nodes. Remarkably, the fiexibility of GRFs sometimes enables them to outperform dense alternatives on Lest negative log probability density and root mean squared error.
We showcase our new techniques by performing Bayesian optimisation on massive graphs, implementing Thompson sampling with > 1M nodes single computer chip."
PRELIMINARIES
Consider an undirected graph 9= = (V,& W), consisting nodes V= {1, N). edges & C V V, and weighted adjacency matrix W R Here Wyy = 0 if (i. $ E. Define the graph Laplacian L = D W, with = diag Σ D Wyy). The normalised graph Laplacian : D-
whose spectrum lies in [0.2] (Chung, 1997).
Graph node kernels. graph node kernel k: 12 V R, mapping pairs of graph nodes to real numbers (Smola and Kondor, 2003). Heuristically, it assigns similarity scores to every pair of nodes, which are assembled into a Gram matrir K [b(i,,)) E NNN Many popular graph node kernels are parame- terised as functions of W, or L, expressed using the power series
symmetric, positive semidefinite function
K.(W) 2 W", R. Vr f0,1,..,,
(1)
The coefficients (a determine the behaviour of the kernel, e.g. whether it upweights
long- or short-range interactions. For instance, the graph diffusion kernel Kar == exp(-3L) takes a, (-8)"/r!, Graph node kernels flexibly capture structural information about 9. providing a natural choice for the GP covariance (Borovitskiy et al., 2021; Reid et al., 2024a).
"Specifically, NVIDIA GeForee RTX 2080 Ti GPU (11 GB memory).
2

---

Published as a conference paper at ICLR 2026
Gaussian processes. Let us now consider modelling functions h: V + R defined on the graph nodes. A common task is to identify the node that maximises h, egg. the most influential social media user, 'patient zcro' in an epidemiological contact network. We may wish to solve x² = arg zax h(z). Suppose we have access to a sequence of T noisy observations of the objective Y/ = h(zz) + E, ~ (V. (0,) & E 21, TT) at distinct nodes T V. on is the noise variance. A common choice of statistical surrogate for the objective function 4t, for which we can perform analytic Bayesian inference given observations Dj- = {(I,, ye 8e) } (also denoted 2 =xx, y}), is the Gaussian process (GP) (Rasmussen, 2003):
1=1
h(x) (()((()),(,,a"")))
(2)
Here, m(x) is the mean function and k(x, the covariance function ('kernel'). In our setting, the imput domain consists of the nodes of a fixed graph. We can 'train' the kernel parameters (==) by maximising the log-marginal likelihood on the training data 2r, and then compute the analytic posterior mean and covariance:
(x)()) = m(x) + k(x, x)) (k x.) + + (g - mm()))
(3)
k,( = k(x, x²) k(r, a) (K(,, + "(,,
(4)
Bayesian oplimisation (BO) (Jones et al., 1998; Močkus, 1974) uses the posterior mean and covariance or related quantities, like samples from the posterior- - to efficiently locate x" in the presence of uncertainty. BO trades off exploration and exploitation in a mathematically principled manner, helping us decide which nodes to query in our attempt to maximise h.
Efficiency and scalability. A core computational challenge with performing Bayesian inference graphs using GPs is that even just evaluating dense graph kernel K, generally incurs (NV time complexity, let alone computing the matrix-vector products and matrix inversions that we will in general require for BO. This is because it involves computing functions like exp(-) or (.)-1 of the N x N weighted adjacency matrix W, which becomes expensive for big graphs. For Euclidean kernels, a common recourse to improve scalability is to use random features (Rahimi and Recht, 2007; Yang et al., 2014): stochastic, finite- dimensional features { E R" whose dot product is equal to the kernel evaluation in expectation, k(i,j) =[(((i) ((j)). In close analogy for discrete domains, researchers recently introduced graph random features (GRFs) (Choromanski, 2023; Reid et al., 2023) sparse random walk-based vectors for unbiased estimation of graph node kernels.
Graph random features: sparse, sharp kernel estimators using random walks. The mathematical details of GRFs are involved and can be safely omitted on a first reading, but their behaviour can be intuitively understood as follows. Consider sequence of scalars (f) satisfying Eg fifr = a, V r, namely, the 'deconvolution' of (a,), We refer to ft as the modulation function. Suppose that the power series v = E f.W converges. Then it is straightforward to see that for symmetric W (undirected graphs), we have vp v = K2, Powers of an adjacency matrix count walks on a graph: for instance, W gives the (weighted) number of walks of length I between nodes and j. Since K, converges, longer walks must uvlnt be discounted, either due to decaying fi or due to multiplication of edge weights that are less than 1. The key insight of GRFs is that we can compute a Monte Carlo estimate & E N that satisfies v E( by importance sampling random walks.
Concretely, we simulate random walks out of every node of the graph. Each random walk of length L consists of number of 'prefix subwalks' namely, for each step i < L, the sequence of the first I nodes visited. We keep track of 1) the weights of edges they traverse, 2) the modulation function f, and 3) thoir marginal probabilities. Using a simple formula, can construct unbiased," sparse N-dimensional vectors that satisly B(p(i)) v()))
'It has been noted that the shared source of randommess actually introduces a (1/nn) bias term for estimates of diagonal kernel entries [K,lyy' This is of little significance for large graphs with many walkers so, following convention (Choromanski, 2023; Reid et al., 2023), we omit further discussion. One could remove this bias by sampling two independent ensembles of random walks and taking K = d d at the cost of losing the positive definiteness guarantee and thus (typically) worse performance.
3

---

Published as a conference paper at ICLR 2026
E [[ K Alg. below provides pscudocode. is deliberately kept high-level
for compactness; the interested reader find more details in App. A.
Algorithm 1: Constructing GRF vector yp(1) E IRN to approximate (((W)
Inputs: Graph 9, modulation function f + IR, random walk sampler p.
2 Output: Set of sparse GRFs ((p(1) RM that satisfy [Kale, = [B(x(()Tx(1)))
3 for E V:
4 initialise ((t) <- 0
5 for walker_idn idx E I,.., n!:
sample random_ walk
for prefix_ subfaik E random_wal walk:
(pi)lprefix srrrefiixssuxaaa([I]]] + (II traversed _eageresddddede__eeeeeeeee _eee__dgeeeeeeeeeee seie__eeehts)
8
[[]ength(preffx___ cenpgsttpprefiisuusabllaae subai/sssbualkk)//ppttftttt (tcpprrrffixssuaakk subwalk)
9
normalise (p(1)/ =n
Remarkably, under mild assumptions on 9 and (a,))== GRFs provide very sharp estimates of K2, In particular, the estimates satisfy exponential concentration bounds, whilst storing only ((1) nonzero entries per feature. See Theorem 1 for formal statement. As will see in Section 3, we can use the sparse kernel estimate K == $ as an efficient alternative to the dense exact kernel K speeding up inference from (NN) (NN/)).
3 SCALABLE POSTERIOR INFERENCE WITH GRFs
Next, we demonstrate how GRFs speed up inference. We begin by proving novel theoretical results (Section 3.1), and then describe our full efficient GP workflow (Section 3.2).
3.1 NOVEL THEORETICAL RESULTS
We first recall the following result for GRFs, proved by Reid et al. (2024b).
Thcorem 1. (GRFs sparse and give sharp kernel estimates (Reid et al., 2024b)). Consider a graph 9 with weighted adjacency matrix W and node degrees {d) s=1" Suppose we sample GRFs { by sampling random walks that terminate with probability p at each timestep, with modulation function f. Suppose also that c = E (ma. W.dd/(1-p) is finite. Then have that
12
P((0(1) ((j)) - [Kllg >1 < exx
(5)
2(2n 1)2c²
Moreover, with probability at least 1- any GRF vii) is guaranteed to be sparse, with at most n log(1 - (1 - 8)m
(027/))og) )log(1 p) nonzero entries.
Proof. The proof, based on McDiarmid's inequality, is reported by Reid et al. (2024b).
Theorem 1 demonstrates that, despite being sparse, GRPs give sharp kernel estimates. In particular, we can use Eq. (5) compute the number of walkers n needed to guarantee an accurate estimate of K., with high probability. Because of the bound, this number is independent of the graph size N. n then determines the number of nonzero entries in the GRF, which also inherits independence of graph size N. We note that Theorem makes the assumption about the graph S that the constant c is finite. This is not controversial; Reid et al. (2024b) provide extensive discussion. Intuitively, it natural that the spectrum of W must lie in some radius of convergence in order for the power series E a, W to converge. The condition for its Monte Carlo estimate to converge is only slightly stronger.
00
4

---

Published as a conference paper at ICLR 2026
For computational reasons we often only sample random walks up to some fixed maximum length mmax: e.g. a fraction of the graph diameter, whereupon fi = 0 [ >> mmxx (discussed in App. C.1). The condition thus trivially holds in any reasonable implementation.
Given Theorem 1, we will henceforth assume that the number of walkers n is constant, confident that this gives sharp kernel estimate. Property (2) of Theorem 2 is novel.
Theorem 2/ (Properties of K). The randomised approximate Gram matrix K == = NNNN has the following properties.
(e(1) p() 1,1=1
1. Property 1. K supports O(N) matrix-vector multiplication;
2. Property 2. The condition number of the approximate Gram matrix [K + I is O(N).
Proof. Property (1) follows trivially from the fact that K has ((N) nonzero entries, whereupon matrix-vector multiplication only requires ((N) operations. Considering (2), since K is positive definite, the smallest possible eigenvalue of K + of is of Then note that
[R S |r =
K R
Σ
( x00) < N maxl((()) +(1))-
(6)
Under the assumptions above 1 S V i, whereupon (x(() ((ti < e² V i,j. Hence, have that x k + T 1++ +1++0 o which is O(N) as claimed.
Theorem 2 immediately implies the following corollary, which is also novel for GRFFs
Lemma 1. Solving the sparse linear system. Consider solving (K + + v, RV This can be achieved with the conjugate gradient method in 0 ( time.
b, where
Proof. Using the conjugate gradient method, it is known that the system can be solved in (.(K+07) iterations' (Shewchuk, 1994), which by property (2) above (N1/2). Each iteration involves matrix-vector multiplication, which is ((N) in our case due to property (1). Combining gives a total time complexity of (N8/22.
We remark that this is substantially less than the ((N3) time complexity of exact GP methods that use K2 rather than K. It is also foorag to see that Theorem 2 and Lemma 1 will continue hold we only consider a subset of the nodes of the graph, c.g. just considering set of training nodes of cardinality Nuran I N.
3.2 FROM PATHWISE CONDITIONING TO CONJUGATE GRADIENTS
We now introduce the three-step 'recipe' of posterior inference using GRFs; kernel initial- isation, hyperparmelerr learning and posterior inference. We will also analyse the overall time and space complexity of this workflow. App. C.1 gives further heuristic guidance for practitioners, including for choosing the number of walkers n.
Kernel initialisation. We compute the Gram matrix using Alg. 1, which involves sampling n random walks for every node on the graph. This yields a sparse kernel approximation:
K TT
(())),) N*N.
In practice, K does not need to be materialised as we can replace the patrix-vecttor product Ke with two fast matrix-vector products ( Te). Each computed in lincar time.
Hyperparaneter learning. Denote the training data Dy = {æ, y}, containing training nodes æ and corresponding noisy observations y. We learn the hy hyperparameters 0, such as observation noise and the modulation function f. by maximising the log marginal likelihood,
N
0()) - H 2 log det det(H)) "y
log(2.))
2
'Following prior work, we use "solving" as shorthand for reducing the error on our approximation of by factor of (Maddox et al., 2021), here set to 10-2, We will see that this is empirically robust.
5

---

Published as a conference paper at ICLR 2026
where Hg (K + 1 We use the Adam optim iser and estimate the gradient,
HHe 00
1
HHg
(Hg" y) - tr H 00
(9)
using iterative methods (Gardner et al., 2018; Lin et al., 2024a). These avoid explicit matrix inverses via iterative linear system solvers such as conjugate gradients (CGs) (Hestenes and Stiefel, 1952; Shewchuk, 1994). Since CGs rely matrix-vector multiplication, this allows
us to leverage the efficient structure of GRFs. Mcanwhile, the trace term is estimated using Hutchinson's ss trace estimator (Hutchinson, 1990),
HH 0 00
HHg 00
OH,
tr H
E Ho
~
H
(10)
S
00
==1
where 2, are random probes satisfying E [z,z5]
= I. This gives a batch of linear systems,
Hg [by " ... Ug] = y 2 21.,,, Eg}
(11)
which can be solved via iterative methods. The solutions allow us to estimate VL.
Posterior inference. We perform posterior inference using pathwise conditioning (Wilson et al., 2020; 2021) and iterative methods- a combination that has attracted recent interest in the literature (Lin et al., 2024b; 2025). This allows us to exploit the efficient structure of GRFs. In particular, pathwise conditioning expresses a sample from the posterior as sample from the prior with an additional correction term,
M = g-) K ())r (Kaa + t (y (g-(g(>)) +e)),
(12)
where (-) is any node of the graph 9, gy a sample from the posterior g is sample from the prior, and e (N, to compute (Kn= + 1) (y (g(x) + E)), which again avoids the explicit inverse and leverages sparse matrix multiplication." Once more, we use CGs (Hestenes and Stiefel, 1952; Shewchuk, 1994) linear system solver, though alternatives have recently been proposed (Lin et al., 2023; 2024c). The structure of the GRFs kernel also admits efficient sampling from the prior via g = &u with w N(O, I)," since (ov((bu))= = 4$ K.
o This facilitates the use of iterative linear system solvers
Algorithm complexity. Kernel initialisation takes ((N) time, since fixed number of random walks are simulated from all N nodes. Training and inference are dominated by CG solvers, with (N332) time complexity (Lemma 1). All stages use sparse matrices (c.g.
K + I) with O(N) nonzero entries, giving overall space complexity O(N).
4 EXPERIMENTAL RESULTS
Here, we present empirical results demonstrating the scalability and practical effectiveness of the GRF-GPs model. In each case, full experimental details are provided in App. C.
4.1 COMPUTATION COMPLEXITY AND ABLATIONS
the importance of efficient implementation. We benchmark posterior inference on synthetic graphs under two GRF implementations. First, consider dense baseline that uses GRPs, but explicitly materialises the N x N kernel approximation and computes its inverse. Second, we take the sparse GRF method described
Dense sparse GRFs:
*An alternative to solving this sparse linear system is to the formation to reduce the dimensionality the features {M E RR, whilst preserving their dot products in expectation (Dasgupta and Gupta, 2003). At the cost of sacrificing sparsity, we then use the Woodbury Identity to efficiently solve a smaller linear system. We describe this App. "g is an N-dimensional vector corresponding to a sample evaluated al all N nodes. One could consider subset of nodes, where the prior sample g(-) now correspouds to the vector's (-)th entry.
trans-

---

Published as a conference paper at ICLR 2026
in Section 3.2, storing the random walk trajectories and solving the corresponding linear systems with CG methods. Table I summarises the results, with full measurements provided in the App. C.2 (Table 2 and Table 3). For graph with 8192 nodes, we observe 50x speedup in total wall-clock time.
Table 1: GRF-GPs have sub-quadratic time scaling and linear memory scaling. Empirical scaling exponents (# s.d.) for memory usage, kernel initialisation, training, and inference with respect to graph size N. In the table, an entry b indicates scaling (NN).
Kernel
Memory
Kernel init. time Training time Inference time
GREs (Dense)
2.00 + 0.00
1.21 + 0.06
1.97 + 0.38
2.16 + 0.33
GRFs (Sparse)
1.00 + 0.00
0.81 + 0.04
1.04 + 0.04
1.04 + 0.05
Figure 2 shows log- log scaling curves. Exponents from the asymptotic regime match those shown in Table 1. As expected, GRPss attain linear memory and initialisation cost, and sub-quadratic training and inference, scaling graphs with 1M nodes. The near-linear runtime trends in training and inference reflect the fixed iteration budget of sparse linear solves; conditioning effects have not yet dominated at these scales.
Gooo0o Dooad
oooooo
Doog
(a) Memory usage
(b) Initialisation time
(c) Training wall-clock time (d) Inference wall-clock time
Figure 2: GRFs scale better (blue curve) when sparsity is leveraged. Scaling experiments for the GRF-GPs. Yellow: brute-force dense implementation. Blue: sparse implementation. Panels (d) correspond to memory footprint, kernel initialisation time, training time and inference time, respectively. The dense model is limited to 8,192 nodes due to its higher memory demands.
Importance sampling ablation. As discussed in Section 2, the key insight of GRFs is that one can replace a function of a weighted adjacency matrix W with a Monte Carlo estimate. This estimate constructed using random walks, weighted by (1) the product of traversed edge weights and (2) the per-walk probability under the sampling mechanism. Following Reid et al. (2024b), one can investigate the signíficance of this principled approach by instead constructing naive random walk-based empirical kernel, without appropriate reweighting. In particular, we replace line 8 of Alg by
p()(prefix_aubxalle Batpss_ubaaaakkkka 1]]
traversed_ edge_degeeeeeeeeee weights) - tanngtthpprefixxuuwaaaaae (13)
removing normalisation by p(prefiixsubwalk) Crucially, this set of features still defines valid kernel on 9, but it is no longer an unbiased estimate of a power series of W. A similar 'ad-hoc" kernel was used in the context of transformer position encodings by Choromanski et al. (2022). Full empirical results are reported in App. C.3, where we find this modification substantially dogrades regression performance. Intuitively, failing to upweight long, unlikely walks by 1/p(prefix ( makes it challenging to model longer-range dependencies.
4.2 REGRESSION TASKS
Next, we apply our method to regression with variety of real-world datasets.
1. Traffic speed prediction To assess predictive capability, we begin with traffic speed forecasting task (Figure 7) on the San Jose freeway sensor network (Chen al., 2001). We follow the setup of Borovitskiy et al. (2021). Experiment details can be found in App. C.5.

---

Published as a conference paper at ICLR 2026
We compare three kernel configurations by measuring the negative log probability density (NLPD) and the root mean squared error (RMSE) of the sanximmm---posterinnn (MAP) predictions. We consider (1) the exact diffusion kernel K Kan: (2) GRF kernel in 'diffusion shape' (namely, the modulation function frozen to approximate Kaim with a learnable lengthscale); and (3) GRF kernel with a flexible, fully learnable modulation function.
Figure 3 (a)-(b) reports the test NLPD and RMSE as a function of the number of random walks per node n. As n increases, the variance of the Monte Carlo approximation K drops. It better captures the underlying graph structure, yielding more accurate predictions. Note that the fully-learnable GRF kernel consistently outperforms the diffusion-shaped variant, highlighting the benefit of implicit kernel learning via a flexible modulation function.
In addition to greater flexibility via learnable f,, another reason GRFs are able to outperform Ka may be that their inbuilt sparsity is actually a sensible inductive bias. Pairs of graph nodes only have nonzero covariance if their respective ensembles of random walks hit, which more likely if they are nearby in 9. This means that a node's predictions depend mostly on information from its local neighborhood, whilst still sampling longer dependencies with lower probability. In contrast, dense kernels can sometimes be prone to the 'oversmoothing" effect they capture spurious long-range correlations driven by noise (Keriven, 2022).
2. Wind interpolation on the globe. Next, we consider the task of interpolating monthly average wind velocities from the ERA5 dataset (Hersbach et al., 2019), from set of locations on the Aeolus satellite track (Reitebuch, 2012). Our problem setup follows that of Wyrwal et al. (2024) and Robert-Nicoud et al. (2023). We discretise the surface of the globe (formally, the manifold 5²) by computing k-nearest neighbours graph from the observation locations. This yields a graph 9 with 10K nodes, with which we can apply our scalable GRF-GPs algorithm. The task is to predict the velocity fields of a held out lest set.
The test NLPD and RMSE of the diffusion-shape and fully-learnable GRF kernels are shown in Figure 3 (c)- Similarly, the predictions improve as n increases. We provide full results and visualisations in App. C.6. This type of implicit manifold GP regression approximating a (possibly unknown) manifold by computing a ncarest neighbour graph 9 and then performing inference therein is a rich area of active research (Borovitskiy et al., 2021; Dunson et al., 2021; Fichera et al., 2023). This is an exciting possible application of GRFs; we hope our initial example will spur future work.
155
303
Naraerr C'Wilen
Nanker Walken
(a) Traffic NLPD
(b) Traffic RMSE
(c) Wind NLPD
(d) Wind RMSE
Figure 3: GREs outperforms diffusion baselines in regression tasks. Panels (a) (d) report test NLPD and RMSE versus the number of random walkers n. Blue: GRF kernel with fully learnable modulation; orange: diffusson-ssape GRF. Shading shows +1 s.d. On Traffic, the learnable GRF surpasses the exact diffusion kernel once n ž 500. On Wind, the exact diffusion kernel is omitted due to ((N)) cost. Again, the fully-learnable GRF kernel consistently achieves lower NLPD and RMSE than the diffusion-shape variant.
4.3 SCALABLE AND RODUST BAYESIAN OPTIMISATION
Having demonstrated the scalability of GRFs (Section 4.1) and their efficacy for GP regression (Section 4.2), we now use them to perform efficient Bayesian optimisation (BO). We consider large graphs with up to 10c nodes, where exact posterior inference becomes prohibitively expensive. For the acquisition strategy we use Thompson sampling, drawing samples from the posterior over the objective function and selecting maximisers as the next
8

---

Published as a conference paper at ICLR 2026
query point (Russo et al., 2018; Thompson, 1933). Posterior sampling is made efficient by pathwise conditioning, given in Equation (12). Alg. in App. C.7 gives full details.
Datasets and baselines. For datasets, consider range of synthetic and real-world graphs. First, we maximise a variety of scalar functions on grids, community and circular graphs, chosen to have different properties, e.g. multimodality and periodicity. Next, identify 'influential' (high node degree) users in a range of social works Enrön, Facebook, Twitch and YouTube. Lastly, we predict the physical location with the greatest windspeed for the ERA5 dataset studied in Section 4.2, considering three different altitudes where the wind behaviour is known to qualitatively different (Wyrwal et al., 2024). In each case, compare efficient BO method with random scarch, breadth first scarch and depth first search policies. In almost all instances, our algorithu achieves lower regret, showing the benefit of uncertainty-awant strategies for large-scale optimisation on graphs.
0.60
04
0.05
(a) Single-modal
(b) Multi-modal
(c) Community
(d) Circle
0.50
0.04
Regret
6.50
0.2
0.03
6.45
0.02
6.40
0.02
0.0
0.00
6.35
0.01
(c) Enron Email
Facebook
(g) Twitch
(h) YouTube
Requet
0.2
0.8
0.0
0.6
0.8
0.4
0.4
0.4
0.5
0.5
0.30
0.5km
( 2km
(k) 5km
Algorithms
0.25
0.4
0.4
Bogget
0.3
0.3
BPS
0.2
DFS
0.2
0.1
0.00
Random
0.1
0.0
-0.05
600 BD Iteratiin
600
1000
BI iteration
steration
Figure 4: GRF-based BO achieves lower regret than search-based baselines in most datasets. Fach panel shows the regret curve of BO for the following datasets: (a)- (d) synthetic datasets, (e)-(h) social networks, and (i)-(k) windspeed in the ERA5 dataset.
4.4 FUTURE WORK: SCALABLE VARIATIONAL GPs FOR CLASSIFICATION
Lastly, we evaluate GRF-GPs on a multi-class node classification task using the Cora citation network benchmark (McCallum et al., 2000). In this non-conjugate inference setting, we handle the non- non-Gaussian likelihood via variational inference (Leibfried et al., 2020). Pathwise conditioning for classification nontrivial (Wilson al., 2021); defer full treatment to a future paper. We can nonetheless assess the performance of GRFs, even without explicit. time complexity guarantees like Lemma Details are provided App. C.8. Once again, sparse GRF kernels achieve very strong performance.
5 CONCLUSION
We demonstrated how graph random features (GRFs), a recently-introdueed Monte Carlo algorithm, can be used to speed training and inference with Gaussian processes on discrete input spaces. Under mild assumptions, GRFs support ((N3/2) time complexity inference much faster than (NN)) for their exact counterpart with probabilistic accuracy guarantecs. This translates to substantial wall-clock time speedups, and unlocks scalable Bayesian optimisation on massive topologies with little or no sacrifice in performance.

---

Published as a conference paper at ICLR 2026
ETHICS AND REPRODUCIBILITY
Ethics. Our work is methodological and does not raise direct ethical concerns. Nenetheless, advances in scalable graph-based ML may amplify risks if misapplied, either by malicious actors or through unforescen downstream consequences.
Reproducibility. All datasets are freely available online, with links to the original sources provided. The code publicly available at http:://ithbb.oom. Brtdhu..co/MMtthewWhaant /Matthe22hangg73/ /rapph-
Trustworthy, scalable and robust machine learning for graphs. Much of graph- based machine learning research to date has focused on graph neural networks. Kernel methods provide an alternative which is inherently more explainable. Our algorithms are applicable to broad range of real world settings, including many tasks in social networks. We demonstrate strong scalability, efficiency and performance whilst retaining the benefits of explainability.
7 FUNDING AND ACKNOWLEDGEMENTS
JAL is supported by the University of Cambridge Harding Distinguished Postgraduate Scholars Programme. AW acknowledges support from EPSRC via Turing AI Fellowship under grant EP GPPVV22279/1 and project FAIR under grant EP EP/V /0///0068///// /11688831..0 AW is also supported by the Leverhulme Trust via CFI. RET is supported by the EPSRC Probabilistic AT Hub (EP/Ye (OPP/02887//1) IR gratefully acknowledges funding from Trinity College, Cam- bridge and from Google. IR was a student researcher al Google DeepMind at the time of submission. This work was funded in part by an unrestricted gift from Huawei.
REFERENCES
Borovitskiy, V., Azangulov, I., Teremin, A., Mostowsky, P., Deisenroth, M., and Durrande, N. Matérn Gaussian processes graphs. International Conference on Artificial Intellgence and Statistics, 2593-2601, 2021.
Carter, L. L., and Cashwell, D. Particle- torticle-ransppott simulation with the Monte Carlo method [Technical report], 1975.
Chen, C., Petty, K., Skabardonis, A., Varaiya, P., and Jia, Z. Freeway Performance Measurement System: Mining Loop Detector Dala. 80th Annual Meeting of the Transporlation Research Board, 2001. https:/ https///pople.eecs.e beople..ee..brreeee.....u
ps.dir/ Mentnng
Choromanski, K. M. Taming graph kernels with random features. International Conference Machine Learning, 5964 5977, 2023.
Choromanski, K., Lin, H., Chen, H., Zhang, T., Sehanobish, A., Likhosherstov, V., Parker-Holdees J., Sarlos, T., Weller, A., and Weingarten, T. From block-toeplitz matrices to differential equations on graphis towards a general theory for scalable masked transformers. International Conference on Machine Learning, 3962-3983, 2022.
Chung, F. R. Spectral graph theory (Vol. 92). American Mathematical Soc., 1997.
contributors, O. OpenStreetMap (Data set], 2024.
Dasgupta, S., and Gupta,
elementary proof theorem Johnson and Lindenstrauss. Rondom Structures & Algorithms,, 22(1), 60-65, 2003.
Dunson, D. B., Wu, H.-T., and Wu, N. Spectral convergence of graph Laplacian and heat kernel reconstruction in Loo from random samples. Applied and Computational Harmonic Analysis, 55, 282-336, 2021.
Fichera, B., Borovitskiy, S., Kranse, A., and Billard, A. G. Implicit manifold gaussian process regression. Advances in Newrul Information Processing Systems, 36, 67701 67720, 2023.
Freksen, C. B. An Introduction to Jolnson-Lindenstrree Lindenstranss Transforms. Corr, 2021. https: https://arxive /aaxiv.oogr/ abs/ 0as/2103.00560
10

---

Published as a conference paper at ICLR 2026
Gardner, J., Pleiss, G., Weinberger, K. Q, Bindel, D., and Wilson, A. G. GPyTorch: Blackbox Matrix- Matrix Gaussian Process Inference with GPU Acceleration. Advances in Neural Information Processing Systems, 2018.
Hendricks, J. S., and Booth, T. E. MCNP variance reduction overview. Monte-Care Carlo Methods and Applications in Neutronies, Photonics and Statistical Physics: Proceedings of the Joint Los Alamos National Taboratory- Commissariat I, I'energie Atomique Meeting Held Cadarache Castle, Provence, France April 22-26, 1985, 83 83-92, 2006.
Hersbach, H., Bell, B., Berristord, P., Biavati, G., Horányi, A., Muñoz Sabater, J., Nicolas, J., Peubey, C., Radu, R., Rozum, I., and others. ERA5 monthly averaged data on single levels from 1979 to present. Copernicus Climate Change Service (C8s) Climote Data Store (CDS), 10, 252-266, 2019.
Hestenes, M. R., and Stiefel, E. Methods of Conjugate Gradients for Solving Linear Systems. Journal of Research of the National Bureau of Standards, 496), 1952.
Hutchinson, M. A Stochastic Estimator of the Trace of the Influence Matrix for Laplacian Smoothing Splines, Communications in Statistics Simulotion and Computation, 19(2), 1990.
Jones, D. R.,, Schonlau, M., and Welch, W. J. Efficient global oplimization of expensive black-box functions. Journal of Global Optimization, 18(4), 455-492, 1998.
Keriven, N. Not too little, not too much: A theoretical analysis of graph (over)smoothing. Advances in Neural Information Processing Systems, 36, 2022.
Leibfried, F., Dutordoir, V., Johm, S. T., and Durrande, N. Tulorial on Sparse Gaussian Processes and Variational Inference, 2020, December. https://doi.org 10.48550/arr /araaaa./X............ 20..2001.13962
Leskover, J., and Krevl, A. SNAP Datasets: Stanford Lorge Network Dataset Collection, 2014.
Lin, A., Ament, S., Balandat, M., and Bakshy, E. Scaling Gaussian Processes for Learning Curve Prediction via Latent Kronecker Structure. Neurips Bayesian Decision-Making and Uncertainty Workshop, 202/bb
Lin, A., Ament, S., Balandat, M., Eriksson, D., Hernández-Labato, J. M., and Bakshy, E. Scalable Gaussian Processes with Latent Kronecker Structure. International Conference on Mochine Learning, 2025.
Lin, A., Antorán, J., Padhy, S., Janz, D.,, Hernández-Lobato,, J. M., and Terenin, A. Sampling from Gaussian Process Posteriors using Stochastic Gradient Descent, Advances in Neural Information Processing Systems, 2023.
Lin, J. A., Padhy, S., Antorán, J., Tripp, A., Terenin, A., Sæepesvári C., Carrnandez-obbato. J. M., and Janz, D. Stochastic Gradient Descent for Gaussian Processes Done Right. International Conference Learning Represenlations, 2024c.
Lin, J. A., Padhy, S., Mlodozeniee, B., Antorán, J., and Hemáández-Lobato, J. M. Improving Linear System Solvers for Hyperparameler Optimisation in Iterative Gaussian Processos. Advonces in Neural Information Processing Systems, 2024a.
Maddox, W. J., Kapoor, S., and Wilson, A. G. When are iterative gaussian processes reliably accurate?. Arzio Preprint Arrio:2112, 15246, 2021.
Matthews, A. G. d. G., van der Wilk, M., Nickson, T., Fujii, K., Boukouvalas, A., León- Villagrá, P., Ghahramani, Z., and Hensman, J. GPflow: A Gaussian process library using TensorFlow. Journal of Machine Learning Research, 18(40), 1-6, 2017. http:// lttp://jmmlr.n org/ papers/ /at/ss///////// 111816-577.hta
McCallum, A., Nigam, K., Rennie, J., and Seymore, J. Automating the Construction of Internet Portals with Machine Learning. Information Retrieval, X2), 127 163, 2000.
Močkus, J. On Bayesian methods for seeking the extremum. IFIP Technical Conference on Optimiza- tion Techniques, 400-404, 1974.
Rahimi, A., and Recht, B. Random features for large-scale kemnel machines. Advances in Neural Information Processing Systems, 20, 2007.
Rasmussen, Gaussian processes machine learning. Sunmer school on machine learning.
pp. 63-71. Springer, 2003.
Reid, I., Choromanski, K., Berger, E., and Weller, A. General graph random features. Arziv Preprint Ar1iv:2310.04850, 2023.
11

---

Published as a conference paper at ICLR 2026
Reid, I,, Dubey, K. A., Jain, D., Whitney, W., Ahmed, A., Ainslie, J., Bewley, A., Jacob, M., Mehta, A., Rendleman, D., and others. Linear transformer topological masking with graph random features. Arriv Preprint Arriu:2410. 03462, 2024b.
Reid, I., Markon, S., Choromanski, K., Turner, R. E., and Weller, A. Variance- Couplings for Random Features. Arrio Preprint Arzie:2405. 16541, 2024a.
Reitebuch, O. The spaceborne wind lidar mission ADM-Acolus, In Atmospherie physics: Background- methods trends, pp. 815-827. Springer, 2012.
Robert-Nicoud, D., Krause, A., and Borovitskiy, V. Intrinsic Gaussian Vector Fields on Manifolds. Arriv Preprint Arriv: 18824, 2023.
Russo, D. J., Van Roy, B., Kazerouni, A., Osband, I, and Wen, Tutorial on Thompson Sampling. Foundations and Trends in Machine Learning, I1(1), 1-96, 2018. https://doi.or 10. 1561 220100200000000
Shewcluk, J. R. An Introduction the Conjugate Gradient Method Wilhout the Agonizing Puin [Technical report], 1994.
Smola, A. J., and Kondor, R. Kernels and regularization graphs. Learning Theory and Kernel Machines: 16th Annual Conference on Learning Theory and 7th Kernel Workshop, Coll/kernu 2003, Washington, DC, USA, August 24-27, 2009. Proceedings, 144-158, 2003.
Thompson, W. R. On the likelihood that one unknown probability exceeds another in view of the evidence of two samples. Biometrika, 25(3-4) -), 285 294, 1933.
Titsias, M. Variational learning of inducing variables in sparse Gaussian processes. Artificial Intelli- gence and Statistics, 567-574, 2009.
van der Wilk, M., Dutordoir, V., John, S., Artemev, A., Adam, V., and Hensman, J. A Framework for Interdomain and Multioutput Gaussian Processes. Arriv:2003.014 2020. https:/ happss//axxvvvv /ttx/v.orgg/////// abs/2003.01115
Wan, X., Osselin, P., Kenlay, H., Ru, B., Osborne, M. A., and Dong, X. Bayesian optimisation of functions on graphs. Advonces in Neural Information Processing Systems, 36, 43012-43040, 2023.
Wilson, J. T., Borovitskiy, V., Terenin, A., Mostowsky, P., and Deisenroth, M. P. Efficiently Sampling Functions from Gaussian Process Posteriors. International Conference on Machine Learning, 2020.
Wilson, J. T., Borovilskiy, V., Terenin, A., Mostowsky, P., and Deisenroth, M. P. Pathwise Coudi- tioning of Gaussian Processes. Journal of Machine Learning Research, 22(1), 2021.
Wyrwal, K., Krause, A., and Borovitskiy, V. Residual Deep Gaussian Processes on Manifolds. Arziv Preprint Arzio:24 1100161, 2024.
Yang, J., Sindhwani, V., Fan, Q.: Avron, H., and Mahoney, M. W. Random laplace feature maps for semigroup kernels on histograms, Proceedings of the IEEE Conference on Compuder Vision and Pattern Recognation, 971-978, 2014.
A FULL GRF ALGORITHM
To complement the pseudocode provided in Alg. 1, Alg. 2 provides a more detailed expla- nation of how one can estimate graph node kernels using graph random features (GRFs). The motivated reader is invited to consult the works of Reid et al. (2023) and Choromanski (2023) for the original accounts, including further intuitions and a proof of unbiasedness.
Algorithm 2: Constructing random feature vector pp(1) RVW to approximate K, (W)
Inputs: weighted adjacency matrix W E NNN for graph 9 with N nodes, vector of unweighted node degrees d E R modulation function f: (N (0)) - R, termination probability Phott E (0, 1), node E N, number of random walks to sample E N.
2 Output: random walk feature vector p(i)= E R.W.
initialise: (p() <- 0 for 1, 1,..., n
12

---

Published as a conference paper at ICLR 2026
A Algorithm 2: Constructing a random feature vector (p(1) RR to approximate K, (W)
initialise: load t-
initialise: current_nood __de <- i
initialise: terminated t- False
initialise: walk_length
while terminated False do
10
(t) (current_node t (t(i) [current_node) + load x f
walk length = walk 1ength+t
12
new new_node (ctrreee
assign to one of neighbours
13
load load dlcurrent node
W[current_non node, new sew_node] update load
P-Pui
14
current_node < _ed_
15
terminated +- (t Unif(0, 1) < Phan))
b draw RV to decide on termination
16
end while
17 end for
18 normalise (p(i) (pii//n
B EFFICIENTLY SOLVING LINEAR SYSTEMS WOODBURY FORMULA
(Kk (K - b wITH THE
In this appendix, provide another algorithm for efficiently solving linear system (K + of = b with the of Woodbury matrir identity formula and Johnson-Tindenstrauss Transform (JLT) (Freksen, 2021). This algorithm has time complexity O(N2m + m²), where < N is the number of the output dimensions (a hyperparameter of the JLT algorithm). While this approach appears promising, we emphasise that our investigation here is prelim- inary. A more thorough evaluation of its empirical performance and potential trade-offs is left to future work.
Take the decomposition of K of the form k= = TΦT Construct random Gaussian matrix G E Nx with entries taken independently at random from the Gaussian distribution with mean == 0 and standard deviation $$ K,K,, where K₁ = - $G E Nxm (and with strong concentration guarantees -1 for m of logarithmic in N order). We can then approximately rewrite (K + b as of (IN + UUT)-"b b, for U K,
1. By the JLT, we can unbiasedly approximate
Now, we can apply the following special case of the celebrated Woodbury Matrix Identity formula:
(IN + = -
(14)
Therefore, conclude that the solution to our linear system be approximated as:
u IN U( (
(15)
The expression on the right side can clearly be computed in time O(Nm + m²) since brute-foree force inversion of X (I,, + UT Rmxm takes O(m²) time and expression U(XU 6)) for U E Nxm be computed in m(N time. Thus, since the computation of K, takes time mN) total time complexity is O 2m + m²).
13

---

Published as a conference paper at ICLR 2026
This approach, using dimensionalityy reduction of GRFs to replace the inverse of an N x N matrix by the inverse of an m X m matrix, provides an interesting alternative to relying on sparse operations to achieve speedups.
C EXPERIMENT DETAILS
We provide experimental details this section. All experiments are conducted on single compute node equipped with an NVIDIA GeForce RIX 2080 Ti GPU (11 GB memory).
C.I CHOOSING GRF HYPERPARAMETERS GUIDANCE FOR PRACTITIONERS
In this appendix, we provide further practical guidance for practitioners when choosing the number of random walks and (if desired) the maximum walk length mmux.
Choosing n. Theorem gives precise formula for choosing the number of walkers n to guarantee an accurate kernel estimate with high probability. In principle, one could use this to derive the minimum n required for sharp estimate, given the constant c, the maximum permissible deviation t, and the maximum permissible probability of deviation P (IP(i)TT (o() - [ > t). However, in practice find this to be unnecessary: choosing to be small multiple of the average node degree already works well. As seen in Figure 3, performance tends to improve as n increases, at the cost of decreasing kernel sparsity and thus slower wall-clock times. We recommend choosing that balances the practitioner's performance and efficiency requirements.
Choosing Imax- In Section 3.1, noted that in implementations it is often convenient and memory-efficient to only sample walks up to some maximum length I mmxx This way, the number of modulation function terms (u)) that must be learned is finite and fixed. We t=0
emphasise that this is not requirement for the time complexity guarantees in Section 3. it is practical (as opposed to mathematical) detail. In principle, one could choose mmax to be sufficiently large that all n walkers will be shorter with high probability, avoiding any e.g. App. A.1 by Reid et al. (2023) for mathematical bound.7 However, in practice we find that choosing Imax to be some modest fraction of the graph diameter sufficient for good performance. In each experiment, we report mmxx the respective appendix.
truncation
C.2 TIME AND SPACE COMPLEXITY MEASUREMENTS
This section reports experimental details for the scaling results in Figure 2 and Table 1.
Synthetic data. We generate synthetic signals on ring graphs of increasing size: = 25 25,2,,... ,..., 220 nodes. The groundtruth functions are smooth periodic functions on the nodes with additive Gaussian noise (o² = 0.1). For graphs with more than 8192 nodes, we only use the sparse GRF implementation, since the dense adjacency matrices exceed the available GPU memory. Random feature matrices & are constructed using 100 random walks per node, with halting probability Phnh = 0.1. Walks longer than hops are truncated.
Measurements taken. For each graph size, and across 5 random seeds, measure:
The memory footprint of the random feature matrices
The random-walk preproccssing time for constructing +
Training wall-clock time, measured as total optimiser runtime over 50 epochs.
Inference wall-clock time, mcasured as posterior mean and covariance evalu- ation time on the test set.
"As the number of walkers grows, the minimum value of max to ensure that all walkers are shorter than Jas with high probability scales logarithmically with
14

---

Published as a conference paper at ICLR 2026
The dense implementation uses the GPilow library for kernels with explicit adjacency materialisation, while the sparse implementation uses a GPyTorch library to implement kernels with customised sparse linear operators maximise efficiency (Gardner et al., 2018; Matthews et al., 2017; van der Wilk et al., 2020). Full empirical measurements are shown in Table 2 and Table 3.
Table 2:
Memory and time measurements for dense implementation: mean + s.d.
Graph Size
Memory (MB)
Kernel init time (s)
Training time (s)
Inference time (s)
32
0.024 + 0.000
0.115 + 0.017
1.726 + 0.336
0.019 + 0.002
64
0.094 + 0.000
0.205 = 0.012
1.430 0.219
0.018 0.002
128
0.375 + 0.000
0.421 + 0.025
1.403 + 0.116
0.017 = 0.002
256
1.500 + 0.000
0.840 = 0.044
1.371 + 0.152
0.016 0.002
512
6.000 + 0.000
1.800 + 0.069
1.370 + 0.288
0.021 = 0.004
1024
24.000 = 0.000
4.189 + 0.204
2.465 + 0.595
0.045 = 0.006
2048
96.000 = 0.000
10.546 + 0.107
7.680 + 1.649
0.173 = 0.001
4096
384.000 = 0.000
31.749 = 1.246
40.376 + 4.080
1.043 = 0.006
8192
1536.000 + 0.000
104.839 + 2.026
307.188 35.938
7.572 0.000
Table 3:
Memory and time measurements for sparse implementation: mean s.d.
Graph Size
Memory (MB)
Kernel init time (s)
Training time (s)
Inference time (s)
32
0.004 + 0.000
0.160 + 0.033
4.103 + 0.216
0.066 + 0.007
64
0.008 + 0.000
0.168 + 0.022
3.823 + 0.136
0.061 + 0.008
128
0.015 + 0.000
0.202 + 0.022
4.036 = 0.191
0.066 + 0.007
256
0.030 + 0.000
0.271 + 0.030
4.369 + 0.349
+
512
0.059 + 0.000
0.379 + 0.021
4.395 + 0.619
0.077 0.019
1024
0.118 0.000
0.552 + 0.024
4.549 + 0.593
0
2048
0.235 + 0.000
0.973 + 0.039
4.416 + 0.320
0
4096
0.470 + 0.000
1.790 + 0.028
4.185 + 0.252
0.078 0.015
8192
0,938 + 0.000
3.481 + 0.074
4.247 + 0.143
0.076 + 0.006
16384
1.876 + 0.000
6.764 + 0.052
5.117 + 0.518
0.100 + 0.016
32768
3.751 + 0.000
13.297 = 0.050
6.623 + 1.048
0.129 + 0.040
65536
7.501 + 0.000
26.569 + 0.063
12.566 + 1.188
0.254 0.061
131072
15.001 + 0.000
53.012 + 0.156
31.534 + 6.376
0.651 0.175
262144
30.000 = 0.000
105.901 + 0.514
60.488 + 17.849
1.216 + 0.443
524288
60.000 = 0.000
212.671 + 0.758
111.672 31.377
2.068 0.775
1048576
120.000 + 0.000
426.074 + 1.562
245.060 + 65.159
4.947 + 1.226
Scaling factor estimation. We estimate empirical complexity exponents by fitting the measured runtime and memory data a power-law model,
y R aNo
using ordinary least squares in log log space, where is the number graph nodes. Uncertainty in the slope b is quantified with 95% confidence intervals derived from the t-
15

---

Published as a conference paper at ICLR 2026
distribution. To capture asymptotic scaling behavior, fits are restricted to sufficiently large graphs: dense GP experiments are fit for N 2 2", while sparse GP experiments are fit for N 2 216. The fitted coefficients and together with confidence intervals and R² values, are summarised in Table 4.
Table 4: Fitted power-law scaling coefficients for memory usage, random-wall initialisation, training, and inference time. Each row reports multiplicative constant a, exponent b with 95% confidence interval, and coefficient of determination R². Fits performed in log log space.
Kernel
a
b
95% CI (b)
R²
Memory (MB)
Sparse
1.37 10 1.00
[1.00, 1.00]
1.00
Dense
2.29 x 10
2.00
[2.00, 2.00]
1.00
Kernel init time (s)
Sparse
3.58 10
0.81
[0.73, 0.88]
0.97
Dense
1.22 x 10-3
1.21
[1.09, 1.33]
0.99
Training time (s)
Sparse
1.32 x 10
1.04
[0.96, 1.12]
1.00
Dense
3.93 10
1.97
[1.20, 2.73]
0.96
Inference time (s)
Sparse
2.79 x 10-6
1.04
[0.93, 1.14]
0.99
Dense 1.92 10-8
2.16
[1.50, 2.81]
0.97
C.3
ABLATION STUDIES
This section reports the results of the ablation experiment described in Section 4.1, where we replace the GRF estimate of a function of a weighted adjacency matrix by an ad-hoc random walk-based kernel. As described in the main body, line of Alg is replaced by
p()lpretix_subeee subvaax([-l]] += (II traversed_er edge__dee_____ee (ecerreed__deee_eehhttse x [[Length(prerxxxeebbag subbalkk)) (16)
removing the normalisation factor ppprefix_subwalee
Data synthesis. We consider synthetic dataset, consisting of regular 30 30 mesh graph (900 nodes). We compute a ground truth diffusion kernel Kir exp(- L) on this mesh graph with known length scale 8 = 10 (hidden from the models), and sample function from the corresponding GP, shown in Figure 5. Noisy observations are made at 10% of the nodes, indicated by black dots. The task is to predict missing measurements.
Kernels comparison. For GP training and inference, we consider three kernels: the exact diffusion kernel Kant = of exp(-3L), a GRF kernel K, and an ad-hoc random walk kernel The Kud.hoc as per Eq. (16). The learned maximum-Bn sxximumna-posterion predictions (posterior mean) are shown in Figure 5 (b)-(d), and the RMSE and NLPD are reported in Table 5. For random walk-based kernels
and K Kadr we sample 10,000 walks per node, truncating any walk exceeding 10 steps. Models trained using the Adam optimiser with learning rate of 0.01 for 1,000 iterations.
Clearly, the ad-hoc kernel fails to capture the underlying structure, producing inaccurate predictions. This shows that a principled importance sampling approach is essential for random walk-based kernels to perform well practice.
16

---

Published as a conference paper at ICLR 2026
80
-32
-64
(a) Ground truth
(b) Diffusion
404
-4
(c) GRFs
(d) Ad-hoc GRFs
Figure 5: The ad-hoc kernel fails to capture longer-range relationships. Panel (a): Ground-truth function on a 30 x 30 mesh graph; black dots mark noisy observations at 10% of the nodes. Panels (b d): Posterior means inferred with the exact diffusion kernel, the GRF kernel, and the ad-hoc kernel, respectively. Unlike the principled GRF estimator, the ad-hoc variant produces poor predictions and misses the underlying structure.
Table 5: The ad-hoc kernel yiclds much worse predictive accuracy. Test RMSE and NLPD for the diffusion kernel, principled GRF kernel, and the ad-hoc variant. The ad- hoc kernel exhibits substantially higher RMSE and NLPD.
Kernel
RMSE
NLPD
Diffusion
0.262
0.090
GRIs
0.339
0.339
Ad-hoc GRFs
0.573
1.265
C.4 CONJUGATE GRADIENTS CONVERGENCE
Here, provide extra discussion to supplement Section 3.1. In particular, we experimentally confirm the (NN) bound on the condition number of the GRF Gram matrix and the (N3322 time complexity for solving the corresponding linear system.
Theorem asserts that the condition number of the GRF Gram matrix K is O(N), assuming that e = Ir (max,,Jel..N) (,,d/(1- - -D))" is constant across the class of graphs considered. We empirically confirm that this is the case below Figure 6. Lemma also uses that the conjugate gradient (CG) method can solve the corresponding linear system in / 1s k + TT iterations (Shewchuk, 1994), which ultimately unlocks our (NN322) scaling. This implicitly assumes that we run CG to fixed error ratio in our case, 10-2 which follows convention for efficient GP methods in the literature (Maddox et al., 2021). We find this CG termination criterion to be empirically robust.
17

---

Published as a conference paper at ICLR 2026
Constrion Namber
Ieccation
10
051 Staeettion
10h
1g1
27
27
27 0mph Stn
27
21
100
132
100
Cooph Size
Ootditio Number
(a) Conditional number scaling
(b) CG iterations scaling
(c) CG iterations vs
Figure 6: Empirical scaling behaviour of CG iterations (a) Condition number vs.
graph size for Erdös-Rényí graphs with edge probability p p= p= 0 NV) The linear fit in log-log space shows that & ~ NO.52 which lies well within the upper bound proved in in Theorem 2, y to fixed relative (c) CG iterations plotted directly against A log-
ie., x (K + I is ((N). (b) CG iterations to solve the (K + o + 10-2 as N increases.
error tolerance log fit gives iter x 0.53 which match the expected (VK)) behaviour.
C.5 REGRESSION TASK: TRAFFIC SPEED PREDICTION
Here we provide further details for the first regression experiment: prodicting traffic speeds the San Jose freeway sensor network (Chen et al., 2001), following the setup of Borovitskiy et al. (2021).
Dataset. We use the San Jose freeway sensor network combined with OpenStreetMap data to construct graph with 1,016 nodes and 1,173 edges (contributors, 2024). Traffic speed measurements (in mph) are available at 325 sensor locations. These values are normalised (zero mean, unit variance), and the data is split into a training set of 250 randomly selectod nodes and test set of the remaining 75 nodes.
Kernel approximation with GRFss We used two variants of GRFs kernels. The first GRF kernel uses a diffusion-shape modulation function f== (-3/2) This is truncated power series expansion of the diffusion kernel, where the Icarnable hyperparameters are length scale 3 and kernel variance of The second kernel directly learns the modulation coefficients (f)))= which are initialised randomly and learned via log marginal likelihood. maximum length of 10, and vary the number of walks per node E (1,2, 4, ..., 8192]. Since the traffic network contains roughly 1,000 nodes, we also include the exact diffusion kernel K as a baseline. The kernel configurations are:
For both GRF variants, we fix Phat = 0.1 and truncate walks
Exact Diffusion: Kaiff oî eX),
Diffusion-shape K; ft (-3/2) 1!
Fully-learnable R Si learned directly.
Regression task. We apply GP inference using the 250 labeled nodes as training data to predict traffic speeds all 1,016 nodes in the network. The kernel hyperparameter and noise variance of are learned by maximising the log marginal likelihood, using Adam. Posterior inference then yields the predicted mean p and covariance Σ of the latent traffic speed function over the graph.
To quantify accuracy, we compute the negative log probability density (NLPD) and root mean squared error (RMSE) the 75 test nodes between the true speeds y ees and the MAP estimate ft:
18

---

Published as a conference paper at ICLR 2026
N
RMSE
(p: v/)
test
f=1
M
1 Nest
NLPD
log (og
Zz Drrai))
The experiment is repeated five times with different random seeds. The results are shown in Figure 3 (a]-(b) in the main text.
Capturing global and local patterns. Using the visualisation toolkits by Borovitskiy et. al. (2021), we illustrate the GRF-GPs posterior inference results on the San Jose traffic network in Figure 7. The left panel provides a global view over the full network, while the right panel zooms in on a specific highway junction. We observe that the global inferred mean (top left) captures large-scale spatial variation across the network speeds are higher main freeway segments and lower in peripheral or downtown regions. Notably, in the zoomed-in in view (top right), the model successfully distinguishes speeds across tightly packed lanes running in opposite directions. Despite spatial proximity, the posterior assigns signif- icantly different mean values to adjacent but directionally distinct segments, demonstrating that GRF-GPs capture connectivity-awd distance. The bottom visualises posterior uncertainty, with standard deviation plotted over the full graph (bottom left) and zoomed in section (bottom right). These results confirm that GRF-GPs respect both global graph structure and local topology, delivering interpretable and spatially coherent predictions on complex, real-world networks.
patterns rather than relying solely on Euclidean
Meen
Mean Eooomed
10
39
20
:0
Standard
Figure 7: Posterior inference using GRF-GPs on the San Jose traffic network. Top left: Mean predictions across the full graph. Top right: Zoomed-in directional differences be- tween closely spaced lanes. Bottom left: Postcrior uncertainty over the network. Bottom right: Zoomed view reveals local variation in confidence. Coloured dots are sensor nodes; white dots indicate training nodes.
C.6 REGRESSION TASK: WIND VELOCITY INTERPOLATION
19

---

Published as a conference paper at ICLR 2026
Here we provide further details about the wind velocity interpolation task from the ERA5 dataset (Hersbach et al., 2019). Our problem setup follows that Wyrwal et al. (2024) and Robert-Nicoud et al. (2023).
Dataset. We use the average wind velocity field from the ERA5 dataset at three altitudes: 0.1 km, 2 km, and 5 km. The surface of the globe (formally, the manifold S²) is discretised al a resolution of 2.5 longitude by 2.5e latitude, yielding a k-nearest, neighbours graph 9 with roughly 10K nodes, on which we apply our scalable GRF-GPs algorithm. The task is to predict the velocity fields on the held-out test nodes. The locations along the Acolus satellite track (1441 nodes) serve as training data, while all remaining nodes are treated as the test set.
Figure 8: Ground-truth wind velocity field from the ERA5 datasct at 0.1 km above sea level. Black vectors show local wind velocities. Red dots mark 1441 Aeolus satellite track locations, used as training data in the interpolation task (Reitebuch, 2012).
Figure 9: Predicted wind velocity field using GRF-GPs. Blue vectors represent. MAP predictions (GP posterior mean).
Figure 10: Prediction uncertainty (GP posterior covariance) using GRF-GPs. Brighter regions indicate higher uncertainty, which is significantly reduced near satellite track.
20

---

Published as a conference paper at ICLR 2026
Figure 11: Absolute error between ground-truth and MAP-predices velocities. GRF-GPs achieve accurate predictions, with error patterns aligned with uncertainty estimates.
Experiment setup. We use the fully-learnable and diffusion-shape GRF kernels, varying the random walk budget, similar to the method described in App. C.5. Note the exact diffusion kernel Kittt cannot be applied on this large graph. We measure the NLPD and RMSE to evaluate the kernel performance. The results are shown in Figure 3 (c-d) in the main text.
Uncertainty-aware wind velocity interpolation. Figures 8 11 visualise GRF-GPs inference on the ERA5 wind dataset at 0.1 km altitude. For visualisation clarity, the k- nearest neighbour graph on the globe downsampled. Figure shows the ground-truth wind field with training node positions marked in red. Figure 9 shows the MAP prediction, and Figure 10 shows the posterior uncertainty, which notably reduced along the Aeolus satellite track. Finally, Figure 11 displays the absolute error field.
C.7 LARGE SCALE BAYESIAN OPTIMISATION ON GRAPHS
Here, we describe our evaluations of the performance and scalability of GRF-GPs on Bayesian optimisation (BO) tasks, detailed in Section 4.3. We test the methodology across three settings: (1) four synthetic graph benchmarks, (2) four real-world social network datasets for identifying influential users, and (3) three wind interpolation datasets. First, let us describe the benchmark datasets.
1. Synthetic benchmarks. We consider four synthetic graph benchmarks.
Unimodal function on grid: a function with a smooth central peak, discretised on 1000 x 1000 grid graph.
Multi-modal function on grid: a function with several randomly placed peaks, discretised on 1000 x 1000 grid graph.
Community graph: a community graph generated via stochastic block model (SBM), with nodes in a community C, assigned score by sampling from N (M,, of).
Circular graph: sinusoidal function defined on ring, discretised into a k-nearest neighbour graph with 10 nodes.
All signals are perturbed with Gaussian noise (o = 0.1). Random features o are computed with 100 walks per node, with halting probability Phatt 0.1. Random walks longer than hops are truncated.
2. Social networks benchmarks: identify the most influential user.
We consider four real-world social network datasets (Table 6) from the Stanford Network Analysis Project (SNAP) (Leskovec and Krevl, 2014), with up to 1.1M nodes. Each node represents user in the network. Following Wan al. (2023), for user influence, and the task is to identify the most 'influential' users in each network.
node degree proxy
21

---

Published as a conference paper at ICLR 2026
Table 6: Summary of four SNAP datasets used for large-scale BO experiments. Fach dataset corresponds to a user-level social network, with node degree used as proxy for influence.
Dataset
Nodes
Edges
Maximum Degree
Description
YouTube
1,134,890
2,987,624
28754
Youtube online social network
Facebook
22,470
171,002
709
Facebook page-page network with page names.
Twitch 168,114
6,797,557
35279
Social network of Twitch users.
Enron
36,652
183,831
1383
Email communication network from Enron
3. ERA5 wind velocity field: predict the location with greatest wind speed.
To demonstrate the utility of GRFs for BO on manifolds, we use the ERA5 wind datasets at three altitudes. Full details of dataset processing are provided in App C.6.
A Igorithm Baselines. We compare GRF-based Thompson sampling against three search heuristics:
Random scarch: uniformly samples nodes without replacement.
Breadth-first search (BFS): sequentially expand observed nodes along the adjacency structure breadth-fisst order.
Depth-first search (DFS): sequentially expand observed nodes along the adja- cency structure in depth-firss first order.
BO setting. In each experiment, algorithms are initialised with up to 1,000 samples and then run for up to 1,000 BO iterations, repeated across five random seeds. At each iteration, we report simple regret, defined as the difference between the global maximum and the best function value observed so far.
Algorithm 3: Graph Thompson Sampling with GREs
Inputs: black-box function h, candidate nodes s_all, initial sample size N_0, number of BO steps T.
Output: augmented dataset (x_obs, y_obs). 1__0
initialise x_ obs {x, ],,] 1 n; Unif(x_aau all)
initialise y_obs <- (h(x;) + + for t 1,...,
model (x y_ y_obs)
7
s_t PosteriorSummle(model, s_all)
8
x_t + ArgMax(s_t)
y_t + h(x_t) +
10
x_obs - x_ obs U x_t
11
y_obs + y obs Uy_t
12 end for
13 return (x_obs, obs
C.8 CLASSIFICATION TASK: CORA CITATION NETWORK
Here we provide more experimental dctails about the classification task on the Cora scientific citation network (McCallum et al., 2000). This experiment highlights the application of GRF-GPs in a more challenging, non-conjugate inference setting.
22

---

Published as a conference paper at ICLR 2026
Dataset and preprocessing. The Cora dataset is standard benchmark in graph-based machine learning. It consists a citation network, where each node corresponds to scientific publication and each edge represents citation. Each publication is labelled with one of seven machine learning topics (Figure 12). While Cora also includes textual features, we focus solely on the graph structure. We extract the largest connected component of the citation graph, resulting in subgraph with 2,485 nodes and 5,069 edges.
(a) Ground-truth labels. Fach color denotes class
(b) Prediction enrors using graph GP with Matérn kernel. Red nodes misclassified
Figure 12: Cora dataset classification with graph GP.
Sparse variational inference for classification. In classification tasks, the likelihood functions are usually non-Gaussian (softmax), so the posterior not analytically tractable. and the M inducing nodes as Define latent
Denote the Numin training nodes function values at the training inputs as h = (h(x) x E x) and function values at induc- ing nodes as u (u(z) : E z). Assume a GP prior p(u) N N0, K and likelihood p(y, h4) (softmax). Choose Gaussian variational postcrior q(u) N(u, Σ) and induce the marginal ((h) p(h u)q(u)) ddu. Under this approximation we maximise the evidence lower bound (ELBO): LELBO E Nenis Eq(h.) (lg.))log p(v, h,)] KL(g(u) p(u)). This varia- tional treatment replaces the intractable posterior with tractable family and supplies principled objective (a lower bound on log p(y)); it yields coherent predictive distributions by integrating over q(h) rather than relying on point approximations, which espocially important when the likelihood breaks conjugacy.
Experiment setup. We compare classification accuracy across exact kernels (diffusion and Matérn) and the GRF kerne.. We use 80/20 train-test split the largest connected component the graph. The goal is to predict the class labels of all nodes based on the graph structure alone. All models are trained using softmax likelihood. Optimisation is performed for up to 1000 iterations using the Adam optimiser. To reduce uncertainty and assess variability, each configuration is repeated five times with different random seeds. We also measure the sparsity of the resulting GRF kernels. Results are reported in Table 7, showing that with a sufficient number random walkers, the flexibility of the GRF kernel allows it to capture the graph structure effectively and perporoor the exact kernel baselines.
C.9 WHY AREN'T SVGPs EFFICIENT ON GRAPHS?
In this section, we briefly compare our method sporse variational Gaussian processes (SVGPs) (Titsias, 2009), a different efficient. GP framework popular for problems in Euclidean space. However, a core difficulty with applying SVGPs to graphs is the fact that, in order to compute a smaller kernel matrix at some subset of inducing points, we still need to compute the entire Gram matrix. This is in general (NN3).
23

---

Published as a conference paper at ICLR 2026
Table 7: The GRF kernel reaches highest accuracy in the Cora benchmark. Classification accuracy on the Cora dataset with different graph kernels. With n = 16384 walks per node (22.17% non-zero entries), the GRF kernel outperforms the diffusion and Matérn kernels (Borovitskiy et al., 2021).
Kernel
Form
Accuracy
Diffusion
K diff exp(-8L)
85.31 + 0.61%
GRFs
87.04 0.53%
Matérn
K Matém = ( L)
86.72 + 0.31%
Comparison with Euclidean setting. Suppose V denotes the full set of N training datapoints, and Vad C 1/ denotes subset of m inducing points, For e.g. the Gaussian kernel, one can compute
(x, x3))
K ind.
exp
E mxm
(17)
20²
i.jE
in (m²) time and space complexity. One need not compute the kernel at points in 1 Vnn Conversely, for graphs, is not straightforward to compute some specific [Kalig without first materialising the full N N matrix
K.(W) Σ W', a, G (0,,,,
10,
(18)
where W is the adjacency matrix for the entire graph 9. One then extract the corre- sponding subset of entries corresponding to the m inducing points. This means that methods like SVGP may not actually provide time complexity gains, since we are still bottlenecked by O(N3) to compule the kernel at the inducing points in the first. place.
Extra baseline. Notwithstanding the above, we can implement SVGPs on graphs to compare performance to GRFs, even if this baseline may not be truly efficient in practice. On the traffic speed prediction task, we trained SVGP models using graph diffusion kernel with 150 inducing points. We used Adam optimiser with learning rate of 0.01 for 1000 steps. The results are shown in Figure 13. As expected, SVGP underperforns compared to the exact kernel. GRFs perform better even with modest number of walkers.
deily stape)
dflly leamable)
1050
16.5
shape)
102.5
Kay louac)
16.0
Kurr (exect)
1000
15.5
NIPO
97.5
RMSE
15.0
95.0
14.5
14.0
13.5
13.0
10P
10 Nunther of Walkers
10*
Number of Walkers
Traffic NLPD with SVGP baseline
(b)Traffic RMSE with SVGP baseline
Figure 13: Extra SVGP baseline. Companion results to Figure 3, with additional SVGP baseline. SVGP as expected performs worse than the Kait baseline both NL.P and RMSE. GRF's achieve better NLPD & RMSE beyond 100 walkers.
24


Reviews:
[{"review_id": "b7F4qrnFT0", "reviewer": "Reviewer_MPVY", "summary": "The authors propose a scalable GP methodology for discrete input spaces using graph random features (GRFs) in combination with conjugate gradient (CG) methods and Hutchinson’s trace estimator. A theoretical analysis is presented showing that the computational complexity of the proposed approach scales as $O(N^{3/2})$, which is a significant improvement over the $O(N^3)$ complexity when using exact kernels. Extensive numerical studies are presented to illustrate the performance of the proposed approach.", "questions": "- Section 3.1 mentions that the condition number bound requires $||\\phi(i)||_1 \\leq c$. Can you clarify if this bound depends on the choice of modulation function?\n- Can the authors derive approximation bounds showing how the error bound on $||\\hat{K} - K_\\alpha||_F$ impacts the error in the posterior predictions?\n- The numerical experiments show empirical scaling exponents of ~1.04 for training/inference (Table 1), which is nearly linear rather than N^{3/2}. What was the convergence tolerance used for CG? Was it an upper bound on the residual and/or maximum number of iterations?\n- Was a preconditioner used for the CG solver? \n- For the wind interpolation task where GRFs outperform exact kernels: (a) is this trend consistent across multiple runs? (b) Can you comment on whether this is due to overfitting of exact kernels or beneficial regularization from GRF sparsity? (c) does the error over the training/validation set show similar patterns?\n- The Woodbury identity approach (Appendix B) trades sparsity for lower-dimensional inversions. Under what conditions is this preferable to the sparse CG approach? Computational cost and approximation quality comparisons would be valuable.", "limitations": "", "rating": 8, "confidence": 4, "soundness": 3, "presentation": 4, "contribution": 3, "strengths": "The paper is well-written and organized. Relevant background and motivation is discussed clearly and details of the setup used for numerical studies are described in sufficient depth in the appendices. The three-step procedure (initialization, hyperparameter learning, inference) is clearly explained. Algorithm pseudocode aids reproducibility.\n\nThe application of GRFs to scalable GP inference is novel, though GRFs themselves were introduced recently (Reid et al. 2023, Choromanski 2023). The key theoretical results are: (1) Theorem 2 proving O(N) condition number for the covariance matrix and (2) Lemma 1 showing $O(N^{3/2})$ complexity for CG solves. The paper presents numerical studies demonstrating 50× speedups on graphs with <10K nodes and  Bayesian optimization on graphs with >1M nodes. \n\nOverall, the paper presents solid methodological work demonstrating that GRFs enable scalable GP inference on large graphs.", "weaknesses": "The theoretical development is largely sound. Theorem 1 (from prior work) provides useful concentration inequalities that the paper builds upon Theorem 2's proof for $O(N)$ condition number follows from standard operator norm arguments. The constant 'c' must be finite for the concentration bounds to hold. While the authors claim this is a “mild assumption,” it effectively requires the spectrum of W to satisfy a specific decay rate. This could be made more explicit.\n\nThe analysis in Lemma 1 that claims $O(N^{3/2})$ complexity for CG is rather crude for a formal theoretical result. In practice, CG is terminated when the residual falls below a specified tolerance ($\\epsilon$), requiring $O(\\sqrt{\\kappa} log(1/\\epsilon))$ iterations. The paper does not specify what tolerance was used, nor does it analyze how the choice of tolerance affects the approximation quality of the posterior. The $O(N^{3/2})$ estimate is equivalent to assuming that CG iterations are terminated after $\\sqrt{\\kappa}$ steps. \n\nThe empirical scaling in Table 1 shows ~1.04 (nearly linear), which the authors attribute to a 'fixed iteration budget’. This trend implicitly suggests that the theoretical O(N^{3/2}) bound may not be achieved in practice without careful tuning of convergence criteria.\n\nI was expecting to see approximation error analysis quantifying how using $\\hat{K}$ instead of $K_\\alpha$ affects the quality of the posterior. Since Theorem 1 provides concentration for individual entries, this should be feasible and would strengthen the contribution.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "PtI5Hje1bR", "reviewer": "Reviewer_sDRg", "summary": "Traditional GPs on graphs suffer from cubic time complexity (O(N^3)), limiting their scalability to small graphs. This paper leverages graph random features, which approximate graph kernels via sparse random walk–based feature maps, enabling sub-quadratic inference with strong theoretical guarantees. GRFs yield unbiased kernel estimates with exponential concentration bounds, allowing the GP covariance matrix to be approximated efficiently while preserving accuracy.", "questions": "1. Clarification on Practical Hyperparameter Sensitivity\n\nThe performance of GRFs depends on parameters such as the number of random walks, the halting probability, and the maximum walk length. Appendix C.1 partially answers Question 1, but is limited to empirical and qualitative descriptions. Could the authors provide more quantitative guidance on how these hyperparameters influence the trade-off between accuracy, sparsity, and computational cost? For example, is there an empirical scaling rule or automatic tuning strategy that generalizes across datasets? \n\n2. Comparison with Other Scalable GP Frameworks\n\nWhile the paper contrasts dense vs. sparse GRFs, it does not directly compare with other established scalable GP techniques (e.g., SVGP, SKI, or Kronecker GPs). Could the authors comment on how GRFs compare empirically or theoretically to these alternatives, especially in terms of scalability and approximation accuracy? \n\n3. Clarification on Bias and Positive-Definiteness\n\nSection 2 mentions a small $O(1/n)$ bias in diagonal kernel entries due to shared randomness among random walks. Could the authors clarify whether this bias ever affects posterior stability or positive definiteness in practice? Have they observed numerical issues during GP training? A short empirical note or ablation (e.g., comparing one vs. two independent random walk ensembles) would help readers understand if this bias is purely theoretical or occasionally relevant.", "limitations": "", "rating": 6, "confidence": 3, "soundness": 3, "presentation": 4, "contribution": 3, "strengths": "1. The authors provide rigorous complexity analysis and concentration bounds, proving that GRF-based GPs achieve ($O(N^{3/2})$) time complexity with probabilistic accuracy guarantees — a solid theoretical contribution that enhances credibility.\n2. The experiments are diverse, covering both synthetic and real-world tasks. Results consistently support the theoretical claims and demonstrate substantial speedups without significant performance loss.\n3. Despite the technical depth, the paper is well-structured and clearly written. Figures, pseudocode, and explanations make the ideas accessible.", "weaknesses": "1. The empirical evaluation focuses mainly on regression and Bayesian optimization; classification and non-conjugate inference are mentioned but not fully explored. This somewhat limits the generality of the method’s demonstrated usefulness.\n2. The performance of GRFs depends on parameters like the number of random walks and the maximum walk length. Although heuristics are provided, the sensitivity to these settings could be further studied or automated.\n3. While the paper compares dense and sparse GRFs, it does not extensively benchmark against other scalable GP frameworks (e.g., inducing point methods, variational GPs, or Kronecker-structured GPs), which would strengthen claims of superiority.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "0eFXf6vDUB", "reviewer": "Reviewer_RsGB", "summary": "This paper proposes Graph Random Features (GRFs) to construct learnable graph node kernels which estimate the sparse GP covariance function reducing the cubic computational complexity of exact GPs and applying these graph-based approximation technique for sparse GP estimation for Bayesian optimization via Thompson sampling. The contributions of this paper are considerably novel, as they build upon the previous work of Reid et. al. on GRFs specifically into GP domain and then provide an application of this technique by applying it into a BO setting showcasing the effectiveness and practicality.", "questions": "I have few concerns.\n-- The assumption that the number of random walkers 𝑛 is independent of 𝑁 is theoretically elegant but practically fragile. In dense or long-range graphs, a small 𝑛 can lead to high-variance kernel estimates and poor posterior calibration.\n-- The sparse approximation $\\hat{K} = \\Phi \\Phi^{\\top}$ is unbiased only in expectation, but in practice, the variance of the estimator is not analyzed in sufficient depth. The paper cites exponential concentration (Theorem~1) but never empirically validates how variance scales with $n$ or the graph topology. As a result, while computational speed is improved, the fidelity of posterior uncertainty (a key factor in Bayesian inference) may degrade unpredictably.\n\n-- The footnote in Section~2 admits an $\\mathcal{O}(1/n)$ bias for diagonal entries and notes  that removing it breaks positive definiteness. Yet, the authors use the biased estimator throughout without quantifying the impact of this bias \non GP likelihood or predictive uncertainty. This is a non-trivial issue --- GP models are sensitive to kernel conditioning, \nand even small biases can affect the posterior mean and variance.", "limitations": "", "rating": 6, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "The paper presents a well-motivated contribution that bridges recent advances in graph kernels and scalable GP inference. The use of GRFs for graph node kernels is an innovative idea that leverages random walks to approximate covariance structures efficiently. The theoretical development is sound and neatly connects GRFs with probabilistic kernel approximation, culminating in a clear argument for the claimed time complexity. The results claimed in Theorem2 and Lemma1 also establish theoretical underpinnings of the proposed work, contributing to the theoretical understanding of GRFs.\n\nEmpirically, the paper demonstrates wall-clock speedups and scalability to graphs with >1M nodes, which is impressive for a sparse GP-based method. The reported experiments, on regression and BO tasks, are diverse enough to illustrate both predictive and computational benefits. \n\nMoreover, the paper is clearly written and well-organized. Overall, the paper is theoretically and empirically solid and easy to follow for community working on GPs and BO.", "weaknesses": "While the core idea of the paper is novel, the experimental validation somewhat lacks depth and breadth. The comparisons are limited mainly to diffusion kernels and do not include strong baselines such as sparse variational GPs, approximation based on inducing-point (support set). Without these benchmarks, it is difficult to judge the real-world competitiveness of GRFs in GP regression and BO tasks. The paper also does not sufficiently explore the sensitivity of the method to hyperparameters such as the number of random walks, halting probability, or walk length.\n\nMoreover, the paper does not provide bounds on the GP posterior approximation error, so while kernel estimation is shown to be accurate, the implications for predictive uncertainty are less clear. Empirical metrics focus on RMSE and NLP, but uncertainty calibration and robustness are not analyzed. The BO results show consistent performance gains, but the lack of statistical testing limits the strength of these claims to some extent.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "E7H5XyWscX", "reviewer": "Reviewer_e8Qw", "summary": "This paper applies the random feature approach in kernel-based learning to graph problems (such as node classification/regression). For context, graph random features are a recently introduced stochastic estimator of a family of kernel matrices defined on graphs. The authors use these features for node regression and Bayesian optimization. They show an $O(N^{3/2})$ time complexity for kernel parameter estimation and demonstrate that a learnable kernel with random features outperforms the graph diffusion kernel, either in the full matrix form or in its random feature representation.", "questions": "See comments above.", "limitations": "", "rating": 6, "confidence": 5, "soundness": 3, "presentation": 4, "contribution": 2, "strengths": "A main bottleneck of Gaussian processes is the $O(N^3)$ scalability barrier. The proposed approach provably reduces the cost to $O(N^{3/2})$. Additionally, the authors demonstrate experiments to show that they can handle 1M graph nodes comfortably by using a single GPU.\n\nThe writing of the paper is excellent and easy to follow, balancing clarity and conciseness for graph random features (GRFs), which contain many subtle details.\n\n(Half strength, half weakness) GRFs are an interesting extension of random features from the Euclidean domain to the graph domain. It is pleasant to see the use of GRFs in machine learning tasks. On the other hand, this paper mostly concerns applying the concept, weakening the contribution. Moreover, as graph neural networks (GNNs) become so widely used for the tasks concerned in this paper, the advantages of GRFs are unclear.", "weaknesses": "On the theoretical side, given the prior work on GRFs and the machinery to perform GP inference and parameter estimation by using random features, the contribution of this work primarily lies in applying the former to the latter. A majority of the background for this application is known (e.g., the concentration of GRF in Theorem 1, the parameter estimation method in Eqn (8)--(11), and the posterior inference in Eqn (12)). To be fair, the authors do provide new results, such as the $O(N^{3/2})$ complexity for inference due to the $O(N)$ spectral norm of the kernel, but this contribution is less significant than the establishment of the GRF itself.\n\nOn the practical side, it is unclear how kernel/GP methods fare with GNNs. There are no such experiments in the paper. GNNs scale better than kernels and they are data-efficient for node-level problems. For example, for the Cora graph and the classification task that the authors discuss in Section 4.4, GNN can achieve over 80% accuracy with 20 training nodes per class (a labeling rate of 5.2% only). Can GRF be as data-efficient?\n\nAnother practical question is the Bayesian optimization over graphs, which has a finite search space. Since predicting the graph signal on each node is cheap and relatively accurate by using a GNN given a small training set, a natural heuristic is that one does not use Thompson sampling for the acquisition but only picks the best GNN predictions. The practical value of the paper can be strengthened by experimenting with such a simple baseline.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}]

Author rebuttal:

Decision: Accept (Poster)