# Continuous Product Graph Neural Networks

Aref Einizade, Fragkiskos D. Malliaros, Jhony H. Giraldo

► To cite this version:

Aref Einizade, Fragkiskos D. Malliaros, Jhony H. Giraldo. Continuous Product Graph Neural Networks. NeurIPS 2024 - 38th Conference on Neural Information Processing Systems, Dec 2024, Vancouver (BC), Canada. hal-04896043

HAL Id: hal-04896043

https://hal.science/hal-04896043v1

Submitted on 19 Jan 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>
Continuous Product Graph Neural Networks

Aref Einizade

LTCI, Télécom Paris
Institut Polytechnique de Paris
aref.einizade@telecom-paris.fr

Fragkiskos D. Malliaros

CentraleSupélec, Inria
Université Paris-Saclay
fragkiskos.malliaros@centralesupelec.fr

Jhony H. Giraldo

LTCI, Télécom Paris
Institut Polytechnique de Paris
jhony.giraldo@telecom-paris.fr

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.

# 1 Introduction

Multidomain (tensorial) data defined on multiple interacting graphs [1-3], referred to as multidomain graph data in this paper, extend the traditional graph machine learning paradigm, which typically deals with single graphs [2, 4]. Tensors, which are multi-dimensional generalizations of matrices (order-2 tensors), appear in various fields like hyperspectral image processing [5], video processing [6], recommendation systems [7], spatiotemporal analysis [8], and brain signal processing [9]. Despite the importance of these applications, learning from multidomain graph data has received little attention in the existing literature [2, 10]. Therefore, developing graph-learning strategies for these tensorial data structures holds significant promise for various practical applications.

The main challenge for learning from multidomain graph data is creating efficient frameworks that model joint interactions across domain-specific graphs [10, 11]. Previous work in this area has utilized discrete graph filtering operations in product graphs (PGs) [10, 12] from the field of graph signal processing (GSP) [13]. However, these methods inherit the well-known issues of over-smoothing and over-squashing from regular graph neural networks (GNNs) [14-16], which restricts the graph's receptive field and hinders long-range interactions [17]. Additionally, these methods often require computationally intensive grid searches to tune hyperparameters and are typically limited to two-domain graph data, such as spatial and temporal dimensions [10, 12, 18, 19].
In regular non-tensorial GNNs, continuous GNNs (CGNNs) emerge as a solution to over-smoothing and over-squashing [20]. CGNNs are neural network solutions to partial differential equations (PDEs) on graphs [20, 21], like the heat or wave equation. The solution to these PDEs introduces the exponential graph filter as a continuous infinity-order generalization of the typical discrete graph convolutional filters [20]. Due to the differentiability of exponential graph filters w.r.t. the graph receptive fields, CGNNs can benefit from both global and local message passing by adaptively learning graph neighborhoods, alleviating the limitations of discrete GNNs [17, 20, 22]. Despite these advantages, CGNNs mostly rely on a *single* graph and lack a principled framework for learning *joint* multi-graph interactions.

In this paper, we introduce tensorial PDE on graphs (TPDEGs) to address the limitations of existing PG-based GNN and CGNN frameworks. TPDEGs provide a principled framework for modeling multidomain data residing on multiple interacting graphs that form a Cartesian product graph heat kernel. We then propose **Continuous Product Graph Neural Networks (CITRUS)** as a continuous solution to TPDEGs. We efficiently implement CITRUS using a small subset of eigenvalue decompositions (EVDs) from the factor graphs. Additionally, we conduct theoretical and experimental analyses of the stability and over-smoothing properties of CITRUS. We evaluate our proposed model on spatiotemporal forecasting tasks, demonstrating that CITRUS achieves state-of-the-art performance compared to previous Temporal-Then-Spatial (TTS), Spatial-Then-Temporal (STT), and Temporal-and-Spatial (T&S) frameworks. Our main contributions are summarized as follows:

* We introduce TPDEGs to handle multidomain graph data. This leads to our proposal of a continuous graph filtering operation in Cartesian product spaces, called CITRUS, which naturally emerges as the solution to TPDEGs.

* We conduct extensive theoretical and experimental analyses to evaluate the stability and 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 a well-established field in mathematics and signal processing [23, 24]. However, to the best of our knowledge, the graph-time convolutional neural network (GTCNN) [10, 18] model is the only GNN-based framework addressing neural networks in PGs. GTCNN defines a discrete filtering operation in PGs to jointly process spatial and temporal data for forecasting applications. Due to the discrete nature of the graph polynomial filters in GTCNN, an expensive grid search is required to correctly tune the graph filter orders, leading to significant computational overhead. Additionally, this discrete aspect restricts node-wise receptive fields, making long-range communication challenging [17, 20]. Furthermore, GTCNN is limited to two factor graphs (space and time), limiting its applicability to general multidomain graph data. In contrast, CITRUS overcomes these limitations by i) employing continuous graph filtering, and ii) providing a general framework for any number of factor graphs.

**Applicability.** We explore spatiotemporal analysis as a specific application of PG processing. Early works in GSP for spatiotemporal forecasting, such as the graph polynomial vector auto-regressive moving average (GP-VARMA) and GP-VAR models [25], do not utilize PGs. Modern 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, 26], where the fusion of temporal and spatial modules is crucial. However, similar to STT networks, T&S frameworks also face challenges with computational complexity [8] and their discretized block-wise nature, which can limit adaptive graph receptive fields [20, 26].

CITRUS leverages *joint* multidomain learning by exploiting continuous separable heat kernels as factor-based functions defined on a Cartesian product graph. The continuity and differentiability
<figure><img src="image_2.png" /><figcaption>Figure 1: Illustration of key concepts of CITRUS. a) Cartesian product between three-factor graphs. b) Continuous product graph function (CITRUS) operating on the multidomain graph data $\underline{U}$.</figcaption></figure>

of these filters allow the graph receptive fields to be learned adaptively during the training process, eliminating the need for a grid search. Additionally, CITRUS maintains low computational complexity by relying on the spectral decomposition of the factor graphs [17, 22]. Furthermore, unlike non-graph approaches, the number of learnable parameters in CITRUS is independent of the factor graphs, ensuring scalability.

## 3 Methodology and Theoretical Analysis

**Preliminaries and notation.** An undirected and weighted graph $\mathcal{G}$ with $N$ vertices can be stated as $\mathcal{G} = \{\mathcal{V}, \mathcal{E}, \mathbf{A}\}$. $\mathcal{V}$ and $\mathcal{E}$ denote the sets of the vertices and edges, respectively. $\mathbf{A} \in \mathbb{R}^{N \times N}$ is the adjacency matrix with $\{\mathbf{A}_{ij} = \mathbf{A}_{ji} \ge 0\}_{i,j=1}^N$ representing the connection between vertices, and $\{\mathbf{A}_{ii} = 0\}_{i=1}^N$ states that there are no self-loops in $\mathcal{G}$. The graph Laplacian $\mathbf{L} \in \mathbb{R}^{N \times N}$ is defined as $\mathbf{L} = \mathbf{D} - \mathbf{A}$, where $\mathbf{D} = \text{diag}(\mathbf{A1})$ is the diagonal degree matrix and $\mathbf{1}$ is the all-one vector. A graph signal is a function that maps a set of nodes to the real values $x: \mathcal{V} \to \mathbb{R}$, and therefore we can represent it as $\mathbf{x} = [x_1, \dots, x_N]^\top$. For a vector $\mathbf{a} = [a_1, \dots, a_N]^\top \in \mathbb{R}^{N \times 1}$, $e^{\mathbf{a}} := [e^{a_1}, \dots, e^{a_N}]^\top$. The vectorization operator is shown as $\text{vec}(\cdot)$. Using the Kronecker product $\otimes$, Kronecker sum (aka Cartesian product) $\oplus$, and the Laplacian factors $\{\mathbf{L}_p \in \mathbb{R}^{N_p \times N_p}\}_{p=1}^P$, we have the following definitions:

$$ \downarrow \otimes_{p=1}^{P} \mathbf{L}_p := \mathbf{L}_P \oplus \dots \oplus \mathbf{L}_1, \quad \otimes_{p=1}^{P} \mathbf{L}_p := \mathbf{L}_1 \otimes \dots \otimes \mathbf{L}_P, \quad \downarrow \otimes_{p=1}^{P} \mathbf{L}_p := \mathbf{L}_P \otimes \dots \otimes \mathbf{L}_1, \quad (1) $$

where the Laplacian matrix of the Cartesian product between two factor graphs can be stated as:

$$ \oplus_{p=1}^{2} \mathbf{L}_p = \mathbf{L}_1 \oplus \mathbf{L}_2 = \mathbf{L}_1 \otimes \mathbf{I}_{N_2} + \mathbf{I}_{N_1} \otimes \mathbf{L}_2, \quad (2) $$

with $\mathbf{I}_n$ the identity matrix of size $n$. See Figure 1a for an example of a Cartesian graph product between three factor graphs. A similar relationship to (2) also holds for adjacency matrices [24, 27]. We define a $D$-dimensional tensor as $\underline{\mathbf{U}} \in \mathbb{R}^{N_1 \times \dots \times N_D}$. The matrix $\underline{\mathbf{U}}_{(i)} \in \mathbb{R}^{N_i \times [\prod_{r=1, r \neq i}^{D} N_r]}$ is the $i$-th mode matricization of $\underline{\mathbf{U}}$ [28]. The mode-$i$ tensorial multiplication of a tensor $\underline{\mathbf{U}}$ by a matrix $\mathbf{X} \in \mathbb{R}^{m \times N_i}$ is denoted by $\underline{\mathbf{G}} = \underline{\mathbf{U}} \times_i \mathbf{X}$, where $\underline{\mathbf{G}} \in \mathbb{R}^{N_1 \times \dots \times N_{i-1} \times m \times N_{i+1} \times \dots \times N_D}$ [28]. For further information about the tensor matricization and mode-$n$ product, please refer to [28, 29] (especially Sections 2.4 and 2.5 in [28]). All the proofs of theorems, propositions, and lemmas of this paper are provided in the Appendix A.

### 3.1 Continuous Product Graph Neural Networks (CITRUS)

We state the generative PDE for CITRUS as follows:

**Definition 3.1** (Tensorial PDE on graphs (TPDEG)). Let $\tilde{\mathbf{U}}_t \in \mathbb{R}^{N_1 \times N_2 \times \dots \times N_P}$ be a multidomain tensor whose elements are functions of time $t$, and let $\{\mathbf{L}_p \in \mathbb{R}^{N_p \times N_p}\}_{p=1}^P$ be the domain-specific factor Laplacians. We can define the TPDEG as follows:

$$ \frac{\partial \tilde{\mathbf{U}}_t}{\partial t} = - \sum_{p=1}^{P} \tilde{\mathbf{U}}_t \times_p \mathbf{L}_p. \quad (3) $$
**Theorem 3.2.** Let $\tilde{\underline{U}}_0$ be the initial value of $\tilde{\underline{U}}_t$. The solution to the TPDEG in (3) is given by:

$$ \tilde{\underline{U}}_t = \tilde{\underline{U}}_0 \times_1 e^{-t\underline{L}_1} \times_2 e^{-t\underline{L}_2} \times_3 \dots \times_P e^{-t\underline{L}_P}. \quad (4) $$

Using Theorem 3.2, we define the core function of CITRUS as follows:

$$ f(\underline{U}_l) = \underline{U}_l \times_1 e^{-t_l \underline{L}_1} \times_2 \dots \times_P e^{-t_l \underline{L}_P} \times_{P+1} \underline{W}_l^\top, \quad (5) $$

where $\underline{U}_l \in \mathbb{R}^{N_1 \times \dots \times N_P \times F_l}$ is the input tensor, $\underline{W}_l \in \mathbb{R}^{F_l \times F_{l+1}}$ is a matrix of learnable parameters, $t_l$ is the learnable graph receptive field, and $f(\underline{U}_l) \in \mathbb{R}^{N_1 \times \dots \times N_P \times F_{l+1}}$ is the output tensor. Figure 1b illustrates a high-level description of $f(\underline{U}_l)$. The next proposition formulates (5) as a graph convolution defined on a product graph.

**Proposition 3.3.** The core function of CITRUS in (5) can be rewritten as:

$$ [f(\underline{U}_l)_{(P+1)}]^\top = e^{-t_l \underline{L}_\circ} [\underline{U}_{l(P+1)}]^\top \underline{W}_l, \quad (6) $$

where $\underline{L}_\circ := \bigoplus_{p=1}^P \underline{L}_p$ is the Laplacian of the Cartesian product graph.

Proposition 3.3 is the main building block for implementing CITRUS. More precisely, we use the spectral decompositions of the factor graphs $\{\underline{L}_p = \underline{V}_p \underline{\Lambda}_p \underline{V}_p^\top\}_{p=1}^P$ and product graph $\{\underline{L}_\circ = \underline{V}_\circ \underline{\Lambda}_\circ \underline{V}_\circ^\top\}$, where $\underline{V}_\circ = \bigotimes_{p=1}^P \underline{V}_p$ and $\underline{\Lambda}_\circ = \bigoplus_{p=1}^P \underline{\Lambda}_p$ [23]. Let $K_p \le N_p$ be the number of selected eigenvalue-eigenvector pairs of the p-th factor Laplacian. When $K_p = N_p$, it can be shown [17, 22] that we can rewrite (6) as follows:

$$ [f(\underline{U}_l)_{(P+1)}]^\top = \underline{V}_\circ^{(K_p)} \left( \underbrace{[\tilde{\underline{\lambda}}_l | \dots | \tilde{\underline{\lambda}}_l]}_{\tilde{\underline{\lambda}}_l} \odot (\underline{V}_\circ^{(K_p)})^\top [\underline{U}_{l(P+1)}]^\top \right) \underline{W}_l, \quad (7) $$

$$ \text{with} \quad \tilde{\underline{\lambda}}_l = \bigotimes_{p=1}^P e^{-t_l \underline{\lambda}_p^{(K_p)}}, \quad \underline{V}_\circ^{(K_p)} = \bigotimes_{p=1}^P \underline{V}_p^{(K_p)}, \quad (8) $$

where $\underline{\lambda}_p^{(K_p)} \in \mathbb{R}^{K_p \times 1}$ and $\underline{V}_p^{(K_p)} \in \mathbb{R}^{N_p \times K_p}$ are the first $K_p \le N_p$ selected eigenvalues and eigenvectors of $\underline{L}_p$ based on largest eigenvalue magnitudes, respectively, and $\odot$ is the point-wise multiplication operation. Finally, we can define the output of the l-th layer of CITRUS as $\underline{U}_{l+1(P+1)} = \sigma(f(\underline{U}_l)_{(P+1)})$, where $\sigma(\cdot)$ is some proper activation function.

**Remark 3.4.** We can consider factor-specific learnable graph receptive fields $\{t_l^{(p)}\}_{p=1}^P$ for formulating $\tilde{\underline{\lambda}}_l$ in (8) as $\tilde{\underline{\lambda}}_l = \bigotimes_{p=1}^P e^{-t_l^{(p)} \underline{\lambda}_p^{(K_p)}}$ to make CITRUS more flexible for each factor graph's receptive field. This technique can be expanded to the channel-wise graph receptive fields $\{t_l^{(p,c)}\}_{p=1,c=1}^{P,F_l}$ such that the c-th column of $\tilde{\underline{\lambda}}_l$ in (7) is obtained by $\tilde{\underline{\lambda}}_{k,c} = \bigotimes_{p=1}^P e^{-t_l^{(p,c)} \underline{\lambda}_p^{(K_p)}}$.

**Remark 3.5.** Computing the EVD in the product graph $\underline{L}_\circ$ in (6) has a general complexity of $O(\left[\prod_{p=1}^P N_p\right]^3)$. However, we obtain a significant reduction in complexity of $O(\left[\sum_{p=1}^P N_p^3\right])$ in (7) by relying upon the properties of product graphs since we perform EVD on each factor graph independently. Besides, if we rely only on the $K_p$ most important eigenvector-eigenvalue pairs [17] of each factor graph, we can reduce the complexity of the spectral decomposition to $O(N_p^2 K_p)$, obtaining a general complexity of $O(\left[\sum_{p=1}^P N_p^2 K_p\right])$.

## 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 to the p-th graph as an addition of an error matrix $\underline{E}_p$ (with upper bound $\varepsilon_p$ for any matrix norm $\|\cdot\|$) to the adjacency matrix $\underline{A}_p$ as:

$$ \tilde{\underline{A}}_p = \underline{A}_p + \underline{E}_p; \quad \|\underline{E}_p\| \le \varepsilon_p. \quad (9) $$

**Proposition 3.6.** Let $\underline{A}$ and $\tilde{\underline{A}}$ be the true and perturbed adjacency matrix of Cartesian product graphs with the perturbation model for each factor graph described as in (9). Then, it holds that:

$$ \tilde{\underline{A}} = \underline{A} + \underline{E}; \quad \|\underline{E}\| \le \sum_{p=1}^{P} \varepsilon_p, \quad (10) $$

where the perturbation matrix $\underline{E}$ also follows the Cartesian structure $\underline{E} = \bigoplus_{p=1}^P \underline{E}_p$.
Finally, let $\varphi(u, t)$ and $\tilde{\varphi}(u, t)$ be the true and perturbed output of the CITRUS model with normalized Laplacian $\hat{\mathbf{L}}_{\diamond}$ being responsible for generating $\mathbf{U}_{t(P+1)}$ in a heat flow PDE with normalized Laplacian $\hat{\mathbf{L}}_{\diamond}$ as $\frac{\partial\varphi(u,t)}{\partial t} = -\hat{\mathbf{L}}_{\diamond}\varphi(u,t)$ in (6) [30], respectively. Then, the following theorem states that the integrated error bound on the stability properties of CITRUS is also separable when dealing with Cartesian product graphs:

**Theorem 3.7.** Consider a Cartesian product graph (without isolated nodes) with the normalized true and perturbed Laplacians $\hat{\mathbf{L}}_{\diamond}$ and $\tilde{\mathbf{L}}_{\diamond}$ in (6). Also, assume we have the perturbation model $\hat{\mathbf{A}}_p = \mathbf{A}_p + \mathbf{E}_p$ with factor error bounds $\{\|\mathbf{E}_p\| \le \varepsilon_p\}_{p=1}^P$ in Proposition 3.6. Then, the stability bound on the true and perturbed outputs of CITRUS, i.e., $\varphi(u, t)$ and $\tilde{\varphi}(u, t)$, respectively, can be described by the summation of the factor stability bounds as:

$$ \|\varphi(u, t) - \tilde{\varphi}(u, t)\| = \sum_{p=1}^{P} \mathcal{O}(\varepsilon_p). \qquad (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 erroneous 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 [20, 31, 32] and can be formulated [20] for the output of a continuous GNN. Precisely, it can be defined on the GNN's output $\mathbf{U}_t = [\mathbf{u}(t)_1, \dots, \mathbf{u}(t)_N]^\top$, with $\mathbf{u}(t)_i$ being the feature vector of the i-th node with the degree $\deg_i$, as:

$$ \lim_{t \to \infty} E(\mathbf{U}_t) \to 0, \qquad (12) $$

$$ \text{where } E(\mathbf{U}) := \frac{1}{2} \sum_{(i,j) \in \mathcal{E}} \left\| \frac{\mathbf{u}_i}{\sqrt{\deg_{i}}} - \frac{\mathbf{u}_j}{\sqrt{\deg_{j}}} \right\|^2 = \text{tr}(\mathbf{U}^{\top} \hat{\mathbf{L}} \mathbf{U}). \qquad (13) $$

Inspired by (13), we extend the definition of Dirichlet energy to the tensorial case as follows:

**Definition 3.8.** (Tensorial Dirichlet energy). By considering the normalized factor Laplacians $\{\hat{\mathbf{L}}_p\}_{p=1}^P$, we define the Tensorial Dirichlet energy for a tensor $\mathbf{U} \in \mathbb{R}^{N_1 \times \dots \times N_P \times F}$ as:

$$ E(\mathbf{U}) := \frac{1}{P} \sum_{f=1}^{F} \sum_{p=1}^{P} \text{tr}(\tilde{\mathbf{U}}_{f(p)}^{\top} \hat{\mathbf{L}}_p \tilde{\mathbf{U}}_{f(p)}), \qquad (14) $$

where $\tilde{\mathbf{U}}_{f(p)} = \mathbf{U}[:, \dots, :; f]_{(p)} \in \mathbb{R}^{N_p \times \prod_{i \neq p}^{P} N_i}$, i.e., the p-th mode matricization of the f-th slice of $\mathbf{U}$ in its $(P+1)$-th dimension.

Based on the previous research in the GSP literature [27, 33, 34], it follows that $E(\mathbf{U})$ in (14) can be rewritten as $E(\mathbf{U}) = \text{tr}(\tilde{\mathbf{U}}^{\top} \hat{\mathbf{L}} \tilde{\mathbf{U}})$, where $\tilde{\mathbf{U}} \in \mathbb{R}^{[\prod_{p=1}^{P} N_p] \times F}$, $\tilde{\mathbf{U}}[:, \dots, :; f] = \text{vec}(\mathbf{U}[:, \dots, :; f])$, and $\hat{\mathbf{L}} := \frac{1}{P} \bigoplus_{p=1}^{P} \hat{\mathbf{L}}_p = \bigoplus_{p=1}^{P} \left(\frac{\hat{\mathbf{L}}_p}{P}\right)$. The next lemma shows interesting and applicable properties of $\hat{\mathbf{L}}$.

**Lemma 3.9.** Consider P factor graphs with normalized adjacencies $\{\hat{\mathbf{A}}_p\}_{p=1}^P$ and Laplacians $\{\hat{\mathbf{L}}_p\}_{p=1}^P$. By constructing the product adjacency as $\hat{\mathbf{A}} := \frac{1}{P} \bigoplus_{p=1}^{P} \hat{\mathbf{A}}_p = \bigoplus_{p=1}^{P} \left(\frac{\hat{\mathbf{A}}_p}{P}\right)$, a similar Cartesian form for the product Laplacian $\hat{\mathbf{L}} := \mathbf{I} - \hat{\mathbf{A}}$ (with $\mathbf{I} \in \mathbb{R}^{[\prod_{p=1}^{P} N_p] \times [\prod_{p=1}^{P} N_p]}$ being the identity matrix) also holds as the following:

$$ \hat{\mathbf{L}} = \frac{1}{P} \bigoplus_{p=1}^{P} \hat{\mathbf{L}}_p = \bigoplus_{p=1}^{P} \left(\frac{\hat{\mathbf{L}}_p}{P}\right). \qquad (15) $$

Besides, the spectrum of $\hat{\mathbf{L}}$ is spanned across the interval of $[0, 2]$.
We observe in (14) that the kernel separability of Cartesian heat kernels also leads to the separability of analyzing over-smoothing in different modes of the data at hand. This is useful in cases of facing factor graphs with different characteristics. Next, we formally analyze over-smoothing in CITRUS. In our case, for the l-th layer with $H_l$ hidden MLP layers, non-linear activations $\sigma(\cdot)$ and learnable weight matrices $\{W_{l,h}\}_{h=1}^{H_l}$, we define:

$$X_{l+1} := f_l(X_l), \quad f_l(\mathbf{X}) := \mathrm{MLP}_l(e^{-\hat{\mathbf{L}}^{(t)}\mathbf{X}}), \quad \mathrm{MLP}_l(\mathbf{X}) := \sigma(\dots\sigma(\sigma(\mathbf{X})\mathbf{W}_{l1})\mathbf{W}_{l2}\dots\mathbf{W}_{lH_l}), \quad (16)$$

with $X_0$ being the initial node feature matrix, $f_l(\cdot)$ is a generalized form of (6), and $e^{-\hat{\mathbf{L}}^{(t)}} := \downarrow \otimes_{p=1}^P e^{-t^{(p)}\frac{\hat{\mathbf{L}}_p}{P}}$, where $t^{(i)}$ is the receptive field of the i-th factor graph. Then, the next theorem describes the over-smoothing criteria.

**Theorem 3.10.** For the product Laplacian $\hat{\mathbf{L}} := \frac{1}{P} \oplus_{p=1}^P \hat{\mathbf{L}}_p$ with domain-specific receptive field $t^{(i)}$ for the i-th factor graph and the activations $\sigma(\cdot)$ in (16) being ReLU or Leaky ReLU, we have:

$$E(\mathbf{X}_l) \le e^{l(\ln s - \frac{2t\tilde{\lambda}}{P})} E(\mathbf{X}_0), \quad (17)$$

where $s := \sup_{l \in \mathbb{N}_+} s_l$ and $s_l := \prod_{h=1}^{H_l} s_{lh}$ with $s_{lh}$ being the square of maximum singular value of $\mathbf{W}_{lh}^\top$. Besides, $\tilde{\lambda} = \lambda^{(m)}$ and $\tilde{t} = t^{(m)}$ with $m = \arg\min_i t^{(i)}\lambda^{(i)}$, where $\lambda^{(i)}$ is the smallest non-zero eigenvalue of $\hat{\mathbf{L}}_i$.

**Corollary 3.11.** When $l \to \infty$, $E(\mathbf{X}_l)$ exponentially converges to 0, when:

$$\ln s - \frac{2}{P} \tilde{t} \tilde{\lambda} < 0. \quad (18)$$

Theorem 3.10 shows that the factor graph with the smallest non-zero eigenvalue (spectral gap) multiplied by its receptive field dominates the overall over-smoothing. The less the spectral gap, the less probability the factor graph is connected [14]. So, we can focus on the (much smaller) factor graphs instead of the resulting massive product graph. These theoretical findings are experimentally validated in Section 4.2.

# 4 Experimental Results

In this section, we experimentally validate the theoretical discussions regarding the stability and over-smoothing analysis and use CITRUS in real-world spatiotemporal forecasting tasks. Similarly, we conduct ablation studies regarding CITRUS in spatiotemporal forecasting. We compare CITRUS against several state-of-the-art methods in forecasting. The results on the case of more than two factor graphs and descriptions of the state-of-the-art are provided in Appendix B and C, respectively, alongside the additional theoretical and experimental discussions in Sections D-I. The implementation codes are available at https://github.com/ArefEinizade2/CITRUS.

## 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_1 = 20$, $N_2 = 30$, and edge probabilities $p_{ER}^{(1)} = p_{ER}^{(2)} = 0.1$ for node regression. We utilize 15% of the nodes in the product graph as the test set, 15% of the remaining nodes as validation, and the rest of the nodes for training. We use a three-layer generative model to generate the outputs based on (6). The continuous parameters for the factor graphs are set as $t^{(1)} = 2$ and $t^{(2)} = 3$. The initial product graph signal $\mathbf{X}^{(0)} \in \mathbb{R}^{N \times F_0}$ is generated from the normal distribution with $F_0 = 6$. Each layer's MLP of the generation model has only one layer with the number of hidden units $\{5, 4, 2\}$. We add noise to the factor adjacency matrices with signal-to-noise ratios

<figure><img src="image_3.png" /><figcaption>Figure 2: Stability analysis vs. different SNR scenarios, related the results in Theorem 3.7.</figcaption></figure>
(SNRs) (in terms of db) of {∞, 20, 10, 0, −10}. Therefore, we construct a two-layer CITRUS and use one-layer MLP for each with $F_1 = F_2 = 4$ to learn the representation in different SNRs to study the effect of each factor graph's noise on the stability performance of the network.

Figure 2 shows the averaged prediction error in MSE between the true and predicted outputs over 10 random realizations. Firstly, we notice that increasing the SNR for each factor graph improves performance, illustrating the stability properties of CITRUS. Therefore, the effect of each factor graph's stability on the overall stability is well depicted in this figure, confirming the theoretical findings in Theorem 3.7. For instance, in SNR=10 for the first graph, the performance is still severely affected by the stability for the second graph, shown by color in Figure 2. Finally, we observe that the standard deviation of the predictions is smaller for larger values of SNR, illustrating CITRUS's robustness.

## 4.2 Experimental Over-smoothing Analysis

For experimentally validating Theorem 3.10, we first generate a initial graph signal $\mathbf{X}_0 \in \mathbb{R}^{N \times F_0}$ with $F_0 = 12$ on a 150-node Cartesian product graph with $N_1 = 10, N_2 = 15$. Therefore, we consider a 10-layer generation process based on (6) with only one-layer MLP for each layer and $\{F_l = 12 - l\}_{l=1}^{10}$, where the weight matrices are generated from normal distribution. We select ReLU as the activation function and $\{t_l = 1\}_{l=1}^{10}$.

The over-smoothing phenomenon strongly depends on the smallest spectral gap of the factor graphs as stated in Theorem 3.10. We consider two scenarios of $\ln s - \frac{2}{P}\tilde{t}\tilde{\lambda} < 0$ and $\ln s - \frac{2}{P}\tilde{t}\tilde{\lambda} > 0$. The first scenario relates to a well-connected product graph obtained from two factor ER graphs with $p_{ER}^{(1)} = 0.05$ and $p_{ER}^{(2)} = 0.95$ and, while in the other, $p_{ER}^{(1)} = p_{ER}^{(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 $\lambda = 0.06$ and $s = 0.004 \pm 10^{-5}$ for the first scenario; and $\lambda = 0.07$ and $s = 7.691 \pm 10^{-6}$ for the other.

We depict the outputs $y(l) = \log \frac{E(\mathbf{X}_l)}{E(\mathbf{X}_0)}$ and the theoretical upper bound $y(l) = l(\ln s_l - \frac{2}{P}\tilde{t}_l\tilde{\lambda})$ across the number layers $l$ in Figure 3. We observe that in the **left** plot where $\{\ln s_l - \frac{2}{P}\tilde{t}_l\tilde{\lambda} < 0\}_{l=1}^{10}$, the theoretical upper bound approximates well the actual outputs during over-smoothing as stated in Theorem 3.10. On the opposite, on the **right** plot, in which $\{\ln s_l - \frac{2}{P}\tilde{t}_l\tilde{\lambda} > 0\}_{l=1}^{10}$, the theoretical upper bound is loose. We leave for future work finding a tighter upper bound.

<figure><img src="image_4.png" /><figcaption>Figure 3: Over-smoothing analysis using Theorem 3.10. **left:** $\ln s - \frac{2}{P}\tilde{t}\tilde{\lambda} < 0$, **right:** $\ln s - \frac{2}{P}\tilde{t}\tilde{\lambda} > 0$.</figcaption></figure>

## 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 training-validation-testing data partitions. All hyperparameters were optimized on the validation set with the details in Section I. We use Adam as optimizer.

**Traffic forecasting.** Here, we work on two well-known datasets for spatiotemporal forecasting: MetrLA [35] and PemsBay [36]. MetrLA consists of recorded traffic data for four months on 207 highways in Los Angeles County with 5-minute resolutions [35]. PemsBay contains 5-minute traffic load data for six months associated with 325 stations in the Bay Area by the same setting with [36]. The spatial and temporal graphs are created by applying Gaussian kernels on the node distances [37] and simple path graphs, respectively. The task for these datasets is, by having the last 30 minutes ($T = 6$ steps) of traffic recordings, we need to predict the following 15-30-60 minutes, i.e., $H = 3, H = 6$, and $H = 12$ steps ahead horizons. Our model first uses a linear encoder network, followed by three CITRUS blocks with 3-layer MLPs for each, where $\{F_l^{\text{MLP}} = 64\}_{l=1}^3$. Therefore, the output node embeddings are concatenated with the initial spatiotemporal data. Finally, we use a linear decoding layer to transform the learned representations to the number of needed horizons. We exploit residual connections within each CITRUS block to stabilize the training process [17]. We use the
<table><caption>Table 1: Traffic forecasting comparison between CITRUS and previous methods.</caption><thead><tr><th rowspan="3">Method</th><th colspan="9">MetrLA</th></tr><tr><th colspan="3">H = 3</th><th colspan="3">H = 6</th><th colspan="3">H = 12</th></tr><tr><th>MAE</th><th>MAPE</th><th>RMSE</th><th>MAE</th><th>MAPE</th><th>RMSE</th><th>MAE</th><th>MAPE</th><th>RMSE</th></tr></thead><tbody><tr><td>ARIMA [35]</td><td>3.99</td><td>9.60%</td><td>8.21</td><td>5.15</td><td>12.70%</td><td>10.45</td><td>6.90</td><td>17.40%</td><td>13.23</td></tr><tr><td>G-VARMA [25]</td><td>3.60</td><td>9.62%</td><td>6.89</td><td>4.05</td><td>11.22%</td><td>7.84</td><td>5.12</td><td>14.00%</td><td>9.58</td></tr><tr><td>GP-VAR [25]</td><td>3.56</td><td>9.55%</td><td>6.54</td><td>3.98</td><td>11.02%</td><td>7.56</td><td>4.87</td><td>13.34%</td><td>9.19</td></tr><tr><td>FC-LSTM [35]</td><td>3.44</td><td>9.60%</td><td>6.30</td><td>3.77</td><td>10.90%</td><td>7.23</td><td>4.37</td><td>13.20%</td><td>8.69</td></tr><tr><td>Graph Wavenet [38]</td><td><u>2.69</u></td><td><u>6.90%</u></td><td><u>5.15</u></td><td><u>3.07</u></td><td><u>8.37%</u></td><td><u>6.22</u></td><td><u>3.53</u></td><td>10.01%</td><td><u>7.37</u></td></tr><tr><td>GMAN [39]</td><td>2.94</td><td>7.51%</td><td>5.89</td><td>3.22</td><td>8.92%</td><td>6.61</td><td>3.68</td><td>10.25%</td><td>7.49</td></tr><tr><td>STGCN [40]</td><td>2.88</td><td>7.62%</td><td>5.74</td><td>3.47</td><td>9.57%</td><td>7.24</td><td>4.59</td><td>12.70%</td><td>9.40</td></tr><tr><td>GGRNN [41]</td><td>2.73</td><td>7.12%</td><td>5.44</td><td>3.31</td><td>8.97%</td><td>6.63</td><td>3.88</td><td>10.59%</td><td>8.14</td></tr><tr><td>GRUGCN [26]</td><td><u>2.69</u></td><td><u>6.61%</u></td><td><u>5.15</u></td><td><u>3.05</u></td><td><u>7.96%</u></td><td><u>6.04</u></td><td><u>3.62</u></td><td>9.92%</td><td><u>7.33</u></td></tr><tr><td>SGP [42]</td><td>3.06</td><td>7.31%</td><td>5.49</td><td>3.43</td><td>8.54%</td><td>6.47</td><td>4.03</td><td>10.53%</td><td>7.81</td></tr><tr><td>GTCNN [10, 18]</td><td><u>2.68</u></td><td><u>6.85%</u></td><td><u>5.17</u></td><td><u>3.02</u></td><td><u>8.30%</u></td><td><u>6.20</u></td><td><u>3.55</u></td><td>10.21%</td><td><u>7.35</u></td></tr><tr><td><strong>CITRUS (Ours)</strong></td><td><strong>2.70</strong></td><td><strong>6.74%</strong></td><td><strong>5.14</strong></td><td><strong>2.98</strong></td><td><strong>7.78%</strong></td><td><strong>5.90</strong></td><td><strong>3.44</strong></td><td><strong>9.28%</strong></td><td><strong>6.85</strong></td></tr></tbody></table><table><thead><tr><th rowspan="3">Method</th><th colspan="9">PemsBay</th></tr><tr><th colspan="3">H = 3</th><th colspan="3">H = 6</th><th colspan="3">H = 12</th></tr><tr><th>MAE</th><th>MAPE</th><th>RMSE</th><th>MAE</th><th>MAPE</th><th>RMSE</th><th>MAE</th><th>MAPE</th><th>RMSE</th></tr></thead><tbody><tr><td>ARIMA [35]</td><td>1.62</td><td>3.50%</td><td>3.30</td><td>2.33</td><td>5.40%</td><td>4.76</td><td>3.38</td><td>8.30%</td><td>6.50</td></tr><tr><td>G-VARMA [25]</td><td>1.88</td><td>4.28%</td><td>3.96</td><td>2.45</td><td>5.42%</td><td>4.70</td><td>3.01</td><td>7.10%</td><td>5.83</td></tr><tr><td>GP-VAR [25]</td><td>1.74</td><td>3.45%</td><td>3.22</td><td>2.16</td><td>5.15%</td><td>4.41</td><td>2.48</td><td>6.18%</td><td>5.04</td></tr><tr><td>FC-LSTM [35]</td><td>2.05</td><td>4.80%</td><td>4.19</td><td>2.20</td><td>5.20%</td><td>4.55</td><td>2.37</td><td>5.70%</td><td>4.74</td></tr><tr><td>Graph Wavenet [38]</td><td><u>1.30</u></td><td><u>2.73%</u></td><td><u>2.74</u></td><td><u>1.63</u></td><td><u>3.67%</u></td><td><u>3.70</u></td><td><u>1.95</u></td><td>4.63%</td><td><u>4.52</u></td></tr><tr><td>GMAN [39]</td><td>1.34</td><td>2.81%</td><td>2.82</td><td>1.62</td><td>3.63%</td><td>3.72</td><td>1.86</td><td>4.31%</td><td><u>4.32</u></td></tr><tr><td>STGCN [40]</td><td>1.36</td><td>2.90%</td><td>2.96</td><td>1.81</td><td>4.17%</td><td>4.27</td><td>2.49</td><td>5.79%</td><td>5.69</td></tr><tr><td>GGRNN [41]</td><td>1.33</td><td>2.83%</td><td>2.81</td><td>1.68</td><td><u>3.79%</u></td><td><u>3.94</u></td><td>2.34</td><td>5.21%</td><td>5.14</td></tr><tr><td>GRUGCN [26]</td><td><u>1.21</u></td><td><u>2.49%</u></td><td><u>2.52</u></td><td><u>1.54</u></td><td><u>3.32%</u></td><td><u>3.46</u></td><td>2.01</td><td>4.72%</td><td>4.65</td></tr><tr><td>SGP [42]</td><td>1.33</td><td>2.71%</td><td>2.84</td><td>1.70</td><td>3.61%</td><td>3.83</td><td>2.24</td><td>5.08%</td><td>5.19</td></tr><tr><td>GTCNN [10, 18]</td><td>1.25</td><td>2.61%</td><td>2.66</td><td><u>1.65</u></td><td>3.82%</td><td>3.68</td><td><u>2.27</u></td><td>5.11%</td><td><u>4.99</u></td></tr><tr><td><strong>CITRUS (Ours)</strong></td><td><strong>1.21</strong></td><td><strong>2.51%</strong></td><td><strong>2.61</strong></td><td><strong>1.48</strong></td><td><strong>3.23%</strong></td><td><strong>3.28</strong></td><td><strong>1.78</strong></td><td><strong>4.08%</strong></td><td><strong>3.99</strong></td></tr></tbody></table>

<table><caption>Table 2: Weather forecasting comparison (by rNMSE) between CITRUS and previous methods.</caption><thead><tr><th rowspan="2">Method</th><th colspan="6">Molene</th><th colspan="6">NOAA</th></tr><tr><th>H = 1</th><th>H = 2</th><th>H = 3</th><th>H = 4</th><th>H = 5</th><th>H = 1</th><th>H = 2</th><th>H = 3</th><th>H = 4</th><th>H = 5</th></tr></thead><tbody><tr><td>GRUGCN [26]</td><td>0.49</td><td>0.56</td><td>0.63</td><td>0.70</td><td>0.75</td><td>0.15</td><td><u>0.18</u></td><td><u>0.24</u></td><td><u>0.30</u></td><td><u>0.37</u></td></tr><tr><td>GGRNN [41]</td><td>0.29</td><td>0.42</td><td>0.54</td><td>0.65</td><td>0.75</td><td>0.19</td><td>0.20</td><td>0.27</td><td>0.38</td><td>0.48</td></tr><tr><td>SGP [42]</td><td><u>0.24</u></td><td><u>0.37</u></td><td><u>0.49</u></td><td><u>0.59</u></td><td>0.69</td><td>0.39</td><td>0.41</td><td>0.43</td><td>0.46</td><td>0.49</td></tr><tr><td>GTCNN [10, 18]</td><td>0.39</td><td>0.45</td><td>0.52</td><td>0.60</td><td><u>0.68</u></td><td><u>0.17</u></td><td>0.19</td><td>0.25</td><td>0.31</td><td><u>0.37</u></td></tr><tr><td><strong>CITRUS (Ours)</strong></td><td><strong>0.23</strong></td><td><strong>0.35</strong></td><td><strong>0.47</strong></td><td><strong>0.58</strong></td><td><strong>0.67</strong></td><td><strong>0.04</strong></td><td><strong>0.12</strong></td><td><strong>0.20</strong></td><td><strong>0.29</strong></td><td><strong>0.36</strong></td></tr></tbody></table>

mean absolute error (MAE) as the loss function and consider a maximum of 300 epochs. The best model on the validation set is applied to the unseen test set. We compare CITRUS against previous methods using MAE, mean absolute percentage error (MAPE), and root mean squared error (RMSE).

Table 1 presents the forecasting results on the MetrLA and PemsBay datasets. We observe that, given the abundant training data in these datasets, the NN-based methods outperform the classic GSP-based methods like ARIMA [35], G-VARMA [25], and GP-VAR [25]. Similarly, the GNN-based models are superior compared to non-graph algorithms like FC-LSTM. More importantly, CITRUS outperforms state-of-the-art baselines in most metrics, especially in larger numbers of horizons ($H > 3$).

**Weather forecasting.** We also test CITRUS for weather forecasting in the Molene [43] and NOAA [44] datasets. The Molene dataset provides recorded temperature measurements for 744 hourly intervals over 32 measurement stations in a region in France. For the NOAA, which is associated with regions in the U.S., the temperature measurements were recorded for 8,579 hourly intervals over 109 stations. The pre-processing and data curation settings were similar to prior works on these datasets [10, 25]. The task here is by having the last 10 hours of measurements, we should predict the next {1 : 5} hours. For the Molene dataset, the embedding dimension and the initial linear layer of
<table><caption>Table 3: Ablation study on comparison between the proposed CITRUS and typically ST pipelines.</caption><thead><tr><th rowspan="3">Method</th><th colspan="9">MetrLA</th></tr><tr><th colspan="3">H = 3</th><th colspan="3">H = 6</th><th colspan="3">H = 12</th></tr><tr><th>MAE</th><th>MAPE</th><th>RMSE</th><th>MAE</th><th>MAPE</th><th>RMSE</th><th>MAE</th><th>MAPE</th><th>RMSE</th></tr></thead><tbody><tr><td>TTS</td><td>2.73</td><td>6.78%</td><td>5.21</td><td>3.10</td><td>8.02%</td><td>6.15</td><td>3.67</td><td>10.05%</td><td>7.46</td></tr><tr><td>STT</td><td><u>2.72</u></td><td>6.74%</td><td><u>5.19</u></td><td>3.07</td><td>7.94%</td><td><u>6.08</u></td><td>3.65</td><td>10.97%</td><td><u>7.37</u></td></tr><tr><td>CTTS</td><td><b><u>2.70</u></b></td><td><u>6.69%</u></td><td><u>5.20</u></td><td><u>3.05</u></td><td>7.93%</td><td>6.18</td><td><u>3.61</u></td><td><u>9.75%</u></td><td>7.51</td></tr><tr><td>CSTT</td><td><b><u>2.70</u></b></td><td><b><u>6.66%</u></b></td><td>5.22</td><td>3.06</td><td>7.92%</td><td>6.19</td><td>3.63</td><td>9.92%</td><td>7.48</td></tr><tr><td><b>CITRUS (Ours)</b></td><td><b><u>2.70</u></b></td><td>6.74%</td><td><b><u>5.14</u></b></td><td><b><u>2.98</u></b></td><td><b>7.78%</b></td><td><b><u>5.90</u></b></td><td><b><u>3.44</u></b></td><td><b>9.28%</b></td><td><b><u>6.85</u></b></td></tr></tbody></table>

<table><caption>Table 4: Training time (per epoch) and forecasting results vs. number of selected eig-eiv on *MetrLA*.</caption><thead><tr><th></th><th>k = 2</th><th>k = 5</th><th>k = 10</th><th>k = 50</th><th>k = 100</th><th>k = 150</th><th>k = 200</th></tr></thead><tbody><tr><td>MAE</td><td>2.80</td><td>2.74</td><td>2.72</td><td>2.71</td><td>2.71</td><td>2.70</td><td>2.70</td></tr><tr><td>Training time (seconds)</td><td>2.52</td><td>2.91</td><td>2.96</td><td>3.86</td><td>5.28</td><td>6.36</td><td>7.78</td></tr></tbody></table>

the CITRUS blocks are 16, while no MLP modules have been used. Regarding NOAA, we set the dimension and the initial linear layer of the CITRUS blocks as 16, and use 3-layer MLPs for channel mixing in each block. The CITRUS blocks have the last activation as a Leaky ReLU. We use MSE as the loss function and consider root-normalized MSE (rNMSE) as the evaluation metric.

Table 2 illustrates the results of the weather forecasting task. We observe that CITRUS shows superior forecasting performance compared to previous methods in all forecasting horizons. We note that the SGP model [42], which is a scalable architecture, is the second-best performing method in Molene, possibly due to the small scale of this dataset. On the contrary, the NOAA dataset is larger, so models with higher parameter budgets perform better.

## 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 a regular GNN; ii) STT, which is exactly the opposite of the TTS; iii) Continuous TTS (CTTS), where we replace the GNN in TTS with a CGNN; and iv) Continuous STT (CSTT), which is the opposite of CTTS. We report the results of this ablation study on the *MetrLA* dataset in Table 3. We observe that CITRUS has superior results probably due to learning joint (and not sequential) spatiotemporal couplings by modeling these dependencies using product graphs with learnable receptive fields. We also observe that CTTS and CSTT outperform their discrete counterparts TTS and STT, probably due to the learning of adaptive graph neighborhoods instead of relying on 1-hop connections in regular GNNs.

**Hyperparameter sensitivity analysis.** The sensitivity analysis is related to the number of selected eigenvector-eigenvalue (eig-eiv) pairs of the factor Laplacians. We select eigenvector-eigenvalue pairs in $k \in \{2, 10, 50, 100, 150, 200\}$ out of 207 components in the *MetrLA* dataset. Table 4 presents the forecasting results in MAE along with the training time (per epoch) of this ablation study. We observe that selecting 50 eigenvector-eigenvalue pairs is enough to get good forecasting performance while having a fast training time. This experiment illustrates the advantages of CITRUS for accelerating training and keeping the performance as consistent as possible.

# 5 Conclusion and Limitations

In this paper, we proposed CITRUS, a novel model for jointly learning multidomain couplings from product graph signals based on tensorial PDEs on graphs (TPDEGs). We modeled these representations as separable continuous heat graph kernels as solutions to the TPDEG. Therefore, we showed that the underlying graph is actually the Cartesian product of the domain-specific factor
graphs. We rigorously studied the stability and over-smoothing aspects of CITRUS theoretically and experimentally. Finally, as a proof of concept, we use CITRUS to tackle the traffic and weather spatiotemporal forecasting tasks on public datasets, illustrating the superior performance of our model compared to state-of-the-art methods.

An interesting future direction for our framework could involve adapting it to handle other types of graph products, such as Kronecker and Strong graph products [23, 24]. Future efforts will also focus on finding tighter and more general upper bounds in the stability and over-smoothing analyses. For example, we will explore the possible relationships between the size of the factor graphs and these properties, either in a general form or for specific well-studied structures, such as ER graphs. Additionally, we plan to investigate more challenging real-world applications of the proposed framework, especially in scenarios involving more than two factor graphs.

## Acknowledgments

This work was supported by the center Hi! PARIS and ANR (French National Research Agency) under the JCJC project GraphIA (ANR-20-CE23-0009-01).

## References

[1] Y. L. Xu, K. Konstantinidis, and D. P. Mandic, "Multi-graph tensor networks," in *Advances in Neural Information Processing Systems - Workshop*, 2020.

[2] Y. L. Xu, K. Konstantinidis, and D. P. Mandic, "Graph tensor networks: An intuitive framework for designing large-scale neural learning systems on multiple domains," *arXiv preprint arXiv:2303.13565*, 2023.

[3] F. Monti, M. Bronstein, and X. Bresson, "Geometric matrix completion with recurrent multi-graph neural networks," in *Advances in Neural Information Processing Systems*, 2017.

[4] T. N. Kipf and M. Welling, "Semi-supervised 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 International 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, "Hierarchical coupled-geometry 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. Sabbaqi and E. Isufi, "Graph-time convolutional neural networks: Architecture and theoretical 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, "ProductGraphSleepNet: Sleep staging using product spatio-temporal graph learning with attentive temporal aggregation," *Neural Networks*, vol. 164, pp. 667-680, 2023.

[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. 5, 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 Y. 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: Representation 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. Klavžar, *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 on 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, "Understanding convolution on graphs via energies," *Transactions on Machine Learning Research*, 2023.

[32] T. K. Rusch, M. Bronstein, and S. Mishra, "A survey on oversmoothing in graph neural networks," *SAM Research Report*, 2023.

[33] M. Indibi and S. Aviyente, "Low-rank and smooth tensor recovery on Cartesian product graphs," in *IEEE International Conference on Sampling Theory and Applications*, 2023.

[34] J. S. Stanley, E. C. Chi, and G. Mishne, "Multiway graph signal processing on tensors: Integrative analysis of irregular geometries," *IEEE Signal Processing Magazine*, vol. 37, no. 6, pp. 160-173, 2020.

[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 P. 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, no. 3, 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-temporal graph modeling," in *International Joint Conference on Artificial Intelligence*, 2019.

[39] C. Zheng, X. Fan, C. Wang, and J. Qi, "GMAN: A graph multi-attention network for traffic prediction," in *AAAI Conference on Artificial Intelligence*, 2020.

[40] B. Yu, H. Yin, and Z. Zhu, "Spatio-temporal graph convolutional networks: a deep learning framework for traffic forecasting," in *International Joint Conference on Artificial Intelligence*, 2018.

[41] L. Ruiz, F. Gama, and A. Ribeiro, "Gated graph recurrent neural networks," *IEEE Transactions on Signal Processing*, vol. 68, pp. 6303-6318, 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 an isometric graph translation," in *IEEE European Signal Processing Conference*, 2015.

[44] A. Arguez, I. Durre, S. Applequist, R. S. Vose, M. F. Squires, X. Yin, R. R. Heim, and T. W. Owen, "NOAA's 1981-2020 US climate normals: an overview," *Bulletin of the American Meteorological Society*, vol. 93, no. 11, pp. 1687-1697, 2012.

[45] K. B. Petersen, M. S. 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 on Learning Representations*, 2023.
# 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 E), the case of inexact Cartesian product graphs (Section F), how to choose the appropriate number of eigenvalue-eigenvector pairs (Section G), the effect of graph receptive fields on the over-smoothing (Section H) and hyperparameter details (Section I), respectively.

## A Proofs

### A.1 Proof of Theorem 3.2

*Proof.* For obtaining the derivative of $\frac{\partial \underline{U}_t}{\partial t}$ from $e^{-t\underline{L}_p}$ (for $p = 1, \dots, P$) in (4), we first express $\frac{\partial \underline{U}_t}{\partial t}$ from the claimed solution in (4) as follows:

$$ \frac{\partial \underline{U}_t}{\partial t} = \sum_{p=1}^{P} \underline{U}_{t,p}, \qquad (19) $$

where

$$ \underline{U}_{t,p} := \underline{U}_0 \times_1 e^{-t\underline{L}_1} \times_2 \dots \times_{p-1} e^{-t\underline{L}_{p-1}} \times_p (-\underline{L}_p e^{-t\underline{L}_p}) \times_{p+1} e^{-t\underline{L}_{p+1}} \times_{p+2} \dots \times_P e^{-t\underline{L}_P}. \qquad (20) $$

Now, by performing mode-$p$ unfolding operation on $\underline{U}_{t,p}$, we have:

$$ \underline{U}_{t,p(p)} = -\underline{L}_p e^{-t\underline{L}_p} \underline{U}_{0(p)} (e^{-t\underline{L}_P} \otimes e^{-t\underline{L}_{P-1}} \otimes \dots \otimes e^{-t\underline{L}_{p+1}} \otimes e^{-t\underline{L}_{p-1}} \otimes \dots \otimes e^{-t\underline{L}_1})^{\top}. \qquad (21) $$

On the other hand:

$$ \underline{U}_{t(p)} = e^{-t\underline{L}_p} \underline{U}_{0(p)} (e^{-t\underline{L}_P} \otimes \dots \otimes e^{-t\underline{L}_{p+1}} \otimes e^{-t\underline{L}_{p-1}} \otimes \dots \otimes e^{-t\underline{L}_1})^{\top}. \qquad (22) $$

Therefore, by combining (21) and (22), one can write:

$$ \underline{U}_{t,p(p)} = -\underline{L}_p \underline{U}_t(p), \qquad (23) $$

which is the mode-$p$ unfolding of the following equation:

$$ \underline{U}_{t,p} = -\underline{U}_t \times_p \underline{L}_p, \qquad (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^{t\underline{L}_\diamond}$ on a Cartesian product graph $\underline{L}_\diamond := \bigoplus_{p=1}^P \underline{L}_p$ can be Kronecker-factored into its factor-based graph kernels as $e^{t\underline{L}_\diamond} = \bigotimes_{p=1}^P e^{t\underline{L}_p}$.

*Proof.* Using the pairwise summation property of the factor eigenvalues in a Cartesian product, one can write:

$$ \begin{aligned} e^{t\underline{L}_\diamond} &= e^{t(\bigoplus_{p=1}^P \underline{L}_p)} \\ &= (\bigoplus_{p=1}^P \underline{V}_p) e^{t(\bigoplus_{p=1}^P \underline{\Lambda}_p)} (\bigoplus_{p=1}^P \underline{V}_p)^{\top} \\ &= (\bigoplus_{p=1}^P \underline{V}_p) (\bigoplus_{p=1}^P e^{t\underline{\Lambda}_p}) (\bigoplus_{p=1}^P \underline{V}_p)^{\top} \\ &= \bigoplus_{p=1}^P e^{t\underline{L}_p}, \end{aligned} \qquad (25) $$

where $\underline{V}_p$, $\underline{\Lambda}_p^{(p)}$, and $\underline{\Lambda}_p$ are the eigenmatrix, $i$th eigenvalue, and the diagonal eigenvalue matrix corresponding to the $p$th Laplacian $\underline{L}_p$. Note that, in the related literature, this property has also been proved from different points of view, e.g., [34]. □
Next, by performing mode-(P + 1) unfolding on both sides of (5), one can write:

$$ (\mathcal{L}_l(\underline{\mathbf{U}}_l))_{(P+1)} = \mathbf{W}_l^\top \underbrace{\underline{\mathbf{U}}_l}_{\substack{\downarrow \\ \otimes_{p=1}^P}} \underbrace{e^{-t_l \underline{\mathbf{L}}_p}}_{\overbrace{\downarrow \otimes_{p=1}^P e^{-t_l \underline{\mathbf{L}}_p}}^{e^{-t_l \tilde{\mathbf{L}}_\circ}}} \top \quad (26) \\ \rightarrow [ (\mathcal{L}_l(\underline{\mathbf{U}}_l))_{(P+1)} ]^\top = \underbrace{\underline{\mathbf{U}}_l}_{\substack{\downarrow \\ \otimes_{p=1}^P}} \top \underbrace{\tilde{\mathbf{L}}_\circ}_{\overbrace{\downarrow \otimes_{p=1}^P e^{-t_l \underline{\mathbf{L}}_p}}^{e^{-t_l \tilde{\mathbf{L}}_\circ}}} \top \mathbf{W}_l. $$

□

### 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:

$$ \begin{aligned} \tilde{\mathbf{A}} &= \tilde{\mathbf{A}}_1 \otimes \mathbf{I}_2 + \mathbf{I}_1 \otimes \tilde{\mathbf{A}}_2 \\ &= (\mathbf{A}_1 + \mathbf{E}_1) \otimes \mathbf{I}_2 + \mathbf{I}_1 \otimes (\mathbf{A}_2 + \mathbf{E}_2) \\ &= \underbrace{(\mathbf{A}_1 \otimes \mathbf{I}_2 + \mathbf{I}_1 \otimes \mathbf{A}_2)}_{\mathbf{A}} + \underbrace{(\mathbf{E}_1 \otimes \mathbf{I}_2 + \mathbf{I}_1 \otimes \mathbf{E}_2)}_{\mathbf{E}} \\ &= \mathbf{A} + \mathbf{E}. \end{aligned} \quad (27) $$

Then, by assuming the theorem holds for the case of P = K and the definitions of $\tilde{\mathbf{A}}'_K := \bigoplus_{p=1}^K \tilde{\mathbf{A}}_p$ and $\tilde{\mathbf{A}}'_K := \bigoplus_{p=1}^K \tilde{\mathbf{A}}_p$, we next show it also holds for P = K + 1 as follows:

$$ \begin{aligned} \tilde{\mathbf{A}} &= \bigoplus_{p=1}^{K+1} \tilde{\mathbf{A}}_p = \underbrace{(\bigoplus_{p=1}^K \tilde{\mathbf{A}}_p)}_{\tilde{\mathbf{A}}'_K} \bigoplus \tilde{\mathbf{A}}_{K+1} = \tilde{\mathbf{A}}'_K \otimes \mathbf{I}_2 + \mathbf{I}_1 \otimes \tilde{\mathbf{A}}_{K+1} \\ &= (\mathbf{A}'_K + \mathbf{E}'_K) \otimes \mathbf{I}_2 + \mathbf{I}_1 \otimes (\mathbf{A}_{K+1} + \mathbf{E}_{K+1}) \\ &= \underbrace{(\mathbf{A}'_K \otimes \mathbf{I}_2 + \mathbf{I}_1 \otimes \mathbf{A}_{K+1})}_{\mathbf{A} = \bigoplus_{p=1}^{K+1} \mathbf{A}_p} + \underbrace{(\mathbf{E}'_K \otimes \mathbf{I}_2 + \mathbf{I}_1 \otimes \mathbf{E}_{K+1})}_{\mathbf{E} = \bigoplus_{p=1}^{K+1} \mathbf{E}_p} \\ &= \mathbf{A} + \mathbf{E}. \end{aligned} \quad (28) $$

For proving the upper bound of the norm of $\mathbf{E}$, in a similar approach to the previous proof, we first prove for the case of P = 2 as:

$$ \begin{aligned} \|\mathbf{E}\| &= \|\mathbf{E}_1 \otimes \mathbf{I}_2 + \mathbf{I}_1 \otimes \mathbf{E}_2\| \le \|\mathbf{E}_1 \otimes \mathbf{I}_2\| + \|\mathbf{I}_1 \otimes \mathbf{E}_2\| \\ &\le \underbrace{\|\mathbf{E}_1\|}_{\mathbf{E}'_1} \cdot \underbrace{\|\mathbf{I}_2\|}_{\mathbf{I}'_1} + \underbrace{\|\mathbf{I}_1\|}_{\mathbf{I}_1} \cdot \underbrace{\|\mathbf{E}_2\|}_{\mathbf{E}_2} = \varepsilon_1 + \varepsilon_2. \end{aligned} \quad (29) $$

Using a similar approach to the proof in (28), we assume the theorem holds for P = K as $\|\mathbf{E}\| \le \sum_{p=1}^K \varepsilon_p$, and, based on this information, then, we prove the theorem for the case of P = K + 1 as follows:

$$ \begin{aligned} \|\mathbf{E}\| &= \|\mathbf{E}'_K \otimes \mathbf{I}_2 + \mathbf{I}_1 \otimes \mathbf{E}_{K+1}\| \le \|\mathbf{E}'_K \otimes \mathbf{I}_2\| + \|\mathbf{I}_1 \otimes \mathbf{E}_{K+1}\| \\ &\le \underbrace{\sum_{p=1}^K \varepsilon_p}_{\|\mathbf{E}'_K\|} \cdot \underbrace{\|\mathbf{I}_2\|}_{\mathbf{I}'_1} + \underbrace{\|\mathbf{I}_1\|}_{\mathbf{I}_1} \cdot \underbrace{\|\mathbf{E}_{K+1}\|}_{\mathbf{E}_{K+1}} = \sum_{p=1}^{K+1} \varepsilon_p, \end{aligned} \quad (30) $$

which concludes the proof. □

### A.4 Proof of Theorem 3.7

*Proof.* The proof can be straightforwardly obtained by applying Proposition 1 in [30] based on the obtained results in Proposition 3.6. □
A.5 Proof of Lemma 3.9

*Proof.* By construction, we first prove the theorem for the case of $P = 2$. The defined Laplacian $\hat{\mathbf{L}}$ in (15) can be rewritten as:

$$
\begin{align*}
\hat{\mathbf{L}} = \mathbf{I} - \hat{\mathbf{A}} &= \mathbf{I} - \left( \frac{\hat{\mathbf{A}}_1}{2} \otimes \mathbf{I}_2 + \mathbf{I}_1 \otimes \frac{\hat{\mathbf{A}}_2}{2} \right) \\
&= \frac{1}{2} \left( 2\mathbf{I} - \left( \hat{\mathbf{A}}_1 \otimes \mathbf{I}_2 + \mathbf{I}_1 \otimes \hat{\mathbf{A}}_2 \right) \right) \\
&= \frac{1}{2} \left( \left[ \mathbf{I} - \left( \hat{\mathbf{A}}_1 \otimes \mathbf{I}_2 \right) \right] + \left[ \mathbf{I} - \left( \mathbf{I}_1 \otimes \hat{\mathbf{A}}_2 \right) \right] \right) \\
&= \frac{1}{2} \left( \left[ \mathbf{I}_1 \otimes \mathbf{I}_2 - \left( \hat{\mathbf{A}}_1 \otimes \mathbf{I}_2 \right) \right] + \left[ \mathbf{I}_1 \otimes \mathbf{I}_2 - \left( \mathbf{I}_1 \otimes \hat{\mathbf{A}}_2 \right) \right] \right) \tag{31} \\
&= \frac{1}{2} \left( \left[ \left( \mathbf{I}_1 - \hat{\mathbf{A}}_1 \right) \otimes \mathbf{I}_2 \right] + \left[ \mathbf{I}_1 \otimes \left( \mathbf{I}_2 - \hat{\mathbf{A}}_2 \right) \right] \right) \\
&= \frac{1}{2} \left( \hat{\mathbf{L}}_1 \otimes \mathbf{I}_2 + \mathbf{I}_1 \otimes \hat{\mathbf{L}}_2 \right) \\
&= \frac{1}{2} \left( \hat{\mathbf{L}}_1 \oplus \hat{\mathbf{L}}_2 \right) \\
&= \left( \frac{\hat{\mathbf{L}}_1}{2} \right) \oplus \left( \frac{\hat{\mathbf{L}}_2}{2} \right),
\end{align*}
$$

where we used the following property of Kronecker products: $\alpha(A \otimes B) = (\alpha A) \otimes B = A \otimes (\alpha B)$ [45]. Next, by assuming that the theorem holds for $P = K$ and the definition of $\hat{\mathbf{A}}'_K := \frac{1}{K}$ and $\oplus_{p=1}^K \hat{\mathbf{A}}_p := \oplus_{p=1}^K (\frac{\hat{\mathbf{A}}_p}{K})$, we prove it for $P = K + 1$:

$$
\begin{align*}
\hat{\mathbf{L}} = \mathbf{I} - \hat{\mathbf{A}} &= \mathbf{I} - \frac{1}{K+1} (K \hat{\mathbf{A}}'_K \oplus \hat{\mathbf{A}}_{K+1}) \\
&= \mathbf{I} - \left( \frac{K \hat{\mathbf{A}}'_K}{K+1} \otimes \mathbf{I}_2 + \mathbf{I}_1 \otimes \frac{\hat{\mathbf{A}}_{K+1}}{K+1} \right) \\
&= \frac{1}{K+1} \left( (K+1)\mathbf{I} - \left( K \hat{\mathbf{A}}'_K \otimes \mathbf{I}_2 + \mathbf{I}_1 \otimes \hat{\mathbf{A}}_{K+1} \right) \right) \\
&= \frac{1}{K+1} \left( K[\mathbf{I} - (\hat{\mathbf{A}}'_K \otimes \mathbf{I}_2)] + [\mathbf{I} - (\mathbf{I}_1 \otimes \hat{\mathbf{A}}'_{K+1})] \right) \\
&= \frac{1}{K+1} \left( \left[ \mathbf{I}_1 \otimes \mathbf{I}_2 - (\hat{\mathbf{A}}'_K \otimes \mathbf{I}_2) \right] + \left[ \mathbf{I}_1 \otimes \mathbf{I}_2 - (\mathbf{I}_1 \otimes \hat{\mathbf{A}}_2) \right] \right) \tag{32} \\
&= \frac{1}{K+1} \left( K \left[ (\mathbf{I}_1 - \hat{\mathbf{A}}'_K) \otimes \mathbf{I}_2 \right] + \left[ \mathbf{I}_1 \otimes (\mathbf{I}_2 - \hat{\mathbf{A}}_2) \right] \right) \\
&= \frac{1}{K+1} \left( K \left[ \oplus_{p=1}^K \left( \frac{\hat{\mathbf{L}}_p}{K} \right) \right] \otimes \mathbf{I}_2 + \mathbf{I}_1 \otimes \hat{\mathbf{L}}_2 \right) \\
&= \frac{1}{K+1} \left( \oplus_{p=1}^{K+1} (\hat{\mathbf{L}}_p) \right) \\
&= \oplus_{p=1}^{K+1} \left( \frac{\hat{\mathbf{L}}_p}{K+1} \right).
\end{align*}
$$

Besides, by the addition rule for the eigenvalues of a product graph adjacency or Laplacian ($\lambda_{ij}^{\diamond} = \lambda_i^{(1)} + \lambda_j^{(2)}$), it can be easily seen that the eigenvalues of the normalized product adjacency ($\tilde{\mathbf{W}}$) and Laplacian ($\hat{\mathbf{L}}$) in (31) are in the intervals of $[-1, 1]$ and $[0, 2]$, respectively, quite similar to the case for the normalized factor graphs, which concludes the proof. $\square$
## A.6 Proof of Theorem 3.10

*Proof.* First, Note that Using $\tilde{x}$ as the Graph Fourier Transform (GFT) [13] of $x$ w.r.t. $\hat{\mathbf{L}}$ (with eigenvalues $\{\lambda_i\}_{i=1}^N$), one can write [15]:

$$E(\mathbf{x}) = \mathbf{x}^{\top}\hat{\mathbf{L}}\mathbf{x} = \sum_{i=1}^{N} \lambda_i \tilde{x}_i^2. \qquad (33)$$

Next, by defining $\lambda$ as the smallest non-zero eigenvalue of the Laplacian $\hat{\mathbf{L}}$, the following lemma characterizes the over-smoothing aspects of applying a heat kernel in the simplest case.

**Lemma A.2.** We have:

$$E(e^{-\hat{\mathbf{L}}\mathbf{x}}) \le e^{-2\lambda} E(\mathbf{x}). \qquad (34)$$

*Proof.* By considering the EVD forms of $\hat{\mathbf{L}} = \mathbf{V}\mathbf{\Lambda}\mathbf{V}^{\top}$ and $e^{-\hat{\mathbf{L}}} = \mathbf{V}e^{-\mathbf{\Lambda}}\mathbf{V}^{\top}$, one can write

$$\begin{aligned} & E(e^{-\hat{\mathbf{L}}\mathbf{x}}) \\ & = \mathbf{x}^{\top} \underbrace{\mathbf{V}e^{-\mathbf{\Lambda}}\mathbf{V}^{\top}}_{\hat{\mathbf{L}}} \underbrace{\mathbf{V}\mathbf{\Lambda}\mathbf{V}^{\top}}_{\hat{\mathbf{L}}} \underbrace{\mathbf{V}e^{-\mathbf{\Lambda}}\mathbf{V}^{\top}}_{e^{-\hat{\mathbf{L}}} \mathbf{x}} = \sum_{i=1}^{N} \lambda_i \tilde{x}_i^2 e^{-2\lambda_i} \le e^{-2\lambda} \left( \sum_{i=1}^{N} \lambda_i \tilde{x}_i^2 \right) = e^{-2\lambda} E(\mathbf{x}). \end{aligned} \qquad (35)$$

Note that in the above proof, we ruled out the zero eigenvalues since they are useless in analyzing Dirichlet energy. □

Then by considering the following lemmas from [15]:

**Lemma A.3.** (Lemma 3.2 in [15]). $E(XW) \le \|W^{\top}\|_2^2 E(X)$.

**Lemma A.4.** (Lemma 3.3 in [15]). For ReLU and Leaky-ReLU nonlinearities $E(\sigma(X)) \le E(X)$.

Therefore, by combining the previous concepts and considering that the minimum non-zero eigenvalue of $\hat{\mathbf{L}}^{(t)}$ in (16) is $\frac{1}{P}t^{(m)}\lambda^{(m)}$ with $m = \arg\min_i t^{(i)}\lambda^{(i)}$, it can be said that:

**Theorem A.5.** For any $l \in \mathbb{N}_+$, we have $E(f_l(\mathbf{X})) \le s_l e^{-\frac{2}{P}\tilde{t}\tilde{\lambda}} E(\mathbf{X})$, where $s_l := \prod_{h=1}^{H_l} s_{lh}$ and $s_{lh}$ is the square of the maximum singular value of $\mathbf{W}_{lh}^{\top}$. Besides, $\tilde{\lambda} = \lambda^{(m)}$ and $\tilde{t} = t^{(m)}$ with $m = \arg\min_i t^{(i)}\lambda^{(i)}$, where $\lambda^{(i)}$ is the smallest non-zero eigenvalue of $\hat{\mathbf{L}}_i$.

Afterward, we can state our final corollary as follows:

**Corollary A.6.** Let $s := \sup_{l \in \mathbb{N}_+} s_l$. We have:

$$E(\mathbf{X}^{(l)}) \le s^l e^{-\frac{2}{P}l\tilde{t}\tilde{\lambda}} E(\mathbf{X}) = e^{l(\ln s - \frac{2}{P}\tilde{t}\tilde{\lambda})} E(\mathbf{X}). \qquad (36)$$

So, $E(\mathbf{X}^{(l)})$ exponentially converges to 0, when $\lim_{l \to \infty} e^{l(\ln s - \frac{2}{P}\tilde{t}\tilde{\lambda})} = 0$, i.e., $\ln s - \frac{2}{P}\tilde{t}\tilde{\lambda} < 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 $p_{ER}^{(1)} = p_{ER}^{(1)} = 0.3$, and varying $p_{ER}^{(3)} \in \{0.1, 0.3\}$ to monitor the performance in 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 5 provides the results across a varying number of layers, which shows that the proposed framework's performance is resistant to adding layers, especially in the case of more sparse graphs, i.e., $p_{ER}^{(3)} = 0.1$. On the other hand, the best results were obtained in $l = 1, 2, 4$, which are the nearest numbers to the actual number of layers in the true generative process which was 3. These observations are validated by comparing the results with the GCN, in which the over-smoothing and performance degradation happen severely faster than CITRUS.
<table><caption>Table 5: Experiments on more than two factor graphs ($p_{ER}^{(1)} = p_{ER}^{(2)} = 0.3$).</caption><thead><tr><th></th><th>l = 1</th><th>l = 2</th><th>l = 4</th><th>l = 8</th><th>l = 16</th><th>l = 32</th></tr></thead><tbody><tr><td>GCN, $p_{ER}^{(3)} = 0.1$</td><td>0.185±0.132</td><td>0.241±0.096</td><td>0.371±0.132</td><td>0.458±0.132</td><td>0.494±0.147</td><td>0.552±0.118</td></tr><tr><td>CITRUS, $p_{ER}^{(3)} = 0.1$</td><td>0.141±0.075</td><td>0.137±0.082</td><td>0.150±0.091</td><td>0.198±0.106</td><td>0.215±0.124</td><td>0.276±0.120</td></tr><tr><td>GCN, $p_{ER}^{(3)} = 0.3$</td><td>0.253±0.211</td><td>0.296±0.113</td><td>0.365±0.123</td><td>0.429±0.097</td><td>0.517±0.143</td><td>0.67±0.153</td></tr><tr><td>CITRUS, $p_{ER}^{(3)} = 0.3$</td><td>0.188±0.094</td><td>0.182±0.095</td><td>0.193±0.096</td><td>0.230±0.089</td><td>0.242±0.096</td><td>0.336±0.107</td></tr></tbody></table>

<table><caption>Table 6: Standard deviation of traffic forecasting comparison in Table 1.</caption><thead><tr><th rowspan="3">Method</th><th colspan="6">MetrLA</th></tr><tr><th colspan="3">H = 3</th><th colspan="3">H = 6</th></tr><tr><th>MAE</th><th>MAPE</th><th>RMSE</th><th>MAE</th><th>MAPE</th><th>RMSE</th></tr></thead><tbody><tr><td>CITRUS (Ours)</td><td>0.015</td><td>0.009</td><td>0.008</td><td>0.007</td><td>0.012</td><td>0.003</td><td>0.002</td><td>0.005</td><td>0.004</td></tr></tbody></table>

<table><thead><tr><th rowspan="3">Method</th><th colspan="6">PemsBay</th></tr><tr><th colspan="3">H = 3</th><th colspan="3">H = 6</th></tr><tr><th>MAE</th><th>MAPE</th><th>RMSE</th><th>MAE</th><th>MAPE</th><th>RMSE</th></tr></thead><tbody><tr><td>CITRUS (Ours)</td><td>0.009</td><td>0.008</td><td>0.003</td><td>0.013</td><td>0.006</td><td>0.005</td><td>0.008</td><td>0.011</td><td>0.008</td></tr></tbody></table>

## C Descriptions of the Compared Methods

* ARIMA [35]: uses Kalman filter to implement auto-regressive integrated moving average framework, but processes data in a timestep-wise manner.

* G-VARMA [25]: generalizes the vector auto-regressive moving average (VARMA) to the case of considering graph-based data by imposing some statistical assumptions on the model and relying on the spatial graph.

* 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 to them.

* Graph WaveNet [38]: a 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]: a hybrid of typical GCNs and 1D convolutional blocks for ST forecasting.

* GGRNN [41]: generalizes the typical RNNs by adapting GCN blocks instead of linear MLPs.

* GRUGCN [26]: a TTS framework, which first processes the temporal information by a GRU cell and then learns spatial representations by applying typical GCN blocks.

* SGP [42]: a 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 1 and weather forecasting regarding the proposed framework in Table 2 are provided in Tables 6 and 7.
<table><caption>Table 7: Standard deviation of weather forecasting comparison in Table 2.</caption><thead><tr><th rowspan="2">Method</th><th colspan="5">Molene</th><th colspan="5">NOAA</th></tr><tr><th>H = 1</th><th>H = 2</th><th>H = 3</th><th>H = 4</th><th>H = 5</th><th>H = 1</th><th>H = 2</th><th>H = 3</th><th>H = 4</th><th>H = 5</th></tr></thead><tbody><tr><td>CITRUS (ours)</td><td>0.0174</td><td>0.0155</td><td>0.0132</td><td>0.0107</td><td>0.0086</td><td>0.0034</td><td>0.0030</td><td>0.0026</td><td>0.0024</td><td>0.0023</td></tr></tbody></table>

<table><caption>Table 8: Intra and Inter-homophily measures on the real-world graphs of MetrLA and PemsBay datasets.</caption><thead><tr><th rowspan="2">Dataset</th><th colspan="3">Intra-graph homophily $\pi_{v_i}^s$</th><th colspan="3">Inter-graph homophily $\pi_{v_i}^T$</th></tr><tr><th>$p_s$</th><th>$q_p$</th><th>$q_n$</th><th>$p_s$</th><th>$q_p$</th><th>$q_n$</th></tr></thead><tbody><tr><td>MetrLA</td><td>0.2273</td><td>0.4732</td><td>0.2995</td><td>0.3325</td><td>0.4920</td><td>0.1755</td></tr><tr><td>PemsBay</td><td>0.1073</td><td>0.5912</td><td>0.3015</td><td>0.2399</td><td>0.6863</td><td>0.0738</td></tr><tr><td>Molene</td><td>0.0148</td><td>0.6248</td><td>0.3602</td><td>0.0152</td><td>0.5545</td><td>0.4301</td></tr><tr><td>NOAA</td><td>0.1249</td><td>0.5345</td><td>0.3405</td><td>0.2862</td><td>0.4351</td><td>0.2786</td></tr></tbody></table>

## E Homophily-Heterophily Trade-off across the Studied Real Datasets

Regarding the homophily-heterophily trade-off, first, we should mention that measuring the homophily-heterophily indices is way more challenging in spatiotemporal settings rather than in single-graph scenarios [46], since we are facing two intra-graph (i.e., the spatial dimension) and 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-heterophily behaviors. We provided the measured homophily-heterophily metrics (full details in [46]) in Table 8. As observed in this table, due to the wide range of variability across the used datasets, we have actually evaluated our performance on various datasets in terms of homophily-heterophily levels. For instance, Molene acts severely heterophilic in both intra and inter-graph domains, while NOAA acts more homophilic than Molene in the inter-graph domain.

## F The Case of Inexact Cartesian Product Modeling

This scenario lies under the general umbrella of inaccurate adjacencies [30], which has been rigorously studied in the stability analysis in Section 3.2. However, in a simplified and insightful case of facing graph products other than Cartesian, we consider two well-known graph products, i.e., Strong and Kronecker products [23, 24]. Consider the true graph is actually a Strong product graph $A = \tilde{A}_1 \otimes A_2 + A_1 \otimes I_2 + I_1 \otimes A_2$, and one mistakenly considers it as Cartesian product $\tilde{A} = A_1 \otimes I_2 + I_1 \otimes A_2$. Then, $E = A - \tilde{A} = A_1 \otimes A_2$, where $||E|| = \lambda_{max}^{(1)} \lambda_{max}^{(2)}$. So, the approximation error depends on the multiplication of the factor graph eigenspaces. Similarly, in the case of Kronecker product graphs $\tilde{A} = A_1 \otimes A_2$, it obtains $||E|| \le \lambda_{max}^{(1)} + \lambda_{max}^{(2)} + \lambda_{max}^{(1)} \lambda_{max}^{(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 an appropriate $k$ using two methodologies: supervised and unsupervised. For the supervised option, we can use cross-validation, as we did in our results. For the unsupervised method, we can analyze the Laplacian matrices for a low-rank approximation task. For example, in Figure 4, we plot the explained variances of the principal components of the spatial Laplacian (i.e., $||v\lambda v^T||_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.
<figure><img src="image_5.png" /><figcaption>Figure 4: Explained variance ratio vs. selected principal components.</figcaption></figure>

<table><caption>Table 9: Details of the training settings and hyperparameters, i.e., T (auto-regressive order), emb (dimension of embedding size in the spatiotemporal encoder), hid (dimension of linear mapping size in the spatiotemporal encoder), F_MLP (dimension of linear mapping size in the MLP layers), n_CITRUS (number of CITRUS blocks), F (second dimension of W_l in (6)), k_i (number of selected eigenvalue-eigenvector pairs in the i-th factor graph), lr (learning rate), n_batch (batch size), and n_epochs (number of epochs).</caption><thead><tr><th>Dataset</th><th>T</th><th>emb</th><th>hid</th><th>F_MLP</th><th>n_CITRUS</th><th>F</th><th>k<sub>1</sub></th><th>k<sub>2</sub></th><th>lr</th><th>n_batch</th><th>n_epochs</th></tr></thead><tbody><tr><td>MetrLA</td><td>6</td><td>16</td><td>32</td><td>64</td><td>3</td><td>64</td><td>205</td><td>4</td><td>0.01</td><td>2048</td><td>300</td></tr><tr><td>PemsBay</td><td>6</td><td>16</td><td>32</td><td>64</td><td>3</td><td>64</td><td>323</td><td>4</td><td>0.01</td><td>2048</td><td>300</td></tr><tr><td>Molene</td><td>10</td><td>16</td><td>16</td><td>16</td><td>3</td><td>16</td><td>30</td><td>8</td><td>0.001</td><td>64</td><td>1500</td></tr><tr><td>NOAA</td><td>10</td><td>4</td><td>4</td><td>4</td><td>3</td><td>4</td><td>107</td><td>8</td><td>0.01</td><td>256</td><td>400</td></tr></tbody></table>

## H The Effect of the Graph Receptive Fields on Controlling Over-smoothing

From a theoretical point of view, based on eq. (18), one intuitive way to alleviate or at least slow down the over-smoothing phenomena is to keep $\ln s - \frac{2}{P} \tilde{t} \lambda$ positive, which tends to $\tilde{t} < \frac{P}{2\lambda} \ln s$. This clarifies why one should not increase the node-neighborhood especially in the case of facing a connected graph or increasing the number of layers. This finding is also compatible with the pioneering research about the over-smoothing concept, e.g., [14]. From an experimental point of view, and if we consider the graph receptive field $t$ as a hyperparameter, we have designed an experiment to validate the claimed statement. Here, we vary $t \in \{0.1, 1, 5, 10, 20\}$ and monitor the over-smoothing process in Figure 5. In this figure, $b = \frac{P}{2\lambda} \ln s$, and as it is observed, for values of $t$ higher than $b$ the over-smoothing phenomena with increasing the number of layers is happening way faster than the other ones, which validates the theoretical discussion, too. Apart from this, the weight normalization technique [14,15] has also been mentioned as a helpful action but it hinders the theoretical foundations of CITRUS. Note that, in practice, only escaping from the over-smoothing is not enough to have good performance. Therefore, for $t < b$, one might need to perform other kinds of analyses, e.g., over-squashing analysis.

## I Hyperparameter Details

The detailed hyperparameters (optimized by cross-validation on the validation data) and/or training settings of CITRUS, i.e., T (auto-regressive order), emb (dimension of embedding size in the spatiotemporal encoder), hid (dimension of linear mapping size in the spatiotemporal encoder), F_MLP (dimension of linear mapping size in the MLP layers), n_CITRUS (number of CITRUS blocks), F (second dimension of W_l in (6)), k_i (number of selected eigenvalue-eigenvector pairs in the i-th factor graph), lr (learning rate), n_batch (batch size), and n_epochs (number of epochs), are provided in Table 9. Full details can be found in the implementation codes in https://github.com/ArefEinizade2/CITRUS.
<figure><img src="image_6.png" /><figcaption>Figure 5: Log relative distance vs. increasing the number of layers for different values of *t* in actual bounds.</figcaption></figure>

# 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 to better traffic management and reduced congestion, while improved weather forecasting can lead to serious urban inconvenience. This research can contribute to sustainable development goals by improving forecasting capabilities and decision-making processes, including sustainable cities and communities, climate action, and industry innovation and infrastructure.