# Quasi-Monte Carlo Graph Random Features

Isaac Reid
University of Cambridge
<EMAIL_ADDRESS>
Krzysztof Choromanski
Google
Columbia University
<EMAIL_ADDRESS>
Adrian Weller
University of Cambridge
Alan Turing Institute
<EMAIL_ADDRESS>Senior lead.

###### Abstract

We present a novel mechanism to improve the accuracy of the recently-
introduced 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 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.111We will
make all code publicly available.

##  1 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 a mathematically principled and elegant way to perform
nonlinear inference using linear learning algorithms. The positive definite
_kernel function_ $k:\mathcal{X}\times\mathcal{X}\to\mathbb{R}$, defined on an
input domain $\mathcal{X}$, measures the ‘similarity’ between two datapoints.
Examples in 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 a 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, 2007,
2008; Avron et al., 2017; Liu et al., 2022) . These randomised mappings
$\phi:\mathbb{R}^{d}\to\mathbb{R}^{s}$ construct low-dimensional feature
vectors whose dot product equals the kernel evaluation in expectation:

$k(\bm{x},\bm{y})=\mathbb{E}\left(\phi(\bm{x})^{\top}\phi(\bm{y})\right).$ (1)

This permits a 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 Euclidean kernels with properties engineered
for the desiderata and symmetries of the particular kernel being approximated
(Dasgupta et al., 2010; Johnson, 1984; Choromanski et al., 2020; Goemans and
Williamson, 2004; 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:\mathcal{V}\times\mathcal{V}\to\mathbb{R}$ on the set of nodes
$\mathcal{V}$ of a graph $\mathrm{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). Substantial
research effort has also been devoted to developing and analysing graph
kernels $k:\mathcal{G}\times\mathcal{G}\rightarrow\mathbb{R}$, now taking
entire graphs $\mathrm{G}\in\mathcal{G}$ from graph spaces $\mathcal{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 $\mathbf{K}$ is typically of at least
cubic time-complexity in the number of nodes $N$. In contrast to kernels
defined on points in $\mathbb{R}^{d}$, random feature methods for fixed graph
kernels (c.f. kernel learning (Fang et al., 2021)) have proved challenging to
construct. Only recently has a viable _graph random feature_ (GRF) mechanism
been proposed, which uses a series of random walkers depositing ‘load’ at each
node they pass through (Choromanski, 2023). GRFs provide a low-rank unbiased
estimate of the matrix $(\mathbf{I}_{N}-\mathbf{U})^{-d}$, where $\mathbf{U}$
is a weighted adjacency matrix of the graph and $d\in\mathbb{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, _orthogonal 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 kernel in
the _orthogonal Johnson-Lindenstrauss transformation_ (Choromanski et al.,
2017). With positive random features, the recently-introduced class of
_simplex random features_ (SimRFs) performs even better (Reid et al., 2023).
This has been used to great effect in estimating the attention mechanism of
Transformers (Vaswani et al., 2017), overcoming its prohibitive quadratic
time-complexity scaling with token sequence length.

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 (2023) first identified the challenge of
constructing _quasi-Monte Carlo GRFs_ (q-GRFs). They suggested a 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 as ‘q-GRFs’ for
brevity) will spur further research on correlations between graph random walks
in machine learning.

We emphasise that, whilst we have presented antithetic termination through the
lens of q-GRFs (an analytically tractable and important use case), it is
fundamentally a procedure to obtain a more diverse ensemble of graph random
walks. It may be of independent interest, e.g. for estimating graphlet
statistics for kernels between graphs (Chen et al., 2016; Wu et al., 2019;
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 it 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.3. We conduct an exhaustive set of experiments in Sec. 4 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 in 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$ is defined by a set of vertices
$\mathcal{V}$ enumerated $1$ to $N$ and a set of edges $E$ given by the
unordered vertex pairs $\\{v_{i},v_{j}\\}$ where $v_{i}$ and $v_{j}$ are
neigbours (denoted $i\sim j$), themselves associated with weights
$W_{ij}\in\mathbb{R}$. The _weighted adjacency matrix_
$\mathbf{W}\in\mathbb{R}^{N\times N}$ has matrix elements $W_{ij}$: that is,
the associated edge weights if $i\sim j$ and $0$ otherwise.

Denote by $\mathbf{D}\in\mathbb{R}^{N\times N}$ the diagonal matrix with
elements $D_{ii}\coloneqq\sum_{j}W_{ij}$, the sum of edge weights connecting a
vertex $i$ to its neighbours. The _Laplacian_ of $G$ is then defined
$\mathbf{L}\coloneqq\mathbf{D}-\mathbf{W}$. The _normalised Laplacian_ is
$\widetilde{\mathbf{L}}\coloneqq\mathbf{D}^{\frac{1}{2}}\mathbf{L}\mathbf{D}^{-\frac{1}{2}}$,
which rescales $\mathbf{L}$ by the (weighted) number of edges per node.
$\mathbf{L}$ and $\widetilde{\mathbf{L}}$ share eigenvectors in the case of a
$d$-regular graph and play a central role in spectral graph theory; their
analytic properties are well-understood (Chung, 1997).

In classical physics, diffusion through continuous media is described by the
equation

$\frac{d\mathbf{u}}{dt}=\nabla^{2}\mathbf{u}$ (2)

where $\nabla^{2}=\frac{\partial^{2}}{\partial
x_{1}^{2}}+\frac{\partial^{2}}{\partial
x_{2}^{2}}+...\frac{\partial^{2}}{\partial x_{N}^{2}}$ is the _Laplacian
operator_ on continuous spaces. The natural analogue on discrete spaces is
$\mathbf{L}$, where we now treat $\mathbf{L}$ as a linear operator on vectors
$\mathbf{u}\in\mathbb{R}^{N}$ (or equivalently functions $u:V\to\mathbb{R}$).
This can be seen by noting that, if we take $\mathbf{W}$ to be the
_unweighted_ adjacency matrix,
$\left<\mathbf{u},\mathbf{L}\mathbf{u}\right>=\mathbf{u}^{\top}\mathbf{L}\mathbf{u}=-\frac{1}{2}\sum_{i\sim
j}(u_{i}-u_{j})^{2}$ so, like its continuous counterpart, $\mathbf{L}$
measures the the local smoothness of its domain. In fact, in this case
$-\mathbf{L}$ is exactly the finite difference discretisation of $\nabla^{2}$
on a square grid in $N$-dimensional Euclidean space. This motivates the
_discrete heat equation_ on $G$,

$\frac{d\mathbf{u}}{dt}=-\widetilde{\mathbf{L}}\mathbf{u},$ (3)

where we followed literature conventions by using the normalised variant whose
spectrum is conveniently contained in $[0,2]$ (Chung, 1997). This has the
solution $\mathbf{u}_{t}=\exp(-\widetilde{\mathbf{L}}t)\mathbf{u}_{0}$, where
$\exp(-\widetilde{\mathbf{L}}t)=\lim_{n\to\infty}\left(1-\frac{t\widetilde{\mathbf{L}}}{n}\right)^{n}$.
The symmetric and positive semi-definite matrix

$\mathbf{K}_{\textrm{diff}}(t)\coloneqq\exp(-\widetilde{\mathbf{L}}t)$ (4)

is referred to as the _heat kernel_ or _diffusion kernel_ (Smola and Kondor,
2003). The exponentiation of the generator $\widetilde{\mathbf{L}}$, which by
construction captures the _local_ structure of $G$, leads to a kernel matrix
$\mathbf{K}_{\textrm{diff}}$ which captures the graph’s _global_ structure.
Upon discretisation of Eq. 3 with the backward Euler step (which is generally
more stable than the forward), we have that

$\mathbf{u}_{t+\delta t}=(\mathbf{I}_{N}+\delta
t\widetilde{\mathbf{L}})^{-1}\mathbf{u}_{t},$ (5)

where the discrete time-evolution operator
$\mathbf{K}_{\textrm{lap}}^{(1)}=(\mathbf{I}_{N}+\delta
t\widetilde{\mathbf{L}})^{-1}$ is referred to as the _$1$ -regularised
Laplacian kernel_. This is a member of the more general family of
$d$-regularised Laplacian kernels,

$\mathbf{K}_{\textrm{lap}}^{(d)}=(\mathbf{I}_{N}+\delta
t\widetilde{\mathbf{L}})^{-d},$ (6)

for which we can construct an unbiased low-rank approximation using GRFs
(Choromanski, 2023). We predominantly consider $d=2$ with the understanding
that estimators for other values of $d$ are straightforward to obtain (App.
8.1). This demonstrates the intimate connection between the graph-diffusion
and Laplacian kernels, and how a QMC scheme that improves the convergence of
our randomised estimator of $\mathbf{K}_{\textrm{lap}}^{(d)}$ will permit more
accurate simulation of diffusion on a graph.

###  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 a full discussion, but for convenience we
provide a cursory summary.

Suppose we would like to estimate the matrix
$\left(\mathbf{I}_{N}-\mathbf{U}\right)^{-2}$, with
$\mathbf{U}\in\mathbb{R}^{N\times N}$ a weighted adjacency matrix of a graph
with $N$ nodes and no loops. (Choromanski, 2023) introduces a novel algorithm
to construct a low-rank decomposition of this matrix. The author uses _graph
random features_ (GRFs) $\phi(i)\in\mathbb{R}^{N}$, with $i$ the index of one
of the nodes, designed such that

$\left(\mathbf{I}_{N}-\mathbf{U}\right)^{-2}_{ij}=\mathbb{E}\left(\phi(i)^{\top}\phi(j)\right).$
(7)

They construct $\phi(i)$ by taking $m$ random walks
$\\{\bar{\Omega}(k,i)\\}_{k=1}^{m}$ on the graph out of node $i$, depositing a
‘load’ at 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 $p$ at every timestep. The $x$th
component of the $i$th random feature is given by

$\bm{\phi}(i)_{x}\coloneqq\frac{1}{m}\sum_{k=1}^{m}\sum_{\omega_{1}\in\Omega_{ix}}\frac{\widetilde{\omega}(\omega_{1})}{p(\omega_{1})}\mathbb{I}(\omega_{1}\in\bar{\Omega}(k,i)),$
(8)

where: $k$ enumerates $m$ random walks we sample out of node $i$; $\omega_{1}$
denotes a particular walk from the set of all walks $\Omega_{ix}$ between
nodes $i$ and $x$; $\widetilde{\omega}(\omega_{1})$ is the product of weights
of the edges traversed by the walk $\omega_{1}$; $p(\omega_{1})$ is the
marginal probability that a walk contains a subwalk $\omega_{1}$, given by
$p(\omega_{1})=((1-p)/d)^{\textrm{len}(\omega_{1})}$ in the simplest case of a
walk of length $\textrm{len}(\omega_{1})$ on a $d$-regular graph;
$\mathbb{I}(\omega_{1}\in\bar{\Omega}(k,i))$ is an indicator function that
evaluates to $1$ when the walk $\omega_{1}$ is a subwalk of the $k$th random
walk sampled from $i$ (itself denoted $\bar{\Omega}(k,i)$) and is $0$
otherwise.

It is simple to see how Eq. 8 satisfies Eq. 7. By construction

$\mathbb{E}\left[\mathbb{I}(\omega_{1}\in\bar{\Omega}(k,i))\mathbb{I}(\omega_{2}\in\bar{\Omega}(l,j))\right]=p(\omega_{1})p(\omega_{2})$
(9)

for independent walks, whereupon

$\begin{multlined}\mathbb{E}\left(\bm{\phi}(i)^{\top}\bm{\phi}(j)\right)=\sum_{x\in\mathcal{V}}\sum_{\omega_{1}\in\Omega_{ix}}\sum_{\omega_{2}\in\Omega_{jx}}\widetilde{\omega}(\omega_{1})\widetilde{\omega}(\omega_{2})=\sum_{\omega\in\Omega_{ij}}(\text{len}(\omega)+1)\widetilde{\omega}(\omega)=\left(\mathbf{I}-\mathbf{U}\right)^{-2}_{ij}.\end{multlined}\mathbb{E}\left(\bm{\phi}(i)^{\top}\bm{\phi}(j)\right)=\sum_{x\in\mathcal{V}}\sum_{\omega_{1}\in\Omega_{ix}}\sum_{\omega_{2}\in\Omega_{jx}}\widetilde{\omega}(\omega_{1})\widetilde{\omega}(\omega_{2})=\sum_{\omega\in\Omega_{ij}}(\text{len}(\omega)+1)\widetilde{\omega}(\omega)=\left(\mathbf{I}-\mathbf{U}\right)^{-2}_{ij}.$
(10)

This shows us that the estimator is unbiased. The central contribution of this
work is a QMC scheme that induces correlations between the $m$ walks out of
each node to suppress the variance of the estimator $\phi(i)^{\top}\phi(j)$
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 _termination random variables_ (TRVs) between $0$ and $1$ from a
uniform distribution, $t_{1,2}\sim\textrm{Unif}\left(0,1\right)$. Each walker
terminates if its respective TRV is less than $p$, $t_{1,2}<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 as $t_{1,2}=\mathrm{Unif}(0,1)$ but are offset by $\frac{1}{2}$,

$t_{2}=\mathrm{mod}_{1}\left(t_{1}+\frac{1}{2}\right),$ (11)

such that we have the conditional distribution

$p(t_{2}|t_{1})=\delta\left(\mathrm{mod}_{1}(t_{2}-t_{1})-\frac{1}{2}\right).$
(12)

Computational cost: we note that the computational cost of generating an
antithetic TRV according to Eq. 11 is no greater than the cost of generating
an independent TRV, so this drop-in replacement in the GRF algorithm is cheap.
We provide a schematic in Fig. 1(a).

Since the marginal distributions over $t_{i}$ are unchanged our estimator
remains unbiased, but the couplings between TRVs lead to statistical
correlations between the walkers’ terminations. Denoting by $s_{1}$ the event
that event that walker $1$ terminates at some timestep, $s_{2}$ the event that
walker $2$ terminates and $\bar{s}_{1,2}$ their complements, it is
straightforward to convince oneself that for $p\leq\frac{1}{2}$

$\displaystyle p(s_{1})=p(s_{2})=p,\hskip
14.22636ptp(\bar{s}_{1})=p(\bar{s}_{2})=1-p,\hskip
14.22636ptp(s_{2}|s_{1})=0,$ (13) $\displaystyle p(\bar{s}_{2}|s_{1})=1,\hskip
14.22636ptp(s_{2}|\bar{s}_{1})=\frac{p}{1-p},\hskip
14.22636ptp(\bar{s}_{2}|\bar{s}_{1})=\frac{1-2p}{1-p}.$

This termination coupling modifies the joint probability over walk lengths. In
the i.i.d. scheme, the walks are independent and are of expected length

$\mathbb{E}(\textrm{len}(\omega))=\frac{1-p}{p}.$ (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

$\mathbb{E}(\textrm{len}(\omega_{2})|\textrm{len}(\omega_{1})=m)=\frac{1-2p}{p}+2\left(\frac{1-2p}{1-p}\right)^{m},$
(15)

which we derive in App. 8.2. It is 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. 1(b) for a schematic. We refer to random features constructed with
antithetic walkers as _quasi-Monte Carlo graph random features_ (q-GRFs).

antithetic_schematic

(a)

graph_figure_2

(b)

Figure 1: _Left_ : schematic of the i.i.d. (GRF) and antithetic (q-GRF)
mechanisms in termination space. $f(t)$ is the probability density of the
termination random variable (TRV) $t$. Vertical arrows represent draws of $t$,
with the walker terminating if they lie in the pink region where $t<p$. With
q-GRFs the TRVs are offset by $\frac{1}{2}$, modifying the joint distribution
over walk lengths. _Right_ : demonstration with $4$ random walks on the karate
graph, beginning at some node labelled $i$. The blue pair of antithetic walks
(q-GRFs) have very different lengths; they cannot terminate simultaneously.
The red pair of i.i.d. walks (GRFs) have similar lengths. We prove that q-GRFs
give lower variance estimators of the $2$-regularised Laplacian kernel.

###  3.1 Theoretical results

In this section, we state and discuss our central theoretical results for the
q-GRFs mechanism. Sec. 3.1.1 provides a sketch, but full proofs are deferred
to App. 8.3. We remind the reader that results for
$(\mathbf{I}_{N}-\mathbf{U})^{-2}$ are trivially applied to
$\mathbf{K}_{\textrm{lap}}^{(2)}$ (see App. 8.1).

######  Theorem 3.2 (Antithetic termination is better than i.i.d.).

For any graph, q-GRFs will give lower variance on estimators of
$\left(\mathbf{I}_{N}-\mathbf{U}\right)^{-2}$ than regular GRFs provided
either i) the termination probability $p$ is sufficiently small or ii) the
spectral radius $\rho(\mathbf{U})$ is sufficiently small.

By ‘sufficiently small’ we mean that for a fixed $\rho(\mathbf{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 $\rho(\mathbf{U})$.
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, it is easy to see that the termination
probabilities in Eq. 13 are not particular to a TRV offset equal to
$\frac{1}{2}$, and in fact hold for any offset $\Delta$ satisfying
$p\leq\Delta\leq 1-p$. An immediate corollary is as follows.

######  Corollary 3.3 (Maximum size of an antithetic ensemble).

For a termination probability $p$, up to $\lfloor p^{-1}\rfloor$ random
walkers can all be conditioned to exhibit mutually antithetic termination.

This is achieved by offsetting their respective TRVs by $\Delta=p$. The
resulting _antithetic ensemble_ will have lower kernel estimator variance than
the equivalent number of i.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)<\Delta<p$ will exhibit
lower variance on estimators of $\left(\mathbf{I}_{N}-\mathbf{U}\right)^{-2}$
than independent walkers, provided either i) the termination probability $p$
is sufficiently small or ii) the spectral radius of the weighted adjacency
matrix $\rho(\mathbf{U})$ is sufficiently small.

This provides an upper limit of $\lfloor(p(1-p))^{-1}\rfloor$ on the number of
walkers we can simultaneously correlate before we can no longer guarantee that
coupling in a further walker’s TRV will be better than sampling it
independently. Intuitively, Theorem 3.4 tells us that we can space TRVs even
more closely than $p$, allowing us to increase the number of simultaneously
anticorrelated random walkers at 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.3.

From its Taylor expansion, the $(ij)$-th element of
$\left(\mathbf{I}_{N}-\mathbf{U}\right)^{-2}$ is nothing other than a sum over
all possible paths between the nodes $i$ and $j$, weighted by their lengths
and the respective products of edge weights. GRFs use a 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 adding
contributions wherever they intersect. In order to sample the space of walks
between nodes $i$ and $j$ more efficiently, it follows that we should make our
ensemble of walks out of each node more diverse. q-GRFs 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 $\phi(i)^{\top}\phi(j)$. Each
term in the resulting sum will take the form

$\small\begin{multlined}\sum_{x,y\in\mathcal{V}}\sum_{\omega_{1}\in\Omega_{ix}}\sum_{\omega_{2}\in\Omega_{jx}}\sum_{\omega_{3}\in\Omega_{iy}}\sum_{\omega_{4}\in\Omega_{jy}}\frac{\widetilde{\omega}(\omega_{1})}{p(\omega_{1})}\frac{\widetilde{\omega}(\omega_{2})}{p(\omega_{2})}\frac{\widetilde{\omega}(\omega_{3})}{p(\omega_{3})}\frac{\widetilde{\omega}(\omega_{4})}{p(\omega_{4})}\\\
\cdot
p(\omega_{1}\in\bar{\Omega}(k_{1},i),\omega_{3}\in\bar{\Omega}(k_{2},i))p(\omega_{2}\in\bar{\Omega}(l_{1},j),\omega_{4}\in\bar{\Omega}(l_{2},j)),\end{multlined}\sum_{x,y\in\mathcal{V}}\sum_{\omega_{1}\in\Omega_{ix}}\sum_{\omega_{2}\in\Omega_{jx}}\sum_{\omega_{3}\in\Omega_{iy}}\sum_{\omega_{4}\in\Omega_{jy}}\frac{\widetilde{\omega}(\omega_{1})}{p(\omega_{1})}\frac{\widetilde{\omega}(\omega_{2})}{p(\omega_{2})}\frac{\widetilde{\omega}(\omega_{3})}{p(\omega_{3})}\frac{\widetilde{\omega}(\omega_{4})}{p(\omega_{4})}\\\
\cdot
p(\omega_{1}\in\bar{\Omega}(k_{1},i),\omega_{3}\in\bar{\Omega}(k_{2},i))p(\omega_{2}\in\bar{\Omega}(l_{1},j),\omega_{4}\in\bar{\Omega}(l_{2},j)),\normalsize$
(16)

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
a function of the _length_ of each of the walks $\omega_{1,2,3,4}$. It is then
natural to write the sum over all walks between nodes $i$ and $x$ as a sum
over walk lengths $m$, with each term weighted by a combinatorial factor that
counts the number of walks of said length. This factor is
$(\mathbf{A}^{m})_{ix}$, with $\mathbf{A}$ the unweighted adjacency matrix.
That is,

$\sum_{\omega_{1}\in\Omega_{ix}}\left(\cdot\right)=\sum_{m=1}^{\infty}(\mathbf{A}^{m})_{ix}\left(\cdot\right).$
(17)

$(\mathbf{A}^{m})_{ix}$ is readily written as its eigendecomposition,
$\sum_{p=1}^{N}\lambda_{p}^{m}k_{pi}k_{px}$, with $\lambda_{p}$ the $p$th
eigenvalue and $k_{pi}$ the $i$-th coordinate of the $p$-th eigenvector
$\bm{k}_{p}$. We put these into Eq. 16 and perform the sums over each path
length $m_{1,2,3,4}$ from $1$ to $\infty$, arriving at

$\sum_{x,y\in\mathcal{V}}\sum_{k_{1},k_{2},k_{3},k_{4}}f(\lambda_{1},\lambda_{2},\lambda_{3},\lambda_{4})k_{1i}k_{1x}k_{2j}k_{2x}k_{3i}k_{3y}k_{4j}k_{4y}$
(18)

where $f(\lambda_{1},\lambda_{2},\lambda_{3},\lambda_{4})$ is a function of
the eigenvalues of $\mathbf{A}$ that depends on whether we correlate the
terminations of the walkers. Since the eigenvectors are orthogonal we have
that $\sum_{x\in\mathcal{V}}k_{1x}k_{2x}=\delta_{12}$, so this reduces to

$\sum_{k_{1},k_{3}}f(\lambda_{1},\lambda_{1},\lambda_{3},\lambda_{3})k_{1i}k_{1j}k_{3i}k_{4j}.$
(19)

Our task becomes to prove that this expression becomes smaller when we induce
antithetic termination. We achieve this by showing that a particular matrix is
negative definite.

Figure 2: Relative Frobenius norm error of estimator of the $2$-regularised
Laplacian kernel with GRFs (red circle) and q-GRFs (green cross). Lower is
better. We also include q-RRW-GRFs (Choromanski, 2023) (blue triangle),
instantiated with $f=\exp$, as a benchmark. Our novel q-GRFs perform the best
on every graph considered. $N$ is the number of nodes and, for the Erdős-Rényi
(ER) graphs, $p$ is the edge-generation probability. One standard deviation is
shaded but it is too small to easily see.

##  4 Experiments

In this section we report on empirical evaluations of q-GRFs. 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 as
described in 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,

$\mathbf{K}_{\textrm{lap}}^{(2)}=(\mathbf{I}_{N}+\sigma^{2}\widetilde{\mathbf{L}})^{-2},$
(20)

where $\widetilde{\mathbf{L}}\in\mathbb{R}^{N\times N}$ is the symmetrically
normalised Laplacian. $0<\sigma<1$ is a regulariser. We use both GRFs and
q-GRFs to generate unbiased estimates
$\widetilde{\mathbf{K}}_{\textrm{lap}}^{(2)}$ (see App. 8.1), then compute the
relative Frobenius norm
$\|\mathbf{K}_{\textrm{lap}}^{(2)}-\widetilde{\mathbf{K}}_{\textrm{lap}}^{(2)}\|_{\textrm{F}}^{2}/\|\mathbf{K}_{\textrm{lap}}^{(2)}\|^{2}_{\textrm{F}}$
between the true and approximated kernel matrices. This enables us to compare
the quality of the estimators. As a benchmark, we also include an
implementation of the high-level reinforced random walk QMC mechanism
suggested (but not tested) by Choromanski (2023). We choose the exponential
mapping as the reinforcement function $f$ (used to downweight the probability
of traversing previously-visited edges), although the optimal choice remains
an open problem. We refer to this mechanism as _q-RRW-GRFs_ to disambiguate
from our instantiation of q-GRFs (which uses antithetic termination).

Fig. 2 presents the results for a broad class of graphs: small Erdős-Rényi,
larger Erdős-Rényi, a binary rooted tree, a ladder, and four real-world
examples available 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 error. We use the regulariser $\sigma=0.1$
and the termination probability $p=0.5$.

The quality of kernel approximation naturally improves with the number of
walkers. Inducing antithetic coupling consistently reduces estimator variance,
with our q-GRF mechanism outperforming regular GRFs in every case. The exact
size of the gain depends on the particular graph (according to the closed
forms derived in App. 8.3), but improvement is always present. It is
intriguing that tree-like and planar graphs tend to enjoy a bigger gap; we
defer a rigorous theoretical analysis to future work. Meanwhile, the q-RRW-GRF
variant with the exponential $f$ is often worse than the regular mechanism and
is substantially more expensive. We do not include it in later experiments.

###  4.2 Scalable and accurate simulation of graph diffusion

In Sec. 2.1, we noted that the Laplacian $\mathbf{L}$ (or
$\widetilde{\mathbf{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-GRFs to provide a lower-variance low-rank decomposition of
$\mathbf{K}_{\textrm{lap}}^{(2)}$, we can simulate graph diffusion in a
scalable and accurate way.

Choosing a finite (even) number of discretisation timesteps $N_{t}$, we can
approximate the final state $\mathbf{u}_{t}$

$\mathbf{u}_{t}\simeq\widetilde{\mathbf{u}}_{t}\coloneqq\left[(\mathbf{I}_{N}+\frac{t}{N_{t}}\widetilde{\mathbf{L}})^{-2}\right]^{\frac{N_{t}}{2}}\mathbf{u}_{0}=\mathbf{K}_{\textrm{lap}}^{(2)\frac{N_{t}}{2}}\mathbf{u}_{0}.$
(21)

We can efficiently compute this using our low-rank GRF or q-GRF decomposition
of $\mathbf{K}_{\textrm{lap}}^{(2)}$ and compare the accuracy of
reconstruction of $\mathbf{u}_{t}$. In particular, we take an initial one-hot
state $\mathbf{u}_{0}=\left(1,0,0,...,0\right)^{\top}$ and simulate diffusion
for $t=1$ divided into $N_{t}=1000$ timesteps, using $m=10$ walkers and a
termination probability $p=0.5$. Fig. 3 gives a schematic. We average the MSE
of $\widetilde{\mathbf{u}}_{t}$ over $1000$ trials. Table 1 reports the
results; q-GRFs approximate the time evolution operator more accurately so
consistently give a lower simulation error. Since the time-evolution operator
is applied repeatedly, even modest improvements in its approximation can lead
to a substantially better reconstruction of the final state. In extreme cases
the simulation error is halved.

Table 1:  MSE of the approximation of the final state $\mathbf{u}_{t}$ by repeated application of a low-rank decomposition of the discrete time-evolution operator, constructed with (q-)GRFs. Brackets give one standard deviation. Graph | $N$ |  Sim error, $\mathrm{MSE}(\widetilde{\mathbf{u}}_{t})$
---|---|---
|  | 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)
$\mathrm{karate}$ | 34 | 0.066(2) | 0.054(1)
$\mathrm{dolphins}$ | 62 | 0.0165(4) | 0.0139(3)
$\mathrm{football}$ | 115 | 0.0170(3) | 0.0160(2)
$\mathrm{eurosis}$ | 1272 | 0.089(1) | 0.084(1)
Figure 3: Schematic of diffusion on a small graph, with heat initially
localised on a single node spreading under the action of the Laplace operator.
We use (q-)GRFs to estimate the quantity on the right.$t=0$$t=1$

###  4.3 Kernelised $k$-means clustering for graph nodes

Table 2: Errors in kernelised $k$-means clustering when approximating the Gram matrix with q-GRFs and GRFs. Lower is better. Graph | $N$ |  Clustering error, $E_{c}$
---|---|---
|  | GRFs | q-GRFs
$\mathrm{karate}$ | 34 | 0.11 | 0.05
$\mathrm{databases}$ | 1046 | 0.17 | 0.11
$\mathrm{polbooks}$ | 105 | 0.28 | 0.28
$\mathrm{dolphins}$ | 62 | 0.40 | 0.38
$\mathrm{football}$ | 115 | 0.10 | 0.09
$\mathrm{citeseer}$ | 3300 | 0.02 | 0.01
$\mathrm{eurosis}$ | 1285 | 0.15 | 0.19
$\mathrm{eu}$-$\mathrm{core}$ | 1005 | 0.04 | 0.05

Next we evaluate the performance of GRFs and q-GRFs on the task of assigning
nodes to clusters using the graph kernel, as described by Dhillon et al.
(2004). We first run the algorithm to obtain $N_{c}=2$ clusters with the exact
$2$-regularised Laplacian kernel, then compare the results when we use its
approximation via GRFs and q-GRFs. In each case, we report the _clustering
error_ , defined by

$E_{c}\coloneqq\frac{\textrm{no. wrong pairs}}{N(N-1)/2}.$ (22)

This is simply the number of misclassified pairs (in the sense of being
assigned to the same cluster when the converse is true or vice versa) divided
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 eight cases
q-GRFs provide equally good or better results.

###  4.4 Kernel regression for node attribute prediction

Table 3: $1-\cos\theta$ between true and predicted node vectors when approximating the Gram matrix with q-GRFs and GRFs. Lower is better. Brackets give one standard deviation. Graph | $N$ |  Pred error, $1-\cos\theta$
---|---|---
|  | GRFs | q-GRFs
$\mathrm{cylinder}$ | 210 | 0.104(1) | 0.101(1)
$\mathrm{teapot}$ | 480 | 0.0531(5) | 0.0493(5)
$\mathrm{idler}$-$\mathrm{riser}$ | 782 | 0.0881(6) | 0.0852(5)
$\mathrm{busted}$ | 1941 | 0.00690(4) | 0.00661(4)
$\mathrm{torus}$ | 4350 | 0.00131(1) |  0.00120(1)

Lastly, we consider the problem of kernel regression on a triangular mesh
graph. Each node is associated with a normal vector $\bm{v}^{(i)}$ (equal to
the mean of the normal vectors of its surrounding faces). We consider five
meshes of different sizes, available from (Dawson-Haggerty, 2023). We predict
a random $5\%$ split of the vectors $\bm{v}^{(i)}$ (‘test’) from the remaining
$95\%$ (‘train’) using

$\widetilde{\bm{v}}^{(i)}\coloneqq\frac{\sum_{j}\mathbf{K}_{\textrm{lap}}^{(2)}(i,j)\bm{v}^{(j)}}{\sum_{j}\mathbf{K}_{\textrm{lap}}^{(2)}(i,j)},$
(23)

where $j$ sums over the training vertices. This is simply a linear combination
of all the training node normal vectors, weighted by their respective kernel
evaluations. We compute the average angular error $1-\cos\theta$ between the
prediction $\widetilde{\bm{v}}^{(i)}$ and groundtruth ${\bm{v}}^{(i)}$ across
the test set, comparing the result when $\mathbf{K}_{\textrm{lap}}^{(2)}$ is
approximated with GRFs and q-GRFs with $m=6$ random walks at a termination
probability $p=0.5$. The regulariser is $\sigma=0.1$. q-GRFs enjoy lower
estimator variance so consistently give better predictions of the missing
vectors.

##  5 Conclusion

We have proposed a novel class of quasi-Monte Carlo graph random features
(q-GRFs) for unbiased and efficient estimation of kernels defined on the nodes
of a 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 in 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.

##  6 Broader impacts and limitations

We envisage several possible impacts of our novel algorithm. First, graphs
provide a 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
a procedure to improve the sampling efficiency of random walks, so it 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 research in this new
domain. Our work is foundational with no immediate direct negative societal
impacts that we can see. However, it is 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 natural directions for future research.
First, although we have derived closed-form expressions for the kernel
estimator variance with GRFs and q-GRFs in App. 8.3, 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.
A 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.

##  7 Relative contributions and acknowledgements

IR devised the antithetic termination scheme, proved all theoretical results
and ran the experiments in Secs 4.1, 4.2 and 4.4. KC provided crucial support
throughout, particularly: running the clustering experiment in 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
4.2 and 4.4. AW gave important guidance and feedback on the manuscript.

IR acknowledges support from a Trinity College External Studentship. AW
acknowledges support from a 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.

## References

  * Albert and Barabási (2002) Réka Albert and Albert-László Barabási.  Statistical mechanics of complex networks.  _Reviews of modern physics_ , 74(1):47, 2002.  URL https://doi.org/10.1103/RevModPhys.74.47.
  * Avron et al. (2017) Haim Avron, Michael Kapralov, Cameron Musco, Christopher Musco, Ameya Velingker, and Amir Zandieh.  Random fourier features for kernel ridge regression: Approximation bounds and statistical guarantees.  In Doina Precup and Yee Whye Teh, editors, _Proceedings of the 34th International Conference on Machine Learning, ICML 2017, Sydney, NSW, Australia, 6-11 August 2017_ , volume 70 of _Proceedings of Machine Learning Research_ , pages 253–262. PMLR, 2017.  URL http://proceedings.mlr.press/v70/avron17a.html.
  * Belkin and Niyogi (2003) Mikhail Belkin and Partha Niyogi.  Laplacian eigenmaps for dimensionality reduction and data representation.  _Neural computation_ , 15(6):1373–1396, 2003\.  URL https://dl.acm.org/doi/10.1162/089976603321780317.
  * Bollacker et al. (1998) 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.  URL https://dl.acm.org/doi/10.1145/280765.280786.
  * Campbell (2002) Colin Campbell.  Kernel methods: a survey of current techniques.  _Neurocomputing_ , 48(1-4):63–84, 2002.  doi: 10.1016/S0925-2312(01)00643-9.  URL https://doi.org/10.1016/S0925-2312(01)00643-9.
  * Canu and Smola (2006) Stéphane Canu and Alexander J. Smola.  Kernel methods and the exponential family.  _Neurocomputing_ , 69(7-9):714–720, 2006.  doi: 10.1016/j.neucom.2005.12.009.  URL https://doi.org/10.1016/j.neucom.2005.12.009.
  * Chen et al. (2016) Xiaowei Chen, Yongkun Li, Pinghui Wang, and John Lui.  A general framework for estimating graphlet statistics via random walk.  _arXiv preprint arXiv:1603.07504_ , 2016.  URL https://doi.org/10.48550/arXiv.1603.07504.
  * Cho and Saul (2011) Youngmin Cho and Lawrence K. Saul.  Analysis and extension of arc-cosine kernels for large margin classification.  _CoRR_ , abs/1112.3712, 2011.  URL https://doi.org/10.48550/arXiv.1112.3712.
  * Choromanski (2023) Krzysztof Choromanski.  Taming graph kernels with random features, 2023.  URL https://doi.org/10.48550/arXiv.2305.00156.
  * Choromanski et al. (2020) 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.14794_ , 2020.  URL https://doi.org/10.48550/arXiv.2009.14794.
  * Choromanski et al. (2017) 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 https://doi.org/10.48550/arXiv.1703.00864.
  * Chung and Yau (1999) Fan R. K. Chung and Shing-Tung Yau.  Coverings, heat kernels and spanning trees.  _Electron. J. Comb._ , 6, 1999.  doi: 10.37236/1444.  URL https://doi.org/10.37236/1444.
  * Chung (1997) Fan RK Chung.  _Spectral graph theory_ , volume 92.  American Mathematical Soc., 1997.
  * Dasgupta et al. (2010) Anirban Dasgupta, Ravi Kumar, and Tamás Sarlós.  A sparse johnson: Lindenstrauss transform.  In Leonard J. Schulman, editor, _Proceedings of the 42nd ACM Symposium on Theory of Computing, STOC 2010, Cambridge, Massachusetts, USA, 5-8 June 2010_ , pages 341–350. ACM, 2010.  doi: 10.1145/1806689.1806737.  URL https://doi.org/10.1145/1806689.1806737.
  * Dawson-Haggerty (2023) Michael Dawson-Haggerty.  Trimesh repository, 2023.  URL https://github.com/mikedh/trimesh.
  * Dhillon et al. (2004) Inderjit S Dhillon, Yuqiang Guan, and Brian Kulis.  Kernel k-means: spectral clustering and normalized cuts.  In _Proceedings of the tenth ACM SIGKDD international conference on Knowledge discovery and data mining_ , pages 551–556, 2004.  URL https://dl.acm.org/doi/10.1145/1014052.1014118.
  * Dick et al. (2013) Josef Dick, Frances Y Kuo, and Ian H Sloan.  High-dimensional integration: the quasi-monte carlo way.  _Acta Numerica_ , 22:133–288, 2013.
  * Fang et al. (2021) 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 https://proceedings.neurips.cc/paper_files/paper/2021/file/93da579a65ce84cd1d4c85c2cbb84fc5-Paper.pdf.
  * Goemans and Williamson (2004) Michel X. Goemans and David P. Williamson.  Approximation algorithms for m${}_{\mbox{ax}}$-3-c${}_{\mbox{ut}}$ and other problems via complex semidefinite programming.  _J. Comput. Syst. Sci._ , 68(2):442–470, 2004\.  doi: 10.1016/j.jcss.2003.07.012.  URL https://doi.org/10.1016/j.jcss.2003.07.012.
  * Hammersley and Morton (1956) JM Hammersley and KW Morton.  A new monte carlo technique: antithetic variates.  In _Mathematical proceedings of the Cambridge philosophical society_ , volume 52, pages 449–475. Cambridge University Press, 1956.
  * Ingraham et al. (2019) 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://www.mit.edu/~vgarg/GenerativeModelsForProteinDesign.pdf.
  * Ivashkin (2023) Vladimir Ivashkin.  Community graphs repository, 2023.  URL https://github.com/vlivashkin/community-graphs.
  * Johnson (1984) William B Johnson.  Extensions of lipschitz mappings into a hilbert space.  _Contemp. Math._ , 26:189–206, 1984.
  * Kondor and Lafferty (2002) Risi Kondor and John D. Lafferty.  Diffusion kernels on graphs and other discrete input spaces.  In Claude Sammut and Achim G. Hoffmann, editors, _Machine Learning, Proceedings of the Nineteenth International Conference (ICML 2002), University of New South Wales, Sydney, Australia, July 8-12, 2002_ , pages 315–322. Morgan Kaufmann, 2002.  URL https://www.ml.cmu.edu/research/dap-papers/kondor-diffusion-kernels.pdf.
  * Kontorovich et al. (2008) Leonid Kontorovich, Corinna Cortes, and Mehryar Mohri.  Kernel methods for learning languages.  _Theor. Comput. Sci._ , 405(3):223–236, 2008\.  doi: 10.1016/j.tcs.2008.06.037.  URL https://doi.org/10.1016/j.tcs.2008.06.037.
  * Leskovec et al. (2007) Jure Leskovec, Jon Kleinberg, and Christos Faloutsos.  Graph evolution: Densification and shrinking diameters.  _ACM transactions on Knowledge Discovery from Data (TKDD)_ , 1(1):2–es, 2007.  URL https://dl.acm.org/doi/10.1145/1217299.1217301.
  * Liu et al. (2022) 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):7128–7148, 2022.  doi: 10.1109/TPAMI.2021.3097011.  URL https://doi.org/10.1109/TPAMI.2021.3097011.
  * Lusseau et al. (2003) David Lusseau, Karsten Schneider, Oliver J 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 Sociobiology_ , 54:396–405, 2003\.
  * Lyu (2017) Yueming Lyu.  Spherical structured feature maps for kernel approximation.  In _International Conference on Machine Learning_ , pages 2256–2264. PMLR, 2017.  URL http://proceedings.mlr.press/v70/lyu17a/lyu17a.pdf.
  * Newman (2006) Mark EJ Newman.  Modularity and community structure in networks.  _Proceedings of the national academy of sciences_ , 103(23):8577–8582, 2006.  URL https://doi.org/10.1073/pnas.0601602103.
  * Nikolentzos and Vazirgiannis (2020) Giannis Nikolentzos and Michalis Vazirgiannis.  Random walk graph neural networks.  _Advances in Neural Information Processing Systems_ , 33:16211–16222, 2020.  URL https://proceedings.neurips.cc/paper/2020/file/ba95d78a7c942571185308775a97a3a0-Paper.pdf.
  * Rahimi and Recht (2007) Ali Rahimi and Benjamin Recht.  Random features for large-scale kernel machines.  _Advances in neural information processing systems_ , 20, 2007.  URL https://people.eecs.berkeley.edu/~brecht/papers/07.rah.rec.nips.pdf.
  * Rahimi and Recht (2008) Ali Rahimi and Benjamin Recht.  Weighted sums of random kitchen sinks: Replacing minimization with randomization in learning.  In Daphne Koller, Dale Schuurmans, Yoshua Bengio, and Léon Bottou, editors, _Advances in Neural Information Processing Systems 21, Proceedings of the Twenty-Second Annual Conference on Neural Information Processing Systems, Vancouver, British Columbia, Canada, December 8-11, 2008_ , pages 1313–1320. Curran Associates, Inc., 2008.  URL https://papers.nips.cc/paper_files/paper/2008/file/0efe32849d230d7f53049ddc4a4b0c60-Paper.pdf.
  * Reid et al. (2023) Isaac Reid, Krzysztof Choromanski, Valerii Likhosherstov, and Adrian Weller.  Simplex random features.  _arXiv preprint arXiv:2301.13856_ , 2023.  URL https://doi.org/10.48550/arXiv.2301.13856.
  * Ribeiro et al. (2021) Pedro Ribeiro, Pedro Paredes, Miguel EP Silva, David Aparicio, and Fernando Silva.  A survey on subgraph counting: concepts, algorithms, and applications to network motifs and graphlets.  _ACM Computing Surveys (CSUR)_ , 54(2):1–36, 2021\.  URL https://doi.org/10.1145/3433652.
  * Roweis and Saul (2000) Sam T Roweis and Lawrence K Saul.  Nonlinear dimensionality reduction by locally linear embedding.  _science_ , 290(5500):2323–2326, 2000.  URL https://www.science.org/doi/10.1126/science.290.5500.2323.
  * Shervashidze and Borgwardt (2009) Nino Shervashidze and Karsten M. Borgwardt.  Fast subtree kernels on graphs.  In Yoshua Bengio, Dale Schuurmans, John D. Lafferty, Christopher K. I. Williams, and Aron Culotta, editors, _Advances in Neural Information Processing Systems 22: 23rd Annual Conference on Neural Information Processing Systems 2009. Proceedings of a meeting held 7-10 December 2009, Vancouver, British Columbia, Canada_ , pages 1660–1668. Curran Associates, Inc., 2009.  URL https://proceedings.neurips.cc/paper/2009/hash/0a49e3c3a03ebde64f85c0bacd8a08e2-Abstract.html.
  * Shervashidze et al. (2009) Nino Shervashidze, S. V. N. Vishwanathan, Tobias Petri, Kurt Mehlhorn, and Karsten M. Borgwardt.  Efficient graphlet kernels for large graph comparison.  In David A. Van Dyk and Max Welling, editors, _Proceedings of the Twelfth International Conference on Artificial Intelligence and Statistics, AISTATS 2009, Clearwater Beach, Florida, USA, April 16-18, 2009_ , volume 5 of _JMLR Proceedings_ , pages 488–495. JMLR.org, 2009.  URL http://proceedings.mlr.press/v5/shervashidze09a.html.
  * Smola and Kondor (2003) Alexander J. Smola and Risi Kondor.  Kernels and regularization on graphs.  In Bernhard Schölkopf and Manfred K. Warmuth, editors, _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.  doi: 10.1007/978-3-540-45167-9\\_12.  URL https://doi.org/10.1007/978-3-540-45167-9_12.
  * Smola and Schölkopf (2002) Alexander J. Smola and Bernhard Schölkopf.  Bayesian kernel methods.  In Shahar Mendelson and Alexander J. Smola, editors, _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. Springer, 2002.  doi: 10.1007/3-540-36434-X\\_3.  URL https://doi.org/10.1007/3-540-36434-X_3.
  * Strubell et al. (2019) Emma Strubell, Ananya Ganesh, and Andrew McCallum.  Energy and policy considerations for deep learning in nlp.  _arXiv preprint arXiv:1906.02243_ , 2019.  URL https://doi.org/10.18653/v1/P19-1355.
  * Vaswani et al. (2017) Ashish Vaswani, Noam Shazeer, Niki Parmar, Jakob Uszkoreit, Llion Jones, Aidan N. Gomez, Lukasz Kaiser, and Illia Polosukhin.  Attention is all you need.  In Isabelle Guyon, Ulrike von Luxburg, Samy Bengio, Hanna M. Wallach, Rob Fergus, S. V. N. Vishwanathan, and Roman Garnett, editors, _Advances in Neural Information Processing Systems 30: Annual Conference on Neural Information Processing Systems 2017, December 4-9, 2017, Long Beach, CA, USA_ , pages 5998–6008, 2017.  URL https://proceedings.neurips.cc/paper/2017/hash/3f5ee243547dee91fbd053c1c4a845aa-Abstract.html.
  * Vishwanathan et al. (2006) S. V. N. Vishwanathan, Karsten M. Borgwardt, and Nicol N. Schraudolph.  Fast computation of graph kernels.  In Bernhard Schölkopf, John C. Platt, and Thomas Hofmann, editors, _Advances in Neural Information Processing Systems 19, Proceedings of the Twentieth Annual Conference on Neural Information Processing Systems, Vancouver, British Columbia, Canada, December 4-7, 2006_ , pages 1449–1456. MIT Press, 2006.  URL https://proceedings.neurips.cc/paper/2006/hash/e37b08dd3015330dcbb5d6663667b8b8-Abstract.html.
  * Williams and Rasmussen (2006) Christopher KI Williams and Carl Edward Rasmussen.  _Gaussian processes for machine learning_ , volume 2.  MIT press Cambridge, MA, 2006.
  * Wu et al. (2019) Lingfei Wu, Ian En-Hsu Yen, Zhen Zhang, Kun Xu, Liang Zhao, Xi Peng, Yinglong Xia, and Charu Aggarwal.  Scalable global alignment graph kernel using random features: From node embedding to graph embedding.  In _Proceedings of the 25th ACM SIGKDD International Conference on Knowledge Discovery & Data Mining_, pages 1418–1428, 2019.  URL https://doi.org/10.48550/arXiv.1911.11119.
  * Yang et al. (2014) Jiyan Yang, Vikas Sindhwani, Haim Avron, and Michael W. Mahoney.  Quasi-monte carlo feature maps for shift-invariant 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_ , pages 485–493. JMLR.org, 2014\.  URL http://proceedings.mlr.press/v32/yangb14.html.
  * Yu et al. (2016) 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 https://doi.org/10.48550/arXiv.1610.09072.
  * Zachary (1977) Wayne W Zachary.  An information flow model for conflict and fission in small groups.  _Journal of anthropological research_ , 33(4):452–473, 1977.  URL https://www.journals.uchicago.edu/doi/10.1086/jar.33.4.3629752.

##  8 Appendix: quasi-Monte Carlo graph random features (q-GRFs)

###  8.1 On the approximation of the $d$-regularised Laplacian using GRFs

In this appendix, we demonstrate how to approximate the $d$-regularised
Laplacian $\mathbf{K}_{\textrm{lap}}^{(d)}$ with GRFs.

Recall that GRFs provide an estimator to the quantity
$\left(\mathbf{I}_{N}-\mathbf{U}\right)^{-2}$ where $\mathbf{U}$ is a weighted
adjacency matrix. Recall also that the matrix elements of the symmetrically
normalised Laplacian $\widetilde{\mathbf{L}}$ are given by

$\widetilde{\mathbf{L}}_{ij}=\begin{cases}1&$if $i=j,\\\
-\frac{\mathbf{W_{ij}}}{\sqrt{\textrm{deg}_{\mathbf{W}}(i)\textrm{deg}_{\mathbf{W}}(j)}}&$if
$i\sim j.\end{cases}$ (24)

where $\textrm{deg}_{\mathbf{W}}(i)=\sum_{j\in\mathcal{V}}\mathbf{W}_{ij}$ is
the weighted degree of the node $i$. We are typically interested in situations
where $\mathbf{W}=\mathbf{A}$, an unweighted adjacency matrix. Now note that

$\mathbf{K}_{\textrm{lap
}ij}^{(2)}=(\mathbf{I}_{N}+\sigma^{2}\widetilde{\mathbf{L}})^{-2}_{ij}=(1+\sigma^{2})^{-2}\left(\mathbf{I}_{N}-\mathbf{U}\right)^{-2}_{ij}$
(25)

where we defined the matrix $\mathbf{U}$ with matrix elements

$\mathbf{U}_{ij}=\frac{\sigma^{2}}{1+\sigma^{2}}\frac{\mathbf{W_{ij}}}{\sqrt{\textrm{deg}_{\mathbf{W}}(i)\textrm{deg}_{\mathbf{W}}(j)}}.$
(26)

This is itself a weighted adjacency matrix, as required. It follows that, by
estimating $\left(\mathbf{I}_{N}-\mathbf{U}\right)^{-2}$ with GRFs, we can
trivially estimate $\mathbf{K}_{\textrm{lap}}^{(2)}$. This was reported in
[Choromanski, 2023].

Supposing that we have constructed a low-rank GRF estimator

$\mathbf{K}_{\textrm{lap}}^{(2)}=\mathbb{E}\left[\mathbf{C}\mathbf{C}^{\top}\right]$
(27)

where the matrix $\mathbf{C}\in\mathbb{R}^{N\times N}$ has rows
$\mathbf{C}_{i}\coloneqq\frac{1}{1+\sigma^{2}}\phi(i)^{\top}$, we note that it
is straightforward to construct the $1$-regularised Laplacian kernel estimator

$\mathbf{K}_{\textrm{lap}}^{(1)}=\mathbb{E}\left[\mathbf{C}\mathbf{D}^{\top}\right]$
(28)

by taking
$\mathbf{D}\coloneqq\left(\mathbf{I}_{N}+\sigma^{2}\widetilde{\mathbf{L}}\right)^{\top}\mathbf{C}$.
It is then trivial to obtain the estimator $\mathbf{K}_{\textrm{lap}}^{(d)}$
for arbitrary $d\in\mathbb{N}$.

###  8.2 Derivation of Eq. 15

In this appendix we derive Eq. 15, which gives the expected length of some
walk $\omega_{2}$ given that its antithetic partner $\omega_{1}$ is of length
$m$: that is,
$\mathbb{E}\left(\textrm{len}(\omega_{2})|\textrm{len}(\omega_{1})=m\right)$.

As a warm-up, consider the simpler marginal expected lengths. Note that

$p\left(\textrm{len}(\omega)=m\right)=(1-p)^{m}p.$ (29)

It follows that

$\mathbb{E}\left(\textrm{len}(\omega)\right)=\sum_{m=0}^{\infty}m(1-p)^{m}p=\frac{1-p}{p}$
(30)

where we computed the arithmetic-geometric series. We reported this result in
Eq. 14. Meanwhile, the probability of a walk being of length $i$ given that
its antithetic partner is of length $m$ is

$p(\textrm{len}(\omega_{2})=i|\textrm{len}(\omega_{1})=m)=\begin{cases}\left(\frac{1-2p}{1-p}\right)^{i}\frac{p}{1-p}&\textrm{if
}i<m,\\\ 0&\textrm{if }i=m,\\\
\left(\frac{1-2p}{1-p}\right)^{m}\left(1-p\right)^{i-m-1}p&\textrm{if
}i>m.\end{cases}$ (31)

The analagous sum then becomes

$\mathbb{E}\left(\textrm{len}(\omega_{2})|\textrm{len}(\omega_{1})=m\right)=\sum_{i=0}^{m}\left(\frac{1-2p}{1-p}\right)^{i}\frac{p}{1-p}i+\sum_{i=m+1}^{\infty}\left(\frac{1-2p}{1-p}\right)^{m}\left(1-p\right)^{i-m-1}pi.$
(32)

After straightforward but tedious algebra, this evaluates to

$\mathbb{E}\left(\textrm{len}(\omega_{2})|\textrm{len}(\omega_{1})=m\right)=\frac{1-2p}{p}+2\left(\frac{1-2p}{1-p}\right)^{m},$
(33)

as stated in Eq. 15. Note that this is greater than
$\mathbb{E}(\textrm{len}(\omega))$ when $m$ is small and smaller than
$\mathbb{E}(\textrm{len}(\omega))$ when $m$ is large; the two walk lengths are
negatively correlated.

###  8.3 On the superiority of q-GRFs (proof of Theorem 3.2)

Here, we provide a proof of the central result of Theorem 3.2: that the
introduction of antithetic termination reduces the variance of estimators of
the matrix $(\mathbf{I}_{N}-\mathbf{U})^{-2}$. From App. 8.1, all our results
will trivially extend to the $2$-regularised Laplacian kernel
$\mathbf{K}_{\textrm{lap}}^{(2)}$.

Notation: to reduce the burden of summation indices, we have used Dirac’s bra-
ket notation from quantum mechanics. $\ket{y}$ can be interpreted as the
vector $\bm{y}$ and $\bra{y}$ as $\bm{y}^{\top}$.

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 of its neighbours with equal probability at each timestep.
We will relax these assumptions in App. 8.4.

We have seen that antithetic termination does not modify the walkers’ marginal
termination behaviour, so the variance of the estimator
$\bm{\phi}(i)^{\top}\bm{\phi}(j)$ is only affected via the second-order term
$\mathbb{E}\left[(\bm{\phi}(i)^{\top}\bm{\phi}(j))^{2}\right]$. Writing out
the sums,

$\small\begin{multlined}(\bm{\phi}(i)^{\top}\bm{\phi}(j))^{2}=\frac{1}{m^{4}}\sum_{x,y\in\mathcal{V}}\sum_{k_{1},l_{1},k_{2},l_{2}=1}^{m}\sum_{\omega_{1}\in\Omega_{ix}}\sum_{\omega_{2}\in\Omega_{jx}}\sum_{\omega_{3}\in\Omega_{iy}}\sum_{\omega_{4}\in\Omega_{jy}}\frac{\widetilde{\omega}(\omega_{1})}{p(\omega_{1})}\frac{\widetilde{\omega}(\omega_{2})}{p(\omega_{2})}\frac{\widetilde{\omega}(\omega_{3})}{p(\omega_{3})}\frac{\widetilde{\omega}(\omega_{4})}{p(\omega_{4})}\\\
\cdot\mathbb{I}(\omega_{1}\in\bar{\Omega}(k_{1},i))\mathbb{I}(\omega_{2}\in\bar{\Omega}(l_{1},j))\mathbb{I}(\omega_{3}\in\bar{\Omega}(k_{2},i))\mathbb{I}(\omega_{4}\in\bar{\Omega}(l_{2},j)).\end{multlined}(\bm{\phi}(i)^{\top}\bm{\phi}(j))^{2}=\frac{1}{m^{4}}\sum_{x,y\in\mathcal{V}}\sum_{k_{1},l_{1},k_{2},l_{2}=1}^{m}\sum_{\omega_{1}\in\Omega_{ix}}\sum_{\omega_{2}\in\Omega_{jx}}\sum_{\omega_{3}\in\Omega_{iy}}\sum_{\omega_{4}\in\Omega_{jy}}\frac{\widetilde{\omega}(\omega_{1})}{p(\omega_{1})}\frac{\widetilde{\omega}(\omega_{2})}{p(\omega_{2})}\frac{\widetilde{\omega}(\omega_{3})}{p(\omega_{3})}\frac{\widetilde{\omega}(\omega_{4})}{p(\omega_{4})}\\\
\cdot\mathbb{I}(\omega_{1}\in\bar{\Omega}(k_{1},i))\mathbb{I}(\omega_{2}\in\bar{\Omega}(l_{1},j))\mathbb{I}(\omega_{3}\in\bar{\Omega}(k_{2},i))\mathbb{I}(\omega_{4}\in\bar{\Omega}(l_{2},j)).\normalsize$
(34)

To remind the reader: the variables $x,y$ sum over the nodes of the graph
$\mathcal{V}$. $k_{1}$ and $l_{1}$ enumerate all the $m$ walks sampled out of
node $i$, whilst $k_{2}$ and $l_{2}$ enumerate walks from $j$. The sum over
$\omega_{1}\in\Omega_{ix}$ is over all possible walks between nodes $i$ and
$x$. $\widetilde{\omega}(\omega_{1})$ evaluates the product of edge weights
traversed by the walk $\omega_{1}$, which is $w^{\textrm{len}(\omega_{1})}$ in
the equal-weights case (with $\textrm{len}(\omega_{1})$ denoting the number of
edges in $\omega_{1}$). $p(\omega_{1})$ is the _marginal_ probability of the
subwalk $\omega_{1}$, which is equal to $((1-p)/d)^{\textrm{len}(\omega_{1})}$
on a $d$-regular graph. Lastly, the indicator function
$\mathbb{I}(\omega_{1}\in\bar{\Omega}(k_{1},i))$ evaluates to $1$ if the
$k_{1}$th walk out of node $i$ (denoted $\bar{\Omega}(k_{1},i)$) contains the
walk $\omega_{1}$ as a subwalk and $0$ otherwise.

We immediately note that our scheme only every correlates walks leaving the
same node, so walks out of different nodes remain independent. Therefore,

$\begin{multlined}\mathbb{E}\left[\mathbb{I}(\omega_{1}\in\bar{\Omega}(k_{1},i))\mathbb{I}(\omega_{2}\in\bar{\Omega}(l_{1},j))\mathbb{I}(\omega_{3}\in\bar{\Omega}(k_{2},i))\mathbb{I}(\omega_{4}\in\bar{\Omega}(l_{2},j))\right]\\\
=p(\omega_{1}\in\bar{\Omega}(k_{1},i),\omega_{3}\in\bar{\Omega}(k_{2},i))p(\omega_{2}\in\bar{\Omega}(l_{1},j),\omega_{4}\in\bar{\Omega}(l_{2},j)).\end{multlined}\mathbb{E}\left[\mathbb{I}(\omega_{1}\in\bar{\Omega}(k_{1},i))\mathbb{I}(\omega_{2}\in\bar{\Omega}(l_{1},j))\mathbb{I}(\omega_{3}\in\bar{\Omega}(k_{2},i))\mathbb{I}(\omega_{4}\in\bar{\Omega}(l_{2},j))\right]\\\
=p(\omega_{1}\in\bar{\Omega}(k_{1},i),\omega_{3}\in\bar{\Omega}(k_{2},i))p(\omega_{2}\in\bar{\Omega}(l_{1},j),\omega_{4}\in\bar{\Omega}(l_{2},j)).$
(35)

Consider the term in the sum corresponding to one particular set of walks
$(k_{1},l_{1},k_{2},l_{2})$,

$\small\begin{multlined}\sum_{x,y\in\mathcal{V}}\sum_{\omega_{1}\in\Omega_{ix}}\sum_{\omega_{2}\in\Omega_{jx}}\sum_{\omega_{3}\in\Omega_{iy}}\sum_{\omega_{4}\in\Omega_{jy}}\frac{\widetilde{\omega}(\omega_{1})}{p(\omega_{1})}\frac{\widetilde{\omega}(\omega_{2})}{p(\omega_{2})}\frac{\widetilde{\omega}(\omega_{3})}{p(\omega_{3})}\frac{\widetilde{\omega}(\omega_{4})}{p(\omega_{4})}\\\
\cdot
p(\omega_{1}\in\bar{\Omega}(k_{1},i),\omega_{3}\in\bar{\Omega}(k_{2},i))p(\omega_{2}\in\bar{\Omega}(l_{1},j),\omega_{4}\in\bar{\Omega}(l_{2},j)).\end{multlined}\sum_{x,y\in\mathcal{V}}\sum_{\omega_{1}\in\Omega_{ix}}\sum_{\omega_{2}\in\Omega_{jx}}\sum_{\omega_{3}\in\Omega_{iy}}\sum_{\omega_{4}\in\Omega_{jy}}\frac{\widetilde{\omega}(\omega_{1})}{p(\omega_{1})}\frac{\widetilde{\omega}(\omega_{2})}{p(\omega_{2})}\frac{\widetilde{\omega}(\omega_{3})}{p(\omega_{3})}\frac{\widetilde{\omega}(\omega_{4})}{p(\omega_{4})}\\\
\cdot
p(\omega_{1}\in\bar{\Omega}(k_{1},i),\omega_{3}\in\bar{\Omega}(k_{2},i))p(\omega_{2}\in\bar{\Omega}(l_{1},j),\omega_{4}\in\bar{\Omega}(l_{2},j)).\normalsize$
(36)

This object will be of central importance and is referred to as the
_correlation term_. In the sum over $k_{1},k_{2},l_{1},l_{2}$, there are three
possibilities to consider. We stress again that $k_{1,2}$ refers to a pair of
walks out of node $i$ and $l_{1,2}$ refers to a pair out of $j$.

  * •

Case 1, same-same, $\bm{k_{1}=k_{2}}$, $\bm{l_{1}=l_{2}}$: the pair of walks
out of $i$ 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, $\bm{k_{1}\neq k_{2}}$, $\bm{l_{1}\neq l_{2}}$:
the walks out of both $i$ and $j$ differ, and each pair may be antithetic or
independent. This term will be modified by the coupling.

  * •

Case 3, same-different. $\bm{k_{1}=k_{2}}$, $\bm{l_{1}\neq l_{2}}$: the walks
out of $i$ 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 $i$ 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. 34 we can conclude that the entire sum and therefore
the Laplacian kernel estimator variance is suppressed. For completeness, we
write out the entire sum from Eq. 34 with the degeneracy factors below:

$\small\begin{multlined}(\bm{\phi}(i)^{\top}\bm{\phi}(j))^{2}=\frac{1}{m^{4}}\sum_{x,y\in\mathcal{V}}\Bigl{\\{}\Bigr{.}\\\
\begin{rcases}m^{2}\sum_{\omega_{1}\in\Omega_{ix}}\sum_{\omega_{2}\in\Omega_{jx}}\sum_{\omega_{3}\in\Omega_{iy}}\sum_{\omega_{4}\in\Omega_{jy}}\frac{\widetilde{\omega}(\omega_{1})}{p(\omega_{1})}\frac{\widetilde{\omega}(\omega_{2})}{p(\omega_{2})}\frac{\widetilde{\omega}(\omega_{3})}{p(\omega_{3})}\frac{\widetilde{\omega}(\omega_{4})}{p(\omega_{4})}\\\
\cdot\mathbb{I}(\omega_{1}\in\bar{\Omega}(k_{1},i))\mathbb{I}(\omega_{2}\in\bar{\Omega}(l_{1},j))\mathbb{I}(\omega_{3}\in\bar{\Omega}(k_{1},i))\mathbb{I}(\omega_{4}\in\bar{\Omega}(l_{1},j))\end{rcases}\text{same-
same (1)}\\\
\begin{rcases}+m^{2}(m-1)^{2}\sum_{\omega_{1}\in\Omega_{ix}}\sum_{\omega_{2}\in\Omega_{jx}}\sum_{\omega_{3}\in\Omega_{iy}}\sum_{\omega_{4}\in\Omega_{jy}}\frac{\widetilde{\omega}(\omega_{1})}{p(\omega_{1})}\frac{\widetilde{\omega}(\omega_{2})}{p(\omega_{2})}\frac{\widetilde{\omega}(\omega_{3})}{p(\omega_{3})}\frac{\widetilde{\omega}(\omega_{4})}{p(\omega_{4})}\\\
\cdot\mathbb{I}(\omega_{1}\in\bar{\Omega}(k_{1},i))\mathbb{I}(\omega_{2}\in\bar{\Omega}(l_{1},j))\mathbb{I}(\omega_{3}\in\bar{\Omega}(k_{2},i))\mathbb{I}(\omega_{4}\in\bar{\Omega}(l_{2},j))\end{rcases}\text{different-
different (2)\hskip 18.49428pt}\\\
\begin{rcases}+m^{2}(m-1)\sum_{\omega_{1}\in\Omega_{ix}}\sum_{\omega_{2}\in\Omega_{jx}}\sum_{\omega_{3}\in\Omega_{iy}}\sum_{\omega_{4}\in\Omega_{jy}}\frac{\widetilde{\omega}(\omega_{1})}{p(\omega_{1})}\frac{\widetilde{\omega}(\omega_{2})}{p(\omega_{2})}\frac{\widetilde{\omega}(\omega_{3})}{p(\omega_{3})}\frac{\widetilde{\omega}(\omega_{4})}{p(\omega_{4})}\\\
\cdot\mathbb{I}(\omega_{1}\in\bar{\Omega}(k_{1},i))\mathbb{I}(\omega_{2}\in\bar{\Omega}(l_{1},j))\mathbb{I}(\omega_{3}\in\bar{\Omega}(k_{2},i))\mathbb{I}(\omega_{4}\in\bar{\Omega}(l_{1},j))\\\
+m^{2}(m-1)\sum_{\omega_{1}\in\Omega_{ix}}\sum_{\omega_{2}\in\Omega_{jx}}\sum_{\omega_{3}\in\Omega_{iy}}\sum_{\omega_{4}\in\Omega_{jy}}\frac{\widetilde{\omega}(\omega_{1})}{p(\omega_{1})}\frac{\widetilde{\omega}(\omega_{2})}{p(\omega_{2})}\frac{\widetilde{\omega}(\omega_{3})}{p(\omega_{3})}\frac{\widetilde{\omega}(\omega_{4})}{p(\omega_{4})}\\\
\cdot\mathbb{I}(\omega_{1}\in\bar{\Omega}(k_{1},i))\mathbb{I}(\omega_{2}\in\bar{\Omega}(l_{1},j))\mathbb{I}(\omega_{3}\in\bar{\Omega}(k_{1},i))\mathbb{I}(\omega_{4}\in\bar{\Omega}(l_{2},j)).\Bigl{.}\Bigr{\\}}\end{rcases}\text{same-
different (3) \hskip
34.1433pt}\end{multlined}(\bm{\phi}(i)^{\top}\bm{\phi}(j))^{2}=\frac{1}{m^{4}}\sum_{x,y\in\mathcal{V}}\Bigl{\\{}\Bigr{.}\\\
\begin{rcases}m^{2}\sum_{\omega_{1}\in\Omega_{ix}}\sum_{\omega_{2}\in\Omega_{jx}}\sum_{\omega_{3}\in\Omega_{iy}}\sum_{\omega_{4}\in\Omega_{jy}}\frac{\widetilde{\omega}(\omega_{1})}{p(\omega_{1})}\frac{\widetilde{\omega}(\omega_{2})}{p(\omega_{2})}\frac{\widetilde{\omega}(\omega_{3})}{p(\omega_{3})}\frac{\widetilde{\omega}(\omega_{4})}{p(\omega_{4})}\\\
\cdot\mathbb{I}(\omega_{1}\in\bar{\Omega}(k_{1},i))\mathbb{I}(\omega_{2}\in\bar{\Omega}(l_{1},j))\mathbb{I}(\omega_{3}\in\bar{\Omega}(k_{1},i))\mathbb{I}(\omega_{4}\in\bar{\Omega}(l_{1},j))\end{rcases}\text{same-
same (1)}\\\
\begin{rcases}+m^{2}(m-1)^{2}\sum_{\omega_{1}\in\Omega_{ix}}\sum_{\omega_{2}\in\Omega_{jx}}\sum_{\omega_{3}\in\Omega_{iy}}\sum_{\omega_{4}\in\Omega_{jy}}\frac{\widetilde{\omega}(\omega_{1})}{p(\omega_{1})}\frac{\widetilde{\omega}(\omega_{2})}{p(\omega_{2})}\frac{\widetilde{\omega}(\omega_{3})}{p(\omega_{3})}\frac{\widetilde{\omega}(\omega_{4})}{p(\omega_{4})}\\\
\cdot\mathbb{I}(\omega_{1}\in\bar{\Omega}(k_{1},i))\mathbb{I}(\omega_{2}\in\bar{\Omega}(l_{1},j))\mathbb{I}(\omega_{3}\in\bar{\Omega}(k_{2},i))\mathbb{I}(\omega_{4}\in\bar{\Omega}(l_{2},j))\end{rcases}\text{different-
different (2)\hskip 18.49428pt}\\\
\begin{rcases}+m^{2}(m-1)\sum_{\omega_{1}\in\Omega_{ix}}\sum_{\omega_{2}\in\Omega_{jx}}\sum_{\omega_{3}\in\Omega_{iy}}\sum_{\omega_{4}\in\Omega_{jy}}\frac{\widetilde{\omega}(\omega_{1})}{p(\omega_{1})}\frac{\widetilde{\omega}(\omega_{2})}{p(\omega_{2})}\frac{\widetilde{\omega}(\omega_{3})}{p(\omega_{3})}\frac{\widetilde{\omega}(\omega_{4})}{p(\omega_{4})}\\\
\cdot\mathbb{I}(\omega_{1}\in\bar{\Omega}(k_{1},i))\mathbb{I}(\omega_{2}\in\bar{\Omega}(l_{1},j))\mathbb{I}(\omega_{3}\in\bar{\Omega}(k_{2},i))\mathbb{I}(\omega_{4}\in\bar{\Omega}(l_{1},j))\\\
+m^{2}(m-1)\sum_{\omega_{1}\in\Omega_{ix}}\sum_{\omega_{2}\in\Omega_{jx}}\sum_{\omega_{3}\in\Omega_{iy}}\sum_{\omega_{4}\in\Omega_{jy}}\frac{\widetilde{\omega}(\omega_{1})}{p(\omega_{1})}\frac{\widetilde{\omega}(\omega_{2})}{p(\omega_{2})}\frac{\widetilde{\omega}(\omega_{3})}{p(\omega_{3})}\frac{\widetilde{\omega}(\omega_{4})}{p(\omega_{4})}\\\
\cdot\mathbb{I}(\omega_{1}\in\bar{\Omega}(k_{1},i))\mathbb{I}(\omega_{2}\in\bar{\Omega}(l_{1},j))\mathbb{I}(\omega_{3}\in\bar{\Omega}(k_{1},i))\mathbb{I}(\omega_{4}\in\bar{\Omega}(l_{2},j)).\Bigl{.}\Bigr{\\}}\end{rcases}\text{same-
different (3) \hskip 34.1433pt}\normalsize$ (37)

We now address each case $1-3$ in turn.

####  8.3.1 Case 1: $k_{1}=k_{2}$, $l_{1}=l_{2}$

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 $i$ and a single walk out of $j$, and all terms of case $1$ are
unchanged by the introduction of antithetic termination.

####  8.3.2 Case 2: $k_{1}\neq k_{2}$, $l_{1}\neq l_{2}$

Now we consider terms where both the walks out of node $i$ 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, we will need to consider the situations where either i)
one or ii) both of the pairs exhibit antithetic termination rather than
i.i.d.. Terms of both kind will appear when we use ensembles of antithetic
pairs. We need to check that in both cases the result is smaller compared to
when both pairs are i.i.d..

To evaluate these terms, we first need to understand how inducing antithetic
termination modifies the joint distribution
$p(\omega_{1}\in\bar{\Omega}(k_{1},i),\omega_{3}\in\bar{\Omega}(k_{2},i))$:
namely, the probability that two randomly sampled walks
$\bar{\Omega}(k_{1},i)$ and $\bar{\Omega}(k_{2},i)$ contain the respective
subwalks $\omega_{1}$ and $\omega_{3}$, given that their termination is either
i.i.d. or antithetic. In the i.i.d. case, it is straightforward to convince
oneself that

$p(\omega_{1}\in\bar{\Omega}(1,i),\omega_{3}\in\bar{\Omega}(3,i))=\left(\frac{1-p}{d}\right)^{m}\left(\frac{1-p}{d}\right)^{n},$
(38)

where $m$ and $n$ denote the lengths of subwalks $\omega_{1}$ and
$\omega_{3}$, respectively. With antithetic termination, from Eq. 13 it
follows that the probability of sampling a walk $\bar{\Omega}_{3}$ of length
$j$ conditioned on sampling an antithetic partner $\bar{\Omega}_{1}$ of length
$i$ is

$p(\textrm{len}(\bar{\Omega}_{3})=j|\textrm{len}(\bar{\Omega}_{1})=i)=\begin{cases}\left(\frac{1-2p}{1-p}\right)^{j}\frac{p}{1-p}&\textrm{if
}j<i,\\\ 0&\textrm{if }j=i,\\\
\left(\frac{1-2p}{1-p}\right)^{i}(1-p)^{j-i-1}p&\textrm{if }j>i.\end{cases}$
(39)

Using these probabilities, it is then straightforward but algebraically
tedious to derive the joint probabilities over subwalks

$p(\omega_{1}\in\bar{\Omega}(1,i),\omega_{3}\in\bar{\Omega}(3,i))=\begin{cases}\frac{1}{d^{m+n}}\left(\frac{1-2p}{1-p}\right)^{n}\left(1-p\right)^{m}&\textrm{if
}n<m,\\\ \frac{1}{d^{2m}}(1-2p)^{m}&\textrm{if }n=m,\\\
\frac{1}{d^{m+n}}\left(\frac{1-2p}{1-p}\right)^{m}\left(1-p\right)^{n}&\textrm{if
}n>m,\end{cases}$ (40)

where $m$ is the length of $\omega_{1}$, $n$ is the length of $\omega_{3}$ and
$i$ is now the index of a particular node.

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 $n<m$ as an example. Using Eq. 39,

$\small\begin{multlined}p(\omega_{1}\in\bar{\Omega}(1,i),\omega_{3}\in\bar{\Omega}(3,i))=\frac{1}{d^{m+n}}\sum_{i=m}^{\infty}\left[\sum_{j=n}^{i-1}\left(\frac{1-2p}{1-p}\right)^{j}\frac{p}{1-p}\left(1-p\right)^{i}p+\right.\\\
+\left.\sum_{j=i+1}^{\infty}\left(\frac{1-2p}{1-p}\right)^{i}(1-p)^{j-i-1}p\left(1-p\right)^{i}p\right],\end{multlined}p(\omega_{1}\in\bar{\Omega}(1,i),\omega_{3}\in\bar{\Omega}(3,i))=\frac{1}{d^{m+n}}\sum_{i=m}^{\infty}\left[\sum_{j=n}^{i-1}\left(\frac{1-2p}{1-p}\right)^{j}\frac{p}{1-p}\left(1-p\right)^{i}p+\right.\\\
+\left.\sum_{j=i+1}^{\infty}\left(\frac{1-2p}{1-p}\right)^{i}(1-p)^{j-i-1}p\left(1-p\right)^{i}p\right],\normalsize$
(41)

where the branching factors of $d$ appeared because at 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, so we stop accumulating these multiplicative
factors. Computing the summations in Eq. 41 (which are all straightforward
geometric series), we quickly arrive at the top line of Eq. 40.

Returning to our main discussion, note that in the $d$-regular, equal-weights
case,

$\small\begin{multlined}\sum_{\omega_{1}\in\Omega_{ix}}\sum_{\omega_{3}\in\Omega_{iy}}\frac{\widetilde{\omega}(\omega_{1})}{p(\omega_{1})}\frac{\widetilde{\omega}(\omega_{3})}{p(\omega_{3})}p(\omega_{1}\in\bar{\Omega}(k_{1},i),\omega_{3}\in\bar{\Omega}(k_{2},i))\\\
=\sum_{\omega_{1}\in\Omega_{ix}}\sum_{\omega_{3}\in\Omega_{iy}}\left(\frac{wd}{1-p}\right)^{m+n}p(\omega_{1}\in\bar{\Omega}(k_{1},i),\omega_{3}\in\bar{\Omega}(k_{2},i)).\end{multlined}\sum_{\omega_{1}\in\Omega_{ix}}\sum_{\omega_{3}\in\Omega_{iy}}\frac{\widetilde{\omega}(\omega_{1})}{p(\omega_{1})}\frac{\widetilde{\omega}(\omega_{3})}{p(\omega_{3})}p(\omega_{1}\in\bar{\Omega}(k_{1},i),\omega_{3}\in\bar{\Omega}(k_{2},i))\\\
=\sum_{\omega_{1}\in\Omega_{ix}}\sum_{\omega_{3}\in\Omega_{iy}}\left(\frac{wd}{1-p}\right)^{m+n}p(\omega_{1}\in\bar{\Omega}(k_{1},i),\omega_{3}\in\bar{\Omega}(k_{2},i)).\normalsize$
(42)

The summand depends only on walk lengths $m,n$ but not direction, which
invites us to decompose the sum
$\sum_{\omega_{1}\in\Omega_{ix}}\left(\cdot\right)$ over paths between nodes
$i$ and $x$ to a sum over path lengths, weighted by the number of paths at
each length. Explicitly,

$\sum_{\omega_{1}\in\Omega_{ix}}\left(\cdot\right)=\sum_{n=1}^{\infty}(\mathbf{A}^{n})_{ix}\left(\cdot\right),$
(43)

with $\mathbf{A}$ the (unweighted) adjacency matrix. We have used the fact
that $(\mathbf{A}^{n})_{ij}$ counts the number of walks of length $n$ between
nodes $i$ and $x$. $\mathbf{A}$ is symmetric so has a convenient decomposition
into orthogonal eigenvectors and real eigenvalues:

$(\mathbf{A}^{n})_{ix}=\sum_{k=1}^{N}\lambda_{k}^{n}\braket{i}{k}\braket{k}{x}$
(44)

where $\ket{k}$ enumerates the $N$ eigenvectors of $\mathbf{A}$ with
corresponding eigenvalues $\lambda_{k}$, and $\bra{i}$ and $\bra{x}$ are unit
vectors in the $i$ and $x$ coordinate axes, respectively. We remind the reader
that we have adopted Dirac’s bra-ket notation; $\ket{y}$ denotes the vector
$\bm{y}$ and $\bra{y}$ denotes $\bm{y}^{\top}$.

Inserting Eqs 44 and 43 into Eq. 42 and using the probability distributions in
Eq. 38 and 40, our all-important variance-determining correlation term from
Eq.36 evaluates to

$\sum_{x,y\in\mathcal{V}}\sum_{k_{1},k_{2},k_{3},k_{4}=1}^{N}B^{(i)}_{k_{1},k_{3}}B^{(j)}_{k_{2},k_{4}}\braket{i}{k_{1}}\braket{k_{1}}{x}\braket{j}{k_{2}}\braket{k_{2}}{x}\braket{i}{k_{3}}\braket{k_{3}}{y}\braket{j}{k_{4}}\braket{k_{4}}{y},$
(45)

where the matrix elements $B^{(i)}_{k_{1},k_{3}}$ and $B^{(j)}_{k_{2},k_{4}}$,
corresponding to the pairs of walkers out of $i$ and $j$ respectively, are
equal to one of the two following expressions:

$B_{k_{1},k_{3}}=\begin{cases}C_{k_{1},k_{3}}\coloneqq\frac{w\lambda_{k_{1}}}{1-w\lambda_{k_{1}}}\frac{w\lambda_{k_{3}}}{1-w\lambda_{k_{3}}}&\textrm{if
i.i.d.}\\\
D_{k_{1},k_{3}}\coloneqq\frac{w\lambda_{k_{1}}}{1-w\lambda_{k_{1}}}\frac{w\lambda_{k_{3}}}{1-w\lambda_{k_{3}}}\frac{c(1-w^{2}\lambda_{k_{1}}\lambda_{k_{3}})}{1-cw^{2}\lambda_{k_{1}}\lambda_{k_{3}}}&\textrm{if
antithetic.}\end{cases}$ (46)

Here, $c$ is a constant defined by $c\coloneqq\frac{1-2p}{(1-p)^{2}}$ with $p$
the termination probability. These forms are straightforward to compute with
good algebraic bookkeeping; we omit details for economy of space.

Eq. 45 can be simplified. Observe that
$\sum_{x\in\mathcal{V}}\ket{x}\bra{x}=\mathbf{I}_{N}$ (‘resolution of the
identity’), and that since the eigenvectors of $\mathbf{A}$ are orthogonal
$\braket{k_{1}}{k_{2}}=\delta_{k_{1},k_{2}}$. Applying this, we can write

$\sum_{k_{1},k_{3}=1}^{N}B^{(i)}_{k_{1},k_{3}}B^{(j)}_{k_{1},k_{3}}\braket{i}{k_{1}}\braket{j}{k_{1}}\braket{i}{k_{3}}\braket{j}{k_{3}}.$
(47)

Our task is then to determine whether 47 is reduced by conditioning that
either one or both of the pairs of walkers are antithetic rather than
independent. That is,

$\sum_{k_{1}=1}^{N}\sum_{k_{3}=1}^{N}\left(C_{k_{1},k_{3}}D_{k_{1},k_{3}}-C_{k_{1},k_{3}}C_{k_{1},k_{3}}\right)\braket{i}{k_{1}}\braket{j}{k_{1}}\braket{i}{k_{3}}\braket{j}{k_{3}}\overset{\textrm{?}}{\leq}0,$
(48)
$\sum_{k_{1}=1}^{N}\sum_{k_{3}=1}^{N}\left(D_{k_{1},k_{3}}D_{k_{1},k_{3}}-C_{k_{1},k_{3}}C_{k_{1},k_{3}}\right)\braket{i}{k_{1}}\braket{j}{k_{1}}\braket{i}{k_{3}}\braket{j}{k_{3}}\overset{\textrm{?}}{\leq}0.$
(49)

Define a vector $\bm{y}\in\mathbb{R}^{N}$ with entries
$y_{p}\coloneqq\braket{i}{k_{p}}\braket{j}{k_{p}}$, such that its $p$th
element is the product of the $i$ and $j$th coordinates of the $p$th
eigenvector $\bm{k}_{p}$. In this notation, Eqs 48 and 49 can be written

$\sum_{p=1}^{N}\sum_{q=1}^{N}\left(C_{pq}D_{pq}-C_{pq}C_{pq}\right)y_{p}y_{q}\overset{\textrm{?}}{\leq}0,$
(50)
$\sum_{p=1}^{N}\sum_{q=1}^{N}\left(D_{pq}D_{pq}-C_{pq}C_{pq}\right)y_{p}y_{q}\overset{\textrm{?}}{\leq}0.$
(51)

For Eqs 50 and 51 to be true for arbitrary graphs, it is sufficient that the
matrices $\mathbf{E}$ and $\mathbf{F}$ with matrix elements $E_{pq}\coloneqq
C_{pq}D_{pq}-C_{pq}C_{pq}$ and $F_{pq}\coloneqq D_{pq}D_{pq}-C_{pq}C_{pq}$ are
_negative definite_. Our next task is to prove that this is the case.

First, consider $\mathbf{E}$, where just one of the two pairs of walkers is
antithetic. Putting in the explicit forms of $C_{pq}$ and $D_{pq}$ from Eq.
46,

$\begin{multlined}E_{pq}=-\left(\frac{\bar{\lambda}_{p}\bar{\lambda}_{q}}{(1-\bar{\lambda}_{p})(1-\bar{\lambda}_{q})}\right)^{2}\frac{p^{2}}{(1-p)^{2}}\frac{1}{1-\frac{1-2p}{(1-p)^{2}}\bar{\lambda}_{p}\bar{\lambda}_{q}}\end{multlined}E_{pq}=-\left(\frac{\bar{\lambda}_{p}\bar{\lambda}_{q}}{(1-\bar{\lambda}_{p})(1-\bar{\lambda}_{q})}\right)^{2}\frac{p^{2}}{(1-p)^{2}}\frac{1}{1-\frac{1-2p}{(1-p)^{2}}\bar{\lambda}_{p}\bar{\lambda}_{q}}$
(52)

where for notational compactness we took $\bar{\lambda}_{p}\coloneqq
w\lambda_{p}$ (the eigenvalues of the _weighted_ adjacency matrix
$\mathbf{U}$). Taylor expanding,

$E_{pq}=-\left(\frac{\bar{\lambda}_{p}\bar{\lambda}_{q}}{(1-\bar{\lambda}_{p})(1-\bar{\lambda}_{q})}\right)^{2}\frac{p^{2}}{(1-p)^{2}}\sum_{m=0}^{\infty}\left(\frac{1-2p}{(1-p)^{2}}\bar{\lambda}_{p}\bar{\lambda}_{q}\right)^{m}.$
(53)

Inserting this into Eq. 50, we get

$\sum_{p=1}^{N}\sum_{q=1}^{N}E_{pq}y_{p}y_{q}=-\frac{p^{2}}{(1-p)^{2}}\sum_{m=0}^{\infty}\left(\sum_{p=1}^{N}\frac{\bar{\lambda}_{p}^{2}}{(1-\bar{\lambda}_{p})^{2}}\left(\frac{\sqrt{1-2p}}{1-p}\bar{\lambda}_{p}\right)^{m}y_{p}\right)^{2}\leq
0,$ (54)

which implies that $\mathbf{E}$ is indeed negative definite. Note that we have
not made any additional assumptions about the values of $p$ and $w$ beyond
those already stipulated: namely, $0<p\leq\frac{1}{2}$ and
$\bar{\lambda}_{\textrm{max}}<1$.

Next, consider $\mathbf{F}$, where both pairs of walkers are antithetic. Again
inserting Eqs 46, we find that

$F_{pq}=\left(\frac{\bar{\lambda}_{p}\bar{\lambda}_{q}}{(1-\bar{\lambda}_{p})(1-\bar{\lambda}_{q})}\right)^{2}\left[\left(\frac{c-c\bar{\lambda}_{p}\bar{\lambda}_{q}}{1-c\bar{\lambda}_{p}\bar{\lambda}_{q}}\right)^{2}-1\right]$
(55)

where we remind the reader that $c=\frac{1-2p}{(1-p)^{2}}$. The Taylor
expansion in $\bar{\lambda}_{p}\bar{\lambda}_{q}$ is

$\begin{multlined}F_{pq}=\left(\frac{\bar{\lambda}_{p}\bar{\lambda}_{q}}{(1-\bar{\lambda}_{p})(1-\bar{\lambda}_{q})}\right)^{2}\left[\sum_{i=0}^{\infty}(\bar{\lambda}_{p}\bar{\lambda}_{q})^{i}(c-1)c^{i}(1+c+i(c-1))\right]\\\
=w^{4}(\lambda_{p}\lambda_{q})^{2}\sum_{i,j,k=0}^{\infty}(\lambda_{p}\lambda_{q})^{i+j+k}w^{2i+j+k}(c-1)c^{i}(1+c+i(c-1))(j+1)(k+1).\end{multlined}F_{pq}=\left(\frac{\bar{\lambda}_{p}\bar{\lambda}_{q}}{(1-\bar{\lambda}_{p})(1-\bar{\lambda}_{q})}\right)^{2}\left[\sum_{i=0}^{\infty}(\bar{\lambda}_{p}\bar{\lambda}_{q})^{i}(c-1)c^{i}(1+c+i(c-1))\right]\\\
=w^{4}(\lambda_{p}\lambda_{q})^{2}\sum_{i,j,k=0}^{\infty}(\lambda_{p}\lambda_{q})^{i+j+k}w^{2i+j+k}(c-1)c^{i}(1+c+i(c-1))(j+1)(k+1).$
(56)

In fact, $\mathbf{F}$ is _not_ generically negative definite, but will be at
sufficiently small $p$ or $w$. Write
$\mathbf{F}=w^{4}(\mathbf{G}+\mathbf{H})$, with

$G_{pq}\coloneqq(\lambda_{p}\lambda_{q})^{2}\left(c^{2}-1\right),$ (57)
$H_{pq}\coloneqq(\lambda_{p}\lambda_{q})^{2}\sum_{i,j,k=0\backslash\\{i=j=k=0\\}}^{\infty}(\lambda_{p}\lambda_{q})^{i+j+k}w^{2i+j+k}(c-1)c^{i}(1+c+i(c-1))(j+1)(k+1).$
(58)

$\mathbf{G}$ is manifestly negative definite because $c<1$ but $\mathbf{H}$
may not be. Treat $\mathbf{H}$ as a perturbation to $\mathbf{G}$.

Recalling that the spectral radius of $\mathbf{H}$ is defined

$\rho(\mathbf{H})\coloneqq\underset{\|\bm{x}\|_{2}=1}{\textrm{max}}\mathbf{H}\bm{x},$
(59)

it is clear that the spectral radius of $\mathbf{H}$ approaches $0$ smoothly
as $w\to 0$ since all its matrix elements vanish. Recall also an important
corollary of Weyl’s perturbation inequality: any perturbed eigenvalue of
$\mathbf{F}+\mathbf{G}$ will be within one spectral radius $\rho(\mathbf{G})$
of the original eigenvalue of $\mathbf{F}$. This means that, by reducing $w$,
we can shrink the spectral radius of $\mathbf{G}$ until
$\rho(\mathbf{G})<(\lambda_{p}\lambda_{q})^{2}\left(1-c^{2}\right)$, at which
point we are guaranteed that $\mathbf{F}$ will be negative definite. Hence, at
sufficiently small $w$, correlation terms with both pairs antithetic are
suppressed as required.

Taylor expanding in $c\to 1$ (which corresponds to $p\to 0$) instead of
$\lambda_{p}\lambda_{q}$, we can make exactly analogous arguments to find that
$\mathbf{F}$ is also guaranteed to be negative definite with when $p$ is
sufficiently small. Briefly: let $c=1-\delta$ with
$\delta=\left(\frac{p}{1-p}\right)^{2}$. Then we have that

$\begin{multlined}F_{pq}=\left(\frac{\bar{\lambda}_{p}\bar{\lambda}_{q}}{(1-\bar{\lambda}_{p})(1-\bar{\lambda}_{q})}\right)^{2}\left(\left(1-\delta\right)^{2}\left(\frac{1-\bar{\lambda}_{p}\bar{\lambda}_{q}}{1-\bar{\lambda}_{p}\bar{\lambda}_{q}+\delta\bar{\lambda}_{p}\bar{\lambda}_{q}}\right)^{2}-1\right)\\\
=\left(\frac{\bar{\lambda}_{p}\bar{\lambda}_{q}}{(1-\bar{\lambda}_{p})(1-\bar{\lambda}_{q})}\right)^{2}\left(\frac{-2\delta}{1-\bar{\lambda}_{p}\bar{\lambda}_{q}}+\mathcal{O}(\delta^{2})\right).\end{multlined}F_{pq}=\left(\frac{\bar{\lambda}_{p}\bar{\lambda}_{q}}{(1-\bar{\lambda}_{p})(1-\bar{\lambda}_{q})}\right)^{2}\left(\left(1-\delta\right)^{2}\left(\frac{1-\bar{\lambda}_{p}\bar{\lambda}_{q}}{1-\bar{\lambda}_{p}\bar{\lambda}_{q}+\delta\bar{\lambda}_{p}\bar{\lambda}_{q}}\right)^{2}-1\right)\\\
=\left(\frac{\bar{\lambda}_{p}\bar{\lambda}_{q}}{(1-\bar{\lambda}_{p})(1-\bar{\lambda}_{q})}\right)^{2}\left(\frac{-2\delta}{1-\bar{\lambda}_{p}\bar{\lambda}_{q}}+\mathcal{O}(\delta^{2})\right).$
(60)

Taylor expanding $\frac{1}{1-\bar{\lambda}_{p}\bar{\lambda}_{q}}$, it is easy
to see that the operator defined by the $\mathcal{O}(\delta)$ term of Eq. 60
is negative definite. This part will dominate over higher order terms (which
are _not_ in general negative definite) when $\delta$ is sufficiently small,
guaranteeing the effectiveness of our mechanism on these terms.

As an aside, we also note that Taylor expanding about $c=0$ (which corresponds
to $p\to\frac{1}{2}$) yields

$F_{pq}=\left(\frac{\bar{\lambda}_{p}\bar{\lambda}_{q}}{(1-\bar{\lambda}_{p})(1-\bar{\lambda}_{q})}\right)^{2}\left(-1+\mathcal{O}(c^{2})\right)$
(61)

which is manifestly negative definite at small enough $c$. Hence,
intriguingly, the $k_{1}\neq k_{2}$ variance contributions are also suppressed
in the $p\to\frac{1}{2}$ limit.

This concludes our study of variance contributions in Eq. 36 where $k_{1}\neq
k_{2}$, $l_{1}\neq l_{2}$. We have found that these correlation terms are
indeed suppressed by antithetic termination when $p$ or $\rho(\mathbf{U})$ is
small enough (or when $p$ is sufficiently close to $\frac{1}{2}$).

####  8.3.3 Case 3: $k_{1}=k_{2}$, $l_{1}\neq l_{2}$

We now consider terms where $k_{1}=k_{2}$ and $l_{1}\neq l_{2}$. We are
considering a total of $3$ walks: just $1$ out of node $i$ but a pair (which
may be antithetic or i.i.d.) out of node $j$. We inspect the term

$\sum_{\omega_{1}\in\Omega_{ix}}\sum_{\omega_{3}\in\Omega_{iy}}\left(\frac{wd}{1-p}\right)^{m+n}p(\omega_{1}\in\bar{\Omega}(k_{1},i),\omega_{3}\in\bar{\Omega}(k_{1},i)),$
(62)

where $m$ denotes the length of $\omega_{1}$ and $n$ denotes the length of
$\omega_{3}$. What is the form of
$p(\omega_{1}\in\bar{\Omega}(k_{1},i),\omega_{3}\in\bar{\Omega}(k_{1},i))$? It
is the probability that a _single_ walk out of node $i$,
$\bar{\Omega}(k_{1},i)$, contains walks $\omega_{1}$ between nodes $i$ and $x$
and $\omega_{3}$ between $i$ and $y$ as subwalks. Such a walk must pass
through all three nodes $i$, $x$ and $y$. After some thought,

$p(\omega_{1},\omega_{3}\in\bar{\Omega}(k_{1},i))=\begin{cases}\left(\frac{1-p}{d}\right)^{m}&\textrm{if
}\omega_{1}=\omega_{3},\\\ \left(\frac{1-p}{d}\right)^{m}&\textrm{if
}\omega_{3}\in\omega_{1},\\\ \left(\frac{1-p}{d}\right)^{n}&\textrm{if
}\omega_{1}\in\omega_{3},\\\ 0&\textrm{otherwise.}\end{cases}$ (63)

Here, $\omega_{1}\in\omega_{3}$ means $\omega_{1}$ is a strict subwalk of
$\omega_{3}$, so the sequence of nodes traversed is $i\to x\to y$. Likewise,
$\omega_{3}\in\omega_{1}$ implies a path $i\to y\to x$. Summing these
contributions,

$\small\begin{multlined}\sum_{\omega_{1}\in\Omega_{ix}}\sum_{\omega_{3}\in\Omega_{iy}}\left(\frac{wd}{1-p}\right)^{m+n}p(\omega_{1}\in\bar{\Omega}(l_{1},i),\omega_{3}\in\bar{\Omega}(l_{1},i))\\\
=\underbrace{\sum_{\omega_{1}\in\Omega_{ix}}\left(\frac{wd}{1-p}\right)^{2\textrm{len}(\omega_{1})}p(\omega_{1}\in\bar{\Omega}(k_{1},i))\delta_{xy}}_{\omega_{1}=\omega_{3},i\to
x=y}\\\
+\underbrace{\sum_{\omega_{1}\in\Omega_{ix}}\left(\frac{wd}{1-p}\right)^{2\textrm{len}(\omega_{1})}p(\omega_{1}\in\bar{\Omega}(k_{1},i))\sum_{\omega_{\delta}\in\Omega_{xy}}\left(\frac{wd}{1-p}\right)^{\textrm{len}(\omega_{\delta})}p(\omega_{\delta}\in\bar{\Omega}(k_{1},x))}_{\omega_{1}\in\omega_{3},i\to
x\to y}\\\
+\underbrace{\sum_{\omega_{3}\in\Omega_{iy}}\left(\frac{wd}{1-p}\right)^{2\textrm{len}(\omega_{3})}p(\omega_{3}\in\bar{\Omega}(k_{1},i))\sum_{\omega_{\delta}\in\Omega_{yx}}\left(\frac{wd}{1-p}\right)^{\textrm{len}(\omega_{\delta})}p(\omega_{\delta}\in\bar{\Omega}(k_{1},y))}_{\omega_{3}\in\omega_{1},i\to
y\to
x}.\end{multlined}\sum_{\omega_{1}\in\Omega_{ix}}\sum_{\omega_{3}\in\Omega_{iy}}\left(\frac{wd}{1-p}\right)^{m+n}p(\omega_{1}\in\bar{\Omega}(l_{1},i),\omega_{3}\in\bar{\Omega}(l_{1},i))\\\
=\underbrace{\sum_{\omega_{1}\in\Omega_{ix}}\left(\frac{wd}{1-p}\right)^{2\textrm{len}(\omega_{1})}p(\omega_{1}\in\bar{\Omega}(k_{1},i))\delta_{xy}}_{\omega_{1}=\omega_{3},i\to
x=y}\\\
+\underbrace{\sum_{\omega_{1}\in\Omega_{ix}}\left(\frac{wd}{1-p}\right)^{2\textrm{len}(\omega_{1})}p(\omega_{1}\in\bar{\Omega}(k_{1},i))\sum_{\omega_{\delta}\in\Omega_{xy}}\left(\frac{wd}{1-p}\right)^{\textrm{len}(\omega_{\delta})}p(\omega_{\delta}\in\bar{\Omega}(k_{1},x))}_{\omega_{1}\in\omega_{3},i\to
x\to y}\\\
+\underbrace{\sum_{\omega_{3}\in\Omega_{iy}}\left(\frac{wd}{1-p}\right)^{2\textrm{len}(\omega_{3})}p(\omega_{3}\in\bar{\Omega}(k_{1},i))\sum_{\omega_{\delta}\in\Omega_{yx}}\left(\frac{wd}{1-p}\right)^{\textrm{len}(\omega_{\delta})}p(\omega_{\delta}\in\bar{\Omega}(k_{1},y))}_{\omega_{3}\in\omega_{1},i\to
y\to x}.\normalsize$ (64)

We introduced $\omega_{\delta}$ for the sum over paths between nodes $x$ and
$y$, and $p(\omega_{\delta}\in\bar{\Omega}(k_{1},x))$ is the probability of
some particular subwalk $x\to y$, equal to
$\left(\frac{1-p}{d}\right)^{\textrm{len}(\omega_{\delta})}$ in the
$d$-regular case. $\omega_{3}$ is a dummy variable so can be relabelled
$\omega_{1}$. The variance-determining correlation term from Eq. 36 becomes

$\small\begin{multlined}\sum_{x,y\in\mathcal{V}}\left[\sum_{\omega_{1}\in\Omega_{ix}}\left(\frac{wd}{1-p}\right)^{2\textrm{len}(\omega_{1})}p(\omega_{1}\in\bar{\Omega}(k_{1},i))\delta_{xy}\right.\\\
+\sum_{\omega_{1}\in\Omega_{ix}}\left(\frac{wd}{1-p}\right)^{2\textrm{len}(\omega_{1})}p(\omega_{1}\in\bar{\Omega}(k_{1},i))\sum_{\omega_{\delta}\in\Omega_{xy}}\left(\frac{wd}{1-p}\right)^{\textrm{len}(\omega_{\delta})}p(\omega_{\delta}\in\bar{\Omega}(k_{1},x))\\\
\left.+\sum_{\omega_{1}\in\Omega_{iy}}\left(\frac{wd}{1-p}\right)^{2\textrm{len}(\omega_{1})}p(\omega_{1}\in\bar{\Omega}(k_{1},i))\sum_{\omega_{\delta}\in\Omega_{yx}}\left(\frac{wd}{1-p}\right)^{\textrm{len}(\omega_{\delta})}p(\omega_{\delta}\in\bar{\Omega}(k_{1},y))\right]\\\
\cdot\sum_{k_{2}=1}^{N}\sum_{k_{4}=1}^{N}B_{k_{2},k_{4}}^{(j)}\braket{j}{k_{2}}\braket{k_{2}}{x}\braket{j}{k_{4}}\braket{k_{4}}{y}.\end{multlined}\sum_{x,y\in\mathcal{V}}\left[\sum_{\omega_{1}\in\Omega_{ix}}\left(\frac{wd}{1-p}\right)^{2\textrm{len}(\omega_{1})}p(\omega_{1}\in\bar{\Omega}(k_{1},i))\delta_{xy}\right.\\\
+\sum_{\omega_{1}\in\Omega_{ix}}\left(\frac{wd}{1-p}\right)^{2\textrm{len}(\omega_{1})}p(\omega_{1}\in\bar{\Omega}(k_{1},i))\sum_{\omega_{\delta}\in\Omega_{xy}}\left(\frac{wd}{1-p}\right)^{\textrm{len}(\omega_{\delta})}p(\omega_{\delta}\in\bar{\Omega}(k_{1},x))\\\
\left.+\sum_{\omega_{1}\in\Omega_{iy}}\left(\frac{wd}{1-p}\right)^{2\textrm{len}(\omega_{1})}p(\omega_{1}\in\bar{\Omega}(k_{1},i))\sum_{\omega_{\delta}\in\Omega_{yx}}\left(\frac{wd}{1-p}\right)^{\textrm{len}(\omega_{\delta})}p(\omega_{\delta}\in\bar{\Omega}(k_{1},y))\right]\\\
\cdot\sum_{k_{2}=1}^{N}\sum_{k_{4}=1}^{N}B_{k_{2},k_{4}}^{(j)}\braket{j}{k_{2}}\braket{k_{2}}{x}\braket{j}{k_{4}}\braket{k_{4}}{y}.\normalsize$
(65)

where $B_{k_{2},k_{4}}^{(j)}$ depends on whether the coupling of the pair of
walkers out of node $j$ is i.i.d. or antithetic, as defined in Eq. 46. $x$ and
$y$ are dummy variables so can also be swapped, and the sum over the paths
$\omega_{\delta}$ is computed via the usual sum over path lengths and
eigendecomposition of $\mathbf{A}$. Using the resolution of the identity and
working through the algebra, we obtain the correlation term

$\begin{multlined}\sum_{x\in\mathcal{V}}\left[\sum_{\omega_{1}\in\Omega_{ix}}\left(\frac{wd}{1-p}\right)^{2\textrm{len}(\omega_{1})}p(\omega_{1}\in\bar{\Omega}(k_{1},i))\right]\\\
\cdot\sum_{k_{2},k_{4}=1}^{N}\left(\frac{1-w^{2}\lambda_{k_{2}}\lambda_{k_{4}}}{(1-w\lambda_{k_{2}})(1-w\lambda_{k_{4}})}\right)B^{(j)}_{k_{2},k_{4}}\braket{x}{k_{2}}\braket{k_{2}}{j}\braket{x}{k_{4}}\braket{k_{4}}{j}.\end{multlined}\sum_{x\in\mathcal{V}}\left[\sum_{\omega_{1}\in\Omega_{ix}}\left(\frac{wd}{1-p}\right)^{2\textrm{len}(\omega_{1})}p(\omega_{1}\in\bar{\Omega}(k_{1},i))\right]\\\
\cdot\sum_{k_{2},k_{4}=1}^{N}\left(\frac{1-w^{2}\lambda_{k_{2}}\lambda_{k_{4}}}{(1-w\lambda_{k_{2}})(1-w\lambda_{k_{4}})}\right)B^{(j)}_{k_{2},k_{4}}\braket{x}{k_{2}}\braket{k_{2}}{j}\braket{x}{k_{4}}\braket{k_{4}}{j}.$
(66)

Now observe that the prefactor in square brackets is positive for any node $x$
since it is the expectation of a squared quantity. This means that, for the
sum in Eq. 66 to be suppressed by antithetic coupling, it is sufficient for
the summation in its lower line to be reduced. Defining a vector
$\bm{y}\in\mathbb{R}^{N}$ with elements
$y_{p}\coloneqq\braket{x}{k_{p}}\braket{k_{p}}{j}$, it becomes clear that we
require that the operator $\mathbf{J}$ with matrix elements

$J_{pq}\coloneqq\left(\frac{1-w^{2}\lambda_{p}\lambda_{q}}{(1-w\lambda_{p})(1-w\lambda_{q})}\right)(D_{pq}-C_{pq})$
(67)

is negative definite. Using the forms in Eq. 46,

$J_{pq}=-\frac{w^{2}\lambda_{p}\lambda_{q}}{(1-w\lambda_{p})^{2}(1-w\lambda_{q})^{2}}\frac{\frac{p^{2}}{(1-p)^{2}}(1-w^{2}\lambda_{p}\lambda_{q})}{1-\frac{1-2p}{(1-p)^{2}}w^{2}\lambda_{p}\lambda_{q}}.$
(68)

Making very similar arguments to in Sec. 8.3.2 (namely, Taylor expanding and
appealing to 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$.

A brief note: Taylor expanding in $c$,

$J_{pq}=-\frac{w^{2}\lambda_{p}\lambda_{q}}{(1-w\lambda_{p})^{2}(1-w\lambda_{q})^{2}}(1-w^{2}\lambda_{p}\lambda_{q})+\mathcal{O}(c),$
(69)

which is only negative definite when we also simultaneously take $w\to 0$.
Interestingly, in contrast to case $2$, these terms are _not_ suppressed by
$p\to\frac{1}{2}$ on its own; we need to control the spectral radius of
$\mathbf{U}$.

This concludes the section of the proof addressing terms $k_{1}=k_{2}$ and
$l_{1}\neq l_{2}$ (case $3$). Again, these variance contributions are always
suppressed by antithetic termination at sufficiently small $p$ or
$\rho(\mathbf{U})$.

Having now considered all the possible variance contributions enumerated by
cases $1-3$ and shown that each is either reduced or unmodified by the
imposition of antithetic termination, we can finally conclude that our novel
mechanism does indeed suppress the $2$-regularised Laplacian kernel estimator
variance for a $d$-regular graph of equal weights at sufficiently small $p$ or
$\rho(\mathbf{U})$. ∎

As mentioned in the main body of the manuscript, these conditions tend not to
be very restrictive in experiments. Intriguingly, small $\rho(\mathbf{U})$
with $p=\frac{1}{2}$ actually works very well.

Our next task is to generalise these results to broader classes of graphs.

###  8.4 Extending the results to arbitrary graphs and sampling strategies
(Theorem 3.2 cont.)

Throughout Sec. 8.3, we considered the simplest setting of a $d$-regular graph
where all edges have equal weight. We have also taken a basic sampling
strategy, with the walker choosing one of its current node’s neighbours at
random at every timestep. Here we relax these assumptions, showing that our
results remain true in more general settings.

####  8.4.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.3.2) are unmodified because taking
$d^{m}\to\prod_{i=1}^{m}d_{i}$ in $p(\omega_{1})$ and
$d^{n}\to\prod_{i=1}^{n}d_{i}$ in $p(\omega_{3})$ is exactly compensated by
the corresponding change in in joint probability
$p(\omega_{1}\in\bar{\Omega}(k_{1},i),\omega_{3}\in\bar{\Omega}(k_{2},i))$.
Our previous arguments all continue to hold.

Case $3$ (Sec. 8.3.3) is only a little harder. Now the prefactor in square
parentheses in the top line of Eq. 66 evaluates to

$\left[\sum_{\omega_{1}\in\Omega_{ix}}\left(\frac{w}{1-p}\right)^{2\textrm{len}(\omega_{1})}\left(\prod_{i=1}^{\textrm{len}(\omega_{1})}d_{i}^{2}\right)p(\omega_{1}\in\bar{\Omega}(k_{1},i))\right]$
(70)

which is still positive for any node $x$. The lower line of Eq. 66 is
unmodified because once again the change $d^{m}\to\prod_{i=1}^{m}d_{i}$
exactly cancels in the marginal and joint probabilities, so $\mathbf{J}$ is
unchanged and our previous conclusions prevail.

####  8.4.2 Weighted graphs

Now we permit edge weights to differ across the graph. Once again, case $2$
(Sec. 8.3.2) is straightforward: instead of Eq. 43, we take

$\sum_{\omega_{1}\in\Omega_{ix}}\widetilde{\omega}(\omega_{1})\left(\cdot\right)=\sum_{n=1}^{\infty}(\mathbf{U}^{n})_{ix}\left(\cdot\right),$
(71)

where $\mathbf{U}$ is the _weighted_ adjacency matrix. We incorporate the
product of each walk’s edge weights into the combinatorial factor, then sum
over path lengths as before. In downstream calculations we drop all instances
of $w$ and reinterpret $\lambda$ as the eigenvalues of the $\mathbf{U}$
instead of $\mathbf{A}$, but our arguments are otherwise unmodified; these
variance contributions will be suppressed if $\rho(\mathbf{U})$ or $p$ is
sufficiently small.

Case $3$ (Sec. 8.3.3) is also easy enough; the bracketed prefactor of 66
becomes

$\left[\sum_{\omega_{1}\in\Omega_{ix}}\left(\frac{1}{1-p}\right)^{2\textrm{len}(\omega_{1})}\left(\prod_{i=1}^{\textrm{len}(\omega_{1})}w_{i\sim
i+1}^{2}d_{i}^{2}\right)p(\omega_{1}\in\bar{\Omega}(k_{1},i))\right]$ (72)

which is again positive. Here, $w_{i\sim i+1}$ denotes the weight associated
with the edge between the $i$ and $i+1$th nodes of the walk. Therefore, it is
sufficient that the matrix $\mathbf{J}$ with matrix elements

$J_{pq}=-\frac{\lambda_{p}\lambda_{q}}{(1-\lambda_{p})^{2}(1-\lambda_{q})^{2}}\frac{\frac{p^{2}}{(1-p)^{2}}(1-\lambda_{p}\lambda_{q})}{1-\frac{1-2p}{(1-p)^{2}}\lambda_{p}\lambda_{q}}$
(73)

is negative definite, with $\lambda_{p}$ now the $p$th eigenvalue of the
_weighted_ adjacency matrix $\mathbf{U}$. Following the same arguments as in
Sec. 8.3.3, this will be the case at small enough $p$ or $\rho(\mathbf{U})$.

####  8.4.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 $\frac{1}{d_{i}}$), but the transition
probability can also be a function of the edge weights. For example, if all
the edge weights are positive, we might take

$p(i\to j|\bar{s})=\frac{w_{ij}}{\sum_{k\sim i}w_{ik}}$ (74)

for the probability of transitioning from node $i$ to $j$ at a given timestep
(with $w_{ij}\coloneqq\mathbf{U}_{ij}$), given that the walker does not
terminate. This strategy increases the probability of taking edges with bigger
weights and which therefore contribute more to
$(\mathbf{I}_{N}-\mathbf{U})^{-2}$ – something that empirically suppresses the
variance on the estimator of the $2$-regularised Laplacian kernel. Does
antithetic termination reduce it further?

Case $2$ (Sec. 8.3.2) is again easy; the $w$-dependent modifications to
$p(\omega_{1})$ and $p(\omega_{3})$ are exactly compensated by adjustments to
$p(\omega_{1}\in\bar{\Omega}(k_{1},i),\omega_{3}\in\bar{\Omega}(k_{2},i))$. To
wit, Eq. 40 becomes

$p(\omega_{3}\in\bar{\Omega}(3,i),\omega_{1}\in\bar{\Omega}(1,i))=\begin{cases}\frac{\widetilde{\omega}(\omega_{1})\widetilde{\omega}(\omega_{3})}{\gamma(\omega_{1})\gamma(\omega_{3})}\left(\frac{1-2p}{1-p}\right)^{n}\left(1-p\right)^{m}&\textrm{if
}n<m\\\
\frac{\widetilde{\omega}(\omega_{1})^{2}}{\gamma(\omega_{1})^{2}}(1-2p)^{m}&\textrm{if
}n=m\\\
\frac{\widetilde{\omega}(\omega_{1})\widetilde{\omega}(\omega_{3})}{\gamma(\omega_{1})\gamma(\omega_{3})}\left(\frac{1-2p}{1-p}\right)^{m}\left(1-p\right)^{n}&\textrm{if
}n>m.\end{cases}$ (75)

where we defined a new function of a a walk,

$\gamma(\omega)\coloneqq\prod_{i\in\omega}\sum_{k\sim i}w_{ik}.$ (76)

$\gamma$ computes the sum of edge weights connected to each node in the walk
$\omega$ (excluding the last), then takes the product of these quantities. It
is straightforward to check that, when all the graph weights are equal,
$\frac{\widetilde{\omega}(\omega)}{\gamma(\omega)}=\frac{1}{d^{m}}$ with $m$
the length of $\omega$. Meanwhile, $p(\omega_{1})$ becomes

$p(\omega_{1})=\frac{(1-p)^{m}\widetilde{\omega}(\omega_{1})}{\gamma(\omega_{1})}$
(77)

such that these modifications cancel out when we evaluate Eq. 36.

Case $3$ (Sec. 49) is also straightforward. The prefactor in square brackets
is equal to 72 and is again positive for any valid sampling strategy
$p(\omega_{1}\in\bar{\Omega}(k_{1},i))$ and $\mathbf{J}$ does not change, so
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.
$\widetilde{\omega}/\gamma$ can be replaced by some more complicated variant
that defines a valid probability distribution
$p(\omega_{1}\in\bar{\Omega}(k_{1},i))$ and antithetic termination will still
prove effective.

####  8.4.4 Summary

In Sec. 8.4, our theoretical results for antithetic termination have proved
robust to generalisations such as 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

$\frac{p(\omega_{1},\omega_{3})}{p(\omega_{1})p(\omega_{3})}$ (78)

is 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_ are 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.3.3 which require a little more
work, but the mathematics conspires that our arguments are again essentially
unmodified, though perhaps without such an intuitive explanation.

###  8.5 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$, $\Delta<p$. Unlike antithetic coupling, it
permits simultaneous termination. Eqs 13 become

$\displaystyle p(s_{1})=p(s_{2})=p,\hskip
14.22636ptp(\bar{s}_{1})=p(\bar{s}_{2})=1-p,\hskip
14.22636ptp(s_{2}|s_{1})=\frac{p-\Delta}{p},$ (79) $\displaystyle
p(\bar{s}_{2}|s_{1})=\frac{\Delta}{p},\hskip
14.22636ptp(s_{2}|\bar{s}_{1})=\frac{\Delta}{1-p},\hskip
14.22636ptp(\bar{s}_{2}|\bar{s}_{1})=\frac{1-p-\Delta}{1-p}.$

The probability of two antithetic walks $\bar{\Omega}(1,i)$ and
$\bar{\Omega}(3,i)$ containing subwalks $\omega_{1}$ and $\omega_{3}$ becomes

$p(\omega_{3}\in\bar{\Omega}(3,i),\omega_{1}\in\bar{\Omega}(1,i))=\begin{cases}\frac{1}{d^{m+n}}\left(\frac{1-p-\Delta}{1-p}\right)^{n}\left(1-p\right)^{m}&\textrm{if
}n<m\\\ \frac{1}{d^{2m}}(1-p-\Delta)^{m}&\textrm{if }n=m\\\
\frac{1}{d^{m+n}}\left(\frac{1-p-\Delta}{1-p}\right)^{m}\left(1-p\right)^{n}&\textrm{if
}n>m,\end{cases}$ (80)

which the reader might compare to Eq. 40. In analogy to Eq. 46, this induces
the matrix

$D^{\Delta}_{k_{1},k_{3}}\coloneqq\frac{w^{2}\lambda_{k_{1}}\lambda_{k_{3}}\frac{1-p-\Delta}{(1-p)^{2}}}{1-w^{2}\lambda_{k_{1}}\lambda_{k_{3}}\frac{1-p-\Delta}{(1-p)^{2}}}\left(\frac{1-w^{2}\lambda_{k_{1}}\lambda_{k_{3}}}{(1-w\lambda_{k_{1}})(1-w\lambda_{k_{3}})}\right).$
(81)

We can immediately observe that this is exactly equal to $C_{k_{1},k_{3}}$
when $\Delta=p(1-p)$, so for a pair of walkers with this TRV offset the
variance will be identical to the i.i.d. result. Replacing $D$ by $D^{\Delta}$
in $E_{pq}$ and $F_{pq}$ and $J_{pq}$ and reasoning about negative
definiteness via their respective Taylor expansions (as well as the new
possible cross-term $D_{k_{1},k_{3}}D^{\Delta}_{k_{1},k_{3}}$), it is
straightforward conclude that variance is suppressed compared to the i.i.d.
case provided $\Delta>p(1-p)$ and $\rho(\mathbf{U})$ or $p$ is sufficiently
small. The $p\to 0$ limit demands a slightly more careful treatment: in order
to stay in the regime $p(1-p)<\Delta<p$ we need to simultaneously take
$\Delta\to 0$, e.g. by defining $\Delta(p)\coloneqq p(1-p)+ap^{2}$ with the
constant $0<a<1$. ∎

This result was reported in Theorem 3.4 of the main text.

###  8.6 What about diagonal terms?

The alert reader might remark that all derivations in Sec. 8.3 have taken
$i\neq j$, considering estimators of the off-diagonal elements of the matrix
$(\mathbf{I}_{N}-\mathbf{U})^{-2}$. In fact, estimators of the diagonal
elements $\phi(i)^{\top}\phi(i)$ will be biased for both GRFs and q-GRFs if
$\phi(i)$ is constructed using the same ensemble of walkers because each
walker is manifestly correlated with, rather than independent of, itself. This
is rectified by taking _two_ ensembles of walkers out of each node, each of
which may exhibit antithetic correlations among itself, then taking the
estimator $\phi_{1}(i)^{\top}\phi_{2}(i)$. It is straightforward to convince
oneself that, in this setup, the estimator is unbiased and q-GRFs will
outperform GRFs. In practice, this technicality has essentially no effect on
(q-)GRF performance and doubles runtime so we omit further discussion.

###  8.7 Further experimental details: compute, datasets and uncertainties

The experiments in Secs. 4.1, 4.2 and 4.4 were carried out on an Intel® Core™
i5-7640X CPU @ 4.00GHz × 4. Each required $\sim 1$ CPU hour. The experiments
in Sec. 4.3 were carried out on a 2-core Xeon 2.2GHz with 13GB RAM and 33GB
HDD. The computations for the largest considered graphs took $\sim 1$ CPU
hour.

The real-world graphs and meshes were accessed from the repositories
[Ivashkin, 2023] and [Dawson-Haggerty, 2023], with further information about
the datasets available therein. Where we were able to locate them, the
original papers presenting the graphs are: [Zachary, 1977, Lusseau et al.,
2003, Newman, 2006, Bollacker et al., 1998, Leskovec et al., 2007].

All our experiments report standard deviations on the means, apart from the
clustering task in Sec. 4.3 because running kernelised $k$-means on large
graphs is expensive.