Paper deep dive
Random Wavelet Features for Graph Kernel Machines
Valentin de Bassompierre, Jean-Charles Delvenne, Laurent Jacques
Intelligence
Status: succeeded | Model: Gemma-4-26B-A4B | Prompt: intel-v1 | Confidence: 93%
Last extracted: 7/21/2026, 2:57:06 AM
Summary
The paper introduces Random Wavelet Features for Graph Kernel Machines, a method to approximate graph kernels using randomized spectral node embeddings. By leveraging the Graph Wavelet Transform on random signals, the authors construct embeddings whose dot products estimate a low-rank approximation of Laplacian-based graph kernels. This approach offers a scalable alternative to direct kernel computation, particularly effective for spectrally localized kernels, and connects to Randomized Singular Value Decomposition (RSVD).
Entities (10)
Relation Signals (7)
Laurent Jacques → affiliatedwith → UCLouvain
confidence 95% · Valentin de Bassompierre, Jean-Charles Delvenne, and Laurent Jacques INMA/ICTEAM, UCLouvain, Belgium
Valentin de Bassompierre → affiliatedwith → UCLouvain
confidence 95% · Valentin de Bassompierre, Jean-Charles Delvenne, and Laurent Jacques INMA/ICTEAM, UCLouvain, Belgium
Jean-Charles Delvenne → affiliatedwith → UCLouvain
confidence 95% · Valentin de Bassompierre, Jean-Charles Delvenne, and Laurent Jacques INMA/ICTEAM, UCLouvain, Belgium
Random Wavelet Features → approximates → Graph Kernel
confidence 95% · we introduce randomized spectral node embeddings whose dot products estimate a low-rank approximation of any specific graph kernel.
Random Wavelet Features → uses → Graph Wavelet Transform
confidence 95% · Specifically, we leverage the graph wavelet transform [5] on random input signals to construct embeddings that provide scalable approximations of graph kernels.
Graph Wavelet Transform → basedon → Graph Laplacian
confidence 90% · GSP leverages the graph Laplacian, and the derived graph Fourier transform... Specifically, we leverage the graph wavelet transform
Random Wavelet Features → relatedto → Randomized Singular Value Decomposition
confidence 85% · Finally, we discuss connections to randomized singular value decomposition (RSVD) [6], highlighting links between our method and classical numerical linear algebra techniques.
Cypher Suggestions (0)
No Cypher suggestions yet.
Abstract
Abstract:Node embeddings map graph vertices into low-dimensional Euclidean spaces while preserving structural information. They are central to tasks such as node classification, link prediction, and signal reconstruction. A key goal is to design node embeddings whose dot products capture meaningful notions of node similarity induced by the graph. Graph kernels offer a principled way to define such similarities, but their direct computation is often prohibitive for large networks. Inspired by random feature methods for kernel approximation in Euclidean spaces, we introduce randomized spectral node embeddings whose dot products estimate a low-rank approximation of any specific graph kernel. We provide theoretical and empirical results showing that our embeddings achieve more accurate kernel approximations than existing methods, particularly for spectrally localized kernels. These results demonstrate the effectiveness of randomized spectral constructions for scalable and principled graph representation learning.
Tags
Links
- Source: https://arxiv.org/abs/2602.15711v1
- Canonical: https://arxiv.org/abs/2602.15711v1
Trouble viewing inline? Open PDF directly →
Full Text
56,818 characters extracted from source content.
Expand or collapse full text
Random Wavelet Features for Graph Kernel Machines (extended version with supplementary material) 111This work was supported by the Fonds de la Recherche Scientifique - FNRS under Grants n° 40038614 and FNRS T.0160.24. During this research, LJ benefited from an INRIA Chair at the Institute for Advanced Study, Collegium of Lyon, working in the OCKHAM team, at ENS Lyon. Valentin de Bassompierre, Jean-Charles Delvenne, and Laurent Jacques INMA/ICTEAM, UCLouvain, Belgium Abstract: Node embeddings map graph vertices into low-dimensional Euclidean spaces while preserving structural information. They are central to tasks such as node classification, link prediction, and signal reconstruction. A key goal is to design node embeddings whose dot products capture meaningful notions of node similarity induced by the graph. Graph kernels offer a principled way to define such similarities, but their direct computation is often prohibitive for large networks. Inspired by random feature methods for kernel approximation in Euclidean spaces, we introduce randomized spectral node embeddings whose dot products estimate a low-rank approximation of any specific graph kernel. We provide theoretical and empirical results showing that our embeddings achieve more accurate kernel approximations than existing methods, particularly for spectrally localized kernels. These results demonstrate the effectiveness of randomized spectral constructions for scalable and principled graph representation learning. 1 Introduction Kernel machines are a class of machine learning algorithms, designed to handle non-linear problems by implicitly mapping data into a high-dimensional feature space. Given some data space X, this is achieved through a positive semi-definite (psd) kernel Γ:(,)∈×↦Γ(,) :( x, y) ×X ( x, y). This kernel can be interpreted, through the kernel trick, as a dot product in a (possibly infinite-dimensional) feature (Hilbert) space ℋH, i.e., there exist a mapping ϕ:→ℋφ:X such that Γ(,)=⟨ϕ(),ϕ()⟩ℋ ( x, y)= φ( x),φ( y) _H. In this context, algorithms perform linear computations in this enriched space without explicitly computing the feature map ϕφ, thereby enabling more expressive linear models. While kernels are often introduced for data in continuous domains such as =ℝdX=R^d, many practical problems involve data supported on graphs, where the underlying relational structure between data samples plays a central role. It is therefore essential to define kernels that respect and exploit this structure. As such, graph kernels [1, 2] provide a way to quantify similarity between pairs of nodes. For a graph G with N nodes, such a kernel can be represented as an N×N× N kernel matrix—most often, a function of the graph Laplacian. In practice, however, computing this matrix is often prohibitive, with a typical cost scaling as (N3)O(N^3), which limits the applicability of many existing graph kernel methods. In ℝdR^d, kernel methods scale poorly on large datasets as computing and storing full kernel matrices is costly. To address this, Rahimi and Recht [3] proposed random (Fourier) features. These explicitly map data into a low-dimensional space via a randomized feature map ϕ:ℝd→ℝD φ:R^d ^D so that, for specific kernels Γ such as the Gaussian or Laplacian kernels, Γ(,)=⟨ϕ(),ϕ()⟩≈⟨ϕ(),ϕ()⟩ ( x, y)=E φ( x), φ( y) ≈ φ( x), φ( y) for all ,∈ℝd x, y ^d. This allows fast linear learning on the transformed points, reducing training and evaluation time as well as memory usage. It is thus natural to ask whether graph kernel computations can be similarly accelerated using random features. We propose to use tools from Graph Signal Processing (GSP) [4] for randomized graph kernel approximation. GSP leverages the graph Laplacian, and the derived graph Fourier transform, to extend the definition of frequency, filtering, and smoothness for graph signals while respecting the graph structure. Specifically, we leverage the graph wavelet transform [5] on random input signals to construct embeddings that provide scalable approximations of graph kernels. Contributions: Building on random feature methods for kernel approximation in Euclidean spaces, we propose randomized spectral embeddings for graphs, where the dot product between embeddings estimates a low-rank approximation of a Laplacian-based graph kernel. We provide both theoretical and empirical evidence that these embeddings yield more accurate kernel approximations than existing approaches, especially when the kernel is spectrally localized. Finally, we discuss connections to randomized singular value decomposition (RSVD) [6], highlighting links between our method and classical numerical linear algebra techniques. 2 Related work In [7], the graph wavelet transform of random signals is used to accelerate the computation of graph matrix functions. Their use case is slightly different from ours, as they approximate correlation matrices based explicitly on wavelets instead of general graph kernels. In [8], the application of graph wavelet transforms to random signals yields embeddings that approximate the distance matrix of classical spectral clustering. We somehow extend these results to general graph kernels. Other works have considered random feature methods for graph kernel approximation, see e.g., [9, 10]. Their random features follow a spatial design relying on random walks to assess the similarities between nodes, and produce unbiased kernel approximations. We show that the convergence to said kernel is slow for bandlimited kernels which are well localized in the spectral domain and spread out in the spatial domain. 3 Mathematical tools 3.A. Spectral graph theory: We consider an undirected weighted graph =(,,w)G=( V, E,w) made of N vertices (or nodes) V identified, up to a correct ordering, with the N=||N=| V| integers [N]:=1,…,N[N]:=\1,…,N\. These nodes are connected with edges in ⊆× E V× V, with a total of E=||E=| E| edges, and w:→ℝ+w: E ^+ is a weight function assigning a positive weight to each edge. Any vector ∈ℝN f ^N can be seen as a function →ℝ V defined on each node of V. The weight matrix ∈ℝN×N W ^N× N of a graph has entries WijW_ij equal to w(i,j)w(i,j) if (i,j)∈(i,j)∈ E, and to 0 otherwise. It is symmetric for undirected graphs. The degree matrix ∈ℝN×N D ^N× N is diagonal with entries Dii≔∑j=1NWijD_i _j=1^NW_ij equal to the degrees of the nodes. Finally, let the normalized graph Laplacian L be ≔−1/2−1/2 L I- D^-1/2 W D^-1/2. The graph Laplacian L plays a central role in GSP and spectral graph theory. Since L is symmetric, it factorizes as =⊤ L= V V , where =(1,…,N) V=( v_1,…, v_N) is an orthogonal matrix whose columns are the eigenvectors k v_k of eigenvalue λk∈[0,2] _k∈[0,2], i.e., k=λkk L v_k= _k v_k, and the diagonal matrix =diag(λ1,…,λN) =diag( _1,…, _N) represents the spectrum of L [11]. The eigenvectors k∈ℝN v_k ^N define a graph spectral domain; as functions over G, they can be interpreted as graph harmonics, discrete analogues of Fourier modes. The eigenvalues act as squared graph frequencies: small (resp. large) λk _k correspond to slowly (resp. rapidly) varying modes on the graph. For any signal ∈ℝN f ^N defined over V, the graph Fourier transform (GFT) projects f onto the eigenvectors of the graph Laplacian, i.e., its forward and inverse operators are ^=ℱ[]:=⊤,=ℱ−1[^]:=^. f 2.79999pt -2.79999pt=F[ f]:= V f, f=F^-1[ f 2.79999pt -2.79999pt]:= V f 2.79999pt -2.79999pt. The representation f 2.79999pt -2.79999pt expresses f in the graph spectral domain. Given a filter g:ℝ→ℝg:R defined in the spectral domain, the filtering g f^g of f by g is defined in the spectral domain as ℱ[g]k:=g(λk)f^k,k∈[N],orℱ[g]=g()^,F [ f^g ]_k:=g( _k) f_k,\ k∈[N],\ or\ F [ f^g ]=g( ) f 2.79999pt -2.79999pt, with g()=diag(g(λ1),…,g(λN))g( )=diag(g( _1),…,g( _N)). The inverse GFT yields g=g()^=g(),withg()=g()⊤. f^g= Vg( ) f 2.79999pt -2.79999pt=g( L) f,\ with\ g( L)= Vg( ) V . (1) Therefore, filtering f amounts to multiplying it by the function of the Laplacian g()g( L). The graph wavelet transform (GWT) implements spectral filtering with specific filters g providing localization in both the frequency and spatial domains. While computing g()g( L) exactly requires a full eigendecomposition (ED) of L, which is prohibitive for large graphs, GWT approximates g by a Chebyshev polynomial. If the Laplacian is sparse, i.e., if the number of edges E=(N)E=O(N), the GWT can be efficiently obtained through a series of fast matrix-vector products with powers of the Laplacian matrix [5]. 3.B. Graph kernels: A graph kernel ∈ℝN×N ^N× N is defined through a psd similarity function Γ:×→ℝ : V× V between nodes. Most graph kernels are defined as a function =h() =h( L) of the graph Laplacian L [1]. Common examples include: =(+σ2)−d, =( I+σ^2 L)^-d, (d-regularized lap.) =(−/σ)2,withσ⩾2, =( I- L/σ)^2,\ with\ σ 2, (22-step random walk) =exp(−σ2), = (-σ^2 L), (diffusion process) =cos(σ2π/4),with|σ|⩽1. = (σ^2 Lπ/4),\ with\ |σ| 1. (inverse cosine) If =h() =h( L), the spectral theorem provides =∑k=1Nh(λk)kk⊤=h()⊤, = _k=1^Nh( _k) v_k v_k = Vh( ) V , (2) i.e., h(λk)h( _k) weights the contribution of each Laplacian eigenvector. The function h thus acts as a spectral filter that emphasizes or suppresses specific frequency components. Since kernels act as smoothing operators, they should preserve low-frequency, smooth components while attenuating high-frequency, rapidly varying ones. Accordingly, h(λ)h(λ) should decrease with λ [1]. For instance, spectral clustering uses a low-pass filter that retains only the smoothest eigenmodes [12]. 4 Proposed Method Our approach provides a random feature map ϕ:i∈↦ϕi=i∈ℝK φ:i∈ V φ_i= δ_i ^K, with a K×NK× N matrix with columns ϕii=1N\ φ_i\_i=1^N, and i∈ℝN δ_i ^N equal to 1 on the i-th node and 0 elsewhere, such that the dot-products ⟨ϕi,ϕj⟩=Γ~ij φ_i, φ_j = _ij define a rank-K approximation ~ of the graph kernel :=h() :=h( L). We assume that this kernel is associated with some known positive, decreasing function h:[0,2]→ℝ+h:[0,2] _+, with h(0)=1h(0)=1. As a reference, the optimal rank-K approximation of is given by its truncated spectral decomposition, (K)=∑k=1Kh(λk)kk⊤=:K:K⊤, ^(K)= _k=1^Kh( _k)\, v_k v_k = V_:K V_:K , (3) which can also be viewed as its projection onto the subspace spanned by :K:=(1,…,K) V_:K:=( v_1,…, v_K) (i.e., the K smoothest graph harmonics). We show that ~ attains an approximation quality comparable to this optimum, ‖−~‖≈‖−(K)‖=h(λK+1)\| - \|≈\| - ^(K)\|=h( _K+1), while bypassing the explicit ED of the Laplacian. The algorithm building follows a two-step procedure: we first estimate span:Kspan V_:K, then is constructed from this estimated subspace and the known filter h (see Alg. 1). 4.A. [Part 1] Range finding: We first estimate span:Kspan V_:K, i.e., we find an orthogonal Q such that span≈span:Kspan Q V_:K. We proceed by filtering K random signals. For simplicity, we take Gaussian random signals ∈ℝN×K G ^N× K (i.e., Gij∼i.i.d.(0,1)G_ij _i.i.d.N(0,1)), but one may also consider using structured (Fourier or Hadamard-based) random matrices as proposed in [6, 13]. From (1), filtering a signal ∈ℝN g ^N with an ideal low-pass filter χK:=[0,λK] 1.42271pt$ χ$ 1.42271pt$ χ$ 1.42271pt$ χ$ 1.42271pt$ χ$_K:= 1_[0, _K] yields a vector =χK() b= 1.42271pt$ χ$ 1.42271pt$ χ$ 1.42271pt$ χ$ 1.42271pt$ χ$_K( L) g lying exactly in the relevant subspace, i.e., ∈span(:K) b ( V_:K). Further, with probability 1, filtering each column of G with χK 1.42271pt$ χ$ 1.42271pt$ χ$ 1.42271pt$ χ$ 1.42271pt$ χ$_K yields a basis =χK() B= 1.42271pt$ χ$ 1.42271pt$ χ$ 1.42271pt$ χ$ 1.42271pt$ χ$_K( L) G for that subspace. The (Gram-Schmidt) orthogonalization ∗=ortho() Q^*=ortho( B) thus provides an exact solution to our problem. In practice, constructing the ideal low-pass filter is infeasible without computing the first K eigenvectors of L (which would defeat the purpose of the range finding). We thus replace the ideal filter with a polynomial approximation222In this paper, we consider pχp_χ to be the Jackson-Chebyshev polynomial approximation to χK 1.42271pt$ χ$ 1.42271pt$ χ$ 1.42271pt$ χ$ 1.42271pt$ χ$_K [14]. This polynomial has a less sharp cut-off than Chebyshev polynomials but doesn’t present Gibbs oscillation phenomenon. pχp_χ, and obtain =ortho(pχ()). Q=ortho(p_χ( L) G). (4) Replacing χK 1.42271pt$ χ$ 1.42271pt$ χ$ 1.42271pt$ χ$ 1.42271pt$ χ$_K with pχp_χ means the filtered vectors do not lie exactly in span:Kspan V_:K. However, if the polynomial approximation is sufficiently accurate, we expect the filtered vectors to lie very close to that subspace. We characterize in Sec. 4.F the effect this has on the approximation error ‖−~‖\| - \|. Notice that to construct pχp_χ one needs access to the value of λK _K. We use [15, Alg. 1], which finds an estimate λ~K λ_K by dichotomy on the interval [0,2][0,2]—estimating the eigenvalue count in [0,λ(i)][0,λ^(i)] at each iteration i. Interestingly, the eigenvalue count estimation is performed by filtering (a very small number of) random Gaussian signals with the same polynomial low-pass filters pχi≈[0,λ(i)]p_ _i≈ 1_[0,λ^(i)] as used in our scheme. 4.B. [Part 2] Kernel approximation: The matrix Q obtained in Sec. 4.A can be viewed alternatively as an orthonormal basis for a subspace close to that spanned by the K foremost eigenvectors of L, or simply as a set of K signals living on the graph. With this second view, it is natural to think that we can filter these signals to yield relevant embeddings for our graph. In particular, defining the linear node embedding :i∈↦ϕi=i∈ℝK :i∈ V φ_i= δ_i ^K with =(ϕ1,…,ϕN)=(h12())⊤, =( φ_1,…, φ_N)=(h 12( L) Q) , (5) then the collection of dot-products (Γ~ij=⟨ϕi,ϕj⟩)i,j=1N( _ij= φ_i, φ_j )_i,j=1^N of the embedding of all node pairs yields a rank-K approximation ~ of the kernel :=h() :=h( L). To show this, let us first consider the ideal case where span=span:Kspan Q=span V_:K, i.e., where ⊤=:K:K⊤ Q Q = V_:K V_:K can be interpreted as the projection onto span:Kspan V_:K, and −⊤ I- Q Q as the projection on its orthogonal complement, the span of K:=(K+1,…,N) V_K:=( v_K+1,…, v_N). Then, Γ~ij _ij =⟨ϕi,ϕj⟩=i⊤h12()⊤h12()j = φ_i, φ_j = δ_i h 12( L) Q h 12( L) δ_j =i⊤h12()⊤:K:K⊤h12()⊤j = δ_i Vh 12( ) V V_:K V_:K Vh 12( ) V δ_j =i⊤:Kh(:K):K⊤j. = δ_i V_:Kh( _:K) V_:K δ_j. with :K=diag(λ1,…,λK) _:K=diag( _1,…, _K). From (2), we thus get ~=⊤=∑k=1Kh(λk)kk⊤. = = _k=1^Kh( _k)\, v_k v_k . The matrix ~ is thus the best rank-K approximation (K) ^(K) of =h()⊤ = Vh( ) V , or its k-truncated SVD version — see (3). Our approach only estimates this ideal case through (4) of Part I. In this case, in (5) yields a kernel ~ that is now an approximation to (K) ^(K): the projector ⊤ Q does not project exactly onto span:Kspan V_:K. As a consequence, some less important graph harmonics diffuse into the kernel approximation, to the detriment of some of the more important graph harmonics. Moreover, an additional error arises in estimating in (5): as with χK 1.42271pt$ χ$ 1.42271pt$ χ$ 1.42271pt$ χ$ 1.42271pt$ χ$_K in (4), the filter h12h 12 is replaced by a polynomial approximation php_h allowing for the fast estimation of =(ph())⊤ =(p_h( L) Q) (6) with the GWT. A more precise characterization of both errors is given in Sec. 4.E and Sec. 4.F. 4.C. Connection with Randomized SVD and oversampling: Part 1 of Alg. 1 bears strong resemblance with the range-finding part of the RSVD algorithm developed in [6]. In RSVD, the range of the K foremost left singular vectors of a matrix A is estimated with =ortho(). Q=ortho( A G). In the graph kernel case, A corresponds to h()h( L). The main difference with (4) is that our algorithm makes use of the knowledge of the structure of L and h()h( L) to design a filter pχ()p_χ( L) better suited than h()h( L) to the range finding task. Oversampling strategy: In both algorithms (RSVD and ours), the range finding part is inexact and propagates errors. RSVD mitigates this by oversampling the Gaussian matrix G; given r≪Kr K (e.g., r=15r=15), it draws K+r>K+r>K Gaussian random vectors in ∈ℝN×(K+r) G ^N×(K+r). This increases the likelihood that the computed subspace aligns well with the one spanned by the top K singular vectors of A. We adopt the same strategy in our algorithm. Algorithm 1 Randomized Low-Rank Kernel Approximation 1:Graph =(,,w)G=( V, E,w), Laplacian ∈ℝN×N L ^N× N, kernel function h:[0,2]→ℝ+h:[0,2] ^+, parameter K, parameter M, over-sampling parameter r 2:Node embedding s.t. ⊤=~≈=h() = ≈ =h( L). 3:Part 1: Range finding 4:Find λK _K using Alg. 1 from [15]. 5:Compute a degree M polynomial approximation pχp_χ of χK 1.42271pt$ χ$ 1.42271pt$ χ$ 1.42271pt$ χ$ 1.42271pt$ χ$_K 6:Generate ∈ℝN×(K+r) G ^N×(K+r) with Gij∼i.i.d.(0,1)∀i,jG_ij _i.i.d.N(0,1)\;∀ i,j 7:Compute =ortho(pχ()) Q=ortho(p_χ( L) G) 8:Part 2: Embedding computation 9:Compute a degree M polynomial approximation php_h of h12h 12 10:Compute ⊤=(ϕ1,…,ϕN)⊤=ph() =( φ_1,…, φ_N) =p_h( L) Q 4.D. Complexity analysis: In Alg. 1, lines 1, 4, and 6 dominate the computational and memory costs. Filtering one signal with a polynomial of order M requires M matrix-vector products with L, yielding a time complexity of (ME)O(ME) and memory (N)O(N). Hence, [15, Alg. 1] runs in (MElogNlog(ϵ−1))O(ME N (ε^-1)) time and (NlogN)O(N N) memory333The algorithm filters (logN)O( N) random signals (log(ϵ−1))O( (ε^-1)) times., where ϵε denotes the tolerance on the estimation error |λ~K−λK| λ_K- _K . Line 4 costs (MEK)O(MEK) to filter G and (NK2)O(NK^2) for Gram–Schmidt, with memory (NK)O(NK) to store G and Q. Similarly, line 6 requires (MEK)O(MEK) time and (NK)O(NK) memory to filter Q and store . Overall, assuming a fixed tolerance ϵε and K=Ω(logN)K= ( N), the total time complexity of Alg. 1 is (MEK+NK2)O(MEK+NK^2), and the memory complexity is (NK)O(NK). If the matrix ~ is explicitly formed as ~=⊤ = , an additional (N2K)O(N^2K) time complexity and (N2)O(N^2) memory cost are incurred. In practice, explicitly forming ~ is rarely necessary, as matrix–vector products can be efficiently computed as ~=⊤() x= ( x)444This further highlights the advantage of low-rank approximations, reducing the cost of matrix–vector products from (N2)O(N^2) to (NK)O(NK).. 4.E. Polynomial approximation error characterization: Let pχp_χ and php_h be Chebyshev polynomials of degree M, approximating χK 1.42271pt$ χ$ 1.42271pt$ χ$ 1.42271pt$ χ$ 1.42271pt$ χ$_K and h12h 12 respectively. We define ϵχ,M _χ,M :=maxk>K|χK(λk)−pχ(λk)|=maxk>K|pχ(λk)|, := _k>K| 1.42271pt$ χ$ 1.42271pt$ χ$ 1.42271pt$ χ$ 1.42271pt$ χ$_K( _k)-p_χ( _k)|= _k>K|p_χ( _k)|, (7) ϵh,M _h,M :=maxk∈[N]|h12(λk)−ph(λk)|. := _k∈[N]|h 12( _k)-p_h( _k)|. (8) Assumption 1. Let ϵχ,M _χ,M be defined as in (7). We suppose λK≠λK+1and 1/2⩽pχ(λk)⩽1+ϵχ,M∀k⩽K. _K≠ _K+1\;\;and\;\;1/2 p_χ( _k) 1+ _χ,M\;∀ k K. From the Weierstrass theorem and since the spectrum of L is discrete, ϵχ,M _χ,M and ϵh,M _h,M vanish for large M. However, we postpone the analysis of their decay rate to a future study. 4.F. Kernel approximation error: Given from (6) and ~:=⊤ := , we can bound the error ℰ:=‖−~‖E:=\| - \| as ℰ =‖h()−ph()⊤ph()‖ =\|h( L)-p_h( L) Q p_h( L)\| ⩽‖h12()(−⊤)h12()‖ \|h 12( L)( I- Q )h 12( L)\| +‖h12()⊤h12()−ph()⊤ph()‖ +\|h 12( L) Q h 12( L)-p_h( L) Q p_h( L)\| =‖(−⊤)h12()‖2 =\|( I- Q )h 12( L)\|^2 +‖h12()⊤h12()−ph()⊤ph()‖ +\|h 12( L) Q h 12( L)-p_h( L) Q p_h( L)\| :=ℰR2+ℰP. :=E_R^2+E_P. Proposition 1. Under Assumption 1, and with all previous notations holding, we have ℰP⩽2ϵh,M+ϵh,M2E_P 2 _h,M+ _h,M^2 (9) and ℰR⩽h12(λK+1)+2Kr−1pχ(λK+1)+2eK+r(∑j>Kpχ(λj)2)1/2 _R h 12( _K+1)+2 Kr-1p_χ( _K+1)\\ +2e K+rr (Σ _j>Kp_χ( _j)^2 )^1/2 (10) The error ℰPE_P is due to replacing the filter h12h 12 with it’s polynomial approximation php_h, while the error ℰRE_R is due to the range finding. The last two terms in the bound (10) indicate the difference with respect to the error of the optimal rank-K projector ‖(−:K:K⊤)h12()‖=h12(λK+1)\|( I- V_:K V_:K )h 12( L)\|=h 12( _K+1). Proof of (9). Given the definition of ϵh,M _h,M in (8) and that |h(λ)|⩽1|h(λ)| 1, we have ℰP _P =∥h12()⊤(h12()−ph()) =\|h 12( L) Q (h 12( L)-p_h( L)) +(h12()−ph())⊤ph()∥ +(h 12( L)-p_h( L)) Q p_h( L)\| ⩽‖h12()−ph()‖(‖h12()‖+‖ph()‖) \|h 12( L)-p_h( L)\|(\|h 12( L)\|+\|p_h( L)\|) ⩽ϵh,M(1+(1+ϵh,M)).∎ _h,M(1+(1+ _h,M)). Proving (10) is trickier, and we defer the proof to Appendix A. 5 Experimental results (a) (b) (c) Figure 1: Average relative spectral norm of the error w.r.t. the ground truth kernel, ‖−~‖/‖\| - \|/\| \|. Shaded areas represent min. and max. values, over 5 trials. Left (a): As a function of σ, and compared to g-GRFs. Center and right (b-c): As a function of the target rank K, for σ=5σ=5, and compared to the best rank-K and rank-(K+r)(K+r) approximations. Center (b): Swiss-Roll graph. Right (c): Community graph. 5.A. Studying kernels with different bandwidths: We compare our method to the general graph random features (g-GRFs) method555We use stopping probability phalt=0.1p_ halt=0.1, and number of random walkers m=8m=8, which are standard parameters for g-GRFs. [10]. As a first experiment, we compare the accuracy of the kernel approximates on a randomly generated Swiss-Roll graph with 5 000 vertices. Both methods are tested on 10 diffusion kernels =exp(−σ2) =exp(-σ^2 L), which controls the bandwidth of the ground-truth kernels (larger σ yields stronger spectral localization). Fig. 1(a) shows the relative approximation error in these settings666The computation time is not influenced by σ. For this experiment, g-GRFs execute on average in 13.57±0.5613.57± 0.56 s, and our method in 11.88±0.5211.88± 0.52 s., for target rank K=800K=800 and oversampling parameter r=K/10=80r=K/10=80. In all experiments, we use M=30M=30 for php_h and M=60M=60 for pχp_χ. Fig. 1(a) highlights that our method is particularly effective for spectrally localized kernels with narrow bandwidths, i.e., kernels for which only a small number of graph Laplacian eigenvectors contribute significantly, and a low-rank approximation is appropriate. Such narrow spectral bandwidth implies spatially widespread interactions, highlighting our methods ability to capture global geometric properties. In contrast, g-GRFs are better suited to kernels with wide spectral bandwidths, corresponding to spatially localized behavior, and struggle to capture long-range dependencies. Consequently, our proposed approach complements g-GRFs by effectively approximating spectrally localized kernels. 5.B. Effect of target rank K: We now examine the influence of the parameter K, for a randomly generated Swiss-Roll graph with 5 000 vertices and the diffusion kernel with σ=5σ=5. We choose r=max(K/10,15)r= (K/10,15) as oversampling parameter. We compare the results of our method with that of the best rank-K approximation, (K) ^(K) — for reminder, our method is constructed to achieve a similar error. Since our method effectively yields a rank-(K+r)(K+r) approximation, we also compare it to the best rank-(K+r)(K+r) approximation, (K+r) ^(K+r), which yields a lower-bound on the achievable error. Fig. 1(b) shows the average relative approximation error for these settings. As expected, the error decreases as the approximation rank increases. We see that our algorithm performs well on the Swiss-Roll graph, closely matching the best rank-K approximation error up to K≈1 000K≈ 1\,000. For values of K greater than 1 000, our method is not as accurate as the optimal rank-K approximation. At this stage, the error stabilizes around 10−610^-6 (which is perfectly acceptable for most applications). A possible explanation is that small errors in the Gram-Schmidt orthogonalization (which increase with K) affect the final error. The error drops again to machine level precision when K+r⩾NK+r N, where = Q= I and the range finding is perfect. We now repeat the experiment but for a randomly generated Community graph with 5 000 nodes and 8 strongly connected communities weakly interconnected to each other. Fig. 1(c) shows the results for this graph. As one can see, our method fails to match the approximation quality of the best rank-K approximation for values greater than K≈80K≈ 80. Moreover, the error increases back to ≈10−1≈ 10^-1 when reaching K≈500K≈ 500. We surmise that the explanation lies in the structural difference between the Swiss-Roll and Community graphs. For the Swiss-Roll graph, the eigenvalues of the Laplacian are evenly spread out, while the first few eigenvalues of the Community graph are well separated from the others. After K≈500K≈ 500, the eigenvalues are very densely packed, such that the eigenvalue count and the estimation of λK _K become inaccurate. The spectral cut-off applied with pχ()p_χ( L) is not sharp enough in these densely packed areas, such that a significant amount of undesirable frequencies contaminate the most important ones. This reveals that the performance of the algorithm depends strongly on the eigenvalue distribution and the chosen approximation rank. 5.C. Experimental complexity: We also time our experiments on the Swiss-Roll graph, for different values of N and K. In practice, we observe that the two most time consuming steps are the estimation of λK _K and the filtering (the computation of pχ()p_χ( L) G and of ph()p_h( L) Q). Fig. 2 shows the time taken by both these steps777Computing pχ()p_χ( L) G takes roughly twice as much time as computing ph()p_h( L) Q, since we chose polynomials of degree 6060 and 3030 respectively. and compares them with the asymptotic complexities derived in Sec. 4.D. It also compares the total time taken by our algorithm with the time taken to compute the ground truth kernel and the best rank-K approximation (K) ^(K) through explicit eigendecomposition. Figure 2: Swiss-Roll graph. Top left: average computation time as a function of K, for N=5 000N=5\,000. Top right, bottom: average computation time as a function of N, for K=800K=800. Shaded areas represent minimum and maximum values, over 5 trials. We observe that the time complexity experimentally matches that derived in Sec. 4.D, with the estimation of λK _K independent of K (except for K⩾NK N, where λK _K is not needed), the filtering growing in (K)O(K) for fixed N (top-left) and in (N)O(N) for fixed K (top-right). The estimation of λK _K seems to be growing only slightly faster than the predicted (NlogN)O(N N) (top-right), and the total time effectively grows in (N)O(N) for fixed K (and K⩽NK N). Our method takes slightly more time than exact decompositions for small N, but the complexity of those is in (N3)O(N^3) (as confirmed by our experiments — bottom), which is why our method is faster for large N. These results indicate that we can hope to achieve significant computational gains for very large graphs. 6 Conclusion We presented a new method for constructing random graph node embeddings that yield accurate low-rank approximations to specified graph kernels. The approach is especially effective for spectrally localized kernels, which are difficult to approximate with existing techniques. Its computational efficiency further enables kernel-based learning on large-scale graphs. Appendix A Proof of Proposition 1: In this section we prove (10) of Proposition 1. Our arguments are largely inspired by the proof of [6, Theorem 9.1]. Before proving (10), let us introduce some useful Lemmas directly originating from [6]. Lemma 1 ( [6], Proposition 8.1). Let ≽ M 0, and ≽M N M. Then, for any real matrix A, ⊤≽and⊤≽⊤. A M A 0 A N A A M A. Lemma 2 ( [6], Proposition 8.2). Let ≽ M 0. Then, −(+)−1≼ I-( I+ M)^-1 M Lemma 3 ( [6], Proposition 8.3). Let ≽ M 0 be a real matrix partitioned as =[⊤]. M= bmatrix A& B\\ B & C bmatrix. Then, ‖⩽‖+‖.\| M\| \| A\|+\| C\|. Lemma 4 ( [6], Proposition 10.1). Fix matrices S, T, and draw a standard Gaussian matrix G (i.e., such that Gij∼i.i.d.(0,1)G_ij _i.i.d.N(0,1) for all i,ji,j). Then, ‖‖≤‖‖F+‖F‖. _ G\| SGT\|≤\| S\|\| T\|_F+\| S\|_F\| T\|. Lemma 5 ( [6], Proposition 10.2). Draw a K×(K+r)K×(K+r) standard Gaussian matrix G with K≥2K≥ 2 and r≥2r≥ 2. Then, (‖†‖F2)1/2=Kr−1 (E\| G \|^2_F)^1/2= Kr-1 and ‖†‖≤eK+r. \| G \|≤ e K+rr. Proof of (10). Let us first introduce some simplifying notations. We define 1=:K _1= _:K, 1=:K V_1= V_:K, 2=K: V_2= V_K: with =(1,2) V=( V_1, V_2), 2=diag(λK+1,…,λN) _2=diag( _K+1,…, _N), 1=1⊤ G_1= V_1 G, 2=2⊤ G_2= V_2 G, and P_ A as the projector onto the column space of A (for an orthogonal matrix Q, =⊤ P_ Q= Q ). Note that because standard Gaussian matrices are invariant under rotations, ⊤ V G is a standard Gaussian matrix, and 1 G_1 and 2 G_2 are also standard Gaussian matrices. Let us further introduce ~ A 2.79999pt -2.79999pt =⊤h12()=[h12(1)1⊤h12(2)2⊤], = V h 12( L)= bmatrixh 12( _1) V_1 \\ h 12( _2) V_2 bmatrix, (11) ~ B =⊤pχ()=[pχ(1)1pχ(2)2]. = V p_χ( L) G= bmatrixp_χ( _1) G_1\\ p_χ( _2) G_2 bmatrix. (12) Recall that =ortho(pχ()) Q=ortho(p_χ( L) G), such that ortho(~)=⊤ ( B)= V Q. Owing to the unitary invariance of the spectral norm and to the fact that ⊤=⊤ V P_ Q V= P_ V Q, we have the identity ℰR:=‖(−)h12()‖ _R:=\|( I- P_ Q)h 12( L)\| =‖⊤(−)~‖ =\| V ( I- P_ Q) V A 2.79999pt -2.79999pt\| =‖(−⊤)~‖ =\|( I- P_ V Q) A 2.79999pt -2.79999pt\| =‖(−~)~‖. =\|( I- P_ B) A 2.79999pt -2.79999pt\|. (13) In view of (13), we can thus execute the proof for the auxiliary matrix ~ A, and the auxiliary projector P~P_ B. The main argument of the proof is based on matrix perturbation theory. We start with a matrix related to ~ B: =[pχ(1)1]. W= bmatrixp_χ( _1) G_1\\ 0 bmatrix. The matrix W has the same range as a related matrix formed by “flattening out” the spectrum of the top block. Indeed, since (with probability 1) pχ(1)1p_χ( _1) G_1 has full row rank, span()=span[].span( W)=span bmatrix I\\ 0 bmatrix. This immediately yields =[],−=[]. P_ W= bmatrix I& 0\\ 0& 0 bmatrix, I- P_ W= bmatrix 0& 0\\ 0& I bmatrix. (14) In other words, the range of W aligns with the first K coordinates, which span the same subspace as the first K left singular vectors of ~ A 2.79999pt -2.79999pt. Therefore, span()span( W) captures the action of ~ A 2.79999pt -2.79999pt, which is what we wanted from span(~)span( B). We treat the matrix ~ B as a perturbation of W, and hope that their ranges are close to each other. Note that ~≈ B≈ W as pχ(2)≈p_χ( _2)≈ 0. To make the comparison rigorous, let us emulate the arguments outlined in the last paragraph. Referring to (12), we flatten out the top block of ~ B to obtain the matrix =~1†pχ(1)−1=[], Z= B\, G_1 p_χ( _1)^-1= bmatrix I\\ F bmatrix, (15) where =pχ(2)21†pχ(1)−1. F=p_χ( _2) G_2 G_1 p_χ( _1)^-1. The construction (15) ensures that span()⊆span(~)span( Z) ( B), therefore the error satisfies ‖(−~)~‖≤‖(−)~‖.\|( I- P_ B) A 2.79999pt -2.79999pt\|≤\|( I- P_ Z) A 2.79999pt -2.79999pt\|. Squaring this relation, we obtain ‖(−~)~‖2 \|( I- P_ B) A 2.79999pt -2.79999pt\|^2 ≤‖(−)~‖2 ≤\|( I- P_ Z) A 2.79999pt -2.79999pt\|^2 =‖~⊤(−)~‖ =\| A 2.79999pt -2.79999pt ( I- P_ Z) A 2.79999pt -2.79999pt\| =‖h12()(−)h12()‖. =\|h 12( )( I- P_ Z)h 12( )\|. (16) The last identity follows from (11) (~=h12()⊤ A 2.79999pt -2.79999pt=h 12( ) V ) and from the unitary invariance of the spectral norm. To continue, we need a detailed representation of the projector − I- P_ Z. The construction (15) ensures that Z has full column rank, so we can express P_ Z explicitly as =(⊤)−1⊤=[](+⊤)−1[⊤]. P_ Z= Z( Z Z)^-1 Z = bmatrix I\\ F bmatrix( I+ F F)^-1 bmatrix I& F bmatrix. Expanding this expression, we find that the complementary projector satisfies −=[−(+⊤)−1−(+⊤)−1⊤−(+⊤)−1−(+⊤)−1⊤]. I- P_ Z= bmatrix I-( I+ F F)^-1&-( I+ F F)^-1 F \\[2.15277pt] - F( I+ F F)^-1& I- F( I+ F F)^-1 F bmatrix. The partitioning here conforms with the partitioning of . This block matrix expression is less fearsome than it looks. Lemma 2 shows that the top-left block verifies −(+⊤)−1≼⊤. I-( I+ F F)^-1 F F. The bottom-right block satisfies −(+⊤)−1⊤≼ I- F( I+ F F)^-1 F I because (+⊤)−1⊤≽ F( I+ F F)^-1 F 0. We abbreviate the off-diagonal blocks with the symbol =−(+⊤)−1⊤ C=-( I+ F F)^-1 F . In summary, −≼[⊤⊤]. I- P_ Z bmatrix F F& C\\ C & I bmatrix. This relation exposes the key structural properties of the projector. Compare this relation with the expression (14) for the “ideal” projector − I- P_ W. We conjugate the last relation by h12()h 12( ), and invoke Lemma 1 to obtain h12()(−)h12()≼[h12(1)⊤h12(1)h12(1)h12(2)h12(2)⊤h12(1)h12(2)h12(2)].h 12( )( I- P_ Z)h 12( ) \\ bmatrixh 12( _1) F Fh 12( _1)&h 12( _1) Ch 12( _2)\\[2.15277pt] h 12( _2) C h 12( _1)&h 12( _2)h 12( _2) bmatrix. Lemma 1 demonstrates that the matrix on the left-hand side is psd, so the matrix on the right-hand side is also psd. Lemma 3 results in the norm bound ‖h12()(−)h12()‖ \|h 12( )( I- P_ Z)h 12( )\| ≤‖h12(1)⊤h12(1)‖ ≤\|h 12( _1) F Fh 12( _1)\| +‖h12(2)h12(2)‖ +\|h 12( _2)h 12( _2)\| =‖h12(1)‖2+‖h12(2)‖2. =\| Fh 12( _1)\|^2+\|h 12( _2)\|^2. Recall that =pχ(2)21†pχ(1)−1 F=p_χ( _2) G_2 G_1 p_χ( _1)^-1, so ‖h12()(−)h12()‖≤‖pχ(2)21†pχ(1)−1h12(1)‖2+‖h12(2)‖2.\|h 12( )( I- P_ Z)h 12( )\|≤\\ \|p_χ( _2) G_2 G_1 p_χ( _1)^-1h 12( _1)\|^2+\|h 12( _2)\|^2. Inserting this into (16) and using (13), we get ℰR2 _R^2 ⩽‖h12(2)‖2+‖pχ(2)21†pχ(1)−1h12(1)‖2 \|h 12( _2)\|^2+\|p_χ( _2) G_2 G_1 p_χ( _1)^-1h 12( _1)\|^2 Now, Assumption 1 yields ‖pχ(1)−1‖⩽2\|p_χ( _1)^-1\| 2, such that ℰR _R ⩽‖h12(2)‖+‖pχ(2)21†‖‖pχ(1)−1‖‖h12(1)‖ \|h 12( _2)\|+\|p_χ( _2) G_2 G_1 \|\|p_χ( _1)^-1\|\|h 12( _1)\| ⩽‖h12(2)‖+2‖pχ(2)21†‖ \|h 12( _2)\|+2\|p_χ( _2) G_2 G_1 \| (17) We are now ready to use Lemma 4 and Lemma 5, which yield (‖pχ(2)21†‖) (\|p_χ( _2) G_2 G_1 \|) =12(‖pχ(2)21†‖) =E_ G_1E_ G_2(\|p_χ( _2) G_2 G_1 \|) ⩽1(‖pχ(2)‖‖1†‖F+‖pχ(2)‖F‖1†‖) _ G_1(\|p_χ( _2)\|\| G_1 \|_F+\|p_χ( _2)\|_F\| G_1 \|) ⩽‖pχ(2)‖1(‖1†‖F2)1/2+‖pχ(2)‖F1(‖1†‖) \|p_χ( _2)\|E_ G_1(\| G_1 \|^2_F)^1/2+\|p_χ( _2)\|_FE_ G_1(\| G_1 \|) ⩽Kr−1‖pχ(2)‖+eK+r‖pχ(2)‖F. Kr-1\|p_χ( _2)\|+e K+rr\|p_χ( _2)\|_F. Inserting back into (17) concludes the proof. ∎ References [1] A. J. Smola and R. Kondor, “Kernels and regularization on graphs,” in Learning theory and kernel machines: 16th annual conference on learning theory and 7th kernel workshop, COLT/kernel 2003, Washington, DC, USA, august 24-27, 2003. Proceedings. Springer, 2003, p. 144–158. [2] R. I. Kondor and J. Lafferty, “Diffusion kernels on graphs and other discrete structures,” in Proceedings of the 19th international conference on machine learning, vol. 2002, 2002, p. 315–322. [3] A. Rahimi and B. Recht, “Random features for large-scale kernel machines,” Advances in neural information processing systems, vol. 20, 2007. [4] D. I. Shuman, S. K. Narang, P. Frossard, A. Ortega, and P. Vandergheynst, “The emerging field of signal processing on graphs: Extending high-dimensional data analysis to networks and other irregular domains,” IEEE signal processing magazine, vol. 30, no. 3, p. 83–98, 2013. [5] D. K. Hammond, P. Vandergheynst, and R. Gribonval, “Wavelets on graphs via spectral graph theory,” Applied and Computational Harmonic Analysis, vol. 30, no. 2, p. 129–150, 2011. [6] N. Halko, P.-G. Martinsson, and J. A. Tropp, “Finding structure with randomness: Probabilistic algorithms for constructing approximate matrix decompositions,” SIAM review, vol. 53, no. 2, p. 217–288, 2011. [7] N. Tremblay and P. Borgnat, “Graph wavelets for multiscale community mining,” IEEE Transactions on Signal Processing, vol. 62, no. 20, p. 5227–5239, 2014. [8] N. Tremblay, G. Puy, R. Gribonval, and P. Vandergheynst, “Compressive spectral clustering,” in International conference on machine learning. PMLR, 2016, p. 1002–1011. [9] K. M. Choromanski, “Taming graph kernels with random features,” in International Conference on Machine Learning. PMLR, 2023, p. 5964–5977. [10] I. Reid, K. Choromanski, E. Berger, and A. Weller, “General graph random features,” arXiv preprint arXiv:2310.04859, 2023. [11] F. R. Chung, Spectral graph theory. American Mathematical Soc., 1997, vol. 92. [12] U. Von Luxburg, “A tutorial on spectral clustering,” Statistics and computing, vol. 17, no. 4, p. 395–416, 2007. [13] F. X. X. Yu, A. T. Suresh, K. M. Choromanski, D. N. Holtmann-Rice, and S. Kumar, “Orthogonal Random Features,” in Advances in Neural Information Processing Systems, vol. 29. Curran Associates, Inc. [14] E. Di Napoli, E. Polizzi, and Y. Saad, “Efficient estimation of eigenvalue counts in an interval,” Numerical Linear Algebra with Applications, vol. 23, no. 4, p. 674–692, 2016. [15] G. Puy, N. Tremblay, R. Gribonval, and P. Vandergheynst, “Random sampling of bandlimited signals on graphs,” vol. 44, no. 2, p. 446–475.