TAMP-S2GCNETS: COUPLING TIME-AWARE MUL-
Year: 2021

Abstract:
Graph Neural Networks (GNNs) are proven to be a powerful machinery for learning complex dependencies in multivariate spatio-temporal processes. However, most existing GNNs have inherently static architectures, and as a result, do not explicitly account for time dependencies of the encoded knowledge and are limited in their ability to simultaneously infer latent time-conditioned relations among entities. We postulate that such hidden time-conditioned properties may be captured by the tools of multipersistence, i.e, a emerging machinery in topological data analysis which allows us to quantify dynamics of the data shape along multiple geometric dimensions. We make the first step toward integrating the two rising research directions, that is, time-aware deep learning and multipersistence, and propose a new model, Time-Aware Multipersistence Spatio-Supra Graph Convolutional Network (TAMP-S2GCNets). We summarize inherent time-conditioned topological properties of the data as time-aware multipersistence Euler-Poincaré surface and prove its stability. We then construct a supragraph convolution module which simultaneously accounts for the extracted intraand inter-dependencies in the data. Our extensive experiments on highway traffic flow, Ethereum token prices, and COVID-19 hospitalizations demonstrate that TAMP-S2GCNets outperforms the state-of-the-art tools in multivariate time series forecasting tasks.

text: Graph Neural Networks (GNNs) are proven to be a powerful machinery for learning complex dependencies in multivariate spatio-temporal processes. However, most existing GNNs have inherently static architectures, and as a result, do not explicitly account for time dependencies of the encoded knowledge and are limited in their ability to simultaneously infer latent time-conditioned relations among entities. We postulate that such hidden time-conditioned properties may be captured by the tools of multipersistence, i.e, a emerging machinery in topological data analysis which allows us to quantify dynamics of the data shape along multiple geometric dimensions. We make the first step toward integrating the two rising research directions, that is, time-aware deep learning and multipersistence, and propose a new model, Time-Aware Multipersistence Spatio-Supra Graph Convolutional Network (TAMP-S2GCNets). We summarize inherent time-conditioned topological properties of the data as time-aware multipersistence Euler-Poincaré surface and prove its stability. We then construct a supragraph convolution module which simultaneously accounts for the extracted intra- and inter-dependencies in the data. Our extensive experiments on highway traffic flow, Ethereum token prices, and COVID-19 hospitalizations demonstrate that TAMP-S2GCNets outperforms the state-of-the-art tools in multivariate time series forecasting tasks.

heading: 1 INTRODUCTION
text: Multivariate time series forecasting plays an integral role in virtually every aspect of societal functioning, from biosurveillance to financial analytics to intelligent transportation solutions. In the last few years, Graph Convolutional Networks (GCNs) have emerged as a powerful alternative to more conventional time series predictive models. Despite their proven success, GCNs tend to be limited in their ability to simultaneously infer latent temporal relations among entities (such as correlations both within a time series and in-between time series or the joint spatio-temporal dependencies). More generally, most existing GCNs architectures are inherently static and as such, do not explicitly account for time-conditioned properties of the encoded knowledge about the complex dynamic phenomena.
At the same time, recent studies on integrating shape properties of the complex data into deep learning (DL) models indicate that topological representations, obtained using the tools of single parameter persistence, can bring an invaluable insight into the system organization and enhance the resulting graph learning mechanisms (Hofer et al., 2019; Carrière et al., 2020; Carlsson & Gabrielsson, 2020; Horn et al., 2021). (By shape here we broadly understand data properties which are invariant under continuous deformations, e.g., stretching, bending, and twisting). However, in many applications, particularly, involving spatio-temporal processes, the data exhibit richer structures which cannot be well encoded with a single parameter persistence. We postulate that many critical hidden time-conditioned interrelations which are inaccessible with other methods can be captured by the emerging machinery of multipersistence. Multipersistence, or multiparameter persistence (MP) generalizes the notion of single parameter persistence to a case when dynamics of the inherent data shape is discerned along multiple geometric dimensions (Carlsson & Zomorodian, 2009). Despite
its premise, applications of MP in any discipline remain nascent at best (Riess & Hansen, 2020; Kerber, 2020).
We make the first step on a path of bridging the two emerging directions, namely, time-aware DL with time-conditioned MP representations of complex dynamic phenomena. By time-conditioned MP representations, we mean the most salient topological properties of the data that manifest themselves over time. To summarize such time-conditioned topological properties, we first introduce a dynamic Euler-Poincaré surface as a new MP summary. We then propose a directed multilayer supra graph abstraction to represent a sequence of time-varying objects and develop a supragraph convolution module which allows us to simultaneously learn co-evolving intra- and inter-dependencies (i.e., spatial and temporal correlations) in the complex high-dimensional data.
The key novelty of this paper are summarized as follows:
• This is the first work to bridge the concepts of MP with the time-aware learning paradigm. Applications of MP in any field of study are currently nascent.
• We introduce a new time-aware multipersistence invariant, a dynamic Euler-Poincaré surface. We prove its stability and show its substantial computational gains and high utility for encoding the time-conditioned knowledge.
• We propose a mathematical abstraction of directed multilayer supra graph for time-conditioned knowledge representation and construct a new Time-Aware Multipersistence Spatio-Supra Graph Convolutional Network (TAMP-S2GCNets) which simultaneously learns latent temporal interand intra-relations among entities in the complex high-dimensional data.
• We perform expansive forecasting experiments, in application to highway traffic flow, Ethereum token prices, and COVID-19 hospitalizations. Our findings demonstrate superior predictive performance, versatility and computational efficiency of TAMP-S2GCNets, compared to the stateof-the-art methods in multivariate time series forecasting.

heading: 2 RELATED WORKS
text: Multipersistence Despite that MP demonstrates very promising results in terms of improving accuracy, tractability and robustness, applications of MP in ML are virtually non-existent (Wright & Zheng, 2020; Riess & Hansen, 2020; Kerber, 2020). Some notable efforts in the direction to develop MP summaries which are suitable for integration with ML models include Multiparameter Persistence Kernel of Corbet et al. (2019), Multiparameter Persistence Landscapes (MP-L) of Vipond (2020), and Multiparameter Persistence Images (MP-I) of Carrière & Blumberg (2020) which are based on the concept of slicing, that is, restricting the MP module to an affine line (or single parameter persistence) (Cerri et al., 2013; Landi, 2014). Such slicing methods enjoy a number of important stability guarantees but tend to be computationally expensive even in static scenarios, which makes them infeasible for time series forecasting tasks. In turn, the most recent results of Beltramo et al. (2021) (i.e., Euler Characteristic Surfaces) and Coskunuzer et al. (2021) (i.e., Multiparameter Persistence Grids) introduce pointwise MP representations, in application to static point clouds and graphs. Such pointwise representations are weaker invariants but are substantially more computationally efficient. Integration of pointwise representation with ML models has not been yet investigated. Here we propose the first time-aware pointwise MP representation, a dynamic Euler-Poincaré surface, derive its theoretical properties and integrate it with GCN in time series forecasting tasks.
Spatial-Temporal Graph Models and Forecasting Recent studies (Li et al., 2018; Yu et al., 2018; Yao et al., 2018) introduce graph convolution methods into spatio-temporal networks for multivariate time series forecasting which, as a result, allows for better modeling of dependencies among entities (Wu et al., 2019; Bai et al., 2020; Cao et al., 2020) and handling data heterogeneity. Despite the GCN successes, designs of the existing spatial-temporal GCNs largely rely on the pre-defined graph structures. As such, GCNs are restricted in their ability to explicitly integrate time dimension into the knowledge representation and learning mechanisms, thereby, limiting model adaptivity to the dynamic environments and requiring more frequent retraining. Most recently, Chen et al. (2021) propose a time-aware GCN, Z-GCNETs, for time series forecasting which integrates zigzag persistence images based on a single filtration, as the primary time-conditioned topological representation. In general, the zigzag concept can be combined with MP, but it requires more fundamental advances in the theory of algebraic topology. As such, our time-aware MP learning approach may
be viewed as complementary to zigzag persistence, while considering time-changing connections of graph structures in dynamic networks.

heading: 3 TIME-AWARE MULTIPERSISTENCE EULER CHARACTERISTIC SURFACES
text: Spatio-Temporal Graph Construction We define a spatial network at time step t as Gt = (Vt, Et, At,Xt), where Vt is a set of nodes and Et is a set of edges. We let |Vt| = N and |Et| = Mt. The adjacency matrix At ∈ RN×N , and Xt = {xt,1, xt,2, . . . , xt,N}> ∈ RN×FN is the node feature matrix with feature dimension FN . To construct the spatial network Gt, we can build the adjacency matrix At based on (i) the prior knowledge of graph structure: first-order neighbours, i.e., At,uv = 1 if the node u and node v have a connection in the dynamic graph at the time step t; and (ii) the Radial Basis Function (RBF): degrees of similarity between instances (i.e., nodes) in Xt, i.e., At,uv = 1exp (−||xt,u−xt,v||2/γ)≤ , where γ denotes the length scale parameter and denotes the threshold parameter filters noisy edges. Let T be the total number of time steps. Given a sequence of observations on a multivariate variable, X = {X1,X2, . . . ,XT } ∈ RN×FN×T with T timestamps and FN node attributes, we construct spatio-temporal networks G = {G1,G2, . . . ,GT } via either prior knowledge of network structure or applying RBF to the node feature matrix.
Single Parameter Persistence Persistent homology (PH) based on one parameter discerns shape of the complex data along a single geometric dimension. The goal is to select some suitable parameter of interest and then to study a graph Gt not as a single object, but as a sequence of nested subgraphs, or a graph filtration G1t ⊆ G2t . . . ⊆ Gmt = Gt, induced by this evolving scale parameter. Armed with such filtration, we can then assess which structural patterns (e.g., loops and cavities) appear/disappear and record their lifespans. To make the counting process more efficient and systematic, we build a simplicial complex Cit from each subgraph Git , resulting in a filtration C1t ⊆ C2t . . . ⊆ Cmt (e.g., clique complexes). For example, to construct such filtration, a common method is to consider a filtering function f : Vt 7→ R and the corresponding increasing set of thresholds {αi}m1 such that Cit = {∆ ∈ Ct : maxv∈∆ f(v) ≤ αi}. The resulting construction is called a sublevel set filtration of f , and f can be selected, for instance, as degree, centrality, or eccentricity function (Hofer et al., 2020; Cai & Wang, 2020). Similarly, f can be defined on the set of edges Et. More details on a single filtration PH can be found in Appendix B.
Time-Aware Multiparameter Persistence (TAMP) Data in many applications, particularly, involving spatio-temporal modeling, might be naturally indexed by multiple parameters, e.g., real time traffic flow and optimal route in urban transportation analytics. Alternatively, the primary focus might be on discerning shape properties of the complex data along multiple dimensions. For instance, to better predict cryptocurrency prices and manage cryptomarket investment performance, we may need to evaluate structural patterns in cryptoasset dynamics not along one dimension but simultaneously along the volume of transactions and transaction graph betweenness, as a measure of joint perception of the cryptomarket volatility among the key investors. Such multidimensional analysis of topological and geometric properties can be addressed using generalization of PH based on a single filtration to a multifiltration case. That is, the MP idea is to simultaneously assess shape characteristics of Gt based on a multivariate filtering function F : Vt 7→ Rd. As a result, e.g., for d = 2 and a set of nondecreasing thresholds {αi}m1 and {βj}n1 , instead of a single filtration of complexes, we get a bifiltration of complexes {Cαi,βjt | 1 ≤ i ≤ m, 1 ≤ j ≤ n} such that if βk < βl, then C αi,βk t ↪→ C αi,βl t and if αi < αj , then C αi,βk t ↪→ C αj ,βk t . Finally, this bifiltration induces a bigraded MP module {Hk(Cαi,βkt )}, where Hk is the kth homology group. Inspired by Beltramo et al. (2021); Coskunuzer et al. (2021), we propose a new time-aware MP summary, namely, a Dynamic Euler-Poincaré Surface. Definition 3.1 (Dynamic Euler-Poincaré Surface). Let {Gt}Tt=1 be a series of time-varying graphs. Let F = (f, g) be a multivariate filtering function F : Vt 7→ R2 with thresholds I = {(αi, βj) | 1 ≤ i ≤ m, 1 ≤ j ≤ n}. Let Cαi,βjt be the clique complex of the induced subgraph G αi,βj t = F−1((−∞, αi] × (−∞, βj ]), t = 1, 2, . . . , T and χ be the Euler–Poincaré characteristic. Then, a sequence of time-evolving m × n-matrices {Et}Tt=1 such that Etij = χ(C αi,βj t ) for 1 ≤ i ≤ m, 1 ≤ j ≤ n is called Dynamic Euler-Poincaré Surface (DEPS). (Figure 2 in Appendix C shows a toy example how DEPS is computed.)
Theoretical Guarantees of DEPS Consider two graphs G+ and G−, where time index t is suppressed for brevity. Let F : V± 7→ R2 be a multivariate filtering function with thresholds
I = {(αi, βj) | 1 ≤ i ≤ m, 1 ≤ j ≤ n}. Let C± be the clique complexes of G±, and let Ĉ± = {C+ij} be the bifiltration induced by (C±, F, I) as before. Let E± be the corresponding Euler-Poincaré Surfaces (i.e., m×n matrices). Then, set ‖E+−E−‖1,1 = ∑m i=1 ∑n j=1 |E + ij −E − ij | as the distance between E+ and E−, where ‖ · ‖1,1 is the vectorized L1 matrix norm. We now introduce an L1-based MP metric instead of using L∞-based metrics like, e.g., matching or interleaving, due to the nature of our summaries E± (see Remark C.1 in Appendix C). Let Dfk (C±) and Dgk(C±) be the kth single parameter persistence diagrams (PDs) of C± for filtrations induced by functions f, g : V± 7→ R, respectively (see Appendix B). Let C±i∗ and C ± ∗j be clique complexes corresponding to G±i∗ = f−1((−∞, αi]) and G ± ∗j = g
−1((−∞, βj ]). Define the ith column distance for the kth PDs as Dki∗(Ĉ +, Ĉ−) =W1(Dgk(C + i∗),D g k(C − i∗)), whereW1 is the Wasserstein-1 distance. Similarly, the jth row distance for kth PDs is Dk∗j(Ĉ +, Ĉ−) =W1(Dfk (C + ∗j),D f k (C − ∗j)).
Definition 3.2 (Weak L1-metric for Multipersistence). The weak L1-metric between Ĉ± is
D(Ĉ+, Ĉ−) = max { Dc(Ĉ +, Ĉ−),Dr(Ĉ +, Ĉ−) } ,
s.t. Dc(Ĉ+, Ĉ−) = ∑M k=0 ∑m i=1 D k i∗(Ĉ +, Ĉ−) and Dr(Ĉ+, Ĉ−) = ∑M k=0 ∑n j=1 D k ∗j(Ĉ +, Ĉ−).
Now, we can state our stability result for Euler-Poincaré Surfaces.
Theorem 3.1. Let G±, F, Ĉ±,E± be as defined above. Then, the Euler-Poincaré Surfaces are stable with respect to the weak L1-metric, i.e., ‖E+ − E−‖1,1 ≤ C ·D(Ĉ+, Ĉ−) for some C > 0.
The proof of the theorem is given in Appendix C. This stability result implies that the distances between multiparameter PDs control the distance between the resulting Euler-Poincaré Surfaces. By combining with the stability result for PDs (Cohen-Steiner et al., 2007), one can conclude that the small changes in the MP filtering function F : Vt 7→ Rd or the small changes in the input data can result only in a small change in DEPS surfaces. For further discussion on implications of the stability result, please refer to Remark C.2.

heading: 4 TIME-AWARE MULTIPERSISTENCE SPATIO-SUPRA GRAPH CONVOLUTIONAL NETWORKS
text: Given T historical observations X T = {Xt−T ,Xt−T+1, . . . ,Xt−1}, the multivariate forecasting model F(·) is learned to predict future observations in the next H + 1 time steps, i.e., {Xt,Xt+1, . . . ,Xt+H} = F(Xt−T ,Xt−T+1, . . . ,Xt−1). Graph Learning Architecture The graph representation learning of our MPS2GCNets is build upon GCN. To learn node representation from graph topology and node features, the input of GCNbased approaches contains the adjacency matrix At of the original input graph Gt and the node feature matrix Xt. However, the prior knowledge of graph structure (1) might restrict the modeling capacity (i.e., graph edges cannot encode the complex relationships between nodes) and (2) leads to the neglect of neighboring information with high diversity. To avoid these limitations, inspired by the adaptive dependency matrix (Wu et al., 2019; Chen et al., 2021), we learn the normalized selfadaptive adjacency matrix S = {suv}N×N with the pre-defined “cost” staying in the same node based on the learnable node embedding Eφ = (e1,φ, e2,φ, . . . , eN,φ)> ∈ RN×dc as
suv =  exp(ReLU(eu,φe>v,φ))∑ v∈V\{u} exp(ReLU(eu,φe > v,φ))+exp (duu) , u 6= v, duu∑
v∈V\{u} exp(ReLU(eu,φe > v,φ))+exp (duu)
, u = v.
Here hyperparameter duu is the “cost” of staying in the same node u and ReLU(·) = max(0, ·) is an activation function, which guarantees suv ≥ 0. Since ∑ v suv = 1 and suv ≥ 0, we can use S as the normalized Laplacian. Towards more effective and robust learning of both spatial and spectral characteristics, we represent the graph diffusion as a matrix power series. That is, let S̃ be an N × N × K Laplacian tensor containing the power series {I,S, . . . ,SK−1} of S, where I ∈ RN×N is the identity matrix and K ≥ 2.
Spatio-Temporal Feature Transformation Note that spatial and temporal domains contain interdependent but highly heterogeneous types of information. Inspired by statistical factor analysis, we facilitate signal extraction from these disparate informational sources, by mapping the original input feature space into the high-level latent feature space, which can be written as follows
H (`+1) i,FT = (X TEφΘ (`) FT ) >U , (1)
where Θ(`)FT ∈ Rdc×P×QFT and U ∈ RT represent the learnable parameters.

heading: 4.1 SUPRAGRAPH CONVOLUTIONAL MODULE IN MULTILAYER SUPRA GRAPH
text: We propose a novel supragraph convolutional module to simultaneously capture spatio-temporal dependencies in dynamic networks. Our key approach is (i) to represent a sequence of time-varying graphs recorded over a sliding time window, as a multilayer supra graph, instead of treating each graph snapshot individually, and (ii) to learn the resulting multilayer supra graph with the random walk exploration which allows us to encode the key details of the time-conditioned relationships among nodes and to boost graph convolutions over multiple edge sets.
Sliding Window Historical Data as Multilayer Supra Graph Network To increase expressive capability of spatio-temporal representation learning, we propose a novel directed multilayer supra graph. Particularly, we treat each graph within a sliding window GT = {Gt−T ,Gt−T +1, . . . ,Gt−1} ∈ RN×N×T as a layer within a directed T -layered network. (Here N is the number of nodes in each graph.) Since information in real world spatio-temporal processes can propagate only from past to present to future, we consider directed multilayer supra graph as abstraction for dynamic knowledge representation. That is, we assume that the information (e.g., spatial features) is shared between layers ta and tb, whenever tb > ta, and time stamps ta, tb ∈ {t − T , . . . , t − 1}. More specifically, we propose a strategy to create such kind of information propagation channels by adding directed virtual edges (i.e., interlayer edges) between every node in the layer ta and its counterparts in other “future” layers tb (tb > ta). Definition 4.1. A directed multilayer supra graph (DMSG) is a tuple, defined as: DMSG =( Gt−T , . . . ,Gt−1, IMtatb ) , where Gta = (Vta , Eta), ta ∈ {t− T , . . . , t− 1} are network layers and IMtatb (Identity Mapping) is an N ×N matrix of node mappings, with IMtatbij : v ta i × v tb j 7→ [0, 1]. That is, IMtatbii is the identity mapping between node vi in the “past” layer ta and the ”future” layer tb (tb > ta): IMtatbii = 1.
The corresponding (r-th power) supra-adjacency matrix W = {wuv}NT ×NT for DMSG is then
wabuv =  (statauv ) r, ta = tb and |u− v| mod N 6= 0, (dtatauu )
r, ta = tb and |u− v| mod N = 0, dtatbuu , ta 6= tb and |u− v| mod N = 0,
(2)
where (statauv ) r is the r-th power of Stata which encodes the r-th step of a random walk among nodes u and v in layer ta, (dtatauu ) r is the “cost” of staying in the same node u and in the same layer ta after r random walk steps, and dtatbuu is the “cost” of jumping from the current node u in layer ta to node u in layer tb. Finally, the generalized (r-th power) supra-Laplacian LSup ∈ RNT ×NT for DMSG is then
LSup =  D11I+(S11)r D12I ··· D1T I
0 D22I+(S22)r ··· D2T I ... ... . . .
... 0 0 ··· DT T I+(ST T )r . (3) Supragraph Diffusion Convolutional Layer Armed with the (r-th power) supra-Laplacian in (3), we now build a supragraph diffusion convolutional layer to encode both the graph structure and node features from the temporal domain. That is, embedded node features are updated by message passing and aggregation via intralayer and interlayer diffusion. The supragraph diffusion convolutional layer is formulated as
H (`+1) i,Sup = (LSupH (`) i,Sup) >EφΘ (`) Sup, (4)
where H(`)i,Sup ∈ RNT ×P and H (`+1) i,Sup ∈ RNT ×QSup are the input and output activations for layer `, H (0) i,Sup = X
T , and Θ(`)Sup ∈ Rdc×P×QSup are the learnable parameters. In this case, the diffusion on layer ta can extend to a fraction of nodes and propagate information of layer tb by interaction.

heading: 4.2 SPATIAL GRAPH CONVOLUTIONAL LAYER
text: To aggregate the features of each node with its multi-hops neighbourhoods to generate node embeddings, we now turn to constructing the spatial graph convolution via
H (`+1) i,Spa = (S̃H (`) i,Spa) >EφΘ (`) Spa, (5)
where Θ(`)Spa ∈ Rdc×K×P×QSpa is the trainable weight, S̃ is the Laplacian tensor, H (`) i,Spa ∈ RN×P and H(`+1)i,Spa ∈ RN×QSpa are the node representations at the `-th layer and (` + 1)-th layer, respectively, and H(0)i,Spa = Xi, i.e., the node features at time step i. To reduce computational costs, we use weight sharing matrix factorization instead of matrix entry assignments. As such, we timely update the latest state of variables with back propagation algorithm and reduce the risk of over-fitting.

heading: 4.3 LEARNING TIME-AWARE MULTIPERSISTENCE REPRESENTATION OF TOPOLOGICAL FEATURES
text: To fully utilize information delivered by the MP representation mechanism, we propose a CNNinspired Time-Aware Multipersistence Euler-Poincaré Surface Representation Learning (DEPSRL) module which employs the CNN base model fθ (including convolutional and pooling layers) and readout layer. In summary, the proposed DEPSRL module offers the following unique innovations: (i) extracting features from the global information, (ii) learning the relationship from DEPS sequences pre and post (i.e., treating the fixed-size sliding window as multi-channels), and (iii) aggregating topological features to make a fixed size representation. The summarized output feature of DEPSRL module is given by
Φ (`+1) i,TAMP = ⊕(fGAP(fθ1({E i}Ti=1)), fGMP(fθ2({Ei}Ti=1))), (6)
where fGAP(·), fGMP(·), and fθj , j = {1, 2}, are global average pooling, global max pooling, and j-th CNN based model, respectively, and ⊕ denotes concatenation, where the fGAP(·) generates summarized feature for each channel. We find that fGMP(·) strengthens the representation learning of the time-conditioned multipersistence features. The matrix Φ(`+1)i,TAMP ∈ RQTAMP is the output of the DEPSRL module. To further learn and fuse multiple latent representations form different views (i.e., spatial information, spatio-temporal correlations, and persistent topological features), we combine these embeddings to obtain the final embedding
Z (`+1) i,out = F (H (`+1) i,Spa ,Υ (`+1) i ,H (`+1) i,FT ). (7)
Here Υ(`+1)i = (1/T ∑ H̃ (`+1) i,Sup )Φ (`+1) i,TAMP integrates global topological and global spatio-temporal information (i.e., we reshape H(`+1)i,Sup ∈ RNT ×QSup to H̃ (`+1) i,Sup ∈ RN×T ×QSup and then average H̃ (`+1) i,Sup over temporal dimension), F (·, ·, ·) is a dimension-wise concatenation function on embed-
dings along the output dimension, and ∑ i∈{Spa,Sup,FT}Qi = Qout. DEPSRL also allows for input of E from different bifiltrations. For more details, please refer to Appendix D.2.

heading: 4.4 MODELING SPATIO-TEMPORAL DYNAMICS
text: To encode spatio-temporal correlations among time series and get hidden state of nodes at a future timestamp, we put the final embedding Z(`+1)i,out into Gated Recurrent Units (GRU). The forward propagation equations of GRU are as follows
<i = ς (W< [Ξi−1,Zi,out] + b<) , =i = ς (W= [Ξi−1,Zi,out] + b=) , Ξ̃i = tanh (WΞ̃ [=i Ξi−1,Zi,out] + bΞ̃), Ξi = <i Ξi−1 + (1−<i) Ξ̃i, (8)
where ς(·) is an activation function which is ReLU in our case, is the elementwise product, <i is the update gate and =i is the reset gate. b<, b=, bΞ̃, W<, W=, and WΞ̃ are learnable parameters. [Ξi−1,Zi,out] and Ξi are the input and output of GRU model, respectively. Now we get Ξ̃i containing the structure, spatio-temporal, and topological information.

heading: 5 EXPERIMENTS
text: Data Description We use three different types of datasets to examine the performance of the proposed TAMP-S2GCNets on dynamic networks: (i) The transportation and traffic data of California state from the freeway Performance Measurement System (PeMS). We use three well-studied traffic networks from literature: PeMSD3, PeMSD4 and PeMSD8 (Chen et al., 2001), (ii) Digital transactions between users in the Ethereum blockchain network. We extract dynamic networks from three token assets: Bytom, Decentraland and Golem (Li et al., 2020), and (iii) The spread of coronavirus disease COVID-19 at county-level in states of California (CA) and Texas (TX). More details on each dataset and parameter settings are provided in Appendix D.1.
Experimental Settings We compare the presented TAMP-S2GCNets with 13 state-of-the-art baselines (see more details of baselines in Appendix D.2) on three evaluation metrics, i.e, Mean Absolute Error (MAE), Root Mean Squared Error (RMSE), and Mean Absolute Percentage Error (MAPE). For transportation networks, following Cao et al. (2020), we (i) split PeMSD3, PeMSD4, and PeMSD8 datasets into a training set (60%), validation set (20%), and test set (20%) in chronological order and (ii) use the one hour historical data to the next 15 minutes data. Following Chen et al. (2021), we (i) split Bytom, Decentraland, and Golem token networks into training set (80%) and test set (20%) in chronological order and (ii) use 7 days historical data to predict future 7 days data. For COVID-19 datasets, in our experiments, we use daily data of 11 months of 2020, from February 1 to December 31, and split the graph signals into training set, first 80% of days (268 days), and test set, last 20% of days (67 days). Further specifics on experimental setup are provided in Appendix D.2. The best results are in bold font and the results with dotted underline are the best performance achieved by the runner-ups. We also perform a one-sided two-sample t-test between the best result and the best performance achieved by the runner-up, where *, **, *** denote significant, statistically significant, highly statistically significant results, respectively. More detailed experiments on all datasets are in Appendix E. The data and code implementation are available at https://www. dropbox.com/sh/n0ajd5l0tdeyb80/AABGn-ejfV1YtRwjf_L0AOsNa?dl=0.
Computational Complexity Although, in general, MP is computationally costly, in our case, to get the DEPS summary, we do not need to compute PDs as in other MP approaches, but only Betti numbers of each filtration cell. In particular, while a computational cost for obtaining kthPD for a graph is O(M3k), where Mk is the number k-simplices (Otter et al., 2017), obtaining Euler Characteristics by sparse matrix methods has computational complexity of O(M0 +M1 + M2) (Edelsbrunner & Parsa, 2014).

heading: 5.1 EXPERIMENTAL RESULTS
text: Transportation Traffic Flow Table 1 shows the performance comparison among 13 state-of-the-art baselines on PeMSD3, PeMSD4, and PeMSD8 for multi-step traffic flow forecasting. Our TAMPS2GCNets consistently outperforms baseline models on all 3 datasets except for PeMSD3, which underscore effectiveness of TAMP-S2GCNets in time series forecasting tasks. On PeMSD3, our TAMP-S2GCNets achieves the best MAE and MAPE, and has an average 2.75% relative gain, compared to the runner-up. The average relative gains of TAMP-S2GCNets over the runner-ups in MAE and MAPE are 2.78% and 2.19% on the 3 datasets, respectively.
Ethereum Blockchain Prices Table 3 (left) summarizes the comparison results in MAPE on 3 Ethereum token networks (i.e., Bytom, Decentrland, and Golen). Table 3 indicates that TAMPS2GCNets is always better than baselines for all dynamic Ethereum networks. We find that, even compared to the baseline which also integrates topological information (i.e., Z-GCNETs), TAMPS2GCNets is highly competitive, implying that the time-aware MP representation with DEPS is capable of better encoding time-conditioned information than the zigzag idea based on one filtration.
COVID-19 Hospitalizations Table 2 shows results on 15-day ahead forecasting of COVID-19 hospitalizations in the U.S. states of California (CA) and Texas (TX) at a county level basis (RMSE is aggregated over each state). For the sake of room, here we only display the top runner models. TAMP-S2GCNets yields significantly better forecasting performance, with relative gains of 1.5%- 24.5% (in CA) and 5.5%-46.7% (in TX).

heading: 5.2 ABLATION STUDIES
text: Contribution of TAMP-S2GCNets components We now conduct the ablation studies on PeMSD4 and PeMSD8 to evaluate contribution of different components of our framework. (Ablation results on Golem and CA are in Table 9 in Appendix E.)
We compare our TAMP-S2GCNets with 4 ablated variants, i.e., (i) TAMP-S2GCNets without DEPSRL module (w/o DEPSRL module), (ii) TAMP-S2GCNets without spatial graph convolutional layer (w/o GCNSpa), (iii) TAMP-S2GCNets without supragraph convolutional module (w/o GCNSup), and (iv) TAMP-S2GCNets without FT (w/o FT). As Table 3 (right panel) shows, ablating each of above causes the
performance drops sharply in comparison with our full TAMP-S2GCNets model. Especially, on PeMSD4, the supragraph convolutional module significantly improves the results as it simultaneously captures the spatial and temporal information. In addition, it is evident that DEPSRL mod-
ule enhances the topological information learning ability in spatio-temporal domain, i.e., TAMPS2GCNets outperforms TAMP-S2GCNets w/o DEPSRL module with an average relative gain 3.92% on RMSE over PeMSD4 and PeMSD8. Besides, as expected, taking the advantages of GCNSpa and FT, we improve the performance via capturing the graph structure in the spatial dimension and consolidating the processing of spatio-temporal correlations between node attributes, respectively.
Contribution of Single vs. Multifiltration Persistence We consider sublevel filtrations based on node degree (Deg), betweenness (Btwns), and graph diameter (power filtration) with edge weights induced by number of transactions (PowerTrns) and volume of transactions (PowerVol). As Table 4 (left) shows, MP always outperforms single filtrations by a significant margin, both in terms of prediction accuracy and variability. The results confirm our premise that MP is able (i) to capture hidden time dependencies of high dimensional time-varying objects which are inaccessible with one-parameter PH and (ii) to yield substantially more stable feature maps in dynamic scenarios.
Contribution of DEPS vs. the Existing MP Summaries Finally, we evaluate performance of the new time-aware MP summary, i.e., DEPS, in comparison to the MP summaries based on the slicing argument, namely, MP-I of Carrière & Blumberg (2020) and MP-L of Vipond (2020) as representation input to TAMP-S2GCNets. As Table 4 (right panel) shows, comparing to MP-I and MP-L, DEPS yields 9%-12% improvement in MAPE, and at least a 10 times decrease in computational time (29.84s for DEPS vs. 401.80s for MP-I and 517.11s for MP-L). Such substantial differences in computational costs are explained by the need of the slicing-based MP summaries to search for the most suitable one-parameter PH representation of MP. In turn, DEPS is based on a point-wise representation argument in linear algebra, and its high computational efficiency makes DEPS the preferred choice for time-conditioned topological representation learning (see the running time in Table 4). Since computational costs remain the major roadblock for MP applications, our DEPS approach shows that MP tools based on scalable sparse matrix algorithms appear to be the most promising direction for integration of MP into machine learning tasks.

heading: 6 CONCLUSION
text: We have explored utility of MP to enhance knowledge representation mechanisms within the timeaware DL paradigm. The developed TAMP-S2GCNets model is shown to yield highly competitive forecasting performance on a wide range of datasets, with much lower computational costs. In the future we plan to explore combination of MP with the zigzag persistence and to investigate asymptotic distributional properties of point-wise MP invariants.
ETHICS STATEMENT
We do not anticipate any negative ethical implications of the proposed methodology. In turn, we believe that introduction of GCNs tools, coupled with multiparameter topological approaches, into biosurveillance opens a new path for more accurate, timely, and robust tracking of infectious diseases with high virulence such as COVID-19. In particular, multiparameter persistence allows for extracting the most salient topological features along multiple geometric dimensions and, as such, can be especially valuable to address disease dynamics as a function of multiple highly heterogeneous variables, e.g., socio-demographic, socio-environmental, and socio-economic factors. In turn, GCNs and, more generally, geometric deep learning allow for capturing sophisticated nonlinear spatio-temporal dependencies among factors which contribute to the disease spread. Hence, we postulate that in the next few years we can expect to see a new set of spatio-temporal biosurveillance artificial intelligence algorithms, based on a combination of geometric deep learning, multiparameter persistence, and more generally, tools of topological data analysis.
REPRODUCIBILITY STATEMENT
To support reproducibility of the results in the paper, we have submitted our datasets and codes (see https://www.dropbox.com/sh/n0ajd5l0tdeyb80/AABGn-ejfV1YtRwjf_ L0AOsNa?dl=0) as the supplementary information.

heading: A NOTATION
text: Frequently used notation is summarized in Table 5.

heading: B SINGLE PARAMETER PERSISTENT HOMOLOGY AND ITS SUMMARIES
text: All extracted topological features can be then summarized as a multiset in R2, called persistence diagram (PD): Dk(C) = {(bi, di) ∈ R2 | bi < di}, where birth bi and death di mark the filtration indexes at which the k-dimensional topological feature ρi appears and disappears, respectively. The farther (bi, di) from the diagonal is (that is, the more persistent ρi is), the likelier ρi is to contain salient information on Gt. Multiplicity of (xb, yd) ∈ D is the number of p-dimensional topological features (p-holes) that are born and die at xb and yb, respectively. Features on the diagonal of D have infinite multiplicity.
Another summary which may be particularly suitable in our context of modeling noisy time-evolving objects is the Euler-Poincaré Characteristics (see Adler et al. (2010); Baryshnikov & Ghrist (2010) on the Euler calculus and its applications, particular, in conjunction with sensor networks).
Definition B.1 (Euler-Poincaré Characteristics). For a given simplicial complex C, Euler-Poincaré Characteristics χ(C) is defined as the alternating sum of the number of k-simplices of C. That is, if nk is the number of k-simplices in C, then χ(C) = ∑M k=0(−1)knk, where M is the dimension of highest dimensional simplex in C. Note that Euler-Poincaré Characteristics is homotopy invariant,
and an alternative definition is χ(C) = ∑M k=0(−1)kBk(C) where Bk(C) is the kth Betti number of C, i.e. rank of the kth homology group Hk(C) (Ghrist, 2018; Hatcher, 2002).

heading: C PROOF OF THEOREM 3.1
text: In this section, we provide the proof of the stability theorem (Theorem 3.1). First, let us recall the key notations we use. Let G+ and G− be two graphs and let C± be the clique complexes of G± (Ghrist, 2018). Let F = (f, g) be a multivariate filter function, i.e. F : V± → R2, where V± is the set of nodes of G±, and let I = {(αi, βj) | 1 ≤ i ≤ m, 1 ≤ j ≤ n} be the corresponding thresholds for F = (f, g).
Next, we define Wasserstein-p distance among PDs, i.e., the important concept we borrow from the theory of single parameter persistence (Edelsbrunner & Harer, 2010).
Definition C.1 (Wasserstein-p distance). Let (C±, f±, I±) be two single parameter filtrations, and Dk(C±) be the corresponding PDs for k-cycles (i.e, k-dimensional topological features). Let q+j = (b+j , d + j ) ∈ Dk(C+) be the birth and death times of a k-dimensional topological feature ρj . Then, Wasserstein-p distance between Dk(C+) and Dk(C−) is defined as
Wp(Dk(C+),Dk(C−)) = inf φ ( ∑ j∈D(C+)∪∆ ‖q+j − φ(q + j )‖ p ∞ ) 1 p , p ∈ Z+, (9)
where φ : Dk(C+) ∪ ∆ → Dk(C−) ∪ ∆ is a bijection (matching) and ∆ is the diagonal set, i.e., ∆ = {(b, b)|b ∈ R}, which contains k-cycles of infinite multiplicity. With ∆ in both sides, we ensure the existence of these bijections even if the cardinalities |{q+j }| and |{q − l }| are different. When p = ∞, (9) corresponds to the bottleneck distance, i.e. W∞. In the proof below, we use p = 1, i.e. W1-distance.
Armed the PD and Wasserstein-p distance concepts, we now turn to the main stability result for the new MP summary. Let Ĉ± = {C±ij} be the bifiltration associated with (C±, F, I), i.e. C ± ij is the clique complex of the induced subgraph G±ij = F−1((−∞, αi] × (−∞, βj ]). (Figure 2 illustrates such bifiltration on graphs with degree and eccentricity functions.) Let E± be the the corresponding Euler-Poincaré Surface, i.e. E± is an m × n-matrix with entries E±ij = χ(C ± ij , ), where χ(.) is the Euler Characteristics (See Figure 2), and D(C+,C−) be the weak L1-distance as defined in Section 3 of the paper (see also Remark C.1 below for the motivation to use this new metric).
Proof of Theorem 3.1: We start from providing the outline of the proof. Notice that Euler-Poincaré Surfaces is the alternating sum of Betti numbers {Bk}. Hence, the distance ‖E+−E−‖1,1 is bounded by the difference of the corresponding Betti numbers for Ĉ±, i.e. Eqs. (10) and (11). By fixing i and k, and by taking column distance Dc into consideration, our main assertion
‖E+ − E−‖1,1 ≤ C ·D(Ĉ+, Ĉ−), C > 0
then reduces to Eqs. (12) and (13). To relate the difference |B(C+j ) − B(C − j )| to W1(D(C+),D(C−)), we observe that the Betti number B(C±j ) is the sum of the number of barcodes inD(C±), containing βj , i.e., #{r | βj ∈ [b±r , d±r )} (Eqs. (13) and (14)). Yet, the differences of the number of barcodes containing βj for C+ and C− can be bounded by the Wasserstein-1 distances of D(C+) and D(C−), i.e. Eqs. (15) and (16). Aggregating these terms on indices i and k finishes the proof.
Now, we turn to derivations in details. Recall that χ(C) = ∑M k=0(−1)kBk(C), where Bk(C) is the kth Betti number of C and M is the maximum homological dimension of C. Then, since E±ij = χ(C ± ij), we have E ± ij = ∑M k=0 Bk(C ± ij). Hence, for any i, j,
|E+ij − E − ij | = | M∑ k=0 (−1)k[Bk(C+ij)− Bk(C − ij)]| ≤ M∑ k=0 |Bk(C+ij)− Bk(C − ij)|. (10)
Since ‖E+ − E−‖1,1 = ∑m i=1 ∑n j=1 |E + ij − E − ij |, we have
‖E+ − E−‖1,1 ≤ m∑ i=1 n∑ j=1 M∑ k=0 |Bk(C+ij)− Bk(C − ij)|. (11)
Now, to complete the proof, we need to bound the latter sum by the MP distance D(Ĉ+, Ĉ−). Recall that by Definition 3.3
Dc(Ĉ +, Ĉ−) = M∑ k=0 m∑ i=1 W1(Dgk(C + i∗),D g k(C − i∗)) ≤ D(Ĉ +, Ĉ−).
Then, by Eq. 11, to prove the theorem, it is enough to show that for any 0 ≤ k0 ≤ M and for any 1 ≤ i0 ≤ m,
n∑ j=1 ∣∣Bk0(C+i0j)− Bk0(C−i0j)∣∣ ≤ C · W1(Dgk0(C+i0∗),Dgk0(C−i0∗)), C > 0. (12) For simplicity, without loss of generality, we fix 1 ≤ i0 ≤ m and 0 ≤ k0 ≤ N . Also, for sake of notations, we drop the subscripts and superscripts, i.e C± = C±i0∗ and C ± j = C ± i0j
, D(C±) = Dgk0(C ± i0∗) and B(.) = Bk0(.). With these simplified notations, Eq. 12 is equivalent to
n∑ j=1 ∣∣B(C+j )− B(C−j )∣∣ ≤ C · W1(D(C+),D(C−)) C > 0. (13) To show (13), we need to evaluate the relation between the Betti numbers B(C±j ) and the persistence diagrams D(C±). Consider D(C) = {(br, dr) | 1 ≤ r ≤ Q}, where br represents the birth time of a k0-cycle ρr in C, while dr represents the death time of ρr. In particular, ρr induces a barcode [br, dr). For each such ρr, define indicator function ψr(t) = I[br,dr)(t), i.e. ψr(t) = 1 for t ∈ [br, dr) and ψr(t) = 0 otherwise. Then, the Betti function of C is defined as B : [β1,∞) → N such that B(t) = ∑Q r=1 ψr(t) for t ∈ [β1,∞). (Here, β1 is the lowest threshold for the filtering function g in the set of thresholds I = {(αi, βj) | 1 ≤ i ≤ m, 1 ≤ j ≤ n}. Similarly, if we consider row sum Dr instead of column sum Dc, then we would use the other filtering function f , and the Betti functions would be defined as B : [α1,∞)→ N.) Hence,
B(C±j ) = Q±∑ r=1 ψ±r (βj), j ∈ Z+.
As a result, we get
n∑ j=1 ∣∣B(C+j )− B(C−j )∣∣ = n∑ j=1 ∣∣∣∣Q +∑ r=1 ψ+r (βj)− Q−∑ r=1 ψ−r (βj) ∣∣∣∣ (14) ≤
n∑ j=1 Q̂∑ r=1 |ψ+r (βj)− ψ−φ(r)(βj)|,
where φ : D(C+) → D(C−) is the best matching in the definition of W1(D(C+),D(C−)) and Q̂ = max{Q+, Q−}.
Notice that the contribution of each |ψ+r (βj) − ψ−φ(r)(βj)| to the sum at βj is less than the number (#r) of intervals [b+r , d + r )4[b+φ(r), b − φ(r)) containing the interval [βj , βj+1), where 4 is the symmetric difference of the sets. Let K = minj(βj+1 − βj). By dividing the length ‖[b+r , d+r )4[b+φ(r), b − φ(r))‖ by K, we ensure that every such interval contributes at least 1 to the Betti number count. Hence, we arrive at the following key inequality:
n∑ j=1 Q̂∑ r=1 |ψ+r (βj)− ψ−φ(r)(βj)| ≤ 1 K · ∑ r |b+r − b−φ(r)|+ |d + r − d−φ(r)|. (15)
Finally, since Q̂∑ r=1 |b+r − b−φ(r)|+ |d + r − d−φ(r)| ≤ 2W1(D(C +),D(C−)),
we obtain n∑ j=1 ∣∣B+(βj)− B−(βj)∣∣≤ 2 K · W1(D(C+),D(C−)). (16)
Since (16) holds for any 1 ≤ i0 ≤ m and 0 ≤ k0 ≤ N , for C = 2/K we get
m∑ i=1 n∑ j=1 M∑ k=0 ∣∣Bk(C+ij)− Bk(C−ij)∣∣≤ C ·Dc(Ĉ+, Ĉ−). By Eq. (11), this implies
‖E+ − E−‖1,1 ≤ C ·D(Ĉ+, Ĉ−), which completes the proof.
Remark C.1 (On the Choice of the Metric). While interleaving and matching distances are very common and useful for multipersistence in various cases (Thomas, 2019), they are not suitable to study the Euler-Poincaré Surfaces. The reason for this is that matching and interleaving distances are based on L∞ metrics, and do not see the quantity of generators of small interleaving distances. However, Euler-Poincaré Surfaces by definition depends on the rank of the homology groups, and the quantity of the generators (Betti numbers) are crucial in their definition even if they exhibit a very short lifespan. In particular, L∞-based metrics find the maximum distance between the essential generators, but such metrics fail to capture the quantity of small generators. That is, while the interleaving or matching distance between the multipersistence modules of the bifiltrations {Ĉ+} and {Ĉ−} is small, there might be a very large distance between the induced Euler-Poincaré Surfaces E+ and E−. In turn, Lp-based distances are successful to capture the rank even if they are small (noise) generators, while interleaving and matching distances are not very sensitive to these. (See Thomas (2019, Section 4.2) for further discussion.) This is why we introduce the weak L1 metric above which naturally suits to study the Euler-Poincaré Surfaces and similar vectorizations in the theory. Furthermore, such weak L1 metric can be considered a straightforward generalization of L1-based norm for multiparameter persistence.
Remark C.2 (Implications of the Stability Result). Note that our stability result (Theorem 3.1) directly implies that the distances between multiparameter PDs control the distance between the resulting Euler-Poincaré Surfaces. By combining with the stability result for PDs (Cohen-Steiner et al., 2007), one can conclude that the small changes in the MP filtering function F : Vt 7→ Rd or the small changes in the input data can result only in a small change in DEPS surfaces. In particular, by Cohen-Steiner et al. (2007), for single parameter persistence, we have
W(Df (C), Dg(C)) ≤ Cḋ(f, g),
where Df is the persistence diagram for filtering function f , and d(f, g) is the distance between the functions f and g. The way we define our weak L1-metric enables us to use this single persistence result in this context as follows. The column distance Dc is the sum of single PD for each column. Each column, (say jth column) is represented by a single variable filtering function F (t, βj). Then, the changes in F (t, βj) reflects in the column distance for PDs corresponding to jth-column. This means if we have a small change in F (t, βj) for each j, then the column distance for MP will be small. So is the row distance induced by F (αi, s).
Similarly, if we have two MP filtering functions F± : V → R2, they will induce two filtrations C±. When d(F+, F−) is small, then by Cohen-Steiner et al. (2007), D(Ĉ+, Ĉ−) is small. When we apply our stability theorem to these filtrations, we see that the distance between the induced DEPS surfaces ‖E+ − E−‖1,1 will be small, too. Since small change in the data can be converted to small change in the filtering functions on the same data, the result can also be interpreted as the small change in the data results in small change in the induced DEPS surfaces.
Remark C.3 (DEPS through the Probabilistic Perspective). Furthermore, the new time-aware MP summary can also be assessed from a probabilistic perspective, if we view a multivariate time series X as real-valued random process indexed by time t and defined on a probability space (Σ,A,P). For 1 ≤ J ≤ L ≤ ∞, define the σ-field zLJ ⊂ A such that zLJ = σ(Xk, J ≤ k ≤ L (k ∈ Z)), i.e. zLJ is a σ-field generated by XJ ,XJ+1, . . . ,XL. Then, Et is a real-valued zt1-measurable random matrix and DEPS {Et}Tt=1 is F-measurable random matrix, where F = (zk1)k>0 is a filtration of σ-fields indexed over time. Given the results of Thoppe et al. (2016) on weak convergence of the Euler characteristic χ(C) for dynamic Erdős-Rényi graphs to the Ornstein-Uhlenbeck process, this opens a potential path to establish asymptotic distribution of Et. Remark C.4 (DEPS for Higher Dimensions). DEPS can be naturally expanded to a case when we consider filtrations along more than two geometric dimensions. In such a case, instead of a m× nmatrix, DEPS will be a tensor. Our stability result (i.e., Theorem 3.1) verbatim extends to such a case.

heading: D REPRODUCIBILITY
text: D.1 DATASETS
We use three different datasets, from transportation, Ethereum blockchain, and biosurveillance domains, to examine the performance of the proposed TAMP-S2GCNets methodology:
1. Transportation: The transportation and traffic data for the U.S. State of California, from the freeway Performance Measurement System (PeMS) data sources (i.e., PeMSD3, PeMSD4, PeMSD8) (Chen et al., 2001), naturally produces a dynamic network where each node is a loop detector and each edge represents a freeway between two nodes. We extract aggregated data to 5 minutes, thus generating dynamic networks with 26,208, 16,992 and 17,856 nets in PeMSD3, PeMSD4 and PeMSD8 datasets, respectively. To capture spatio-temporal dependencies, we rebuild the traffic graph structure, at time t, and compute edge weights wt,uv , between pair nodes (u, v), via the Radial Basis Function wt,uv = e−||xt,u−xt,v||
2/γ ; where γ = 1.0 and each node xt,· has features speed, occupancy and flow rate. To investigate the dynamic of the traffic graph topology, we restrict the final graphs to only keep edges whose are below or equal to threshold α = 0.01.
2. Blockchain: Ethereum blockchain data contain token assets which naturally represents network layers of the Ethereum network (di Angelo & Salzer, 2020; Web, e), in which nodes and edges are addresses of users and digital transactions, respectively. Using the publicly available blockchain, we extract dynamic networks from three tokens with market value above $100,000,000: Bytom, Decentraland and Golem. In our experiments, we extract dynamic networks by daily transactions and use daily closed prices (Web, d), hence, creating a graph/net for each day. Since each token is created at different date, in all cases last day is 05/07/2018, the dynamic networks vary in the number of nets: Bytom (285 nets), Decentraland (206 nets) and Golem (443 nets). For the sake of reasonable computation time, we reduce the original net size using the maximum weight subgraph approximation method of Vassilevska et al. (2006) from more than 400,000 nodes, average of 120,000 edges, to 100 nodes using number of transactions and price-volume as edge weights.
3. Biosurveillance: SARS-CoV-2, which is the virus that causes COVID-19, has spread around every place in the world. To analyze the progression of coronavirus disease we create a dynamic network based on daily records on COVID-19 cases and hospitalizations from the CovidActNow project (Web, c) and Johns Hopkins University (Web, f; Dong et al., 2020), and population numbers from the U.S. Census Bureau (Web, b). We choose data at county level from California and Texas states as our cases of study, and focus on forecasting the number of hospitalized patients. Similarly to (i) we rebuild the coronavirus-spread graph structure, at day-time t, based on the county adjacency graph (Web, a) and Radial Basis Function between pair of nodes. We build a dynamic network for each state, each dynamic network contains 335 nets. The number of confirmed COVID-19 cases and county population size serves as node features, and final graphs only keep edges whose are below or equal to the threshold α = 0.01.
More details on experimental settings and filtering functions are provided in Section D.2. Table 6 summarize the characteristics of all datasets used in our experiments.
D.2 EXPERIMENTAL SETTINGS
We implement our TAMP-S2GCNets with Pytorch framework on NVIDIA GeForce RTX 3090 GPU. Further, for all datasets, TAMP-S2GCNets is trained end-to-end by using the Adam optimizer with a L1 loss function. The tuning of our proposed TAMP-S2GCNets on each dataset is done via grid search over a fixed set of choices and the same cross validation setup is used to tune baselines. We compare TAMP-S2GCNets with 13 types of state-of-the-art methods, including
FC-LSTM (Sutskever et al., 2014), SFM (Zhang et al., 2017), N-BEATS (Oreshkin et al., 2019), DCRNN (Li et al., 2018), LSTNet (Lai et al., 2018), STGCN (Yu et al., 2018), TCN (Bai et al., 2018), DeepState (Rangapuram et al., 2018), GraphWaveNet (Wu et al., 2019), DeepGLP (Sen et al., 2019), AGCRN (Bai et al., 2020), Z-GCNETs (Chen et al., 2021), and StemGNN (Cao et al., 2020).
MODEL SPECIFICATIONS For PeMSD3, TAMP-S2GCNets consists of 2 layers whose hidden feature dimension is 64. The learning rate is 5e-3 and the batch size is set as 64. The “cost” dtatauu of staying in the same node u (in the layer ta) and the “cost” dtatbuu of jumping from the current node u in layer ta to node u in layer tb are set to be 0.1 and 0.2, respectively.
For PeMSD4, TAMP-S2GCNets consists of 2 layers whose hidden feature dimension is 128. The learning rate is 5e-3 and the batch size is set as 64. The “cost” dtatauu of staying in the same node u (in the layer ta) and the “cost” dtatbuu of jumping from the current node u in layer ta to node u in layer tb are set to be 0.1 and 0.2, respectively.
For PeMSD8, TAMP-S2GCNets consists of 2 layers whose hidden feature dimension is 64. The learning rate is 2e-3 and the batch size is set to 16. The “cost” dtatauu of staying in the same node u (in the layer ta) and the “cost” dtatbuu of jumping from the current node u in layer ta to node u in layer tb are set to be 0.1 and 0, respectively. The dimension of node embedding dc (in Eφ) are 10 (PeMSD3), 10 (PeMSD4), and 20 (PeMSD8), respectively. We set the dimension K of Laplacian tensor S̃ as 2 for all transportation networks.
For all Ethereum token networks (i.e., Bytom, Decentraland, and Golem), TAMP-S2GCNets consists of 2 layers whose hidden feature dimension is 32. The final hyperparameter setting is learning rate of 1e-2 and batch size of 8. The “cost” dtatauu of staying in the same node u (in the layer ta) and the “cost” dtatbuu of jumping from the current node u in layer ta to node u in layer tb are set to be 0.01 and 0.01, respectively. In addition, we set the embedding dimension dc in Eφ as 1 for all Ethereum token networks. We set the dimension K of Laplacian tensor S̃ as 2, 3, 3 for Bytom, Decentraland, and Golem, respectively.
For COVID-19 hospitalizations in CA and TX, TAMP-S2GCNets consists of 2 layers whose hidden feature dimension is 128. The learning rate is set to 1e-1 and 1e-2 for CA and TX, respectively. The batch size is set to 8, and the embedding dimension dc in Eφ is 10. The “cost” dtatauu of staying in the same node u (in the layer ta) and the “cost” dtatbuu of jumping from the current node u in layer ta to node u in layer tb are set to be 0.1 and 0.2, respectively. We set the dimension K of Laplacian tensor S̃ as 2 for both CA and TX.
ALTERNATIVE OPTIONS FOR DEPS WITHIN DEPSRL MODULE In our experiments, for all datasets, the CNN based model (i.e., fθj where j = {1, 2} in Time-Aware Multipersistence Euler-Poincaré Surface Representation Learning (DEPSRL) module) consists of 2 CNN layers. For transportation networks and COVID, The filter size, kernel size, and stride is set to be 8, 2, 2 respectively. For Ethereum token networks (i.e., Bytom, Decentraland, and Golem), the filter size, kernel size, and stride is set to be 8, 1, 2 respectively. We set the size of global average pooling and global max pooling as 3 × 3. Table 7 shows the optimal output dimensions of spatial graph convolutional layer (QSpa), supragraph diffusion convolutional layer (QSup), and spatio-temporal feature transformation (QFT) on all datasets.
Table 7: The output dimensions of spatial graph convolutional layer (QSpa), supragraph diffusion convolutional layer (QSup), and spatio-temporal feature transformation (QFT) on different datasets.
Output dimension PeMSD3 PeMSD4 PeMSD8 Bytom Decentraland Golem CA TX
QSpa Qout 4 Qout 4 3Qout 4 Qout 2 Qout 2 3Qout 4 Qout 4 Qout 4 QSup Qout 4 Qout 4 Qout 8 Qout 4 Qout 4 Qout 8 Qout 4 Qout 4 QFT Qout 2 Qout 2 Qout 8 Qout 4 Qout 4 Qout 8 Qout 2 Qout 2
Furthermore, in addition to the Time-Aware Multipersistence Euler-Poincaré Surface Representation Learning (DEPSRL) module (see Section 4.3 in our main paper):
Φ (`+1) i,TAMP = ⊕(fGAP(fθ1({E i}Ti=1)), fGMP(fθ2({Ei}Ti=1))),
we also have considered passing the aggregated DEPS as input to DEPSRL, that is, fθj (Ē), j = 1, 2, where Ē = 1/T ∑T i=1 ∈ Rm×nEi.
As expected, DEPS {Ei}Ti=1 tends to yield better performance within the DEPSRL module than its aggregated counterpart Ē, which can be explained by lower loss of time-conditioned topological information within DEPS. However, in the case of PeMSD4 dataset, we find that Ē slightly improves the performance of TAMP-S2GCNets (i.e., the MAE of TAMP-S2GCNets based on the aggregated DEPS is 17.58 and the MAE of TAMP-S2GCNets based on DEPS is 17.69), which is likely due to higher data heterogeneity of PeMSD4 vs. other transportation datasets. In all other datasets, we report results based on DEPS.
COMPARISON TO OTHER MP SUMMARIES For fair comparison, computation of all MP representations (as shown in Table 11 and Table 4 (right panel) in our main paper) have been run on an 8-cores DO droplet machine with Intel Xeon Scalable processors at base frequency of 2.5 Ghz. Our source codes are available for revision1. Source codes for computation of MP-I and MP-L are freely available in Github2. We follow the filtering function recommendations for graphs, on both MP-I and MP-L, from supplementary material of Carrière & Blumberg (2020). Further technical details on the 8-core machine are publicly available3.
ON TYPES OF MULTIFILTRATIONS Good filtrations should have rich range to get a fine resolution, and it should reflect an important property relevant to the question at hand. For multipersistence, for a good combination of two such filtering functions, one should choose functions uncorrelated. One approach is to choose one function from the domain of the question, and the other from strictly graph properties (e.g. degree, betweenness, etc.). We assess shape characteristics of Gt via filtering functions based on node degree (Deg), betweenness centrality (Btwns), and edge weights (Power). Edge weights for each dataset are computed as described in D.1; i.e using transactions between addresses (Trns), volume of transactions (Vol), traffic measurements (Traffic) and COVID-19 incidence (Incidence). Our analysis focuses on multivariate filtering function F : Vt 7→ Rd with d = 2. For Deg and Btwns the set of nondecreasing thresholds runs on values over nodes, whilst for filtering via power filtration (Power) of Gt nondecreasing thresholds runs on weighted paths between nodes. Notice that each pair-combination of filtering functions produces a specific MP representation in a form of the Euler-Poincaré surface, which impacts the algorithm performance, as shown in our multiple experiments. Table 12 presents the optimal multifiltration for each dataset.

heading: E ADDITIONAL EXPERIMENTS AND RUNNING TIME
text: Why MP Helps? More Experiments To validate gains (if any) delivered by multifiltration compared to the single filtration persistence and ensemble of multiple stacked together single parameter
1https://www.dropbox.com/sh/n0ajd5l0tdeyb80/AABGn-ejfV1YtRwjf_L0AOsNa? dl=0
2https://github.com/MathieuCarriere/multipers 3https://www.digitalocean.com/blog/premium-droplets-intel-cascade-lake-amd-epyc-rome/
filtrations, we also report the results of single filtration, filtration ensemble, and multifiltration on Decentraland (see Table 10). (In case of a filtration ensemble, we aggregate two single filtrations and feed the ensemble representation into the model.) We find that MP substantially outperforms both the single filtration persistence and the filtration ensemble. In addition, the MP results are generally more stable than ones based on the single filtration and filtration ensemble. This result can be expected, as by construction of DEPS, each individual single parameter persistence corresponds either to a row or a column of DEPS, that is, we have only a one-dimensional fingerprint of the data. In contrast, in MP and its DEPS summary, we obtain a two-dimensional fingerprint of the data, meaning any cell in the output simultaneously related to the all 8 neighboring cells (East, West, Northeast, etc.). As such, MP provides a 2D resolution of the data space and delivers much finer information on the underlyng data topology.
Table 11 shows the computational costs and performance of different MP summaries on Decentraland and Golem datasets. We observe that our new time-aware MP summary DEPS not only outperforms the existing MP summaries (i.e., MP-I and MP-L) but also induces substantially lower computational costs. These findings demonstrate that our DEPS can effectively capture salient timeaware topological information with lower time and memory cost.
COVID-19 Biosurveillance: More Experiments Finally, Figure 3 depicts the results on 15-day ahead forecasting of hospitalizations in California on a county level basis, using the top 3 models (i.e., AGCRN (Bai et al., 2020), StemGNN (Cao et al., 2020), and our TAMP-S2GCNets). (Here, to differentiate among RMSEs on a county level basis, we omit the map for DCRNN (Li et al.,
2018) due to its lower predictive performance, relative to RMSE scale of AGCRN, StemGNN, and our TAMP-S2GCNets.) Overall, TAMP-S2GCNets and StemGNN tend to perform similarly for less populated counties in Northern California, with TAMP-S2GCNets being slightly better than StemGNN over the Pacific Coast and StemGNN yeilding some edge in the North East. In turn, AGCRN tends to deliver more competitive performance in North-Eastern California, except of Shasta county.
However, performance among models starts to differ drastically in more populated counties and, in particular, in Southern California where it is much harder to forecast COVID-19 hospitalizations due to higher population density and more disparate socio-economic factors. Here, TAMP-S2GCNets
substantially outperforms both StemGNN and AGCRN almost in all counties of Southern California, yielding predictive gains from 3.8% vs. StemGNN and 29.8% vs. AGCRN in San Bernardino to 57.7% vs. StemGNN and 55.7% vs. AGCRN in San Diego. While all models perform very poorly in Los Angeles county (i.e., one of the most populated counties), TAMP-S2GCNets still yields the lowest RMSE. One important phenomenon we observe is that TAMP-S2GCNets systematically tends to substantially outperform StemGNN and AGCRN in more highly populated counties, including both Southern and Northern California. For instance, the predictive gains of TAMP-S2GCNets in Santa Clara are 39.2% vs. StemGNN and 67.9% vs. AGCRN. These findings indicate that timeaware MP topological information contained in TAMP-S2GCNets can address hidden relationships among spatio-temporal factors contributing to COVID-19 dynamics which other models cannot.

heading: F ADDITIONAL DETAILS OF TAMP-S2GCNETS ARCHITECTURE AND EXPERIMENTS
text: Remark F.1 (Spatio-Temporal Feature Transformation). We utilize and share the learnable node embedding Eφ to (i) construct the normalized self-adaptive adjacency matrix S and (ii) project the target node from the input feature space (in spatial and temporal dimensions) into a high-level latent space (together with ΘFT ). Hence, this shared-update weights can not only mitigate the risk of overfitting but allow us to capture shared patterns among spatial-temporal features and graph structures over the time axis. The learnable node embedding Eφ in TAMP-S2GCNets can be obtained from the backward propagation through time.
We use Einstein summation convention in Eq. 1 for the multiplication of XT ∈ Rb×T ×N×P (where b is the batch size) and EφΘ (`) FT . The implementation of matrix multiplication in Eq. 1 is as follows: (i) we first perform the matrix multiplication between Eφ ∈ RN×dc and Θ(`)FT ∈ Rdc×P×QFT to get an intermediate weight term denoted by Θ̃(`)FT ∈ RN×P×QFT ; (ii) then we compute the matrix multiplication between XT and Θ̃(`)FT through Einstein summation convention, i.e., in Pytorch, output = torch.einsum(‘btni, nio → btno’, X, intermediate weight term) (where X is XT and intermediate weight term is Θ̃(`)FT ).
Remark F.2 (Normalized Self-Adaptive Adjacency Matrix). We utilize the node embedding dictionaries Eφ ∈ RN×dc to construct the normalized self-adaptive adjacency matrix S via multiplying Eφ and E>φ . As such, construction of the normalized self-adaptive adjacency matrix S does not impose any requirements on density of the input graph. That is, our graph learning architecture (i.e., the self-adaptive adjacency matrix construction mechanism) overcomes this problem of the dense input graph met by GAT.
To explore more computationally efficient mechanisms of learning the normalized self-adaptive adjacency matrix, in the future work, we will (i) apply the sparse sampling strategies to the graph learning architecture; (ii) learn the graph based on the original graph structure information, i.e., we only infer the dependencies between the adjacent nodes (i.e., there exists an edge between them in the original graph) via multiplying eu,φ and e>v,φ; (iii) construct hypergraph representation instead of graph which not only alleviate the computational complexity but capture heterogeneous higher-order structures information.
For S̃, in our experiments, we performed an extensive grid search for the hyperparameter K (i.e., truncation order). We find the optimal K of 2 for the transportation network, the optimal K of 2 for COVID-19 network dataset, and the optimal K of 3 for Ethereum blockchain network. That is, for larger networks the truncation order K tends to be lower so that the computational costs are mitigated.
Furthermore, we have performed ablation studies on PeMSD4 and PeMSD8 to assess utility of the normalized self-adaptive adjacency matrix S in TAMP-S2GCNets. The results in Tables 13 and 14 indicate that TAMP-S2GCNets outperforms TAMP-S2GCNets with the plain adjacency matrix (instead of the normalized self-adaptive adjacency matrix) for both PeMSD4 and PeMSD8 across all the evaluation metrics. The results illustrate the importance of learning the node-specific patterns/embeddings over the spatial-temporal time series dataset.
Remark F.3 (Global Average Pooling and Global Max Pooling). We utilize the global average pooling (GAP) to preserve information about the Euler-Poincaré surfaces in a fixed-size representation through averaging of the learnt topological embeddings. In our experiments, we find that adding global max pooling (GMP) tends to strengthen the Euler-Poincaré surfaces representation and to improve the overall performance. Before GAP and GMP, we first use CNN base models (i.e., fθ1(·) and fθ2(·)) to learn the topological features of Euler-Poincaré surfaces (i.e., fθi({Ei}Ti=1), where i = {1, 2}); then we employ GAP and GMP to the outputs of CNN base models, respectively. We have conducted the ablation study to evaluate the performance of TAMP-S2GCNets without (w/o) GAP or GMP on Golem. From Table 15, we find that applying both GAP and GMP enhances prediction performance.
Remark F.4 (Supra-Laplacian). In our experiments, for large-scale networks, we apply sparse sampling strategy for the supra-Laplacian construction. For example, on transportation network such as PeMSD3, we randomly select τ = 3 timestamps out of the T = 12 timestamps in the sliding window; as a result, we feed the sparse supra-Laplacian (i.e., 3N × 3N ; where N is the number of nodes) into the supragraph diffusion convolutional layer. For smaller networks such as Ethereum token networks, we utilize the T timestamps (i.e., the whole sliding window information) to construct the supra-Laplacian. In the future we will explore various graph sparsification techniques, e.g., subsampling, to reduce computational costs for larger networks. Remark F.5 (More Details on the Framework of TAMP-S2GCNets in Figure 1). In Figure 1, the input is the target network over past 3 time slices ({Gt−2,Gt−1,Gt}) and our TAMP-S2GCNets consists of 5 components. That is, (i) equipped with node feature matrix Xt, we apply spatial graph convolutional layer (i.e., GCN Layer in Figure 1) on Gt extracts spatial information at time t; (ii) feature transformation (i.e., FT in Figure 1) learns representation of the spatio-temporal data XT over a sliding window of size T (in Figure 1, for better understanding the operation procedure, we
set T = 3, i.e., XT = {Xt−2, Xt−1, Xt}); (iii) we create a multiplex network based on the input networks (i.e., {Gt−2,Gt−1,Gt}) and construct the corresponding supra-Laplacian; lastly, we feed the supra-Laplacian into supragraph convolutional module (i.e., Supra GCN Layer in Figure 1); (iv) in Figure 1, we consider two types of Euler-Poincaré Surface Representation (DEPS), i.e., DEPSf,g and DEPSf,h are Euler-Poincaré surface representations under different multifiltrations, where f, h, g represent different types of filtrations; in (iv), we feed Euler-Poincaré surface representations into CNN base models and employ global average/max pooling to the outputs of CNN base models. Finally, we combine these embeddings to obtain the final embedding, where is fed into GRU module for multi-step forecasting.

heading: G ADDITIONAL COMPUTATIONAL COMPLEXITY AND TRAINING TIME COMPARISON
text: Remark G.1 (Training Time Comparison). We report the average training time (per epoch in seconds) for our TAMP-S2GCNets and 4 baselines on PeMSD4 and Decentraland datasets in Tables 16 and 17. In terms of training time, TAMP-S2GCNets runs slightly slower than the baselines but the demonstrated running time is comparable to the baselines. Although the time per epoch is slightly higher for TAMP-S2GCNets, we find that the performance of TAMP-S2GCNets is significantly better than the state-of-the-art baselines. As such, this extra time appears to be worth it.
Remark G.2 (More Details on Computational Complexity of Topological Summaries). Computational complexity of Dynamic Euler-Poincaré Surfaces (DEPS) is much lower than that of persistence diagrams. This is because Euler Characteristics is an alternating sum of the Betti numbers, and one does not need to compute the persistence diagrams to compute the Betti numbers. The computational complexity of a single persistence diagram PDk is O(N 3), where N is the number of k-simplices (Otter et al., 2017). On the other hand, by using sparse matrix methods, computational complexity for the Euler Characteristics of a simplicial complex withM simplices isO(M) (Edelsbrunner & Parsa, 2014). If p is the resolution size of the multipersistence grid, the resulting complexity for DEPS is O(p2M). However, since one does not need to compute the persistence diagrams to obtain the Betti numbers, faster algorithms are possible. For instance, Lesnick & Wright (2019) developed a much faster algorithm (RIVET software) to compute Betti numbers in multipersistence setting. This reduction approach leads to computational complexity of O(R3), where R is the size of the filtered complex giving the minimal representation for the bipersistence module and R is far less than the original number of simplices. Furthermore, by utilizing sparse matrix reductions, recent results of Kerber & Rolle (2021) significantly improved the RIVET approach. In particular, computational complexity for multipersistence landscapes with RIVET would be O(R5), where R is again the size of the reduced filtered complex (Vipond, 2020).
We have provided both computational complexity and running time of 8 topological summaries on Decentraland and Golem (see Table 18).
The computational complexity of the overall approach is: O(N2 + NT FNQsup + NT 2Qsup/2 + Qsup ∑t−1 `=t−T M (`) +R3 +WGRU ). That is, we have: (i) spatial graph convolution module: O(N2)
aDeg: node degree; Btwns: node betweenness centrality; Power: edge weight; Trns: transaction between addresses. Please refer to on types of filtrations (Appendix D.2) for a more detailed description of types of multifiltration.
bWhere p is the resolution of MP-I, Q is the number of lines, B is the maximum number of bars of the fibered barcodes.
cWhere R is the size of the reduced filtered complex (see previous paragraph).M is the number of simplices in the simplicial complex.
dThe reported time for DEPS is given with the O(p2 ×M) algorithm. eThis computational complexity is for the time after the computation of persistence diagrams.
(where N is the number of nodes); (ii) supragraph diffusion convolution module: O(NT FNQsup + NT 2Qsup/2 +Qsup ∑t−1 `=t−T M
(`)) (where FN is the number of node attribute features, T denotes the sliding window size, Qsup is the output dimension of the supragraph diffusion convolutional layer, and M (`) is the number of edges at the `-th layer); (iii) DEPS: O(R3), where R is the size of the reduced filtered complex; (iv) GRU:O(WGRU ), whereWGRU is the total number of parameters in the GRU.

Review summaries:
This paper proposes the use of multi-parameter persistence (an emergent research topic in topological data analysis) to capture latent time-conditioned relations among nodes in a GNN. To do so, the authors introduce a dynamic Euler-Poincaré surface as a new multi-parameter persistence summary, and prove its stability and empirically show its computational efficiency. A "supra" graph convolution module is then proposed in order to allow for simultaneously learning co-evolving spatial and temporal correlations in the complex multivariate time-series data. The superior predictive performance of the proposed approach is shown empirically based on 3 different multi-variate time-series datasets from 3 different domains: highway traffic flow, Ethereum token prices, and COVID-19 hospitalizations.