Efficient Graph Continual Learning via Lightweight Graph Neural Tangent Kernels-based Dataset Distillation

Q1: Why use Bernoulli distribution than largest r eigenvalues?

We adopt Bernoulli sampling because it effectively balances high-frequency and low-frequency components, preserving both local and global graph structures. In contrast, directly selecting the largest $r$ eigenvalues mainly captures local variations, while selecting the smallest $r$ eigenvalues emphasizes global structures.

1. Theoretical Justification

Under Lipschitz gradient assumption [1], we assume that the GNN $f_\theta$ satisfies

$$\|\nabla_\theta f_\theta(A) - \nabla_\theta f_\theta(B)\| \leq L_\nabla \|A - B\|,$$

$L_\nabla$ is the constant. Therefore, when approximating $L$ with $\tilde{L}$, the induced error in the gradients $\| \nabla f_\theta(LX) - \nabla f_\theta(\tilde{L} X) \|$ is controlled by $\| f_\theta(LX) - f_\theta(\tilde{L} X) \|$. Given eigen-decomposition: $L = U \Lambda U^\top, \tilde{L} = \tilde{U} \tilde{\Lambda} \tilde{U}^\top,$ where $\tilde{L}$ is a low-rank approximation, we introduce a selection indicator $\xi_i$ for each eigenvalue $\lambda_i$, leading to an error formulation:

$$(L-\tilde{L})X = \sum_{i=1}^{n}(1-\xi_i)\lambda_i u_i (u_i^\top X).$$

To ensure a fair comparison, we analyze the normalized error bound:

$$\mathbb{E}_\xi\left[\frac{\|(L-\tilde{L})X\|}{\|LX\|}\right] = \frac{\sum \mathbb{E}[1-\xi_i]\ \lambda_i^2\|u_i^\top X\|^2}{\sum \lambda_i^2\|u_i^\top X\|^2}$$

For Bernoulli, where $\xi_i \sim \text{Bernoulli}(p)$ with $p = \frac{r}{n}$, $\mathbb{E}[1-\xi_i] = 1-p \approx 1-\frac{r}{n}.$

Thus, the normalized error bound is: $\text{Error}_{\text{Bern}} \approx 1 - \frac{r}{n}.$ This result is independent of the eigenvalue distribution, ensuring robustness across different graphs.

In comparison,

$$\frac{\sum_{i=r+1}^{n}\lambda_i^2\|u_i^\top X\|^2}{\sum_{i=1}^{n}\lambda_i^2\|u_i^\top X\|^2} = \text{Error}_{\text{Largest}}.$$ $$\frac{\sum_{i=1}^{n-r}\lambda_i^2\|u_i^\top X\|^2}{\sum_{i=1}^{n}\lambda_i^2\|u_i^\top X\|^2} = \text{Error}_{\text{Smallest}} .$$

These methods depend on the eigenvalue spectrum; simply dropping local or global signals will lead to significantly larger errors [2, 3].

2. Experimental Validation

We conducted a comparative analysis of various sampling strategies:

Experimental results demonstrate that Bernoulli sampling achieves the best performance in preserving both global and local graph information.

[1]LECTURES ON LIPSCHITZ ANALYSIS.

[2]Self-supervised graph-level representation learning with local and global structure. PMLR, 2021.

[3]From local to global: Spectral-inspired graph neural networks. NeurIPS 2022.

---

Q2: Explanation of low-rank optimization

We acknowledge that our initial complexity analysis contained an error. Below, we provide a corrected and more detailed explanation.

The efficiency improvement of optimization primarily stems from reducing the computational complexity of matrix operations. Instead of directly computing $N \times N$ dense matrix $L$, we use a sampled low-rank approximation: $L \approx U_r \Lambda_r U_r^T,$ where: $U_r \in \mathbb{R}^{N \times r}, \Lambda_r \in \mathbb{R}^{r \times r}.$

To compute $\tilde{L} X = U_r \Lambda_r U_r^T X$ efficiently, we decompose into 3 sequential steps:

• 1. $Z = U_r^T X, \quad Z \in \mathbb{R}^{r \times d}, \quad \text{cost: } O(rNd)$

$Y = \Lambda_r Z = \Lambda_r (U_r^T X), \quad Y \in \mathbb{R}^{r \times d}, \quad \text{cost: } O(r^2 d)$

• 3. $\tilde{L} X = U_r Y = U_r (\Lambda_r (U_r^T X)), \quad \text{cost: } O(Nrd)$

Thus, overall computational complexity: $O(rNd) + O(r^2 d) + O(Nrd) \approx O(Nrd),$ where $r \ll N$. This significantly reduces the complexity compared to the naive approach $O(N^2 d)$. This low-rank approximation approach is conceptually related to linear attention in Transformers [4].

[4]Transformers are RNNs: Fast Autoregressive Transformers with Linear Attention. ICML 2020.

---

Q3: Unclear definition of Θ(Gi)

We will correct this inconsistency in the notation of Θ in the final version of the paper by carefully distinguishing between the kernel function and the matrix and updating notation table to ensure clarity.

---

Q4: Definition of equation (5)

Equation (5) is a combination of two terms:

• Cross-entropy loss: $\frac{1}{M} \sum_{i=1}^{M} - o_i \log(\hat{o}_i)$ measures the difference between true labels $o_i$ and predicted probabilities $\hat{o}_i$.
• Regularization: $\eta \|\Theta\|^2$, is a regularization term that adds a penalty proportional to the squared $L_2$ norm of model parameters $\Theta$ to prevent overfitting.

---

S1&S2: Typo Errors

We will carefully proofread and improve the manuscript. Thanks.

We sincerely appreciate your comments, suggestions, and every effort spent on reviewing our work. Here we attempt to address all your remaining concerns. In the following, we quote your comments and then give our detailed response point-by-point.

---

W1: Ablation studies on key components are missing.

We acknowledge that our initial submission lacked comprehensive ablation studies. To rigorously validate this, we have performed extensive ablation experiments focusing on the sampling probability and the use of layer-specific gradients.

1. Sampling Probability

We investigated the impact of dataset distillation quality and computation cost under different sampling rates for Bernoulli sampling:

Table 1. Performance Comparison among Different Sampling Rates

Table 2. Time consumption (seconds) among Different Sampling Rates

The results indicate that a sampling probability of 0.1 achieves the best balance between training efficiency and model accuracy. Lower probabilities leed to performance decay due to inadequate training data, while higher probabilities increase computational cost without significant gains in performance.

2. Layer-Specific Gradients

We analyzed the effect of employing layer-specific gradients by comparing it with LIGHTGNTK that uses gradients across all layers, first layer and last layer.

Table 3. Performance of gradient computation using layer-specific gradients

The results indicates using gradients across all layers can provide best performance, while choosing gradients from the last layer offers balance between accuracy and efficiency.

---

W2: Provide computational efficiency of LIGHTGNTK.

We compared the training time of GNTK and LIGHTGNTK on more datasets. The results, detailed in the table, indicate that LIGHTGNTK achieves an average of 30% reduction in training time.

Table 4. Time consumption (seconds) between GNTK and LIGHTGNTK

---

W3: More discussion of Graph Representation Learning models

To comprehensively evaluate the impact of different pretrain methods, we conducted additional experiments using graph contrastive learning (GraphCL [1], CI-GCL [2]) to pretrain our GNN backbone, comparing it with GPT-GNN-based pretraining approach.

Table 5. Performance Comparison under Different Pretrain Methods

The experimental results indicate the generality of our LIGHTGNTK, which is effective under different GNN backbones pretrain methods, and the stronger the ability of the base model, the better the quality of data selecting.

[1] Community-invariant graph contrastive learning. ICML.

[2] Graph contrastive learning with augmentations. NeurIPS.

We sincerely appreciate your comments, suggestions, and the time spent reviewing our work. Below, we address each of your concerns in detail.

---

W1: Confusion in using symbols $\theta$ and $\Theta$ throughout the text

We sincerely apologize for the inconsistency in the notation of $\theta$ and $\Theta$ in our manuscript. This was an oversight on our part, and we greatly appreciate your feedback on this issue.

To address this, we will carefully distinguish between the kernel function and the matrix in the final version of the paper. Additionally, we will update the notation table to ensure clarity and maintain consistency throughout the manuscript.

---

W2: The absence of parameter sensitivity analyses renders the experimental section incomplete

We acknowledge that our initial submission lacked comprehensive ablation studies. To address this, we have conducted extensive experiments analyzing the impact of sampling probability and layer-specific gradients.

1. Sampling Probability Analysis

We evaluated dataset distillation quality and computational cost under different sampling probabilities in Bernoulli sampling:

Table 1. Performance Comparison Across Different Sampling Rates

Table 2. Computation Time (seconds) Across Different Sampling Rates

Our findings indicate that a sampling probability of 0.1 achieves the optimal balance between training efficiency and model accuracy.
A lower sampling rate leads to a decline in performance due to insufficient training data, whereas a higher rate increases computational costs without significant performance gains.

2. Layer-Specific Gradient Analysis

We further examined the impact of employing layer-specific gradients by comparing LIGHTGNTK's performance when using gradients across all layers, only the first layer, and only the last layer.

Table 3. Performance Comparison of Layer-Specific Gradient Computation

The results demonstrate that utilizing gradients across all layers yields the best performance, while using gradients from only the last layer provides a favorable trade-off between accuracy and efficiency.

---

W3: The validation-test distribution assumption constitutes a strong prior that lacks sufficient empirical support

We acknowledge the concern regarding Assumption 1 (validation-test distribution consistency) and would like to clarify this assumption both theoretically and empirically.

1. Theoretical Justification

From a theoretical perspective, assuming an identical distribution across training, validation, and test sets is a widely adopted principle in machine learning, particularly when data splitting is performed randomly under the i.i.d. (independent and identically distributed) assumption. This practice aligns with foundational theories established in various works, such as [1][2].

2. Empirical Verification

To further substantiate this assumption, we conducted quantitative analyses using Maximum Mean Discrepancy (MMD) and Kullback-Leibler (KL) divergence on the learned GNN embeddings after training.

Table 4. MMD and KL Divergence between Validation and Test Sets

Following [3][4], the MMD and KL divergence values are below the critical thresholds commonly used to assess distributional shifts. These results provide strong empirical support that the validation and test sets follow statistically similar distributions, thereby validating the i.i.d. assumption.

[1] Sampling: Design and Analysis. Chapman and Hall/CRC, 2021.
[2] Asymptotic Properties of Random Restricted Partitions. Mathematics, 2023.
[3] Prompt-based Distribution Alignment for Unsupervised Domain Adaptation. AAAI, 2024.
[4] KL Guided Domain Adaptation. ICLR, 2022.

We sincerely appreciate your comments, suggestions, and the time spent reviewing our work. Below, we address each of your concerns in detail.

---

W1: Concepts of GNTK have been previously discussed

Thank you for your valuable comment. While the concepts of GNTK have indeed been explored in KIDD [1], our approach in LIGHTGNTK differs significantly. KIDD is specifically designed for graph classification and employs GNTK for kernel ridge regression-based dataset distillation. However, the kernel regression paradigm in KIDD does not guarantee that training graph samples remain relevant during the inference stage.

In contrast, our LIGHTGNTK framework unifies graph classification, node classification, and link prediction tasks. We introduce a Bernoulli-sampled low-rank Laplacian-approximated GNTK to enhance similarity-based dataset distillation. This novel formulation enables efficient, scalable, and task-adaptive distillation beyond the scope of KIDD.

[1] Kernel Ridge Regression-Based Graph Dataset Distillation

---

W2: GNTK similarity matrix may become computationally intensive

Thank you for raising this concern. While constructing the GNTK similarity matrix has a theoretical complexity of $O(N^2)$, we mitigate computational costs by selecting only the most informative samples and reducing redundant information, following the LESS [2] framework to maintain distillation performance.

Furthermore, to further optimize computational efficiency, we explore pooling-based methods (e.g., mean pooling) to derive a coarse-grained validation set and a smaller similarity matrix. The experimental results (on node classification tasks using CiteSeer and Cora) are presented below:

Mean-Pooling Results:

Pairwise Calculation Results:

The results indicate that mean pooling retains a considerable portion of the performance achieved by pairwise computation while significantly reducing computational costs. Although pairwise similarity remains superior, mean pooling presents a viable alternative in scenarios where computational efficiency is a priority.

[2] LESS: Selecting Influential Data for Targeted Instruction Tuning

---

W3: Why could GNTK capture both structural and feature relationships?

Thank you for your insightful question. The ability of GNTK to capture both structural and feature relationships is rooted in its gradient-based dataset distillation mechanism. In our LIGHTGNTK framework, structural information is propagated through the Laplacian matrix $L$, while feature information is encoded via the weight matrix $W$.

Mathematically, in standard (non-graph) dataset distillation, the gradient of the loss function with respect to the weight matrix at layer $l$ is given by:

$$\frac{\partial\mathcal{L}}{\partial W^{(l)}} = \frac{\partial\mathcal{L}}{\partial Z^{(l)}} \cdot a^{(l-1)T}$$

where $a^{(l-1)}$ represents the activations from the previous layer, and $T$ denotes the matrix transpose to ensure correct dimensional alignment.

However, in the GNTK-based formulation, the gradient incorporates graph structural dependencies via the Laplacian matrix $L$ and is expressed as:

$$\frac{\partial\mathcal{L}}{\partial W^{(l)}} = H^{(l-1)T} L \left( \frac{\partial\mathcal{L}}{\partial Z^{(l)}} \odot \sigma^{\prime} (Z^{(l)}) \right)$$

where:

• $H^{(l-1)}$ is the node representation at layer $(l-1)$.
• $L$ is the graph Laplacian, encoding structural relationships.
• $\sigma^{\prime} (Z^{(l)})$ represents the element-wise derivative of the activation function.

From this formulation, it is evident that GNTK-based gradients inherently embed structural information through $L$ while simultaneously encoding feature dependencies via $H^{(l-1)}$ and the weight parameters. This allows GNTK to effectively distill graph-structured datasets while preserving essential structural and feature relationships.