General Graph Random Features

We thank the reviewer for reading the manuscript and are pleased that they report that the paper reads well. We address their questions and points of minor misunderstanding below.

1. Difference from Choromanski (2023): The reviewer correctly notes that this work extends Choromanski's paper (https://arxiv.org/abs/2305.00156, ICML 2023, oral presentation). Whilst this previous work presents a random feature mechanism for unbiased approximation of the $d$-regularised Laplace kernel, our extension allows unbiased approximation of any function of a weighted adjacency matrix. It is therefore significantly more general; it contains Choromanski's scheme as a special case. Moreover, unlike its predecessor, our method also allows for implicit kernel learning via a neural modulation function. Hence, our work is substantially broader in scope.
2. Accuracy vs compute: Our algorithm has the convenient feature that walkers are independent and straightforwardly distributed across machines. This means that the run time to construct features is directly proportional to the number of walkers, so the plots in e.g. Fig 2 would look identical if the $x$-axis variable were replaced with simulation time. Moreover, from Alg. 1 it is clear that the time-complexity to construct features is subquadratic so scales better than exact methods (cubic). That said, we agree that adding some wall-clock times to compare the exact and approximated speeds might be interesting for practitioners. We are measuring these and will add the results to the appendix.
3. Classes of kernels being approximated: The reviewer is correct that our algorithm can approximate arbitrary functions of a weighted adjacency matrix, which is not the same as any kernel. This is the sense in which they are 'universal', as we state in the first sentence of the abstract. We chose this name to disambiguate from 'graph random features' (https://arxiv.org/abs/2305.00156, ICML 2023) which can only approximate the $d$-regularised Laplacian function. We intended to be unambiguous and welcome suggestions to make our contribution clearer -- does the reviewer think 'general graph random features' would avoid confusion? We would be happy to make any reasonable changes to the name.
4. Normalisation factor: the $\frac{1}{m}$ normalisation factor is actually correct since we are interested in the distribution over subwalks of the average walk. We invite the reviewer to read Sec. A.2 for more detail. This might seem strange compared to e.g. random Fourier features (https://people.eecs.berkeley.edu/~brecht/papers/07.rah.rec.nips.pdf, NIPS 2007) which are normalised by the square root of the feature length, but upon careful inspection this is indeed the right factor.
5. Shaded area and standard deviation: The shaded area is indeed the standard deviation on the estimate of the mean, i.e. the standard deviation of one trial divided by the square root of the number of trials. We will clarify this in the manuscript. We will also clarify bracket notation for standard deviation on the last digit.
6. Table 3 and kernel learning: To produce Table 3, the modulation function was trained for better quality approximation on the smallest graph ('on the small Erdős-Rényi graph ($N = 20$) with $m = 16$ walks, we minimise the loss'). We then freeze the modulation function and test the quality of kernel approximation on the other bigger graphs. Even though $f$ was only trained on to be optimal on the smallest graph (top row of table), it also improves approximation quality for all the remaining graphs (remaining rows), even when $N$ is bigger by a factor of $60$. We will amend the text to make this important point clearer.
7. When the kernel is not defined: The reviewer is correct to observe that if the kernel diverges then we cannot hope to approximate it using finite random features. In the literature, this is avoided by multiplying the weighted adjacency matrix by a regulariser $\sigma^2$ which downscales its spectral radius. We will add further discussion to the manuscript.
8. Symmetric vs asymmetric approaches: We sincerely thank the reviewer for their interesting question. Eq. 4 encodes is the condition for an unbiased estimator; symmetric functions is the special case when we constrain $f_1 = f_2$. This has the benefit of requiring half as many parameters in the neural modulation functions, and empirically during training functions tend to become close to symmetric anyway (see Fig. 6). A detailed study of the effect of different choices among the possible modulation function pairs, including whether the symmetric case is best, is an open and exciting research direction. We will add this discussion and look forward to tackling these questions fully in the future.

We again warmly thank the reviewer and invite any further questions. We hope that, if satisfied, they will consider raising their score.

The reviewer might be interested to see the updated main and supplementary pdfs. Changes are indicated in red. Prompted by their helpful suggestion, we have now added a further experiment in Sec. A.6 plotting the wall-clock times for exact and u-GRF kernel evaluation as the number of nodes grows. For this experiment, exact methods become slower for graphs with more than $\sim 1600$ nodes, and by the time we reach the biggest graph considered ($12800$ nodes) our algorithm is already about $8$ times faster. We also plot the kernel approximation error as $N$ grows to check that the quality of the estimate is not degrading. We hope this clarifies the rough speedups that can be expected from using u-GRFs in different regimes. We have also explained the standard deviation shading and notation, added comments on regularisation and the requirement that the target kernel does not diverge (Sec. A.7), and added discussion of the strengths and weaknesses of symmetric and asymmetric modulation functions. We thank the reviewer again for suggesting these improvements.

We thank the reviewer for their reading of the manuscript and detailed feedback. We are happy that they note the clarity of the writing and wide applicability of u-GRFs. We answer their interesting questions below.

1. Graph size: We consider graph node clustering with $N=3300$ and node attribute prediction with $N=4350$ which are already sufficiently large that exact methods, penalised by cubic time-complexity scaling, become slow (for clustering, by a factor of $> 5$ compared to u-GRFs). Nonetheless, we agree that it would be interesting to continue to push the methods to yet larger graphs. Therefore, following the reviewer's suggestion, we also ran the node attribute prediction experiment on the 'cycloidal' mesh which has $21384$ nodes. In this regime exact kernel evaluation is extremely slow, but once again u-GRFs perform well and the learned kernel (trained on the cylinder graph) does best. Specifically, the normalised differences in prediction error compared to the learned kernel are: $0.065(4)$, $0.129(6)$ and $0.011(1)$ for the $1$-regularised Laplacian, $2$-regularised Laplacian and diffusion kernels respectively. Indeed, throughout the work we see that performance is similar across graph size (see e.g. Fig. 2, where the plots are essentially identical even as $N$ differs by a factor of $60$). However, we will also add a plot to the appendix explicitly comparing number of $N$ and the kernel approximation error; we thank the reviewer for their excellent suggestion.

2. Clustering experiment and metric: We are unsure what the reviewer means by 'groundtruth clusters' for these graphs, but will be happy to respond if they clarify. The pairwise clustering metric in Eq. 15 is a standard way to capture the quality of clustering (the accuracy becomes low if lots of pairs that should be in the same cluster are in different clusters or vice versa), but our scheme would perform well under any reasonable way of comparing the clusterings with the exact and approximated kernels.

3. Modulation function symmetry: The reviewer is correct that Eq. 4 is sufficient for an unbiased estimator of a graph kernel, so $f_1$ and $f_2$ need not in general be equal. The symmetric special case is useful because it needs half the number of learnable parameters compared to the general case when we use neural modulation functions. Moreover, modulation functions initialised differently empirically tend to converge towards being close to symmetric anyway (see Fig. 6 in Sec. A5), suggesting that these give lower-variance estimates. The reviewer is right to identify this as an exciting and challenging direction for future theoretical work; we thank them for their interest.

4. Standard deviation notation: the bracketed number gives the standard deviation of the final digit of a long number more compactly, e.g. $0.123(4)$ means $0.123 \pm 0.004$. We will clarify this notation in the manuscript.

We again thank the reviewer for their time and detailed review. We welcome any further questions.

As a further update, we have now uploaded new versions of the main and supplementary pdfs. Additions and changes are indicated in red. The reviewer might be particularly interested to see Sec. A.6 of the Appendix, which includes a new plot of wall-clock time and approximation error vs the number of nodes in an Erdős-Rényi graph. As suggested, we scale graphs to as big as 12800 nodes where exact computation becomes prohibitively slow. Our method continues to perform well in this regime, with the relative Frobenius norm error between the true and approximated Gram matrices remaining very small. We have also added discussion comparing symmetric and asymmetric modulation functions and clarified the standard deviation notation. We warmly thank the reviewer for suggesting these interesting additions and important clarifications.

We thank the reviewer for their reading of the manuscript and thoughtful comments. We are pleased that they identify the interesting theoretical and numerical results and clear writing. We address all questions and points of concern below.

1. Efficiency gains and the kernel inverse: The reviewer correctly notes that the large cost of kernel inversion is a big motivating factor for constructing random feature mechanisms. In particular, in the presence of explicit random features for a kernel $K$ so that $K \simeq \phi^\top \phi$, we can avoid needing to compute its inverse $K^{-1}$. This is well-established in the random features literature. However, graph kernel methods have the additional challenge that even computing $K$ scales cubically with the number of nodes $N$. Our predominant focus is on overcoming this (using our sub-quadratic algorithm), with the understanding that the other benefits of explicitly-manifested features also follow. We will make this clearer in the manuscript.

2. Bethe Hessian and community detection: The reviewer rightly identifies one of the novel contributions of the paper as a scalable new technique for implicit kernel learning in random feature space. Though this only constitutes one short section (the paper also introduces an unbiased mechanism for kernel estimation and learns the modulation function to improve kernel approximation quality), we sincerely thank them for their interest in this exciting research direction. We agree that there is plenty of scope for interesting work, including their excellent suggestion to apply it to community detection problems and see whether the implicitly learned kernel approximates any popular choices from the literature. We are excited to explore this aspect to u-GRFs further, but for now must defer it as future work.

3. Choice of parameters: When approximating a fixed graph kernel, the only free choice of parameter is $p_\text{halt}$. It encodes a simple trade-off between speed and approximation quality; if $p_\text{halt}$ is lower, walks last longer and sample more subwalks so give a better approximation of graph kernels. In practice, any small value works well (as also reported in https://arxiv.org/abs/2305.00156 which introduced the GRF mechanism -- ICML 2023) so we typically take $p_\text{halt} \sim 0.1$. We will add this discussion to the manuscript.

4. Sparsity and performance: We thank the reviewer for their interesting suggestions. We have included experiments on both synthetic and real-world graphs which naturally have a range of sparsities, and in all cases the performance is similar (see e.g. the plots in Fig. 2 or the clustering errors in Table 2 -- they do not change much between graphs). That said, we agree that it would be interesting to systematically investigate how approximation error changes with sparsity so we are running this extra experiment and will add the results to the Appendix.

5. Additional comments and typos: We thank the reviewer for their observation about Eq. 9 and agree that their suggestion would make the notation clearer, so we have added an $\alpha$ subscript. We have also added the word 'for' to the sentence at the start of Sec. 3.1. In Eq. 14, $\boldsymbol{\phi}_j$ actually refers to the $N \times N$ random feature matrix used to approximate the particular time evolution operator $\exp(\mathbf{W}(t - \tau_j))$ (i.e. the diffusion kernel corresponding to one particular timestep indexed by $j$) rather than $j$ indexing a row. We will clarify this notation and thank the reviewer for flagging this.

We again thank the reviewer for their careful reading of the manuscript and invite them to respond with any further questions. If satisfied, we hope they will consider raising their score.

To provide a further update, we have now uploaded new versions of the main and supplementary pdfs. Changes and additions are indicated in red. Prompted by the reviewer's interesting suggestion, we have included an additional experiment in the appendix (Sec. A.6) investigating the effect of sparsity on the quality of kernel approximation. We fix the number of nodes of an Erdős-Rényi graph and vary the edge-generation probability from 0.1 to 0.9, and measure the relative Frobenius norm error between the true and approximated Gram matrices with the diffusion kernel in every case. The error actually remains very small throughout, but does increase slightly as $p_\text{edge}$ grows. Indeed, it is smallest for sparse graphs which the reviewer rightly identifies as corresponding to 'real data'. We have also added a section (A.7) to discuss the choice of parameters and corrected typos and unclear notation. We thank the reviewer for prompting these additions and improvements.

We thank the reviewer for reading the manuscript. We address questions, concerns and points of minor misunderstanding below.

1. Proof clarification: We are grateful that the reviewer checked the supplementary material, but we respectfully maintain that there is no mistake: Eq. 19 is in fact correct and the theorem follows. We invite the reviewer to consider the paragraph directly after Eq. 19 -- in particular, '$\mathbb{I}$ is the indicator function which evaluates to 1 if its argument is true (namely, $\omega_{iv}$ is a subwalk of $\bar{\Omega}_k^{(i)}$, the $k$-th walk sampled out of [node] $i$), and $0$ otherwise'. To be clear, \omega_i$$_v \in $\bar{\Omega}_k^{(i)}$ is the event that \omega_i$$_v is a 'prefix' of $\bar{\Omega}_k^{(i)}$, which has probability: p(\omega) $$= (1-p)^{\text{len}(\omega)}\prod_1^{\text{len}(\omega)} \frac{1}{d_i} (we dropped the $iv$ subscript because of markdown rendering errors and the product is over $i$). This expression does not need an extra factor of $p$ because we just require that the walk we sample includes \omega_i$$_v as a subwalk, not that it terminates at $v$ -- the walk is free to continue afterwards. Note then that, for the $k$th walk sampled out of node $i$ (denoted $\bar{\Omega}_k^{(i)}$), considering all the $v$ coordinates of the random feature, we get contributions to the vector $\phi(i)$ for every walk \omega_i$$_v that is a prefix of $\bar{\Omega}_k^{(i)}$. This happens $\text{len}(\bar{\Omega}_k^{(i)})$+1 times (all the prefix walks, including the walk of length $0$) at different coordinates $v$, which corresponds exactly to the updates at every timestep the reviewer identifies in line 8 of Alg. 1. The crux here is the relationship between \omega_i$$_v (the variable in the sum over all the possible walks between $i$ and $v$ on the graph) and $\bar{\Omega}_k^{(i)}$ (a particular finite-length walk we have sampled out of $i$) via the event \omega_i$$_v \in \bar{\Omega}_k^{(i)} (\omega_i$$_v is a prefix subwalk of $\bar{\Omega}_k^{(i)}$). On careful inspection, Eq. 19 then corresponds exactly to the estimator generated by the algorithm. We invite the reviewer to respond with further questions.

2. Experiments: The reviewer is correct that many of our experiments compare results with the approximated kernel to the groundtruth. This is because our central contribution is a novel algorithm for unbiased and efficient approximation of arbitrary functions of a weighted adjacency matrix. A natural way to check whether the approximation is good by comparing the exact and approximated results. This includes on big 'real' graphs where exact methods are already slow due to the cubic-time complexity scaling (e.g. 'citeseer' has $3300$ nodes and 'torus' has $4350$ nodes for the clustering and kernel regression experiments respectively). Here, the speedup from using u-GRFs is already a factor of $> 5$. Following the reviewer's comment we ran the node attribute prediction task again on the 'cycloidal' graph which has $21384$ nodes. The kernel we trained on the small 'cylinder' graph again performed best, with normalised differences in prediction error of: $0.065(4)$, $0.129(6)$ and $0.011(1)$ for the $1$-regularised Laplacian, $2$-regularised Laplacian and diffusion kernels respectively. This is a very big graph for which exact kernel computation is prohibitively slow, but our scaleable u-GRF method continues to perform well.

3. 'Universal' graph random features: We added the word 'universal' to disambiguate from 'graph random features' (https://arxiv.org/abs/2305.00156, ICML 2023, oral presentation) which can only approximate the $d$-regularised Laplace function. In the first sentence of the abstract we state that this refers to the ability to approximate arbitrary functions of a weighted adjacency matrix. In the second sentence we clarify that this includes 'many' (not all) popular graph node kernels. We intended to be unambiguous about this fact and strongly agree that it is important to be clear about contributions, so we welcome feedback and suggestions for how to make this more obvious. Does the reviewer think that e.g. 'general graph random features' would fit better? We are happy to rename them to any sensible suggestion.

We thank the reviewer again and invite them to respond with further queries. We hope the clarifications will be sufficient for them to consider raising their score.