Structured Matrix Basis for Multivariate Time Series Forecasting with Interpretable Dynamics

Authors: xiaodan  chen
xiucheng  li
xinyang  chen
zhijun  li
Venue: NeurIPS_cc_2024_Conference
Published: 2024-05-02 10:17:56

Abstract: Multivariate time series forecasting is of central importance in modern intelligent decision systems. The dynamics of multivariate time series are jointly characterized by temporal dependencies and spatial correlations. Hence, it is equally important to build the forecasting models from both perspectives. The real-world multivariate time series data often presents spatial correlations that show structures and evolve dynamically. To capture such dynamic spatial structures, the existing forecasting approaches often rely on a two-stage learning process (learning dynamic series representations and then generating spatial structures), which is sensitive to the small time-window input data and has high variance. To address this, we propose a novel forecasting model with a structured matrix basis. At its core is a dynamic spatial structure generation function whose output space is well-constrained and the generated structures have lower variance, meanwhile, it is more expressive and can offer interpretable dynamics. This is achieved via a novel structured parameterization and imposing structure regularization on the matrix basis. The resulting forecasting model can achieve up to $8.5\%$ improvements over the existing methods on six benchmark datasets, and meanwhile, it enables us to gain insights into the dynamics of underlying systems.

Main paper: 
 introduction: 
 Multivariate time series forecasting plays a pivotal role in a wide range of fields, such as traffic flow management, electricity consumption, and weather prediction. Multivariate time series data records quantities of interest from N series spanning over T time steps, and its underlying dynamics are jointly characterized by the temporal correlations (intra-series dependencies) and spatial structures (inter-series dependencies). Inspired by the advancements in Natural Language Processing, substantial research has been proposed to apply RNNs and Transformers to capture the underlying temporal dependencies [41, 32, 44, 45, 30] . Besides, the convolution paradigm also exhibits promising temporal correlation modeling capability and excels in long-term forecasting [33, 36, 29] .As the underlying dynamics are jointly described by the intra-and inter-correlations, it is equally important to explore the spatial structures for an ideal forecasting model design. Many proposals employ dense connection to capture the spatial correlations implicitly [31, 24] . However, the dense connection lacks clear structures and is prone to introduce noise from uncorrelated spatial dimensions. The development of Graph Neural Networks (GNNs) [18, 17, 7] offers an effective solution to model non-Euclidean structure data. DCRNN [22] builds graphs based on spatial proximity and conducts graph convolution to capture spatial correlations for traffic forecasting. To apply GNNs to more general scenarios where graph structures are unavailable, the forecasting methods propose to learn the graph adaptively through learnable node embeddings [37, 1, 38, 15] , which significantly enhances the forecasting performance.Despite the progress achieved, the spatial structures remain static across time steps for the aforementioned forecasting methods, which may not reflect the actual inter-series correlation. Because in many scenarios, the spatial correlations are also changing dynamically, for example, the traffic speeds of certain road segments manifest correlation only in peak hours. To relax this restriction, many dynamic graph-based methods have been proposed [43, 39, 42, 28, 34] to learn the spatial structures within a short time window. The intuition is that the spatial correlations in real-world applications often evolve continuously and tend to be stable over a short period of time. The general idea is to learn each series a dynamic representation via a nonlinear transformation f dynm by taking as input the current time window data, and it then generates the spatial structures by pairwise interacting the dynamic representations via a transformation f pair , the composition of two transformations forms the spatial structure generation function f spatial = f pair • f dynm . The dynamic representation function f dynm is often implemented as MLP ( [43, 42, 28, 34] ) or RNN ( [39] ) whereas f pair is mostly instantiated by attention mechanism or inner product. One severe issue of these methods is that the output space of f dynm is not well constrained and unbounded, this makes the learned dynamic representations very sensitive to the short time window input data and the unboundedness will be exaggerated by the inner product operation in f pair , which will lead to the outputs of f spatial fluctuate drastically and have high variance. The issue becomes even more severe in the presence of anomaly patterns. To reduce the variance, TPGNN [26] first learns a static spatial structure A (adjacency matrix) and then generates the dynamic spatial structures with a matrix polynomial M m=1 α m A m where α m is determined by the timestamps of the current time window. However, it has two drawbacks: 1) the coefficient α m solely depends on the timestamps and cannot adapt to the current window data, and 2) the matrix power basis A, A 2 , . . . , A M is overly restricted and has limited expressive capability. In addition, the existing forecasting methods often lack interpretable dynamics.In this paper, we propose a dynamic multivariate time series forecasting model with a structured matrix basis. Instead of relying on the two-stage spatial structure learning process (learning dynamic series representations and then generating spatial structures), we directly parameterize f spatial with a learnable matrix basis B 1 , B 2 , . . . , B M and represent any spatial structure with a convex combination M m=1 α m B m , α m ≥ 0 and M m=1 α m = 1. To learn the matrix basis effectively, we propose a novel structured matrix parameterization method and impose structure regularization on the basis to enhance parameter efficiency and reduce complexity. In contrast to the two-stage spatial structure learning methods, the output space of our proposed f spatial is well constrained. Consequently, the generated spatial structures have lower variance and the resulting model is easier to learn. In comparison to TPGNN, our matrix basis is more expressive since it is not limited by the matrix power constraint; the coefficient α can also be computed adaptively via the interaction of current time-window data and the basis. In addition, the coefficient α offers a fashion to track the spatial structure evolution and enables us to gain insights into the underlying dynamics. Thus, the resulting model is more interpretable.In summary, our proposed f spatial has the following appealing properties: 1) lower variance and easier to learn, 2) it is more expressive, and 3) it can yield more interpretable results. This is achieved through a novel structured matrix parameterization and structure regularization. By integrating f spatial into the forecasting framework, we evaluate the efficacy of the proposed method on six benchmark datasets, it achieves up to 8.5% improvements over existing forecasting methods across various prediction lengths and can also offer interpretable dynamics. 
related work: 
 Temporal Dependency Modeling Early deep sequential methods adopt recurrent neural networks to capture nonlinear temporal dynamics [41, 32, 5] . Motivated by the wide receptive fields of attention mechanism, various Transformer-based methods have been developed to capture the longterm temporal dependencies in forecasting. To reduce the quadratic complexity of the vanilla attention, LogTrans [21] , Informer [44] , Autoformer [35] , and FEDformer [45] have been proposed successively. Non-stationary Transformer [27] attempts to mitigate the difficulties caused by nonstationarity in modeling temporal correlation. PatchTST [30] explores the strategies of patch-level semantic modeling and channel-independence. VQ-MTM [13] explores the well-defined semantic units for the Transformer architecture in time series modeling. In state-space models, the transition matrix is employed to model long-term temporal dependencies, and their recent representative works include Hippo [10] , LSSL [11] . To reduce the computational complexity, S4 [12] and Mamba [9] have been proposed successively. The convolution-based methods have also shown promising results in time series modeling. MICN [33] explores isometric convolution to capture non-local temporal patterns. TimesNet [36] reshapes 1D signals into 2D by aligning according to inherent multiple periodicity and it then employs the 2D kernels to capture both intraperiod-and interperiod-variations. ModernTCN [29] utilizes large kernels to model long-term dependencies. In addition, the Fourier basis parameterization has also been proposed to model long-term dependencies in time series imputation [25] .The spatial correlation plays an equally important role in time series modeling. DeepAR [31] and Pyraformer [24] employ dense connection to model the spatial dependencies. However, the dense connection fails to explore the underlying structure and may introduce noises from unrelated series. The advancement of GNNs offers an effective way to model non-Euclidean structure data. DCRNN [22] constructs the graph by using spatial proximity and proposes to fuse the spatial information via graph convolution operation, and thus it is only applicable when the underlying graph structures are easily accessible. To sidestep this limitation, the adaptive GNNs-based methods [1, 38, 15, 14, 16] propose to learn each series a representation and then generate the spatial structure via the interaction of the representations. BiTGraph [4] further develops the method to account for missing patterns in message passing to handle the time series with missing values. However, these approaches implicitly assume the series representations are shared across the entire time steps, and hence the underlying spatial graphs remain static over time.In real-world applications, the inter-series correlations or spatial structures are often evolving dynamically. To adapt to these scenarios, many dynamic spatial structure methods are proposed, which share a similar two-stage spatial structure learning process as the static adaptive GNNs-based methods (as discussed in Section 1). The difference is that the series representations are generated by conditioning a small time window rather than the entire timeline data, and hence, the underlying spatial structures can change dynamically over time. GMAN [43] , iTransformer [28] , Crossformer [42] , and Card [34] implement the dynamic node representation function via MLP, whereas ESG [39] adopts the RNN. As the output spaces of their f spatial are not well constrained, the learned spatial structures are very sensitive to the change of time-window data and hence have high variance. To reduce the variance, TPGNN [26] proposes to represent the dynamic graphs with matrix polynomial. As the polynomial coefficients are solely determined by the timestamps, such a method fails to utilize the current time-window data. Moreover, the matrix power basis is overly restricted and has weak expressive capability. The proposals [40, 8] attempt to build dynamic graphs by merging temporal and spatial dimensions. Nevertheless, the entanglement of temporal correlations and spatial dependencies makes the models hard to optimize. In addition, the existing multivariate time series forecasting methods also lack interpretability regarding the underlying dynamics.Question: Does the paper discuss the limitations of the work performed by the authors?Answer: [Yes]Justification: See Appendix F.Guidelines:• The answer NA means that the paper has no limitation while the answer No means that the paper has limitations, but those are not discussed in the paper. • The authors are encouraged to create a separate "Limitations" section in their paper.• The paper should point out any strong assumptions and how robust the results are to violations of these assumptions (e.g., independence assumptions, noiseless settings, model well-specification, asymptotic approximations only holding locally). The authors should reflect on how these assumptions might be violated in practice and what the implications would be. • The authors should reflect on the scope of the claims made, e.g., if the approach was only tested on a few datasets or with a few runs. In general, empirical results often depend on implicit assumptions, which should be articulated. • The authors should reflect on the factors that influence the performance of the approach.For example, a facial recognition algorithm may perform poorly when image resolution is low or images are taken in low lighting. Or a speech-to-text system might not be used reliably to provide closed captions for online lectures because it fails to handle technical jargon. • The authors should discuss the computational efficiency of the proposed algorithms and how they scale with dataset size. • If applicable, the authors should discuss possible limitations of their approach to address problems of privacy and fairness. • While the authors might fear that complete honesty about limitations might be used by reviewers as grounds for rejection, a worse outcome might be that reviewers discover limitations that aren't acknowledged in the paper. The authors should use their best judgment and recognize that individual actions in favor of transparency play an important role in developing norms that preserve the integrity of the community. Reviewers will be specifically instructed to not penalize honesty concerning limitations. 
methodology: 
 Notation From a generative perspective, the multivariate time series X ∈ R N ×T ×D records a D-dimensional physical quantities of interest generated by N series (i.e., sensors or instances) over T time steps. We use X (n) ∈ R T ×D to represent the observations from n-th sensor and X t ∈ R N ×D to indicate the observations at the t-th timestamp. The slice notation X t-H:t ∈ R N ×H×D denotes the values within a window spanning from the time interval [t-H, t). The operator diag : R N ×N → R N takes the diagonal elements of a square matrix and returns it as a vector, the operator vec reshapes a matrix or tensor into a vector. i . The Multi-Scale TCN captures the temporal dependencies by performing multi-scale temporal convolution operation in a channel-independent manner, we choose the kernel sizes 1 × 2, 1 × 3, 1 × 6, and 1 × 7 in this paper. The Dynamic GCN module comprises two functions, namely, the spatial structure generation function f spatial and graph convolution function f gcn . Our proposed f spatial generates the dynamic spatial structure A t ∈ R N ×N (adjacency matrix) t by performing the graph convolution operation. Note that the spatial and temporal dimensions are kept unchanged during the entire transformation in the L blocks.The core layer of our proposed forecasting method is the spatial structure generation function f spatial . It serves to infer the optimal graph structure A t that best characterizes the present spatial correlations from the intermediate representation Z ′ t (we drop the layer index ℓ to keep the notation uncluttered in this subsection). As mentioned in Section 1, the existing methods all adopt the two-stage learning process, which results in unconstrained output function space and high graph structure variance. To address this, we propose to directly parameterize a learnable matrix basis B = {B m } M m=1 of dimension M where B m ∈ R N ×N and define the spatial structure generation function asEQUATIONHere, we choose the convex combination instead of the linear combination to better control the output space. This is reasonable since the basis B is free to optimize in the training stage. The intuition behind Eq. 1 is that the basis B can be shared and optimized globally across different time windows, and to infer the spatial structure dynamically, we only need to adaptively compute the coefficient α t ∈ R M by conditioning on Z ′ t . However, two challenges remain in adopting the spatial structure generation function in Eq. 1. 1) The number of parameters to be learned equals M N 2 , which grows quadratically to N . As these M matrices lack connection and constraints, learning by brute force will sooner become infeasible even for medium-size N . This actually is the reason why the existing methods resort to the two-stage learning process, i.e., learning each series an embedding whose learnable parameters are N D with D being the embedding dimension. 2) Intuitively, the best α m should be simultaneously determined by Z ′ and B m , but identifying B m with N 2 parameters will make α m hard to compute when N is large.To circumvent the two challenges, we propose to parameterize the basis matrices in a structured manner and impose additional structure regularization on the basis. The idea is to represent each basis matrix B m in its SVD (singular value decomposition) factor product form B m = U m Σ m V ⊤ m and then parameterize the factors U m , Σ m , V m , where U m , V m ∈ R N ×N are orthogonal matrices, and Σ m is a diagonal matrix consisting of the singular values of B m . One benefit of such structured parameterization is that it permits us to establish connections between the basis matrices and impose constraints, and consequently, enhance the parameter efficiency and ease of the model learning. To be specific, we impose the constraint that all B m for m = 1, 2, . . . , M share the same parameterized orthogonal matrices U, V and each matrix has its unique Σ m , i.e.,EQUATIONThe rationality behind such a choice stems from the geometry interpretation of orthogonal matrixvector multiplication, i.e., left multiplying a vector by an orthogonal matrix is equivalent to coordinate transformation, which is provided in Appendix A.1. Hence, we implicitly require that all basis matrices share the same pair of coordinate transformations, which can be considered a sort of implicit regularization since it reduces extra freedom and guides the model to find the coordinate transformations (U and V) that best suit the basis.Dynamic Coefficient Generation Such a parameter sharing regularization mechanism along with the structured parameterization also brings another benefit, it allows us to treat diag(Σ m ) ∈ R N as a fingerprint to identify each B m . Hence, we can compute the α t to infer the dynamic spatial structure by simultaneously conditioning on Z ′ t and diag(Σ m ). To this end, we design an adaptive matching module that takes as input Z ′ t ∈ R N ×Do and [diag(Σ 1 ), diag(Σ 2 ), . . . , diag(Σ M )], and it yields the coefficient α t at the t time step, as follows.EQUATIONα t = softmax (W d D) ⊤ W z z/ √ d ∈ R Mwhere TimeEncoding is the timestamp encoding function, W d and W z are used to match the dimension, i.e., mapping D and z into the R d space. Given the structured basis in Eq. 2, each α t represents a dynamic spatial structure or graph at step t asEQUATIONGiven A t at each time step, we can perform graph convolution operation to aggregate information from the spatial dimension. The process is illustrated in Figure 1 -(b). The output space of f spatial is well constrained on the premise that Σ m are well bounded, which is given by the following theorem. The proof is presented in Appendix A.2 and Theorem 3.1 states that the variance of the learned structures is controllable by restricting the maximum value of Σ m , which is easy to achieve since Σ m can simply be parameterized by a vector with nonnegative values. Besides, by tracking the change of α t over time, our proposed method enables us to gain insight into the underlying dynamics of the system, and thus offers additional interpretability, as we will show in Section 4.4.Orthogonality To impose the orthogonal constraint, one may attempt to apply an orthogonality penalty, i.e., by adding the penalty term ∥U ⊤ U -I∥ + ∥V ⊤ V -I∥ to the optimized objective. However, such a hard penalty cannot guarantee genuine orthogonality and the extra penalty term also increases the learning difficulty. Hence, rather than relying on the hard penalty, we opt for the orthogonal parameterization. In particular, we restrict our attention to the special orthogonal groupEQUATIONwhich are flexible enough to represent the coordinate transformation. A nice property of SO(N ) is that it is both a compact Lie group and a smooth manifold [6] . The tangent space of the Lie group at the identity forms a vector space equipped with a Lie bracket operation, namely, its Lie algebra. The Lie algebra of SO(N ) is the set of all skew-symmetric matrices, denoted by so(N ), i.e.,EQUATIONIn differentiable manifold and Lie group, a well-known result is that the matrix exponential exp : so → SO establishes the connection between a Lie group and its Lie algebra, i.e., for any g ∈ so(N ) we have exp(g) ∈ SO(N ). The matrix exponential is defined asEQUATIONThe computation of the matrix exponential map is costly but for the special orthogonal group it has a cheap first order approximation, also known as the Clay map [3, 2] ,EQUATIONThe Clay map can be implemented in parallel by the Gaussian elimination algorithm in a numerically stable way.Low Rank Approximation By using the Clay map, the number of parameters required to parameterize U and V is N (N -1) (two skew-symmetric matrices N (N -1)/2 + N (N -1)/2). Thus the total number of parameters required by the basis in Eq. 2 is N (N -1) + M N . To further reduce the parameter count, we can apply the low rank approximation. In real-world applications, the spatial structures (adjacency matrices) of multivariate time series are often low-ranked. Suppose the rank is K, we can only preserve the first K columns of U and V as well as parameterize each Σ m with a length K vector, which leads to a parameter count upper bound N K + M K. Since K is often much less than N in practice, the low rank approximation can enhance both the parameter and computation efficiency when N is large.By stacking L blocks, the Sumba significantly enhances its capability to model temporal correlations and spatial dependencies effectively. We initializeZ (0) with the original input X ∈ R N ×H×D . The output of L-th block Z (L) ∈ R N ×H×D (L)produces the multi-step prediction Ŷt0:t0+F through a linear transformation. The model is optimized by minimizing mean absolute error (MAE):EQUATIONQuestion: For each theoretical result, does the paper provide the full set of assumptions and a complete (and correct) proof?Answer Guidelines:• The answer NA means that there is no societal impact of the work performed.• If the authors answer NA or No, they should explain why their work has no societal impact or why the paper does not address societal impact. • Examples of negative societal impacts include potential malicious or unintended uses (e.g., disinformation, generating fake profiles, surveillance), fairness considerations (e.g., deployment of technologies that could make decisions that unfairly impact specific groups), privacy considerations, and security considerations. • The conference expects that many papers will be foundational research and not tied to particular applications, let alone deployments. However, if there is a direct path to any negative applications, the authors should point it out. For example, it is legitimate to point out that an improvement in the quality of generative models could be used to generate deepfakes for disinformation. On the other hand, it is not needed to point out that a generic algorithm for optimizing neural could enable people to train models that generate Deepfakes faster. • The authors should consider possible harms that could arise when the technology is being used as intended and functioning correctly, harms that could arise when the technology is being used as intended but gives incorrect results, and harms following from (intentional or unintentional) misuse of the technology. • If there are negative societal impacts, the authors could also discuss possible mitigation strategies (e.g., gated release of models, providing defenses in addition to attacks, mechanisms for monitoring misuse, mechanisms to monitor how a system learns from feedback over time, improving the efficiency and accessibility of ML). 
experiments: 
 In this section, we evaluate our approach Sumba against 15 time series forecasting methods on six benchmark datasets (Section 4.2 and Appendix C). The ablation studies are presented in Section 4.3. We demonstrate the interpretability of our method with case studies in Section 4.4. The sensitivity of hyperparameters is provided in Section 4.5 and Appendix D; the computational cost is empirically studied in Appendix E. The code of Sumba is available at: https://github.com/chenxiaodanhit/Sumba/.Datasets We conduct experiments on six commonly adopted public datasets including: (1) Electricity [44] contains hourly electricity consumption of 321 clients. (2) Weather [35] includes 21 meteorological factors collected every 10 minutes from the weather station of the Max Planck Biogeochemistry Institute. (3) PEMS [23] records traffic data of 358 variates in California sampled every 5 minutes. (4) ETTh2 [44] contains hourly data from 7 electricity transformers. ( 5) Traffic [35] measures the hourly road occupancy rates of 862 sensors on San Francisco Bay area freeways. ( 6) Solar-Energy [20] records the solar power production, which is sampled every 10 minutes from 137 PV plants.Baselines We compare our method with the following baselines: ( 0.324 ± 0.002 0.452 ± 0.002 2.921 ± 0.010 0.0821 ± 0.001 0.220 ± 0.002 1.533 ± 0.022 S4 [12] 0.331 ± 0.002 0.479 ± 0.002 2.792 ± 0.015 0.0921 ± 0.001 0.292 ± 0.003 4.782 ± 0.055 MTGNN [38] 0.153 ± 0.001 0.245 ± 0.002 1.340 ± 0.004 0.0693 ± 0.000 0.219 ± 0.001 1.624 ± 0.018 MegaCRN [15] 0.272 ± 0.001 0.394 ± 0.002 2.429 ± 0.026 0.0674 ± 0.000 0.215 ± 0.001 1.458 ± 0.032 iTransformer [28] 0.158 ± 0.001 0.248 ± 0.000 1.472 ± 0.007 0.0644 ± 0.000 0.216 ± 0.001 1.628 ± 0.026 CrossFormer [42] 0.156 ± 0.000 0.246 ± 0.001 1.370 ± 0.008 0.0643 ± 0.002 0.210 ± 0.003 1.417 ± 0.035 Card [34] 0.152 ± 0.001 0.244 ± 0.001 1.460 ± 0.013 0.0641 ± 0.000 0.215 ± 0.001 1.523 ± 0.007 ESG [39] 0.153 ± 0.002 0.247 ± 0.001 1.434 ± 0.005 0.0753 ± 0.000 0.215 ± 0.001 1.837 ± 0.072 FourierGNN [40] 0 0.191 ± 0.002 0.286 ± 0.001 1.719 ± 0.005 0.0934 ± 0.000 0.241 ± 0.001 1.712 ± 0.020 PatchTST [30] 0.188 ± 0.001 0.290 ± 0.001 1.547 ± 0.009 0.0746 ± 0.000 0.239 ± 0.001 3.138 ± 0.096 FEDformer [45] 0.271 ± 0.002 0.384 ± 0.001 2.541 ± 0.008 0.230 ± 0.001 0.364 ± 0.003 4.118 ± 0.027 Autoformer [35] 0.269 ± 0.004 0.378 ± 0.002 2.643 ± 0.007 0.246 ± 0.002 0.399 ± 0.003 12.782 ± 1.039 Reformer [19] 0.348 ± 0.001 0.483 ± 0.001 3.189 ± 0.006 0.117 ± 0.001 0.263 ± 0.002 3.611 ± 0.025 S4 [12] 0.341 ± 0.003 0.489 ± 0.002 3.265 ± 0.009 0.0926 ± 0.001 0.285 ± 0.003 3.932 ± 0.037 MTGNN [38] 0.169 ± 0.001 0.268 ± 0.003 1.636 ± 0.003 0.0817 ± 0.000 0.256 ± 0.001 1.987 ± 0.017 MegaCRN [15] 0.292 ± 0.002 0.419 ± 0.002 2.743 ± 0.011 0.0819 ± 0.000 0.242 ± 0.001 1.972 ± 0.019 iTransformer [28] 0.171 ± 0.002 0.269 ± 0.001 1.578 ± 0.005 0.0804 ± 0.000 0.243 ± 0.002 1.774 ± 0.023 CrossFormer [42] 0.169 ± 0.001 0.270 ± 0.001 1.541 ± 0.006 0.0882 ± 0.000 0.237 ± 0.002 2.931 ± 0.059 Card [34] 0.166 ± 0.001 0.267 ± 0.001 1.576 ± 0.002 0.0778 ± 0.000 0.242 ± 0.002 1.748 ± 0.009 ESG [39] 0.168 ± 0.002 0.275 ± 0.003 1.549 ± 0.003 0.0932 ± 0.005 0.242 ± 0.001 2.322 ± 0.026 FourierGNN [40] 0.237 ± 0.003 0.369 ± 0.002 1.978 ± 0.007 0.0935 ± 0.004 0.257 ± 0.002 2.869 ± 0.038 TPGNN [26] 0.234 ± 0.002 0.367 ± 0.002 2.042 ± 0.007 0.0782 ± 0.000 0.246 ± 0.002 1.731 ± 0.030 Sumba 0.163 ± 0.001 0.260 ± 0.001 1.530 ± 0.005 0.0704 ± 0.000 0.234 ± 0.002 1.663 ± 0.013Method ETTh2 Traffic (F = 6) MAE RMSE MAPE(%) MAE RMSE MAPE(%) MICN [33]0.281 ± 0.001 0.405 ± 0.001 1.692 ± 0.006 0.275 ± 0.002 0.571 ± 0.001 3.369 ± 0.019 ModernTCN [29] 0.272 ± 0.002 0.380 ± 0.003 1.819 ± 0.015 0.261 ± 0.002 0.597 ± 0.001 3.219 ± 0.028 PatchTST [30] 0.210 ± 0.002 0.328 ± 0.002 0.956 ± 0.005 0.231 ± 0.001 0.554 ± 0.002 2.465 ± 0.012 FEDformer [45] 0.349 ± 0.006 0.508 ± 0.007 1.523 ± 0.021 0.344 ± 0.003 0.718 ± 0.007 3.845 ± 0.033 Autoformer [35] 0.322 ± 0.002 0.473 ± 0.003 1.350 ± 0.010 0.360 ± 0.002 0.730 ± 0.003 3.891 ± 0.018 Reformer [19] 0.488 ± 0.012 0.690 ± 0.007 2.604 ± 0.031 0.338 ± 0.002 0.775 ± 0.006 3.613 ± 0.025 S4 [12] 0.562 ± 0.003 0.747 ± 0.009 2.755 ± 0.021 0.301 ± 0.002 0.738 ± 0.003 2.952 ± 0.020 MTGNN [38] 0.212 ± 0.001 0.323 ± 0.002 0.994 ± 0.006 0.227 ± 0.002 0.551 ± 0.004 2.562 ± 0.007 MegaCRN [15] 0.287 ± 0.004 0.411 ± 0.003 1.437 ± 0.009 0.338 ± 0.002 0.754 ± 0.002 3.994 ± 0.015 iTransformer [28] 0.212 ± 0.001 0.329 ± 0.001 1.029 ± 0.003 0.211 ± 0.002 0.524 ± 0.001 2.269 ± 0.005 CrossFormer [42] 0.249 ± 0.002 0.387 ± 0.001 1.097 ± 0.006 0.216 ± 0.002 0.634 ± 0.003 2.140 ± 0.014 Card [34] 0.209 ± 0.001 0.319 ± 0.001 0.951 ± 0.004 0.215 ± 0.001 0.545 ± 0.002 2.341 ± 0.006 ESG [39] 0.254 ± 0.006 0.371 ± 0.013 1.027 ± 0.012 0.221 ± 0.002 0.626 ± 0.002 2.248 ± 0.003 FourierGNN [40] 0.265 ± 0.002 0.393 ± 0.006 1.933 ± 0.008 0.291 ± 0.001 0.665 ± 0.002 3.608 ± 0.028 TPGNN [26] 0.236 ± 0.001 0.379 ± 0.003 0.155 ± 0.001 0.228 ± 0.002 1.180 ± 0.004 0.077 ± 0.000 0.136 ± 0.001 0.782 ± 0.003 ModernTCN [29] 0.143 ± 0.001 0.211 ± 0.000 1.102 ± 0.003 0.065 ± 0.000 0.125 ± 0.001 0.608 ± 0.002 PatchTST [30] 0.179 ± 0.001 0.242 ± 0.001 1.137 ± 0.003 0.082 ± 0.000 0.142 ± 0.001 0.647 ± 0.002 FEDformer [45] 0.246 ± 0.002 0.349 ± 0.003 1.765 ± 0.005 0.169 ± 0.001 0.247 ± 0.001 1.224 ± 0.006 Autoformer [35] 0.462 ± 0.003 0.670 ± 0.002 2.446 ± 0.008 0.213 ± 0.001 0.290 ± 0.002 1.340 ± 0.004 Reformer [19] 0.190 ± 0.001 0.300 ± 0.002 1.506 ± 0.003 0.068 ± 0.000 0.147 ± 0.001 0.818 ± 0.004 S4 [12] 0.232 ± 0.002 0.331 ± 0.003 2.029 ± 0.007 0.092 ± 0.001 0.139 ± 0.001 0.974 ± 0.003 MTGNN [38] 0.145 ± 0.000 0.215 ± 0.000 1.127 ± 0.004 0.060 ± 0.000 0.135 ± 0.001 0.629 ± 0.002 MegaCRN [15] 0.154 ± 0.001 0.232 ± 0.003 1.209 ± 0.004 0.057 ± 0.000 0.129 ± 0.001 0.658 ± 0.002 iTransformer [28] 0.146 ± 0.001 0.215 ± 0.001 1.147 ± 0.007 0.060 ± 0.000 0.126 ± 0.001 0.609 ± 0.002 CrossFormer [42] 0.144 ± 0.001 0.219 ± 0.001 1.139 ± 0.002 0.047 ± 0.000 0.116 ± 0.001 0.553 ± 0.003 Card [34] 0.155 ± 0.001 0.226 ± 0.001 1.139 ± 0.002 0.057 ± 0.000 0.134 ± 0.001 0.684 ± 0.002 ESG [39] 0.142 ± 0.001 0.212 ± 0.001 1.089 ± 0.003 0.050 ± 0.000 0.122 ± 0.001 0.606 ± 0.002 FourierGNN [40] 0.154 ± 0.001 0.236 ± 0.002 1.292 ± 0.003 0.063 ± 0.000 0.129 ± 0.002 0.647 ± 0.003 TPGNN [26] 0 0.202 ± 0.001 0.316 ± 0.002 1.655 ± 0.003 0.083 ± 0.000 0.175 ± 0.001 0.919 ± 0.004 S4 [12] 0.244 ± 0.003 0.344 ± 0.005 2.177 ± 0.006 0.101 ± 0.001 0.158 ± 0.002 1.156 ± 0.004 MTGNN [38] 0.158 ± 0.001 0.234 ± 0.001 1.229 ± 0.003 0.083 ± 0.000 0.169 ± 0.001 0.810 ± 0.002 MegaCRN [15] 0.167 ± 0.002 0.254 ± 0.003 1.356 ± 0.005 0.087 ± 0.000 0.172 ± 0.002 0.924 ± 0.005 iTransformer [28] 0.157 ± 0.001 0.233 ± 0.002 1.241 ± 0.003 0.071 ± 0.000 0.163 ± 0.001 0.823 ± 0.002 CrossFormer [42] 0.151 ± 0.001 0.233 ± 0.002 1.234 ± 0.002 0.063 ± 0.001 0.156 ± 0.002 0.739 ± 0.003 Card [34] 0.152 ± 0.002 0.228 ± 0.003 1.182 ± 0.004 0.062 ± 0.000 0.160 ± 0.002 0.775 ± 0.003 ESG [39] 0.151 ± 0.001 0.226 ± 0.001 1.159 ± 0.004 0.068 ± 0.000 0.158 ± 0.001 0.757 ± 0.005 FourierGNN [40] 0.182 ± 0.003 0.267 ± 0.002 1.327 ± 0.003 0.096 ± 0.001 0.177 ± 0.002 0.872 ± 0.002 TPGNN [26] 0.169 ± 0.001 0.256 ± 0.002 1.292 ± 0.003 0.090 ± 0.000 0.196 ± 0.002 0.872 ± 0.002 Sumba 0.147 ± 0.001 0.222 ± 0.002 1.155 ± 0.002 0.062 ± 0.000 0.147 ± 0.001 0.734 ± 0.002 methods, and the history window size H is set to 168. All methods are trained on Nvidia V100 GPUs.Our method is implemented with PyTorch 2.0 and we use the source codes publicly released by the authors for all baseline methods. We adjust the hyperparameters of baseline methods to obtain the best performance on each dataset, and evaluate the performance of different methods in terms of Mean Absolute Error (MAE), Root Mean Square Error (RMSE), and Mean Absolute Percentage Error (MAPE). We conduct the ablation study on the Electricity, Weather, and PEMS datasets with prediction horizons F of 3 and 6 to verify the efficacy of the proposed modules. In particular, we consider the following variants of our proposed model: 3 , the performance of w/o. dynamic declines rapidly by up to 13.83%, 8.59%, and 16.31% on three metrics, which verifies the efficacy of the proposed dynamic GCN model. The absence of structure regularization (w/. U m ,V m ) and the elimination of orthogonal parameterization (w/o. orthogonality) on structured matrix basis lead to averaged performance drops of 4.95% and 3.05%, respectively. This proves the effectiveness of the proposed structure regularization and orthogonal constraint strategy.As discussed in Section 3.1, one appealing feature of our proposed Sumba is that it enables us to gain insights into the underlying time series dynamics by tracking the change of the matrix basis coefficient α over time. To verify this, we present the heatmap of the change of α on two datasets-Solar-Energy and Electricity in Figure 2 , in which the x-axis denotes the time (hourly), the y-axis represents the no. of 5 basis matrices, and the color indicates the weight of α. It can be observed that most of the weights are concentrated on two basis matrices, no. 1 and no. 4, which implies that there are two dominant spatial structures on the two datasets. Furthermore, the two dominant spatial structures appear alternatively and regularly, which actually corresponds to the day and night for the Solar-Energy dataset. This aligns well with our intuition that solar energy observations should manifest different spatial correlations as the light intensity varies. For the Electricity dataset, there is a spatial correlation (no. 4) that spans two hours and only emerges in midnight, this reveals an interesting electricity consumption pattern, which is unnoticed by previous methods. Therefore, our proposed method is able to offer more interpretable dynamics of the underlying systems through tracking α. We evaluate the impact of M on performance of our model, where the M ranges from 1, 2, 5, to 10. M = 1 corresponds to the static graph. As shown in Figure 3 , the performance of our model improves significantly when M > 1. When M = 2, our method achieves the best performance on Electricity and Solar-Energy datasets, this coincides with the fact that the two datasets have two dominant spatial structure patterns, i.e., day and night. The experiments empirically show that the best setting of M is 5 on ETTh2, Traffic, and PEMS datasets and 10 on the Weather dataset. We hypothesize this is due to the Weather dataset having more complex dynamics, and thus it requires more basis matrices to cover the various patterns. The sensitivity analysis of L, K and H is presented in Appendix D. 
conclusion: 
 In this paper, we propose a time series forecasting method with the structured matrix basis, Sumba, to capture dynamic spatial structures. To this end, we propose a novel structured parameterization and impose structure regularization on the basis to enhance parameter efficiency, the output space of the spatial structure function is thus well constrained and the generated spatial structures have lower variance. Our proposed method offers us a manner to gain insights into the dynamics of the underlying systems, and thus it is more interpretable. The experiments on six benchmark datasets verify the superiority of our proposed method. Our extensive ablation studies prove the effectiveness of each proposed component. In addition, the case study shows that our method can offer desirable interpretability. In the future, we would like to explore how to better regularize the learned matrix vector space and set the dimension of basis M in a data-driven manner. We will also explore the possibility of integrating our proposed spatial structure modeling with other temporal encoders such as Transformers, and Structured State Space models, and apply our interpretable dynamics into more datasets to discover more interesting and hidden patterns. 
safeguards: 
 Question: Does the paper describe safeguards that have been put in place for responsible release of data or models that have a high risk for misuse (e.g., pretrained language models, image generators, or scraped datasets)?Answer: [NA]Justification: The paper poses no such risks.Guidelines:• The answer NA means that the paper poses no such risks.• Released models that have a high risk for misuse or dual-use should be released with necessary safeguards to allow for controlled use of the model, for example by requiring that users adhere to usage guidelines or restrictions to access the model or implementing safety filters. • Datasets that have been scraped from the Internet could pose safety risks. The authors should describe how they avoided releasing unsafe images. • We recognize that providing effective safeguards is challenging, and many papers do not require this, but we encourage authors to take this into account and make a best faith effort.12. Licenses for existing assets Question: Are the creators or original owners of assets (e.g., code, data, models), used in the paper, properly credited, and are the license and terms of use explicitly mentioned and properly respected?Answer: [Yes] Justification: We cite the original paper that produced the code package or dataset. See Section 2 and Section 4.Guidelines:• The answer NA means that the paper does not use existing assets.• The authors should cite the original paper that produced the code package or dataset.• The authors should state which version of the asset is used and, if possible, include a URL. • The name of the license (e.g., CC-BY 4.0) should be included for each asset.• For scraped data from a particular source (e.g., website), the copyright and terms of service of that source should be provided. • If assets are released, the license, copyright information, and terms of use in the package should be provided. For popular datasets, paperswithcode.com/datasets has curated licenses for some datasets. Their licensing guide can help determine the license of a dataset. • For existing datasets that are re-packaged, both the original license and the license of the derived asset (if it has changed) should be provided. • If this information is not available online, the authors are encouraged to reach out to the asset's creators. 
new assets: 
 Question: Are new assets introduced in the paper well documented and is the documentation provided alongside assets? Answer: [NA] Justification: The paper does not release new assets. Guidelines:• The answer NA means that the paper does not release new assets.• Researchers should communicate the details of the dataset/code/model as part of their submissions via structured templates. This includes details about training, license, limitations, etc. • The paper should discuss whether and how consent was obtained from people whose asset is used. • We recognize that the procedures for this may vary significantly between institutions and locations, and we expect authors to adhere to the NeurIPS Code of Ethics and the guidelines for their institution. • For initial submissions, do not include any information that would break anonymity (if applicable), such as the institution conducting the review. 
 
 Appendix: Recall that given a matrix A = UΣV ⊤ , the multiplication Ax can be interpreted as 1) transforming x into the new coordinate system x ′ ≜ V ⊤ x, 2) doing the elementwise product in the new coordinate y ′ = Σx ′ , and 3) transforming the coordinate system again y = Uy ′ . Proof.EQUATIONEQUATIONEQUATIONEQUATION)The proof is completed.In our experiments, we compare our proposed Sumba with state-of-the-art (SOTA) baseline models.• • MTGNN: It builds an adaptive directed graph using learnable node embedding and aggregates information along spatial dimensions through mix-hop propagation.• MegaCRN: It exploits memory-enhanced node embedding to build the graph structure.• iTransformer: It embeds the whole time series into a spatial token and captures spatial correlations using the self-attention mechanism.• Crossformer: It constructs both temporal and spatial dependencies using the attention mechanism and introduces a router mechanism to decrease the complexity of the spatial attention module.• Card: It introduces summarized spatial tokens to decrease the complexity of the spatial attention module and employs a token blend mechanism in the temporal attention module to extract local temporal correlations.• ESG: It designs multi-scale evolving node embedding based on Gated Recurrent Unit to capture dynamic spatial dependencies.• FourierGNN: It builds dynamic graphs by merging spatial and temporal dimensions and designs the Fourier Graph Operator to decrease complexity.• TPGNN: It regards the dynamic spatial correlations as the matrix polynomial graph where the time-varying coefficients are determined by timestamps.The evaluation metrics including MAE, RMSE, and MAPE are as follows:EQUATIONwhere Ω denotes the index set along temporal and spatial dimensions.Table 4 provides the results with prediction lengths of 12 and 24. The results underscore the significant role of explicitly capturing spatial correlations in improving forecasting accuracy. Furthermore, dynamic graphs exhibit a greater ability to enhance the forecasting performance than modeling static spatial correlations. Among the baseline methods, Card, Crossformer, and ESG demonstrate superior forecasting accuracy. Notably, our model achieves state-of-the-art performance, yielding improvements more than 20% in terms of MAE and 25% in terms of MAPE over the existing methods on the Solar-Energy dataset. In addition, we conduct experiments with a forecasting horizon of 96, and as Table 5 shows, Sumba gives rise to favorable performance in comparison to PatchTST, FEDformer, and iTransformer which are designed specifically for long-term forecasting. Table 6 : The complexity of different methods.Method ComplexityMTGNN [38] O(N 2 D) MegaCRN [15]O(N 2 D) iTransformer [28] O(N 2 D) CrossFormer [42] O(N rD) Card [34] O(N rD) ESG [39] O(N 2 D) FourierGNN [40] O(N Hlog(N H)) TPGNN [26] O(N 2 D)w/o. Common U-V basis O(M N KD) Sumba O(N KD)Convergence Speed Figure 7 illustrates the training curves for the first 10 epochs on the Electricity and Weather datasets. Sumba, ModernTCN, and ESG exhibit fast convergence, achieving it within 5 epochs. In contrast, methods such as Card, FourierGNN, and TPGNN take several tens of epochs to reach convergence. Although our proposed Sumba demonstrates superior performance in multivariate time series forecasting, the dimension of matrix basis M is empirically tuned in our experiments. Besides, the long-term forecasting capability of the proposed model requires further enhancement.In this paper, we propose a novel time series forecasting method to capture dynamic spatial correlations with a structured matrix basis. Our research aims to contribute to the advancement of the relevant community without any negative social impact.NeurIPS Paper Checklist Justification:Guidelines:• The answer NA means that the abstract and introduction do not include the claims made in the paper. • The abstract and/or introduction should clearly state the claims made, including the contributions made in the paper and important assumptions and limitations. A No or NA answer to this question will not be perceived well by the reviewers. • The claims made should match theoretical and experimental results, and reflect how much the results can be expected to generalize to other settings. • It is fine to include aspirational goals as motivation as long as it is clear that these goals are not attained by the paper.38th Conference on Neural Information Processing Systems (NeurIPS 2024). 


Decision: Thanks for your submission to NeurIPS! The reviewers have reach a
consensus, and it is clear that the reviewers are positive with its
quality. The recommended decision from the AC is accept. The authors may
take into account the feedbacks of the reviewers to further improve the
paper. There could be open or extension questions that may not be easy
to answer while possible comments (probably in the supplemental) could
help.

Bests

The AC

Reviews: review: paper_summary: paper_summary: This paper presents a method to learn dynamic spatial structures in spatio-temporal forecasting tasks. Specifically, it proposes to parametrize the dynamic structures with a convex combination of fixed matrix bases, and the bases are further confined to be in the same coordinate system. Beyond that, it also imposes low-rank assumption on the coordinates to further reduce complexity. Empirical evidences show the proposed method achieve impressive accuracy in a number of spatio-temporal forecasting benchmarks, and the found spatial structures are highly interpretable.
strength_weakness: strength_weakness: The forecasting horizon in the experiments is quite different from the setting in some baselines, such as PatchTST and FEDformer. While long-term forecasting is not major claim of the paper, I wonder if the method scales well with prediction length.
questions: questions: 1. Why is TCN selected for temporal encoding?
2. Is there any issue that the learned bases degenerate to be trivial or converge to be close?
limitations: limitations: Limitations are discussed, and I agree that the selection of $M$, i.e. the number of bases, is highly empirical and ad-hoc.
paper_summary: paper_summary: This paper presents a multivariate forecasting model underlined by a dynamic spatial structure generation function, enabled by SVD-based parameterization and theoretically-bounded output space. Experiments on six commonly benchmark datasets in comparison to several existing baselines demonstrated the overall improved forecasting performance of the presented method as measured by MAE/RMSE/MAPE. Additional ablation studies further demonstrated the benefits of dynamic coefficient generation and the structured parameterization.
strength_weakness: strength_weakness: The obtained performance gain (Table 1) was overall marginal (up to 2-3 decimal points). The practical implication of such margin of improvements is not clear and need to be clarified. The effect of hyperparameters and random initialization on such margin of improvements also needs to be examined — adding statistics to the results over different random seeds of experiments is important.

In switching dynamic systems, it is also common to model the transition matrix over time as a linear combination of several global matrices, with time-varying mixing coefficients (without considering the proposed SVD decomposition and parameterization). It’d strengthen the paper to add discussion about the relation with this line of works, as well as experimental comparison.


Based on the results in Fig 3, it appears that increasing M from 1 to 2 in general introduced small difference in performance except in the Electricity dataset. This again raises some question on the significance of the dynamic spatial structure introduced in this paper. Stronger clarification on performance improvement and adding error bars to the results in Fig 3 will be important for verifying the contribution of the paper.
questions: questions: Clarification on the observed margin of performance improvements and statistics/error bars to the results will be appreciated. 

Better clarification and empirical results on relation with switching dynamics systems will be appreciated. 

What is the computational cost difference in the ablation with or without common U-V basis (i.e., vs. w/ Um,Vm)?
limitations: limitations: The authors did not include sufficient discussion about the limitation or potential negative societal impact of the presented work.
paper_summary: paper_summary: This paper proposes a novel approach for effectively capturing spatial and temporal correlations in multivariate time series forecasting and enhancing the interpretability of prediction models. The aim is to address the shortcomings of existing methods, which often fail to adequately reflect dynamic spatial correlations and exhibit high variability. To tackle this issue, the paper introduces a method that bypasses the traditional two-stage learning process and directly generates dynamic structures using a structured matrix basis. Furthermore, the basis matrices are parameterized through singular value decomposition (SVD), and all basis matrices share the same orthogonal matrices to improve the efficiency of model training.
strength_weakness: strength_weakness: 1. The proposed model may be complex to implement due to the use of structured matrix bases and singular value decomposition, potentially leading to a steep initial learning curve in practical applications. Does the appendix indicate similar results for datasets other than Electricity?
2. While the model has been validated on various datasets, it remains to be seen how well it performs on very specialized domain data or data with unclear characteristics (i.e., weak temporal dependencies or spatial correlations). How does the model address these scenarios?
questions: questions: 1. The proposed model may be complex to implement due to the use of structured matrix bases and singular value decomposition, potentially leading to a steep initial learning curve in practical applications. Does the appendix indicate similar results for datasets other than Electricity?
2. While the model has been validated on various datasets, it remains to be seen how well it performs on very specialized domain data or data with unclear characteristics (i.e., weak temporal dependencies or spatial correlations). How does the model address these scenarios?
limitations: limitations: Including more analysis and comparison with graph-based spatio-temporal models would enhance the paper.

Comments: comment: Thank you for your timely response and the suggestions for improving this paper. We will further revise the paper according to your suggestions and the results in the rebuttal, including adding more learning curves, the corresponding experiments to verify the performance on datasets with independent channels, and more analyses and comparisons with graph-based spatio-temporal models.
Thanks for your response, after reading your response and other reviews, I am updating my score.
We sincerely thank the reviewer for raising the rating! We will include error bars in our revised manuscript.
Thanks for the rebuttal effort -- I appreciate the new baseline and the error bars added to the pdf. It'd be highly recommended for the authors to add the error bar (or their numerical values) to the main text of the paper if accepted, as that's the best evidence for the significance of the margins of improvements (instead of stating that the range is common in literature). I will raise my rating.