General Graph Random Features
Authors: Isaac Reid
Krzysztof Marcin Choromanski
Eli Berger
Adrian Weller
Venue: ICLR 2024 poster (2024)
Keywords: Graphs
kernels
random walks
Laplacian
adjacency matrix
kernel learning
ordinary differential equation
neural network
Primary area: learning on graphs and other geometries & topologies
TL;DR: A novel random walk-based algorithm for unbiased estimation of arbitrary functions of a weighted adjacency matrix

Abstract: We propose a novel random walk-based algorithm for unbiased estimation of arbitrary functions of a weighted adjacency matrix, coined 
general graph random features (g-GRFs). This includes many of the most popular examples of kernels defined on the nodes of a graph. Our algorithm enjoys subquadratic time complexity with respect to the number of nodes, overcoming the notoriously prohibitive cubic scaling of exact graph kernel evaluation. It can also be trivially distributed across machines, permitting learning on much larger networks. At the heart of the algorithm is a modulation function which upweights or downweights the contribution from different random walks depending on their lengths. We show that by parameterising it with a neural network we can obtain g-GRFs that give higher-quality kernel estimates or perform efficient, scalable kernel learning. We provide robust theoretical analysis and support our findings with experiments including pointwise estimation of fixed graph kernels, solving non-homogeneous graph ordinary differential equations, node clustering and kernel regression on triangular meshes.

Published conference paper at ICIR 2024
GENERAL GRAPH RANDOM FEATURES
Isaac Reid', Krzysztof Choromanski Eli Berger",, Adrian Weller Uniiveriity of Cambridge, "Google DeepMind, *Columbia University, UUniversity of Haifa, Alan Turing Institute ir3370cam. ac. uk, kchoro@google.com
ABSTRACT
We propose novel random walk- based algorithm for unbiased estimation of arbitrary functions of a weighted adjacency matrix, coined general graph random features (g-GRFss. This includes many of the most popular ex- amples of kernels defined on the nodes of graph. Our algorithm en- joys subquadratic time complexity with respect to the number of nodes, overcoming the notoriously prohibitive cubic scaling of exact graph kernel evaluation. It can also be trivially distributed across machines, permit- ting learning on much larger networks. At the heart of the algorithm is a modulation function which upweights or downweights the contribution from different random walks depending on their lengths. We show that by parameterising it with a neural network we can obtain g-GREss that give higher-quality kernel estimates or perform efficient, scalable kernel learn- ing. We provide robust theoretical analysis and support our findings with experiments including pointwise estimation of fixed graph kernels, solving non non-homen non-homogeneous graph ordinary differential equations, node clustering and kernel regression on triangular meshes]
INTRODUCTION AND RELATED WORK
The kernel trick is a powerful technique to perform nonlinear inference using linear learn-
ing algorithms (Campbell 2002, Kontorovich et al.) 2008 Canu and Smola 2006 Smola have set of N datapoints X { F it replaccs Euclidean dot products x; x; with evaluations of a kernel function K : X' x X + R, cap-
and Schölkopf 2002 Supposing
turing the 'similarity' of the datapoints by instcad taking inner product between implicit (possibly simenie-dimessional) feature vectors in some Hilbert space Hk.
An object of key importance is the Gram malric K E NRN whose entries emumerate the pairwise kernel evaluations, K = [K(a;; my , Despite the theoretical rigour and empirical success enjoyed by kernel-based learning algorithms, the requirement to manifest and invert this matrix leads to notoriously poor O(N3) time-complexity scaling. This has spurred research dedicated to efficiently approximating K, the chief example of which is random features (Rahimi and Recht 2007): a Monte-Carlo approach which gives explicitly manifested, finite dimensional vectors o(x;) E R whose Euclidean dot product is equal to the kernel evaluation in expectation,
= [o(x,)
(1)
This allows one to construct a low-rank decomposition of K which provides much belter seal- ability. Testament to its utility, a rich taxonomy of random features exists to approximate many different Euclidean kernels including the Gaussian, softmax, and angular and lincar
kernels (Johnson] 1984 Dasgupta et al. 2010 Goemans and Williamson 2004 Choromanski
et al. 2020).
Kernels defined discrete input spaces, e.g. K A x W E with N the set of nodes of a graph G Smola and Kondor 2003. Kondor and Lafferty 2002 Chung and Yau 1999
"Equal contribution.
"Code available at httpe///cchiib.ooo//saaaar/d//er graph random. features

---

Published as a conference paper at ICIR 2024
enjoy widespread applications including in bioinformatics (Borgwardt et al. 2005 commu- uily detection (Kloster and Gleich 2014) and recommender systems (Yaimaa 2006]. More recently, they have been used in applications as diverse as manifold learning for deep gener- ative modelling (Zhou et al. 2020) and for solving single- and multiple-souse shortest path problems (Crane et al, 2017 p. However, for these graph-baso kernel methods the problem of poor scalability is particularly acute. This is because even computing the corresponding Gram matrix K is typically of at least cubic time complexity in the number of nodes N, requiring e.g. the inversion of an N x N malrix or computation of multiple matrix-matrix products. Despite the presence of this computational bottleneck, random feature methods for graph kernels have proved clusive. Indood, only recently was a viable graph random feo- ture (GRF) mochanism proposed by Choromanski (2023). Their algorithm uses an ensemble of random walkers which deposit a 'load" at every vertex they pass through that depends on i) the product of weights of edges traversed by the walker and ii) the marginal probability of the subwalk. Using this scheme, it is possible to construct random features {o(1)] C R such that ((())(()) gives an unbiased approvimation to the ij-th matrix element of the 2-regularised Laplacian kernel. Multiple independent approximations can be combined to estimate the d-regularised I Laplacian kernel with d & 2 or the diffusion kernel (although the latter is only asymptotically unbiased). The GRFs algorithm enjoys both subquadratic time complexity and strong empirical performance on tasks like k-means node clustering, and it is trivial to distribule across machines when working with massive graphs.
However, key limitation of GRFs is that they only address limited family of graph kernels which may not be suitable for the task at hand. Our central contribution is a simple modification which generalises the algorithm to arbitrary functions of a weighted adjacency matria, allowing efficient and unbiased approximation much broader class of graph node kernels. We achieve this by introducing an extra modulation function f that controls each walker's load as it traverses the graph. As well as empowering practitioners to approximate many more fixed kernels, we demonstrate that f can also be parameterised by neural network and learned. We use this powerful approach to optimisc g- g-GRFs for higher-quality approximation of fixed kernels and for scalable implicit kernel learning.
The remainder of the manuscript is organised as follows. In Sec. 2 we introduce gen- eral graph random features (g-GRF's) and prove that they enable scalable approximation of arbitrary functions of a weighted adjacency matrix, including many of the most popular examples of kernels defined on the nodes of graph. We also extend the core algorithm with neural modulation functions, replacing one component of the g-GRF mechanism with a neural network, and derive generalisation bounds for the corresponding class of learnable graph kernels (Sec. 2.11. In Sec. 3 we run extensive experiments, including: pointwise estimation of a variety of popular graph kernels (Sec. 3.1); sinulation of time evolution un- der non-homoge homogeneous graph ordinary differential equations (Sec. 3.2) kernelised k-means node clustering including on large graphs (Sec. 3.3); training neural modulation function to suppress the mean square error of fixed kernel estimates (See 3.4 and training neu- ral modulation function to learn a kernel for node altribule prediction on triangular mesh graphs (Sec. 3.5).
2 GENERAL GRAPH RANDOM FEATURES
Consider a directed weighted graph G(N, E, W == [wug]sen where: i) N == (1.....N) N) is the set of nodes; ii) & is the set of edges, with (t.j) E E if there is a directed edge from to j in G; and iii) W is the weighted adjacency matris, with wij the weight of the directed edge from i to j (equal to 0 if no such edge exists). Note that undirected graph can be described as dirocted with the symmetric weighted adjacency matrix W. Now consider the matrices (a(W) E R RNXN where a = =aa)) and ag E R:
Ka
Σ ag W k=O
(2)
We assume that the sum above converges for all W under consideration, which can be ensured with regulariser W oW, a E R4. Without loss of generality, we also assume

---

Published as a conference paper at ICIR 2024
that a is normalised such that op = 1. The matrix Ka(W) can be associated with graph function Kg :N x N + IR mapping from a pair of graph nodes to a real number.
Note that if G is an undirected graph then (W( automatically inherits the symmetry of W. In this case, it follows from Weyl's perturbation inequality (Bai et al, 2000) that (W( positive semidefinite for any given a provided the spectral radius p(W) = maxxea(W) ((A)) is sufficiently small (with A(W) the set of eigenvalues of W). This can again be ensured by multiplying the weight matrix W by regulariser E R4. It then follows that (W) can be considered the Gram malrix of a graph kernel function Kg.
With suitably chosen a = (a the class described by Eq. 2 includes many popular examples of graph node kernels in the literature (Smola and Kondor 2003 Chapelle et al. 2002 They measure connectivity between nodes and are typically functions of the graph Laplacian matrir, defined by L = I- W with W = -wws// /ddd Here, d; := E, way the weighted degree of node i such that W is the normalised weighted adjacency matrix. For reference, Table I gives the kernel definitions and normalised cocfficients ag (correspouding to powers W) to be considered later in the manuscript. In practice, factors in ax equal to a quantity raisod to the power are absorbed into the normalisation of W.
Name
Form
d-regularised Laplacian p-step random walk
(IN + +
(1+09)
(aIN 22
(a-
Diffusion
exp(-o'L//2) cos (LT/1)
d(
# (g)* {(-)) if lo even, (-1) if odd]
Inverse Cosine
Table 1: Different graph functions/keree kkennlss Kg X R. The exp and cos mappings are defined via Taylor scries expansions rather than clement-wisc, c.g. exp(M) == lim (IN + M/n)" and cos(M) == Re(exp(iM)) a and are regularisers. Note that the diffusion kernel is sometimes instead defined by exp(o>([x- -)) but these forms are equivalent up to normalisation. Also note that the p-step random walk kernel is closely related to the graph Matérn kernel (Borovitskiy et al. 2021
The chief goal of this work is to construct a random feature mop Φ(1) : + R with e N that provides unbiased approximation of (.(W) as in Eq. 1 To do so, we consider the following algorithm.
Algorithm 1 Constructing a random feature vector (1(1) E R to approximate (W)
Input: weighted adjacency matrix W e NNN, vector of unweighted node degrees (no. neighbours) E R& modulation function f (N U [0}) + R, termination probability Phal E (0, 1), node E N,, number of random walks to sample N.
Output: random feature vector (i( R
1: initialise: (i) <
2: for w = 1,.... do
3: initialise: load 6-1
4: initialise: current_nodee _odeeenoodee-eeeeee <- i
5: initialise: terminated +- False
6: initialise: malk_hength_t +
7: while terminated False do
8: ey(1) ((fr)))curettt_noe
+ ( +loadxf (walk_ length)
9:
walk Backs_legth__h +- walk length+1
10:
new node new_node
+ Unif [A' (urrentt
b assign to one of neighbours
11:
load <- loadx d[current_node] W [current_noo neweeeeewwm__ new__ode]]
& update load
1- Phals
12:
purrent_node t- new_nodw _ode
13: terminated +- (t Unif(O, 1) < Phatt)
draw RV to decide on termination
14: end while
15:
end for
16: normalise: 0y(1) m1mmm

---

Published conference paper at ICLR 2024
f(1)
Figure 1: Schematic for random walk graph (solid red) and an accompanying modu- lation lunction J (dashed blue) used to approximate an arbiirary graph node functiin KS.
This is identical to the algorithm presented by Choromanski (2023) for constructing features to approximate the 2-rogularised Laplacian kernel, apart from the presence of the extra modulation function f (NU(()) - line that upweights or downweights contributions from walks depending on their length (see Fig, []. We refer to dy as general graph rondom features (g-GRl's), where the subscript f identilies the modulation function. Crucially, the time complerity of Alg. 7 subquadratic the namber of nodes N, contrast to exact methods which are ((N))
We now state the following central result, proved in App. A.2
Theorem 2.1 (Unbiased approximation of Ka via convolutions). For two modulation fune- tions: f2 (N U (0)) + R, g g-GRFs (ø,, (1)) (6), constructed according to Alg. 1 give unbiased approrimation of K..
Kaly = [0 (1) [1)()))
(3)
for kernels with on arbitrory Taylor expansion a = (aw) provided that a = f Here, * is the discrete convolution of the modulation functions f1, f2; that is, for all k (NU(0)),
Σ Si(k - p)j((p) ax.
Clearly the class of pairs of modulation functions f1, S2 that satisfy Eq. 4 is highly degen- erate. Indeed, it is possible to solve for fi given any f2 and a provided f2(0) o 0. For instance, trivial solution is given by: fi( 0g f2(1) = T(i 0) with I(-)) the indicator function. In this case, the walkers corresponding to f2 are 'lazy', depositing all their load al the node al which they begin. Contribulions ions to the estimator Φy (i) (f) only come from walkers initialised at traversing all the way to rather than two walkers both passing through an intermediate node. Also great interest the case of symmetric modulation functions fi f2, where now intersections do contribute. In this case, the following is true (proved in App. A.3
Theorem 2.2 (Computing symmetric modulation functious). Supposing Ji S2 S, Bg. is solued by a function which is unique (up to a sign) and given by
Σ (mrtette..) ky kg
(
f(0)
(5)
Moreover, f() can be efficiently computed with the iterative formula.
f(0)= = =/200
- (i-p)f(++1)
(6)
f(i+1)
for 2 0.
2/10)
"Concretely, Alg. 1 yields a pair matrioes Φ12 (0()) NXN such that E(d db.
subquadratic time. Of course, explicitly multiplying the matrices evalnate every element of would be (LV but avoid this since in applications just evaluate Φi (1 where R vector and the brackets give the order of computation. This O

---

Published as a conference paper at ICLR 2024
For symmetric modulation functions, the random features ( and (% are identical apart from the particular sample of random walks used to construct them. They cannot share the same sample estimates of diagonal kernel elements [Ko will be biased.
Computational cost: Note that when running Alg. one only needs evaluate the modulation functions (i,2()) up to the length of the longest walk one samples. A batch of size b, (J1,2(0) can be pre-computed in time O(b2) and reused for random features corresponding to different nodes and even different graphs. Further values of f can be computed at runtime if is too small and also reused for later computations. Morcover, the minimum length b required to ensure that all m walks are shorter than b with high probability len(()) s b) 1-6, 5 & 1) scales only logarithmically with m (see App. [A.]). This means that, despite estimating a much more general family of graph functions, g-GRFs are essentially more expensive than the original GRF algorithm. Moreover, any techniques used for dimensionality reduction of regular GRFs (e.g. applying the transform (Dasgupta et al., 2010) or using 'anchor points' (Choromanski 2023 can also
be used with g- GRGFs, providing lurther elficiency gains.
Generating functions: Inserting the constraint for unbiasedness in Eq. back into the definition of Ka (W), we immodiately have that
(a(W) = K (W)K,,(W)
(7)
where K,, (W)- (W) = ft)) is the generating function corresponding to the sequence (0) This is natural because the (discrete) Fourier transform (discrete) convolu- tion returns the product of the (discrete) Fourier transforms of the respective functions. In the symmetric case f1 f2. it follows that
(W (Ka (W))
(8)
If the RHS has a simple Taylor expansion (c.g. (W) exp( W) so ( = exp( W this enables us obtain ʃ without recourse to the conditional sum in Eq. 5 or the iterative
2 )),
expression in Eq. 0 This is the case for many popular graph kernels; we provide some promi- nent examples in the table left. A notable ex- ception is the inverse cosine kernel.
Name
f(i)
d-regularised Laplacian
(d-24-25J
(2))(()22)
p-step random walk
(i)
As an interesting corollary, by considering the dill'usion kernel we have also proved that Σ 2npt 2 p(k-p)) M!
Diffusion
27
2.1 NEURAL MODULATION FUNCTIONS, KERNEL LEARNING AND GENERALISATION
Instead of using fixed modulation function f to estimate a fixed kernel, it is possible to parameterise it more flexibly. For example, ((N) (N U {0}) - by neural network (with a restricted domain) whose input and output. dimensionalitinn are equal to 1. During training, we can choose the loss function to target particular desiderata of g-GRFs: for example, suppress the mean square error of estimates of some particular fixed kernel (Sec. (.4)) or to learn kernel which performs better in a downstream task (Sec. 3.5 Implicitly learning Ka via p(a)) is more scalable than learning Ka directly because obviates the need to repeatedly compute the exact kernel, which is typically of O(N²) time complexity. Since any modulation function f maps to a unique by Eq. 4 it is also always straightforward to recover the exact kernel which the g-GRI's estimate, e.g. once the training is complete.
define neural modulation function
Supposing we have (implicitly) learned a, how can the learned kernel Ky be expected to generalise? Let 1hx, x + Ii,,, denote the feature mapping from the input space to the reproducing kernel Hilbert space Hk., induced by the kernel KG. Define the hypothesis set
H = (z + w T "kk (z) |as] (M [|wllz 1 1).
(9)
where we restricted our family of kernels so that the absolute value of each Taylor coefficient a; is smaller than some maximum value a; (M) Following very similar arguments to Cortes al. 2010), the following true.

---

Published conference paper at ICLR 2024
Theorem 2.3 (Empirical Rademacher complexity bound), For a fized sample (=:)==> the emppriiall Rademacher complexity R(H) bounded by
R(H) S 5 (M) a; P(W),
(10)
where p( (W) the spectral radius of the weighted adjacency matrix W.
Naturally, the bound on R(H) increases monotonically with p (W) Following standard ar- guments the literature, this immediately yiclds gencralisation bounds for learning kernels K2
We discuss this in detail, including proving Theorem 2.3 in App. A.4
3
EXPERIMENTS
Here
Lest the empirical performance of g-GREs, both for approximating fixed kernels (Secs 3.1 3.3) and with learnable neural modulation functions (Secs 3.4
3.1
UNBIASED POINTWISE ESTIMATION OF FIXED KERNELS
We begin by confirming that GRFs do indeed give unbiased estimates of the graph kernels listed in Table taking regularisers 0.25 and a 20 with Phalt 0.1. We use symmetric modulation functions f. computed with the closed forms where available and using the iterative scheme in Eq. D if not. Fig. 2 plots the relative Frobenius norm error between the true kernels K and their approximations with g- g-GRFs K (that is, K - - F/|| against the number of random wa walkers m. We consider 8 different graphs: small random Brdös-Rényí graph, larger Erdös-Rényí graph, binary tree, d-regular graph and 4 real world examples (karate, dolphins, football and eurosis) (Ivashkin 2023) They vary substantially in size. For every graph and for all kernels, the quality of the estimate improvee modest number of walkers.
grows and becomes very small with
ER IN= p0 02)
ER (N = 100. p= p=0.041
Binary tree UN 127)
d-regular (V 100. d= 101
norm aror
Freh norm Sererr
Frob. nori METOT
30
Terg
random
random Dolphins
13
No. random walks
10
No. random
22
Karabe 34)
Football (W- 115)
(N- 1272)
Frel. marm artor
Frob. METUT
2
No: random walls
No. random walks
No. random walks
No. random walks
Figure 2: Unbiased estimation of popular kernels on different graphs using g- g-GRss The approximation error (y-axis) improves with the number of walkers (x-axis). We repeat 10 times; standard deviation of the mean is shaded.
3.2 SOLVING DIFFERENTIAL EQUATIONS ON GRAPHS
An intriguing application of g- g-GR for fixed kernels is efficiently computing approximate solutions of timo-invariant non-homogencous ordinary differential oquations (ODEs) on graphs. Consider the following ODE defined on the nodes N of the graph G: (11)
dx(t) = Wx(l) + y(),
dt

---

Published as a conference paper at ICLR 2024
where æt)) E RN the state of the graph at time t, W R is a weighted adjacency matrix and y()) is a (known) driving term. Assuming the null initial condition x=0) = 0, Eq- T is solved by the convolution
x: (t)
exp(W(t ((r = EEER exp(W(t- ())(=)))((7))
(12)
(P(T)
where P is probability distribution the interval [0, equipped with (preferably efficient) sampling mechanism and probability density function p(r). Taking n E N Monte Carlo samples (Ti) ~ P. we can construct the unbiased estimator:
2(1) = 1 Σ p(Tz)) TI
exp(W(L - -,)))(T))
(13)
Note that exp(W(t - r,)) is nothing other than the diffusion kernel, which is expensive to compute exactly for large N but can be efficiently approximated with g-GRFs. Take exp(W((L (0(--77 7y)) = $ with +, = (0(1)) N x N matrix whose rows g-GRIs constructed to approximate the kernel at particular Ty. Then we have that
1
2 B (T(z) a(t) = P(T.)
(14)
3=1
which can be compuled in quadratic time (c.f. cubic if the heat kernel is computed exactly). Further speed gains are possible if dimensionnality reduction techniques are applied to the
g. g-GRFs (Choromanski] 2023 Dasgupta et al. 2010).
As an example, we consider diffusion on three real-world graphs with a fixed source at one node, taking W L (the graph Laplacian) and y()) = y = (1,0,0, 0, ...) The steady slate a(x)) = (-')- We simulate evolution un- 1 with
Diffusion simulation
walkers
dolshine
footbal
0.38-
der the ODE for sation timesteps and P uniform, approximating (xp(W(t- - T))) with different numbers of walk- m. As grows, the quality of approxima- tion improves and the (normalised) error on the final state [](()) (1))l/2 [1](1)] drops for ev- cry graph. We take 100 repcats for statistics and plot the results in Fig. 3 One standard devia- tion of the mean error is shaded. Phalt = 0.1 and
10 discrcti-
Sim.
0.34-
6.32
10
12
14
No.
Figure 3: ODE simulation error decreases as the number of walkers grows.
the regulariser is given by =
3.3 EFFICIENT KERNELISED GRAPH NODE CLUSTERING
As further demonstration of the utility of our new mechanism, show how estimates of the kernel K e))( can be used to assign nodes to l = 3 clusters. Here, we choose W to be the (unweighted) adjacency matrix and the regulariser o² 0.2. We follow the algorithm proposed by Dhillon et al.) (2004), comparing the clusters when we use exact and g-GRR-approximaatee ker- nels. For all graphs, m S 80. Table 2 re- ports the clustering error, defined by
Table 2: Errors kernelised k-means clustering when approximating the kernel exp(o" W) with g- g-GRFs.
Graph
N
Clustering error, E.
karate
31
0.08
dolphins
62
0.16
polbooks
105
0.12
football
115
0.02
databases
1046
0.10
eurosis
1272
0.09
E, = no. wrong pairs
cora
2485
0.01
(15)
citeseer
3300
0.04
N(N - 1)/2
This the number of misclassified pairs of nodes (assigned to the same cluster when the converse is true or vice versa) divided by the total number of pairs. The error is small even with a modest number of walkers and on large graphs; kernel estimates efficiently constructed using g-GRFs be readily deployed on downstream tasks where exact methods are slow.

---

Published as a conference paper at ICLR 2024
3.4 LEARNING f(N) FOR BETTER KERNEL APPROXIMATION
Following the discussion in Sec. 2. we now replaced fixed f with a neural modulation
function ((N) parameterised by a simple neural network with hidden layer of size 1:
(z) = sooftllu (W2OReLU + b1) +
(16)
where wi, b1, w2. by R and TnLU and ssoftplus are the ReL.U and softplus activation fune- tious, respectively. Bigger, more expressive architectures (including allowing f f(X) take negative values) can be used but this is found to be sufficient for our purposes.
We define our loss function to be the Frobenius norm error between target Gram matrix and our g- aRF-approxdm Gram matrix the small Erdös- Rényi graph (N = 20) with m = 16 walks; For the target, we choose the 2-regularised Laplacian kernel. We train symmetric f1 = f2 but provide a brief discussion of the asymmetric case (including its respective strengths and weaknessos) in App. A.5 On this graph, we minimisc the loss with the Adam optimiser and decaying learning rate (LR = 0.01, 0.975, 1000 epochs). We make the following striking observation: f(h)) does not generically converge to the unique unbiased (symmetric) modulation function implied by a. but instead to some different f that though biased gives a smaller mean squared error (MSE). This is possible because by downweighting long walks the learned f(N) gives estimators with a smaller variance, which is sufficient to suppress the MSE the kernel matrix elements even though it no longer gives the larget value in expectation. We then fix ((V) and use it for kernel approximation on the remaining graphs. The learned, biased f(,)) still provides better kernel estimates, including for graphs with very different topologies and much greater number of nodes: eurosis bigger by a factor of over 60. See Table for the results. Phalt = 0.5 and a = 0.8.
Naturally, the learned f(N) is dependent on the the number of random walks m; as m grows, the variance on the kernel approximation drops so it is intuitive that the learned f(V) will approach the unbiased f. Fig. 4 empirically confirms this is the case, showing the learned f(8) for different numbers of walkers. The line labelled co is the unbiased modulation function, which for the 2-rogularised Laplacian kernel is constant.
Table 3: Kernel approximation error with 16 walks and unbiased or learned mod- ulation functions. Lower is better. Brackets give one standard deviation of the last digit.
Figure 4: Learned modulation function with different numbers of random walkers m. It approaches the unbiased f(X) as m + 00
Medulation functions for different no, walkers
Graph
N
Frob. norm error K Uubiased
1.2
Learned
Small ER
20
0.0188(9)
0.0437(9)
Larger ER
100
0.0503(4)
0.0448(4)
16
Binary tree
127
0.0453(4)
0.0410(4)
0,5-
12
64
d-regular
100
0.0490(2)
0.0434(2)
128
karate
34
0.0492(6)
0.0439(6)
dolphins
62
0,0505(5)
0.0449(5)
02-
football
115
0.0520(2)
0.0459(2)
eurosis
1272
0.0551(2)
0.0481(2)
0.3
These learned modulation functions might guide the more principled construction of biased, low-MSE GRFs in the future. An analogue in Euclidean space is provided by structured orthogonal random features (SORFs), which replace a random Gaussian matrix with a HD- product to estimate the Gaussian kernel (Choromanski et al. 2017 Yu et al. 2016). This likewise improves kernel approximation quality despite breaking estimator unbiasedness.
3.5 IMPLICIT KERNEL LEARNING FOR NODE ATTRIBUTE PREDICTION
As suggested in Sec. 2.1 it is also possible to train the neural modulation function f(N) directly using performance on a downstream task, performing implicit kernel learning. We have argued that this is much more scalable than optimising K. directly.

---

Published as a conference paper at ICIR 2024
In this spirit, we now address the problem of kornel regression on triangular mesh graphs, as previously considered by Reid et al. (2023). For graphs in this dataset (Dawson- Haggorty] 2023), every node associated with a normal vector e(6) R3 equal to the mean of the normal vectors of its 3 surrounding faces. The task is to predict the direc- tions of missing vectors (a random 5% split) from the remainder. Our (unnormalised) predictions are given by 70) = E; ((
where K(N) is a kernel estimate con- structed using g-GRFs with a neural modulation function f(x) (see Eq. 3 The an- gular prediction error is
cos with 0g the angle between the true a and and
approximate g6) normals, averagod over the missing veclors. We train a symmelric pair f(M) using this angular prediction error the small cylinder graph (N 210) as the loss function. Then freeze f(x) and compute the angular prediction error for other larger graphs. Fig. 5 shows the learned f(l) as well as some other modulation functions corresponding to popular fixed kernels. Note also that learning ((8) al- ready includes (but is not limited to) optimis- ing the lengthscale of given kernel: taking W - W is identical to f(i) + f())@' V i.
Modulation functions ons for regracsion tadk
3.6-
reg 2 rog Lap dittuiion
3:24
Figure 5: Fixed and learned modulation functions for kernel regression
The prediction errors highly correlated be- tween the different modulation functions for a given random draw of walks; ensembles which 'explore' the graph poorly, terminating quickly repeating edges, will give worse g- g-GRF estimators for every For this reason, pute the prediction errors as the average normalised difference compared to the learned kernel result. Table | for every graph and every fixed kernel, meaning the learned kernel always performs best.
reports the results. Crucially, this differrnce is found to be positive
Table 4: Normalised difference in angular prediction error compared to the learned kernel defined by f(N) (trained on the cylinder graph). All entrics are positive since the learned kernel always performs best. We take 100 repeats (but only 10 for the very large cycloidal graph). The brackets give one standard deviation on the final digit.
Graph
Normalised A(pred error) c.f. learned
1-reg Lap
2-reg Lap
Diffusion
cylinder
210
+0.40(2)
+0.85(4)
-0.02((3)
teapot
480
+0)81(5)
+1.78(8)
-0.059(3)
idler-riser
782
+0.52(2)
+1.12(1)
-0.042(2)
busted
1941
+0.81(2)
+1.60(4)
-0.063(2)
4350
+2.13(5)
+5.3(1)
-0.067
cycloidal
21381
+0.065(1)
+0.13(1)
-0.01
It is remarkable that the learned ((0) gives the smallest error for all the graphs even though it was just trained on cylinder, the smallest one. We have implicitly learned a good kernel for this task which generalises well across topologies. It is also intriguing that the diffusion kernel performs only slightly worse. This is to be expected because their modulation functions are similar (see Fig. 5) so they encode very similar kernels, but this will not always be the case depending on the task at hand.
4 CONCLUSION
We have introduced 'general graph random features' (g-GRFs), novel random walk-based algorithm for time-efficient estimation of arbitrary functions of weighted adjacency matrix. The mechanism is conceptually simple and trivially distributed across machines, unlocking kernel-based machine learning on very large graphs. By parameterising the random features with simple neural network, we can futher suppress the mean square error of estionao and perform scalable implicit kernel learning.
component of

---

Published conference paper at ICLR 2024
5 ETHICS AND REPRODUCIBILITY
Ethics: Our work foundational. There are no direct ethical concerns that we can see, though of course increases in scalability afforded by graph random features might amplify risks of graph-based machine learning, from bad actors or as unintended consequences.
Reproducibility: To foster reproducibility, we clearly state the central algorithm Alg. Source code is available at httpss///gtuub.omo/ httpg/aaw........eeriiiillll /isaaiiariiddaggnneeeee /general graph random features, All theoretical results accompanied by proofs in Appendices A.I.A.4 where any assump- tions are made clear. The datasets we use correspond to standard graphs and are freely available online. We link suitable repositories in every instance. Except where prohibitively computationally expensive, results are reported with uncertainties help comparison.
RELATIVE CONTRIBUTIONS AND ACKNOWLEDGEMMETS
EB initially proposed using a modulation function to generalise GRFs to estimate the dif- fusion kernel and derived its mathematical expression. IR and KC then developed the full g g-GRFs algorithm for general functions of weighted adjacency matrix, proving all theoret- ical results, running the experiments and preparing the manuscript. AW provided helpful discussions and advice.
IR acknowledges support from Trinity College External Studentship. AW acknowledges support from a Turing AI fellowship under grant EP/VO25279/ and the Leverhulme Trust via CFI.
We thank Kenza Tazi and Austin Tripp for their careful readings of the manuscript. Richard Turner provided valuable suggestions and support throughout the project.
REFERENCES
Zhaojun Bai, James Demmel, Jack Dongarta, Axel Ruhe, and Henk van der Vorst. Templates for the solution of algebraic eigenvalue problems: pruclical guide. SIAM, 2000. URL html
https://wwx cs ucdavis. seddavis.edu/-ballle
Peter L Bartlett and Shahar Mendelson. Rademacher and gaussian complexities: Risk bounds and structural results. Journal of Machine Learning Reseorch, 4(No 2002. URL https://dl, acm. org/ doi http///pppp//fttt/// 995.//////11.65555499999 944944.
5555/944919
Karsten Borgwardt, Cheng Soon Ong, Stefan Schönauer, SVN Vishwanathan, Alex J Smola, and Hans-Peter Kriegel. Protein function prediction via graph kernels. Bioinfor- matics, 21(suppl 1):147 156, 2005. URL
bti1007
Viacheslav Borovitskiy, Iskander Azangulov, Alexander Terenin, Peter Mostowsky, Marc Deisenroth, and Nicolas Durrande. Matérn gaussian procosses on graphs. In International Conference on Artificial Intelligence and Slatislics, pages 2593-2601, PMLR, 2021. URL
httpps://doi. org 485444855//// aat:s0//XXiiv..o 2010. 15538
Colin Campbell. Kernel methods: survey of current techniques. Nearocompating, 48(1-4): 63 84, 2002. doi:
URL https:///ooi.org/10..00
20925-23120011
Stéphane Canu and Alexander J. Smola. Kernel methods and the exponential fam- ily. Neurocomputing, (9(7-9) 2006. doi: 10.1016/j.neucom.2DD 20w.omm.00050.0.000000000 URL
https:://oii orridiiooor//001 1016// neucom. 2005. 12. 009
Olivier Chapelle, Jason Weston, and Bernhard Schölkopf. Cluster kernels for semi-supervised
learning. Advonces in neurøl information processing systems, 15, 2002. URL https:
7/dl. acm. wtwp.////oi.........9.......1...
296866
10

---

Published as a conference paper at ICLR 2024
Krzysztof Choromanski, Valerii Likhoshersoov, David Dohan, Xingyou Song, Androca Gane, Tamas Sarlos, Peter Hawkins, Jared Davis, Afroz Mohiuddin, Lukasz Kaiser, et al. Rethinking atlention with performers. Xri preprind Xrv:::2099.14 14794, 2020. URL
httpps///oo. org/10 httg/10.48555 48550/arXiv...o. 2009. .14994
Krzysztof M Choromanski, Mark Rowland, and Adrian Weller. The unrcasonable effee- liveness of structured random orthogonal embeddiings. Advances in newrul information processing systems, 30, 2017. URL https::/tpeet ptttpp::/doo..oom org/10 10. 48550/ 48550/arXiv..../////////// arxiv. 1703. 008864
Krzysatof Marcin Choromanski. Taming graph kernels with random features. In- lernational Conference on Machine Learning, pages 5964-5977. PMLR, 2023. URL
Satps://doos org/10. 10. 48550/ arrivv 2305..305000 00156
Fan R. K. Chung and Shing-TTung Yau. Coverings, heat kernels and spanning trees. Electron. J. Comb., 1999. doi: 10.37236/1 / 1444. URL httphtttt///oo.ooo httr::::/600.org//13//33333 37736//143373
Corinna Cortes, Mehryar Mohri, and Afshin Rostamizadeh. Generalization bounds for learning kernels. In Proceedings of the 27th Annual International Conference on Machine Learning (ICMI, 2010), 2010. URL atttp://ww.. ca. nyu.
Keenan Cranc, Clarisse Weischodel, and Max Wardetzky. The heat method for distance computation. Communications of the ACM, 60(11):90 99, 2017. URL https://dl.aco acm. org/do1/100 1 1145/3131280
Anirban Dasgupta, Ravi Kumar, and Tamás Sarlós. A sparse johnson: Lindenstrauss trans- form. In Leonard J. Schulman, editor, Proceedings of the 42nd ACM Symposium on The- ory of Computing. STOC 2010, Cambridge, Massachusetts. USA, 5-8 June 2010, pages 341-350. ACM, 2010. doi: 10.1145/ 00.1145/1806699.0 1806689. 1806689. 1806737
URL hattps://do http:///ooooorg/1/00 org/10. 1145/
Michael Jawson-Haggeryy Trimesh repository, 2023. URL trimesh
Inderjit S Dhillon, Yuqiang Guan, and Brian Kulis. Kernel k-means: spectral clustering and normalized cuts. In Proceedings of the tenth ACM SIGKDD international conference on Knowledge discovery and dala mining, pages 556-556, 2004. URL https://ddl acm.
1111 1145/1014052..0 1014118
Michel X. Goemans and David P. Williamson. Approximation algorithms for max-3-cut and other problems via complex semidefinite programming. J. Comput. Syst. Sci., 68 doi: 10.1016/j.jcss (oe.://:..sss.00....... 2003.07.012. URL /tttps://doi.oo http://oo...rg/1.11111//////////
(2):442 470, 2004.
j css. 20003007001.
Vladimir Ivashkin. Community graphs repository, 2023. URL https://github..co vlivashkinn / community-ggraphhe
William B Johnson. Extensions of lipschitz mappings into a hilbert space, Contemp. Math.,
26:189 206, 1984. URL http://stanford. www.efodddddddddeeeeeeeeeeeeeeee cwwss/caasll///// /aaaaa//sss////aaaaaa//L-Looooe pdf.
Kyle Kloster and David F Gleich. Heat kernel based community detection. In Proceedings of the 20th ACM SIGKDD international conference on Knowledge discovery and data mining, pages 1386 1395, 2014. URL. httpp:/tts://aoioooo org/ 10. 48550/ arXiv. 1408.3148
Vladimir Koltchinskii and Dmitry Panchenko. Empirical margin distributions and bounding the generalization error of combined classifiers. The Annals of Statistics, 30(1):1-50, 2002. URL ptttps:///ot. org/10 10. 48550/ 38550/arXiv www.mattthpt... /4th/4445333)
Risi Kondor and John D. Lafferty. Diffusion kernels on graphs and other discrete input spaces. In Claude Sammut and Achim G. Hoffmann, editors, Machine Learning, Proceed- ings of the Nineteenth International Conference (ICML 2002), University of New South Wales, Sydney, Australia, July 8 12, 2002, pages 315-322, Morgan Kaufmann, 2002, URL -on.oeeeeelll..eee
httpa:/ h/ppa:///www.. ml. http:///ww....aectt//////pppe//a

---

Published conference paper at ICL. 2024
Leonid Kontorovich, Corinna Cortes, and Mchryar Mohri. Kernel methods for learning languages. Theor. Compul. Sci., 40:(3):223 236, 2008. doi: 10.1016/j..cs.00000 200...00000.............
URL https://doo: http//t....ggggg//// org/10..//.. 1016/ tcs. 2008. 06.037.
Ali Rahimi and Benjamin Recht. Random features for large-scale kernel machines. Advonces in neural information processing systems, 20, 2007. URL htttp://dd. org/doi/10..
2981710.
Isaac Reid, Krzysztof Choromanski, and Adrian Weller. Quasi-monte carlo graph random 2023. URL https://doi. org/10. 488500
features. arXiv preprint
arXiv. 1247.12470
Alexander Smola and Risi Kondor. Kernels and regularization on graphs. In Learning The- ory and Kernel Machines: 16th Annual Conference on Learning Theory and 7th Kernel Workshop, COLT/Kernel 2003, Washington, DC, USA, August 24-27, 2003. Proceed- ings, pages 144 158. Springer, 2003. URL hittps://peopps cs. www.ccccccc....oo www.ccaa...uuuuaeeeeeee edu/-risi/ Sapers/SSolaaoondoo .df.
Alexander J. Smola and Bernhard Schölkopf. Bayesian kernel methods. In Shahar Mendel- and Alexander J. Smola, editors, Advanced Lectures on Machine Leorning, Machine Learning Summer School 2002, Canberro, Australia, February 11-22, 2002, Revised Lec- lures, volume 2600 Lectare Notes in Compuler Science, pages 65 157.17.2 Springer, 2002. URL http://ppg.///.o.o////////00/2/3 1110//0000000----------------0
Yasutoshi Yajima. One-class support vector machines for recommendation tasks. In Pacific- Springer,
Asia Conference Knowledge Discovery and Dato Mining, pages 2006. URL https://dll www.daa.cc...roooooo..ooo http:///a.aa..rrr//////./.....00 /ttp://000//..000//73319999
Felix Ximnan X Yu, Ananda Theertha Suresh, Krzysztof M Choromanski, Daniel N Holtmann-Rice, and Sanjiv Kumar. Orthogonal random features. Advances newrdl in- formation processing systems, 29, 2016. URL httph:tt////e org/ 4.00.4855//aaa att.://axvvv..110000
09072
Yufan Zhou, Changyou Chen, and Jinhui Xu. Learning manifold implicitly via explicit heat- kernel learning. Advances in Newral Information Processing Systems, 33:477-487, 2020. URL http/:ttt:///t org/10. 48550/855K5i 2010.01761
12

---

Published conference paper at ICI. 2024
A APPENDICES
A.1 MINIMUM BATCH SIZE SCALES LOGARITHIICALLY WITH NUMBER OF WALKS
Consider an ensemble of m random walks S = (wi]m=) whose lengths are sampled from a
geometric distribution with termination probability Trivially, Pr(Ien(wx < 6) = - (1
p)º. Given m independent walkers,
Pr (Ilex)))s(een(ss) < b) Pr (e(())
<) = -1-(11 (1-))"²
(17)
Take 'with high probability" to mean with probability at least 1 1-6 6, with < fixed. Then
log(1 - (1 )) log(d) - log(m)
(18)
log p)
log(1 - p)
As the number of walkers m grows, the minimum value of b to ensure that all walkers are shorter than b with high probability scales logarithmically with m.
A.2 PROOF OF THEOREM 2.1 (UNBIASED APPROXIMATION OF Ka VIA CONVOLUTIONS)
We begin by proving that g- G-GF constructed according to Alg, [ give unbiased estimation of graph functions with Taylor coefficients (ax) provided the discrete convolution relation Eq. is fulfilled.
Denote the set of m walks sampled out of node by (5(())
|==1- Carefully considering Alg. I it is straightforward to convince oneself that the g-GItF' estimator takes the form
Σ
Mawio E 52())
(19)
p(w)))
Here, Skt denotes the set of all graph walks between nodes and v, of which each wie a member. len(ww) is function that gives the length of walk Who and (wtr)) evaluates the products of edge weights it traverses. p(wa)) = (1 - p) p)lmt())))) [[ T. denotes the walk's marginal probability. I is the indicator function which evaluates to if its argument is true (namely, walk Wae is prefix subwalk of ((9) the kth walk sampled out of i) and 0 otherwise. Trivially, we have that [Maso E (2(2)] p(ute) (by construction to make the estimator unbiased) Now note that:
E To01)
EE Σ 2i H07
LVEN
Σ Σ
( ((le)(((((((())))))))))
Win Elty Eyg
20
x W wh (44)((2)
(20)
1==0(==0
x
wh W's Wh (la))(1a)
1==0==
8
Σ Wh Σ F( 1==0
"To flag a subtle point: in the sum we take wir E Qi to mean that wi is a member of the set of all possible walks of any length between nodes and v, P... On the other hand, inside the indicator function by wis (( we mean that wio is a prefix subwalk of one particular walk 5(0) sampled of node i. In these two cases the interpretation of the symbol should be slightly different.
13

---

Published conference paper at ICL. 2024
From the first the second line we used the definition of GRFs Eq. 19 the sum of the product of edge weights over all possible paths as powers of the weighted adjacency matrix W, with 4,2 corresponding to walk lengths. To get to the fourth line, we changed indices in the infinite sums, then the final line followed simply.
We then rewrote
The final expression Eq. have that
20 is exactly equal to Ka = En=oan Wh provided we
h
Σ 114-4))(4)) = on:
(21)
1==0
as stated in Eq. the main text (with variables renamed to lt and p).
A.3 PROOF OF THEOREM 2.2
(COMPUTING SYMMETRIC MODULATION FUNCTIONS)
Here, we show how compute under the constraint that the modulation functions are identical, fi = f2 = We will use the relationship in Eq. g reproduced below for the reader's convenience:
Σ f(i- = a.. 1=
(22)
The iterative form in Eq. is easy to show. For i = have that f(0)² = og, so f(0) Vai #l (where we used the normalisation condition on 1). Now note that
Σ P=0
2f(0)f(i+1)
ait1
(23)
so clearly
f(i+1) 1) 0++1 p= f(i p)f(p+11
(24)
2f(
This enables us to compute f(i+ 1) given 0441 and (f(p))
The analytic form in Eq,
5 is only little harder. Inserting the discrete convolution relation
in Eq. back into Eq. 2 have that
(o(W) K, (W)
(25)
where ( = E W is the generating function corresponding the sequence (f) Constraining fi = fz.
= + ( (K.
(26)
We also discuss this the "gencrating functions' scction of Sec. 2 where we use it derive simple closed forms for f for some special kernels. Now we have that
Σ foo =
Σ
aaa WD
= (1+ W + 2
(a +..))
(27)
(a W + +...)"
n=0
We need to equate powers of W between the generating functions. Consider the terms proportional W Clearly no such terms will feature when n > i, so we can restrict the sum on the RHS to 0 S & i. Meanwhile, the term in (as W + a2 W + +...) proportional to W is nothing other than
Σ Any Boggttgg.....
bs ksnks.
(28)
kikzka....
14

---

Published conference paper at ICLR 2024
where (2nk2t the multinomial coefficient, Combining,
(N n
(Miket k112ks...)
(()) =
2
Σ
(afiozaag
ka
(29)
shown in Eq. 5 Though this expression gives f purely in terms of the presence of the conditional sum limits practical utility compared to the iterative form.
A.4 PROOF OF THEOREM 2.3 (EMPIRICAL RADEMACHER COMPLEXITY BOUND)
In this appendix derive the bound the empirical Rademacher complexity stated in Theorem 2.3 and show the consequent generalisation bounds. The early stages closely follow arguments made by Cortes et al. (2010). Recall that we have defined the hypothesis set
H (c w 4K. (2) |asl a
||wl|a 1}.
(30)
with v0,, x + K., the feature mapping from the input space to the reproducing kernel Hilbert space Hk,, induced by the kernel Kg.
The empirical Rademacher complexity R(H) for an arbilrary fixed sample $ = (m)) is defined by
R(H)
- Er m
sup
(31)
keHf
the expectation taken over (01, -0 Fm) with 03 Unif(+11) i.d. Rademacher random variables.
Begin by noting that
M(x) = w uk (x) = Σ (,K,(
(32)
where 3 are the coordinates of the orthogonal projections of on Hs
(pan, (Z))..., thk (Im)), where 3 K.B s 1. Then we have that
1 E. sup o K R(H)
(33)
The supremum supg or K.3 is reached when K B is collinear with K/ a, and making || as large as possible gives
R(H) E. sup o Kod
(34)
Now note that
o K.o=
Σ ao W'
(35)
W'a may take either sign, and the sum maxinised by laking or = a, for positive terms and O; -a (M) for negative terms. Observe that
supsup V KoT
m i=0
a
(M), P(W)
(36)
whereupon from Eq. 31 it follows that
8
R(H) &
(M) pl
(W)
(37)
i=0

---

Published conference paper at ICLR 2024
as stated in Thm 2.3 This bound not tight for general graphs, but will be for specific ex- amples: for example, when W is proportional to the identity so the only edges are self-loops. Nonetheleess, it provides some intuition for how the learned kernel's complexity depends on W.
As stated the main text, Eq. 37 immediately yields generalisation bounds for learning kernels. Again closely following Cortes et al. 2010), consider the application of binary clas- sification where nodes assigned label yi EI. Denote by R(k) the true generalisation error of E H;
R(h) = Pr(yh(x) < 0).
(38)
Consider training sample S ((x4, Vs (s)) and define the p-empirical margin loss by
(h) = (in uh((x./pp.
(39)
where p > 0. For any 8 > 0, with probability at least 8, the following bound holds for 2002 Koltchinskii and Panchenko 2002):
any h E H |Barlett and Mendelson
R(h) & (b) 2 RR(H) +3
log
(40)
2m
Inserting the bound on the empirical Rademacher complexity in Eq. 37 have that
immediately
R(h) - (b) +
2
(M) P
(W)
+31
log
(41)
2m
which shows how the generalisation error can be controlled via (M)) or p(W).
A.5 LEARNING ASYMMETRIC ((N)
Recalling from Eq. that the modulation func- tions (f1, f2) do not necessarily need to be equal for unbiased estimation of some kernel Ka natural extension of Sec. 3.4 to rndrud two separate neural modulation functions fi, and train them both following the scheme in Sec. 3.4 Intriguingly, even with an initialisa- tion where f2 encodes 'lazy' behaviour (de- posits almost all its load at the starting node see Sec. 2 and f(a) flat, upon training the neural modulation functions quickly become very similar (though not identical). Seo Fig. 0 We use the samne optimisation hyperparaneterr and network architectures as in Sec. 3.4 Rigor- ously proving the best possible choice of (fi. f2).
Learned asymmetric modulsti functions
intllds trained
14-
Figure 6: Modulation functions (f1, f2) parameterised by separate neural net- works before and after training target the 2-regularised Laplacian kernel.
including whether e.g. symmetric pair is opti-
mal, is left as an exciting open theoretical prob- lem. We note that, whilst parameterising two separate neural modulation functions gives a more general mechanism (with the symmetaic pair as a special case), also doubles the mumber of parameters required so slows training and evaluation.
A.6 FURTHER APPROXIMATION ERROR RESULTS FOR ERDÖS-RÉNYI GRAPHS
Here we further investigate the behaviour of the GRFs kernel approximation as the prop- erties of the graph change. To do this, we generate random graphs using the Erdös-Réngt model, where each of the (N) possible edges are present with probability Pedge or absent with probability Podge- Edges independent.
16

---

Published conference paper at ICIR 2024
Firstly, we investigate how the approximation error varies with graph sparsity. We take a fixed number of nodes (N = 100) and control the sparsily by varying Pedge between 0.1 and 0.9. We then approximate the diffusion kernel using g-GRP' with lermination probability Pbalt 0.1 and walkers. For each graph we compute the relative Frobenius norm error between the true (K) and approximated (K) kernels: namely, ||- || We repeat 100 times to obtain the standard deviation of the mean error estimate. Fig. 7 shows the results; approximation quality degrades slightly as Pedge grows, but remains very good throughout (error 0.00275).
Next, we investigate the scalability of the method by comparing the time for exact and approximate evaluation of the diffusion kernel as the size of the graph grows. We fix Pedge 0.5 and vary the number of nodes N between 100 and 12800. For every graph, measure the wall-clock time for exact computation of K by computing the matrix exponential and ii) approximate computation of K using g- g-GREs (8 walkers, Phah = 0.5, regulariscr a 0.5). Fig. [ shows the results. Naturally, the exact method scales worse, becoming slower for graphs bigger than a few thousand nodes, and by the largest. graph g-GRF's are already faster by factor of 7.8. We also measure the relative Frobenius norm between the true and approximated Gram matrices to check that the quality of estimation remains good. This is shown by the green line, which indeed remains almost constant and takes very small values for the error 0.005).
Figure 7: Approximation error vs. edge- generation probability for the diffuusio ker- nel graph of size N = 100. The quality remains good as sparsity changes, varying within a narrow range.
Figure 8: Wall clock time for exact and ap- proximate kernel evaluation as the number of nodes varies. g gRGs scale better and are faster with a few thousand nodes. The ap- proximation error remains low as N grows.
Apgroximation error edge probability
Wall-clock time grapp
100000
exatt ex. time
approximation enor
angenx hime
6.3923
click time
00005
-2224--
00/50
10-1
6.30250
nodes
A.7 FURTHER EXPERIMENTAL DETAILS
In this short section, provide further experimental details and discussion to supplement
the main text.
Choice of Phote: The termination probability encodes simple trade-off between approximation quality and speed: if Push is lower, walks tend to last longer and sample more subwalks so give better approximation of graph kernels. In prac- tice, any reasonably small value works well, also reported for the original GRTss mechanism [Choromanski 2023). In experiments we typically choose Phalt ~ 0.1.
2. Choice of a and kernel convergence: After Eq. 2 we noted that, when ap- proximating some fixed kernel (W) Ex= Og W we assume that the sum does not diverge. It would not be possible to construct a finite random feature estinator if this were not the case. We require that E (gp((W)* is finite, with p( (W) == maxxex(W (lAl) the spectral radius of the weighted adjacency matrix (A(W) is the set of eigenvalues of W). When considering e.g. the diffusion kernel in the literature, this ensured by multiplying W by some regulariser R,,, tak-
17

---

Published conference paper at ICL. 2024
ing c.g. W - a W whereupon X; -i= I..,,N This the reason for extra parameters {a, ax in the kernel expressions in Table [ As we noted in Sec. 3.5 this is exactly equivalent to transforming the modulation function f(i) + [(i)a* 7 Where space allows, we report the exact choice of a with each experiment (though empirically it docs not modify our conclusions). In Sec. 3.5 we take the weighted adjacency matrix = [a.j /dd// with d;j if node i connected to j and 0 otherwise, and dj the degree of node i. We then regularise by taking W - a WW with 0.025, which is small enough for convergence even with the largest graph considered. We use m 16 walkers and Phalt 0.5.
18


Reviews:
[{"review_id": "5chWjkKCY9", "reviewer": "Reviewer_44jF", "summary": "This paper extends a result of Choromanski (2023) in order to compute efficiently an approximation of any functions of a weighted adjacency matrix. The main idea is to add a so-called modulation function in an algorithm based on random walks in order to compute the feature projector.\n\n###\nAfter a discussion with the authors, I changed my ratings about soundness and contribution and my general rating as the authors clarified my mistake about the proof of their Theorem.", "questions": "I think equation (19) in the appendix is not correct. This makes the proof of the main result Th 2.1 incorrect. Equation (19) involves only the modulation function evaluated on the total length of the walk, but in Algorithm 1 on line 8, the update for the feature vector uses all values of the modulation function evaluated at $k\\leq$ length(walk).\nYou probably mean in the indicator function that $\\omega_{iv}$ is a 'prefix' of the total walk?\nBut then, when you consider the full walk, you need to have a factor $p$ and not only $(1-p)^len$?\n\nI find the term universal quite misleading. You are approximating arbitrary functions of the weighted adjacency matrix, but these are not all possible graph kernels. Indeed, it would be nice to have the expressive power of the kernel you are approximating if possible.", "limitations": "", "rating": 8, "confidence": 5, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "This result is a nice extension of the previous result due to Choromanski (2023) allowing us to compute approximations of a variety of graph kernels.", "weaknesses": "I do not understand the proof of the main Theorem 2.1 see below.\n\nThe experiments could be improved. The authors mainly compare the results given by their algorithm to the real kernel, i.e., measure the approximation error due to their algorithm. It would be more convincing to have used their algorithm on a 'real' graph problem with a large number of nodes and where standard kernel methods are not possible.", "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": "mxHakhlwaY", "reviewer": "Reviewer_jWMX", "summary": "This paper introduces universal graph random features (u-GRFs) based on random walks in the graph. The authors show that u-GRFs can be used to unbiasedly approximate functions of the weighted adjacency matrix of a graph, which include popular graph kernels such as d-regularized Laplacian kernel, p-step random walk kernel and the heat diffusion kernel. In addition, the authors extend their algorithm for estimating fixed graph kernels to a learnable approach, parametrizing a fixed modulation function with a two-layer neural network. The authors provide a simple generalization bound for the corresponding class of learnable graph kernels. Finally, extensive experiments over various application settings are provided to demonstrate the good performance and potentially wide applicability of u-GRFs.", "questions": "- Given Theorem 2.1 and Equation 4, as long as f_1 and f_2 satisfy Equation 4 we will have an unbiased estimator. In that case, what is particularly interesting/special about symmetric modulation functions? Does using symmetric modulation functions reduce estimator variance (either empirically or theoretically)?\n\n- Why did you use the pairwise metric in Equation 15 as opposed to node-wise accuracy compared with the clustering result using the exact kernel?\n\n- If I understood correctly, the numbers in the parenthesis in Table 3 and Table 4 are standard deviations. Why are they integers?", "limitations": "", "rating": 8, "confidence": 4, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "- This paper extends the recent work of Choromanski (2023) on random walk-based graph random features. Compared to the previous work, this paper has two notable innovations. The first is an introduction of a modulation function which enables estimation of arbitrary functions of the weighted adjacency matrix, the second is an extension to a parametrized setting with a generalization bound. While this work is not the first to study graph random features, it has an easy-to-identify contribution which I think makes it a nice addition to the community.\n\n- The overall presentation and writing are clear and easy to follow. I enjoyed reading this paper.\n\n- The experiments, though limited to a few small graphs, are reasonably extensive and showcase the potentially wide applicability of the proposed universal graph random features.", "weaknesses": "- The experiments are carried over small graphs. It would be nice to include some experiments on larger graphs, e.g. N=10,000. For example, fix the Erdos-Renyi model but vary the number of nodes. I wonder if and how the approximation error would scale with N, while keeping everything else fixed.\n\n- For experiments on node clustering, Table 2 shows the difference between u-GRFs and the exact kernel. It would be nice to have another table that shows the actual clustering accuracy when compared against the ground-truth information, given the ground-truth number of clusters in each graph.", "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": "a91hVYWNRz", "reviewer": "Reviewer_9UzR", "summary": "A large span of interesting problems involving graph data require to compute a function of the weighted adjacency matrix. In most cases, using a naive procedure to compute this function scales (at least) cubicly with the number of nodes in the graph. \n\nThe authors propose a new algorithm to circumvent this computational limitation. Inspired by a recent work by Krzysztof Marcin Choromanski (KMC), the authors introduce a new algorithm to approximate the value of arbitrary functions of the weighted matrix of the graph by using graph random features. Graph random features rely on the possibility to write the function of the weighted adjacency matrix as the expectation of finite dimensional vectors. The former vectors are obtained by using random walks on the graph. Contrary to the KMC's work, the author introduce the notion of modulation function that controls each walker’s load as it traverses the graph. This extra modulation function allows the algorithm to approximate arbitrary functions of the weighted adjacency matrix.\n\nThe authors also show that this modulation function can be learned using neural network to solve some given task of interest.", "questions": "I thank the authors for their nice work. Here are my questions and comments.\n\n- The authors did a good job in providing several experiments to illustrate their method. Nevertheless, I think it would have been good to propose an experiment on real data that tackles one of the main applications mentioned in the introduction (such as community detection or recommender systems). \n\nIn particular, it is still an ongoing research question to find the right representation of a graph to use spectral methods for community detection. Indeed, using the adjacency matrix is known to be suboptimal in particular when dealing with sparse graphs. That's why people have proposed to use the non-backtracking matrix or the Bethe-Hessian. Since the authors show the interesting possibility to learn the modularity function to automatically selects a graph representation well suited for some downstream task, I think that a very interesting experiment would be to use their method to automatically select a good kernel for community detection. In particular, does the algorithm converge to a representation close to the non backtracking matrix ?\n\n\n- The authors do not discuss really how some parameters of the method are chosen in practice (such as the parameter $p_{halt}$), and how these choices influence the performance of the algorithm.\n\n- Real world graphs are known to be sparse. It would be interesting to comment of the influence on the algorithm when dealing with sparse graphs. A possibility would be to conduct simulated experiments by varying the sparsity level of the graph and see the impact of the performance.\n\nHere are some typos and additional comments: \n\n- In Eq.(9), I think a notation is missing since $\\alpha$ is not appearing in the mapping $x\\mapsto w^{\\top} \\psi_{K}(x)$. I think the correction would be $x\\mapsto w^{\\top} \\psi_{K_{\\alpha}}(x)$.\n\n- At the first line in section 3.1, I think a \"for\" is missing. It should be \"do indeed give unbiased estimates for the graph kernels\".\n\n- Between Eq.(13) and Eq.(14), I think there is an issue with the notations. Should $\\Phi$ depend on $\\tau_j$ ? Based on what I could understand, $\\Phi$ is a rectangular matrix of size $n\\times N$ where the $j$-th row is $\\Phi_j = ((t-\\tau_j) \\phi(i))_{i=1}^N$.", "limitations": "", "rating": 8, "confidence": 4, "soundness": 2, "presentation": 3, "contribution": 3, "strengths": "- The paper brings two main contributions. First the introduction of the notion of modulation function and its use in a random walk algorithm to estimate arbitrary functions of the weighted adjacency matrix of a graph. Second, the possibility to learn using neural network this modularity function to automatically selects a graph representation well suited for some downstream task. \n\nTherefore, this work not only allows to reduce the computational burden to solve some important graph related questions using graph random features and the proposed algorithm, but the method can be used to automatically find the relevant graph representation for a task of interest.\n\n- The paper includes both interesting theoretical results and several numerical experiments (such as node attribute prediction on triangular mesh graphs).\n\n- The paper is well written and the authors pan to make their code publicly available in case their work is accepted.", "weaknesses": "- The authors highlight several times in the paper that one major issue regarding in terms of computational complexity is to inverse the matrix K(W) (where W is the weighted adjacency matrix of the graph). However, in the numerical experiments presented in the paper, the authors do not propose an example of application of their method to inverse such matrix K. It might be relevant to add such experiment in the paper or at least to explain how the method can be used to efficiently compute the inverse.\n\n- My other comments or suggestions for improvements are mainly related to the experiment section and I refer to the list of my questions.", "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": "KXeR3ra5XV", "reviewer": "Reviewer_TqUg", "summary": "The paper introduces techniques to approximate a family of kernels for comparing the nodes of a graph. The methods rely on modulation functions that are either analytically derived from the original kernel, or learned from data by optimizing the downstream task.\nExperiments show that as the representation grows the performance improves and that the learned kernels perform well in inductive tasks too.", "questions": "- Alg. 1: shouldn't the normalization be $1/\\sqrt{m}$ instead of $1/m$?\n- Fig. 3: please specify if the extremely narrow shaded area around the lines refers to +/- one standard deviation.\n- Sec. 3.4: it is claimed that the modulation function learned on one data set can be effectively applied to other graphs, referring to Tab. 3. However, I don't see specified how the table has been generated. Please provide more details.\n- Tab. 3 and 4: I suggest introducing the notation used for the standard deviation.\n- I suggest discussing the impact of the fact that the sum in Eq. 2 is not always well-defined - an aspect that is only briefly mentioned in the paper.\n- I suggest discussing the pros and cons of symmetric vs asymmetric approaches.", "limitations": "", "rating": 8, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 2, "strengths": "- The proposed method allows approximating the entire family of kernels based on weighted adjacency matrices.\n- Learning the modulation function shows how to successfully retrieve effective kernels within the considered family.\n- Examples show how to apply the proposed methods in different scenarios and tasks.\n- The paper reads well, and the contribution w.r.t. prior work is fair and clearly stated.", "weaknesses": "- The proposed method is an extension of previous works, in particular, Choromansky (2023).\n- The experiments do not clearly show the price that one has to pay in terms of compute time when requesting better accuracy. Experiments show only relations between the number of walks and the accuracy/error. Linking accuracy and execution time (e.g., in standard machines) would have given a better sense of the relevance of the proposed methods.\n- The kernels apply only to weighted graphs without node and edge attributes.\n- The use of the term 'universal' is therefore misleading; as the authors say, the provided random features can approximate only kernels from a specific family.", "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)