Continuous Product Graph Neural Networks
Authors: Aref Einizade
Fragkiskos D. Malliaros
Jhony H. Giraldo
Venue: NeurIPS 2024 poster (2024)
Keywords: Continuous graph neural networks
product graphs
Cartesian product
Primary area: graph_neural_networks

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.

Continuous Product Graph Neural Networks
Aref Einizade
Fragkiskos D. Malliaros
LTCI, Télécom Paris
CentraleSupélec, Inria
Institut Polytechnique de Paris
Université Paris-Saclay
aref, ein.iiaddeeeeeeeeomppriis....
fragkiskos mal.llarresseenttreeepeee...
Jhony H. Giraldo
LTCI, Télécom Paris
Institut Polytechnique de Paris
jhony. gww.lllddteleeecmmprrr.aa......e
Abstract
Processing multidomain data defined on multiple graphs holds significant potential in various practical applications in computer science. However, current meth- ods are mostly limited to discrete graph filtering operations. Tensorial partial differential equations graphs (TPDEGs) provide principled framework for modeling structured data across multiple interacting graphs, addressing the limita- tions 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 the performance. We evaluate CITRUS well-known traffic and weather spatiotemporal forecasting datasets, demonstrating superior performance over existing approaches. The implementation codes available http:///ahhbb.oomm/Aeeenmddddddt
Introduction
Multidomain (tensorial) data defined on multiple interacting graphs (1-3). referred as multidomain graph data in this paper, extend the traditional graph machine learning paradigm, which typically deals with single graphs 2 Tensors, which 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 upplications. 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 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) [33. However, these methods inherit the well-known issues of over-
smoothing and wer-squashing from regular graph neural networks (GNNs) [14 16). which restricts the graph's 's receptive field and hinders long-range interactions [17]. Additionally, these methods often
require computationally intensive grid searches tune hyperparameters and are typically limited to
two-domain graph data, such as spatial and temporal dimensions 10
12
18
19
38th Conference on Neural Information Processing Systems (NcurIPS 2024).

---

In regular non-tensonee non-tensorial GNNs, continuous GNNs (CGNNs) emerge as solution to over-smoothing and over-squashing [20]. CGNNs are neural network solutions to partial differential equations (PDEs) on graphs [ 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 mul- tidomain data residing on multiple interacting graphs that form Cartesian product graph heat kernel. We then propose Continuous Product Graph Neural Networks (CITRUS) as 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 spatiotem- poral forecasting tasks, demonstrating that CITRUS achieves state-of-the-art performance compared to previous Temporal-Thee Seepoorl-Thee-Sppee Then-Spatial (TTS), Spatial- Spatial-Then- Temporal (STT), and Temporal-ann and-Spatial (T&S) frameworks. Our main contributions 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 over- smoothing properties of CITRUS.
We test CITRUS on spatiotemporal forecasting tasks, demonstrating that our model achieves state-of-the-art performance.
2
Related Work
We divide the related work into two parts covering the fundamental and applicability aspects.
Fundamentals. The study of PGs is well-established field in mathematics and signal processing 24.. However, to the best of our knowledge, the graph-time convolutional neural network GICCNN [8) model is the only GNN-based framework addressing neural networks in PGs. GTCNN defines a discrete filtering operation 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 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 [7 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 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 GAR models (25). do not utilize PGs. Modern GNN- based 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], 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 bastdr-based functions defined on Cartesian product graph. The continuity and differentiability

---

Factor Graphs
Product Graph
Continuous Product Graph Function (CITRUS)
L + L2 Lg)
U
U
L 2
L
tL,
L2
-tg
x X1e
x2
x3
(a) Cartesian product graph.
(b) Continous product graph function in
Figure I: Illustration of key concepts of CITRUS. a) Cartesian product between three-factor graphs. b) Continous product graph function (CITRUS) operating on the multidomain graph data U.
of these filters allow the graph receptive fields to be learned adaptively during the training process, climinating the need for grid search. Additionally, CITRUS maintains low computational complexity by relying on the spectral decomposition of the factor graphs [2] Furthermore, unlike non- non-graph approaches, the number of learnable parameters in CITRUS independent of the factor graphs, ensuring scalability.
Methodology and Theoretical Analysis
Preliminaries and notation. An undirected and weighted graph 9 with N vertices can be stated as g E Aj. J. V and denote the sets of the vertices and edges, respectively. A E NxN is the adjacency matrix with {A. Ajt 0] representing the connection between vertices, and {A, = 0} states that there are no self-loops in G. The graph Laplacian E NxN is defined as L graph signal is function that maps a set of nodes to the real values x V - R. and therefore we can represent it The vectorization operator is shown as vec(.). Using the Kronecker product @, Kronecker sum (aka Cartesian product) @ and the Laplacian factors (L, E RRN, )=== we have the following definitions:
A, where D diagí A1) 11 is the diagonal degree matrix and 1 is the all-one vector. A
[æ1...., .N For vector
[e.......N]]]
Nx1, ea == [
4 p L, == Lp .... ① L1, Lpp .... Lp, 4 & Lp = Lp & ... L1, (1)
where the Laplacian matrix of the Cartesian product between two factor graphs can be stated as:
LL = L + L2 = L 6 IN,
INj X L2,
with I, the identity matrix of size n. See Figure la 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 D-dimensional tensor as U R xp.Nxx..xNN. The matrix U E R.N. [
the i-th mode matricization of U 28 The mode-i tensorial multiplication of tensor
matrix X E mm*XN N; is denoted by G = U x X, where G E R.Nxx...x x Noj x No [28
For further information about the tensor matricization and mode-n product, please refer (28 29
(especially Sections 2.4 and 2.5 in (28). All the proofs of theorems, propositions, and lemmas of this
paper are provided the Appendix A
3.1 Continuous Product Graph Neural Networks (CITRUS)
We state the generative PDE for CITRUS as follows:
Definition 3.1 (Tensorial PDE graphs (TPDEG)), Let U, N.N............. be multidomain tensor whose elements are functions of time t, and let (Lp E R. xN,,,,,, xN 1p=1 be the domain specific factor Laplacians. We can define the TPDEG follows:
aú, 01
Ú Xp Lp Σ p=1
(3)
3

---

Theorem 3.2. Let Ug be the initial value of U. The solution to the TPDEG in 3 is given by:
U. To X1 the X2 the X3 xp
Using Theorem 3.2 we define the core function of CITRUS as follows:
(W) U X1 444 X2... Xpe t
&.p XP+1 WT,
(5)
where U, E the learnable graph receptive field, and f(U,) I x..xxxx..p.... x P the oulput tensor. Figure
Npx Fi is the input tensor, W, E
E R is matrix of learnable parameters,
[b illustrate a high-level description of f(U,). The next proposition formulates D as a graph convolution defined on product graph.
Proposition 3.3. The core function of CITRUS in (5) can be rewritten as:
[((L))(++)] =e e e ((((4)) Wt,
(6)
where L. : Op= 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 [Lp Ap V
V p- and product graph {L.
V.A. V ) where V. =4 o sp= Vp and A, =4 A 23 Let K2 S Np be the number of selected aggenvalue-eigenved that we can rewrite (6) as follows:
pairs of the p-th factor Laplacian. When Kp Np. it can be
shown 17 22
Fi times
[)(U))(PD)]] v(Ka)
X A
(K())
U ((PH)))) Wi,
A
with X =1 10
A1 (Kp) V((K)) =1 P v(Ks),
(8)
where 2 (K)p) E K and (((k) N, are the first Kp I Np selected eigenvalues and eigenvectors of L, based on largest eigenvalue magnitudes, respectively, and @ is the point- wise multiplication operation. Finally, we can define the output of the 1-th layer of CITRUS as U1+1(P41) (f(U))(++1)) where (-)) is some proper activation function.
Remark 3.4. We can consider factor-specife learnable graph receptive fields (t/²) for formulating
d 8 X =1 p ((p) make CITRUS more flexible for each factor graph's 's receptive
field. This technique can be expanded to the channel-wist graph receptive fields (t) (p.c) P.F.
(
1p=1,011 such
that the e-th column of Å₁ in 7 is obtained by År. =1 co -
((Pi6)) x (Ky)
Remark 3.5. Computing the EVD in the product graph L, in (6) has a general complexity of
P 77 by relying upon the properties of product graphs since we perform EVD on each factor graph independently. Besides, if we rely only on the Kp most important eigenvector-eigenvalue pairs [17 reduce the complexity of the spectral decomposition to O(N? )
O p= N21²). However, we obtain a significant reduction in complexity of O(E NN,1) in
of each factor graph,
obtaining a general complexity of O(EP)) N3 Kpl).
3.2 Stability Analysis
We study the stability of CITRUS against possible perturbations in the factor graphs. First, as defined in the literature [30]. we model the perturbation the p-th graph an addition of an error matrix E, (with upper bound Ep for any matrix norm III (I)) to the adjacency matrix Ap as:
Å A, + Ep; || Ep.
(9)
Proposition 3.6. Let and À be the true and perturbed udjacency matrix of Cartesian product graphs with the perturbation model for each factor graph described as in (9). Then, it holds that:
- A
EEI
Eps
(10)
where the perturbation matrix E also follows the Cartesian structure E p Ep.

---

Finally, let p(u, t) and (fu, t) be the true and perturbed output of the CITRUS model with normalized Laplacian Lo being responsible for generating (P+1) in heat flow PDE with normalized Laplacian Lg @p(a.t) ((, in (6) (30), respectively. Then, the following theorem states that the 26 integrated error bound on the stability properties of CITRUS also separable when dealing with Cartesian product graphs:
Theorem 3.7. Consider Cartesian product graph (without isolated nodes) with the normalized
true and perturbed Laplacians Lo and Lo in (6). Also, assume we have the perturbation model Ap - Ap + Ep with factor error bounds Ell S En in Proposition 3.6 Then, the stability bound on the true and perturbed outputs of CITRUS, i.e., p(u, 1) and (u,, 1}, respectively, can be described by the summation of the factor stability bounds as:
l,y( ((u, t)ll = O((p).
(11)
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 erroncous factor graphs, and will be numerically validated in Section 4.1
3.3
Over-smoothing Analysis
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 (1 and can be formulated (20) for the output of continuous GNN. Precisely, [(())]] with u(!), being the feature vector
can be defined on the GNN's output U, of the i-th node with the degree deg,, as:
lim E(U,) - 0,
(12)
where E(U) : =
Σ
u
tr (U(LU
(13)
deg,
deg,
(KJ)EE
Inspired by (3), 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
{[p) we define the Tensorial Dirichlet energy for tensor U E R XN.xx..XNxxx as:
1 F
EU
tr(U fm L, (,
(14)
f=1
where Ui. Ul....., in fle N, xII iee, the p-th mode matricization of the f-th slice of in its (P 1)-th dimension.
Based on the previous research in the GSP literature 27
33 34
it follows that E(U) in can Ül: f] = Ul. fl), and
be rewritten as E(U) = tr( Ü LŪ), where Ü R!T
:: # fp i, 3/- (4z) The next lemma shows interesting and applicable properties L. Lemma 3.9. Consider P factor graphs with normalized adjacencies { }pm and Laplacians (Lp)yer By constructing the product adjacency as À = + S Å = 0f P
a similar Cartesian form for the product Laplacian L == I À (with I being the identity matrix) also holds as the following:
i
p=
L.
p=
P
(15)
Besides, the spectrum of L is spanned across the interval of [0, 2].
5

---

We observe in (14) that the kernel separability of Cartesian heat kernels also leads to the separabilityy of analyzing wer-smotthin different modes of the data at hand. This useful in cases of facing factor graphs with different characteristics. Next, we formally analyze over-smoothing in CITRUS. In our case, for the 1-th layer with H, hidden MIP layers, non-linear activations a()) and learnable weight matrices {WiJee WiJker we define:
X = Ji(X)), Ji(X == (L LG X), MLP == ((((..
Wir) )
(16)
with Xo being the initial node feature matrix, fi(-) generalized form of ( and t. :=4
8 e 4, where {(() 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 ::= to of=1 Lp with domain sppaif- receptive field ( for the i-th factor graph and the activations o(.) in (16) being ReLU or Leaky ReLU, we have:
E(XI) S In E(Xo)
(17)
where s == -.UIIEN, 81 and 81 == I 8th with Sth being the square of maximum singular value of WT. Besides, (mm) and = {(() with m = arg min, { (6), where X is the smallest non-zero eigenvalue L,,
Corollary 3.11. When + co, E(xi) exponentially converges to 0, when:
2 Ins P i
(18)
Theorem 3. 3.
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
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, conduct ablation studies regarding CITRUS in spatiotemporal forecasting. We compare CIT- RUS against several state-of-the--e methods in forecasting. The results on the case of more than two factor graphs and descriptions of the state-of-the-art provided in Appendix and C respectively, alongside the additional theoretical and experimental discussions in Sections D The implementation codes available at https///gihubb.oom//effffiiaaaat
4.1 Experimental Stability Analysis
We generate a 600-node product graph consisting of two connected Erdös-Rényi (ER) factor graphs with N₁ 20, N₂ 30, and edge probabilities (1) (2) 0.1 for node regression. We utilize 0.15
Synthetic Data
0.20
PER PFR
15% of the nodes in the product graph as the test set, 15% of the remaining nodes validation, and the rest of the nodes for training. We use a three-layer generative model to generate the outputs based (6). The continuous parameters for the factor graphs are set as (1) = and (2) = 3. The initial product graph signal (0) E N is generated from the normal distribution with Fo = 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 -o-ise ratios
MSE
SNRgminf
SNRg=20
0.10-
SNR==10
SNR,-0
SNR
0.08
10
10
SNR:
Figure 2: Stability analysis vs. different SNR scenarios, related the results in Theorem 3.7

---

(SNRs) (in terms of db) of foo, 20, 10, -10). Therefore, we construct two-layer CITRUS and use one-layer MLP for each with Fi = F2 = learn the representation in different SNRs 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 is 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.
4.2 Experimental Over-smoothing Analysis
For experimentally validating Theo- rem 3.10 we first generate initial graph signal Xy E R Nx F with Fu = 12 on 150-node Cartesian product graph with N, = 10, N₂ = 15. Therefore, we consider 10-layer gen- eration process based on (6) with only one-layer MI.P for each layer and {Fi 12 - 1 where the weight matrices are generated from normal distribution. We select ReLU as the activation function and (t, 1)
Pos: 10.05, 0.95), 0.06, s: 0.004
Pri: [0.1, 0.1 A: 0.07, 7.691
20.
Log Rebaivv Distance
Lag DistaiveeDstttnn
-15-
10-
-50-
-10
-75
Actual
Actual
Theorem
-00-
Theorem
-100
0
10
10
Namber of layers
Number of layers
Figure 3: Over-smoothing analysis using Theorem 3.10 left:
In 3id < 0. right: In = > 0.
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 In 3 < 0 and In - _ >> 0. The first scenario relates to a well-connected product graph obtained from two factor ER graphs with
(2) 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
(1) = 0.05 and PER
(2) = 0.95 and, while in the other, PER PeR
and = 0.004 = 10 5 for the first scenario; and A = 0.07 and = 7.691 + 10 6 for the other.
E(X) and the theoretical upper bound y() = 1(n 81 ( We depict the outputs yl) = log El X
across the number layers I in Figure We observe that in the left plot where {In st 3 < 0)
the theoretical upper bound approximates well the actual outputs during over-smoothing stated in Theorem 3.10 On the opposite, the right plot, in which (In 8 ξ 0) the theoretical
upper bound is loose. We leave for future work finding a tighter upper hound.
4.3 Experiments on Real-world Data
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 vaaiinn-valddatio--testtt data partitions. All hyperparameters were optimized on the validation set with the details in Section We use Adam as optimizer.
Traffic forecasting. Here, we work on two well-known datasets for spatiotemporal forecasting: MetrLA (35) and PemsBay (36). MetrI.A consists of recorded traffic data for four months on 207
highways 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 = 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 {FMMPP 64) 64) Therefore, the output node embeddings are concatenated with the initial spatiotemporal data. Finally, we use línear 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 [7]. We use the
7

---

Table 1: Traffic forecasting comparison between CITRUS and previous methods.
MenrL
H
H
H
Method
MAE
MAPE
RMSE
MAE
MAPE
RMSE
MAE
MAPE
RMSE
ARIMA
35
3,99
9.60%
8.21
5.15
12.70%
10.45
6.90
17.40%
13.23
G-RARMA
25
3.60
6.89
11.22%
7.84
5.12
14.00%
GP-VAR
3.56
9.55%
6.54
11.02%
7.56
4.87
13.34%
FC-LSTM
3.44
9.60%
6.30
3.77
10.90%
7.23
4.37
13.20%
8.69
Graph Wavenet
38
2.69
6.90%
5.15
3.07
8.37%
6.22
3.53
10.01%
7.37
GMAN
39
2.94
7.51%
5.89
3.22
8.92%
6.61
3.68
10.25%
7.49
STGCN
40
2.88
7.62%
5.74
3.47
9.57%
7.24
4.59
12.70%
9.40
GGRNN
TT
2.73
7.12%
5.44
3.31
8.97%
6.63
3.88
10.59%
8.14
GRUGCN
2.69
6.61%
5.15
7.96%
6.04
SGP
3,06
7.31%
5.49
3.43
8.54%
6.47
10.53%
7,81
GTCNN
10 18
2,68
6.85%
5.17
3.02
8.30%
6.20
3.55
10.21%
7.35
CITRUS (Ours)
2.70
6.74%
5.14
2.98
7.78%
5.90
3.44
9.28%
6.85
PemsBay
H
6
Method
MAE
MAPE
RMSE
MAPE
RMSE
MAE
MAPE
RMSE
ARIMA
35
1.62
3.50%
3.30
2.33
5.40%
4.76
3.38
8.30%
G-VARMA
25
1.88
4.28%
3.96
2.45
5.42%
4.70
7.10%
GP-VAR
1.74
3.45%
3.22
2.16
5.15%
4.41
2.48
6.18%
5.04
FC-LSTM
35
2.05
4.80%
4.19
2.20
5.20%
4.55
2.37
5.70%
4.74
Graph Wavenet
38
1.30
2.73%
2.74
1.63
3.67%
3,70
1.95
4.63%
4.52
GMAN
39
1.34
2.81%
2.82
1.62
3.63%
3.72
1.86
4.31%
4.32
STGCN
40
1.36
2.90%
2.96
1.81
4.17%
4.27
2.49
5,79%
5.69
GGRNN
4T
1.33
2.83%
2.81
3.79%
3.94
5.21%
GRUGCN
1,21
2.49%
2.52
3,32%
3.46
4.72%
SGP
1.33
2.71%
2.84
1.70
3.61%
3.83
2.24
5.08%
5.19
GTCNN 10
18
1,25
2.61%
2.66
1.65
3.82%
3.68
2.27
5.11%
4.99
CITRUS (Our))
1,21
2.51%
2.61
1.48
3.23%
3.28
1.78
4,08%
3.99
Table 2: Weather forecasting comparison (by rNMSH) between CITRUS and previous methods.
Molene
NOAA
Method
H
H
3
H
H
H
H =
3
H
H
GRUGCN
26
0.49
0.56
0.63
0.70
0.75
0.15
0.18
0.24
0.30
0.37
GGRNN NT
0.29
0.42
0.54
0.65
0.75
0.19
0.20
0.27
0.38
0.48
SGP
0.24
0.37
0.49
0.59
0.69
0.39
0.41
0.43
0.46
0.49
GTCNN 100 18
0.39
0..45
0.52
0.60
0.68
0.17
0.19
0.25
0.31
0.37
CITRUS (Ours)
0.23
0.35
0.47
0.58
0.67
0.04
0.12
0.20
0.29
0,36
mean absolute error (MAE) as the loss function and consider 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).
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- 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-an art baselines in most metrics, especially in larger numbers of horizons (H > 3).
Table
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 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 these datasets [1025] The task here by having the 10 hours measurements, we should predict the next {1 5) hours. For the Molene dataset, the embedding dimension and the initial linear layer of

---

Table 3: Ablation study comparison between the proposed CITRUS and typically ST pipelines.
MenLL
3
H =
12
Method
MAE
МАРE
RMSE
MAE
MAPE
RMSE
MAE
MAPE
RMSE
TTS
2.73
6.78%
5.21
3.10
8.02%
6.15
3.67
10.05%
7.46
STT
2.72
6.74%
5.19
3.07
7.94%
6.08
3.65
10.97%
7.37
CTTS
2.70
6.69%
5.20
3.05
7.93%
6.18
3.61
9.75%
7.51
CSTT
2.70
6.66%
5.22
3.06
7.92%
6.19
3.63
9.92%
7.48
CITRUS (Ours)
2.70
6.74%
5.14
2.98
7.78%
9.28%
6.85
Table 4: Training time (per epoch) and forecasting results vs. number of selected eig-eiv MetrLA.
2
k
10
k
50
100
150
200
MAE
2.80
2.74
2.72
2.71
2.71
2.70
2.70
Training time (seconds)
2.52
2.91
2.96
3.86
5.28
6.36
7.78
the CITRUS blocks are 16, while MLP modules have been used. Regarding NOAA, set the dimension and the initial linear layer of the CITRUS blocks as 16, and use 3-layer MIPs for channel mixing in each block. The CITRUS blocks have the last activation the loss function and consider root-normalized MSE (rNMSE) as the evaluation metric.
Leaky ReLU. We MSE as
Table illustrates the results of the weather forecasting task. We observe that CITRUS show superior forecasting performance compared to previous methods all forecasting horizons. We note that SGP model [42). which is scalable architecture, is the second-bes 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.
4.4 Ablation Study and Hyperparameter Sensitivity Analysis
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 regular GNN: ii) STT, which exactly the opposite of the TTS; iii) Continuous TTS (CTTS), where we replace the GNN in TTS with CGNN; and iv) Continous STT (CSTT), which the opposite of CTTS. We report the results of this ablation study on the MetrLA dataset in Table 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.
Hy Hyperparameter sensitivity analysis. The sensitivity analysis is related to the number of selected digenvector-eegenvale (cig-civ) pairs of the factor Laplacians. We select eigenvector-eigenaale pairs in k (2, 10, 50, 100, 150, 200) out of 207 components in the MetrLA dataset. Table presents the forecasting results in MAE along with the training time (per epoch) of this ablation study. We observe that selecting 50 weweeveeooo--igevvllle having fast training time. This experiment illustrates the advantages of CITRUS for accelerating training and keeping the performance as consistent as possible.
pairs is enough to get good forecasting performance while
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 separable continuous heat graph kernels solutions to the TPDEG. Therefore. 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 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-an 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 general form or for specific well-studied structures, such as ER graphs. Additionally, we plan investigate more challenging real- wolld applications of the proposed framework, especially in scenarios involving more than two factor graphs.
Acknowledgments
This work has been carried out the 4Cimmae Interdisciplinary Center (E4C) of IP Paris and Ecole Nationale des Ponts et Chaussées, which is in part supported by 3rd Programme Investissements d' Avenir VANR-18-EUR-OUA nique with the private sponsor Fonds Ifker pour le Climat, financed by Stéphane Agnès Ifker. This work was also supported by the center Hi! PARIS and ANR (French National Research Agency) under the JCJC project GraphIA (ANR-20-CE23-0099-01]
and by the Foundation of Ecole Polytech-
References
[1] Y. L. Xu, K. Konstantinidis, and P. Mandic, "Multi-graph tensor networks," Advances in Neural Information Processing Systems Workshop, 2020.
[2] Y. L. Xu, K. Konstantinidis, and D. P. Mandic, "Graph tensor networks: An intuitive frame- work for designing large-scale neural learning systems on multiple domains," arXiv preprint arXiv:2303. 13565, 2023.
[3] Monti, M. Bronstein, and X. Bresson, "Geometric matrix completion with recurrent multi- graph neural networks," in Advances in Neural Information Processing Systems, 2017.
[4] N. Kipf and M. Welling, "Semi-superrised classification with graph convolutional networks," in International Conference on Learning Representations, 2017.
[5] X. Song, K. Wu, and L. Chai, "Brain network analysis of schizophrenia patients based on hypergraph signal processing," " IEEE Transactions on Image Processing, vol. 32, pp. 4964- 4976, 2023.
[6] N. Shahid, F. Grassi, and P. Vandergheynst, "Tensor robust PCA on graphs," in IEEE Interna- tional Conference on Acoustics, Speech and Signal Processing, 2019.
[7] G. Ortiz-Jiménez, M. Coutino, S. P. Chepuri, and G. Leus, "Sparse sampling for inverse problems with tensors," IEEE Transactions on Signal Processing, vol. 67, no. 12, pp. 3272-3286, 2019.
[8] A. Cini, I. Marisca, D. Zambon, and C. Alippi, "Graph deep learning for time series forecasting." arXiv preprint arXiv:2310. 15978, 2023.
[9] G. Mishne, R. Talmon, R. Meir, J. Schiller, M. Lavzin, U. Dubin, and R. R. Coifman, "Hierar- chical
analysis for neuronal structure and activity pattern discovery." IEEE Journal of Selected Topics in Signal Processing, vol. 10, no. 7, pp. 1238-1253. 2016.
[10) M. Sabbagi and E. Isufi, "Graph-time convolutional neural networks: Architecture and theoreti- cal analysis," " IEEE Transactions on Pattern Analysis and Machine Intelligence, vol. 45, no. 12, pp. 14625 14638, 2023.
[11] A. Einizade and S. H. Sardouie, "Learning product graphs from spectral templates." IEEE Transactions on Signal and Information Processing over Networks, vol. 9, pp. 357-372, 2023.
[12] A. Einizade, S. Nasiri, S. H. Sardouie, and G. D. Clifford, "ProduetGraphhSleepNer Sleep staging using product spatio-temporal graph learning with attentive temporal aggregation," Neural Networks, vol. 164, pp. 667-680, 2023.
10

---

[13] A. Ortega, P. Frossard, J. Kovačević, J. M. Moura, and P. Vandergheynst. "Graph signal processing: Overview, challenges, and applications." Proceedings of the IEEE, vol. 106, no. pp. 808-828, 2018.
[14] K. Oono and T. Suzuki, "Graph neural networks exponentially lose expressive power for node classification," in International Conference on Learning Representations, 2020.
[15] C. Cai and Wang. "A note on over-smoothing for graph neural networks," in International Conference on Machine Learning Workshop, 2020.
[16] J. H. Giraldo, K. Skianis, T. Bouwmans, and F. D. Malliaros, "On the trade-off between over-smoothing and over-squashing in deep graph neural networks," " in ACM International Conference on Information and Knowledge Management, 2023.
[17] M. Behmanesh, M. Krahn, and M. Ovsjanikov, "TIDE: Time derivative diffusion for deep learning on graphs," in International Conference on Machine Learning, 2023,
[18] E. Isufi and G. Mazzola, "Graph-time convolutional neural networks," in IEEE Data Science and Learning Workshop, 2021.
[19] J. A. Castro-Correa, J. H. Giraldo, M. Badiey, and F. D. Malliaros, "Gegenbauer graph neural networks for time-varying signal reconstruction." IEEE Transactions on Neural Networks and Learning Systems, vol. 35, no. 9, pp. 11734-11745, 2024.
[20] A. Han, D. Shi, L. Lin, and J. Gao, "From continuous dynamics to graph neural networks: Neural diffusion and beyond," Transactions on Machine Learning Research, 2024.
[21] L.-P. Xhonneux, M. Qu, and J. Tang, "Continuous graph neural networks," in International Conference on Machine Learning, 2020.
[22] N. Sharp, S. Attaiki, K. Crane, and M. Ovsjanikov, "DiffusionNet: Discretization agnostic learning on surfaces," ACM Transactions on Graphics, vol. 41, no. 3, pp. 1-16, 2022.
[23] A. Sandryhaila and J. M. Moura, "Big data analysis with signal processing on graphs: Repre- sentation and processing of massive data sets with irregular structure," " IEEE Signal Processing Magazine, vol. 31, no. 5, pp. 80-90, 2014.
[24] R. Hammack, W. Imrich, and S. Klavvar, Handbook of product graphs. CRC press, 2011.
[25] E. Isufi, A. Loukas, N. Perraudin, and G. Leus, "Forecasting time series with varma recursions on graphs," IEEE Transactions on Signal Processing, vol. 67, no. 18, pp. 4870-4885, 2019.
[26] J. Gao and B. Ribeiro, "On the equivalence between temporal and static equivariant graph representations," in International Conference on Machine Learning. 2022.
[27] S. K. Kadambari and S. P. Chepuri, "Product graph learning from multi-domain data with sparsity and rank constraints," IEEE Transactions on Signal Processing, vol. 69, pp. 5665-5680, 2021.
[28] T. G. Kolda and B. W. Bader, "Tensor decompositions and applications." " SIAM Review, vol. 51, no. 3, pp. 455-500, 2009.
[29] N. D. Sidiropoulos, L. De Lathauwer, X. Fu, K. Huang, E. E. Papalexakis, and C. Faloutsos, "Tensor decomposition for signal processing and machine learning," IEEE Transactions Signal Processing, vol. 65, no. 13, pp. 3551-3582, 2017.
[30] Y. Song, Q. Kang, S. Wang, K. Zhao, and W. P. Tay, "On the robustness of graph neural diffusion to topology perturbations," in Advances in Neural Information Processing Systems, 2022.
[31] F. D. Giovanni, J. Rowbottom, B. P. Chamberlain, T. Markovich, and M. M. Bronstein, "Un- derstanding convolution on graphs via energies," Transactions on Machine Learning Research, 2023.
[32] T. K. Rusch, M. Bronstein, and S. Mishra, "A survey oversmoothing in graph neural networks." SAM Research Report, 2023.
[33] M. Indibi and Aviyente, "Low-rank and smooth tensor recovery on Cartesian product graphs,"
in IEEE International Conference Sampling Theory and Applications, 2023.
[34] J. S. Stanley, E. C. Chi, and G. Mishne, "Multiway graph signal processing on tensors: Inte- grative analysis of irregular geometries," IEFE Signal Processing Magazine, vol. 37, no. pp. 160-173, 2020.
11

---

[35] Y. Li, R. Yu, C. Shahabi, and Y. Liu, "Diffusion convolutional recurrent neural network: Data- driven traffic forecasting." in International Conference on Learning Representations. 2018.
[36] D. I. Shuman, S. K. Narang, P. Frossard, A. Ortega, and Vandergheynst, "The emerging field of signal processing on graphs: Extending high-dimensional data analysis to networks and other irregular domains," IEEE Signal Processing Magazine, vol. 30, pp. 83-98, 2013.
[37] A. Cini and I. Marisca, "Torch Spatiotemporal," 3 2022.
[38] Z. Wu, S. Pan, G. Long, J. Jiang, and C. Zhang. "Graph WaveNet for deep spatial-tempord graph modeling," " in International Joint Conference on Artificial Intelligence, 2019.
[39] C. Zheng, X. Fan, C. Wang, and J. Qi, "GMAN: graph multi-attention network for traffic prediction," in AAAI Conference on Artificial Intelligence, 2020.
[40] B. Yu, H. Yin, and Z. Zhu, "Spatio-temport"" graph convolutional networks: a deep learning framework for traffic forecasting." in International Joint Conference on Artificial Intelligence, 2018.
[41] L. Ruiz, Gama, and A. Ribeiro, "Gated graph recurrent neural networks." " IEEE Transactions on Signal Processing, vol. 68, pp. 6303-66 6318, 88, 2020.
[42] A. Cini, I. Marisca, F. M. Bianchi, and C. Alippi, "Scalable spatiotemporal graph neural networks," in AAAI Conference on Artificial Intelligence, 2023.
[43] B. Girault, "Stationary graph signals using isometric graph translation," in IEEE European Signal Processing Conference, 2015.
[44] A. Arguez, Durre, S. Applequist, R. Vose, M. Squires, X. Yin, R. R. Heim, and T. W. Owen, "NOAA's ss 1981-2020 US climate normals: an overview," Bulletin of the American Meteorological Society, vol. 93, no. 11, pp. 1687-1697, 2012.
[45] K. Petersen, M. Pedersen, et al., "The matrix cookbook." " Technical University of Denmark, vol. 7, no. 15, p. 510, 2008.
[46] Z. Zhou, Q. Huang, G. Lin, K. Yang. L. Bai, and Y. Wang, "GReTo: Remedying dynamic graph
topology-task discordance via target homophily," in International Conference Learning Representations, 2023.
12

---

Appendix
The appendix contains the proof of the claimed theoretical statements (Section A), descriptions of the compared methods (Section C the standard deviation of the results (Section D) homophily-
heterophily trade-off (Section Eb. the case of inexact Cartesian product graphs (Section F. how to choose the appropriate number of eigenvalue-eigenveede tive fields on the over-smoothing (Section and hyperparameter details (Section respectively.
pairs (Section G), the effect of graph recep-
A Proofs
A.1 Proof of Theorem 3.2
Proof. For obtaining the derivative of au, from e thp (for p I..... P) in 44 we first express 8t
au, from the claimed solution in (4) as follows:
OU,
Σ p=1
(19)
Ōt
where
U U. == Up x1 e -te x Xp-1 6 -pp- Xp ( (~Ipe =bp) xpt1 e -lepts xpt2... Xpe ebbr
(20)
Now, by performing mode-p unfolding operation on Upp we have:
U (.p)))) - L dpp Ln(g) (e-tkr -lr - @ e...e chpt 1e.tp-i....000 (21)
On the other hand:
(( c (p) (e (e . thp+ 1 @ 8... 0..
(22)
Therefore, by combining 21 and [22 one can write:
U UP = ..))
(23)
which is the mode-p unfolding of the following equation:
U = -- Xp Lp.
(24)
and the proof is completed.
A.2 Proof of Theorem 3.3
Proof. First, we prove the following lemma about the separability the heat kernels define on product graphs w.r.t. the factor graphs:
Lemma A.1. The exponential graph kernel e.e on a Cartesian product graph L. @pmIIp can be Kratecker-factored into its factor-based graph kernels as e e
Proof. Using the pairwise summation property of the factor eigenvalues in Cartesian product, one can write:
cL =((
Vp) ( Vp) = (
(25)
= (of Vp) (8f (n) (8fm
e
where Vp. X(2), and Ap are the eigenmatrix, ith eigenvalue, and the diagonal eigenvalue matrix corresponding to the pth Laplacian Lp. Note that, the related literature, this property has also been proved from different points of view, e.g., (34]
D
13

---

Next, by performing mode-(P 1) unfolding on both sides of ( one can write:
(C((U)))(+4) = W, U (P)+)
0 0
ubp
e
(26)
((G(U)))))
[ C b 4p
U
U(441) Wi.
A.3 Proof of Proposition 3.6
Proof. The proof is presented by construction. Therefore, we first show it is true for the case of P = 2 as:
Á = A @ 1 1 A2
(A1 E1 1 2 (A2 + E2)
E
(27)
AA & +1 I @ A2) + E1 X I I 8
A E.
Then, by assuming the theorem holds for the case of P = K and the definitions of A's = 0b Ap
and Axx Bb Ap. we next show it also holds for p K as follows:
À = KKIA A, (0p GAKHI A @ 12 + & AKtI
(A's EE & 12 +1 I @ (A + Ek-1)
(28)
A= Ap
E= EEE Ex
(A & L2 + I @ AR+1) + (Ex @ I2 + I & EK+1)
A.
For proving the upper bound of the norm of E, in a similar approach to the previous proof, first prove for the case P 2 as:
E|l
E
I
Eall S E @ L
I & Eg
(29)
E2 = E + E2.
Using a similar approach to the proof in 28 we assume the theorem holds for P = K as E <
ΣK- Ep. and, based on this information, then, prove the theorem for the of P K + 1 as follows:
B|| E|
I
EK+1
<
E
Iall
Ex+ill |||
ΣK Ep
K+1
(30)
Ex
III
1.
L
E
Ex
Σ
Ep:
p=1
which concludes the proof.
A.4 Proof of Theorem
Proof. The poof can straightforwardly obtained by applying Proposition (30] based on the obtained results in Proposition 3.6
D
14

---

A.5 Proof of Lemma 3.9
Proof. By construction, we first prove the theorem for the case of P 2. The defined Laplacian Ż in 15) can be rewritten as:
Å
A
-=I-A=I-
12
+1
@
2
2
A 1 11841100,,
21
[T - (A,
& + (I s A2ll)
2
(
+ 1 - (I1 @ Aa)))
=
2
(31)
(
(I2 (12-)2))
A2))
L
+ L2
L 2
L2
+
where used the following property of Kronecker products: a<A & B) (aA) @ B A (aB)
45 Next, by assuming that the theorem holds for P K and the definition of A 1 and
p= A
= op
we prove for P K + 1:
K
i=I-A=I
(KÅ' + AKT1)
K
K+1
KAK
AK+1
& II I 8
K+I
K K+1
1
(KK +[]++11---- - (KÅ's @ 12 @ AktI)
K+1
1
[K[I-(A) @ 12) + [ (1 @ Ak++)))
K+1
1
(11 - (A& @ 12)
+ I 8 l2
(I 8 A2)
K+1
(32)
1
[K]
(I2 A Lg
+ 1 8 (12 Az)
K K+1
1
K
12+ho 12+ho
K +1
K
1
(t.p)
K +
p=1
Lp
+ p=1
K +11
Besides, by the addition rule for the eigenvalues of a product graph adjacency or Laplacian (A ( + x(2)) ). it can be easily seen that the cigenvalues of the normalized product adjacency (W) and Laplacian (L) (3) are in the intervals of [-1, and [0, 2]. respectively, quite similar to the case for the normalized factor graphs, which concludes the proof.
15

---

A.6 Proof of Theorem 3.10
Proof. First, Note that Using as the Graph Fourier Transform (GFT) 13 of x w.r.t É (with eigenvalues Jem)) one can write 15
E(x) = x Lx Σ A
(33)
Next, by defining À as the smallest non-zero eigenvalue of the Laplacian L, the following lemma characterizes the over-smoothing aspects of applying heat kernel in the simplest case.
Lemma A.2. We have:
Ele xx) S 5x
(34)
Proof. By considering the EVD forms of VAVT and -- Ve VVT, one can write
E(e xx)
Ve Ve
VAV
T
(35)
N Σ
=x
L
x Σ A
Se-2
Exx).
Note that in the above proof, we ruled out the zero eigenvalues since they are useless analyzing Dirichlet energy.
Then by considering the following lemmas from 15
Lemma A.3. (Lemma 3.2 in (15).. ((XW) s | | (X).
Lemma A.4. (Lemma 3.3 in 15 ). For ReLU and Leaky-ReLU nonlinearities E(o(X)) S E(X).
Therefore, by combining the previous concepts and considering that the minimum non-zero eigenvalue of L (4) in T6 is (((m) (m) with m arg min, {() ( it can be said that:
Theorem A.5. For any I N4, we have (X): S sie -iA where 81 := II 8th and 8th is the square of the maximum singular value of WT. Besides, X X(m) and = ((m) with arg min, ((0) (1), where (3) is the smallest non-zero eigenvalue of L.
Afterward, we can state our final corollary as follows:
Corollary A.6. Let == SUPIEN, 8t. We have:
E(X") () < s'e E(X) = (1ns- E(X).
(36)
So, F)X exponentially converges to 0, when limy e (1us-2i8) = 0, ie., Ins -75 < 0.
B Experiments on more than Two Factor Graphs
We consider the node regression task similar to the settings described in Section 4.1 but for three
ER graphs with PiR (1) (1) 0.3, and varying PER e {0.1, 0.3) to monitor the performance in Pix
(3)
two different connectivity scenarios. Note that here we do not perturb factor graphs but, to make the experiment more challenging, the SNR on the graph data is set to SNR=0. Table Provide provides the results across a varying number of layers, which shows that the proposed framework's 's performance is (3) 0.1. On the other resistant to adding layers, especially in the case of more sparse graphs, i.e., PER hand, the best results were obtained in 1, 2, 4, which are the nearest numbers to the actual number of layers in the true generative process which was 3. These observations validated by comparing the results with the GCN, which the over-smoothing and performance degradation happen severely faster than CITRUS.
16

---

Table 5: Experiments on more than two factor graphs (PER
(1)
PPx
(2)
0.3).
=
GCN, PeR
0.185::0.12
0.241
0.371:10.32
0.458:0.132
0.49410.147
0.552::0.118
CITRUS. Pm
(5) 0.1
0.141=00075
0.137=0.082
0.150=0.091
0.198 106
0.215=0.124
0.276=0.120
GCN, Pur 0.3
0.253-+0.21
0.296+0.113
0.365+0.123
0.429+0.097
0.517+0.143
0.67+00.53
CITRUS. PER
(3)
0.3
0.18840.094
0.182=0.095
0.19310.096
0.23010.089
0.242:10.096
0.336-0.107
Table 6: Standard deviation of traffic forecasting comparison in Table
MetLA
H 3
6
H 12
Method
MAE
MAPE
RMSE
MAE
MAPE
RMSE
MAE
MAPE
RMSE
CITRUS (Ours)
0.015
0.009
0.008
0.007
0.012
0.003
0.002
0.005
0.004
PemsBay
H = 3
6
H = 12
Method
MАE
МAPE
RMSE
MAE
MAPE
RMSE
MAE
MAPE
RMSE
CITRUS (Ours)
0.009
0.008
0.003
0.013
0.006
0.005
0.008
0.011
0.008
C Descriptions of the Compared Methods
uses Kalman filter to implement auto-regressive integrated moving average framework, but processes data in a timestep-wise manner.
ARIMA 35
generalizes the vector auto-regressive moving average (VARMA) the case of considering graph-based data by imposing some statistical assumptions on the model and relying on the spatial graph.
G-VARMA (25
GP-VAR 25 : adapts AR model by spatial graph polynomial filters, with fewer number of learnable parameters.
FC-LSTM 35);
a hybrid of fully connected MLPs and LSTM cells that considers each time step separately by assigning one LSTM cell them.
Graph WaveNet (38 hybrid adapted convolutional and attentional network with GNNs for ST forecasting.
GMAN 39
exploits multiple ST attention layers within an encoder-decoder framework.
STGCN
40
hybrid typical GCNs and ID convolutional blocks for ST forecasting.
GGRNN
generalizes the typical RNNs by adapting GCN blocks instead linear
MLPs.
TTS framework, which first processes the temporal information by GRU cell and then learns spatial representations by applying typical GCN blocks.
GRUGCN 26
SGP scalable graph-based ST framework for considerably reducing the number of learnable parameters, Note that we used the version without positional embedding blocks.
D Standard Deviation of the Results
The standard deviation of traffic forecasting in Table and weather forecasting regarding the proposed framework in Table 2 are provided in Tables 6 and
17

---

Table 7: Standard deviation of weather forecasting comparison in Table
Molene
NOAA
Method
H =
H = 2
H =
3
H =4
H =
H
H =
H =
3
H = 4
H = 5
CITRUS (ours)
0.0174
0.0155
0.0132
0.0107
0.0086
0.0034
0.0030
0.0026
0.0024
0.0023
Table 8: Intra and nnter-hompphily measures on the real-world graphs of MetrLA and PemsBay datasets.
Intra-graph homophily क.;
Inter-graph homophily T..
Dataset
Ps
4p
9n
Ps
qn
MetrI.A
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
E Homophily-eteroohilly Trade-off across the Studied Real Datasets
Regarding the honophily-heteroohily trade-off, first, we should mention that measuring the homophily-heterophhly 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 inter- graph (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 nomophily-heterophily heterophily behaviors. We metrics (full details in (46]) in Table X As observed
provided the measured this table, due the wide range of variability across the used datasets, we have actually evaluated our performance on various datasets in terms of homophily-heterophin 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.
cemoopily-heteropee
The Case of Inexact Cartesian Product Modeling
This scenario lies under the general umbrella of inaccurate adjacencies 33), 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
24). Consider the true graph is actually Strong product graph A A & A2+ A @ 1 & A2, and one mistakenly considers it as Cartesian product À A & I Ø Ag.
Then, E = À = A & A2, where ||E|| (I)
X x12) So, the approximation error depends on the multiplication of the factor graph eigenspaces. Similarly, in the case of Kronecker product graphs
À A₁ Ag, it obtains E|| (1) (2) + X (2) So, in this case, the summation of the factor eigenspaces matters too. Briefly, the approximation error bound depends on the factor graph spectrums.
G Choosing Appropriate k
We can choose appropriate i 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., ||vAv 7 for and À being the eigenvectors and eigenvalues) on the MetrLA dataset. We observe 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.
18

---

Explained variance ratio
Individual explained variance
0.015
0.010
0.005
0.000
50
100
200
Principal components
Figure 4: Explained variance ratio vs. selected principal components.
Table 9: Details of the training settings and hyperparameters, i.e., T Gauto-regressive order), emb (dimension of embedding size in the spatiotemporal encoder), hid (dimension of linear mapping size in the spatiotemporal encoder), PMPP (dimension of linear mapping size in the MLP layers), Victtrus (number of CITRUS blocks), F (second dimension of W, in (6)) k; (number of selected sigenvalue-eigenvectoo pairs in the i-th factor graph), Ir (learning rate), Ybbachh (batch size), and repochs (number of epochs).
Dataset
T
emb
hid
FMPr
ACITRUS
F
h1
Ir
Thbach
Reepocas
MetrLA
6
16
32
64
3
64
205
4
0.01
2048
300
PemsBay
6
16
32
64
3
64
323
0.01
2048
300
Molene
10
16
16
3
16
30
0.001
64
1500
NOAA
10
4
4
3
4
107
8
0.01
256
400
H
The Effect of the Graph Receptive Fields on Controlling Over-smoothing
From a theoretical point of view, based on eq. (88 one intuitive way to alleviate or at least slow down the over-smoothing phenomena is to keep In ži positive, which tends to f < P. In s. 2x
This clarifies why one should not increase the node-neighte seade--eigbohhood especially in the case of facing connected graph or increasing the number of layers. This finding is also compatible with the pioneering research about the over-smoothing concept, egg., (14). From an experimental point of view, and if we consider the graph receptive field as hyperparameter, we have designed an experiment to validate the claimed statement. Here, we vary & E (0.1, 5, 10, 20) and monitor the over-smoothing process in Figure 5 In this figure, b = P. sns, and as it is observed, for values of higher than b the 25,
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 15415) has also been mentioned as 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 < b, one might need to perform other kinds of analyses, e.g., over-squashing analysis.
Hyperparameter Details
The detailed hyperparameters (optimized by cross- validatoii on the validation data) and/or training settings of CITRUS i.e., T reuts-regr order), emb (dimension of embedding size in the spatiotemporal encoder), hid (dimension of linear mapping size in the spatiotemporal encoder), FML (dimension linear mapping size in the MLP layers), VeCTrUS (number of CITRUS blocks), F (second dimension of W, in (6). k; (number of selected eigenvalue-eigenvetto pairs in the i-th factor graph), Ir (learning rate), Tibatch (batch size), and Nepochs (number of epochs), are provided in Table Full details can be found in the implementation codes in https//githbbcom////tttaaaaa///a
19

---

b 7.08
0
Log Relative Distanee
-20
40
880
100
0.0
2.5
5.0
7.5
10.0
Number of layers
0.1
20 f=
10
Figure 5: Log relative distance vs. increasing the number of layers for different values of t actual bounds.
J
Experiments Compute Resources
All the experiments were run on one GTX A100 GPU device with 40 GB of RAM.
K
Broader Impacts
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 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.
20

---

NeurIPS Paper Checklist
1. Claims
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 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 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 long is clear that these goals are not attained by the paper.
2. Limitations
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 separate paragraph for planning future work.
Guidelines:
The answer NA means that the paper has 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 (c.g.. independence assumptions, noiseless settings, model sellsspeciiiationn 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 few datasets or with 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, facial recognition algorithm may perform poorly when image resolution low or images are taken in low lighting. Or speech-to-tex system might not be used reliably to provide closed captions for online lectures because it fails 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 faimesss.
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 't acknowledged in the paper. The authors should use their best judgment and recognize that individual actions in favor of transparency play an impor- tant role in developing norms that preserve the integrity of the community, Reviewers will be specifically instructed to not penalize honesty concerning limitations.
3. Theory Assumptions and Proofs
Question: For each theoretical result, does the paper provide the full set of assumptions and complete (and correct) proof?
21

---

Answer: [Yes]
Justification:
The Appendix provides all the needed proofs.
Guidelines:
The answer NA means that the paper does not include theoretical results.
All the theorems, formulas, and proofs in the paper should be numbered and cross- referenced.
All assumptions should be clearly stated or referenced in the statement of any theorems.
The proof's can either appear in the main paper or the supplemental material, but if they appear in the supplemental material, the authors are encouraged to provide a short proof sketch provide intuition.
Inversely, any informal proof provided in the core of the paper should complemented by formal proofs provided in appendix or supplemental material.
Theorems and Lemmas that the proof relies upon should be properly referenced.
4. Experimental Result Reproducibility
Does the paper fully disclose all the information needed to reproduce the main ex- perimental results the paper to the extent that affects the main claims and/or conclusions of the paper (regardless of whether the code and data are provided not)?
Question:
Answer: [Yes]
Justification: Section 4 and the Appendix contain any needed detail for reproducing the results the proposed approach.
Guidelines:
The answer NA means that the paper does not include experiments.
If the paper includes experiments, a No answer to this question will not be perceived well by the reviewers: Making the paper reproducible is important, regardless of whether the code and data are provided or not.
If the contribution is dataset and/or model, the authors should describe the steps taken to make their results reproducible or verifiable.
Depending on the contribution, reproducibility can be accomplished in various ways. For example, if the contribution is a novel architecture, descríbing the architecture fully might suffice, or if the contribution is specific model and empirical evaluation, it may be necessary to either make possible for others to replicate the model with the same dataset, provide access to the model. In general. releasing code and data is often one good way to accomplish this, but reproducibility can also be provided via detailed instructions for how to replicate the results, access to a hosted model (e.g.. g., in the case of a large language model), releasing of a model checkpoint, or other means that are appropriate to the research performed.
While NeurIPS does not require releasing code, the conference does require submis- sions to provide some reasonable avenue for reproducibility, which may depend on the nature of the contribution. For example
(a) the contribution is primarily a new algorithm, the paper should make clear how reproduce that algorithm. (b) If the contribution is primarily new model architecture, the paper should describe the architecture clearly and fully.
(c) the contribution is new model (e-g., a large language model), then there should either be way to access this model for reproducing the results or way to reproduce the model (e.g., with an open-source dataset or instructions for how to construct the dataset).
(d) We recognize that reproducibility may be tricky in some cases, in which case authors are welcome describe the particular way they provide for reproducibility. In the case of closed-source models, it may be that access to the model is limited in some way (e.g., to registered users), but it should possible for other researchers to have some path to reproducing verifying the results.
5. Open access to data and code
22

---

Does the paper provide open access to the data and code, with sufficient instruc- tions faithfully reproduce the main experimental results, as described in supplemental material?
Question:
Answer: [Yes]
Justification: We have uploaded the reproducible codes alongside paper and 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.
Please see the NeurIPS code and data submission guidelines (https://nips.co cc/
for more details.
While we encourage the release of code and data, we understand that this might not be possible, "No" is an acceptable answer. Papers cannot be rejected simply for not including code, unless this is central to the contribution (e.g., for a new open-source benchmark).
The instructions should contain exact command and environment needed to run to reproduce the results. See the NeurlPS code and data submission guidelines (https
y /nips. /nips.cc//ppppp https//cccccccguuddeeeeeee/ee//e huabbiiccguidde/odeeeuubiiii
for more details.
The authors should provide instructions on data access and preparation, including how access the raw data, preprocessed data, intermediate data, and generated data, etc.
The authors should provide scripts to reproduce all experimental results for the new proposed method and baselines. only subset of experiments are reproducible, they should state which ones are omitted from the script and why.
At submission time, to preserve anonymity, the authors should release anonymized versions (if applicable).
Providing as much information as possible in supplemental material (appended to the paper) recommended, but including URLs to data and code permitted.
6. Experimental Setting/ettails
Question:
Does the paper specify all the training and test details (e.g., data splits, hyper- parameters, how they were chosen, type of optimizer, etc.) necessary 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 the core the paper to 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, as supplemental material.
7. Experiment Statistical Significance
Does the paper report error bars suitably and correctly defined or other appropriate information about the statistical significance of the experiments?
Question:
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 emor bars, confi- dence intervals, or statistical significance tests, least for the experiments that support the main claims of the paper.
23

---

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 the standard error of the mean.
It is OK to report 1-sigma error bars, but one should state it. The authors should preferably report a 2-sigma error bar than state that they have 96% CI, if the hypothesis of Normality of errors is not verified.
For asymmetric distributions, the authors should be careful not show in tables figures symmetric error bars that would yield results that are out range (e.g. negative error rates).
If error bars are reported in tables or plots, The authors should explain the text how they were calculated and reference the corresponding figures or tables in the text.
8. Experiments Compute Resources
For each experiment, does the paper provide sufficient information on the com- puter resources (type of compute workers, memory, time of execution) needed reproduce the experiments?
Question:
Answer: [Yes]
Justification:
We have reported the used compute resources in the Appendix, Section E.
Guidelines:
The answer NA means that the paper does not include experiments.
The paper should indicate the type of compute workers CPU or GPU, internal cluster, or cloud provider, including relevant memory and storage.
The paper should provide the amount compute required for each of the individual experimental runs as well estimate the total compute.
The paper should disclose whether the full research project required more compute than the experiments reported in the paper (e.g., preliminary or failed experiments that didn't make into the paper).
9. Code Of Ethics
Question: Does the research conducted the paper conform, in every respect, with the NeurlPS Code of Ethics https: : /ttee://eeuuiieeeee
Answer: [Yes]
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 NeurlPS Code of Ethics.
If the authors answer No, they should explain the special circumstances that require deviation from the Code of Ethics.
The authors should make sure to preserve anonymity (e.g., there a special consid- cration due to laws or regulations their jurisdiction).
10. 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 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), faimess 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 tied to particular applications, let alone deployments. However, if there is 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 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 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 system learns from feedback over time, improving the efficiency and accessibility of ML).
11. Safeguards
Does the paper describe safeguards that have been put in place for responsible release of data models that have high risk for misuse (e.g., pretrained language models, image generators, or scraped datasets)?
Question:
Answer: [NA]
Justification: Our paper poses no such risks.
Guidelines:
The answer NA means that the paper poses such risks.
Released models that have 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 usage guidelines or restrictions 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 challenging, and many papers do not require this, but encourage authors to take this into account and make 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
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 used and, 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.
25

---

If assets are released, the license, copyright information, and terms of use in the thcode. has curated licenses for some datasets. Their licensing guide can help determine the license of dataset.
package should be provided. For popular datasets, paperswi
batasets.papersstt
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 encouraged to reach out the asset's s creators.
13. New Assets
Question: Are assets introduced in the paper well documented and is the documentation provided alongside the assets?
Answer: [Yes]
Justification: The only assets of our work the implementation codes with detailed Readme files.
Guidelines:
The answer NA means that the paper does not release new assets.
Researchers should communicate the details of the dataset/code/moodl as part of their submissions via structured templates. This includes details about training, license, limitations, etc.
The paper should discuss whether and how consent was obtained from people whose asset is used.
At submission time, remember to anonymize your assets (if applicable). You can either create an anonymized URL or include an anonymized zip file.
14. Crowdsoureing and Research with Human Subjects
Question: For crowdsourcing experiments and research with human subjects, does the paper include the full text of instructions given to participants and screenshots, if applicable, as well details about compensation (if any)?
Answer: [NA]
Justification: Our paper does not involve crowdsourcing nor research with human subjects.
Guidelines:
The answer NA means that the paper does not involve crowdsourcing nor research with human subjects.
Including this information in the supplemental material is fine, but if the main contribu- tion of the paper involves human subjects, then much detail as possible should be included the main paper.
According to the NeurlPS Code of Ethics, workers involved in data collection, curation, other labor should be paid at least the minimum wage in the country of the data collector.
15. Institutional Review Board (IRB) Approvals or Equivalent for Research with Human Subjects
Question: Does the paper describe potential risks incurred by study participants, whether such risks were disclosed to the subjects, and whether Institutional Review Board (IRB) approvals (or an equivalent approval/review based on the requirements of your country or institution) were obtained?
Answer: [NA]
Justification: Our paper does not involve crowdsourcing nor research with human subjects.
Guidelines:
The answer NA means that the paper does involve crowdsourcing nor research with human subjects.
Depending the country in which research is conducted, IRB approval (or equivalent) may be required for any human subjects research. If you obtained IRB approval, you should clearly state this in the paper.
26

---

We recognize that the procedures for this may vary significantly between institutions and locations, and we expect authors adhere 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.
27


Reviews:
[{"review_id": "Ut7OjJAuqV", "reviewer": "Reviewer_wQ9d", "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.", "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": "The authors address the limitations of their work as the last paragraph of their paper.", "rating": 6, "confidence": 2, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "* The paper is well written.\n* The proposed methods are backed by theory.\n* The experimental results are SOTA.\n* Ablation study is provided.", "weaknesses": "See questions.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "One example where real-world graphs cannot be represented as a Cartesian product of factor graphs is the case where we have time-varying dynamic graphs, which might lead to some approximation errors by assuming the resulting graph is a Cartesian product graph. Indeed, this comment relates to the general issue of inaccurate adjacencies, denoted as $\\tilde{A}=A+E$, which we studied in the stability analysis of CITRUS in Section 3.2. Therefore, the approximation error is bounded in Theorem 3.7.\n\nAnother case is when the real-world graphs follow a different product graph operation. For example, let's consider two well-known graph products: the Strong product and the Kronecker product. Suppose the true graph is actually a Strong product graph, $A=A_1\\otimes A_2+A_1\\otimes I_2+I_1\\otimes A_2$, but it is mistakenly treated as a Cartesian product $\\tilde{A}=A_1\\otimes I_2+I_1\\otimes A_2$. In this case, the error $E$ is given by $E=A-\\tilde{A}=A_1\\otimes A_2$, where $||E||=\\lambda^{(1)}\\_{max}\\lambda^{(2)}\\_{max}$. Therefore, the approximation error depends on the multiplication of the maximum eigenvalues of the factor graphs.\n\nSimilarly, if we consider Kronecker product graphs, $\\tilde{A}=A_1\\otimes A_2$, we find that $||E||\\le\\lambda^{(1)}\\_{max}+\\lambda^{(2)}\\_{max}+\\lambda^{(1)}\\_{max}\\lambda^{(2)}\\_{max}$.  In this case, the summation of the maximum eigenvalues of the factor graphs matters too.\n\nIn summary, the bound on the approximation error depends on the spectrum of the factor graphs.", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "66u2mJOeSj", "reviewer": "Reviewer_gs8b", "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.", "questions": "1. Could the authors provide the standard deviation for the experimental results reported in table 1-4?\n2. What are typical problems involving more than two factor graphs?\n3. 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.\n4. It is also interesting that CITRUS seems to outperform especially in longer horizon predictions. Any investigation why it is the case?\n5. Any guided way to choose k? Also any idea why a very small k seems to be a good enough approximation?", "limitations": "n/a", "rating": 6, "confidence": 3, "soundness": 3, "presentation": 2, "contribution": 3, "strengths": "1. The paper is overall well written and easy to follow (except some notations could be better defined).\n2. Handling multi-domain graph data is an interesting and challenging problem in practice.\n3. Learning joint (and not sequential) spatiotemporal couplings by modeling these dependencies using product graphs with learnable receptive fields is an interesting idea that is shown to be effective.\n4. The evaluation is done with respect to a diverse set of baselines.", "weaknesses": "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.\n- For example, what is purpose of having tilde in (3) and (4)? Then it disappeared in (5).\n- Is lower bar for $U$ necessary? why not just go with $U$?\n- is $\\times_i$ mode-i tensorial multiplication? what is the precise definition?\n- One minor thing --- $:=$ is used in line 120 but not line 123. Please check for consistency of other notations as well. \n2. 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?", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "**W1**. We’ve significantly simplified the notations in our paper. For example, we removed the unnecessary tildes from the notation in the revised manuscript.\n- Regarding the lower bar for tensors, we find it important to differentiate tensors (higher-order data) from regular matrices as in [9], but, if you have found it hard to read, we will think of alternating notations, like in [1].\n- For the tensor product, this is an extension of the matrix multiplication operation. For example, on 3D tensors, $\\underline{X}=\\underline{G}\\times\\_1 A$ is equal to $\\underline{X}\\_{(1)}=A\\underline{G}\\_{(1)}$, where the matrix $\\underline{G}_{(1)}$ is obtained by concatenating the mode-1 slices of the tensor $\\underline{G}$. We’ve included a detailed description in the Appendix of the revised manuscript. For more information please refer to Sections 2.4-5 in [1].\n- In line 120, we have the definition of the product Laplacian, but the equations in line 123 are EVD forms of Laplacians, not definitions themselves.\n\n**W2**.\n- In Fig. 3, we forgot to mention that the y-axis is on a logarithmic scale. Indeed, both lines are very close in non-logarithmic scales. We’ve corrected the label on the y-axis. Regarding the divergence, please note that our theorem states that $\\ln{(\\frac{E(X_l)}{E(X_0)})}\\le l(\\ln{s}-\\frac{2}{P}t\\lambda)$. Therefore, when $\\ln{s}-\\frac{2}{P}t\\lambda$ is positive, the RHS here grows with an increase in the number of layers. On the other hand, there is no problem with $E(X_l)$ going to zero because the theorem holds and is validated by the experimental results. However, in this case, the bound is looser than the case of $\\ln{s}-\\frac{2}{P}t\\lambda<0$.\n- Theorem 3.9 aims to characterize and provide insights about the over-smoothing aspects of CITRUS. To see the effect of all factor graphs, by using Lemmas A.2-A.4, we can also obtain $E(X_{l})\\le e^{l\\left(\\ln{s}-\\frac{2\\sum_{p=1}^{P}{t^{(p)}\\lambda^{(p)}}}{P}\\right)}E(X_{0})$. Therefore, we observe that over-smoothing is affected by the weighted average of the factor graph spectra, weighted by their receptive fields. We have provided further over-smoothing diagrams for varying sparsity levels for factor ER graphs in the **uploaded PDF** on the main rebuttal (Fig. R1). We observe that the denser the factor graphs (i.e., higher values for edge probability $p$), the faster the convergence to the overall over-smoothing state. This addresses the reviewer’s concern about the factor graph effects on over-smoothing. This aligns with results in the literature, e.g., [2], about a higher possibility of over-smoothing with denser graphs. We’ve added this discussion as a Remark in the Appendix of the revised manuscript.\n\n**Q1**. The standard deviations (STDs) have been provided in the **uploaded PDF** (Tab. R1) on the main rebuttal. Additionally, please note that we included some STDs in Tables 5, 6, and 7 in the Appendix.\n\n**Q2**. There are several possible examples of product graphs with more than two-factor graphs. For example, video data can be represented as a three-factor graph (width, height, and time). Similarly, in sleep staging with brain signals, there are multiple dimensions, such as temporal sleep windows, the spatial dimension of electrodes on the head, the frequency domain, and the time domain [3].\n\nAnother possible application could be in recommendation systems, where we have an item graph, a user graph, and a feature space of the user-item elements.\n\n**Q3**. Cartesian product graphs are useful for processing spatiotemporal data because we can model it as the Cartesian product of a spatial graph and a simple path graph (time dimension).  By copying the spatial graph through time, one can model time-varying graph signals. This type of modeling has a rich literature in graph signal processing. The reviewer is kindly referred to [4-7], and also to Section 5 in the first chapter of [8].\n\n**Q4**. The main difference between CITRUS and other baselines lies in its ability to learn the graph receptive fields $t$, which is difficult for other baselines (with fixed or non-adaptive receptive fields) to estimate. We've plotted the learned $t$ for different numbers of horizons in the **uploaded PDF** on the main rebuttal (Fig. R3). We observe that the estimation of $t$ is more robust for longer horizons since $t<b=P\\ln{s}/2\\lambda$ (here, $b=7.08$). According to our theoretical analysis, this alleviates the over-smoothing phenomenon for longer horizons.\n\n**Q5**. 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 manuscript. For the unsupervised method, we can analyze the Laplacian matrices for a low-rank approximation task. For example, in the **uploaded PDF** on the main rebuttal (Fig. R2), we plot the explained variances of the principal components of the spatial Laplacian (i.e., $||v\\lambda v^\\top ||_F^2$ for $v$ and $\\lambda$ 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.\n\nReferences:\n-----\n[1] “Tensor decompositions and applications,” SIAM Review, 2009\n\n[2] “Graph neural networks exponentially lose expressive power for node classification,” ICLR, 2020\n\n[3] “Learning product graphs from spectral templates,” IEEE TSIPN, 2023\n\n[4] “Big data analysis with signal processing on graphs”, IEEE SPM, 2014\n\n[5] \"Product graph learning from multi-domain data with sparsity and rank constraints,” IEEE TSP, 2021\n\n[6] \"Learning product graphs from multidomain signals,” IEEE ICASSP 2020\n\n[7] “Product graph Gaussian processes for multi-domain data imputation and active learning,” IEEE EUSIPCO, 2023\n\n[8] \"Vertex-frequency analysis of graph signals\", Springer, 2019\n\n[9] \"Era of big data processing: A new approach via tensor networks and tensor decompositions\". 2014", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "TkMY4WXQQS", "reviewer": "Reviewer_QFxx", "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.", "questions": "1) what are the training times for the CITRUS and its comparison to baseline methods (or the most relevant baseline)?\n\n2) 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? \n\n3) 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? \n\n4) 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": "yes", "rating": 5, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 2, "strengths": "1) improvements on prediction accuracy fro the proposed method compared to existing methods.\n2) Overall well-analyzed methods with some theoretical support on the behavior of the proposed method", "weaknesses": "1) There is no clear strategy to overcome oversmoothing.\n2) 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.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "**Q1**. The training time comparison with the most relevant baseline, GTCNN, has been provided in the **uploaded PDF** on the main rebuttal (Tab. R2) for the NOAA and MetrLA datasets. We observe that our model requires less time per epoch to be trained.\n\n**Q2**. \n- Based on Eqn. (18), we can expect to slow down the over-smoothing phenomenon by keeping $\\tilde{t}<\\frac{P}{2\\tilde{\\lambda}}\\ln{s}$. This explains why we should not increase the receptive field parameter $t$, especially when dealing with strongly connected graphs or deep GNNs. This finding is consistent with previous research on over-smoothing in discrete GNNs [1]. In the **uploaded PDF** on the main rebuttal (Fig. R3), we show the values of the learned receptive fields $t$ by CITRUS across different numbers of horizons. For this experiment, we have $b=7.08$, where $b=\\frac{P}{2\\tilde{\\lambda}}\\ln{s}$. We observe in the figure that CITRUS always learns $t<b$ to slow down the over-smoothing for all horizons.\n- To further exemplify the implications of our theoretical analysis of over-smoothing, we now consider the graph receptive field $t$ as a hyperparameter in CITRUS and design an experiment for different values of $t$. In this experiment, we vary $t\\in\\{0.1,1,5,10,20\\}$ and monitor the convergence to the over-smoothing state. We include the results of this experiment in the **uploaded PDF** on the main rebuttal (Fig. R4). We observe that for values of $t$ higher than $b$, CITRUS converges faster to the over-smoothing state for a larger number of layers. In practice, alleviating over-smoothing might not be enough to achieve good performance. Therefore, for $t<b$, we might need to perform other kinds of analyses like over-squashing, which we leave for future work.\n\n**Q3**. CITRUS has several hyperparameters related to typical architectural and optimization choices in neural networks, such as the number of layers, learning rate, weight decay, etc.  We’ve included all hyperparameters in the Appendix of the revised manuscript. Apart from these typical hyperparameters, the number of selected eigenvalue-eigenvector pairs $k$ is unique to our framework.\n\nWe 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 manuscript (Table 4). For the unsupervised method, we can analyze the Laplacian matrices for a low-rank approximation task. For example, in the **uploaded PDF** on the main rebuttal (Fig. R2), we plot the explained variances of the principal components (i.e., $||v\\lambda v^\\top ||_F^2$ for $v$ and $\\lambda$ being the eigenvectors and eigenvalues) on the MetrLA dataset. We observe a strong concentration of variance in a few components. We also observe that we can capture almost 80% of the explained variance by only relying on 50 eigenvalue-eigenvector pairs. In this case, the recommended range of $k$ is about 10-25% of the number of nodes.\n\nChoosing an appropriate $k$ is indeed related to an efficient rank approximation on the factor Laplacians.\n\n**Q4**. Our framework does not make any assumption about homophily or heterophily properties of the underlying graphs. We have reported the intra-homophily (spatial dimension) and inter-homophily (temporal dimension) of the MetrLA and PemsBay datasets in Tab. R3 in the **uploaded PDF** on the main rebuttal, using the proposed algorithms in [2,3]. We observe in Tab. R3 that these datasets have a mix of low and high measures, indicating a combination of homophily and heterophily behaviors. Therefore, our framework can efficiently handle both homophily and heterophily cases.\n\nReferences:\n-----\n[1] “Graph neural networks exponentially lose expressive power for node classification,” ICLR, 2020\n\n[2] “Greto: Remedying dynamic graph topology-task discordance via target homophily,” ICLR, 2023\n\n[3] “Graph neural networks for graphs with heterophily: A survey”, 2022", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "myA3sbMMPK", "reviewer": "Reviewer_Y9UP", "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.", "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?\n\n2. 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?\n\n[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": "Future directions to improve the theoretical results are discussed.", "rating": 5, "confidence": 4, "soundness": 2, "presentation": 3, "contribution": 2, "strengths": "The authors provide a rather thorough analysis of the proposed methods. Theoretical results are provided. Experiments are detailed. The paper is well-written and easy to follow.", "weaknesses": "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.\n\nAnother 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.\n\nMinor points: figure 2 is not very intuitive; figure 3 is also not very convincing.\n\n[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.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "**We’ve expanded some points of our rebuttal in additional comments. We kindly invite the reviewer to read them if required.**\n\n**W1**. We should have mentioned that, first, we need to define the Dirichlet energy for the tensorial data, $D_T(.)$, as the summation of the factor-wise Dirichlet energies on the unfolded metricized forms of the tensorial data $D_T(x):=\\frac{1}{P}\\sum_{p=1}^{P}{tr(\\underline{X}\\_{(p)}^\\top\\hat{L}\\_p\\underline{X}\\_{(p)})}$, where $x=vec(\\underline{X})$ and $\\underline{X}_{(p)}$ is the $p$-th mode matricization of $\\underline{X}$. The main intuition for this definition of $D_T(.)$ is that multi-domain data might be smooth on one of the factor graphs and non-smooth on the other factor graphs. Therefore, we cannot consider such cases as “pure” smooth signals on the product graph. This intuition aligns with relevant works in this direction [1,9], where these previous studies defined Diricihelet energy on the combinatorial Laplacian of the product graphs as the summation of the factor-wise Diricihelet energies. Note that the over-smoothing analysis in our paper holds in this case too. One can easily show that $D_T(x)$ is equivalent to writing the total Dirichlet energy using the product Laplacian as defined in Lemma 3.8 as $D_T(x)=x^\\top \\hat{L} x$, where $x=vec(\\underline{X})$. \n\nOnce we have this definition for the Dirichlet energy, our definition of the Laplacian matrix in Eqn. (14) emerges naturally. This matrix has two favorable properties: i) it’s positive semidefinite, and ii) its spectrum lies in $[0,2]$ (as in regular normalized Laplacians) having a bounded spectrum for numerical instability. We agree with the reviewer that this Laplacian might not be generally a valid normalized Laplacian for the product graph and deserves further exploration in future work. We’ve corrected some sentences and added a remark regarding this discussion in the Appendix section of the revised manuscript.\n\n**W2**. Let’s consider $h(L)=e^{-L\\_{1}}\\otimes \\sum_{i=1}^{K}{L^{i}\\_2}$ as a separable graph polynomial filter. Since $h(L)=(V_1\\otimes V_2)(e^{-\\Lambda_1}\\otimes\\sum_{i=1}^{K}{\\Lambda^i_2})(V_1\\otimes V_2)^\\top$, the KP structure of the eigenvectors holds but we are not facing a valid well-defined graph product. However, instead of computing direct EVD of $h$ with complexity $(N_1N_2)^3$, we can first compute the factor graphs EVD with complexity $N_1^3+N_2^3$, and then $V_1\\otimes V_2$ and $e^{-\\Lambda_1}\\otimes\\sum_{i=1}^{K}{\\Lambda^i_2}$. Although we agree about the relationships between having a heat diffusion filter defined on a Cartesian Laplacian and KP of separable heat filters, please note that computational benefits come from the possibility of selective eig-eiv pairs. For example, consider a heat diffusion defined on a Kronecker product graph adjacency. In this case, we have $h(A)=e^{-A_1\\otimes A_2}=(V_1\\otimes V_2)(e^{-\\Lambda_1\\otimes\\Lambda_2})(V_1\\otimes V_2)^\\top$, and since we cannot write it as $h_1(A_1)\\otimes h_2(A_2)$, we cannot use a selective strategy. While, in the Cartesian products, we can write $e^{-L_1}\\otimes e^{-L_2}$ and also we have $e^{-L_1}=\\sum_{i=1}^{N_1}{e^{-\\lambda^{(1)}_i} v^{(1)}_1 {v^{(1)}_1}^\\top}$. Therefore, selecting the most important pairs (as a low-rank approximation problem) has a theoretical justification, which is not generally true for Kronecker graphs.\n\n**W3**. These minor points are answered in the comments.\n\n**Q1**.\n- It has been shown that computing the forward propagation in a spectral GNN using the Chebyshev polynomials is $\\mathcal{O}(K|\\mathcal{E}|F\\_{in}F\\_{out})$, with $K$, $|\\mathcal{E}|$, $F\\_{in}$, and $F\\_{out}$ being graph filter order, number of edges of the product graph, input and output features, respectively [4]. For our product graph model, the number of edges can be expressed as $|\\mathcal{E}|=\\sum_{p=1}^{P}{(\\prod_{i=1\\ne p}^{P}{N\\_p})|\\mathcal{E}\\_p|}$ [5]. Therefore, the complexity of our forward propagation in Eqn. (7) is $\\mathcal{O}(K\\_P+K\\_PF\\_l+K\\_PNF\\_l+K\\_PF\\_lF\\_{l+1}+NK\\_PF\\_{l+1})$, where $K_P=\\prod_{p=1}^{P}{K_p}$. Assuming often $K_P\\ll N$, our complexity approximately takes the form of $\\mathcal{O}(K_PN(F_l+F_{l+1}))$. So, we can say that the complexity of our forward function is linear in terms of the number of nodes of the product graph, and can be reduced by selecting the most important eig-eiv pairs in factor graphs $K_P$. Note that, we also need one preprocessing step of the factor graph’s EVD with a complexity $\\mathcal{O}(N^2_pK_p)$, but it's only needed once.\n\n- The PDEs on graphs are defined using Laplacian and not adjacency matrices, so using the adjacency leads to a lack of theory.\n\n**Q2**.\n- The full details about the baselines are provided in Appendix B.\n- The most general method that relies on product graphs, to the best of our knowledge, is GTCNN which is already compared in our paper. However, GTCNN is limited to only two-factor graphs.\n- Please note that we have already compared with the traditional (graph) SP methods GP-VAR [6] and GP-VARMA [6]. Indeed, GP-VAR and GP-VARMA are built upon [7].\n\nReferences:\n-----\n[1] “Low-rank and smooth tensor recovery on Cartesian product graphs” IEEE International Conference on Sampling Theory and Applications, 2023\n\n[2] “Graph neural networks exponentially lose expressive power for node classification” ICLR, 2020\n\n[3] “A note on over-smoothing for graph neural networks” 2020\n\n[4] “Convolutional neural networks on graphs with fast localized spectral filtering” NeurIPS, 2016\n\n[5] “Learning product graphs underlying smooth graph signals” 2020\n\n[6] “Forecasting time series with Varma recursions on graphs” IEEE TSP, 2019\n\n[7] “A time-vertex signal processing framework: Scalable processing and meaningful representations for time-series on graphs” IEEE TSP, 2017\n\n[8] \"The emerging field of signal processing on graphs\" IEEE SPM 2013\n\n[9] “Product graph learning from multi-domain data with sparsity and rank constraints.” IEEE TSP, 2021", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}]

Author rebuttal:
We express our gratitude to the reviewers for their thoughtful and constructive feedback. We are encouraged by their recognition of several strengths in our work. In particular, the reviewers found that: "*The proposed methods are backed by theory*" (reviewers wQ9d, QFxx, and Y9UP), "*The experimental results are SOTA*" (reviewers wQ9d and QFxx), "*The paper is overall well written and easy to follow*" (reviewers wQ9d, gs8b, and Y9UP), and the paper presents a "*thorough analysis of the proposed methods*" (reviewers QFxx and Y9UP). Besides, they found that “*handling multi-domain graph data is an interesting and challenging problem in practice*” (reviewer gs8b).

We have addressed each of the reviewers' comments in a detailed manner, focusing on clarifying any points of uncertainty and resolving any misunderstandings that may have arisen. Our responses are organized in a point-by-point format, as outlined in the subsequent sections of this rebuttal. A summary of these responses is listed as follows:

* We have provided additional analyses on the over-smoothing phenomenon and discussed how to alleviate it.

* We discussed the approximation error that occurs due to the Cartesian product assumption when a graph cannot be fully represented as a Cartesian product of factor graphs.

* We clarified the methodology used for tuning and analyzing the hyperparameter of CITRUS.

* We added additional theoretical and experimental comparisons with the most relevant baselines and spectral GNNs to gain more insights into the advantages of our method.

* We clarified the main intuition behind using Cartesian product graphs to model time-varying spatiotemporal time series.

* We outlined the ability of the proposed method to handle both homophily and heterophily in graphs.

* We highlighted the computational advantages of the separability of the graph filters.

**We have uploaded a one-page PDF file and included additional results as tables and figures. We kindly invite reviewers to refer to them whenever required.**

Meta-review: 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.
Decision: Accept (poster)