On Time Series Clustering with Graph Neural Networks
Authors: Jonas Berg Hansen
Andrea Cini
Filippo Maria Bianchi
Venue: Accepted by TMLR ()

Abstract: Graph clustering and pooling operators have been adopted in graph-based architectures to capture meaningful patterns in time series data by leveraging both temporal and relational structures. However, the contribution of each design choice and the behavior of different operators remain underexplored. This work introduces a streamlined deep learning framework based on a spatio-temporal graph neural network (STGNN) for clustering time series, which can leverage prior knowledge on the spatial structure of the data. The STGNN-based model flexibly identifies clusters in various data settings through an encoder-decoder architecture with a bottleneck, showing that a spatio-temporal approach can identify meaningful clusters even in datasets that do not explicitly include spatial relations. We validate the framework’s qualitative performance through experiments on synthetic and real-world data, showing its effectiveness in different scenarios. We also provide a heuristic for model selection in unsupervised settings via a self-supervised forecasting loss. Code is available at https://github.com/NGMLGroup/Time-Series-Clustering-with-GNNs

Published in Transactions on Machine Learning Research (08/2025)
On Time Series Clustering with Graph Neural Networks
Jonas Berg Hansen
wmmas.bhansen@@hi...o
Deportment of Mathematics and Statistics, UiT The Arctic University of Norway
Andrea Cini
andrea. middra.ciin@@usi.dee
Swiss National Science Foundation Postdoc Fellow
IDSIA USI-SUPSI, Università della Suizzera italiona
Filippo Maria Bianchi
filippo. m. bianchi@ust..rr
Deportment of Mathematics and Statistics, UTT The Aretic University of Norway NORCE Norwegion Research Centre AS
Reviewed on OpenReview: https: // openreview. porm? 6d=MHQX iXer3
Abstract
Graph clustering and pooling operators have been adopted in graph-based architectures to capture meaningful patterns in time series data by leveraging both temporal and relational structures. However, the contribution of each design choice and the behavior of different operators remain underexplored. This work introduces a streamlined deep learning frame- work based on spatio-temporal graph neural network (STGNN) for clustering time series, which can leverage prior knowledge on the spatial structure of the data. The STGNN-based model flexibly identifies clusters in various data settings through an encoder-decoder archi- tecture with a bottleneck, showing that a spatio-temporal approach can identify meaningful clusters even in datasots that do not explicitly include spatial relations. We validate the framework's qualitative performance through experiments on synthetic and real-world data, showing its effectiveness in different scenarios. We also provide heuristic for model lection in unsupervised settings via self-supervised forecasting loss. Code is available at https://gite thhttbb/ com hon//GGMGGroup/Timm -errsrrrooiriimm---riiessssss Clustering-wiie with-GNNa
1 Introduction
In recent years, graph-based time series processing archilectures have emerged as powerful tools for handling multivariate time series by modeling with a graph the (spatial) relationships between different sequences (Jin et al., 2024; Cini et al., 2023a). Such relationships can either be given part of the data derived from temporal characteristies of the time series (Kipf et al., 2018; Cini et al., 2023c). Clustering and graph pooling have been employed both as preprocessing techniques and as integral components of graph-based architectures to enhance performance in tasks such as forecasting (Cini et al., 2020; 2024) and missing data imputation (Marisca et al., 2024). These methods capture dependencies among time scries, enabling the learning of representations that account for multiple scales and dynamics at different resolutions. Despite the effectiveness of clustering and graph pooling components, previous work mainly focused on the performance of downstream tasks and put less emphasis on analyzing the generated clustering partitions. In this sense, a focused study about the capability of graph deep learning models to learn meaningful clusters by integrating both spatial and temporal information is still lacking.
To fill this gap, we analyze a streamlinnd graph-based model for clustering time series with sparse spatial relations represented by a graph. This model is implemented as a spatio-temporal graph neural network (STGNN) with an encodor-decoder architecture. The STGNN incorporates temporal processing units (Jin et al., 2024), message passing (MP) layers (Battaglia et al., 2018), and graph clustering (pooling) mod- ule (Grattarola et al., 2022). The latter creates bottleneck positioned between the encoder and the decoder

---

Published in Transactions on Machine Learning Research (08/2025)
of the STGNN and learns end-to-end clustering partition of the time series collection. The pooled repre- sentation is then mapped back into the original graph. The full model is trained end-to-end to perform time series forecasting as a (self-supervised) pretext task.
Given the unsupervised nature of the clustering problem, standard cross-validation procedures cannot be followed to select the optimal hyperparameterss of the STGNN model. However, we show that, with the proposed encoder-decoder architecture, the forecasting accuracy correlates well with the quality of the ex- tracted clusters. This observation allows for adopting a model selection strategy based on the self-supervised forecasting loss, thus avoiding reliance supervised information.
To evaluate the effectiveness and limitations of the STGNN model in clustering time series, we perform two types of experiments. The first consists of clustering a set of synthetic data, where both spatial ure and temporal dependencies need to be accounted for to cluster time series correctly. The second is a qualitative analysis of the clusters obtained on a real-world dataset of time scries, where the spatial structure is not available a prior. In the latter case, we show that, even when an explicit graph structure describing the relationships across the time series is not provided, an STGNN can be used to obtain meaningful clusters if the underlying graph is built based on the similarities of the time series.
We summarize our main contributions below.
We systematically analyze the chustering performance of reference STGNN model in different settings, both from a quantitative and qualitative perspective.
We show that, within our framework, forecasting loss can be used to perform model selection in a completely unsupervised fashion.
show that STGNNs can be used to learn meaningful clusters from collections of time series, even when spatial relationships are not given pre-defined.
The paper is structured as follows. Section 2 introduces the background and the problem definition. Sec- tion 3 describes the deep learning framework used to perform clustering, along with its training and model selection. In Section 4 we discuss the related work. Section describes the synthetic and real-world datasets. In Section 6, we analyze the performance and the clustering results obtained with the proposed STGNN framework compared to traditional approaches for time series and graph node clustering, deep learning ap- proaches for node clustering of static graphs, and few variants of the proposed framework. Finally, we draw our conclusions in Section 7.
2 Preliminaries
2.1 Problem definition
Consider a group of N multivariate time series with C 2 variables. We define the observations at time as X, E NXC, Furthermore, a sequence of observations from time t to t T is denoted as X E TxNNNC Analogously, we use x; E RCC to denote observations w.r.t. the i-th time series at time step t. All time series are assumed to be homogenous, i.e,, to have the same C' channels referring to the same variables. The time series can be associated with M covariates that, at time I, are described by a matrix U, NNM We assume that each time series generated by time-invariant stochastic process and the presence of Granger causality among them (Granger, 1969).
Time series clustering The goal of time series clustering is to group time series that exhibit similar temporal patterns. Standard approaches typically rely on similarity scores among time series or representa- tions thereof. Clustering assignments can be described by a (soft) assignment matrix S E R X where K denotes the number of clusters, the entry $1 is the membership of time series i to cluster j, and Σ, St = 1. Hard assignments are obtained by taking the row-wise argmax of S. The quality of a clustering partition is generally measured in terms of the compactness and separation of the clusters (Maulik & Bandyopadhyay, 2002). Information-theoretia measures be used to quantify how much the elements in cluster similar
2

---

Published in Transactions on Machine Learning Research (08/2025)
to each other and different from those of other clusters (Gokeay & Principe, 2002). Finally, when supervised information is available, the performance can be measured as the agreement between the class and the cluster labels (Rand, 1971). For example, commonly used metrics include normalized mutual information (NMI), homogencity score (HS), and completeness score (CS). The NMI measures the mutual dependence between the ground-truth labels and the cluster assignments, normalizee by their entropy. Formally, it is defined as:
21(Y:C) NMT(Y,C) H(Y) H(C)'
where Y denotes the true labels, C the prodicted cluster labels, I(Y; C) is the mutual information between Y and C, and (-) represents the entropy. The HS evaluates whether each cluster contains only data points that belong to a single class. It is given by:
HS H(YC)) H(Y)
(2)
where HY is the conditional entropy of the class labels given the cluster assignments. The CS assesses whether all data points belonging to the same class are grouped into the same cluster, and is defined as:
H(D CS = 1 H(C)
(3)
where H(C))) is the conditional entropy of the cluster assignments given the true labels. All these metrics take values in the interval [0, 1], where a value of indicates perfect clustering performance.
Relational priors In many real-world scenarios, time scries are not independent but show a rich de- pendency structure, e.g., in sensor networks, traffic data, or biological signals. Specifically, we assume the presence of Granger predictive causality (Granger, 1969) among time series, i.e., we assume that observa- tions at related time series can be used to make more accurate predictions. Further, we assume dependencies among time series to be sparse and encode them in an adjacency matrix A E RNN, The entries of A can be binary or real, in which case weight is associated with each edge. These dependencies, represented graph, can be used to inform the clustering process, c.g., as regularization. Note that we use the term spatial in a broad sense to refer to the dimension spanning the time series collection. In this regard, the adjacency matrix is used to encode functional dependencies and sparsity priors on the dependency structure.
2.2 Graph Neural Networks with Pooling
This section introduces the reference graph neural network (GNN) operators used in Section 3 as components of the proposed framework. Hicrarchical GNN architectures are usually built by alternating message passing (MP) layers, which use structural relations to update the node representations, with graph pooling layers, which gradually reduce the size of the graph and the node representations (Grattarola et al., 2022).
Graph convolutional operators As MP layer, we use as reference a simple graph convolutional operator with parametrized skip connection (Bianchi et al,, 2021). Let H E NNF be the matrix of the node features associated with the nodes in the graph. Given adjacency matrix A, the Graph Conv Skip (GCS) convolutional operator is defined as:
H' = D--1/2 AD 1/2 HOMP + HOsc + b).
(4)
where D = diag(A1), and OMP E FXF Osc E R FXF' and b R² are learnable parameters.
Graph pooling Several families of hierarchical graph pooling operators exist, which differ in how they compute the topology and features a coarsened graph (Grattarola et al., 2022; Liu et al., 2023a). In this work, we focus on pooling operators that cluster the nodes of the original graph, and then let each cluster become a node in the pooled graph. Specifically, the graph pooling operator learns a soft assignment matrix S E NNK, which maps each node of the original graph into one more of the K supernodes of the pooled
3

---

Published in Transactions on Machine Learning Research (08/2025)
graph. The assignment matrix S is typically obtained by passing node features learned by more MP layers into a learnable function such as multilayer perceptron (MLP). The structure of S is influenced by the graph topology, the node features, the loss of the downstream task, and a set of auxiliary losses that encourage the formation well-behavod clusters (Grattarola et al., 2022). The clustering- based pooling operators we consider compute the node features of the pooled graph as:
pool = S H E RxF,
(5)
where the k-th row Hoool is the combination node features assigned to supernode k.
We also consider a lifting (or unpooling) operation used to map the pooled nodes back to the original node space:
Hig Sf NX
(6)
where f()) represents processing (e.g,, the application of one or more MP layers) between the pooling and unpooling operations.
Graph clustering regularizations One of the main differences between the existing clustering- based graph pooling methods, such as MinCutPool (Bianchi ct al., 2020a), DiffPool (Ying et al., 2018), DMoN (Tsit- sulin et al., 2020), and TVGNN (Hansen & Bianchi, 2023), in how they regularize the learning of the cluster assignment matrix S; usually, this happens through the use of an auxiliary clustering loss. The mathematical formulation of the different losses is deferred to Appendix A.
While different methods exist, the loss usually consists of two components (Bianchi, 2022):
Laux Liopo + fquality.
(7)
The first component, Lrope; ensures that the cluster assignments consistent with the graph topology, while Lquality encourages clusters to be balanced and well-formed. These auxiliary losses are generally inspired by graph-theoretical properties and are typically formulated based on the static structure of the graph. However, when graph nodes represent time series, it not immediately clear whether and how strongly these regularizations impact the formation of meaningful temporal clusters. The impact of these regularizations may vary depending on how accuralely the adjacency matrix A captures relationships relevant to the clustering task. In this study, we explore whether such regularization contributes meaningfully the clustering process and investigate under which conditions incorporating these losses provides practical advantages.
2.3 Spatio-temporal Graph Neural Networks
STGNNs integrate GNNs and temporal processing units to capture dependencies across both the spatial and temporal dimensions. The core component of an STGNN is spatio-temporal MP operator (Cini et al., 2023a; Marisca et al., 2024). Let hid be the sequence of representations up to of the i-th node at the l-th layer. Following Cini el al. (2023a), generic spatio-temporal MP operator updates representations at each layer as:
h hi441 = ACGR Φ A JEN(i)
(n h24 (gh.))))))))))))) a:)):
(8)
where ays is entry (i.j) of A, M(i) is the neighborhood of the i-th node, AGGR is a differentiable and permutation-invarinn aggregation function (e.g,, sum or mean), and 26, oi are differentiable update and message functions, respectively.
STGNN architectures generally follow two main paradigms depending on the order in which temporal and spatial information processed (Gao & Ribeiro, 2022; Cini et al., 2023a). The time-then-space (TTS) paradigm first encodes temporal information and subsoquently processes the resulting node embeddings spatially. This improves computational elliciency, as spatial propagation is performed only with respect to the last time step, thus reducing computational overhead and memory consumption. Conversely, the time- and-space (T&S) paradigm interleaves temporal and spatial propagation. While the T&S approach is more flexible, the streamlined architecture of TTS models makes them valid choice in many scenarios, e.g., when scalability is requirement.
4

---

Published in Transactions on Machine Learning Research (08/2025)
INPUT
OUTPUT
A
X, T-iit Ut-Iwa
H
O2
X Ki+TT
Conditional Block
Sequence Encoder
Pool Operator
Temporal Decoder
Lift Operator
MP Layers
Decoder
R
H
Zet-To Tu
S
Encoder
Figure 1: Architectural diagram of the model.
3 Time series clustering with STGNNs
In the following, present reference architecture for time series clustering with STGNNs. In particular, propose time series clustering pipeline where the input time series are first encoded in static graph rep- resentation, which is then pooled into a smaller graph where each node represents cluster. Representations corresponding to each node are mapped into predictions and then mapped back to the original nodes. this way, the predictions for each node depend only on their cluster membership. The architecture is trained end-to-ent minimize forecasting errors at each node in a self-supervised fashion, clustering together time series with similar dynamics. The training loss is optimized jointly with a graph clustering loss, which regu- larizes the learned assignments based the predefined relational structure. As result, the graph structure influences the learned model both through the loss function and the message passing operators.
Section 3.1 introduces each processing block, while Section 3.2 discusses the training strategy, the loss function, and model sclection.
3.1 Model architecture
We start by giving a general overview of the architecture of the STGNN model. Then, we detail the structure of the encoder and decoder. diagram of the whole architecture shown in Figure
3.1.1 Overview
The proposed architecture is an STGNN with coder block of the model follows a TTS paradigm, i.c., the temporal dimension is processed first and encoded into a vector embedding for each node. The time embeddings are then processed by GNN block consisting of stacked MP layers followed by a pooling operator that results a reduced set of N nodes. More details are discussed in Section 3.1.2.
encoder-decoder structure. The en-
After pooling, the coarsened representation is processed by decoder that first maps the representation each node into different encoding for cach target forecasting time step, reintroducing a temporal dimension. We refer lo the resulting K time series as latent factors (Wang el al., 2019b). Such latent factors can be interpreted as a set of bases, each corresponding to one of the clusters that the model has learned. Since the decoder generates predictions for all the original nodes through lift operation from such latent factors, their formation directly influences the performance of the downstream task. A detailed description of the decoder is deferred to Section 3.1.3.
5

---

Published in Transactions on Machine Learning Research (08/2025)
3.1.2 Encoder
Let Tw be the length of the temporal window of the input time series. At the beginning of the encoder, a combination of two MLPs used to combine the input time series Xt-Twst and covariates Ut- as
r = CondBlock(az ut) = M=;) + MLPu
where rf E RR and H is the hidden size, i.e., the embedding dimension. If exogenous variables are absent, MLPy is not utilized.
Sequence encoder These initial representations are then processed by a sequence encoder (SeqEnc) which embeds the temporal information of each node into a static vector:
h SeqEnc(rhw:))
(10)
where h RR"; we use H. Nx to denote the matrix obtained by stacking the node representations. The sequence encoder can be implemented by any sequence modeling block, e.g., a temporal convolutional network recurrent neural network, which produces vector encoding of the entire sequence.
In practice, as the temporal encoder, we use stack of dilated temporal convolutional networks (TCNs) (Lea et al., 2016) with increasing dilation factors, to capture at each layer different input dynamics. To obtain a vectorial representation, the Tw representations generated by the TCN layers are combined with an allention-based aggregator.
Spatial message passing After the temporal information has been encoded, the spatial dependencies captured by the graph topology embedded into node representations through GNN encoder consisting of a stack of I. MP layers:
H = ,))
(11)
where H. E NxH and we use GNN( A) to indicate the stack of MP layers applicd w.r.t. to A. While, in principle, any MP operator can be considered, we use the graph convolution operator of Eq. 4 in our implementation.
Clustering step There a strong connection between node clustering and the family of graph pooling operators mentioned in Section 2.2, particularly those that optimize losses inspired by graph theoretical objectives, such graph cuts. Indeed, one use the graph pooling framework to model the assignment of nodes to supernodes, i.c., clusters. In practice, similarly to Cini et al. (2024), we rely on a static cluster assignment matrix S E N by parameterizing it directly with a table of learnable parameters 4 E N By doing so, we ensure consistency in the cluster assignments across time steps. Note that this is different from what typically done in standard pooling operators, where the cluster assignment matrix is obtained as function of input node features (Grattarola et al., 2022).
To obtain normalized assignments, apply row-wisc softmax to the unnormalized scores क The output of the GNN encoder is then pooled to obtain a reduced node set with K nodes as:
S.k
O = ST ST S.A
(12)
ΣM S.N
where O E K and Š is the rescaled cluster assignment matrix. By scaling S before pooling, the magnitude of the resulting node representations O, becomes independent of cluster sizes. This prevents the magnitude of each representation from scaling with the number of nodes it is derived from. This mitigates numerical instabilitics for very small clusters and facilitates the formation of latent factors that better capture average patterns in the clusters time step t.
3.1.3 Decoder
The docoder first processes the pooled representations with an MLP that maps cach node represcntation into a sequence of length Th: ZEATT MT.PTD(O.), Z, E Rxc
(13)
6

---

Published in Transactions on Machine Learning Research (08/2025)
The time series Zu are the latent factors, which serve as cluster representatives and are instrumental to the reconstruction of the input time series while encouraging meaningful clusters. As described in Eq. 6, the soft assignments are used to lift the features of the pooled graph back to the input node space:
X, = SZz, T t.....,,
(14)
Here, Th is the horizon of the time series forecasting used to train the model, and Xe-+TT are the final predictions. By design, the latent factors ZEETTH act as bases for the forecasts, which are obtained as convex combinations as shown in Eq. 14. In the case of sharp cluster assignments, i.e. when the rows of the cluster assignment matrix S are close to being one-hot vectors, each node is almost exclusively mapped to a single cluster. In this scenario, the lift operation directly maps the corresponding latent factor to all nodes in that cluster, which will more explicitly connect forecasting performance to the groupings made during pooling, Conversely, if the cluster assignments are smooth, different latent factors are combined to build a prediction. Therefore, enforcing sharpness in S can act as further regularization. We propose a training strategy that encourages each latent factor to serve as a distinctive (and representative) average predictor for its corresponding cluster, thereby leading to the formation of more meaningful clustors.
3.2 Training
The proposed model is trained in a self-supervised fashion to forecast the next Tи steps by minimizing a mean absolute error (MAE) loss Fpred- The optimization procedure will also directly affect the parameters in S and, thus, the formation of the clusters.
Graph-based regularization A key aspect of optimizing GNNs with pooling layers is to promote the formation of informative clusters while avoiding the degenerate solutions that often arise at the start of training. Those typically include assigning all nodes to the same cluster or memberships being near-uniform and uninformative. To this end, we consider auxiliary losses provided by the existing graph pooling opera- tors (Grattarola et al., 2022). These losses are used to modify the learnable parameters & that are used to compute the soft assignments S. The balance between the forecasting loss and the auxiliary pooling losses is key to capturing both spatial and temporal information, which, in turn, enables obtaining meaningful clustering result. To control this balance, the pooling losses are weighted with coefficients:
Liotal = Cpred + C1 Cropo + +2Fquality
(15)
where Ltopo and Lquatity are the two terms the auxiliary loss from Eq. 7 associated with the pooling method of choice. As discussed at the end of Section 3.1.3, our goal is to produce sharp assignments to obtain latent factors that represent each cluster well. To this end, during training, we gradually lower the temperature T in the softmas used to compute the soft cluster assignments S.
3.3
Model selection
As there are no universal optimal criteria. for asscssing clustering performance and measuring clustcring quality, model selection largely depends prior assumptions about the relevant structures in the data and priorities regarding the similarities that one wishes to model. For example, we might want to group time series that exhibit similar dynamics or to identify patterns that appear more often in a certain class. Class labels, however, are not available during training and validation in a completely unsupervised setting. This implies that we cannot. leverage metrics such as the NMT, HS, and CS to perform model selection, as these motrics quantify the agroement between identified clusters and class labels.
In our architecture, predictions are obtained directly from the latent factors associated with each cluster. As result, forecasting accuracy will depend on how well each latent factor represents the dynamics of the time scries in the associated cluster. As shown in Section 6, our approach highlights positive correlation between cluster-class correspondence and forecast performance. Specifically, we show a significant correlation bet ween the forecasting loss and the clustering quality evaluated in terms of NMI and HS. Under the assumption that clusters contain time scries with similar patterns and that those patterns characterize, to some extent, the different classes, we can perform model selection by heuristically monitoring the forecasting loss Epred proxy of clustering quality.
7

---

Published in Transactions on Machine Learning Research (08/2025)
4
Related Work
There is a vast literature on time series clustering methods based on standard statistical approaches (Liao, 2005; Izakian et al., 2015), machine learning methods (Aghabozorgi et al., 2015; Bianchi et al., 2020b), and deep learning approaches (Ma et al., 2019). Similarly, node clustering is a well-known problem in graph processing research (Schaeffer, 2007) that can be addressed with traditional techniques such as spectral clustering (von Luxburg, 2007), traditional deep learning architectures (Wang et al., 2019a; Pan & Kang, 2021), and modern GNNs (Bianchi et al., 2020a; Tsitsulin et. al., 2020; Bianchi, 2022). Methods to cluster time series based modeling spatio-temporal relationships have traditionally relied on density estimation techniques that leverage various forms of kernel methods (Nakaya & Yano, 2010; Hu et al., 2018; Jiang et al., 2012) and ST-DBSCAN (Birant & Kut, 2007; Gomide et al., 2011). A more recent example is the work of Deb & Karmakar (2023), which proposed to adopt a convex combination of spatial and temporal distance matrices. In this work, we mainly focus on spatio-temporal deep learning models for correlated time series.
Graph clustering and pooling in spatio-temporal graph neural networks (STGNNss Recent advancements in apatio-temporal forecasting have leveraged GNNs to handle complex dependencies in spatial and temporal data (Cini et al., 2023a). Similar to our proposed approach, spatio-temporal processing in STGNNs is typically done by alternating MP with temporal processing layers. Several works also leverage clustering techniques to pre-group data points or enhance model performance (Caggiani et al., 2017: Cini et al., 2020). However, the way clustcring is applicd in these works differs significantly from our proposed graph pooling strategy. For example, Ji el al. (2023) proposes a framework for traffic flow prediction, which leverages adaptive data augmentation and clustering- based self-supervised tasks to model spatial and temporal heterogeneity to enhance the representation of traffic patterns. To enhance traffic flow forecasting accuracy, Chen et al. (2022) leverages spatial-temporal clustering to group traffic nodes based on recent and periodic traffic states, enabling the STGNN model to capture the spatio-temporal dependencies better. Graph pooling has been increasingly adopted in STGNN architectures to learn richer representations, but mostly focusing on improving the performance on downstream tasks rather than on the quality of the clustering partition. Cini et al. (2024) and Ma et al. (2022) incorporate graph pooling operators, such as MinCutPool (Bianchi et al., 2020a) and DifiPool (Ying et al., 2018), to create hicrarchical representations of time series data. Marisca et al. (2024) also extract hierarchical representation but utilize a non-trainable k- MIS pooling operator. Other works adopted more conventional pooling operators like global max pooling (Liu et al., 2023b). To our knowledge, the closest related work to our framework is found in Brenner et al. (2023), which presents a GNN model with an encoder-decoder structure employing graph pooling applied to capacity expansion problems. However, their model is trained to perform reconstruction rather than prediction. In addition, they employ a temporal window in processing, whereas our approach integrates both spatial and temporal dynamics more comprchensively. Also, while Brenner et al. (2023) puts some focus on clustcring, their work is still mainly performance-driver and does not perform a thorough analysis of the clustering results. In addition, the cluster assignment matrix for graph pooling is computed directly from the output of a block of GNN layers rather than being learned directly, and, thus, the cluster assignments will vary with time rather than being static.
To summarize, all these works demonstrate the utility of graph pooling in learning multi-scale representations, but limited focus is put on the quality and the interpretability of the clustering results, and the extent to which spatial and temporal information is accounted for in producing the groupings.
5 Dataset description
This section presents the datasets considered our experimental evaluation.
5.1 Synthetic datasets
We have designed a series of synthetic datasets to assess the capability of the proposed STGNN architecture to leverage both temporal and spatial dependencics to form clusters. Specifically, each datasct consists of a graph with nodes associated with a time series and class label. The class of node is determined both by its role in the graph and by the associated time series, in proportion that varies across datasets. In other
8

---

Published in Transactions on Machine Learning Research (08/2025)
words, in some datasets, the class of node depends more on the graph topology, and in others it depends more on the characteristios of the time series.
The graph of each dataset is constructed by connecting subgraphs obtained from the Barabási-Albes (BA) model (Barabási & Albert, 1999; Barabási & Frangos, 2002). The nodes of each subgraph are associated with one of the classes, and the subgraph is constructed using a unique choice of parameters the BA model. The full graph is obtained by connecting all the subgraphs through a few random edges. The BA model was chosen because, compared to other models such as the stochastic block model, its edge sampling algorithm offers a reasonable amount of flexibility to generate subgraphs that differ in both edge density and structure. For example, one can obtain range of topologies spanning from more sparsely connected tree-like structures to more dense structures with heavily linked hub nodes. Each subgraph is also guaranteed to be connected (no disjoint regions), meaning that only a few edges between the subgraphs are necessary to ensure that the full graph is connected.
The time series are generated using a graph polynomial vector autoregressive (GPVAR) (Zambon & Alippi, 2022) system with local effects. The procedure is inspired by Cini et al. (2023b), with the difference that we use simusoids to model the local effects
H =
(16)
9=1
X = tan(H +a sin a pprt/ppa 6) m.
(17)
where @ E RQXL are the autoregressive coefficients, a, Φ E RN denote the sine wave parameters, and n E @N(00, denote a vector of N random values generated from a Gaussian distribution centered at zero with the variance 6² being a predefined parameter. To inject the class information within the structure of the time series, all the time series associated with nodes of class k share the set of parameters ak; Pk, and ΦR. The intra-class differences arise from the Gaussian noise n. and the structure of A, which combines each node with different neighbors.
5.2 Smart meter dataset
Real-world collections of time series with node-level labels and associated relational information are rare, which complicates computing quantitative evaluations of the partitions performance metrics such as NMI rely on supervised information. On the other hand, lime series classification datasets are common and, even if the graph topology is not part of the original inputs, an STGNN can still effectively perform time series clustering. In this case, a suitable adjacency matrix must be derived from the data in a way that reflects meaningful relational information across the time series.
Specifically, test the effectiveness of the proposed STGNN model to perform clustering in a real-world set- ting, we consider the Commission for Energy Regulation (CER) dataset (Commission for Energy Regulation (CER), 2012). This dataset contains time series of energy load consumption from over 5,000 Irish homes and businesses, which were collected between 2009 and 2010. Each time series is associated with an electricity customer labeled as "residential", "SME" (small/medie medium enterprise), or "other". Since the "other" class primarily consists of unlabeled data from the other two categories, we focus only on clustering time series from the "residential" and "SME" groups Furthermore, due to the large amount of missing values in the datasct, we filter out customers with more than 5% missing load values. And finally, to somewhat mitigate the imbalance between "residential" and "SME", we randomly sample 25% of the residential customers. These steps yield final dataset of 1,541 time series, with 70/30 split between "residential" and "SME".
Because no explicit relational dependencies are provided in the original data, we construct a graph based on the level of association between the lime series. In this instance, the prior spatial information that establishes a dependency structure between the time series is based on the assumption that only series exhibiting a degree of similarity should be considered as related. To build the graph, we first compute the similarities among the time series and then apply a k-nearest neighbor (k = 3) approach to create a sparse connectivity matrix. Then, we binarize and symmetrize the connectivity matrix, obtaining the final
9

---

Published in Transactions on Machine Learning Research (08/2025)
(a) Balanced-U
(b) Balanced-NU
(c) OnlyGraph
(d) OnlySeries
Figure 2: Visual representation of the graph of each synthetic dataset. The colors/shapes correspond to the different underlying classes. Node positions were generated using the algorithm.
force-directed
adjacency matrix A that describes an undirected graph with all edges weighted equally. In Section 6.2, we perform a study to investigate the impact of taking different approaches for building the relational graph, and we evaluate the overall impact of the spatial component of the STGNN. Additional details regarding the CER subset and the constructed graphs are given in Appendix B.2.
6 Results and discussion
In this section, we present and discuss the experimental results. Section 6.1 compares the proposed framework against several baselines and provides a qualitative cluster analysis on the syuthetic datasets; Section 6.2 investigates the impact of using different approaches to build the adjacency matrix and provides a qualitative cluster analysis on the CER dataset; Scction 6.3 reports ablation studies to assess the contribution of individual components; and finally, Section 6.4 evaluates the performance of few variants of the proposed framework.
For each dataset, hyperparameter selection was made through grid sweep where only the pooling loss coefficients were searched, with candidate values consisting of 6 linearly spaced points spanning from 0.1 to 2.5. Furthermore, each configuration was run times, and the configuration that achieved the lowest mean validation mean absolute error (MAE) was selected for evaluation.
We split the datasets along the temporal dimension, where the first 80% of the time series was used for training and the next 10% for validation. The remaining 10% kept aside for making qualitative plots. For those baselines that did not require hyperparameter selection, we also used the validation set for training,
For the synthetic datasets, the window size Tw was set to 16, and for CER it was set to 72 (hours). For both the synthotic and real-world datasets, the proposed STGNN was set to do a one-step prodiction, i.c., Th 1. Further details about model configuration and additional results are deferred Appendix C.
6.1 Synthetic data clustering
Four different datasets were created following the procedure described in Scction 5.1. Each dataset consists of graph with 1500 nodes divided into 5 classes, and there a time series of 2000 time steps associated with each node. Fig. 2 shows the graph topology and the class of the nodes in cach datasct.
In the first two datasets, named Balanced-U U and Balanced-NU, the parameters of the graph and the time series were chosen such that the class information influences both the spatial and temporal components of the datasets. The only difference between the two datascts is that in Balancod-U the cardinality of the elements in each class is the same, while in Balanced-NU the class distribution is non-uniform. In the latter two datasets, named OnlyGraph and OnlySeries, the class information characterizes either the graph topology or the time scries, respectively, but not both. In the OnlyScrics dataset, the underlying graph is a single Erdös Erdös-Rényí graph (Erdös & Rényi, 1959) rather than being a composition of different subgraphs, one for each class. On the other hand, in the GnlyGraph dataset, all series generated from the same AR(2)
10

---

Published in Transactions on Machine Learning Research (08/2025)
model. For OnlyGraph and OnlySeries, the class distribution is uniform. The details about the generations of the different synthetic datasets are reported in Appendix B.1.
As baselines, we compare against methods that either utilize only the spatial or the temporal information. We also consider as basclines architectures for static graphs where the node attributes are a set of temporal and spectral features of the time series. As purely spatial methods, we consider spectral clustering (SC), Walktrap (Pons & Latapy, 2005), and Leiden (Traag et al., 2019), which cluster the nodes by looking at the graph topology. On the other hand, the chosen temporal methods apply spectral clustering to similarity matrix computed from the time series and disregard the original graph. As time series similarity measures, we use the Euclidean similarity (obtained by taking the inverse of the distance scaled to [0,1]), the cosine similarity, and the dynamic time warping (DTW) algorithm (Berndt & Clifford, 1994). We also consider the cosine similarity between the vectorial embeddings obtained by processing the time series with an echo state network (ESN), as described in Bianchi et al. (2020b).
As discussed in Section 4, there is a lack of STGNN-based approaches specifically designed for clustering time series. To provide comparison with methods that can consider both time series data and an underlying graph structure, we extract a set of static features from the time series and apply methods for clustering graph nodes that can account for the node attributes. We rely on the TSFEL (Barandas et al., 2020) package to extract the time series features. By using the default parameters, we obtain a total of 165 features consisting of temporal, statistical, and spectral properties of the time series, To process the static graphs we consider a standard graph autoencoder (GAE) and a variational graph autoencoder (VGAE) (Kipf & Welling, 2016), which are trained to perform link reconstruction. After training, we take the latent representations generated by the encoder of the GAE and VGAE, and we cluster them using standard K-means algorithm. As an additional baseline, we consider a GNN followed by graph pooling layer, implemented as MinCutPool (Bianchi et al., 2020a), which outputs a soft assignment matrix S. This latter method is named StaticPool and is trained end-to-end by minimizing the MinCutPoo losses. The loss coofficients and softmax temperature used in StaticPool are all fixed at 1. In all these basclines (GAE, VGAE, and StaticPool), the GNN encoder is implemented according to equation 11 and, thus, is similar to the one in our proposed STGNN architecture.
For the proposed STGNN model, report the results obtained using the auxiliary pooling losses of MinCut- Pool. Also, we report results for a configuration trained only with the forecasting loss named "NoPoolLoss The results with the other pooling operators are deferred to Appendix C.2. For both the baselines and the proposed STGNN model, the number of clusters K was set equal to the number of node classes.
Numerical results displaying the cluster-class correspondence in terms of NMT, HS, and CS are displayed in Table 1. For each configuration, we report the averages and standard deviations of 5 independent runs. For the balanced datasets, the proposed model with MinCutPool outperforms the other methods in all metrics for both datasets, which shows that the STGNN manages to capture both temporal and spatial information in the data when forming the partition. Additionally, while the NoPoolLoss variant performs reasonably well on average, it fails to leverage the graph information effectively and shows greater inconsistencies in the clustering partitions oblained across the different runs. In the dataset where the class information is only contained in the time series, MinCutPool resulted in more inconsistency than NoPoolLoss. On the other hand, with only graph information, the model with MinCutPool excels, while the NoPoolLoss model fails. It should be noted, however, that in the OnlyGraph dataset, the time series do not contain features that can be leveraged to guide the model selection via the forecasting loss and, thus, the optimal hyperparameter configuration is essentially chosen at random. A more detailed discussion about this point is provided in Appendix C.2.
The performance difference between the proposed model and the methods processing the graphs with static, extracted features can be explained by the main differences in their optimization. GAE, VGAE, and Stat- icPool are optimized with loss functions that primarily focus on the structural properties of the graph, specifically link reconstruction and graph cut regularization. At the same time, these methods rely on dis- tinetive node attributes to produce meaningful latent representations and subsequent clusters. So, when both the graph and the series have a balanced set of distinctive features, they perform very well. However,
11

---

Published in Transactions on Machine Learning Research (08/2025)
Table I: Synthetic data experiment results. SC: spectral clustering, DTW: dynamic time warping, ESN: echo state network. For each metric and dataset, the highest value in bold, and the second highest is underlined. The group Static feats refers the methods processing graphs with static node features extracted from the time scries. Values in percentages (average deviation).
Method
Dataset
Balanced-U HS
Balanced-NU
NM1
CS
NM1
HS
CS
47.7
40.2
58.5+0.0
52.6+0.0
65.920.0
Spatial
Walktrap
61.7+00.
88.020.0
47.6=0.0
60.9
87.7+0.0
46.6=0.0
Leiden
77.710.0
91.820.0
67.4
62.8+0.0
86.020.0
49.4
Euclidean
92.3++.0
92.1+
92.3+ta.0
87.2+0.0
88.8+0.0
85.6#u.0
Cosine
91.8+o.0
91.6*#.0
91.9+0.0
86.8+0.0
88.5+0.0
85.3=0.0
Temporal
DTW
67.0+0.1
63.9
70,2+0.1
66.4+0.4
65.2
67. 67.7
ESN
86.3+1.5
82.4
90.6#2.5
81.241.2
77.8+0.9
84.8
GAE VGAE
94.116.2
92.0
96,4
93.4
91.5
95.4+3.2
Static feats
98.210.9
98.220.9
98.2:1..9
94.842.0
94.8+2.0
94.71220
StaticPool
94.7*4..0
94.524.2
94.943.8
69.8+4.6
72.2+4.9
67.524.2
NoPoolLoss
74.2+37.1
73.9+37.D
94.642.9
83.6#1.8
85.0+1.6
82.2
Proposed
MinCutPool
100.010.0
100.010.0
100.010.0
95.520.2
95.811.2
95.310.2
OnlySeries
Only Graph
NMI
HS
CS
NMI
HS
CS
SC
0.5+0.0
0.4=0.0
0.6
100.0#0.0
100.0z00
100.0=0.0
Spatial
Walktrap
7.940.0
14.520.0
5.420.0
85.620.0
85.620.0
85.620.0
Leiden
0,8+0.0
1.1+0.0
0.740.0
99.7+..0
99.7+0.0
99.7+0.0
Euclidean
100.0
100.020.0
100.0+0.0
0.3+0.1
0.3+n.1
0.3_0.1
Cosine
100.0#0.0
100.020.0
100.0+0.0
0.3+0.0
0.3+n.0
0.3=o.0
Temporal
DTW
99.510.0
99.5
99.5
0.210.0
0.210.0
0.2-0.0
ESN
99.410.6
99.410.6
99.4.10.6
0.310.1
0.3
0.3=0.1
GAE
8.7+2.3
8.3
9.1+2.3
21.5#9.
20.1
23.1+10.4
Static feats
VGAE
11.7+0.9
11.7+t09
11.8#t..9
41.4+s.0
41.3+8.1
41.5#8.0
StaticPool
13.742.6
13.7+2.6
13.722.6
6.3
6.3+1.
6.3=1.6
NoPoolLoss
100.0:0.0
100.020.0
100.0+0.0
0.040.0
0.0#0.0
0.0=0.0
Proposed
MinCutPool
75. 1438.8
74.74-38.6
75.5+38.6
100.0+0.0
100.0
100.0=o.a
Cluster 1
Cluster
Cluster
Cluster
Cluster
MMM o
0
100
100
100
100
100
Figure 3: Quantile plot for a snippet of the time series of each cluster for the Balanced-U dataset. The median is given as a solid line, and the corresponding latent factor from the model given as a dashed line. The shaded areas correspond to the 25% and 75% quantiles.
compared to the proposed approach, they struggle when there is less information either in the spatial or the temporal domain. See Appendix C.2 for additional details and further experiments on this setting.
The latent factors produced by the STGNN model with MinCutPool for a segment of the series are displayed in Fig. 3, together with the quantile plots of the time series in the corresponding clusters. As we can see, the latent factors closely match the cluster medians, and we argue that they can be considered valid representatives of the clusters.
12

---

Published in Transactions on Machine Learning Research (08/2025)
DiffPool
MinCutPool
DMoN
TVGNN
1.0
1.0
1.0
1.0
mean NMI
0.5
0.5
0.5
0.5
0.0
0.0
0.0
0.0
0.75
0.70
0.75
0,70
0.70
0.68
0.75
0.70
Mean Validation MAE
Mean Validation MAE
Mean Validation MAE
Mean Validation MAE
Figure 4: Plot of linear regression for mean NMI as a function of mean validation MAE for the Balanced-U dataset (the x-axis is inverted).
Unsupervised model selection As discussed in Section 3.3, we rely solely on the unsupervised forecasting loss for model selection. Fig. 4 plots the NMT against the MAF for the configurations considered during the hyperparameter sweep of the pooling coefficients. Notably, the configurations that yield low prediction MAE generally also achieve high NMI, which supports hypothesis that the forccasting performance and the quality of the clustering partition are closely related in our proposed architecture. Moreover, this correlation holds even when other losses different from those of MinCutPool are used. We underline that relying on the MAE to perform model selection is a noteworthy result, given the well-known difficulties in tuning the hyperparameters a clustering model in the absence of supervision.
6.2 Qualitative cluster analysis on real-world time series data
The main goal of this experiment is to study the quality of the partitions generated by the proposed spatio- temporal clustering approach when applied to a purely temporal dataset, i.e., where graph encoding the relationships among the time series is not given. In particular, we consider the real-world time series dataset presented in Section 5.2 and we manually construct an adjacency matrix before applying the given framework.
We consider different methods to construct the adjacency matrix. First, as described Section 5.2, we build k-nn graph by looking at the similarities between the time series. As the similarities, we consider the Euclidean distance, Pearson correlation, and correntropy (Liu et al., 2007) across the different time series. In addition, perform an ablation study by considering three adjaccncy matrices that do not account for the temporal information: Erdös-Rént Rényi graph. These latter configurations represent baselines with non-informative graph structures to explore whether more informative spatial features are beneficial for time series clustering. In each sctting, we consider the NoPoolLoss and MinCutPool variants of the proposed framework. Notably, NoPoolLoss, configured with the identity matrix as adjacency, represents a purely temporal deep learning framework optimized solely with forecasting loss.
identily matrix, fully connected matrix, and the adjacency matrix of an
We also consider two different settings by letting the number of clusters be either equal (K = 2) greater (K' = 5) than the number of classes. The latter setting allows for deeper qualitative assessment of the latent factors as the clusters are smaller, more compact, and should show more diversity. Therefore, in the 5-cluster setting, our quantitative comparisons focus exclusively on the homogeneity score, which measures the class purity within each cluster.
The results comparing the clustering performance relative to the "residential" and "SME" classes are reported in Table 2. In the 2-cluster setting, utilizing uninformative graph that disregards the relationships among time series results in significantly poor performance. Furthermore, if no pooling losses are applied, the clustering performance does not improve even when more informative graphs are utilized. On the other hand, when the Mool losses are minimized during training, the performance is noticeably enhanced when using an informative graph that describes the relationship between time series. The most significant improvement is observed when the graph is built using the Pearson similarity across the time scries. When adding the pooling losses to the configurations with uninformative graphs, the performance did, however, not improve but rather decreased. Additionally, in some of the runs with the model utilizing the fully connected
13

---

Published in Transactions on Machine Learning Research (08/2025)
Table 2: Performance on subsampled CER dataset for the different adjacency construction methods.
Method
2 clusters
5 clusters
NMT
HS
CS
HS
Identity
8.1:2.9
8.5
7.7
29.31:3.1
Fully Connected
8.3
8.7
7.912.3
28. 289.7
Random
8.4+2.1
8.8
8.0+2.0
21.1+11.2
NoPoolLoss
Euclidean
7.2+1.3
7.6+1.3
6.9+1.2
32.4
Pearson
4.9
5.2+3.3
4.743.0
27.1
Correntropy
5.822.0
6.1
5.521.9
29.424.6
Identity
6.942.3
7.3+2.5
6,5+2.2
15.1
Fully Connected
0.3
0.2+n.3
64.7243.4
6.3
Random
2.242.6
2.3+2.7
2.1+2.6
0.2+0.2
MinCutPool
Euclidean
26.8:0.9
28,010.9
25.7:0.9
58 58.61114
Pearson
45.3
44.511.
46.142.5
50.3
Correntropy
22.5+1.1
23.4+1.2
21.7+11
58.7
Cluster
Cluster 2
Cluster 3
Cluster 4
Cluster 5
0.50
1.0
7.5
0.25
0.5
2.5
50
100
50
100
50
100
50
100
0
50
100
Figure 5: Quantile plot for a portion of the CER test set for the case with 5 clusters obtained with the proposed model applied with the correntropy k-nn graph. The cluster median is given as solid line, 25% and 75% quantiles are given as shaded areas, and the corresponding latent factor is given as a dashed line.
graph, we observed collapse into a single cluster, which resulted in a high average and standard deviation in CS.
In the 5-cluster scenario, with some exceptions, the different configurations improve in terms of HS. This improvement is expected since a finer partition allows for splitting larger clusters that combine samples from different classes. However, the extent of performance enhancement in terms of HS varies significantly. When the proposed STGNN is configured with MinCutPool, the configurations utilizing Euclidean and correntropy-based graphs exhibit the greatest improvement, with correntropy achieving the highest, mean HS value overall. The NoPoolLoss variant also achieves better overall performance when using more clusters and, expected, outperforms the model with the MinCutPool losses when using the uninformative graphs. This latter result is not surprising because if the graph structure is completely unrelated to the actual time series, optimizing the clustering objective of the MinCutPPol losses has an adversarial effect. To further qualitatively analyze the clustering results, in Fig. 5 we compare the latent factors identified by the model using the correntropy-basee graph to the quantiles of the time series assigned to each cluster. As we can see, the latent factors closely match the median of the time series distribution within each group, meaning that they can effectively be used as representatives of the different clusters. In Fig. 6 we show the class medians and the class distribution within each cluster. All clusters exhibit relatively high degree of class purily (i.e., high HS value), as they mostly contain time series from the same class. This also means that the latent factors closely match the energy consumption profiles for the most frequent class in each cluster. We also note that in clusters like the fifth one, there is a notable resemblance between the median profiles of the two classes.
It is also worth noting that the difference magnitude between the time series significantly affects the clustering results. This is direct consoquence of how we normalized the data, i.e., we standardized each time series with the global mean and standard deviation computed across the whole dataset. To put more
14

---

Published in Transactions on Machine Learning Research (08/2025)
Cluster I
Cluster
500
500
Value
Count
Value
1.0
Count
0.2
250
250
0.1
0.5
0
50
100
R
S
50
100
R
S
Time
Time
Cluster
Cluster
500
500
Value
Count
Value
Count
250
250
0
100
R
100
Time
Time
Cluster
500
Value
Count
250
0.5
50
100
R
8
Time
Figure 6: Class distribution within the same clusters in Fig. 5. For cach cluster, the line plot shows the median for each class with "residential" represented as dashed blue line and "SME" solid orange line. The bar plot to the right of each line plot shows the class distribution with R: "residential" and S: "SME"
Identity
Euclidean
Pearson
Correntropy
1.0
1.0
1.0
1.0
mean NMI
0.5
0.5
0.5
0.5
0.0
0.0
0.0
0.0
0.38
0.36
0.37
0.36
0.35
0.38
0.36 Mean Validation MAE
0.37
0.36
0.35
Mean Validation MAE
Mean Validation MAE
Mean Validation MAE
(a)
Identity
Euclidean
Pearson
Correntropy
1.0
1.0
1.0
mean HS
0.5
0.5
0.5
0.5
0.0
0.0
0.0
0.375
0.350
0.325
0.35
0.30
0.38
0.36
0.35
0.30
Mean Validation MAE
Mean Validation MAE
Mean Validation MAE
Mean Validation MAE
(b)
Figure Linear regression plot for CER (the x-axis is inverted) with (a): 2 clusters and (b): clusters.
emphasis on the shape rather than on the magnitude, one should normalize each time series by own mean and standard deviation.
Unsupervised model selection As the synthetic dataset experiment, Fig. 7 we display the relation- ship between NMI HS and the validation MAE the settings where the graph is built using the time series similarities (Fuclidean, Pearson, and Correntropy). For comparison, we also include the setting where the
15

---

Published in Transactions on Machine Learning Research (08/2025)
Table 3: Ablation results based on hyperparameter search over the loss coefficients. For each configuration, the first three columns contain the NMI, HS, and val MAE obtained by selecting the model setup with the lowest val MAE (based on mean values from 3 runs). The next two columns give the highest observed NMI and HS values in the scarch. The last two columns present correlation metrics between NMI and val MAE, where Tp: Pearson correlation coefficient and Spearman's rank correlation coefficient.
Config
Obtained from search
Best NMI and HS
NMI-MAE
NMI
HS
MAE
NMI
HS
Tp
Ts
AB1
85.2+2.1
86.2:1.9
0.757=0.001
88 88.221.2
88. 18.
-0.959
-0.968
AB2
95.1+0.2
95.4=0.1
0.68511..001
95.4=0.1
95.1+0.2
-0.922
-0.886
AB3
51.0+33.8
52.7+3.9
+0.00s
57.0+28.1
57.2+-312
-0.583
-0.570
AB4
75.2
77.6#6.8
0..01
75.2+7.3
77.6#66.
-0.743
-0.638
Proposed
95.520.2
95.820.2
0.68220.000
95.520.2
95.820.2
-0.915
-0.863
graph is represented by the identity matrix, which corresponds to configuration where graph information is given. Also in this case, there is pronounced correlation between the two metrics, which empirically supports our hypothesis of leveraging MAE for doing model selection in the absence of supervision.
The low dispersion of the points around the linear fit indicates that there is a strong relationship between NMI/HS and MAE, no matter how the graph is constructed. However, the numerical values achieved by the different graph construction methods cannot be compared directly. For example, the configuration that achieves the lowest mean MAF does not necessarily imply the best NMI, as evidenced by the comparison between Euclidean and Pearson in Fig. 7a.
6.3 Additional ablation studies
The experiments on the synthetic data and CER have shown the effect of different design choices, including testing different loss functions and coefficient values, utilizing different adjacency construction methods, and training with only the forecasting loss. In this section, we report further ablation experiments aimed at studying different designs of the sequence encoder and different types of scheduling for the softmax temperature. The results are limited to the Balanced-NU dataset, which is the more complex of the two balanced datasets. The tested configurations are as follows:
AB1: A configuration where the TCN in the encoder is removed and basic mean aggregator is used to reduce the temporal dimension before message passing.
A configuration that uses a mean aggregator after TCN instead of an attention-based aggre- gator.
AB2:
AB3 and AB4: Configurations with static softmax temperature during training, fixed at 1.0 and 0.01, respectively.
The results of each ablation configuration are presented in Table 3. The performance of AB1 and AB2 shows that utilizing deep learning based sequence encoder is beneficial to the clustering, but in this case, there is no significant difference between using a TCN with mean or attention-based aggregation. AB4 and AB5 show the benefit the temperature scheduling in the softmax used to produce the cluster assignments and subsequently perform pooling. In particular, fixing the temperature 1 results in significantly worse cluster-class correspondence and worse NMI-MAE correlation.
6.4 Additional configuration experiments
Here, we report additional experiments whose goal is investigate further the importance and impact of the chosen sequence encoder, MP layers, and temporal decoder by replacing them with different components. In particular, consider the following variants:
16

---

Published in Transactions on Machine Learning Research (08/2025)
Table 4: Additional experimental results.
The experimental setup is the same as for Table
Config
Obtained from search
Best NMI and HS
NMI-MAE corr
NMT
HS
MAE
NMT
HS
Tp
Ts
AAN-STGNNEBCCON
95.7+0.1
96.0=0.1
0.68311.000
95.710.1
96.010.1
-0.979
-0.993
EAT-steGNEEccoer
95.1
95.5=0.1
0.683
95.310.1
95.610.1
-0.940
-0.974
EucDecSkip
72.5+0.5
72.9=a.5
0.563+40.001
73.3
73.7
-0.838
-0.759
EmbeediingDecoder
14.6+0.9
15.2=0.9
0.714++0.004
86.9+
87.2+10.7
-0.303
-0.272
Proposed
95.5+0.2
95.820.2
0.68210.000
95.5+0.2
95.8
-0.915
-0.863
RNN-STGNNEnooder: In this variant, the TCN is replaced with a GRU RNN (Cho et al., 2014) with the same hidden size and number of layers, where the contraction the temporal dimension is accomplished by using only the last state of the RNN. While we have already established that using a deep learning based sequence encoding improves the performance, it is important to check if using RNNs improve both forecasting and clustering performance. Indeed, RNNs outperform TCNs in some cases and learn more impactful latent embeddings.
This variant replaces the GCS layers in the encoder with GAT lay- ers (Veličković et al., 2018). The purpose of testing this configuration is to ascertain to what extent adopting a more sophisticated attention-based MP approach influences the clustering outcome.
GAT-STCNNEncoder:
EncDecSkip: In this variant, we add a skip connection that bypasses the bottleneck between the encoder and docoder. Specifically, in addition being passed to the pooling operator, the embed- dings from the MP layers are passed to different MLP that reintroduces the temporal dimension. The output of this second decoder is added to the output of the lift operator to form the final prediction. This variant investigates the role and the importance of the bottleneck by incorporating a computational path that bypasses it entirely.
EmbbeddingDeeoder:: This is a variant that relains the bottleneck without any bypass option, but increases the capacity of the components processing the pooled representations in the decoder. Spocifically, the number of temporal channels the latent factors is increased from to 8, and, after the lift operation, a learned embedding with dimensions N 8 is used to weight the different channels and produce final prediction for each node. The rationale behind this architecture to test whether having different predictions for each time series (obtained by using local node embeddings) would improve or degrade clustering performance.
The experimental setup is otherwise the same as with the additional ablations of Section 6.3, and the results are reported in Table Starting with the variant using an RNN sequence encoder (RNN-SIGNNEnnoder)) we see comparable clustering performance and an improvement in the correlation metrics. This latter im- provement indicates more monotonic relationship between NMI and MAE, and a possible higher guarantee that basing model selection on MAE performance will result in higher cluster quality. However, as we already mentioned, the resulting cluster-class correspondence is not significantly improved, and the RNN- based model requires a 35% longer training time on the same machine with an RTX 4090 GPU. So, for this particular dataset, we argue that the TCN-based model is preferable. Replacing GCS with GAT (GAT- STGNNEncoder) does not significantly change the clustering performance, although the correlation between NMI and MAE becomes stronger, while keeping the computational complexity comparable with GCS.
Regarding EncDecSkip, introducing skip connection to bypass the bottleneck resulted in improved forecast- ing performance, which is expected. However, the clustering performance is degraded, and the NMI-MAE correlation is weaker. On the other hand, for Embedding and additional local embeddings did not enhance the forecasting performance. In addition to the lack of improved predictive performance, this variant produced significant detrimental effect on the NMI-MAE correlation, resulting in a model selection that yields poor clustering performance. In particular, selecting the model with the lowest MAF also resulted in much lower NMT and HS than the highest observed in the
the added capacity via more channels
17

---

Published in Transactions on Machine Learning Research (08/2025)
hyperparameter search, further illustrating that the model selection process becomes unreliable. In both of the latter two variants, plausible reason for the weakened correlation is that the clustering and forecasting objectives become less aligned when the architecture becomes more suitable for forecasting.
7 Conclusions
In this work, we have introduced a framework for clustering time series using spatio-temporal graph neural networks (STGNNs). Our findings highlight the importance of considering not only the temporal features of the individual time series but also the spatial axis describing the relat ionships across the time series to find meaningful clustering partitions. By employing an encoder-decoder architecture with intermediate and interpretable latent factors, the proposed framework adjusts the underlying data structures and manages to form meaningful clusters in different scenarios.
Through experimentation with both synthetic and real-world datasets, we have validated the robustness and flexibility of our model. A series of quantitative and qualitative analyses underscores the model's capability to leverage structural and temporal data, thereby forming informative clusters that show good correspondence with the class labels.
The optimization of the loss coefficients vital component in learning a meaningful clustering partition. particular, several of the existing pooling losses promote balanced partition, and adjusting the weights of the different components of the loss can lead to the formation of clusters with non-uniform size. However, finding the optimal hyperparameter configurations without relying supervised cross-validation is well- known challenge. To address this issue, we showed that the self-supervised forecasting loss function can act as an efficient tool for guiding the clustering partition. This important connection allows us to perform model selection in a completely unsupervised setting.
Clustering spatio-temporal data using STGNNs remains a largely unexplored area, and our study represents an important initial step toward systematically addressing this gap. Through extensive quantitative and qualitative analyses, we have demonstrated the fcasibility and potential of STGNN-basod clustcring. To continue the research in this direction, there is an urgent need for standardized benchmarking frameworks and comprehensive datasets to more effectively evaluate and compare emerging methodologies and STGNN architectures. While further research and additional real-world datasets are essential, we believe our proposed evaluation framework and the introduced synthetic dataset provide a solid foundation, facilitating future advancements and inspiring continued exploration in this promising research direction.
Acknowledgments
This work was supported by the Norwegian Research Council grant no. 345017 (RELAY: Relational Deep Learning for Energy Analytics). We wish to thank Nvidia Corporation for donating some of the GPUs used in this project.
References
Saeed Aghabozorgi, Ali Seyed Shirkhorshidi, and Teh Ying Wah. Time-series clustering a decade review. Information systems, 53:16-3: 38, 2015.
Albert-László Barabási and Réka Albert. Emergence of scaling in random networks. science, 286(5439): 509-512, 1999.
Albert-László Barabási and Jennifer Frangos. Linked: the new science of networks science of networks. Basic books, 2002.
Marília Barandas, Duarte Folgado, Letícia Fernandes, Sara Santos, Mariana Abreu, Patrícia Bota, Hui Liu, Tanja Schultz, and Hugo Gamboa. Tsfel: Time series feature extraction library. SoftwareX, 11:100456, 2020.
18

---

Published in Transactions on Machine Learning Research (08/2025)
Peter W. Battaglia, Jessica B. Hamrick, Victor Bapst, Alvaro Sanchez-Gonzale Vinicius Zambaldi, Mateusz Malinowski, Andrea Tacchetti, David Raposo, Adam Santoro, Ryan Faulkner, et al. Relational inductive biases, deep learning, and graph networks. arXic preprint 0rXiv:8866.01661 2018.
Donald J. Berndt and James Clifford. Using dynamic time warping to find patterns in time series. In Proceedings of the 3rd international conference on knowledge discovery and dala mining, pp. 359-377 370, 1994.
Filippo Maria Bianchi. Simplifying clustering with graph neural networks. arv preprint orXiv:2007.0877) 2022.
Filippo Maria Bianchi, Daniele Graltarola, and Cesare Alippi. Spectral clustering with graph neural networks for graph pooling. In International Conference on Machine Learning, pp. 874-883. PMLR, 2020a.
Filippo Maria Bianchi, Simone Scardapane, Sigurd Løkse, and Robert Jenssen. Reservoir computing ap- proaches for representation and classification of multivariate time series. IEEE transactions on neural networks and learning systems, 32(5):2169-211 2020b.
Filippo Maria Bianchi, Daniele Grattarola, Lorenzo Livi, and Cesare Alippi. Graph neural networks with convolutional arma filters. IEEE transactions on pattern analysis and machine intelligence, 44(7):3496 3507, 2021.
Derya Birant and Alp Kut. St-dbscan: An algorithm for clustering engineering, 60(1):208 221, 2007.
data. Data & knowledge
Aron Brenner, Rahman Khorramfar, and Saurabh Amin. Learning spatio-temporal aggregations for large- scale capacity expansion problems. In Proceedings of the ACM/IEEE 14th International Conference on Cyber-Physical Systems (with CPS-loT Week 2023), pp. 68- 77, 2023.
Leonardo Caggiani, Michele Ottomanelli, Rosalia Camporeale, and Mario Binctti. Spatio-temporal clustering and forecasting method for free-foating bike sharing systems. In Advances in Systems Science: Proceedings of the International Conference on Systems Science 2016 (ICSS 2016) 19, pp. 244 254. Springer, 2017.
Yan Chen, Tian Shu, Xiaokang Zhou, Xuzhe Zheng, Akira Kawai, Kaoru Fueda, Zheng Yan, Wei Liang, Kevin, and Kai Wang. Graph attention network with spatial-temporal clustering for traffic flow fore- casting in intelligent transportation system. IEEE Transactions on Intelligent Transportation Systems, 24 (8):8727-8B8) 2022.
Kyunghyun Cho, Bart Van Merriënboer, Dzmitry Bahdanau, and Yoshua Bengio. On the properties of neural machine translation: Encoder-decoder approaches. arXiv preprint 4rXiv:1409.1259 2014.
Andrea Cini, Slobodan Lukovic, and Cesare Alippi. Cluster-based aggregate load forecasting with deep neural networks. In 2020 International Joint Conference on Neural Networks (LJCNN), pp. 1-8. IEEE, 2020.
Andrea Cini, Ivan Marisca, Daniele Zambon, and Cesare Alippi. Graph deep learning for time series fore- casting. arXiv preprint orXiv:2810.15978 2023a.
Andrea Cini, Ivan Marisca, Danicle Zambon, and Cesare Alippi. Taming local effects in graph-basod spa- tiotemporal forecasting. In Advances Neural Information Processing Systems, volume 36, pp. 55375 55393. Curran Associates, Inc., 2023b.
Andrea Cini, Daniele Zambon, and Cesare Alippi. Sparse graph learning from spatiotemporal time series. Journal of Machine Learning Research, 24(242):1 36, 2023c.
Andrea Cini, Danilo Mandic, and Cesare Alippi. Graph- based time series clustering for end-to-end hierar- chical forecasting. In International Conference on Machine Learning. PMLR, 2024.
19

---

Published in Transactions on Machine Learning Research (08/2025)
Commission for Energy Regulation (CER). Cer smart metering project electricity customer behaviour trial, 2009-2010 [dataset]. Trish Social Science Dala Archive, SN: 0012-00, 2012. URL https://www. ucd. ie/
Soudeep Deb and Sayar Karmakar. A novel spatio-temporal clustering algorithm with applications on covid- 19 data from the united states. Compulational Slalistics & Dala Analysis, 188:107810, 2023.
Paul Erdös and Alfréd Rényi. On random graphs i. In Publicationes Mathematicne Debrecen, volume 6, pp. 290, 1959.
Jianfei Gao and Bruno Ribeiro. On the equivalence between temporal and static equivariant graph repre- sentations. In International Conference on Machine Learning, pp. 7052-7076. PMLR, 2022.
Erhan Gokcay and Jose C. Principe. Information theoretic clustering. IEEE transactions on pottern analysis and machine intelligence, 24(2):158 171, 2002.
Janaína Gomide, Adriano Veloso, Wagner Meira Jr, Virgilio Almeida, Fabrício Benevenuto, Fernanda Ferraz, and Mauro Teixeira. Dengue surveillance based on a computational model of spatio-temporal locality of twitter. In Proceedings of the 3rd international web science conference, pp. 1-8, 2011.
Clive W.J. Granger. Invest investigatigg causal relations by econometric models and cross-spectral methods. Econo- metrica: journal of the Econometric Society, pp.
1969.
Danicle Grattarola, Danicle Zambon, Filippo Maria Bianchi, and Cesare Alippi. Understanding pooling in graph neural networks. IEEE Transactions on Neural Networks and Learning Systems, 2022.
Jonas Berg Hansen and Filippo Maria Bianchi. Total variation graph neural networks. In Proceedings of the 40th international conference on Machine learning. ACM, 2023.
Yujic HU, Fahui Wang, Cocile Guin, and Haojic Zhu. spatio-temporal kernel density estimation framework for predictive crime hotspot mapping and evaluation. Applied geography, 99:89-97, 2018.
Hesam Izakian, Witold Pedrycz, and Iqbal Jamal. Fuzzy clustering of time series data using dynamic time warping distance. Engineering Applications of Artificial Intelligence, 39:235 244, 2015.
Jiahao Ji, Jingyuan Wang, Chao Huang, Junjic Wu, Boren Xu, Zhenhe Wu, Junbo Zhang, and Yu Zheng. self-super learning for traffic flow prediction. In Proceedings of the AAAI conference on artificial intelligence, volume 37, pp. 4356-4364, 2023.
Shan Jiang, Joseph Ferreira Jr, and Marta C. Gonzalez. Discovering urban spatial-temporal structure from human activity patterns. In Proceedings of the ACM SIGKDD international workshop on urban computing, pp. 95-102, 2012.
Ming Jin, Huan Yee Koh, Qingsong Wen, Danicle Zambon, Cesare Alippi, Geoffrey I. Webb, Irwin King, and Shirui Pan. A survey on graph neural networks for time series: Forecasting, classification, imputation, and anomaly detection. IEEE Thansactions Pattern Analysis and Machine Intelligence, 2024.
Thomas Kipf, Ethan Fetaya, Kuan-Chieh Wang, Max Welling, and Richard Zemel. Neural relational inference for interacting systems. In International conference machine learning, pp. 2688-2699 Pmlr, 2018.
Thomas N. Kipf and Max Welling. Variational graph auto-encoders. ar preprint ar 16l::1661.0008, 2016.
Colin Lea, Rene Vidal, Austin Reiter, and Gregory D. Hager. Temporal convolutional networks: A uni- fied approach to action segmentation. In Computer Vision-ECCV 2016 Workshops: Amsterdam, The Netherlands, October 8-10 and 15-16, 2016, Proceedings, Part III 14, pp. 47-54. Springer, 2016.
T. Warren Liao. Clustering of time series data survey. Patlern recognition, 38(11):1857-1535 2005.
20

---

Published in Transactions on Machine Learning Research (08/2025)
Chuang Liu, Yibing Zhan, Jia Wu, Chang Li, Bo Du, Wenbin Hu, Tongliang Liu, and Dacheng Tao. Graph pooling for graph neural networks: progress, challenges, and opportunities. In Proceedings of the Thirly- Second International Joint Conference on Artificial Intelligence, IJCAI '23, 2023a. ISBN 978-1-956792- 03-4. doi: 10.24963 224463/jjcai../ 20c3 URL https: Stopss:/en /ttp::///ooooo.../o/ooo 24963/ijcal 2023/752.
Weifeng Liu, Puskal P. Pokharel, and Jose C Principe. Correntropy: Properties and applications in non- gaussian signal processing. IEEE Transactions on signal processing, :5(11) 5298, 2007.
Zhening Liu, Li He, Liang Yuan, Kai Lv, Runhao Zhong, and Yaohua Chen. Stagp: Spatio-temporal adaptive graph pooling network for pedestrian trajectory prediction. IEEE Robotics and Automation Letters, 2023b,
Qianli Ma, Jiawei Zheng, Sen Li, and Gary W Cottrell. Learning representations for time series clustering. Advances in neural information processing systems, 32, 2019.
Yihong Ma, Patrick Gerard, Yijun Tian, Zhichun Guo, and Nitesh V. Chawla. Hierarchical spatio-temporal graph neural networks for pandemic forecasting. Tn Proceedings of the 31st A ACM International Conference on Information & Knowledge Management, pp. 1481-1490, 2022.
Ivan Marisca, Cesare Alippi, and Filippo Maria Bianchi. Graph-based forecasting with missing data through spatiotemporal downsampling. In International Conference on Machine Learning. PMLR, 2024.
Uijwal Maulik and Sanghamitra Bandyopadhyay. Performance evaluation of some clustering algorithms and validity indices. IEEE Transactions on pattern analysis and machine intelligence, 24(12): 1650-1644, 2002.
Tomoki Nakaya and Keiji Yano. Visualising crime clusters a space-time cube: An exploratory data- analysis approach using space-time kernel density estimation and scan statistics. Transactions in GIS, 14 (3):223-22) 239, 2010.
Erlin Pan and Zhao Kang. Multi-view contrastive graph clustering. Advances in neural information processing systems, 34:2148 2159, 2021.
Pascal Pons and Malthieu Latapy. Computing communities in large networks using random walks. In Com- puter and Information Sciences-ISCIS 2005: 20th International Symposium, Istanbul, Turkey, October 26-28, 2005. Proceedings 20, pp. 284 293. Springer, 2005.
William M. Rand. Objective criteria for the evaluation of clustering methods. Journal of the American Statistical association, 66(336):846 850, 1971.
Satu Elisa Schaeffer. Graph clustering. Computer science review, 1(1):27-64, -4, 2007.
Vincent A. Traag, Ludo Waltman, and Nees J. van Eck. From louvain to leiden: guaranteeing well-connected communities. Scientific Reports, 9:5233, 2019.
Anton Tsitsulin, John Palowitch, Bryan Perozzi, and Emmanuel Müller. Graph clustering with graph neural networks. arXiv preprinl arXiv:2006. 16904, 2020.
Petar Veličkovié, Guillem Cucurull, Arantxa Casanova, Adriana Romero, Pictro Liò, and Yoshua Bengio. Graph attention networks. In International Conference on Learning Representations, 2018.
Clrike von Luxburg. A tutorial spectral clustering. Statistics and computing, 17:395 416, 2007.
Chun Wang, Shirui Pan, Ruiqi Hu, Guodong Long, Jing Jiang, and Chenggi Zhang. Auributed graph clustering: A deep attentional embedding approach. arXiv preprint arXiv:1906.063: 2019a.
Yuyang Wang, Alex Smola, Danielle Maddix, Jan Gasthaus, Dean Foster, and Tim Januschowski. Deep factors for forecasting. In International conference on machine learning, pp. 0607-6617 PMLR, 2019b.
Zhitao Ying, Jiaxuan You, Christopher Morris, Xiang Ren, Will Hamilton, and Jure Leskovec. Hicrarchi- cal graph representation learning with differentiable pooling. Advances in neural informalion processing systems, 31, 2018.
Daniele Zambon and Cesare Alippi. wzwhitnnsss test: test for signal uncorrelation on spatio-temporal graphs. Advances in Neural Information Processing Systems, 35:11975- 11966, 2022.
21

---

Published in Transactions on Machine Learning Research (08/2025)
A Additional background material
Table summarizes the components of the pooling losses considered in this work. The implementation of DMoN in Pytorch Geometric also uses the orthogonality loss of MinCutPool, but we consider only the modularity and collapse regularization losses that are proposed in the original paper. We refer the interested readers to the reference papers for a more in-depth explanation about the different pooling losses.
Table Pooling loss components.
Method
Topelegy-oriented loss
Quality-orintted loss
Tr(s" AS)
S
Ic
MinCutPool (Bianchi et al., 2020a)
Tr(S DS)
|| Sll vc F
DiffPool (Ying et al., 2018)
A softmax (S)softmax(()) II
1 ΣN H(S.)S
N
VC
DMoN (Tsitsulin et al., 2020)
Tr(S" BS)
|> C 1|
n
- 1) -(K
S -uqquuutt k (S.
TVGNN (Hansen & Bianchi, 2023)
2E Ensr Acy1Sen SS,al
N(K 1)
B
Dataset details
In this section, additional details about the applied datasets are given.
(a) MostlyGraph
(b) MostlySeries
Figure 8: Visual representation of the graphs of MostlyGraph and MostlyScries. The colors/shapes corre- spond to the different underlying classes. Node positions were generated using the Seanteemma--Riingeee force-directed algorithm.
B.1 Synthetic data
In addition the four synthetic datasets presented in Section 6.1, two more variants were constructed: MosthyGraph and MosthySeries. As the names suggest, these two datasets incorporate more class-related information into either the series or the graph compared to the two balanced datasets (Balanced-U and Balanced-NU). The graphs of MostlySeries and MostlyGraph are depicted in Fig. 8.
The parameters utilized to generate the six datasets are presented Table 6. Each dataset comprises 1500 nodes and 2000 time steps. In the table, Caistr represents the class distribution, m and ma denote the BA parameters, O, a, Ø, p, and are the series parameters presented in Eqs. 16 and 17. For the OnlyGraph dataset, 0 contains the AR parameters of the shared AR(2) model employed to generate the series. For
DNoMPooling.
22

---

Published in Transactions on Machine Learning Research (08/2025)
OnlySeries, the parameter prob denotes the probability of connection between two nodes in the Erdös- Rényi graph, and connected True indicates that after all initial edges are sampled, additional edges are randomly added until the graph is connected.
Table 6: Paramcters for the synthetic datasets
Balanced-U
Balanced-NU
MostlyGraph
Cauas: Uniform
Caiar [.4, -,,1,-15,-11
Caar: Uniform
m
[2, 4, 2
Rest the same as Balanced-U
m [2, 2,3, 1]
mg
2,
mg 3,
.6
0
c: Same Balanced-U
2
[.4, 3, 6, 4,.4,
a [.5,.5, .5.
[.00, 15, -., 20, .20]
Φ: [.2, 2, 2, 2,
[12, 14, 10, 16, 12]
p: [12, 12, 12, 12, 12]
.6 MostlySeries
.6
OnlyGraph
OnlySeries
Caias: Uniform
Catar:: Uniform
Ciaar:: Uniform
m: 1,111
m 2,2,
prob .005 connected True
mo 2,2,2,2,2
mo 3,
1
0 (.5,-3
@: Same MestlySeries
.3
2
[8, 8, 8888888 .8]
N/A
[1.2, 1.5, 8, 2.0, 1.4] Φ [[0, -2, 220] [12, 1,,8, 16, 20) a=
6 [.0, 15, -15 200
N/A
76 12, 14, 8, 16, 20
N/A
B.2 Real world load data
The sub-sampled and filtered CER dataset consists of 1,541 nodes, where 1,056 "residential" and 485 are "SME". The time series of each node consists of 12,865 values after an hourly resampling procedure. The constructed graphs using Euclidean distance, Pearson correlation, correntropy, and Erdös-Rényí are displayed in Fig. 9.
(a) Euclidean
(b) Pearson
Correntropy
(d) Random
Figure Visual representation of the graphs computed for the CER dataset when using four different adjacency construction methods. The colors/shapes correspond the different underlying classes: blue round nodes denote "residential" and squared orange nodes the "SME" classes, respectively. Node positions were generated using the Fruchtermaan-Reingol. force-directed algorithm.
C Experimental details and additional results
This section presents the model and training parameters used in the experiments and some additional results.
23

---

Published in Transactions on Machine Learning Research (08/2025)
Table 7: Pooling loss coefficients for the model with MinCutPool obtained from hyperparameter searches. Cj: topology oriented, C2: regularization oriented.
Synthetic
CER
Dataset
5 clusters
Adjacency
2 clusters
clusters
02
C2
el
C2
Balanced-U
0.10
0.10
Identity
1.06
0.10
2.02
0.10
Balanced-NU
1.06
0.10
Full
2.50
0.58
1.06
2.50
MostlySeries
0.10
0.10
Random
0.10
0.10
0.10
1.54
MostlyGraph
0.58
0.10
Euclidean
2.02
0.10
1.06
0.10
OnlySeries
0.10
0.10
Pearson
2.50
0.10
0.58
0.10
OnlyGraph
1.54
0.58
Correntropy
2.50
0.10
2.50
0.58
C.1 Model and training parameters
Model configuration In each layer of the model, the hidden size H fixed at 16. The MLP components of the conditional block comprise a single hidden layer and an output layer, both with the same size H and ReLU activation functions. The bias term is omitted from the output layer of the ML. processing the covariates. The TCN used as the temporal encoder consists of two temporal layers, with the dilation of the second layer set to 2. The ReLU function is used the activation in cach layer of the TCN. The hidden and output sizes are set to H, and the kernel size of the temporal convolution is set to 3. The TCN is also configured with a skip connection between the temporal layers. Finally, the temporal dimension is aggregated using a linear attention layer. The GNN encoder comprises two message passing layers, implemented as described in Eq. 4, with hidden and output sizes H. The temporal decoder is an MLP with a single hidden layer of size and an output size Tи =1
Training procedure In all experiments and configurations, the model was trained for 250 epochs with Adam using multi-step learning rate scheduler reducing the learning rate by factor of 0.5 every 50 epochs, with the initial learning rate set to 0.001. The optimizer was also configured with weight decay regularization (I.2 penalty) set to 0.0001.
The softmax temperature used when computing the cluster assignment matrix was initially set to 1.0 and then decreased linearly by 0.0099 per epoch down to a minimum 0.01.
Gradient clipping was applied during backpropagation with clip value of 5.
Hyperparameter search In all experiments, grid search conducted to find suitable values for the pooling loss coefficients where Cı, c2 E {0.1, 0.58, 1.06, 1.54, 2.02, 2.5) (six linearly spaced points between 0.1 and 2.5). The resulting coefficients from the searches are reported in Table 7
C.2 Synthetic data experiment
For the synthetic datasets, the proposed model was trained with a batch size of 16. The model was not given covariates, so the CondBlock of Eq. 9 was comprised of a single MLP layer with an output size of 16 and ReLU activation.
The obtained numeric results for MostlyGraph and MostlySeries are presented in Table The results are aligned with those presented Section 6.1 for the other dataset configurations, where the proposed model better handles the combination of spatial and temporal features, even compared to the methods applied to the graph with extracted features, as they rely more heavily the presence of relevant features in both the graph structure and node attributes to produce meaningful clusters.
The results obtained for the different pooling losses are reported in Table 9. With the exception of Diff- Pool, the different losses produce similar results. TVGNN differs significantly from the top performers for OnlyGraph. However, for this dataset, the time scries do, by intent, not contain any distinctive features help guide the clustering. All distinctive features are in the graph, and therefore, the correlation between forecasting performance and cluster-class correspondence does not hold, as be seen in Fig. 10. The model
24

---

Published in Transactions on Machine Learning Research (08/2025)
Table 8: Additional synthetic data experiment results.
Method
Dataset
MostlySeries HS
MostlyGraph HS
NMI
CS
NMI
CS
SC
16,910.0
13.6=0.0
22.410.0
81.510.0
80.5=0.0
82.610.0
Spatial
Walktrap
34.610.0
65.1=0.0
23.6
74.610.0
90,9=0.0
63.210.0
Leiden
32.240.0
52.3=0.0
23.2+0.0
90.8+0.0
99.540.0
83.5+0.0
Euclidean
100.0=0.0
100.0
100.0+0.0
43.522.0
42.5=0.0
44.7
Cosine
100.0=0.0
100.0
100.010.0
49.110.11
49.0=0.1
49.240.1
Temporal
DTW
74.2+0.0
71.4=0.0
77.2
23.7+0.1
22.6=0.1
25.1
ESN
99.8+0.2
99.8=0.2
99.8=0.2
6.241.4
1.8
6.841.5
GAE 95.240.9
95.2=0.9
0.
39.8
34.8=3.8
46,6144.
Static feats
VGAE
92,3+-14
92.2+144
92.3+1.4
47.0+3.8
44.7=5.4
49.842.5
StaticPool
84.6++4.T
84.3=0.0
85.0+4.4
39.7++2.2
39.6=2.1
39.8+2.3
NoPoolLoss MinCutPool
100.0=0.0
100.0:0.0
100.010.0
46.412.8
45.3=2.7
47.542.
Proposed
100.0:0.0
100.0:0.0
100.010.0
86.8
86.8-6.7
86.9
Table 9: Synthetic data results for different pooling methods.
Method
Dataset
Balanced-U HS
Balanced-NU HS
OnlySeries HS
NMT
CS
NMT
CS
NMT
CS
NoPoolLoss
74.2437.1
73.9237.0
94.612.9
83.611..
85.021.5
82.221.2
100.020.0
100.020.0
100.0=0.0
DifiPool
64.01332.2
59.8230.7
69.213344
58.6 28.6229..2
58.7128.0
61.7130.9
97.315.4
96.616.9
98.113.8
MinCutPool
100.0:0.0
100.0:0.0
100.0+0.0
95.5+0.2
95.8+0.2
95.3++.22
75.1
74.7238.6
75.5=38.6
DMoN
96.5+a.8
96.546.8
96.5+6.8
86.741.
88.3+1.4
85.141.6
100.0#0.0
100.0=0.
100.0=0.0
TVGNN
93.3=8.2
93.318.2
93.318.2
99.9
99.940.1
99.110.1
100.0#0.0
100.0=0.0
100.0=0.0
OnlyGraph HS
MostlySeries HS
MostlyGraph HS
NMT
CS
NMT
CS
NMI
CS
NoPoolLoss
0.0=0.0
0.020.0
0.020.0
100.020.0
100.0:10.0
100.0:0:0
46.412.8
45.512.7
47.512.
DifiPool
0.4:0..1
0.410.1
0.410.1
93.245.7
89.716..4
97.413.5
17.013.8
16.813.8
17.213.8
MinCutPool
100.0:0.0
100.020.0
100.0+0.0
100.0=0.6
100.0#0.0
100.0:0.0
86.8+6.7
86.846.7
86.946.7
DMoN
100.010.0
100.0200
100.040.0
100.0=0.0
100.0:0.0
100.020.0
13.9
93.9+10.8
93.9-10.8
TVGNN
1.340.3
1.3+0.3
1.3+0.3
100.0=0.0
100.0:0.0
100.020.0
98.5422
98.512.4
98.5+2.3
DiftPool
MinCutPool
DMoN
TVGNN
1.0
mean NMI
0.5
0.5
0.5
0.5
0.0
eat
5.0
2.5
3
2
0
2
1o-5-7.985e-
1e-647.9851e-1
7.98-17 Mean Validation MAE
10-6-7.98516
Mean Validation MAE
Mean Validation MAE
Mean Validation MAE
Figure 10: NMI-MAE correlation for OnlyGraph.
selection based on forecasting MAE thus ineffective, as expected. It is also noteworthy that the different losses exhibit significantly different behaviors within the selected scarch space for the loss coefficients, with DMoN consistently delivering less erratic outcomes than the other two graph-cut related pooling losses, and DiffPool consistently underperforms.
25

---

Published in Transactions on Machine Learning Research (08/2025)
C.3 Real world data experiment
There are a few differences in the training selup and model configuration used in the real-world data exper- iment. Firstly, a set of exogenous variables, Uta+T, based on sinusoidal functions with daily, weekly, and yearly periods was utilized to help the model better capture the scasonal patterns in the data. Secondly, the batch size was set to &
Lastly, in the lift operation of the model, the softmax with temperature is replaced with a hardmax using the straight-through trick. More specifically, during forward propagation, the hardmax is applied, while the softmax with temperature is applied during backpropagation. This change was made due to the presence of multiple outliers, which produced low validation MAE when near-hard cluster assignment matrix, which had also collapsed into a single cluster, was combined with high-maguitude latent factors corresponding to the near-zero assignment values. This phenomenon is illustrated in Fig. 11.
0.6
1.0
30
Cluster
Assignment value
Mean NMI
0.4
Magniuuee
Cluster 2
20
Cluster
Cluster
0.2
10
0.0
0.0
0.35
0.30
0.25
20
50 Node
100
Mean Validation MAB
Time
(a)
(b)
(c)
Figure 11: The effect of not using the straight-through trick on CER with Pearson adjacency. (a): NMI- MAE correlation plot, (b): Latent factor snippet of each cluster for the rightmost outlier in (a), (c): Soft assignments for the first 100 nodes for the same outlier.
26


Reviews:
[{"review_id": "zkRibNoYWq", "reviewer": "Reviewer_gKd4", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "Strengths:\n\n1.) The fact that time series can be clustered using both spatial and temporal information in a way that is useful for downstream processing is useful when access to labeled data is limited.\n\n2.) The proposed architecture is a natural generalization of GNNs to spatiotemporal data.\n\nWeaknesses:\n\n1.) Regarding comparison with other methods, it seems odd to compare only with methods that use either spatial or temporal information, but not both.  It seems intuitively obvious to me that having access to both types of data can improve clustering quality.  The performed analysis leaves open the possibility that higher-quality clusters could be produced with a simpler architecture using hand-engineered time series features.\n\n2.)  I do not see much discussion of the reasoning behind the choice of synthetic dataset.  For example, why is the Barabasi-Albert model a good one for the structure of each cluster?  Would the empirical results be the same if, for example, each cluster were an Erdos-Renyi graph with sufficiently large edge probability compared to inter-class probabilities (so the entire graph would be a stochastic block model)?   Potentially, competing methods might perform better if one changes the intra-class graph structure distribution.", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "No", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "This paper proposes a spatiotemporal graph neural network (ST-GNN) architecture and training method for clustering time series in an unsupervised fashion.  That is, the method produces clusters without respect to a downstream task.  The authors show that metrics of clustering quality on their produced clusters correlate well with performance of models trained on these clusters on downstream forecasting tasks.  The authors compare clustering metrics (NMI, HS, CS) of their method with those of various other spatial and temporal time series clustering methods in the literature.", "claims_and_evidence_yn": "No", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "Both weaknesses 1 and 2 should be addressed in some fashion.  I feel that both are important for a comprehensive evaluation.  While the results are promising, I do not quite feel that they are comprehensive enough.", "broader_impact_concerns": "N/A", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "Xd0ziOBNdn", "reviewer": "Reviewer_VMxE", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "**Strengths**\n\n- The paper is well-written and easy to follow.\n- This work is largely diagnostic and asks questions about capabilities of the existing STGNN paradigm, which is a valuable contribution besides designing new SOTA models.\n- This work conducts comprehensive studies and experiments that generally support the author's conclusions, especially the connection between forecasting MAE and NMI.\n- Code is provided.\n\n**Weaknesses**\n\n- This work is mostly focused on a specific proposed STGNN architecture. This makes the conclusions of this work largely specific to their model, and the results would be considerably stronger if other STGNN models were evaluated.\n- The central claims of this work are relatively weak. For example, it is well-expected that better forecasting performance should lead to better clustering performance. It is also expected that spatial/temporal information should be heavily exploited when both are provided per the methods' OnlySeries and OnlyGraph performance.\n- The real-world dataset validation is rather weak with the max NMI with 2 clusters being less than 50. This is also on a relatively small-size dataset.\n- The choice to use sinusoids to model local effects is unjustified. Results would be more robust if other types of noise were benchmarked.", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "This work investigates the effectiveness of STGNNs for time series clustering, discussing the usefulness of spatial relationships and evaluating STGNNs on tasks with and without an explicit graph structure. This work also evaluates the roles and effects of forecasting and topological loss functions commonly used in temporal graph learning.", "claims_and_evidence_yn": "Yes", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "See weaknesses.", "broader_impact_concerns": "N/A", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "d3pCWJRR4u", "reviewer": "Reviewer_jB4K", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "**Strengths**\n\n- The paper provides empirical evidence that forecasting MAE can serve as a proxy for clustering quality and offers a practical solution for hyperparameter tuning without ground truth labels.\n\n-  The synthetic datasets vary the importance of spatial vs temporal information, thereby enabling controlled evaluation of when graph structure helps clustering performance.\n\n**Weaknesses**\n\n- The actual contribution amounts to applying existing graph pooling methods to time series data. The paper provides no new insights into *how* or *why* these methods are effective for spatio-temporal clustering.\n\n- Only one real-world dataset (CER smart meter) is used, and crucially, this dataset has no natural spatial relationships. The authors artificially construct a graph using similarity measures, which contradicts their claim in the abstract: *\"can leverage prior knowledge on the spatial structure of the data\"*\n\n- Section 3.1 fails to provide clear reasoning or motivation for why the proposed model uses pooling operators in this specific configuration.\n\n- The authors provide no ablation studies on their architectural choices, leaving key design decisions unjustified.", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "- The paper demonstrates how existing graph pooling methods (MinCutPool, DiffPool, DMoN, TVGNN) can be incorporated into a spatio-temporal GNN architecture to cluster multivariate time series.\nThe authors show that spatial graph structure can improve clustering performance even when artificially constructed from temporal similarities.\n\n- The authors empirically show that forecasting loss correlates with clustering quality metrics in their experimental settings.\n\n- The paper introduces a set of synthetic datasets where class membership is determined by varying combinations of graph topology and temporal dynamics. The authors provide controlled test cases for evaluating spatio-temporal clustering methods.", "claims_and_evidence_yn": "Yes", "claims_explanation": "", "audience_yn": "No", "audience_explanation": "", "requested_changes": "- Q1. The authors mention on page 7 that \"during training, we gradually lower the temperature $\\tau$ in the softmax used to compute the soft cluster assignments $S$,\" which suggests the use of Gumbel softmax in Eq.(12), similar to Cini et al. (2024). However, the authors write as if they are not using Gumbel softmax. Could you clarify whether Gumbel softmax is actually being used? \n\n- Q2. The authors claim in the introduction that \"the forecasting accuracy correlates well with the quality of the extracted clusters.\" However, Figure 10 shows that this correlation does not appear for  OnlyGraph dataset. This requires further analysis, or the claim in the introduction should be revised to reflect this limitation. \n\n- Q3. The claim that prediction loss correlates with clustering quality for unsupervised model selection relies solely on empirical evidence. Can the authors provide theoretical justification or cite related work that supports why this correlation should hold in general? \n\n- Q4. Can the authors provide ablation studies on their architectural choices to justify key design decisions?\n\n---\n\n- Minor comment: The explanation of $\\eta_t$ in Eq.(17) is not sufficiently clear and should be expanded for better readability.", "broader_impact_concerns": "N/A", "additional_comments": "", "extra_scores": {}, "extra_text": {}}]

Author rebuttal:

Decision: Accept (poster)