LINOCS: Lookahead Inference of Networked Operators for Continuous Stability
Authors: Noga Mudrik
Eva Yezerets
Yenho Chen
Christopher John Rozell
Adam Shabti Charles
Venue: Accepted by TMLR ()

Abstract: Identifying latent interactions within complex systems is key to unlocking deeper insights into their operational dynamics, including how their elements affect each other and contribute to the overall system behavior. For instance, in neuroscience, discovering neuron-to-neuron interactions is essential for understanding brain function; in ecology, recognizing the interactions among populations is key for understanding complex ecosystems. Such systems, often modeled as dynamical systems, typically exhibit noisy high-dimensional and non-stationary temporal behavior that renders their identification challenging. Existing dynamical system identification methods often yield operators that accurately capture short-term behavior  but  fail to predict long-term trends, suggesting an incomplete capture of the underlying process. Methods that consider extended forecasts (e.g., recurrent neural networks) lack explicit representations of element-wise interactions and require substantial training data, thereby failing to capture interpretable network operators. Here we introduce Lookahead-driven Inference of Networked Operators for Continuous Stability (LINOCS), a robust learning procedure for identifying hidden dynamical interactions in noisy time-series data. LINOCS integrates several multi-step predictions with adaptive weights during training to recover dynamical operators that can yield accurate long-term predictions. We demonstrate LINOCS' ability to recover the ground truth dynamical operators underlying synthetic time-series data for multiple dynamical systems models (including  linear, piece-wise linear, time-changing linear systems' decomposition, and regularized linear time-varying systems) as well as its capability to produce meaningful operators with robust reconstructions through various real-world examples.

Published in Transactions on Machine Learning Research (09/2024)
LINOCS: Lookahead Inference of Networked Operators for Continuous Stability
Noga Mudrik
nmudrik @rjhu. edu
Biomedical Engineering, Kavli NDI, CIS The Johns Hopkins Universily Baltimore, MD, 21218.
Eva Yezerets
eyezerel @fjhu edu
Biomedical Engineering, CIS
The Johns Hopkins University Baltimore, MD, 21218
Yenho Chen
yenho&ggatech.co edu
Department of Biomedical Engineering Georgia Institute of Technology Atlanta, GA 30392.
Christopher J. Rozell
crozell@gatech. edu
School of Electrical and Computer Engineering Georgia Institute of Technology
Atlanta, GA 30332.
Adam S. Charles
adamsc edu
Biomedical Engineering, Kaoli NDI, CIS The Johns Hopkins Universily Baltimore, MD, 21218.
Reviewed on OpenReview: https: // openreview. net/ forum? id- A6D3PYSygJ
Abstract
Identifying hidden interactions within complex systems is key to unlocking deeper insights into their operational dynamics, including how their elements aflect each other and con- tribute to the overall system behavior. For instance, neuroscience, discovering neuron- to-neuron interactions is essential for understanding brain function; in ecology, recognizing interactions among populations is key to understanding complex ecosystems. Such sys- lems, often modeled dynamical systems, typically exhibit noisy high-dimensional and non-stationary temporal behavior that renders their identification challenging. Existing dynamical system identification methods typically yield operators that accurately capture short-term behavior but fail to predict long-term trends, suggesting an incomplete capture of the underlying process. Methods that consider extended forecasts (e.g., recurrent neural networks) lack explicit representations of element-wise interactions and require substantial training data, thereby failing to capture interpretable network operators. Here introduce Lookahead-driven Inference of Networked Operators for Continuous Stability (LLNOCS), robust lcarning procedure for identifying hidden dynamical interactions in noisy time- series data. LINOCS integrates several multi-step predictions with adaptive weights during training to recover dynamical operators that can yield accurate long-term predictions, We demonstrate LINOCS' ability recover the ground truth dynamical operators underlying synthetic time-series data for multiple dynamical systems models (including linear, pieco- wise linear, time-changing linear systems' decomposition, and regularized linear time-varying

---

Published in Transactions on Machine Learning Research (09/2024)
systems) as well as its capability to produce meaningful operators with robust reconstruc- tions through various real-world examples.
1 Introduction
Uncovering the dynamics underlying high-dimensional time-series data is crucial for deciphering the fun- damental principles that govern temporally evolving systems. This is apparent across significant scientific domains, including neuroscience (where neurons or ensembles interact over time, e.g., Vaadia et al. (1995):
D'Aleo et al. (2019); Mudrik et al. (2024ch), immunology (where cells regulate immune responses, e-g-, [Svvll] et al. (002)] and ecology (where understanding population interactions yiclds insights into ecosystem dy- namics, e.g. Stein et al. (2013)). Hence, scientific research necessitales the development of procedures adept at learning dynamic operators that can accurately capture the non-linear and non-stationary evolution of real systems.
Existing approaches for dynamical systems identification; though, often rely on either "black-box" deep learning methods, which while powerful, yield uninterpretable representations, or on simple learning pro- cedures that maximize reconstruction between consccutive samplos, and thus fail to accurately predict the system's behavior for longer scales. Specifically, common dynamical system identification models regularly rely on optimizing dynamics by minimizing the prediction error for each time point based on projecting the preceding one through dynamics. Consequently, when using such procedures to learn the operators, post-learning long-term predictions of the system's values (by iteratively estimating the system's state at the next time point) usually result in undesired divergence away from the real system's 's values. This difficulty in long-term predictions, importantly, implics that operators recovered by these models that are based on local cost functions may not capture the underlying system correctly. The challenge in identifying such un- derlying operators, therefore, lies in the need to incorporate long-term predictions directly into the learning procedure, which can be especially challenging in cases where the dynamics are non-stationary, non-linear, or otherwise constrained in ways that reflect real-world system behavior.
To address this challenge, we present a learning procedure that introduces Lookahead Inference of Networked Operators for Continuous Stability (LINOCS). LINOCS bridges the gap between minimizing reconstruction costs based on single time-step projections (which typically result operators that quickly diverge in long- term forecasts), and optimizing multi-step training that relies on reconstructions from past time points (which can lead to unstable predictions). Particularly, LINOCS achieves this by integrating adaptive re- re-weighted multi-step reconstructions into the dynamics inference, progressively building up the cost over training iterations while simultaneously considering several multi-step reconstruction terms to identify operators that enable stable, long-term reconstruction post-training. Through this process, LINOCS also avoids relying on massive amounts of data (like other methods, including RNNs,, for example, require). We demonstrate the effectiveness and adaptability of LINOCS in a variety of dynamical systems models, including linear, switching-linear, decomposed systems, and smoothly linear time-varying systems, achieving significantly improved accuracy in operator identification and long-term predictions compared to 1-step optimization approaches (Fig. 4 34).
The main advantages of LINOCS over existing approaches is its integration of adaptive weights for dif- ferent reconstruction orders. This adaptive nature allows LINOCS to stably regulate noisy measurements and extract highly accurate dynamics that can predict far into the future. Moreover, while other meth primarily focus on specific dynamical system architectures (e.g., RNNs), LINOCS can be implemented into various dynamic structures, as we demonstrate through applications to switching linear systems and decom- posed linear systems. Finally, LINOCS improves the identification of operators that expose pairwise element interactions, which is important for scientific interpretation.
Our contributions in this paper notably include:
We propose LINOCS, a novel learning procedure that incorporates re-weighting multi-step predic- tions into the cost for operator identification.
We demonstrate that applying LINOCS improves the ability to recognize ground-truth operators.
2

---

Published in Transactions on Machine Learning Research (09/2024)
We show LINOCS' efficacy across diverse range of dynamical systems, including linear, periodically linear, linear time-varying (ITTV), and decomposed linear.
Finally, we demonstrate LINOCS' ability to work real-world brain recordings, resulting in better long-term reconstruction compared to baselines.
2 Background and Terminology
Consider system with p interacting elements (e.g.; neurons in the brain) whose time-changing state X E RexT evolves over discrete time points f == 1... as I++1 = g( b,, t), where I3 E I2 refers to the state at time t, b, E IP represents an offsct at time t, and g is a function g IR x R x Z + Rr. For example, neuroscience, de can represent the time-evolving activily of p recorded neurons during a record- ing session with T time points, or 24 can represent the activation levels of p immune cells when modeling the immune system.
In this paper, we focus on linear, piece-wise linear, and locally linear time-varying systems. Specifically, we limit our analysis to functions g(-) that can be written as:
a++1 g( b., () == A.T: + b.,
(1)
where A, E RPxP represent the transition matrices at each time t. Our focus on locally linear dynamics is supported by the fact that even highly nonlinear functions can be well approximated over small time intervals using local linearization (Khalil, 2001 Sastry 2013). Importantly, this formulation's advantage lies in its "network interprotability" the retention of the ability to casily extract the system's pairwise interactions, including non-stationary changes in A, and b, over time. Specifically, any operator entry [Alky for = =..=1...p in every t 1...f can be interpreted as the effect of element on element at time t.
In practice, however, robustly recovering operators that can accurately describe the system's evolution in non-stationary and non-linear settings, faces numerous computational challenges. Chiefly, if we adopt a naive approach to identify the operators as: A,, b, arg min | A,æ; - bell for every = 1....." this problem is statistically unidentifable. Specifically, the problem has p² +p unknowns for each time point, but only p equations.
One approach improve inference in these settings is introduce additional structure via a prior over A, and b,, that constrains the solution space and is commonly grounded in application-driven assumptions. For instance, in many scientific settings, it is reasonable to assume that interactions change smoothly over time. Therefore, adding a temporal smoothness constraint on A₁ and b, (c.g., ||A, - A < bb 112 < es) can be beneficial for both interpretability and accuracy. In addition, the inclusion of such constraints can be crucial, particularly in noisy settings, to prevent overfitting. The addition such constraints transforms the problem to:
and
Å, b = arg min |0 - A,x; + A,
+ RRA,, b.).
(2)
where R(A,, b,) can represent regularization on the dynamic operators.
Although these solutions may offer good short-term predictions (i.e., predicting æ; from xx-1) with minor errors, they often struggle to fully reconstruet the dynamics over longer time-scales due to the build-up of estimation errors over multiple iterative predictions.
The challenge of operator identification is further complicated by the practical issue that we can typically only observe noisy measurements of æ;;
き = xt + E, where e E T represents some noise, e.g.,, ~ N(
Particularly, if A. is obtained from arg mig || A, the distances between the real and the estimated operators {I[A, AT increase with the noise furthr hindering operator recov- ery. For instance, given 6
o2(0,22) i.dd Gaussian noise, ie., 2, has covariance of
3

---

Published in Transactions on Machine Learning Research (09/2024)
I, then T441_ AT, + c A + c) + A Ae + meaning the new noise term is A + e ((0, (AAT + T)). For higher orders, the accumulated noise variance over K steps will result in: E o² (AT)), which can quickly increase based on the eigenspectra of A.
As the accuracy of the fit of a dynamical system to data often needs to be evaluated based on its ability to prodict future values accurately (Tabar) 2019), an inability to capture long-term prodictions suggests that the learned operators demonstrate limited capacity to fully describe the system.
In this context, we define below three types of prediction styles that will be used to evaluate and compare methods throughout this paper:
1-Step Prediction (xt-1(æ²); 1-Step prediction involves using the state each time point estimate the state at the next time point (t (1+11 1).
Iterative Multi-Step Prediction (IMS) of Order K E N (0++a(æ)): IMS involves iteratively, for K times, forecasting 1-step ahead values and using these forecasts as inputs for further 1-step ahead forecasts. Namely, predicting Ee+x|0e+k- vk= = 0... - 1, where ãt-1 := _t-1- Here, we will notate an IMS prediction of order K' by E We chose to name this prediction style "IMS" as to be consistent with the literature (Chevillon 2007..
Full Lookahead Prediction (w This method enhances IMS by forecasting the state at each time point Dt starting from the initial observations (ão). It achieves this by sequentially applying transition matrices to the estimation from the previous time point @ starting from æo, resulting in: E = A Agão =1
(Note: the formula above is presented without offsets {be)T=o for simplicity, though they may be included).
Importantly, IMS Prediction and Full Lookahead Prediction typically exhibit instability due the accumu- lation of errors in the sequential reconstructing process (Fig. []].
3 Prior relevant approaches:
Thcoretical literature long-term prediction instability traces back to Cox 1961 and Klein (2019). who, respectively, introduced exponential smoothing and direct estimation of distant future states. Subsequent studies, including Findley (1983 1985); Weiss (1991); Tiao and Xu (1993 Lin and Granger (1994): Kang (2003), evaluated the effectiveness of dynamical system identification methods in producing long- term predic- tions. These approaches, however, build on either "I-step training", which focuses on ideutifying dynamical operators by minimizing the reconstruction error when projecting the state from one time point to the next, or on "direct forecasting", which aims to identify a mapping function F&t : ER + R² that predicts future states by lzt* (x)( bypassing explicit identification of intermediate dynamical operators.
While the "direct estimation" approach naturally results in more stable long-term predictions compared to 1-step optimization, such approach fails to provide an interpretable "network" meaning to the operators, i.e. the ability to interpret each entry (i.j) of the operator (A.,(i.j) as the effect of clement j on element at time 4. This important e.g., in neuroscience, where understanding the brain's 's interactions entails discerning the time-changing fast interactions of neurons (Sussillo 2014), or epidemiology where tracking disease spread dynamics matters (Aguiar et al., (020)). In addition, when Marcellino al. (2006) compared between iterated and direct estimates using macroeconomic data, they found that in contrast to previous assumptions, iterated forecast methods outperform direct forecast methods, especially when models can auto-select long-lag specifications raising questions about which approach is more suitable for learning dynamical operators.
Markov and Hidden Markov Models (Florian et al. 2011 Ou et al. 2013 Bilmes 2006 are widely used to model time series data and capture temporal dependencies in dynamical systems. However, HMMs struggle with long-term predictions due their reliance on the Markovian assumption, and require extensions
4

---

Published in Transactions on Machine Learning Research (09/2024)
for modeling non-stationary systems (e.g. Sin and Kim (1995) ). Other models, including low-rank auto-
aggressors (Harris et al. 2021 Basu et al. 2019) and other low-dimensional linear models (e.g.DDM) and its variants Tu (2013); Askham and Kutz (2018); Sashidhar and Kutz (2022); Ferré et al. (2023); Kutz et al. (2016) are widely used for capturing and predicting complex dynamical systems in various fields such as fluid dynamics, neuroscience, and financial modeling. However, they do not specifically address the issue of long-term divergence.
More recently, Venkatraman et al. (2015) proposed general approach called DAD that reuses training data to build no-regret learner with multi-step prediction that includes "fixing" and updating the model itself based on every step within the multi-step prediction. However, the authors presented it as a general, nonspecific, approach to consider without specific implementation details. More importantly, DAD does not discuss the possibility for including priors over the operators (e.g., temporal smoothness between consecutive operators) during the training, nor did they consider the need to find operators that are not only expressive but also interpretable.
Other approaches are based on multiple shooting, include, e.g. Jordana et al 2021 who proposed proposed learning dynamics with multiple shooting and multi-step training under each sub-tra sub-trajectory (i.e., between shooting nodes). However, their method requires matching boundary conditions as they start only from the shooting nodes, which can be prone to noise sensitivity. Additionally, they give similar weights to different multi-step orders, which may be hard to tune some orders present high errors at different periods during training.
Generalized Linear Model (GLM) with multi-step inference were proposed in the context of learning neural spikes from spiking data (Hocker and Park 2017). Their Poisson GLM model maximizes log-likelihood cost that incorporate the spiking multiple steps ahead in the future, thus addressing runaway self-excitation of neuron activity. However, their idea is specific to the Poisson GLM model and is not designed for non- stationary systems, thus limiting its applicability under varying applications.
Full forward and/or backward passes through e.g., Backpropagation Through Time (BPTT) as in Recurrent Neural Networks (RNNs), can partially handle long-term prediction instability. Nevertheless, approaches relying on this process may prone the vanishing or exploding gradient. More advanced RNN models that incorporate future steps integration to improve RNN forecasting include (Hess et al. 2023), that addresses the vanishing/ hwtpl//xxppnggggeee gradients through modification in RNN training to promote bounded gradients in training on chaotic systems and in piecewise linear RNNs Pal et al. (2021). on the other hand, proposed using particle flow to approximate the posterior distribution of future RNN states using a spatial graph of
element similaritiies Other RNN models, including Unni et al. (2023); | Xiao el al. (2023 Yeung et al. (2019)
leverage the Koopman operator for improved long-term prediction, yet, while powerful, they yield operators that cannot be directly interpreted as pairwise interactions between our elements of interest. Moreover, all
above RNN models are currently limited by their architecture and do not support handling more maanced dynamics (e.g., decomposed systems as in Mudrik et al. (2024a)), or effectively capturing smoothly-changing non-stationary behaviors. Another limitation of these methods is that they do not incorporate constraints on the operators (e.g., sparsity), which can be important for enhancing the operators' interpretability.
Hafner et al. 2019) proposed a model of reinforcement learning with future steps learning that learns the environment dynamics and chooses actions through online planning in a latent space. For the dynamics modeling, they use both stochastic and deterministic components based generalized variational objective that encourages multi-step predictions. However, their latent model is action- based and does not explicitly limit the dynamics to locally lincar functions, potentially complicating the understanding of evolution in the latent space. Additionally, they use gradient descent to learn the parameters, which may be prone to convergence local minima.
An additional approach to understanding dynamical systems involves identifying a sparse set of functions that jointly decompose the observations. For example, SINDy (Sparse Identification of Nonlinear Dynam- ics, Brunton et al. (2016) employs a data-driven approach to discover governing fundamental equations from data using sparse regression. Although SINDy and its extensions (e.g., Kaheman et al. (2020) are promising for discovering governing data equations, such representation does not provide explicit insight into the sime-changing interactions between the state elements.
5

---

Published in Transactions on Machine Learning Research (09/2024)
4 Specific models considered in this work:
Of particular interest in this work is improving the model fit of a core set of linear dynamical systems with different temporal constraints on the system cvolution, including 1) Inve-invariant Linear Dynamical
Systems (LDS)) 2) Switching Linear Dynamical Systems (SLDS) (Ackerson and Fu
1970:
Bar-Shalom and
Li 1990; Hamilton 1990; Ghahramani and Hinton [996] Murphy 1998 Fox et al.
2008
Linderman et al,
2017 3) decomposed Linear Dynamical Systems (dLDS, Mudrik et al. (2024a]; Chen ct al. 2024 : Mudrik
et al. 2024b)), and 4) regularized Linear Time- Varying (LTV) Dynamical Systems.
Time- Invariav Linear Dynamical Systems (LDS). In linear systems analysis, the evolution of a general state X over T ++1 1 time points can be typically represented as Lr+1 = Ac, + b, where A E Rexp is the time-invariant dynamics matrix and b E R remain constant over time. One common method for determining A (and if it is assumed that an unknown offset exists) involves a 1-step optimization approach that includes applying least squares across all time points. This entails solving
A, b arg min A.b
[(.TT) ( - [+]xTT @ bll?.
(3)
where X 1.r. and X [10:T-1] represents the noisy observations of the state from the second time point (t = 1) up to the last time point (T) and from the first time point (t = 0) up to T 1, respoctively. [1]xx @ b represents the horizontal concatenation of the column vector b, for T times. Here, A captures the average influence from - to 24, for all & I.
While such linear (time-invariant) systems are simple and therefore interpretable in terms of capturing referred to here "network" interpretability their overly simplistic nature often prevents them from adequately representing the complexities of real-world time scries.
element interactions- interactios-rreferree
Bar-Shalom and Li 1990; Hamilton 1990; Ghahramani and Hinton 1996. Murphy 1998] Fox et al. 2008 Linderman et al. 2017 along with other piece-wise stationary models (c.g., Song et al. (2021)) aim to provide interpretable representations of dynamics by identifying operators that govern periods of stationary behavior, with the system transitioning between these operators over time. Variations of SLDS include, e.g., recurrent SLDS (rSLDS), which introduces additional dependency between discrete switches and the previous state's location in space (Linderman et al.) 2017); ) and treo-structured recurrent SLDS, which extends rSLDS by incorporating a generalized stick-breaking procedure (Nassar el. al. 2018).
Switching Linear Dynamical Systemss (SLDS). SLDS models (Ackersn and Fu, 1970
While SLDS models usually involve transitioning from an observed to a latent low-dimensional space, here we chose focus on the case where switches occur within the observation space, essentially enforcing the transition to the latent space to be the identity operator. Tf we denote X E RNT as the noisy observations subjected to i.i.d Gaussian noise, SLDS models the evolution of T. using set of J discrete states (j f= I...),, where each state is associated with its own linear dynamical system fj. These discrete states switch between them abruptly at certain time points following an HMM model. During each inter- switch period, if the system is in the j-th discrete state, SLDS models the evolution of the state linearly as x, = f + bj. where fú represents the lincar transition matrix for the j-th discrete state and b; denotes a constant offset term for that discrete state. SLDS can be trained by an alternating set of steps between dynamic learning and the HMM update of the operators. This way, SLDS tackles the crucial task of capturing non-stationarities while preserving interpretability; however it inherently lacks the capability to distinguish between multiple co-occurring processes or overlapping subsystems. More information about the model assumptions, limitations, and parameter selection can be found in [Lnnderman et al. 2016 2017
dccomposed Lincar Dynamical Systems (dLDS). The Decomposed Lincar Dynamical Systems (dLDS, Mudrik et al. (2024a) model relaxed the time-invariant or piece-wise linear limitation of LDSs and SLDSs to support the discovery of co-occurring processes while maintaining interpretability. To em- phasize, here, for simplicity, we focus on the case where the dynamics evolution described directly in the observation space, while the full model presented in (Mudrik et al. 2024a) supports learning the dynamics within an identified latent state. Specifically, dLDS models the dynamics evolution (A,ã,- - w.) using a sparse timo-changing decomposition of lincar dynamical operators such that A. = (E)=1) resulting in (2)=1)93)) T -> Tt. These dynamical operators ((f))/z) are global, i.e., not time dependent, and hence
6

---

Published in Transactions on Machine Learning Research (09/2024)
are interpretable globally. However, their time-changing weights (ct) enable modeling non-stationary and non-linear dynamics (Fig. 2 right). Notably, dLDS is trained through an Eapectation-Maxiiizttion (EM) procedure where the global dynamic operators {fo and their time-changing coefficients {(s))=)))=
are updated iteratively to maximize the posteriors as:
(G arg max P(fc X (f,})
(4)
(exr)
(B)= arg max X,
(5)
(f)
Interostingly, dLDS can also capture linear or switching behaviors described carlier, by fixing the dLDS coefficients over time (for linear behavior, Fig. 2 far left) and by allowing abrupt change of coefficients
in specific time points (for switching behaviors, Fig. 2 middle left). More information about the model assumptions, limitations, and parameter selection can be found in (Mudrik et al. 2024a).
While dLDS presents several dynamic modeling advantages (e.g., captures non-stationarities, promotes in- terpretability), as it estimates the parameters for each time t solely based on the values of the precoding state at time t - 1, it does not address the issue of inaccurate long-term predictions" divergence.
Smooth or Sparse Linear Time- Varying Systems (LTV). In this paper, we refer to LTV systems that can be described by: 3t+1 = A, for all = 1..... We further assume that a regularization R(A,) may be applied the operators [A] This regularization can be inspired by the application (e.g., smoothness of operators time, || A,- || 62 or operator sparsity [|)ec(A,)]lo <1) and can mitigate the ill-posed
nature of finding A, separately for each time point.
5 LINOCS
In LINOCS we aim to learn the unknown dynamic operators (A,)== by integrating several multi-step predictions simultaneously into the inference procedure. This approach yields not only more accurate full-lookahead post-learning reconstruction but also operators that are more closely aligned with the ground truth. Particularly, for every t = 1... T, LINOCS finds the most likely estimate of (A,, by} given K+1 (K E Zzo hyper-parameter) multi-step reconstructions of orders O...K with different weights {wm)
Å. = arg min A.
wk 1 A
(6)
where Ti41 are the observations at time t ++1 and E is the k-order lookahead estimation, defined by the recursive rule for predicting the state at time 1, starting from the observations
Namely,
a A A, A
The weights {wwJ k associated with the orders 0...K K are dynamically adjusted throughout the infer- ence process (Fig. [3. This adjustment considers both the order mumber (k)) and the current reconstruction error related to that order, Eb (e.g., the & norm, eμ = ||, of ll)). Unlike other multi-step methods Venkatraman et al. (2015)), LINOCS adapts the weights of the reconstruction orders to prioritize the minimization of large errors in lower orders before considering higher orders (Fig. 3 Specifically, it grad- ually increases the weight of the best lookahead reconstruction until convergence conditions are satisfied. In our implementations, the weights can be chosen from a list of built-in choices such as uniform, linearly decreasing, and exponentially decreasing weights. Additionally, our framework allows custom weight func- tions that suit their specific needs. In the experiments presented in this paper, focus showcasing three specific options for the weights:
Adapting the weights to soquentially introduce higher multi-step reconstruction orders once the error for each preceding order falls below designated threshold, while continuing maintain the activation of lower orders. Specifically, in the initial iterations, only wi 0, with all other weights

---

being special case in Published in Transactions on Machine Learning Research (09/2024)
A: Problem
B: LINOCS (for K=3)
C: Expected Effect
D: Weight Function
x
11 F22
7 K=T
High single-stg ence mabty 20 sy Small single-sepp erce 01)
Donning licuienn
0.1-
Jennion IS
w
8
22:
2:
2:
22
Iteration
C-
e
955
toenntion
e-
1:
1::
1::
11:
0.0
4
k
80
Figure 1: Schematic of the problem and approach. A: Models that perform well 1-step predic- tion (i..., prediction order 0, in red), commonly fail in higher orders reconstruction (e.g., order-T, in green). B: LINOCS (e.g., for training order K 3) integrates weighted multi-step reconstructions for all k = 0... K orders. It adapts the weights of these reconstruction orders to prioritize minimizing large crrors at lower orders before addressing higher orders. The system gradually increases the weight of the most effective lookahead reconstruction until convergence conditions are met. C: LINOCS improves long-term reconstruc- tion during training iterations, such that the predictions during training are expected to increasingly align with the real system. Green color indicates full-lookahcad predictions across different iterations. D: Weights of different lookahead training orders (k). we R² + R is function of the order & and the k-th order reconstruction error e. Top: Illustration of an exemplary effect of on w when fixing k. Bottom: Illustration of an exemplary effect of on w when fixing e.
my 0 for E [], As the error of each step's reconstruction falls below this threshold, the subsequent weight Wj+1 is activated. For instance, if wo is the last activated weight and the error for the first (1-step) reconstruction falls below the threshold, w1 becomes active (w) > 0), and this sequence of activation continues for higher-order weights as each subsequent step achieves the required accuracy.
Constant weights over iterations with an exponential decay over k defined as wi = exp ob for some scalar Ry0-
exp-ok
weight function that considers both and c exhibiting a monotonic decrease in k and an increase e, with & decreasing faster than increases (Fig. DD).
Importantly, throughout this paper, we distinguish between two concepts: training order and prediction order. We denote "training order" (Kuain) as the maximum order considered during LINOCS training. In line with this, "K-step optimization" "K-step training" specifically refer to the use of the K-step cost (i.e. estimating [ during training, with "1-step optimization/traime which the training considers only consecutive time points for inference.
In contrast, prediction order refers to post-training predictions that involve iteratively propagating the iden- tified operators for Kpred steps into the future.
Here, we demonstrate the contribution of LINOCS for accurate long-scale predictions in four types of systems: 1) time-invariant linear; 2) switching linear; 3) decomposed linear; and 4) LTV systems. Importantly, in our experiments, assume that observe the underlying system under additive i.id Gaussian noisc conditions, however LINOCS can be easily adjusted to other noise stalistics.
5.1 LINOCS for Linear Dynamics
We first present the LINOCS learning rule for the simplest case time- invariant linear dynamical systems (Fig. leftmost subplot). Let X E RexT be the observations of state X, such that X X 7), with 7 being an i. Gaussian noise (n N(0,a2)) In this time-invariant linear case, X evolves linearly #++1 Ax. for all = 1... where b E Rex1 is constant offset.
8

---

Published in Transactions on Machine Learning Research (09/2024)
Linear
Switching
Sreudd-Switching
Sparse Decomposition
Coefficients
C2
C3
Time
Time
Time
Time
Figure 2: Example time-varying LDS systems A: The baseline linear time-invariant dynamical system
will present dynamical operator constant over time. B: Switching linear dynamical systems (SLDS) jump between different linear operators that are time-invariant between jumps. C: "Pseudo-Switching dynamics is similar to SLDS with the inclusion of smoother transitions between periods of constant linear dynamics. D: The decomposed linear dynamical systems (dLDS) model is generalization of SLDS to sparse time- changing linear combinations of linear operators. dLDS can also model negative coefficients.
In this case, LINOCS estimates A and b by
K Σ
k A 4bbl
A, = arg min A,b k=0
wv
T++1
A
A+
Tt
(8)
j=0
where A ++1 is taking the transition matrix to the power of k + 1, is an hyperpameter that dictates the maximum reconstruction order, and the set {wy] can be either predefined or automatically adapted over training basod on each order error. Please refer Algorithm and Appendix A.5 for further details the inference of the operator and offset for the linear case.
Algorithm Linear-LINOCS
Inputs: Observations X and maximun order K. Build K Multi-step reconstructions to get the set of {
described in Appendix D. via Eq. as described in Sec. A.5
A.5.
Infer Transition A
Infer Offset b
D via Eq.
16
5.2 LINOCS for Switching Linear Dynamical Systems (SLDS)
For SLDS (Fig. 2 middle-left), we integrate LINOCS into the SLDS operator inference stage (to infer
(f))=1, {b} see Soc. 4 using the SSM framework proposed by Linderman ct al. (2020) which is the current framework for rumning SLDS/1SLDS as described in e.g. Linderman et al. (2016] 2017 We maintain the existing SLDS approach to estimating switch times that delineate the boundaries of the linear periods al.)
between switches (i.e., we kept the switching times inference step as implemented by Linderman () To recover the operators within these identified linear periods, integrate the learning rule for the time-invariant linear case described above (Sec. 5.1 Please refer to Algorithm 2 for the procedural steps.
Algorithm 2 SLDS-LINOCS
Input: Observations X, maximum order K R+, number of itcrations. Initialize: (f))=1 Uniform
Initial paramoters
Start with uniform distribution
for each iteration until defined number of iterations do find switching time infer operalor in each period
& Based on current SSM framework Linderman et al. 2020
b Use linear inference as described in Alg
end for
9

---

Published in Transactions on Machine Learning Research (09/2024)
Roal
dbserved
LINCSS
stop
DAO
DAD reweigh
Rowdigh
A
Aquureed
Ausos
Anure
Paway
Avun
Oneo
Avelyy
0.98
4.2
0.55
2.98
-0.2
0.56
0.92
-0.2
0.53
D.7
-0.86
2.97
-0.
0.53
0.:97
+b.2
0:51
0.98
0.73
B
17.15
55
0.017
-30
-18.55
37.15
E
Tain Order
30
10
Prediction Order
LINOCS
Stop
Ai2
2i0 pereeggeiionnnnng Prediction Order (time)
Figure 3: Linear System Experiment. A: Real vs. identified operators and oflsets. B: Quiver plots of real and identified operators. C: Highlighted differences in effects between real operators and inferred operators showing how small differences in dynamic operators gain prominence during lookahead reconstruction (cal- culation details in Section A.II. D: Full lookahead reconstructions (ground truth vs. baselines) show swift convergence to the origin the 1-step optimization (yellow) and divergence for DAD-based results (three most right subplots). E: Frobenius norm of the differences between the ground truth operators (A) and the identified operators (A). F: MSE under increasing prediction orders. For all orders, LINOCS achieves better (lower) MSE compared to 1-step optimization. G: Full lookahead post-training predictions using operators identified by 1-step optimization (yellow) vs. the predictions using those identified by LINOCS (bluc) under various training orders (rows) and prediction orders (columns). H, I: LINOCS reconstruction compared to 1-step optimization under increasing noise levels (a) demonstrates its robustness. J: Propagating the identified operators until reaching relative reconstruction error of ~ 1. LINOCS identifics operators that can accurately predict ~ 35,000 time points (~ 1380 full rotations), much higher than 1-step training that decay immediately.
10

---

Published in Transactions on Machine Learning Research (09/2024)
5.3 LINOCS for decomposed Linear Dynamical Systems (dLDS)
For dLDS (Fig. 2 two rightmost subplots), in SLDS, we incorporate LINOCS into the dynamical systems update step. Let ãi+1 denote the k-th order reconstruction of #++1, calculated by iteratively propagating the dLDS reconstruction k + times, starting from æt-k- Furthermore, let I++1 Se=1eeffan. where represents the current estimate of the dLDS coefficients and (f, 1/=1 denotes the current estimate of the basis operators, We can now write the k-th order reconstruction (3/_1) as
Σ Ceff
Σ
1) TM.
Σ 2e ki Z.- - x
as follows with the recursive update rule
aff1: where T.. Tr.
To effectively integrate LINOCS into dLDS, we incorporate multi-step predictions into the training procedure of dLDS itself. Specifically, in each iteration, we start with the update the dynamics coefficients C₁. For this, we first define Fat E RPx as the horizontal concatenation
F = If a f2øt ................................ faat]] for some E O...K. K.
Next, we define new matrix FK that extends Fr, to all K + reconstructions stacked vertically, resulting in FK 2t e {(
Wi Sy F x
FK 2t
(10)
WK- F Fag
where we is the weight of the k-th multi-step order. This matrix can then be used to infer the coefficients (c:) while simultaneously considering different reconstruction orders with varying weights.
To mirror this concatenated matrix of dynamics with the observations (ã++1), we further define a matching concatenated state vector (#z+1)wert E [g((++)xx1, This vector (#.41) obtained by vertically stacking K + times the observations Ti+1 Rex1 at time t weighted by their corresponding wy values, re- sulting in (#++1) [wn W1 (1]. In simplified terms, (@t+1) == w & ē where w = [wo: wii- ; E [K++1)x1 are the training orders' currents weights, and @ denotes the Kronecker product.
The coefficients, et are thus updated every iteration by minimizing the squared l norm
& = arg min [I( pK cxlll xellcall1-
(11)
where de the weight of the &
regularization on the coefficients.
Note that the multiplication in the first. term (F* e E (((++)xx1) produces vector of estimates of the observations at + 1 (x4+) computed from all different K + past states. This way, the estimator in Equation 11 secks the vector that best predicts #++1 considering all K ++1 lookahcads.
The next step within every iteration includes updating the dynamic operators (fj)- One addit ional modification we make (compared to the original learning of dLDS as presented by Mudrik et al. (2024a) is that rather than updating cach f, using gradient descent, we infer the dLDS' basis dynamics operators (5al-n by fully and directly minimizing the cost. Specifically, let.
Fau [fi, f2, fa,...,,.
.2,...,,f.. fr] DPxpJ
11

---

Published in Transactions on Machine Learning Research (09/2024)
be the concatenated matrix of all (fol)=1-
Also, let (x.), = ([]) ( (cr @1]xx)) Rexx1 where denoted the Kronecker product and denotes
element-wise multiplication, and let. X. ReJxT be the horizontal concatenation of all {(ar))))=0) )
With these definitions, the dLDS operators (fi))=1 are updated by
Faul arg min X Faul (X Fall
(12)
with each {fyly= then being extracted from Fall and normalized to Frobenius norm of 1. Please see Algorithm 3 for the procedural steps.
Algorithm 3 dLDS-LINOCS
Input: Observations X, maximum order K R+ Initialize: Ce R/ and (f) Uniform f) while not converged do Compute & lookaheads of ž using current estimates of {aH and (f)
> Normalizing each f, to unit spectral norm
D via Eq.
9
Construct FK #s
via Eq.
10
Update e using the LASSO method.
> via Eq.
Update Fall.
via Eq.
Split Fall RPxpJ into {f,}]=1 (cach f, E ROXP).
end while
5.4 LINOCS for regularized Linear Time Hone-Varying Systems (LTV)
Finally, we focus on the more general case of regularized linear time varying systems that are not necessarily switching or decomposed.
particular, we focus on two types of regularizations, 1) the operators change smoothly over time, i.., ||A, - A.- small, and 2) the operators are sparse, i.c., ||A is small.
For the LTV case, we apply LINOCS to find the time-changing operators {A.)= by iteratively integrating multi-step reconstruction with the appropriate regularization. The operators are initialized with a regularized 1-step optimization
Å, arg min p A R(A.) = arg min A. |.
A, + AlIA; A,-1
(13)
where À is the regularization weight. While in this paper we chose to focus on R(4.) = AI|A, 1 other regularization terms can be used in a more general sense.
We integrate LINOCS into the estimation process by iteratively updating the operator estimates one at a time. Specifically, during each round of updates, we loop over every time point 0... T and hold all operators of times + ((A)) fixed their former estimates. We then update A, by
K
k+1
A, arg min A:
wk
Σ
... A... A A A
+ R(AA)
k=0
1
where & = 0...A K denotes the order of the reconstruction and 1 k; denotes the starting point of the reconstruction. The weights " are set as in the previous models. Please see Algorithm 4 for the procedural steps.
6 Results
To showcase LINOCS' ability capture the dynamics in multiple models, we applied LINOCS the above systems under diverse settings. The hyper-parameter used in each experiment are summarized in Section A.3
12

---

Published in Transactions on Machine Learning Research (09/2024)
Algorithm 4 LTV-LINOCS
Input: Observations X, maximum order K Rt,, number of iterations M. Initialize: (At)== Uniform
Initialize from uniform distribution
for each iteration = 1...M do for each time point -= 1 do
To infer A,, consider K + windows with shifts up to K
Define the following quantities:
x = diag( -k E RoK**.*
1141 18+1 1K
diag( ( | k==1k==0)
k+1
TROK" xk* xK"
B == concat((1]]] k1
j= Ar
12 k;=1kk0p axis 1) E RexpK"
k+1 1K
b horizontal concatenation of A
estimates for windows
C = :==oncat(([[I] j-1
- - k+1 18 axis = 0) E RRP xp
vertical concatenation of A
estimates for shifts
V = Bæ E ROOK
A. arg minA, A, 1| CAN A-1
Solve for A, through Ridge regression
optimization
end for end for
6.1 LINOCS more accurately identifies ground truth linear systems under noisy observations
We first. test LINOCS' ability to robustly learn time-invariant linear dynamical systems from noisy obser- vations. We then simulate the dynamics E R2x2 as a rotational transition operator and a random offsct b E R2x1, where each b, (niform(O, 1). We build the synthetic state æ, R2x1 for T = 500 time points, V, Ax + b starting from random initial value zo E Uniform(0, 1)2x1, such that the noisy observations are H x + 7). where the noise (() = (O, 0.3²) (Fig. 3AB)..
We compare the learned operators using LINOCS against four baselines. First we compare to traditional 1-step optimization (Eqn. (3)) We further compare the linear LINOCS to our implementations of the DAD approach (Venkatraman al. 2015), as it the approach closest LINOCS terms of integrating multi-step predictions into model training.
Our implementation of DAD integrates expert and non-expert demonstrations for model training, inspired by the Dataset- Aggregattnn (DAgger) approach (Ross et al. 2011). Specifically test three implementations of DAD. For each implementation we initialized the transition matrix (Ainit) and the offset (bima) using the optimal estimate from 1-step optimization. We then, in each DAD implementation, train the model through 100 iterations where at each iteration, we used our last estimates of A and b to perform full lookahead reconstruction, starting from time O. We then update our estimates of A and using the optimal 1-step optimization while considering both the observations and the above full lookahead reconstruction.
Particularly, for the comparisons to DAD, we tested all three options outlined below:
DAD with full model update:
At each iteration, tion of the state (w.) calculated based the last operators estimate. Namely, (Aner+I, ( = arg arggminA.b + ET-1 || - (Aâ. + b)l..
update A and b based on the lookahead reconstruc-
DAgger-inspired DAD (reweighed DAD):
For reweighted DAD, we estimate and at each iteration using both the observations and the lookahead reconstruction from the last estimates of A and b. In particular, let. [R,; क.. E Rox2 be horizontal concatenation of the lookahead reconstruction and of the observations time t. Then iteratively solve: {Aiter41, b
Diter+ }
arg min (A.b) + - [. E..1 (A | Ø. b)ll
DAgger with lz constraint (reweighted DAD with lg):
This option solved similarly to the reweighted DAD, with the addition of the Frobenius norm (II []r) on the operators (A) and (b) during training.
13

---

Published in Transactions on Machine Learning Research (09/2024)
We find that LINOCS identifies operators that yield accurate dynamics in long-time scale predictions (Fig. 3D). The other methods we tested, including 1-step optimization (Fig. 3C,E, cyan) and DAD-based implementations (Fig. BD, reds), instead indeed decay to zero (away from the real system), indicating a less accurate ostimation the dynamics. The improved accuracy of the operators identified by LINOCS (Fig. 3E) becomes apparent when examining the effects of the operators' estimation errors (Fig. 3C). These errors are larger for the other methods, showcasing that those methods accumulate more errors over shorter time spans than LINOCS does.
We next investigated the effect of the training order in LINOCS on long-term reconstruction. We trained LINOCS on the noisy observations with incrcasing training orders and then tested the performance un- der increasing prediction orders (Fig. 3G). Compared to the baselines tested, LINOCS exhibits increased performance even with very low training orders (e.g. 5), with higher orders resulting in almost perfect roconstruction (Fig. 3G bottom-right. subplot).
Additionally, exploring LINOCS's robustness to noise reveals that, unlike one-step reconstruction, LINOCS is robust even under very high levels of noise (Fig, 3G, H, blue). The resulting MSE compared to the ground truth dynamics is much lower in LINOCS, even under very high noise levels (Fig. 31,, 1 blue vs. orange-red).
When examining the duration for which LINOCS remains robust without converging, we observe that our approach accurately predicts approximately 35,000 time points into the future before deviating from the real system and decaying to 0 demonstraaiin stability over exceptionally long time scales (Fig. BT).
We further tested LINOCS on linear systems with structured noise (Fig. 13 well on simulation of 3-dimensional cylinder (Fig. 14p, yielding similar results. For structured noise, we modeled the observation ã. 2, + a syt) with o = 0.5 and 3 for = 1 501. Unlike other methods, LINOCS found operators that led to accurate long-term predictions (Fig. 13 D). Moreover, when examined under increasing training and prediction orders, we found that LINOCS is robust for long-term prodictions, even for full lookahead reconstructtions (Nipred = 501, Fig. 13 E,H). When evaluating its robustness to increasing structured noise levels (o)), we found that even for very high noise levels (a = 0.9), LINOCS achieved much more robust results than 1-step optimization (Fig. 13 F, F,G). Additionally, when exploring how far into the future it enables robust reconstruction before converging, we found that it is capable of full lookahead for approximately 70,000 time points a testament. to its ability to find more robust operators that can adequately describe the system (Fig. 13 1).
For the 3D cylinder case (Fig. 14 with Gaussian noise (a= 0.4), we similarly demonstrate that LINOCS recovers more accurate operators, leading to significantly more robust long-term predictions and enabling full recovery of the process (Fig. 14 C,D,C, both under increasing prediction orders (Fig. 14 E, G) and noise (a) levels (Fig. 14 H). It further exhibits an impressive ability to reconstruct lookahead predictions (starting from 2o) for very long periods (approximately 70,000 time points) before converging to similar error of 1-step (Fig. 14 Ij. In contrast, 1-step optimization yields high-error within few prediction orders.
6.2 LINOCS identifies accurate interactions in switching systems
We next tested LLNOCS-driven SLDS as detailed in Section 5.2 on simulated data comprising of J = 3 discrete states. The transition operators for each of the distinct states was set to a 3 x 3 rotational matrix oriented in a different direction. Additionally, the offset for each state (b, E R3x1) was set to be the same random vector drawn from a uniform distribution between 0 and 1 (Fig. 15D).
Notably, since the method is invariant to the order of the operators, to compare the identified operators to the ground truth operators, sorted the operators using the "linear sum assignment" problem (SciPy's implementation, by Crouse (2016)), with the cost function being the Frobenius norm between each pair of fs (ground truth vs. estimated for each model). As baselines, we compare the results of LINOCSSaagmenteed SLDS with standard SLDS and recurrent SL.DS (Linderman et al. 2016) with varying numbers of iterations.
When comparing LINOCS-SLDS to the basclines (Fig. 4 LINOCS consistently outperformed the other approaches across multiple metrics inchuding operator recovery (Fig. [C,E), switching times recovery (Fig. (A,)), and dynamics reconstruction (Fig. 43, Fig. 21 In particular, L.NOCS-SLDD accurately
14

---

Published in Transactions on Machine Learning Research (09/2024)
UNOCSE
EMOCES
LNDOO
5 4:22
0317
CONT
COO
LOUT
3-5
4.07
333
340
.COM
421
400
428
327
2000
LOOD
2003
A00A
C007
LCM
2--
4.90
-O0MSS
4000
<K
32
4766
Todt
8-
XXI
411
0M0
6.00
3000
612
2331
100
026
180
230
220033
A4600
14000
LC00)
1c00
12w00
40003
80
30
60
Figure Results on switching systems. A: Active discrete states for LINOCS (blue) compared baselines, including SLDS and rSLDS with 10 100 training epochs. B: Correlation (p) between the ground truth dynamics (a) and the full-lookahead reconstructed dynamics (ø). C: Correlation (p) between the ground truth operators (ff) and the identified operators (f). D: Correlation (p) between the ground truth coefficients (c) and the identified coefficients (c). E: Difference between the ground-truth sub-dynamics (F) and reconstructed basis dynamics by different models. LINOCS was able to achieve sub-dynamics that are much closer to the ground truth than the other basclines.
identified switching times, whereas classical SLDS and rSLDS tended to introduce additional redundant switches (Fig. AA). Moreover, the discrepancies between the ground truth operators and those identified by LINOCS (Fig. 4 right-most four columns) were substantially smaller than the differences observed with classical SLDS SLDS/ISLDS ISLDS (Fig, | E, left columns), as evidenced by the higher correlations between LINOCS operators and the ground truth (Fig. 4C). Furthermore, when examining the cigenvalues of the identified operators compared to the ground truth (Fig. 16), the eigenspectrum derived from the LINOCS-driven solver closely resembled the ground truth eigenspectrum more than the classical SLDS and rSLDS cases, highlighting the effectiveness of LINOCS in capturing the underlying dynamics.
6.3 LINOCS finds dLDS operators that yield accurate dLDS lookahead predictions
Next, applied LINOCS to dLDS, as described in Section 5.3 First we generated ground-truth data that represent a "pseudo-switching (Fig. 2 process process-i.e. linear dynamics that switch more smoothly (in our case between J = 3 systems) compared to SLDS where operators switch abruptly. This creates overlap periods where two dynamical systems are active at once as they trade off (Fig. 17). LINOCS-dLDS demonstratod significantly improved slability full lookahead reconstruction compared to single-step dLDS (Fig. 5 No- tably, training with orders approximately greater than 35 (Kirain > 35) on our synthetic dataset (containing 1000 time points) resulted in highly accurate full reconstruction (Fig, 5 Additionally, when comparing MSF and correlation of the time-evolving operator F, = E Caff the ground truth, we observed a monotonic decrease in MSE with increasing maximal LINOCS training orders (Krrai), while the correlation showed monotonic increase (Fig. 5 3).
Interestingly, although the 1-step prediction (post-training) is seemingly good also for non-LINOCS dLDS (or low-order LINOCS-ALDS) (Fig. P left), the advantage of LINOCS is revealed in Figure 5D right, under the full lookahead reconstruction. This implies that evaluating dynamical models not only based on their ability to predict immediate steps but also on their performance in further steps (i.e., under multistep roconstructions) is critical for more nuanced evaluations, as reconstruction errors could be obscured in 1- step predictions. Additionally, this comparison also highlights the importance of integrating multiple orders simultaneously during training.
15

---

Published in Transactions on Machine Learning Research (09/2024)
A
Ground Truth
Avas
Peay
Rosa
Keay 15
Pyuk 23
Roou 36
1-0000 prediction
B5
R5
05
4.5
9.
05
6P55
9.
60
65
b.
6.00
Lockateod predition
25
50
6.0 6
-10
-5
50
0
00
0:30
30
50
B
C
D
1-step (MSE)
full looksheed (MSE)
1Ix 2012
Time
Time
Figure Decomposed linear dynamical systems results. A: Ground truth dynamics compared to 1-step (top) and full lookahead (bottom) reconstructions for non-LINOCS dLDS (Kirain 1) and LINOCS- dLDS with different training orders (Kerain E [2,5, 7, 15, 23, 36)) B: MSE (pink) and correlation (green) between ground truth operators and the operators identified by LINOCS under different orders. C: MSE (pink) and correlation (green) between ground truth dynamics and full lookahead reconstructions using the different LINOCS training orders. D: Local MSE for 1-step (left) and for full lookahead reconstruction (right) over the time points of the dynamics.
We further extended our study encompass more nuanced dLDS settings, exhibiting prolonged time scales and recurring patterns of identical active operators across distinct intervals (Fig. 18 We found analogous enhancements of LINOCS over the traditional 1-step dLDS implementation. Specifically, LINOCS demon- strates robust accurate long-term predictions, including full lookahcad prediction (Fig. 6B), in contrast to 1-step optimization, which yield high lookahead error (Fig. 6C, last three subplots). Furthermore, also for - cafs
this more complex example, upon comparing the identified time-varying transition operators F; = to the ground truth, LINOCS revealed operators with cigenvalues significantly more correlated with the real operalors'' evaluations (Fig. GD, E) compared to the 1-step optimization results. Additionally, when comparing the operators' values themselves against the ground truth, those identified by LLNOCS exhibited higher correlation and smaller MSE with the ground truth compared to these identified by 1-step dLDS (Fig. F.G.H).
6.4 LINOCS finds interactions that yield robust lookahead predictions in Linear Time-Varying (LTV) systems
To test the applicability of LINOCS to more general LTV systems, implemented LINOCS-LTV to capture the chaotic behavior of the Lorenz attractor (Sec. A.4b through a smoothly changing L.V approximation (Fig. 7 We compared LINOCS-LTV with several other LTV solvers with varying constraints, including sness and sparsity (r = 6, 7,8 and smoothness with weights A = 2, 20, refer to Sec. A.2 for details). Unlike methods relying on 1-step optimization, LINOCS, despite similar regularization constraints, achieved superior full lookahead reconstruction (Fig. bottom).
Also here, while different methods performed satisfactorily in the 1-step (post-training) prediction (Fig. 7 top, B red, C red), disparities emerged in higher-orders lookahead predictions where alternative meth- ods failed. While all methods, including LINOCS, achieved commendable 1-step reconstruction, LINOCS demonstrated a markedly lower full lookahead error (Fig. 7 green, 5 most right bar pairs) and superior full-lookahadd reconstruction correlation with the ground truth (Fig. 7 green, five most right bar pairs).
In addition, we analyzed operators identified across various training iterations of LINOCS to assess their proficiency in achieving lookahead reconstruction (Fig. 19 For this analysis, we used the Lorenz attractor with 900 time points with intervals of 0.1/9 arbitrary units (a.u.), and applied a smoothness constraint with
16

---

Published in Transactions on Machine Learning Research (09/2024)
A
B
Ground Truth
Kyrev
Kyedd
Kavod 10
Kyrca
Kped
Kyrdd 2000
10
1Y
C4
13
-30
04
11
04
6
-1
-29
-10
C
REAL
LINOCS 1 STEP
14
CY
-2
G
H
10
c0m6)) (mxA))
p(Ra(A), RolA))
40-
0:1-
18 u5e a00
Figure 6: Additional decomposed linear dynamical systems results (Kursin = 50). A: Ground truth dynamics compared to (B) LINOCS results and (C) 1-step optimization results, for increasing prediction orders. D,E: Correlation between the eigenvalues of the ground truth transition matrix (F. = Ei carf;) and the eigenvalues of the one identified by LINOCS. D: imaginary part; E: real part. Results display the average correlation over time. Eigenvalues were matched using the "lincar sum assignment problem" (Scipy's Crouse (2016) F,G: Comparing the identified time-changing transition matrix (Fi) identified by LINOCS vs. 1-step optimization in terms of correlation (F) and MSE (G). H: Comparing the MSE of the identified F; over time.
weight of == 0.1. We observed that over training iterations, LINOCS adaptively influenced the predicted lookahcad dynamics to gradually converge towards the ground truth dynamics (Fig. 19 C), with monotonic decrease MSE (Fig. 194, B).
When analyzing which time points of the dynamics contributed to higher MSEs in the full post-training lookahead prediction, we noticed that carly training itcrations tended to produce higher full lookahead prediction errors later time points of the dynamics (Fig. 19 top right). However, over subsequent iterations, the effect of LINOCS managed to mitigate the accumulation of errors these late time points bottom right).
(Fig. 19 193,
6.5 LINOCS finds robust interactions in real-world neural data
Finally, we applied LINOCS to real-world dataset described by Kyzar et al. (2024), which consists of high density electrode array of populations of single units in the human medial temporal and medial frontal lobes while subjects were engaged a screening task. We applied linear LINOCS, SLDS, dI.DS-LINOCS, and LTV-LINOCS to a single recording session that includes recordings from five brain arcas (amygdala left and right, cingulate cortex, hippocampus, pre-supplementary motor area). All the dynamical systems models were trained on the firing rate data, which we inferred from the spike-sorted electrophysiology via Gaussian kernel convolution (sce Sec. C.6 for details).
We investigated several LINOCS models to showcase their distinct characteristies. First, we examined the linear case for each brain area individually and explored the mean field interactions between areas (Fig. S Importantly, while typical real-world brain dynamics are assumed our aim in starting with the linear model was to demonstrate how LINOCS can identify the fundamental background neural interactions under linear assumptions and check how its identified interactions defer from
be non-linear and non-stationary,
17

---

Published in Transactions on Machine Learning Research (09/2024)
LNNCS
LINOCS
Full lockahead prodiction
Stop
Soout
LIMCS
C
Ras
14-300
60
Figure 7: LTV approximation of the Lorenz attractor. A: 1-step post-training prediction (top, pink- red) vs. full lookahead prediction (bottom, green) for different baselines, including LINOCS with various smoothness levels or to regularization. Color indicates the local MSE. refers the regularization weight. B: Squared lg of the error between the ground truth and 1-step prediction (red) vs. full lookahead predictions (starry green), for the different. methods. C: Correlations between the ground truth and 1-step prediction (red) vs. full lookahead predictions (starry green), for the different methods.
these identified by the 1-step approach. We first applied the linear LINOCS on the firing rale activity from all neurons within each region to identify between-region interactions. We observed that LINOCS identified different linear interactions within areas compared to the 1-step optimization. Drawing from our conclusions based on synthetic data linear results, this suggests that LINOCS may provide more muanced linear approximation of brain activity compared to the common 1-step optimization (Fig. S A), though further comparisons with established neuroscience methods are needed to fully validate this advantage.
Then, we also applied the linear LINOCS on the mean activity of each region, and found that when examining the full lookahead reconstructions (Fig. g the 1-step optimization, in contrast to LINOCS, decayed to zero activity due to small accumulated deviations in operator values. In contrast, LINOCS managed to maintain activity closer to the average values of the dynamics. However, due to linear enforcement, neither approach could capture fluctuations in dynamics. Moreover, the full lookahead reconstruction error for LLNOCS-lineas was overall much smaller compared to the classical 1-step (Fig. 3)
We next applied LINOCSSSLDS with three discrete states and compared it with regular SLDS using the same number of iterations. LINOCS identified operators that exhibit slight differences compared to those found by classical SLDS (Fig. 9
B vs Fig. 22 Fig. 23 A vs. B) well as slightly different switching patterns between the two approaches Figg. 23 C,D). These operators resulted in significantly more robust lookahead predictions. Specifically, differences are evident in both connection presence, weights, and distribution among global operators. For example, in the " "Amygdala left" region, both classical SLDS and LLNOCS-driven SLDS identify a connection from neuron 10 to neuron as part of fa, albeit with varying weights. Additionally, both methods identify connections from 5 to 3 (in f1 for classical SL.DS and in f2 for LNNOCS-SLDDS) as well as from to 3 (in fa for classical SLDS and in fi for LINOCS-SLDS), but with differing weights. Similar discrepancies are observed in other regions. Furthermore, LINOCS-SLDS and classical-SLDS each identify connections that the other overlooks; for instance, in the "Amygdala right" region, LINOCS-SLDS identifies strong connection from 11 to 6, whereas classical SLDS does not. Conversely, classical SLDS identifies a connection from 3 to 6 (in fa), which LINOCS-SLDS does not recognize.
Importantly, the operators found by LINOCS enable full lookahead reconstruction without diverging, in contrast to regular classical SLDS that diverge to extreme values in full lookahead prediction (Fig. 9C). Moreover, the reconstruction error for the full lookahead prediction was overall much smaller for LINOCS- SLDS compared to the classical SLDS (Fig. [A). These observations suggest that if the real neural process
18

---

Published in Transactions on Machine Learning Research (09/2024)
A
Haporaa Pro supiometir moter 3r05 left
1 step
LINOCS
Datiaal on aatittor cingolaie catioa left
Antrgeoia right
Any let
B
U50C5
Ground Truth
stop
dameal cingurate cofTeN
phe supphementary orea
amygdale
amygdala
Nppocampus
Figure 8: Application of linear-LINOCS to multi-region neural recordings. A: The per-region and between-region mcan field lincar dynamics operators identified by 1-step lincar optimization vs. LINOCS with a linear lime-invariant system. Each within-region network describes the linear operator Aregion derived by applying LINOCS the firing rate matrix constrained to include only the neurons from that region. Fach node represents a unit of the spike sorted data, ideally corresponding to single neuron. Colors of nodes represent their estimated brain region. Interactions between areas are found by applying LINOCS to the mean activity per area. Edges width indicate connection ion strength. B: Full-lokkahead reconstruction error for 1-step linear optimization vs. LINOCS- approach. C: Ground truth mean activity per region compared to the lookahead prediction trace of the mean field activity by 1-step linear optimization vs. LINOCS-lined
follows switching dynamics, LINOCS may capture the underlying dynamics more effectively, as inferred from our analysis of the synthetic case.
We observed similar patterns using ALDS-LINOCS, which revealed underlying global brain interactions potentially fundamental to brain function (Fig. 10 A). When examining their dynamic activations (c:), we noted "background" interaction consistently active, with slight modulations over time (Fig. 10 B, brown), alongside gradually changing activities of other interactions (Fig, 10 B, gray-blue-purple) Importantly, these results provided lookahead predictions that did not decay and maintained a high correlation with the observations (Fig. 10 C).
Finally, employed the LTV-LINOCS on all neurons from all regions simultaneously while imposing a smooth- constraint consecutive operators (with regularization of X 0.1 on ||A, 11. Fig. [.p. Our findings reveal that LINOCS identifies operators capable of producing full lookahead reconstructions with- out divergence, closcly approximating observed data. Comparative analysis against 1-step optimization with various smoothhess levels (Fig. CI,,C, underscores LINOCS' ability to achieve better reconstructions than other approaches. Additionally, examination of error evolution over time suggests a monotonic increase in error for non-LINOCS approaches (Fig. (1C). Morcover, we obscrved notable discrepancies between the operators identified by by LINOCS and the baselines (Fig. ITD). These results highlight the efficacy of LTV-LINOCS in capturing complex temporal dynamics in real world data while maintaining data fidelity.
19

---

Published in Transactions on Machine Learning Research (09/2024)
A
Original SLDS
LINOSS SLDS
B
Attygdbe
Hippocapi
10a
10
||x 211
10
10
10
107
10a
109
10-3
Dorfa Aommon
Hppocampus kitt
metior
3r03 kift
Orggin
3300
20-
1300
12
<16 -17-
300-
300
Timo,
TE090
Time,
Figure 9: SLDS results on real data. A: Relative error the full lookahcad prediction compared to the ground truth, for both LINOCS-SLDS vs. classical SLDS,for each brain area. B: The networks identified by LINOCS-SLDS (see Fig. 22 for the networks identified by the classical-ST.DS). C: Lookahead reconstruction using LINOCS-SLUD (solid curve) vs. classical SLDS (dashed curve). Classical-SLDS diverge to extreme values in the full lookahead reconstruction.
A
A
f
Amygdala
Hippocampus left
C1
C3
C5
C7
Amygdala right
Pre supplementary motor area left
C2
C4
C6
C8
Dorsal anterior cingulate cortex left
B
D
225
06-
610-
02
100ms
Time x
100ms
Time x 100ms
Figure 10: dLDS-LINOCS results on the real neural data. A: The identified operators { B: The identified sparse coefficients ct. C: Correlation between the hull-lookahead reconstruction results and the observations. D: The two largest eigenvalues (A) of the time-varying transition operator F, ΣS= cjffj. Black (top): real part. Red (bottom): imaginary part.
Overall, we showed that in all these real-world neural versions, LINOCS was able to recover more robust descriptions of the dynamic evolution for the long run, which, based on our synthetic results, may imply that these closer to the real unknown interactions.
20

---

Published in Transactions on Machine Learning Research (09/2024)
A
B
C
D
A=0.1
12ah
A=02
UNOCS
47-
64
-2
371
15"
E
-07
2900
1.0
30
70-
23-
2000
1000
600
-04
02
-01
-02
70
10-
F
1.0
-2002
1.000
1.75
-25
15
10
1200
3030
13
-10
04
04
004
35
200
02
200
000
Time,
Time,
lime,
Figure 11: LTV-LINOCS results real neural data. A: The ground truth data. B: MSE between the ground truth data and the full lookahead reconstructions. C: MSE of the reconstruction over time points, compared to 1-step optimization with different smoothness levels. D: Frobenius norm of the differ- ences between As identified by the different models, normalized by the magnitude of the operator identified by LINOCS. B: Full-lookahead predictions produced by LINOCS-LTV (left) and 1-step optimization with increasing smoothness constraints (A). F: Difference between the observations and the full lookahead recon- struction for LINOCS-LTV vs. 1-step optimization with increasing smoothness levels.
7 Discussion
In this paper, we introduced LINOCS (Lookahead Inference of Networked Operators for Continuous Stabil- ity), a learning procedure to improve stability and accuracy of dynamical system inference that leverages multiple lookahead estimations. By iteratively integrating re-weighted multi-step reconstructions with addi- tional constraints on the operators, LINOCS enables robust inference of worked operators in dynamical systems, even in the presence of noise and nonlinearity.
Our experimental results highlight LINOCS offectiveness across various dynamical systems, including Linear Systems (LDSs), Switching Linear Systems (SLDS), decomposed LDSs (dLDSS (Mudrik al., 2024a) ; and Lincar Time- Varying Systems (LTV) in both simulated and real-world neural data. LINOCS achieves more precise full lookahead reconstruction and more accurately retrieves ground truth operators in synthetic data compared to baseline methods, highlighting its ability to better capture muanes in the underlying system. These findings suggest that LINOCS holds greater potential than alternative approaches for accurately identifying unknown hidden interactions also in real-world data, where the real underlying interactions are typically obscured but pivotal for robust. scientific interpretation.
Looking ahcad, several promising avenues exist for future work. These include applying LINOCS a broader range of datasets uncover new scientific discoveries about component interactions (e.g.- using LINOCS. dlds to identify functional interaction motifs in C. elegans, done for 1-step inference in Yezerets et. al.
21

---

Published in Transactions on Machine Learning Research (09/2024)
2024))- Additionally, from computational standpoint, LINOCS could be extended to improve operator inference and reconstruction robustness in models like TVART Harris et al. (2021), or be integrated into RNN training to enhance learning and prediction stability. Particularly, if extending LINOCS to deep networks, the integration of multi-step reconstructions into the networks' training, may be used to address issucs such as vanishing or exploding gradients. Another future step would be to integrate an offset term identification into the dLDS application (e.g., as addressed in Chen al. (2024) for 1-step inference). Additionally, extending LINOCS handle non-linear local transformations (i.e., 22 -(A,3t=1) i), where g is a non-linear activation), well as non- non-Gaussian noise, could enhance its applicability to wider range of real-world scenarios.
8 Code and Data Availability
The code is shared on GitHub and available at http:://ggtttb..co////ggaatttait notebooks available here (for the Linear-LINOCS example) and at https://githut. http::///ggccuu....oooo/ogggaaaa com/NogaMudrik/ LINocs/bloo/mmain/ru aLDS Lorenn_ caammpee.pynnt for the ALDS-LINOCS example. The human neural
with tutorial
recordings data we used to excmplify LINOCS available at (Kyzar et al. 2024
A video briefly explaining the paper is available
9 Acknowledgments
We thank our funding sources for supporting us during this project. N.M. was funded by The Kavli Founda- tion NeuroData Discovery award. Y.C. and C.R. were funded by James McDonnell Foundation grant num- ber 22002039, with Y.C. being further funded by National Institutes of Health grant number 2T32EB025816, and C.R. being further funded by the Julian T. Hightower Chair. A.S.C. and E.Y. were partially supported Award 2340338 and a Johns Hopkins Bridge Grant.
by the NSF
22

---

Published in Transactions on Machine Learning Research (09/2024)
References
Guy A. Ackerson and K. S. Fu. On State Estimation in Switching Environments. IEEE Transactions Automatic Control, AC-15(1):10-17, 1970. ISSN 15582523. doi: N0.1109/TAC.1970.100
Maira Aguiar, Carlos Braumann, Bob Kooi, Andrea Puglicse, Nico Stollenwerk, Ezio Venturino, et al.
Current Trends in Dynomical Systems in Biology and Natural Sciences:(Prefacc)) Springer, 2020.
Travis Askham and J Nathan Kutz. Variable projection methods for an optimized dynamic mode decompo- sition. STAM Journal on Applied Dynamical Systems, 17(1):380-44 -166 2018
Yaakov Bar-Shalom and Xiao-Rong Li. Estimation and tracking. Artech House, Boston, MA, 1990.
Sumanta Basu, Xianqi Li, and George Michailidis. Low rank and structured modeling of high nigh-dimensional vector autoregressions. IEEE Transactions on Signal Processing, 67(5):1207-1222,, 2019.
Jeff A Bilmes. What hmms can do. IEICE TRANSACTIONS Information and Systems, 89(3):869 891, 2006.
Steven Brunton, Joshua Proctor, and J Nathan Kutz. Discovering governing equations from data by sparse identification of nonlinear dynamical systems. Proceedings of the national academyy of sciences, 113 (15):3932 3937, 2016.
Yenho Chen, Noga Mudrik, Kyle A Johnsen, Sankaraleengan Alagapan, Adam S Charles, and Christopher J Rozell. Probabilistic decomposed linear dynamical systems for robust discovery of latent neural dynamics. arXiv preprint 2rv: 16862, 2024.
Guillaume Chevillon. Direct multi-step estimation and forccasting. Journal of Economic Surveys, 21 (4): 746 785, 2007.
David R Cox. Prediction by exponentially weighted moving averages and related methods. Journal of the Royal Statislical Sociely: Series B (Methodologicaa), 23(2):414-422 422, 1961.
David F Crouse. On implementing 2d rectangular assignment algorithms. IEEE Transaclions Aerospace
and Electronic Systems, 52(4):1679-166) 1696, 2016.
Raina D'Aleo, Adam Rouse, Mare Schieber, and Sridevi Sarma. Quantifying interactions between neul ral populations during behavior using dynamical systems models. In 2019 41st Annual International Conference of the IEEE Engineering in Medicine and Biology Society (EMBC), pages 1960-1964. IEEE, 2019.
John Ferré, Ariel Rokem, Elizabeth Buffalo, Nathan Kutz, and Adrienne Fairhall. Non-stationary dynamic mode decomposition. IEEE Access, 2023.
David F Findley. On the use of multiple models for multi-period forecasting. In Proceedings of Business and Economic Statistics, American Statistical Association, pages 528 531, 1983.
David F Findley. Model selection for multi-step-ahcaa forccasting. IFAC Proceedings Volumes, 18(5):1039- 1044, 1985.
Blaettler Florian, Kollmorgen Sepp, Herbst Joshua, and Hahuloser Richard. Hidden markov models in the neurosciences. Hidden Markov Models, Theory and Applications, page 169, 2011.
Emily Fox, Erik Sudderth, Michael Jordan, and Alan Willsky. Nonparametric bayesian learning of switching linear dynamical systems. In D. Koller, D. Schuurmans, Y. Bengio, and L. Bottou, editors, Advances in Newral Information Processing Systems, volume 21. Curran Associates, Inc., 2008. URL https://
Paper www.t/9--appppp... .app
Zoubin Ghahramani and Geoffrey E. Hinton. Switching state-space models. Technical report, King's College Road, Toronto M5S 3H5, 1996.
23

---

Published in Transactions on Machine Learning Research (09/2024)
Danijar Hafner, Timothy Lillicrap, Ian Fischer, Ruben Villegas, David Ha, Honglak Lee, and James Davidson. Learning latent dynamics for planning from pixels. In International conference on machine learning, pages 2555 2563. PMLR, 2019.
James D. Hamilton. Analysis of time series subject to changes regime. Journal of Econometrics, 45(1-2): 39 70, 1990. ISSN 03044076. doi: 00.10160030444066999000000
Kameron Decker Harris, Aleksandr Aravkin, Rajesh Rao, and Bingni Wen Brunton. Timo-varying autore- gression with low-rank tensors. SIAM Journal on Applied Dynamical Systems, 20(4):2335 2358, 2021.
Florian Hess, Zahra Monfared, Manuel Brenner, and Daniel Durstewitz. Generalized teacher forcing for learning chaotic dynamics. In Andreas Krause, Emma Brunskill, Kyunghyun Cho, Barbara Engelhardt, Sivan Sabato, and Jonathan Scarlett, editors, Proceedings of the 40th International Conference Machine Learning, volume 202 of Proceedings of Machine Learning Research, pages 13017-13049. PMLR, 23-29 Jul 2023. URL https://proceeedingg mml. wtt.://eeeiaa.......rre//e////// httpe///v222/heess220 html
David Hocker and Il Memming Park. Multistep inference for generalized linear spiking models curbs runaway excitation. In 2017 8th International IEEE/EMBS Conference on Neural Engineering (NER), pages 613 616. IEEE, 2017.
Armand Jordana, Justin Carpentier, and Ludovic Righetti. Learning dynamical systems from noisy sensor measurements using multiple shooting. arXiv preprint orXiv:2106.11122 2021.
Kadierdan Kaheman, J Nathan Kulz, and Steven L Brunton. Sindy-pi: robust algorithm for parallel im- plicit sparse identification of nonlinear dynamics. Proceedings of the Royal Society A, 476(2222):202000279 2020.
Jan Kamiñikii Shannon Sullivan, Jeffrey M Chung, lan B Ross, Adamn N Mamelak, and Ueli Rutishauser. Persistently active neurons in human medial frontal and medial temporal lobe support working memory. Nalure neuroscience, 20(4):590 601, 2017.
Jan Kaminssii, Aneta Brzezicka, Adam N Mamelak, and Ueli Rutishauser. Combined phase-rate coding by persistently active neurons as a mechanism for maintaining multiple items in working memory in humans. Newron, 106()):256-2444 2020.
In-Bong Kang. Multi-period forecasting using different models for different horizons: an application to economic time series data. International Journal of Forecosting, 19(3))387-400) 2003.
Hassan K Khalil. Nonlinear systems, 2001.
Lawrence R Klein. An essay on the accuracy of economic prediction. In Pwul Samuelson and the Foundalions of Modern Economics, pages 29 69. Routledge, 2019.
Nathan Kutz, Xing Fu, and Steven Brunton. Multiresolution dynamic mode decomposition. SIAM Journal on Applied Dynamical Systems, 15(2):713 735,,,17733--11 2016.
Michael Kyzar, Jan Kamiński, Aneta Brzezicka, Chrystal M Reed, Jeffrey M Chung, Adam N Mamelak, and Ueli Rutishauser. Dataset of human-single neuron activity during sternberg working memory task. Scientifie Data, (1())889, 2024.
Jin-Lung Lin and Clive WJ Granger. Forecasting from non-linear models in practice. Journal of Forecasling, 13(1):1-9, 1994.
Scott Linderman, Matthew Johnson, Andrew Miller, Ryan Adams, David Blei, and Liam Paninski. Bayesian lcarning and inference in recurrent switching lincar dynamical systems. In Artificial Intelligence and Statistics, pages 914 922. PMLR, 2017.
Scott Linderman, Benjamin Antin, David Zoltowski, and Joshua Glaser, SSM: Bayesian Learning and Inference for State Models. httpps//ggithue..cc//aaaaaaaarr
2020.
24

---

Published in Transactions on Machine Learning Research (09/2024)
Scott W Linderman, Andrew C Miller, Ryan P Adams, David M Blei, Liam Paninski, and Matthew J Johnson. Recurrent switching linear dynamical systems. arXie preprint arXiv:1610000400 2016.
Massimiliano Marcellino, James H Stock, and Mark Watson. comparison direct and iterated multistep ar methods for forecasting macrocconomic time scries. Journal of econometrics, (5(12):4499-555 -56,, 2006.
Noga Mudrik, Yenho Chen, Eva Yezerets, Christopher J Rozell, and Adam S Charles. Decomposed linear dynamical systems (dlds) for lcarning the latent components of neural dynamics. Journal of Machine Learning Research, 25(59):1 -4,, 2024a.
Noga Mudrik, Ryan Ly, Oliver Ruebel, and Adam S Charles. Creimbo: Cross ensemble interactions in multi-view brain observations. arXiv preprint ar arxi:22405. 17395, 2024b.
Noga Mudrik, Gal Mishne, and Adam Shabti Charles. SiBBIInGS: similarity-driven building-bloch inference using graphs across states. In Ruslan Salakhutdinov, Zico Kolter, Katherine Hcller, Adrian Weller, Nuria Oliver, Jonathan Scarlett, and Felix Berkenkamp, editors, Proceedings of the 41st International Conference on Machine Learning, volume 235 of Proceedings of Machine Leorning Reseorch, pages 36504-36530, PMLR, 21- 27 Jul 2024c. URL httph: httpps//pproceedinggs mlr. http:///ww22mmmrmr22222
K. P. Murphy, Switching Kalman filters. Technical report, Technical Report, UC Berkeley, 1998.
Josue Nassar, Scott W Linderman, Monica Bugallo, and II Memming Park. Tree-structured recurrent switching linear dynamical systems for multi-scale modeling. arXiv preprint arXiv:1811. 2018.
Jinli Ou, Li Xic, Peng Wang, Xiang Li, Dajiang Zhu, Rongxin Jiang, Yufeng Wang, Yaowu Chen, Jing Zhang, and Tianming Liu. Modeling brain functional dynamics via hidden markov models. In 2013 6th International IEFE/EME /MMBS Conference on Newral Engineering (NER), pages 569-572. IEFE, 2013.
Soumyasundar Pal, Liheng Ma, Yingxue Zhang, and Mark Coates. Rnn with particle flow for probabilistic spatio-temporal forocasting. In International Conference Machine Learning, pages 8336-8348. PMLR,
2021.
Stéphane Ross, Geoffrey Gordon, and Drew Bagnell. A reduction of imitation learning and structured pre- diction to no-regret online learning. In Proceedings of the fourteenth international conference on artificial intelligence and statistics, pages 627 637-655.44 JMLR Workshop and Conference Proceedings, 2011.
Oliver Rübel, Andrew Tritt, Ryan Ly. Benjamin K Dichter, Satrajit Ghosh, Lawrence Niu, Pamela Baker, Ivan Soltesz, Lydia Ng, Karel Svoboda, et al. The neurodata without borders ecosystem for neurophysio- logical data science. Elife, 11:e78362, 2022.
Diya Sashidhar and J Nathan Kutz. Bagging, optimized dynamic mode decomposition for robust, stable forecasting with spatial and temporal uncertainty quantification. Philosophical Transactions of the Royal Society A, (800(229):2021100))
2022.
Shankar Sastry. Nonlinear systems: analysis, stability, and control, volume 10. Springer Science & Business Media, 2013.
John Savill, Ian Dransfield, Chris Gregory, and Chris Haslett. A blast from the past: clearance of apoptotic cells regulates immune responses. Nature Reviews Immmnnlogy, 2(12):965-99)) 975, 2002.
Bongkee Sin and Jin H Kim. Nonstationary hidden markov model. Signal Processing, 16(1):31-465 1995.
Andrew H. Song, Demba Ba, and Emery N Brown. Plso: A generative framework for decomposing nonstation- ary time-series into piecewise stationary oscillatory components. In Uncertainly in Artificial Intelligence, pages 1371-1381. PMLR, 2021.
Richard R Stein, Vanni Bucci, Nora C Toussaint, Charlic Buffie, Gunnar Rätsch, Eric Pamer, Chris Sander, and Joao Xavier. Ecological modeling from time-series inference: insight into dynamics and stability of intestinal microbiota. PLoS compulational biology, ((12)=e100333) 2013.
25

---

Published in Transactions on Machine Learning Research (09/2024)
David Sussillo. Neural circuits as computational dynamical systems. Current opinion in newrobioloyy, 25: 156-163, 2014.
Rahimi Tabar. Analysis and dota-based reconstruction of complez nonlinear dynamical systems, volume 730. Springer, 2019.
George C Tiao and Daming Xu. Robustness of maximum likelihood estimates for multi-step predictions:
the exponential smoothing case. Biometrika, 80(3):623 641, 1993.
Jonathan H Tu. Dynamic mode decomposition: Theory and applications. PhD thesis, Princeton University, 2013.
Vinit S Unni, Ashish Mittal, Preethi Jyothi, and Sunita Sarawagi. Improving rnn-transducers with acoustic lookahead. arXiv preprint urXin:2307.05006 2023.
E Vaadia, I Haalman, Abeles, Hagit Bergman, Prut, Hi Slovin, and AMHJ Aertsen. Dynamics of neuronal interactions in monkey corlex relation to behavioural events. Nalure, 373(6514):555 518, 1995.
Arun Venkatraman, Martial Hebert, and Bagnell. Improving multi-step prediction of learned time series models. In Proceedings of the AAAI Conference on Artificial Intelligence, volume 29, 2015.
Andrew A Weiss. Multi-step estimation and forecasting in dynamic models. Journal of Econometrics, 48 (1-2):135-149, 1991.
Yongqian Xiao, Zixin Tang, Xin Xu, Xinglong Zhang, and Yifei Shi. A deep koopman operator-based mod- elling approach for long-term prediction of dynamics with pixcl-level measurements. CAAI Transactions on Intelligence Technology, 2023.
Enoch Yeung, Soumya Kundu, and Nathan Hodas. Learning deep neural network representations for koop- man operators of nonlinear dynamical systems. In 2019 Americon Control Conference (ACC), pages 4832-4839. IEEE, 2019.
Eva Yezerets, Noga Mudrik, and Adam Charles. Decomposed lincar dynamical systems (dlds) models reveal context-dependent dynamic connectivity c. elegans. bioRziv, pages 2024 05, 2024.
26

---

Published in Transactions on Machine Learning Research (09/2024)
A Appendix
A.1 calculation details for operators differences (Fig. 3C)
We computed the operator differences illustrated in Figure 3C using the expression:
((A 1) (A 1)w) x factor
Here, À represents the operators identified by the different methods, A denotes the ground-truth operators, and 'factor' is a scalar used for visualization purposes only (identical for all methods).
A.2 Regularization options in LTV system experiment we applied the following regularization options in Sec. 6.4
Smoothness: The objective function becomes: A, - arg min |60 A,Til1 AI|A, AL-11
2...T
Sparsity Regularization: A, = arg min |0 - A, s.1. ||Aello <=T
A.3 Hyperparamters used in experiments
Table 1: Hyporparameter settings for DAD bascline in linear experiment
Parameter
Value
Additional Info
seed
random seed
500
number of time points
Wig
weight of la regularization on dynamics
Wdecay
docay of regularization cocfficient over itora- tions
N Niterations
100
number of iterations
init type
step
initialize with 1-step optimization
reweight
False
whether to reweight the observations and the
lookahead during training.
Table 2: Hyperparameter settings for "DAD reweigh" baseline in linear experiment
Parameter
Value
Additional Info
seed
random seed
T
500
number of time points
Wia
weight of & regularization on dynamics
Wdecay
docay of regularization coofficient over itera- tions
Naterations
100
number of iterations
init type
'step'
initialization type for matrix A
reweight
True
whether to reweight the observations and the lookahead during training.
27

---

Published in Transactions on Machine Learning Research (09/2024)
Table 3: Hyperparameter settings for "reweigh fz" baseline in linear experiment
Parameter
Value
Additional Info
seed
0
random seed
T
500
mumber of time points
W82
1
weight of &g regularization on dynamics
Wdecay
1
decay of regularization coefficient over itera- tions
Niterations
100
number of iterations
A init type
'step'
initialization type for matrix A
reweight
True
whether to reweight the observations and the lookahead during training.
Table 4: Hyperparameter settings for LINOCS in linear experiment
Parameter
Value
Additional Info
K
80
maximum training order
with offset
True
whether to look for an offset term b
cal offset
True
whether to calculate offset
weights_sste
exponential
weight style for orders
Fu
exponential param-
0.01
eter for the weights
infer b way
'each'
approach to infer the onset. Based on each order.
K₂
maximum lookahead training order for the off- sets
weights style b
'exponential'
weights style for the offsets
A.4 Lorenz equations
The Lorenz attractor follows:
dæ oly x), dl
dy
dt x(p 2) y. dz
dt zy $2,
where y, and represent the state variables, and 0, p. and are the system parameters, set to 10, 28, and B = 3
A.5 Calculation of Operator and Offset for Linear Experiment
We bogin by assuming that X represents the observations in this subsection. Our analysis involves two main objectives: estimating the transition operator A and the offset b.
We define set of matrices { where each th is found by solving:
4/ arg min || X tb v/x (xl [1]xx-=+]]]
4
Here, [ is row vector of ones, and [X. - ki [1]xxT=]] is the vertical concatenation. Each thy captures 1+1 + +K A b,,, b, with the last column of th specifically capturing Σi, A Alb. b.
We extract all but the last column from each tha
À arg min A
(zbb)...pl?
A
28

---

Published in Transactions on Machine Learning Research (09/2024)
Table 5: Hyperparameter settings for dLDS experiment
Parameter
Value
Additional Info
K
50
maximum lookahead order
additional update
True
whether to include additional update step value of & the coefficients for the 1-st iter- ation
4 init
0
max iters
200
maximum number of iterations
||
0
&g norm on the basis of dLDS dynamics frequency of updating {f,)
frequpdater
Ozdecay
0.99
decay of the l norm the coefficients
decay
0.9999,
decay of the & norm on the coefficients
Wemoothdeccy
1.01,
decay of smoothing weight on coefficients
1 smnooth time
1.1,
regularization weight on coefficients smooth- ing
L2
0,
Eg on the coefficients smoothness on C₁ over time std of the noise to add to coefficients during training.
Wimoothiime
0.1,
nooigge
0.05
2.5,
& regularization on C₁ wether to decorrelate the basis dynamics maximum interval to increase l during train-
decor
False
max interval
ing
to norm
True
whether to normalize the basis dynamics
370
(0.2160895,
0.97627445,
ground truth initial state (at = 0)
0.00623026]
number of basis dynamics operators (fs)
Table 6: Hyperparameter setting for the LTV-LLrenz experiment
Parameter
Value
Additional Info
K
maximum lookahead order
Wsmooth
2 (or 20)
smoothness weight
eslimalestres
2
(hreshold to increase value
with future costr or
True
consider both fast and future reconstructions
with future
Truc
wether to apply time smoothness to future weights
W Wsmooth future weighlsstyle
0.05
weights future smoothness
'uni'
Uniform weight style.
|wllF
Frobenius norm weights
init. style
'step'
initialize based on optimal 1-step
maFseers
maximum number iterations
errorthres
threshold to increase error whether to include reverse update how many zeros (relevant only to the sparse networks)
withieeerre
False
N Nzeros
or 2 or
This operation isolates the transition operator by minimizing the Frobenius norm of the differences across all k levels.
29

---

Published in Transactions on Machine Learning Research (09/2024)
Using the last columns of the original &b matrices, estimate b, considering previously estimated A. We model each observation ae as:
Lt ~ Ax +
(14)
Expanding recursively for & steps, we derive:
k-1 Dt R x Σ
b
(15)
k,=0
To organize our data for analysis, we define shifted observation matrices X+* = x and corresponding
estimated matrices X+ = A X These concatenated vertically for all is = 1. K, adjusting the dimensions as necessary with NaN padding.
Additionally, we construct A Rxxx a vertical concatenation of (2 A K
Jmml
We solve for the offset b by minimizing:
b = arg min |X
X+ AbI
(16)
This minimization is achieved via simple least squares, and b is then averaged across times.
B Explanation of Baselines for Comparing LTV-LINOCS
We compared the LTV-LINOCS system against two other multi-step approaches, each with various parameter
settings:
SSM-RNN Jordana et al. (2021): : We used the code from the GitHub repository linked to the original paper Jordana al. (2021). We experimented with different numbers of epochs, network architectures (fully connected vs. locally linear), seeds, and numbers of shooting points. Details are provided in Table g
GTF-shPLRNN Hess et al. (2023): : We utilized the code available the GitHub repository here,
which is linked to the original paper Hess et al. 2023) and written in Julia. We based the parameters on the values defined in the file ShPLANNP http://aaaeeperrreeeeee http://pprraaaaaa/////////////// __CTF. jl but modified the observation model to the idenlity.
It important note that while we did not perform exhaustive hyperparameter tuning, we examined 12- 16 hyperparameter combinations around the default settings. Details of the hyperparameters are listed in Tables 7 and 9
For reference, used the ground truth version of the Lorenz attractor, shown in Figure and present the results in Figure 20
C
Computational Complexity
C.1 Linear System
Given that æ E NXT and A E x computing the highest power K of A (A*) through repeated squaring involves matrix multiplications, each with a complexity of O(N3), resulting in O(KN3).
If the optimizer is set maximum of M iterations, then the total complexity O(MKN3). In the lincar case, the addition of the offset term integrated into the inference by extending the size of A by which does not affect the computational complexity scale (assuming N >> 1).
30

---

Published in Transactions on Machine Learning Research (09/2024)
Table 7: Parameters for GTF baseline options
option 3
option.
option.
option
12
forcing
lafest_
16
16
16
16
scalar saving Interval
500
500
500
500
500
500
500
500
500
500)
50H
500
Luncher luecing Interval
16
10
16
15
16
16
16
16
10
16
16
16
guussion noise level
B
glf
TRUE
TRUR
TRUE
TRUE
TRUE
TRUE
TRUE
TRUE BADAM
TRUR
TRUE RADAMM
TRUE
TRUE
RADAM
RADAMI
RADAM
RADAM
RADAM
RADAM
HADAM
RADAM
balck sise
16
16
16
16
start
800
0001
0.001
0/001
600
0.001
0.00
0000
(0001
00001
0000
path to Impuls
gff alpha
0.8
088
688
689
0.9
hustehow
10
50
50
90
50
Alplo derety
model
Reutily
Menttyy
Mentty
Fentity
Keatlly
Mentes
Identity
Meatity
Mentley
Beneity
Mentity
las. modial regularbaatin
0.3001
0.009
0.000
(10001
010001
0100/1
(10001
D slop biss
30
10
30
30
30
PE B
20
20
30
20
20
20
30
201
20
30
30
end Ir
Te(l
Tedii
Te4ll
Te:li
Te-00
Te466
1e:06
Te:1i
Texiii
Textl
Te:0i
Texti
device
apu
gradieut clipping Dorm
slap scaling
image soving interval
500
500
500
500
500
500
500
500
500
300
500
500
fum Inooe
W
50
50
90
50
30
50
301
14)
hinklon alm
M1
512
112
612
512
512
612
512
512
312
1122
512
50
ula mantal regalarizntion
1o-06
1040
10006
alpha updatt Interval
oportos
M00
50.0
5000
9:00
3000
3000
5000
9000
50.00
30000
5000
10000
MAR lambida
11
amoothing
opoch
Table 8: GTF Model of Each option model
option 1
clppeeSSballo.PLRNN
option 2
clipeddShallow.PLNN
option 3
ClipeeSSballow.PLRRNN
option 4
clipeeSShallowPPLRNN
option 5
clippedShallowPLLNNN
option 6
BiypeeSShlllowPIRNN
option 7
shallowPLRNN
option 8
shallowPLRNN
option 9
option 10
shallowPLRN
option 11
shallowPLRN
option 12
shallowPLRNN
C.2 dLDS
The direct dLDS version we propose in the paper works directly the observations rather than a latent low-dimensional space. This process involves both identifying the dynamic operators (fi)/= and their coefficients ((c:)T=))
We will start with the complexity for the coefficients inference. Particularly, the inference of e₁ for each 1... .T involves: C. arg mitte | F c. l7 where FK E R(K+1lpxJ and (#t+1)oer: R Hence, assuming {K+1)p >>> J, performing the above least squares to infer each requires ((J Op) assuming K > 1 for the pscudo-inverse and O(J(K + 1)p) = O(JKp) for the matrix-vector multiplication, resulting in overall Ky) for each time point, and O(MTJ'Kp) for T time points and M overall model iterations.
To infer the elements of the set (fz)
Faul = arg min X Faul (x. Fall
31

---

Published in Transactions on Machine Learning Research (09/2024)
Table Parameters for RNN Multiple-Shootiae baseline options
num shooting
len
filt.
max
it
argT
alpha
log name
bs
seed
num epochs
Option 1
100
50
10000
100
locally linear
40
5
Option 2
100
50
210
100
locally linear
40
4
1000
Option 8
100
50
210
100
locally lincar
40
1000
Option 4
200
100
50
10000
100
40
4
5
Option 5
100
50
210
100
locally lincar
40
4
1000
Option 6
100
50
100
locally linear
40
1000
Option 7
100
50
210
100
locally linear
40
4
100
Option 8
100
50
100
locally linear
40
Option 9
100
50
210
100
locally linear
40
480
100
Option 10
100
50
100
locally linear
40
480
Option 11
100
50
210
20
locally linear
40
480
100
Option 12
100
50
locally linear
40
480
100
Option 13
100
50
210
locally linear
40
480
100
Option 14
100
10
210
fully connected
40
480
100
Option 15
100
10
210
fully connected
40
480
1000
Option 16
100
10
210
100
fully connected
40
480
1000
Table 10: Parameters for LINOCS baseline options for comparison figure
Option 1
Option 2
Option 3
Option 4
Option 5
Option 6
Option 7
Agmooth
0
0
0
0
0
2
20
num zeros
0
0
1
2
3
0
0
max K
8
3
4
4
4
5
5
max iters
15
15
25
25
25
40
40
where X E (T+()+I) and X. E RDJxT Assuming pJ < T, the overall complexity will be J(P for the the useudo-inverse step and O(TpJy) = O(Ty²J. Hence, for each iteration O(p2 and overall
(2T.
Hence, assuming K > p. the overall complexity of the dLDS direct LINOCS extension would be O(MT maxly, K)) pMTJ2pk).
C.3 Locally Linear System (Time Invariant)
Assuming we limit the number of iterations to M. For each iteration 1. M, we iterate over 1...2 T time points to infer A, For each A, we consider K + 1 windows (k = 0... K) with ką 1...k + shifts. Hence, under each combination of time and iteration m, we get an overall
K k+1 (K+ [X++)[1K++)
(K + 1)(K + 2)
K
equations for each A, For simplicity, we will neglect the weights wy in the following computational com- plexity analysis since they do not add complexity and do not nee to be considered here.
Let:
k+1 x =+ diag (( 11
limo TRKK xk*) (each æ, RPx1).
x diag( (( 4+1
Jk= e Bok" xoK"
C' Bpk xp be vertical concatenation of the set: (1T1/=1 ki A, k.+ 1k+1 1
B p2xpk' be an horizontal concatenation of the set: i A 1
V = -B RoxK"
32

---

Published in Transactions on Machine Learning Research (09/2024)
csserved
c5o
rewelsh::ss
reweigh far:50)
Figure 12: Zoom in for the roconstruction of the lincar experiment.
v E (px((+p)) be V Lxppl, namely, horizontal consternation of v and Lxxp where X is the weight of smoothness regularization on A.
at be an horizontal consternation of x+ and x @ Ar-1 where @ is the Kronecker product, 11 a vector of ones with K' elements, and R+ is the weight of smoothness regularization on A,.
Each A, can thus be solved via extended least squares:
A = arg min A;
x
CA,Be 11 A||A, - A1-1ll?
arg min
CA, V + +1IA, Ar111
A;
= arg min A
x
CA.TE
Assuming K > p, the overall process thus involves the following computational steps:
Calculating (K*)))+)))
2. Pseudo-Inverse of v: O(1) +)) = OWPK*.
3. Pseudo Inverse of C: Opp² Kp) O(Kp').
4. Matrix multiplication of pinv aC) pinv (V) (1P((*)3)
Overall, the most computational complexity in each iteration and time step ((p((K+))) Since we have T time points and assuming that we limit the number of iterations to M, the overall complexity would be:
OMTTpR3
C.4 Observation vs. Latent Spaces in Dynamical Systems
In many dynamical systems models, the observation space, denoted as S R.N contains the directly measured variables. For a system observed over time, we represent these observations as æ; X, where denotes the time index. The latent (low-dimensional) space, denoted Z C RM with M s N, contains hidden variables that are not directly observed but inferred from the observations. Let E Z represent the latent state at time t. The relationship between the obscrvation space and the latent space is often modeled by an observation function h Z + X, such that:
æc = 1122)))) +
where & ~ _N(00 021) represents observation noise. Several models, including SLDS and dLDS, include such transition to a latent state to capture the underlying dynamics low-dimensional space. In this paper, by fixing the observation function the identity operator, in essence learn the dynamic evolution directly the observations.
33

---

Published in Transactions on Machine Learning Research (09/2024)
C.5 More Information about Neural Data
The data we used was collected by the Rutishauser lab at Cedars-Sinai Medical Center Kyzar el al. (2024) with detailed descriptions in Kyzar et al. (2024); Kamiński et al. (2017) 2020). The dataset includes electro- physiological recordings from 21 epileptic patients who were implanted with depth electrodes and Behnke- Fried microwires in the human medial temporal lobe and medial frontal cortex. The recordings were obtained while the subjects performed a memory task (*Sternberg task"). In the Sternberg task, participants were required to memorize set of 1 images. These images were pseudo-randomly chosen from group of five that elicited the strongest selective responses during the screening task. Following maintenance period, participants were shown probe image and asked to identify whether it was included in the initial set. The images used in these tasks represented broad array of subjects and complex natural scenes. Please see Kyzar et al. (2024) for more details.
We loaded the data from the DANDI Archive an NWB (Neurodata Without Borders) format Rübel ct al. (2022)), ), and used single session from it. This session includes recordings 63-year-old male subject (Subject 10 in the data) recorded in June 2023 while performing the Sternberg task.
C.6 Pre-processing Neural Data
To process the data, we took the spike times of the p = neurons within the chosen session. We then convolved the spike times with 100-ms-width Gaussian kernel to get firing rate estimation. We then normalized the data by dividing cach neurons' estimated firing rate by its top 99% firing rate, and used it the initial 850 samples.
34

---

Published in Transactions on Machine Learning Research (09/202))
A
8:22
0.72
6:75
B
E
Fredicttoo Order
220
zido iloo ald Picdiction (Smo)
Figure 13: Linear System under Structured Noise. A: Real vs. identified operators and oflsets. B: Quiver plots of real and identified operators (lower 2x2 sub-matrix, ie., focusing the rotational part of v2-wa) present patterns that appear similar, rendering it challenging to discern differences when examined in isolation.
C: The differences in effects between real operators and inferred operators highlight how minor distinctions in dynamic operators gain prominence during lookahead reconstruction (calculation details in A.Ij. D: Full lookahead reconstruction (ground truth vs. baselines) shows swift convergence to the circle's center for the one-step optimization results due to small differences in dynamic values (mid-yellow subplot) and divergence for DAD-based results (three most-right subplots). E, H: MSE under increasing prediction orders. LINOCS achicves better (lower) MSE compared to 1-step optimization. F, G: LINOCS reconstruction compared to 1-step optinization under increasing noise values reveals that LINOCS maintains good reconstruction even under extreme noise conditions. I: By propagating identified operators until a relative reconstruction error of ~ 1, LINOCS enables future predictions full rotations), contrasting with immediate convergence in one-step optimization. Black indicates error differences bet ween one-step optimization and LINOCS.
70,000 time points (~ 4280
35

---

Published in Transactions on Machine Learning Research (09/2024)
B
-18.0
125.0
17.35
38
A22
52g
575
118.00
17.35
Arg
Akg
bxg
F
Prediction Order
Prodiction Order (Uiee)
Figure 14: Linear System for 3D Cylinder. A: Real vs. identified operators and offsets. B: Quiver plots of real and identified operators present patterns that appear similar, rendering it challenging to discern differences when examined in isolation. C: The differences in effects between real operators and inferred operators highlight how minor distinctions in dynamic operators gain prominence during lookahead recon- struction (calculation dotails in A. D: Full lookahead reconstruction (ground truth operators vs. basclines) shows swift convergence to the cylinder's center for the 1-step optimization results due to small differences in dynamic values (yellow) and divergence for DAD-based results (three most-right subplots). E, G: MSE un- der increasing prediction orders. LINOCS achieves better (lower) MSE compared 1-step optimization with perfect full lookahead reconstruction under high-enough training order (E right bottom). F F,H: LINOCS reconstruction com compared to 1-step optimization under increasing noise values reveals that LINOCS main- tains good reconstruction even under extreme noise conditions. I: Propagating the identified operators until reaching a relative reconstruction error of ~ shows that LINOCS identifies operators that enable future prediction 70,000 time points (~ 2680 full rotations) before converging, unlike one-step optimization that converges immodiately. black: error difference between one-step optimization and LINOCS.
36

---

Published in Transactions on Machine Learning Research (09/2024)
A
B
D
5
X1
X1
0.61
0.59
K1
0:81
0.59
-05
40
0.61
0.59
K2
Kgm
0.51
000
2.0
8.
X1
0
0.39
0.31
X5-
5/
0.81
Xs-
0
X1
23 0X
X2
X1
X5
X3
X1
X2
Xg
5X
8X
Time
Figure 15: Ground Truth operators and coefficients for the SLDS experiment. A: The ground truth basis dynamics operators {ff 1 consist of rotational matrices oriented in various directions. B: Ground truth operators' coefficients (c). C: Ground truth state æ. D: offset applied to the operators.
SLUSS Cround Truth
{Y} (for f1)
{Y} (for f2)
{Y} (for fz)
Imagt()
0.5
(magii) (moz(A)
Ré(A)
Re(A)
Re(A)
M Mleew
Imaguh)
0.5
0.0
-0.5
ReiA)
Re(A)
Re(A)
SLUSSO
ImagiAl
(mogli)
Re(A)
RelA)
Re(A)
Imag(A)
05
(moglA)
-0.5
Reli)
Re(A)
Re((A)
imagiAl
0.5
-0.5
Re(A)
Re(A)
RE(A)
L00
ImaplA)
Imag/A)
Rdix)
RclA)
Re(A)
Imag
0.5
(mogtA)
0.0
-0.5
Re(A)
Re(A)
Re(A)
LINOCS LINOCSSS
(magiA)
(msolA)
ReiA)
RC(A)
Rd(A)
imag(A)
(mapli)
(moplA)
0.0
-0.5
-02
0.0
0.
1.0
-0.2
0.4
0.6
-02
0.0
0.2
1.0
Re(A)
Re(A)
Re(A)
Figure 16: Eigenvalues of identified operators by the different systems compared to the eigen- values of the ground truth operators. Rows represent different methods, with LINOCS (with different parameter combination) in the four last rows. Columns represent the three eigenvalues of each of the three different x 3 linear operators. LINOCS enabled the identification of almost perfect eigenvalues while the other methods found at least one wrong eigenvalue per operator, explaining the decaying docaying/divorgenoe of their reconstruction.
37

---

Published in Transactions on Machine Learning Research (09/2024)
A
B
f
10-
X5
0.99
X2
0:31
X2
Kg
0.11
0P
-22
X5
024
-03
-02
X2
X3
X3
X1
200 0/9 300
1804
0X
Figure 17: Ground Truth operators and coefficients for the "pseudo-switching"" dLDS experi- ment. A: The ground truth basis dynamics operators {fyly= consist of rotational matrices oriented in various directions. B: Ground truth operators' coefficients (c). C: Ground truth state a. The colors of the markers are a time-changing re-weighted average of the active operators.
GLI
41
X2n
83
M2-
01
-45
X3
43
X3+
-12
X1
X2 6X
X5
250 193 190 100
6X
6X
Figure 18: Ground Truth operators and coefficients for the second dLDS experiment. A: The
ground truth basis dynamics operators {f))g=] consist of rotational matrices oriented in various directions.
B: Ground truth operators' coefficients (c). C: Ground truth state æ.
A
B
20
10a
10a
log MSE)
Reration
10
-10
10-1
10
10.5-
10
10
0
38
76
114
152
190
223
266
304
342
380
250
284
552
609
846
772
760
790
636
674
Iterations
C
fime
Iteration
Iteration
Iteration
Iteration
Iteration 5
500
600
500
000
2000
400
100
20000
3000
200
-10
0
10
-20
-50000
10000
5002
2000
4000
2000
20 e0 000
750
3309000
0000
-100
-200
Iteration
Iteration
Iteration 8
Iteration
Iteration 10
X
89
X
20 o 22
70
20
-10
-20
-30
-20
Figure 19: Demonstration of LINOCS' effect over iterations for Lorenz attractor with 900 time points (Kurain = 19). A: MSE over itcrations (curve corresponds to median values; shade represents 25%-75% percentiles over lime). B: MSE over lime-points (horizontal) and training iterations (vertical). C: Full lookahead reconstruction based on operators identified under different iterations. Dotted black: (as exemplified in the top left) ground truth; Solid Cyan: Full lookahead LINOCS reconstruction.
38

---

Published in Transactions on Machine Learning Research (09/2024)
A
Ground Truth
y 10-
SSM Multiple Shooting
LINOCS
37
B
LINOCS Option
LINOCS Option
LINOCS Option
LINOCS Option 4
LINOCS Option 5
LINOCS Option 6
LINOCS Option 7
X
x
50
50
50
50
25
0
25
25
25
25
25
-20
X1
25
20
X1
25
20
1
25
20
1
254
20
X1
20
1
0
25
43
20
1
0
C
D
Shooting Option
Shooting Option
Shooting Option 3
Shooting Option
GTF Option
GTF Option
GTF Option
Shoating Option
Shaoting Option
Shoeting Option 7
Shooting Option
GTF Option
GTF Option
GTF Option
Rtors
10
Shoating Option
Shoating Option 10
Shosting Option 11
Shaoting Option 12
Option
Option
Ailer 202
Shooting Option 13
Shooting Option 14
Shooting Option 15
Shooting Option 16
Option
GTF Option
Option
Anas
Figure 20: Comparison of LINOCS Lorenz locally linear time-varying results to other multi- step approaches. A: rMSE of full lookahead reconstruction by LINOCS (blue) compared to the other baselines (Jordana et al. 2021 (Shooting, red), and (Hess et al. 2023 (GTF, orange). B: LINOCS full lookahead reconstruction under different hyper-parameter settings, with the options described in Table IO C: Shooting (Jordana et. al. 2021 full lookahead reconstruction ion under different. hyper-parameter settings, with the options described in Table g D: GTF (Hess et al. 2023 full lookahead reconstruction under different hyper-parameter settings, with the options described Tables and
39

---

Published in Transactions on Machine Learning Research (09/2024)
A
LINOCS
B
Original SLDS (1 step)
Arrysaain PUGH SFFFER
Dorsal Srttrro
Hippocampia EFFT
PCC suppp parplemett
C
antertor
Dorsal
supplemettary mator area
D
Amygdala
Amyodala
cingulate
cortex
Hoppocampus
LINOCS SLDS
Origiall SLDS
1970 199.000
Figure 21: Operators and reconstructions by the LINOCS-driven SLDS compared to classical
SLDS. A: Dynamical operators identified by LINOCS-driven SLDS. B: Full lookahead reconstruction.
40

---

Published in Transactions on Machine Learning Research (09/2024)
Dorsal
Pre
anterior
supplementar
cingulate
motor
Amygdala left
Amygdala right
cortex
Hippocampus left
area
left
left
10
8
14
Figure 22: Classical SLDS results (non LINOCS) the real world data. The identified net- works ((f))/)) by the non-LINOCS SLDS code [Lnddemaan et al. 2020 for each region in the real world data Kyzar el al. 2024h.
41

---

Published in Transactions on Machine Learning Research (09/2024)
A
LINOCS
B
Classical SLDS
p Pasa 3 B2aSE
COIOOO enterne
mooor
C
supplementary
D
cingulate
metter
M0M SOOPP $LDS
Amyyddala
Amyydalaa
cortex
Hoppocampus
UNNOS
SOTO COOKTT
30
40
Figure 23: SLDS results real data. A: Identified operators per region by LINOCS. B: Identified operators per region by classical SLDS. C: Switch times by LINOCS-SLDS classical SLDS. D: Number of switches for LINOCS-SLDS vs. classical SLDS.
DAD
DAD reweigh
reweigh t
1023
7.22666 <00
7.2116
M MM W
x.34x100
7.3
722156
722156 10
728
733180
124 x10
26
722196 x301
85 d:
lerations
lterations
Figure 24: Lookahead prediction error x of DAD-algorithms over training iterations (for the white noise linear experiment).
42


Reviews:
[{"review_id": "UkjicfOXTc", "reviewer": "Reviewer_nQ8x", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "# Strengths\n\nIn general, I think the paper is well-written and very easy to follow. Moreover, I'm *very* impressed the breadth of the experimental section which I think speaks to the broad applicability of LINCOS.\n\n# Weaknesses\n\nI think the two biggest weaknesses is the novelty and positioning of the paper; I will first focus on the novelty aspect. Training dynamical systems---and state-space models in general---using Iterative Multi-Step Prediction (IMS) is not novel and has been done before. For instance, in [1] IMS was used to train deterministic dynamical systems while in [2] and [3] IMS was used to train state-space models. Thus, it is imperative that the authors discuss the differences between the LINCOS and the studies I have linked as I think that there is substantial overlap between them. Moreover, I think it is imperative that the some/all of the methods should be compared against LINCOS in the experiments.\n\nNext, I think the positioning of the paper needs to be revamped. LINCOS is a general learning algorithm that can be applied to arbitrary dynamical systems, though the authors focus on linear dynamical systems and its extensions. This is important to note because the authors state in the introduction that \"LINOCS also avoids relying on massive amounts of data (like RNNs, for example, require)\", where the massive amount of data is really a comment about the parameterization of the dynamical system (RNN) and NOT about the learning algorithm used. For instance, LINCOS can be applied to RNNs but still might require a massive amount of data due to the sheer number of parameters in the RNN. I am also making this point as methods proposed in [2] and [3] can also be used for linear dynamical systems as they are, just like LINCOS, learning algorithms that can be applied to arbitrary dynamical systems. Thus, what needs to be made clear in the paper is what is LINCOS solving that the previous papers have not? For instance, it might be the case that previous approaches don't empirically perform well on linear dynamical systems and that is what LINCOS solves. Or it might be the case that adaptive weighting in LINCOS allows for more stable training compared to previous works.\n\n\n# References\n[1] *Learning Dynamical Systems from Noisy Sensor Measurements using Multiple Shooting*. Armand Jordana, Justin Carpentier, Ludovic Righetti. https://arxiv.org/abs/2106.11712\n\n[2] *Learning Latent Dynamics for Planning from Pixels*. Learning Latent Dynamics for Planning from Pixels.* Danijar Hafner, Timothy Lillicrap, Ian Fischer, Ruben Villegas, David Ha, Honglak Lee, James Davidson. https://arxiv.org/abs/1811.04551\n\n[3] *RNN with Particle Flow for Probabilistic Spatio-temporal Forecasting.* Soumyasundar Pal, Liheng Ma, Yingxue Zhang, Mark Coates. https://arxiv.org/abs/2106.06064", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "In this work the authors introduce LINCOS, a multi-step ahead training procedure for training dynamical systems. While in principle LINCOS could be applied to any arbitrary dynamical system, the authors focus on linear dynamical systems, switching linear dynamical systems and time-varying linear dynamical systems.", "claims_and_evidence_yn": "Yes", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "1. Discussing similarities and differences between LINCOS and previous IMS training procedures. [Critical]\n2. Including comparisons to previous IMS training procedures. [Critical]", "broader_impact_concerns": "N/A", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "Y1OANRWbcr", "reviewer": "Reviewer_MKNA", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "Strengths:\\\nThey provide a whole host of evidence and figures demonstrating the utility of their method. They apply their method on a wide variety of simulations as well as one real dataset. They provide code for the community to recreate results.\n\nWeaknesses:\\\nCould be more concise. The figures especially are quite small and hard to read. It is also lacking in some references. Specifically, there is a whole host of methods analyzing dynamical systems: DMD (OPT-DMD, BOP-DMD, NS-DMD, MR-DMD, etc.), TVART, PLSO, to name a few.", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "The authors present a new method, LINOCS, to analyze and model non-stationary and non-linear dynamical systems. They demonstrate their method on a variety of simulations, which outperforms previous methods in prediction accuracy. Particularly, it focuses on combining multiple multi-step predictions to create robust estimates of the states, allowing for robust modeling and forecasting. They further demonstrate it's utility on a published neural dataset, where it identifies possible underlying neural connections.", "claims_and_evidence_yn": "Yes", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "Critical:\n1. Add citations to \"common\" methods. For example, there is a whole host of methods, such as DMD (OPT-DMD, BOP-DMD, NS-DMD, MR-DMD, etc.), TVART, PLSO, to name a few.\n2. Background \"(e.g., if $x_t | x_{t−1}$ has a variance of $v_t$, then $x_t|x_{t−K}$ will have a variance of $\\sum_{k=0}^{K-1} v_{t-k}$ due to the property of variance summation).\" Doesn't this assume that all samples are independent from the previous ones? In which case this isn't quite right.\n3. Sec 2.1, Mention SLDS and others, but this should be cited\n4. Sec 2.1.1 \"...Koopman operator for improved long-term prediction. However, these methods often require extensive learning data and remain uninterpretable in the “network” sense, as non-linear state measurements in Koopman operators can obscure understanding of pairwise state interactions.\" Koopman operator methods (e.g. something like data-driven DMD approaches) specifically target finding interpretability and understanding of pairwise state interactions. This should be clarified.\n5. Sec 2.1.1 \"...such representation does not provide explicit insight into the time-changing interactions between the state elements.\" It does if you choose the correct equations, but this is hard to do. Again, this should be clarified.\n6. Sec 2.1: \"...“network” interpretability. Expand upon this. What does network interpretability mean here?\n7. Sec 2.1: \"...observation space, essentially enforcing the transition to the latent space to be the identity operator.\" First time talking about observation and latent spaces; should expand on what these are.\n8. Fig 1. Is this a cartoon? Or how is this generated?\n9. Sec 3. \"[LINOCS]...yields not only a more accurate...\" This model is \"better\" compared to what?\n10. \"Specifically, it gradually increases the weight of the best lookahead reconstruction until convergence conditions are satisfied.\" From the text, it looks like $w$ is being learned while simultaneously finding A? Please add an algorithm to help aid with these descriptions.\n11. Fig 1. Could use more description. E.g. what are x1 and x2? Which models are you using that perform well on 1-step predictions? B, An algorithm approach could be extremely helpful to explain this. C, What is being shown? What is \"Iter?\"\n12. Eq 9: When do you update c vs f vs w (if included)? Since $F_{all}$ depends on $c_t$, while $c_t$ depends on $F_{all}$? An algorithm may be useful here too.\n13. Figures are extremely small, such as figures 4, 7, 11, 12, etc. These should either be enlarged or condensed for better readability. E.g. in Fig 7A, it's extremely hard to tell by eye if one trace is better than another.\n14. Sec 4.3 \"This underscores the necessity of multi-step reconstructions for accurately estimating dynamical system dynamics.\" This just shows that this method works while some others don't, but this doesn't prove that multi-step reconstructions are \"necessary.\" Might want to reword.\n15. Sec 4.4 \"We compared LINOCS-LTV with several other LTV solvers with varying constraints, including smoothness and sparsity.\" These should be listed.\n16. Sec 4.5 Please provide a little bit more description of the task for the monkeys. E.g. \"screening task?\" Also, the results presented should be task aligned; what alignment are you using? E.g. in Fig 8, when does this timing occur? Is there more data that is excluded from analysis?\n17. Sec 4.5 \"..this suggests that LINOCS may better capture the linear approximation of brain activity than the common 1-step optimization.\" This is really big claim without the necessary evidence to back it up. Ideally, would want to compare to standard neuro methods and show that your method can recover standard results (e.g. from granger causality) while something like 1-step optimization cannot. Since this dataset was already published, is there some result that can be recovered?\n18. Fig 8A: Might want to make it clear in the text that each node is a neuron. And how does the region to region connectivity work? Is it averaged across neurons? Is there any normalization applied? \n19. \"LINOCS not only achieves more precise full lookahead reconstruction compared to baseline methods but also successfully retrieves ground truth operators in synthetic data...\" This seems to imply that other methods do NOT recover ground truth operators, which isn't necessarily true. Might want to rewrite this to clarify this point.\n20. \"[LINOCS]...may help mitigate issues such as vanishing or exploding gradients.\" This is a fairly big claim, what evidence is there to support this so far?\n21. Fig 16/17/19: Given that all three rotational matrices (all in different directions) occur, shouldn't there be THREE \"rings\"?  Why are there only two? As a slightly more general test, the three matrices could be rotations in three non-orthogonal directions, to show that the transitions are found regardless if there's an overlap in transition matrices.\n22. Fig 20: What is the prediction and training order?\n23. The following items should be added/discussed in the discussion section: \n24. \"Additionally, our framework allows custom weight functions that suit their specific needs. In the experiments presented in this paper, we concentrate on showcasing three specific options for the weights...\" How does one choose a weight function?\n25. You mention 4 different methods in Sec 3. (pure linear, switching linear, etc.). How does one choose which method to choose? What happens if you choose a more general method on a simpler case?\n26. For the neuro applications, which method should be used? You applied all four methods, but does one seem like the most appropriate for neuroscience? Does each one provide different information?\n27.For the SLDS like methods and applications, how does one choose the number of modes? Especially in the neural application, wouldn't this number affect results?\n28. How does one choose hyperparameters in general. For example, how should one choose K?\n\nTo strengthen work:\n1. The term \"often\" (and other similar terms) is overused and citations should be added to help the reader. E.g. Intro P2: \"often rely on either “black-box.\" Intro P2: \"often result in undesired divergence away from the system.\" Background P6 \"they often struggle to fully reconstruct the dynamics over longer time-scales.\"\n2. Intro \"...in cases where the dynamics are non-stationary, non-linear, or otherwise constrained in ways that reflect real-world system behavior.\" What other ways exist?\n3. Intro \"...achieving significantly improved accuracy in operator identification and long-term predictions.\" Should say compared to what. Also, it may help to reference later sections (e.g. Sec. X) so one can easily find this information.\n4. Sec 2.1: \"One common method for determining A (and b if it is assumed that an unknown offset exists).\" In a lot of methods (e.g. PCA), the data is normalized such that \"b\" doesn't exist. Could say something about this decision. Also, this was referred to earlier in the manuscript as a \"forcing term,\" but now this is just an offset?\n5. Eq 6. Is $w_k$ normalized?\n6. Fig 2. Far right panel, do the coefficients have to be roughly constant? Or can they fluctuate in more general forms?\n7. \"...however LINOCS can be easily adjusted to other noise statistics.\" How?\n8. Eq in Sec. 3.1: A is raised to an exponent right? Not an index, could be helpful to specify this.\n9. Sec 3.2: SSM framework? And this could use a little bit more description. This first uses SLDS to find the boundaries of linear periods, and only then finds the linear dynamics? Can you summarize the method here? How do you deal with Kth order predictions, which may extend past the boundary?\n10. Comment on Fig 3: not to add another model, but OPT-DMD should also work extremely well here.\n11. Fig 3D right 3 panels: Can you interpret this? From Fig 3B, I would expect the values to still \"rotate.\" Why does it diverge?\n12. Fig 3E: What is theoretically best considering the level of noise?\n13. Sec 4.1: you mention using DAD with 100 iterations. Is this enough to converge? What would happen if you tested with more iterations?\n14. Comment on the statement: \"...instead decay to zero (away from the real system).\" In fact, this is known in DMD (Dawson 2016), and is fixed with OPT-DMD.\n15. Sec 4.1: \"...approximately 35,000 time points.\" What is the inherent time scale? E.g. how many rotations is this for in the plots, or what is delta t / the sampling rate?\n16. \"We further tested LINOCS on linear systems with structured noise (Fig. 12).\" This may be better included as SI material instead of the main text? Also, the cylinder case. What about using a much larger system? E.g. around 100 dimensions?\n17. Fig 4:  It looks like SLDS and rSLDS improve with training epochs. If training for much longer, would they find the right transitions? And if they do find the right transitions, how do the algorithm speeds compare?\n18. End of section 4.2. This is the first point when the importance of the eigenvalue spectrum is mentioned, which feels out of place at this point. Could specify that this is an objective earlier? \n19. Sec 4.3 \"Training with orders...resulted in high-quality full reconstruction\". While impressive, I'm not sure this is a \"perfect\" full reconstruction. May want to reword.\n20. \"In addition, we analyzed operators identified across various training iterations of LINOCS to assess their proficiency in achieving lookahead reconstruction (Fig. 20).\" What are the prediction and training orders? Hard to tell how impressive the result is without this information.\n21. Fig 6H, which prediction order is used? Also, a similar plot would be interesting for looking at predictions when switches occur, since a switch could be impossible to predict. Or, perhaps a plot of the prediction order vs training order matrix. I would expect to see the error increase with lower training orders and larger prediction orders.\n22. Fig 6A, Why does K pred of 50 look worse than 800 or 1000?\n23. Fig 8C: Why does the 1-step approach fail to capture the average activity? Shouldn't this be an area where the 1-step approach works well? Also, typically in neuro applications, the data is normalized; how does that affect things here?\n24. Fig 9C: Doesn't this figure say that SLDS does well during parts of the time period? In which case, how do the results compare during these times?\n25. Fig 11: What new information is D telling us?\n26. Fig 18: Is there a different way to show the same info? This doesn't seem accurate by eye, but how good is it actually?\n27. Fig 22: What does this plot tell us in terms of the differences between classical SLDS and LINOCS SLDS? E.g. what does the number of switches tell us?\n28. Possible discussion items:\n29. For the SLDS model, why does your setup lead to the correct boundary transitions in Fig 4 while SLDS does not? I would assume that for a simple setup this should be doable with SLDS.\n30.The simulated cases here are very small. What would happen for large matricies, which is fairly common for a lot of systems?\n\nTypos or small fixes:\n1. Sec 2.1: typo: Switching, not switched.\n2. HMM should be expanded at least once and cited.\n3. In descriptions of methods (3.1, 3.2...) it could be nice to reference associated figures. This is done for 3.2, but not any others.\n4. Fig 3A: Titles have \"F\" and \"b,\" but shouldn't \"F\" be \"A\" according to the text?\n5. Fig 3A-D: Ylabels are confusing (Values? Quiver? y hat is never used?).\n6. Fig 3F/I are not referenced in the text and the colors aren't labeled.\n7. Fig 3 text \"...increasing noise levels...\" put sigma here.\n8. First eq in 3.4: Write out the regularization term at least once. \n9. \"...The operators are initialized with a regularized 1-step optimization (Fig 3 cyan)...\" Cyan should be orange I think? Also, check which figures are referenced to which others.\n10. \"...even under very high sigma noise levels...\" Sigma should be labeled inside the Gaussian noise definition? Otherwise, this is the first time seeing this.\n11. Sec 4.1 \"...even extreme Gaussian noise levels, sigma = 0.9 (Fig. 3D, blue).\" This should be Fig 3H I think\n12. Sec 4.2 \"Additionally, the offset for each state was set to a random...\" This offset is not shown in Fig 16.\n13. Sec 4.3 \"...the high-order LINOCS-dLDS achieved a structure much closer to the ground truth (Fig. 18).\" There is nothing but LINOCS? So there is no comparison?\n14. Fig 4: Mention that correlation is rho in plots vector\n15. Fig 5: A: Specify in the figure itself for which ones are made by non-LINOCS and LINOCS? B: The colors could be aligned with other figures.\n16. Fig 7: What is lambda? Some hyperparameter? Also should specify somewhere how l0 regularization is implemented. Also, maybe use a log scale or something else for B? With the current scale, it's hard to tell how much better one is compared to another (e.g. the red ones).\n17. Fig 8: units for time. Also, it might make sense to align the patterns between B and C (e.g. LINOCS be solid for both instead of solid for B and dashed for C).\n18. Fig 9C: Maybe do logscale? Can't tell how good or bad the solid line is. Also \"(full curve)\" should be \"solid curve.\"\n19. Fig 10B: Colors are hard to distinguish. Units needed.\n20. Fig 10D: The real and imag part should be on different plots given that the scale is completely different. Should also say lambda in the text for easier readability.\n21. Fig 11: A: wasn't \"w\" called \"lambda\" before? Either way, labels should be consistent. B/E, consider fixing the color scale across plots.\n22. App B was previously said in the main text and doesn't add new information.\n23. Fig 12: Similar comments as fig 4. Also, it may help to order the plots the same way.\n24. Fig 13: Similar comments to fig 4. Some numbers don't fit in their respective boxes. Also, for D, x1 / x2 (which is shown as the circular bit in B), is not the circular bit here; so the axes are not consistent.\n25. Fig 16: \"b\" is not shown. Inconsistent axes (x,y,z notation for C and $x_0$,$x_1$,$x_2$ notation for A). Time is cutoff in B. It would be handy in C to show the temporal evolution, to see the initial conditions and how it evolves. Initial conditions are not specified.\n26. Fig 18: What do the colors mean? And in the text you mention a comparison, but it doesn't look like the comparison is there.\n27. Fig 20: Ground truth doesn't exist, despite being labeled?\n28. Fig 22: White space could be avoided. Also (C), isn't SLDS used, not 1 step?", "broader_impact_concerns": "N/A", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "Pgfilju6lG", "reviewer": "Reviewer_GGPj", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "Strengths\n- Introducing multi-step predictions cost improves long-term reconstruction\n- LINOCS is applicable to various existing dynamical system models\n- LDS with accurate system identification produce reliable and interpretable interaction between state dimensions (element)\n- The real-world neural recoding example demonstrates the practicality of LINOCS.\n\nWeaknesses\n- The scope is limited to linear systems. Though the motivation of LDS is partially interpretability, nonlinear systems would also benefit from multiple-step cost. In addition, the expressive power of LDS is limited to such as fixed points and spiral.\n- The idea of multiple-step prediction was also used by other studies, e.g. Hocker 2017 (10.1109/NER.2017.8008426).\n- Lack of complexity analysis, i.e., how much more time the multiple-step prediction costs the training.\n- Lack of comparison or discussion with other methods aiming for long-term predictions.", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "The paper introduces a learning procedure, LINOCS, for linear dynamical systems to produce accurate long-term predictions and identify interactions in noisy time-series data. LINOCS integrates multiple-step predictions with adaptive weights to accurately capture long-term trends. The authors demonstrate the effectiveness of LINOCS in recovering ground-truth dynamic operators in various linear dynamical system models, including time-invariant linear, switching linear, decomposed linear, and regularized linear systems. They also demonstrated the application on real-world neural recordings, identifying meaningful interactions between multiple regions.", "claims_and_evidence_yn": "Yes", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "- discussion and citations to existing methods with similar ideas\n- presentation of computational cost at different orders\n- compare or discuss with other methods aiming for long-term predictions, e.g. generalized teacher forcing or any other methods to the authors knowledge.", "broader_impact_concerns": "No broader impact concerns.", "additional_comments": "", "extra_scores": {}, "extra_text": {}}]

Author rebuttal:

Meta-review: All reviewers were initially positive about the submission, but had questions of clarity on how LINOCS was different from other methods and how specific technical details were conveyed. The authors provided a detailed and thorough rebuttal and revision, greatly enhancing the quality of their original submission. Due to this lengthy set of revisions, and due to the fact that the reviewers seemed satisfied with them, I am recommending accept as is. 

Note that one of the three reviewers was not able to change their original recommendation, which was made prematurely before the authors provided their responses. However, given that both of the other two reviewers recommend acceptance, and because I believe the authors have sufficiently addressed the comments of that specific reviewer, I am making my decision now (so as to not unnecessarily prolong the process).
Decision: Accept (poster)