Continuous Simplicial Neural Networks

We appreciate the reviewer for providing constructive comments on our paper. We'll include the feedback in the camera-ready paper. In the following, we have addressed the reviewer's concerns:

Q1: Firstly, our framework and its theoretical claims are also valid and straightforwardly transferable to the case of weighted structures, not just binarized. Therefore, the incidence matrices can be non-integer in general. Also, the proof of Theorem 4.3 does not need the error matrices to have small norm. Theorem 4.3 holds in general, and having a small norm is only needed in the proof of Corollary 4.4. We will clarify these points in the camera-ready.

We'd also like to clarify the meaning of having a small norm for $\mathbf{E}_k$. From the proof of Corollary 4.4, we observe that being small is w.r.t the scale of the maximum eigenvalue of the Hodge Laplacians, i.e., $\lambda\_{max}(\mathbf{L}\_{k,d})$. Therefore, it is not necessarily considered as being small in an absolute manner.

Finally, regarding the nature of error matrices $\mathbf{E}_k$, we clarify that the removal, addition, and weight changing of edges can be modeled by the additive model $\hat{\mathbf{B}}_1=\mathbf{B}_1+\mathbf{E}_1$.

To simplify the discussion, we start with a simple toy example. Let's consider an undirected three-node graph in which $\mathcal{V}=${$1,2,3$} and $\mathcal{E}=${$(1,2),(1,3)$}. In this case, we have $\mathbf{B}_1=[[1 \\:\\:\\: 1][-1 \\:\\:\\: 0][0 \\:\\:\\: -1]]\in\mathbb{R}^{3\times2}$. Then, one special form of an incidence error matrix corresponding to changes in the current edge weights can be stated as $\mathbf{E}\_1=[[-\epsilon\_{12} \\:\\:\\: -\epsilon\_{13}][\epsilon\_{12} \\:\\:\\: 0][0 \\:\\:\\: \epsilon\_{13}]]\in\mathbb{R}^{3\times2}$. Therefore, after structural perturbation, we have $\hat{\mathbf{B}}_1=\mathbf{B}_1+\mathbf{E}_1=[[1-\epsilon\_{12} \\:\\:\\: 1-\epsilon\_{13}][-1+\epsilon\_{12} \\:\\:\\: 0][0 \\:\\:\\: -1+\epsilon\_{13}]]\in\mathbb{R}^{3\times2}$. So, we have $\\|\mathbf{E}\_1\\|=\sqrt{2}\sqrt{\epsilon^2\_{12}+\epsilon^2\_{13}}$, and, given the scale of $\epsilon\_{12}$ and $\epsilon\_{13}$, this norm can be small or not. As a conclusion, since our definition of simplicial perturbation is not limited to only the addition or removal of edges (or higher-order connections), the norm of the error matrices can take different scales. On the other hand, and in the case of weighted graphs and in connection with the previously mentioned examples, assume that we have the edge weight $e\_{12}=0.0625$ (instead of 1). This implies that $\mathbf{B}_1=[[0.25 \\:\\:\\: 1][-0.25 \\:\\:\\: 0][0 \\:\\:\\: -1]]\in\mathbb{R}^{3\times2}$. Now, for only modeling the removal of this edge, we just need to have $\epsilon\_{12}=0.25$ in $\mathbf{E}_1$, leading to a small norm for $\mathbf{E}\_1$. Therefore, the scale of the norm for $\mathbf{E}$ heavily depends on the scale of edge weights.

Q2: In connection with the response to the previous question, the perturbation of $\mathbf{B}_k$ could model potential erroneous structures where we have the addition of non-existing simplicial connections, the removal of existing simplices, or changes in simplicial connection weights in a principled manner. For further clarification, one can state the toy graph from the response to Q1 as $\mathcal{V}=${$1,2,3$} and $\mathcal{E}=${$(1,2),(1,3),(2,3)$} in which edge $(2,3)$ takes zero value because of disconnectivity. Then, the incidence matrix is stated as $\mathbf{B}_1=[[1 \\:\\:\\: 1 \\:\\:\\: 0][-1 \\:\\:\\: 0 \\:\\:\\: 0][0 \\:\\:\\: -1 \\:\\:\\: 0]]\in\mathbb{R}^{3\times3}$. Now, as an example, for modeling removing edge $(1,2)$ and adding edge $(2,3)$, the error matrix can take the form of $\mathbf{E}_1=[[-1\\:\\:\\: 0 \\:\\:\\: 0][1 \\:\\:\\: 0 \\:\\:\\: 1][0 \\:\\:\\: 0 \\:\\:\\: -1]]\in\mathbb{R}^{3\times3}$. Again, if we consider the case of changing edge weights (not removal or addition) or a weighted graph, the norm of the error matrix can take a very small value. Therefore, the stability of COSIMO against perturbations shows its robustness to such cases.

Q3: We'll clarify in the answer to Q4 that $t$ is indeed a parameter of the model, not a hyperparameter. Previous studies typically analyze the over-smoothing phenomenon by stacking several graph or simplicial layers. Therefore, we provide comprehensive results on the $`high-school`$ node classification dataset averaged over three random seeds in the following table regarding the tradeoff between over-smoothing and accuracy. We observe that the proposed COSIMO can maintain the performance even with an increasing number of layers up to 16 in terms of accuracy and Dirichlet energy. We also observe that COSIMO is more robust against over-smoothing than SCCNN, where extremely smoothness happens after $n_l=4$ and leads to severely reduced performance. The key observation and trend here are also fairly consistent with the results in Figure 3.

Acc:

Dirichlet Energy:

Q4: As stated in Remark 5.6, and due to the differentiability of our framework w.r.t. the simplicial receptive fields, in all the experimental results (except the over-smoothing analysis in Section 6.2), the values of $t$ are learned during the training process, i.e., $t$ are parameters, not hyperparameters. We considered $t$ as a hyperparameter only in the over-smoothing analysis in Section 6.2 to illustrate the effect of different values for $t$ on the over-smoothing phenomenon.

Q5: Eq. (8) implies a definition. In fact, it integrates joint and independent dynamics that act as source terms; this has been shown to improve robustness to over-smoothing and stability as outlined in [1]. The analytical form of the solution to the set of PDEs in eqs. (5), (6), and (7) take the form of eq. (9), which can be stated in a general form as eq. (8). Note that $\mathbf{x}\_{k,d}(t\_d) + \mathbf{x}\_{k,u}(t\_u)$ in eq.(8) is actually the general symbolic form of $\mathbf{x}\_{k}(t\_d,t\_u)$ in eq. (9). We will modify and clarify these points in the camera-ready version.

---

References:

[1] A. Han, et al, “From continuous dynamics to graph neural networks: Neural diffusion and beyond,” TMLR 2025

We appreciate the reviewer for providing constructive comments on our paper. We'll include the feedback in the camera-ready paper. In the following, we have addressed the reviewer's concerns:

W1 and Q1: First, we emphasize that performing EVD on the Hodge Laplacians is done once and as a preprocessing step, so they can be precomputed. Additionally, as will be presented in the following, one can rely on a small subset of eigenvalue-eigenvector pairs, and then, the computational complexity of performing EVD is reduced from $\mathcal{O}(N^3)$ to $\mathcal{O}(N^2K)$. In our framework, one can select a proper $K$ using supervised (cross-validation) or unsupervised methodologies. We used the cross-validation approach in the paper. Regarding the unsupervised methodologies, one can exploit the following alternative strategies:

1. Eigengap heuristic [1]: The number $K$ is often found using the eigengap heuristic:

$K = \text{argmax}\_{i}{(\lambda\_{i+1}-\lambda\_i)}$.

A large gap indicates a natural cutoff point.

2. Energy-based Criterion [2]: The Laplacian EVD can be viewed like PCA, such that one can choose enough eigenvectors to explain a certain percentage of spectral energy. Therefore, the total spectral energy is defined as:

$E\_{total}:=\sum\_{i=1}^{n}{\lambda\_i}$.

One can thus choose the smallest $K$ such that $\frac{\sum_{i=1}^{K}{\lambda_i}}{E_{total}}\ge\eta$, where $\eta\in[0,1]$ is a threshold like 90%.

3. Spectral Gap Ratio or Knee Detection [3]: Instead of just picking the largest gap, use a normalized eigengap ratio:

$\text{gap ratio}(K):=\frac{\lambda_{K+1}-\lambda_{K}}{\lambda_{K}}$.

This method accounts for the scaling of eigenvalues and can find relative jumps.

4. Spectral Entropy-Information-Theoretic Criteria [4]: The spectral entropy can be defined as:

$H:=-\sum\_{i=1}^{n}{p\_i\log{p\_i}}; \text{ where } p\_i=\frac{\lambda\_i}{\sum\_{j}{\lambda\_j}}$.

Therefore, one can choose $K$ where the entropy contribution of the next eigenvalues becomes negligible.

Apart from introducing and listing the relevant methods, we have also experimentally applied these strategies on a subset of the $`high-school`$ node classification dataset in the following table, accompanied by the accuracy performance. Note the supervised methodology picked $\sim5.1\\%$ of eigenvalues for $**L**_1$, which is consistent with the unsupervised methods ($\sim5.5\\%$).

Table: Selected optimal values of eigenvalue-eigenvector $(K)$ pairs from $\sim$ 5800 ones for $\mathbf{L}_1$ by different relevant methods.

Regarding the effect of varying $K$ on the performance, we provide comprehensive experimental results on a subset of the $`high-school`$ node classification dataset in the following table. From these results, we can conclude that choosing the best $K$ heavily relies on the specific data distribution, but generally, only a small fraction of the eigenvalues is enough to reach acceptable performance. We also observe that increasing $K$ does not necessarily lead to better performance, as high-frequency components might contain more noisy or non-relevant content [5].

Table: Performance (accuracy) results across a selected set of eigenvalue-eigenvector $(K)$ pairs from $\sim$ 300 ones for $\mathbf{L}_0$.

---

References:

[1] J. Shi, et al, “Normalized cuts and image segmentation,” IEEE TPAMI, 2000

[2] M. Belkin, et al, “Laplacian eigenmaps for dimensionality reduction and data representation,” Neural computation, 2003

[3] I. M. Johnstone, “On the distribution of the largest eigenvalue in principal components analysis,” The Annals of statistics, 2001

[4] M. De Domenico, et al, “Spectral entropies as information-theoretic tools for complex network comparison,” Physical Review X, 2016

[5] Dong, X., et al. Learning Laplacian matrix in smooth graph signal representations. IEEE TSP 2016.

We appreciate the reviewer for providing constructive comments on our paper. We'll include the feedback in the camera-ready paper. In the following, we have addressed the reviewer's concerns:

W1: To complement the experimental analysis, we compare COSIMO against GRAND and Graph-coupled oscillator networks (GraphCON) [1] in the following table. We observe the superiority of COSIMO over GNNs, which only rely on the node space. This is a general observation in the comparison of SNNs vs. GNNs when higher-order information is available.

W2: We provide comprehensive results on the $`high-school`$ node classification dataset averaged over three random seeds in the following table regarding the tradeoff between over-smoothing and accuracy. We observe that the proposed COSIMO can maintain the performance even with an increasing number of layers up to 16 in terms of accuracy and Dirichlet energy. We also observe that COSIMO is more robust against over-smoothing than SCCNN, where extremely smoothness happens after $n_l=4$ and leads to severely reduced performance. The key observation and trend here are also fairly consistent with the results in Figure 3.

Acc:

Dirichlet Energy:

W3: We have repeated the experiment in Section 6.3, including comparisons with the baseline methods SCCNN and SAN, and provided the results in the following table in terms of prediction error. We observe that the baselines are vulnerable in moderate to high amounts of noise (-5 or 0 db SNRs). In comparison, COSIMO has superior and more robust performance compared to the baselines even in the low-noisy regime (SNRs of 10 or 20 db).

SNR$_2$=-5:

SNR$_2$=0:

SNR$_2$=10:

Q1: Based on the provided results regarding over-smoothing, stability, and performance, we argue that relying on simplicial complexes is meaningful when higher-order information is available. This is because SNNs, like COSIMO, can benefit from the higher-order dynamics. One practical example is the $`high-school`$ dataset, where friendship groups of multiple people (maximum cliques) are modeled as simplices. The description of the dataset has been explicitly mentioned on the official website as follows:

"This is a temporal higher-order network dataset, which here means a sequence of timestamped simplices where each simplex is a set of nodes. The dataset is constructed from interactions recorded by wearable sensors by people at a high school. The sensors record interactions at a resolution of 20 seconds (recording all interactions from the previous 20 seconds). Nodes are the people, and simplices are maximal cliques of interacting individuals from an interval."

It is clear from our experimental results that relying on the typical graph structure (only node space) cannot leverage the information dynamics lying on the simplices in this dataset, so the performance of classical GNNs is not competitive against SNNs.

To explore more diverse real-world applications from computer graphics to drug discovery and biotechnology, we kindly refer the reviewer to this recent ICML position paper [2].

---

References:

[1] T. K. Rusch, et al, “Graph-coupled oscillator networks,” ICML 2022

[2] Papamarkou, et al. "Position: Topological deep learning is the new frontier for relational learning". ICML 2024.

We appreciate the reviewer for providing constructive comments on our paper. We'll include the feedback in the camera-ready paper. In the following, we have addressed the reviewer's concerns:

Q1: Firstly, we clarify that we used cross-validation for tuning the hyperparameter. For experimental results on $`synthetic`$ and $`ocean-drifts`$, we followed the experimental settings from reference papers [1,2]. In addition to adding a comprehensive table of details like the following table for the camera-ready, we will publicly release the codes (which have already been uploaded alongside our submission) and refer to them to make these settings as clear as possible.

Q2: In order to make the connection between the proposed COSIMO and reference Graph-PDE counterparts [3,4,10], we first highlight the generalization aspect of simplicial complexes over graphs. For instance, let's consider a simple case of a simplicial complex containing nodes, edges, and triangles. In this case, by relying on graph structure, we neglect edge and triangle dynamics. Therefore, we only deal with $\mathbf{L}_0$ (in the nodes space with $k=0$), and equations (6) and (7) in our framework are canceled. On the other hand, since $k=0$, only the lower node space is present and equation (5) in our framework turns to $\frac{\partial \mathbf{x}(t)}{\partial t}=-\mathbf{L}_0\mathbf{x}(t)$. Then, by assuming the usage of the normalized version of this Laplacian, i.e., $\hat{\mathbf{L}}_0=\mathbf{I}-\mathbf{A}$ with $\mathbf{A}$ being the normalized adjacency matrix, this equation turns to $\frac{\partial \mathbf{x}(t)}{\partial t}=(\mathbf{A}-\mathbf{I})\mathbf{x}(t)$, which is equivalent to eq. (2) in GRAND [4] and eq. (7) in CGNN [3] frameworks in the isotropic linear cases. Therefore, COSIMO generalizes the CGNN and GRAND models by considering edge and triangular dynamics. Apart from these theoretical/architectural differences, we have also experimentally compared COSIMO against GRAND [4] and Graph-coupled oscillator networks (GraphCON) [10] in the following table, showcasing the superior performance by relying on higher-order dynamics.

Q3: Apart from the outlined possible limitations in the conclusion section, we also point out the case of i) eventually facing over-smoothing in Figure 3, and ii) an extremely noisy regime (for example, $SNR_1=-5$ and $SNR_2=-5$ in Figure 4). To mitigate these limitations in these cases, it has been shown that equipping the underlying set of PDEs with source terms [5] can effectively improve the stability and robustness to over-smoothing. Therefore, the corresponding equations of our framework, i.e., independent and joint dynamics, turn to:

$\frac{\partial \mathbf{x}\_{k,d}(t\_d)}{\partial t\_d} = -\mathbf{L}\_{k,d}\mathbf{x}\_{k,d}(t\_d)\rightarrow\frac{\partial \mathbf{x}\_{k,d}(t\_d)}{\partial t\_d} = -\mathbf{L}\_{k,d}\mathbf{x}\_{k,d}(t\_d)+\mathbf{E}\_d,$

$\frac{\partial \mathbf{x}\_{k,u}(t\_u)}{\partial t\_u} = -\mathbf{L}\_{k,u}\mathbf{x}\_{k,u}(t\_u)\rightarrow\frac{\partial \mathbf{x}\_{k,u}(t\_u)}{\partial t\_u} = -\mathbf{L}\_{k,u}\mathbf{x}\_{k,u}(t\_u)+\mathbf{E}\_u,$

$\frac{\partial \mathbf{x}\_{k}(t\_d, t\_u)}{\partial t\_d} + \frac{\partial \mathbf{x}\_{k}(t\_d, t\_u)}{\partial t\_u} = -\mathbf{L}\_{k,d}\mathbf{x}\_{k}(t\_d, \infty) - \mathbf{L}\_{k,u}\mathbf{x}\_{k}(\infty, t\_u)\rightarrow\frac{\partial \mathbf{x}\_{k}(t\_d, t\_u)}{\partial t\_d} + \frac{\partial \mathbf{x}\_{k}(t\_d, t\_u)}{\partial t\_u} = -\mathbf{L}\_{k,d}\mathbf{x}\_{k}(t\_d, \infty) - \mathbf{L}\_{k,u}\mathbf{x}\_{k}(\infty, t\_u)+\mathbf{E}\_{d,u},$

where $\mathbf{E}\_d$, $\mathbf{E}\_u$, and $\mathbf{E}\_{d,u}$ are source terms. Although this modification might bring the above-mentioned advantages, the solutions pose more computational burden, as detailed in different frameworks in [5].

Q4: As stated in Remark 5.6, and due to the differentiability of our framework w.r.t. the simplicial receptive fields, in all the experimental results (except the over-smoothing analysis in Section 6.2), the values of $t$ are learned during the training process, i.e., $t$ are parameters, not hyperparameters. We considered $t$ as hyperparameters only in the over-smoothing analysis in Section 6.2 to illustrate the effect of different values for $t$ on the over-smoothing phenomenon. In general, we use cross-validation over a range of possible values to tune the hyperparameters of COSIMO. Regarding the hyperparameter $M$, we also exploited a multi-branch architecture in trajectory prediction tasks in which we have already provided its sensitivity and accuracy trade-off in the Appendix, i.e., Figure 7, by using a cross-validation scheme on the corresponding range for $M$. Furthermore, one can select a proper $K$ using supervised (cross-validation) or unsupervised methodologies. We used the cross-validation approach in the paper. Regarding the unsupervised methodologies, one can exploit the following alternative strategies:

1. Eigengap heuristic [6]: The number $K$ is often found using the eigengap heuristic:

$K = \text{argmax}\_{i}{(\lambda\_{i+1}-\lambda\_i)}$.

A large gap indicates a natural cutoff point.

2. Energy-based Criterion [7]: The Laplacian EVD can be viewed like PCA, such that one can choose enough eigenvectors to explain a certain percentage of spectral energy. Therefore, the total spectral energy is defined as:

$E\_{total}:=\sum\_{i=1}^{n}{\lambda\_i}$.

One can thus choose the smallest $K$ such that $\frac{\sum_{i=1}^{K}{\lambda_i}}{E_{total}}\ge\eta$, where $\eta\in[0,1]$ is a threshold like 90%.

3. Spectral Gap Ratio or Knee Detection [8]: Instead of just picking the largest gap, use a normalized eigengap ratio:

$\text{gap ratio}(K):=\frac{\lambda_{K+1}-\lambda_{K}}{\lambda_{K}}$.

This method accounts for the scaling of eigenvalues and can find relative jumps.

4. Spectral Entropy-Information-Theoretic Criteria [9]: The spectral entropy can be defined as:

$H:=-\sum\_{i=1}^{n}{p\_i\log{p\_i}}; \text{ where } p\_i=\frac{\lambda\_i}{\sum\_{j}{\lambda\_j}}$.

Therefore, one can choose $K$ where the entropy contribution of the next eigenvalues becomes negligible.

Apart from introducing and listing the relevant methods, we have also experimentally applied these strategies on a subset of the $`high-school`$ node classification dataset in the following table, accompanied by the accuracy performance. Note the supervised methodology picked $\sim5.1\\%$ of eigenvalues for $**L**_1$, which is consistent with the unsupervised methods ($\sim5.5\\%$).

Table: Selected optimal values of eigenvalue-eigenvector $(K)$ pairs from $\sim$ 5800 ones for $\mathbf{L}_1$ by different relevant methods.

In conclusion and by observing results on various datasets, choosing 5-10$\\%$ of eigenvalue-eigenvector pairs and number of branches ($M$) from 1-3 seems to provide a fair compromise of performance on downstream tasks.

---

References:

[1] M. Roddenberry, et al, “Principled simplicial neural networks for trajectory prediction,” ICML 2021

[2] M. Yang, et al, “Convolutional learning on simplicial complexes,” arXiv 2025

[3] L.-P. Xhonneux, et al, “Continuous graph neural networks,” ICML 2020

[4] B. Chamberlain, et al, “Grand: Graph neural diffusion,” ICML 2021

[5] A. Han, et al, “From continuous dynamics to graph neural networks: Neural diffusion and beyond,” TMLR 2024

[6] J. Shi, et al, “Normalized cuts and image segmentation,” IEEE TPAMI, 2000

[7] M. Belkin, et al, “Laplacian eigenmaps for dimensionality reduction and data representation,” Neural computation, 2003

[8] I. M. Johnstone, “On the distribution of the largest eigenvalue in principal components analysis,” The Annals of statistics, 2001

[9] M. De Domenico, et al, “Spectral entropies as information-theoretic tools for complex network comparison,” Physical Review X, 2016

[10] T. K. Rusch, et al, “Graph-coupled oscillator networks,” ICML 2022