A Unifying Framework for Parallelizing Sequential Models with Linear Dynamical Systems
Authors: Xavier Gonzalez
E. Kelly Buchanan
Hyun Dong Lee
Jerry Weihong Liu
Ke Alexander Wang
David M. Zoltowski
Leo Kozachkov
Christopher Re
Scott Linderman
Venue: Accepted by TMLR ()

Abstract: Harnessing parallelism in seemingly sequential models is a central challenge for modern machine learning. Several approaches have been proposed for evaluating sequential processes in parallel using iterative fixed-point methods, like Newton, Picard, and Jacobi iterations. In this work, we show that these methods can be understood within a common framework based on linear dynamical systems (LDSs), where different iteration schemes arise naturally as approximate linearizations of a nonlinear recursion. Moreover, we theoretically analyze the rates of convergence of these methods, and we verify the predictions of this theory with several case studies. This unifying framework highlights shared principles behind these techniques and clarifies when particular fixed-point methods are most likely to be effective. By bridging diverse algorithms through the language of LDSs, the framework provides a clearer theoretical foundation for parallelizing sequential models and points toward new opportunities for efficient and scalable computation.

Published in Transactions on Machine Learning Research (4/2026)
A Unifying Framework for Parallelizing Sequential Models with Linear Dynamical Systems
Xavier Gonzalez*
zavier 88stinoodd edu
Statistics Deportment and Wu. Tsai Neurosciences Institute
Stanford University
E. Kelly Buchanan*
kelly. wuchannnn@ttnnforrrrorooo edu
Statistics Deportment and Wu. Tsai Neurosciences Institute Slanford University
Hyun Dong Lee
hdlle @stanford.ccmm edu
Computer Science Deportment end Wu Tsoi Neurosciences Institute Stanford University
Jerry Weihong Liu
jul500stanford edu
Institute for Computational and Mathematical Engineering Stanford University
Ke Alexander Wang
alcwang@cs.stanforee edu
Computer Science Deportment
Stanford University
David M. Zoltowski
dzoltowskii@stannordd edu
Statistics Department and Wu. Tsai Neurosciences Institute
Stanford University
Leo Kozachkov
leokoz8@broune edu
School of Engineering and Carney Institute for Brain Science Brown University
Christopher Ré
chrismre whrismre@sstanforddc edu
Computer Science Department
Stanford University
Scott W. Linderman
scott. BinttermamnSsttnforrr edu
Statistics Department and Wu Tsai Neurosciences Institute
Stanford University
Reviewed on OpenReview: https: // openreview. net/ forum? id= fugGAAIGrr
Abstract
Harnessing parallelism in seemingly sequential models is a central challenge for modern machine learning. Several approaches have been proposed for evaluating sequential processes in parallel using itcrative fixed-point methods, like Newton, Picard, and Jacobi itcrations. In this work, show that these methods can be understood within common framework based linear dynamical systems (L.Sss), where different iteration schemes arise naturally as approximate linearizations of nonlinear recursion. Moreover, theoretically analyze the rates of convergence of these methods, and we verify the predictions of this theory with several case studies. This unifying framework highlights shared principles behind these techniques and clarifies when particular fixed-point methods are most likely to effective. By bridging diverse algorithms through the language LDSs, the framework provides clearer theoretical foundation for parallelizing sequential models and points toward new opportunities for efficient and scalable computation.
indicates lead authors. Author contribution statement in Appendix A.

---

Published in Transactions on Machine Learning Research (4/2026)
1 Introduction
Sequential processes are ubiquitous in machine learning models. Evaluating a recurrent neural net- work (Goodfellow et al., 2016), sampling a diffusion model (Sohl-Dickstein et al., 2015; Ho et al., 2020; Song et al., 2021b), generating from a deep state space model (Gu et al., 2022; Smith et al., 2023; Orvieto et al., 2023; Gu & Dao, 2024), and unrolling layers of a deep neural network (He al., 2016; Vaswani et al., 2017) all involve sequential computations. Naively, these sequential computations require time proportional to the sequence length the network depth and do not take full advantage of hardware accelerators like GPUs and TPUs (Lim et al., 2024). However, for one important class of computations-linee recursions or linear dynamical systems (LDSs) this bottleneck has been overcome using techniques like the parallel scan" (Blelloch, 1990; Fatahalian & Olukotun, 2024). Indeed, the parallelizability of linear recursions is key to scalably evaluating many deep state space models (Smith et al., 2023; Gu & Dao, 2024). A natural question whether other sequential processes in machine learning, such as those with nonlinear recursions, can be similarly accelerated.
On first inspection, the parallel scan algorithm does not seem to generalize to noulinear recursions. The parallel scan applics to linear recursions because the composition of two lincar functions remains lincar, whereas the composition of two nonlinear functions is generally more complicated. For example, the com- position of two quadratic functions is quartic. Nevertheless, recent works have proposed several techniques to parallelize nonlinear recursions, including Jacobi (Song et al., 2021a), Picard (Shih et al., 2023; Lu et al., 2025), Newton (Danieli et al., 2023; Lim et al., 2024) and quasi-Newton (Tang et al., 2024; Gonzalez et al., 2024) iterations. These techniques were originally applied to different machine learning problems, includ- ing parallclizing the training of nonlincar RNNs and sampling from diffusion modcls, and follow different notations and intuitions, which obscures their underlying similarities.
In this paper, we show that Jacobi, Picard, and Newton iterations all parallelize nonlinear recursions over the sequence length by iteratively linearizing the nonlinear recurrence and evaluating the resulting lincar dynamical system parallel. While the comnections between Picard and Newton iterations and their con- vergence rates are well-known in applied mathematics (Ortega & Rheinboldt, 1970), we make them explicit for nonlincar recursions and characterize when each method works and why.
Our contributions include the following:
A unifying framework casting all four mothods (Newton, quasi-Newton, Picard, and Jacobi) as iterative LDS evaluations (Section 2);
A theoretical analysis showing that the rate of convergence depends on the stability of the resulting LDS and the fidelity of its approximation to the true linearized dynamics (Section 3); and
Three case studies the group word problem, evaluating a nonlinear RNN, and sampling from a dis- crotizod Langevin cattusion-valdaattiinn that this analysis prodicts which method suits which problem (Section 4).
These contributions unify and clarify several recently proposed methods in the machine learning literature, and they highlight the rauirality of linear dynamical systems for parallelizing seemingly sequential processes.
2 Unifying fixedpoiitt iterations using linear dynamical systems
We first introduce notation for evaluating a generic sequence model. Let F, E R denote the state at time t, and let ft denote the corresponding transition function at that time point. Throughout this paper, we will use D to denote the dimension of the hidden state and T' to denote the sequence length.
Problem Statement (Sequential Evaluation): Evaluate the sequence X1.7 (x1, T2, ... TT) starting
from To via the recursion,
Ft+1 = Se+1 (xz).
(1)
I See Appendix B for introduction to the parallel scan.
2

---

Published in Transactions on Machine Learning Research (4/2026)
We omit input dependencies for simplicity, but note that an imput uy can be incorporated into the definition of the transition function by letting f++1(1)) == [( ue).
The recurrence described in eq. (1) cannot be evaluated in parallel in its original form because I++1 depends directly on 2t, creating a chain of dependencies. As a result, the computation of each state must wait for the previous state to be computed. This approach takes ((T) time to evaluate the sequence. Moreover, this inherently sequential approach prevents us from fully leveraging modern hardware accelerators, which can dramatically accclerate parallelizable computations. These nonlinear recursions are ubiquitous, appearing, for example, in the denoising pass of diffusion model, in the forward pass of a nonlinear RNN, or in the recurrence relations in implicit layers and deep equilibrium models.
Fixed-point methods ofler a promising alternative: rather than computing the sequence step by step, make (0) initial guess for the entire trajectory. We denote this initial guess by xi.T. We then iteratively refine this guess, operating over the entire sequence length in parallel, denoting the guess after i fixed-point iterations x In particular, the current guess iteration is further refined by applying fixed-point operator A: R + R.DD (which depends on the functions fe and the initial condition zo) according to
(++1 I: = A (
(2)
In order to be a fixed-point operator, A should have a fixed point xi.r- That is, should satisfy that xi.r = (xtrr) where xi.r is the unique root of the system of equations given by
T++1 fxt1 O, V ... TO
(3)
A stopping criterion is used to determine when the fixed point iterations in eq. (2) have converged up to some level of desired numerical accuracy.
Many fixed-point operators can be constructed to satisfy this constraint. However, in the context parallel evaluation of sequences, we can often be much more specific than a generic operator A. In fact, for the four fixed-point methods we consider in this work, the fixed-point operators A solve a linear varying system the sequence length, with the common form,
(f+1) T Fi+1 = + Aeti (I (6+11 (1((tt)), ((9)
where the transition matrix dott E RRD is determined by the dynamics functions fe+1 and the current guess for the state 2,(") Importantly, eq. (4) is lincar dynamical system for the yet unknown x xi- (i+1) of the currently known X1,7, as #++1 is linear function of x (5+1) plus bias 6,+1 = f(+1(1(") (1)) Ac+14() which only depends on the states from the previous fixed-point itcration.
(++1) in terms
The transition matrix Åtς1 can be thought of as an approximation to the Jacobian of the dynamics function 8/4 1//z// Different fixed-point methods simply use different approaches to linearizing the dynamics function, as shown in Table 1.
Parallel Scan
Init
A1
A₂
As
A
000o 7) dipph
Step
A
AzAr
As
A.A.
Importantly, because the recursion in eq. (4) is an LDS, it can be evaluated with a parallel scan (Blel- loch, 1990; Nguyen, 2007). The parallel sean core primitive that allows the unrolling of T steps from an LDS, 2441 14 + b in O"logT) time sors. The logarithmic time in the sequence length comes from a divide-and-comquue divide-and-comquee approach made possible from the fact that composition affine functions closed, i.c. another affine function. To show the closure of composition of affine func- tions more explicitly, note that if fe(x:) Aşa + by
A
AzAr
AsAzA
A,A
Endt
Zg e
A
3
A,A:A2A
34
machine with O(T') proces-
T processors
inactive
Figure 1: Parallel scan for LDS. We illustrate the parallel sean for sequence of length T = 4 for the simple LDS x = A
3

---

Published in Transactions on Machine Learning Research (4/2026)
Table 1: Summary of fixed-point iteration schemes as linear dynamical systems. We list the methods by the order of their approximation. While higher order methods may converge in fewer iterations, each iteration may be more costly. For example, the prefix sum and parallel scan have Olog T') depth, while single Jacobi iteration has constant depth. For all the methods, each iteration is an LDS, i.e. they can be written in the form of eq. (4) where A++1 is the transition matrix. These methods are guaranteed to converge in at most T iterations (Gonzalez et al., 2024, Proposition 1). "Order" mcans the highest number of derivatives laken: Newton and quasi-Newlon methods use first. derivatives, while Picard and Jacobi methods do use derivatives of ft.
Fixed-point method
Order
Transition matrix Å++1
Parallelization
0f+t (1) dz afft1 Ovs
Parallel Scan
Newton
first-order
(at
(dense matrix multiplication)
(x) (1))
Parallel Scan
Quasi-Newton
quasi first-order
diag
(elementwise vector multiplication)
Prefix Sum
Picard
zeroth-order
(vector addition)
Map
Jacobi
zeroth-order
(embarrassingly parallel)
Algorithm Fixed-point methods for evaluating sequences using LDSs and parallel scan
procedure PARALLELFIXEDDOME
To, initial guess (0)
x1.7> tolerance
for i i= 0, 1,.,,..7 T do
Å LNBARRZDDYAAMCSS(,, 1:T) (++1) 1:7
For all / in parallel
+ EVALUATELDS((o, AET f. 1
Using parallel scan
if COMPUTEERROR(so) x (i-1) ff) < then
break
return x xiT
(1+1)
and (x) A x + b141, then fe1(f((x)) A411A + (At42b, + bet1). Therefore, as shown in Figure 1, by combining neighboring affine functions together, we oblain transition maps from the initial condition zo to all other states in O(log T) time. We provide further discussion of the parallel scan algorithm for the interested reader in Appendix B.
Next, we discuss in more detail how the four prominent fixed-point methods discussed in Table 1 reduce to iterative application of L.Ss when used to solve a recursion. Fundamentally, all of these methods parallelize nonlinear recursions by iteratively linearizing and evaluating them, as we indicate Algorithm However, the different methods use different approximations of the Jacobian 0fe-s/aw, of the dynamics function for their transition matrices A(+1-
2.1 Newton iterations
Bellen & Zennaro (1989), Gander & Vandewalle (2007), Danieli et al. (2023), and Lim et al. (2024) demon- strated that when the fixed-point operator A is constructed as an appropriately designed LDS, built from a linearization of original noulinear recursion f, then each application of equivalent to an iteration of New- ton's root-finding method on the system of equations given in eq. (3). Specifically, each Newton fixed-point (f+1) Ax (x)) is defined by the lincar recursion.
iteration, X1.T
(++1) (++1)= f1+1(()))
0h+t dz,
(o( (z (i+1) -1((0) =
(5)
4

---

Published in Transactions on Machine Learning Research (4/2026)
which we recognize the first-order Taylor expansion of the nonlinear recursion. Because this fixed-point iteration uses first derivative, we refer to it as a first-order fixed-point iteration. We note that for Newton iterations, the transition matrix is exactly the Jacobian of the dynamies function, Á4t 0f++ Ox. )
Crucially, because eq. (5) is a linear dynamical system, it can be evaluated in parallel over the sequence length using associative scan (Biclloch, 1990). However, each iteration requires O(TD²) memory to store the T Jacobian matrices, and OTT work to compute the matrix-matrix multiplications in the parallel scan. This expense is prohibitive for large state size or sequence length, which motivates the quasi-Newton iterations we discuss next.
2.2 Quasi-Newtor iterations
The compute and memory costs of full Newton iterations have motivated wide literature on quasi-Newton methods (Nocedal & Wright, 2006). A particularly simple way to turn the Newton iteration in eq. (5) into a quasi-Newton iteration that uses a parallel associative scan was proposed by Gonzalez et al. (2024): just use the diagonal of the Jacobian of the dynamics function. Specifically, each of these quasi-Newtor fixed-poine iterations AqN is given by
(1+1) x "+4 = f1+1(x + diag
05++1 Oz.
(x
(6+1) xx =(?),
(6)
which we recognize as an LDS with transition matrix b.t1 = diag Ofet1 Ox
(I) (x(²))
With this transition matrix, the parallel scan requires only O(TD) memory and O(TD) work, and is therefore more computationally efficient than the full Newton iteration eq. (5). However, quasi-Newton methods may take more fixed-point iterations to converge. As we will show in Sections 3 and 4, whether Newton and quasi-Newton methods are faster for evaluating given sequence depends accuracy of the approximate Jacobian, the stability of the linearized system, the choice of hardware, and the scale of the problem. In many cases, the large memory consumption of full Newton iterations renders it infeasible for large-scale problems.
many factors including the
Another consideration is the manner in which A++1 computed. A primary goal of quasi-Newton methods (Broyden, 1970; Dennis & Schnabel, 1996) is avoid computation of 0fes1/aa, Simply computing these Jacobians with autodifferentiation and then taking their diagonal still requires D) function calls, which can be costly. However, in sequence modeling with recurrent neural networks (RNNs), the diagonal elements of the Jacobian can often be written in closed form, allowing for its evaluation in a single function call. This approach is taken in Gonzalez et al. (2024) and Danieli et al. (2026). Another approach proposed in Zoltowski et al. (2025) is to stochastically estimate the diagonal using the Hutchinson estimator. This stochastic approach also requires only one function call but generally applicable even when closed forms of Art1 are not available. We employ this efficient stochastic estimator in all of our experiments. Finally, third approach comes from the multisecant or Broyden family of methods (Fang & Saad, 2009; Walker & Ni, 2011), including Anderson accelcration (Anderson, 1965; Bai et al., 2019; Tang et al., 2024), which build up estimates of derivative information over the fixed-point iterations.
In general, any structured approximation of the Jacobian matrix that remains closed under matrix multipli- cation2 could be used to form a quasi-Newton fixed-point iteration amenable to a parallel scan. For example, Zoltowski et al. (2025) introduces a parallel quasi-Newlon method where each block of the matrix is diago- nal. The broader development of quasi-Newton methods that fit into the unifying LDS framework discussed in this paper an important direction for future work, which we elaborate on in Section 6. However, for simplicity, henceforth in this paper when we refer to "quasi-Newton iterations," " we restrict ourselves to the simple diagonal approximation shown in eq. (6).
2Closed in the sense that the product, Ar+1 At, has the same structure as Ar+1 and Ai, with diagonal matrices. See Ap- pendix B for more detail.
5

---

Published in Transactions on Machine Learning Research (4/2026)
2.3 Picard iterations
A seemingly different approach to using fixed-point iterations is Picard ions, which were used by Shih et al. (2023) to parallelize sampling in diffusion models. Picard iterations are often used in the context of evaluating ODEs, where
& = g(, ()-
After Euler discretization with step size A, we obtain the discrete-time recursion,
Le+1 It g( () A.
(7)
The Picard fixed-point iteration, x x1.7"
Ap(x( ; then given by,
(++1) 1+1
rD + Σ
g(xxx))) A.
(8)
==0
Because Picard iterations do not use any derivatives of the discrete-time recursion, we call them zeroth-order fixed-point iterations.
Shih et al. (2023) proves by induction that for any dynamical system given by eq. (7), the fixed-point iterations given by eq. (8) will converge to the true trajectory in at most T iterations. Similarly, Gonzalez et al. (2024) proved that for any dynamical system given by eq. (1), both the Newton fixed-point iterations given by eq- (5) and the quasi-Newton fixed-point iterations given by eq. (6) will converge to the true trajectory in at most T iterations. The similarity of these results and techniques begs the question as to how Picard and Newton iterations relate to each other. Our first result shows that Picard iterations are in fact a type of quasi-Newton iteration, where we approximate the Jacobian of the dynamics function by the identity matrix.
The Picard iteration operator Ap given by eq. (8) is a special case of an LDS, eq. (4), where the transition matriz the identity,
Proposition 1.
A41 ID.
Proof. Define fi+1(Z)) = It + g(x,, t) A. Then, from eq. (8) it follows that
(i+1)
zt (i+1)
.t)
(6)
#i+1
((+1) (1) +x( (f)
+ g(x(() ,t) A
xt Itt
= S (1) + (ft (n((tt)
(i+1) - (
This is exactly of the form of the generic linear recursion shown in eq. (4), with A4t1 = Lp-
An important consequence of Proposition 1 is that like Newton iterations and quasi-Newton iterations, Picard iterations can also be cast as an LDS. In Newton iterations, the full Jacobian 0f-n/aw is used in LDS; in quasi- Newton iterations, the diagonal approximation diag("f--1aax.] used; and in Picard iterations, the identity In is used. The Picard iteration is more compute and memory efficient. than even quasi-Newton, requiring only vector additions via the prefix sum algorithm. However, it is generally a less faithful approximation and takes more iterations to converge, unless the Jacobian is well-approximated by the identity.
2.4 Jacobi iterations
Yel another seemingly different class of fixed-point methods are Jacobi iterations (Orlega & Rheinboldt, 1970), which were used by Song et al. (2021a) to accelerate computation variety of settings in machine learning, such feedforward networks with skip connections. Jacobi iterations are also zeroth-order fixed- point method, and are commonly used to solve systems of multivariate nonlinear equations of the form,
hx(x1:)) 00 Vt (1,..., T).

---

Published in Transactions on Machine Learning Research (4/2026)
(i+1) (c) Instead, the Jacobi fixed-point operator, X As(x);r), equations in parallel to obtain XIT
solves the following system of T univariate
(6+1)
( (6) by (x1))
(3) Lt, 4t1 1.1.... o WE (I....,T))
(9)
Song et al. (2021a) considers in particular the problem of solving recurrence relations of the form #++1 (x,e)) and proves that, for such system, Jacobí iterations converge in at most. 7 iterations. This result is directly analogous to similar such results and proofs in Gonzalez al. (2024) and Tang et al. (2024) for Newton and quasi-Newton iterations and Shih et al. (2023) for Picard iterations. In fact, in the context of iteratively applying LDSs to parallelize Markovian state space models, we prove that Jacobi iterations are a type of degencrate quasi-Newton iterations, where we "approximate" the Jacobian of the dynamics function by zero.
Proposition 2. When applied to a Markovian state space model as in eq. (1), the Jacobi iteration operator As specified by eq. (9) special case of the common form, eq. (4), where,
A+t1= =
Proof. a Markovian state space model, the recurrence relation always takes the form specified eq. (1), i.c. Z441 fit1 (2+). Thus, Jacobi iterations take the simple form
(i-1) Fit1 = S++1
(++1) (f+1) does not depend on xt= we see that the transition matrix is zero. Because 2(41
We have shown that evaluating nonlinear recursions via Newton, quasi-Newton, Picard, and Jacobi fixed- point iterations reduces to itcratively solving lincar dynamical systems. An important corollary is that these methods are guaranteed converge all problem seltings in most. iterations (Gonzalez et al., 2024, Prop. 1). However, the precise convergence rates of the different fixed- fixed-point methods will depend on the accuracy of the Jacobian approximation,
show next.
3 Analysis of convergence rates
In this section, analyze the convergence properties of the fixed-point methods above. We show that the convergence rate of these fixed-point methods depends on two factors: (i) how well the transition matrix 44t1 approximates the true dynamics Jacobian A == 0f41/002; and (ii) the stability of the LDS with transition matrices A
3.1 Derivation of fixed- trentpoint iterations from high-dimensions root-finding.
We first review the derivation of the Newton fixed-point iterations in eq. (5) from the system of equations given by eq. (3), highlighting how the other fixed-point iterations naturally arise from different approxima- tions of 8f41//o.. This background will be useful for proving convergence rates in Section 3.2.
Returning to the system of nonlinear equations in eq. (3), define the residual function r: PDRTD RTD
r(xiT)= = [x] fzo), x2 - fe(x1), F3 - fs(a2)....,77 - fr
(10)
We define J(x1r) E TDXTD to be the Jacobian of the residual r(x1:r).
Id
0
0
8f2
(x1)
ID
0
De
Òr (xx.T) = 0xiT
8fa
(r2)
Ip
0
0
(xI:T)
(11)
In
0
0
(TT-1)
Ip/
7

---

Published in Transactions on Machine Learning Research (4/2026)
Therefore, in this notation, every Newton's method update takes the form
(++1) AN( = x ==
Using matrix manipulations, it can be shown (Danieli et al., 2023; Lim et al., 2024; Gonzalez et al., 2024) that each Newlon fixed-point iteration is given by the linear recursion in eq. (5).
Analogously, for the other fixed point iterations considered in this paper (ef. Table 1), which use transition matrices A to approximate 0ft-//m1, we can define a block-bidiagonal E RDDxTD simply by substituting (x4) for 0fe1s/Oz, in eq. (11). The corresponding fixed point iteration A takes the form
A(x = x -r((11)
(12)
Different fixed-point methods A give risc to different matrices J, which impacts their convergence rates.
3.2 Convergence rates of fixed-poing poil-poont iterations
We derive convergence rates for all fixed-point operators discussod in this paper. The analysis follows standard techniques in the literature (cf. Kelley (1995, Thm 5.4.1)), which Lu et al. (2025, Prop 4.) applied to parallelizing sequential models with Picard iterations. Our unifying framework allows us to that their analysis generalizes straightforwardly.
(6) x xi.r RTD denote the error of trajectory Xi7- Note that the error depends the current trajectory and the dynamics functions, but suppress these dependencies to keep notation lightweight. For any of the fixed-point methods discussed in this paper, can bound the convergence rate of this error.
Let ei.
Proposition 3 (c.f. Proposition 4 of Lu et al. (2025)). Consider fired-point solver with updates given by eg. (12) for some matrie [(x with form specified by eq. (19). Let L, be the marimumm of the Lipschilz constants of aff/am.... Then |/e /ef(+) 112 satisfies
|
(ftt S
(1
L l + 2 (1lef2
(13)
where || |l2 denoles the spectral norm of a matria and the &a norm of a vector.
Proof. Starting from eq. (12), we subtract xix from both sides to obtain
e e
Next, we Taylor expand r(-) around ((0 to obtain
r(xi.r) -[x R(e
where ((( (1) is the sccond-order remainder function and has norm bounded by 1/// times the Lipschitz constant of [x which Theorem 3 of Gonzalez et al. (2025) shows is bounded by L. Since r(xi.)) = 0, it follows that
(++1)
- -xxx ( (x
e R(e(.))
e
(14)
The result follows by taking norms on both sides and using the triangle inequality.
3.3 Intuitions about rates of convergence
Proposition 3 shows that, to first order, the error is determined by the discrepancy between the chosen lincar operator J and the true Jacobian J of the residual. Moreover, that as [+) approaches zero, the
8

---

Published in Transactions on Machine Learning Research (4/2026)
first-order lerm, which is linear in [| must eventually (under strong enough continuity assumptions) dominate the sccond-order term, which is quadratic in 1er.rll.. Typically, we would say the rate of decrease in | approaches an asymptotic linear rate y given by
7 = [ 1 -(xt:r)- J(xi:r)
(15)
Discussions of asymptotic linear rate subtle in our setting, where all fixed-point methods guaranteed to converge in T iterations (see our discussion in Appendix F.1). Nonetheless, the functional form of y provides useful intuition about the convergence rates of different fixed-point methods.
Let's study the two factors the asymptotic linear rate y, which must be less than one to be meaningful. Lu et al. (2025) left quantifying these two factors for future work, which we now bring to fruition. The first factor is the spectral norm of the inverse of the approxinate Jacobian. Because (x::r) is block bidiagonal matrix (see eq. (11) and the definition that follows), its inverse has a block lower triangular structure of the form
ID
0
0
0
A2
Ip
0
0
-(x1=a)===
(16)
AzA2
Aa
In
A4
A4A3
A4
1p/
shown above for T = 4.
From eq- (16), we see that the blocks of J(x1r))-1 are products of the transition matrices A from the chosen fixed-point method. In particular, the chosen fixed point method results in an unstable LDS with ||A++1ll2 > for many time points, the spectral norm of ((x:r)) can be much larger than one. Intuitively, fixed-point methods resulting in unstable LDSs will have slower rates of convergence since instability leads to mumerical divergence. We elaborate this point in Appendix C.3.2.
The second factor is the spectral norm of the difference between the approximate and true residual Jacobian. The following lemma shows that we can control this factor in terms of the spectral norms of the differences between the approximate and true dynamics Jacobians. Intuitively, the rate of convergence will be faster if the transition matrices Å4t1 are closer to the true dynamics Jacobian A == 0f11//zz, in spectral norm.
Lemma 1. If ((x:T) is given by eq. (19) and J(x1:) is given by eg. (11), then
diff(AA) = [((1:)) J(x
[(44) A441(2)))22)
max 151ST-1
where A++1
0f (Iz).
Owr
Proof. Appendix C.3.1.
Proposition 3 is an upper bound on the norm of the error of each fixed-point iterate. As an upper bound, this result, can not always fully predict the precise trajectory of the norm of the error (see Appendix C.2). Nevertheless, this proposition provides helpful intuitions that are borne out in the case studies below.
4 Case studies: performance of the different fixed-point methods
In this section, we consider three empirical case studies that illustrate how the unifying framework and convergence analysis presented in this paper provides guidance about which fixed-point schemes will excel in which settings. This concordance based on the structure of the Jacobian of f++1 and the relative computational cost of different fixed-point methods. In nutshell, we advise:
Use the simplest approrimate Jacobian as possible, but no simpler.
To elaborate: simpler approximate Jacobians are less computationally expensive, meaning that each fixed- point iteration is more efficient. So, if lower-order fixed-point. method (e.g., Picard Jacobi) converges in

---

Published in Transactions on Machine Learning Research (4/2026)
small number of iterations, it may reach x" in less time than higher-order methods. However, if the higher- order fixed-point method (e.g. Newton or quasi-Newton) converges in far fewer iterations than lower-order methods, then the increased computation of the higher-order fixed-point method may be worthwhile. As supported by the thcoretical analysis in Section 3, the number of iterations needed for a fixed-point method to converge is related to the difference in spectral norm between A, and A, == 0ft+s/am. We support this intuition with the following case studies, with further experimental details provided in Appendix E. Our code is available at https:// httpa:///gthhb....//e http://gttuub.cooo/l.a/ httpc////.ddeeaaaaaaaaaaaallaaaa
Baapiing_wtthh_heeeeeeeee _ds
4.1 Case study #1: Solving the group word problem with Newton iterations
Newton iterations should outperform quasi-Newton and Picard iterations in settings where the Jacobian of the recursion, ft+1, is not well approximated by its diagonal, the identity matrix, or the zero matrix. One example of such a recursion is the group word problem, which has been used to theoretically and empirically assess the limits of sequential modeling architectures for state-tracking tasks (Kim & Schuster, 2023; Liu et al., 2023; Merrill et al., 2024; Grazzi al., 2025; Schöne et al., 2025). In the sequence-modeling community, the term "group word problem" is defined follows.
Definition 1 (Group Word Problom). Let G be finite group and let 91, 12.99..... 9r be sequence of group elements. The group word problem is to evaluate the producl 9ı 92 ..2-- Since each gr E G, the producl of these group elements belongs to G as well.
Merrill al. (2024) emphasizes that nonlinear RNNs in both theory and practice are able to learn the group word problem in arbitrary groups to high accuracy with only single layer, whereas compositions of popular linear RNNs linked by nonlinearities, such as $4 (Gu et al., 2022) and Mamba3 (Gu & Dao, 2024), require number of layers that grows with T. Merrill et al. (2024) emphasizes that recurrent architectures with nonlinear transitions are well-suited for solving the group word problem, because in theory and practice, such architectures can learn the group word problem to high accuracy with a single layer. Other literature has explored the value of matrix-valued states (Beck et al., 2024; Grazzi et al., 2025). However, in Proposition 4 below, we show that neither nonlinearity nor matrix-valued states are needed to understand or solve the group word problem. Instead, the problem can be formulated as LDS with vector-valued states and input-dependent transition matrices.
Proposition 4. Let G be a finite group. Then there erists some D S GG] for which we can represent the group word problem a lime-varying LDS, (e+1(x)) A4t13T, with slates in RDD denoling the running product of group elements and transition matrices Act1 E DxD that depend on the input 9t+1-
Proof. By Cayley's theorem, any finite group G can be embedded in a symmetric group Sp, for some D < [G]. Therefore, by choosing the initial state xp E PR to have D distinct entries (a "vocabulary" of size D), we can use the tabular representation of permutations (Artin, 2011, eq. 1.3.2) to represent an element of SD as xe (by a permutation of the elements of zo). We can also choose A411 E TR. DxD to be the permutation matrix corresponding to the embedding of 9t+1 in Sp. since any element of Sp can be represented as a D D permulation matrix (e.g., see Figure 2B). Consequently It. AAA A2 is an embedding of an element of G in Sp in the tabular representation. In fact, I; DR represents the running product 9192. ...2...99 91-19 which is preciscly the goal of the group word problem.
Though we have cast the group word problem as a time-varying LDS with (xt) = A we can still evaluate this recursion with any of the fixed-point methods described above. Since the dynamics are linear, the Newton iteration corresponds to evaluating the LDS with a parallel scan, and it converges in one iteration. While other methods would require more iterations to converge, they could still be more efficient in wall-clock time, since they use less memory and compute per iteration. However, we can use the Jacobian approximation error of the different fixed-point methods (Lemma 1) to get sense if the other fixed-poit methods are likely to excel in this setting.
Mamba allows input-dependent dynamics matrices but they must be diagonal, which prevents single Mamba layer from implementing the particular LDS Proposition 4, which uses permutation matrices. Merrill al. (2024) also demonstrate that linear varying system with dense transition matrix can learn the group word problem.
10

---

Published in Transactions on Machine Learning Research (4/2026)
A) Group Word Prohlem
C)) Convergence Rates
Jacobi
Quasi-Newton
Sequential
Picard
Newton
Number of Fixed-Point Iters
21+1 A+t1 (in)( A++1 E RDD
10a
10a
10'
10a
10a
104
B) 1
Permutation 51 2 3)
Wallelock Time (8)
100
10
10-2
10-3
0ft:1/0r=A111
10a
10a
Sequence Length (T)
Figure 2: A single Newton iteration solves the Ss group word problem, whereas the number of iterations required for the other methods increases with sequence length. We consider the task of evaluating the product S5 group elements. A: The group word problem can be expressed as an LDS with input-dependent state- taae-transition matrices. B: An example input-dependent transition matrix A for permutation 5 4 3), in cycle notation. C: For each fixed-point method and range of sequence lengths, T, we compute the median (over ten random seeds) number of fixed-point iterations to converge (top) and the median wall-clock time (bottom). While a single Newton iteration is sufficient to solve the S problem, the number of iterations required for the other methods increases with the sequence length.
The state transition matrices of the group word problem are permutation matrices with spectral norm one, and (A,) 1. Furthermore, since with high probability there will be state transition matrix with no non-zero entries on the diagonal, it follows that (AAg)) = 1 while diff(Ar) = 2. Since we would expect to need (AA < 1 for a fixed-point method A to be effective, our theoretical analysis in Section 3 suggests that none of the fixed-point methods other than Newton will be effective on the group word problem.
We test this hypothesis with simple experiment simulating the Ss word problem, a standard problem the sequence modeling literature (Merrill al., 2024; Grazzi et al., 2025). In this setting, Figure 2 shows that quasi-Newton, Picard, and Jacobi iterations require nearly T iterations converge. On the other hand, see that Newton's method solves the S₂ word problem with just one fixed-point iteration, as expected since the true dynamics are linear. The speed-up is also apparent in the wall-clock time comparison, where we see that Newton is faster than other methods, regardless of T.
4.2 Case Study #2: Picard iterations struggle to parallelize RNNs
Next, we consider a task where Picard iterations struggle and other fixed-point methods excel, namely, parallelizing recurrent neural networks (RNNs) like the Gated Recurrent Unit (GRU; Cho et al., 2014). We evaluate GRUs with random paramcter initialization for different hidden dimension sizes D and sequence lengths T using sequential evaluation as well as fixed-point iterations. Figure 3B shows that at initialization the Jacobian has small entries (on the order of 0.1). In this regime, diff(A) and diff( are less than one, while diff (AAp) is greater than one. (Sce also, Figure Appendix E.3.) Thus, as expected, Newton, quasi- Newton, and Jacobi iterations excel in this setting, while Picard iterations converge prohibitively slowly, shown in Figure 3C.
11

---

Published in Transactions on Machine Learning Research (4/2026)
A) GRU dynamics
C) Convergence Rates
Jacobi Picard
Quasi-Newton
Sequential
Newton
2++1 (L 24) C +
D
[(()))
Number of Fixod- Point Iters
10
2, = (Linnear[[u; ([4 æ:l))
10a
10a
10a
F = or (LLinear[[,, mx]))
10a
10a
10a
Z Lanb (aube(Lnuaar((u (u,, rx. xz
101
B) Jacobian of GRU dynamics
10a
10a
10a
(D 8)
101
10.
Wallelock Time (8)
100
0.1
10
10"1
10
0.0
10-2
10
-0.1
10-3
10
10a
Sequence Length (T)
Figure 3: Picard itcrations struggle to parallelize RNNs. We evaluate GRUs with random parameter initialization for different sequence lengths T and hidden state sizes D. A: The nonlinear dynamics of GRU, following Feng et al. (2024), where x, is the hidden state, u, is the input, and the notation Linear[:-]] indicates a linear readout from the concatenation of two vectors. B: A representative Jacobian matrix 0f./oz from GRU trajectory, which is not well approximated by the identity matrix. C: For each fixed-point method and a range of sequence lengths, T, and state sizes, D, we compute the median (over ten random seeds) number of fixed-point iterations to converge (top row) and the median wall-clock time (bottom row). Picard iterations take nearly T iterations to converge, while the other fixed point methods yield order-of- magnitude speed-ups over sequential evaluation
4.3 Case Study #3: Jacobi iterations struggle to parallelize discretized Langevin diffusion
Based on the theoretical analysis presented in Proposition 3, we expect that if the Jacobian of the dynamics function is well-approximated by the identity matrix, then Picard should converge relatively quickly and at considerably lower cost, especially when compared to the other zeroth-order method of Jacobi flerations A canonical example such a system where the dynamics are close to identity comes from a discretization of Langevin dynamics (Langevin, 1908; Friedman, 2022). Langevin dynamics are a workhorse for MCMC (Besag, 1994) and motivated the development of score-matching methods (Song & Ermon, 2019), which are closely related to diffusion models DSihl-Dickstein et al., 2015; Ho et al., 2020; Song et al., 2021b).
In Langevin dynamics for potential é, the state evolves according to a nonlincar recursion with dynamics,
fe-1 = = eVo(z)) + V2ew,,
where e is the step size and w, sid In) Since the Jacobian is, (x)) In -V2 it follows that the identity approximation should work well with small step-sizes More generally, the identity approximation tends to be well-suited to problems where differential equation is discretized with small step sizes, such as when sampling from diffusion models (Holderrieth & Erives, 2025; Lai et al., 2025). In fact, simply by observing the structure of the Jacobian in Panel B of Figure 4, we observe that diff(.) operator for Newton, quasi-Newton, and Picard iterations in this setting will be close to zero, while diff (AA) will be close to one. Therefore, based on our analysis in Proposition 3, we hypothesize that the other fixed-point methods should dramatically outperform Jacobi iterations in this setting.
12

---

Published in Transactions on Machine Learning Research (4/2026)
Langevin Dynamics
C? Convergence Rates
Jacobi Picard
Quasi-Newton
Sequential
Newton
250
2141 = d4 (Vol))) + 2 un
Poind Iters
10
10
10a
Number
10a
10
0f
In (())
10a
107
102
102
Oz
10a
101
10a
10
Jacobian of Langevin dynsmiis
D 32)
102
10a
107
10
10'
10
10a
101
1.0
109
10'
100
0.5
Time
10-1
109
0.0
10-11
-0.5
10-2
10"
10
10
10-3
81.41 1114//2-Ac41
Sequenace Length (T)
Figure 4: Jacobi iterations struggle when the dynamics Jacobian is close to the identity. We evaluate Langevin dynamics for a potential 6. A: The nonlinear dynamics of Langevin dynamics for a potential & and step size c, where Ia is the state and wy is Gaussian noise. B: The Jacobian for Langevin dynamics well-approximated by the identity matrix, especially for small step size e = 1 x 10-5, C: We random seeds. Jacobi iteration
evaluate Langevin dynamics for larger dimensions, plotting the median consistently take T steps and are always slower than sequential, while the other fixed-point methods converge in fewer T steps and can be faster than soquential. The missing Newton iteration points indicate the GPU ran out of memory.
We test this hypothesis with a simple experiment shown in Figure We simulate Langevin dynamics on a potential ø given by the negative log probability of the mixture of two anisotropic Gaussians. In this sctting, Picard iterations take far fewer than T iterations to converge and can be faster than sequential evaluation. We note that quasi-Newton iterations, which include information only about the diagonal of the Jacobian of the dynamics, appear to have comparable wall-clock time, by virtue of taking fewer iterations to converge (though each fixed-point iteration involves more work).
Whether fixed-point. iterations are faster than sequential evaluation also depends on memory utilization. For example, Shih et al. (2023) and Lu et al. (2025) demonstrated wall-clock speed-ups when using Picard iterations for sampling from a diffusion model using a "sliding window" to only evaluate chunks of the sequence length where the parallel scan algorithm can fit. in memory. We discuss these considerations further in Appendix E.5.
5 Related Work
In this paper we unify prominent fixed-point methods for the parallel evaluation of sequences in the language of linear dynamical systems. While many papers have employed different fixed-point iterations for different problems in machine learning Lim et al. (2024), Danieli al. (2023), and Danieli et al. (2026) using Newton iterations, Tang et al. (2024) and Gonzalez et al. (2024) using quasi-Newton iterations, Shih et al. (2023) using Picard iterations, and Song et al. (2021a) using Jacobi iterations, among other works- to the best of our knowledge no one in the machine learning community has explicitly unified these different methods in the language of linear dynamical systems. It is important to note, however, that Gander & Vandewalle (2007) propose an expansive notion of parareal iterations which include these LDSs as a special case. However, Gander & Vandewalle (2007) focus on solving continuous time differential equations and do not discuss the parallel scan that has become workhorse for the application of parallel-in -nn----i---ie time methods on GPUs in machine learning.
13

---

Published in Transactions on Machine Learning Research (4/2026)
General unification of fixed-point methods: parallel-chord methods While connections between Newton's method and Picard iterations have been made before outside of the machine learning literature, our contribution is the tight coupling of these methods to LDSs in the context of parallel evaluation of nonlinear sequences. Ortega & Rheinboldt (1970, Ch. 7) considered the problem of solving a nonlinear equation F(x) = O. They showed that. Newton and Picard iterations are special cases of general iterative methods where each iterate given by
(1+1) =x(0) -x())_((x -[(x(0)-1F(xx)
(17)
for some matrix ))). We discuss the relationship between the unifying frameworks put forward in Ortega & Rheinboldt (1970) and in our paper at greater length in Appendix F.1. The primary difference is that by locusing on the setting of nonlinear sequence evaluation, we bring into greater locus the role of the Jacobian of the dynamics function. Moreover, by unifying fixed-point iterations in the language of LDSs, we emphasize their parallelizability over the sequence length using the parallel scan (Blelloch, 1990).
Convergence rates of fixed-point methods In the context of analysis of fixed-point methods in general, there is broad literature (Ortega & Rheinboldt, 1970; Young, 2014) on the convergence rates of different fixed-point methods. For example, Ortega & Rheinboldt. (1970) also proved convergence rates for iterative methods of the form in eq. (17). Though their methods have much in common with the proof techniques used to prove Proposition 3 of this paper, their provided results are actually trivial in the setting considered in this paper. Part of the reason for the inapplicability of the convergence results from Ortega & Rheinboldt (1970) to our paper is that Ortega & Rheinboldt (1970) consider the asymptotic setting, while has been firmly established that in the particular setting considered in this paper, Jacobi, Picard, quasi-Newton, and Newton iterations all globally converge in at most T iterations (Bellen & Zennaro, 1989; Shih et al., 2023; Tang et al., 2024; Gonzalez et al., 2024). For moving beyond this worst-case analysis, Gonzalez et al. (2025) uses an optimization perspective to show that the difficulty of parallelizing a dynamical system is directly related to the stability of the system, which can be thought of as the "average" spectral norm of 0f++//0m. Proposition 4 of Lu et al. (2025) nicely presents the foundations of the convergence analysis we present in Proposition 3 of this paper. However, we extend their work by applying it a wider varicty of fixod-point methods, explicitly bounding many quantities of interest, and demonstrating its relevance in simulation.
Other fixed-point methods: mixing sequential with parallel In this note, we focus on Jacobi, Picard, and Newton iterations because of their prominence (Song et al., 2019; 2021a; Shih et al., 2023; Danicli et al., 2023; Lim et al., 2024; Gonzalez et al., 2024; Grazzi & Zanella, 2025; Zoltowski et al., 2025; Farsang & Grosu, 2025; Danieli et al., 2026; Iacob et al., 2025) and their relationship to LDSs, as listed in Table 1. However, there is a wide literature on iterative solvers (Ortega & Rheinboldt, 1970; Young, 2014). Many of these other methods can also be parallelized over the sequence length, or provide a mixture of parallel and sequential computation. For example, Horton et al. (1995) apply the parallel associative scan (which they call parallel cyclic reduction) for parallel-in-time differential equation solvers, though they focus Jacobi and Gauss- Seidel iterations. In discrete time, Naumov (2017) shows how evaluating Markov chains can be cast as a system of nonlincar cquations and discusscs many techniques from numerical analysis for solving them, again focusing on Jacobi and Gauss--eidel,, Song et al. (2021a) extend this program with many deep learning experiments. Although Gauss-Seidel iterations reduce to sequential evaluation when applied to Markovian processes, Song et al. (2021a) also emphasize how the structure of the problem and hardware considerations dictate the optimal mixture of parallel and sequential computation. Parareal iterations (Lions et al., 2001) mix parallel and sequential computation by applying parallelization at multiple length scales, and have also been used to parallclize diffusion models (Selvam et al., 2024). Tang et al. (2024) also parallelize diffusion models using both generalization of Jacobi iterations, well as Anderson acceleration (Anderson, 1965; Walker & Ni, 2011), which they modify to be form of quasi-Newton. We discuss other fixed-point methods based on trust-region techniques in Appendix F.2.
We elaborate Appendix F.I.
14

---

Published in Transactions on Machine Learning Research (4/2026)
6
Discussion
This work unified a variety of approaches for parallelizing recursions via fixed-point iterations-inchuuee zeroth-order methods like Jacobi and Picard iterations as well as first-order methods like Newton and quasi- Newton crrationssunded a common framework. In each case, the iterates reduce to evaluating an appro- priately constructed linear dynamical system, which approximates the nonlinear recursion interest. This unifying framework provides insight into which different problems in machine learning are likely to benefit from which types of fixed-point iterations. In particular, we demonstrate that the structure of the Jacobian matrix of the dynamics function plays key role in determining which fixed-point method to use.
For this reason, understanding the structure of the Jacobian of the dynamics function is important for using our framework. Fortunately, there are many problems where the structure of the Jacobian matrix is known in advancc. As we showed in Section 4.1, the group word problem can always be simulated with permutation matrices for its dynamics. As we showed used in sampling from diffusion models and rolling out neural ODEs, have 0f/0= equal to the identity matrix plus correction term proportional to the discretization step-size. Morcover, shown in Zoltowski et al. (2025), the dynamics of position and momenta variables in Hamiltonian Monte Carlo (HMC) results in banded matrices. Furthermore, in sequence modeling, one can design recurrent neural network to have Jacobians with desired structure (Farsang & Grosu, 2025; Zattra et al., 2025; Danieli et al., 2026). Finally, if thore is truly no analytic information about the Jacobian in advance, its structure can be probed with fiite-differrnce methods.
Section 4.3, discretized roll-outs from differential equations,
Future directions Clarifying the relationships and properties of these approaches through the lens of linear dynamical systems also suggests promising areas for future study. One clear direction of future work is to explore additional approaches for exploiting problem-specific structure, using our unifying framework to develop new fixed-point iterations. For example, intermediate between Picard and quasi-Newton methods is a scaled identity approximation, A++1 a++1/p. If we had prior knowledge on the appropriate scaling factors, 0++1 E TR, we could avoid computing any Jacobian-vector product evaluations. More generally, there exist other groups of structured matrices with compact representations that are closed under composition such that parallel evaluation of the LDS would be computationall) efficient. Examples include permutation matrices, block block-diagonal matrices, and block matrices where each sub-block diagonal, among others. Future work should enumerate these cases and investigate problem-speeific applications where they are appropriate. One example application is for more efficient parallelization of the group word problem using compact representation of permutation matrices, as was done in concurrent work by Terzić et al. (2025).
Understanding the shared backbonc of these fixed-point methods can give practitioners guidance about which methods to use for which problems. As parallel evaluation of seemingly sequential processes becomes increasingly important in machine learning, these insights may provide valuable guidance to the field.
15

---

Published in Transactions on Machine Learning Research (4/2026)
References
Donald G Anderson. Iterative procedures for nonlinear integral equations. Journal of the ACM (JACM), 12 (4):547 560, 1965.
Michael Artin. Abstract Algebra. Pearson, 2nd edition, 2011. ISBN 9780132413770.
Shaojie Bai, J Zico Kolter, and Vladlen Koltun. Deep equilibrium models. In Advances in Neural Information Processing Systems (NeurIPS), 2019.
Maximilian Beck, Korbinian Pöppel, Markus Spanring, Andreas Auer, Oleksandra Prudnikova, Michael Kopp, Günter Klambauer, Johannes Brandstetter, and Sepp Hochreiter. xLSTM: Extended Long Short- Term Memory, In Advances in Neural Information Processing Systems (NeurIPS), 2024.
Alfredo Bellen and Marino Zemnaro. Parallel algorithms for initial-value problems for difference and differ- ential equations. Journal of Computational and Applied Mathematics, 2333 1989.
Julian Besag. Comment on "Representations of knowledge in complex systems" " by Grenander and Miller. Journal of the Royal Statisticnl Sociely: Series B (Metholologioal) 56(4):549-55) 581, 1994.
Guy E. Blelloch. Prefix Sums and Their Applications. Technical Report CMU-CS-90-199, Carnegie Mellon University, School of Computer Science, 1990.
James Bradbury, Roy Frostig, Peter Hawkins, Malthew James Johnson, Chris Leary, Dougal Maclaurin, George Necula, Adam Paszke, Jake VanderPlas, Skye Wanderman-Milne, and Qiao Zhang. JAX: compos- able transformations of Python-NumPP programs, 2018.
C.G. Broyden. The convergence of class of double-rank minimization algorithms. IMA Journal of Applied Mathematics, 6(1):76-90, 1970.
Kyunghyun Cho, Bart van Merriënboer, Caglar Gulcehre, Dzmitry Bahdanau, Fethi Bougares, Holger Schwenk, and Yoshua Bengio. Learning Phrase Representations using RNN Encoder-Decodd Decoder for Sta- tistical Machine Translation. In Proceedings of the 2014 Conference on Empirieal Methods in Natural Language Processing (EMNLP), 2014.
C.M. da Fonseca. On the eigenvalues of some tridiagonal matrices. Journal of Computational and Applied Mathematics, 200(1):283 286, 2007.
Federico Danieli, Miguel Sarabia, Xavier Suau, Pau Rodríguez, and Luca Zappella. DeepPCR: Paralleliz- ing Sequential Operations in Neural Networks. In Advances in Neural Information Processing Systems (NeurIPS), 2023.
Federico Danicli, Pau Rodriguez, Miguel Sarabia, Xavier Suan, and Luca Zappella. ParaRNN: Unlocking parallel training of nonlinear RNNs for large language models. In International Conference on Learning Representations (ICLR), 2026.
Tri Dao and Albert Gu. Transformers are SSMs: Generalized models and efficient algorithms through structured state space duality, In International Conference on Machine Learning (ICML), 2024.
Tri Dao, Daniel Y. Fu, Stefano Ermon, Atri Rudra, and Christopher Ré. FlashAttention: Fast and memory- efficient exact attention with 10-awareness. In Advances in Newral Information Processing Systems (NeurIPS), 2022.
John E Dennis Jr and Robert B Schnabel. Numerical methods for unconstrained optimization and nonlinear equations. SIAM, 1996.
Haw-ren Fang and Yousef Saad. Two classes of multisccant methods for nonlinear acceleration. Numerical linear algebra with applications, 16(3):197 221, 2009.
Mónika Farsang and Radu Grosu. Scaling up liquid-resistance liquid-capacitance networks for efficient sequence modeling. In Advances in Neural Information Processing Systems (NewrIPS), 2025.
16

---

Published in Transactions on Machine Learning Research (4/2026)
Kayvon Fatahalian and Kunle Olukotu.. Lecture 8: Data-parallel thinking. Lecture slides, CS149: Parallel Computing, Stanford University, 2024.
Leo Feng, Frederick Tung, Mohamed Osama Ahmed, Yoshua Bengio, and Hossein Hajimirsadeghi. Were
RNNs All We Needed? orXiv, 2024.
Roy Friedman. A Simplified Overview of Langevin Dynamics. Blog post, 2022.
Martin J. Gander and Stefan Vandewalle. Analysis of the parareal time-parallel time-integration method. SIAM Journal Scientific Computing, 29(2):556-51) 2007.
V. A. Gasilov, V. F. Tishkin, A. P. Favorskii, and M. Yu. Shashkov. The use of the parallel-chord method to solve hydrodynamic difference equations. U.S.S.R. Computational Mathematics and Mathematical Physics, 21(3):178 192, 1981. ISSN 0041-5553. doi: (5111/0000155538188888188888
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. arXiv preprint "xXi:2552200177 2025.
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.
Xavier Gonzalez, Leo Kozachkov, David M. Zoltowski, Kenneth L. Clarkson, and Scott W. Linderman. Predictability enables parallelization of nonlinear state space models. In Neural Information Processing Systems (NeurIPS), 2025.
Tan Goodfellow, Yoshua Bengio, and Aaron Courville. Deep learning, volume 1. MIT Press, 2016.
Riccardo Grazzi, Julien Siems, Jörg KH Franke, Arber Zela, Frank Hutter, and Massimiliano Pontil. Un- locking State-Tracking in Linear RNNs Through Negative Eigenvalues. In International Conference Learning Represenlations (ICLR), 2025.
Sebastiano Grazzi and Giacomo Zanella. Parallel computations for Metropolis Markov chains with Picard maps, 2025.
Albert Gu and Tri Dao. Mamba: Linear-Tine Sequence Modeling with Selective State Spaces. In Conference on Language Modeling (COLM), 2024.
Albert Gu, Isys Johnson, Karan Goel, Khaled Saab, Tri Dao, Atri Rudra, and Christopher Ré. Combining Recurrent, Convolutional, and Continuous-time Models with Linear State-Space Layers. In Advances in Newral Information Processing Systems (NeurIPS), 2021.
Albert Gu, Karan Goel, and Christopher Ré. Efficiently Modeling Long Sequences with Structured State Spaces. In The International Conference on Learning Representations (ICLR), 2022.
Mark Harris, Shubhabrata Sengupta, and John D. Owens. Parallel prefix sum (scan) with CUDA. In Hubert Nguyen (ed.), GPU Gems 3, chapter 39, pp. 851-876 Addison-Wesley Professional, Upper Saddle River, NJ, August 2007.
Syeda Sakira Hassan, Simo Särkkä, and Ángel F Farcia-Fernández. Temporal parallelization of inference in hidden markov models. IEEE Transactions on Signal Processing, 69:4875-4887, 2021.
Kaiming He, Xiangyu Zhang, Shaoqing Ren, and Jian Sun. Deep residual learning for image recognition. In Proceedings of the IEEE conference on computer vision and pattern recognition, pp. 770-778, 2016.
Franz A. Heinsen. Efficient parallelization of ubiquitous sequential computation, 2023.
W Daniel Hillis and Guy L Steele Jr. Data parallel algorithms. Communications of the ACM, 29(12):
11 1986.
17

---

Published in Transactions on Machine Learning Research (4/2026)
Jonathan Ho, Ajay Jain, and Picter Abbeel. Denoising diffusion probabilistic models. In Advances in Neural Information. Processing Systems (NewrIPS), 2020.
Peter Holderrieth and Ezra Erives. Introduction flow matching and diffusion models, 2025. MIT course.
Graham Horton, Stefan Vandewalle, and P Worley. An algorithm with polylog parallel complexity for solving parabolic partial differential equations. SIAM Journal on Scientific Computing, 16(3):5311544 5411, 1995.
Amber Hu, Henry Smith, and Scott Linderman. SING: SDE Inference via Natural Gradients. In Aduances in Neural Information Processing Systems (NeurIPS), 2025.
Michael F Hutchinson. A stochastic estimator of the trace of the influence matrix for Laplacian smoothing splines. Communications in SatisliicsSimmulation and Computation, 18(3):1059- 1001-0999--66,,, 1989.
Casian Tacob, Hassan Razavi, and Simo Särkkä. A parallel-in-time newton's 's method-based ode solver. arXiv preprint trXfo:2511.01465, 2025.
Matthew J Johnson, David K Duvenaud, Alex Wiltschko, Ryan P Adams, and Sandeep R Datta. Composing graphical models with neural networks for structured representations and fast inference. Advances in Newrul Information Processing Systems, 2016.
R. E. Kalman. A new approach to linear filtering and prediction problems. Journal of Basic Engineering, 82(1):35-45 45, 1960.
Craig Kapfer, Kurt Stine, Balasubramanian Narasimhan, Christopher Mentzel, and Emmanuel Candès. Marlowe: Stanford's GPU-based Computational Instrument. Zenodo, 2025. Version 0.1.
Carl T Kelley. Iterative methods for linear and nonlinear equations. SIAM, 1995.
Patrick Kidger and Cristian Garcia. Equinox: neural networks in JAX via callable PyTrees and filtered transformations. CoRR, abs/2111.0025 2021.
Najoung Kim and Sebastian Schuster. Entity tracking in language models. arXiv preprint orXiv:2305.02363, 2023.
Volodymyr Kyrylov. Accelerated scan, 2024. URL http:://gttttoocooopprggg///aaee GitHub repository.
Richard E Ladner and Michael J Fischer. Parallel prefix computation. Journal of the ACM (JACM), 27(4): 831-838, 1980.
Chieh-Hsin Lai, Yang Song, Dongjun Kim, Yuki Mitsufuji, and Stefano Ermon. The principles of diffusion models. arXiv preprint arXiv:2510.21890 2025.
Sivaramarrishnan Lakshmivarahan and Sudarshan K Dhall. Parallel computing using the prefic problem. Oxford University Press, 1994.
Paul Langevin. On the Theory of Brownian Motion. American Journal of Physics, 651):1107 1081, 1908. English translation, introduced by D. Lemons and translated by A. Gythiel. Original: C. R. Acad. Sci. 146, 530-534 (1908).
Kenneth Levenberg. A method for the solution of certain non-linear problems in least squares. Quarterly of Applied Mathematics, 2:164-168, 1944.
Yi Heng Lim, Qi Zhu, Joshua Selfridge, and Muhammad Firmansyah Kasim. Parallelizing non-lineat sequen- tial models over the sequence length. In International Conference on Learning Representations (ICLR), 2024.
Scott W Linderman, Peter Chang, Giles Harper-Donnelly, Alcyna Kara, Xinglong Li, Gerardo Duran-Martin, and Kevin Murphy. Dynamax: A Python package for probabilistic state space modeling with JAX. Journal of Open Source Software, 10110)):7069, 2025.
18

---

Published in Transactions on Machine Learning Research (4/2026)
Jacques-Louis Lions, Yvon Maday, and Gabriel Turinici. A "parareal" in time discretization of PDE's. Comples Rendus de I' l'Acudémie des Sciences Series T- Malhemalics, 332(():6611660) 2001.
Bingbin Liu, Jordan T. Ash, Surbhi Goel, Akshay Krishnamurthy, and Cyril Zhang. Transformers learn short- cuts to automata. In Proceedings of the International Conference on Learning Representations (ICLR), 2023.
Jianrong Lu, Zhiyu Zhu, and Junhui Hou. Parasolver: A hierarchical parallel integral solver for diffusion models. In International Conference on Learning Representalions (TCLR), 2025.
Donald W. Marquardt. An algorithm for least-squares estimation of nonlinear parameters. Journal of the Society for Industrial and Applied Mathematics, 11(2):431 441, 1963.
Eric Martin and Chris Cundy, Parallelizing Linear Recurrent Neural Nets Over Sequence Length. In International Conference on Learning Representations (ICLR), 2018.
William Merrill, Jackson Petty, and Ashish Sabharwal. The Illusion of State in State-Space Models. In International Conference on Machine Learning (ICML), 2024.
Kevin P Murphy. Prolabilistic machine learning: Advanced lopics. MIT Press, 2023.
Maxim Naumov. Parallel complexity of forward and backward propagation. arXiv preprint orXi:: 2017.
Hubert Nguyen (ed.). GPU Gems 3. Aidison-Wesley, Professional, August 2007. ISBN 978-0321515261.
Jorge Nocedal and Stephen Wright. Numerical Optimization. Springer, 2 edition, 2006.
James M Ortega and Werner C Rheinboldt. Iterative Solution of Nontinear Equations in Several Variables. Academic Press, New York and London, 1970. Republished by SIAM in 2000.
Antonio Orvieto, Samuel L 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 (ICML), 2023.
Adam Paszke, Sam Gross, Francisco Massa, Adam Lerer, James Bradbury, Gregory Chanan, Trevor Killeen, Zeming Lin, Natalia Gimelshein, Luca Antiga, Alban Desmaison, Andreas Köpf, Edward Z. Yang, Zach DeVito, Martin Raison, Alykhan Tejani, Sasank Chilamkurthy, Benoit Steiner, Lu Fang, Junjie Bai, and Soumith Chintala. PyTorch: An Imperative Style, High-Performant Deep Learning Library. Advances in Neural Information Processing Systems (NewrIPS), 2019.
E. Rauch, Tung, and C. T. Striebel. Maximum likelihood estimates linear dynamic systems. AIAA Journal, 3(8):1445-144 1965.
Simo Särkkä and Ángel F. actions on Automatic Control, 66(1):299-366, 2021. doi: 10.109//AA.2200.297666.
Temporal parallelization of bayesian smoothers. IEEE Trans-
Simo Särkkä and Lonnart Svensson. Bayesion filtering and smoothing, volume 17. Cambridge University Press, 2023.
Felix Sarnthein. Linear recurrences accessible everyone. In ICLR Blogposts, 2025.
Mark Schöne, Babak Rahmani, Heiner Kremer, Fabian Falck, Hitesh Ballani, and Jannes Gladrow. Implicit Language Models are RNNs: Balancing Parallelization and Expressivity. In International Conference on Machine Learning (ICML), 2025.
Nikil Roashan Selvam, Amil Merchant, and Stefano Ermon. Self-Refining Diffusion Samplers: Enabling Parallelization via Parareal Iterations. In Advances in Neural Information Processing Systems (NeurIPS), 2024.
19

---

Published in Transactions on Machine Learning Research (4/2026)
Andy Shih, Suneel Belkhale, Stefano Ermon, Dorsa Sadigh, and Nima Anari. Parallel Sampling of Diffusion Models. In Advances in Neural Information Processing Systems (NewrIPS), 2023.
Jimmy T.H. Smith, Andrew Warrington, and Scott W. Linderman. Simplified State Space Layers for Se- quence Modeling. In International Conference on Learning Representations (ICLR), 2023.
Jascha Sohl-Dickstein, Eric Weiss, Niru Maheswaranathann, and Surya Ganguli. Deep unsupervised learning using nonequilibrium thermodynamics. In International Conference on Machine Learning (ICML), 2015.
Yang Song and Stefano Ermon. Generative Modeling by Estimating Gradients of the Data Distribution. In Advances in Newral Information Processing Systems (NeurIPS), 2019.
Yang Song, Chenlin Mong, and Stefano Ermon. Mintnet: Building invertible neural networks with masked convolutions. Advances in Newral Information Processing Systems (NeurIPS), 2019.
Yang Song, Chenlin Meng, Renjie Liao, and Stefano Ermon. Accelerating Feedforward Computation via Parallel Nonlinear Equation Solving. In International Conference on Machine Learning (ICML), 2021a.
Yang Song, Jascha Sohl-Dickstein, Diederik P Kingma, Abhishek Kumar, Stefano Ermon, and Ben Poole. Score-Based Generative Modeling through Stochastic Differential Equations. In International Conference on Learning Representations (ICLR), 2021b.
Harold S. Stone. An efficient parallel algorithm for the solution of a tridiagonal linear system of equations. Journal of the ACM, 20(1):27-33 38, 1973.
Simo Särkkä and Lennart Svensson.
and line-search extended Kalman smoothers. In ICASSP 2020 2020 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), pp. 5875 5879. IEEE, 2020.
Zhiwei Tang, Jiasheng Tang, Hao Luo, Fan Wang, and Tsung- Chang. Accelerating Parallel Sampling of Diffusion Models. In International Conference on Machine Learning (ICML), 2024.
Aleksandar Terzić, Nicolas Menet, Michael Hersche, Thomas Hofmman, and Abbas Rahimi. Structure sparse transition matrices enable state tracking in state-space models. In Advances in Neural Information Processing Systems (NeurIPS), 2025.
Ashish Vaswani, Noam Shazeer, Niki Parmar, Jakob Uszkoreit, Llion Jones, Aidan N Gomez, Lukasz Kaiser, and Illia Polosukhin. Attention is All You Need. In Advances in Neural Information Processing Systems (NeurIPS), 2017.
Homer F Walker and Peng Ni. Anderson acceleration for fixed-point iterations. SIAM Journal on Numerical Analysis, 49(4):1715 1735, 2011.
Fatemeh Yaghoobi, Adrien Corenflos, Sakira Hassan, and Simo Särkkä. Parallel square-root statistical lincar regression for inference in nonlinear state space models. SIAM Journal on Scientific Computing, 47(2): B454, 2025.
Songlin Yang, Bailin Wang, Yikang Shon, Rameswar Panda, and Yoon Kim. Gated Linear Attention Trans- formers with Hardware--Efiiient Training. In International Conference on Machine Learning (ICML), 2024.
David M. Young. Ilerative Solution of Large Linear Systems. Elsevier, 2014. ISBN 975-0-12-773050-..
Riccardo Zattra, Giacomo Baggio, Umberto Casti, Augusto Ferrante, and Francesco Ticozzi. Context- selective state space models: Feedback all you need. arXiv preprint arXiw:2510.11 140077 2025.
Yixiu Zhao and Scott Linderman. Revisiting structured variational autoencoders. In International Confer- ence on Machine Learning (ICML), 2023.
David M. Zoltowski, Skyler Wu, Xavier Gonzalez, Leo Kozachkov, and Scott W. Linderman. Parallelizing MCMC Across the Sequence Length. In Advances in Neural Information Processing Systems (NeurIPS), 2025.
20

---

Published in Transactions on Machine Learning Research (4/2026)
A Author Contributions and Acknowledgments
Author Contributions XG, EKB, and SWL conceived the project. XG and SWL proved Proposition 1 and Proposition 2. XG provided most of the content of Section 3, with imput from all authors. XG and JWL proved Proposition 4. HDL developed, ran, and wrote up the state tracking experiments. XG, EKB, HDL, JWL, and KAW all wrote substantial code in various exploratory aspects of the project. DZ and CR contributed important insights about the relationships between algorithmu and computational efficiency. LK contributed important insights about rates of convergence and the importance of stability. All authors made important contributions conceptualization, brainstorming, exploration, writing, and editing. XG, EKB, and SWL wrote the manuscript with input from all authors. SWL supervised the project.
Acknowledgments We thank Nicolas Zucchet, Francois Chaubard, Libby Zhang, Henry Smith, and lan Christopher Tanoh, well the anonymous reviewers on Open Review, for helpful feedback. X.G. acknowl- edges support from the Walter Byers Graduate Scholarship from the NCAA. J.W.I. was supported by the Department of Energy Computational Science Graduate Fellowship under Award Number DE-SC0023112. S.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 Com- puting Center for providing computational resources and support. Additional computations were performed Marlowe (Kapfer et al., 2025), Stanford University's GPU-based Computational Instrument, supported by Stanford Data Science and Stanford Research Computing. The authors have no competing interests.
B The Parallel Scan: A Gentle Introduction
B.1 very simple example: multiplying matrices
The parallel scon (Stone, 1973; Blelloch, 1990), also known as the associative scan and, colloquially, pscan, is a well-known primitive in the parallel computing literature (Hillis Steele Jr, 1986; Ladner Fischer, 1980; Lakshmivarahan Dhall, 1994). The core idea of the parallel scan is a andoo--and-oonquer algorithm. We illustrate this point in the simple example of multiplying scries of matricos together.
Simple Example (Multiplying Matrices): Given a series of square matrices A,, A2. ... Ar....,AA-LL. 1, Ar, compute their produet AT A2A1. The simplest way to carry out the matrix multiplication is sequentially: first compute A1, then compute A2A,, then compute AsA2A1, and so on. Such approach takes (TT) time.
A core insight of the parallel sean is that matrix multiplication closed; that is, if A, E RDxD and A, e DxD, then ArAg DDD, Thus, matrix products can be computed recursively in pairs, as illustrated in Figure 5.
A1
A2
Aa
A
A2A1
A4
AsAsAzA,
Figure Parallel Scan for Matrix Multiplication. We illustrate a divide-and- aivide-and-conqque approach compute the product A442A2A1 Note that this divide-and-conquer approach naturally leads to O(log T) depth.
Note that we have the matrices act left-multiplication over sequence length, because this the most common way write matrix vector products.
21

---

Published in Transactions on Machine Learning Research (4/2026)
Because of the divide-and-comquae (binary-tree-ikke) nature of this approach to multiplying matrices, with O(T) processors, the time needed to get the matrix product is only Ollog 77). This simple example illustrates the core intuition behind the parallel scan: a closed operation leading to a anc-ouu approach that parallelizes a computation so that it takes sublincar time. However, there are two additional details of the parallel associalive sean that we should address: arbitrary binary associative operators and closure; and getting intermediate products.
B.2 Detail #1: Parallel scans for arbitrary binary associative operators
Matrix multiplication an associative operator, as A1(A2A1) = (A1 In general, consider a binary associative operator &, which would satisfy 9s @ (g2 2 q1) = (qs @ 92) @ q1 Now, let us further assume that this binary associative operator is closed: Definition 2 (Closure). A binory associative operator is closed over a set S if it satisfies the property:
q1 E S2 92 E + 92 9 ES.
(18)
If is closed, then we can again use a parallel scan compute the cumulative product of the operands. A wide range of binary associative operators are closed, and can thus be parallelized with the parallel scan. Some examples include:
Scalar addition:
The fact that addition of scalars (and vectors) is closed allows cumulative sums to be computed with the parallel scan algorithm. In this instance, it is also known as the prefir sun algorithm. Clearly, addition is associative and closed, and so summing a series of scalars can be done with divide-and-comquer approach.
Composition of affine functions, in LDSs. Consider the affine function f((x) A bi- Notice that the composition of affine functions is also affine, as fy(f(z)) A, + (6, + A,6.). Thus, if we represent the operands as ordered pairs (A,, b.) and (Aj.by), we can write the associative operator & for the composition of affine functions as
(A,, bi) @ (A),b,) = (AjA,, + A,bi).
Thus, we observe that in this setting, & is closed. We also should check that @ is associative: we can do so with either elementary algebra, or by observing that function composition is associative.
This observation that composition of affine functions can be parallelized with the associative scan is what lets us parallelize LDSs. The parallelization of LDSs is what allows for parallel computation in many important architectures in sequence modeling, such as linear RNNs (Martin & Cundy, 2018; Orvieto al., 2023), deep SSMs (Smith et al., 2023; Gu & Dao, 2024), and nonlinear RNNS (Lim et al., 2024; Gonzalez el al., 2024; Farsang & Grosu, 2025; Danieli et al., 2026). This parallel scan for LDSs is the core primitive uniting the fixed-point methods discussed in this paper.
Kalman filtering and smoothing: Parallel scans can also be utilized in probabilistic modeling. A standard probabilistic model is the linear Gaussian state space model (LGSSM), where the latent variables x; follow linear dynamics with Gaussian noise, and emit obscrvations ye with linear readouts with Gaussian noise (Murphy, 2023; Särkkä & Svensson, 2023). See Figure 6.
Two canonical inferential targets in the LGSSMM are the filtering distributions, p(ax. y):)), and the smoothing distributions, p(æ, B1.T). The Kalman filter (Kalman, 1960) and Rauch- Strieeet (RTS) smoother (Rauch et al., 1965) obtain the filtering and smoothing distributions (respectively) in an LGSSM. The Kalman filter makes a single pass forward in time to get the filtering distributions, while the RTS smoother then makes an additional pass backwards in time to get the smoothing distributions. Both the Kalman filter and RTS smoother would seem to be inherently sequential algorithms, requiring O(T) time. However, Särkkä & Carcia-Fernááne (2021) demonstrated that
6The parallel scan is used in nonlinear RNNs via the iterative fixed-point methods discussed in this paper, Newton and Newton iterations, see Algorithm 1.
22

---

Published in Transactions on Machine Learning Research (4/2026)
30
J/2
F3
1/1
1/2
y3
Figure 6: A linear Gaussian state space model (LGSSM): The LGSSM consists of latent variables and observed variables ye- The generative model of the LGSSM consists of dynamics Z++1 ~ N (Ar, Q) and emissions 1141 ~ N (Hæ141, Rt.
the Kalman filter and RTS smoother can also be parallelized over the sequence length via the construction of custom binary associalive operators and a parallel scan. While leave the details (2021), we note that it is intuitively plausible to
of this construction to Särkkä & be able to parallelize filtering and smoothing in an LGSSM with a parallel scan because
the dynamical backbone is an LDS, for which we have a parallel scan;
since everything is lincar and Gaussian, all distributions remain Gaussian, hinting at closure; and
we can combine p(a,, (o B1ss) with p 1s) to obtain p( (2n Vrze). suggesting divide- and-conquer strategy.
These parallel filtering and smoothing algorithms are useful in machine learning, allowing for par- allelization of structured variational autoencoders (Johnson et al., 2016; Zhao & Linderman, 2023). Similar approaches also work for Hidden Markov Models (Hassan et al., 2021) and for computing log-normalizing constants (Hu et al., 2025).
Parallelizing fixed point iterations: applicable to the types of parallel fixed-point itcrations that are the focus of this paper. In partic- ular, the ELK algorithm (Gonzalez et al., 2024) Evaluating Levenberg-Marquardde via Kalman stabilizes the Newton iterations (Lim et al., 2024) we discuss this paper using the Levenberg- Marquardt (trust-region) method (Levenberg, 1944; Marquardt, 1963). The trust-region updates are able to be computed using a parallel scan because they are equivalent to Kalman smoother in an appropriately constructed L.SSMM (Särkkä & Svensson, 2020).
Finally, these parallel filtering and smoothing are directly
What about arbitrary function composition? The astute reader might note that the composition of functions, i.e. f2, always binary associative operator. So, why do we have all these special cases of parallel scans, and not simply one parallel scan for the composition of arbitrary functions f.?
The reason to have many different parallel scans is preciscly the importance of having the binary associative operator be closed. In all the previous examples, the binary associative operator @ satisfies Definition 2, letting us easily store combinations of operands q & and so employ a divide-and-conquer technique.
While we could consider some gigantic function space F, for which function composition would be closod, the practical question then becomes: how would we store the combinations of operands? do not have some compact representation for elements of F, then cannot parallel scan in practice, even though the parallel scan may seem applicable in theory.
B.3 Detail #2: Obtaining the intermediate terms in the product
When we evaluate lincar dynamical system, we often do not want only the final term TT, but also all the intermediate terms X i.c the full trajectory. So far, the parallel scan as presented would only yield the final term Fr as well intermediate terms that happened to be powers of 2, i.e. 21, F., 38, etc.
23

---

Published in Transactions on Machine Learning Research (4/2026)
However, the parallel scan can be easily adjusted to obtain all the intermediate terms well. Let return to our very simple example of matrix multiplication to illustrate, in particular the setting where T = 8. We will denote the individual matrices as A1, A2, As, and their products as Asa, i.e. A = AgA5.
The first phase of the parallel scan the up-sweep, and takes (T) iterations and O(T) memory. We start multiplying adjacent pairs of matrices together, going, for example", from As to A78 to As:8 to A
Then, in the down-sweep, we fill in the missing products obtain all the cumulative products Aue for stt S Intuitively, the down-sweep also takes Ollo T) iterations, for the same reason that any natural number 7 can be represented using + (og2(7) digits in binary.
Table 2: Up-sweep for multiplying AI, .g....AA. Position
Position
Position
Position
Position
Position
Position
Position
A1
A2
Ag
A4
As
Ag
A7
As
A
A1.2
Ag
A3:4
A5
A5:6
A7
A7:8
A1
A1
A3
A14
A5
A5.6
A7
Asis
A
A1.
Ag
A14
A5
A5.6
A7
A
Table 3: Down-sweep for multiplying A1, Ag...... Ag.
Position 1
Position 2
Position
Position 4 Position
Position
Position
Position
A1
A1.2
A3
A14
As
A1.
A7
A
A
A1.2
A1.3
A14
A
A
Aur
A
Thus, together, the up-sweep and the down-sweep of the parallel scan run in Ollog I)) time O(T) proces- sors, and at the end of this algorithm, we get all of the intermediate products (the "prefix sums").
B.4 Implementation considerations
This gentle introduction provides the main ideas and intuition of the parallel scan: closed binary associa- tive operators leading a divide-and-cooqque algorithm to leverage parallel processors to compute sequential compositions sublinear time. However, there are also a host of implementation details for using the par- allel scan when programming on accelerated hardware like GPUs (Harris et al., 2007; Yang el al., 2024; Sarnthein, 2025). For example, the presence of general-purpose parallel scan is, as of the time of writing, major difference between JAX (Bradbury et al., 2018) and PyTorch (Paszke et al., 2019), two leading Python libraries for deep learning. JAX has a general purpose parallel scan (jax. lax. associative scan) as fundamental primitive, which allows for implementation of a wide range of parallel scans. For example, dynamax, a JAX library for probabilistic state space modeling (Linderman et al., 2025), implements the parallel filtering and smoothing algorithms from Särkkä García- Ferfándezz (2021). In contrast, PyTorch currently has only torch. cumsum, which the parallel scan where the binary associative operator addi- tion', and torch. cumprod (for scalar multiplication). This difference is why we implement the experiments in this paper in JAX. This lack of general purpose parallel scan in PyTorc has also led to the custom devel- opment of highly-optimized, hardware-awarr custom CUDA kernels for parallel scans. These custom parallel scans appear most prominently in Mamba (Gu & Dao, 2024), a leading SSM for language modeling, and ParaRNN (Danieli et al., 2026), which applies the Newton iterations discussed in this paper to 7B parameter nonlinear RNNs to achieve strong language modeling performance. There also exist useful implementations of parallel scans for scalar /diagonal L.DSs in PyTorch such as Kyrylov (2024), which we used to implement quasi-Newer iterations Pororch this repo: https: htpps:///itthule htt.h/....coo//dddmmmaaaaaaa/a/a
See Position 8 of Table
88cc the last of Table 3.
9 Althooggh Heinsen (2023) shows that clever uses cuneum can parallelize scalar/dlagonal LDSs, of the type that used in quasi- Newton iterations (eq. (6)).
24

---

Published in Transactions on Machine Learning Research (4/2026)
C Further Discussion of Convergence Analysis
In this appendix, we provide a deeper discussion of the details, intuitions, and limitations of the convergence analysis presented in Section
C.1 Formula for approximate Jacobian
For completeness, we provide the formula for the block-bidiagona approximate Jacobian J(xIT) E RIDxID that the ixietpoinnt methods considered this paper give rise to.
ID
0
0
0
A2(m1)
ID
0
0
0
ID
[xx:)) =
(19)
0
ID
0
0
(r(rr-1)
ID/
For Newton's method AN, the resulting (x(x)))) = (x::), where I(x1:) is defined in eq. (11). For Picard iterations, Jx(x)r)) takes the form
ID
0
0
0
-ID
ID
0
0
0
0
- In
In
0
0
JP(XIT)=
(20)
1000 0
ID
0
0
Ip
Tp
For Jacobi iterations, J)(x).)) is always the identity matrix LTD.
C.2 Limitations of Proposition 3
Proposition 3 only guarantees decrease in the error when the iterate (1)) already in basin of decrease Bd given by
1 [(x 11, -(x - (( L [(xir))
Bp = XLT |e(x < 2-
2
However, since we know from Proposition of Gonzalez et al. (2024) that all the fixed-point algorithms considered in this paper must eventually converge, we know that the iterates x must all eventually enter this basin of decrease Bn if B, * 0. For this reason, Proposition 3 provides helpful intuition about which fixed-point algorithms are useful for which dynamical systems.
For example, let us define the basin of linear rate B2 to comprise those [(xx.r) - ((x [/ (1 L.. the expression linear in [|o right side of equation 13 dominates the expression quadratic in [|e(s) ll.. It follows that B,, is given by
X where
2 2[[(xx))) J(x1x)
By
X1.r [|((x))))22
L
17 Bp n Br, it follows that the norm of the error guaranteed to decrease by factor of J [ -1 J [xx Moreover, [ approaches zero, the guaranteed factor of decrease approaches the value given by eq. (15).
Therefore, when x (1)
(1)
25

---

Published in Transactions on Machine Learning Research (4/2026)
C.3 Extended discussion about intuitions from the rate of convergence
In this appendix, we elaborate the inluitions about the rate of convergence presented in Section 3.3.
C.3.1 Intuitions from Jacobian approximation error
First, elaborating on the section discussing [[xxrr) (x1:T) we provide a proof of Lemma 1.
Proof of Lemma 1. Plugging in the functional forms of J() and J(), if we define E = Act1(2t) Ac+1(
then
E2
E3
0
[xIT)) - J(xIT))
Fr
0
The spectral norm of matrix M is equal to the square root of the largest eigenvalue of M M. Defining M = ((x.T) - J(xIT), we see that
0
0
0
Eg
E2
0
0
0
Fgg
F3
0
M M
0
0
0
E
ET
0
0
ET ET
Since M T M block-diagonal matrix, its eigenvalues are equal to the union of the eigenvalues of each of the blocks E, F4. Thus, it follows that the maximum eigenvalue of M M is equal to the maximum of all the eigenvalues of all the matrices E, E,, and so the maximum singular value of ((xir) J (x.r)) given by maxigesr-i|lAuqde (z.)) Ax41(2)))12
C.3.2 Intuitions from the norm of the inverse of the approximate Jacobian
Next, we claborate on our bounds for [(x1:r)
Recall that j(x1r) takes the form (Dao & Gu, 2024; Gonzalez et al., 2025)
Tn
0
0
0
-'x1a))===
In
0
0
Åg
Tp
0
A443A2
ÁqA3
A4
Tp)
shown above for T = 4.
Theorem of Gonzalez et al. (2025) controls [(x terms of p == sUP251STwER? [x)((2) and p :: infgsIsT.ERR [x)]( as
1 1- p" max(1, (((i))))111 <
1-p"
for 7 F We would therefore expect that fixed-point methods that give rise to unstable L.Dss with transition matrices having spectral norms much larger than one should have slower rates of convergence. Methods that give rise to unstable LDSs suffer from numerical blowup, especially for large T.
26

---

Published in Transactions on Machine Learning Research (4/2026)
Jacobi
Picard
Convergence Speed
100
10
0.8
0..7
80
10 1.
0.6
Iterations to reach |ee
|lc80
60-
Jacobi
0.55
10
Picard
0.4
40
10
0.3
0.2
10
0
25
100
0
25
75
100
0.2
0.4
0.6
0.8
Iteration (5)
Iteration (4)
Figure 7: Comparing Picard and Jacobi iterations on diagonal LDS. For the underlying dynamical (i) system I(+1 az; with T 100, plot the norm of the error c1.7 for Jacobi (Left) and Picard iterations
(Center). In the (Right) panel, we also show the number of iterations needed for the norm of the error to below 1 x 10-5
analytically, For Jacobi iterations, Jı = For Picard iterations, the expression for [ is more complicated, but it scales as O(T) (see Lemma 2 in Appendix D.1).
Morcover, in the special cases of Jacobi and Picard itcrations, we can compute [
T
Because
J ll > الار|| for large T, the formula for given by eq. (15) yields the following expectation:
In settings where the A4t1 from Picard vs. Jacobi iterations approximates the true dynamics Jacobian 0f+t1/oz, equally well, we expect Jacobi iterations to converge more quickly because | [Jp'1 2--
We test this hypothesis in the next section with simple simulation designed to show how Proposition 3 provides helpful intuition about the convergence rates of different fixed-point methods.
C.4 simulation using Proposition 3 to distinguish between Jacobi and Picard iterations
We demonstrate the helpfulness of the intuitions stemming from Proposition in simple simulation. We consider the LDS æz-1 az2, for 2, IRa Because this is an LDS with diagonal dynamics, both the Newton and quasi-Newton itcrations considered in this paper converge in one iteration. However, this simulation is useful for comparing Jacobi versus Picard iterations. This comparison is particularly fruitful in light of the formula for given by eq. (15) and Lemma. because, in this setting,
[J J l2 [ 2a 1
However, || = 1, while [ scales lincarly with T. Therefore, when comparing the number of Jacobi iterations needed to converge when the dynamics are multiplication by a to the number of Picard iterations needed to converge when the dynamics are multiplication by - aa, we expect fewer Jacobi iterations should needed than Picard iterations, as 7J < 7p.
We observe precisely this behavior in Figure 7. For a 0.5, when رقا ر p we see that Jacobi iterations converge in far fewer iterations than Picard iterations. Moreover, when comparing the behavior of
27

---

Published in Transactions on Machine Learning Research (4/2026)
Jacobi for simulating f(+1(x)) az, with Picard for simulating f(+1(1z) (1-a)x,, observe that Jacobi iterations always converge laster. However, when comparing for the same value of a, we see that Picard can be faster than Jacobi when a is closer to one, This behavior makes sense, because in those settings the true Jacobian 0f01 is closer to ID than to O.
Moreover, observe that in this setting, the error eI:T for Jacobi iteration shows clear linear convergence rale, as predicted by Proposition 3. The slope of the norm of the errors of the Jacobi iterates should be la) by eq. (15) and Lemma 1, and in fact those values are exactly the slopes of the lines in Figure 7 (Left pancl).
(8)
Theoretical Details: Additional Proofs
D.1 Lemma
Lemma 2. Let Jp be as in eq. (20). Then
|J sin
202++1)
By the small-angle opprorimation for sine, [J p ll2 scales as (T)).
Proof. Consider
ID
1a
ID
la
ID
2Ip
2ID
2ID
K : Jp
ID
2
3ID
3ID
ID
2ID
3ID
TID
We know that Amax(K))/ |Jp' 1-- Since the spectrum of K is equivalent to the spectrum of M (just with all eigenvalues having multiplicity D). Therefore, seek to find the spectrum of M E RTTT
Kronecker product M @ ID. where Ma, min
However, the spectrum of M is known in the literature. For example, Theorem 2.1 of da Fonseca (2007)
shows that if T 2 3, then the cigenvalues {A=
T-1 of M are given by
2k + 1
-1
A
008
77 27 +
2
1
2k +
-2
sin
2(2T 1) m)))
where the second equality comes from the half-angle formula. We observe that the largest eigenvalue is therefore Ao and so the result follows after we take a square root.
E Experimental Details and Additional Experiments
We implemented our experiments using the Equinox library (Kidger & Garcia, 2021) in JAX (Bradbury al., 2018). In our experiments that report wall-clock time, iterations (7oltowski et al., 2025) that estimates the diagonal using the Hutchinson estimator (Hutchinson, 1989); see Section 3.4 of Zoltowski et al. (2025) for details.
stochastic implementation of quasi-Newton
The stopping criterion we use for deciding when the fixed-point iterations in eq. (2) have converged is based on the merit function ((xIr), which is defined as:
r(x1:T) = [ ff)=)(x)) x2 f(11), f(z2). ... FT f(xr-1)1
[(xim) = (000)=]ll(xxxx (1((
28

---

Published in Transactions on Machine Learning Research (4/2026)
In our experiments use tolerance of 5 x 10-4, that is we terminate the fixed-point iterations when iterate
i satisfies
C(x 1.T) 55x10
E.1
Experimental Details for Appendix C.4
For each scalar LDS with scalar multiplication a, we use state size of D =2 length of 7= = 100. We consider 10 random seeds, where the randomness controls the initial state Tp- We initialize x i.T at all zeros. We plot the median number of steps to convergence in Figure 7.
2, i... E R2, and a sequence
(0)
E.2 Experimental Details for Section 4.1
We use 10 random seeds, where the randomness controls the sequence of the S₁ group elements, Each seed uses a batch size of 16, i.c. 16 different Ss word problems are evaluated for cach run. For cach scod, we time how long it lakes for runs of the fixed point solver (sequential, Picard, quasi-Newlon, or Newton) to evaluate, and record this mean wall-clock time. In Figure 2, then plot the median of these 10 mean wall-clock times. We run an H100 with 80GB of onboard memory.
E.3 Experimental Details for Section 4.2
The experiment depicted in Figure 3 closcly follows the experiments in Section 4.1 of Lim et al. (2024) and Section 6.1 of Gonzalez et al. (2024, Section 6.1). We use 10 random seeds, where the randomness cont the initialization of the GRUs, the random inputs to the GRUs, and Rademacher variables used in the stochastic implementation of quasi-Newton iterations (Zoltowski et al., 2025). The GRUs were initialized following the standard initialization practice in Equinox. Each seed uses batch size of 16, i.e. 16 different GRU trajectories are evaluated for each run. For each seed, we time how long it takes for 5 runs of the fixed point solver (sequential, Picard, quasi-Newton, or Newton) to evaluate, and record this mean wall-clock time. Figure 3, we then plot the median of these 10 mean wall-clock times. We run on an H100 with 80GB of onboard memory.
To demonstrate the different values of the diff(-) operator for quasi-Newton, Jacobi, and Picard iterations in this setting, we consider the setting D and T 1000. In Figure 8, we plot variety of quantities relevant for the 7 (ef. eq. (15)) for ten random seeds. We observe that lower values of (i.e., faster rates of asymptotic linear convergence) coincide with fewer fixed-point iterations needed in Figure 3. We observe that both (A,, and (ifl() are both below one always, which corresponds to their fast rates of convergence demonstrated in Figure 3. In contrast, (AAp) is always greater than one, which corresponds to the slow rates of convergence of Picard iteration the experiment depicted in Figure 3.
E.4 Experimental Details for Section 4.3
The potential used for the experiment depicted in Figure 4 is the negative log probability of a mixture of Gaussians. Fach Gaussian is D-dimensional, and has a random covariance matrix drawn from a Wishart distribution. For each fixed-point method, we ran 10 random seeds, where the randomness controls the randomly chosen covariance matrices for the mixture of Gaussians and the random inputs for Langevin dynamics. Fach seed uses a batch size of 16, i.e. 16 different Langevin trajectories are evaluated for each run. We use step size = 1 10-5 for the discrete Langevin steps. For each sced, we time how long takes for 5 runs of the fixed point solver (sequential, Picard, quasi-Newton, or Newton) evaluate, and record this mean wall-clock time. In Figure 4, we then plot the median of these 10 mean wall-clock times. We run an H100 with 80GB of onboard memory. In order to get convergence for Newton iterations for longer sequence lengths, we had to run with increased precision, using the highest option for the jax_djxeaaaa fox__dffalt__atteee ctattt_aamm___peiiiists pprrilsppeeiioonn flag. See Appendix F.5 for further discussion of the importance of numerical precision for implementation of these LDS-based fixed-point methods.
We also include a small additional experiment in Figure 9: instead of varying the state size dimension, we instead vary K, the number of Gaussians that make up the potential & determining the Langevin dynamics.
29

---

Published in Transactions on Machine Learning Research (4/2026)
2.0
3000
1.5
2000
1.0
1000
0.5
0.0
0.1
0.2
0.3
0.4
1.26
1.28
1.30
1.32
1.34
1.36
||
1|2
[ [JaNlla
3000
4
2000
2
1000
0.25
0.50
0.75
1.00
1.25
1.50
1.75
0.8
0.9
1.0
1.1
||
A411l2
Jacobi
Picard
Quasi-Newton
Figure 8: Understanding the convergence rates in Figure 3. In the setting of the GRU experiment for D = 8 and T = 1000, we plot relevant quantities for understanding the convergence rates of different methods over 10 random seeds. (Top left.) For all T = 1000 time steps and over all 10 random seeds, we plot the spectral norm of the approximate Jacobian for the quasi-Newton iterations we consider in this paper, iee diag(of++/o<<x;;)) (Top right.) For each of the 10 random seeds, we plot [[JoN(xt,)) observe that they are always larger than one. (Bottom left.) We plot the difference between approximate Jacobians and true dynamics Jacobians over all time steps and seeds for quasi-Newton, Jacobi, and Picard iterations. We observe that this difference for Picard iterations is always larger than one, and so would intuitively expect Picard iteration to be very slow for parallelizing GRUs, This behavior is precisely what we see Figure (Bottom right.) Across the 10 random seeds, we plot the value of y for Jacobi and quasi-Newton iterations (Picard would be O(T) and so is not shown). Because [T]( PJ's are equivalent to the maximum values from the differences in (bottom left) over the 10 random seeds. However, since (top right) shows that [|Jox(xi.r) than the values of ?u in (bottom right). In summary, because the values of y are smaller than the values of TQN, we would intuitively expect Jacobi to converge in fewer fixed-point iterations, which is exactly what we observe in Figure 3.
We
= 1, the 10
1, observe that the values of "QN
larger
We observe qualitatively similar results to the experiment shown in Figure 4: viz, that Picard and quasi- Newton iterations enjoy similar convergence rates and wall-clock time in settings where the Jacobian is well-approximated by the identity matrix.
Finally, we include an experiment designed to show how the density of the Jacobian favors Newton iterations Picard iterations. In the experiment depicted in Figure 10, again parallelize discretized Langevin
30

---

Published in Transactions on Machine Learning Research (4/2026)
Sequential
Picard
Quasi-Newton
K 2
Number of Fixed-Point Iters
10a
10
104
103
103
103
10a
10a
102
10
Wallclock (s)
10a
10
10-1
103
10a
10a
10
10
10*
Sequence Length (T), Experiment is Langevin Dynamics on K-well potential
Figure 9: In similar setting as Figure 4, we instead consider evaluating Langevin dynamics on a potential
defined by the negative log probability of mixture of K anisotropic Gaussians. Here keep D 128 throughout. We observe qualitatively similar behavior to that shown in Figure 4, where both Picard and quasi-Newton iterations enjoy similar convergence and wall-clock speed.
diffusion two-well potentials. However, this time we use three different discretization step-sizes e: we use x 1x10-5,11 1 10-1, and 10-3 We showed that for step size 1 10-5 (the smallest considered, Jacobian most diagonal) that Newton and Picard both converged in a small number of fixed-point iterations, with Picard therefore running faster overall. In contrast, for e 1 x 10-3 (the largest considered, Jacobian closer to dense), Newton continued to converge in small number of iterations whereas Picard required considcrably more iterations to converge. In this setting, Newton remained faster than sequential evaluation on wallclock time, whereas Picard did not.
E.5 Implementation Considerations
In this section, we note that while the presented fixed-point methods are parallelizable, their real-world efficiency depends on the compute environment. Choosing the appropriate method requires balancing con- vergence speed, compulational intensily, and resource availability.
Parallel Associative Scan is Hardware-Sensitive As we have shown, the seiet-point-oo methods dis- cussed can be cast L.DSs and therefore be parallelizod over the sequence length using a parallel associative scan. However, the practical performance of these operations depends strongly on the hardware and low-level implementation details.
For example, modern GPUs (e.g., NVIDIA A100, H100) are highly optimized for tensor operations and large batch matrix multiplications, which favor methods like Newton and quasi-Newton iterations that perform fewer, more intensive steps.
The performance gains from using more iterations of lighter-weight updates (e.g., Picard iterations) are hardware-dependent. For example, Sarnthein (2025) showed that writing custom kernels for linear recur- rences tensor operations yield near-optimal memory-bound performance.
31

---

Published in Transactions on Machine Learning Research (4/2026)
<= 10-3
10
0-22
2-4
D-8
12-22
Dm4
D-4
22
Numnber Frens
Thtee 60
2270
10
Sequential
Picard
Quasi-Newter
13
Newton
Tixe 61
107 107 107
107
Figure 10: Varying the step size In this variant of discretized Langevin diffusion two-well potential, vary the step size For the smallest step size. 1 10-5 the Jacobian is most diagonal dominant of these three settings, and Picard iterations are faster on a wallclock time than Newton iterations. For the largest step size e = x 10-3 the Jacobian is least diagonal dominant of these three settings, and Picard iterations are slower on wallclock time than Newton iterations. We do not show Jacobi iterations because they already were demonstrated struggle discretizations of Langevin diffusions in Figure 4.
32

---

Published in Transactions on Machine Learning Research (4/2026)
Memory Usage and Tradcoffs stantially. In particular:
The memory requirements of different fixed-point iterations vary sub-
Newton iterations require storing and manipulating full Jacobians, which scales O(D2T) in space. This can become prohibitive for long sequences high-dimensional hidden states.
Quasi-Newton iterations reduce memory cost by using diagonal Jacobians, bringing the com- plexity down to O(DT). This is often a sweet spot for balancing memory usage and convergence rate.
Picard and Jacobi iterations are the most memory-effccient, requiring only the storage of current and previous state estimates (O(DT)), and no Jacobian-related storage.
In practice, when parallelizing over long sequences (T' > D), the memory cost is often dominated by the size of intermediate state reprosentations and the need to unroll computations over multiple fixed-point iterations. Chunking (dividing the sequence into smaller windows) and truncation (limiting the number of fixed-point iterations) are useful strategies to reduce memory usage in these settings (Dao et al., 2022; Shih et al., 2023; Selvam et al., 2024; Geiping et al., 2025; Zoltowski et al., 2025)
Numerical Stability For all fixed-point methods, numerical stability is a concern (Yaghoobi et al., 2025). In particular, LDS matrices with spectral norm close to or greater than one can cause numerical instabilitics in the parallel scan operation (Gonzalez el. al., 2024; 2025). This is especially critical high-precision tasks or long sequences, and practitioners should monitor for numerical divergence or the accumulation of floating-point error.
F
Extended related work
F.1 Discussion of parallel chord methods
Ortega & Rheinboldt (1970) discuss at length ilerative methods for solving arbitrary systems of nonlinear equations F(x) 0 using iterations of the form
x (++1) = -j((())) (F)(())
(21)
for some matrix (x). In general, can be a function of the current iterate x(( a fixed and constant matrix. Newton's method corresponds to
((x OF
8x ())
When J is fixed and constant, Ortega & Rheinboldt (1970) describe the resulting family of fixed-point iterations as parallel-chord methods. However, we will use this term for all iterative methods with updates of the form in eq. (21), which includes both Newton and Picard itcrations.
The term "parallel" in this context does not have to do with applying a parallel scan the sequence length (which we discuss at length in this paper). Instead, "parallel" in "parallel-chord methods" refers to the way in which Newton's method finds the zero of a function by making a guess for the zero, and then forming a chord that is parallel to the function at the current guess (see Figure 11). In one-dimension, the lincarization is a line (a chord), while in higher-dimensions the lincarization is in general a hyperplane. In Newton's is method, the chord/ chord/hyperplan /hyperplaee tangent to the function at the current guess, while for other parallel-chord methods the approximate linearization is in general not tangent.
The equation F(x) is fully general way to represent system of nonlincar equations. However, in this paper, focus parallelizing Markovian state space models. Thus, consider the special setting where
F(x)= r(xi:) [m] f(00), x2 {(x1), f(x2), ,T ----------
(22)
33

---

Published in Transactions on Machine Learning Research (4/2026)
Newton's Method
[())
(00
C
Figure 11:
The term "parallel" in parallel-chord methods. Here we illustrate iterations of Newton's
method for root-finding on the one-dimensional cubic function ((x) = (x - 0.4) + 0.45(x 0.4). We observe that each iteration of Newton's method involves forming a "parallel chord" to the function (shown in color).
where f the nonlinear transition function, as defined in eq. (1), and the residual function defined in eq. (10). Thus, in the context of parallelizing sequences, it follows that
ID
0
f2
(x1)
ID
0
de
OF
0 (I2)
ID
0
0
(x)
8x
(23)
0x
0
TD
0
0
0
ofr
(IT-1)
Ipl
dr
When plug this form of the Jacobian into J in eq. (21) and simplify, we obtain the linear dynamical system in eq. (5), ie. Newton's method.
In their treatment of Picard iterations, Ortega & Rheinboldt (1970) consider a more general formulation than that presented Shih et al. (2023) or in eq. (8). Instead, similar to the definition presented Appendix C.2.3 of Gu al. (2021), Ortega & Rheinboldt (1970) define Picard iterations in the setting where we have removed linear component of F, namely we have written
F(x) = Jx 1 G((),
(24)
for some constant, nonsingular matrix and nonlinear function G(-). Note that such redefinition of F(.) in terms of and G(.) is always possible and not uniquely determined. After making such a redefinition, Ortega & Rheinboldt (1970) define a Picard iteration as an update of the form
(1+1) = -j-(G(x0))
(25)
However, by multiplying both sides of eq. (24) by j it follows that
J j- = x (i) - --(FF(x()),
showing that the Picard iterations as defined in eq. (25) fit into the parallel-chord framework set eq. (21). Note that Picard iterations defined by Shih et al. (2023) in eq. (8) of this paper also fit into
34

---

Published in Transactions on Machine Learning Research (4/2026)
the framework of eq. (24): in the context of evaluating discretized ODEs, the residual defined in eq. (22) becomes
(x+1(x) = I++1 Iz -g((
Thus, in the context of eq. (24); we have that the resulting Gx(x) = (9t--((t-11), while the resulting J operator is given by eq. (20).
When we plug this J into eq. (21) and simplify, we obtain the linear dynamical system in the "Picard" row of Table 1. In general, the fixed-point methods of the common form given by eq. (4) all give rise to J RTDXTD matrices of the form shown in eq. (19).
Thus, Ortega & Rheinboldt (1970) unites Newton and Picard iterations for the general root finding problem F(x) 0 0 under the umbrella of parallel-chord methods, which are iterative updates of the form of eq. (21). The framework we provide in Table 1 can be understood as a specialization of parallel-chord methods for the particular problem of sequential evaluation discussed in eq. (1). Nonetheless, we focus on how in the specific problem of sequential evaluation, which is of great interest in many arcas of machine learning, a wide variely of fixed-point methods become iterative application of LDSs, allowing them to be parallelized over the sequence length with an associative scan. This important perspective about parallelizability, which is of great interest in machine learning, is not discussed in Ortega & Rheinboldt (1970) because they are considering more general problem.
Ortega & Rheinboldt (1970) also discuss in their Chapters 7 and 10 how the closeness the "parallel chord" (in general and higher dimensions, the "approximating hyperplane") to the true linearization of the function (Newton's method) affects the number of iterations needed for the parallel-chord method converge. This analysis is directly analogous to our study of the effect of [](x):)) J J(x1r)
|| the rate of convergence of fixed-point methods, see Lemma 1 In particular, in Chapter 10 of Ortega & Rheinboldt (1970), they consider the rates of convergence of fixed-point methods with updates taking the form of
x (i+1) = =x(x
(26)
for some function U(... Ortega & Rheinboldt (1970) use the name one-step stationary methods such fixed- point methods with updates with the form eq. (26).
For parallel-chord methods of the form given in eq. (21), it follows that
U(x()) = x - -
(27)
In particular, in their Chapters 7 and 10, Ortega & Rheinboldt (1970) introduce and study (U, F, x,). which determines the rate of convergence of iterative methods with updates of the form given by eq. (26) the solution x* of F(x)= = 0. Thoy define as
a ,U,,F = ( OU
(xxx(xx))) (xx)
(28)
where p(M) denotes the spectral radius of matrix M. In the context of parallel-chord methods where U(.) is given by eq. (27), it follows that
OU 0x
OF
= J(xx)= --x(x)
Ox (x"),
because F(x*) Thus follows that if J /x/ax(x*), then o Thus, lower values of indicate that J is a good approximation of the Jacobian matrix 07/8x evaluated at the zero x* of F, while higher values of a indicate that J is poor approximation for 8F/2x. Ortega & Rheinboldt (1970) then use a in their Chapter 10 (in particular, thoir Thcorem 10.1.4) to prove linear rates of convergence 10 for one-step stationary methods within a neighborhood of the solution x"
10 where the rate is given by
35

---

Published in Transactions on Machine Learning Research (4/2026)
Thus, takeaway from Ortega & Rheinboldt (1970) (as paraphrased from Gasilov et al. (1981)) is that the closer J is to 0a/2x, the fewer iterations are needed for convergence to x". This takeaway is extremely similar to our guidance, though we specialize to the partieular system of equations given by eq. (3) that results from the goal of rolling out the Markov process given by eq. (1).
However, in the setting we consider in this paper paper-using using fixed-point iterations of the form eq. (4) to solve nonlinear equations of the form eq. (3) Theorem 10.1.4 of Ortega & Rheinboldt (1970) is actually trivial. By "trivial," " we mean that it does not distinguish between the convergence rates of any of the exed-point iterations we focus on in this paper.
To make this point more precisely, we review the notion of rool-convergence, more commonly known R- R-convergence.
Definition 3 (R-convergencee) Let be a fized-point operator with fized-point x'. Let C(A,x*) be the set of all sequences generated by A which converge x*. Then the Ri-factors of A at x* are given by
R1(A,x*) = sup lim |xx(9 x*1 b20 E 0(Axx))
(29)
1+000
Intuitively, Rí (AA, x') gives the rate of linear convergence of fixed-point. operator A to its fixed-point x Theorem 10.1.4 of Ortega & Rheinboldt (1970) implies that if A is a one-step stationary method with update given by U(.,, then R1 (A4,x)) = (U, (U,,F,**), Therefore, if a > 0, then o is the rate of R-linear convergence of A to x*, while if a = 0 we say that A converges R-superlinearly. However, it is important to note that these definitions are asymptotic in nature.
The fixed-point iterations considered in this paper, i.c. following the common form eq. (4), all have a 0. and therefore can be said to converge R-superlinearly.
Proposition Let F(x)= = be a nonlinear equation of the form eq. (3) with solution x*. Let be parallel-chord method with stred-point>2 x*. Then
a (U,F,x') =
Proof. Both ox))x(x*) and J(x*) are lower-triangular matrices with all D x D identity matrices on their main biagoadia In particular, J-1 is also a lower-triangular matrix with all D x D identity matrices on its main block-diagonal. Consequently, the product J- j-1 OF is also lower-triangular matrix with all D x D Õx
identity matrices its main block-diagonal. As result, I- j- OF zeros on its main block-diagonal, and so has all its cigenvalues equal to 0. Consequently, its spectral radius is equual to zero.
is a lower-triangular matrix with all
dx
It may scem counterintuitive that even Jacobi iterations technically enjoy R-superlincar convergence in the context of parallelizing Markov processes. However, this seemingly strange result stems from the asymptotic nature of Definition 3 of R-convergence, and the fact that Proposition 1 of Gonzalez et al. (2024) guarantees that all fixed-point iterations of the form given by eq. (4) will converge to x" in a finite number of iterations
(T, to be exact). Therefore, for any LDS fixed-point scheme, we always have lim [ x1 =
However, in both Proposition 4 of Lu et al. (2025) and Proposition 3 of this paper, we effectively get around this difficulty by considering the spectral norm instead of the spectral radius. The spectral norm always bounds the spectral radius, and so by focusing on spectral radius, Ortega & Rheinboldt (1970) could get tighter bounds (faster rates of convergence). However, in our setting the spectral radius cannot distinguish between any of the fixed-point methods, and so we instead use the looser bound provided by the spectral norm, which can distinguish between the different fixed- point methods. Note that the core entities are effectively the same, as y defined in eq. (15) is equal to /xx///x))ll>
We fallow the presentation of Chapter 9 of Ortega & Rheinboldt (1970), particular Definition 9.2.1.
36

---

Published in Transactions on Machine Learning Research (4/2026)
Finally, again, because all of our fixed-point methods converge in most T iterations, asymptotic notions of linear convergence are not suitable to fully capture the behavior of these fixed point methods. For this reason, we use empirical case studies in Section 4 to show the efficacy of the intuition, inspired by Proposition 3, that the closeness of A, to 0f/0z impacts the number of iterations needed for A to converge. This empirical approach also highlights how the increased compulational cost. of higher-order fixed-point methods affects wall-clock time on GPUs.
F.2 Trust-region fixed-point methods
Many other fixed-point techniques can fit into the general framework we propose in Table and eq. (4) (see the general algorithm in Algorithm 1). We focus on Newton, quasi-Newton, Picard, and Jacobi in the main text because of their prominence and canonicity. However, other fixed-point iterations also fit into this framework.
For example, Gonzalez al. (2024) introduces the scale-ELK algorithm which for some k E [0, 1] sets the transition matrix A, to be
A, (1 k) ff++
Ox
Scale-ELK can also be applied to the diagonal approximation for quasi-Newton methods. However, scale- ELK introduces an additional hyperparameter
Therefore, we propose clip-ELK, which is a hyperparameter--F free approach to achieve the same goal of a stable LDS. Clip-ELK applies to the diagonal approximation only, and simply clips each element of Å, (which in this sctting is diagonal matrix) to be between [-1, Clip-ELK is immediately also part of the framework set out in Table of fixed-point methods that parallelize the evaluation of sequences via iterative application of LDS; moreover, it too must also converge globally by Proposition 1 of Gonzalez et al. (2024). An interesting and important direction for future work includes both developing and interpreting more lixed-point methods the context of this framework that we propose.
37


Reviews:
[{"review_id": "tQFHGAhG7r", "reviewer": "Reviewer_N66A", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "## Contributions: \n\nThis paper proposes a new framework to unify four different algorithms for solving fixed-point iteration in the language of linear dynamical systems. The resulting LDS is clean, showing that the four different algorithms—Jacobi, Picard, Quasi-Newton, and Newton—differ only in the transition matrix in the iteration rule. The resulting practical guidance, \"Use the simplest approximate Jacobian as possible, but no simpler,\" is clean and precisely captures the essence: to use as much full Jacobian information as possible at minimal cost.\n\n## Strengths: \n\nThe strengths are clear; the proposed framework and summarization of ideas are clean and indeed provide a novel, unified view of these algorithms.\n\n## Weaknesses: \n\nDespite generally enjoying this work, I feel its depth is insufficient in both algorithmic and practical aspects.\n\nFirstly, from an algorithmic aspect, the current version only discusses overly-simplified algorithms (Jacobi, Picard) or overly-expensive algorithms (diagonalized Newton and Newton). These algorithms, most of the time, are not optimal in practical scenarios. A majority of scenarios in machine learning involve a more general form of quasi-Newton method, where the transition matrix \nis estimated in an economical way; for example, Anderson acceleration considers using history residuals to estimate the full Jacobian. Therefore, I think the current version misses a discussion on more general forms of quasi-Newton methods, which are more widely used and are more practical choices.\n\nSecondly, from an experimental aspect, the depth is also not enough to support the claim \"Use the simplest approximate Jacobian as possible, but no simpler.\" The current experiments use three different settings, where they have different Jacobian structures specific to their problem setting. However, this is not very instructive guidance, because most of the time a practitioner cannot know the Jacobian structure very well before running these algorithms, especially when the full Jacobian information is extremely hard to evaluate.", "claims_and_evidence_yn": "Yes", "claims_explanation": "The main claim is that the authors propose a framework that unifies several classical fixed-point iteration algorithms, which is evident in the summarization in Table 1.\n\nThe experimental section (Section 3) also validates how the problem structure favors different algorithms.", "audience_yn": "Yes", "audience_explanation": "Yes, because fixed-point iteration is a powerful and ubiquitous modeling tool in the machine learning community.", "requested_changes": "Overall, I have three suggestions for the authors:\n\nIn Table 1 and the discussion, I suggest the authors include more practical forms of the quasi-Newton method, where the transition matrix is estimated rather than being exactly computed. This will make the framework cover more widely used and practical algorithms.\n\nIf possible, I suggest the authors discuss ways for practitioners to probe the black-boxed Jacobian structure before deciding on an algorithm. For example, via finite differences or other economical ways.\n\nThirdly, I suggest the authors use a more controlled experiment to rigorously demonstrate the differences between the algorithms. For example, in the Langevin diffusion section, you could set epsilon to various values to toggle the Jacobian matrix from more diagonalized to more dense, and then validate how the performance of different algorithms varies among these settings. This would provide a more precise view of how the Jacobian structure affects the algorithmic choice.", "broader_impact_concerns": "This paper focuses on a unified algorithmic framework for fixed-point iteration, and I did not see any ethical issues.", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "BDOFxcQpOD", "reviewer": "Reviewer_LPFg", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "No", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "This paper studies the topic of solving the fixed-point problem of sequential processes in parallel. In particular, the paper investigates a unified underlying framework for Newton, Quasi-Newton, Picard, and Jacobi iterations. By contrasting the difference of the algorithms above under the same framework, the paper identifies conditions where each algorithm is most suitable for. The paper accompanied the observation with examples in the real-world applications, and supported the claim with experimental results.", "claims_and_evidence_yn": "No", "claims_explanation": "First, the paper claims that the iterates of Newton, Quasi Newton, Picard, and Jacobi iterations for nonlinear recursions can be modeled as linear dynamical systems. However, in the formulation of Eq. (4), the iterate-wise update still depends on a nonlinear function $f_{t+1}(x_t^{(i)})$. It is not clear why Eq. (4) falls into the category of a linear dynamical system or why it can be solved using parallel scan algorithms. Moreover, the paper claims that how the structure of the problems informs which method one should use. The definition of the structure of the problem is vague here, as in the examples of Section 3, the paper simply discussed the dependency on the property of the per-iteration Jacobian. I believe this is a dependency more on the iterates rather than the general structure of the problem itself.", "audience_yn": "Yes", "audience_explanation": "The general problem of sequential process indeed appears in various machine learning settings, including unrolled networks, RNNs, and diffusion models. Speeding up this process can easy the computation during the inference of these models. The paper discusses the situations where different existing algorithms are most suitable for, which can be a guideline for practitioner to choose the correct algorithm.", "requested_changes": "1. I suggest the paper to include a discussion of the parallel scan algorithm, and how it applies to the particular case of solving linear dynamical systems.\n\n\n2. The paper should clarify the assumption of $f_{t+1}(x_t^{(i)})$ in Eq. (4), and, if it is a nonlinear function, explain why the resuling updates falls into the category of a linear dynamical systems.\n\n3. The paper should be more specific on the \"structure\" of the problem instance, and explain that the choice of the algorithms depends on the per-iterate Jacobian matrix.", "broader_impact_concerns": "None.", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "8aa5g2jGZU", "reviewer": "Reviewer_Hz5N", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "The paper discusses techniques to parallelize sequential models. The authors propose a general framework based on linear dynamical systems (LDS) that unifies several recent fixed-point methods that parallelize sequential models, including Newton, Picard, and Jacobi iterations. A parallel scan algorithm is then used to solve the dynamical system in log time. Case studies on solving the group word problem, RNN’s, and Langevin diffusions are given to highlight the differences between the three iteration methods.\n\nStrengths:\n\n1. The proposed LDS framework is elegant and effectively captures the three discussed iteration methods. Clear derivations of the correspondences between the proposed framework and the three methods are also given.\n\n2. The paper verifies the methods empirically and conducted ablation studies between different methods in different experiment settings. Through the experiments, the paper also gives an empirical rule for choosing an iteration method.\n\nWeaknesses or Suggestions:\n\n1. The contributions are limited. The Newton, Picard, and Jacobi iterations are all previous works, and turning these methods into the form of $x_{t+1}^{(i+1)} = f_{t+1}\\left(x_t^{(i)}\\right) + A_{t+1}\\left(x_t^{(i+1)}-x_t^{(i)}\\right)$ is rather straightforward.\n\n2. The main heuristic for choosing between different iteration methods is based on the Jacobian $\\frac{\\partial f_{t+1}}{\\partial x_t}$, where the Newton iteration uses the full Jacobian for the $A$ matrix. While the proposed framework uses a general $A$ matrix, the discussion in the paper seems to be still centered around the Newton method. Is it possible to discuss a whole new family of iteration methods based on the general framework and analyze the new family both theoretically and empirically?\n\n3. Is it possible to give an analytical bound for the number of fixed point iterations based on some metric of distance between $A_{t+1}$ and the Jacobian $\\frac{\\partial f_{t+1}}{\\partial x_t}$?\n\n4. The experiments are mostly synthetic. It would be stronger to show speed-up results on trained real tasks, e.g. a real diffusion sampling with SDE or ODE samplers.\n\n5. Although an introduction of the parallel scan algorithm is given in the appendix, it would be better to outline the exact parallel scan algorithm to solve the LDS’s mentioned in the main paper.", "claims_and_evidence_yn": "Yes", "claims_explanation": "The claims made are generally supported by evidence. However, as mentioned above, since the experiments are mostly synthetic, it would be more convincing to include more real-world experiments with larger scales.", "audience_yn": "Yes", "audience_explanation": "Since sequential sampling is prevail in LLMs, developing efficient parallel sampling techniques is definitely valuable and of broad interest to the community.", "requested_changes": "Please see my suggestions in the Summary Of Contributions section.", "broader_impact_concerns": "N/A", "additional_comments": "", "extra_scores": {}, "extra_text": {}}]

Author rebuttal:

Decision: Accept (poster)