Graph Random Features (GRFs)

• Graph Random Features (GRFs) are randomized embeddings that approximate graph kernels using Monte Carlo random walks and sampling techniques.
• The approach leverages random-walk estimators and convolution identities to replace cubic kernel evaluations with subquadratic, scalable computation.
• GRFs enable diverse applications from node-level predictions to graph-level embeddings, offering efficient performance with controlled variance.

Graph Random Features (GRFs) are explicit randomized embeddings that replace expensive exact graph-kernel computations by inner products of sampled feature vectors. In the literature, the term has been used for several related constructions: random-walk estimators for graph-node kernels such as regularized Laplacian, diffusion, and other power-series kernels; graphlet-sampling methods that replace graph-isomorphism tests by random features, including optical random features; global-alignment graph embeddings obtained from random graphs; and closely related random-neural-feature embeddings built from permutation-invariant graph neural networks (Choromanski, 2023, Ghanem et al., 2020, Wu et al., 2019, Zambon et al., 2019).

1. Problem setting and main formulations

A central GRF setting starts from an undirected weighted graph $G=(V,E,W)$ with $|V|=N$ and a graph-node kernel written as a power series

$K_{\alpha}(W)=\sum_{k=0}^{\infty}\alpha_k W^k.$

Under mild spectral-radius conditions this converges and defines a positive-definite kernel on $V\times V$. The computational bottleneck is that exact evaluation of dense graph kernels typically requires $O(N^3)$ time through eigendecomposition, matrix products, or linear solves. GRFs replace those exact computations by Monte Carlo feature maps $\phi(v)$ such that

$$\mathbb E[\phi(u)^\top \phi(v)] = K_{\alpha}(W)[u,v],$$

or, in a transition-matrix formulation with $P=D^{-1}A$,

$$\mathbb E[\phi_{f_1}(i)^\top \phi_{f_2}(j)] = \Bigl[\sum_{k=0}^{\infty}\alpha_k P^k\Bigr]_{ij}.$$

This establishes GRFs as graph analogues of classical random features, but with random walks and graph-dependent load deposition replacing Fourier sampling (Choromanski, 2023, Reid et al., 2023).

The term also appears in graph-level settings, where the target object is not a node-kernel matrix but an explicit embedding of an entire graph. In those cases the randomization may act on sampled graphlets, on random reference graphs, or on random GNN parameters. The common denominator is the same: a high- or infinite-dimensional graph similarity is replaced by an explicit Monte Carlo embedding whose inner products approximate the desired kernel or metric.

Family	Random object	Output
Random-walk GRFs	random walks with modulation $f$	node-kernel estimator
Graphlet GRFs	sampled $|V|=N$0-node subgraphs	graph embedding $|V|=N$1
Random graph embeddings	random graphs $|V|=N$2	graph feature vector $|V|=N$3
Graph random neural features	random GNN parameters $|V|=N$4	embedding $|V|=N$5

A common misconception is that GRFs denote a single fixed algorithm. The literature instead presents a family of randomized graph embeddings sharing the same computational motif—explicit features in place of exact kernel evaluation—but differing substantially in the graph object being randomized and in the target similarity being approximated.

2. Core random-walk mechanism on graph nodes

The canonical GRF construction uses a modulation sequence $|V|=N$6 satisfying a convolution identity. In the basic two-fold case,

$|V|=N$7

If $|V|=N$8, then formally $|V|=N$9. GRFs estimate rows of $K_{\alpha}(W)=\sum_{k=0}^{\infty}\alpha_k W^k.$0 by tracing random walks. For a target node $K_{\alpha}(W)=\sum_{k=0}^{\infty}\alpha_k W^k.$1, one repeats $K_{\alpha}(W)=\sum_{k=0}^{\infty}\alpha_k W^k.$2 independent random walks, maintains a running load $K_{\alpha}(W)=\sum_{k=0}^{\infty}\alpha_k W^k.$3, deposits $K_{\alpha}(W)=\sum_{k=0}^{\infty}\alpha_k W^k.$4 at the current node at step $K_{\alpha}(W)=\sum_{k=0}^{\infty}\alpha_k W^k.$5, updates the load by the importance-corrected edge weight, and halts according to a prescribed length mechanism. After averaging over walks, one obtains a sparse feature vector $K_{\alpha}(W)=\sum_{k=0}^{\infty}\alpha_k W^k.$6 whose inner products unbiasedly estimate the target kernel (Choromanski, 2023, Choromanski et al., 9 Oct 2025).

A transition-matrix formulation makes the same idea explicit for arbitrary functions of a weighted adjacency matrix. Let $K_{\alpha}(W)=\sum_{k=0}^{\infty}\alpha_k W^k.$7 and suppose the desired kernel is

$K_{\alpha}(W)=\sum_{k=0}^{\infty}\alpha_k W^k.$8

Choose modulators $K_{\alpha}(W)=\sum_{k=0}^{\infty}\alpha_k W^k.$9 satisfying

$V\times V$0

For each node $V\times V$1, one performs $V\times V$2 independent random walks with length $V\times V$3, deposits load $V\times V$4 along visited nodes, and averages the resulting vectors. The unbiasedness identity

$V\times V$5

follows from the discrete convolution constraint. This framework covers diffusion, regularized Laplacian, and $V\times V$6-step random-walk kernels (Reid et al., 2023).

Several papers modify the walk-length mechanism without changing unbiasedness. Instead of Bernoulli halting with fixed $V\times V$7, one may pre-sample a walk length $V\times V$8, define $V\times V$9, add $O(N^3)$0 at step $O(N^3)$1, and stop exactly when $O(N^3)$2. This generalized termination remains unbiased if the $O(N^3)$3 correction is used. The same line of work also introduces walk stitching: rather than sampling one long walk, GRFs++ form products of estimators built from several independent sets of short walks, using a $O(N^3)$4-fold convolution constraint on $O(N^3)$5 to preserve unbiasedness (Choromanski et al., 9 Oct 2025).

3. Approximation guarantees, sparsity, and variance control

The original motivation for GRFs is algorithmic: exact graph-kernel evaluation is cubic in the number of nodes, whereas random-walk estimators can be computed from local graph access. For u-GRFs, each walk has expected length $O(N^3)$6, so with $O(N^3)$7 walks per node and $O(N^3)$8 nodes the total time is $O(N^3)$9; if $\phi(v)$0 is constant and $\phi(v)$1, this is subquadratic in $\phi(v)$2, and with high probability one can show $\phi(v)$3. Parallelism is immediate because walks from different nodes and different walkers are independent (Reid et al., 2023).

Theoretical control of the Monte Carlo estimator takes several forms. For general u-GRFs, the kernel estimate is unbiased and standard concentration yields error $\phi(v)$4 with high probability. For topological masking in linear transformers, the first exponential concentration bound for GRFs states that if

$\phi(v)$5

then for all $\phi(v)$6 and any $\phi(v)$7,

$\phi(v)$8

The same work proves a sparsity lemma: with $\phi(v)$9 walks of Geometric$$\mathbb E[\phi(u)^\top \phi(v)] = K_{\alpha}(W)[u,v],$$0 length, the number of nonzero entries in $$\mathbb E[\phi(u)^\top \phi(v)] = K_{\alpha}(W)[u,v],$$1 is at most $$\mathbb E[\phi(u)^\top \phi(v)] = K_{\alpha}(W)[u,v],$$2 with probability at least $$\mathbb E[\phi(u)^\top \phi(v)] = K_{\alpha}(W)[u,v],$$3, where

$$\mathbb E[\phi(u)^\top \phi(v)] = K_{\alpha}(W)[u,v],$$4

These results formalize why GRFs can be used inside $$\mathbb E[\phi(u)^\top \phi(v)] = K_{\alpha}(W)[u,v],$$5 linear-attention pipelines (Reid et al., 2024).

Variance reduction has become a distinct GRF subtopic. Quasi-Monte Carlo Graph Random Features (q-GRFs) couple walkers by antithetic termination, preserving the marginals of the walk-length distribution while inducing negative correlations between paired walk lengths. For the $$\mathbb E[\phi(u)^\top \phi(v)] = K_{\alpha}(W)[u,v],$$6-regularised Laplacian kernel, this yields lower-variance estimators under mild conditions, and the method is described as the first rigorously studied quasi-Monte Carlo scheme for kernels defined on combinatorial objects (Reid et al., 2023). Follow-up work on the diffusion kernel reports that q-GRFs achieve lower variance estimators on Ladder graphs, with the number of rungs affecting performance; by contrast, binary trees did not benefit in those experiments (Feinberg et al., 2024).

GRFs++ address a different limitation: long-range interactions. Standard GRFs rely on two independent walks intersecting, and the probability of capturing long-distance similarity decays exponentially in distance. Walk stitching replaces one long sequential walk by several shorter walks combined through sparse matrix products, and Theorem 3.4 states that under the standard Bernoulli-halting scheme the mean-squared error decreases monotonically for degree $$\mathbb E[\phi(u)^\top \phi(v)] = K_{\alpha}(W)[u,v],$$7 (Choromanski et al., 9 Oct 2025).

4. Graph-level GRFs: graphlets, random graphs, and random GNNs

A separate GRF line targets entire graphs rather than nodes. In graphlet kernels, the classical construction enumerates induced $$\mathbb E[\phi(u)^\top \phi(v)] = K_{\alpha}(W)[u,v],$$8-node subgraphs and compares them through the one-hot isomorphism map

$$\mathbb E[\phi(u)^\top \phi(v)] = K_{\alpha}(W)[u,v],$$9

The normalized histogram

$P=D^{-1}A$0

defines the graphlet kernel $P=D^{-1}A$1, but the cost

$P=D^{-1}A$2

is intractable for moderate $P=D^{-1}A$3. The GRF replacement samples subgraphs $P=D^{-1}A$4, computes

$P=D^{-1}A$5

and links $P=D^{-1}A$6 to the Maximum Mean Discrepancy between subgraph distributions. With random features on the base subgraph kernel, the deviation from the MMD is $P=D^{-1}A$7. Optical random features further replace the usual $P=D^{-1}A$8 random projection by an Optical Processing Unit readout

$P=D^{-1}A$9

with cost $$\mathbb E[\phi_{f_1}(i)^\top \phi_{f_2}(j)] = \Bigl[\sum_{k=0}^{\infty}\alpha_k P^k\Bigr]_{ij}.$$0 in the compute-heavy part (Ghanem et al., 2020).

Global-alignment graph kernels provide another graph-level formulation. Here each graph is first mapped to geometric node embeddings using the $$\mathbb E[\phi_{f_1}(i)^\top \phi_{f_2}(j)] = \Bigl[\sum_{k=0}^{\infty}\alpha_k P^k\Bigr]_{ij}.$$1 smallest eigenvectors of the normalized Laplacian, together with normalized degree weights $$\mathbb E[\phi_{f_1}(i)^\top \phi_{f_2}(j)] = \Bigl[\sum_{k=0}^{\infty}\alpha_k P^k\Bigr]_{ij}.$$2. Similarity between graphs is measured through Earth Mover’s Distance, and the kernel is defined by integrating

$$\mathbb E[\phi_{f_1}(i)^\top \phi_{f_2}(j)] = \Bigl[\sum_{k=0}^{\infty}\alpha_k P^k\Bigr]_{ij}.$$3

over a distribution $$\mathbb E[\phi_{f_1}(i)^\top \phi_{f_2}(j)] = \Bigl[\sum_{k=0}^{\infty}\alpha_k P^k\Bigr]_{ij}.$$4 on random graphs. Monte Carlo sampling of $$\mathbb E[\phi_{f_1}(i)^\top \phi_{f_2}(j)] = \Bigl[\sum_{k=0}^{\infty}\alpha_k P^k\Bigr]_{ij}.$$5 random graphs yields the explicit feature vector

$$\mathbb E[\phi_{f_1}(i)^\top \phi_{f_2}(j)] = \Bigl[\sum_{k=0}^{\infty}\alpha_k P^k\Bigr]_{ij}.$$6

with a uniform approximation bound from Theorem 3.1 and $$\mathbb E[\phi_{f_1}(i)^\top \phi_{f_2}(j)] = \Bigl[\sum_{k=0}^{\infty}\alpha_k P^k\Bigr]_{ij}.$$7quasi-$$\mathbb E[\phi_{f_1}(i)^\top \phi_{f_2}(j)] = \Bigl[\sum_{k=0}^{\infty}\alpha_k P^k\Bigr]_{ij}.$$8linear scalability with respect to the number and the size of the graphs (Wu et al., 2019).

A related but distinct construction is Graph Random Neural Features (GRNF). One fixes a family $$\mathbb E[\phi_{f_1}(i)^\top \phi_{f_2}(j)] = \Bigl[\sum_{k=0}^{\infty}\alpha_k P^k\Bigr]_{ij}.$$9 of node-permutation-invariant one-hidden-layer GNNs $f$0, places a probability distribution $f$1 on the parameter space $f$2, samples $f$3, and defines

$f$4

The induced kernel $f$5 is complete, the induced distance $f$6 is a metric, and for any two fixed graphs $f$7 one has

$f$8

GRNF can be used as a standalone embedding or as a training-free input layer of a GNN (Zambon et al., 2019).

5. Applications and systems integration

GRFs have been integrated into several large-scale learning systems because they expose graph kernels as explicit sparse features rather than dense $f$9 matrices. In scalable Gaussian processes on graphs, GRFs produce a sparse feature matrix $|V|=N$00 with $|V|=N$01 nonzeros when $|V|=N$02 walkers per node are used. This allows $|V|=N$03-time matrix-vector products with $|V|=N$04, and under mild assumptions conjugate gradients yields $|V|=N$05 time for Bayesian inference, versus $|V|=N$06 for exact kernels. The same framework reports $|V|=N$07 memory and Bayesian optimisation on graphs with over $|V|=N$08 nodes on a single computer chip (Zhang et al., 3 Sep 2025).

In graph-aware transformers, GRFs are used to approximate topological masks $|V|=N$09 while preserving linear attention. The construction tensorizes ordinary query/key random features with sparse GRF topological features, leading to an $|V|=N$10 masking algorithm for arbitrary weighted graphs. Reported experiments cover image graphs, point-cloud graphs with $|V|=N$11k nodes, and a preliminary video setting; the stated comparison is that the fastest previous alternative was $|V|=N$12 and only suitable for specific graphs (Reid et al., 2024).

GRFs have also become building blocks in continuous and implicit settings. Manifold Random Features discretize a manifold, compute GRF signatures on the induced graph, and train a bounded nonnegative field $|V|=N$13 so that

$|V|=N$14

On refined grids, the resulting rescaled GRF signatures converge to Gaussian random features, and the paper states a deep asymptotic connection between GRFs and continuous random features used for regular kernels. Empirically, the method reports $|V|=N$15–$|V|=N$16 speed-up against spectral baselines on synthetic manifolds (Parashar et al., 3 Feb 2026).

For implicitly defined graphs, SWING replaces discrete walks on materialized adjacency matrices by walks in the continuous feature space in which the graph is embedded. Edge weights are linearized with random features, Gumbel-softmax sampling is factorized, and the resulting system scales as $|V|=N$17 time and $|V|=N$18 memory, without input graph materialization (Manenti et al., 13 Feb 2026).

6. Limitations, misconceptions, and active directions

The GRF literature repeatedly emphasizes that the method’s advantages are inseparable from Monte Carlo trade-offs. Variance decreases as the number of sampled walks or features increases, but sufficiently large $|V|=N$19, $|V|=N$20, or $|V|=N$21 may still be required for accurate estimates; this is explicit in the $|V|=N$22 or $|V|=N$23 approximation rates and in the dimension rule $|V|=N$24 for GRNF (Reid et al., 2023, Zambon et al., 2019). Standard GRFs also struggle with long-range relationships because intersection probabilities between independent walks decay exponentially in distance, which motivated walk stitching and generalized length distributions in GRFs++ (Choromanski et al., 9 Oct 2025).

Several implementation-specific limitations are likewise documented. In transformer masking, the symmetric construction has a high ambient feature dimension $|V|=N$25 and benefits from sparse linear-algebra support; asymmetric versions incur larger Monte Carlo variance (Reid et al., 2024). In GRF-based Gaussian processes, shared-walk randomness induces $|V|=N$26 bias in diagonal entries unless separate ensembles are used, and preconditioning of $|V|=N$27 remains an open opportunity (Zhang et al., 3 Sep 2025). For very large graphs, the cost of a single forward pass in graph-random-neural-feature models may still be expensive, even though the embedding is linear in $|V|=N$28 (Zambon et al., 2019).

Research directions in the cited papers are diverse but coherent. They include learned modulation functions $|V|=N$29 for variance reduction or end-to-end kernel learning, dynamic graphs, higher-order random walks on edges or subgraphs, Johnson–Lindenstrauss compression of graph features, quasi-Monte Carlo couplings beyond antithetic termination, and continuous-space GRFs for implicit graphs and manifolds (Reid et al., 2023, Reid et al., 2024, Manenti et al., 13 Feb 2026). Taken together, these developments suggest that GRFs have evolved from a specialized Monte Carlo estimator for regularized Laplacian kernels into a broader design pattern for scalable kernelization on discrete, graph-structured, and even discretized continuous domains.