Predictability Enables Parallelization of Nonlinear State Space Models
Authors: Xavier Gonzalez
Leo Kozachkov
David M. Zoltowski
Kenneth L. Clarkson
Scott Linderman
Venue: NeurIPS 2025 poster ()
Keywords: dynamical systems
optimization
largest Lyapunov exponent
contraction analysis
RNNs
Newton's method
Parallel algorithms
Scalability
Numerical Stability
state-space models
DEER
Primary area: optimization
TL;DR: Nonlinear systems whose future behavior is not overly sensitive to small perturbations can be efficiently parallelized; whereas unpredictable dynamical systems cannot be efficiently parallelized.

Abstract: The rise of parallel computing hardware has made it increasingly important to understand which nonlinear state space models can be efficiently parallelized. Recent advances like DEER  and DeepPCR recast sequential evaluation as a parallelizable optimization problem, sometimes yielding dramatic speedups. However, the factors governing the difficulty of these optimization problems remained unclear, limiting broader adoption. In this work, we establish a precise relationship between a system's dynamics and the conditioning of its corresponding optimization problem, as measured by its Polyak-Łojasiewicz (PL) constant. We show that the predictability of a system, defined as the degree to which small perturbations in state influence future behavior and quantified by the largest Lyapunov exponent (LLE), impacts the number of optimization steps required for evaluation. For predictable systems, the state trajectory can be computed in at worst $\mathcal{O}((\log T)^2)$ time, where $T$ is the sequence length: a major improvement over the conventional sequential approach. In contrast, chaotic or unpredictable systems exhibit poor conditioning, with the consequence that parallel evaluation converges too slowly to be useful. Importantly, our theoretical analysis shows that predictable systems always yield well-conditioned optimization problems, whereas unpredictable systems lead to severe conditioning degradation. We validate our claims through extensive experiments, providing practical guidance on when nonlinear dynamical systems can be efficiently parallelized. We highlight predictability as a key design principle for parallelizable models.

Predictability Enables Parallelization of Nonlinear State Space Models
Xavier Gonzalez" Stanford University vavier18@stanfrrd. edu
Leo Kozachkov IBM Research
David M. Zoltowski
Stanford University
1eokoz8@broom.edu
tzoltow@stanforz edu
Kenneth L. Clarkson IBM Research
Scott W. Linderman
Stanford University
klclarks@ur imm com
scott www.llidmmantttannoooooooo edu
Abstract
The rise of parallel computing hardware has made it increasingly important understand which nonlinear state space models can be efficiently parallelized. Recent advances like DEER [1] and DeepPCR [2] recast sequential evaluation as a parallelizable optimization problem, sometimes yielding dramatic specdups. However, the factors governing the difficulty of these optimization problems re- mained unclear, limiting broader adoption. In this work, we establish a precise relationship between system's 's dynamics and the conditioning of its correspond- ing optimization problem, as measured by its We show that the predictability of system, defined the degree to which small perturbations in state influence future behavior and quantified by the largest Lya- punov exponent (LLE), impacts the number of optimization steps required for evaluation. For predictable systems, the state trajectory can be computed in at worst ((log 2)2) time, where T is the sequence length: major improvement over the conventional sequential approach. In contrast, chaotic or unpredictable systems exhibit poor conditioning, with the consequence that parallel evaluation converges too slowly to be useful. Importantly, our theoretical analysis shows that predictable systems always yield well-conditioned optimization problems, whereas unpredictable systems lead to severe conditioning degradation. We val- idate our claims through extensive experiments, providing practical guidance on when nonlinear dynamical systems can be efficiently parallelized. We highlight predictability as a key design principle for parallelizable models.
(PL) constant.
Introduction
Parallelization has been central to breakthroughs in deep learning, with GPUs enabling the fast training of large neural networks. In contrast, nonlinear state space models like recurrent neural networks (RNNs) have resisted efficient parallelization on GPUs due to their sequential nature.
Recent work addresses this mismatch by reformulating sequential dynamics into parallelizable op- timization problems. Notably, the DEER/DeepPCR algorithm [1, 2] evaluates nonlinear state space dynamics by minimizing residual-based merit function, facilitating efficient parallel computation
"Equal contribution.
Now at Brown University.
39th Conference on Neural Information Processing Systems (NcurIPS 2025).

---

via the Gauss-Newire method. Gonzalez al. [3] further developed these methods, including quasi-Newton methods and trust- region methods for parallel evaluation of non- linear dynamical systems. These methods evaluate nonlinear dynamical systems by it- eratively linearizing the nonlinear system and evaluating the resulting linear dynamical sys- tem (LDS) with parallel (a.k.a. associative) sean [4, 5]. Each parallel evaluation of an LDS implements optimization step [1-3, 6].
Predictable
Unpredictable
State Space
State Space
Merit Function
Merit Function
St
The usefulness of this optimization- re- formulation depends on two key factors: (a) the computational time per optimization step, and (b) the number of optimization steps re- quired. The computational time per optimiza- tion step is only logarithmic in the sequence length, thanks to its parallel structure. How- ever, the number of steps is governed by the conditioning of the merit function, and that remains poorly understood. In this puper. we charac- terize the merit function's 's conditioning, allowing to draw sharp distinction between systems that are amenable to efficient parallelization via merit function minimization and those that are not (see Figure 1, which is generated from trajectories of an RNN). Geometrically, show that unpre- dictable systems lead to merit functions that have regions extreme flatness, which can lead to very slow convergence.
8t-1
8
Figure I: Predictable nonlinear state space mod- els can be recast as well-conditioned, parallelizable optimization problems.
Drawing from nonlinear dynamical systems theory- particclaal contraction analysis [7] and Lya- punov exponent methods [8]-we formalize the relationship between system predictability and the conditioning of the merit function. Unpredictable systems are dynamical systems whose future behavior is highly sensitive small perturbations.. A common example is chaotic system, like the weather: a butterfly flapping its wings in Tokyo today can lead to a thunderstorm in Manhattan next month [9, 10). By contrast, predictable systems [1 1, 12] those in which small perturbations "forgotten." A familiar example is aviation: a patch of choppy air rarely makes an airplane land at the wrong airport. A more formal definition of (un) pundiiicictyyllityy is given in Definition 1. Our results establish key theoretical principles, make connections between optimization theory and dynamical systems, and demonstrate the practical applicability of parallel computations across wide range of nonlinear state space modeling tasks.
Contributions Outline Our central finding is that predictable systems give rise to well- conditioned merit functions, making them amenable to efficient parallelization. Unpredictable (e.g., chaotic) systems produce poorly conditioned merit functions and are not easily parallelizable.
The paper is organized as follows. Section 2 provides background, with formal definitions of pre- dictable and unpredictable nonlinear state space models. Section 3 presents two key theoretical sults that characterize the conditioning of the merit function, showing that the Polyak-LLojssiwwiez (PL) constant şu of the merit function is controlled by the predictability of the dynamics (Theorem 2), and that the Lipschitz constant of the residual function Jacobian is governed by the nonlinearity of the dynamics (Theorem 3). Section then uses the results about the conditioning of the merit fune- tion to prove results about for Gauss-Newton, with the precise rate scaling with the unpredictability of the problem (Theo- rem 4), and we characterize the basin of quadratic convergence in terms of the predictability and nonlinearity of the underlying dynamics (Theorem 5). In Section we illustrate our results with experiments, and in Section 6 we conclude by summarizing context, implications, limitations, and future directions.
in particular. We prove global linear rates of convergence
'DEER [1] and DeepPCR [2] were concurrent works that both proposed to use the Gauss-Newton method for optimizing nonlinear sum of squares to parallelize sequential processes. In this paper, we therefore use interchangeably.
DEER, DeepPCR, and

---

2 Problem Statement & Background
Notation Throughout the paper, we use to denote the length of sequence and D represent the dimensionality of a nonlinear state space model. Elements in IR, RD or RDx are written using non-bold symbols, while elements R IR. TDxTD are denoted with bold symbols.
Sequential Evaluation vs. Merit Function Optimization We consider the D-dimensional non- linear state space model
8, = (((t--) E ERD.
(I) -, Bu,), where W and are
simple example is an input-driven nonlinear RNN, 8t weight matrices and is the input into the network time t. We want to compute the state trajec- tory, (S1,..., 8T), starting from an initial condition 8p, for a given sequence of functions fi, Ji,.... //r
Systems of the form (1) are widespread across science and engineering, Examples include physics (numerical weather prediction, molecular dynamics), biology (gene regulatory networks, popula- tion dynamics), engineering (control, robotics), and economics (macroeconomic forecasting, asset pricing). In machine learning, sequential operations arise in recurrent neural networks, iterative op- timization, and the sampling pass of diffusion model [2. 13J. Sequential operations even appear probabilistic modeling, se- quential operations arise in Markov Chain Monte Carlo [20]. In all of these cases, the state evolves through nonlinear transformations that capture the system's 's underlying dynamics.
in the problem of evaluating transformer blocks over depth [14-19].
The obvious approach is to sequentially compute the states according to eq. (1), taking T steps. Alternatively, can cast state evaluation advantage of this approach is that it admits parallel computation [1-3]. Depending the properties of the nonlinear state space model, the optimization algorithm, and the available hardware, the latter approach can be significantly faster than sequential evaluation.
optimization problem. While less intuitive, an
We define the residual and corresponding merit function L by stacking the elements s E R of trajectory into a T dimeessiona vector s and considering the vector of temporal differences,
1 C(s) 2 |r(s)]
r(s) vec ([s1 - f(()) ST fr(ST-1))) R TD
where vec(-) denotes the flattening of a sequence of vectors into a single column vector. The true trajectory s* is then obtained by minimizing L(s). Note that the residual is zero only at the true tra- jectory, i.e. when $1, S2, ". ,T satisfy (1) at every time point, so s" is the unique global minimum of C(s).
DeepPCR [2] and DEER [1] minimize the merit function using Gauss-Newton updates. Each update takes the form
(1-1) = (
(3)
where (s) denotes the Jacobian of the residual function, evaluated at the current iterate (()) The Jacobian a TD x TD matrix with D x block bidiagonal structure
ID
0
I
0
0
dr J(s(0) = Os
(s())=
=
:
where 5
Of. Os.-1
(8 (8(2). (4)
0
TD
0
ID/
Due to this block bidiagonal structure, solving (s) sion, which can be done in Olog T') time with parallel scan [1, 3, 5, 22, 23]. Further details are given in Appendix A.
r(s amounts solving linear recur-
This sublinear time complexity per step only useful if the number of optimization steps required to minimize the merit function is small, otherwise it would be more efficient to evaluate the recursion sequentially. Thus, we seek to characterize the conditioning of the merit function - determining
'While minimizing "merit function" is admittedly counterintuitive, we follow Nocedal and Wright [21, eq. 11.35] in this convention.
3

---

when it is well-conditioned and when it is not since this affects the difficulty of finding its min- imum. Equation (4) already offers an important clue. The presence the nonlinear state-space model Jacobians Jr, which measure the local stability and predictability of the nonlinear dynamics, foreshadows our central finding: the system's predictability dictates the conditioning of the merit function.
Predictable Systems: Lyapunov Exponents and Contraction Predictability usually defined through its antonym: wipredictability [9, 10]. In an unpredictable system, the system's 's intrinsic sensitivity amplifies small perturbations and leads massive divergence of trajectories. Predictable systems show the opposite behavior: small perturbations are diminished over time, rather than am- plified. The notion of (un)predictability can be formalized through various routes such as chaos theory [24, 25] and contraction analysis 26].
The definition of predictability comes from the Largest Lyapunov Exponent (LLE) [8, 10]:
Definition 1 (Predictability and Unpredictability). Consider sequence of Jacobians, J1, J2, ... Jr. We define the associated Largest Lyapunov Exponent (LLE) to be
1 T- T
LLE
lim
log J ... Jill)
X,
(5)
where I| an induced operator norm. If X < O, we say that the nonlinear state space model is predictable so- Otherwise, we say is unpredictable.
Suppose wish to evaluate the nonlinear state space model (1) from an initial condition so. but only have access to an approximate measurement s that differs slightly from the true initial state. If the system is unpredictable (A > 0), then the distance between nearby trajectories grows
|| s 11 -61l-
(6)
Letting A denote the maximum acceptable deviation beyond which we consider the prediction to have failed, the time horizon over which the prediction remains reliable scales as
1
A
Time to degrade to A prediction error
A log
50
This relationship highlights key limitation in unpredictable systems: even significant improve- ments in the accuracy of the initial state estimate yield only logarithmic gains in prediction time. The system's 's inherent sensitivity to initial conditions overwhelms any such improvements. Predictable systems, such as contracting systems, have the opposite property: trajectories initially separated by some distance will eventually converge towards another (Figure 1), improving prediction accuracy over time.
3 Conditioning of Merit Function Depends on Predictability of Model
The number of optimization steps required minimize the merit function (2) is impacted by its conditioning, which in our setting is determined by the smallest singular value of the residual func- tion Jacobian. As will see, what determines the smallest singular value of the residual function Jacobian is the stability, or predictability, of the underlying nonlinear state space model (1).
3.1 The Merit Function is PL
To begin, show that the merit function (2) satisfies the Polyak-L.ojasiewicz (PL) condition [27, 28], also known as the gradient dominance condition [29]. A function L(s) is p-PL if it satisfies, for 0,
2 [| 2 P (C(s) L(s"))
(8)
for all s. The largest μ for which eq. (8) holds for all is called the PL constant of C(s). Proposition 1. The merit function L(s) defined in eq. (2) satisfies eq. for
P. inf o2 min (J(s).
(9)

---

Proof. See Appendix B. This result, known in the literature for general sum-of-square, [30], is included here for context and completeness.
Proposition is important as characterizes the flainess the merit function. If μ is very small in a certain region, this indicates that the norm of the gradient can be very small in that region, which can make gradient- based optimization inefficient. Proposition also links immpootnnn for characterizing the conditioning of Joto the geometry of the merit function landscape.
3.2 Merit Function PL Constant is Controlled by the Largest Lyapunov Exponent of Dynamics
As stated earlier, the Largest Lyapunov Exponent is a commonly used way to define the (un)predictability of nonlinear state space model. In order to proceed, we need control more carefully how the product of Jacobian matrices in (5) behaves for finite-time products. We will as- sume that there exists a "burn-in" period where the norm of Jacobian products can transiently differ from the LLE. In particular, we assume that
Vt 1, Ak 2 0, Vs,
b ee S 1| Jell I e
(10)
where 2 and b < 1. The constant quantifies the potential for transient growth or over- shoot-in in the norm of Jacobian products before their long-term behavior emerges, while quantifies the potential for undershoot.
Theorem 2. Assume that the LLE regularity condition (10) holds. Then the PL constant p satisfies
1 -1
(1
a exT -1 S V S min
(11)
Proof. See Appendix C for the full proof and discussion. We provide a brief sketch. Because (min(J) = a/ (J 1), it suffices to control ||- We can write J = I N where N is a nilpotent matrix. Thus, it follows that J- ET N As we discuss further in Appendix C, the matrix powers N are intimately related to the dynamics of the system. The upper bound on | ||2 follows after applying the triangle inequality and the formula for geometric sum. The lower bound follows from considering ||NT- 22-
k=
Theorem 2 is our main result, offering a novel connection between the predictability À of non- linear state space model and the conditioning μ the corresponding merit function, which affects whether the system can be effectively parallelized. If the underlying dynamics are unpredictable (A > 0), then the merit function quickly becomes poorly conditioned with increasing T. because the denominators of both the lower and upper bounds explode due the exponentially growing factor. Predictable dynamics < lead to good conditioning of the optimization problem, and parallel methods based on merit function minimization can be expected to perform well in these cases.
The proof mechanism we have sketched upper and lower bounds ||J- terms of norms of Jacobian products. We only use the assumption in eq. (10) to express those bounds in terms of X. As we discuss at length in Appendix C, we can use different assumptions from eq. (10) to get similar results. Theorem 2 and its proof should be thought of as framework, where different assumptions (which may be more or less relevant in different settings) can be plugged in to yield specific results.
Why Unpredictable Systems have Excessively Flat Merit Functions Theorem 2 demonstrates that the merit function becomes extremely flat for unpredictable systems and long trajectories. This flatness poses a fundamental challenge for any method that seeks to compute state trajectories by minimizing the merit function. We now provide further intuition to explain why unpredictability in the system naturally leads to flat merit landscape.
Suppose that we use an optimizer to minimize the merit function (2) for unpredictable system until it halts with some precision. Let us further assume that the first state of the output of this optimizer following the initial condition is e-close the true first state, |ss - - e. Suppose also that the residuals for all times greater than one are precisely vero-in other words, the optimizer starts with "true" trajectory starting from initial condition $1. Then the overall residual norm is at most
[| = ||s| - f(s)))2 s (l)s s7ll + || f() = |ls - = = 2.
5

---

However, since and 5i are by construction both trajectories of an unpredictable system starting from slightly different initial conditions s₁ and si the distance between them will grow exponen- tially a consequence of eq. (7). By contrast, predictable systems will have errors that shrink exponentially. This shows that changing the initial state 81 by small amount can lead to massive change in the trajectory of an unpredictable system, but a tiny change in the merit function. Geomet- rically, this corresponds to the merit function landscape for unpredictable systems having excessive flatness around the true solution (Figure 1, bottom right panel). Predictable systems do not exhibit such flatness, since small residuals imply small errors. Theorem formalizes this idea.
3.3 Residual function Jacobian Inherits the Lipschitzness of the Nonlinear State Space Model
In addition the parameter 14, which measures the conditioning of the merit function, the difficulty of minimizing the merit function is also influenced by the Lipschitz continuity of its Jacobian J. The following theorem establishes how the Lipschitz continuity of the underlying sequence model induces Lipschitz continuity in J.
Theorem 3. If the dynamics of the underlying nonlinear state space model have L-Lipschitz Jaco- bians, i.e.
Vt>1, 1, s,'' E
[ Je('')l - L|| - s'll.
then the residual function Jacobian also L-Lipschitz, with the same L.
Proof. See Appendix D.
Theorem 3 will be important for the analysis in Section 4, where we consider convergence rates. Because Gauss-Newtor methods rely on iteratively linearizing the dynamics (or equivalently the residual), they converge in single step for linear dynamics L = 0, and converge more quickly if system is close to linear (L is closer to 0).
4 Rates of Convergence for Optimizing the Merit Function
In Section 3, established that the predictability of the nonlinear state space model directly influ- ences the conditioning of the merit function. This insight is critical for analyzing any optimization method used to compute trajectories via minimization of the merit function.
In this section, we apply those results to study the convergence behavior of the Gauss-Newton (DEER) algorithm for the merit function defined in eq. (2). See Appendix A for brief overview of DEER. We derive worst-case bounds on the number of optimization steps required for conver- gence. In addition, we present average-case analysis of DEER that is less conservative than the worst-case bounds and more consistent with empirical observations.
DEER Always Converges Globally at Linear Rate Although DEER is based on the Gauss- Newton method, which generally lacks global convergence guarantees, we prove that DEER always converges globally a linear rate. This result relies on the problem's 's specific hierarchical structure, which ensures that both the residual function Jacobian and its inverse are lower Mock-triangular. In particular we prove the following theorem
Theorem 4. Let the DEER (Gauss-Newtond updates be given by eq. (3), and let s denote the i-th iterate. Let (+)) == (()) s" denote the error at iteration i, and assume the regularity condition in eq. (10). Then the error converges zero al a linear rate:
| ll2 S Xw 3 (
for some constant Xw. independent if and a convergence rale 0 5
Proof. See Appendix E.
Theorem is unexpected since, general, Gauss-Newton methods do not enjoy global convergence. The key caveat of this theorem is the multiplicative factor Xw. which can grow exponentially with the sequence length T. This factor governs the extent of transient error growth before the decay term 3r eventually dominates.

---

Theorem 4 has several useful, practical consequences. First, when the nonlinear state space model sufficiently contracting (A is sufficiently negative), then X Theorem 4 can be made small, implying that in this case DEER converges with little-to-no overshoot (Appendix F).
Theorem 4 also lets us establish key worst-case and average-case bounds the number of steps needed for Gauss-Newton to converge to within a given distance of the solution. In particular, when Xw does not depend on the sequence length T' then Theorem 4 implies Gauss-Newton will only require O (llo ?) total computational time, with one log factor coming from the parallel scan at each optimization step and the other coming from the total number of optimization steps needed. We elaborate on these points in Appendix G.
Size of DEER Basin of Quadratic Convergence It is natural that DEER depends on the Lips- chitzness of J since Neuus- converges in one step for linear problems, where L = 0. In Section 3, we showed that the conditioning of the merit function, as measured by the PL-constant μ4, depends on the stability, or predictability, of the nonlinear dynamics. Thus, the performance of DEER depends on the ratio of the nonlinearity and stability of the underlying nonlinear state space model. Note that once is inside the basin of quadratic convergence, it takes Ollog log(1/))) steps to reach e residual (effectively a constant number of steps).
Theorem 5. Let ş denote the PL-constant of the merit function, which Theorem 2 relates to the LLE A. Let L denote the Lipschitz constant of the Jacobian of the dynamics function Js). Then, 2n/L lower bounds the radius of the basin of quadratic convergence of DEER; that is, if
2p
||rss00|l2 <
(12)
L
then ((0) is inside the basin of quadratic convergence. In terms of the LLE X, it follows that if
2
e^11
2
[|(s <
a²L
XT
then (00) is inside the basin of quadratic convergence.
Proof. See Appendix H. We make no claim about the originality of lower bounding the size of the basin of quadratic convergence in Gauss-Newton the convergence analysis of Newton's 's method in Section 9.5.3 of Boyd and Vandenberghe [31]. Our contribution is we highlight the elegant way the predictability À and nonlinearity L of dynamical system influence an important feature of its merit function's 's landscape.
In fact, our proof of Theorem 5 closely follows
Experiments
We conduct experiments to support the theory developed above. demonstrating that predictability enables parallelization of nonlinear SSMs. To illustrate this point, we use Gauss-Newton opti- mization (aka DEER [1, 2)). We provide more experimental details in Appendix K. Our code is com wiinddemanaabprrdeieeete
at /attts::/gittee http::////...aaaooooo/o//aaaaaaa sewrratcttailltttyeblleesesee wnwales_ppaaaleelle.
The Convergence Rate Exhibits a Threshold between Predictable and Chaotic Dynamies Theorem 2 predicts a sharp phase transition in the conditioning of the merit function at = 0, which should be reflected in the number of optimization steps required for convergence. To em- pirically validate this prediction, we vary both the LLE and sequence length T within a parametric family of recurrent neural networks (RNNs). and measure the number of steps DEER takes con- verge. We generate mean-field RNNs following Engelken et al. [32], scaling standard normal weight matrices by a single parameter that controls their variance and therefore the expected LLE. In Fig- ure 2, we observe a striking correspondence between the conditioning of the optimization problem (represented by log P. where ğ is the lower bound for μ from Theorem 2) and the number of steps DEER takes to converge. This relationship holds across the range of LLEs, X, and sequence lengths, T. There a rapid threshold phenomenon around X = O, which divides predictable from unpredictable dynamics, precisely as expected from Theorem 2. As we discuss in Appendix K. the correspondence between log ji and the number of optimization steps needed for convergence can be explained by DEER iterates approaching the basin of quadratic convergence with linear rate.
7

---

Theory log(())
Experiment (steps to convergence)
Experiment: 7 1000
10000
2800
10000
steps to converge
1000
2100
7500
5000
1400
5000
500
2500
0
0
0.6-0.4 11-00.44... 0.2 0.0
0.6 06...-0..2 0.0
-0.6-0. 1..-0.4-0..0 0.0
LLE (A)
LLE (A)
LLE (A)
Figure 2: Threshold phenomenon in DEER convergence based on system predictability. In family of RNNs, DEER has fast convergence for predictable systems and prohibitively slow conver- gence for chaotic systems, Left (Theory): We depict Theorem 2, illustrating how the conditioning of the optimization problem degrades as T and the LLE (A) increase. Center (Experiment): We vary À across the family of RNNs, and observe a striking concordance in the number of DEER optimization steps empirically needed for convergence with our theoretical characterization of the conditioning of the optimization problem. Right: For 20 seeds, each with 50 different values of X, we plot the relationship between and the number of DEER steps needed for convergence for sequence length 7 = 1000 (gray line in left and center panels). We observe sharp increase in the number of optimization steps at precisely transition between predictability and unpredictability.
In Appendix K.3, we provide additional experiments in this setting. We parallelize the sequential rollout with other optimizers like quasi-Newtor and gradient descent, and observe that the number of steps these optimizers take to converge also scales with the LLE. We also record wallclock times on an H100, and observe that DEER is faster than sequential by order of magnitude in predictable settings, but slower by an order of magnitude in unpredictable settings.
DEER can converge quickly for predictable trajectories passing through unpredictable regions DEER may still converge quickly even if the system is unpredictable in certain regions. As long as the system is predictable on average, as indicated by a negative LLE, DEER can still converge quickly. This phenomenon is why we framed Theorem in terms of the LLE À and burn-in constants a, as opposed to weaker result that assumes the system Jacobians have singular values less than one over the entire state space (see discussion of condition (10) vs. condition (22) in Appendix C).
To illustrate, we apply DEER to Langevin dynamies ure 3 for D = 2). The dynamics are stable within each well but unstable in the region between them. Despite this local instability, the system's overall behavior governed by time spent in the wells, resulting in negative LLE and sublinear growth in DEER's convergence steps with sequence length (Figure 3, right subplot). Additional details and discussion are Appendix K.4.
two-well potential (visualized in Fig-
Notably, prior works such as Lim et al. [1] and Gonzalez et al. [3] initialized optimization from (0)) 0, which lies entirely in the unstable region. Thus, our theoretical insights into predictability and parallelizability suggest practical improvements for initialization.
Application: Chaotic Observers Finally, demonstrate practical application of our theory in the efficient parallelization of chaotic observers. Observers are commonly used to reconstruct the full state of system from partial measurements [33, 34]. On nine chaotic flows from the dysts benchmark dataset [35], Table shows that while DEER converges prohibitively slowly on chaotic systems, it converges rapidly stable observers of these systems, in accordance with our theory that predictability implies parallelizabilit.. For more details, see Appendix K.5.
6 Discussion
Recent work demonstrated that parallel computing hardware like GPUs can be used to rapidly com- pute state trajectories of nonlinear state space models (nSSMs) by recasting the trajectory as the solution to an optimization problem. In this work, provide the first precise characterization of

---

Two-well state space
Y coordinate over time
Steps needed convergencee
1.050
Spectral Norm of J(s)
Number DEER steps
20
1.025
15
0.975
0.950
0
2500 5000
7500
10000
5000
10000
Time
Sequence Length T
Figure 3: DEER converges quickly for Langevin dynamics in two-well potential. (Left) An illustration of the two-well potential state space in 2. We superimpose a contour plot of the potential color scheme showing the spectral norm of the dynamics Jacobian (blue indicates stability, red instability). (Center) A trace plot for the y-coordinate. The LLE the system is -0.0145. (Right) We observe that this system, which has negative LLE, enjoys sublinear scaling in the sequence length T in the number of DEER iterations needed to converge. We plot the median number of DEER steps convergence over 20 random seeds.
Table 1:
Comparison of system and observer LLEs and number of DEER steps for T 30, 000 and Euler discretization step size At 0.01.
LLE (System)
LLE (Observer)
DEER Steps (System)
DEER Steps
System
(Observer)
ABC
0.16
-0.08
4243
Chua's Circuit
0.02
-1.37
Kawczynski-Strizhan
0.01
-3.08
29396
Lorenz
1.02
-6.28
30000
Nosé Hoover Thermostat
0,02
-0.13
29765
Rössler
0.01
-0.07
29288
Sprott
0.20
-0.39
29486
Thomas
0.01
-3.07
12747
Vallis EI Niño
0.58
-2.48
30000
the optimization problem's 's inherent difficulty, which determines if parallelization will be faster in practice than sequential evaluation. We show that the conditioning of the optimization problem is governed by the predictability of the underlying dynamics. We translate this insight into worst- case performance guarantees for specific optimizers, including Gauss-Newton (DEER). Our main takeaway is: Predictable dynamics yield well-conditioned merit functions, enabling rapid conver- gence. Unpredictable dynamics produce flat or ill-conditioned merit landscapes, resulting in slow convergence or numerical failure.
While Lim et al. [1] and Danieli et al. [2] introduced parallel Newton methods, they did prove their global convergence. Gonzalez et al. [3] proved global convergence, though only with worst-case bounds of optimization steps. These prior works did not address the relationship between system dynamics and conditioning, establish global linear convergence rates.
Related Work
Global convergence rates for Gauss-Newton are rare, despite the breadth of optimization literature [21, 31, 36, 37]. Theorem 4 establishes global convergence with linear rate for Gauss-Newton by leveraging our specific problem structure, though similar results have existed for local linear convergence [38], most famously the DeotonnKKattrovine theorem [39].
Fifty years ago, Hyafil and Kung [40] and Kung [41] showed that linear recursions enjoy speedups from parallel processors while nonlinear recursions of rational functions with degree larger than one cannot. These prescient works set the stage for our more general findings, which explicitly link the dynamical properties of the recursion to its parallelizability. dtrallel-in-time methods for continuous

---

systems also have long history [42-44]. with Chartier and Philippe [45] showing that dissipative systems can be parallelized using multiple shooting. Furthermore. Danieli and MacLachlan [46] and De Sterck et al. [47] study the CFI. number for determining the usefulness of multigrid systems. Connecting this work with our paper is an interesting direction for future research.
More recently, several works have parallelized diffusion models via fixed-point iteration, including worst-case guarantees of I' steps [13, 48, 49] as well as polylogarithmic rates in T [50, 51J. Lu et al. [52] develops quasi-Newton methods for sampling from diffusion models and, like us, shows two-phase model of linear followed by quadratic convergence. Crucially, prior work has not focused on the merit function, which we can define for any discrete-time dynamical system and optimizer.
To our knowledge, no prior work connects the LLE of dynamical system to the conditioning of the corresponding optimization landscape, as established in Theorem 2. In particular, we showed that systems with high unpredictability yield poorly conditioned (i.e., flat) merit functions, linking dynamical instability to optimization difficulty in a geometrically appealing way.
The centrality of parallel sequence modeling architectures like transformers [53]. deep SSMs [54, 23, 55], and linear RNNs [56] in modern machine learning underscores the need for our theo- retical work. Merrill et al. [57] explored the question of parallelizability through the lens of circuit complexity, analyzing when deep learning models can solve structured tasks in constant depth. Their focus complements ours, and suggests an opportunity for synthesis in future work [58].
Implications Our work unlocks three key implications for nonlinear state space models.
First, it provides a principled way to determine, a priori, whether optimization- based parallelization of given model is practical. In many robotic control systems, particularly ones that are strongly dissipative, this insight can enable orders-of -faanitu-- speed-ups on GPUs [59-67].
For example, the concurrent work of Zoltowski et al. [20] developed and leveraged quasi-Newton methods to parallelize Markov Chain Monte Carlo over the sequence length, attaining order-of- magnitude speed-ups. These speed-ups occurred because the quasi-Newton methods converged quickly in the settings considered. Suggestively, MCMC chains are contractive in many settings [68-70]. A precise characterization of what makes MCMC algorithm and target distribution predictable would provide useful guidance for when one should aim to parallelize MCMC over the sequence length. Providing precise theoretical justification for parallelizing MCMC over the sequence length is an exciting avenue for future work.
Second, our results impact architecture design. When constructing nonlinear dynamical systems in machine learning-such as novel RNNs-parallelization benefits are maximized when the system is made predictable. Given the large body of work on training stable RNNs [71-82], many effective techniques already exist for enforcing stability or predictability during training. A common approach is to parameterize the model's weights so that the model is always stable (see Appendix I).
Notably, the concurrent work of Farsang et al. [82] and Danicli al. [83] develop nonlinear SSMs and train them with DEER, with Danieli et al. [83] scaling to very strong performance as a 7B pa- rameter language model. Both highlight the fast convergence of DEER, which result of the contractivity of their architectures: Farsang et al. [82] parameterizes their LreSSM to be contractive, while Danieli et al. [83] clip the norms of their weight matrices. Ensuring a negative largest Lya- punov exponent through parameterization guarantees parallelizability for the entire training process, enabling faster and more scalable learning. Our contribution provides theoretical foundation for why stability is essential in designing efficiently parallelizable nonlinear SSMs.
Finally, our results have implications for the interpretation of stable nSSMs. Because each Gauss- Newton step in DEER is a linear dynamical system (LDS), and because we prove in Theorem 4 that DEER converges in (log T') steps for stable nSSM, we can interpret a stable nSSM as being equivalent to "stack" of Oloo T) LDSs coupled by nonlinearities (cf. Appendix J).
Limitations and Future Work While this work focuses establishing the fundamental concepts and theoretical foundations, several practical considerations arise for scaling to large systems. No- tably, DEER incurs significant memory footprint. While this issue can be alleviated through quasi- Newton methods [3, 20]. these approaches require more optimization steps to converge. Studying quasi-Newton methods with our theory could provide new insight into the efficacy of these methods.
Overall, the theoretical tools developed here have immediate implications for parallelizing nonlinear systems, and they open several exciting avenues for future work.
10

---

Acknowledgments
We thank members of the Linderman Lab for helpful feedback, particularly Noah Cowan, Henry Smith, Skyler Wu, and Etaash Katiyar. We also thank Federico Danieli, Will Merrill, Julien Siems, and Riccardo Grazzi for helpful discussions. We thank the anonymous NeurIPS reviewers whose feedback improved this paper.
X.G. acknowledges support from the Walter Byers Graduate Scholarship from the NCAA. L.K. was a Goldstine Fellow at IBM Research while conducting this research. S. W.W.L. was supported by fellowships from the Simons Collaboration on the Global Brain, the Alfred P. Sloan Foundation, and the McKnight Foundation.
We thank Stanford University and the Stanford Research Computing Center for providing computa- tional resources and support. Additional computations were performed on Marlowe [84], Stanford University's GPU-based Computational Instrument, supported by Stanford Data Science and Stan- ford Research Computing.
The authors have no competing interests to declare.
References
[1] Yi Heng Lim, Qi Zhu, Joshua Selfridge, and Muhammad Firmansyah Kasim. Parallelizing non-linear sequential models over the sequence length. The Twelfth International Confer- ence Learning Representations, 2024.
[2] Federico Danieli, Miguel Sarabia, Xavier Suau Cuadros, Pau Rodriguez, and Luca Zappella. DeepPCR: Parallelizing sequential operations in neural networks, Advances in Neural Infor- mation Processing Systems, 36:47598-47625, 2023.
[3] Xavier Gonzalez, Andrew Warrington, Jimmy T.H. Smith, and Scott W. Linderman. Towards Scalable and Stable Parallelization of Nonlinear RNNs. In Advances in Neural Information Processing Systems (NeurIPS), 2024.
[4] Harold S. Stone. An efficient parallel algorithm for the solution of tridiagonal linear system of equations. Journal of the ACM, 20(1):27-38, 1973.
[5] Guy E. Blelloch. Prefix sums and their applications. Technical Report CMU-CS-90-199 School of Computer Science, Carnegie Mellon University, November 1990.
[6] Xavier Gonzalez, E. Kelly Buchanan, Hyun Dong Lee, Jerry Weihong Liu, Ke Alexander Wang, David M. Zoltowski, Christopher Ré, and Scott W. Linderman. A unifying framework for parallelizing sequential models with linear dynamical systems. Transaction on Machine Learning Research (TMLR), 2026.
[7] Winfried Lohmiller and Jean-Jacques E Slotine. On contraction analysis for non-linear sys- tems. Automatica, 34(6):683-69)) 1998.
[8] Arkady Pikovsky and Antonio Politi. Lyapunov exponents: tool to explore complex dynam- ics. Cambridge University Press, 2016.
[9] Michael James Lighthill. The recently recognized failure of predictability in Newtonian dy- namics. Proceedings the Royal Society of London. A. Mathematical and Physical Sciences, 407(1833)335500 1986.
[10] Steven H Strogatz. Nonlinear dynamics and chaos with student solutions manual: With ap- plications to physics, biology, chemistry, and engineering. CRC press, 2018.
[11] Philip Duncan Thompson. Uncertainty of initial state factor in the predictability of large scale atmospheric flow patterns. Tellus. 93)):755-255, 1957.
[12] Edward N Lorenz. Predictability: problem partly solved. In Proceedings of the Seminar Predictabilityy volume 1, pages 1-18. ECMWF Reading, UK, 1996.

---

[13] Andy Shih, Suneel Belkhale, Stefano Ermon, Dorsa Sadigh, and Nima Anari. Parallel sam- pling of diffusion models. 37th Conference on Neural Information Processing Systems, 2023. 37th Conference on Neural Information Processing Systems.
[14] Mostafa Dehghani, Stephan Gouws, Oriol Vinyals, Jakob Uszkoreit, and Lukasz Kaiser. Uni- versal transformers. In International Conference on Learning Representations (ICLR), 2019.
URL https:// arxiv. org/ abs /1807 03819.
[15] Mark Schöne, Babak Rahmani, Heiner Kremer, Fabian Falek, Hitesh Ballani, and Jannes Gladrow. Implicit language models are mas: Balancing parallelization and expressivity. In ICML, 2025. doi: (048550aaXxxv22022.7777 URL https://arxiv. www.aaiwiw...iaaaaaaaeaaaaeeeeee herggabs/25000 07827.
[16] Jonas Geiping, Sean McLeish, Neel Jain, John Kirchenbauer, Siddharth Singh, Brian R. Bartoldson, Bhavya Kailkhura, Abhinav Bhatele, and Tom Goldstein. Scaling up test- time compute with latent reasoning: A recurrent depth approach, 2025. URL https: 05171.
[17] Ramón Calvo-González, Daniele Paliotta, Matteo Pagliardini, Martin Jaggi, and François Fleuret. Leveraging the true depth of llms, 2025. URL atttps///arie oog.oog/aaao 2b22222.. 02790. Introduces Layer Parallelism for parallelizing adjacent Transformer layers.
[18] Guan Wang, Jin Li, Yuhao Sun, Xing Chen, Changling Liu, Yue Wu, Meng Lu, Sen Song, and Yasin Abbasi Yadkori. Hierarchical reasoning model, 2025. URL https: //axxiv. org/ ht5/abs/25066 21734.
[19] Alexia Jolicoeur Less is more: Recursive reasoning with tiny networks, 2025.
[20] David M. Zoltowski, Skyler Wu, Xavier Gonzalez, Leo Kozachkov, and Scott W. Linder- man. Parallelizing MCMC Across the Sequence Length. In Advances in Neural Information Processing Systems (NeurIPS), 2025.
[21] Jorge Nocedal and Stephen J. Wright. Numerical Optimization. Springer, edition, 2006.
[22] Eric Martin and Chris Cundy. Parallelizing linear recurrent neural nets over sequence length. In International Conference on Learning Representations, 2018.
[23] Jimmy T.H. Smith, Andrew Warrington, and Scott W. Linderman. Simplified state space lay- ers for sequence modeling. International Conference on Learning Representations (ICLR), 2023.
[24] James Gleick. Chaos: Making a new science. Penguin, 2008
[25] Heinz Georg Schuster and Wolfram Just. Deterministic chaos: an introduction. John Wiley & Sons, 2006.
[26] F. Bullo. Contraction Theory for Dynamical Systems. Kindle Direct Publishing, 1.2 edition, 2024. ISBN 979-8836646806.
[27] Boris Polyak. Gradient methods for the minimisation of functionals. Zh. Vychisl. Mat. Mat. Fiz, 3(4))644-663, 1963.
[28] Hamed Karimi, Julie Nutini, and Mark Schmidt. Linear convergence of gradient and proximal-gradient methods under the Polyak-t.ojasiewicr, condition. In Machine Learning and Knowledge Discovery in Databases: European Conference, ECML PKDD 2016, Riva del Garda, Italy, September 19-23, 2016, Proceedings, Part 16, pages 795-811. Springer, 2016.
[29] Maryam Fazel, Rong Ge, Sham Kakade, and Mehran Mesbahi. Global convergence of policy gradient methods for the linear quadratic regulator. In International Conference on Machine Learning (ICML), 2018.
[30] Yurii Nesterov and B. T. Polyak. Cubic regularization of Newton method and its global performance. Mathematical Programming. Series A, 178(1): 2006.
12

---

[31] Stephen Boyd and Lieven Vandenberghe. Conver Optimization. Cambridge University Press, Cambridge, UK, 2004. ISBN 9780521833783
[32] Rainer Engelken, Fred Wolf, and Larry F Abbott. Lyapunov spectra of chaotic recurrent neural networks. Physical Review Reseurch, 5(4):043044, 2023.
[33] David G Luenberger. Introduction dynamic systems: theory, models, and applications. John Wiley & Sons, 1979.
[34] Dan Simon. Optimal state estimation: Kalman, H infinity, and nonlinear approaches. John Wiley & Sons, 2006.
[35] William Gilpin. Chaos as an interpretable benchmark for forecasting and data- In Joaquin Vanschoren and Sai-Kit Yeung, editors, Pro-
driven modelling. ceedings of the Neural Information Processing Systems Track on Datasets and Benchmarks 1 (NeurIPS Datasets and Benchmarks 2021). December 2021, virtual, 2021. URL ttaap://ddatasete-eenchi Srraaaaceeeeeccccccccceeeeeeeeee newggip.......... cc/paper/ 2021/hash/ h htta:///eeccaaaaaooa//e/eeeaaaaa -eee---bsttcttruuteaaee
html. htmml.
[36] Jim Zhao, Sidak Pal Singh, and Aurelien Lucchi. Theoretical characterisation of the gauss-newton conditioning in neural networks. In Neural Information Processing Systems (NeurlPS), 2024.
[37] Yurii Nesterov. Lectures on Convex Optimization, volume 137 of Springer Optimization and Its Applications. Springer, 2nd edition, 2018. ISBN 978-3-319-91577-4. doi: 10.1007/ 978-3-319-91578-1
[38] James Ortega and Werner C Rheinboldt. Iterative Solution of Nonlinear Equations in Severul Variables. Academic Press, New York and London, 1970. Republished by SIAM in 2000.
[39] L. V. Kantorovich. Functional analysis and applied mathematics. Uspekhi Matematicheskiich Nauk, 366):89-185, 1948. In Russian. English translation in: NBS Report 1509, Washington D.C., 1952.
[40] L Hyafil and HT Kung. Bounds on the speed-up of parallel evaluation of recurrences. Camegie Mellon University, Department of Computer Science, 1975.
[41] HT Kung. New algorithms and lower bounds for the parallel evaluation of certain rational expressions and recurrences. Journal of the ACM (JACM), 23(2):252-261, 1976.
[42] J. Nievergelt. Parallel methods for integrating ordinary differrnnial equations. Communica- tions of the ACM, 7(12):331733, December 1964.
[43] Martin J Gander. 50 years of time parallel time integration. In Multiple Shooting and Time Domain Decomposition Methods: MuS-TDD, Heidelberg, May 6-8, 2013, pages 69-1133 Springer, 2015.
[44] Benjamin W Ong and Jacob B Schroder. Applications of time parallelization. and Visualization Science, 13:1-15, 2020.
Computing
[45] Philippe Chartier and Bernard Philippe. Eine parallele "shooting" technik zur lösung dissipa- tiver gewöhnlicher differentialgleichungen Computing, 51:209-236. 1993.
[46] Federico Danieli and Scott Maclach Multigrid reduction in time for non-linear hyperbolic equations. arXiv preprint 2eXiv:21040944/ 2021.
[47] Hans De Sterek, Stephanie Friedhoff, Oliver A Krzysik, and Scott P MacLachlan. Multigrid Reduction- In-mee Convergence for Advection Problems: A Fourier Analysis Perspective. Numerical Linear Algebra with Applications, 32(1):e2593, 2025.
[48] Zhiwei Tang, Jiasheng Tang, Hao Luo, Fan Wang, and Tsung-Hui Chang. Accelerating paral- lel sampling of diffusion models. In Forty-first International Conference on Machine Learn-
ing, 2024.
13

---

[49] Nikil Selvam, Amil Merchant, and Stefano Ermon. Self-Refinin Diffusion Samplers: En- abling Parallelization via Parareal Iterations. In Advances in Neural Information Processing Systems, volume 37, pages 5429-5453, 2024.
[50] Nima Anari, Sinho Chewi, and Thuy-Duong Vuong. Fast parallel sampling under isoperime- try. In The Thirty Seventh Annual Conference on Learning Theory, pages 161-185. PMLR, 2024.
[51] Haoxuan Chen, Yinuo Ren, Lexing Ying, and Grant M. Rotskoff. Accelerating diffusion mod- els with parallel sampling: Inference sub-linear time complexity. In Neural Information Processing Systems (NeurIPS), 2024.
[52] Jianrong Lu, Zhiyu Zhu, and Junhui Hou. Parasolver: A hierarchical parallel integral solver for diffusion models. In International Conference on Learning Representations (ICLR), 2025.
[53] Ashish Vaswani, Noam Shazeer, Niki Parmar, Jakob Uszkoreit, Llion Jones, Aidan Gomez, tukasz Kaiser, and Illia Polosukhin. Attention is All You Need. In Advances in Neural Information Processing Systems (NeurIPS), 2017.
[54] Albert Gu, Karan Goel, and Christopher Ré. Efficiently modeling long sequences with struc- tured state spaces. In International Conference on Leurning Representations (ICLR), 2021.
[55] Albert Gu and Tri Dao. Mamba: Linear-Time Sequence Modeling with Selective State Spaces. In Conference on Language Modeling (COLM), 2024.
[56] Songlin Yang. Bailin Wang, Yu Zhang, Yikang Shen, and Yoon Kim. Parallelizing linear transformers with the delta rule over sequence length. In Proceedings of NeurIPS, 2024.
[57] William Merrill, Jackson Petty, and Ashish Sabharwal. The illusion of state in state-space models. Forty-first International Conference on Machine Learning, 2024.
[58] Yuxi Liu, Konpat Preechakul, Kananart Kuwaranancharoen, and Yutong Bai. The serial scal- ing hypothesis. In International Conference on Learning Representations (ICLR), 2026.
[59] J Zico Kolter and Gaurav Manek. Learning stable deep dynamics models. Advances in neural information processing systems, 32, 2019.
[60] Hadi Beik-Mohammadi, Søren Hauberg, Georgios Arvanitidis, Nadia Figueroa, Gerhard Neumann, and Leonel Rozo. Neural contractive dynamical systems. arXiv preprint arXiv:2401.09352, 2024.
[61] Sean Jaffe, Alexander Davydov, Deniz Lapsekili, Ambuj K Singh, and Francesco Bullo. Learning neural contracting dynamics: Extended linearization and global guarantees. Ad- vances in Neural Information Processing Systems, 37:66204-66225, 2024.
[62] Fletcher Fan, Bowen Yi, David Rye, Guodong Shi, and Ian R Manchester. Learning sta- ble koopman embeddings. In 2022 American Control Conference (ACC), pages 2742-2747 IEEE, 2022.
[63] Vikas Sindhwani, Stephen Tu, and Mohi Khansari. Learning contracting vector fields for stable imitation learning. arXiv preprint arXiv: 1804.04878, 2018.
[64] Dawei Sun, Susmit Jha, and Chuchu Fan. Learning certified control using contraction metric. conference on Robot Leurning, pages 1519-1539. PMLR, 2021.
[65] Hiroyasu Tsukamoto, Soon-Jo Chung, and Jean-Jacques E Slotine. Contraction theory for nonlinear stability analysis and learning-based control: A tutorial overview. Annual Reviews in Control, 52: 52:135-169, 2021.
[66] Max Revay, Ruigang Wang, and Ian R Manchester. Recurrent equilibrium networks: Flexible dynamic models with guaranteed stability and robustness. IEEE Transactions Automatic Control, 2855528555 2023.
[67] Alexander Davydov and Francesco Bullo. Perspectives on contractivity in control, optimiza- tion, and learning. IEFE Control Systems Letters, 2024.
14

---

[68] Nawaf Bou-Rabee, Andreas Eberle, and Raphael Zimmer. Coupling and con-
vergence for Hamiltonian Monte Carlo. The Annals of Applied Probabil-
iy. 209 June 2020. doi: 10.1214/19-AAP5528 URL https:
//projecteuclid.
issue-
1214/
A9-528.
[69] Oren Mangoubi and Aaron Smith. Mixing of Hamiltonian Monte Carlo on strongly log-concave distributions: Continuous dynamics. The Annals of Applied Probability. 31(5):
2019-2045, October 2021. doi: 10. 10.1214/20-AA11440 URL https://projecteucclll org/ wwwpurunnnnnaaaaalaeoaaeeaaaaaae
of Hww.-Haaamtttna-ttaaaaaa
10. 122144/0--AP1644
[70] Persi Diaconis and David Freedman. Iterated random functions. SIAM Review, 41(1):45-76, 1999.
[71] Sepp Hochreiter. Untersuchungen dynamischen neuronalen Netzen. Diploma thesis, Tech- nische Universität München, Munich, Germany, 1991.
[72] Sepp Hochreiter and Jürgen Schmidhuber. Long short-term memory. Neural computation, (8):1735-1780) 1997.
[73] Ilya Sutskever. Training Recurrent Neural Networks. PhD thesis, University of Toronto, Graduate Department of Computer Science, Toronto, Canada, 2013. PhD thesis.
[74] John Miller and Moritz Hardt. Stable recurrent models. In International Conference Learning Representations, 2019.
[75] N Benjamin Erichson, Omri Azencot, Alejandro Queiruga, Liam Hodgkinson, and Michael W Mahoney. Lipschitz recurrent neural networks. arXiv preprint arXiv:2006. 12070, 2020.
[76] Leo Kozachkov, Michaela Ennis, and Jean-Jacques Slotine, Rons mns: Recursive con- struction of stable assemblies of recurrent neural networks. In Advances in Neural Informa- tion Processing Systems (NeurIPS), 2022.
[77] Karan Goel, Albert Gu, Chris Donahue, and Christopher Re. It's raw! audio generation with state-space models. In International Conference on Machine Learning, 2022.
[78] Dmitry Krotov. new frontier for Hopfield networks. Nature Reviews Physics, 5(7):366-367, 2023.
[79] Rainer Engelken. Gradient flossing: Improving gradient descent through dynamic control of jacobians. Advances in Neural Information Processing Systems (NeurIPS), 2023.
[80] Antonio Orvieto, Samuel Smith, Albert Gu, Anushan Fernando, Caglar Gulcehre, Razvan Pascanu, and Soham De. Resurrecting recurrent neural networks for long sequences. In International Conference on Machine Learning, pages 26670-26698. PMLR, 2023.
[81] Nicolas Zucchet and Antonio Orvieto. Recurrent neural networks: vanishing and exploding gradients are not the end of the story. In Neural Information Processing Systems (NeurIPS), 2024.
[82] Mónika Farsang, Ramin Hasani, Daniela Rus, and Radu Grosu. Scaling Up SiquiddResiistnce. Siquid-Capacitanee Networks for Efficient Sequence Modeling. In Advances in Neural In- formation Processing Systems (NeurIPS), 2025.
[83] Federico Danicli, Pau Rodriguez, Miguel Sarabia, Xavier Suau, and Luca Zappella. Paramn: Unlocking parallel training of nonlinear mns for large language models. In International Conference on Learning Representations (ICLR), 2026.
[84] Craig Kapfer, Kurt Stine, Balasubramanian Narasimhan, Christopher Mentzel, and Em- manuel Candès. Marlowe: Stanford's GPU-based Computational Instrument, Zenodo, 2025. Version O.
15

---

[85] Xavier Gonzalez. Parallelizing Nonlinear RNNs with the Ungulates: DEER and ELK, De- cember 2024. URL https://lindermanlaa github. io/hackatthnss/..... Linderman
Lab Blog.
[86] Sinho Chewi and Austin J. Stromme. The ballistic limit of the log-soboley constant equals the polyak-tojasiewiez constant. Annales de l'Institut Henri Poincaré (B) Probabilités Statistiques, 2025. URL httph:tt////e //txxv/axivvv///a///// org/ org/abs/// abbs/411... 11415.
[87] Tri Dao and Albert Gu. Transformers are SSMs: Generalized models and efficient algorithms through structured state space duality. In International Conference Machine Learning (ICML), 2024.
[88] Chaoyue Liu, Libin Zhu, and Mikhail Belkin. Loss landscapes and optimization in over- parameterized non-linear systems and neural networks. Applied and Computational Har- monic Analysis, 59:85-116, 2022.
[89] Mark Konstantinovich Gavurin. Nonlinear functional equations and continuous analogues of iteration methods. Izvestiya Vysshikh Uchebnykh Zavedenii. Matemalika, pages 18-31, 1958.
[90] Diederik P. Kingma and Jimmy Ba. Adam: A method for stochastic optimization. In 3rd International Conference on Learning Representations (ICLR), 2015, URL https: / a/raiiiv. org/ aba/ 1412. 6980.
[91] Vineet Gupta, Tomer Koren, and Yoram Singer. Shampoo: Preconditioned stochastic tensor optimization. In Jennifer Dy and Andreas Krause, editors, Proceedings of the 35th Inter- national Conference on Machine Learning, volume 80 of Proceedings of Machine Learning Research, pages 1842-1850. PMLR, 10-15 Jul 2018. URL https: https://pproceedings .ml.. mlr. press/v80/ gupta18a. html.
[92] Nikhil Vyas, Depen Morwani, Rosie Zhao, Itai Shapira. David Brandfonbrener, Lucas Janson, and Sham M. Kakade. SOAP: Improving and Stabilizing Shampoo using Adam for Language Modeling. In International Conference on Learning Representations (ICLR), 2025.
[93] Paul Langevin. On the Theory of Brownian Motion. American Journal of Physics, 65 (1 (119:177,10881 1908. English translation, introduced by D. Lemons and translated by A. Gythiel. Original: C. R. Acad. Sei. 146, 530-533 (1908).
[94] Roy Friedman. A simplified overview of Langevin dynamics. Blog post, https:// fr friedmanroy. wewedmanoo..gihuu..iie.. http:///ooo//////aaaa/e///////// 2022.
[95] Desmond J. Higham. An algorithmic introduction to numerical simulation of stochastic dif- ferential equations. SIAM Review, 43(3):525-546, 2001.
[96] Matthew MacKay. Stats 305c lecture 17: Stochastic differential equations. Course slides, STATS 305C: Applied Statistics III, Stanford University, 2022.
[97] Peter Holderrieth. Langevin dynamies: An introduction for machine learn- http://www.pterhhoerrrreeee com/blog/2025
ing engineers, 2023. URL. worooo-foo---aaaaae iatroductttooo-------hhhheeee -egg Blog Langevin-Dynamics An- inntroductiinnnnooo post.
[98] William Gilpin. Chaos an interpretable benchmark for forecasting and data-driven mod- elling. In Thirty-fifth Conference on Neural Information Processing Systems Datasets and Benchmarks Track (Round 2), 2021.
[99] Louis M Pecora and Thomas L Carroll. Synchronization in chaotic systems. Physical review letters, 64(8):821, 1990.
[100] Ali Zemouche and Mohamed Boutayeb. Observer design for Lipschitz nonlinear systems: the discrete-time case. IEEE Transactions on Circuits and Systems II: Express Briefs, 53(8): 777-781, 2006.
[101] Bertrand Legras and Robert Vautard. A guide to Liapunov vectors. In Proceedings of the ECMWF Seminar on Predictability, 4-8 September 1995, volume 1, pages 143-156, Reading, UK, 1996. FCMWE.
16

---

Appendix
Brief Overview DEER/DeePPCR
18
B Merit Function is PL
18
Merit Function PL Constant Controlled by Largest Lyapunov Exponent of Model
19
C.1 A Generalized Proof that the Largest Lyapunov Exponent Controls the PL Constant
24
DEER Merit Function Inherits Lipschitzness of Dynamics Jacobians
27
E
DEER Always Converges Linearly
27
DEER Converges Globally with Small Overshoot for Sufficiently Strongly Contracting Systems
29
Alternative Descent Techniques Worst/Average Complexity
H Proof of Size of Basin of Quadratic Convergence
31
Parameterizing nonlinear SSMs be contractive
33
Interpreting nonlinear SSMs as stacks of linear dynamical systems
33
Experimental Details and Discussion
34
K.1 Deriving the Empirical Scaling of DEER
34
K.2 Details and Discussion for mecan-fiele RNN experiment
34
K.3 Additional experiment for the mean-field RNN: other optimizers and wallclock time
36
K.4 Additional details for the two-well potential
37
K.5
Building Stable Observers for Chaotic Systems
38
K.6
Numerical computation of the discrete-time LLE
39
Discrete and Continuous Time LLE
39
17

---

A Brief Overview of DEER/DeePPCR
This section provides background on DEER/DeepPCR needed to support section of the main text. Other options for further background on DEER are sections 2-4 of Gonzalez et al. [3] and the corresponding blog post [85].
We begin with a brief review of DEER/DeepPCR [2. 1, 3]. As mentioned in the introduction, the choice of optimizer erucial for this procedure outperform sequential evaluation in terms of wall clock time. Indeed, for this reason DEER uses the Gauss-Newton method (GN) to minimize the residual loss, since GN exhibits quadratic convergence rates near the optimum [21]. Recall from eq. (3) that the i-th step of the DEER algorithm is,
g(++) s() -(s(1)) -(((())
This step requires inverting the TD x TD matrix, J(s( Rather than explicitly inverting it, which generally infeasible, DEER solves for the updates by running linear time-varying recursion [3]:
As (6-1) A 1 As rr(s))),
where
As (6+1)
(i+1)
(
(13)
Unlike the standard sequential rollout, this recursion can be parallelized and computed in O T') time using parallel associative scan [5]. When the number of optimization steps needed for DEER converge the true trajectory is relatively small, DEER can yield faster overall evaluation than the sequential approach. Since -Nwusn converges quadratically when the initial guess is sufficiently close to the true optimum [1, 21J, DEER potentially only requires tiny number of iterations to converge. Our first key result is to prove that DEER always converges globally with linear rate, and will thus always reach this basin of quadratic convergence after sufficient time.
A note about notation The DEER quantities:
residual r(s) RTD
Jacobian (s) E RDDTDD
merit function C(s) E
are functions of the current guess for the trajectory = vec(s1,..., ST) E RTD. As much as possible, try to emphasize the dependence the current guess for the trajectory, but sometimes we will drop the dependence for notational compactness.
B Merit Function is PL
This section provides a proof main text Proposition 1. We first note that Proposition applies to optimizing any nonlinear sum of squares problem where L(s) |r(s)]|, not just the consider this paper (defined in eq. (2)).
Proposition (Proposition 1). The merit function L(s) defined in eg. (2) satisfies eq. (8)
;= inf o (J()).
Proof. Observe that
and C(s") = 0. Substituting these expressions into the PL inequality in eq. (8) we obtain
VL(s) =
J() pr' r.
Therefore, if is full rank, then the merit function L is PL, where
# inf Amin JJs))
inf 22 (s)))
18

---

To be precise, we must have > O for L to satisfy the definition of PL. Therefore, condition that must apply for be PL that we must have inf, min (J(s)) O. We note that the proof strategy of Theorem 2 ensures that se3"" mmin (J(s) > 0 if we assume eq. (22), which holds for dynamical systems that are globally contracting.
By the chain rule, eq. (22) also holds for functions of the form f(s) = (WWs), where W E D and é is a scalar function with bounded derivative that is applied elementwise. In particular, such a function o(W's) satisfies eq. (22) whether or not it is globally contracting. This function class extremely common in deep learning (nonlinearities with bounded derivatives include tanh, the logistic function and ReLU).
In our statement and proof of Proposition 1, we deliberately do not specify the over which take the infimum. The result is true regardless of what this set is taken to be. The largest such set would be RTD, but other sets that could be interest are the optimization trajectory {s()), N), or alternatively neighborhood of the solution s", We discuss further in Appendix C.
Some more general notes on the PL inequality The PL inequality or gradient dominance con- dition is stated differently in different texts [30, 29, 86]. We follow the presentation of Karimi et al. [28]. Karimi al. [28] emphasizes that PI. is often weaker than many other conditions that had been assumed in the literature to prove linear convergence rates.
We note that the PL inequality stated in eq. (8) is not invariant the scaling of C. However, in Definition 3 of Nesterov and Polyak [30], they broaden the definition to be gradient dominant of degree p [1, 2]. The PL inequality we state in eq. (8) corresponds gradient dominance of degree 2. Note that gradient dominance of degree is scale-invariant.
C Merit Function PL Constant is Controlled by Largest Lyapunov Exponent of Model
This section provides the proof of main text Theorem 2.
Theorem (Theorem 2). Assume that the LLE regularity condition from eq. (10) holds. Then if X + the PL constant of the meril function
functinn in (8) satisfies
1
11
1 (xx(-11111)
S
Ve s min
(14)
ET -1
If X then the bounds are instead
1 aT
1 2D bV 1
V min
Proof. We present two proofs. A shorter, direct proof of (14) assuming ||. is the standard Euclidean norm, and then more general version in Appendix C.I. which will be useful later on.
Notice that the residual function Jacobian (4) can be written as the difference of the identity and a T-nilpotent matrix N, as
J LTD - N with N OLD
Because N is nilpotent, the Neumann series for J-1 finite sum:
T-1
J- = (IDD
Σ N
(15)
R=D
Straightforward linear algebra also shows that the norms of the powers of this nilpotent matrix are bounded, which enables one to upper bound the inverse of the Jacobian
T-1
T-1
1- AT 1-
12 S a e Ak and therefore |J <
Σ K=O
N
2
S
Σ ae
(16)
k=
The powers of N are closely related to the dynamics of the nonlinear state space model. We pro- vide dynamical interpretation below, in the paragraph "The dynamical interpretation of and its powers".
19

---

To lower bound ||3"| Ila. we observe that by the SVD, property of the spectral norm is that
112
sup 1= 1
I
JJJy.. y.
(17)
We pick two unit vectors u and v, both in IR
TDD
that are zero everywhere other than where they need
to be to pull out the bottom-left block of J (i.e., the only non-zero block in N which is equal to JIJr-I....d5 J2). Doing so, we get JJ-'v= (J J2)õ,
where ū and e are unit vectors in R and are equal to the nonzero entries of u and v.
Note, therefore, that because of eq. (17), it follows that
/ JJrJT-1.... J2) s || ||
(18)
Le. we also have lower bound on 112.
Furthermore, choosing ū and ö to make
ãJ (JTJT-1.... J2) ē = | 1 J2ll2.
we can plug in this choice of ū and ē into eq. (18), to obtain
|| J...
J2I12 2 _ || 1
Applying the regularity conditions (10) for k T and = 2 obtain
b M(T-1) < ||-
(19)
Because
Amin (JJ")
JJ-1 B
the result for À * O follows by applying eq. (16) and eq. (19) all s(0) along the optimization trajectory.
Note that any choice of å and 0 results in a lower bound, i.e. we could also have targetted the block identity matrices. So, also follows that S ||J- and so
max (b A(T 1), 1) < ||-
Finally, let conclude by considering the case x 0. In this setting, the lower bound on v/u follows from L'Hôpital's rule. For the upper bound, we again must lower bound ||J 11|2- To do so, we leverage the relationship between spectral and Frobenius norms, namely that for an n X n matrix A.
| Tz
& ||Al12 I ||A||e.
(20)
We can find the squared Frobenius norm, i.c. the entries. The squared Frobenius the sum of the squared Frobenius norms of the blocks. We know that each block has spectral norm lower bounded by so each block also has Frobenius norm lower bounded by b. Therefore, summing up over all of the blocks, it follows that
||- which the sum the squares of all of
factors over the block structure of the matrix, i.e. ||J-1|||
T(T + 1)
6²
1
2
and
|J |J'1|E
VTD|J -1
112.
Putting these equations together, it follows that
T(T 11)
bv
TTD J J- 1|2
or
T
by
2
2D
and the upper bound on / when follows from taking reciprocals.
The above proof sheds light on how many dynamical system properties fall out of the structure of J (s), which now discuss further.
20

---

Discussion of why small min (J(s) leads to ill-conditioned optimization Recall that our goal find lower bound on the smallest singular value of [(s,, which denote by m (J))). This quantity controls the difficulty of optimizing L. For example, the Gauss-Newton update is given by (ss)-r(ss) Recall that
maxx (J(s) = 1/omin (J(=)) [[](s) |12.
Recall that an interpretation of the spectral norm [ is how much multiplication by J(s) can increase the length of vector. Therefore,
that very small values of min (J))) result in large values of [J(s I|2, which means that [J(s)) (J(s)='r(s)le can become extremely large as well, and small perturbations in
lead to very different Newttn updates (i.e. the problem is ill-conditioned, cf. Nocedal and Wright [21] Appendix A.I).
Furthermore, we observe that in the X > 0 (unpredictable) setting and the large T limit, the upper and lower bounds in (14) are tight, as they are both (ee Thus, the upper and lower bounds together ensure that unpredictable dynamics will suffer from degrading conditioning.
In contrast, in the À < 0 (predictable) setting, the lower bound on vp converges to -= which is bounded away from zero and independent of the sequence length. Thus, in predictable dynamics, there is lower bound on min or, equivalently, an upper bound on F (J (J-1).
The dynamical interpretation of N and its powers
As shown in the above proof,
T-1 EN N(s)*
J(s) - -IrD-N(s))-
K=
It is worth noting explicitly that
0
0
J2
Of.
((s)
:
where J.
(st-1);
(21)
:8.-1
0
0
Jr
i.e. ((s) collects the Jacobians of the dynamics function along the first lower diagonal. Each matrix power 2 therefore collects length & products along the kth lower diagonal. Thus, multiplication by (ss) Σ N(s)* recovers running forward linearized form of the dynamics, which is one of the core insights of DeepPCR and DEER [2, 1].
21

---

Concretely, in the setting where T = 4, we have
lD
000
0
lD
N
00000
Ip
ID.
0
0
0
0
J2
0
0
N
0
J3
0
00
0
J4
0
0
0
0
2
J3J2
0
00
JAJ3
0
0
N
0
0
J4J3J2
0
0
LD
0
0
1D
0
0
J-1=
J3J2
J3
ID
0
JaJzJ2
Ja
J4
Ip/
Connection to semiseparable matrices and Mamba2 Having depicted the structure of J we note the connection between J- in this paper and the attention or sequence mixer matrix M in Dao and Gu [87], which introduced the Mamba2 architecture (see equation 6 or Figure 2 of Dao and Gu [87] for the form M, and compare with J-1 above).
Mamba2 is deep learning sequence modeling architecture. Its sequence mixer in each layer has at its core linear dynamical system. Dao and Gu [87] observe that while linear dynamical system (LDS) can be evaluated recurrently (sequentially) or in parallel (for example, with parallel scan), it can also be evaluated multiplying the inputs to the LDS by the matrix M. Since each DEER iteration is also a linear dynamical system, with the transition matrices given by {Ja H- it follows that M in Dao and Gu [87] and J- in our paper are the same object, and so results about these objects from these two papers transfer.
In particular, we observe that, in the language from Dao and Gu [87], the J-1 consider in this paper is D-semiseparuble (see Definition 3.1 from Dao and Gu [87]). Thus, any efficient, hardware- algorithms and implementations developed for D-semiseparable matrices could also be ap- plied to accelcrating each iteration of DEER, though we note that Dao and Gu [87] focus on the I-semiseparable setting, which they call state space dual or SSD layer. any case, using these connections to accelerate each iteration of DEER and related parallel Newton algorithms from systems implementation perspective would be an interesting direction for future work.
framing of Theorem 2 based on global bounds on ||Je||2 We chose to prove Theorem 2 using condition (10) in order to highlight the natural connection between the smallest singular value of and system stability (as measured by its LLE). However, an assumption with a different framing would be to impose uniform bound the spectral norm of the Jacobian the entire state space:
sup [|]() S p.
(22)
=ERD
For p 1, this assumption corresponds global contraction of the dynamics [7].
we replace the LLE regularity condition (10) with the global spectral norm bound (22) in the proof of Theorem obtain that the PI. constant bounded away from zero, i.e.
0 p" -
inf SERTD
of (J(s)).
In particular, if the dynamics are contracting everywhere (i.e., p 1), the condition (22) guarantees good conditioning of J throughout the entire state space.
22

---

Discussion of the LLE regularity conditions The LLE regularity conditions in eq. (10) highlight the more natural "average case" behavior experienced along actual trajectories s RTD This "average case" behavior is highlighted, for example, by our experiments with the two-well system (cf. Section 5 and Appendix K.4), where even though global upper bound on [J+(s)) over all of state space would be greater than 1 (i.e., there are unstable regions of state space), we observe fast convergence DEER because the system as a whole has negative LLE (its trajectories are stable on average).
We also note the pleasing relationship the LLE regularity conditions have with the definition of the LLE given in eq. (5). Note that in the LLE regularity conditions in eq. (10), the variable &o denotes the sequence length under consideration. Taking logs and dividing by k. we therefore obtain
log b k
1
log
+ +
log (1Jc4R-1Je44-- -.ll S K.
A
Therefore, as + T. and T + oo (i.c., we consider longer and longer sequences). we observe that the finite-time estimates of the LLE converge to the true LLE A.
We observe that as g! approaches the true solution s", the regularity conditions eq. (10) become increasingly reasonable. Since any successful optimization trajectory must eventually enter neigh- borhood of s", is natural to expect these conditions to hold there. In fact, rather than requiring the regularity conditions over all of state space along the entire optimization trajectory, one could alternatively assume that they hold within a neighborhood of s*, and prove a corresponding version of Theorem 2.
We now do so, using the additional assumption that L-LLpschitz.
Theorem 6. I is L-Lipschitz, then there exists a ball of radius R. around the solution s", denoted B(s", R), such that
Vs E B(s",R)
J [mii (J(s)) - -miin (J(s*)) S LR
Proof. The argument parallels the proof of Theorem 2 in Liu et al. [88].
fact stemming from the reverse triangle inequality is that for any two matrices A and B,
mmin (A) T (B) ||A B|
Applying this with A J(s) and B J(s"), obtain
Omin (J(s) > amin (J(")) - [[](s) J(s"))
If the Jacobian J-) is L-Lipschitz, then
[J(s) - J(s")] < T. -s'll-
Combining, we get
mmin (J()) 2 mmin (J(s")) )- - L||s - =*1
and
Omin (J(s") 2
mmiu (J(s)) L|ss s*11,
which gives
mmin (J(s)) L|ls s*|| S amin (J(s) < Tmin (J(s") )) + Llls - s*
Ensuring that |s s*1| < R completes the proof.
A consequence of Theorem 6 is that if the system is unpredictable, then there exists finite ball around s* where the conditioning of the merit function landscape is provably bad.
As concrete example, suppose that miin (J(*)) the loss function inside the ball B(s", R) is e + If is small (bad conditioning) then R can be chosen such that the PL constant inside the ball (ss", RR) is also small.
and L = 1. Then at best, the PL constant of
23

---

Controlling maxx In our proof of Theorem 2. we proved upper and lower bounds for mmin (m)n(((s) that depended on the sequence length T. We can also prove upper and lower bounds for ( (s)), but these do not depend on the sequence length.
Assuming condition (22), upper bound max (J) is straightforward compute via the triangle inequality,
maxx (J) = ||
= | N
1 ||
Recalling the definition of N in (21), observe that it is composed of {Je) along its lower block diagonal, and so we have
|[N(s)]l2 = sup [|J{(
sup SERTD
[|(s) sup [s](
SERD
Elaborating. for particular choice of trajectory s E RTD, |[[(s)]2 controlled by the maximum spectral norm of the Jacobians (s,) along this trajectory. Analogously, suPseR79 |[N(s)[2- i.c. the supremum of the spectral norm of N(s) over all possible trajectories s E RTD,, i.c. the optimization space-is upper bounded by sup .UD.RRD [s)(s) i.e. the supremum of the spectral norm of the system Jacobians over the state space RD
Thus, it follows that
max (J) S 11+p p.
(23)
Importantly, the upper bound on maax (J) does not scale with the sequence length T.
To obtain the lower bound on max (J), we notice that it has all ones along its main diagonal, and so simply by using the unit vector e1, we obtain
e Jei = S mmax
(24)
Condition number of J Note that the condition number k of a matrix is defined as the ratio of its maximum and minimum singular values, i.e.
T max (J)
s(J)
mmin (J)
However, because our bounds in eq. (23) and eq. (24) on max (J) do not scale with the sequence length 7' it follows that the scaling with T' of an upper bound on x(J)-the conditioning of the optimization problem-is controlled solely by the bounds on Tmin (J) that we provided in Theorem The importance of studying how the conditioning scales with T stems from the fact that we would like to understand if there are regimes particularly involving large sequence lengths and parallel computers-wiie where parallel evaluation can be faster than sequential evaluation.
C.1 Generalized Proof that the Largest Lyapunov Exponent Controls the PL Constant
Lower Singular Value Bound Recall the following sequence of observations.
Amin (JJT) iiin (J)
J-
max
Thus, to lower bound the eigenvalues of JJ" desired, we can upper bound the spectral norm of
General Bound As discussed in the main text, the predictability of the nonlinear state space model characterized by the products of its Jacobians along a trajectory. We will need to control how this product behaves. To reduce notational burden, we will drop the DEER iteration superscript i. In particular, we will assume that there exists a function gJ No + such that
Jr Jk-1 Jalle S 9 i)
24

---

holds for all products J 1........ Jj with k > i, where || lle is the matrix operator norm induced by the vector norm I| Ils. Intuitively, the function 9J measures the stability of the nonlinear state space model. For example, suppose the model is contracting with rate 1. Then the product of Jacobians exponentially decreases, which we can write as
g)(i) = ap',
for some 2 1. The larger the value of a, the larger the potential "overshoot", before exponential shrinkage begins.
Lemma 1. Let lle be matrix operator norm induced by the vector norm -le. Suppose there is a function y. No I such that
JR-1 Jille s gs(k
holds for all products JR-1 Je with k> Define
GT(T) Σ 95(i).
OSiET
Then
|J S GT(T).
Proof. Let y= = =J-1x, By backward substitution for the blockwise entries of J-1x, we have
= (J æz.
SE[K]
Omitting the subscript in the norms for brevity and applying the triangle inequality and the induced-norm property,
|lyl 1
1l S Σ 1 Jill |
kEl2
ME|A
By assumption, -.l. S g kk i). Hence,
|||| S Σ Σ
[s - 1) |2i||
Σ KE(T]
|
gi = Σ
|لlll (T
RE(T] RE(M]
Am
iE[T]
Since G)(()) is nondecreasing in I, the largest multiplier in these sums is GT(T)) In the worst case, |x| = || |. Thus,
| |J-'x||
ly
||
<
GT(T).
x
x
This completes the proof.
Remark (Contraction in The Identity Metric). Recall that system is contracting in the identity
metric when the system Jacobians have singular values less than one:
Vi,
|| S = JT < 11
In this case, take
g)(s) p!.
Then, by Lemma 1,
D -1 p-
p p 1,
|- S
p
(25)
1=
p 1=
where in the case 1, there are T' summands and each term equals 1.
25

---

Remark (Contraction in Time- Varying. State-Deeendent Metries), Recall that system is con- tracting in metric M, (s,, i) if the following linear matrix inequality satisfied
Vi E [T - 11, J MM41 J, 3 22 M..
Equivalently, this condition can be written
constraint
E [T - 1|, IMM ' S p
Using these metrics, define the block-diagonal, symmetric, positive-definite matrix
M diag(M , Mr)
as well the similarity transform of the residual function Jacobian, based on this matrix
JM = 1/2 J JM-1/2
Then the off-diagonal block entries of JM
M J, M, 1/2 for i -1-
while its diagonal block entries are the identity matrix. If the off-diagonal blocks of JM satisfy product bound function 9Jm (j) as in Lemma 1, then JM has norm bounded by GJM (T). Hence,
1| = M 1/2 - -M 1/2 M
I | 1
M M 1/2
|M 1311
|M--1/2
1/21| |JM|| |
S 1/2| GJM(()) |
KM GJM (T),
where
Amax (M)
RM
Amin (M)
In this case, may again take IJ2c (1) = p², and we obtain the bound
p²
||ー S KM
Σ
KM P
P 1# 1,
0Sj<T
KA
p p=
Remark 3 (Contraction After Burn-In). Suppose that
gr(ر) S Aj Aj
where a 2 and measures the degree "overshoot" the system can undergo before eventually converging, and À > 0. In partieular, assume for concreteness that
|| S
Then, the product of two Jacobians can grow, if a > e^, since
|| Jell S ae a> In general, the product of Jacobians can transiently grow (i.e., overshoot) for
koversboot
log a
time steps, at which point the product of & > hiowerstoue Jacobians will remain less than 1, and will in fact decay exponentially with rate X.
In this case, by Lemma
-x -
|ー
<
Σ j=0
c ~x
e
26

---

D DEER Merit Function Inherits Lipschitzness of Dynamics Jacobians
This section provides proof of main text Theorem 3.
Theorem (Theorem 3). If the dynamics of the underlying nonlinear state space model have L- Lipsehitz Jacobians, Jacobiaan, L... vt>1, ERD: [J((8) Je(s')Il Li|s s'11,
then the residual function Jacobian J is also L-Lipschitz, with the same L.
Proof. By assumption, for each 6, Vs, E ER::
[[ - Je(s-)l & L |1st s/ll2.
Define D, = J+(o)) - J(( and
D = J(s')
- J(s).
Since D places the blocks D, along one subdiagonal, we have
|DI2 max Delll.
But each block D, satisfies the Lipschitz bound
| S L ||s - soll2.
|D|22 max |D2ll2 s L max || sell2
L ||s -s|l2.
Hence, it follows that
[J((') J(s)ll2 = |DI|2 I L ||s s||2-
Thus J L-Lipschitz.
E DEER Always Converges Linearly
This section provides proof of Theorem 4.
While proofs of global convergence and this can be exploited to provide a global proof of convergence. particular, we will exploit the hierarchical nature of DEER, which is reflected in the fact that J and J are lower block-triangular
challenging in general for GN, DEER is highly structured,
Theorem (Theorem 4). Let the DEER
updates be given by eq. (3), and let s(i) denote the i-th iterate. Let e == s()) s* denote the error at iteration i, and assume the regularity condition eq. (10). Then the error converges to zero at a linear rate:
|| S X 8 (0)) 12
for some constant Xw independent i and a convergence rate 0 < 3 <
Proof. Our general strategy for deriving DEER convergence bounds will be to fix some weighted norm || |!w = W//2 W-1/2 || such that each DEER step is contraction, with contraction factor 3 E [0, 1). This will imply that the DEER emor iterates decay to zero with linear rate, as
| www S 8 (e I|w-
To convert this bound back to standard Euclidean space, we incur an additional multiplicative factor that depends on the conditioning of W:
Amaxe (W)
|ee(0) 1 S Xw (0
where =
(26)
Amin (W)
27

---

DEER as Contraction Mapping Recall that the DEER (Gauss-Newton) updates are given by
((+ =s( --
Recalling that r(s*) and subtracting the fixed point from both sides, we have that
(f++)) el J (s( + +J (s() e - J (s(0) (r(s!() r(s))
This equation can be written using the mean value theorem as
e (f+1) II J (- ( el
B = / J(s" + re dr
where
From this, we can conclude that the DEER iterates will converge (i.e., the error shrinks to zero) if
| - J
|J (s(0) (J(s(0) B()) 11w S 1.
(27)
Constructing the Weighted Norm
We will choose a diagonal weighted norm, given by
W = Diag(ID, wº Under the norm induced by (28) we have
wIID, ID. .... w (D) RDDTDD
0.
(28)
(s(()))) - ll I 2wpp
(29)
J ( < 1 (we*)T
(30)
wed
where o upper bounds || la over all states in the DEER optimization trajectory.
Multiplying (29) and (30) yields
1 (weA)T
' (s ( llw (s(()))
B 2awp
(31)
we
To ensure the right-hand hand side of (31) does not exceed prescribed | E [0, 1), choose
w
(32)
2pa + Be
With this choice,
2amp
we 1,
and
3,
(33)
1- we
so the geometric series in (30) is convergent and the bound in (31) holds for all T, because
(we>)T
J |- (J](s(0) - B llll <
2awpp
(we>)²) B.
wed
This shows that we can always pick a weighted norm so that DEER converges with linear rate that norm. Converting back into the standard Euclidean norm using (26) and substituting in the condition number of 1/12 one finds that
(2pa + Be
| 125
B (0e 2
(34)
3
Thus, the DEER converges with linear rate towards
Remark 4. The multiplicative overshoot factor arising from the conditioning of W grows expo- nentially in the sequence length T, leading potentially to long convergence times. Indeed, quick calculation shows that the number of steps needed to bring the DEER error is upper bounded as O(T)) because of this multiplicative constant.
28

---

Remark 5. One can ask under what conditions choosing w = 1 in (32) is possible, which eliminates the overshoot. We will address this more detail the next section. To provide a simple result here, we can assume that the system is contracting at every time step so that
and 1. Then we have that
B
1
2pa + Bed
2ea + Be
Solving for X, we have that
3
A S log
log(3),
2
then w can be chosen to be equal to one, meaning the DEER converges globally with rate 8 and no overshoot,
DEER Converges Globally with Small Overshoot for Sufficiently Strongly Contracting Systems
In this section we show that DEER converges globally to the optimum s* when the nonlinear state space model (1) is sufficiently strongly contracting. To do so, we first briefly recall the assumptions of Lemma 1. Let IIs be the matrix operator norm induced by the vector norm || le. Suppose there function gu No R such that
Ja-1Jk-2- ... Jalls S gu(k - i)
holds for all products JR-1 Jr with k > i. Define
GT(T) Σ 9s(1).
0Si<T
Then
|J s GT(T).
For example, if there is no structure which can be exploited in the products of Jacobians J., we may consider the "one-step" growth/decay factor
Vt,
| s
which yields
Gj(T) E 9r(i) 1- XT
9r(5) = A2
USi<T l-
Theorem. DEER exhibits linear, global convergence the optimum s* with rate B E 10,1) in the matrix operator norm lle if
(g) 8
Proof. Recall that the DEER (Gauss-Newtood updates are given by
(++1) (() - J (s (((())
Define the error at DEER iteration (i) as (4) (()) - s" Recalling that r(s") = =0 and subtracting the fixed point s* from both sides, have that
(++1) = ( - J- +)))(s())))))))))))))))))))))))) r(s*) = e - J-)((s()) (r(s(0) r(s*))
This equation can written in terms the mean value theorem as
where B = / J(s" + (e( dr
(f+1) I - J-J-'(sf) B (9) (1)
29

---

This follows from the identity:
r(s (0) r(()) [( J(Ts +1- dr (s (s) (s-
(f,
J(s J(s"
roll)) dr e(
This identity can be proven by starting from the fundamental theorem of calculus, by letting s(r) s + reli) =s + T(s(() - =s), E [0, 1],
which defines a straight-line path from s()), The fundamental theorem of calculus then says that
r(s") -(sss((s// -/h d -न -r(((()) dr. r(s -
Applying the chain rule inside the integral gives the result, because
d
d dr s()) J(s + re!0) e0)
dr r(s(r)) = J(((T)
From this, we can conclude that the DEER iterates will converge (i.e., the error shrinks zero) if
J- (s()) 1
|J (s!0 (J(S!)) - (()) lle <
(35)
By Lemma we have that
| [e < (|J- lle (s() B lelle s 22gr(1)G>(() |1
Thus, if there exists some 3 E [0, 1) such that
(g)( < 8,
then the DEER error converges globally to zero in the weighted norm:
||el S (le
Corollary. Suppose the state space model is contracting in constant metric M, ie.,
J,M-1/2 /
|/2
|| S e <1
If is sufficiently small, in particular if
2 2 B
3
then the DEER errors converge to zero with rate B.
Proof. Suppose the state space model is contracting in constant metric M, so that
|| J,M-1/21 = ||J.|| S e
for all 6, Then, by Lemma we have that
1- 1-e2T XxT
| (e((+1) a < J |J J|
B lIa(le I| s
T- 26 ||
Thus, in order to achieve linear convergence of the DEER iterates with rate B E (0,1),
||e(+ 5 1I
||e(9) ||M S (811e(0)
we require that
XT
2e
SB<1.
1 e
A simple sufficient condition for satisfying this inequality is
8
2+ 3
3'
or,
A < - log(3).
30

---

Number of Steps to reach basin of quadratic convergence [0, 1) such that
Let us assume that there exists B
8
2 2 B
then the number of steps to reach the basin of quadratic convergence is upper bounded as
e e>
ki S log
log
G Alternative Descent Techniques & Worst/Average Complexity
DEER uses the Gauss-Newton algorithm, which converges quadratically near the optimum but can be slow outside this basin. This motivates inexact GN methods that guarantee a certain loss decrease per step, such as line-search and trust-region techniques. These trade increased computation and possibly more iterations for faster convergence guarantees.
In practice, found that plain GN reliably converged quickly the global optimum in contracting systems, so such safeguards were unnecessary. Still, it is useful to analyze DEER's worst-case path to the quadratic basin.
Many inexact GN variants achieve global convergence from any starting point. These include step- size schemes that approximate continuous flow [89], trust-regions that bound update size (yielding ELK when applied to DEER [3]), and backtracking line search ensuring loss reduction at each step [21].
One can also use simpler algorithm outside of the basin of quadratic convergence, and then switch to GN when needed. We will consider this latter option, and choose gradient descent our simpler algorithm. Because the merit function is PL (see section 3.1), the number of steps required for gradient descent reach the quadratic convergence region scales as:
|
k2
log
(36)
/
where ||r(0]] is the residual at initialization. For unpredictable systems, μ may shrink arbitrarily with increasing sequence length T. leading to unbounded growth in the number of optimization steps ky. By contrast, for predictable systems, μ remains bounded, implying that the number of optimization steps does not increase with sequence length. Since cost of sequential evaluation always increases with T, DEER can, even in the worst case, compute the true rollout faster than sequential evaluation for predictable systems-cspeciaalll for long sequences. Indeed, assuming the system is contracting with rate ed 1, then the number of steps needed the reach the basin of quadratic convergence is O (log [[))
Thus, if the initial error grows polynoaial in T, i.e., ||r(0)|| TP, then this implies that the number of gradient descent steps needed to reach the basin of quadratic convergence is only Olog T). and thus the total computational time is O(log In practice, for randomly initialized DEER, observe p
In practice we observe that DEER converges much faster than the worst-case analysis (36) would suggest. In particular, we observe that DEER converges in roughly log 1, steps, even for unpre- dictable systems. This behavior can be explained with simple "two-phase" model, wherein the DEER iterates move towards the basin of quadratic convergence at rate which is independent of the PL-constant (see Appendix K.I).
H Proof of Size of Basin of Quadratic Convergence
This section provides proof of Theorem 5:
Theorem (Theorem 5). Let μ denote the PL-constant of the meril function, which Theorem 2 relates to the LLE X. Let L denote the Lipschitz constant of the Jacobian of the dynamics function (s). Then, MLL lower bounds the radius of the basin of quadratic convergence of DEER; that if
(r(s
p
31

---

then g00) is inside the basin of quadratic convergence. In terms of the LLE X, it follows that if
|r(s - ( e²
1
a²
edT -
then s( is inside the basin of quadratic convergence.
Suppose we are point s(1) E IR. change in the trajectory obtained from eq. (3) is,
TD (i.e. DEER iterate i), and we want to get s(++1), The
(+) == -
(where the iteration number will hopefully be clear from context). The merit function is L(s) [|r[s)]] so if we can get some control over ||r(s(()]]2,, we will be well on our way to proving quadratic rate of convergence.
First, leveraging the form of the Gauss-Newton update, we can simply "add zero" to write
r(s ((s()+))) = rs((() As(())
= =r((()) + As (1) - r(s(1) J )
Next, we can write the difference r(s(() + As(") -r(s()) as the integral of the Jacobian, i.c.
r(s(0) + As())) - r(s = //¹ J (s + = As dr.
Therefore,
r(s r(s(++1)) = L// (J (I(S0) + + -((s(1)) As dr
Taking &2 norms and using the triangle inequality, it follows that
(s ||r(s(+D))ll - [. [[ (s( + = J(s( As || dr.
Now, if we assume that J is L-Lipschitz and use the definition of spectral norm, it follows that
(J (0)) + + ((s As( S ST.
and so taking the integral we obtain
L
[[((s
- |ss
2
L
2 r(s
(1)
-
J(
By definition, v is a lower bound on all singular values of J(s()), Therefore, ||(((())~11 < 1/v for all i, and it follows that
|r(t(+
-
L (1))
(37)
24
which is the direct analogy of Boyd and Vandenberghe [31, 9.33]. To reiterate, here L is the Lips- chitz constant of J, while μ infien 2 (J(s())).
While this is quadratic convergence result for GN, this result is not useful unless r(s(i+1))|l2 S [ (i.c would backtracking line search accept this update). However, if we have (1(s(() I| < بي then every step guarantees reduction in r because in this case
[|((s(+)))]22 [|)
Therefore, we have [|r(s()) < 2u for all > i. Thus, we have related the size of the basin of
L quadratic convergence of GN on the DEER objective to the properties of J. Note that with linear dynamics, each Jr is constant in and so each J, is O-Lipschitz. Thus, the basin of quadratic
convergence becomes infinite. Intuitively, if J, doesn't 't change too quickly with 8, then DEER becomes more and more potent method.
32

---

log T
LDS
nSSM
LDS +
log T
T
LDS
Figure 4: Equivalence between a contractive nSSM and an Olog T) stack of linear state-space models. Contractivity implies that the nonlinear dynamics can be decomposed into a hierarchy of Olog 7') layers of linear SSMs, or linear dynamics systems (LDS), each of which can be evaluated in Oloog T) time by a parallel sean.
I Parameterizing nonlinear SSMs to be contractive
In this section, we highlight practical strategy for speeding up the training of nonlinear state space models (SSMs) based on our theoretical findings.
Our results indicate that nonlinear SSMs with negative largest Lyapunov exponents (LLEs) are ef- ficiently parallelizable. To exploit this during training, one must ensure that the model maintains negative LLEs throughout optimization. One straightforward and effective method to achieve this is by design, through parameterization. In particular, by introducing an auxiliary variable enforce the desired constraint (in this case, negative LLE), and then performing unconstrained optimization on this variable.
This strategy is particularly well-suited to neural network-based SSMs. For example, consider the scalar nonlinear SSM:
Ft = tanh(wr;-1 + u4)
To guarantee negative LLE, it suffiices ensure that the Jacobian norm is strictly less than one:
[Je] w sech2 u)] S [w]
Thus, enforcing |w| < is sufficient. This can achieved by reparameterizing w tanh(b), where b is trainable, unconstrained auxiliary variable. This guarantees that E (-1, 1) for all finite b, ensuring contractivity and, hence, negative LLE, A similar argument holds in the multivariate case, using the spectral norm.
J Interpreting nonlinear SSMs as stacks of linear dynamical systems
As mentioned in our Discussion in Section 6, an important implication of our results that a con- tractive nSSM can be interpreted as a hierarchical composition of linear state-space layers (SSMs), or equivalently, linear dynamical system (LDS) layers. Each layer can be evaluated in Olog T)) time with a parallel scan, and the total number of layers required scales as O(log T). This perspective shows that nonlinear temporal dependencies can be captured through a logarithmic-depth stacking of linear dynamics. Figure 4 provides schematic illustration of this equivalence.
More explicitly, each iteration of DEER is given by the linear dynamical system
(++1) = [()(0)
(0)) (s (1+2)
(38)
33

---

Therefore, we can interpret each "iteration" (2) of DEER as sequence mixing "layer" (i), where the sequence- mixing layer is an imput-eependent switching linear dynamical system, like Mamba [55]. The inputs to "layer" (i ++1] 1) is the state trajectory of the immediately preceding "iteration" or "layer" (i). Because we prove that DEER converges linearly in Theorem 4, it follows that contractive nSSM can be simulated in Oloo 7') LDS layers of the form shown in eq. (38), assuming the initial error grows polynomially in the sequence length.
K
Experimental Details and Discussion
All of our experiments use FP64 to, as much as possible, focus on algorithmic factors controlling the rate of convergence of DEER, opposed to numerical factors. As noted in [3], DEER can be prone to numerical overflow in lower precision. While such numerical overflow can be overcome by resetting NaÑs to their initialized value, such approach resets the optimization and leads to rates that are slower than what Gauss-Newton would achieve in infinite precision (exact values in R).
K.1 Deriving the Empirical Scaling of DEER
In our experiments, we observed that DEER typically converges in (log(1/H)) steps (see, for example, Figure 2). To understand this scaling behavior, we propose simple two-phase model of DEER convergence. In the first phase, the iterates approach the basin of quadratic convergence at linear rate, guaranteed by Theorem 4. In the second phase, rapid quadratic convergence occurs, typically requiring only one or two steps to reach the true solution (up floating point precision).
Although Theorem 4 shows that, in unpredictable systems, the overshoot factor may be exponen- tially large in the sequence length T. this reflects a worst-case analysis. In practice, DEER behaves as though the overshoot factor is negligible. To formalize this observation, recall from Theorem that the residuals satisfy the linear convergence bound
|ra/l S KmrB"|roll
for some 3 E [0, 1) and X 2 where is always independent of T. our two-phase model, we assume that X is also independent of T, even when the largest Lyapunov exponent À is positive. We upper-bound the number of steps required enter the basin of quadratic convergence, whose size scales as /ı/ L (as given by (12)). Solving
Xw (Ioll
L =
X |roll
k == log B
log
(39)
we recover the empirically observed logarithmic scaling.
K.2 Details and Discussion for mean-field RNN experiment
We rolled out trajectories from a mean-field RNN with step size 1 for 20 different random seeds. The dynamics equations follow the form
1 Wtanh(st) + uy
for mild sinusoidal inputs u, We have s, RD, where in our experiments 100. Note that because of the placement of the saturating nonlinearity, here 8, represents current, not voltage.
draw each entry Wa ~ N00,, s²/D), where is a scalar parameter. We then set Wi: = 0 for all i (no self- coupling of the neurons). A key point of Engelken al. [32] is that by scaling the single parameter the resulting RNN goes from predictable to chaotic behavior. While Engelken et al. [32] computes the full Lyapunov spectrum in the limit D - co, for finite D we can compute a very accurate numerical approximation to the LLE (cf. Appendix K.6). In Figure 5, we verify numerically that there is monotonic relationship between g and the LLE the resulting system, and that the min- max range for 20 seeds is small. Accordingly, when making Figure 2 (Center), we use the monotonic relationship between y and the L.E from Figure to map the average number of DEER steps (over 20 different seeds) needed for convergence for different values of g to the appropriate value of the
In the design of the weight matrix W, we follow Engelken et al. [32]. In particular,
tid
34

---

Median LGE
XinoMax rxage
11.03
-01
3.5 0.55
Figure 5: Robust relationship in mean field RNN between variance parameter y and LLE of the system. For 20 seeds, we observe a robust and non-decreasing relationship between the sealar parameter g and the LLE of the resulting mean-field RNN. The plot above is made for 50 different values of g from 0.5 to 2.0 (linearly spaced). We estimate the LLE over a sequence length of T = 9999.
LLE. We use 50 values of T from 9 to 9999 (log spaced) to make Figure 2 (Center). We highlight T 1000 Figure (Right).
For the purposes of Figure 2, we define
e² XT
# =
ic. the lower bound μ from Theorem 2, with 1.
In Figure we observe that around g 1.2, the RNNs have LLE around 0, which is the threshold between predictability and chaos. Working with chaotic dynamics in finite precision for long time series led some interesting difficulties.
First, as discussed in Gonzalez et al. [3], DEER can experience numerical overflow when deployed on unstable systems. While we reset the initialization (in this experiment we initialized S with iid draws from [0, I)), doing so slows convergence. Thus, many of our runs for À >> 0 and large T take the maximum number of DEER iterations we allow (we do not allow more than T iterations, as this is the theoretical upper bound for number of DEER iterations before convergence, cf. Proposition of [3]), which helps explain the slight increase in red space for experiment (center plot of Figure 2) vs. theory (left plot of Figure 2). Note, however, that for T = 1000 (the sequence length shown in the right plot), there is no numerical overflow for the DEER trajectories for any of the 20 random seeds or 50 values of g tried.
Second, observe that for many values of À in the chaotic range, even after the maximum number of DEER steps (T)) was taken, there was still a large discrepancy between the true sequential rollout and the converged DEER iteration, even though the converged DEER iteration had numerically zero merit function. For example, in Figure (Right), there are a series of points in the top right of the graph that all sit on the line T 1000, and while they have numerically zero merit function value, the converged DEER trajectories are quite different from the true sequential trajectories. The reason for this behavior precisely stems from the fact that for large values g (equivalently A), these mean-field RNNs are chaotic. Even working in FP64, if slight numerical errors are introduced at any time point in the sequence (say = 1), then over the sequence length we can observe exponential divergence from the true trajectories, as illustrated in Figure 6. This experimental observation is complemented by our discussion of why unpredictable systems have excessively flat merit functions in Section 3.2, and provides numerical perspective on why ill-conditioned landscapes are hard to optimize: if the landscape extremely flat, many potential trajectories $1:7 can have numerically zero merit function, even in extremely high precision.
35

---

10
1st coordinate of state
10-2
10-5
MAD
10-3
10
2
DEER sequential
10-14
0
200
400
600
800
1000
200 400
600
800
1000
T
Figure 6: Chaotic behavior means numerically zero merit function can still be far from se- quential trajectory. For g = 1.85 and T 1000. we show the final DEER vs sequential trajectory. The DEER trajectory has merit function (2) numerically equal to zero. However: (Left) the mean absolute deviation (MAD) at each time point between the final DEER iteration s(?) and the quential rollout st grows exponentially. This exponential growth of error is a signature of chaos: compare, for example, with Figure 9.3.5 of Strogatz [10]. The saturation of the error eventually occurs because of the saturating nonlinearity present in the RNN. (Right) We visualize the first coordinate both the final DEER iteration and the sequential trajectory, showing that while they initially coincide, they diverge around / 300.
K.3 Additional experiment for the mean-field RNN: other optimizers and wallelock time
In this section, provide further experiments in the setting of the mean-field RNN (Figure 2). In particular, we showcase the generality of our theory beyond DEER (Gauss-Newton optimization), and the practicality of our theory by reporting wallclock times. We consider the setting in the right most panel of Figure 2, where we evaluate a mean field RNN over sequence length of length T 1000.
Quasi-Newton and Gradient Descent Instead of only using Gauss-Newton optimization (DEER) to parallelize the sequence length, also consider other optimization algorithms (quasi-Newton and gradient descent) to showcase the generality of our theory.
We include a quasi-Newton algorithm proposed in Gonzalez et al. [3] called quasi-DEER. Quasi- DEER simply replaces the J, defined eq. (4) with diag(J:), and so is also parallelizable over sequence length with parallel scan. Furthermore, function, which embarrassingly parallel over the sequence length. the top panel of Figure 7, we observe that the number of steps for gradient descent and quasi-DEER to converge also scales mono- tonically with the LLE, as we expect from Thcorem 2. DEER (Gauss-Newton) converges in a small number of steps all the way up to the threshold between predictability and unpredictabilitty (A 0). Intuitively, the performance of the other optimizers degrades more quickly as unpredictability in- creases because quasi-Newton and gradient descent less information about the curvature of the loss landscape.
also include gradient descent on the merit
Even though gradient descent was slower to converge in this setting, we only tried gradient descent with a fixed step size. An advantage first-order method like gradient descent over second- order method like Gauss-Newton (DEER) is that the first-order method is embarrassingly parallel (and so with sufficient parallel processors, the update runs in constant time), while DEER and quasi- DEER use parallel scans (and so the update runs in Ollog T') time). Exploring accelerated first-order methods like Adam [90], or particularly Shampoo [91] or SOAP [92] (which are often preferred in recurrent settings like eq. (I))-or in general trying to remove the parallel scan-are therefore very interesting directions for future work.
Sequential evaluation of eq. (1) can also be thought of as block coordinate descent on the merit function C(s). where the block 8, E IRP is optimized at optimization step (1). The optimization of each block is convex problem: simply minmmize ||x: - ( or equivalently set = f(
As sequential evaluation will always take T steps to converge, we do not include the top panel of Figure?? 7.
36

---

Wallclock time In the bottom panel of Figure 7, we also report the wallclock times for these algorithms to run (our experiments are run on an H100 with 80 GB onboard memory). We observe that the run time of sequential evaluation (green) is effectively constant with respect to X. We observe that in the predictable setting, DEER an order of magnitude faster than sequential evaluation, while in the unpredictable regime, DEER is 1-2 orders of magnitude slower than sequential evaluation. This importance using parallel evaluation only in predictable settings is a core practical takeaway from our theoretical contributions.
the experiment in Figure 7 on a smaller scale than the ex- periment in Figure 2 (Right). In Figure we consider 5 random seeds for 16 values of g equispaced between 0.5 and 2.0. Each wallclock time reported is the average of runs for the same seed. We use a batch size of 1. While DEER (Gauss-Newton) and quasi-DEER effectively do not have a
Further details We
Comvrggnce rales optimizatio algotitim
Quae-Daik
Gradiad Doseat
stge musvErree
step size (they use step size of 1 always). For each value of g, we ran gradient de- scent with the following set of step sizes a: 0.01, 0.1, 0.25, 0.5, 0.6, 0.7, 0.8, 0.9, and 1.0. For each value of g. we then pick the step size a that results in the fastest convergence of gradient descent. For the smallest value of g = 0.5, we use a = 0.6; for g = 0.6, we use 0.5; and for all other values of 9, use = 0.25. Future work may investigate more the step size a, or use a learning rate schedule.
walblodk (co))
adaptive ways
We use larger tolerance of C(+// S 10-4 to declare convergence than in the rest of the paper (where we use tolerance of 10-10) be- cause gradient descent often did not converge to the same degree of numerical precision sequential, quasi- DEER, or DEER. However, this is per sime average eπror the or- der of 10-1, in a system where D = 100 and each state has current the order of 1. Nonetheless, it is an interesting direction for future work to investigate how get gradient descent to converge to greater degrees of nu- merical precision in these settings; and, in gen- eral, how to improve the performance of all of these parallel sequence evaluators in lower numerical
Figure 7: Convergence rates and wallclock time for many optimizers. We supplement the mean- field RNN experiment by also considering quasi- Newton and gradient descent methods (top), and recording wallelock time, including for sequential evaluation (bottom).
precision.
K.4 Additional details for the two-well potential
We form the two-well potential for our experiment in Section as a sum of two quadratic potentials. Concretely, we define the potential é the negative log probability of the mixture of two Gaussians, where one is centered at (0, -1.4) and the other is centered at (0, 1.6), and they both have diagonal covariance. In Langevin dynamics [93, 94] for potential 6, the state s; evolves according to
St+1 St eVo(s))) + V2ew, ldd
where is the step size and we ~ (0, Ip). ) In our experiments, we use e 0.01. $ Accordingly. the Jacobians of the dynamics (those used in DEER) take the form
J. = Ip - <V2 ((st).
'Ntice that this is a discretization (with time step <) of the Langevin Diffusion SDE ds(() -o S (dw(()) where w(t) is Brownian motion [95-97].
37

---

LLE Iterations
Mecit Function
Stepe
0.0
Morit function
1047
Number of DEER
30
-0.5
1010
LLE
10a
10-11
10-11
-20
100
30
100
20
30
5000 10000
DEER Iteration
DEFR
Iroration
Sequenon Length T
Median
Min-Max (20 aeeds)
Figure 8: In this plot, we provide additional information about the behavior of DEER when rolling out Langevin dynamics on two-well potential. (Left) We observe that across 20 random seeds (including different Langevin dynamics trajectories), the LLE for intermediate DEER iterations be- comes negative after the first iteration. Consequently, we observe that the merit function (Center) experiences spike on the very first DEER iteration (following initialization, which was the only trajectory with positive LLE), before trending towards convergence. As the system spends most of its time contracting regions, we observe (Right) that the number of DEER iterations needed for convergence scales sublinearly with the sequence length T. We plot the min-max range for 20 seeds, and observe that even out of 20 seeds, the maximum number of DEER iterations needed to converge on sequence length of ' 10, 000 is around 35.
As a result, the dynamics are contracting in regions where & has positive curvature (inside of the wells, where the dynamics are robustly oriented towards one of the two basins) and unstable in regions where ø has negative curvature (in the region between the two wells, where the stochastic inputs strongly influence which basin the trajectory heads towards). We observe that there are regions in state space where the dynamics are not contracting, the resulting trajectories have negative LLE. Accordingly, in Figure 3 (Right), observe that the number of DEER iterations needed for convergence scales sublinearly. as the LLE of all the intermediate DEER trajectories after initialization negative. These results demonstrate that if the DEER optimization path remains in contractive regions on average, we can still attain fast convergence rates as the sequence length grows.
though
Moreover, a further added benefit of our theory is demonstrated by our choice of initialization of
DEER. Both [1] and [3] exclusively initialized all entries of ( to zero. However, such an initial- ization can be extremely pathological if the region of state space containing is unstable, as is the case for the particular two well potential we consider. For this reason, we initialize (()) at random (as iid standard normals).
An important consequence of this experiment is that it shows that there are systems that are not globally contracting that nonetheless enjoy fast rates of convergence with DEER. This fact is impor- tant because globally contractive neural network may not be so interesting/useful for classification, while a locally contracting network could be.
Futhermore, in this experiment we show empirically that Langevin dynamics can have negative LLE (cf. Figure 3). This results suggest that the Metropolis-adjuster Langevin algorithm (MALA), a workhorse MCMC, may also be predictable in settings of interest, including multimodal distri- butions.
K.S Building Stable Observers for Chaotic Systems
To further demonstrate the applicability of our results and to validate them in the context of non- autonomous systems- we construct nonlinear observers. Observers are commonly used in science and engineering to reconstruct the full state of system from partial measurements [33, 34]. As benchmark, we consider nine chaotic flows from the dysts dataset [98]. According to Theorem (2), these systems exhibit poorly conditioned merit function landscapes and are thus not well-suited for parallelization via DEER, the corresponding observers are stable, then they should be suitable for DEER.
38

---

We design observers for these systems using two standard approaches: (1) by directly substituting the observation into the observer dynamics, following Pecora and Carroll [99]. or (2) by incorpo- rating the observation as feedback through a gain matrix, as in Zemouche and Boutayeb [100]. We then apply DEER to compute the trajectories of both the original chaotic systems and their corre- sponding stable observers. As anticipated by Theorem (2), the chaotic systems exhibit slow con- vergence- -otee requiring the full sequence length-wherean the stable observers converge rapidly (Figure 9).
As with the two-well experiment, initialize our guess for 8 (0) iid standard normals.
Comparison of DEER Steps, Chaotic Systems vs. Stable Observer
Lng(DFER Steps)
:01 101 10+
Observer System
104
35
Rissler
Larenz
Thomas
ABC
Sprottti
Kawczyyski Strizahk
Vallis El Mñña
Chus's Circult
Nosé Hoover Thermosist
Figure 9: Comparison of DEER convergence behavior for original chaotic systems (red) and cor- responding stable observers (blue) across nine flows taken from the dysts dataset. As predicted by Theorem (2), the chaotic systems converge slowly-often taking the whole sequence length T, denoted by the horizontal dashed dune-dee poorly conditioned merit landscapes, while the stable observers achieve rapid convergence
K.6 Numerical computation of the discrete-time LLE
The Largest Lyapunov Exponent (LLE), which often denote by X, is defined in Definition 1. However, for long sequences T, naively computing it would be numerically unstable. Thus, we use Algorithm to compute the LLE in numerically stable way. Note that the algorithm nominally depends on the initial unit vector 2gg. For this reason, we choose 3 different unit vectors (initialized random on the unit sphere) and average over the 3 stochastic estimates. However, in practice observe that the estimate is very stable with respect to choice un, and agrees with systems for which the true LLE is known, such as the Henon and logistics maps.
Algorithm Numerically Stable Computation of Largest Lyapunov Exponent (LLE)
1: Input: Initial unit vector un, total iterations T
2: Initialize: LLE
3: for trtt= to 7 do
4: Compute evolved vector: u4 <- J,04-1
5: Compute stretch factor: de f- |/well
6: Normalize vector: u, t u,/>,
7: Accumulate logarithmic stretch: LLE +- LLE + log X
8: Output: Estimated LLE À <- LLE/T
L Discrete and Continuous Time LLE
We provide the definition of the LLE of discrete-tima dynamical system (often called map) in Definition our paper studies discrete-time SSMs as (1), this discrete-time definition LLE makes sense for our setting. However, many of our experiments involve the discretization of continuous time systems, we want review how the LLESs of discrete and continuous time systems relate to each other. Helpful references on this topic include [7, 101].
39

---

The LLE quantifies what happens to perturbation" dz over time. Does its magnitude || grow or shrink over time? The LLE is that value that makes eq. (6) hold, i..
||6x(()]| ~ e (0))
This notion of the change in the size of perturbation 8x(() over time-as quantified by the LLE- makes sense for both a discrete-ting SSM 2, = fi ( as well as a continuous time system F(x, t).
Let us consider how a perturbation dz evolves in both a discrete time system 24 = f(( and continuous time system á F(x,t). An infinitesimal perturbation óx intimately related to derivatives with respect to æ. Therefore, their variational equations are:
af((x+1) 5x
discrete-time:
ox; = (x(1).
dx
OF
continuous-time:
5x(()) Fx (z((),
Case study: discretizing a continuous-time system
All of our experiments are the forward Euler discretizations of continuous-time systems. So, in this section, we work out what happens to the LLE in such a setting. Let's consider running our continuous-time setting for length of time T.
In this setting, we discretize by timestep Al, our resulting discrete-time map f is given by
f++A1(0:) x4 + AtF( t).
Therefore, it follows that
0f (xx) In Al OF Or
J4
(2x).
Oz
We naturally define the continuous-time LLE the limit of the products of these as At O, i.e.
>e = lim
log
II JKN k=
At-
T
where the number of discrete time-steps N is given by N= T/At.
Notice that this extremely similar to the definition of the discrete-time LLE A we gave Defi nition 1, except division is occurring by the length of the time window T' instead of the number of discrete steps taken
Tt (Of course,
T' At 1.)
Therefore, if we have a discretization of a continuous- time system, and naively plug into our discrete-ting time LLE defined in Definition (i.e., divide by number of discrete-time steps N instead of the length of the time window 7), the resulting >a(At) will satisfy
>a(At)
A
At
Note that naively plugging in this discrete-time estimator would therefore result in a different esti- mate of the LLE depending on the size of the time-step At. However, in most of this paper we still report the naive discrete-time LLE Aa At) as this is the quantity relates to μ via J. The exception our convention in Table 1, where we report our estimates X to better coincide with the intu- itions and expectations of readers with backgrounds in continuous time systems. For all systems in Table we use step size of 0.01, and so one can translate from the reported continuous time LLEs to discrete time LLEs by dividing by 100.
Ultimately, dividing by At, which positive, does not change the sign of X, i.e. whether the system is predictable or unpredictable.
'technically virtual displacement, cf. Lohmiller and Slotine [7].
40

---

NeurIPS Paper Checklist
Claims
Question: Do the main claims made in the abstract and introduction accurately reflect the paper's contributions and scope?
Answer: [Yes]
Justification: Our main claim is that predictability implies parallclizability. We demonstrate this theoretically by analyzing the condition of the merit loss function, and illustrate these findings experimentally well.
Guidelines:
The answer NA means that the abstract and introduction do not include the claims made
in the paper.
The abstract and/or introduction should clearly state the claims made, including the con- tributions made in the paper and important assumptions and limitations. A No NA answer to this question will not be perceived well by the reviewers.
The claims made should match theoretical and experimental results, and reflect how much the results can be expected to generalize to other settings.
It is fine to include aspirational goals as motivation as long is clear that these goals are not attained by the paper.
2. Limitations
Question: Does the paper discuss the limitations of the work performed by the authors?
Answer: [Yes]
Justification:
See Limitations section in the conclusion.
Guidelines:
The answer NA means that the paper has limitation while the answer No means that the paper has limitations, but those are not discussed in the paper.
The authors are encouraged create separate "Limitations" section in their paper.
The paper should point out any strong assumptions and how robust the results are to vi- olations of these assumptions (e.g.. independence assumptions, noiseless settings, model seelfspeciiicatton, asymptotic approximations only holding locally). The authors should reflect on how these assumptions might be violated in practice and what the implications would be.
The authors should reflect on the scope of the claims made, e.g., if the approach was only tested on few datasets or with few runs. In general, empirical results often depend on implicit assumptions, which should be articulated.
The authors should reflect the factors that influence the performance of the approach. For example, a facial recognition algorithm may perform poorly when image resolution is low or images are taken in low lighting. Or a speech-to-text system might not be used reliably to provide closed captions for online lectures because it fails to handle technical jargon.
The authors should discuss the computational efficiency the proposed algorithms and how they scale with dataset size.
If applicable, the authors should discuss possible limitations of their approach to address problems of privacy and faimesss.
While the authors might fear that complete honesty about limitations might be used by reviewers grounds for rejection, worse outcome might be that reviewers discover limitations that aren't I acknowledged in the paper. The authors should use their best judgment and recognize that individual actions in favor of transparency play important role in developing norms that preserve the integrity of the community. Reviewers will be specifically instructed to not penalize honesty concerning limitations.
3. Theory Assumptions and Proofs
Question: For each theoretical result, does the paper provide the full set of assumptions and a complete (and correct) proof?
Answer: [es]]
41

---

Justification:
We provide detailed, complete, and correct proofs of all of our claims in the
Appendix.
Guidelines:
The answer NA means that the paper does not include theoretical results. All the theorems, formulas, and proof' in the paper should be numbered and cross- referenced.
All assumptions should be clearly stated or referenced in the statement of any theorems.
The proofs can either appear in the main paper or the supplemental material, but if they appear in the supplemental material, the authors are encouraged to provide short proof sketch to provide intuition.
Inversely, any informal proof provided in the core of the paper should be complemented by formal proofs provided appendix or supplemental material. Theorems and Lemmas that the proof relies upon should be properly referenced.
4. Experimental Result Reproducibility
Question: Does the paper fully disclose all the information needed to reproduce the main experimental results of the paper to the extent that it affects the main claims and/or conclu- sions of the paper (regardless of whether the code and data are provided or not)?
Answer: [Yes]
Justification: We provide code httprraaaaatttoooooooeerarraaaaa
https///gthab..co//nndeaaaaa
Guidelines:
The answer NA means that the paper does not include experiments.
If the paper includes experiments, No answer this question will not be perceived well by the reviewers: Making the paper reproducible is important, regardless of whether the code and data provided not.
If the contribution is a dataset and/or model, the authors should describe the steps taken to make their results reproducible verifiable.
Depending on the contribution, reproducibility can be accomplished in various ways. For example, if the contribution is novel architecture, describing the architecture fully might suffice, or if the contribution a specific model and empirical evaluation, it may be necessary either make it possible for others replicate the model with the same dataset, or provide access to the model. In general. releasing code and data is often one good way accomplish this, but reproducibility can also be provided via detailed instructions for how to replicate the results, access to hosted model (e.g., in the case of large language model), releasing of a model checkpoint, or other means that are appropriate to the research performed.
While NeurIPS does not require releasing code, the conference does require all submis- sions to provide some reasonable avenue for reproducibility, which may depend on the nature of the contribution. For example
(a) If the contribution is primarily new algorithm, the paper should make it clear how to reproduce that algorithm.
(b) If the contribution is primarily a new model architecture, the paper should describe the architecture clearly and fully.
(c) If the contribution is a new model (e.g., large language model), then there should either be a way to access this model for reproducing the results or a way reproduce the model (e.g., with an open-source dataset or instructions for how to construct the dataset).
(d) We recognize that reproducibility may be tricky in some cases, in which case authors are welcome to describe the particular way they provide for reproducibility. In the case of dosed-source models, may be that access to the model is limited in some way (e.g., to registered users), but it should be possible for other researchers to have some path to reproducing or verifying the results.
5. Open access to data and code
Question: Does the paper provide open access to the data and code, with sufficient instruc- tions to faithfully reproduce the main experimental results, as described in supplemental material?
42

---

Answer: [Yes]
Justification:
We provide code the in supplement.
Guidelines:
The answer NA means that paper does not include experiments requiring code.
Please the NeurlPS code and data submission guidelines (https://nips. cc/ public/ Betee//iidee/odeSumiiiii oolliiccy))))o for more details.
While encourage the release of code and data, understand that this might not be possible, so "No" is an acceptable answer. Papers cannot be rejected simply for not including code, unless this central to the contribution (e.g., for new open-source benchmark).
The instructions should contain the exact command and environment needed to run to reproduce the results. See the NeurIPS code and data submission guidelines (https:
assionoliccy myy for more details.
The authors should provide instructions on data access and preparation, including how to access the raw data, preprocessed data, intermediate data, and generated data, etc.
The authors should provide scripts to reproduce all experimental results for the new pro- posed method and baselines. If only a subset of experiments are reproducible, they should state which ones omitted from the script and why.
At submission time, to preserve anonymity, the authors should release anonymized ver- sions (if applicable).
Providing as much information as possible in supplemental material (appended to the paper) is recommended, but including URLs to data and code is permitted.
6. Experimental Setting/Details
Does the paper specify all the training and test details (c.g.. data splits, hyper- parameters, how they were chosen, type of optimizer, etc.) necessary to understand the results?
Question:
Answer: [Yes]
Justification:
We provide experimental descriptions in the Appendix.
Guidelines:
The answer NA means that the paper does not include experiments.
The experimental setting should be presented in the core of the paper level of detail that necessary to appreciate the results and make sense of them.
The full details can be provided either with the code, in appendix, or supplemental material.
7. Experiment Statistical Significance
Question: Does the paper report error bars suitably and correctly defined or other appropri- ate information about the statistical significance of the experiments?
Answer: [es]
Justification: Our experiments reported in the Appendix.
with 20 random seeds, with the full range being
Guidelines:
The answer NA means that the paper does not include experiments.
The authors should answer "Yes" if the results are accompanied by error bars, confidence intervals, or statistical significance tests, at least for the experiments that support the main
claims of the paper.
The factors of variability that the error bars are capturing should be clearly stated (for example, train/test split, initialization, random drawing of some parameter, overall run with given experimental conditions).
The method for calculating the error bars should be explained (closed form formula, call to a library function, bootstrap, etc.)
The assumptions made should be given (e.g.. Normally distributed errors).
It should be clear whether the error bar the standard deviation or the standard error of the mean.
43

---

It is OK to report 1 1-sigma error bars, but one should state it. The authors should prefer- ably report 2-sigma error bar than state that they have 96% CI, if the hypothesis of Normality of errors is not verified.
For asymmetric distributions, the authors should be careful not to show in tables or fig- ures symmetric error bars that would yield results that are out of range (e.g. negative error rates).
If error bars reported in tables plots. The authors should explain in the text how they were calculated and reference the corresponding figures or tables in the text.
8. Experiments Compute Resources
For each experiment, does the paper provide sufficient information on the com- puter resources (type of compute workers, memory, time of execution) needed to reproduce the experiments?
Question:
Answer: [Yes]
Justification:
We provide experimental details in the Appendix.
Guidelines:
The answer NA means that the paper does not include experiments.
The paper should indicate the type of compute workers CPU or GPU, internal cluster, or cloud provider, including relevant memory and storage.
The paper should provide the amount of compute required for each of the individual experimental runs as well as estimate the total compute.
The paper should disclose whether the full research project required more compute than the experiments reported in the paper (e.g., preliminary or failed experiments that didn't make it into the paper).
9. Code Of Ethics
Question: Does the research conducted in the paper conform, in every respect, with the NeurIPS Code of Ethics https: ftttps://neurip:: cwterripscccppuueee wccppblic//thieeecceeeeeeeceeeee CEthicsuuiidellneeee
Answer: [Yes]
Justification: Our paper is theoretical and doesn't 't involve human subjects or proprietary data.
Guidelines:
The answer NA means that the authors have not reviewed the NeurIPS Code of Ethics.
If the authors answer No, they should explain the special circumstances that require deviation from the Code of Ethics.
The authors should make sure preserve anonymity (e.g., there a special consider- ation due to laws regulations in their jurisdiction).
10. Broader Impacts
Question: Does the paper discuss both potential positive societal impacts and negative societal impacts of the work performed?
Answer: [Yes]
Justification: We discuss the wide-ranging applicability of our theoretical contribution: that understanding the predictability implies parallelizability allows for the parallelizability of nonlinear systems across many domains.
Guidelines:
The answer NA means that there is no societal impact of the work performed.
If the authors answer NA or No, they should explain why their work has no societal
impact why the paper does address societal impact.
Examples of negative societal impacts include potential malicious or unintended uses (e.g., disinformation, generating fake profiles, surveillance). fairness considerations (e.g., deployment of technologies that could make decisions that unfairly impact specific groups), privacy considerations, and security considerations.
The conference expects that many papers will be foundational research and not tied to particular applications, let alone deployments. However, if there is direct path any negative applications, the authors should point it out. For example, it is legitimate to point
44

---

out that an improvement in the quality of generative models could be used to generate deepfakes for disinformation. On the other hand, it is needed point out that generic algorithm for optimizing neural networks could enable people to train models that generate Deepfakes faster.
The authors should consider possible harms that could arise when the technology is being used as intended and functioning correctly, harms that could arise when the technology is being used as intended but gives incorrect results, and harms following from (intentional or unintentional) misuse of the technology.
If there are negative societal impacts, the authors could also discuss possible mitigation strategies (e.g., gated release of models, providing defenses in addition to attacks, mecha- nisms for monitoring misuse, mechanisms to monitor how system learns from feedback over time, improving the efficiency and accessibility of MI.).
11. Safeguards
Does the paper describe safeguards that have been put in place for responsible release of data models that have high risk for misuse (e.g., pretrained language models, image generators, or scraped datasets)?
Question:
Answer: [NA]
Justification:
We do not create any language models or generators or scraped datasets.
Guidelines:
The answer NA means that the paper poses such risks.
Released models that have a high risk for misuse or dual-use should be released with necessary safeguards to allow for controlled use of the model, for example by requiring that users adhere to usage guidelines or restrictions to access the model or implementing safety filters.
Datasets that have been scraped from the Internet could pose safety risks. The authors should describe how they avoided releasing unsafe images. We recognize that providing effective safeguards is challenging, and many papers do not require this, but we encourage authors to take this into account and make best faith effort.
12. Licenses for existing assets
Question: Are the creators or original owners of assets (e.g., code, data, models), used in the paper, properly credited and are the license and terms of use explicitly mentioned and properly respected?
Answer: [Yes]
Justification:
We cite original creators, such as the creator of the dysts benchmark dataset.
Guidelines:
The answer NA means that the paper does not use existing assets.
The authors should cite the original paper that produced the code package or dataset.
The authors should state which version of the asset is used and, if possible, include URL.
The name of the license (e.g., CC-BY 4.0) should be included for each asset.
For scraped data from a particular source (e.g., website), the copyright and terms of service of that source should be provided.
If assets are released, the license, copyright information, and terms of use in the pack- age should be provided. For popular datasets, paperswi saperswithcon thcode com/dataatte datasets has curated licenses for some datasets. Their licensing guide can help determine the license of a dataset.
For existing datasets that are re-packaeed both the original license and the license of the derived asset (if it has changed) should be provided. If this information is not available online, the authors are encouraged to reach out to the asset's creators.
13. New Assets
Question: Are new assets introduced in the paper well documented and the documenta- tion provided alongside the assets?
45

---

Answer: [Yes]
Justification: We provide our code and documentation for it at htt//gtttubb.oom com/ www.beess.aaaallliiaaioo..
Guidelines:
The answer NA means that the paper does not release new assets.
Researchers should communicate the details the dataset/code/model as part of their submissions via structured templates. This includes details about training, license, limi- tations, etc.
The paper should discuss whether and how consent was obtained from people whose asset is used.
At submission time, remember to anonymize your assets (if applicable). You either create an anonymized URL or include an anonymized zip file.
14. Crowdsourcing and Research with Human Subjects
Question: For crowdsourcing experiments and research with human subjects, does the pa- per include the full text of instructions given participants and sereenshots, applicable, well as details about compensation (if any)?
Answer: [NA]
Justification: Our paper does not involve research with human subjects.
Guidelines:
The answer NA means that the paper does not involve crowdsourcing nor research with human subjects.
Including this information in the supplemental material is fine, but if the main contri- bution of the paper involves human subjects, then as much detail as possible should be included in the main paper.
According the NeurIPS Code of Ethies, workers involved in data collection, curation, or other labor should be paid at least the minimum wage in the country of the data col- lector.
15. Institutional Review Board (IRB) Approvals or Equivalent for Research with Human Subjects
Question: Does the paper describe potential risks incurred by study participants, whether such risks were disclosed to the subjects, and whether Institutional Review Board (IRB) approvals (or equivalent approval/review based the requirements of your country or institution) were obtained?
Answer: [NA]
Justification: the paper does not involve crowdsourcing nor research with human subjects.
Guidelines:
The answer NA means that the paper does not involve crowdsoureing nor research with human subjects.
Depending on the country in which research is conducted, IRB approval (or equivalent)
may be required for any human subjects research. If you obtained IRB approval, you should clearly state this in the paper.
We recognize that the procedures for this may vary significantly between institutions and locations, and we expect authors to adhere to the NeurIPS Code of Ethics and the guidelines for their institution.
For initial submissions, do not include any information that would break anonymity (if applicable), such as the institution conducting the review.
16. Declaration of LLM usage
Question: Does the paper describe the usage of LLMs if it is an important, original, or non-standard component of the core methods in this research? Note that if the LLM is used only for writing, editing, or formatting purposes and does not impact the core methodology, scientific rigorousness, or originality of the research, declaration is not required.
Answer: [NA]
46

---

Justification:
The core methods introduced in this paper do not involve LLMs as any im- portant, original, or non-standard components.
Guidelines:
The answer NA means that the core method development in this research does not involve
LLMs as any important, original, or non-standard components.
for Please refer to our LLM policy httpg:/ss/aaauuiee.cc//o/ooeeees what should or should not be described.
47


Reviews:
[{"review_id": "DNUX6D9lSL", "reviewer": "Reviewer_HYQS", "summary": "The paper presents a theoretical framework demonstrating that a system’s predictability governs the number of optimization steps needed for evaluation. This relationship is important for assessing the efficiency of parallel computing. The authors validate their approach by analyzing the convergence behavior of the Gauss-Newton DEER algorithm using a worst-case bound on optimization steps.", "questions": "- How does the degree of predictability (i.e., sensitivity to perturbations) affect the results?\n\n- How does the system’s dimension impact the results?\n\n- Could you clarify the example in lines 177–182? You state residuals for later times are assumed to be zero, but later mention residuals grow over time. What exactly is assumed in that example?", "limitations": "yes", "rating": 4, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 3, "significance": 3, "originality": 2, "strengths_and_weaknesses": "## Strengths\n- Clear presentation of the problem, motivation, and overall approach (though the organization could be improved, see below).\n- The motivation behind the theory is well explained.\n- The authors provided proofs where appropriate.\n- Demonstrated the model using the DEER algorithm.\n\n## Weaknesses\n- The evaluation is limited to the DEER algorithm. Including comparisons to additional algorithms will help determine if the theory applies only to DEER or to a wider range of algorithms.\n- It would be helpful to include some comparisons to other methods; i.e., what theoretical guarantees exist about the link between certain system quantities and the corresponding optimization landscape. I see you mentioned in lines 272–273 that there are no prior works on LLE, but I assume there may be some works that link it to other system quantities. Without such comparisons, it’s hard to evaluate the strength of the results.\n- The paper feels organized weirdly. Putting the related work section in the middle throws off the flow and makes it harder to follow.", "was_revised": true, "final_justification": "I thank the authors for their response and believe that the clarifications, additional comparisons, and discussion points will strengthen the paper. I increased my score to 4.", "rebuttal": "We thank the reviewer for their constructive comments. We address each point, and have incorporated all changes in the revised manuscript. We believe we have significantly strengthened our submission in doing so.\n\n## Generality beyond DEER\n\n> will help determine if the theory applies only to DEER or to a wider range of algorithms\n\nWe thank the reviewer for emphasizing the importance of highlighting the generality of our results. \n\nTheorem 2 establishes a relationship between the conditioning of the merit function and the Largest Lyapunov Exponent (LLE) $\\lambda$. Specifically, we show that the Polyak–Łojasiewicz (PL) constant $\\mu$ of the merit function increases monotonically with $\\lambda$. The PL constant $\\mu$ quantifies the conditioning of the optimization problem—i.e., how \"well-conditioned\" the merit function landscape is—while the LLE $\\lambda$ measures the predictability of the underlying dynamical system.\n\nTheorem 2 thus demonstrates that systems with low predictability ($\\lambda > 0$) exhibit poorly conditioned merit function landscapes, whereas more predictable systems ($\\lambda < 0$) yield well-conditioned landscapes.\n\nImportantly, Theorem 2 concerns *the landscape of the merit function itself*, independent of any specific optimization algorithm. Therefore, our theory is not limited to Newton’s method (DEER), but extends to a broad class of optimization methods.\n\nFor instance, as discussed on line 131, the PL condition guarantees that gradient descent achieves a linear convergence rate. However, the rate constant depends on the value of $\\mu$, which we becomes prohibitively small in the case of unpredictable systems—leading to slow convergence. This highlights that our theoretical results are not limited to Gauss-Newton method (DEER) but also extend to gradient descent. We have clarified this point in the revised manuscript to better emphasize the broad applicability of our contributions across optimization algorithms.\n\n## The degree of predictability has a monotonic relationship with the conditioning of the merit function landscape\n\n> How does the degree of predictability (i.e., sensitivity to perturbations) affect the results?\n\nTheorem 2 (equation 9) precisely indicates how the degree of predictability $\\lambda$ influences the difficulty of the optimization problem, as measured by $\\mu$. The relationship is monotonic: the higher the degree of unpredictability, the more difficult the corresponding optimization problem. \n\nIn the paper we particularly emphasized the sharp change that occurs at the transition between predictability and unpredictability (at $\\lambda=0$). We did so because the $\\exp(\\lambda T)$ in the denominator of equation 9 grows exponentially for positive $\\lambda$. We have revised our paper to more heavily emphasize that the degree of predictability affects the difficulty of the optimization problem, and that the particular functional form of this relationship is what makes the predictability/unpredictability transition so stark.\n\nThe effect of the degree of predictability also becomes more clear in our experimental evaluations with other methods that we have added, thanks to the reviewer’s feedback.\n\n## Strengthened Experimental Evaluation\n\nWe have strengthened our experimental evaluations by including results for both **gradient descent** and a **quasi-Newton method**. Specifically, we revisited our experiment involving the mean-field RNN (see Figure 2 in the paper). In the setting depicted in the rightmost plot of this figure (state dimension $D = 100$, sequence length $T = 954$), we also measured the rates of convergence of a quasi-Newton method and gradient descent. The results are summarized in the following table:\n\n###  Iterations until convergence\n| LLE               |   Newton |   quasi-Newton |   gradient descent |\n|:------------------|---------:|---------------:|-------------------:|\n| -0.7              |        3 |             10 |                 31 |\n| -0.5              |        4 |             13 |                 58 |\n| -0.35             |        4 |             17 |                206 |\n| -0.2              |        4 |             26 |                435 |\n| -0.1              |        4 |             66 |                954 |\n| -0.1 to 0         |       18 |            674 |                954 |\n| 0 to 0.05         |      147 |            954 |                954 |\n| 0.05 to 0.1       |      642 |            954 |                954 |\n| 0.1 to 0.14       |      895 |            954 |                954 |\n| greater than 0.14 |      954 |            954 |                954 |\n\n\nThis table clarifies two points raised in the review:\n\n1. The effect of the degree of predictability (the LLE) on rate of convergence.\n\nFor all 3 optimizers, the number of steps to convergence increases monotonically with the value of the LLE. This is the case because, as we proved in Theorem 2, the conditioning $\\mu$ of the merit function landscape degrades monotonically as the LLE $\\lambda$ increases. \n\n2. The applicability of our theory across optimizers.\n\nWe observe that this monotonic relationship between LLE and numbers of steps to converge holds across optimizers, thus strengthening our evaluation as per the reviewer’s recommendation. Newton’s method converges in a small number of steps all the way up to the threshold between predictability and unpredictability (LLE= 0). Intuitively, the performance of the other optimizers degrades more quickly as unpredictability increases because quasi-Newton and gradient descent use less information about the curvature of the loss landscape. \n\n## Comparisons with Related Work\n\n>  the related works section\n\nWe agree about the location of related works. We have moved it to the conclusion to improve the flow.\n\n> comparison to other methods\n\nTo the best of our knowledge, our work is novel in connecting the LLE of a dynamical system with the conditioning of its merit function landscape. Nonetheless, we thank the reviewer for highlighting the need to expand the comparison of our contributions with the literature. We have done so in our revision, better contextualizing the strength of our results.\n\nOther system-theoretic quantities have been proposed to assess the “parallelizability” of dynamical systems. Particularly relevant to our work is the study we cited by Chartier and Philippe (reference [21] in the paper), which investigates the parallelizability of input-free continuous-time systems using a multiple-shooting approach. Their analysis centers on the *logarithmic norm* of the system Jacobian, a quantity that can be roughly interpreted as a “one-step” analogue of the LLE. Specifically, while the LLE captures the *asymptotic* growth rate of perturbations over time, the logarithmic norm quantifies the *instantaneous* contractivity of the system’s Jacobian. As such, a negative logarithmic norm implies a negative LLE, but the converse does not necessarily hold. In this way, our work is a significant generalization of Chartier and Philippe.  \n\nWe also surveyed literature in theoretical computer science relating to the parallelization of nonlinear recursions. In particular, we cited Hyafil and Kung (ref [17]) and Kung (ref [18]). They showed that when the nonlinear recursion is a rational function (i.e., the ratio of two polynomials), and if the degree of that rational function is greater than one, then such a nonlinear recursion does not benefit from parallelization as the number of processors grows (cf. Theorem 4 of [17]). However, our work greatly generalizes this foundational contribution, because our formulation in terms of the LLE applies to *any* recursion, not just rational functions. \n\n## Additional clarifications\n\n> How does the system’s dimension impact the results?\n\nTheorem 2 does not explicitly depend on the system dimension $D$, just the LLE $\\lambda$. However, the LLE itself can depend on the system dimension. For example, consider a linear dynamical system with transition matrix $A$. If the entries of $A$ are drawn i.i.d from a normal distribution with mean 0 and variance $g$, then classical results from random matrix theory (e.g., Sompolinsky, et al. \"Chaos in random neural networks”) show that, as the dimension $D$ grows, the largest eigenvalue scales like $g\\sqrt{D}$. Consequently, if $g$ is not appropriately scaled with dimension, the LLE may increase and eventually become positive as $D$ increases.\n> the example in lines 177–182\n\nIn this example, we provide intuition as to why unpredictable systems have flat merit landscapes. We used two different ideas: the *residuals* and the *errors*.\n\nWe define the residual $\\mathbf{r}(\\mathbf{s})$ in equation 2. The residual at time $t$ is $s_t - f(s_{t-1})$. The residual is a function only of our current guess for the sequence trajectory $\\mathbf{s}$.\n\nWe have clarified the definition of “errors” in our revision. The error is the distance between our current guess for the sequence $\\mathbf{s}$ and the true sequence trajectory $\\mathbf{s*}$. More precisely, we can define the error at time $t$ as $e_t(s_t, s_t*) = | s_t - s_t* |$. In contrast to the residual, the error is a function both of our current guess $\\mathbf{s}$ and of the true trajectory $\\mathbf{s*}$. \n\nThis example shows that for an unpredictable system, small residuals do *not* necessarily imply small errors. \n\nIn this example, we do assume that the residual for all times greater than 1 is 0, that is $s_t = f(s_{t-1})$ for $t > 1$. However, we do not state that the “residuals grow over time.” Instead, we note that the *error* (that is, the difference between $s_t*$ and $s_t$) grows over time (see lines 180-182). This exponential growth in the error—in spite of the small and bounded value of the residual—highlights why many trajectories can have small residual for unpredictable systems, even though they may be very different from the true trajectory $\\mathbf{s*}$. This is why unpredictable systems have flat merit landscapes."}, {"review_id": "1LF4KXIsSi", "reviewer": "Reviewer_XTsi", "summary": "The paper investigates the condition in which nonlinear state space models can be efficiently\nparallelized. Although recent methods like DEER (non-linear Differential Equation as fixed point\nitERation) where the model was able to apply parallel processes to efficiently train sequential\nmodels, what factors determine the convergence speed of this optimization remains unclear.\nTherefore, the authors demonstrate that a predictable system, measured by the Largest\nLyapunov Exponent, corresponds to well-conditional optimization problems that can be solved in\nparallel process, wherares unpredictable systems lead to ill-conditioned problems that converge\ntoo slowly.", "questions": "NA", "limitations": "Yes", "rating": 5, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 4, "significance": 3, "originality": 3, "strengths_and_weaknesses": "The paper’s strength lies in its main contribution of connecting the Largest Lyapunov Exponent\nto PL conditioning in optimization problems. The formalization of predictable and unpredictable\nsystems was also creative and paved the way for the author’s main work. The proof was also\nclear, and the author used a strong intuition, such as the \"flatness\" of the merit landscape for\nchaotic systems, to explain complex mathematical results. The work has practical significance\nfor state-space models and recurrent neural networks to operate on parallel hardware.\n\nHowever, the paper’s theoretical convergence results and analysis is closely tight with the\nDEER algorithm. The authors did not fully explore other methods. \nAdditionally, the work has a niche relevance since a lot of the deep learning community has been shifted to inherently\nparallel architectures like Transformers.", "was_revised": false, "final_justification": "", "rebuttal": "Firstly, we thank the reviewer for taking the time to review our submission and for their positive comments! We were particularly pleased by the comments on the creativity, clarity, and practicality of our work. \n\nThe reviewer provided two main areas for improvement: an exploration of other optimization algorithms, and a better discussion of the relevance of our theory in the context of transformers’ popularity. By incorporating this feedback, we believe that we have strengthened our submission. We discuss in detail below.\n\n## Our theory also applies to other optimizers like gradient descent and quasi-Newton methods\n\nWe agree that it is important to emphasize the broad applicability of our theoretical contribution to optimization techniques beyond the Gauss-Newton method (DEER).\n\nTo recapitulate, Theorem 2 establishes a relationship between the conditioning of the merit function and the Largest Lyapunov Exponent (LLE) $\\lambda$. Specifically, we show that the Polyak–Łojasiewicz (PL) constant $\\mu$ of the merit function increases monotonically with $\\lambda$. The PL constant $\\mu$ quantifies the conditioning of the optimization problem—i.e., how \"well-conditioned\" the merit function landscape is—while the LLE $\\lambda$ measures the predictability of the underlying dynamical system.\nTheorem 2 thus demonstrates that systems with low predictability (high $\\lambda$) exhibit poorly conditioned merit function landscapes, whereas more predictable systems (low $\\lambda$) yield well-conditioned landscapes. Importantly, this result concerns *the landscape of the merit function itself*, independent of any specific optimization algorithm. Therefore, our theory is not limited to the Gauss-Newton method (DEER), but extends to a broad class of optimization methods and settings.\nFor instance, as discussed on line 131, the PL condition guarantees that gradient descent achieves a linear convergence rate on the merit function landscape. However, the rate constant depends on the value of $\\mu$, which we observe to become prohibitively small in the case of unpredictable systems—leading to extremely slow overall convergence. This highlights that our theoretical results are not limited to Gauss-Newton method (DEER) but also extend to gradient descent. We have clarified and expanded on this point in the revised manuscript to better emphasize the broad applicability of our theoretical contributions across optimization algorithms.\nWe have also strengthened our experimental evaluation by including results for both **gradient descent** and a **quasi-Newton method**. Specifically, we revisited our first experiment involving the mean-field RNN (see Figure 2 in the paper). In the setting depicted in the rightmost plot of this figure (state dimension $D = 100$, sequence length $T = 954$), we conducted an additional experiment comparing the convergence rates of Newton’s method, a quasi-Newton method, and gradient descent. The results are summarized in the following table (unfortunately, images are not permitted in this year’s NeurIPS rebuttals):\n\n###  Iterations until convergence\n| LLE               |   Newton |   quasi-Newton |   gradient descent |\n|:------------------|---------:|---------------:|-------------------:|\n| -0.7              |        3 |             10 |                 31 |\n| -0.5              |        4 |             13 |                 58 |\n| -0.35             |        4 |             17 |                206 |\n| -0.2              |        4 |             26 |                435 |\n| -0.1              |        4 |             66 |                954 |\n| -0.1 to 0         |       18 |            674 |                954 |\n| 0 to 0.05         |      147 |            954 |                954 |\n| 0.05 to 0.1       |      642 |            954 |                954 |\n| 0.1 to 0.14       |      895 |            954 |                954 |\n| greater than 0.14 |      954 |            954 |                954 |\n\nFor all three optimizers, we observe that the number of steps to convergence increases monotonically with the value of the LLE. This occurs because, as we proved in Theorem 2, the conditioning $\\mu$ of the merit function landscape degrades monotonically as the LLE $\\lambda$ increases.\n\n## Parallelizing sequential operations is broadly important\n\nWe thank the reviewer for emphasizing the need to highlight the importance and relevance of our work in the context of transformers. While the original DEER paper (Lim et al) focused exclusively on RNNs and Neural ODEs for sequence modeling, we should have emphasized that DEER can apply to any Markovian system. In our revisions, we have more heavily emphasized the importance of parallelizing sequential computation generally, and not just in sequence modeling. The below text indicates our reframing of the contextualization:\n\nMarkovian systems involving sequential computation abound across deep learning, probabilistic modeling, and scientific computing. In particular:\n\n### 1. Parallelization over *depth* is relevant for transformers\n\nWe agree that transformers are an extremely important architecture, and that the attention mechanism is parallelized *over the sequence length*. However, almost all transformers consist of multiple transformer blocks that must be applied sequentially *in depth*. Therefore, sequential operations also occur in transformers. As transformers become increasingly deep, parallelizing Transformers over depth will become even more important. In fact, preliminary efforts to this effect have already been attempted (Calvo Gonzalez et al, “Leveraging the true depths of LLMs”, ‘25). Another important trend is “recurrent depth” in transformers. “Recurrent depth” treats test-time-compute/latent space “thinking” as an RNN where the time step function is a transformer block (see Schone et al, “Implicit Language Models are LLMs”, ICML ‘25 Spotlight;  and Geipeng et al, “Scaling up Test-Time Compute with Latent Reasoning: A Recurrent Depth Approach”, ‘25). We believe DEER could play a key role in enabling efficient, complete parallelization for such models, especially those that involve iterative or recurrent computation over depth.\n\n### 2. Stateful models are relevant across ML\n\nA huge range of important models in deep learning and machine learning involve stateful computation. These include the sampling pass of a diffusion model, the forward pass of a feedforward network, and the state of a reinforcement learning agent. Moreover, stateful operations are important in probabilistic modelling, such as Markov Chain Monte Carlo (MCMC). In fact, our results in Figure 3 about the parallelizability of Langevin dynamics is highly suggestive about applications to MCMC because MALA (Metropolis-adjusted Langevin algorithm) is an important MCMC workhorse. We have better emphasized this connection in our revision.\n\n### 3. Stateful models are relevant in scientific computing\n\nFinally, the ability to parallelize seemingly sequential computation is important in scientific computing broadly. Notable examples include the solving of differential equations and computing the trajectories of systems under various forms of control. In fact, our results in Table 1 about the difficulties of parallelizing chaotic observers, and the ability to parallelize stable observers, relate to both these fields (solving ODEs and controlling systems). We have more heavily emphasized these connections in our revision.\n\n### In summary\n\nWe have expanded our introduction to appropriately motivate the importance of sequential operations as a primitive across machine learning and scientific computing. We situate our contribution as providing an important theoretical basis for determining in which settings sequential computation can be efficiently parallelized."}, {"review_id": "oXGCZF9Yhi", "reviewer": "Reviewer_XHnm", "summary": "The paper establishes a direct link between how predictable the system is, measured by its largest Lyapunov exponent, and the minimum singular value of the merit function, which is a set of sequential estimation errors. Specifically, the authors show that if small perturbations are diminished over time, which is defined by the largest Lyapunov exponent, the authors also show that the merit function is well-conditioned, which allows for parallel evaluation and utilizing algorithms such as DEER. On the contrary, if perturbations are not diminished over time, the authors also show that it will take an exponential number of steps for convergence. The authors showcase the theoretical results in mean-field RNNs and Langevin dynamics to show how DEER scales with stable and chaotic systems.", "questions": "The theoretical conditions developed by the paper are well-defined and well-presented. However, I would like to ask some questions about the algorithm.\n\nThis does not affect the theoretical results, but I was curious whether there are existing bounds between the finite-time estimates of the Lyapunov exponent versus the infinite-horizon Lyapunov exponent used in the theoretical results.\n\nI believe the authors require a global Lipschitz condition to hold for global convergence of the results, which is defined in Theorem 4. I wonder if it's possible to establish similar results with local Lipschitz assumptions while understanding that the convergence would depend on the initialization.\n\nIf the model is only locally Lipschitz, does this imply that the Lyapunov exponent is positive, and hence, it's computationally expensive to perform parallel estimation? The authors consider systems that satisfy Lipschitz conditions globally, but I was curious if these conditions do not hold with locally Lipschitz systems, which is very common in non-linear systems.", "limitations": "yes", "rating": 5, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 4, "clarity": 4, "significance": 4, "originality": 3, "strengths_and_weaknesses": "The paper provides a well-defined notion of stability of the non-linear systems and the conditioning of the merit functions, which provides a clear signal on the parallelizability of the inference of non-linear systems, such as using algorithms like DEER.\n\nThe only weakness that I can find is, the paper assumes that the non-linear systems satisfy global Lipschitz conditions, which may not hold for non-linear systems that often appear in practice. However, the global Lipschitz conditions may be necessary for parallel inference with the metric functions, and the authors may clarify this in the rebuttal phase.", "was_revised": true, "final_justification": "I think the paper can be accepted for publication. The authors addressed my comments and the other reviewers comments well, and I initially recommended the acceptance of the paper.", "rebuttal": "We thank the reviewer for their thoughtful and positive review, and for their interesting questions about Lipschitzness and the Largest Lyapunov Exponent. We discuss the points they raise below.\n\n> The only weakness that I can find is, the paper assumes that the non-linear systems satisfy global Lipschitz conditions, which may not hold for non-linear systems that often appear in practice.\n\nMany of the results in our paper do not depend on Lipschitz conditions (Theorem 1, Theorem 2, and Theorem 4). In particular, we emphasize that our main result, Theorem 2—which establishes a direct link between how predictable the system is, measured by its largest Lyapunov exponent, and the minimum singular value of the merit function Jacobian—does not have Lipschitzness as an assumption. \n\nHowever, the reviewer is correct that Theorem 5 does require a Lipschitz assumption. We agree that an interesting direction for future work could be to investigate if Lipchitzness conditions for this theorem could be relaxed.\n\nFinally, while we agree that there are important non-linear systems that are not Lipschitz, we note that many systems in deep learning are Lipschitz by design. In particular, if the transition function is some combination of multiplication by a finite dimensional matrix and the application of a Lipschitz nonlinearity (such as a ReLU or saturating nonlinearities), then that transition function is also Lipschitz (cf. eg. Szegedy et al, “Intriguing properties of neural networks”). \n\n> This does not affect the theoretical results, but I was curious whether there are existing bounds between the finite-time estimates of the Lyapunov exponent versus the infinite-horizon Lyapunov exponent used in the theoretical results.\n\nWe thank the reviewer for raising this interesting question. \n\nThis topic has been addressed in the literature. For example, see Sections 4 (\"Error Analysis\") and 5 of [1] below, where bounds are derived between the finite-time Lyapunov exponents and their infinite-time limits. These results apply to all Lyapunov exponents, not just the largest one.\n\nAn alternative approach to this question is to consider it probabilistically. In particular, where one can assert—with high probability—that the finite-time estimate deviates from the infinite-time exponent by no more than a specified amount. These bounds are usually found in the context of ergodic or random dynamical systems and are framed using large deviation principles. See reference [2] below and references therein for examples of this second approach. \n\n[1] Dieci, L., Russell, R. D., & Van Vleck, E. S. (1997). On the computation of Lyapunov exponents for continuous dynamical systems. SIAM journal on numerical analysis, 34(1), 402-423.\n\n[2] Prasad, A., & Ramaswamy, R. (1999). Characteristic distributions of finite-time Lyapunov exponents. Physical Review E, 60(3), 2761.\n\n> However, the global Lipschitz conditions may be necessary for parallel inference with the metric functions, and the authors may clarify this in the rebuttal phase.\nand\n> I believe the authors require a global Lipschitz condition to hold for global convergence of the results, which is defined in Theorem 4. I wonder if it's possible to establish similar results with local Lipschitz assumptions while understanding that the convergence would depend on the initialization.\n\nWe thank the reviewer for the opportunity to clarify this point. Importantly, the proof of Theorem 4 does not require a global Lipschitz condition. However, a global Lipschitz condition is used elsewhere in the paper (e.g., in Theorem 5). We have revised the manuscript to clarify this distinction.\nWe agree that relaxing the global Lipschitz requirement in favor of a local one, with convergence guarantees that depend explicitly on the initialization, is an important (and promising!) direction. While we have not pursued this in the current work, we are optimistic that our results could be extended in this way. \n> If the model is only locally Lipschitz, does this imply that the Lyapunov exponent is positive, and hence, it's computationally expensive to perform parallel estimation? The authors consider systems that satisfy Lipschitz conditions globally, but I was curious if these conditions do not hold with locally Lipschitz systems, which is very common in non-linear systems.\nWe thank the reviewer for this interesting question. The Lipschitzness of a model and its Lyapunov exponents can in fact be decoupled, in the sense that a model only being locally Lipschitz does not imply a positive Lyapunov exponent, and vice versa.\n\nThis can be seen with an example. Consider the system:\n\n$x_t = (1-h) x_{t-1} - h x^3_{t-1}$.\n\nWhere $0 < h \\ll 1$. Note that this system is the Euler discretization of the stable, continuous-time dynamics $dx/dt = -x - x^3$, with Euler step-size h. \n\nThis system is not globally Lipschitz, because its Jacobian is \n$J(x_{t-1}) = (1- h) - 2 h x^2_{t-1}$\n\nSince $J’(x) = -4hx$ can become unbounded with increasing $x$, the Jacobian is not Lipschitz. \n\nLet us now evaluate the LLE around the origin. When $x = 0$, we have that \n\n$|J(0)| = |1- h| < 1 $.\n\nThus, the origin is a locally stable trajectory and thus has negative LLE. Indeed, for this system, one can show that if the initial condition $x_0$ is drawn from a finite ball of radius $R$ around the origin, then there exists a sufficiently small $h(R)$ such that: 1) the ball is forward-invariant—that is, all trajectories remain within it for all time; and 2) the Jacobian along all trajectories has absolute value less than one, implying that the system is contracting and the largest Lyapunov exponent (LLE) is negative for all trajectories.\nThus, our results (in particular equation 9 and Theorem 4) hold for systems which are not globally Lipschitz."}, {"review_id": "sNf3Pyymer", "reviewer": "Reviewer_HUX6", "summary": "The paper presents a theoretical framework relating the convergence speed of parallel evaluation in nonlinear state space models to the predictability of the system dynamics through the conditioning of a merit function. The authors support the theory with several experiments showing strong agreement between empirical results and theoretical predictions.", "questions": "1. What are the axis labels in Figure 1? What do the colors of the contour plot mean? Its unclear what is being shown in each panel.\n2. In the experiments, it is only shown the number of steps to converge for parallel evaluation. How do these results compare to the speed of sequential evaluation? \n3. What is the cost of computing the LLE? What is the memory footprint and runtime complexity? Is it cheap enough to compute during a training loop when training an RNN?\n4. Is the magnitude of the LLE important for interpretation? It seems from Table 1 that the same number of newton steps required can correspond to very different LLE scores. \n5. How can this theory be exploited to improve the speed of training nonlinear state space models?", "limitations": "yes", "rating": 5, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 4, "significance": 3, "originality": 3, "strengths_and_weaknesses": "Strength\n1. Analytic results are theoretically sound and a very nice contribution to the current literature on parallelizing nonlinear SSMs.\n2. Empirical experiments are well thought out and directly validate the proposed theory.\n3. The use tools from dynamical systems to characterize parallel evaluation of RNNs is a novel contribution.\n4. The text is well-written and explanations are intuitive and easy to follow.\n\n\nWeakness\n1. **No comparison with sequential speed.** The experiments report only the number of parallel steps to convergence, without comparing against sequential evaluation. It is unclear whether the theory offers a practical criterion for deciding between sequential and parallel evaluation. Reporting wall-clock times for sequential and parallel approaches for each experiment would help clarify this connection between LLE and the decision to use parallelism.\n2. **No benchmarking of LLE computation cost.** While the LLE is proposed as a criterion for when parallel evaluation is feasible, the practicality of computing the LLE itself is unclear. What is the computational complexity of estimating the LLE, and can it be efficiently computed in the context of large-scale models?\n3. **Unclear connection to RNN training efficiency.** Although the experiments empirically validate aspects of the theory, it remains unclear how this theory can be used to accelerate training of nonlinear sequence models such as RNNs. I would like to see either a concrete demonstration of improved training speed using the proposed approach, or at least a discussion/outline of how the theory could inform better training strategies for recurrent models.", "was_revised": true, "final_justification": "The authors provide additional runtime experiments and a more detailed discussion on the cost and practicality of computing LLE, which should be included in the main paper or supplementary material. I also appreciate the discussion on how the proposed theory can be applied to improve RNN training efficiency. My questions and concerns have been fully addressed, and I have increased my score to accept.", "rebuttal": "We thank the reviewer for taking the time to review our paper. We are especially pleased by the comments that our contribution is “very nice” and “novel”; our analysis is \"theoretically sound”;  our experiments are “well thought out” and “directly validate” our theory; and our text is “well-written”. We also appreciate the thoughtful feedback, which we have addressed and incorporated in our revision.\n\n## Comparison with sequential speed\n\nWe have added wall-clock times for both parallel and sequential evaluation in our experiments. Our experiments use an H100 with 80GB onboard memory, and we average our results over five runs.\n\nFor example, in the setting of Experiment 1 (Figure 2, right) in our paper, we have added wallclock time comparisons, as follows.\n\n### Wall‑clock time (ms) by LLE for Newton & sequential evaluation\n| LLE               |   Newton |   sequential evaluation |\n|:------------------|---------:|------------------------:|\n| -0.7              |     1.22 |                    9.56 |\n| -0.5              |     1.59 |                    9.65 |\n| -0.35             |     1.59 |                    9.17 |\n| -0.2              |     1.61 |                    9.34 |\n| -0.1              |     1.73 |                    9.46 |\n| -0.1 to 0         |     6.84 |                    9.82 |\n| 0 to 0.05         |    53.59 |                    9.92 |\n| 0.05 to 0.1       |   233.33 |                    9.97 |\n| 0.1 to 0.14       |   325.71 |                   10.08 |\n| greater than 0.14 |   345.47 |                   10.09 |\n\nWe observe that in the predictable setting, Newton is an order of magnitude faster than sequential evaluation, while in the unpredictable regime, Newton is 1-2 orders of magnitude slower than sequential evaluation. This importance of using parallel evaluation only in predictable settings is a core practical takeaway from our theoretical contributions.\n\n## Discussion of Largest Lyapunov Exponent\n\n### LLE computation cost\n\nWe use a standard algorithm (Algorithm 1 in Appendix J.5) for computing the largest Lyapunov exponent (LLE). This algorithm propagates a unit vector through a sequence of $T$ Jacobians, normalizing at each step for numerical stability of the estimation (see e.g., Section 3.1 of Politi, A., & Pikovsky, A. (2010). *Lyapunov Exponents: A Tool to Explore Complex Dynamics.*).\n\nAlgorithm 1 has a memory footprint of only $O(D^2)$. We can also use Jacobian vector products (JVPs) to reduce the memory footprint to $O(D)$. Using these JVPs, the computational complexity is $O(TD)$. \n\n### Can LLE be computed and used during training? \n\nThe LLE is indeed cheap enough to compute during an RNN training loop, as demonstrated in [1] where _all_ Lyapunov exponents—not just the top one—are computed during training and used to regularize the RNN away from chaotic solutions. \n\nFinally, we emphasize that in many scenarios, computing the LLE during training is unnecessary. For instance, if a system is known or constrained to be contracting/stable *a priori*, our results imply that it is parallelizable regardless of the exact LLE value. This situation arises, for example, in modern state space models and many RNN architectures (see [2-3] below). \n\n[1] Engelken, Rainer. \"Gradient flossing: Improving gradient descent through dynamic control of Jacobians.\" Advances in Neural Information Processing Systems 36 (2023): 10412-10439.\n\n[2] Gu, Albert, and Tri Dao. \"Mamba: Linear-time sequence modeling with selective state spaces.\" arXiv:2312.00752 (2023).\n\n[3] Krotov, D. A new frontier for Hopfield networks. Nat Rev Phys 5, 366–367 (2023).\n\n### Interpreting the magnitude of the LLE\nWe thank the reviewer for raising this interesting and subtle point. In addition to the Lyapunov Exponent (LLE), our theory (specifically equation 9) shows that the “overshoot” parameter $a$ also plays a role in determining the flatness of the merit function. In principle, two systems could have the same LLE value but different overshoot parameters, resulting in different numbers of optimization steps. We focus on LLE rather than overshoot because the dependence on LLE is exponential, whereas the dependence on overshoot is only linear. Therefore, for large sequences, it is primarily the LLE that determines whether parallelization via merit function optimization is feasible. We have revised the text to better highlight this important point.\n## Application to RNN training efficiency\n\nWe thank the reviewer for the insightful question. Below, we highlight a practical strategy for speeding up the training of nonlinear state space models (SSMs) based on our theoretical findings.\n\nOur results indicate that nonlinear SSMs with negative largest Lyapunov exponents (LLEs) are efficiently parallelizable. To exploit this during training, one must ensure that the model maintains negative LLEs throughout optimization. One straightforward and effective method to achieve this is by design, through **parameterization**. In particular, by introducing an auxiliary variable to enforce the desired constraint (in this case, negative LLE), and then performing unconstrained optimization on this variable.\n\nThis strategy is particularly well-suited to neural network-based SSMs. For example, consider the scalar nonlinear SSM:\n\n$x_t = \\tanh(w x_{t-1} + u_t).$\n\nTo guarantee negative LLE, it suffices to ensure that the Jacobian norm is strictly less than one:\n\n$ |J_t| = |w \\cdot \\text{sech}^2(w x_{t-1} + u_t)| \\leq |w|.$\n\nThus, enforcing $|w| < 1$ is sufficient. This can be achieved by reparameterizing $w = \\tanh(b)$, where $b$ is a trainable, unconstrained auxiliary variable. This guarantees that $w \\in (-1, 1)$ for all finite $b$, ensuring contractivity and, hence, negative LLE.\n\nThis approach has precedent in the literature on training stable RNNs [e.g., Rainer et al., Manchester et al., and Farsang et al.]. We particularly highlight the concurrent work of Farsang et al. (arXiv:2505.21717), which scales up a nonlinear SSM using parallel training with DEER. Farsang et al explicitly parameterizes their SSM to be contractive (see their Appendix A.1). However, they do not make the connection that this contractivity is precisely why DEER enjoys fast convergence throughout training. Ensuring negative LLE via parameterization enables models to remain parallelizable throughout training, providing a pathway to faster and more scalable learning.\n\nWe have more heavily emphasized how enforcing stability in the design of nonlinear RNNs allows for their efficient parallelization and scalability. Our work plays an important role from moving from heuristics to firm theoretical grounding for the importance of stability in the design of efficiently parallelizable nonlinear SSMs.\n\n\n## Clarifying Figure 1\n\nFigure 1 is our banner schematic, so we appreciate feedback on it. While this year’s NeurIPS does not allow us to upload images, we provide answers to the reviewer’s first question and have adjusted the schematic to be more clear (especially adding axis labels).\n\nFigure 1 is a schematic illustrating Theorem 2, which proves that the predictability $\\lambda$ of a dynamical system (illustrated by the top row, labelled “state space”) determines the conditioning $\\mu$ of its merit function landscape (illustrated by the bottom row, labelled “merit function”). The left column (blue font), shows that predictable/contracting systems result in well-conditioned loss landscapes. The right column (red font) shows that unpredictable systems result in ill-conditioned (flat) landscapes.\n\nThe system we used to generate this banner schematic is a vanilla RNN with a hidden state of dimension 1; we can adjust the scalar weight to result in either predictable or unpredictable dynamics. \n\nThe top row shows the trajectories over time of these RNNs. Thus, the x-axis is time, while the y-axis is the value in state space. The gray arrows show flow fields based on the particular value in state space, while the solid lines (blue for predictable, red for unpredictable) indicate the trajectories resulting from two different initial conditions. Trajectories from predictable systems converge over time, while trajectories from unpredictable systems diverge over time.\n\nThe bottom row (the contour plots) show the resulting merit function landscapes (see equation 2 for the definition of the merit function). While the merit function is a function of $\\mathbf{s} \\in \\mathbb{R}^T$, for visualization we pick two time steps for the axis labels. In this particular visualization, the x-axis is $s_{50}$ and the y-axis is $s_{51}$. The color of the contour plot is the value of the merit function for particular values of $\\mathbf{s}$, with darker color indicating a lower value of the merit function (pitch black correspond to the minimum $\\mathcal{L}(\\mathbf{s}^*) = 0$). The blue/red dotted lines indicate the optimization trajectory of DEER for these two RNNs. In the predictable setting, the merit function landscape is well-conditioned (bowl shaped), and DEER quickly moves to the minima (the blue star). In the unpredictable setting, the merit function landscape is ill-conditioned (has a long, flat, narrow valley), and DEER oscillates around the minima (the red star). \n\nWe have added axis labels and expanded the caption to improve the readability of this banner schematic."}]

Author rebuttal:

Meta-review: This paper presents a theoretical framework demonstrating that a system’s predictability governs the number of optimization steps required for evaluation. This relationship is particularly significant for assessing the efficiency of parallel computing. A key strength of the work lies in its central contribution of linking the Largest Lyapunov Exponent to PL conditioning in optimization problems. The formalization of predictable and unpredictable systems is also creative and provides a strong foundation for the authors’ main results. As most of the reviewers’ comments have been adequately addressed, this paper is suitable for publication.
Decision: Accept (poster)