Attractor Memory for Long-Term Time Series Forecasting: A Chaos Perspective

Authors: jiaxi  hu
yuehong  hu
wei  chen
ming  jin
shirui  pan
qingsong  wen
yuxuan  liang
Venue: NeurIPS_cc_2024_Conference
Published: 2024-05-09 17:00:23

Abstract: In long-term time series forecasting (LTSF) tasks, an increasing number of works have acknowledged that discrete time series originate from continuous dynamic systems and have attempted to model their underlying dynamics. Recognizing the chaotic nature of real-world data, our model, Attraos, incorporates chaos theory into LTSF, perceiving real-world time series as low-dimensional observations from unknown high-dimensional chaotic dynamical systems. Under the concept of attractor invariance, Attraos utilizes non-parametric Phase Space Reconstruction embedding along with a novel multi-resolution dynamic memory unit to memorize historical dynamical structures, and evolves by a frequency-enhanced local evolution strategy. Detailed theoretical analysis and abundant empirical evidence consistently show that Attraos outperforms various LTSF methods on mainstream LTSF datasets and chaotic datasets with only one-twelfth of the parameters compared to PatchTST.

Main paper: 
 introduction: 
 In the intricate dance of time, time series unfold. Emerged from continuous dynamical systems [36, 11, 39] in the physical world, these series are meticulously collected at specific sampling frequencies. Like musical notes in a composition, they harmonize, revealing patterns that resonate through the symphony of temporal evolution. In this realm, Long-term Time Series Forecasting (LTSF) stands as one of the enduring focal points within the machine learning community, achieving widespread recognition in real-world applications, such as weather forecasting, financial risk assessment, and traffic prediction [28, 22, 30, 27, 18] .Building on the success of various deep LTSF models [49, 54, 52, 50, 29, 22] , which primarily leverage neural networks to learn temporal dependencies from discretely sampled data. Currently, researchers [46, 32] have been investigating the application of Koopman theory [48, 25] in recovering continuous dynamical systems, which applies linear evolution operators to analyze dynamical system characteristics in a sufficiently high-dimensional Koopman function space. Nevertheless, the existence of Koopman space relies on the deterministic system, posing challenges given the chaotic nature of real-world time series data, evidenced through the Maximal Lyapunov Exponent in Appendix E.2.In this paper, inspired by chaos theory [8] , we revisit LSTF tasks from a chaos perspective: Linear or complex nonlinear dynamical systems exhibit stable patterns in their trajectories after sufficient evolution, known as attractors. As illustrated in Figure 1 (a), attractors can be classified into four types: Fixed Point, indicating stable, invariant systems; Limited cycle, representing periodic behavior; Limited Toroidal, exhibiting quasi-periodic behavior with non-intersecting rings in a 2D plane, reflecting temporal distribution shifts; and Strange Attractor, characterized by nonlinear behavior and complex, non-topological shapes. Supported by chaos theory, we can transcend the limitations of deterministic dynamical systems to construct generalized dynamical system models. Figure 1(b-c ) showcases various classical chaotic dynamical systems and the dynamical trajectories of real-world LTSF datasets using the phase space reconstruction method [9] . Notably, the dynamical system trajectories in these LTSF datasets exhibit fixed structures akin to those in typical chaotic systems.Given this chaos perspective, we consider real-world time series as stemming from an unidentified high-dimensional underlying chaotic system, broadly encompassing nonlinear behaviors beyond periodicity. Our focus centers on recovering continuous chaotic dynamical systems from discretely sampled data for LTSF tasks, with the goal of predicting future time steps through the lens of attractor invariance. Specifically, this problem can be decomposed into two key questions: (i) how to model the underlying continuous chaotic dynamical system based on discretely sampled time series data; (ii) how to enhance forecasting performance by utilizing the attractors within the system.In this context, Attraos emerges with the goal of capturing the underlying order within the seeming chaos via attractors. For tackling the first question, we employ a non-parametric phase space reconstruction method to recover the temporal dynamics and propose a Multi-resolution Dynamic Memory Unit (MDMU) to memorize the structural dynamics within historical sampled data. Specifically, as polynomials have been proven to be universal approximators for dynamical systems [4] , MDMU expands upon the work of the State Space Model (SSM) [12, 10, 18] to different orthogonal polynomial subspaces. This allows for memorizing diverse dynamical structures that encompass various attractor patterns, while theoretically minimizing the boundary of attractor evolution error.To address the second question, we devise a frequency-enhanced local evolution strategy, which is built upon the recognition that attractor differences are amplified in the frequency domain, as observed in the field of neuroscience [5, 14, 6] . Concretely, for dynamical system components that belong to the same attractor, we apply a consistent evolution operator to derive their future states in the frequency domain. Our contributions can be summarized as follows:• A Chaos Lens on LTSF. We incorporate chaos theory into LTSF tasks by leveraging the concept of attractor invariance, leading to a principal way to model the underlying continuous dynamics• Efficient Dynamic Modeling. Our model Attraos employs a non-parametric embedding to obtain high-dimensional dynamical representations, leverages MDMU to capture the multi-scale dynamical structure, and performs the evolution in the frequency domain. Remarkably, Attraos achieves this with only about one-twelfth the parameter count of PatchTST. Furthermore, we utilize the Blelloch scan algorithm [3] to enable efficient computation of the MDMU.• Empirical Evidence. Various experiments validate the superior performance of Attraos. Besides Leveraging the properties of chaotic dynamical systems, we explore their extended applications in LTSF tasks, focusing on chaotic evolution, modulation, representation, and reconstruction. 
preliminary: 
 Attractor in Chaos Theory. In chaos theory, the interaction of three or more variables exhibiting periodic behavior gives rise to a complex dynamical system characterized by chaos. According to Takens's theorem [43, 34] , assuming an ideal dynamical system F : M → M that "lives" on attractor A in manifold space M which locally C N (N-times differentiable), time series data {z i } ∈ R can be interpreted as the observation of it by an unknown observation function h. To explore the properties of the unknown ideal dynamical system, we can employ the phase space reconstruction (PSR) method to establish an approximation K : R m → R m which lives in differential homomorphism attractor Ã in the Euclidean space with suitable dimension m [9] . The whole process is illustrated in Equation (1) , where {z i }, {u i } are the sampled data from two dynamical systems. Strictly speaking, in our paper, the chaotic attractor structure Ã = { Ãi } we focused on is in phase space R m . To facilitate understanding, we further provide a visual example of the Lorenz96 system in Appendix E. 3. In the forecasting stage, the local prediction method emerges as a prominent one: u i+1 = K (i) (u i ), where the local evolution K (i) can be either linear or nonlinear neural network [45, 2, 40] , with the parameter being shared among the points in the neighborhood of u i or belong to the same local attractor. Considering the universal approximation capabilities of polynomials for dynamical systems, we leverage the polynomial to describe the chaotic dynamical structures.ai ∈ A ⊂ M F -→ ai+1 ∈ A ⊂ M ↓ h ↓ h zi ∈ R zi+1 ∈ R ↓P SR ↓P SR ui ∈ Ã ⊂ R m K -→ ui+1 ∈ Ã ⊂ R m (1)x ′ (t) = Ax(t) + Bu(t) (2a)EQUATIONEQUATIONPolynomial Projection with Measure Window. We only consider the first part of the SSM (2a), which is a parameterized map that transforms the input u(t) into an N -dimensional latent space. According to Hippo [11] , it is mathematically equivalent to: given an input u(s), a set of orthogonal polynomial basisϕ n (t, s) that t -∞ ϕ m (t, s)ϕ n (t, s)ds = δ m,n, and an inner product probability measure µ(t, s). This enables us to project the input u(s) onto the polynomial basis along time dimension (3), and we can combine ϕ n (t, s)ω(t, s) as a kernel K n (t, s) (4). When ω(t, s) is defined in a time window I[t, t + θ], it represents approximating the input over each window θ.⟨u, ϕ n ⟩ µ = t -∞ u(s)ϕ n (t, s)ω(t, s)ds (3)x n (t) = u(s)K n (t, s)I(t, s)dsWhen the basis and measure are solely dependent on time t, it can be expressed in a convolution form (2b). In this paper, we will utilize this property to project the dynamical trajectories {u n } in phase space onto the polynomial spectral domain with kernel e tA B (2c) for characterization.Question: Does the paper discuss the limitations of the work performed by the authors?Answer: [Yes] .Justification: In Section 5, we discuss the limitations of this work.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. 
theoretical analysis & methods: 
 The overall structure of Attraos is illustrated in Figure 2 . In this section, we provide a comprehensive description of its components, including the Phase Space Reconstruction embedding, the Multi-Resolution Dynamic Memory Unit (MDMU), and the frequency-enhanced local evolution, as well as the efficient computational methods employed.According to chaos theory, the initial step involves constructing a topologically equivalent dynamical structure through the PSR. The preferred embedding method is typically the Coordinate Delay Reconstruction [34] , which does not rely on any prior knowledge of the underlying dynamical system. By utilizing the discretely sampled data {z i } and incorporating two hyper-parameters, namely, embedding dimension m and time delay τ , a high-dimensional dynamical trajectory {u i } in phase space can be constructed by Eq. (5) .EQUATIONFigure 2 : Overall architecture of Attraos. Initially, the PSR technique is employed to restore the underlying dynamical structures from historical data {zi}. Subsequently, the dynamical system trajectory is fed into MDMU, projected onto polynomial space G N θ using a time window θ and polynomial order N . Gradually, a hierarchical projection is performed to obtain more macroscopic memories of the dynamical system structure. Finally, local evolution operator K (i) in the frequency domain is employed to obtain future state, thereby for the prediction.For multivariate time series data with C variables, we have observed considerable variations in the Lyapunov exponents of each variable, hence a channel-independent strategy [33] is employed to construct a unified dynamical system K ⊂ R m . To accelerate model convergence and reduce complexity, we apply non-overlapping patching to obtain K patch (6) . We denote the number of patches as L, using u ∈ R B×L×D to represent the tensor used for computing, where D = mp. The determination of m and τ is achieved by applying the CC method [23] as shown in Appendix C.1.Remark 3.1. This represents the pioneering non-parametric embedding in LTSF, effectively reducing the model parameters in the embedding and output projection process (m is typically single-digit).Remark 3.2. In a large body of dynamical literature [44, 41, 20, 25] , local linear approximation serves as an effective method for modeling dynamical systems, providing a basis for the effectiveness of the patching operations in this paper.Proposition 1. A = diag{-1, -1, . . . } is a rough approximation of normal Hippo-LegT [11] matrix, which utilizes polynomial projection under a finite measure window (Lebesgue measure). Remark 3.3. All proofs in this section can be found in Appendix B.Adhere to Mamba [10] , we generate B and measure window θ by linear layer, and propose a novel parameterized method (Proposition 1) for A to instantiate Eq. (2a):EQUATIONwhere A ∈ R D×N represents the dynamical characteristics of the system's forward evolution in the polynomial space. Due to its diagonal nature, its representational capacity is comparable to R D×N ×N ; Matrix B ∈ R B×L×N controls the process of projecting u onto the polynomial domain like a gate mechanism; The learnable approximation window ∆ ∈ R B×L×D , similar to an attention mechanism, enables adaptive focus on specific dynamical structures (attractors).Remark 3.4. In the Hippo theory, the measure window is denoted by θ, while in SSMs [10, 12] , the discrete step size is represented by ∆. These two terms can be considered approximately equivalent.As shown in Figure 3 (a), in practical computations, we need to discretize Eq. (2a) to fit the discrete dynamical trajectories. We apply zero-order hold (ZOH) discretization [11] to matrix A, while opting for a combination of Forward Euler discretization for B (instead of the commonly used EQUATIONUp to this point, the construction of the underlying continuous dynamical system K and its application to discretely sampled data z n have been established. However, in this scenario, we are still limited to a single representation of the dynamical structures with measure window θ. The strange attractors, on the other hand, are often composed of multiple fundamental topological structures. Therefore, we require a multi-scale hierarchical representation of dynamical structures to capture their complexity.To address this, as illustrated in Figure 3 (b), we progressively increase the length of the window θ by powers of 2. The region previously approximated by g θ1 ∈ G N θ (left half) and g θ2 ∈ G N θ (right half) will now be approximated by g 2θ ∈ G N 2θ . Since the piecewise polynomial function space can be defined as the following form:EQUATIONwith polynomial order k ∈ N, piecewise scale r ∈ Z + ∪ {0}, and piecewise internal indexl ∈ {0, 1, ..., 2 r -1}, it is evident that dim(G k (r) ) = 2 r k, implying that G k θ possesses a superior function capacity compared to G k 2θ .All functions in G k 2θ are encompassed within the domain of G k θ . Moreover, since G k θ and G k 2θ can be represented as space spanned by basis functions {ϕ θ i (x)} and {ϕ 2θ i (x)}, any function including the basis function within the G k 2θ space can be precisely expressed as a linear combination of basis functions from the G k θ space with a proper tilted measure µ 2θ :EQUATIONThe projection coefficients on these two spaces can be mutually transformed using the linear projection matrix H and its inverse matrix H † based on the odd or even positions along the L dimension in x.EQUATIONRemark 3.5. Although ∆ is akin to an attention mechanism leads to different measure windows, ie., θ 1 ̸ = θ 2 , it still maintains the linear projection property for up and down projection:G θ1 , G θ2 ↔ G θ1+θ2 . In our illustration, we have used a unified measure window for simplicity. Theorem 2. (Approximation Error Bound) Function f : [0, 1] ∈ R is k times continuously differentiable, the piecewise polynomial g ∈ G N r approximates f with mean error bounded as:∥f -g∥ ≤ 2 -rN 2 4 N N ! sup x∈[0,1] f (N ) (x) .Iteratively repeating this process enables us to model the dynamical structure from a more macroscopic perspective. Theorem 2 indicates that under this approach, the polynomial projection error has convergence of order N . Hyper-parameter analysis can be found in Figure 5 . Remark 3.6. When the weight function is uniformly equal across the dynamical structure, H θ1 and H θ2 are shared in each projection level as the projection matrix for the left and right interval.Theorem 3. The mean attractor evolution errorK • Ã -Ã of evolution operator K = {K (i) } is bounded by ∥K • Ã -Ã∥(N -1)E (-β∇ i + 1), with the number of random patterns N ≥ √ pc d-14 stored in the system by interaction paradigm E in an ideal spherical retrieval paradigm.Based on this hierarchical projection, we additionally want to uphold the constancy of attractor patterns throughout the evolution, which is equivalent to minimizing attractor evolution errors. According to Theorem 3, it is imperative to ensure the separation between attractors, denoted as ∇ i := min j,j̸ =i ÃT i Ãi -ÃT i Ãj , is sufficiently large. While there is an intersection between the G N θ and G N 2θ spaces, which limits the attainment of a sufficiently large ∇. To address this issue, we define an orthogonal complement space asG N θ = S N θ G N 2θ , to establish a series of orthogonal function spaces {S N θ , S N 2θ , ..., S N 2 L θ , G N 2 L θ }.We can extend Eq. ( 12) as:EQUATIONEQUATIONTheorem 4 states that the coarsest-grained G space, along with a series of orthogonal complement S spaces, can approximate any given dynamical system structure with finite error bounds in Theorem 2. Theorem 4. (Completeness in L 2 space) The orthonormal system B N = {ϕ j : j = 1, . . . , N } ∪ {ψ rl j : j = 1, . . . , N ; r = 0, 1, 2, . . . ; l = 0, . . . , 2 r -1} spans L 2 [0, 1]. Remark 3.7. H, G, H † , G † ∈ R N ×N are obtained by applying Gaussian Quadrature to Legendre polynomials [13] . The gradients of these matrices are subsequently utilized for adaptive optimization.The hierarchical projection can be implemented explicitly using iterative display or can be efficiently computed through an implicit implementation proposed in the following Section 3.3.Following chaos theory, we employ a local evolution method x i+1 = K (i) (x n ) to forecast future state. Given dynamics representations s/x, the subsequent step involves partitioning the attractor and utilizing the operator K i belonging to A i for system evolution. We present three evolution strategies:• Direct Evolution: We employ each representation x t /s t ∈ R D×N as the feature to partition adjacent points in the dynamical system trajectory into the same attractor using the K-means method. Subsequently, we evolve these points using a local operator K (i) . • Frequency-enhanced Evolution: Inspired by neuroscience [6] , where attractor structures are amplified in the frequency domain, we first obtain the frequency domain representation of the dynamical structure through Fourier transformation. Considering the dominant modes as attractors, we employ K (i) to drive the system's evolution in the frequency domain.• Hopfield Evolution: Hopfield networks [16] are designed specifically for attractor memory retrieval (see Appendix A.2). In our approach, we utilize a modern version [38] of the Hopfield network for the evolution of dynamical systems, employing cross-attention operations. We treat the trainable attractor library as the Key and Value, while different scales of dynamical structure representations {s θ , s 2θ , ..., s 2 L θ , x 2 L θ } serve as the Query, enabling sequence-to-sequence evolution.Our experimental results in Table 4 demonstrate that the frequency-enhanced evolution strategy outperforms others comprehensively, and we introduce two implementation approaches:Explicit Evolution. Initially, the finest dynamical representation x θ is obtained using K, and then expanded to multiple scales s θ , s 2θ , ..., s 2 L θ , x 2 L θ . By applying the Fourier Transform, we select M low-frequency components as the primary modes. Each mode undergoes linear evolutionK (i) = W i ∈ C N ×N, followed by back projection to the original scale.Implicit Evolution. However, explicit evolution methods inevitably increase the time complexity.To address this issue, inspired by the Blelloch algorithm and hierarchical projection, which both utilize tree-like computational graphs, we propose an implicit evolution method.A Blelloch scan [3, 42] (Figure 3 (d)) defines a binary operator q i •q j := (q j,a ⊙q i,a , q j,a ⊗q i,b +q j,b ) used to compute the linear recurrence x k = Ax k-1 + Bu k . We take L = 4 as an example:Up sweep for even positionr2 = c1 • c2 = A, Bu1 • A, Bu2 = A 2 , ABu1 + Bu2 q4 = c3 • c4 = A, Bu3 • A, Bu4 = A 2 , ABu3 + Bu4 r4 = r2 • q4 = A 2 , ABu1 + Bu2 • A 2 , ABu3 + Bu4 = A 4 , A 3 Bu1 + A 2 Bu2 + ABu3 + Bu4 .r1 = r0 • c1 = (I, 0) • A, Bu1 = A, Bu1 r3 = r2 • c3 = A 2 , ABu1 + Bu2 • A, Bu3 = A 3 , A 2 Bu1 + ABu2 + Bu3 x1, x2, x3, x4 = r1[R], r2[R], r3[R], r4[R].The process commences by computing the values of variable x at even positions through an upward sweep. Subsequently, these even position values are employed during a downward sweep to calculate the values at odd positions. The corresponding value of x resides in the right node of r. Thus, we can modify the binary operator as q i • q j := (q j,a ⊙ q i,a , H i ⊗ (q j,a ⊗ q i,b + q j,b )), thereby implicitly integrating hierarchical projection into the scanning operation. This leads to the scales of x i :scale(x i ) =        0 if i = 0 scale(x i-1 ) + 1 if i is odd log 2 (i)if i is even and a power of 21if i is even and not a power of 2We simplify the original H θ1/θ2 , G θ1/θ2 , H †θ1/θ2 , G †θ1/θ2 with just H ∈ R B×L×N , which is generated directly through a linear layer, and omit the reconstruction process. This approach directly sparsifies the kernels e At B in different subspaces and using the linear layer to generate hierarchical space projection matrix. Afterwards, we learn the evolution using the data x in the frequency domain.EQUATIONIn this paper, we utilize this indirect and efficient hierarchical projection as the default setting. Ultimately, Attraos projects from the polynomial spectral space back to the phase space by employing another gating projection W out ∈ R B×L×N ×1 (15) to x ∈ R B×L×D×N , and derives the prediction results using an observation function parameterized by W h ∈ R LD×H (flattening the patches).Question: For each theoretical result, does the paper provide the full set of assumptions and a complete (and correct) proof? • Providing as much information as possible in supplemental material (appended to the paper) is recommended, but including URLs to data and code is permitted. 
experiments: 
 In this section, we commence by conducting a comprehensive performance comparison of Attraos against other state-of-the-art models in seven mainstream LTSF datasets along with two typical chaotic datasets, namely Lorenz96-3d and Air-Convection, followed by ablation experiments pertaining to model architectures. Furthermore, leveraging the properties of chaotic dynamical systems, we explore their extended applications, including experiments on chaotic evolution, representation, reconstruction, and modulation (Appendix E). Finally, we provide the complexity analysis and robustness analysis. For detailed information regarding baseline models, dataset descriptions, experimental settings, and hyper-parameter analysis, please refer to Appendix D. 1 , we can observe that: (a) Attraos consistently exhibits the best performance, closely followed by RWKV-TS and Koopa, which underscores the crucial role of modeling temporal dynamics in LTSF tasks. (b) The models based on state space models (Mamba4TS, S-Mamba) generally outperform the Transformer-based models (PatchTST, InvTrm), indicating the potential superiority of state space models as fundamental frameworks for temporal modeling. (c) The performance of the GPT-TS model, which relies on a pre-trained large language model, is relatively average, suggesting the inherent challenge in directly capturing the dynamics of temporal data using such models. A promising avenue for future research lies in training a dynamical foundational model from scratch on large-scale physical datasets or leveraging pre-training to obtain an attractor tokenizer that is better suited for inputs to the large language model. Chaotic Datasets. Table 2 presents the results on both artificial and real-world chaotic datasets. It can be observed that: (a) Attraos exhibits superior performance on both datasets, thanks to the utilization of the PSR and MDMU modules, which effectively capture the multi-scale attractor structures. (b) In Lorenz96 dataset, where there is a prior knowledge about the phase space dimension, Attraos outperforms other models by a significant margin. This highlights the importance of PSR in recovering the complete temporal dynamics. (c) Apart from Attraos, the linear model (DLinear) demonstrate the best predictive results due to their robustness. Conversely, deep learning models based on transformers exhibit weaker performance and fail to model chaotic dynamics effectively.Ablation Studies. Next, we turn off each module of Attraos to assess their individual effects. As shown in Table 3 , we observe consistent performance decline of Attraos when deleting each module: (a) The removal of Phase Space Reconstruction exhibits the most severe performance degradation, indicating that the process of reconstructing the dynamical structure through PSR forms the foundation for Attraos' efficient capture of attractor structures. (b) Multi-scale hierarchical projection effectively captures the complex topological structure of the singular attractor, leading to improved performance. However, overfitting may occur in certain prediction lengths. (c) Time-varying B and W out acts as a gating attention mechanism, allowing for a more focused emphasis on dynamical structure segments that potentially contain attractor structures, thereby enhancing performance. (d) The initialization method proposed for A matrix demonstrates marginal yet consistent improvements, underscoring the importance of prior inductive bias for machine learning models. (e) Frequency domain evolution methods significantly reduce temporal noise information and amplify attractor structures. We will further analyze the importance of frequency domain evolution in subsequent analysis.Chaotic Evolution Strategy. We further compare various dynamical system evolution strategies mentioned in Section 3.3. From Table 4 , it is evident that the frequency-enhanced evolution strategy outperforms the others. Moreover, our proposed efficient implicit evolution method can adaptively explore multi-scale dynamical structure information, avoiding redundant cyclic computations and mitigating overfitting. (b) An inherent characteristic of time series data is significant noise, making it challenging to capture the underlying dynamical structures in the time domain. Direct evolution strategies, whether linear or non-linear neural network-based, do not yield satisfactory results. Moreover, according to the theorem 2 in FiLM [53] , the recursion computation for dynamic projection further accumulates noise information. (c) Applying Hopfield networks in the time domain also proves to be unsatisfactory, and even adding more patterns (Key and Value) can have adverse effects. A potential solution is to apply Hopfield networks in the frequency domain instead.Complexity Analysis. As depicted in Figure 4 , we present a comprehensive visual analysis comparing Attraos with various baseline models in terms of their average performance on the ETTh1 dataset. The x-axis represents the training time, the y-axis represents the test loss, and the circle radius corresponds to the model parameters. In this analysis, we substituted GPT-TS with FiLM due to its limited relevance to this specific evaluation. The results clearly demonstrate that Attraos surpasses other models in both time and space complexity, maintaining a significant advantage. Notably, when compared to the PatchTST model with a hidden dimension of 256 (2.4M parameters), Attraos (0.2M parameters) possesses only one-twelfth of its parameter count.Robustness Analysis. We add a 0.1 * N (0, 1) Gaussian noise to the training dataset to test the robustness of Attraos. As shown in Table 5 , it can be observed that Attraos exhibits strong robustness against noisy data, and increasing the level of noise can even lead to further performance improvement. This is attributed to the frequency domain evolution strategy, where we retain only the dominant modes as attractor structures, effectively removing the noise information. Furthermore, an interesting phenomenon has been observed in our experiments: as noise is introduced, the model's convergence speed increases. This discovery warrants further exploration in future studies. Chaotic Reconstruction. As illustrated in the left of Figure 5 , we visualize the phase space forecasting results of Attraos on the Lorenz96 system. It can be observed that: Although some minor details may be missing due to the sparsity introduced by the frequency domain evolution, Attraos successfully reconstructs the chaotic dynamic structure of Lorenz96. Moreover, modeling time series based on the dynamics structure of the phase space can be viewed as a form of data augmentation, e.g., two-dimensional time figuring [50] or seasonal decomposition [54] .Chaotic Representation & Hyper-parameter Analysis. In the right of Figure 5 , we validated the impact of polynomial dimensions on the model performance. Noteworthy observations include:As noted in substantial literature on SSMs, a polynomial dimension of 256 is generally required to approximate input time series signals with sufficient accuracy. However, we found that the patch operation effectively reduces this threshold in a linear fashion. For instance, in the ETT dataset with a phase space dimension of 4, the required polynomial dimension drops to 256/4 = 64 dimensions. Discussion. Recently, across various fields of machine learning, an increasing number of works have focused on the underlying physical properties hidden within real-world observational data [20, 21, 26, 51] . By leveraging physical priors, these models achieve significant improvements in generalization, accuracy, and interpretability. We hope Attraos introduces a fresh perspective to the LTSF community and encourages the emergence of physics-guided time series analysis models. 
conclusion and future work: 
 LTSF tasks have long been a focal point of research in the machine-learning community. However, mainstream deep learning models currently overlook the crucial aspect that time series data is derived from discretely sampling underlying continuous dynamical systems. Inspired by chaotic theory, our model, Attraos, considers time series as generated by a generalized chaotic dynamical system. By leveraging the invariance of attractors, Attraos enhances predictive performance and provides empirical and theoretical explanations. In the future, we will verify our proposed Attraos on largescale chaotic datasets and utilize the implicit neural representations for the phase space coordinates, enabling a trainable embedding process to address the instability of PSR techniques. 
experimental setting/details: 
 Question: Does the paper specify all the training and test details (e.g., data splits, hyperparameters, how they were chosen, type of optimizer, etc.) necessary to understand the results?Answer: [Yes] .Justification: We have provided detailed experimental settings in Appendix D to facilitate reproducibility.Guidelines:• The answer NA means that the paper does not include experiments.• The experimental setting should be presented in the core of the paper to a level of detail that is necessary to appreciate the results and make sense of them. • The full details can be provided either with the code, in appendix, or as supplemental material. 
experiment statistical significance: 
 Question: Does the paper report error bars suitably and correctly defined or other appropriate information about the statistical significance of the experiments?Answer: [No] .Justification: In time series forecasting (TSF) tasks, it is common not to provide error bars but instead directly calculate the average by conducting multiple experiments.Guidelines:• The answer NA means that the paper does not include experiments.• The authors should answer "Yes" if the results are accompanied by error bars, confidence intervals, or statistical significance tests, at least for the experiments that support the main claims of the paper.• The factors of variability that the error bars are capturing should be clearly stated (for example, train/test split, initialization, random drawing of some parameter, or overall run with given experimental conditions). • The method for calculating the error bars should be explained (closed form formula, call to a library function, bootstrap, etc.) • The assumptions made should be given (e.g., Normally distributed errors). • It should be clear whether the error bar is the standard deviation or the standard error of the mean. Justification: In Conclusion, we state that our goal is to offer the machine-learning community a fresh perspective and inspire further research on the essence of time series dynamics. 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 networks 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). 
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 datasets chosen in this paper are commonly used benchmark datasets for time series forecasting tasks.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: Yes, we have.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 the assets?Answer: [Yes] .Justification: This paper follows CC 4.0, and the code is in an anonymized URL.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. • At submission time, remember to anonymize your assets (if applicable). You can either create an anonymized URL or include an anonymized zip file. 
crowdsourcing and research with human subjects: 
 Question: For crowdsourcing experiments and research with human subjects, does the paper include the full text of instructions given to participants and screenshots, if applicable, as well as details about compensation (if any)?Answer: [NA] .Justification: This paper does not involve crowdsourcing nor research with human subjects.Guidelines:• The answer NA means that the paper does not involve crowdsourcing nor research with human subjects. • Including this information in the supplemental material is fine, but if the main contribution of the paper involves human subjects, then as much detail as possible should be included in the main paper. • According to the NeurIPS Code of Ethics, workers involved in data collection, curation, or other labor should be paid at least the minimum wage in the country of the data collector. 
institutional review board (irb) approvals or equivalent for research with human subjects: 
 Question: Does the paper describe potential risks incurred by study participants, whether such risks were disclosed to the subjects, and whether Institutional Review Board (IRB) approvals (or an equivalent approval/review based on the requirements of your country or institution) were obtained?Answer: [NA] .Justification: This paper does not involve crowdsourcing nor research with human subjects.Guidelines:• The answer NA means that the paper does not involve crowdsourcing nor research with human subjects. • Depending on the country in which research is conducted, IRB approval (or equivalent) may be required for any human subjects research. If you obtained IRB approval, you should clearly state this in the paper.• 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: A.1 Takens TheoremTakens' theorem, introduced by Floris Takens in 1981, provides a framework for reconstructing the dynamics of a chaotic system based on a series of observations. The theorem establishes conditions under which the state space of a dynamical system can be reconstructed from time-delay measurements of a single observable. The reconstructed system retains properties that are invariant under smooth transformations, known as diffeomorphisms.In the context of discrete-time dynamical systems, Takens' theorem is commonly applied. Consider a dynamical system whose state space is a v-dimensional manifold M, with the evolution of the system governed by a smooth mapF : M → M.Assume there exists a strange attractor A ⊂ M with a box-counting dimension d A . According to Whitney's embedding theorem [37] , the attractor A can be embedded into an m-dimensionalEuclidean space if m > 2d A .This implies that there is a diffeomorphism φ mapping the attractor A into R N , such that the Jacobian matrix of φ is of full rank. In the delay embedding process, the embedding function is constructed using an observation function. This observation function h : M → R needs to be twice continuously differentiable and must assign a real number to each point on the attractor A, ensuring that it is typical, meaning its derivative is of full rank without exhibiting any special symmetries.The delay embedding theorem states that the mappingφ T (x) = h(x), h(F(x)), . . . , h F k-1 (x)constitutes an embedding of the strange attractor A into R N . Now, consider a d-dimensional state vector x t , which evolves according to an unknown but deterministic and continuous dynamic process. Suppose a one-dimensional observable y exists, which is a smooth function of x and is coupled to all the components of x. At any given time, we can observe not only the current measurement y(t) but also measurements from past times separated by a lag τ : y t+τ , y t+2τ , and so on. Using m such lags results in an m-dimensional vector. As the number of lags increases, the motion in the reconstructed space becomes more predictable, and in the limit as m → ∞, the dynamics could become deterministic. In reality, the deterministic nature of the dynamics is achieved at a finite dimension, with the reconstructed dynamics being equivalent to the original system's dynamics, related by a smooth, invertible change of coordinates (a diffeomorphism).The theorem specifically asserts that deterministic behavior emerges once the dimension reaches 2d + 1, with the minimal embedding dimension often being lower.A Hopfield network is a form of recurrent artificial neural network with binary neurons. It is characterized by:• Binary neurons with states +1 or -1.• Symmetric weight matrix with zero diagonal (no self-connections).• Energy function that is minimized at stable states.• Asynchronous update of neuron states.The state of each neuron is updated according to:EQUATIONWhere s i (t + 1) is the state of neuron i at time t + 1, w ij is the weight between neurons i and j, and s j (t) is the state of neuron j at time t.The energy of the network is defined as:EQUATIONMemories are stored in the network by adjusting the weights to minimize the network energy, often using the Hebbian learning rule. EQUATIONwith lse(log -sum -exp) interaction function:EQUATIONSering 0E2M 2 . Using p = softmax βX T ξ , The novel update rule is:EQUATION)The new update rule can be viewed as the attention in Transformers. Firstly, N stored (key) patterns y i and S state (query) patterns r i that are mapped to the Hopfield space of dimension d k . Then set T and R = (r 1 , . . . , r S ) T combine the y i and r i as row vectors. By defining the matricesx i = W T K y i , ξ i = W T Q r i ,X T = K = Y W K , Ξ T = Q = RW Q , and V = Y W K W V = X T W V , where W K ∈ R dy×d k , W Q ∈ R dr×d k , W V ∈ R d k ×dv , β = 1/ √ d k and softmax ∈ R N ischanged to a row vector, the update rule multiplied by W V is:Z = softmax 1/ d k QK T V = softmax βRW Q W T K Y T Y W K W V .In the Hopfield Evolution strategy of Attraos, The Query is settled as the dynamic structures and the Key/Value is settled as trainable vectors.Note: In this section, we have selectively extracted key content from the appendix of Hippo [11] that is pertinent to our work, for the convenience of the reader. We take Legendre polynomials as an example.Under the usual definition of the canonical Legendre polynomial P n , they are orthogonal with respect to the measureω leg = 1 [-1,1] : 2n + 1 2 1 -1 P n (x)P m (x)dx = δ nmAlso, they satisfyP n (1) = 1 P n (-1) = (-1) n .With respect to the measure 1 θ I[t -θ, t], the normalized orthogonal polynomials are(2n + 1) 1/2 P n 2 x -t θ + 1In general, the orthonormal basis for any uniform measure consists of (2n + 1)1 2 times the corresponding linearly shifted version of P n .We note the following recurrence relations on Legendre polynomials:(2n + 1)P n = P ′ n+1 -P ′ n-1P ′ n+1 = (n + 1)P n + xP ′ nThe first equation yields P ′ n+1 = (2n + 1)P n + (2n -3)P n-2 + . . . , where the sum stops at P 0 or P 1 . These equations directly implyP ′ n = (2n -1)P n-1 + (2n -5)P n-3 + . . . and (x + 1)P ′ n (x) = P ′ n+1 + P ′ n -(n + 1)P n = nP n + (2n -1)P n-1 + (2n -3)P n-2 + . . .To sum up, The Legendre polynomials are in closed-recursive form.A = diag{-1, -1, . . . } is a rough approximation of shifted Hippo-LegT [11] matrix.Proof. Hippo 3 provides a mathematical framework for deriving the AB matrix for polynomial projection. The mainstream SSMs [10] typically initialize matrix A = diag{-1, -2, -3, -4, . . . }, representing the negative real diagonal elements of the normalized Hippo-LegS matrix (21) . Since the LegS matrix is a mathematical approximation of exponentially decaying Legendre polynomials, initializing A = diag{-1, -2, -3, -4, . . . } can be seen as a rough approximation of exponential decay. Similarly, for the normalized Hippo-LegT matrix (22) , which approximates a uniform measure, we can consider A = diag{-1, -1, . . . } as a rough approximation of a finite window of Legendre polynomials, which better for non-stationary time series as well [19] . The use of negative values for the elements is to ensure gradient stability during training.EQUATIONEQUATIONTheorem 6. (Approximation Error Bound) Suppose that the function f : [0, 1] ∈ R is k times continuously differentiable, the piecewise polynomial g approximates f with mean error bounded as follows:∥f -g∥ ≤ 2 -rk 2 4 N k! sup x∈[0,1] f (k) (x) .Proof. Similar to [1] , we divide the interval [0, 1] into subintervals on which g is a polynomial; the restriction of g to one such subinterval I r,l is the polynomial of degree less than k that approximates f with minimum mean error. Also, the optimal g can be regarded as the orthonormal projection Q N r f onto G N (r) . We then use the maximum error estimate for the polynomial, which interpolates f at Chebyshev nodes of order k on I r,l . We define I r,l = [2 -r l, 2 -r (l + 1)] for l = 0, 1, . . . , 2 r -1, and obtainQ N r f -f 2 = 1 0 Q N r f (x) -f (x) 2 dx = l I r,l Q N r f (x) -f (x) 2 dx ≤ l I r,l C N r,l f (x) -f (x) 2 dx ≤ l I r,l 2 1-rk 4 N k! sup x∈I r,l f (k) (x) 2 dx ≤ 2 1-rk 4 N k! sup x∈[0,1] f (k) (x) 2and by taking square roots we have bound (7) . Here C N r,l f denotes the polynomial of degree k which agrees with f at the Chebyshev nodes of order k on I r,l , and we have used the well-known maximum error bound for Chebyshev interpolation.The error of the approximation Q N r f of f therefore decays like 2 -rk and, since S N r has a basis of 2 r k elements, we have convergence of order k. For the generalization to m dimensions in the dynamic structure modeling, a similar argument shows that the rate of convergence is of order k/m. Theorem 7. (Evolution Error Bound) By Jacobian value, the mean attractor evolution error K • Ã -Ã of evolution operator K = {K (i) } is bounded byK • Ã -Ã (N -1)E (-β∇ i + 1)with the number of random patterns N ≥ √ pc d-1 4 stored in the system by interaction paradigm E in an ideal spherical retrieval paradigm.Proof. Due to the numerous assumptions and extensive lemmas involved in the proof of this theorem, we will provide a brief exposition of its main ideas. For a detailed proof, please refer to [38] .Firstly, the theorem defines the matching between patterns as: Definition 1. (Pattern match). Assuming that around every pattern x i a sphere S i is given. We say x i is matched with ξ if there is a single fixed point x * i ∈ S i to which all points ξ ∈ S i converge.As shown in Theorem 3 in [38] , according to the upper branch of the Lambert W function [35] , we can obtain the number of random patterns stored in a system is N ≥ √ pc d-1 4 . In our paper, we define the pattern1 as the Attractor Ã in phase space, pattern2 as the future state evaluated by operator K, noted as K • Ã. From Lemma A4 in [38] , when the radius of the pattern matching sphere is M , we can describe the evolution process by jacobian value and get the matching error as:K • Ã -Ã 2ϵMwhere in the Equation (179) of [38] :ϵ = (N -1) exp -β ∇ i -2 max K • Ã -Ã , Ã * i -Ãi M .The Equation (404) of [38] says Ã * i -Ãi1 2βM and K • Ã -Ã 1 2βM, so we can get: ϵe(N -1)M exp (-β∇ i ) .In our paper, we replace the exponential interaction function with an unknown function E(•) to finally get the version in Theorem 3. Theorem 8. (Completeness in L 2 space) The orthonormal system B k = {ϕ j : j = 1, . . . , k}∪{ψ rl j : j = 1, . . . , k; r = 0, 1, 2, . . . ; l = 0, . . . , 2 r -1} spans L 2 [0, 1].Proof. We define the space G N to be the union of the G N (r) , given by the formula:EQUATIONand observe thatG N = L 2 [0, 1].In particular, G N contains the Haar basis for L 2 [0, 1], consisting of functions piecewise constant on each of the subintervals (2 -r l, 2 -r (l + 1)). Here the closure G N is defined with respect to the L 2 -norm,∥f ∥ = ⟨f, f ⟩ 1/2where the inner product ⟨f, g⟩ is defined by the formula⟨f, g⟩ = 1 0 f (x)g(x)dx.Also, we have:EQUATIONand S N r = linear span {ψ l j,r : ψ l j,r (x) = 2 r/2 ψ j (2 r x -l) , j = 1, . . . , k; n = 0, . . . , 2 r -1} . (25)We let {ϕ 1 , . . . , ϕ k } denote an orthonormal basis for G N 0 ; in view of Equation 23, 24, and 25, the orthonormal system B k = {ϕ j : j = 1, . . . , k} ∪ h l j,r : j = 1, . . . , k; r = 0, 1, 2, . . . ; l = 0, . . . , 2 r -1} spans L 2 [0, 1]. Now we construct a basis for L 2 (R) by defining, for r ∈ Z, the space GN r by the formula GN r = {f : the restriction of f to the interval (2 -r n, 2 -r (n + 1)) is a polynomial of degree less than k, for n ∈ Z } and observing that the space GN r+1 \ GN r is spanned by the orthonormal setψ l j,r : ψ l j,r (x) = 2 r/2 ψ j (2 r x -l) , j = 1, . . . , k; l ∈ Z .Thus L 2 (R), which is contained in r GN r , has an orthonormal basis ψ n j,m : j = 1, . . . , k; r, l ∈ ZPhase space reconstruction is a crucial technique in the analysis of dynamical systems, particularly in the study of time series data. This method transforms a one-dimensional time series into a multidimensional phase space, revealing the underlying dynamics of the system. The cross-correlation mutual information (CC mutual information) method is a statistical tool used to analyze the dependencies between two different time series. We provide a detailed mathematical description of these methods.Phase space reconstruction involves transforming a single-variable time series into a multidimensional space to unveil the dynamics of the system that generated the data. The technique is based on Takens' Embedding Theorem.Takens' Embedding Theorem allows the reconstruction of a dynamical system's phase space from a sequence of observations. The theorem states that under generic conditions, a map:EQUATIONcan be constructed, where n is the dimension of the original phase space, and m is the embedding dimension, typically m <= 2n + 1 is sufficient to recover dynamics.The most common approach for phase space reconstruction is the time delay embedding method. Given a time series {x(t)}, the reconstructed phase space is:EQUATION)where τ is the time delay, and m is the embedding dimension. In the context of coordinate delay reconstruction in a complete scenario, it involves a parameter called time window, which selectively uses discrete one-dimensional data for phase space reconstruction. In order to maximize the preservation of time series information, the time window parameter is typically set to 1. As a result, this process is reversible, meaning that the original data can be reconstructed without loss of information.The CC mutual information method is used to analyze the dependencies between two time series. It is a measure of the amount of information obtained about one time series through the other.Mutual information I(X; Y ) between two random variables X and Y is defined as:EQUATIONwhere p(x, y) is the joint probability distribution function of X and Y , and p(x) and p(y) are the marginal probability distribution functions.For time series analysis, the mutual information is extended to account for the time-lagged relationships:EQUATION)where τ is the time lag, and y(τ ) represents the time series Y shifted by τ .Time delay (τ ) is an interval used in reconstructing the phase space of a time series. The right choice of τ is essential for revealing the dynamic properties of the system.For a given time series {x t }, calculate the mutual information between the time series and its time-shifted version for different delays τ :EQUATIONwhere p(x t , x t+τ ) is the joint probability distribution, and p(x t ) and p(x t+τ ) are the marginal probability distributions.Plot the mutual information I(τ ) against τ . Choose the τ at the first local minimum of this plot. This represents the delay where the series' points provide maximal mutual information.Embedding dimension (m) is the dimension of the reconstructed phase space. The correct m ensures that trajectories in the phase space do not intersect each other.For each dimension m, calculate the false nearest neighbors (FNN) F (m), which reflects the complexity of the trajectories reconstructed in m-dimensional space:EQUATIONwhere C i (m, r) is the count of points within a distance r from point i in m-dimensional space, and N is the length of the time series.As m increases, F (m) typically increases and reaches a saturation point. Choose the smallest m for which F (m) does not significantly increase, indicating that increasing the dimension does not reveal more information about the dynamics.Our experiments are carried out on Five real-world datasets and two chaos datasets as described below:• ETT2 dataset are procured from two electricity substations over two years. They provide a summary of load and oil temperature data across seven variables. For ETTm1 and ETTm2, the "m" signifies that data was recorded every 15 minutes, yielding 69,680 time steps. ETTh1 and ETTh2 represent the hourly equivalents of ETTm1 and ETTm2, each containing 17,420 time steps. • Exchange3 dataset details the daily foreign exchange rates of eight countries, including Australia, British, Canada, Switzerland, China, Japan, New Zealand, and Singapore from 1990 to 2016. • Weather4 dataset is a meteorological collection featuring 21 variates, gathered in the United States over a four-year span. • Lorenz96 dataset is an artificial dataset. We simulate the data of 30,000 time steps with an initial dimension of 40d according to the Lorenz96 equation and map to 3D by a randomly initialized Linear network to simulate the realistic chaotic time series generated from an unknown underlying chaotic system. More details are in Appendix E.3. • Crypots comprises historical transaction data for various cryptocurrencies, including Bitcoin and Ethereum. We select samples with Asset_ID set to 0, and remove the column Count. We download the data from https://www.kaggle.com/competitions/ g-research-crypto-forecasting/data?select=supplemental_train.csv. • Air Convection5 dataset is obtained by scraping the data from NOAA, and it includes the data for the entire year of 2023. The dataset consists of 20 variables, including air humidity, pressure, convection characteristics, and others. The data was sampled at intervals of 15 minutes and averaged over the course of the entire year.• Lorenz96-3d: See Appendix E.3.Our baseline models include: Mamba4TS: Mamba4TS is a novel SSM architecture tailored for TSF tasks, featuring a parallel scan (https://github.com/alxndrTL/mamba.py/tree/main). Additionally, this model adopts a patching operation with both patch length and stride set to 16. We use the recommended configuration as our experimental settings with a batch size of 32, and the learning rate is 0.0001. [47] : S-Mamba utilizes a linear tokenization of variates and a bidirectional Mamba layer to efficiently capture inter-variate correlations and temporal dependencies. This approach underscores its potential as a scalable alternative to Transformer technologies in TSF. We download the source code from: https://github.com/wzhwzhwzh0921/S-D-Mamba and adopt the recommended setting as its experimental configuration.RWKV-TS [17] : RWKV-TS is an innovative RNN-based architecture for TSF that offers linear time and memory efficiency. We download the source code from: https://github.com/howard-hou/RWKV-TS. We follow the recommended settings as experimental configuration.Koopa [32] : Koopa is a novel forecasting model that tackles non-stationary time series using Koopman theory to differentiate time-variant dynamics. It features a Fourier Filter and Koopman Predictors within a stackable block architecture, optimizing hierarchical dynamics learning. We download the source code from: https://github.com/thuml/Koopa. We set the lookback window to fixed values of {96, 192, 336, 720} instead of twice the output length as in the original experimental settings.iTransformer [31] : iTransformer modifies traditional Transformer models for time series forecasting by inverting dimensions and applying attention and feed-forward networks across variate tokens. We download the source code from: https://github.com/thuml/iTransformer. We follow the recommended settings as experimental configuration.PatchTST [33] : PatchTST introduces a novel design for Transformer-based models tailored to time series forecasting. It incorporates two essential components: patching and channel-independent structure. We obtain the source code from: https://github.com/PatchTST. This code serves as our baseline for long-term forecasting, and we follow the recommended settings for our experiments.LTSF-Linear [52] : In LTSF-Linear family, DLinear decomposes raw data into trend and seasonal components and NLinear is just a single linear models to capture temporal relationships between input and output sequences. We obtain the source code from: https://github.com/cure-lab/LTSF-Linear, using it as our long-term forecasting baseline and adhering to recommended settings for experimental configuration.This study explores the application of pre-trained language models to time series analysis tasks, demonstrating that the Frozen Pretrained Transformer (FPT), without modifications to its core architecture, achieves state-of-the-art results across various tasks. We download the source code from: https://github.com/DAMO-DI-ML/NeurIPS2023-One-Fits-All. We follow the recommended settings as experimental configuration.All experiments are conducted on the NVIDIA RTX3090-24G and A6000-48G GPUs. The Adam optimizer is chosen. A grid search is performed to determine the optimal hyperparameters, including the learning rate from {0.0001, 0.0005, 0.001}, model layer from {1, 2} (typically 1), polynomial order from {32, 64, 256}, patch length from {8, 16}, projection level based on logL and no more than 3, primary modes is 16, and the input length is 96. Due to the sensitivity of the CC method to numerical calculations, we take the average result from three iterations of the calculation. As shown in Figure 6 , we present the phase space structures of various real-time series. Please note that due to our visualization limitations in three dimensions, the shapes of many attractors for these time series are manifested in higher dimensions. We can only display slices of the attractors in the first three dimensions.Lyapunov exponents are a set of quantities that characterize the rate of separation of infinitesimally close trajectories in a dynamical system. They play a crucial role in understanding the stability and chaotic behavior of the system [24] . The formal definition of the Lyapunov exponent for a trajectory of a dynamical system is given by:EQUATIONwhere δX(t) is the deviation vector at time t, and λ is the Lyapunov exponent. The set of Lyapunov exponents for a system provides a spectrum, indicating the behavior of trajectories in each dimension of the system's phase space.When the direction of the initial separation vector is different, the separation rate is also different. Thus, the spectrum of Lyapunov exponents exists, which have the same number of dimensions as the phase space. The largest of these is often referred to as the Maximal Lyapunov exponent (MLE). We reconstruct dimension m according to different phase spaces to get the maximum Lyapunov index. A positive MLE is often considered an indicator of chaotic behavior in the system, while a negative value indicates convergence, implying stability in the system's behavior.As shown in Table 6 , we calculate the maximum Lyapunov exponent for all mainstream LTSF datasets. Surprisingly, we find that their MLEs are all positive, indicating the presence of chaos to varying degrees in these datasets. This directly supports the motivation proposed by Attraos. Among them, the weather dataset exhibits the strongest chaotic behavior. Note that the existence of at least one positive Lyapunov exponent is sufficient to determine the presence of chaos. For example, the classical Lorenz63 system has three Lyapunov exponents with values negative, zero, and positive respectively.Remark E.1. For multivariate time series, we take the average. Lorenz63 SystemThe Lorenz63 system, introduced by Edward Lorenz in 1963, is a simplified mathematical model for atmospheric convection. The model is a system of three ordinary differential equations now known as the Lorenz equations:EQUATIONEQUATIONEQUATIONHere, x, y, and z make up the state of the system. The parameters σ, ρ, and β represent the Prandtl number, Rayleigh number, and certain physical dimensions of the layer, respectively. Lorenz derived these equations to model the way air moves around in the atmosphere. This model is famous for exhibiting chaotic behavior for certain parameter values and initial conditions.The Lorenz96 system is a mathematical model that was introduced by Edward N. Lorenz in 1996 as an extension of the original Lorenz model. It is commonly used to study chaotic dynamics in systems with multiple interacting variables.The Lorenz96 system consists of a set of ordinary differential equations that describe the time evolution of a set of variables. In its simplest form, the system is defined as follows:dx i dt = (x i+1 -x i-2 )x i-1 -x i + FHere, x i represents the state variable at position i, and F is a forcing term that controls the overall behavior of the system. The nonlinear term (x i+1 -x i-2 )x i-1 captures the interactions between neighboring variables, leading to the emergence of chaotic behavior. The Lorenz96 system exhibits a range of fascinating phenomena, including intermittent chaos, phase transitions, and the presence of multiple stable and unstable regimes. Its dynamics have been extensively studied to gain insights into the behavior of complex systems and to explore the limits of predictability. To simulate time series generated from an unknown chaotic system, we first construct a 40dimensional Lorenz96 system to represent the underlying chaotic system. We then use a randomly initialized linear neural network to simulate the observation function h. Through h, we map the Lorenz96-40d system, which resides in the manifold space, to a 3-dimensional Euclidean space to obtain a multivariate time series dataset with three variables. As shown in the middle of Figure 7 , in the observation domain, the dynamical system structure of Lorenz96-40d is difficult to discern specific shapes, so we need to reconstruct it to 40 dimensions using the Phase Space Reconstruction (PSR) technique to study its dynamic characteristics in a topologically equivalent structure. It is worth noting that on the right side of Figure 7 , we visualize the last variable of the system and find striking similarities with real-world time series, even exhibiting some periodic behaviors. This further supports our hypothesis that real-world time series are generated by underlying chaotic systems. In LTSF tasks, it is commonly observed that predictive performance deteriorates as the length of the forecasting window increases. This phenomenon bears a striking resemblance to the inherent challenge in long-term predictions of chaotic dynamical systems, which are highly sensitive to initial conditions. Recent studies [15] in the field of chaotic dynamical systems have highlighted that, to address the issue of gradient divergence caused by positive maximal Lyapunov exponents indicative of chaos, the implementation of teaching forcing as a method to incrementally constrain the trajectory of the dynamical system presents a straightforward yet effective framework. To enhance Hippo-RNN, we have implemented the following modifications: zt := (1 -α)z t + αz t , z t = F θ ( zt-1 ) , with 0 ≤ α ≤ 1, where we intervene the evolution state z by utilizing the ground truth hidden state z. As shown in Figure 8 , this modification leads to notable improvements in both mainstream LTSF datasets and real-world chaotic datasets. However, the autoregressive nature of the teaching forcing methods introduces additional computational overhead and may potentially reduce its generalization capabilities, which presents a challenge for integrating it into sequence-mapping models.As shown in Table 7 , we showcase the full experiment results for all mainstream LTSF datasets. NeurIPS Paper Checklist Justification: Attraos, considers time series as generated by a generalized chaotic dynamic system. By leveraging the invariance of attractors, Attraos enhances time series prediction performance and provides empirical and theoretical explainations.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.https://github.com/zhouhaoyi/ETDatasethttps://github.com/laiguokun/multivariate-time-series-datahttps://www.bgc-jena.mpg.de/wetterhttps://www.psl.noaa.gov/ 


Decision: The paper presents a novel model for time-series forecasting building on a chaotic dynamical systems view: the approach is original and robust, both empirically and theoretically. This has been widely recognized by the Reviewers. 
Issues were mainly associated with presentation aspects and they have been appropriately addressed and clarified through the rebuttal and the discussion period. 
The work can therefore be accepted pending that the Authors thoroughly revise and rewrite Section A.1 which appears to be a verbatim (and unacknowledged) copy of the Wikipedia page on Takens' theorem.

Reviews: review: paper_summary: paper_summary: The paper introduces a novel approach, named Attraos, for Long-term Time Series Forecasting (LTSF) based on treating the observed time series as high dimensional chaotic dynamical system. The model first estimates an embedding of the data through Phase Space Reconstruction (Takens embedding) and then utilizes a memory unit through state space model that represents the dynamical system through polynomials which allows to evolve the dynamics into future steps and forecast the data. The model is evaluated on multiple LTSF datasets and additional experiments such as ablations and robustness under noise are performed.
strength_weakness: strength_weakness: W1. There seems to be a detachment between the theorems and properties proved and the proposed system. There seems to be no discussion about how these properties lead to the particular setup in the paper. For example why this particular SSM was used? How the hyperparameters were chosen in light of the propositions? What are the cases that the model is limited and not expected to be effective? Also interpretation of the results is lacking.

W2. For non toy experimental results exposition and interpretation of the resulting dynamical system is missing. What is the dynamical system/s that are being obtained for the benchmark and what is the distribution of polynomials?
questions: questions: See questions in W1 and W2 and also:

Q1. The Phase Space Reduction approach seems to be the well known Takens embedding and authors mention Takens thm, however call the method differently. Is there something different that I'm missing in PSR that is not in the classical computation of Takens embedding?

Q2. With local evolution I assumed that authors mean single step "autoregressive" forecasting. Is that the case? How would that change for multi-step?

Q3. It would be informative if the benefit/limitation/motivation of frequency evol vs. direct evol be more elaborately explained.
limitations: limitations: See W1-W2 and Q1-Q3. More extended discussion of limitations would contribute to better evaluation of the work vs generalization.
paper_summary: paper_summary: The paper introduces chaos theory into a long-term time series forcasting (LTSF) model called Attraos (a play on the words "attractor" and "chaos"). They propose a Multi-resolution Dynamic Memory Unit (MDMU) which is inspired by (and looks a lot like) the State Space Models (SSMs) used in the Mamba family of models. However, unlike SSMs, Attraos assumes additional structure on top of time-series signals by utilizing attractor behavior from chaos theory in the phase space. Attraos outperforms many existing architectures (e.g., Mamba, RWKV-TS, PatchTST) on time-series tasks at a fraction of the computational complexity and cost.
strength_weakness: strength_weakness: - W1: **Presentation of results could be clearer**. I was often confused by the choice of coloring in the tables. See the following for my suggestions:
    1. [Table 1,2,4]. The "Red/blue" colorscheme is unfortunate, because most brains (including mine) associate the color "red" with "poor performance". I suggest using **bold** for best, <u>underline</u> for second best, as is often done in AI papers. 
    2. [Table 3,5]. For "improvement/decline" in performance, I suggest one of two options: (1) show "better performance" with blue, and "worse performance" with bold and red to emphasize decreased performance from an ablation, or (2) showing only the "deltas" with a "stock ticker arrow" (🔺 or 🔻, colored green for improvement, red for decline).
- W2: **Complete architectural description is missing**. Most of the paper focuses on describing how to implement the special components of Attraos, but the complete picture of the architecture is missing. Additionally, there is no section describing the hyperparameters for the training setup for each experiment. See Q1 for specific questions.
- W3: **Missing error bars**. Table 1 represents average results of long term forecasting across a swath of model classes, but these averages are not accompanied by reports of standard deviation. Thus, it is difficult to tell whether improved average performance is actually signficant. Improve by including error bars for all numerical results, which may require re-running some experiments under different random seeds.
- W4: **Grandiose, non-academic language**. The paper starts with "In the intricate dance of time" [L15] 😂 and makes statements like "we can transcend the limitations of deterministic dynamical systems" [L34]. I suggest reworking these lines to preserve the professionalism of the rest of the work.

These weakenesses are admittedly small and easily fixable by the authors during the review process. My overall attitude towards the paper is that it is of high quality and should be accepted. However, a lot of the theory was beyond my ability to evaluate, and I am willing to adjust my score positively/negatively as I become more familiar with this work throughout the review process.
questions: questions: - Q1: The complete architecture for Attraos is not described. Do you just take Mamba and modify the SSM component? What is the complete set of learnable parameters? Where does the Hopfield Network come into play and how are its memories trained? What are the training hyperparameters for each experiment?
limitations: limitations: Limitations are adequately described in the conclusion.
paper_summary: paper_summary: This paper introduces Attraos, a new model for long-term time series forecasting (LTSF) that incorporates chaos theory and views time series data as observations from high-dimensional chaotic dynamic systems. Attraos utilizes attractor invariance, non-parametric Phase Space Reconstruction, and a multi-scale dynamic memory unit to effectively capture historical dynamics and forecast future states with substantially fewer parameters than existing models such as PatchTST. Empirical evidence demonstrates Attraos' superior performance across various LTSF and chaotic datasets, providing a new perspective on the underlying dynamics of time series data.
strength_weakness: strength_weakness: -	Most of the descriptions for problem setup (section 2), methods, and theoretical analysis (section 3) are very confusing and hard to follow. There are always notations suddenly coming out without much information. I suspect there are also several typos that cause difficulties in understanding and assessments. See **Questions** for details.
-	If I understood correctly, the title is a little misleading since seemingly the key insight/perspective comes from the embedding theorem by Takens, a counterpart of Whitney theorem for attractors. Are there any other insights related to chaotic systems? If not, I would suggest concretizing this point in the title, abstract, introduction, and conclusion.
-	Section A.1 in appendix is entirely copied from Wikipedia, https://en.wikipedia.org/wiki/Takens%27s_theorem , except for some changes on notations. However the notations in this section are not coherent, e.g. $k$ and $N$ in L414.
-	Similar to the first point, the paper does not provide enough background introduction.

If the authors worry containing all the necessary details will exceed the page limit, I would suggest presenting a detailed version of problem setup, backgrounds, analysis and methods in appendix.
questions: questions: 1.	L90 Is it a weighted integral $\mu(ds)$ instead of $ds$? For polynomials $\phi$, the integral with Lebesgue measure does not converge. 
2.	Is $\mu$ same as $\omega$ in L94?
3.	Isn’t L94 giving the formula of projection to polynomial subspace and isn’t $K_n$ the kernel? How are they connected to $e^{tA}B$, where $A, B$ only appear in the induced dynamics of $x(t)$?
4.	Does $\mathcal{A}_i$ refer that there are multiple attractors in the dynamics? $\mathcal{A}$ are sets, then what does transpose in L158 stand for?
5.	L112, 120, if $A$ is a (D,N) tensor, $B$ is a (B,L,N) tensor and $u$ is a $(B,L,D)$ tensor, how is $Ax(t)+Bu(t)$ in (2a) defined? 
6.	Could the authors give more interpretation on what they intend to do $\Delta$ and $\theta$? I have no idea what the authors mean by ‘$\Delta$ is similar to an attention mechanism’ and how this is related to the previous context.
7.	Again, given the shape of $B$ and $\Delta$, it is confusing what $\Delta B$ refers to.
8.	How is $H$ in L144 related to previous text?
limitations: limitations: Yes.

Comments: title: Authors Rebuttal Response
Kindly Request for Reviewer's Feedback
comment: Thank you for your response and support! We would be delighted to provide you with further clarification on the theory and approach choices.

* **Motivation and Significance**: Machine learning has historically been bifurcated into physics-driven approaches (such as PDE numerical solvers, PINN, etc.) and data-driven methods (like large models, various neural networks, neural operators, etc.). While data-driven techniques have excelled in time series modeling (as well as in computer vision, natural language processing, etc.), their black-box nature poses challenges for sustained advancement within the machine learning community. Recent publications in esteemed journals like Nature, Science, PNAS, etc., have showcased a fusion of physics-driven and data-driven methodologies, hinting at the emergence of a prominent research trend involving the integration of physics priors into deep learning models. **Our research marks the inaugural incorporation of PDE dynamical system knowledge into the realm of time series prediction, promising to facilitate interdisciplinary amalgamation and the creation of interpretable deep time series models.**
* **Dynamic System**: also known as PDE dynamical systems. Let the domain $S$ be an open subset of $\mathbb{R}^d$ and set an integer $k \geq 1$. Define the system state as $\boldsymbol{x}: S \mapsto \mathbb{R}^m$ where $\boldsymbol{x}=\left(x^1, \ldots, x^m\right)$. Then, an expression of the form:$\mathcal{F}\left(D^k \boldsymbol{x}(s), D^{k-1} \boldsymbol{x}(s), \ldots, D \boldsymbol{x}(s), \boldsymbol{x}(s), s\right)=0$ is called a $k^{\text {th }}$-order system of partial differential equation, where $\mathcal{F}: \mathbb{R}^{m d^k} \times \mathbb{R}^{m d^{k-1}} \times \ldots \times \mathbb{R}^{m d} \times \mathbb{R}^m \times S \mapsto \mathbb{R}^m$ and $s \in S$.

real-time series are observed values obtained from the system through a function h. By starting from dynamical systems instead of the original value, it becomes possible to have a better understanding of the fundamental temporal behavior. Moreover, this approach offers the advantages of interpretability and visualization.

* **Chaos Theory**: The study of chaos can be considered a branch of dynamical systems, focusing on the fact that both linear and nonlinear dynamical systems tend to exhibit certain fixed shapes in their trajectories, known as attractors. What may appear as irregular behavior in the time domain often reveals stable structures in the dynamical trajectories. The concept of attractors can readily be correlated with current deep-learning pattern recognition technologies.

* **Why SSM?**: The reason for choosing SSM to encode dynamical structures is that polynomials are widely recognized as universal approximators in dynamical systems research. The mathematical interpretation of SSM conveniently aligns with polynomial projection, making SSM well-suited for encoding dynamical structures. And we have introduced a novel SSM evolution matrix {-1, -1, -1, ...} to describe finite window approximations.

* **Why SSM + Multiple Orthogonal Subspaces**: Building upon SSM, we drew inspiration from chaos-related research (**Theorem 3 Attractor evolution error, line 163**) and proposed an enhanced version of SSM that utilizes distinct orthogonal subspaces to store various attractor structures. This innovation was validated through experimental results.

* **Why Frequency Evolution**: In the field of neuroscience research (where EEG and ECG are recognized as chaos datasets), it has been established that attractors are amplified in the frequency domain. Therefore, this paper adopts a frequency-domain evolution strategy. Besides, by discarding high-frequency components, noise effects can be mitigated, reducing computational complexity. Additionally, we explored the effects of Hopfield evolution and KNN clustering evolution, and through experimentation, **we found that these two methods were not as effective as frequency-domain evolution. We provided explanations for these findings in line 257**.

* **Follow-up Work**: In our recent research, we have discovered that **incorporating dynamical system priors leads to performance enhancements across all four types of time-series tasks (prediction, classification, interpolation, anomaly detection) and all four neural network architectures (convolutional, attention-based, linear, SSM). These improvements include a tenfold reduction in parameter count, as well as more stable gradients.** This work is set to be released soon—stay tuned for more updates.
I would like to thank the authors for their concise yet informative rebuttal that have clarified most of my questions. After reading authors rebuttal to my review and other reviewers I remain positive and supportive of the work and novelty of the approach.

One issue that I am still wondering about and would like to know if authors are able to further clarify is the detachment of the theory and the approach choices. Could the authors summarize the outcomes of the properties developed in the work and how they led/constrained model choices?
Dear Reviewer zUbZ,

Since the End of author/reviewer discussions is coming in one day, may we know if our response addresses your main concerns? If so, we kindly ask for your reconsideration of the score. Should you have any further advice on the paper and/or our rebuttal, please let us know and we will be more than happy to engage in more discussion and paper improvements.

Thank you so much for devoting time to improving our paper!
Thank you for your feedback and validation of our paper! We greatly appreciate your help in improving our scores! Wishing you all the best!
Thank you for your response! We greatly appreciate your support in recommending our paper!  Wishing you all the best!
Thanks for the clarification. Hope the authors will polish their writing in the final version.
In this year's review cycle I am unable to see your paper's revisions, but I trust the authors that my suggested changes have been made. Similarly, I thank the authors for answering my questions. After reading the other reviews and responses, I see no reason to decrease my score and keep my vote to accept this paper.