Year: 2021

Abstract:
We introduce space-time graph neural network (ST-GNN), a novel GNN architecture, tailored to jointly process the underlying space-time topology of time-varying network data. The cornerstone of our proposed architecture is the composition of time and graph convolutional filters followed by pointwise nonlinear activation functions. We introduce a generic definition of convolution operators that mimic the diffusion process of signals over its underlying support. On top of this definition, we propose space-time graph convolutions that are built upon a composition of time and graph shift operators. We prove that ST-GNNs with multivariate integral Lipschitz filters are stable to small perturbations in the underlying graphs as well as small perturbations in the time domain caused by time warping. Our analysis shows that small variations in the network topology and time evolution of a system does not significantly affect the performance of ST-GNNs. Numerical experiments with decentralized control systems showcase the effectiveness and stability of the proposed ST-GNNs.

text: We introduce space-time graph neural network (ST-GNN), a novel GNN architecture, tailored to jointly process the underlying space-time topology of time-varying network data. The cornerstone of our proposed architecture is the composition of time and graph convolutional filters followed by pointwise nonlinear activation functions. We introduce a generic definition of convolution operators that mimic the diffusion process of signals over its underlying support. On top of this definition, we propose space-time graph convolutions that are built upon a composition of time and graph shift operators. We prove that ST-GNNs with multivariate integral Lipschitz filters are stable to small perturbations in the underlying graphs as well as small perturbations in the time domain caused by time warping. Our analysis shows that small variations in the network topology and time evolution of a system does not significantly affect the performance of ST-GNNs. Numerical experiments with decentralized control systems showcase the effectiveness and stability of the proposed ST-GNNs.

heading: 1 INTRODUCTION
text: Graph Neural Networks (GNNs) are powerful convolutional architectures designed for network data. GNNs inherit all the favorable properties convolutional neural networks (CNNs) admit, while they also exploit the graph structure. An important feature of GNNs, germane to their success, is that the number of learnable parameters is independent of the size of the underlying networks. GNNs have manifested remarkable performance in a plethora of applications, e.g., recommendation systems (Ying et al., 2018; Wu et al., 2020), drug discovery and biology (Gainza et al., 2020; Strokach et al., 2020; Wu et al., 2021; Jiang et al., 2021), resource allocation in autonomous systems (Lima et al., 2020; Cranmer et al., 2021), to name a few.
Recently, there has been an increased interest in time-varying network data, as they appear in various systems and carry valuable dynamical information. This interest is mostly prominent in applications as decentralized controllers (Tolstaya et al., 2020; Gama et al., 2020b; Yang and Matni, 2021; Gama and Sojoudi, 2021), traffic-flow forecasting (Yu et al., 2018; Li et al., 2018; Fang et al., 2021), and skeleton-based action detection (Yan et al., 2018; Cheng et al., 2020; Pan et al., 2021). The state-ofthe-art (SOTA) usually deploys an additional architecture side by side with a GNN so the latter learns patterns from the graph domain while the former works on the time sequences. One choice could be a CNN as in (Li et al., 2020; Isufi and Mazzola, 2021; Wang et al., 2021) or a recurrent neural network (RNN) as in (Seo et al., 2018; Nicolicioiu et al., 2019; Ruiz et al., 2020). However, these joint architectures are performed in a centralized manner in the sense that the up-to-date data of all nodes are given at any given time. While this is well suited for, but not limited to, the case of social networks and recommendation systems, many physical-network applications are decentralized in nature and suffer from time delays in delivering messages.
In this paper, we close the gap by developing a causal space-time convolutional architecture that jointly processes the graph-time underlying structure. That is, the convolutional layers preserve time delays in message passing. Our work is motivated by the following question. Is it possible to transfer learning between signals and datasets defined over different space-time underlying structures? This is a well motivated question since in practice we execute these architectures on graphs that are different from the graphs used in training and signals are sampled at different sampling rates between training and execution. The answer to the above question was provided for the case of static graph signals in (Gama et al., 2020a), where the stability of traditional GNNs to graph perturbations was proved. In
this work we give an affirmative answer to the above question in the case when time-varying graph signals are considered and space-time convolutional architectures are employed.
The contribution of this paper is twofold. First we introduce a novel convolutional architecture for time-varying graph signals, and second we prove its stability. Specifically, we provide a general definition of convolutions for any arbitrary shift operator and define a space-time shift operator (STSO) as the linear composition of the graph shift operator (GSO) and time-shift operator (TSO). We then introduce space-time graph neural networks (ST-GNNs), a cascade of layers that consist of space-time graph filters followed by point-wise nonlinear activation functions. The proposed ST-GNN allows processing continuous-time graph signals, which is pivotal in the stability analysis. Furthermore, we study the effect of relative perturbations on ST-GNNs and prove that small variations in the graph and/or irregularities in the sampling process of time-varying graph signals do not essentially affect the performance of the proposed ST-GNN architecture. Our theoretical findings are also supported by thorough experimental analysis based on decentralized control applications.
The rest of this paper is structured as follows. The related work is summarized in Section 2. Sections 3 and 4 present our contributions listed above. Numerical experiments and conclusions are presented in Sections 5 and 6, respectively. The proofs and extended experiments are provided in the appendices.
Notation: Bold small and large symbols, i.e. x and X, denote vectors and matrices, respectively. Calligraphic letters mainly represent time-varying graph signals unless otherwise is stated.

heading: 2 RELATED WORK
text: GNNs for Time-varying Graph Signals. One of the early architectures was (Yan et al., 2018), which introduced a convolutional filter that aggregates information only from 1-hop neighbors. Relying on product graphs, (Isufi and Mazzola, 2021) introduced a convolutional architecture, where each node has access to the present and past data of the other nodes in the graph. A similar convolutional layer was studied in (Pan et al., 2021; Loukas and Foucard, 2016), and while restricted to a domain that preserves time delays, (Pan et al., 2021) shows that it is not necessarily less expressive. Meanwhile, (Wang et al., 2021) performs GNNs over graph signals at each time instance separately before a temporal convolution is performed at the GNN outputs to capture the evolution of graph embeddings. Graph (RNNs) are another architecture developed to deal with time-varying graph signals. For example, (Pareja et al., 2020) uses an RNN to evolve the GNN parameters over time. Moreover, (Hajiramezanali et al., 2019) combined the GRNN with a variational autoencoder (VGAE) to improve the former’s expressive power. However, all these architectures did not take into account the physical restrictions in the form of time delays, associated with decentralized applications.
Similar to our work, other architectures considered the diffusion equation to form message-passing layers, e.g., (Xhonneux et al., 2020; Poli et al., 2021; Fang et al., 2021). These architectures parameterize the dynamics of graph signals with GNNs. In simple words, they learn a parameterization that helps find the current state variables from the previous ones. The architecture in (Chamberlain et al., 2021) is another example of these architectures, but learns the graph weights instead of a limited number of filter coefficients (i.e., it resembles the parameterization of graph attention networks in (Veličković et al., 2018)).
Stability. Deformation stability of CNNs was studied in (Bruna and Mallat, 2013; Bietti and Mairal, 2019). The notion of stability was then introduced to graph scattering transforms in (Gama et al., 2019; Zou and Lerman, 2020). In a following work, Gama et al. (2020a) presented a study of GNN stability to graph absolute and relative perturbations. Graphon neural networks was also analyzed in terms of its stability in (Ruiz et al., 2021). Moreover, (Pan et al., 2021) proves the stability to absolute perturbations of space-time graph scattering transforms.

heading: 3 SPACE-TIME GRAPH NEURAL NETWORKS
text: In this section, we present the proposed ST-GNN architecture for time-varying graph signals. First, we provide a general definition of convolutions and then we develop the space-time graph filters, which are the cornerstone of the ST-GNN architecture.
Our analysis starts with the homogeneous diffusion equation, which is defined with respect to the Laplacian differential operator L
∂X (t) ∂t
= −LX (t). (1) Equation (1) describes processes of signals that evolve across the diffusion dimension, t. The diffused signal, X (t), is modeled as an abstract vector in a vector space V. The vector space is associated with an inner product ⟨., .⟩V, and is closed under addition and scalar multiplication, e.g., n-dimensional vector spaces and function spaces. The solution of (1) is X (t) = e−tLX (0) = e−tLX0 = e−tLX (the subscript is omitted in what follows), and describes the signal after time t from the initial state.

heading: 3.1 CONVOLUTION
text: With the diffusion equation in our hands, we can now define the convolution operation as the linear combination of the diffusion sequence with a filter h(t). Definition 1 (Convolution Operator). For a linear shift-invariant filter with coefficients h(t), t ≥ 0, the convolution between the filter and an input signal X ∈ V, with respect to a linear differential operator L : V → V, is defined as
h ∗D X = ∫ ∞ 0 h(t)e−tLXdt, (2)
where ∗D denotes the convolution operator applied on signals with underlying structure D. Definition 1 is a valid convolution operator. With abuse of notation, we refer to L as the shift operator but it should be understood that the shift operation is executed using e−L. Definition 1 establishes a generalized form of convolutions for a wide variety of shift operators and signals. Next we discuss two convolution types of practical interest that follow directly from (2), i.e., time convolutions and graph convolutions.
Time Convolutions: They involve continuous-time signals, denoted by x(τ), that are elements of square-integrable function spaces, i.e, x(τ) ∈ L2(R). The underlying structure of x(τ) is the real line, since the time variable τ ∈ R. In order to generate the traditional definition of time convolution, the time shift operator (TSO) is chosen as the differential operator, Lτ = ∂/∂τ . One can show that e−t∂/∂τx(τ) = x(τ − t), using the Taylor series,
e−t ∂/∂τx(τ) = ∞∑ n=0 (−t)n n! · ∂ n ∂τn x(τ) = ∞∑ n=0 (u− τ)n n! · ∂ n ∂τn x(τ) = x(u),
with u = τ − t. In other words, the TSO performs a translation operation. Then (2) reduces to x(τ) ∗R h(τ) = ∫ ∞ 0 h(t)e−t ∂/∂τx(τ)dt = ∫ ∞ 0 h(t)x(τ − t)dt, (3)
which is indeed the traditional definition of convolutions between time signals. It is also worth noting that ejωτ , ∀ω ∈ [0,∞) are eigenfunctions of the TSO Lτ = ∂/∂τ with associated eigenvalues jω. This observation is pivotal in the stability analysis of section 4.
Graph Convolutions: They involve graph signals, x ∈ L2(RN ), that are N -dimensional squaresummable vectors containing information associated with N different nodes. The underlying structure of x is represented by a graph G = (V, E ,W), where V is the set of nodes with cardinality |V| = N , E ⊆ V × V is the set of edges, and W : E → R is a map assigning weights to the edges (Shuman et al., 2012). The graph shift operator (GSO), S ∈ RN×N , is a matrix representation of the graph sparsity, e.g., graph adjacency or Laplacian (Ortega et al., 2018). Applying S to a signal x results in a signal Sx, where [Sx]i represents aggregated data from all nodes j that satisfy that (i, j) ∈ E or (j, i) ∈ E . We focus on undirected graphs, and therefore, S is a real symmetric and diagonalizable matrix with a set of N eigenvectors {vi}Ni=1, each associated with a real eigenvalue λi. The GSO is a linear operator S : L2(RN ) → L2(RN ) that deals with discrete signals. Therefore, it is more common to discretize the diffusion process over graphs at times kTs, k ∈ Z+, where Ts is the sampling period. The convolution is then defined as
x ∗G h = K−1∑ k=0 hke −kSx, (4)
where {hk}k is the filter coefficients and K is the number of filter taps. The finite number of taps is a direct consequence of Cayley-Hamilton theorem, and h is a finite-impulse response (FIR) filter. Equation (4) generalizes graph convolutions in the literature, where convolutions are defined as polynomials in the shift operator, i.e., x ∗G h = ∑K−1 k=0 hkS
kx, (Sandryhaila and Moura, 2013). This definition executes the shift operation using the matrix S directly compared to e−S in (4).

heading: 3.2 SPACE-TIME CONVOLUTION AND GRAPH FILTERS
text: The generalized definition of convolution in 1, as well as the definitions of time and graph convolutions in (3), (4) respectively, set the ground for the space-time convolution. Space-time convolutions involve signals that are time-varying and are also supported on a graph. They are represented by X ∈ L2(RN ) ⊗ L2(R), where ⊗ is the tensor product [see (Kadison and Ringrose, 1983; Grillet, 2007)] between the vector spaces of graph signals and continuous-time signals. In the following, calligraphic symbols (except L) refer to time-varying graph signals unless otherwise is stated. The space-time shift operator (STSO) is the linear operator that jointly shifts the signals over the space-time underlying support. Therefore, we choose the STSO to be the linear composition of the GSO and TSO, i.e., L = S ◦ Lτ . The convolution operator in Definition 1 can be specified with respect to the STSO for time-varying graph signals as
h ∗G×R X = ∫ ∞ 0 h(t)e−tS◦LτXdt =: H(S,Lτ )X , (5)
where H(S,Lτ ) = ∫∞ 0
h(t)e−tS◦Lτ dt represents the space-time graph filter. The frequency response of the filter is then derived as [see Appendix A]
h̃(λ, jω) = ∫ ∞ 0 h(t)e−t(λ+jω)dt. (6)
Remark 1. The space-time graph filters can be implemented distributively. For each node, the operator ( e−S◦Lτ )t aggregates data from the t-hop neighbors after they are shifted in time. This process can be done locally after each node acquires information from their neighbors. The time shift t represents the time required for the data to arrive from the t-hop neighbors.

heading: 3.3 SPACE-TIME GRAPH NEURAL NETWORKS
text: The natural extension of our previous analysis is to define the convolutional layer of ST-GNNs based on the space-time graph filters in (5). However, learning infinite-impulse response (IIR) filters is impractical and thus we employ FIR space-time graph filters, defined as Hd(S,Lτ ) =∑K−1
k=0 hke −kTsS◦Lτ . The ST-GNN architecture employs a cascade of L layers, each of which consists of a bank of space-time graph filters followed by a pointwise nonlinear activation functions, denoted by σ. The input to layer ℓ is the output of the previous layer, Xℓ−1, and the output of the ℓth layer is written as
X fℓ = σ Fℓ−1∑ g=1 K−1∑ k=0 hfgkℓ e −kTsS◦LτX gℓ−1  , (7) for each feature X fℓ , f = 1, . . . , Fℓ. The number of features at the output of each layer is denoted by Fℓ, ℓ = 1, . . . , L. A concise representation of the ST-GNN can be written as Φ(X ;H,S◦Lτ ), where H is the set of all learnable parameters {hfgkℓ } f,g k,ℓ . An interesting remark is that the time-varying graph signals processed in (7) are continuous in time, but learning a finite number of parameters. Remark 2. The ST-GNNs are causal architectures, that is, the architecture does not allow the nodes to process data that are still unavailable to them. This is achieved by delaying the data coming from the k-hop neighbors by k time steps [c.f. (7)]. Therefore, each node has access to its own up-to-date data and outdated data from their neighbors.

heading: 4 STABILITY TO PERTURBATIONS
text: In this section, we perform a stability analysis of our proposed ST-GNN architecture. In particular, we characterize the difference between the filters H(S,Lτ ) and H(Ŝ, L̂τ ), where Ŝ and L̂τ are perturbed versions of the graph and time shift operators, respectively.

heading: 4.1 PERTURBATION MODELS
text: Before we study the effect of graph and time perturbations on the stability of our proposed ST-GNN, we first need to introduce a perturbation model in space and time. For graph perturbations, we follow the relative perturbation model proposed in (Gama et al., 2020a):
PT0 ŜP0 = S+ SE+ES, (8)
where Ŝ is the perturbed GSO, E is the error matrix and P0 is a permutation matrix. Taking a closer look in equation (8) we observe that the entries of S and Ŝ become more dissimilar as the norm ∥E∥ increases. The dissimilarity is measured by the difference [S]ij − [PT0 ŜP0]ij = [SE]ij + [ES]ij , where each summand is a weighted sum of the edges that are connected to nodes i and j, respectively. For nodes i and j with high node degrees, the sum includes a high number of non-zero entries from S, which leads to perturbations of high volume. Consequently, the considered graph perturbations are relative to node degrees. It is also worth noticing that equation (8) is invariant to node permutations of the GSO. This is due to the use of P0 which reconciles for node relabelling of perturbed graphs (see Section 4.2).
While there has been plenty of work regarding graph perturbation models, time perturbation in the context of graph neural networks is not well studied. In this paper we fill this gap and propose a practical model for time perturbations. In particular, we study irregularities that appear in the process of sampling continuous-time signals, x(τ). It is often the case that sampling is imperfect and not exactly equispaced, i.e., the discrete signals are captured at times kTs±z(kTs), where k ∈ Z, Ts ∈ R is the sampling period and z : R → R is a differentiable function. To model this phenomenon, we use time-warping operations. In particular, we consider a perturbed timeline τ̂ = τ + z(τ), and we observe the signal x(τ̂) = x(τ + z(τ)) instead of x(τ). The diffusion process in (1) for a continuous-time signal, over the perturbed timeline, can then be cast as
∂ ∂t x(t, τ̂) = − ∂ ∂τ̂ x(t, τ̂) ⇒ ∂ ∂t x(t, τ + z(τ)) = − (1 + z′(τ)) ∂ ∂τ x(t, τ + z(τ)). (9)
Now, let ξ(τ) = z′(τ), Lτ = ∂∂τ , and L̂τ = ∂∂τ̂ . Then the perturbed TSO takes the form L̂τ = (1 + ξ(τ)) Lτ . (10)
The effect of time perturbation is controlled by ξ(τ). As the norm ∥ξ(τ)∥2 grows, the dissimilarity between the original and perturbed time shift operators increases. This shows that the considered time perturbation model is relative to the rate of change of z(τ); faster changes in z(τ) produce more significant perturbations (i.e., longer time shifts), whereas smoother choices of z(τ) yield lower values of ∥ξ(τ)∥2 and lower perturbation levels. The time convolution in (3) with a perturbed TSO becomes x(τ + z(τ)) ∗R h(τ) = ∫ ∞ 0 h(t)e−t(1+ξ(τ))Lτx(τ + z(τ))dt
= ∫ ∞ 0 h(t)e−t(1+ξ(τ))Lτ ez(τ)Lτx(τ)dt, (11)
since x(τ + z(τ)) = ez(τ)Lτx(τ), which models a signal translation ez(τ)Lτ caused by time perturbations. Moreover, the perturbed TSO e−t(1+ξ(τ))Lτ models a time-varying shift of t(1 + ξ(τ)) when ξ(τ) is not a constant (which is the case in our analysis).

heading: 4.2 JOINT OPERATOR-DISTANCE MODULO
text: In order to proceed with the stability analysis, we need to define the operators used to measure the effect of perturbations. We start with the measure of graph perturbations as we utilize the operator distance modulo permutation, ∥.∥P , introduced in (Gama et al., 2020a). For any two graph operators A and Â : L2(RN ) → L2(RN ), the distance between them can be calculated as
∥A− Â∥P = min P∈P max x:∥x∥2=1 ∥Ax−PT ÂPx∥2, (12)
where P is the set of all permutation matrices. This distance measures how close an operator Â to being a permuted version of A. The distance in (12) is minimum (i.e., 0) when the two graphs are permuted versions of each other.
Next, we define a new operator distance modulo translation for time operations:
Definition 2 (Operator Distance Modulo Translation). Given the linear time operators B and B̂ : L2(R) → L2(R), define the operator distance modulo translation as
∥B − B̂∥T = min s∈R max x:∥x∥2=1 ∥∥∥Bx(τ)− (e−sLτ ◦ B̂)x(τ)∥∥∥ 2 . (13)
The expression in (13) measures how far or close two time operators are, to being translated versions of each other. In the case where L̂τ = (1 + ξ(τ)) Lτ as in (10), the distance between B = e−tLτ and B̂ = e−tL̂τ reduces to
∥e−tLτ − e−tL̂τ ∥T = min s∈R max x:∥x∥2=1 ∥x(τ − t)− x(τ − s− t− tξ(τ − s))∥2 . (14)
Since we are dealing with time-varying graph signals and joint time and graph perturbations, we combine the two previously defined distances and introduce the operator distance modulo permutationtranslation. Definition 3 (Operator Distance Modulo Permutation-Translation). Given the graph and time operators A, Â : L2(RN ) → L2(RN ), B, B̂ : L2(R) → L2(R), define the operator distance modulo permutation-translation as
∥A ◦ B − Â ◦ B̂∥P,T = min P∈P min s∈R max X :∥X∥F=1 ∥∥∥(A ◦ B)X − (PT ÂP) ◦ (e−sLτ ◦ B̂)X∥∥∥ F . (15)

heading: 4.3 STABILITY PROPERTIES OF SPACE-TIME GRAPH FILTERS
text: The last step before we present our stability results is to discuss integral Lipschitz filters, which are defined as: Definition 4 (Multivariate Integral Lipschitz Filters). Let the frequency response of a multivariate filter be h̃(w) : Rn → R. The filter is said to be integral Lipschitz if there exists a constant C > 0 such that for all w1 and w2, ∣∣∣h̃(w2)− h̃(w1)∣∣∣ ≤ 2C ∥w2 −w1∥2∥w2 +w1∥2 . (16) The condition in (16) requires the frequency response of the filter to be Lipschitz over line segments defined by arbitrary w1 and w2 ∈ Rn, with a Lipschitz constant that is inversely proportional to their midpoint. For differentiable filters the above condition reduces to
∣∣∣ ∂∂ζ h̃(w1)∣∣∣ ≤ C 1∥w1∥2 for every ζ that is a component of w1. The proposed space-time graph filters admit a frequency response which is bivariate with variables λ and jω [cf. (6)]. Following Definition 4, a space-time graph filter h̃(λ, jω) is integral Lipschitz, if
|λ+ jω| · ∣∣∣∣ ∂∂ζ h̃(λ, jω) ∣∣∣∣ ≤ C, ∀ζ ∈ {λ, jω}. (17) Note that the vector w1 in this case is given as [λ, ω]T with ∥w1∥2 = |λ+jω|. The conditions in (17) indicate that the frequency response of an integral Lipschitz space-time graph filter can vary rapidly at low frequencies close to 0. Therefore, the filter can discriminate between close low-frequency spectral components, but is more flat at high frequencies, prohibiting discriminability between the spectral features in these bands. Note that when we mention low frequencies in the context of space time graph filters (or ST-GNNs later), we refer to low values of λ and ω, whereas high frequencies correspond to high values of λ or ω. Proposition 1 (Stability to Graph Perturbations). Consider a space-time graph filter H along with graph shift operators S and Ŝ, and a time shift operator Lτ . If it holds that (A1) the GSOs are related by PT0 ŜP0 = S+ SE+ES with P0 being a permutation matrix, (A2) the error matrix E has a norm ∥E∥ ≤ ϵs, and an eigenvector misalignment δ relative to S, and (A3) the filter H(S,Lτ ) is an integral Lipschitz filter with a Lipschitz constant C > 0, then the difference between the space-time graph filters H(S,Lτ ) and H(Ŝ,Lτ ) is characterized by
∥H(S,Lτ )−H(Ŝ,Lτ )∥P ≤ 2Cϵs ( 1 + δ √ N ) +O(ϵ2s), (18)
where δ = (∥U−V∥+ 1)2 − 1, and N is the number of nodes in S.
Proposition 1 states that space-time graph filters are stable to graph perturbations of order O(ϵs) with a stability constant 2C(1 + δ √ N). The constant is uniform for all graphs of the same size and depends on the Lipschitz constant of the filter. The bound is proportional to the eigenvector misalignment between the GSO and the error matrix. The misalignment is bounded and does not grow with the graph size as δ ≤ (∥U∥+ ∥V∥+ 1)2 − 1 = 8 with U and V being unitary matrices. Now, consider a type of perturbations known as graph dilation, where the edge weights stretch by a value ϵs, i.e., the GSO becomes Ŝ = (1 + ϵs)S. In this case, the eigenvalues of the GSO also stretch by ϵs, but the eigenvectors are no longer misaligned, i.e., δ = 0. The bound for graph dilation then reduces to 2Cϵ and does not depend on the structure of the underlying graph. Therefore, it is inferred that the bound in (18) is split into (i) a term that reflects the difference in eigenvalues between S and Ŝ, and (ii) a term that arises from the eigenvector misalignment. Note that the stability constant is not affected by any property of the TSO, and consequently has the same form of the one derived by Gama et al. (2020a) for traditional graph filters. Proposition 2 (Stability to Time Perturbations). Consider a space-time graph filter H along with time shift operators Lτ and L̂τ , and a graph shift operator S. If it holds that (A1) the TSOs are related by L̂τ = (1 + ξ(τ)) Lτ , (A2) the error function ξ(τ) is infinitely differentiable with a norm ∥ξ(τ)∥2 ≤ κϵτ , and the norm of the mth-order derivative∥ξ(m)(τ)∥2 is of order O(ϵm+1τ ) where κ is an absolute constant, and (A3) the filter H(S,Lτ ) is an integral Lipschitz filter with a Lipschitz constant C > 0, then the difference between the space-time graph filters H(S,Lτ ) and H(S, L̂τ ) satisfies
∥H(S,Lτ )−H(S, L̂τ )∥T ≤ Cκϵτ +O(ϵ2τ ). (19) Similarly to Proposition 1, the filter difference is bounded by a constant proportional to the size of perturbation, i.e., ∥ξ(τ)∥2, and the bound is also affected by the Lipschitz constant of the filter. The perturbed TSO shares the same eigenfunctions with the original TSO since L̂τejωτ = jω(1 + ξ(τ))ejωτ . Therefore, the difference between the filters only arises from the difference in eigenvalues. It is worth mentioning that the assumption (A2) of Proposition 2 is rather reasonable. One example of a time-warping function z(τ) that satisfies this assumption is z(τ) = √ ϵτ cos(ϵττ)e
−ϵττ . The error function is then ξ(τ) = z′(τ) = −ϵ3/2τ (sin ϵττ + cos ϵττ)e−ϵττ with a norm of √ 3/4ϵτ . The first derivative of the error function is ξ′(τ) = −2ϵ5/2τ cos(ϵττ)e−ϵττ and its norm is of order O(ϵ2τ ), which increases exponentially with the order of the derivatives.
In Theorem 1, we consider joint time and graph perturbations, and characterize the difference between the filters under the STSO, S ◦ Lτ , and its perturbed version Ŝ ◦ L̂τ . Theorem 1 (Space-Time Graph Filter Stability). Under the assumptions of Propositions 1 and 2, the distance between the space-time graph filters H(S,Lτ ) and H(Ŝ, L̂τ ) satisfies
∥H(S,Lτ )−H(Ŝ, L̂τ )∥P,T ≤ 2Cϵs ( 1 + δ √ N ) + Cκϵτ +O(ϵ2), (20)
where ϵ = max{ϵs, ϵτ}. In order to keep this stability bound small, space-time graph filters should be designed with small Lipschitz constant C. This comes at the cost of filter discriminability at high frequencies [cf. (17)], whereas at low frequencies, the integral Lipschitz filters are allowed to vary rapidly. Thus Theorem 1 shows that space-time graph filters are stable and discriminative at low frequencies but cannot be both stable and discriminative at higher frequencies.

heading: 4.4 STABILITY PROPERTIES OF SPACE-TIME GRAPH NEURAL NETWORKS
text: Theorem 2 establishes the stability of the proposed ST-GNN architecture. Theorem 2 (ST-GNNs Stability). Consider an L single-feature ST-GNN, Φ(·;H,S ◦ Lτ ) along with the GSOs S and Ŝ, and the TSOs Lτ and L̂τ . If (A1) the filters at each layer are integral Lipschitz with a Lipschitz constant C > 0 and have unit operator norms, i.e., ∥Hℓ(S,Lτ )∥ = 1,∀ℓ = 1, . . . , L, (A2) the nonlinearirties σ are Lipschitz-continuous with a Lipschitz constant of 1, i.e., ∥σ(x2) −
σ(x1)∥2 ≤ ∥x2 − x1∥2,
(A3) the GSOs satisfy (A1) and (A2) of Proposition 1, and the TSOs satisfy (A1) and (A2) of Proposition 2, then
∥Φ(·;H,S ◦ Lτ )−Φ(·;H, Ŝ ◦ L̂τ )∥P,T ≤ 2CLϵs ( 1 + δ √ N ) + CLκϵτ +O(ϵ2), (21)
where ϵ = max{ϵs, ϵτ}. Theorem 2 shows that the stability bound of ST-GNNs depends on the number of layers L, in addition to the factors affecting the stability of space-time graph filters. Compared to space-time graph filters, ST-GNNs can be both stable and discriminative. This is due to the nonlinearities applied at the end of each layer. The effect of the pointwise nonlinearity is that it demodulates (i.e., spills) the energy of the higher-frequency components into lower frequencies [see (Gama et al., 2020a)]. Once the frequency content gets spread over the whole spectrum by the consecutive application of the nonlinearities, the integral Lipschiz filters in the deeper layers become discriminative. The stability/discriminability property of ST-GNNs is a pivotal reason why they outperform space-time graph filters.

heading: 5 NUMERICAL EXPERIMENTS
text: We examine two decentralized-controller applications: flocking, and motion planning. Specifically, we are given a network of N agents, where agent i has position pi,n ∈ R2, velocity vi,n ∈ R2 and acceleration ui,n ∈ R2 at time steps n < T . The agents collaborate to learn controller actions and complete a specific task. For each task, the goal is to learn a controller U that imitates an optimal centralized controller U∗. Therefore, we parameterize U with ST-GNNs and aim to find H∗ that satisfies
H∗ = arg min H∈H
1
M M∑ m=0 ℓ ( Φ ( Xm;H,Lm ) −U∗m ) , (22)
where H is the set of all possible parameterizations, ℓ(.) is the mean-squared loss function, Lm = (S ◦ Lτ )m is the STSO of the m-th example, and M is the size of the training dataset. The input signal Xm ∈ Rq×N×T contains q state variables of N agents, e.g. the position and velocity, over T time steps. Detailed description of the ST-GNN implementation can be found in Appendix G.
A. Network Consensus and Flocking. A swarm of N agents collaborates to learn a reference velocity rn ∈ R2,∀n, according to which the agents move to avoid collisions. Ideally the reference velocity and the agent velocity vi,n, should be the same, which is not the case in practice since each agent observes a noisy version of the reference velocity r̃i,n. Therefore the state variables (inputs) are the estimated velocities {vi,n}i,n, the observed velocities {r̃i,n}i,n and the relative masses {qi,n|qi,n = ∑ j∈Ni(pi,n − pj,n)}i,n, where Ni is the set of direct neighbors of the i-th agent. We train an ST-GNN to predict agent accelerations {ui,n}i,n according to which the agents will move to a new position with certain velocity. Further details can be found in Appendix G.
Experiment 1. First we organize 100 agents in a mesh grid, but the agents are not allowed to move. The reason is that we want to test the ST-GNN on a setting where the graph is not changing over time. Note that although the agents do not move they still estimate an acceleration and therefore a velocity according to which they should move. In simple words the goal is to estimate the optimal acceleration of each agent in the case they would be required to move. Consequently, the training and testing is executed over the same graph but with different state variables and outputs. Fig. 1 (Left) illustrates the mean of agent velocities, defined as vn = ∑N i=1 vi,n, along with the reference velocity rn. We observe that the agents succeed to reach consensus and follow the reference velocity. We also execute the previously trained ST-GNN on a setting with perturbed graph and time shift operators according to (8) and (10), respectively. Specifically, the error matrix E is diagonal with ∥E∥ ≤ ϵ, and z(τ) = √ ϵ cos(ϵτ)e−ϵτ . Fig. 1 (Middle) shows that the relative RMSE of the outputs {ui,n}i,n, defined as ∥Φ(X;H,S ◦ Lτ ) − Φ(X;H, Ŝ ◦ L̂τ )∥F /∥Φ(X;H,S ◦ Lτ )∥F . We observe that the relative RMSE is linearly proportional to the perturbation size ϵ, following the results of Theorem 2.
Experiment 2. We now allow 50 agents to move and form a communication network that changes with their movement, i.e., the graph of the agents varies over time. Consequently, the ST-GNN is trained on a dynamic graph [see Appendix G]. Fig. 1 (Right) shows the difference between the agent and reference velocities, vi,n and rn, averaged over N agents. The small standard deviation in the velocity difference (represented in red) indicates that the agents reached a consensus on their velocities. The figure also shows that the ST-GNN outperforms the decentralized policy, described in Appendix G. This experiment demonstrates the effectiveness of the proposed ST-GNN architecture
on dynamic graphs, even though our theoretical analysis is derived for fixed graphs. Further results are presented in Appendix G.
B. Unlabeled Motion Planning. The goal of this task is to assign N unlabeled agents to N target locations. The agents collaborate to find their free-collision trajectories to the assigned targets. We train an ST-GNN to learn the control accelerations {ui,n}i,n according to (22). The state variables for agent i are position {pi,n}n, velocity {vi,n}n, the position and velocity of agent’s M -nearest neighbors, and the position of the M -nearest targets at each time step n. The underlying graph is constructed as (i, j) ∈ En if and only if j ∈ Ni,n or i ∈ Nj,n, where Ni,n is the set of the M -nearest neighbors of agent i at time step n. The rest of the training parameters are shown in Appendix G.
We executed the learned ST-GNN on a test dataset that is designed according to the parameter M and sampling time Ts. A snapshot of the predicted trajectories are shown in Fig. 2 (Left). We observe that the agents approach their target goals without collisions. The average distance between the agent’s final position and the desired target is, d̂pg = 0.524 with variance equal to 0.367. We now test the sensitivity of ST-GNNs to the right choice of M and Ts in the testing phase. Fig. 2 (Middle) shows the relative error, which is calculated as (d̂pg,pert − d̂pg,org)/d̂pg,org, when we change M . The error increases with ∆M because changing the neighborhood size induces a change in the underlying graph. According to our theoretical analysis, the error increases with the size of graph perturbations, which matches the results in Fig. 2 (Middle). The same remark can be observed in Fig. 2 (Right) when we change the sampling time, which induces a change in the underlying time structure.

heading: 6 CONCLUSIONS
text: In this paper we developed a novel ST-GNN architecture, tailored for time-varying signals that are also supported on a graph. We showed that the proposed architecture is stable under both graph and time perturbations under certain conditions. Our conditions are practical and provide guidelines on how to design space-time graph filters to achieve desirable stability. Our theoretical analysis is supported by strong experimental results. Simulations on the tasks of decentralized flocking and unlabeled motion planning corroborated our theoretical results and demonstrated the effectiveness of the proposed ST-GNN architecture.

heading: A (PROOF OF) FREQUENCY RESPONSE OF SPACE-TIME GRAPH FILTERS
text: We start from the abstract representation of signal X ∈ V, on which we apply a shift operator L : V → V. Then we can write the signal X as
X = ∫ i ⟨X ,U(i)⟩V U(i)di, (23)
where {U(i)}i is the set of the orthonormal eigenfunctions of L that spans the space V. The form in (23) is a reminiscence of the inverse Fourier transform and the scalar values ⟨X ,U(i)⟩V are the Fourier coefficients. In classical signal processing, the Fourier coefficients represent the spectral components of the signal X . Now, we consider the example of time-varying graph signals, X ∈ L2(RN ) ⊗ L2(R), which are constructed as the tensor product of a continuous-time signal and a graph signal. Shifting the signal X is executed using the STSO, S ◦ Lτ , which is the linear composition of the GSO and TSO. The eigenfunctions of the STSO can then be defined as the tensor product vi ⊗ ejωτ , where vi is an eigenvector of the GSO and ejωτ is an eigenfunction of the TSO. Recall that the eigenfunction is itself a vector in the space L2(RN )⊗ L2(R), and thus applying the STSO on vi ⊗ ejωτ means to jointly shift the the graph vector and the continuous-time function. In other words, we have
e−S◦Lτ (vi ⊗ ejωτ ) = e−Svi ⊗ e−Lτ ejωτ = e−(λi+jω)(vi ⊗ ejωτ ). (24) Thus we conclude that every eigenfunction vi ⊗ ejωτ of the STSO is associated with an eigenvalue λi + jω.
Defining the eigenfunctions of the STSO allows us to re-write (23) for time-varying graph signals as
X = 1 2π N∑ i=1 ∫ ∞ 0 X̃ (λi, jω)(vi ⊗ ejωτ )dω, (25)
where X̃ (λi, jω) is the spectral component of the signal X at the frequency pair (λi, ω), and the constant 12π is to normalize the eigenfunctions. It is inferred that the frequency domain of timevarying graph signals is bivariate. Considering the space-time graph filter H(S,Lτ ), the filter output can be expressed as
Y = H(S,Lτ )X = ∫ ∞ 0 h(t)e−tS◦LτXdt
(a) =
1
2π N∑ i=1 ∫ ∞ 0 X̃ (λi, jω) ∫ ∞ 0 h(t)e−tS◦Lτ (vi ⊗ ejωτ )dtdω
(b) =
1
2π N∑ i=1 ∫ ∞ 0 X̃ (λi, jω) ∫ ∞ 0 h(t)e−t(λi+jω)dt(vi ⊗ ejωτ )dω,
(26)
where (a) results from expressing X by its spectral representation [cf. (25)] and (b) from recalling (24). As in (25), the filter output can be also written as Y = 12π ∑N i=1 ∫∞ 0
Ỹ(λi, jω)(vi ⊗ ejωτ )dω, and thus we have Ỹ(λi, jω) = X̃ (λi, jω) ∫∞ 0
h(t)e−t(λi+jω)dt, following from (26). From the convolution theorem, the convolution operator implies a multiplication in the spectral domain. Therefore, the frequency response of the space-time graph filter is
h̃(λ, jω) = ∫ ∞ 0 h(t)e−t(λ+jω)dt, (27)
which is the Laplace transform of the impulse response h(t). It is worth noting that the filter response depends on the filter coefficients, which is irrespective of the graph. The spectral coefficients of the filter applied to one specific graph are, however, obtained by instantiating the frequency response h̃(λ, jω) on its eigenvalues {λi}Ni=1.

heading: B PROOF OF PROPOSITION 1
text: First, we define the eigenvector misalignment between S and E in the following lemma (Gama et al., 2020a).
Lemma 1. Let S = VΛVH and E = UMUH with ∥E∥ ≤ ϵs. For any eigenvector vi of S, it holds that
Evi = mivi +Eivi, (28)
with ∥Ei∥ ≤ ϵδ, where δ = (∥U−V∥+ 1)2 − 1.
Now we proceed with the proof of Proposition 1.
Proof of Proposition 1. Proposition 1 bounds the distance between the space-time graph filters before and after adding graph perturbations. From (12), this distance can be evaluated as
∥H(S,Lτ )−H(Ŝ,Lτ )∥P = min P∈P max X :∥X∥F=1 ∥H(S,Lτ )X −PTH(Ŝ,Lτ )PX∥F
= min P∈P max X :∥X∥F=1
∥H(S,Lτ )X −H(PT ŜP,Lτ )X∥F , (29)
where the last step is due to H being permutation equivariant, i.e., PTH(S,Lτ )P = ∫ ∞ 0 h(t)PT e−tS◦LτPdt (a) = ∫ ∞ 0 h(t)e−tP TSP◦Lτ = H(PTSP,Lτ ). (30)
Equality (a) results from the fact that matrix P commutes with the exponential operator and als with the TSO. From the minimum operator in (29), there exists a matrix P0 ∈ P that satisfies
∥H(S,Lτ )−H(Ŝ,Lτ )∥P ≤ max X :∥X∥F=1 ∥H(S,Lτ )X −H(PT0 ŜP0,Lτ )X∥F
= ∥H(S,Lτ )−H(PT0 ŜP0,Lτ )∥ (b) = ∥H(S,Lτ )−H(S+ES+ SE,Lτ )∥,
(31)
where ∥.∥ is the operator norm, and (b) follows from the assumption (A1). Define the the filter difference as ∆s := H(S+ES+ SE,Lτ )−H(S,Lτ ). The next step is to find the norm of ∆s. ∆s can be expressed [cf. (5)] as
∆s = ∫ ∞ 0 h(t)e−t(S+ES+SE)◦Lτ dt− ∫ ∞ 0 h(t)e−tS◦Lτ dt. (32)
From Taylor series, e−tS◦Lτ = ∑∞
n=0 1 n! (−t)n(S ◦ Lτ )n and the difference is then written as
∆s = ∫ ∞ 0 h(t) ∞∑ n=0 (−t)n n! ( ((S+ES+ SE) ◦ Lτ )n − (S ◦ Lτ )n ) dt. (33)
By induction, we expand the term ((S+ES+ SE) ◦ Lτ )n to the first order on E, so we have
((S+ES+ SE) ◦ Lτ )n (a) = (S ◦ Lτ + (ES+ SE) ◦ Lτ )n
= (S ◦ Lτ )n + n−1∑ r=0 (S ◦ Lτ )r ( (ES+ SE) ◦ Lτ ) (S ◦ Lτ )n−r−1 +O2(E),
(34)
where O2(E) combines the other higher order terms, and (a) follows from ◦ being a linear composition. Hence, (33) reduces to
∆s = ∫ ∞ 0 h(t) ∞∑ n=0 (−t)n n! n−1∑ r=0 (S ◦ Lτ )r ( (ES+ SE) ◦ Lτ ) (S ◦ Lτ )n−r−1dt+O(E), (35)
with O(E) = ∫∞ 0 h(t) ∑∞ n=0 (−t)n n! O2(E)dt. Note that the frequency response of the filter is an analytic function. Thus the norm of O(E) is of order O(∥E∥2) since
0 < lim ∥E∥→0 ∥O(E)∥ ∥E∥2 < ∞. (36)
Denote the first term in the right hand side of (35) by ∆(S). After splitting the summands (ES + SE) ◦ Lτ in (35) into ES ◦ Lτ + SE ◦ Lτ , we get
∆(S) = ∫ ∞ 0 h(t) ∞∑ n=0 (−t)n n! n−1∑ r=0 (S ◦ Lτ )rE(S ◦ Lτ )n−rdt
+ ∫ ∞ 0 h(t) ∞∑ n=0 (−t)n n! n−1∑ r=0 (S ◦ Lτ )r+1E(S ◦ Lτ )n−r−1dt. (37)
The first line is straightforward and the second follows from commuting the matrix E with the TSO, Lτ . The required norm becomes ∥∆s∥ ≤ ∥∆(S)∥ + ∥O(E)∥ from the triangle inequality, with ∥O(E)∥ being of order O(ϵ2s) by the assumption (A2) and (36). Thus the aim of the proof reduces to bound the norm ∥∆(S)∥ by 2Cϵs(1 + δ √ N).
The norm ∥∆(S)∥ is defined as ∥∆(S)∥ = maxX :∥X∥F=1 ∥∆(S)X∥F . We express ∆(S)X as
∆(S)X = 1 2π N∑ i=1 ∫ ∞ 0 X̃ (λi, jω) ∫ ∞ 0 h(t) ∞∑ n=0 (−t)n n! n−1∑ r=0 (S ◦ Lτ )rE(S ◦ Lτ )n−r(vi ⊗ ejωτ )dtdω
+ 1
2π N∑ i=1 ∫ ∞ 0 X̃ (λi, jω) ∫ ∞ 0 h(t) ∞∑ n=0 (−t)n n! n−1∑ r=0 (S ◦ Lτ )r+1E(S ◦ Lτ )n−r−1(vi ⊗ ejωτ )dtdω,
(38)
which follows from expressing ∆(S) as in (37) and X as in (25). Then
∆(S)X = 1 2π N∑ i=1 ∫ ∞ 0 X̃ (λi, jω) ∫ ∞ 0 h(t) ∞∑ n=0 (−t)n n! n−1∑ r=0 (λi + jω) n−r(S ◦ Lτ )rE(vi ⊗ ejωτ )dtdω
+ 1
2π N∑ i=1 ∫ ∞ 0 X̃ (λi, jω) ∫ ∞ 0 h(t) ∞∑ n=0 (−t)n n! n−1∑ r=0 (λi + jω) n−r−1(S ◦ Lτ )r+1E(vi ⊗ ejωτ )dtdω,
(39)
since (S ◦ Lτ )q(vi ⊗ ejωτ ) = (λi + jω)q(vi ⊗ ejωτ ). The assumption (A2) indicates that matrix E has eigenvectors that are not aligned with vi ∀i. From Lemma 1, it follows that E(vi ⊗ ejωτ ) = Evi ⊗ ejωτ = mivi ⊗ ejωτ +Eivi ⊗ ejωτ . The terms (S ◦ Lτ )pE(vi ⊗ ejωτ ) of (39) can then be simplified as
(S ◦ Lτ )pE(vi ⊗ ejωτ ) = mi(S ◦ Lτ )p(vi ⊗ ejωτ ) + (S ◦ Lτ )pEi(vi ⊗ ejωτ ) = mi(λi + jω) p(vi ⊗ ejωτ ) + (S ◦ Lτ )pEi(vi ⊗ ejωτ ). (40)
Substituting the first term of (40) in (39) results in
∆1(S)X := 1
2π N∑ i=1 ∫ ∞ 0 X̃ (λi, jω) ∫ ∞ 0 h(t) ∞∑ n=0 (−t)n n! n−1∑ r=0 2mi(λi + jω) n(vi ⊗ ejωτ )dtdω.
(41)
Doing the same for the second term of (40) leads to
∆2(S)X := 1
2π N∑ i=1 ∫ ∞ 0 X̃ (λi, jω) ∫ ∞ 0 h(t) ∞∑ n=0 (−t)n n! n−1∑ r=0 (λi + jω) n−r(S ◦ Lτ )rEi(vi ⊗ ejωτ )dtdω
+ 1
2π N∑ i=1 ∫ ∞ 0 X̃ (λi, jω) ∫ ∞ 0 h(t) ∞∑ n=0 (−t)n n! n−1∑ r=0 (λi + jω) n−r−1(S ◦ Lτ )r+1Ei(vi ⊗ ejωτ )dtdω,
(42)
where ∆(S)X = ∆1(S)X +∆2(S)X . The next step is to find the norm of ∆1(S) and ∆2(S). For the term ∆1(S)X in (41), we notice that the inner summation no longer depends on r and can be written as 2nmi(λi + jω)n(vi ⊗ ejωτ ). Therefore, we get
2mi ∫ ∞ 0 h(t) ∞∑ n=0 n (−t)n n! (λi + jω) ndt = 2mi(λi + jω) ∂ ∂λ h̃(λi, jω), (43)
where ∂∂λ h̃(λi, jω) = ∫∞ 0 h(t) ∑∞ n=0 n n! (−t)n(λi + jω)n−1dt is the partial derivative of the frequency response of the filter computed at λ = λi. Thus (41) is written as
∆1(S)X = 1
2π N∑ i=1 ∫ ∞ 0 2miX̃ (λi, jω) (λi + jω) ∂ ∂λ h̃(λi, jω) (vi ⊗ ejωτ )dω. (44)
From Pythagoras’ theorem we know that the squared norm of a sum of orthogonal terms is the sum of the squared norm of the individual summands. Then
∥∆1(S)X∥2F = N∑ i=1 ∫ ∞ 0 4m2i |X̃ (λi, jω)|2 · |λi + jω|2 · ∣∣∣∣ ∂∂λh̃(λi, jω) ∣∣∣∣2 dω ≤ 4ϵ2sC2∥X∥2F , (45) since |mi| ≤ ∥E∥ ≤ ϵs, the filter is integral Lipschitz [cf. (17)], and the eigenfunction basis vectors are orthonormal, i.e., ∥vk ⊗ 12π ejω∗τ∥2F = 1. The norm ∥X∥2F is defined as∑N
i=1 ∫∞ 0
|X̃ (λi, jω)|2dω. We now rewritte ∆2(S)X in (42) as
∆2(S)X = 1
2π N∑ i=1 ∫ ∞ 0 X̃ (λi, jω)KiEi(vi ⊗ ejωτ )dω, (46)
with Ki = ∫ ∞ 0 h(t) ∞∑ n=0 (−t)n n! n−1∑ r=0 ( (λi + jω) n−r(S ◦ Lτ )r + (λi + jω)n−r−1(S ◦ Lτ )r+1 ) dt. (47)
We next find the norm of Ki in order to find ∥∆2(S)∥. Note that Ki : L2(RN ) ⊗ L2(R) → L2(RN ) ⊗ L2(R) is a linear operator and so is the matrix S : L2(RN ) → L2(RN ). A matrix S can be expressed using the outer product of its eigenvectors as S = ∑N i=1 λiviv H i . Since the tensor product generalizes the outer product of n-dimensional vectors, we can draw a parallel and have
(S ◦ Lτ )q = 1
(2π)2 N∑ k=1 ∫ ∞ 0 (λk + jω∗) q(vk ⊗ ejω∗τ )⊗ (vk ⊗ ejω∗τ )Hdω∗. (48)
Calculating (47) requires evaluating the term (λi + jω)n−r(S ◦ Lτ )r + (λi + jω)n−r−1(S ◦ Lτ )r+1 first, which can be simplified as
(λi + jω) n−r(S ◦ Lτ )r + (λi + jω)n−r−1(S ◦ Lτ )r+1
= (λi + jω) n ( (S ◦ Lτ )r (λi + jω)r + (S ◦ Lτ )r+1 (λi + jω)r+1 ) .
(49)
When we substitute (48) in (49), we get
(λi + jω) n−r(S ◦ Lτ )r + (λi + jω)n−r−1(S ◦ Lτ )r+1
= (λi + jω)
n
(2π)2 N∑ k=1 ∫ ∞ 0 ( (λk + jω∗) r (λi + jω)r + (λk + jω∗) r+1 (λi + jω)r+1 ) (vk ⊗ ejω∗τ )⊗ (vk ⊗ ejω∗τ )Hdω∗.
(50)
The inner summation in (47) reduces to the sum of two geometric series, each of which has the form
n−1∑ r=0 ( λk + jω∗ λi + jω )r =
1 (λi + jω)n−1 (λi + jω) n − (λk + jω∗)n λi + jω − (λk + jω∗) , (51)
and the reader can confirm with a simple algebraic manipulation that the right hand side of (51) follows from the geometric sum ∑n−1 r=0 a
r = (1 − an)/(1 − a). It is also straightforward to show that
(λi+jω) n n−1∑ r=0 ( λk + jω∗ λi + jω )r + ( λk + jω∗ λi + jω )r+1 = λi + λk + j(ω + ω∗) λi − λk + j(ω − ω∗) ( (λi+jω) n−(λk+jω∗)n ) . (52) with some algebraic manipulations. We eventually can write Ki in (47) as
Ki = ∫ ∞ 0 h(t) ∞∑ n=0 (−t)n (2π)2 n! N∑ k=1 ∫ ∞ 0 λi + λk + j(ω + ω∗) λi − λk + j(ω − ω∗) ( (λi + jω) n − (λk + jω∗)n )
(vk ⊗ ejω∗τ )⊗ (vk ⊗ ejω∗τ )Hdω∗dt. (53)
Since we have h̃(λi, jω) = ∫∞ 0 h(t)e−t(λi+jω)dt = ∫∞ 0 h(t) ∑∞ n=0 1 n! (−t)n(λi + jω)ndt from Taylor series, Ki reduces to
Ki = 1
(2π)2 N∑ k=1 ∫ ∞ 0 (λi + jω) + (λk + jω∗) (λi + jω)− (λk − jω∗) ( h̃(λi, jω)− h̃(λk, jω∗) ) (vk ⊗ ejω∗τ )⊗ (vk ⊗ ejω∗τ )Hdω∗.
(54)
Note that (54) has the form (w2 +w1)(h̃(w2)− h̃(w1))/(w2 −w1) of (16) with w being replaced by its complex form, λ+ jω. This remark will help bound the norm of Ki. The norm ∥Ki∥ is defined as maxX :∥X∥F=1 ∥KiX∥F for all i, and we have
KiX (a) = Ki
( 1
2π N∑ m=1 ∫ ∞ 0 X̃ (λm, jω̇)(vm ⊗ ejω̇τ )dω̇ )
(b) =
1
2π N∑ k=1 ∫ ∞ 0 X̃ (λk, jω∗) (λi + jω) + (λk + jω∗) (λi + jω)− (λk − jω∗) ( h̃(λi, jω)− h̃(λk, jω∗) ) (vk ⊗ ejω∗τ )dω∗.
(55)
The RHS in (a) involves the inner products 1(2π)2 (vk ⊗ ejω∗τ )H(vm ⊗ ejωτ ) = 1(2π)2 ⟨vk ⊗ ejω∗τ ,vm ⊗ ejωτ ⟩, ∀m, k. These inner products are between the orthogonal eigenvectors and are nonzero if and only if m = k and ω̇ = ω∗. Therefore, only the terms that have m = k and ω̇ = ω∗ appeared in (b), and the inner products are 1.
From Pythagoras’ theorem, the squared norm of KiX is then
∥KiX∥2F = N∑
k=1
∫ ∞ 0 |X̃ (λk, jω∗)|2 ∣∣∣∣ (λi + jω) + (λk + jω∗)(λi + jω)− (λk − jω∗) ( h̃(λi, jω)− h̃(λk, jω∗) )∣∣∣∣2 dω∗ ≤ 4C2∥X∥2F , ∀i,
(56)
followed from the assumption (A3) that the filter is Lipschitz filter [cf. Definition 4]. Therefore, we get ∥Ki∥ ≤ 2C for a unit-norm signal X . Now, it is left to take the norm of (46), which yields
∥∆2(S)X∥F = ∥∥∥∥∥ 12π N∑ i=1 ∫ ∞ 0 X̃ (λi, jω)KiEi(vi ⊗ ejωτ )dω ∥∥∥∥∥ F
≤ N∑ i=1 ∫ ∞ 0 |X̃(λi, jω)| · ∥Ki∥ · ∥Ei∥dω ≤ 2Cϵδ N∑ i=1 ∫ ∞ 0 |X̃(λi, jω)|dω, (57)
where we have ∥Ei∥ ≤ ϵδ, ∀i from Lemma 1. To bound the summation in (57), we can write N∑ i=1 ∫ ∞ 0 |X̃ (λi, jω)| · 1 dω (a) ≤ ( N∑ i=1 ∫ ∞ 0 |X̃ (λi, jω)|2dω ) 1 2 ( N∑ i=1 ∫ ∞ 0 1 dω ) 1 2
(b) = ( N∑ i=1 ∫ ∞ 0 |X̃ (λi, jω)|2dω ) 1 2 ( N∑ i=1 ∫ ∞ 0 δ̂(τ)dτ ) 1 2 = √ N∥X∥F ,
(58)
where (a) follows from Cauchy-Schwartz inequality, and (b) from Parseval’s theorem. Note that δ̂(τ) is the Dirac delta function. Eventually, we obtain
∥∆2(S)X∥F ≤ 2Cϵsδ √ N∥X∥F . (59)
Recall that ∆(S)X = ∆1(S)X +∆2(S)X . From the triangle inequality, we then obtain ∥∆(S)X∥F ≤ ∥∆1(S)X∥F + ∥∆2(S)X∥F ≤ 2Cϵs(1 + δ √ N)∥X∥F . (60)
following from and using (45) and (59). Hence, ∥∆(S)∥ = maxX :∥X∥F=1 |∆(S)X∥F ≤ 2Cϵs(1 + δ √ N), which completes the proof.

heading: C PROOF OF PROPOSITION 2
text: The aim of the proof is to bound the distance between the space-time graph filters before and after adding time perturbation. This distance is calculated using the operator distance modulo translation. From Definition 3, it holds that there exists a real value s such that
∥H(S,Lτ )−H(S, L̂τ )∥T = min s∈R max x:∥x∥F=1 ∥∥∥H(S,Lτ )X − e−sLτH(S, L̂τ )X∥∥∥ F
(a) ≤ max X :∥X∥F=1
∥H(S,Lτ )X −H(S, L̂τ )X∥F (b) = ∥H(S,Lτ )−H(S, (1 + ξ(τ))Lτ )∥,
(61)
where in (a) we let the minimum norm value be lower than or equal to the norm when s = 0, and in (b) we use assumption (A1). The difference between the filters is evaluated as
H(S, L̂τ )−H(S,Lτ ) = ∫ ∞ 0 h(t)e−tS◦(Lτ+ξ(τ)Lτ )dt− ∫ ∞ 0 h(t)e−tS◦Lτ dt. (62)
Let ∆τ denote H(S, L̂τ )−H(S,Lτ ), which can be written as
∆τ = ∫ ∞ 0 h(t) ∞∑ n=0 (−t)n n! ( (S ◦ (Lτ + ξ(τ)Lτ ))n − (S ◦ Lτ )n ) dt, (63)
following from expanding the exponentials to e−tS◦Lτ = ∑∞
n=0 1 n! (−t)n(S ◦ Lτ )n using Taylor
series. We also expand (S ◦ Lτ + S ◦ ξ(τ)Lτ )n to the first order of ξ(τ):
(S ◦ Lτ + S ◦ ξ(τ)Lτ )n = (S ◦ Lτ )n + n−1∑ r=0 (S ◦ Lτ )rξ(τ)(S ◦ Lτ )n−r +O2(ξ(τ)), (64)
where O2(ξ(τ)) is a polynomial of the higher powers of ξ(τ). By substituting (64) in (63), we can re-write the latter as
∆τ = ∫ ∞ 0 h(t) ∞∑ n=0 (−t)n n! n−1∑ r=0 (S ◦ Lτ )rξ(τ)(S ◦ Lτ )n−rdt+O(ξ(τ)), (65)
where O(ξ(τ)) = ∫∞ 0 h(t) ∑∞ n=0 (−t)n n! O2(ξ(τ))dt. The quantity O(ξ(τ)) has all the higher power terms and satisfies
0 < lim ∥ξ(τ)∥→0 ∥O(ξ(τ))∥2 ∥ξ(τ)∥22 < ∞. (66)
Therefore, the norm ∥O(ξ(τ))∥2 is of order O(ϵ2τ ) following from assumption (A2). Remember that if the norm is of order ϵτ , it means that ∥ξ(τ)∥2 ≤ κϵτ with κ being an absolute constant. Our goal reduces to bound the norm of the first term of the right hand side of (65).
Since the TSO and GSO can commute, it holds that
(S ◦ Lτ )rξ(τ)(S ◦ Lτ )n−r = (Sr ◦ Lrτ )ξ(τ)(Sn−r ◦ Ln−rτ ) = Sn ◦ Lrτξ(τ)Ln−rτ . (67) However, Lτ and ξ(τ) do not commute due to their dependence on τ . Recalling that Lτ is a differential operator and applying the chain rule, we obtain
Lrτ ( ξ(τ)Ln−rτ ) (a) = r∑ m=0 ( r m ) ξ(m)(τ)Ln−mτ = ξ(τ)Lnτ +G2(ξ′(τ)), (68)
where (a) is valid by induction, ξ(m)(τ) = ∂mξ(τ)/∂τm and the term G2(ξ′(τ)) contains the higherorder derivatives starting from the first derivative. Substituting (67) and (68) in (65), we get
∆τ = ∫ ∞ 0 h(t) ∞∑ n=0 (−t)n n! n−1∑ r=0 Sn ◦ ξ(τ)Lnτ dt+G(ξ′(τ)) +O(ξ(τ))
(a) = ∫ ∞ 0 h(t) ∞∑ n=0 (−t)n n! nξ(τ)(S ◦ Lτ )ndt+G(ξ′(τ)) +O(ξ(τ)). (69)
The summands of the inner summation in (69) no longer depend on r leading to the form in (a). We also have G(ξ′(τ)) = ∫∞ 0 h(t) ∑∞ n=0 (−t)n n! ∑n−1 r=0 G2(ξ
′(τ))dt. It can now be shown that the norm ∥G(ξ′(τ))∥2 is of order O(ϵ2τ ) since
0 < lim ∥ξ′(τ)∥2→0 ∥G(ξ′(τ))∥2 ∥ξ′(τ)∥2 < ∞, (70)
and ∥ξ′(τ)∥2 is of order ϵ2τ according to assumption (A2). Denote the first term in the right hand side of (69) by ∆(Lτ ). Since both ∥O(ξ(τ))∥ and ∥G(ξ′(τ))∥ are of order ϵ2τ , the rest of the proof aims to show that the norm of ∆(Lτ ) is bounded by Cκϵτ . We have
∆(Lτ )X = (∫ ∞
0
h(t) ∞∑ n=0 (−t)n n!
nξ(τ)(S ◦ Lτ )ndt )( 1
2π N∑ i=1 ∫ ∞ 0 X̃ (λi, jω)(vi ⊗ ejωτ )dω )
= 1
2π N∑ i=1 ∫ ∞ 0 X̃ (λi, jω) ∫ ∞ 0 h(t) ∞∑ n=0 (−t)n n! nξ(τ)(S ◦ Lτ )n(vi ⊗ ejωτ )dtdω.
(71)
Recalling that (vi ⊗ ejωτ ) is an eigenfunction of (S ◦ Lτ ), we can re-write (71) as
∆(Lτ )X = 1
2π N∑ i=1 ∫ ∞ 0 X̃ (λi, jω)ξ(τ) ∫ ∞ 0 h(t) ∞∑ n=0 n n! (−t)n(λi + jω)ndt(vi ⊗ ejωτ )dω
= 1
2π N∑ i=1 ∫ ∞ 0 X̃ (λi, jω)ξ(τ)(λi + jω) ∂ ∂jω h̃(λi, jω)(vi ⊗ ejωτ )dω,
(72)
where ∂∂jω h̃(λi, jω) = ∫∞ 0 h(t) ∑∞ n=0 n n! (−t)n(λi + jω)n−1dt. The the squared norm of (72) can then be evaluated as
∥∆(Lτ )X∥2F (a) ≤ ∫ ∞ 0 N∑ i=1 ∫ ∞ 0 ∣∣∣X̃ (λi, jω)∣∣∣2 |ξ(τ)|2|λi + jω|2 ∣∣∣∣ ∂∂jω h̃(λi, jω) ∣∣∣∣2 dωdτ
= |λi + jω|2 ∣∣∣∣ ∂∂jω h̃(λi, jω) ∣∣∣∣2(∫ ∞ 0 |ξ(τ)|2dτ )(∫ ∞ 0 N∑ i=1 ∣∣∣X̃ (λi, jω)∣∣∣2 dω)
= |λi + jω|2 ∣∣∣∣ ∂∂jω h̃(λi, jω) ∣∣∣∣2 ∥ξ(τ)∥22 ∥X∥2F (b)≤ C2κ2ϵ2τ∥X∥2F . (73)
In (a), we use the inequality | ∫∞ 0 xdx|2 ≤ ∫∞ 0
|x|2dx. In (b), we use assumption (A3), which states that the filter is integral Lipschitz with a constant C, and assumption (A2), which has the norm ∥ξ(τ)∥2 bounded by κϵτ . Finally, the required norm can be calculated as ∥∆(Lτ)∥ = maxX :∥X∥F=1 ∥∆(Lτ)X∥F ≤ Cκϵτ , which completes the proof.

heading: D PROOF OF THEOREM 1
text: The aim of the proof is to bound the distance between the space-time graph filters before and after adding joint perturbations. From (15), this distance can be evaluated as
∥H(S,Lτ )−H(Ŝ, L̂τ )∥P,T = min P∈P min s∈R max X :∥X∥F=1 ∥∥∥H(S,Lτ )X −H(PT ŜP, e−sLτ L̂τ)X∥∥∥ F
(a) ≤ max X :∥X∥F=1 ∥∥∥H(S,Lτ )X −H(S+ SE+ES, L̂τ )X∥∥∥ F
= ∥∥∥H(S,Lτ )−H(S+ SE+ES, L̂τ )∥∥∥ ,
(74)
where (a) is true for any specified values of P and s, and we choose P0 and 0 as in (31) and (61), respectively. Then, we add and subtract H(S+ SE+ES,Lτ ) from the filter difference to get
H(S,Lτ )−H(S+ SE+ES, L̂τ ) = H(S,Lτ )−H(S+ SE+ES,Lτ ) +H(S+ SE+ES,Lτ )−H(S+ SE+ES, L̂τ ). (75)
Using the triangular inequality, the norm of the filter difference is bounded by
∥H(S,Lτ )−H(S+ SE+ES, L̂τ )∥ ≤ ∥H(S,Lτ )−H(S+ SE+ES,Lτ )∥ + ∥H(S+ SE+ES,Lτ )−H(S+ SE+ES, L̂τ )∥ = ∥∆s∥+ ∥∆τ∥.
(76)
Note that in Proposition 2, ∆τ is defined for the GSO S but the proposition is valid for any GSO. The bounds of ∥∆s∥ and ∥∆τ∥ are obtained in Propositions 1 and 2. Therefore, we get
∥H(S,Lτ )−H(S+ SE+ES, L̂τ )∥ ≤ 2Cϵs(1 + δ √ N) + Cκϵτ +O(ϵ2), (77)
where O(ϵ2) is the highest among O(ϵ2s) and O(ϵ2τ ). This completes the proof.
D.1 FINITE-IMPULSE RESPONSE FILTERS
While Theorem 1 handles space-time graph filters with continuous-time impulse response h(t), its results can be extended to finite-impulse response filters following the same proofs in Sections B, C, and D. In this context, FIR filters can be thought of as a special case of the filter h(t), where we have a finite number of filter coefficients. We summarize this remark in the following lemma. Lemma 2. Define a space-time graph filter with a finite impulse response as
Hd(S,Lτ ) = K∑
k=0
hke −kTsS◦Lτ . (78)
Under the GSOs defined in Proposition 1 and TSOs in Proposition 2, the distance between the FIR space-time graph filters Hd(S,Lτ ) and Hd(Ŝ, L̂τ ) satisfies
∥Hd(S,Lτ )−Hd(Ŝ, L̂τ )∥P,T ≤ 2Cϵs ( 1 + δ √ N ) + Cκϵτ +O(ϵ2). (79)

heading: E PROOF OF THEOREM 2
text: The goal of the proof is to show the difference between the GNN output before and after adding perturbations. The GNN output is the L layer’s output, i.e., XL. From (7), we can write the difference between the two GNNs as
∥Φ(X0,H,S ◦ Lτ )−Φ(X0,H, Ŝ ◦ L̂τ )∥F = ∥XL − X̂L∥F = ∥∥∥σ(HL(S,Lτ )XL−1)− σ (HL(Ŝ, L̂τ )X̂L−1)∥∥∥
F ,
(80)
where HL(S,Lτ ) is the filter in (78) at Layer L, Xℓ is the output of layer ℓ, and X̂ℓ is the corresponding output after adding the perturbations.
From assumption (A2), we get ∥σ(X2) − σ(X1)∥2 ≤ ∥X2 − X1∥2. Accordingly, the norm of the output difference at any layer ℓ ≤ L is bounded by
∥Xℓ − X̂ℓ∥F ≤ ∥∥∥Hℓ(S,Lτ )Xℓ−1 −Hℓ(Ŝ, L̂τ )X̂ℓ−1∥∥∥
F . (81)
Add and subtract Hℓ(Ŝ, L̂τ )Xℓ−1 inside the norm in the right hand side. The norm of (81) can then be wriien as
∥Xℓ − X̂ℓ∥F ≤ ∥∥∥Hℓ(S,Lτ )Xℓ−1 −Hℓ(Ŝ, L̂τ )Xℓ−1∥∥∥+ ∥∥∥Hℓ(Ŝ, L̂τ )Xℓ−1 −Hℓ(Ŝ, L̂τ )X̂ℓ−1∥∥∥
F ≤ ∥∥∥Hℓ(S,Lτ )−Hℓ(Ŝ, L̂τ )∥∥∥ ∥Xℓ−1∥F + ∥∥∥Hℓ(Ŝ, L̂τ )∥∥∥∥∥∥Xℓ−1 − X̂ℓ−1∥∥∥
F ,
(82) using the triangle inequality in (a) and Cauchy-Schwarz inequality in (b). From assumption (A1), the filters have unit operator norm, i.e., ∥Hℓ(S,Lτ )∥ = 1,∀ℓ = 1, . . . , L, and hence, ∥Xℓ∥F ≤ ∥Xℓ−1∥F ≤ ∥X0∥F . From the stability of graph filters in Theorem 1 (with Lemma 2), the difference becomes
∥XL − X̂L∥F ≤ ( 2Cϵs ( 1 + δ √ N ) + Cκϵτ +O(ϵ2) ) ∥X0∥F + ∥∥∥X̂ℓ−1 −Xℓ−1∥∥∥ F . (83)
Substituting (82) in (83) recursively, we obtain ∥XL − X̂L∥F ≤ 2CLϵ ( 1 + δ √ N ) ∥X0∥F + CLκϵτ∥X0∥F +O(ϵ2). (84)
Note that the base case of the recursion is calculated as ∥X1 − X̂1∥F ≤ ∥∥∥H1(S,Lτ )X0 −H1(Ŝ, L̂τ )X0∥∥∥
F ≤ ∥∥∥H1(S,Lτ )−H1(Ŝ, L̂τ )∥∥∥ ∥X0∥F
≤ ( 2Cϵs ( 1 + δ √ N ) + Cκϵτ +O(ϵ2) ) ∥X0∥F
(85)
Eventually, as in (74), it follows that the joint operator distance modulo of the output of ST-GNNs is ∥Φ(.;H,S ◦ Lτ )−Φ(.;H, Ŝ ◦ L̂τ )∥P,T ≤ 2CLϵs ( 1 + δ √ N ) + CLκϵτ +O(ϵ2), (86)
completing the proof.

heading: F STABILITY OF MULTIPLE-INPUT MULTIPLE-OUTPUT ST-GNNS
text: Theorem 2 only considers the case of single-feature layers. In this section, we extend the stability analysis to the case of multiple features. We derive a bound for the difference in the output of ST-GNNs in Corollary 1, where each layer as well as the input signals have multiple features. We let the features at the hidden layers to be the same for simplicity. Corollary 1 (ST-GNNs Stability). Consider the assumptions of Theorem 2 but let Fℓ = F be the number of features per each layer for 1 ≤ ℓ ≤ L− 1. Let F0 and FL be the number of the features of the input and output signals, respectively. Then,
∥Φ(.;H,S ◦ Lτ )−Φ(.;H, Ŝ ◦ L̂τ )∥P,T ≤√ FL ( FL−1F0 +
L−1∑ l=1 F l
)( 2Cϵs(1 + δ √ N) + Cκϵτ ) +O(ϵ2).
(87)
Proof. We aim to find the difference between the output of the Lth layer before and after adding joint perturbations. For a layer ℓ and feature f , we can define its output as
X fℓ = σ ( F∑ g=1 Hfgℓ (S,Lτ )X g ℓ−1 ) , (88)
where Hfgℓ is a single-input signle-output filter as the one used in Theorem 2. From assumption (A2) of Theorem 2, we have the identity ∥σ(x2)− σ(x1)∥2 ≤ ∥x2 − x1∥2. We can then bound the output difference for all 1 < ℓ < L:
∥X fℓ − X̂ f ℓ ∥F
(a) ≤ F∑
g=1
∥∥∥Hfgℓ (S,Lτ )X gℓ−1 −Hfgℓ (Ŝ, L̂τ )X̂ gℓ−1∥∥∥ F
(b) ≤ F∑
g=1
(∥∥∥Hfgℓ (S,Lτ )−Hfgℓ (Ŝ, L̂τ )∥∥∥∥∥X gℓ−1∥∥F + ∥∥∥Hfgℓ (Ŝ, L̂τ )∥∥∥∥∥∥X̂ gℓ−1 −X gℓ−1∥∥∥F ) ,
(c) ≤ F ( 2Cϵ ( 1 + δ √ N ) + Cκϵτ +O(ϵ2) )∥∥X 1ℓ−1∥∥F + F ∥∥∥X̂ 1ℓ−1 −X 1ℓ−1∥∥∥F , (89)
where (a) follows from the triangle inequality, (b) from (82), and (c) from Theorem 1. The summation in (b) has equivalent F summands since we assume that all the filters satisfies assumption (A1) in Theorem 2. With applying ∥∥∥X̂ 1ℓ−1 −X 1ℓ−1∥∥∥ F recursively, we can re-write (89) as
∥X fℓ − X̂ f ℓ ∥F ≤ ( F ℓ−1F0 + ℓ−1∑ l=1 F l )( 2Cϵ ( 1 + δ √ N + Cκϵτ ) +O(ϵ2) ) . (90)
Note that the base case of the recursion is calculated as ∥X f1 − X̂ f1 ∥F ≤ F0∑ g=1 ∥∥∥Hfg1 (S,Lτ )X g0 −Hfg1 (Ŝ, L̂τ )X g0 ∥∥∥ F
≤ F0 ( 2Cϵs ( 1 + δ √ N ) + Cκϵτ +O(ϵ2) ) ∥X 10 ∥F .
(91)
Finally, we can express the difference between the output at layer L. It has multiple features, and therefore, the norm of the difference is expressed as∥∥∥XL − X̂L∥∥∥2
F = ∥∥∥∥∥ FL∑ g=1 ( X fL − X̂ fL )∥∥∥∥∥ 2
F
≤ FL ∥∥∥X 1L − X̂ 1L∥∥∥2
F
(a) ≤ FL ( FL−1F0 +
L−1∑ l=1 F l
)2 ( 2Cϵs ( 1 + δ √ N ) + Cκϵτ +O(ϵ2) )2 , (92)
where (a) is calculated from (90) for ℓ = L. Taking the square root yields the inequality in (87), which completes the proof.

heading: G EXTENDED NUMERICAL RESULTS
text: We consider the problem of decentralized controllers, where we are given a team of N agents, each of which has a position pi,n ∈ R2, a velocity vi,n ∈ R2 and an acceleration ui,n ∈ R2, that are captured at times nTs, n ∈ Z+. The goal is to learn controller actions that allow the agents to move together and complete a specific task. Two tasks are considered in this paper: flocking and unlabeled motion planing. Optimal centralized controllers for the two tasks are derived in the literature. However, centralized controllers require access to the information at all the agents, and therefore, the computation complexity scales fast with the number of agents. With the help of ST-GNNs, we can find decentralized controllers that imitate the centralized solutions according to (22). In this section, we aim to provide a detailed description of the experiments in Section 5 along with further experiments.
Communication networks. The underlying graphs represent the communication networks between the agents. Two criteria can be used to assemble the graph. First, each agent is connected to its M -nearest neighbors. The graph at a time step n is then represented by a binary graph Gn such that (i, j) ∈ En if and only if j ∈ Ni,n or i ∈ Nj,n, where Ni,n is the set of the M -nearest neighbors of the agent i. The second is to connect each agent to the neighbors within a communication range R. The graph in this case is also represented by a binary graph such that (i, j) ∈ En if and only if ∥pi,n − pj,n∥ < R. In both cases, when the agents move, the graph Gn changes with n. ST-GNN Implementation. As indicated above, the communication networks (i.e., graphs) change with the movement of the agents (i.e., with time). However, ST-GNNs in (7) are designed for fixed graphs. Moreover, ST-GNNs are designed for continuous-time signals, whereas the signals in our experiments are discrete. To turn around these challenges, we implement the FIR space-time graph filter recursively, where at every time step we use the corresponding underlying graph. The output of the filter at a time step n is expressed as
yn = h0xn + K−1∑ k=1 hk
( k∏
m=1
Sn−m ) xn−k, (93)
where xn and Sn are the graph signal and the GSO at time step n. The graph diffusion is implemented in (93) with the GSO directly instead of the operator exponential e−Sn . The latter can be expressed as
e−Sn = ∞∑ l=0 (−1)l l! Sln = I− Sn +O(Sn), (94)
where O(Sn) contains the higher-order terms. Equation (94) shows that the GSO is a first-order approximation of e−Sn , and the sign difference is absorbed in the learning parameters {hk}Kk=0.
G.1 APPLICATION I: FLOCKING AND NETWORK CONSENSUS
In this application, the agents collaborate to avoid collisions and learn to move according to a reference velocity rn ∈ R2. The reference velocity is generated randomly as rn+1 = rn + Ts∆rn, ∀n ∈ Z+, n < T . The initial value r0 is sampled from a Gaussian distribution and so is ∆rn, with zero mean and expected norms E[∥r0∥] and E[∥r0∥], respectively. At each time step n, each agent i observes a biased reference velocity r̃i,n such that r̃i,n = rn +∆r̃i with ∆r̃i being white Gaussian with independent and identically-distributed (i.i.d.) components. The goal is to learn acceleration actions {ui,n}i,n that allow the agents to form a swarm moving with the same velocity. Mobility model. Given acceleration actions {ui,n}i,n, each agent moves according to the equations of motion:
vi,n+1 = vi,n + Tsui,n, pi,n+1 = pi,n + Tsvi,n + T 2 s/2 ui,n. (95)
The initial velocity is assumed to be vi,0 = r0 +∆v, for all i, where ∆v is white Gaussian with i.i.d. components. In our experiments, we aim to learn accelerations {ui,n}i,n that make vi,n for all i be as close as possible to rn and prevent the position differences pij,n = pi,n − pj,n from being zero for all i ̸= j at any time step n. We solve this problem under the constraints ∥ui,n∥2 ≤ µ,∀i, n.
Optimal centralized controllers. The problem described above has an optimal solution, which, for all i = 1, . . . , N and n ≤ T , is given as
u∗i,n = −1 2Ts vi,n − 1 N N∑ j=1 r̃j,n − 1 2Ts N∑ j=1 ∇pi,nC(pi,n − pj,n), (96)
if ∥u∗i,n∥2 ≤ µ, and otherwise u∗i,n = µ. The collision avoidance potential C(.) is defined as C(pi,pj) = {
1/∥pij∥22 − log(∥pij∥22), ∥pij∥2 ≤ γ, 1/γ2 − log(γ2), otherwise, (97)
where pij = pi − pj , and γ = 1 (Tanner et al., 2003; Gama et al., 2020b). Note that (96) requires access to data from all the agents, and therefore, this solution is known to be the optimal centralized solution. However, in a decentralized setting, the agents exchange information only within their K-neighborhood only. Therefore, our goal is to find acceleration controls that imitate the centralized solution in (96).
Decentralized Controllers. We compare our results to a decentralized controller that only has access to the information shared within the K-neighborhood of the agents. Unlike the centralized policy where all the information is available to the central unit, the agents only have access to the data that they receive from their neighbors. More importantly, the data get delayed through the communication network, and the agents have to make their predictions based on outdated data. Therefore, the estimated accelerations in (96) are approximated by
udeci,n = −1 2Ts vi,n − K∑ k=0
1 K|N ki,n| ∑
j∈Nki,n
r̃j,(n−k) − 1 2Ts K∑ k=1 ∑ j∈Nki,n ∇pi,nC ( pi,n − pj,(n−k) ) ,
(98) where |.| is the set cardinality, and N ki,n = {j′ ∈ N k−1j,(n−1) | j ∈ Ni,n} is the set of the k-hop neighbors.
G.1.1 EXPERIMENT 1: NETWORK CONSENSUS
In this section, we present a detailed description of the training process. Remember that the agents in this experiment are not moving and the graphs are fixed mesh grids. The state variables are the estimated and observed velocities, vi,n and r̃i,n,∀i, n. Training. The dataset is generated according to the mobility model in (95) and (96). The dataset consists of 500 time-varying graph signals {Xm}500m=1 that are calculated under optimal centralized policies {U∗m}500m=1. We split the data into 460 examples for training, 20 for validation and 20 for testing. We train a 2-layer ST-GNN on the training data and optimize the mean squared loss using ADAM algorithm with learning rate 0.01 and decaying factors β1 = 0.9 and β2 = 0.999.1 We then keep the model with the lowest cost (across the validation data) among 30 epochs while the cost is averaged over the T steps. For a single time step, the cost is calculated as
c(un) = 1
2N N∑ i=1 ∥∥∥∥∥∥vi,n − 1N N∑ j=1 r̃j,n ∥∥∥∥∥∥ 2
2
+ 1
2N N∑ i=1 ∥Tsui,n∥22, ∀n < T. (99)
All the training and simulation parameters are shown in Table 1.
Execution. In addition to the 20 test examples generated with the training set, we generate another 20 examples under joint graph and time perturbations for each value of ϵ. The perturbation size of both graph and time topologies is chosen the same. The results shown in Fig. 1 (Middle) suggest that the closer the space-time structures, the closer the outputs of the trained ST-GNN under the two cases are. This shows that an ST-GNN, which is trained with signals defined on one particular graph and sampled at a certain sampling rate, can be generalized to signals with a different underlying topology as long as the two topologies are close.
1We used the GNN library at https://github.com/alelab-upenn/ graph-neural-networks
G.1.2 EXPERIMENT 2: FLOCKING
In this experiment, we have 50 agents that are spread uniformly with initial density ρ0. The agents are allowed to move with their velocities {vi,n}i,n. At each time step n, the graph Gn represents their communication network based on a communication range R. We train a 2-layer ST-GNN that is implemented in (93), and the training procedure is exactly as in the previous experiment. The dataset consists of 800 training, 100 validation, and 100 test examples. We execute the trained ST-GNNs on signals defined over graphs constructed with the same initial agent density ρ0 and sampled at the same sampling rate 1/Ts. Fig. 3 depicts three paradigms of the learned velocities compared to the reference velocity. We notice that the agents follow the trend in the reference velocity. The mean and variance of the difference between the agent and reference velocities were shown before in Fig. 1 (Right). Fig. 3d illustrates the positions of the swarm after 9s and shows that the swarm moved together with same velocity in the same direction. Fig. 1 (Right) also shows that the proposed architecture outperforms the decentralized policy in (98). The difference between the estimated and reference velocities is lower under ST-GNNs than the difference under the decentralized policy.
Examples of the test dataset are shown in Fig. 4. It is clear that ST-GNNs help the agents to mitigate the delays in the received information and learn accelerations that follow the reference velocity.
Graph and time perturbations appear in real applications as a change in the underlying graph-time topology. In the following, we provide extended experiments to show the effect of perturbations. In particular, we repeat the same experiment either with graphs constructed with different agent densities or under different values of sampling time Ts.
Experiment 2.I. In this experiment, we execute the trained ST-GNN on graphs constructed under different initial agent densities (i.e., a source of graph perturbations). We execute the trained STGNNs on 50 signals generated with a different initial agent density. We repeat the experiments under densities 2, 12 , 1 8 , 1 32 , 1 128 , 1 512 agents/m
2, respectively, and plot the relative costs after 10s in Figure 5 (Left). The cost represents the mean difference between the agent velocity and the mean of the observed velocities, which is calculated as
cost = 1
2N N∑ i=1 ∥∥∥∥∥∥vi,T − 1N N∑ j=1 r̃j,T ∥∥∥∥∥∥ 2
2
. (100)
The relative cost is then calculated as a relative deviation from the cost under a density of 2 agents/m2. Figure 5 (Left) shows that for the small changes in the density ∆ρ (e.g., from the original density 2 to 0.5 agents/m2), the relative cost remains small. However, the higher values of ∆ρ result in higher relative costs, i.e., degradation in the performance. Note that using smaller agent densities (i.e., higher ∆ρ) at the same communication range R results in sparser graphs. Therefore, their distances to the graphs generated at the original density (2 agents/s) increases leading to the aforementioned performance degradation.
Experiment 2.II. In this experiment, we execute the trained ST-GNN on signals sampled at different sampling rates (i.e., a source of time perturbations). For several values of ∆Ts, we replace the sampling time Ts with Ts +∆Ts and calculate the relative costs as in Experiment 2. Figure 5 (Right) shows that the smaller the value of ∆Ts, the smaller the relative cost is. This matches our theoretical results.
G.2 APPLICATION II: UNLABELED MOTION PLANNING
In this problem, we aim to assign N unlabeled agents to N target goals, {gj ∈ R2}Nj=1 through planning their free-collision trajectories. The term unlabeled implies that the assignment is not pre-determined and it is executed online. At the start of the experiment, the agents are spread with minimum inter-agent distance d and so do the goal targets. A centralized solution to this problem was introduced in (Turpin et al., 2014), which we refer to as the CAPT solution. This solution gives the agent trajectories, {Qi ∈ R2×T }Ni=1, while the agent velocities and accelerations can be calculated as vi,n = (pi,n+1 − pi,n)/Ts, u∗i,n = (pi,n+1 − pi,n)/T 2s . (101) Note that pi,n is the n-th column of matrix Qi. We use the optimal accelerations in (101) to learn an ST-GNN parameterization that predicts the agent accelerations {ui,n}i,n according to (22). The input Xm consists of 6M + 4 state variables for each agent, which are the agent position {pi,n}n and velocity {vi,n}n, the position of the nearest M neighbors {Pi,n ∈ RM×2 | [Pi,n]j = pj,n ∀j ∈ Ni,n}n along with their velocities, and the position of the nearest M target goals. The CAPT accelerations and the corresponding state variables constitutes together one pair in the dataset.
The dataset consists of 55000 training examples and 125 validation examples. We train an ST-GNN on the training data and optimize the mean squared loss using ADAM algorithm with learning rate 0.0005 and decaying factors β1 = 0.9 and β2 = 0.999. We then keep the model with the lowest mean distance between the agent final position and its desired target among 60 epochs. The other training parameters are shown in Table 2. The trained ST-GNN is then executed on a test dataset that consists of 1000 examples. The ST-GNN output for each example is the estimated accelerations, and the planned trajectories are then calculated using (95). One example was shown in Fig. 2 (Left), which depicts free-collision trajectories.
Similar to the flocking experiment, we test the trained ST-GNNs with different graph-time topologies, generated with either different neighborhood sizes or different sampling times. Figures 2 (Middle) and (Right) uphold the conviction that the ST-GNN output difference increases with the distance between the underlying topologies.

Review summaries:
This paper proposes a new deep learning architecture called ST-GNN that learns representations on graphs that evolve over time. Their work focuses on developing an interesting time-varying convolutional architecture, which exploits the graph-time underlying structure of the signals, processing across both the graph and time domains. The authors introduce new graph shift operators (GSO), time shift operators (TSO), followed by a linear combination of those which results in a space-time shift operator termed STSO. The advantage of this paper is that it can handle continuous-time graph signals and the stability of this proposed architecture is studied theoretically (Section 4) with a significant result outlined in Proposition 1 of the paper where the difference between the space-time graph filters of the original graph and the perturbed graph are upper bounded by the order of error introduced.    

This paper introduces a new spatio-temporal Graph Neural Network, ST-GNN, for making predictions on temporal network. Its proposed space-time convolution operator is a composition of temporal convolution and graph diffusion. The paper further proves that under practical conditions their ST-GNN with Integral Lipschitz filters is stale to small perturbation in both time and graph domain.