Quasi-Monte Carlo Graph Random Features
Authors: Isaac Reid
Krzysztof Marcin Choromanski
Adrian Weller
Venue: NeurIPS 2023 spotlight (2023)
Keywords: Graph
discrete mathematics
quasi-Monte Carlo
kernel
scalability
Laplacian
clustering
random walks
TL;DR: A novel QMC mechanism that induces correlations between the terminations of random walkers on a graph to improve kernel estimator variance.

Abstract: We present a novel mechanism to improve the accuracy of the recently-introduced class of graph random features (GRFs). Our method induces negative correlations between the lengths of the algorithm's random walks by imposing antithetic termination: a procedure to sample more diverse random walks which may be of independent interest. It has a trivial drop-in implementation. We derive strong theoretical guarantees on the properties of these quasi-Monte Carlo GRFs (q-GRFs), proving that they yield lower-variance estimators of the $2$-regularised Laplacian kernel under mild conditions. Remarkably, our results hold for any graph topology. We demonstrate empirical accuracy improvements on a variety of tasks including a new practical application: time-efficient approximation of the graph diffusion process. To our knowledge, q-GRFs constitute the first rigorously studied quasi-Monte Carlo scheme for kernels defined on combinatorial objects, inviting new research on correlations between graph random walks.

Quasi-Monte Carlo Graph Random Features
Isaac Reid University of Cambridge ir3370cam.ar ac uk
Krzysztof Choromanski
Adrian Weller University of Cambridge Alan Turing Institute aw6650cas camm. ac. uk
Google DeepMind
Columbia University
kcchroogoogge... com
Abstract
We present novel mechanism improve the accuracy of the recently-introducee class of graph random features (GRFs) [Choromanski 2023 Our method induces negative correlations between the lengths of the algorithm's random walks by imposing antithetic fermination: procedure to sample more diverse random walks which may be independent interest. has a trivial drop-in implementation. We derive strong theoretical guarantees on the properties of these quasi-Monte Carlo GRFs (q-GRFs), proving that they yield lower-variance estimators the 2-regularised Laplacian kernel under mild conditions. Remarkably, results hold for any graph topology. We demonstrate empirical accuracy improvements on a variety of tasks including the graph diffusion process. To our knowledge, q-GRFs constitute the first rigorously studied quasi-Monte Carlo scheme for kernels defined combinatorial objects, inviting new research on correlations between graph random walks [
practical application: time-efficient approximation
Introduction and related work
Kernel methods are ubiquitous in machine learning (Canu and Smola 2006 Smola and Schölkopf 2002 Kontorovich et al. 2008 Campbell 2002). Via the kernel trick, they provide mathematically principled and elegant way to perform nonlinear inference using linear learning algorithms. The positive definite kernel fianttion X' X IR, defined on input domain X', measures the "similarity' between two datapoints. Examples Euclidean space include the Gaussian, linear, Matérn, angular and arc-cosine kernels (Williams and Rasmussen 2006 Cho and Saul 2011].
Though very effective on small datasets, kernel methods suffer from poor scalability. The need to materialise and invert the kernel matrix typically leads to time-complexity cubic in the size of the dataset. Substantial research has been dedicated to improving scalability by approximating this matrix. notably including random features (RFs) [Rahimi and Recht] [007 2008 [Avron et al.) 2017
Liu et al. [0222] These randomised mappings o R + R construct low-dimensional feature vectors whose dot product equals the kernel evaluation in expectation:
k(æ, y) E (Φ(x) To()))
This permits low-rank decomposition of the kernel matrix which enables better time- and space- complexity than exact kernel methods. Random feature methods exist for a variety of kernels with properties engineered for the desiderata and symmetries of the particular kernel being
Euclidean
"Senior lead.
'Source code available atthgs///ggthhbbccmmmmscrreealaa http:///weee.eeeaiitiitieeeeeeen termination
37th Conference on Neural Information Processing Systems (NcurIPS 2023).

---

approximated Dasgupta et al. 2010 Johnson 1984 Williamson 2001
Choromanski et al. 2020 Goemans and
Rahimi and Recht 2007
Kernels can also be defined on discrete input spaces such as graphs, which are the natural way to represent data characterised by local relationships (e.g. social networks or interacting chemicals [Albert and Barabási, 2002)) or when data is restricted to a lower-dimensional manifold than the original space [Roweis and Saul, 2000 Belkin and Niyogi 2003]. We consider graph kernels k: N x N + I on the set of nodes N of a graph G. Examples include the diffusion, regularised
Laplacian, p-step random walk and cosine kernels [Smola and Kondor 2003 Kondor and Lafferty
2002 Chung and Yau 1999], which have found applications including in bioinformatics [Borgwardt al. 2005]. community detection [Kloster and Gleich| 2014) and recommender systems (Yajmma
2006). More recently. these kernels have been used in mamifold learning for deep generative modelling
[Zhou et al. [020] and for solving shortest path problems [Crane et al., 2017 7 Substantial research
effort has also been devoted to developing and analysing graph kernels k: G x g + IR, now taking entire graphs G E 9 from graph spaces G as inputs rather than their nodes (Shervashidze et al.) 2009 Vishwanathan et al. 2006 Shervashidze and Borgwardt 2009], but we stress that these are not the subject of this paper.
The problem of poor kernel scalability is exacerbated in the graph domain because even computing the kernel matrix K is typically of least cubic complexxy in the number of nodes N. In contrast to kernels defined on points in Rd. random feature methods for fixed graph kernels (c.f. kernel learning (Fang et al. [ have proved challenging to construct. Only recently has a viable graph random feature ( (GRF) mechanism been proposed, which uses series of random walkers depositing 'load' at each node they pass through [Choromanski] 2023). GRFs provide an unbiased estimate of the matrix (IN~ U)- where U is weighted adjacency matrix of the graph and d e N. The decomposition supports subquadratic time-complexity (again with respect to the number of nodes N) in downstream algorithms applying regularised Laplacian kernels. Moreover, the computation of GRFs admits a simple distributed algorithm that can be applied if a large graph needs to be split across machines. The author demonstrates the strong empirical performance of GRFs in speed tests, Frobenius relative error analysis and a k-means graph clustering task.
In the Euclidean setting, significant research has been dedicated to developing quasi-Monte Carlo (QMC) variants to RF methods that enjoy better convergence properties (Yang et al. 2014 Lyu 2017 Dick et al.] 2013). By using correlated ensembles rather than i.i.d. random variables in the feature maps, one can suppress the mean squared error (MSE) of the kernel estimator. For example, orthogomal random features (ORFs) (Yu et al.) 2016 Choromanski et al.) 2020) improve the quality of approximation of the Gaussian kernel when using trigonometric or positive random features, and of the linear kereel in the orthogonal 2017]. With positive random features, the recently-introducee class of simplex random features (SimRFs) performs even better [Reid et 20.3). J. This has been used to great effect in estimating the attention mechanism of Transformers |Vaswani et al. [017]. overcoming its prohibitive quadratic time-complexity scaling with token sequence length.
transformation [Choromanski et al.,
This invites the central question of this work: how can we implement a QMC mechanism for random walks on a graph? What do we mean by a 'diverse' sample in this context? Choromanski [2033] first identified the challenge of constructing quasi-Monte Carlo GRFs (R-GRss. They suggested high-level approach of reinforced random walks, but left to future work its theoretical and empirical analysis. In this paper, we provide a first concrete implementation of q-GRFs, proposing an unbiased scheme that correlates the length of random walks by imposing antithetic termination. We derive strong theoretical guarantees on the properties of this new class, proving that the correlations reduce the variance of estimators of the 2-regularised Laplacian kernel under mild conditions. Our results hold for any graph topology. We also demonstrate empirical accuracy improvements on a variety of tasks. We hope our new algorithm (hereafter referred to 'q-GRFs' for brevity) will spur further research on correlations between graph random walks in machine learning.
We emphasise that, although we present antithetic termination through the lens of q-GRFs (an analytically tractable and important use case), is fundamentally a procedure to obtain a more diverse ensemble of graph random walks and may be of independent interest. To illustrate this, in App 8. we show how it can alternatively be used to improve approximations the PageRank vector [Page et al. 1998 It could also be relevant e.g. for estimating graphlet statistics for kernels between graphs
2

---

Chenn et al. 2016 Ribeiro et al. 2021 or in some GNN architectures Nikolentzos and Vazirgiannis 2020
The remainder of the manuscript is organised as follows. In Sec. 2 we introduce the mathematical concepts and existing algorithms to be used in the paper, including the d-regularised Laplacian and diffusion kernels and the GRF mechanism. Sec. 3 presents our novel q-GRFs mechanism and discusses its strong theoretical guarantees in particular, that provides lower kernel estimator variance than its regular predecessor (GRFs) under mild conditions. We provide a brief proof-sketch
for intuition but defer full technical details to App. 8.4 We conduct an exhaustive set of experiments in Sec. to compare q GRFs to GRFs, including: (a) quality of kernel approximation via computation of the relative Frobenius norm; (b) simulation of the graph diffusion process; (c) kernelised k-means node clustering; and (d) kernel regression for node attribute prediction. q-GRFs nearly always perform better and some applications the difference is substantial.
2 Graph kernels and GRFs
2.1 The Laplacian, heat kernels and diffusion on graphs
An undirected, weighted graph G(N, E) is defined by set of vertices W = (l, N) and a set of edges & given by the unordered vertex pairs (i.j) E E where i and j are neighbours, themselves associated with weights Way R. The weighted adjacency matrix W RNx has matrix elements Wy: that is, the associated edge weights if (i,j) E E and 0 otherwise. Let (ii) = {j E (i. e E} denote the set of neighbours of node
Denote by D E N the diagonal matrix with elements Du == Ej Wij. the sum of edge weights connecting a vertex i to its neighbours. The Laplacian of G is then defined L == D W. The normalised Laplacian is L == D -- LDDLDD 1, which rescales L by the (weighted) number of edges per node. L and L share eigenvectors in the case of an equal-weights d-regular graph and play a central role spectral graph theory: their analytic properties are well-understood [Chung 1997..
In classical physics, diffusion through continuous media is described by the equation
du v'u
(2)
dt
where V
8² is the Laplacian operator continuous spaces. The natural 0m3 Day analogue discrete spaces is L, where can be seen by noting that, if we take W to be the unweighted adjacency matrix, (u, Lu) = u Lu = )² so, like its continuous counterpart, L measures the local smoothness of its domain. In fact, in this case L is exactly the finite difference discretisation of 2 on square grid in N-dimensional Euclidean space. This motivates the discrete heat equation on G,
02 82
au?
treat L as linear operator on vectors RN This
- "(Cj)
du dt Lu,,
where we followed literature conventions by using the normalised variant whose spectrum is con- veniently contained in [0,2] Chung 1997]. This has the solution u, =ex((-L))mm where exp(-L/) L/)= = -im
I.. The symmetric and positive semi-definite matrix
(tn == exp(- Lt)
(4)
is referred to the heat kernel or diffusion kernel (Smola and Kondor 2003 The exponentiation
of the generator L, which by construction captures the local structure of G, leads to kernel matrix Kann which captures the graph's 's global structure. Upon discretisation of Eq. 3 with the backward Fuler step (which is generally more stable than the forward), we have that
Up+-6t (IN + 8tL) ues
(5)
where the discrete time-evolution operator K (1) = (IN + stī) is referred to the 1-regularised
Laplacian kernel. This is member of the more general family of d-regularised Laplacian kemnlss
K
(d) lap
(IN + StL)
(6)
3

---

for which we can construct an unbiased approximation using GRFs [Choromanski] 2023]. We predominantly consider d = with the understanding that estimators for other values (and the diffusion kernel) are straightforward to obtain see App. 8.2 This demonstrates the intimate connection between the graph-diffusion and Laplacian kernels, and how a QMC scheme that improves the convergence of K (2(2) will be of broad interest.
2.2 Graph random features (GRFs)
Here we recall the GRF mechanism, which offers a rich playground for our novel antithetic termination QMC scheme. The reader should consult Choromanski [2023] (especially Algorithm 1.1) for full discussion, but for convenience we provide cursory summary.
Suppose we would like to estimate the matrix (Iy -Iv-U) with U RNN a weighted adjacency matrix of a graph with N nodes and no loops. Choromanski 2023] proposed a novel algorithm to
construct a randomised estimator of this matrix. The author uses graph random features (GRFs)
((1) E R with the index of one of the nodes, designed such that
(IN -
(+(i) T(1)).
(7)
They construct (b)) by taking m random walks {$2(k; i) )= on the graph out of node i, depositing
a 'load' every node that depends on i) the product of edge weights traversed by the subwalk and ii) the marginal probability of the subwalk. Importantly, each walk terminates with probability at every timestep. The æth component the ith random feature is given by
1 Σ Σ
((w1) M(w1 E (((,,))) p(w1)
<(i)= =
(8)
Wi
where: l enumerates mn random walks we sample out of node w1 denotes a particular walk from the set of all walks Du between nodes and z; ((w1) is the product of weights of the edges traversed by the walk wi: p(w1) the marginal probability that a walk contains a prefix subwalk w1, given by p(w1) = ((1 p)ddd)//ddda /lm(e))))) in the simplest case of walk of length len(wi) on d-regular graph; I[wa E (2kks i)) is an indicator function that evaluates to 1 when the walk w1 is a prefix subwalk of the kith random walk sampled from i (itself denoted (2k, i)) and is 0 otherwise.
is simple to see how Eq. satisfies Eq. 7 By construction
E [[(n E 12Kks ()) E (2(1,))) p(w))p(w2))
(9)
for independent walks, whereupon
E ( = ΣΣΣ Σ
(Ww)) =
Σ (len(w) + 1) = (I - U)
-I-U),²
DEN we EQte w2El
wElly
(10)
This shows us that the estimator is unbiased. The central contribution of this work is a QMC scheme that induces correlations between the walks out of each node to suppress the variance of the estimator (()) (()) without breaking this unbiasedness.
3 q-GRFs and antithetic termination
We will now present our novel antithetic termination mechanism. It generalises the notion of antithetic variates a common, computationally cheap variance-reduction technique when sampling in Euclidean space [Hammersley and Morton 1956) to the termination behaviour of random walks.
We have seen that, in the i.i.d. implementation of the GRF algorithm, each walker terminates independently with probability p at every timestep. For a pair of i.i.d. walkers out of node i, this is implemented by independently sampling two fermination random variables (TRVs) between 0 and 1 from a uniform distribution, t1,2 Unif (0, 1). Each walker terminates if its respective TRV is less than p. t1 < p. In contrast, we define the antithetic walker as follows.
Definition 3.1 (Antithetic walkers). We refer to a pair of walkers as antithetic if their TRVs are marginally distributed 112 Unif(0, 1) but are offset by 2.
t2 = modi
t1 + 2
(11)

---

00
1
((8)
((1)
i.i.d.
antithetic
(a)
(b)
Figure I: Left: schematic of the i.i.d. (GRF) and antithetic (q-GRF) mechanisms in termination space. [() is the probability density of the termination random variable (TRV) L. Vertical arrows represent draws of with the walker terminating if they lie in the pink region where f p. With q-GRFs the TRVs are offset by 1. modifying the joint distribution over walk lengths. Right: demonstration with 4 random walks on the karate graph, beginning some node labelled The blue pair of antithetic walks (q-GRFs) have very different lengths; they cannot terminate simultaneously. The red pair of i..d walks (GRFs) have similar lengths. We prove that q- q-GRFs give lower variance estimators the 2 2-regularised Laplacian kernel.
such that have the conditional probability density
p(t2lt1) = (modd - -1)
(12)
In both schemes, the TRVs are resampled at every timestep until the corresponding walker terminates.
Computational cost: We note that the computational cost of generating antithetic TRV according to Eq. iis is no greater than the cost of generating an independent TRV, so this drop-in replacement in the GRF algorithm is cheap. We provide schematic in Fig. 1a
Since the marginal distributions over t; are unchanged our estimator remains unbiased, but the couplings between TRVs lead statistical correlations between the walkers' terminations, Denoting by 81 the event that walker terminates some timestep, 82 event that walker 2 terminates and $1.2 their complements, it is straightforward to convince oneself that for
p(s)) p(22) = P. p(5)) ((52) = P.
p(52 (1) 0,
1 1 2p
(13)
(( 1, ((22(51) 1-p
31 = 1 1-P
This termination coupling modifies the joint probability distribution over walk lengths. In the i.i.d. scheme, the walks are independent and are of expected length
1 1-P
[len(w)) =
(14)
These marginal expectations are preserved in the antithetic scheme, but now the expected length of one walk conditioned on the length of the other is
1 2p + (1 m (1-p)
B(len(w>) (ee(w())))))))) m)
(15)
p
which derive in App. 8.3 It straightforward to see that the two lengths are negatively correlated. Antithetic termination 'diversifies' the lengths of random walks we sample, preventing them from clustering together, and in the spirit of QMC this turns out to suppress the kernel estimator variance. See Fig. 1b for a schematic. We refer to random features constructed with antithetic walkers as Munse Carlo graph random features (q-GRFs).
Though this paper focuses on antithetic termination for GRFs, a more general QMC scheme applicable in algorithms that sample walks with geometrically distributed lengths. To illustrate this, we again direct the reader to App. 8. where we show how can used to improve numerical estimates of the PageRank vector [Page et al. 1998).
5

---

3.1 Theoretical results
In this section, we state and discuss our central theoretical results for the q-GRFs mechanism. Sec. 3.1.I provides a sketch, but full proofs are deferred to App. 8.4 We remind the reader that results for (IN- U)-² are trivially applied K (see App. 8.2
Theorem 3.2 (Antithetic termination is better than i.i.d.). For any graph, q-GRFs will give lower vari- ance on estimators of (IN - -U) than regular GRFss provided either i) the termination probability p is sufficiently small or ii) the spectral radius P(U) is sufficiently small.
By 'sufficiently small' we mean that for a fixed p(U) there exists some value of p below which antithetic termination will outperform i.i.d. Likewise, for fixed p there exists some value of p UU) These conditions turn out to not be too restrictive in our experiments; antithetic termination is actually very effective at p = 0.5 which we use for practical applications.
Considering Eq. 11 carefully, is easy see that the termination probabilities in Eq. 13 are not particular to a TRV offset equal to 2. and in fact hold for any offset A satisfying S A P. An immediate corollary is as follows.
Corollary 3.3 (Maximum size of an antithetic ensemble). For a termination probability p. up to
! random walkers can all be conditioned to exhibit mutually antithetic termination.
This is achieved by offsetting their respective TRVs by A = The resulting antithetic ensemble will have lower kernel estimator variance than the equivalent number of i...d. walkers or a set of mutually independent antithetic pairs. We make one further interesting remark.
Theorem 3.4 (Termination correlations beyond antithetic). A pair of random walkers with TRVs offset by p(1 p) < A p will exhibit lower variance on estimators of (IN U)~ than independent walkers, provided either i) the termination probability is sufficiently small or ii) the spectral radius of the weighted adjacency matrix p(U) is sufficiently small.
This provides an upper limit of L(x(1 p}}- on the number of walkers can simultaneously correlate before we can no longer guarantee that coupling in a further walker's TRV will be better than sampling independently. Intuitively, Theorem tel]tells us that we can space TRVs even more closely than p. allowing us to increase the number of simultaneously anticorrelated random walkers the cost of the strength of negative correlations between walkers with neighbouring TRVs.
3.1.1 Proof sketch
In this section, we will outline a proof strategy for the results reported earlier in this section. Full technical details are reported in App. 8.4
From its Taylor expansion, the (ih)-th element of (IN U)-² is nothing other than a sum over all possible walks between the nodes and j. weighted by their lengths and the respective products of edge weights. GRFs use Monte Carlo scheme to approximate this sum by sampling such walks at random - concretely, by first sampling separate random walks out of nodes i and j and summing contributions wherever they intersect. In order to sample the space of walks between nodes i and more efficiently, it follows that we should make our ensemble of walks out each node more diverse. q-GRT's achieve this by inducing negative correlations such that they are different lengths.
In more detail, it is clear that the variance of the kernel estimator will depend upon the expectation of the square of ((1) Φ(j). Each term in the resulting sum will take the form ((w)) (wwz ((wa) (wa)) p(wn) p(ws) p(ws) p(ws)
Σ Σ3 Σ
R. T.VENiEERL wg Wggltt dig Engy
(16)
p(wn S2(ki, ()), wa E S(kk, ())p(s E (1(,. wa E (((,,))),
where we direct the reader to Sec. 2.2 for the symbol definitions. Supposing that all edge weights are equal (an assumption we relax later), the summand is function of the length of each of the walks w1,2.3,4 It is then natural to write the sum over all walks between nodes and x as a sum over walk lengths m, with each term weighted by combinatorial factor that counts the number of walks of said length. This factor is
(A m )iz, with A the unweighted adjacency matrix. That is,
Σ Wy E0
Σ m==
(A
m2 iz (-).
(17)

---

p- 0.2)
100, p 0.04)
Binary (N 127)
Ladder (N
GRFS
Frah norm eror
Trob. porm cror
Frob, norm error
Q-GRIs lourst
norm
10
Cx 10-1
4x10
10-
10
12
10
Karate (N 34)
Dol phins (V = 62)
Football (N 115)
Eurosis (N 1272)
Freb norm eror
Freb. norm encr
Freb. norm errer
1C-
20
k000m
Ex10-
10-1
No. random walks
No. random walks
No. random walks
No. random walks
Figure 2: Relative Frobenius norm error of estimator of the 2-regularised Laplacian kernel with GRFs (red circle) and q-GRF's (green cross). Lower is better. We also include q 9-RRW-GRIT [Choromanski 2023] (blue triangle), instantiated with f = exp, as benchmark. Our novel q-GRFs perform the best every graph considered. N is the number of nodes and, for the Erdōs-Réng (ER) graphs, p is the edge-generation probability. One standard deviation is shaded but is too small to easily see.
(A' AA") the is readily written as its eigendecomposition, Ep=1 Xo kipskpas with Ap the pth eigenvalue and kps the i-th coordinate of the p-th eigenvector kp. We put these into Eq. 16 and perform the sums over each walk length 2,,2,3,4 from to oc, arriving at
Σ E.ENN Ao2
Σ ΣE
A(A1, X2, Aa, da)kiskiz kan kis
(18)
where f(A1, A2, λ3 A4) is function of the eigenvalues of A that depends whether we correlate the terminations of the walkers. Since the eigenvectors are orthogonal we have that ΣzEN kizk2. = 612. so this reduces to
f(A1,A1, As. Sa
(19)
Ag.
Our task becomes prove that this expression becomes smaller when we induce antithetic termination. We achieve this by showing that particular matrix (with elements given by the difference in between the schemes) is negative definite: a task we achieve by appealing to Weyl's perturbation inequality [Bai al. 2000].
Experiments
In this section report empirical evaluations of G-GRss We confirm that they give lower kernel estimator variance than regular GRFs and show that this often leads to substantially better performance in downstream tasks, including simulation of graph diffusion, k-means node clustering and kernel regression for node attribute prediction. We use ensembles of antithetic pairs (Def. 3.1
4.1 Estimation of the 2-regularised Laplacian kernel
We begin with the simplest of tasks: estimation of the 2-regularised Laplacian kernel,
[I + o'
(20)
lap
where L E NNN, the symmetrically normalised Laplacian. < o < is a regulariser. We use both GRFs and q-GRFs to generate unbiased estimates ((() (see App. 8.2). then compute the relative Frobenius norm | 1
(2)
K || (2) ( between the true and approximated kernel
matrices. This enables to compare the quality of the estimators. As benchmark, also include
7

---

an implementation of the high-level reinforced random walk QMC mechanism suggested (but not tested) by CChrrmmansi 2023]. We choose the exponential mapping as the reinforcement function (used to downweight the probability of traversing previously-visited edges), although the optimal choice remains an open problem. We refer to this mechanism as g-RRW-GRF's to disambiguate from our instantiation of q-GRFs (which uses antithetic termination).
Fig. 2 presents the results for broad class of graphs: small Erdös-Rényi, larger Erdös-Rényí.
binary rooted tree, a ladder, and four real-world examples from Ivashkin 2023 (karate, dolphins,
football and eurosis). We consider 2, 4, 8 and 16 walks, taking 100 repeats for the variance of the approximation enror. We use the regulariser o = 0.1 and the termination probability = 0.5.
The quality of kernel approximation naturally improves with the number of walkers. Inducing anti- thetic termination consistently reduces estimator variance, with our q-GRF mechanism outperforming the particular graph (according
regular GRFs in every case. The exact size of the gain depends to the closed forms derived in App. 8.44, but improvement is always present. It is intriguing that tree-like and planar graphs tend enjoy bigger gap; we defer a rigorous theoretical analysis to future work. Meanwhile, the q-RRW-GRF variant with exponential f is often worse than the regular mechanism and substantially more expensive. We do not include in later experiments.
4.2
Scalable and accurate simulation of graph diffusion
In Sec. 2.1 we noted that the Laplacian L (or L) is the natural operator to describe diffusion on discrete spaces and that the 1-regularised Laplacian kernel constitutes the corresponding discrete time-evolution operator. Here we will show how, by leveraging q g-GFss provide a lower-variane, decomposition kpp we can simulate graph diffusion in a scalable and accurate way.
Choosing a finite (even) number of discretisation timesteps No. we can approximate the final state by
(2)
NIN
N₂
ub =
lap
Ug
(21)
This is equivalent to approximating the action of the diffusion kernel Kair (see App. 8.2p. We can efficiently compute this using our GRF q-GRF decomposition of K( and compare the accuracy of reconstruction of u,
| In particular, we take an initial one-hot state uq = (1, O..... 0) and simulate diffusion for 1 1 divided into N₂ = 1000 timesteps, using = 10 walkers and termination probability p = 0.5. Fig. 3 gives schematic. We average the MSE of i, over 1000 trials. Table reports the results; q q-GRFs approximate the time evolution operator more accurately so consistently give lower simulation error. Since the time-evolution operator is applied repeatedly, even modest improvements in its approximation can lead to substantially better reconstruction of the final state. In extreme cases the simulation error is halved.
Table I: MSE of the approximation of the final state u, by repeated application of an approxi- mation of the discrete time-evolution operator, constructed with GRFs Brackets give one standard deviation.
Figure 3: Schematic of diffusion on a small graph, with heat initially localised node spreading under the action of the Laplace operator. We use (q-)GRF's estimate the quan- tity on the right.
single
Graph
N
Sim error, MSE(Ü.) )
GRFs
q-GRFs
Small ER
20
0.0210(5)
0.0160(3)
Larger ER
100
0.0179(9)
0.0085(3)
Binary tree
127
0.0161(6)
0.0106(3)
Ladder
100
0.0190(8)
0.0105(3)
karate
34
0.066(2)
0.054(1)
dolphins
62
0.0165(4)
0.0139(3)
football
115
0.0170(3)
0.0160(2)
eurosis
1272
0.089(1)
0.084(1)
t t=0 0
"To be asymptotically unbiased, we should construct N independent estimates of k(2) (see App. 8.2 However, here we observe that the error is small when we apply the same approximation repeatedly. T.e.
estimate (2 once and use to simulate diffusion at every timestep, despite this remaining biased N.
8

---

4.3 Kernelised k-means clustering for graph nodes
Next we evaluate the performance of GRFs and q-GRFs on the task of assigning nodes to clus- ters using the graph kernel, following the ap- proach described by Dhillon et al. [2004. We first run the algorithm to obtain N: = 2 clusters with the exact 2-regularised Laplacian kernel, then compare the results when we use its ap- proximation via GRFs and q-GRFs. In each case, we report the clustering error, defined by
Table 2: Errors in kernliised k-means clustering when approximating the Gram matrix with q- GRFs and GRFs. Lower is better.
Graph
N
Clustering enrrr, E
GRFs
q-GRF's
karate
34
0.11
0.05
dolphins
62
0.40
0.38
polbooks
105
0.28
0.28
football
115
0.10
0.09
Ee N(N 1)/2
wrong pairs
(22)
eu-core
1005
0.04
0.05
databases
1046
0.17
0.11
eurosis
1285
0.15
0.19
This is simply the number of misclassified pairs (in the sense being assigned the same clus- ter when the converse is true or vice versa) di- vided by the total number of pairs. The results are less straightforward because curiously the clustering error does not generally vary monotonically with the variance on the kernel estimator, but nonetheless in six out of cight cases q-GRFs provide equally good or better results.
citeseer
3300
0.02
0.01
4.4 Kernel regression for node attribute prediction
Lastly, we consider the problem of kernel regres- sion on a triangular mesh graph. Each node is associated with normal vector ((),, equal to the mean of the normal vectors of its surround- ing faces. The task is to predict the directions of missing vectors (a random 5% split) from the remainder. We consider five meshes of different sizes, available from Dawson- -Dawsn--aggertty (2023) Our prediction for the direction the normal vector corresponding to the ith node is
Table 3: - cos & between true and predicted node vectors when approximating the Gram ma- trix with q-GRFs and GRFs. Lower is better. Brackets give standard deviation.
Graph
N
Pred error, 1 ons 0
GRFs
q-GRFs
cylinder
210
0.104(1)
0.101(1)
leapot
480
0,0531(5)
0.0493(5)
idler-riser
782
0.0881(6)
0.0852(5)
busted
1941
0.00690(4)
0.00661(4)
E, E, K (1, ()))
torus
4350
0,00131(1)
0.00120(1)
(23)
E (i,
where sums over the vertices with known vectors (a 95% split). This is simply a linear combination of all the training node normal vectors, weighted by their respective kernel evaluations. We compute the average angular eлror 1 - cos 0 between the prediction ((() and groundtruth (8) across the missing vectors, comparing the result when K is approximated with GRFs and g-GRFs with m random walks at termination probability p 0.5. The regulariser is 0.1. q-GRFs enjoy lower kernel estimator variance so consistently give better predictions for the missing vectors.
Conclusion
We have proposed novel class of quasi-Monte Carlo graph random features (q-GRFs) for unbiased and efficient estimation of kernels defined on the nodes of graph. We have proved that our new algorithm, which induces negative statistical correlations between the lengths of graph random walks via antithetic termination, enjoys better convergence properties than its regular predecessor (GRFs). This very often permits better performance empirical tasks and for some applications the improvement is substantial. Our work ushers in further research in this new domain of quasi-Monte Carlo methods for kernels on combinatorial objects. It may be of broader interest including for algorithms that sample random walks.

---

Broader impacts and limitations
We envisage several possible impacts of our novel algorithm. First, graphs provide natural way to describe systems characterised by complex biological interactions such as the proteins at the proteome scale or drugs in the body [Ingraham et al. 2019]. Our novel QMC algorithm might be of translational impact in this bioinformatics setting. Antithetic termination is at its heart procedure to improve the sampling efficiency
random walks, so could also help mitigate the notoriously high energy and carbon cost of large models (Strubell et al.) 2019]. Lastly, we believe our results are of intrinsic interest as the first (to our knowledge) rigorously studied QMC scheme defined on a combinatorial object. They might spur further rescarch in this new domain. Our work is foundational with immediate direct negative socictal impacts that we can see. However, it important to note that increases in scalability afforded by GRF and q-GRF algorithms could amplify risks of graph-based machine learning, either from bad actors or as unintended consequences.
The work also has some limitations and directions for future research. First, although we have derived closed-fore expressions for the kernel estimator variance with GRFs and q- q-GRFs in App. 8.4] they are still complicated functions of the spectra of the respective graph Laplacians. Understanding what characterises graphs that particularly benefit from antithetic termination (empirically, tree- like and planar graphs) is an important future direction. Moreover, our scheme only correlates walk lengths. more sophisticated mechanism that couples walk directions might do even better. Lastly, further work is needed to fully understand the applicability of antithetic termination beyond the GRF setting.
Relative contributions and acknowledgements
IR devised the antithetic termination scheme, proved all theoretical results and ran the experiments in
Secs [.1 4.2 and 4.4, KC provided crucial support throughout, particularly: running the clustering
experiment Sec. 4.3 showing how d-regularised Laplacian kernels can be used to approximate the graph diffusion process, and proposing to apply q-GRFs to the problems in Secs and 4.4 AW gave important guidance and feedback the manuscript.
IR acknowledges support from a Trinity College External Studentship. AW acknowledges support from Turing AI fellowship under grant EP/V025279/1 and the Leverhulme Trust via CFI.
We thank Austin Tripp and Kenza Tazi for their thoughtful feedback on earlier versions of the text, and Michael Ren for his excellent suggestion to treat the negative definite property perturbatively. We are also grateful to the anonymous reviewers.
10

---

References
Réka Albert and Albert-László Barabási. Statistical mechanics of complex networks. Reviews of modemn physics, 74(1):47, 2002. htt.httpp:///oi.rrrg////// 1tt///100.110///ee1 1e110/////oodPyysseee..e 74.47
Haim Avron, Michael Kapralov, Cameron Musco, Christopher Musco, Ameya Velingker, and Amir Zandieh. Random fourier features for kernel ridge regression: Approximation bounds and sta- tistical guarantees. In Proceedings of the 34th International Conference Machine Learning, ICML 2017, Sydney, NSW, Australia, 6-11 August 2017, volume 70 of Proceedings of Machine .ww.eng.............m........... press
Learning Research, pages 253-262. PMLR, 2017. URL
v70/avron17a html
1mmm
Zhaojun Bai, James Demmel, Jack Dongarra, Axel Ruhe, and Henk van der Vorst. Templates for the solution of algebraic eigenvalue problems: a practical guide. SIAM, 2000. URL https 7/ /www..cm cs. ucdavis. edu/" "////bai/EE1 Ttts/ET//ooteetts..tst contents. html
Mikhail Belkin and Partha Niyogi. Laplacian eigenmaps for dimensionality reduction and data representation. Neural computation, 15(6):1333-136, 2003. URL https: acm. org/doi/ 10. 00.1662/089976603378880113
Kurt D Bollacker, Steve Lawrence, and C Lee Giles. Citeseer: An autonomous web agent for automatic retrieval and identification of interesting publications. In Proceedings of the second international conference on Autonomous agents, pages 116-123, 1998. UtL.https://// acm
org/doi/10 100oo/110..00 1e1111111118811 280765.0770 280786
Karsten M Borgwardt, Cheng Soon Ong. Stefan Schönauer, SVN Vishwanathan, Alex Smola, and Hans-Peter Kriegel. Protein function prediction via graph kernels. Bioinformatics, 21(suppl I): i47-i56, 2005. URL https::/doi..rrr http://///ooooooooaaaact////////
1093/
Colin Campbell. Kernel methods: survey of current techniques. Neurocomputiing, 48(1-4):63 -3,4 2002.
Stéphane Canu and Alex Smola. Kernel methods and the exponential family. Neurocomputing, 69
(7-9):714-720, 2006. URL http:::////::ooorr////// 1016/ neucom. 2005 12. .00
Xiaowei Chen, Yongkun Li, Pinghui Wang, and John Lui. A general framework for estimating graphlet statistics via random walk. arXiv preprint arXiv: 2016. URL https
//doi.org/10.. 48550//555/////
1603.07502
Youngmin Cho and Lawrence K Saul. Analysis and extension of arc-cosine kernels for large margin classification. arXiv preprint arXiv:1112.3712, 2011. URL ortt::://doi..or/10...
1112.31112.372
Krzysztof Choromanski, Valerii Likhosherstov, David Dohan, Xingyou Song, Andreea Gane, Tamas Sarlos, Peter Hawkins, Jared Davis, Afroz Mohiuddin, Lukasz Kaiser, et al. Rethinking attention with performers. arXiv preprint arXiv:2009. /4794, 2020. htRpssttt:://att http::::///o....//1././/// 48550/
arXiv arXiv.2009.11 2009.14794
Krzysztof M Choromanski, Mark Rowland, and Adrian Weller. The unreasonable effectiveness of structured random orthogonal embeddings. Advances in neural information processing systems, 30, 2017. URL httpp:://o// org/ 48550/arXis 1703. 00864
Krzysztof Marcin Choromanski. Taming graph kemnels with random features. In International Conference on Machine Learning, pages 5964-5977 PMLR, 2023. URL https:///o// org/10 48550/ arXiv
Fan R. K. Chung and Shing-Tung Yau. Coverings, heat kernels and spanning trees. Electron. J.
Comb., 6, 1999, doi: 10.37236/1444, URL
Fan RK Chung. Spectral graph theory, volume 92. American Mathematical Soc., 1997.
Keenan Crane, Clarisse Weischedel, and Max Wardetzky. The heat method for distance computation. Communications of the ACM, 60(11):90-99) 2017. URI.https: saattps://att //dl. acm. org/doi/10. 1145/
3131280

---

Anirban Dasgupta, Ravi Kumar, and Tamás Sarlós. A sparse johnson: Lindenstrauss transform. In Proceedings of the forty-second ACM symposium on Theory of computing. pages 341-350, 2010.
URL https: //dois org/10. 1145/1806683 1806737
Michael Dawson-Haggeryy
Trimesh repository, 2023. URL https://gttuu..com/m/kkk...
trimesh
Inderjit S Dhillon, Yuqiang Guan, and Brian Kulis. Kernel k-means: spectral clustering and nor- malized cuts. In Proceedings of the tenth ACM SIGKDD international conference on Knowledge discovery and data mining, pages 551-556, 2004. (tL/httts:///// acm. www..rgg/o.//1.....//// 1145/
1014052. 1014118
Josef Dick, Frances Y Kuo, and Ian H Sloan. High-dimensional integration: the quasi-monte carlo way. Acta Numerica, 22:133-288, 2013.
Jinyuan Fang, Qiang Zhang, Zaiqiao Meng. and Shangsong Liang. Structure-aware random
fourier kernel for graphs. Advances in Neural Information Processing Systems, 34:17681- 17694, 2021. URL ptpp//ttt////rcceeningg... http://oooeeddinn.eeererrrpppppp cc/ http://ar....lleppppp/////////// pdf
Dániel Fogaras, Balázs Rácz, Károly Csalogány, and Tamás Sarlós. Towards scaling fully personalized pagerank: Algorithms, lower bounds, and experiments. Internet Mathematics, 2(3):333- 2005.
URL https:://dooo org/ 10. 1007
Michel X Goemans and David Williamson. Approximation algorithms for max-3-cut and other problems via complex semidefinite programming. In Proceedings of the thirty-thire anmual ACM symposium Theory of computing, pages 443 452, 2001. URL hittps://dot org. 10//10..11 1016/
j. 5.jcas cas. 2003.000000000000000 .1.007..1111 .012.
JM Hammersley and KW Morton. A new monte carlo technique: antithetic variates. In Mathematical proceedings of the Cambridge philosophical society, volume 52, pages 449 4794775. Cambridge University Press, 1956.
John Ingraham, Vikas Garg, Regina Barzilay, and Tommi Jaakkola. Generative models for graph- based protein design. Advances in neural information processing systems, 32, 2019. URL https:/ htt://wwwmttt wtt.:mitttttttmmemmmmeemmm edu/"vgarg/ atgaar//eeeerttttaaloodeeeeorrro
.dff
Vladimir Ivashkin. Community graphs repository, 2023.
URL. https://githi: com/
wlivashkin/ /Livashhin/ccomuntt--rrrnn http:/ootoonnttgaaasss..
William B Johnson. Extensions of lipschitz mappings into hilbert space. Contemp. Math., 26:
189-206, 1984.
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 2014. URL huttps://oii. org/10. wtw./..../..///////////////// aar//v.V... 1403. 3148
Risi Kondor and John D. Lafferty. Diffusion kernels on graphs and other discrete input spaces. In Machine Learning, Proceedings of the Nineteenth International Conference (ICML 2002), University of New South Wales, Sydney, Australia, July 8- 12. 2002, pages 315-322 2002. URL cmu wwweeeeeeerceeaaaeeeeeeeeee httpeesaarrccaaappppsrsssss
attps://wwy,
diffusionnn--o kerneleeeelll.... pdf
Leonid Kontorovich, Corinna Cortes, and Mchryar Mohri. Kernel methods for learning languages. 06.033..66... URL //doi
Theor Comput. Sci., 203(3):223-22) 236, 2008. doi: 10.
org/10 10 2og/10.1016/j.1 tcs. 2008.00000000 037
Jure Leskovec, Jon Kleinberg, and Christos Faloutsos. Graph evolution: Densification and shrinking diameters. ACM transactions on Knowledge Discovery from Data (TKDD), 1(1):2ees, 2007. URL
https: acm. 1tt.//////.111111111199999999990 121131
Fanghui Liu, Xiaolin Huang, Yudong Chen, and Johan A. K. Suykens. Random features for kernel approximation: A survey on algorithms, theory, and beyond. IEEE Trans. Pattern Anal. Mach. Intell, 44(10)) 2022. doi: 009//PAAM..022.330000000 URI. https
//doi fooo.orrg/100 10. 1109/ /TPMMI 2021 3097011
12

---

David Lusseau, Karsten Schneider, Oliver Boisseau, Patti Haase, Elisabeth Slooten, and Steve M Dawson. The bottlenose dolphin community of doubtful sound features a large proportion of long- lasting associations: can geographic isolation explain this unique trait? Behavioral Ecology and URL htlphhttp://iiin http::::///innsssiinneeee springer... cow httm//aticle/1111000
Sociobiology, 54:396-405
2003.
500265-003-0651-1
Yueming Lyu. Spherical structured feature maps for kernel approximation. In International Confer- ence on Machine Learning, pages 2256-2264. PMLR, 2017. URL http://procceedings mlr
press//70//Ayuu7a///ppll1
Mark EJ Newman. Modularity and community structure in networks. Proceedings of the na-
fional academy of sciences, 103(23):5577-8882, 2006. htt.httts://do..rrg// wt.//////...00 pna
0601602103
Giannis Nikolentzos and Michalis Vazirgiannis, Random walk graph neural networks. Advances in Neural Information Processing Systems, 33:16211-16222, 2020. URL htppt:tt///pooeeedinng Paperr pdf
Lawrence Page. Sergey Brin, Rajeev Motwani, and Terry Winograd. The pagerank citation ranking:
Bring order to the web. Technical report, Technical report, stanford University, 1998. URL
httpp://www. cis. upenn. http://ddeerrsseecccaaal/lleeeee www.rr/eee....eeee.aaaaaaaaa pdd
Ali Rahimi and Benjamin Recht. Random features for large-scale kernel machines. Advances in neural information processing systems, 20, 2007. httpps:p:///p eecs. berkeleyy
edu/"brechte wwwpechttpppers///eee 07 rrah.rec..r nic....s..ddf
Ali Rahimi and Benjamin Recht. Weighted sums of random kitchen sinks: Replacing min- imization with randomization in learning. In Advances in Neural Information Process- ing Systems 21, Proceedings of the Twenty-Second Annual Conference on Neural Infor- mation Processing Systems, Vancouver, British Columbia, Canada, December 8-11, 2008, 2008. https:tttt://paaaata www..iipppppaaaeeee paappr_filesssss pdf
pages 00re372422301133054444a000000
Isaac Reid, Krzysztof Marcin Choromanski, Valerü Likhosherstov, and Adrian Weller. Simplex random features. In International Conference on Machine Learning, pages 28864-28888. PMLR, 2023. URL httpt: otg.dooo.ogg////oo 10. 48550//5Xiii 2301 13856
Pedro Ribeiro, Pedro Paredes, Miguel EP Silva, David Aparicio, and Fernando Silva. survey on subgraph counting:
concepts, algorithms, and applications to network motifs and graphlets. ACM Computing Surveys (CSUR), 54(2):1-36, 2021. http:ttpps:/.ii.rrr///// 1145/34336523
Sam T Roweis and Lawrence K Saul. Nonlinear dimensionality reduction by locally linear embed- ding. science, 29065500) 2000. URL http::://ww..tt http://wwwwwsiinncccco ort. 1126/ science 5500.. 2323
Nino Shervashidze and Karsten M. Borgwardt. Fast subtree kemness graphs. In Advances in Neural Information Processing Systems 22: 23rd Annual Conference on Neural Information Processing 2009, Vancouver, British Columbia, Canada, pages 1660-1668, 1668. Curran Associates, Inc., 2009. URL neurips.
Systems 2009. Proceedings of meeting held 7-10 December
waw.ceaaaa00ee--bbaaaattet Abstract httml
Nino Shervashidze, S. V. N. Vishwanathan, Tobias Petri, Kurt Mehlhorn, and Karsten M. Borgwardt. Efficient graphlet kernels for large graph comparison. In Proceedings of the Twelfth International Conference on Artificial Intelligence and Statistics, AISTATS 2009, Clearwater Beach, Florida, USA, April 18-18, 2009, volume 5 of JMLR Proceedings, pages 488-495, 2009. URL http: Utthttp! \/proceedings. mlr. wwwpppppsseeepppppppp https////bbbSrrraaaade00.0. html
Alexander J. Smola and Risi Kondor. Kernels and regularization on graphs. In Computational Learning Theory and Kernel Machines, 16th Annual Conference on Computational Learning
Theory and 7th Kernel Workshop, COLT/Kernel 2003, Washington, DC, USA, August 24-27, 2003, Proceedings, volume 2777 of Lecture Notes in Computer Science, pages 144-158. Springer, 2003. URL https: //doii org/ 1007 0977/997-000 540 544-540-----9-1222
13

---

Alexander J. Smola and Bernhard Schölkopf. Bayesian kernel methods. In Advanced Lectures on Machine Learning. Machine Learning Summer School 2002, Canberra, Australia, February 11-22, 2002, Revised Lectures, volume 2600 of Lecture Notes in Computer Science, pages 65-117. 2002. URI. htt.:/tttss///aa...o////0001/1/1
Emma Strubell, Ananya Ganesh, and Andrew McCallum. Energy and policy considerations for deep learning in nlp. arXiv preprint arXiv:1906.022244 2019. URL https:///oi. org/ wttg//..// 18653/
v1/P19-1355
Ashish Vaswani, Noam Shazeer, Niki Parmar, Jakob Uszkoreit, Llion Jones, Aidan N. Gomez, Lukasz Kaiser, and Illia Polosukhin. Attention is all you need. In Advances in Neural Information Processing Systems 30: Annual Conference on Neural Information Processing Systems 2017, De- cember 4-9, 2017, Long Beach, CA, USA, pages 5998-6008, 2017. / www.aaaaaaaa--aaattttt... .htll
neurips. cc/ppprr/////////////////
S. V. N. Vishwanathan, Karsten M. Borgwardt, and Nicol N. Schraudolph. Fast computation of graph kernels. In Advances Neural Information Processing Systems 19, Pruceedings of the Twentieth Annual Conference on Neural Information Processing Systems, Vancouver, British Columbia, Canada, December 4-7, 2006, pages 1449-1456. MIT Press, 2006. www.hhtsss:/ppcceedeeeee htwp://///bbB---bbbbaaaeea Abstract html
neurips. cc/ wwtppcppppcccppppccccccccc.
Christopher KI Williams and Carl Edward Rasmussen. Gaussian processes for machine learning, volume 2. MIT press Cambridge, MA, 2006.
Yasutoshi Yajima. One-class support vector machines for recommendation tasks. In Pacific-Asia Conference on Knowledge Discovery and Data Mining, pages 230-239, Springer, 2006. URL
httpps//ddl. acm. wwwo/og///.....o /tt/:0.i111111/ 1007/1173113 org//bb0/10..100//133333 28
Jiyan Yang. Vikas Sindhwani, Haim Avron, and Michael W. Mahoney. Quasi-monte carlo feature maps for shift- invariatt kernels. In Proceedings of the 31th International Conference on Machine Learning, ICML 2014, Beijing, China, 21-26 June 2014, volume 32 of JMLR Workshop and Conference Proceedings, 2014. http http:///poeeeeen...eeee mlr.
html
Felix Xinnan X Yu, Ananda Theertha Suresh, Krzysztof M Choromanski, Daniel N Holtmann-Rice and Sanjiv Kumar. Orthogonal random features. Advances in neural information processing systems, 29, 2016. URL htttstss: http:///oo.oog//..//// 48550/ arXiv. 1610. 09072
Wayne W Zachary. An information flow model for conflict and fission in small groups Journal of anthropological research, 33(4):452-47,, 1977. URL https://t:/// www.www....rrall..o.m uchicago
33. 3629752
Yufan Zhou, Changyou Chen, and Jinhui Xu. Learning manifold implicitly via explicit heat-kernel learning. Advances in Neural Information Processing Systems, 33:477-487, 2020. URL https
//doi. (eddi.or//////
org/ 10tg//1.... 48550/ arxiv. 2010.0100 01761
14

---

8 Appendices
Notation:
reduce the burden of summation indices, in places we use Dirac's bra-ket notation from quantum mechanics. (1) can be interpreted as the vector and (y) as y
8.1 Antithetic termination for PageRank
Here demonstrate the broader applicability of antithetic termination beyond the GRF setting by showing how it can be used to improve the convergence of numerical estimators of the PageRank vector. The treatment is similar to App. 8.4 the reader may benefit from reviewing this section first or in parallel.
The PageRank vector is a measure of the relative importance of the nodes in graph, proposed by Page et al.] [1998] to rank websites in search engine results. It is the steady state of Markov chain with N X N transition matrix
P E.
P= (1 p)T
(24)
where p E (0, 1) is scalar and N = LV is the number of graph nodes. The matrix E == []]ee has entries that are all equal to 1, and the matrix given by
1 Túj
E otherwise
(25)
O
with d, = Ejex [((i,)) E E), the degree (number of neighbours) of the ith node. The state space of the Markov chain is the nodes N. T encodes the behaviour of a simple random walk.
Formally, we are solving for the vector π E IR
such that
TP=- 11=1
(26)
where 1 == 1,.,,..., 1). Physically, πy is the stationary probability that surfer is located at node j. supposing that at every timestep it either i) moves to of its neighhouring nodes (with probability
1 - p) or ii) teleports to random node (with probability p).
For large graphs, becomes very expensive to compute exactly, motivating a host of research dedicated to approximating it. Solving Eq. and Taylor expanding (1 - (1 - -p)T)- -1 is given by
the solution
8
P N JEN K=0
p) TT
(27)
Note that the term
P >
(1-p)
(28)
is nothing other than a sum over all walks between nodes j and i, weighted by p(1 f with & the walk length and d, the degree of the ith node visited on the walk. Supposing a walker terminates with probability p at each timestep and picks neighbour to visit at random, this weight equal to the probability of sampling a particular walk between nodes j and i, allowing to write
II
Σ JEN
p (walker starting at terminates i).
(29)
This invites the Monte Carlo estimate
[(walker starting at j terminates at
(30)
with m the number of walks sampled out of each node. We have arrived the algorithm proposed by Fogaras et al. [2005]: simulate m walks out of each node j of the graph, and compute the estimator π, as the proportion of walkers terminating node i.
Since we are sampling random walks that terminate with probability p at every timestep, this is natural setting to try antithetic termination. The following is true.
15

---

Theorem 8.1 (Superiority of antithetic termination for PageRank estimation). Supposing an ensemble of random walkers exhibits either i..dd or ii) antithetic termination and p S 2
antiiiiteei T::) S
Variad (T:)-
(31)
Proof: It is sufficient to consider the behaviour of single pair of walkers out of node j that are either i..d. or exhibit antithetic termination. Since in both cases the marginal probability of walkers terminating at node is identical, the estimators are unbiased. It is sufficient to determine whether the expectation value
TE (T[22 terminates at node ] [[22 terminates at node [} = p(both terminate at i) is suppressed by antithetic termination, where l1 denotes two walks simulated out of j.
(32)
If the walkers independent, Eq. 32 evaluates to
N
2
pboth terminate at t) = p² Σ
p) (ilbaa (Na
(33)
a=1
where Xa and (ka) denote the a-th eigenvalue and eigenvector of the transition matrix T respectively, Meanwhile, for walkers with antithetic termination, it is straightforward to derive that
1-3p
lf < 4,
I-P
I-P
((2) l2 (e1)) (1)
if l2 4
(34)
(1-30)
(1 - 2
if l2 > 4
where the walk lengths are no longer independent. It follows that
11
1 2p
p TT
p(l2 terminates at i len(S),) = (1) =
p
1==0
(35)
2p
)
4
(1
111))~--111p
p
4=1++1
We then sum over all the possible walks Q2 at each length 41 After tedious but straightforward algebra,
p²
(1- p)daa
p(both terminate at i)
Σ
(ilka) (kial)) (ilkg) (kgli)
T- (1 - -p)>a
p 0,8m1
(36)
(- - (1 -1110 p) To assess whether antithetic termination suppresses the variance, we need to consider the difference between Eqs 33 and 36 This evaluates to
1
1- (1 Au
1
(Aa 1) 1)
Ap(both terminate at i) =
Σ a,3=1
-p² 1- - p)>a - (1 p)d 1 - (1 Ap)
(37)
(ilk:a) (baalt) (j) (i)ka)) (kell))
For any pair of nodes i, the products (jlka) (ka () and (ilkg) (k can be treated the and (9- coordinates of an N-dimensional vector. Then we simply require that the matrix whose o, B-th element is given by
1
(2. - 1)(A 1)
(38)
1- 11-1-))A, p) - (1 2p)
negative definite, which is manifestly the case provided S (e.g. by Taylor expanding the denominator and observing that it is sum of terms symmetric in Aa. with negative coefficients). It follows that Ap(both terminate at i) S O, and the claim in Eq. 31 follows.
We run simple numerical test to check Thm 8.1 for number of graphs (see Sec. 4.1 of the main text for details), and find that using antithetic termination does indeed improve the quality of estimation of the PageRank vector. Fig. 4 reports the results, showing the error on the estimator || क ll2 for both schemes with 2, 4, and walkers. We choose a termination probability p = 0.3. As per the theoretical guarantees, antithetic random walks consistently yield more accurate estimators #.
16

---

20, p-
100, p 0.04)
Binary
IV
lid.
Pager eror
Pagger ant ArTOr
Paper ank ast, eror
Pagerant
10
13
12
52
random
random
Karatt (W = 34)
Dolphins (N 62)
Football (N 115)
Eurosis (N 1272)
Pugarank est. crrer
Pegerank
Po Pegerank
No. random walks
No. random walks
No. random walks
No. random walks
Figure 4: Error the estimator the PageRank vector | + when approximated using i) i..dd and ii) antithetic random walks. Lower is better. Antithetic termination yields more accurate estimate for every graph considered. N is the number of nodes and, for the Brdōs-Rényí (ER) graphs, p is the edge-generation probability. One standard deviation is shaded but it is too small easily see.
8.2 On the approximation of the d-regularised Laplacian using GRFs
In this appendix, we demonstrate how to approximate the d-regularised Laplacian k with GRFs.
Recall that GRFs provide an estimator to the quantity (IN U) where U is a weighted adjacency matrix. Recall also that the matrix elements of the symmetrically normalised Laplacian L given by
i J if i
Lij
Wi
(39)
deww de
where degw (t) EJEN JEN Wiy is the weighted degree of the node i. We are often interested in situations where W A, an unweighted adjacency matrix. Now note that
K (2) lup = NIN + o oDL); (1 0 (IN U),²²
(40)
where we defined the matrix U with matrix elements
W
Ug 1 ddegw (i)
(41)
This is itself a weighted adjacency matrix. It follows that, by estimating (IN - U) with GRFs, we
estimate K(2) This reported by Choromanski 2023
Supposing that we have constructed a GRF estimator
k(() lap
= |CC
(42)
where the matrix C E N has construct the 1-regularised Laplacian kernel estimator
C; 1 (1)
1
note that it is straightforward to
K
(1) lap
= E CD
(43)
by taking D (IN o²L,)
C. By multiplying independent copies of the estimators
K and
(1) it is then trivial to obtain unbiased estimator K(O for arbitrary d E N. For instance,
K (5) C,
(44)
TTechiically, for unbiasedness the random walks used for C and CT should be different, though this rarely makes any observable difference. See App. 8.7 further discussion.
17

---

Note also that
Kaif exp(-oOLt) lim
L' 1-
1
n
1
o Lt
(45)
n
with >> 1 the number of discretisation timesteps. We again included a regulariser a. This shows that the ability to approximate the d-regularised Laplacian kernel for arbitrary d permits us approximate the diffusion kernel, with the corresponding estimator asymptotically unbiased as d -> oc. In Sec. 4.2 see that the quality of approximation (in this case measured by the accuracy of simulated diffusion) is good even with modest d.
We also note that the variance of estimates [ and Kain constructed in this manner will vary mono- tonically with the variance of the approximation of the constituent estimators k so improvements provided by our QMC scheme (and any theoretical guarantees) will follow for this more general class of kernels.
8.3 Walk lengths with antithetic termination are anticorrelated (derivation of Eq. 15
In this appendix derive Eq. [5| which gives the expected length of some walk w2 given that its antithetic partner w/ is of length m: that is, IE (Ien(w>) (en((w)) = m))
As warm-up, consider the simpler marginal expected lengths. Note that
p (len(w) m) (1 pp
(46)
follows that
= Σ m(1 p) 1 1-P
E (len(w))
(47)
P
m=0
where we evaluated the cechheetic-gemmettii series. We reported this result in Eq. 14 Meanwhile, the probability walk being of length given that its antithetic partner is length m is
1-2p T-P
if < m,
T-p
p(len(w2) illen(wi) mm)
if i = m,
(48)
1p
(1 - p)
4-m-1
if i >> m.
T-p
The analagous sum then becomes
m-1 1 2p
1 1- 2p m
E (len(wy) [(e(())) m)
p
Σ
(1 p) i-m pi. (49)
1
1-p
1 - p
After straightforward but tedious algebra, this evaluates to
1
2p
2p m
E (en(w)))len((()
m)
+
(50)
p
1
Note that this is greater than [[[en(w)) when m is small and smaller than [[le((w)) when is large; the two walk lengths negatively correlated.
as stated in Eq. 15
8.4 On the superiority of q-GRFs (proof of Theorem 3.2
Here, provide a proof of the central result of Theorem 3.2 termination reduces the variance of estimators of the matrix [IN U)-2. From App. 8.2 all our
that the introduction of antithetic
results will trivially extend to the 2-regularised Laplacian kernel 1ap
We will begin by assuming that the graph is d-regular, that all edges have equal weights denoted w, and that our sampling strategy involves the random walker choosing one its neighbours with equal probability each timestep. We will relax these assumptions in App. 8.5
18

---

We have seen that antithetic termination does not modify the walkers' marginal termination be- haviour, so the variance of the estimator ((i) T (x() is only affected via the second-order term [()(1] Φ()))].. Writing out the sums,
(0(1) (
Φ())
Σ Σ
Σ E
Σ
x(wi) ((w)) ((سا)) x(wا))
m
wi Eft, w2 wy 44 p(wi) plws) p(wa) p(ww)
2.yEN 11.11. A2-
I[w E (2 ()) e [(( likk, i)) (2((2. j)). (51)
To remind the reader: the variables z, sum over the nodes of the graph M. 1ı and 4 enumerate all the m walks sampled out of node whilst ky and la enumerate walks from j. The sum over wi E Oiz is over all possible walks between nodes and x. (ww) evaluates the product of edge weights traversed by the walk w1, which is (e((()) in the equal-weights case (with len(wr) denoting the number of edges in w1). p(w1) is the marginal probability of the prefix subwalk w1, which is equal ((1 - p)/ dd)len(ar) on d-regular graph. Lastly, the indicator function I(w1 E ((k), evaluates to if the ki th walk out of node i (denoted (2 i)) contains the walk w1 prefix subwalk and O otherwise.
We immediately note that our scheme only ever correlates walks leaving the same node, so walks out of different nodes remain independent. Therefore.
E [[wi] E ((kk, E ((h:s)) E O1k,, i)) (2(2,1))
(52)
= =p(wi E Oiki, ws 02k2, i))p(us 5) ()2,,
Consider the the sum corresponding to one particular walks (k1, 1 k2 12),
Σ3ΣΣN Σ
((wا) xحwا) [سwاة) (مسهت)
p(wn) p(wz) p(ws) p(w))
a.VEN En wzED wgE0 4e En
(53)
-p(wi (2(k),)) ws E f(kk, 1))s E (((,.j) E ((2,j)))-
This object will be of central importance and referred to as correlation term. the sum over k1, k2, l1, 12, there are three possibilities to consider. We stress again that 1 refers to pair of walks out of node and 2 refers to a pair out of j.
Case 1, same-same, k; k2, 4 lz: the pair of walks out of are identical and the pair of walks out of j are identical. This term will not be modified by antithetic coupling since the marginal walk behaviour is unmodified and walks out of different nodes remain independent.
Case 2, different-different, k * 12, 4 * l2: the walks out of both i and differ, and each pair may be antithetic or independent. This term will be modified by the coupling.
Case 3. same-different. &1 K2. 4 * lz: the walks out of differ and may exhibit antithetic or independent termination but the walks out of j are the same. This term will be modified by the coupling. Note that the and j labels are arbitrary so we have chosen one ordering for concreteness.
If we can reason that the contributions from each of these possibilities 1 3 either remains the same or is reduced by the introduction of antithetic coupling, then from Eq. 5T we can conclude that the entire sum and therefore the Laplacian kernel estimator variance is suppressed. For completeness,
19

---

write out the entire sum from Eq. 51 with the degeneracy factors below:
(0(1) (0(1) -0)))² = m
Σ
Z.NEN
n²² Ewiern eRie EwzeR EuaeAry Ewiesty plun) plwa) p(wa) p(=a))
(ifnn) ((tar) (lses) (iwes))
same-same (1)
II(n (1k,, ()) E R(h,, f)) E (1ki, ()) E ())
+m²(m 10-11] 1) En)) EngERfE Enger, Emeeg p(wi) p(wz) p(ws) p(wa)
afen) ((wa)) l(wa) ((w))
different-different (2)
1(wi (2(k ()) MN) E OKk ()) E (l2,,
D(ws) D(ws) ((wn) plan) -m²(m 1) Ewienti Emseajo Emserit Emeetty plaa) plwa) plwa) p(ma)
T(w1 E (2(k1, ()) E ((11))) E l1kk, ()) E (11,)))
B(w1) ((wa) ((wn) ((m)) +m²(m Enjen Ewyea Engent Emeen plwi) plwz) plws) plua)
same-different (3)
T(w1 f2(k,, (ww ((11,)))T(( (2(1,, T)) ((2.)))).)
(54)
We now address each case 1 - 3 in turn.
8.4.1 Case 1: k1
k2, 1 = l2
As argued above, case 1 is trivial. By design, antithetic termination does not affect the marginal walk behaviour (a sufficient condition for the estimator to remain unbiased). This means that it cannot affect terms that consider a single walk out of node and a single walk out of j. and all terms of case 1 are unchanged by the introduction of antithetic termination.
8.4.2 Case 2:
k2, l1 + l2
Now we consider terms where both the walks out of node and the walks out of node j differ. To emphasise, we are considering 4 different random walks: 2 out of i and 2 out of j.
Within this setting, will need consider the situations where either i) pairs exhibit antithetic termination rather than i.i.d.. Terms of both kind will appear when we use ensembles of antithetic pairs We need check that in both cases the result is smaller compared to when both pairs are i.i.d..
ii) both of the
To evaluate these terms, we first need to understand how inducing antithetic termination modifies the joint distribution p(w1 E (2kk, i), wa E ((kk,, i)): namely, the probability that two randomly sampled walks Q(ki, i) and U(k2, i) contain the respective subwalks W1 and wa. given that their termination is either iid, or antithetic. In the i.d.d case, it is straightforward to convince oneself that
p(wi (2(1,)),us i(3,)) (12²) m
(55)
d
where m and n denote the lengths subwalks w1 and w3. respectively. With antithetic termination, from Eq. 13 follows that the probability of sampling walk Zg of length conditioned sampling an antithetic partner 2g of length is
T-P
I-p
if <i
v(len(Aa)) []ee((Q)==))
if =
(56)
(122) (1 1ーー))ーーーーp
if >
Using these probabilities, it is then straightforward but algebraically tedious derive the joint probabilities over subwalks
1-2p
(1 - p) if m,
T-p
p(wr ((1,1), ()3
2m
2p) m
m, if n > m,
(57)
2p
(1 (1-2) (1-p)"
Fm
where is the length of wi, is the length of wa and i is now the index of a particular node.
20

---

To be explicit, we have integrated over the conditional probabilities of walks of particular lengths
(i. j) to obtain the joint probabilities of sampled walks containing subwalks of particular lengths (m, n)). Let us consider the case of < m as example. Using Eq. 56
plwi E (1(11,). E (2(3. 0) d
1 2p p
P_ (1
p) p+
I-p p
(58)
2p
(1 p) p p)²
j=++1
where the branching factors of d appeared because every timestep the subwalks have d possible edges to choose from. After we have completed the particular subwalks of lengths m and n we no longer care about where the walks go, just their lengths,
stop accumulating these multiplicative factors. Computing the summations in Eq. 588 (which are all straightforward geometric series), we quickly arrive at the top line of Eq. 57
Returning to our main discussion, note that in the d-regular, equal- case,
((w)) ((wa) p(ur) p(wr) p(wr (l(k,, w) wa (2(,, i))
Σ
Σ
W1 El
wig Eugy
(59)
Σ
Σ
wd
p(wi E (.(k, ws E ((k2,
I-p)
The summand depends only on walk lengths m. n but not direction, which invites us to decompose E (-) over walks between nodes i and x sum over walk lengths, weighted by the number of walks at each length. Explicitly.
the sum Ew, EQ..
- ())
(60)
n=0
with A the (unweighted) adjacency matrix. We have used the fact that (A (A"), counts the number of walks of length between nodes i and x. A is symmetric so has convenient decomposition into orthogonal eigenvectors and real eigenvalues:
N
(A Ji
Σ X (ilk) (k/²)
(61)
Am
where [k) enumerates the N eigenvectors of A with corresponding eigenvalues λs and and (ar are unit vectors in the i and æ coordinate axes, respectively. We remind the reader that we have adopted Dirac's bra-ket notation; (y) denotes the vector y and {y/ denotes T.
Inserting Eqs 61 and 60 into Eq. 599 and using the probability distributions in Eq. 55 and 57 our all-important variance- variance-determining correlation term from Eq 53 evvlluates to
N
Σ
Σ
B
kg Blgia (ilks) (kalz) (jlk)) (k2)z) (il ka) (kalv) (ikk) (kulg).
(62)
Z.BEN lo kg.
where the matrix elements Byg k2 and B(ú) corresponding to the pairs of walkers out of i and respectively, are equal to of the two following expressions:
Cki,As = 1- whe, 1-whes
ifi if antithetic.
B ky [ = whr,, 1- I-ow Ant
1- = Xe, Auea
(63)
1-2p Here, is constant defined by [I-p]² with the termination probability. These forms straightforward to compute with good algebraic bookkeeping; we omit details for economy of space.
Eq. 62 can be simplified. Observe that EPEN (ax) (ar IN ('resolution of the identity"), and that since the eigenvectors of are orthogonal (ki [k2) = Ō Applying this, we can write
N
B
B ks (ilki) (1k)) (álka) (ika))
(64)
=1
21

---

Our task is then to determine whether 64 is reduced by conditioning that either one or both of the pairs of walkers are antithetic rather than independent. That
(K Dki ,ky Cky (f(k1) (dlba) (i)bs) (j) kis)) S 0,
(65)
===1
E((sss Dky,ks Cky,ka Rgg Cki, (ilk)) (ilk)) (ilka) (ilks) & 0.
(66)
k1=1ka=1
Define vector y E TRN with entries Up == (ilkp) (j) (jlkp)), such that its pth element is the product of the and jth coordinates of the pth eigenvector kpp In this notation, Eqs 65 and 66 can be written
N
N
(Cpg - Cm) Dplla S 0.
(67)
Di - Cpe) Nylla 2
(68)
For Eqs 67 and 68 to be true for arbitrary graphs. it is sufficient that the matrices E and F with matrix elements Eng CpgDp - Cp and Fm == Dm Cp are negative definite. Our next task is to prove that this is the case.
First, consider E, where just one of the two pairs of walkers is antithetic. Putting in the explicit forms of Cpg and Dyg from Eq. 63
Xp) Prg (1- (1->))->>>) 2n))
2 p²
Xp 1-2p
(69)
(1 p)2
Ap
where for notational compactness we took Xp = wdp (the eigenvalues of the weighted adjacency matrix U). Taylor expanding,
Epg
[1 Ap)
2
Ap1n
p
OC:
(1 - p)2 Apr,y)
1 2p
(70)
Aa)
(1
p)²
Inserting this into Eq.
67
we get
N
N
p2
N
X
VI- 2p
2
E
yp
S
0,
(71)
(1 p)2
Xp
1
p=1q=1
m=0
which implies that E is indeed negative definite. Note that we have not made any additional assumptions about the values of p and beyond those already stipulated: namely, < 2 and
Xmax 1.
Next, consider F, where both pairs of walkers are antithetic. Again inserting Eqs 63
we find that
(1= (
2 -
Fpz =
(72)
where we remind the reader that
1-22. The Taylor expansion in Xpd, is
Fpg = (11 x,)(1 (1 Any) (1 (x--(11X111>)))))
)
2
1)c² 2ic 1 + + 1)c
(73)
not generically negative definite, but it will be sufficiently small p or w. Taylor expanding in w,
Fpg = 2 >> + Oww))
(74)
The (ww2) contribution manifestly negative definite because e < We cannot guarantee this for the higher order terms, but we can take the limit 0 and treat them perturbation. It is
22

---

clear that the spectral radius of the perturbation will go to 0 faster than the spectral radius of the leading term (all its contributions are proportional higher powers of w). follows from Weyl's perturbation inequality [Bai et al. 2000] that Fp will indeed by negative definite provided w is small enough. Therefore, at sufficiently small w correlation terms with both pairs antithetic are suppressed as required.
Taylor expanding - (which corresponds p - 0) instead of ApAq, we can make exactly analogous arguments to find that F is also guaranteed to be negative definite with when p is sufficiently
small. Briefly: let
with å (I-x)
Then have that
2 >
Pm ( (1- Ap) (1--))(1XX)) -xp) m)))
1 (11
1- Ap (I--
-1)
Xp22 + 6xpp1
(75)
( (1- -1ーー>>a>>)
1
2 25 25 + ((82)
Xpd2
Taylor expanding 1- it is easy that the operator defined by the (6) term of Eq. 75 negative definite, This part will dominate over higher order terms (which are not general negative definite) when 5 is sufficiently small, guaranteeing the effectiveness of our mechanism these terms.
This concludes our study of variance contributions in Eq. 53 where 1ą * 12, * 12- We have found that these correlation terms are indeed suppressed by antithetic termination when p(U) is small enough.
8.4.3 Case 3: k: = k2 11 o 12
We now consider terms where 1ą = l2 and 1 * l2. We are considering total of walks: just 1 out of node but pair (which may be antithetic or i.i.d.) out of node j. We inspect the term
wd
ΣΣ WIER WaER.
p(wi E 1(k,, 1), wa E (2(k),
(76)
1-p
where m denotes the length of wi and denotes the length of wa- What is the form of p(wi ((k,, i), wa E ()){ It is the probability that a single walk out node i, f2(ka, i), contains walks w1 between nodes and x and w3 between and y as subwalks. Such a walk must pass through all three nodes i, and y. After some thought,
if Ww W3, if wa E Wl if if E Wa:
p(u1,w3 E ((k),
(77)
otherwise.
Here, W1 E wa means wi is a strict subwalk of wa. so the sequence of nodes traversed i + x > Likewise, Wg E wi implies walk y x. Summing these contributions,
Σ Σ
wd
m+n
plun U1h, ws E f2(4.
I-p)
My El wny Ey
wd 2len(w1)
p(ws E 0(ki, 1)
wd 2len(w1)
ud ka(ws)
(78)
Σ
p(wi E (2(k),0)
Σ
plws (2(k), x))
I p
wi Ei
W WyER
wa i
+
wdd
2kenwwg)
plws E (2(k),
Σ
wd kn(wrg) plws E (())..9)..
I
ageRt
WAEft
wi WyG
23

---

We introduced Wi for the sum over walks between nodes z and y, and p(ws E (2(k1, a)) is the probability of some particular subwalk a + equal to (--2) len( in the d-regular case. Wi is dummy variable so can be relabelled w1. The -ettrmin correlation term from Eq. 53 becomes
Zlen(w1)
Σ
p(ur E Q(k,, 1))aay
1-p/
wd
(w1)(
plw1 (1(1)
Σ
wd
kun(ws)
plws E (1(11,
wn ERii
Hg ubeRta
(79)
wd
2len(w>)
p(wr ()k),
Σ
wd
len(ws)
plws E O2(11, v))
1-p/
1-p
ERg
Ma
B K2.44 (ilk2) (122 (ilks) (kaly).
where B depends on whether the coupling of the pair of walkers out of node , is i.i.d. or antithetic, as defined Eq. [3] x and y are dummy variables so can also be swapped, and the sum over the walks Wy computed via the usual sum over walk lengths and eigendecomposition of A. Using the resolution of the identity and working through the algebra, we obtain the correlation term
2le((w1)
wd
Σ
p(ui (2(), ())
-
sEN
Lui Eflir
(80)
N
1 - 2 Akg Aee
B() Kzika (2))ke) (Kel)) (z/(k) (kalj).
(1- wAk )1 wdkx)
k2.
Now observe that the prefactor in square brackets is positive for any node æ since it the expectation of squared quantity. This means that, for the sum in Eq. 80 to be suppressed by antithetic coupling, is sufficient for the summation in its lower line to be reduced. Defining a vector y E R with elements Up kx/)kp) (kpl)),, ), it becomes clear that we require that the operator J with matrix elements
w² Xp 1 wAp)(1 Jpg = (1L (1 why) ) (Dpg Cp4)
(81)
is negative definite. Using the forms in Eq.
63
2 Xp
(1 (1-p)²
p²
w² Xp
Jpq
(82)
(1 - ) wxy))
2p
2 Xp
Making very similar arguments to in Sec. 8.4.2 (namely, Taylor expanding and appealing Weyl's perturbation inequality), we can show that, whilst this operator is not generically negative definite, it will be at sufficiently small p or w.
This concludes the section of the proof addressing terms k = k2 and / * l2 (case 3). Again, these variance contributions are always suppressed by antithetic termination at sufficiently small p or p(U).
Having now considered all the possible variance contributions enumerated by cases 3 and shown that each is either reduced unmodified by the imposition of antithetic termination, can finally conclude that our novel mechanism does indeed suppress the 2-regularised Laplacian kernel estimator variance for d-regular graph of equal weights at sufficiently small p(U).
As mentioned in the main body of the manuscript, these conditions tend not to very restrictive in experiments. Intriguingly, small p(U) with p 1/ actually works very well.
Our next task generalise these results broader classes of graphs.
8.5 Extending the results to arbitrary graphs and sampling strategies (Theorem 3.2 cont.)
Throughout Sec, 8.4 we considered the simplest setting of d-regular graph where all edges have equal weight. We have also taken a basic sampling strategy, with the walker choosing one of its
24

---

current node's neighbours at random at every timestep. Here we relax these assumptions, showing that our results remain in more general settings.
8.5.1 Relaxing d-regularity
First, we consider graphs whose vertex degrees differ. It is straightforward to see that the terms in case
2 (Sec. 8.4.2) are unmodified because taking dtd -> II -,=_1 d, p(w1) and dm - [[."= | di in plus)
exactly compensated by the corresponding change in in joint probability p(w1 E ( 1), Wa (2k2, i)). Our previous arguments all continue to hold.
Case (Sec. (8.43) is only little harder. Now the prefactor in square parentheses in the top line of Eq. 80 evaluates to
2len(ws)
Ien((or) II df
Σ W1E
w
2) plus E 12(k1, i))
(83)
1
-
which is still positive for any node z. The lower line Eq. 80 unmodified because once again the change dt + I di exactly cancels in the marginal and joint probabilities, J is unchanged and our previous conclusions prevail.
8.5.2 Weighted graphs
Now we permit edge weights differ across the graph. Once again, case 2 (Sec. 8.4.2 is straightfor- ward: instead of Eq. [60 we take
Σ (ww1) Σ (U)) jar (),
(84)
W1 Wi
where the weighted adjacency matrix. We incorporate the product of each walk's edge weights into the combinatorial factor, then sum over walk lengths as before. In downstream calculations we drop all instances of ur and reinterpret the eigenvalues of the instead of A, but our arguments otherwise unmodified; these variance contributions will be suppressed if p(U) p is sufficiently small.
Case 3 (Sec. 4.3
is also easy enough; the bracketed prefactor of 80 becomes
2len(wr)
/len(w))
Σ
II
w wf
plws E ())
(85)
1=1
which is again positive. Here, W/-4+1 denotes the weight associated with the edge between the and i 1th nodes of the walk. Therefore, it is sufficient that the matrix J with matrix elements
Xp
[-
Ap
Jpg
(86)
X,) )
1- 2p
Apd
negative definite, with Ap now the pth eigenvalue the weighted adjacency matrix U. Following the same arguments as in Sec. 8.4.3 this will be the case small enough p or ((U).
8.5.3 Different sampling strategies
Finally, we consider modifying the sampling strategy for random walks on the graph. We have previously assumed that the walker takes successive edges at random (i.e. with probability J, ). but the transition probability also be function of the edge weights. For example, if all the edge weights are positive, we might take
pii
j|8)
Wsj
(87)
E
Wik
for the probability of transitioning from node to j at a given timestep (with waj Uz)) given that the walker does not terminate. This strategy increases the probability of taking edges with bigger
25

---

weights and which therefore contribute more to (IN - -U) 2 something that empirically suppresses the variance on the estimator of the 2-regularised Laplacian kernel. Does antithetic termination reduce it further?
Case (Sec. 8.4.2) is again easy; the w-dependent modifications to p(w1) and p(wa) are exactly compensated by adjustments to p(w1 E 92(k), i), wz E S2(k2, 2)). To wit, Eq. 57 becomes
1-2p
(1 p) m
if < m
(ww)
T-p
p(ws E 1) E ((1,t))
(1
2p) m
if n m
(88)
r( )a
Wan [á((w))
1p
(1 (}
if n > m.
((w))7(W))
1-p
where we defined new function of a walk,
7W)
II
Σ
Wik-
(89)
NEW kEN(()
computes the sum of edge weights connected to each node in the walk w (excluding the last), then takes the product of these quantities. It is straightforward check that, when all the graph weights are equal, ([w)) 7(4) am with m the length of w. Meanwhile, p(w1) becomes
p(w1) = (1 - (w1)
(90)
7(11)
such that these modifications cancel when evaluate Eq. 53
Case 3 (Sec. (66) also straightforward. The prefactor in square brackets is equal to 85 and again positive for any valid sampling strategy p(w1 e (2(k), i)) and J does not change, our arguments still hold and these variance contributions are reduced by antithetic coupling.
We note that these arguments will generalise straightforwardly to any weight-dependent sampling strategy and are not particular to the linear case. w/y can be replaced by some more complicated variant that defines valid probability distribution p(w1 (2(k1, i)) and antithetic termination will still prove effective.
8.5.4
Summary
In Sec. 8.5 our theoretical results for antithetic termination have proved robust to generalisations such relaxing d-regularity and changing the walk sampling strategy. A qualitative explanation for this is as follows: upon making the changes, the ratio of the joint to marginal probablities
p(w1, wa)
(91)
p(wi (y(ws)
unmodified, This is because we know how we are modifying the probability over walks and construct the estimator to compensate for it. Meanwhile, the correlations between walk lengths insensitive to the walk directions, so in every case they continue to suppress the kernel estimator variance. The only kink is the terms described in Sec. 8.4.3 which require a little work, but the mathematics conspires that our arguments are again essentially unmodified, though perhaps without such an intuitive explanation.
8.6 Beyond antithetic coupling (proof of Theorem 3.4
Our final theoretical contribution is to consider random walk behaviour when TRVs are offset by less than p, A < p. Unlike antithetic coupling, it permits simultaneous termination. Eqs 13 become
p A
p(s1) = p(s2) p.
p(51) p(52) P.
p(s2(s1)
(92)
1 1-
p(32(51)=
p(s2|51) --
p(32(31)
1
26

---

The probability of two antithetic walks S2(1, i) and 12(3, i) containing subwalks wi and wy becomes
1- p--- A
"
(1 p)
if
dm
I-p
p(ws E () E ()1,))
p
A )m
if
(93)
1-p-A
m
(1 p)
if > m,
which the reader might compare to Eq. 57
In analogy to Eq. 63 this induces the matrix
1 w2 Aky ABs
D
(94)
u Ak, Aks (1-p)=
1-p-A A (1 wd [11 w>
We can immediately observe that this exactly equal to Cky, when A = p(1 11pp, p), so for a pair of walkers with this TRV offset the variance will be identical the i..d. result. Replacing D by D in Egg and Fpg and Jpg and reasoning about negative definiteness via their respective Taylor expansions (as well as the new possible cross-term Du, Akg D kn Ay ), is straightforward conclude that variance is suppressed compared to the i.i.d. case provided A > p(1 p) and p(U) or p is sufficiently small. The + 0 limit demands slightly more careful treatment: in order to stay in the regime
p(1 p) < A < we need to simultaneously take A - 0, e.g. by defining A (p) = p(1 - p) + ap? with the constant 0 < a < 1.
This result was reported in Theorem 3.4 of the main text.
8.7 What about diagonal terms?
The alert reader might remark that all derivations in Sec. 8.4 have taken i * j, considering estimators of the off-diagonal elements of the matrix (IN - U)
fact, estimators of the diagonal elements o(i) (ii) will be biased for both GRFs and q-GRFs if d(i) is constructed using the same ensemble of walkers because each walker is manifestly correlated with, rather than independent of, itself. This rectified by taking two ensembles of walkers out of each node, each of which may exhibit antithetic correlations among itself, then taking the estimator Φ1 (i) (2()). It is straightforward to convince oneself that, in this setup, the estimator is unbiased and q- q-GRFs will outperform GRFs. practice, this technicality has essentially no effect on (q- (Q-)GRF performance and doubles runtime so we omit further discussion.
8.8 Further experimental details: compute, datasets and uncertainties
The experiments in Secs. 4.1 4.2 and 4.4 were carried out on an Intel CMre i5-7640X CPU @ 4.00GHz 4. Each required I CPU hour. The experiments in Sec. 4.3 were carried 2-core
Xeon 2.2GHz with 13GB RAM and 33GB HDD. The computations for the largest considered graphs
took 1 CPU hour.
The real-world graphs and meshes were accessed from Ivashkin [2023] and Dwwson-Haggenn [0033 )
with further information about the datasets available therein. Where were able to locate them. the
original papers presenting the graphs are: /Zachary [1977] Lusseau et al.) (0033
Newman 2006
Bollacker et al. [9988] Leskovee et al. [2007].
All our experiments report standard deviations on the means, apart from the clustering task in Sec. 4.3 because running kemellised k-means large graphs is expensive.
27


Reviews:
[{"review_id": "G4Swhk2Tkv", "reviewer": "Reviewer_wGeX", "summary": "The paper introduces a method for graph random features (GRFs). Graph random features is the graph equivalent (i.e. kernels defined on graphs) of kernel matrix approximation using Random Features (Rahimi & Recht). Recently Chromanski et al. had proposed a method for GRF which uses a series of random walks. \n\nIn the Euclidean setting, the quasi-Monte Carlo variants of the Random Features methods have been developed since they provide better convergence properties. So the goal of this paper is to devise an algorithm for the quasi-Monte Carlo GRFs setting. \n\nThe idea is to correlates the length of the random walks by imposing 'antithetic termination' which is effectively a procedure to obtain a more diverse ensemble of graph random walks. \n\nThe paper provides several theoretic results related to the error guarantees as well as a brief empirical study.", "questions": "* It seems that all the experiments are on synthetic datasets. Can the authors confirm that or clarify why they didn't test this on at least on real dataset?\n\n* Sec 4.3 cites the Dhillon paper regarding the kernel matrix they used but that paper has experiments on both handwritten image data and a gene dataset. It is not clear what dataset this experiment is based on. Please clarify this point.\n\n* Additionally,  I would suggest adding a sentence in the text describing the nature of the datasets so that the reader doesn't have to go back to the original papers to understand the eval/test data. Similarly for Sec 4.4 what are the matrices used by Dawson-Haggerty, 2023? Please add a 1 sentence clarification. ", "limitations": "The authors address the limitations and potential downstream impact in a dedicated section at the end.", "rating": 9, "confidence": 4, "soundness": 3, "presentation": 3, "contribution": 4, "strengths": "+ Original idea with strong practical impact and wide ranging downstream applications as time-efficient approximation of the graph diffusion process is used in many scientific fields. \n+ The paper is of very high quality: Well written and clear.  \n+ Good literature review and problem justification, as well as empirical analysis. \n\n", "weaknesses": "- Some typos/grammatical mistakes (esp. definite/indefinite articles), please proof read before publishing\n- Some questions below as well.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the reviewer for their reading of the manuscript and positive comments. We are delighted that they recognise the scheme’s originality, strong practical impact and broad downstream applications. They note that the work is well-motivated and the manuscript is clear. We address their questions in detail below and will be happy to follow up with any further clarifications.\n\n1. *Typos*: we thank the reviewer for flagging minor errors and will be sure to check the text again prior to submission of the final draft. \n\n2. *Synthetic and real datasets*: we refer to a graph as ‘synthetic’ if it is generated according to some mathematical scheme (deterministic or stochastic) and ‘real-world’ if it is derived from actual data. By this definition, the experiments in sections 4.1 (direct estimation of the $2$-regularised Laplacian kernel) and 4.2 (simulation of graph diffusion) use an equal split of synthetic (2x ER, binary tree and ladder) and real-world (karate, dolphins, football, eurosis) graphs. Section 4.3 (clustering) uses exclusively real-world graphs, all of which are standard for benchmarking. Section 4.4 (kernel regression) uses triangular mesh graphs from a popular online repository, corresponding to a variety of 3D shapes such as a cylinder and a teapot. It is unclear whether such computer graphics should be considered synthetic or real but they are certainly standard. We hope this clarifies the split of synthetic and real data but will be happy to discuss further if we misinterpret the reviewer’s question.\n\n3. *Relationship to Dhillon et al. (2004)*: we cite Dhillon’s paper because it provides a clear description of kernelised $k$-means clustering (Alg. 1 in Sec. 2.1, page 552). In our work, we implement this algorithm to cluster the nodes of $8$ real world graphs, comparing the results when we use the exact kernel (yielding the groundtruth clusters) and its approximation with i) GRFs and ii) q-GRFs. We find that, since q-GRFs provide a provably better approximation of the Gram matrix, they tend to return clusters closer to the groundtruth. We stress that all graphs considered correspond to real-world data and are of standard use in benchmarking. The reviewer correctly notes that Dhillon’s paper uses different datasets to us in experiments designed to elucidate the relationship between spectral clustering and kernelised $k$-means clustering, but these are not of direct relevance to our work. \n\n4. *Clarification of nature of datasets, especially Dawson-Haggerty (2023)*: we thank the reviewer for their comment and direct them back to point 2 of our rebuttal. Dawson-Haggerty (2023) is a popular repository of triangular mesh graphs corresponding to 3D computer graphic objects, and here the task is to use kernel regression to predict the missing normal vector directions of the nodes. We will be sure to clarify this important point in the manuscript to save the reader from needing to follow the link to the repository. \n\nWe again thank the reviewer for their attentive reading of the draft and excellent review.  "}, {"review_id": "n5tGWUia7y", "reviewer": "Reviewer_SAEb", "summary": "This work proposes an efficient random-walk sampling approach that accelerates the estimation of the feature mapping for 2-regularized Laplacian kernels. The random walk sampling is based on antithetic termination, which is a variance reduction technique and sounds novel  when being applied to the estimation of graph random features. The technique sounds solid. The authors also conducted experiments from different angles to verify the effectiveness of the proposed approach in practice. ", "questions": "1. Why is the 2-regularized Laplacian kernel important? Adding experiments to compare with other graph kernels (in the aspects of both inference accuracy and computation) may improve the draft. \n2. Is the random walk approach only applied to the estimation of 2-regularized Laplacian kernel?", "limitations": "N/A", "rating": 5, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "1. Using random walk with antithetic termination to estimate graph random features is novel, to the best knowledge of me. \n2. Although I did not check the appendix, the technique and the benefits it may bring look reasonable to me.\n3. The authors did a great job to explain the key insights behind. The paper is very well written. \n4. The experiments are extensive. The benefits of the proposed approach is justified from various angles. ", "weaknesses": "1. The biggest weakness of the work is the limited scope of this work. It is unclear why the 2-regularized Laplacian kernel is important and why we need to compute it. It seems that 1-regularized Lap kernel can be efficient computed by Pagerank algorithms, why we may need 2-regularized Lap kernel. \n2. It is unclear how much the estimation of the 2-regularized Laplacian kernel relies on the random walk with antithetic termination. Can this random walk approach be only applied to 2-regularized Laplacian kernels or a broader range of graph kernels? This question is kind of related to whether the paper has over-claimed the contributions. Neither the title nor the introduction reflect that the improvement is only for  2-regularized Laplacian kernels, and they claim a broader range of contributions. I expect the authors to make the actual contributions more clear in the revised version.  ", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We sincerely thank the reviewer for their encouraging comments, especially pertaining to the novelty of the mechanism, the quality of the writing and the exhaustive nature of the experimental section. Below, we address minor points of misunderstanding and answer the reviewer’s questions. \n\n1. *Scope of the work*: though it is our predominant focus, **the scope of antithetic termination is not limited to estimation of the $2$-regularised Laplacian**. As we show explicitly in Sec. 8.1, GRFs (and, by extension, our more efficient q-GRF variant) can be used to construct scalable estimates to $d$-regularised Laplacians for arbitrary integer $d$, covering a broad family of very popular graph kernels. Moreover, in Sec. 2.1 we show how these techniques can be used to estimate the heat kernel $K_\\textrm{diff} = \\textrm{exp}(- \\widetilde{\\textbf{L}} t)$ – a different graph kernel which is also widely-used – and in Table 1 we see that here the gains from employing antithetic termination are substantial. This covers what are often considered the two most popular classes of graph kernels, listed as the top examples in Eqs 17 and 18 of Smola and Kondor’s seminal paper (https://people.cs.uchicago.edu/~risi/papers/SmolaKondor.pdf). With strong physical motivation and practical efficacy, they are well-documented to enjoy widespread applications in any algorithms employing kernel-based algorithms on graphs, e.g. SVMs, Gaussian Processes, and kernelised PCA (e.g. https://academic.oup.com/bioinformatics/article/21/suppl_1/i47/202991, https://arxiv.org/abs/1403.3148, https://link.springer.com/chapter/10.1007/11731139_28). Exciting recent work also suggests that they might also be used to modulate regular attention to provide a structural inductive bias in graph-based Transformers (https://arxiv.org/abs/2106.05234, https://arxiv.org/abs/1710.10903, https://arxiv.org/abs/2107.07999). We hope the reviewer will agree that this demonstrates the broad scope of q-GRFs for kernel estimation, and will be sure to make this point clearer in the manuscript. What’s more, we would like to stress that **antithetic termination is fundamentally a procedure to sample more diverse discrete random walks which can be straightforwardly applied in Monte Carlo schemes beyond graph kernel estimation**. This further extends our work’s scope. We discuss this in detail later in the rebuttal (see point 4).\n\n2. *Experimental comparison of graph kernels*: we again stress that the manuscript already includes experimental and theoretical discussion of different types of kernel beyond the $2$-regularised Laplacian (see Secs 2.1, 4.1 and 8.1). However, the purpose of the paper is not to compare how they perform on particular learning tasks, but rather present a general QMC scheme that can better the estimation of any one of them. Concretely, we are less interested in whether e.g. the $1$- or $2$-regularised Laplace kernel is better for node attribute prediction on trimesh graphs, but rather in whether approximating either of them with our new QMC scheme can improve efficiency and performance compared to the vanilla MC mechanism.\n\n3. *Contributions, title and introduction*: as explained above, GRFs can be used for scalable estimation of a broad range  of kernels beyond the $2$-regularised Laplacian, and our theoretical and experimental results show that in every case our new q-GRFs variant converges better. Whilst we anticipate that the popular $d=2$ case might be of most direct interest, we agree that this is an important point and will emphasise it. \nOur paper’s title uses standard terminology: Choromanski coined ‘GRFs’ to describe their new technique for estimation of graph-based kernels, and ‘q-GRFs’ is a natural choice for our more efficient QMC variant. As early as the third sentence of the abstract we state that our chief interest will be the $2$-regularised Laplacian. We are transparent about this fact. That said, we again stress that the method’s applicability is more general: it is at its heart a technique to improve sampling on graphs.   \n4. *How are (q-)GRFs related to PageRank*? We sincerely thank the reviewer for referring us to PageRank, a popular measure of the importance of the nodes of a network. The PageRank vector is the steady state distribution of a particular random walk process, or equivalently the top eigenvector of a $1$-regularised Laplacian kernel. The reviewer correctly points out that efficient MC algorithms exist to numerically approximate this vector, but the scope of GRFs (and the q-GRFs variant we present) is broader: they provide a means to approximate the entire $d$-regularised Laplace kernel for arbitrary integer $d$, rather than just its largest eigenvector when $d=1$.\nThat said, we agree that the literature on PageRank presents interesting connections to our work. Especially relevant is Alg. 1 of Avrachenkov’s ‘Monte Carlo Methods in PageRank Computation’ (https://www.jstor.org/stable/40232890), which samples $K$ random walks out of each node with probability of termination $p$ at every timestep and records where each walk ends. It is straightforward to incorporate antithetic termination into Avrachenkov’s algorithm, pairing up RWs and using our mechanism to anticorrelate their lengths. **Following the reviewer’s comment, we implemented this experiment and observed large gains: for a fixed number of walkers, antithetic termination  suppresses the $L_2$-norm error on the approximation of the PageRank vector by a factor of approximately 2** (see the attached PDF). We have added this remarkable result to the manuscript and will include accompanying theoretical discussion in the final draft. This further demonstrates the scope of antithetic termination beyond graph-based kernels. \n\nWe again thank the reviewer for their helpful comments and questions. We hope the above clarifications and additions will be sufficient to allay any concerns and invite them to reach out with further questions."}, {"review_id": "8rr4ix2bxX", "reviewer": "Reviewer_mBGP", "summary": "The author(s) extend the work of Choromanski [2023] and introduce quasi monte carlo graph random features.\nTo estimate the 2-regularized Laplacian kernel, they now use random walk features for each node which have coupled probabilities of terminating after a certain number of steps.\nThis increases diversity in the features and interestingly leads to provable reduction of variance of the global estimates.\nAn empirical evaluation supports the theoretical findings by improved practical variance, as well as in a clustering and supervised node learning scenario.\n", "questions": "\nThe last sentence of the proof sketch seems to encompass quite some work. I suggest it to be extended to mention at least what kind of matrix is supposed to be negative definite.\n\nSuggestions \nl. 107: 'the the'\nl. 164: 'the event that event that'\nl. 207 and l. 222: 'path' should probably be 'walk'\nReferences: I suggest the author(s) revise the references for consistent formatting. E.g. Arxiv vs Corr.Abs., as another example, Choromanski [2023] does not list any venue, at all.\n", "limitations": "The authors provide a dedicated section discussing potential (far away) societal impacts, as well as open questions with the current approach. \nThe discussion is adequate.\n", "rating": 7, "confidence": 4, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "\nThe paper is well-written and mostly easy to follow. \nWhile proofs are moved to the appendix, a sketch of proof of the main result is included in the paper. I like that.\n\nThe introduced way of increasing diversity in random walks is (up to my knowledge) novel and very easy to implement. While it is evaluated in the context of estimating the 2-regularized laplacian of a graph, I agree with the authors that this method is of independent interest. \nAn interesting question is if this kind of sampling improves results of other random walk based approaches, as well, or if the parctical impact is restricted to just a few scenarios (GRF being one of them).\nAs a result, I see potential impact for further research beyond this immediate application.\n\n", "weaknesses": "The introduction is structurally very similar to that of Choromanski [2023]. I'm assuming now that I know one of the authors of this submission. Otherwise, I would recommend the authors to revise the introduction to decrease similarity. \n", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the reviewer for their careful reading of the manuscript and thoughtful feedback. We are delighted that they note this novel mechanism to increase the diversity of random walks is simple to implement and offers strong theoretical guarantees and empirical performance – in fact not only for the $2$-regularised Laplacian but also the $d$-regularised Laplacian (Sec. 8.1) and the heat kernel (Secs 2.1 and 4.1). They note that the method might be of independent interest for other random walk-based algorithms,  a possibility that we expect to dedicate substantial future work to. Prompted by comments from another reviewer, we have already implemented a new experiment applying antithetic termination to numerical estimation of the PageRank vector (a popular measure of a node’s importance in a graph – https://www.jstor.org/stable/40232890). We observe very large gains, showing the generality of the method beyond graph kernel estimation. Please see the attached PDF. It seems this simple scheme to improve random walk sampling efficiency is of very broad scope.  \n\nWe address the reviewer’s questions and suggestions in detail below.\n\n1. *Structure of introduction*: we thank the reviewer for their observation, but will refrain from further comment in the interests of preserving anonymity.\n2. *Final sentence of proof sketch*: the reviewer correctly observes that a key ingredient of our theoretical work is identifying a particular matrix that must be negative definite in order that antithetic termination suppresses the kernel estimator variance, and using Weyl’s perturbation inequality to identify regimes in which this is the case. As they say, this belies substantial technical work and, whilst the proof sketch is not intended to reproduce every detail, we agree that supplementing the discussion might benefit the motivated reader. We have updated the manuscript to include this.  \n3. *Typos and reference formatting*: we thank the reviewer for their careful observations and have corrected these minor mistakes in the manuscript.\n\nWe again thank the reviewer for their close reading of the manuscript and helpful suggestions. \n"}, {"review_id": "LsS2J2Niew", "reviewer": "Reviewer_qCjY", "summary": "Kernel methods are important in many ML applications, but they suffer from scalability issues. In the Euclidean setting, random features methods (Rahimi & Recht NeurIPS 2017) can be used to \"sketch\" the kernel matrices and speed up computations. In the graph setting, however, no such random feature seems to have been available until a 2023 preprint by Choromanski. \n\nThe present paper present an improvement over Choromanski's method. To explain it, note that the basic idea in the 2023 preprint is to consider random features for the regularized graph Laplacian, and to obtain an approximation of an operator of the form $(I-U)^{-d}$ via importance-weighted random walks. These random walks have geometrically distributed lengths. \n \nThe main theoretical contribution of the present paper can now be explained: a simple modification of Choromanski's method, where the the geometrics of pairs of random walks become correlated, has the same zero bias and a smaller variance than the above method. \nThis requires an elementary, but long proof with a long case analysis.\n\nExperimentally, the gains observed were:\n\n* visible, but not too big for approximation of the kernel;\n* impressive for graph diffusions (which is natural because this requires exponentiating the Laplacian);\n* fairly significant for clustering and node attribute prediction methods.", "questions": "Can you point to more specific and realistic applications where this methodology would be useful? (A somewhat speculative discussion would be nice.)\n\n--\n\nEdit (08/11) -- addressed in rebuttal. ", "limitations": "These are properly discussed. ", "rating": 7, "confidence": 3, "soundness": 4, "presentation": 3, "contribution": 3, "strengths": "The paper is clear and seemingly original. In general it is fairly well-written. Moreover, the \"anthitetic\" coupling is a simple little idea that could be applicable more broadly.\n\n", "weaknesses": "The present paper has been submitted before Choromanski's preprint appeared. As a result, applications of graph kernel random sketches are quite incipient still. It is hard to say whether the present paper will have a significant direct impact. _[Edit on 08/11: this was partly addressed in the rebuttal.]_\n\nThe proof of variance reduction is very long and difficult to check. (Though I believe it.)  (An \"actionable version\" of this point is that it would pay off to make this paper a bit less dependent on Chorumovski's preprint by eg. adding more experiments or theory.)  \n\n\n_Less important comments_\n\nThe notation $t_{1,2}$ is a bit weird as it represents a pair of RVs. Moreover, if I understand correctly, these r.v.s are resampled at each step of the walk. \n\nIn terms of exposition, it seems that all the authors need to do is to correlate the coin flips for the two walks at each step. Could you write their joint distribution purely in terms of the probabilities of 00, 01, 10 and 11?  ", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the reviewer for reading the manuscript, and are delighted that they note the work’s novelty, clarity and strong experimental contributions. We agree that it is remarkable that such a conceptually simple technique to correlate walker lengths can lead to such rich behaviour and demand such a long proof. We hope the brief proof sketch will help explain the high-level strategy to the interested reader. \n\nWe also agree that antithetic coupling is of broader interest for algorithms sampling random walks with geometrically distributed lengths. In fact, prompted by such comments from the reviewers, we have since implemented antithetic termination for numerical approximation of the PageRank vector, a popular measure of the importance of a node in a graph (https://www.jstor.org/stable/40232890). We find very large gains from our QMC scheme, reducing the $L_2$-norm error on the vector’s estimation by a factor of approximately 2. This further demonstrates the power and generality of our new sampling strategy. Plots are included in the attached PDF.   \n\nWe now address the reviewer’s questions and comments in detail. \n\n1. *Relationship with Choromanski’s paper*: the reviewer correctly notes that the field of graph kernel random sketches is young and the impact of these mechanisms will become clearer with time. That said, we note that since submitting this draft Choromanski’s original preprint has been accepted to ICML with an oral presentation, suggesting that these scalable GRF techniques — and by extension their more efficient QMC variants — might be of interest to the community. The present manuscript also goes well beyond Choromanski’s work in many ways, including: considering a *quasi*-Monte Carlo scheme, using (q-)GRFs to approximate the heat kernel (Secs 2.1 and 4.1) and providing detailed theoretical analysis of both the old and new schemes (Appendices 8.3-8.6). Indeed, the $d$-regularised Laplacian and heat kernels are possibly the two most popular classes of kernels defined on graph nodes, listed as top examples in Eqs 17 and 18 Smola and Kondor’s seminal paper (https://people.cs.uchicago.edu/~risi/papers/SmolaKondor.pdf). We therefore anticipate that these new schemes to efficiently approximate them will prove very impactful. We will emphasise this in the manuscript.\nMoreover, we stress that antithetic termination is at its heart a scheme to sample more diverse random walks. As such, its applicability extends beyond the (nonetheless important) application of q-GRFs. This is made clear by our new PageRank results. We are eager to hear of any further suggestions for possible applications.\n2. *Notation for termination random variables (TRVs)*: we thank the reviewer for their stylistic comment and will clarify the notation in the manuscript. They are correct that the TRVs are resampled at every timestep — an important point that we will highlight.\n3. *Joint distribution over coin flips*: we thank the reviewer for their interesting comment. In Eq. 13 we opted to present the conditional probabilities. From here it is trivial to obtain the joint probabilities, but we agree that explicitly including them might make the exposition clearer. We have updated the manuscript.\n4. *Applications where this methodology might be useful*: to quote Smola and Kondor’s seminal paper  (https://people.cs.uchicago.edu/~risi/papers/SmolaKondor.pdf), ‘since kernel-based algorithms, such as Support Vector Machines, Gaussian Processes, Kernel PCA, etc. capture the structure of $\\mathcal{X}$ via the kernel $K : \\mathcal{X} \\times \\mathcal{X} \\to \\mathbb{R}$, as long as we can define an appropriate kernel on our discrete input space, these algorithms can be imported wholesale, together with their error analysis, theoretical guarantees and empirical success’. Since any graph-based learning algorithm using the $d$-regularised and heat kernels can benefit from our efficient q-GRF construction, it follows that our techniques will enjoy applications in a very wide variety of downstream graph-based learning tasks. To give some concrete examples: inferring protein function (https://academic.oup.com/bioinformatics/article/21/suppl_1/i47/202991), community detection (https://arxiv.org/abs/1403.3148) and recommender systems (https://link.springer.com/chapter/10.1007/11731139_28). More recently, the heat kernel has been used in settings as diverse as manifold learning for deep generative modelling (https://arxiv.org/abs/2010.01761) and for solving single- and multiple-source shortest path problems (https://arxiv.org/abs/2010.01761). Speculatively, these kernels might also be used to modulate the regular attention mechanism in graph-based Transformers, incorporating topological signal and endowing a structural inductive bias (https://arxiv.org/abs/2106.05234, https://arxiv.org/abs/1710.10903, https://arxiv.org/abs/2107.07999). For big graphs where exact kernel evaluation is expensive, we expect that our efficient algorithm will prove especially impactful. \nBeyond kernel estimation, we have argued that antithetic termination is likely to be effective in any Monte Carlo algorithm sampling random walks with geometrically distributed lengths. Our new results confirm that this is the case for numerical estimation of the PageRank vector, immediately opening up diverse applications from ranking websites to analysing social networks. \nWe agree that it is important to properly motivate the work and will add this discussion to the draft.    \n\nWe again sincerely thank the reviewer and invite them to respond with any further questions. \n"}]

Author rebuttal:
We are grateful for the reviewers’ careful readings of the manuscript and thoughtful feedback. We are delighted to receive such positive comments about the elegance of the mechanism and our exhaustive theoretical and empirical analysis. They also make several helpful suggestions which have further improved the draft, in particular prompting us to include a new section applying antithetic termination to PageRank (please see our individual responses for details and the attached PDF for results). This further showcases the strength and generality of the QMC approach even beyond graph kernel estimation. We address concerns and points of minor misunderstanding in the rebuttal and warmly encourage the reviewers to respond with any other queries. 

Meta-review: The reviewers were in consensus that the antithetic termination idea introduced in this paper is simple, novel, empirically effective, and may have further applications beyond those explored in the work. All felt that the paper was also well-written. Thus, the paper was judged to be a clear accept.
Decision: Accept (spotlight)