Continuous Product Graph Neural Networks

Authors: einizade  aref
télécom  ltci
  paris
fragkiskos d malliaros
jhony h giraldo ltci
télécom  paris
Venue: NeurIPS_cc_2024_Conference
Published: 2024-05-02 17:44:20

Abstract: Processing multidomain data defined on multiple graphs holds significant potential in various practical applications in computer science. However, current methods are mostly limited to discrete graph filtering operations. Tensorial partial differential equations on graphs (TPDEGs) provide a principled framework for modeling structured data across multiple interacting graphs, addressing the limitations of the existing discrete methodologies. In this paper, we introduce Continuous Product Graph Neural Networks (CITRUS) that emerge as a natural solution to the TPDEG. CITRUS leverages the separability of continuous heat kernels from Cartesian graph products to efficiently implement graph spectral decomposition. We conduct thorough theoretical analyses of the stability and over-smoothing properties of CITRUS in response to domain-specific graph perturbations and graph spectra effects on the performance. We evaluate CITRUS on well-known traffic and weather spatiotemporal forecasting datasets, demonstrating superior performance over existing approaches. The implementation codes are available at https://github.com/ArefEinizade2/CITRUS.

Main paper: 
 introduction: 
 Multidomain (tensorial) data defined on multiple interacting graphs [1] [2] [3] , referred to as multidomain graph data in this paper, extend the traditional graph machine learning paradigm, which typically deals with single graphs [2, 4] . Tensors, which are multi-dimensional generalizations of matrices (order-2 tensors), appear in various fields like hyperspectral image processing [5] , video processing [6] , recommendation systems [7] , spatiotemporal analysis [8] , and brain signal processing [9] . Despite the importance of these applications, learning from multidomain graph data has received little attention in the existing literature [2, 10] . Therefore, developing graph-learning strategies for these tensorial data structures holds significant promise for various practical applications.The main challenge for learning from multidomain graph data is creating efficient frameworks that model joint interactions across domain-specific graphs [10, 11] . Previous work in this area has utilized discrete graph filtering operations in product graphs (PGs) [10, 12] from the field of graph signal processing (GSP) [13] . However, these methods inherit the well-known issues of oversmoothing and over-squashing from regular graph neural networks (GNNs) [14] [15] [16] , which restricts the graph's receptive field and hinders long-range interactions [17] . Additionally, these methods often require computationally intensive grid searches to tune hyperparameters and are typically limited to two-domain graph data, such as spatial and temporal dimensions [10, 12, 18, 19] .In regular non-tensorial GNNs, continuous GNNs (CGNNs) emerge as a solution to over-smoothing and over-squashing [20] . CGNNs are neural network solutions to partial differential equations (PDEs) on graphs [20, 21] , like the heat or wave equation. The solution to these PDEs introduces the exponential graph filter as a continuous infinity-order generalization of the typical discrete graph convolutional filters [20] . Due to the differentiability of exponential graph filters w.r.t. the graph receptive fields, CGNNs can benefit from both global and local message passing by adaptively learning graph neighborhoods, alleviating the limitations of discrete GNNs [17, 20, 22] . Despite these advantages, CGNNs mostly rely on a single graph and lack a principled framework for learning joint multi-graph interactions.In this paper, we introduce tensorial PDE on graphs (TPDEGs) to address the limitations of existing PG-based GNN and CGNN frameworks. TPDEGs provide a principled framework for modeling multidomain data residing on multiple interacting graphs that form a Cartesian product graph heat kernel. We then propose Continuous Product Graph Neural Networks (CITRUS) as a continuous solution to TPDEGs. We efficiently implement CITRUS using a small subset of eigenvalue decompositions (EVDs) from the factor graphs. Additionally, we conduct theoretical and experimental analyses of the stability and over-smoothing properties of CITRUS. We evaluate our proposed model on spatiotemporal forecasting tasks, demonstrating that CITRUS achieves state-of-the-art performance compared to previous Temporal-Then-Spatial (TTS), Spatial-Then-Temporal (STT), and Temporal-and-Spatial (T&S) frameworks. Our main contributions are summarized as follows:• We introduce TPDEGs to handle multidomain graph data. This leads to our proposal of a continuous graph filtering operation in Cartesian product spaces, called CITRUS, which naturally emerges as the solution to TPDEGs.• We conduct extensive theoretical and experimental analyses to evaluate the stability and oversmoothing properties of CITRUS.• We test CITRUS on spatiotemporal forecasting tasks, demonstrating that our model achieves state-of-the-art performance.Question: Do the main claims made in the abstract and introduction accurately reflect the paper's contributions and scope?Answer: [Yes] Justification: We have tried to include all the paper's contributions and scope in the abstract and introduction.Guidelines:• The answer NA means that the abstract and introduction do not include the claims made in the paper. • The abstract and/or introduction should clearly state the claims made, including the contributions made in the paper and important assumptions and limitations. A No or NA answer to this question will not be perceived well by the reviewers. • The claims made should match theoretical and experimental results, and reflect how much the results can be expected to generalize to other settings. • It is fine to include aspirational goals as motivation as long as it is clear that these goals are not attained by the paper. 
related work: 
 We divide the related work into two parts covering the fundamental and applicability aspects.Fundamentals. The study of PGs is a well-established field in mathematics and signal processing [23, 24] . However, to the best of our knowledge, the graph-time convolutional neural network (GTCNN) [10, 18] model is the only GNN-based framework addressing neural networks in PGs. GTCNN defines a discrete filtering operation in PGs to jointly process spatial and temporal data for forecasting applications. Due to the discrete nature of the graph polynomial filters in GTCNN, an expensive grid search is required to correctly tune the graph filter orders, leading to significant computational overhead. Additionally, this discrete aspect restricts node-wise receptive fields, making long-range communication challenging [17, 20] . Furthermore, GTCNN is limited to two factor graphs (space and time), limiting its applicability to general multidomain graph data. In contrast, CITRUS overcomes these limitations by i) employing continuous graph filtering, and ii) providing a general framework for any number of factor graphs.Applicability. We explore spatiotemporal analysis as a specific application of PG processing. Early works in GSP for spatiotemporal forecasting, such as the graph polynomial vector auto-regressive moving average (GP-VARMA) and GP-VAR models [25] , do not utilize PGs. Modern GNNbased architectures can be broadly classified into TTS and STT frameworks. TTS networks, due to their sequential nature, often encounter propagation information bottlenecks during training [8] . Conversely, STT architectures enhance node representations first and then integrate temporal information, but the use of recurrent neural networks (RNNs) in the aggregation stages leads to high computational complexity. Additionally, there are T&S frameworks designed for simultaneous learning of spatiotemporal representations [8, 26] , where the fusion of temporal and spatial modules is crucial. However, similar to STT networks, T&S frameworks also face challenges with computational complexity [8] and their discretized block-wise nature, which can limit adaptive graph receptive fields [20, 26] .CITRUS leverages joint multidomain learning by exploiting continuous separable heat kernels as factor-based functions defined on a Cartesian product graph. The continuity and differentiability of these filters allow the graph receptive fields to be learned adaptively during the training process, eliminating the need for a grid search. Additionally, CITRUS maintains low computational complexity by relying on the spectral decomposition of the factor graphs [17, 22] . Furthermore, unlike non-graph approaches, the number of learnable parameters in CITRUS is independent of the factor graphs, ensuring scalability.Question: Does the paper discuss the limitations of the work performed by the authors?Answer: [Yes]Justification: We have discussed it at the end of the Conclusion and Limitations section as a separate paragraph for planning future work.Guidelines:• The answer NA means that the paper has no limitation while the answer No means that the paper has limitations, but those are not discussed in the paper. • The authors are encouraged to create a separate "Limitations" section in their paper.• The paper should point out any strong assumptions and how robust the results are to violations of these assumptions (e.g., independence assumptions, noiseless settings, model well-specification, asymptotic approximations only holding locally). The authors should reflect on how these assumptions might be violated in practice and what the implications would be. • The authors should reflect on the scope of the claims made, e.g., if the approach was only tested on a few datasets or with a few runs. In general, empirical results often depend on implicit assumptions, which should be articulated. • The authors should reflect on the factors that influence the performance of the approach.For example, a facial recognition algorithm may perform poorly when image resolution is low or images are taken in low lighting. Or a speech-to-text system might not be used reliably to provide closed captions for online lectures because it fails to handle technical jargon. • The authors should discuss the computational efficiency of the proposed algorithms and how they scale with dataset size. • If applicable, the authors should discuss possible limitations of their approach to address problems of privacy and fairness. • While the authors might fear that complete honesty about limitations might be used by reviewers as grounds for rejection, a worse outcome might be that reviewers discover limitations that aren't acknowledged in the paper. The authors should use their best judgment and recognize that individual actions in favor of transparency play an important role in developing norms that preserve the integrity of the community. Reviewers will be specifically instructed to not penalize honesty concerning limitations. 
methodology and theoretical analysis: 
 Preliminaries and notation. An undirected and weighted graph G with N vertices can be stated as G = {V, E, A}. V and E denote the sets of the vertices and edges, respectively. A ∈ R N ×N is the adjacency matrix with {A ij = A ji ≥ 0} N i,j=1 representing the connection between vertices, and {A ii = 0} N i=1 states that there are no self-loops in G. The graph Laplacian L ∈ R N ×N is defined as L = D -A, where D = diag(A1) is the diagonal degree matrix and 1 is the all-one vector. A graph signal is a function that maps a set of nodes to the real values x : V → R, and therefore we can represent it as x = [x 1 , . . . , x N ] ⊤ . For a vector a = [a 1 , . . . , a N ] ⊤ ∈ R N ×1 , e a := [e a1 , . . . , e a N ] ⊤ . The vectorization operator is shown as vec(.). Using the Kronecker product ⊗, Kronecker sum (aka Cartesian product) ⊕, and the Laplacian factors {L p ∈ R Np×Np } P p=1 , we have the following definitions:EQUATION) where the Laplacian matrix of the Cartesian product between two factor graphs can be stated as:EQUATIONwith I n the identity matrix of size n. See Figure 1a for an example of a Cartesian graph product between three factor graphs. A similar relationship to (2) also holds for adjacency matrices [24, 27] . We define a D-dimensional tensor as U ∈ R N1×...×N D . The matrixU (i) ∈ R Ni×[ D r=1,r̸ =i Nr]is the i-th mode matricization of U [28] . The mode-i tensorial multiplication of a tensor U by a matrixX ∈ R m×Ni is denoted by G = U × i X, where G ∈ R N1×...×Ni-1×m×Ni+1×...×N D [28].For further information about the tensor matricization and mode-n product, please refer to [28, 29] (especially Sections 2.4 and 2.5 in [28] ). All the proofs of theorems, propositions, and lemmas of this paper are provided in the Appendix A.We state the generative PDE for CITRUS as follows:Definition 3.1 (Tensorial PDE on graphs (TPDEG)). Let Ũt ∈ R N1×N2×...×N P be a multidomain tensor whose elements are functions of time t, and let {L p ∈ R Np×Np } P p=1 be the domain-specific factor Laplacians. We can define the TPDEG as follows:EQUATIONTheorem 3.2. Let Ũ0 be the the initial value of Ũt . The solution to the TPDEG in (3) is given by: Ũt = Ũ0 × 1 e -tL1 × 2 e -tL2 × 3 . . . × P e -tL P .Using Theorem 3.2, we define the core function of CITRUS as follows:EQUATIONwhere U l ∈ R N1×...×N P ×F l is the input tensor, W l ∈ R F l ×F l+1 is a matrix of learnable parameters, t l is the learnable graph receptive field, and f (U l ) ∈ R N1×...×N P ×F l+1 is the output tensor. Figure 1b illustrates a high-level description of f (U l ). The next proposition formulates (5) as a graph convolution defined on a product graph.Proposition 3.3. The core function of CITRUS in (5) can be rewritten as:EQUATIONwhere L ⋄ :=↓ ⊕ P p=1 L p is the Laplacian of the Cartesian product graph. Proposition 3.3 is the main building block for implementing CITRUS. More precisely, we use the spectral decompositions of the factor graphs {L p = V p Λ p V ⊤ p } P p=1 and product graph [23] . Let K p ≤ N p be the number of selected eigenvalue-eigenvector pairs of the p-th factor Laplacian. When K p = N p , it can be shown [17, 22] that we can rewrite (6) as follows:{L ⋄ = V ⋄ Λ ⋄ V ⊤ ⋄ }, where V ⋄ =↓ ⊗ P p=1 V p and Λ ⋄ =↓ ⊕ P p=1 Λ pEQUATIONwithEQUATIONwhereλ (Kp) p ∈ R Kp×1 and V (Kp) p∈ R Np×Kp are the first K p ≤ N p selected eigenvalues and eigenvectors of L p based on largest eigenvalue magnitudes, respectively, and ⊙ is the pointwise multiplication operation. Finally, we can define the output of the l-th layer of CITRUS as U l+1 (P +1) = σ f (U l ) (P +1) , where σ(•) is some proper activation function. P p=1 N 3 p ]) in (7) by relying upon the properties of product graphs since we perform EVD on each factor graph independently. Besides, if we rely only on the K p most important eigenvector-eigenvalue pairs [17] of each factor graph, we can reduce the complexity of the spectral decomposition to O(N 2 p K p ), obtaining a general complexity of O([P p=1 N 2 p K p ]).We study the stability of CITRUS against possible perturbations in the factor graphs. First, as defined in the literature [30] , we model the perturbation to the p-th graph as an addition of an error matrix E p (with upper bound ε p for any matrix norm |||•|||) to the adjacency matrix A p as:EQUATIONProposition 3.6. Let A and Ã be the true and perturbed adjacency matrix of Cartesian product graphs with the perturbation model for each factor graph described as in (9) . Then, it holds that:EQUATIONwhere the perturbation matrix E also follows the Cartesian structure E = ⊕ P p=1 E p .Finally, let φ(u, t) and φ(u, t) be the true and perturbed output of the CITRUS model with normalized Laplacian L⋄ being responsible for generating U t(P +1) in a heat flow PDE with normalized Laplacian L⋄ as ∂φ(u,t) ∂t = -L⋄ φ(u, t) in (6) [30] , respectively. Then, the following theorem states that the integrated error bound on the stability properties of CITRUS is also separable when dealing with Cartesian product graphs: Theorem 3.7. Consider a Cartesian product graph (without isolated nodes) with the normalized true and perturbed Laplacians L⋄ and L⋄ in (6) . Also, assume we have the perturbation model Ãp = A p + E p with factor error bounds {|||E p ||| ≤ ε p } P p=1 in Proposition 3.6. Then, the stability bound on the true and perturbed outputs of CITRUS , i.e., φ(u, t) and φ(u, t), respectively, can be described by the summation of the factor stability bounds as:EQUATION)Theorem 3.7 states that the solution of the TPDEG ( 4) is robust w.r.t. the scale of the factor graph perturbations. This is a desired property when dealing with erroneous factor graphs, and will be numerically validated in Section 4.1.There is a prevalent notion of describing over-smoothing in GNNs based on decaying Dirichlet energy against increasing the number of layers. This has been theoretically analyzed in related literature [20, 31, 32] and can be formulated [20] for the output of a continuous GNN. Precisely, it can be defined on the GNN's output U t = [u(t) 1 , . . . , u(t) N ] ⊤ , with u(t) i being the feature vector of the i-th node with the degree deg i , as:EQUATIONwhere E(U) := 1 2EQUATIONInspired by (13) , we extend the definition of Dirichlet energy to the tensorial case as follows:Definition 3.8. (Tensorial Dirichlet energy). By considering the normalized factor Laplacians { Lp } P p=1 , we define the Tensorial Dirichlet energy for a tensor U ∈ R N1×...×N P ×F as:EQUATIONwhereŨf (p) = U[:, . . . , :, f ] (p) ∈ R Np× P i̸ =p Ni , i.e., the p-th mode matricization of the f -th slice of U in its (P + 1)-th dimension.Based on the previous research in the GSP literature [27, 33, 34] , it follows that E(U) in ( 14) can be rewritten as p=1 Np] being the identity matrix) also holds as the following:E(U) = tr( Ũ⊤ L Ũ), where Ũ ∈ R [Π P p=1 Np]×F , Ũ[:, f ] = vec(EQUATIONBesides, the spectrum of L is spanned across the interval of [0, 2].We observe in ( 14) that the kernel separability of Cartesian heat kernels also leads to the separability of analyzing over-smoothing in different modes of the data at hand. This is useful in cases of facing factor graphs with different characteristics. Next, we formally analyze over-smoothing in CITRUS.In our case, for the l-th layer with H l hidden MLP layers, non-linear activations σ(•) and learnable weight matrices {W lh } H l h=1 , we define:EQUATIONwith X 0 being the initial node feature matrix, f l (•) is a generalized form of (6) , and e -L(t) :=↓ ⊗ P p=1 e -t (p) Lp P , where t (i) is the receptive field of the i-th factor graph. Then, the next theorem describes the over-smoothing criteria.Theorem 3.10. For the product Laplacian L := 1 P ⊕ P p=1 Lp with domain-specific receptive field t (i) for the i-th factor graph and the activations σ(•) in (16) being ReLU or Leaky ReLU, we have:EQUATIONwhere s := sup l∈N+ s l and s l := H l h=1 s lh with s lh being the square of maximum singular value of W ⊤ lh . Besides, λ = λ (m) and t = t (m) with m = arg min i t (i) λ (i) , where λ (i) is the smallest non-zero eigenvalue of Li .Corollary 3.11. When l → ∞, E(X l ) exponentially converges to 0, when:EQUATIONTheorem 3.10 shows that the factor graph with the smallest non-zero eigenvalue (spectral gap) multiplied by its receptive field dominates the overall over-smoothing. The less the spectral gap, the less probability the factor graph is connected [14] . So, we can focus on the (much smaller) factor graphs instead of the resulting massive product graph. These theoretical findings are experimentally validated in Section 4.2.Question: For each theoretical result, does the paper provide the full set of assumptions and a complete (and correct) proof?Question: Does the paper provide open access to the data and code, with sufficient instructions to faithfully reproduce the main experimental results, as described in supplemental material?Answer: [Yes] Justification: We have uploaded the reproducible codes alongside our paper and a Readme file to guide the main results' reproductions. Besides, all of the real datasets we used are open-access benchmark.Guidelines:• The answer NA means that paper does not include experiments requiring code. • Providing as much information as possible in supplemental material (appended to the paper) is recommended, but including URLs to data and code is permitted. 
experimental results: 
 In this section, we experimentally validate the theoretical discussions regarding the stability and over-smoothing analysis and use CITRUS in real-world spatiotemporal forecasting tasks. Similarly, we conduct ablation studies regarding CITRUS in spatiotemporal forecasting. We compare CIT-RUS against several state-of-the-art methods in forecasting. The results on the case of more than two factor graphs and descriptions of the state-of-the-art are provided in Appendix B and C, respectively, alongside the additional theoretical and experimental discussions in Sections D-I. The implementation codes are available at https://github.com/ArefEinizade2/CITRUS. We generate a 600-node product graph consisting of two connected Erdős-Rényi (ER) factor graphs with N 1 = 20, N 2 = 30, and edge probabilitiesp (1) ER = p (2)ER = 0.1 for node regression. We utilize 15% of the nodes in the product graph as the test set, 15% of the remaining nodes as validation, and the rest of the nodes for training. We use a three-layer generative model to generate the outputs based on (6) . The continuous parameters for the factor graphs are set as t (1) = 2 and t (2) = 3. The initial product graph signal X (0) ∈ R N ×F0 is generated from the normal distribution with F 0 = 6. Each layer's MLP of the generation model has only one layer with the number of hidden units {5, 4, 2}. We add noise to the factor adjacency matrices with signal-to-noise ratios (SNRs) (in terms of db) of {∞, 20, 10, 0, -10}. Therefore, we construct a two-layer CITRUS and use one-layer MLP for each with F 1 = F 2 = 4 to learn the representation in different SNRs to study the effect of each factor graph's noise on the stability performance of the network.Figure 2 shows the averaged prediction error in MSE between the true and predicted outputs over 10 random realizations. Firstly, we notice that increasing the SNR for each factor graph improves performance, illustrating the stability properties of CITRUS. Therefore, the effect of each factor graph's stability on the overall stability is well depicted in this figure, confirming the theoretical findings in Theorem 3.7. For instance, in SNR=10 for the first graph, the performance is still severely affected by the stability for the second graph, shown by color in Figure 2 . Finally, we observe that the standard deviation of the predictions is smaller for larger values of SNR, illustrating CITRUS's robustness. For experimentally validating Theorem 3.10, we first generate a initial graph signal X 0 ∈ R N ×F0 with F 0 = 12 on a 150-node Cartesian product graph with N 1 = 10, N 2 = 15. Therefore, we consider a 10-layer generation process based on ( 6) with only one-layer MLP for each layer and {F l = 12 -l} 10 l=1 , where the weight matrices are generated from normal distribution. We select ReLU as the activation function and {t l = 1} 10 l=1 . The over-smoothing phenomenon strongly depends on the smallest spectral gap of the factor graphs as stated in Theorem 3.10. We consider two scenarios of ln s - (2) ER = 0.1. So, we scale the MLP weight matrices by 100 and 2.5 to limit their maximum singular values. These settings result in λ = 0.06 and s = 0.004 ± 10 -5 for the first scenario; and λ = 0.07 and s = 7.691 ± 10 -6 for the other.We depict the outputs y(l) = log E(X l ) E(X0) and the theoretical upper bound y(l) = l(ln s l -2 P tl λ) across the number layers l in Figure 3 . We observe that in the left plot where {ln s l -2 P tl λ < 0} 10 l=1 , the theoretical upper bound approximates well the actual outputs during over-smoothing as stated in Theorem 3.10. On the opposite, on the right plot, in which {ln s l -2 P tl λ > 0} 10 l=1 , the theoretical upper bound is loose. We leave for future work finding a tighter upper bound.We evaluate and compare CITRUS on real-world traffic and weather spatiotemporal forecasting tasks, where we follow the settings in [10] for pre-processing and training-validation-testing data partitions. All hyperparameters were optimized on the validation set with the details in Section I. We use Adam as optimizer.Traffic forecasting. Here, we work on two well-known datasets for spatiotemporal forecasting: MetrLA [35] and PemsBay [36] . MetrLA consists of recorded traffic data for four months on 207 highways in Los Angeles County with 5-minute resolutions [35] . PemsBay contains 5-minute traffic load data for six months associated with 325 stations in the Bay Area by the same setting with [36] . The spatial and temporal graphs are created by applying Gaussian kernels on the node distances [37] and simple path graphs, respectively. The task for these datasets is, by having the last 30 minutes (T = 6 steps) of traffic recordings, we need to predict the following 15-30-60 minutes, i.e., H = 3, H = 6, and H = 12 steps ahead horizons. Our model first uses a linear encoder network, followed by three CITRUS blocks with 3-layer MLPs for each, where {F MLP l = 64} 3 l=1 . Therefore, the output node embeddings are concatenated with the initial spatiotemporal data. Finally, we use a linear decoding layer to transform the learned representations to the number of needed horizons. We exploit residual connections within each CITRUS block to stabilize the training process [17] . We use the mean absolute error (MAE) as the loss function and consider a maximum of 300 epochs. The best model on the validation set is applied to the unseen test set. We compare CITRUS against previous methods using MAE, mean absolute percentage error (MAPE), and root mean squared error (RMSE).Table 1 presents the forecasting results on the MetrLA and PemsBay datasets. We observe that, given the abundant training data in these datasets, the NN-based methods outperform the classic GSP-based methods like ARIMA [35] , G-VARMA [25] , and GP-VAR [25] . Similarly, the GNN-based models are superior compared to non-graph algorithms like FC-LSTM. More importantly, CITRUS outperforms state-of-the-art baselines in most metrics, especially in larger numbers of horizons (H > 3).Weather forecasting. We also test CITRUS for weather forecasting in the Molene [43] and NOAA [44] datasets. The Molene dataset provides recorded temperature measurements for 744 hourly intervals over 32 measurement stations in a region in France. For the NOAA, which is associated with regions in the U.S., the temperature measurements were recorded for 8, 579 hourly intervals over 109 stations. The pre-processing and data curation settings were similar to prior works on these datasets [10, 25] . The task here is by having the last 10 hours of measurements, we should predict the next {1 : 5} hours. For the Molene dataset, the embedding dimension and the initial linear layer of Table 2 illustrates the results of the weather forecasting task. We observe that CITRUS show superior forecasting performance compared to previous methods in all forecasting horizons. We note that the SGP model [42] , which is a scalable architecture, is the second-best performing method in Molene, possibly due to the small scale of this dataset. On the contrary, the NOAA dataset is larger, so models with higher parameter budgets perform better.Ablation study. The ablation study concerns the CITRUS modules responsible for jointly learning spatiotemporal couplings. To this end, we use two well-known architectures: TTS and STT. Therefore, we evaluate four possible configurations: i) TTS, where we first apply an RNN-based network (here, GRU), and then the outputs are fed into a regular GNN; ii) STT, which is exactly the opposite of the TTS; iii) Continuous TTS (CTTS), where we replace the GNN in TTS with a CGNN; and iv) Continous STT (CSTT), which is the opposite of CTTS. We report the results of this ablation study on the MetrLA dataset in Table 3 . We observe that CITRUS has superior results probably due to learning joint (and not sequential) spatiotemporal couplings by modeling these dependencies using product graphs with learnable receptive fields. We also observe that CTTS and CSTT outperform their discrete counterparts TTS and STT, probably due to the learning of adaptive graph neighborhoods instead of relying on 1-hop connections in regular GNNs.Hyperparameter sensitivity analysis. The sensitivity analysis is related to the number of selected eigenvector-eigenvalue (eig-eiv) pairs of the factor Laplacians. We select eigenvector-eigenvalue pairs in k ∈ {2, 10, 50, 100, 150, 200} out of 207 components in the MetrLA dataset. Table 4 presents the forecasting results in MAE along with the training time (per epoch) of this ablation study. We observe that selecting 50 eigenvector-eigenvalue pairs is enough to get good forecasting performance while having a fast training time. This experiment illustrates the advantages of CITRUS for accelerating training and keeping the performance as consistent as possible. 
conclusion and limitations: 
 In this paper, we proposed CITRUS, a novel model for jointly learning multidomain couplings from product graph signals based on tensorial PDEs on graphs (TPDEGs). We modeled these representations as separable continuous heat graph kernels as solutions to the TPDEG. Therefore, we showed that the underlying graph is actually the Cartesian product of the domain-specific factor graphs. We rigorously studied the stability and over-smoothing aspects of CITRUS theoretically and experimentally. Finally, as a proof of concept, we use CITRUS to tackle the traffic and weather spatiotemporal forecasting tasks on public datasets, illustrating the superior performance of our model compared to state-of-the-art methods.An interesting future direction for our framework could involve adapting it to handle other types of graph products, such as Kronecker and Strong graph products [23, 24] . Future efforts will also focus on finding tighter and more general upper bounds in the stability and over-smoothing analyses. For example, we will explore the possible relationships between the size of the factor graphs and these properties, either in a general form or for specific well-studied structures, such as ER graphs. Additionally, we plan to investigate more challenging real-world applications of the proposed framework, especially in scenarios involving more than two factor graphs.Next, by performing mode-(P + 1) unfolding on both sides of (5), one can write:EQUATIONA.3 Proof of Proposition 3.6Proof. The proof is presented by construction. Therefore, we first show it is true for the case of P = 2 as:EQUATIONThen, by assuming the theorem holds for the case of P = K and the definitions of Ã′ K := ⊕ K p=1 Ãp and Ã′ K := ⊕ K p=1 Ãp , we next show it also holds for P = K + 1 as follows:EQUATIONFor proving the upper bound of the norm of E, in a similar approach to the previous proof, we first prove for the case of P = 2 as:EQUATION)Using a similar approach to the proof in (28) , we assume the theorem holds for P = K as |||E||| ≤ K p=1 ε p , and, based on this information, then, we prove the theorem for the case of P = K + 1 as follows:EQUATIONwhich concludes the proof.A.4 Proof of Theorem 3.7Proof. The poof can be straightforwardly obtained by applying Proposition 1 in [30] based on the obtained results in Proposition 3.6.A.5 Proof of Lemma 3.9Proof. By construction, we first prove the theorem for the case of P = 2. The defined Laplacian L in (15) can be rewritten as:EQUATION)where we used the following property of Kronecker products: α(A ⊗ B) = (αA) ⊗ B = A ⊗ (αB) [45] . Next, by assuming that the theorem holds for P = K and the definition ofÂ′ K := 1 K and ⊕ K p=1 Âp := ⊕ K p=1Âp K , we prove it for P = K + 1:EQUATIONBesides, by the addition rule for the eigenvalues of a product graph adjacency or Laplacian (λ⋄ ij = λ (1) i + λ (2)j ), it can be easily seen that the eigenvalues of the normalized product adjacency ( Ŵ) and Laplacian ( L) in (31) are in the intervals of [-1, 1] and [0, 2], respectively, quite similar to the case for the normalized factor graphs, which concludes the proof. = 1 H = 2 H = 3 H = 4 H = 5 H = 1 H = 2 H = 3 H = 4 H = 5CITRUS (ours) 0.0174 0.0155 0.0132 0.0107 0.0086 0.0034 0.0030 0.0026 0.0024 0.0023  
experimental setting/details: 
 Question: Does the paper specify all the training and test details (e.g., data splits, hyperparameters, how they were chosen, type of optimizer, etc.) necessary to understand the results?Answer: [Yes] Justification: Section 4 and the Appendix mentioned all the needed details.Guidelines:• The answer NA means that the paper does not include experiments.• The experimental setting should be presented in the core of the paper to a level of detail that is necessary to appreciate the results and make sense of them. • The full details can be provided either with the code, in appendix, or as supplemental material. 
experiment statistical significance: 
 Question: Does the paper report error bars suitably and correctly defined or other appropriate information about the statistical significance of the experiments?Answer: [Yes] Justification: Our results are averaged over multiple random realizations and we have reported the standard deviations at the end of the Appendix (Section D).Guidelines:• The answer NA means that the paper does not include experiments.• The authors should answer "Yes" if the results are accompanied by error bars, confidence intervals, or statistical significance tests, at least for the experiments that support the main claims of the paper.• The factors of variability that the error bars are capturing should be clearly stated (for example, train/test split, initialization, random drawing of some parameter, or overall run with given experimental conditions). • The method for calculating the error bars should be explained (closed form formula, call to a library function, bootstrap, etc.) • The assumptions made should be given (e.g., Normally distributed errors). • It should be clear whether the error bar is the standard deviation or the standard error of the mean. Justification: We confirm that our paper, in every respect, meets the NeurIPS Code of Ethics.Guidelines:• The answer NA means that the authors have not reviewed the NeurIPS Code of Ethics.• If the authors answer No, they should explain the special circumstances that require a deviation from the Code of Ethics. • The authors should make sure to preserve anonymity (e.g., if there is a special consideration due to laws or regulations in their jurisdiction). 
broader impacts: 
 Question: Does the paper discuss both potential positive societal impacts and negative societal impacts of the work performed?Answer: [Yes]Justification: We have discussed about Broader Impacts of our work in the Appendix, i.e., Section F.Guidelines:• The answer NA means that there is no societal impact of the work performed.• If the authors answer NA or No, they should explain why their work has no societal impact or why the paper does not address societal impact. • Examples of negative societal impacts include potential malicious or unintended uses (e.g., disinformation, generating fake profiles, surveillance), fairness considerations (e.g., deployment of technologies that could make decisions that unfairly impact specific groups), privacy considerations, and security considerations. • The conference expects that many papers will be foundational research and not tied to particular applications, let alone deployments. However, if there is a direct path to any negative applications, the authors should point it out. For example, it is legitimate to point out that an improvement in the quality of generative models could be used to generate deepfakes for disinformation. On the other hand, it is not needed to point out that a generic algorithm for optimizing neural networks could enable people to train models that generate Deepfakes faster. • The authors should consider possible harms that could arise when the technology is being used as intended and functioning correctly, harms that could arise when the technology is being used as intended but gives incorrect results, and harms following from (intentional or unintentional) misuse of the technology. • If there are negative societal impacts, the authors could also discuss possible mitigation strategies (e.g., gated release of models, providing defenses in addition to attacks, mechanisms for monitoring misuse, mechanisms to monitor how a system learns from feedback over time, improving the efficiency and accessibility of ML). 
safeguards: 
 Question: Does the paper describe safeguards that have been put in place for responsible release of data or models that have a high risk for misuse (e.g., pretrained language models, image generators, or scraped datasets)?Answer: [NA] Justification: Our paper poses no such risks.Guidelines:• The answer NA means that the paper poses no such risks.• Released models that have a high risk for misuse or dual-use should be released with necessary safeguards to allow for controlled use of the model, for example by requiring that users adhere to usage guidelines or restrictions to access the model or implementing safety filters. • Datasets that have been scraped from the Internet could pose safety risks. The authors should describe how they avoided releasing unsafe images. • We recognize that providing effective safeguards is challenging, and many papers do not require this, but we encourage authors to take this into account and make a best faith effort.12. Licenses for existing assets Question: Are the creators or original owners of assets (e.g., code, data, models), used in the paper, properly credited and are the license and terms of use explicitly mentioned and properly respected?Answer: [Yes]Justification: All the needed references and credits have been explicitly mentioned in our paper. Since the used datasets are open-access, no permission was needed.Guidelines:• The answer NA means that the paper does not use existing assets.• The authors should cite the original paper that produced the code package or dataset.• The authors should state which version of the asset is used and, if possible, include a URL. • The name of the license (e.g., CC-BY 4.0) should be included for each asset.• For scraped data from a particular source (e.g., website), the copyright and terms of service of that source should be provided.• If assets are released, the license, copyright information, and terms of use in the package should be provided. For popular datasets, paperswithcode.com/datasets has curated licenses for some datasets. Their licensing guide can help determine the license of a dataset. • For existing datasets that are re-packaged, both the original license and the license of the derived asset (if it has changed) should be provided. • If this information is not available online, the authors are encouraged to reach out to the asset's creators. • We recognize that the procedures for this may vary significantly between institutions and locations, and we expect authors to adhere to the NeurIPS Code of Ethics and the guidelines for their institution. • For initial submissions, do not include any information that would break anonymity (if applicable), such as the institution conducting the review. 
 
 Appendix: Regarding the homophily-heterophily trade-off, first, we should mention that measuring the homophily-heterophily indices is way more challenging in spatiotemporal settings rather than in single-graph scenarios [46] , since we are facing two intra-graph (i.e., the spatial dimension) and intergraph (i.e., the temporal domain of changing time series characteristics) aspects. In this way, in the recent pioneer literature, they considered MetrLA, PemsBay, and temperature datasets (like Molene and NOAA used in our paper) as complicated in terms of homophily-heterophily behaviors. We provided the measured homophily-heterophily metrics (full details in [46] ) in Table 8 . As observed in this table, due to the wide range of variability across the used datasets, we have actually evaluated our performance on various datasets in terms of homophily-heterophily levels. For instance, Molene acts severely heterophilic in both intra and inter-graph domains, while NOAA acts more homophilic than Molene in the inter-graph domain.This scenario lies under the general umbrella of inaccurate adjacencies [30] , which has been rigorously studied in the stability analysis in Section 3.2. However, in a simplified and insightful case of facing graph products other than Cartesian, we consider two well-known graph products, i.e., Strong and Kronecker products [23, 24] . Consider the true graph is actually a Strong product graphA = A 1 ⊗ A 2 +A 1 ⊗I 2 +I 1 ⊗A 2 ,Ã = A 1 ⊗I 2 +I 1 ⊗A 2 . Then, E = A -Ã = A 1 ⊗ A 2 , where ∥E∥ = λ (1) max λ (2)max . So, the approximation error depends on the multiplication of the factor graph eigenspaces. Similarly, in the case of Kronecker product graphsÃ = A 1 ⊗ A 2 , it obtains ∥E∥ ≤ λ (1) max + λ (2) max + λ (1) max λ (2)max . So, in this case, the summation of the factor eigenspaces matters too. Briefly, the approximation error bound depends on the factor graph spectrums.We can choose an appropriate k using two methodologies: supervised and unsupervised. For the supervised option, we can use cross-validation, as we did in our results. For the unsupervised method, we can analyze the Laplacian matrices for a low-rank approximation task. For example, in Figure 4 , we plot the explained variances of the principal components of the spatial Laplacian (i.e., ∥vλv ⊤ ∥ 2 F for v and λ being the eigenvectors and eigenvalues) on the MetrLA dataset. We observe a strong concentration of variance in a few eigenvalues. We also observe that we can capture almost 80% of the explained variance by only relying on 50 components. Table 9 : Details of the training settings and hyperparameters, i.e., T (auto-regressive order), emb (dimension of embedding size in the spatiotemporal encoder), hid (dimension of linear mapping size in the spatiotemporal encoder), F MLP (dimension of linear mapping size in the MLP layers), n CITRUS (number of CITRUS blocks), F (second dimension of W l in ( 6)), k i (number of selected eigenvalue-eigenvector pairs in the i-th factor graph), lr (learning rate), n batch (batch size), and n epochs (number of epochs). Dataset T emb hid F MLP n CITRUS F k 1 k 2 lrFrom a theoretical point of view, based on eq. ( 18), one intuitive way to alleviate or at least slow down the over-smoothing phenomena is to keep ln s -2 P tλ positive, which tends to t < P 2 λ ln s. This clarifies why one should not increase the node-neighborhood especially in the case of facing a connected graph or increasing the number of layers. This finding is also compatible with the pioneering research about the over-smoothing concept, e.g., [14] . From an experimental point of view, and if we consider the graph receptive field t as a hyperparameter, we have designed an experiment to validate the claimed statement. Here, we vary t ∈ {0.1, 1, 5, 10, 20} and monitor the over-smoothing process in Figure 5 . In this figure, b = P 2 λ ln s, and as it is observed, for values of t higher than b the over-smoothing phenomena with increasing the number of layers is happening way faster than the other ones, which validates the theoretical discussion, too. Apart from this, the weight normalization technique [14, 15] has also been mentioned as a helpful action but it hinders the theoretical foundations of CITRUS. Note that, in practice, only escaping from the over-smoothing is not enough to have good performance. Therefore, for t < b, one might need to perform other kinds of analyses, e.g., over-squashing analysis.The detailed hyperparameters (optimized by cross-validation on the validation data) and/or training settings of CITRUS , i.e., T (auto-regressive order), emb (dimension of embedding size in the spatiotemporal encoder), hid (dimension of linear mapping size in the spatiotemporal encoder), F MLP (dimension of linear mapping size in the MLP layers), n CITRUS (number of CITRUS blocks), F (second dimension of W l in ( 6)), k i (number of selected eigenvalue-eigenvector pairs in the i-th factor graph), lr (learning rate), n batch (batch size), and n epochs (number of epochs), are provided in Table 9 . Full details can be found in the implementation codes in https://github.com/ArefEinizade2/CITRUS. All the experiments were run on one GTX A100 GPU device with 40 GB of RAM.The demonstrated superior performance of CITRUS in spatiotemporal forecasting tasks, such as traffic and weather prediction, can significantly benefit urban planning and public safety. Accurate traffic predictions can lead to better traffic management and reduced congestion, while improved weather forecasting can lead to serious urban inconvenience. This research can contribute to sustainable development goals by improving forecasting capabilities and decision-making processes, including sustainable cities and communities, climate action, and industry innovation and infrastructure. 


Decision: Overall, the reviewers agree that this is a good publication in terms of the proposed approach, clarity, theoretical analysis, and experimental results, and based on these merits, the paper should be accepted for publication. The feedback provided to the authors by the reviews as well as the additional experiments the authors conducted in the rebuttal should be included to improve the final version of their paper.

Reviews: review: paper_summary: paper_summary: The authors propose Tensorial Partial Differential Equations (TPDEG) to model multidomain data. Then they propose Continuous Product Graph Neural Networks (CITRUS) as a continuous solution. They provide theoretical and experimental analysis of their proposed approaches. They test their approach on spatiotemporal forecasting tasks and show SOTA performance.
strength_weakness: strength_weakness: See questions.
questions: questions: 1. Is it possible that real world graphs won't be representable as a cartesian product of factor graphs? How much approximation error is observed due to this assumption?
limitations: limitations: The authors address the limitations of their work as the last paragraph of their paper.
paper_summary: paper_summary: The paper proposes CITRUS, a novel model for jointly learning multidomain couplings from product graph signals based on tensorial PDEs on graphs (TPDEGs). By modelling these representations as separable continuous heat graph kernels as solutions to the TPDEG, the paper shows that the underlying graph is actually the Cartesian product of the domain-specific factor graphs. Then the paper studied the stability and over-smoothing aspects of CITRUS theoretically and experimentally. Finally, the paper applies CITRUS to tackle the traffic and weather spatiotemporal forecasting tasks on public datasets, illustrating its effectiveness compared to state-of-the-art methods.
strength_weakness: strength_weakness: 1. Bad notations: too many tilde, lower bar, subscript that seems unnecessary to me. They make it hard to read the math. Math is hard to communicate so good writing is essential.
- For example, what is purpose of having tilde in (3) and (4)? Then it disappeared in (5).
- Is lower bar for $U$ necessary? why not just go with $U$?
- is $\times_i$ mode-i tensorial multiplication? what is the precise definition?
- One minor thing --- $:=$ is used in line 120 but not line 123. Please check for consistency of other notations as well. 
2. The bound in Theorem 3.9 seems to be too loose to be meaningful in some cases, as shown in Figure 3, right plot that while the bound suggests divergence, the actual $E(X_l)$ converges to zero. Also, what is the implication of Theorem 3.9 on over-smoothing besides that we can focus on one factor graph for the phenomenon?
questions: questions: 1. Could the authors provide the standard deviation for the experimental results reported in table 1-4?
2. What are typical problems involving more than two factor graphs?
3. I see that Cartesian product arises in solutions to the TPDEG, that is why the CITRUS uses Cartesian product. I wonder if the authors have any intuition why this product is a good for modelling spatiotemporal couplings in practice.
4. It is also interesting that CITRUS seems to outperform especially in longer horizon predictions. Any investigation why it is the case?
5. Any guided way to choose k? Also any idea why a very small k seems to be a good enough approximation?
limitations: limitations: n/a
paper_summary: paper_summary: This paper proposes tensor PDE for graph neural networks temporal graph prediction. The construction is seems to be good and theoretical analysis for oversmoothing and stability are provided. Experiments show good performance with the proposed method compared to several baseline methods.
strength_weakness: strength_weakness: 1) There is no clear strategy to overcome oversmoothing.
2) Comprehensive details of hyperparameters is required. The lack of nature of hyperparameters and their ranges leads to difficulties in understanding the overall computational efficiency of the method.
questions: questions: 1) what are the training times for the CITRUS and its comparison to baseline methods (or the most relevant baseline)?

2) Over-smoothing effect is provided with a theorem and experimental validation, however, is there any suggestion or improvement in the proposed methods to overcome the over-smoothing? 

3) what are the hyperparameters of CITRUS?  how would you tune them? Do you hyparameter tune for eigenvector-eigenvalue pairs? Does some quantity such as rank is required for the tensor construction? 

4) Is there a specific structure of graph such as homophily and hitherophily in temporal graph that the proposed method is in more favor of?
limitations: limitations: yes
paper_summary: paper_summary: The authors propose a new spectral GNN for Cartesian product graphs. The graph filter is chosen as Laplacian heat diffusion, which is separable across factor graphs. The stability and over-smoothing of the proposed model are studied. The authors perform synthetic experiments to validate the analysis. Experiments on real datasets also show the advantages of the proposed methods.
strength_weakness: strength_weakness: My main concern is with the over-smoothing analysis. The normalized product Laplacian as defined in Lemma 3.8. is generally not the normalized Laplacian of the product graph. It is more of a manufacture, so that the smoothness of the product graph decomposes to factor-wise smoothness terms. This also means that Eq. 15 is not true, if the LHS is defined as the 'actual' smoothness on the normalized product graph as in Eq. 13. I think the authors need to revise the analysis here.

Another concern is that the advantages of a separable graph filter are not made explicit. This is somewhat discussed in Remark 3.5., where the authors argue that factor-wise EVD is more efficient than EVD of the product graph. However, note that this is merely a consequence of the KP eigenvector structure of the Cartesian graph product [1], not the separability of the heat diffusion filter. You don't need a full EVD of the product Laplacian even if the filter is not separable. How a separable filter further benefits the computation is unclear to me.

Minor points: figure 2 is not very intuitive; figure 3 is also not very convincing.

[1] Stanley, Jay S., Eric C. Chi, and Gal Mishne. "Multiway graph signal processing on tensors: Integrative analysis of irregular geometries." IEEE signal processing magazine 37.6 (2020): 160-173.
questions: questions: 1. Spectral GNNs that use polynomials of graph adjacency matrices only require local aggregation and do not need to compute global EVD. How does the complexity of the proposed method compare to that? Are there good reasons for choosing Laplacian-based filters over adjacency-based filters in spectral GNNs?

2. Can the authors provide more details on the competing methods for real graph data? What about the methods that also explicitly consider the product structure? The real datasets all have a temporal dimension. How do traditional SP methods such as [2] perform?

[2] Grassi, Francesco, et al. "A time-vertex signal processing framework: Scalable processing and meaningful representations for time-series on graphs." IEEE Transactions on Signal Processing 66.3 (2017): 817-829.
limitations: limitations: Future directions to improve the theoretical results are discussed.

Comments: title: Thank you
Response to authors
Further Responses
Response to author's rebuttal
More details about the responses
comment: We thank you again for the insightful discussion.

Regarding your mentioned point, actually, this is almost what we expected based on our previously provided detailed response, in which we showed both GTCNN and CITRUS naturally embed a graph filter smoothness term on the product graphs in their loss functions. Based on the deep study of the behavior of the smoothness regularization terms [1,2,3], they essentially do not allow the response to have high variations (around the mean). But, as we already outlined in our previous response, this term has a pure accumulating nature in GTCNN leading to (most probably) higher importance to gain a smooth response, while, in CITRUS, this term is getting adapted and modified with learnable importance (due to learnable graph receptive fields in our approach). In summary, we expected smooth results from both GTCNN and CITRUS, but more accurate and optimized for CITRUS. 

Thanks again.

References:
---
[1] "Graph regularized nonnegative matrix factorization for data representation." IEEE TPAMI, 2010

[2] "Fast robust PCA on graphs." IEEE JSTSP, 2016

[3] "Robust principal component analysis on graphs." ICCV, 2015
I thank the authors for their detailed answers, which help me understand the work better. For Q4, one thing I am not sure is that whether the smaller standard deviations in prediction performance shown in Tab. R1 are due to more accurate estimation of $t$ in longer horizons --- as the trend for standard deviations seems similar for GTCNN. 

I will keep my score and stay on the positive side.
Thank you for your insightful comments and for engaging in the discussion.
- As suggested, we will list the range and precise values of our hyperparameters in the Appendix.
- Regarding the heterophilic datasets, we should first mention that measuring the homophily indices in spatiotemporal settings is more challenging than in regular node classification scenarios. In our case, we are facing two intra-graph (the spatial dimension) and inter-graph (the temporal domain) aspects. The study in [1] considers Metr-LA, PemsBay, and temperature datasets (similar to Molene and NOAA used in our paper) as highly heterophilic.
To complement our rebuttal, we’ve studied the heterophilic measurements on the Molene and NOAA datasets, which we provide in the following table (${\pi}^{s}\_{v\_i}$ and ${\pi}^{T}\_{v\_i}$ corresponds to intra and inter-graph homophily, respectively). Due to the wide range of variability across these datasets, we have already evaluated our performance on various datasets in terms of homophily-heterophily levels. For instance, Molene is highly heterophilic in both the intra and inter-graph domains, while NOAA acts more homophilic than Molene in the inter-graph domain. We can also compare the provided metrics for NOAA with the Temperature or KnowAir datasets used in [1], where we observe similar dynamics. We will include these discussions related to homophily-heterophily metrics in the Appendix of the revised manuscript.
| Datasest |  ${\pi}^{s}\_{v\_i}:p\_s$ | ${\pi}^{s}\_{v\_i}:q\_p$ | ${\pi}^{s}\_{v\_i}:q\_n$ | ${\pi}^{T}\_{v\_i}:p\_s$ |  ${\pi}^{T}\_{v\_i}:q\_p$ | ${\pi}^{T}\_{v\_i}:q\_n$ |
|----------|----------|----------|----------|----------|----------|----------|
| MetrLA   | 0.2273 | 0.4732 | 0.2995 |  0.3325 | 0.4920 | 0.1755 |
| PemsBay    | 0.1073 | 0.5912 | 0.3015 |  0.2399 | 0.6863 | 0.0738 |
| Molene    | 0.0148 | 0.6248 | 0.3602 | 0.0152 | 0.5545 | 0.4301 |
| NOAA    | 0.1249 | 0.5345 | 0.3405 | 0.2862 | 0.4351 | 0.2786 |

References:
---
[1] "Greto: Remedying dynamic graph topology-task discordance via target homophily”, ICLR, 2023
Thank you for the response.

I think the listing of the hyper-parameters and their ranges are important to have in the revised version of the paper. Furthermore, I think more studies should be helpful on commonly used hetherophilic graph datasets in the revised paper.
Thank you for engaging in the discussion and for your helpful and insightful comments.


**W1**. We will revise all that's required to align with the correct conclusions from our discussion.


We agree that for the *general* normalized Laplacian, the decomposability might not be true. However, please notice that for our definition in Eqn. (14) in our paper, the decomposability holds. For example, if we consider two (general, and not even graph-related) factor matrices, the following relationship is valid for any general matrices $A,B,C$ with consistent dimensions [1]:
- $AX+XB=C \leftrightarrow (A\oplus B^\top)vec(X)=vec(C).$


Now, by considering $x=vec(X)$ and using this relationship, one can write:
- $vec(X)^\top\left[(A\oplus B^\top)vec(X)\right]=vec(X)^\top\left[vec(AX+XB)\right]=vec(X)^\top vec(AX)+vec(X)^\top vec(XB)=tr(X^\top AX) + tr(X^\top XB)=tr(X^\top AX) + tr(XBX^\top).$


If $A=\hat{L}_1$ and $B=\hat{L}_2$, the proof is completed. Please note that our over-smoothing analysis is easily adaptable for general combinatorial Laplacian matrices as well.


**W2**. We agree that even in the non-separable filters, one can select a graph frequency from the factor graphs. However, we don't find a strong theoretical justification for this. There might exist a logic in specific cases, like in [2], but this requires the general study of the product graph's spectrum. For example, consider the non-separable filters are $h(\lambda_1,\lambda_2)=e^{\lambda_1\lambda^2_2-\lambda^2_1\lambda_2}$, or even $h(\lambda_1,\lambda_2)=\cos{(\lambda_1\lambda_2)}$. In these cases, the importance of $\lambda_1$ on the spectrum of the product graphs is also tied with the importance of $\lambda_2$. Therefore, the most important eigenvalue in the first-factor graph for performing low-pass filtering is not necessarily the one with the highest importance in the product graph spectrum. This is because the importance of the factor graphs' spectra is hindered by the non-separable function. Please notice that the theoretical justification is clear in our formulation since, based on our response in the previous rebuttal stage, the selection procedure in our paper comes from factor-wise low-rank approximation subproblems.


As a final remark, please note that this whole discussion about the separability boils down to just one remark in our paper, and not to one of our main contributions. Therefore, we can easily make the relevant modifications without an important repercussion in our main takeaways.


References
----
[1] "The matrix cookbook", 2008

[2] "Learning Cartesian Product Graphs with Laplacian Constraints", AISTATS, 2024.
Thank you for engaging in the discussion.


**W2**. While it's true the right plot of Fig. 3 is not tight enough, we should emphasize this is an open research question. Our findings align with previous work in over-smoothing, please refer to Fig. 2 in [1] for example, where a similar loose bound was found.


We now provide insights for possible future directions to find a tighter bound:
- The RHS in Theorem 3.9 also depends on the maximum singular value of the learnable weights matrices $W\_{l1},...,W\_{ll}$, *i.e.*, $s$, which solely depends on the training process. We can use techniques like weight normalization [1] to bound $s$.
- When there's a big gap between the arguments of the Lipschitz functions, we might use integral Lipschitz functions [2] since in that case $|f(x\_2)-f(x\_1)|\le C\frac{|x\_2-x\_1|}{|x\_2+x\_1|/2}$. This is tighter than an absolute Lipschitz function, but its effect on the factor graph spectrums needs further theoretical exploration.
- Previous works have found that the Frobenius norm might not be very tight as an upper bound for the spectral norm [3]. Precisely, $\\|T\\|\le\\|T^K\\|^{\frac{1}{K}}=\sup\_{x\ne 0}{\left(\frac{\\|T^Kx\\|}{\\|x\\|}\right)^{\frac{1}{K}}}=\sup\_{x:\\|x\\|=1}{(\\|T^Kx\\|)^{\frac{1}{K}}}$, for some finite $K$ [3]. Since $E(x)=x^\top L x=\\|L^{1/2}x\\|\_F^2=\sum_{i=1}^{N}{\lambda\_i\tilde{x}^2\_i}$ with $\tilde{x}$ being the GFT of $x$ (which is a weighted nuclear norm with direct relationship with the spectral norm), we might use other norms like Sobolev [4], gradient-adapted [5], or other $p-$norms [4,6].
- Recently, some works have explored alternatives to the regular triangular inequalities [7], which might not be tight for the over-smoothing analysis.


In summary, there are some future directions to make our bounds tighter, but they need to be carefully studied in the case of continuous filters in product graphs.


**Q4**. By robustness, we mean smaller standard deviations in prediction performance as shown in Tab. R1 in the **uploaded PDF file of the main rebuttal**.


Regarding the longer horizons, the authors from the GTCNN model [8] (which can be considered as a particular case of CITRUS for discrete, two-factor graphs) said: "The GTCNN outperforms the other models in a short horizon while Graph WaveNet and GMAN work better for longer horizons. The benefits in the short term are due to high-order spatiotemporal aggregation in the GTCNN which allows capturing efficiently the spatiotemporal patterns in the data". Therefore, we state the GTCNN model doesn't perform well in longer horizons due to the non-optimal receptive field of the discrete case. More formally, consider a one-layer CITRUS network without non-linearity for predicting one horizon as $\tilde{y}=e^{-tL\_{\diamond}}Xw$. The loss function is therefore given by $\min\_{w}{\\|y-e^{-tL\_{\diamond}}Xw\\|\_2^2}=\min\_{w}{w^\top X^\top e^{-2tL\_{\diamond}}Xw-2y^\top e^{-tL\_{\diamond}}Xw}$. Here, the second term tries to maximize the inner product (the similarity) between the target $y$ and predicted output $\tilde{y}=e^{-tL\_{\diamond}}Xw$. The first term $f_{CITRUS}=w^\top X^\top e^{-2tL\_{\diamond}}Xw$ enforces the output to be as smooth as possible on the heat diffusion filters $e^{-tL\_{\diamond}}$. For a simple GTCNN [8], this smoothness term takes the form of $f\_{GTCNN}=w^\top X^\top L\_{\diamond}^2Xw$. If we consider longer horizons, these terms will accumulate leading to bigger smoothness terms in the loss function, likely leading to over-smoothing. The key difference between $f\_{CITRUS}$ and $f\_{GTCNN}$ is that our model learns $t$ to alleviate this accumulation issue as observed in Fig. R3 in the **uploaded PDF file of the main rebuttal**. Please also note that we can also have better control of the spectrum of $f\_{CITRUS}$ because of the learnable heat kernel.


Apart from the role of $t$ in $f_{CITRUS}$, other aspects besides over-smoothing could also play an important role, like over-squashing and over-fitting, which we leave for future work.


References
----
[1] “Graph neural networks exponentially lose expressive power for node classification”, ICLR, 2020

[2] “Stability properties of graph neural networks”, IEEE TSP, 2020

[3] “Learning interface conditions in domain decomposition solvers”, NeurIPS, 2022

[4] “Reconstruction of time-varying graph signals via sobolev smoothness”, IEEE TSIPN, 2022

[5] “Time-varying graph signal reconstruction”, IEEE JSTSP, 2017

[6] “Flow smoothing and denoising: Graph signal processing in the edge-space”, IEEE GlobalSIP, 2018

[7] “Studying the effect of gnn spatial convolutions on the embedding space’s geometry”, UAI, 2023

[8] "Graph-time convolutional neural networks: Architecture and theoretical analysis", IEEE TPAMI, 2023.
Thank you for reply. I acknowledge reading the rebuttal.
I thank the authors for their detailed reply.

W1. Please make sure to revise the paragraph at line 173. Also the author should realize that [9] used combinatorial Laplacian and that's why the smoothness is naturally decomposable. Normalized Laplacian is a whole different story. The authors should also pay attention to Eq.13 and Eq.15, and make clear distinction between them. Also make sure you don't accidentally use $E(U_t)$ as the form in Eq.13 in the proof.

W2. I think one certainly can make the filter 'selective' even without separability. For example, consider a low-pass filter $f_i$ on some factor $i$ and a non-separable filter $g$ on the product graph. Their composition $f_i \circ g$ is not separable either, but selective on factor $i$. To be more specific, you can compute the EVD of factor graphs, transform the data to spectral domain, apply a non-separable filter, then select a subset of frequencies for each factor. I don't see how this is not as efficient as using a separable filter like the heat diffusion.
Thank the authors for their rebuttal. There are a few further comments and questions I want to discuss:

W2: I guess my concern was not that the Theorem is wrong but not tight enough. It seems that the authors admit that the current bound might be vacuous in certain cases. Is there a way to make the results tighter? 

Q4: What do you mean when saying "the estimation of  is more robust for longer horizons" in Fig R3 --- do you mean the variance is smaller? I can see that this learned t is not under the scenario of oversmoothing established in the theoretical results and hence it helps alleviate oversmoothing. But I wonder if alleviating oversmoothing itself is the only reason behind the good long horizon performance (so other methods are bad because they oversmooth?)
**W1**. 
- As a simple toy example, consider a $4\times3$ matrix $U=[1,2,3;1,2,3;4,5,6;4,5,6]$, and also assume the factor graph as 4 and 3-node simple path graphs. So, $U$ is smooth on the 4-node graph but it is not the case with the other 3-node one, and, so, can not be considered a pure smooth signal on the product graph.

**W2**.
- For the sake of completing the discussion, consider we could write the resulting multi-way graph filter as $h(L)=h\_1(L_1)\otimes h\_2(L\_2)$, where $h\_{1(2)}()$ is a polynomial graph filter (with possible infinite order like heat graph filters). Then, one can use low-rank expansion $h\_1(L\_1)=\sum_{i=1}^{K\_1}{h\_1(\lambda\_i) v\_i v^T\_i}$. Therefore, the computation can be improved by choosing an appropriate small enough $K_1$.

**W3**.
- Figure 2 illustrates the effect of factor stabilities on the overall stability, as Theorem 3.7 states. By varying adjacency SNR for each factor (related to $\mathcal{O}(\epsilon_p)$), this figure validates the presence of factor stability on the overall performance. We've included a 3D plot in the **uploaded PDF** on the main rebuttal (Fig. R5) with a more intuitive view of the same experiment for better illustrating the interpretation of Theorem 3.7.

- Regarding Figure 3, please first note that it is on a logarithmic scale, which we forgot to mention and has been corrected in the revised version of the manuscript. The left plot is associated with the case of $\ln{s}-\frac{2\tilde{t}\tilde{\lambda}}{P}<0$, which is prone to over-smoothing and converges to zero very fast (exponentially). That's why the difference between $\ln{(\frac{E(X_l)}{E(X_0)})}$ and $l(\ln{s}-\frac{2\tilde{t}\tilde{\lambda}}{P})$ is extremely small. For the right plot in Figure 3, we have that  $\ln{s}-\frac{2\tilde{t}\tilde{\lambda}}{P}>0$. Here, again the theorem is validated by the RHS bound in Eqn. (17). However, this bound is not tight and we leave further analysis as future work. We kindly refer the reviewer to [2] to get more familiarized with this kind of analysis and plots regarding over-smoothing in discrete GNNs.