Temporal Graph Neural Tangent Kernel with Graphon-Guaranteed

Thanks very much for the review! We greatly appreciate your enjoy on our paper's insights, theoretical analysis, and experiments. We address your questions in the form of Q&A below.

Moreover, for the length limit, we distillate our answer and provide brief but important answers, more detailed intermediate steps, reasonings, and derivations can be provided during the upcoming author-reviewer discussion phase. All the detailed reply and supplementary experiments will be added to our updated paper.

Q1. How do you deal with the node features in Equation (1)?

The reason we omit node features in Equation (1) is that prevalent benchmarks only support time-aware edge features.

But if node features are given, our framework still works. Theoretically, Equation (1) can be extended as follows: $h_v(t) = c \cdot \sum_{(u, \bar{t})}[\mathbf{t}\_{enc}(t - \bar{t}) || \mathbf{e}\_{uv}(t) || \mathbf{x}_{u}(\bar{t})]$, where $\mathbf{x}_u(\bar{t})$ is the node features of $u$ at time $\bar{t}$.

Adding node features will not affect our conclusions and framework, as the construction of our kernel depends on $h_v(t)$ as a general term, instead of its components, such as node / edge features. Moreover, Theorem 5.1 and Theorem 5.2 still hold if node features are given.

Q2. How weights W in Equation (6) are initialized

For Equation (6), the weights $\mathbf{W}$ are initialized i.i.d from the normal distribution $\mathcal{N}(0, 1)$.

Q3. To obtain (6) and (13), aren't you actually training a different GNN for each time $t$?

We want to clarify that our method does not involve training GNNs.

The kernel function in Equation 6 is just a starter and illustrated by a temporal graph neural network for easy interpretation, because the temporal graph neural network exists earlier than the proposed temporal graph neural tangent kernel (i.e., ours).

The derivation of our kernel function is given in Sec. 4, especially in Equation (12), which just relies on the covariance matrix computation through Equation (11), to Equation (10), to Equation (9), to Equation (8), and finally to Equation (7).

A typo is found during the rebuttal. For the last term in Eq (11), $l-1$ should be $l$.

We will strengthen the above derivation in the updated paper.

Q4. The kernel (12) does not capture the temporal relationships between different time instances. I am not sure how this models or compares with typical temporal graph learning approaches ...

Temporal dependencies are captured in Equation (1), where we construct the temporal node representations for a node $v$ at time $t$ by aggregating information (node, edge features, and time difference) from its previous temporal neighbors $(u, \bar{t}), \bar{t} < t$. Intuitively, this models the current state of $v$'s neighborhood at time $t$ based on its own neighborhood at the previous time $\bar{t}$. Since we would want a vector representation for $v$, we perform vector summation over all information of temporal neighbors $(u, \bar{t})$.

It is true that there is an established line of works in temporal graph learning using recurrent architecture and attention to learn temporal dependencies like TGN and TGAT, but there are works employing simpler architectures, such as GraphMixer, which aggregates information from recent temporal neighbors and processes them with MLP-Mixer layers.

Q5. Why do you not consider the kernel that takes the whole temporal graph sequence and as inputs? I can understand that this is more sophisticated ...

If we understand your question correctly, your concern is about Equation (13) that current time prediction only depends on other samples at the current time.

We admit that taking the historical behavior as you suggested is a promising approach, and actually this is what we have done in the paper. Instead of taking each past snapshot individually, recall our Equation (1), temporal dependencies are captured, where we construct the temporal node representations for a node $v$ at time $t$ by aggregating information (node, edge features, and time difference) from its previous temporal neighbors $(u, \bar{t}), \bar{t} < t$. Hence, the current time representation is the cumulation of the past. In this angle, the sequence of the temporal graph is observed by our method.

Q6. Table 1: I fail to see why Temp-G3NTK is computationally cheaper than GraphMixer. (Also, you have not defined $K$ for GraphMixer.) ...

First, $K$ for GraphMixer is the padding length, i.e., the maximum number of recent interactions a node can sample.

The theoretical complexity compassion is not that straightforward by only looking at the power to number of nodes. As we stated in the paper, GraphMixer needs the complex neural computation like gradient descent and back propagation, but our method does not. That’s why we claim our method is faster.

To in-depth verify our claim, we added empirical running time comparison of all baselines across all datasets. The results are reported as Table 1 in the 1-page pdf in “general rebuttal”. According to the result, our method is roughly 10x-20x faster than GraphMixer, and also faster than most baselines.

We will include these new results in our paper to in-depth verify my hypothesis.

Q7. GraphMixer is not included in Figure 1 because it is "designed to directly process continuous-time dynamic graphs, which require observing the entire length of the input temporal graph instead of processing snapshots". I don't think this claim is accurate. You can easily implement GraphMixer by restricting it to snapshots.

Technically, GraphMixer can process discrete time dynamic graphs, by only considering part of the past snapshots as the continuous dynamic graphs, although it was originally designed for continuous graphs.

Since we can use snapshot timestamp as edge timestamp for each edge in this snapshot, then we are now running GraphMixer and will add its plot during the author-reviewer discussion phase.

We are happy to report that the additional running experiments required for GraphMixer with the rest 7 baselines are now finished.

Due to the rebuttal policy, we can not add an external link to direct you to check the figure. However, we describe the overall pattern and provide the important data points in table format as follows. The new figure and corresponding analysis are promised to be added to our paper.

• After involving GraphMixer, the general pattern in Figure 1 is not changed, i.e., our method still outperforms baselines, and the gap is obvious.

• For providing specific data points, we extract the important comparison, i.e., GraphMixer vs. Ours, and made the tabular representation as below. According to the table, we can see our method outperformed GraphMixer in most cases.

Table. For involving GraphMixer in Figure 1 of the paper, the temporal graph classification accuracy on different stages of temporal graphs.

Thanks for your appreciation, and we are glad to answer your follow-up question!

Performance of complex neural architectures and simple neural architectures

According to the most recent temporal graph learning work, e.g., MLP-Mixer [1] and GraphMixer [2], there is a trend to replace the complex neural architecture in the representation learning process, but using simplified neural architecture to maintain or even boost the performance.

According to GraphMixer [2], learning positional encoding plus simple MLPs can achieve very competitive performance. For example, in the link prediction task on Reddit datasets, GraphMixer achieved 99.93% precision and outperforms the 99.30% of JODIE [3] (using recurrent architecture) and the 99.80% of TGN [4] (also using recurrent architecture).

Motivated by this direction, we designed our temporal graph neural tangent kernel. On the one hand, we even simplified the complexity of SOTA temporal graph neural networks like GraphMixer [2], TGN [4] (using recurrent arch), and DyRep [5] (using recurrent arch), and outperformed no matter

• in the theoretical complexity (as shown in Table 1 of the paper);
• or in the empirical complexity (as shown in Table 1 of the 1-page PDF file in “general rebuttal”);
• or in the theoretical effectiveness (as shown in Theorem 5.1 and entire Appendix B for proof);
• or in the empirical effectiveness (as shown in Table 2 of the paper for graph-level tasks and in Table 3 of the paper for node-level tasks).

Moreover, besides our two existing baselines that already used recurrent architectures (TGN [4] and DyRep [5]), we also added your recommended work to our comparison, and the result is shown below.

Table. Temporal Graph Classification Accuracy on Infectious Dataset

From the above table, we can see that EvolveGCN performs similarly with another recurrent architecture baseline TGN, which aligns with the analysis [1,2] that complex neural architectures can be replaced by simple architectures in the temporal graph learning community. Since it is replaceable, we next expand the possible trade-off.

Possible Trade-off

Here, we would like to share the trade-off of this evolution trend in the temporal graph learning community. To simplify the neural architecture but most importantly maintain (or even boost) the performance, it requires the researchers to devote time and efforts to propose new paradigms, just like what we did in the paper to replace the neural computing with exact mathematical expressions and ensure the theoretical and empirical effectiveness.

On the other hand, our method still triggers the future direction for the temporal graph learning community. For example, even though we already achieved theoretical and empirical effectiveness and efficiency, we discerned that the room still exists to do better for future research. For example, in our Equation (13), making the decision for the testing data, it usually needs a bunch of observed training data. Even though our computation is simple and not related to neural computation like gradient descent and back propagation, we still wish to decrease the number of training data to meet the specific case if the training data is inadequate.

Here, we would like to share a possible solution. We can involve the prototype learning in a non-neural way, such that a bunch of training data can be presented by a single prototype representation. Hence, the number items in Equation (13) can be decreased and the efficiency can be largely increased. This is an interesting topic, and we would like explore it as our future research direction.

Reference

[1] A Generalization of ViT/MLP-Mixer to Graphs, ICML 2023

[2] Do We Really Need Complicated Model Architectures For Temporal Networks, ICLR 2023

[3] Predicting Dynamic Embedding Trajectory in Temporal Interaction Networks, KDD 2019

[4] Temporal graph networks for deep learning on dynamic graphs, arXiv 2020

[5] Dyrep: Learning representations over dynamic graphs, ICLR 2019

Thanks very much for the review! We greatly appreciate your acknowledgment of our paper's motivation, theoretical analysis, and experiments. We address your questions in the form of Q&A below.

Moreover, for the length limit, we distillate our answer and provide brief but important answers, more detailed intermediate steps, reasonings, and derivations can be provided during the upcoming author-reviewer discussion phase. All the detailed reply and supplementary experiments will be added to our updated paper.

W1. Expanding the related work section

On one hand, we expand the detailed introduction to Graphons, Graphon Neural Network, and Graphon Neural Tangent Kernel.

For example, following details will be added to Section C.1 of our Appendix.

Graphons. A graphon is a symmetric, bounded, and measurable function $W : [0, 1]^2 \rightarrow [0, 1]$ that acts as the limit object of dense graph sequences and defines a family of similar graphs. Specifically, given an undirected graph $\mathbf{F}$ and a graph sequence $\{G_n\}$, whose $i$-th graph in the sequence, $G_i$, has $i$ nodes, then $\{G_n\}$ converges to $W$ if $\lim_{n \rightarrow \infty} t(\mathbf{F}, G_n) = t(\mathbf{F}, W)$, where $t(\mathbf{F}, G_n)$ defines the density of homomorphisms between $\mathbf{F}$ and $G_n$, and $t(\mathbf{F}, W)$ can be defined similarly. In this way, $W$ can also be regarded as a generative models for a certain graph family.

Graphon Neural Networks (WNNs) [26]. Similar to how GNNs perform graph convolutions on a graph given the adjacency matrix, $A$, and node / edge features $H$, Graphon Neural Networks (WNNs) also perform graphon convolutions using $W$ and $X$, where $X: [0, 1] \rightarrow \mathbb{R}$ (corresponds to node-level) or $X: [0, 1]^2 \rightarrow \mathbb{R}$ (corresponds to edge-level) is graphon signals. Intuitively, we could consider WNNs as the counterpart of GNNs that operate on graphons. As proven in [26], as the graph's size of the graph sequence $\{(G_n, H_n)\}$ generated by $W, X$ grows, i.e $n \rightarrow \infty$, the output of GNNs operating on $\{(G_n, H_n)\}$ converges to the output of WNNs operating on $W, X$. More details can be found in Section C.1 in our Appendix.

On the other hand, we include more temporal graph representation learning methods discussion in our paper, like CAWN [1], and DySAT [2].

For example, following details will be added to Related Work of our paper.

* CAWN obtains the causal walks for each node, anonymize the node indices on each walk, further process these walks with recurrent architectures, and the final temporal node representation is determined by aggregating walks encodings.

* DySAT first transforms the given CTDG into discrete temporal graph snapshots, then applies Graph Attention Network (GAT) on each snapshots to obtain the spatial information, $\mathbf{X}(t)$, where $t$-th corresponds to the graph snapshot at time $t$, and denote the spatial information for node $u$ in the $t-$th graph snapshot as $[\mathbf{X}(t)]_u$. Then $u$'s temporal representation is obtained by applying Transformer on its spatial information at different timestamps: $h_u(t_k), h_u(t_{k - 1}), \dots, h_u(t_0) = \text{Transformer}([\mathbf{X}(t_k)]_u, [\mathbf{X}(t_{k - 1})]_u, \dots, [\mathbf{X}(t_0)]_u)$.

Finally, we would be more than willing to discuss any related work the reviewer recommend during the upcoming author-reviewer discussion phase.

Reference:

[1] Inductive Representation Learning in Temporal Networks via Causal Anonymous Walks, ICLR 2021

[2] DySAT: Deep Neural Representation Learning on Dynamic Graphs via Self-Attention Networks, WSDM 2020

W2. More explanation regrading the simplicity and interpretation of the proposed method

For simplicity, as stated in the paper, our kernel method does not involve complex neural architecture or training algorithms. Moreover, all terms appearing in the method can be computed in an iterative manner as shown from Equation 7 to Equation 11.

Additionally, we added the empirical running time for comparing the efficiency of our method and 8 baselines across 4 datasets, as shown in Table 1 in the 1-page pdf in "general rebuttal". Overall, the empirical running time aligns with our theoretical analysis of time complexity in Table 1 of the paper. That is, our method belongs to the graph kernel category, where the node-wise comparison is usually inevitable, and our time complexity is lower; Compared to the neural baselines, since our method does not rely on complex neural training like backpropagation, our method is still efficient.

In brief, the interpretation of our graph neural tangent kernel relates to the training dynamics of the temporal graph neural architecture that is defined in Section 3.1. Specifically, in the limit that the layer width of MLPs goes towards infinity, then the graph prediction using the architecture defined in Section 3.1 is equivalent to the prediction of using kernel regression, as stated in Equation 13.

Q1. This is good work. Is it possible to use more graph neural architecture-based methods in experiments for two-level tasks? Also the datasets selection can be more inclusive.

As stated in Section 6.1, for graph-level datasets, we cover a variety of different domains from DBLP (citation graph), Facebook (social graph), Infectious (human-contact graph), and Tumblr (social graph). Moreover, in Appendix D.4, we also show how our methods can scale on larger datasets (Wiki, Reddit, MOOC, LastFM), in terms of number of nodes and edges.

Our baselines include static graph kernel methods, static GNN methods, temporal GNN methods, and ours stand for the temporal graph kernel method.

During the rebuttal, we also include all baselines into a new graph dataset from https://chrsmrrs.github.io/datasets/docs/datasets/, the performance is shown in Table 2 in 1-page pdf in “general rebuttal”. According to the new results, our method still achieves the best performance.

Thanks very much for the review! We greatly appreciate your acknowledgment of our paper's novelty, presentation, and model performance. We address your questions in the form of Q&A below.

Moreover, for the length limit, we distillate our answer and provide brief but important answers, more detailed intermediate steps, reasonings, and derivations can be provided during the upcoming author-reviewer discussion phase. All the detailed reply and supplementary experiments will be added to our updated paper.

Q1. Might be good to add plots or tables showing the compute time for each dataset

We have added the empirical running time of our method and 8 other baselines on 4 benchmark datasets. The table is plotted as Table 1 in the 1-page pdf file in the “general rebuttal”.

Overall, the empirical running time aligns with our theoretical analysis of time complexity in Table 1 in the paper. That is, our method belongs to the graph kernel category, where the node-wise comparison is usually inevitable, and our time complexity is lower. Compared to the neural baselines, since our method does not rely on complex neural training like backpropagation, our method is still efficient.

Given our method achieved the best classification accuracy, as shown in Table 2 in the paper, according to the corresponding running time reported in Table 1 of the uploaded 1-page pdf during the rebuttal, our method is (1) more than 10x - 20x faster then complex temporal graph neural network methods like GraphMixer [6] and TGN [24]; (2) similarly efficient as simple kernel methods like WL-Subtree [27] and Shortest Path [4] and embedding methods like NetLSD [30] and GL2Vec [5]; and only Graph2Vec [21] is running faster than our method, but our performance is roughly 1.4x better.

Q2. Is it possible to adapt our method to continuous time dynamic graphs (CTDGs)?

The answer is yes.

In the paper, we choose discrete time just for a better illustration of the proposed method. The transformation between discrete time and continuous time is possible, as described in lines 63-65 of the paper.

In other words, given 2 CTDGs $G, G'$, suppose that we already computed the Temp-G$^{3}$NTK value between $G, G'$ at time $t$. If new links occurring at time $t', t' > t$ are added into $G, G'$, then we update the neighborhood of each node in $G, G'$, and re-compute the Temp-G$^{3}$NTK value for $(G, G')$ at time $t'$.

An extreme case needs more caution, i.e., spare continuous time dynamic graph. For example, only one edge is updated per timestamp. This scenario will make our method run slower, because more snapshots are emerging. But multiple edges updated per timestamps can be adapted seamlessly. Moreover, one edge update per timestamp is very rare no matter in the real-world setting.

Q3. How would the model convergence be if the network both grows and shrinks over time. This is quite common in real world temporal graphs.

In brief, our method handles the case if some graph shrinks at the intermediate timestamp, and our theoretical findings hold as long as the life-long pattern is growing.

In other words, for a graph that has some growing stages and some shrinking stages in its life, if its number of nodes follow an overall growing trend, then Theorem 5.2 holds.

To be more specific, let the number of nodes of a temporal graph snapshot at time $t$, $G^{(t)}$, be $n(t)$, then Theorem 5.2 holds if $n(t) \rightarrow \infty$ as $t \rightarrow \infty$. Intuitively, the term $||K_W(W^{(t)}, W'^{(t)}) - K_W(W, W')||$ could be decomposed and bounded by terms that depend on $||W^{(t)} - W||, ||W'^{(t)} - W'||$ and $||X^{(t)} - X||, ||X'^{(t)} - X'||$ (where $X$ is $W$'s signal function and $X^{(t)}, W^{(t)}$ denote the induced graphon and induced graphon corresponds to temporal snapshot $G^{(t)}$). Then, $||W^{(t)} - W||, ||X^{(t)} - X||$ are bounded by terms that depend on the number of nodes $1 / n(t)$, so $||W^{(t)} - W||, ||X^{(t)} - X|| \rightarrow 0$ as $n(t) \rightarrow \infty$.

Moreover, according to the real-world temporal graph dataset benchmark like (https://snap.stanford.edu/data/index.html), especially social graphs and citation graphs, a growing pattern is usually more common.

Q4. Possible negative societal impact

In general, graph deep learning often involves analyzing complex networks, which can include sensitive personal information. If not managed properly, this can lead to privacy violations, especially when dealing with social networks, financial systems, or medical data. Also, if the data used to train graph models contains biases, these biases can be propagated and even amplified by the algorithms.

To the best of the authors’ knowledge, we do not discern specific negative societal impacts caused by this proposed method. However, we will add the above general discussions regarding the risks of graph deep learning in our updated paper.

Thanks very much for the review! We greatly appreciate your acknowledgment of our paper's theoretical analysis. We address your questions in the form of Q&A below.

Moreover, for the length limit, we distillate our answer and provide brief but important answers, more detailed intermediate steps, reasonings, and derivations can be provided during the upcoming author-reviewer discussion phase. All the detailed reply and supplementary experiments will be added to our updated paper.

Q1. Explain using $\mathbf{t}_{enc}$ instead of $t - \bar{t}$ directly? What does $\cos(\Delta t \mathbf{w})$ indicate?

Theoretically, mapping a discrete scalar to a continuous high-dimensional vector space has two unique advantages. First, the embedding vector can capture more complex information like non-linear relationships. Moreover, the embedding does not hurt the original scalar representation, i.e., close time scalars still have similar embedding vectors, and the embedding vectors can also capture periodic behaviors [1]. Second, the embedding vector can jointly work with deep learning models, i.e., embedding vectors for time can be potentially used within any neural architecture’s representation learning process, which is incapable of using discrete time scalars.

In our case, $\mathbf{t}\_{enc}(\Delta t) = \cos(\Delta t \mathbf{w})$ maps a scalar $\Delta t$ to an $d_t-$dimensional vector, whose $i-$th entry is $\cos(\Delta t \cdot \alpha^{-(i-1) / \beta})$. Given two real positive numbers $p, q > 0$ such that $p > q$. As $i \rightarrow d\_t$, the $i-$th entries of $\mathbf{t}\_{enc}(p), \mathbf{t}\_{enc}(q)$ go towards $0$. Since $p < q$, it is likely that $\mathbf{t}\_{enc}(p)$ contains more $0$ entries than $\mathbf{t}\_{enc}(q)$.

Moreover, there are many choices for time encoding. However, most functions require trainable parameters that add extra computational burden. But our $\mathbf{t}\_{enc}(.)$ relies on mathematical properties and only needs fixed parameters ($\alpha, \beta$), which is suitable for kernel-based methods.

Empirically, as shown in Table 4 of our paper, using relative difference $t - \bar{t}$ (3rd line of Table 4) directly yields a lower accuracy than using $\mathbf{t}\_{enc}(t - \bar{t})$, justifying the effectiveness of $\mathbf{t}\_{enc}(.)$.

Reference: [1] Time2Vec: Learning a Vector Representation of Time

Q2. Does the whole framework still work well if the node features are given by the dataset?

The reason we omit node features in formula (1) is that prevalent benchmarks only support time-aware edge features.

But if node features are given, our framework still works. Theoretically, formula (1) can be extended as follows: $h_v(t) = c \cdot \sum_{(u, \bar{t})}[\mathbf{t}\_{enc}(t - \bar{t}) || \mathbf{e}\_{uv}(t) || \mathbf{x}\_{u}(\bar{t})]$, where $\mathbf{x}\_{u}(\bar{t})$ is the node features of $u$ at time $\bar{t}$.

Adding node features will not affect our conclusions and framework, as the construction of our kernel depends on $h_v(t)$ as a general term, instead of its components, such as node / edge features. Moreover, Theorem 5.1 and Theorem 5.2 still hold if node features are given.

Q3. Proof of symmetric and semi-definite properties.

We have proved the symmetric and semi-definite properties of our method. Due to the rebuttal length limit, the detailed proof will be uploaded during the author-reviewer discussion phase via the “official comment” channel.

Q4. However, there should be theoretical proof that the proposed method (1) and (2) indeed converge to the Gaussian process or kernel method, like in [12].

In short, our proposed method naturally enjoys the convergence in [12].

First, we provide some preliminary.

Let $f_{nn}(\mathbf{x}, \theta)$ be the output of a fully connected neural network, with parameters $\theta \in \mathbb{R}^{p}$ and $\mathbf{x} \in \mathbb{R}^d$ as the input.

Then the Neural Tangent Kernel (NTK) corresponds to $f$, and any arbitrary input $\mathbf{x}, \mathbf{x'}$ would be $\nabla_{\theta}f_{nn}(\mathbf{x}, \theta)^ \top \nabla_{\theta}f_{nn}(\mathbf{x}', \theta)$

Suppose that $f_{nn}$ has $L$ layers, the $l-$th layer's width is denoted as $d_l$. Theorem 1 from [12] states that in the limit of $d_1, \dots, d_L \rightarrow \infty$, the NTK converges in probability to a deterministic limiting kernel.

Next, we analyze these terms in our paper's context.

If we apply Theorem 1 from [12], then the NTK (defined in [12]) associating with a deep neural network consisting of $L$ layers applying on temporal node representation $v, v'$ converges, with high probability, to the deterministic limiting kernel $\boldsymbol{\Theta}^{(L)}(G^{(t)}, G^{(t)})\_{vv'} \otimes \text{Id}\_{dL}$, (where $\text{Id}\_{d_L}$ is the identity matrix with dimension $d_L$ and $d_L$ is the width of the $L$-layer MLP), if $d_1, \dots, d_L \rightarrow \infty$, and $\boldsymbol{\Theta}^{(L)}$ is denoted as a scalar kernel. To be clear, the kernel value between two graphs in Equation (13) is defined based on the summation of limiting scalar kernel $\boldsymbol{\Theta}^{(L)}$ taken over all node pairs $v, v'$.

Regarding Gaussian process, based on Proposition 1 from [12], the output of the deep neural network consisting of $L$ layers operating on temporal node representations defined in Section 3.1 converges to Gaussian process of covariance $\boldsymbol{\Sigma}^{(L)}$, with $d_1, \dots, d_{L - 1} \rightarrow \infty$ (as defined in Equation 9).

To better sum up for understanding, our method is constructed based on limiting objects of the NTK convergence theories.

We will add the above derivation to our updated paper.

[Symmetric proof]

We need to prove $K(G^{(t)}, G'^{(t)}) = K(G'^{(t)}, G^{(t)})$.

Given our proposed kernel function, $K(G^{(t)}, G'^{(t)}) = \sum\_{v \in V(t)} \sum\_{v' \in V'(t)} \boldsymbol{\Theta}^{(L)}(G^{(t)}, G'^{(t)})\_{vv'}$, we first write down another equation, where the internal order is flipped, i.e., $K(G'^{(t)}, G^{(t)}) = \sum\_{v' \in V'(t)} \sum\_{v \in V(t)} \boldsymbol{\Theta}^{(L)}(G'^{(t)}, G^{(t)})\_{v'v}$

We first prove that $\boldsymbol{\Theta}^{(l)}(G^{(t)}, G'^{(t)})\_{vv'} = \boldsymbol{\Theta}^{(l)}(G'^{(t)}, G^{(t)})\_{v'v}, \forall l, 1 \leq l \leq L.$

For $l = 1$, we have

$\boldsymbol{\Theta}^{(1)}(G^{(t)}, G'^{(t)})\_{vv'} = \boldsymbol{\Theta}^{(0)}(G^{(t)}, G'^{(t)})\_{vv'} \cdot \dot{\boldsymbol{\Sigma}}^{(1)}(G^{(t)}, G'^{(t)})\_{vv'} + \boldsymbol{\Sigma}^{(1)}(G^{(t)}, G'^{(t)})\_{vv'}$

= $(h\_v(t)^\top h\_{v'}(t)) \cdot \frac{\pi - \arccos(\boldsymbol{\Sigma}^{(0)}(G^{(t)}, G'^{(t)})\_{vv'})}{2\pi} + \frac{\pi - \arccos(\boldsymbol{\Sigma}^{(0)}(G^{(t)}, G'^{(t)})\_{vv'})}{2\pi} + \frac{\sqrt{1 - \boldsymbol{\Sigma}^{(0)}(G^{(t)}, G'^{(t)})^2\_{vv'}}}{2\pi}$

= $(h_v(t)^\top h_{v'}(t)) \cdot \frac{\pi - \arccos(h_v(t)^\top h_{v'}(t))}{2\pi} + \frac{\pi - \arccos(h_v(t)^\top h_{v'}(t))}{2\pi} + \frac{\sqrt{1 - (h_v(t)^\top h_{v'}(t))^2}}{2\pi}$

= $(h_{v'}(t)^\top h_{v}(t)) \cdot \frac{\pi - \arccos(h\_{v'}(t)^\top h\_{v}(t))}{2\pi} + \frac{\pi - \arccos(h\_{v'}(t)^\top h\_{v}(t))}{2\pi} + \frac{\sqrt{1 - (h\_{v'}(t)^\top h\_{v}(t))^2}}{2\pi}$

= $(h_{v'}(t)^\top h\_{v}(t)) \cdot \frac{\pi - \arccos(\boldsymbol{\Sigma}^{(0)}(G'^{(t)}, G^{(t)})_{v'v})}{2\pi} + \frac{\pi - \arccos(\boldsymbol{\Sigma}^{(0)}(G'^{(t)}, G^{(t)})\_{v'v})}{2\pi} + \frac{\sqrt{1 - \boldsymbol{\Sigma}^{(0)}(G'^{(t)}, G^{(t)})^2\_{v'v}}}{2\pi}$

= $\boldsymbol{\Theta}^{(0)}(G'^{(t)}, G^{(t)})\_{v'v} \cdot \dot{\boldsymbol{\Sigma}}^{(1)}(G'^{(t)}, G^{(t)})\_{v'v} + \boldsymbol{\Sigma}^{(1)}(G'^{(t)}, G^{(t)})\_{v'v}$

= $\boldsymbol{\Theta}^{(1)}(G'^{(t)}, G^{(t)})_{v'v}$

Then, suppose $\exists k, 1 \leq k \leq L$, we can have $\boldsymbol{\Theta}^{(k)}(G^{(t)}, G'^{(t)})\_{vv'} = \boldsymbol{\Theta}^{(k)}(G'^{(t)}, G^{(t)})\_{v'v}$

Therefore,

$\boldsymbol{\Theta}^{(k + 1)}(G^{(t)}, G'^{(t)})\_{vv'} = \boldsymbol{\Theta}^{(k)}(G^{(t)}, G'^{(t)})\_{vv'} \cdot \dot{\boldsymbol{\Sigma}}^{(k + 1)}(G^{(t)}, G'^{(t)})\_{vv'} + \boldsymbol{\Sigma}^{(k + 1)}(G^{(t)}, G'^{(t)})\_{vv'}$

= $\boldsymbol{\Theta}^{(k)}(G^{(t)}, G'^{(t)})\_{vv'} \cdot \frac{\pi - \arccos(\boldsymbol{\Sigma}^{(k)}(G^{(t)}, G'^{(t)})\_{vv'})}{2\pi} + \frac{\pi - \arccos(\boldsymbol{\Sigma}^{(k)}(G^{(t)}, G'^{(t)})\_{vv'})}{2\pi} + \frac{\sqrt{1 - \boldsymbol{\Sigma}^{(k)}(G^{(t)}, G'^{(t)})^2\_{vv'}}}{2\pi}$

= $\boldsymbol{\Theta}^{(k)}(G'^{(t)}, G^{(t)})\_{v'v} \cdot \frac{\pi - \arccos(\boldsymbol{\Sigma}^{(k)}(G'^{(t)}, G^{(t)})\_{v'v})}{2\pi} + \frac{\pi - \arccos(\boldsymbol{\Sigma}^{(k)}(G'^{(t)}, G^{(t)})\_{v'v})}{2\pi} + \frac{\sqrt{1 - \boldsymbol{\Sigma}^{(k)}(G'^{(t)}, G^{(t)})^2\_{v'v}}}{2\pi}$

= $\boldsymbol{\Theta}^{(k + 1)}(G^{(t)}, G'^{(t)})\_{vv'}$

Hence, if $\boldsymbol{\Theta}^{(k)}(G^{(t)}, G'^{(t)})\_{vv'} = \boldsymbol{\Theta}^{(k)}(G'^{(t)}, G^{(t)})\_{v'v}$, then $\boldsymbol{\Theta}^{(k + 1)}(G^{(t)}, G'^{(t)})\_{vv'} = \boldsymbol{\Theta}^{(k + 1)}(G'^{(t)}, G^{(t)})\_{v'v}$.

Moreover, we have proven that $\boldsymbol{\Theta}^{(1)}(G^{(t)}, G'^{(t)})\_{vv'} = \boldsymbol{\Theta}^{(1)}(G'^{(t)}, G^{(t)})\_{v'v}$.

Thus, by induction, we have $\boldsymbol{\Theta}^{(l)}(G^{(t)}, G'^{(t)})\_{vv'} = \boldsymbol{\Theta}^{(l)}(G'^{(t)}, G^{(t)})\_{v'v}, \forall l, 1 \leq l \leq L.$

Finally, we have $K(G(t), G'(t)) = \sum_{v \in V(t)} \sum_{v' \in V'(t)} \boldsymbol{\Theta}^{(L)}(G^{(t)}, G'^{(t)})\_{vv'} = \sum_{v' \in V'(t)} \sum_{v \in V(t)} \boldsymbol{\Theta}^{(L)}(G'^{(t)}, G^{(t)})\_{v'v} = K(G'(t), G(t))$

The proof is completed.

[Semi-definite proof]

We need to prove $K(G^{(t)}, G^{(t)}) \geq 0$

Given our proposed kernel function $K(G^{(t)}, G'^{(t)}) = \sum_{v \in V(t)} \sum_{v' \in V(t)} \boldsymbol{\Theta}^{(L)}(G^{(t)}, G'^{(t)})_{vv'}$, we first prove that $\boldsymbol{\Theta}^{(L)}(G^{(t)}, G'^{(t)})\_{vv'} \geq 0, \forall l, 1 \leq l \leq L$.

For $l = 1$,

$\boldsymbol{\Theta}^{(1)}(G^{(t)}, G'^{(t)})\_{vv'} = \boldsymbol{\Theta}^{(0)}(G^{(t)}, G'^{(t)})\_{vv'} \cdot \dot{\boldsymbol{\Sigma}}^{(1)}(G^{(t)}, G'^{(t)})\_{vv'} + \boldsymbol{\Sigma}^{(1)}(G^{(t)}, G'^{(t)})\_{vv'}$

= $(h_v(t)^\top h_{v'}(t)) \cdot \frac{\pi - \arccos(\boldsymbol{\Sigma}^{(0)}(G^{(t)}, G'^{(t)})\_{vv'})}{2\pi} + \frac{\pi - \arccos(\boldsymbol{\Sigma}^{(0)}(G^{(t)}, G'^{(t)})\_{vv'})}{2\pi} + \frac{\sqrt{1 - \boldsymbol{\Sigma}^{(0)}(G^{(t)}, G'^{(t)})^2\_{vv'}}}{2\pi}$

= $(h_v(t)^\top h_{v'}(t)) \cdot \frac{\pi - \arccos(h_v(t)^\top h_{v'}(t))}{2\pi} + \frac{\pi - \arccos(h_v(t)^\top h_{v'}(t))}{2\pi} + \frac{\sqrt{1 - (h_v(t)^\top h_{v'}(t))^2}}{2\pi}$

= $x \cdot \frac{\pi - \arccos(x)}{2\pi} + \frac{\pi - \arccos(x)}{2\pi} + \frac{\sqrt{1 - x^2}}{2\pi} ~(\text{Let } x = (h\_v(t)^\top h\_{v'}(t)))$

As visualized in Figure 1 in the uploaded 1-page PDF file in "general rebuttal", we know that $x \cdot \frac{\pi - \arccos(x)}{2\pi} + \frac{\pi - \arccos(x)}{2\pi} + \frac{\sqrt{1 - x^2}}{2\pi} \geq 0, \forall x \in [-1, 1]$

As $x$ is the input of $\arccos(.)$ function, so it must be constrained to interval $[-1, 1]$. Also, we can achieve the constraint $x \in [-1, 1]$ by normalizing $h_v(t), h_v'(t)$.

Therefore, $\boldsymbol{\Theta}^{(1)}(G^{(t)}, G'^{(t)})\_{vv'} = \boldsymbol{\Theta}^{(0)}(G^{(t)}, G'^{(t)})\_{vv'} \cdot \dot{\boldsymbol{\Sigma}}^{(1)}(G^{(t)}, G'^{(t)})\_{vv'} + \boldsymbol{\Sigma}^{(1)}(G^{(t)}, G'^{(t)})\_{vv'} \geq 0, \forall v, v'$

Suppose $\exists k, 1 \leq k \leq L$, we can have $\boldsymbol{\Theta}^{(k)}(G^{(t)}, G'^{(t)})_{vv'} \geq 0$

Then, we can derive the three following equations that $\boldsymbol{\Theta}^{(k + 1)}(G^{(t)}, G'^{(t)})\_{vv'} = \boldsymbol{\Theta}^{(k)}(G^{(t)}, G'^{(t)})\_{vv'} \cdot \dot{\boldsymbol{\Sigma}}^{(k + 1)}(G^{(t)}, G'^{(t)})\_{vv'} + \boldsymbol{\Sigma}^{(k + 1)}(G^{(t)}, G'^{(t)})\_{vv'}$,

$\dot{\boldsymbol{\Sigma}}^{(k + 1)}(G^{(t)}, G'^{(t)})\_{vv'} = \mathbb{E}_{(a, b) \sim \mathcal{N}(0, \boldsymbol{\Lambda}^{(k+1)}(G^{(t)}, G'^{(t)})\_{vv'})} [\dot{\sigma}(a) \cdot \dot{\sigma}(b)] \geq 0$, as $\dot{\sigma(x)} \geq 0, \forall x$,

and $\boldsymbol{\Sigma}^{(k + 1)}(G^{(t)}, G'^{(t)})\_{vv'} = \mathbb{E}_{(a, b) \sim \mathcal{N}(0, \boldsymbol{\Lambda}^{(k+1)}(G^{(t)}, G'^{(t)})\_{vv'})} [\sigma(a) \cdot \sigma(b)] \geq 0$, as $\sigma(x) \geq 0, \forall x$.

Moreover, as $\boldsymbol{\Theta}^{(k)}(G^{(t)}, G'^{(t)})\_{vv'} \geq 0$, we now have $\boldsymbol{\Theta}^{(k + 1)}(G^{(t)}, G'^{(t)})\_{vv'} \geq 0$.

We also have $\boldsymbol{\Theta}^{(1)}(G^{(t)}, G'^{(t)})\_{vv'} \geq 0$, so by induction, $\boldsymbol{\Theta}^{(l)}(G^{(t)}, G'^{(t)})\_{vv'} \geq 0, \forall l \leq 1 \leq L$.

Finally, we have $K(G^{(t)}, G'^{(t)}) = \sum_{v \in V(t)} \sum_{v' \in V(t)} \boldsymbol{\Theta}^{(L)}(G^{(t)}, G'^{(t)})_{vv'} \geq 0$

The proof is completed.