Detecting Invariant Manifolds in ReLU-Based RNNs
Authors: Lukas Eisenmann
Alena Brändle
Zahra Monfared
Daniel Durstewitz
Venue: ICLR 2026 Poster (2026)
Keywords: Dynamical Systems
Recurrent Neural Networks
Attractors
Nonlinear Dynamics
Analysis
Primary area: learning on time series and dynamical systems
TL;DR: A novel semi-analytical algorithm for detecting the stable and unstable manifolds in RNN state spaces, which separate basins of attraction and lead to chaos if they intersect.

Abstract: Recurrent Neural Networks (RNNs) have found widespread applications in machine learning for time series prediction and dynamical systems reconstruction, and experienced a recent renaissance with improved training algorithms and architectural designs. Understanding why and how trained RNNs produce their behavior is important for scientific and medical applications, and explainable AI more generally. An RNN's dynamical repertoire depends on the topological and geometrical properties of its state space. Stable and unstable manifolds of periodic points play a particularly important role: They dissect a dynamical system's state space into different basins of attraction, and their intersections lead to chaotic dynamics with fractal geometry. Here we introduce a novel algorithm for detecting these manifolds, with a focus on piecewise-linear RNNs (PLRNNs) employing rectified linear units (ReLUs) as their activation function. We demonstrate how the algorithm can be used to trace the boundaries between different basins of attraction, and hence to characterize multistability, a computationally important property. We further show its utility in finding so-called homoclinic points, the intersections between stable and unstable manifolds, and thus establish the existence of chaos in PLRNNs. Finally we show for an empirical example, electrophysiological recordings from a cortical neuron, how insights into the underlying dynamics could be gained through our method.

Published as a conference paper at ICLR 2026
DETECTING INVARIANT MANIFOLDS IN RELU-BASED RNNS
Lukas Eisenmann '2"', Alena Brändle 1.2" Zahra Monfared and Daniel Durstewitz
www.ddaannni...diii.t.... @zi-mannheim.de
Dept. of Theoretical Neuroscience, Central Institute of Mental Health, Mannheim, Germany Faculty of Physics and Astronomy, Heidelberg University, Heidelberg, Germany
è IInterdisciplinaryy Center for Scientific Computing, Heidelberg University, Germany
*Faculty of Mathematics and Computer Science, Heidelberg University, Heidelherg, Germany
"These authors contributed equally
ABSTRACT
Recurrent Neural Networks (RNNs) have found widespread applications in ma- chine learning for time series prediction and dynamical systems reconstruction, and experienced a recent renaissance with improved training algorithms and ar- chitectural designs. Understanding why and how trained RNNs produce their behavior is important for scientific and medical applications, and explainable Al more generally. An RNN's dynamical repertoire depends on the topological and geometrical properties of its state space. Stable and unstable manifolds of periodic points play a particularly important role: They dissect dynamical system's 's state space into different basins of attraction, and their intersections lead to chaotic dynamics with fractal geometry. Here introduce novel algorithm for detecting these manifolds, with a focus on piecewise-lineas RNNs (PLRNNs) employing rectified linear units (ReLUs) as their activation function. We demonstrate how the algorithm can be used to trace the boundaries between different basins of attraction, and hence to characterize multistability, a computationally important property. We further show its utility in finding so-called homoclinic points, the intersections between stable and unstable manifolds, and thus establish the existence of chaos PLRNNs. Finally we show for an empirical example, electrophysiologica record- ings from a cortical neuron, how insights into the underlying dynamics could be gained through our method.
INTRODUCTION
Recurrent neural networks (RNNs) are widely employed for time series forecasting (Park et al. 2018 Gu and Dao 2023) and dynamical systems (DS) reconstruction (Hess et al. 2023) Brenner al. 2022] 2024b and Ou (023) or neuroscience (Durstewitz al.] 2023), as well as in medical domains. RNNs experienced a recent revival due to the advance of new powerful training algorithms that avoid vanishing or exploding gradients (Hess et al.
Platt et al. 2023), especially in scientific applications like climate modeling (Patel
2023). and novel network architectures (Schmidt et al. 2021; Rusch and Mishra 2020 Peng et al.
2023 Gu and Dao. 2023 Rusch et al.), in particular the context of state space models (which are essentially linear RNNs with nonlinear readouts and input gating: (Gu and Dao 2023 Orvieto et al. (023)). ). What lags behind, like in many other areas of deep learning, a thorough theoretical understanding of the behavior of these systems and how they achieve the tasks they were trained on. Yet, such an understanding is crucial especially in scientific and medical areas where we are often interested in using trained RNNs as surrogate models for the underlying DS, providing mechanistic insight into the dynamical processes that generated the observed time series.
Formally, RNNs recursive maps, hence discrete-time dynamical systems (Mikhaeil et al. (022). Dynamical systems theory (DST) offers a rich repertoire of mathematical tools for analyzing the behavior of such systems (Guckenheimer and Holmes| 2013). Yet, exploration of high-dimensional

---

Published as a conference paper at ICLR 2026
RNN dynamics remains challenging due to limitations in current numerical methods, which struggle to scale and often only yield approximate results (Katz and Reggia 2018 Golub and Sussillo 2018).
The dynamical behavior of system is governed by its state space topology and geometry (Hasselblatt and Katok, 2002). most prominently topological objects like attractors, such as stable fixed points, cycles, chaotic sets, which determine the system's long-term behavior. Similarly important are the stable and unstable manifolds of fixed points and cycles, although they received much less attention in the scientific and ML communities. Stable manifolds are the sets of points in a system's 's state space that converge towards an equilibrium or periodic orbit in forward-time, while unstable manifolds, conversely, are the sets of points that converge points delineate the boundaries between basins of attraction in multistable systems, that is, systems that harbor multiple attractor objects between which it can be driven back and forth by perturbations like external inputs or noise (Feudel et al., 2018). Multistability has been hypothesized to be an important computational property of RNNs, most prominently in computational neuroscience where it has been linked to working memory (Durstewitz et al., 2000 or decision making (Wang) 2008).
in backward timee Stable manifolds of saddle
Tracing out stable and unstable manifolds is also important for finding homo- and heteroclinic orbits, which connect a cyclic point to itself or to another such point, respectively. These orbits provide skeleton for the dynamics, forming structures like separatrix cycles (Perko (00T) and heteroclinic channels (Rabinovich et al. 2008) which have been implied in flexible sequence generation. Intersec- tions between stable and unstable manifolds further give rise to so-called homoclinic points which create sensitive regions in a system's state space associated with chaos (Wiggins) 1988). Chaotic behavior, in turn, or regimes at the edge-of-chaos, have been associated with increased expressivity (larger function classes that can be emulated) and computational power in RNNs (Siegelmann and Sontag 1992 Bertschinger and Natschläger 2004 Pereira-Obilinovy et al. 2022
2023
(Mikhaeil et al.
ess et al. 2023).
Here we introduce a novel algorithm for computing the stable and unstable manifolds of fixed and cyclic points for the class of piecewise-lineat (PL) RNNs (PLRNNs), which use PL nonlinearities like the ReLU as their activation function. Unlike traditional techniques designed for smooth dynamical systems- -ucch as numerical continuation methods (Krauskopf et al. 2007) -the proposed method explicitly exploits the piecewise- linear structure of PLRNNs to enable the exact location of stable and unstable manifolds. In contrast to ODE systems, methods for discrete-time systems are much more scarce, especially if these involve discontinuities in their Jacobians like ReLU- RNNs. PLRNNs have emerged as one of the most powerful classes of models for DS reconstruction (Hess)
et al.
2023 Brenner al. 2022 2024b Nassar et al.) 2019), partly because the linear subspaces of such models support the indefinite online retention of memory contents without exploding or vanishing gradients (Schmidt et al. [021 Orvieto et al. 2023 Gu and Dao 2023), and partly because their PL structure makes them more tractable (Coombes et al.) 2024). In fact, PL systems have been popular in engineering (Bemporad et al. 2000), Carmona et al.) 2002) Storace and De Feo 2004)
and mathematics (Alligood et al.) for decades for these reasons, and variety PL models similar in structure to PLRNNS exist, like threshold-linear networks (TLNs; Curto and Morrison (2023); Parmelee et al.) (0022); Morrison et al.
et al. 2019 Coombes et al.] Siipson 2023)
(2016): Hahnloser and Seung (2000)), switching linear DS (SLDS: Linderman et al.
(2016) 2017)), 2016). with
or Lur' systems in control theory (Su et al., 22023] Waitman et al. Biil et al.
substantial body of mathematical theory built around them (Parmelce et al. 2022 Morrison et al.
2016 Brenner al. [202aa] Avrutin et al. 2019 Amaral et al. |006 SSimpon (2023 Coombes et al.
2024 Casdagli]] 1989; Ives and Dakos [012 Costa et al.) 2019). Further, efficient algorithms for exactly localizing fixed and cyclic points in polynomial time exist (Eisenmann et al.) 2023), which we will build here. We illustrate how our algorithm can be employed to delineate the boundaries of basins of attraction and to detect homoclinic intersections, thus establishing the existence of chaos, and show how it can be used on empirical data to gain insight into dynamical mechanisms.
2 RELATED WORK
While there is an extensive literature by now on DS reconstruction with RNNs (renner et al. 2022 2024a Hess al. 2023 Patel and Oti 2023 Platt et al. 2023; Cestmik and Abel 2019). work on algorithms for dissecting topological properties of trained systems is much more scarce. Existing research is almost exclusively focused on algorithms for locating stable and unstable fixed points
2

---

Published as a conference paper at ICLR 2026
and cycles (Golub and Sussillo, 2018: Eisenmann et al. 2023), but not the invariant manifolds associated with them, despite their crucial role in structuring the state space and dynamics. More generally in the DST literature, related class of methods are so-called continuation methods which numerically trace back important curves and manifolds in dynamical systems, mostly for systems of ordinary differential equations (Krauskopf et al., 2007; Osinga 2014). Recently methods have been developed specifically for PL maps (Simpson (023). but, similar to continuation methods, they seriously suffer from the curse of dimensionality. Thus, current methods effectively work only for very cow-dimensional (5 5d) systems, in contrast to the size of many modern RNNs used for time series forecasting or dynamical systems reconstruction. The specific structure of PL maps, of which all types of PLRNNs are specific examples, considerably eases the computation of certain topological properties. This also provided the basis for the SCYFI algorithm for locating fixed points and cycles, which scales pol polynomially, in fact often linearly, with the dimensionality of the PLRNN's latent space (Eisenmann et al., 2023). To our knowledge, no existing method efficiently detects stable and unstable manifolds in discrete- timeeet-tmeeeee RNNs. Our approach fills this gap by leveraging the PL nature of ReLUUbased RNNs to construct manifolds directly, bypassing the limitations of traditional techniques.
3 METHODS
3.1 PLRNN ARCHITECTURES
PLRNN (Durstewitz 2017 is simply a ReLU-based RNN
2e Az-1 + W W + h.
where A E RMXM is a diagonal matrix, W RMxxM an off-diagonal matrix, the element-wise ReLU non- con-linearity ((-) max(0. -), and h E IR M bias term. Several variants of the PLRNN have been introduced to enhance its expressivity or reduce its dimensionality (Brenner al. 2022). of which we picked for demonstration here the shallow PLRNN (shPLRNN: (Hess et al. 2023))
A A + W₁ Φ 22 + 122) + h, with Wi E R M W₂ E R and the recently proposed almost-linear RNN (ALRNN; (Brenner et al. 2024a)). which attempts to use as few nonlincarities as needed for the problem at hand, (zz)) = 1,,
P,r,, max 2 -MM)) max (0, #M.,))]
The PL structure of these PLRNNS can be exposed by rewriting the ReLU in Equation (1 as
2: = 1 = (A + W 1)-1) + h,
Dj == diag(d;) (d)) with d. = -(d,,.
dM.e) and di,t = 0 if zi,e S O and di,, otherwise Monfared and Durstewitz. (0200) For the ALRNN, we have dı P.t = 1 Hence, we have 2
M linear subregions for the standard PLRNN, S En (#) for the shPLRNN (Pals et al. 2024 and 2f for the ALRNN.
3.2 MATHEMATICAL PRELIMINARIES
Recall that fixed point of recursive map F(xx-1) is point xx for which xx = F(x*), and a cyclic point with period m is a fixed point of the m-times iterated map, i.c. such that at+mm = 1. (xxx) = 2, (hence, a fixed point is period-1 point). An m-cycle is a periodic sequence {x) ...,,,,, of such points with all points distinct, x; A x V1 s j S m. Definition (Un-/stable manifold). Let F: M + RM be a map and p be a hyperbolic period-m cyclic point of F. The local stable manifold of p, w (p). is defined as
W-+1) (p) = fx E RM Funm (x) + p as n + 00
The local unstable manifold of p. w (p), is defined
w-(") (p) == {a E TR F-"m(a) +
n x
The global stable manifold is the union of all preimages of the local stable manifold, and the global unstable manifold is created by the union of all images (forward iterations) of the local unstable manifold (Patra and Banerjee 2018):
p)(+)((p) = U F (=((((t))
Loc (//tt!)
(p)).
w(-11(p) = U Fm loc
(w))
(3)
n=l
3

---

Published as a conference paper at ICLR 2026
If the map F' is noninvertible (i.e., does not have unique inverse), the (global) stable manifold could be disconnected, making its computation hard as we need to trace back disconnected sets of points in time to determine it. However, most commonly continuous-time DS å = (() in DS reconstruction through our RNN map Fe. For such systems, the Picard-Lindelöf theorem guarantees uniqueness of solutions and the flow (solution) operator étt, zn) 25 + So f()) du will be diffeomorphism (Perko 2001). hence invertible, such that is reasonable to assume (or enforce) invertibility for Pa as well.
aim to approximate smooth
For PLRNN, since it is PL map, in each linear subregion we have the following explicit (non- recursive) expression for the dynamics along trajectories (z+ W 2t-11 +h, W A+ W (Zz-1))) 1)):
Zt Σ CoAfus + X r=l
1)! wr
Σw*h.
(4)
k=
where the o; E are determined from the initial condition, v, are the eigenvectors, w; are general- ized eigenvectors (associated with defective eigenvalues with geometric multiplicity less than the algebraic). This formula describes the exponential evolution for all non-defective eigenvalues X, and polynomially modulated exponentials for the defective eigenvalue X Appx. B for details (Perko 2001, Hirsch et al. 2013). While this expression is for single orbits, it descrihes how curvature in the manifolds can arise
the boundaries of linear subregions (ef. sect. 3.3
Stable manifolds of saddle objects segregate the state space into different basins of attraction (Ganatra and Banerjee, 2022), which are sets of points from which the state evolves toward one or the other attractor. Formally, they are given by (Alligood et al. 1996):
Definition 2 (Basin of attraction). Let F' IRM attractor A the largest open set ((A) S IM' (containing A) such that for every point a E B(A), the iterates of under the map F converge to A in the forward limit:
R. be a map. The basin of attraction of an
B(A) = x RMM, d(F*(()))))) (x). + 0 as k + 00)
{
where F denotes the k-times iterated map F.
Stable and unstable manifolds are crucially important for the system dynamics not only because they dissect the state space into different regions of flow, but also because their intersections can give rise to complex types of dynamics like heteroclinic channels or chaos (Rabinovich et al., 2008 Perko 2001}.
Definition 3 (Homoelinic orbit). Let p be saddle fixed point (or saddle cycle) of the map F; R + R A homoclinic orbit is a trajectory Ohom = {xn JneN that connects p to itself, i.e., 2n + p as - =o0. Hence,
Ohom C cW(1)
(p) n w(-))(p)
If x: * p is point where the stable and unstable manifolds of intersect, then x is referred to as homoclinic point or intersection (Patra and Banerjee 2018 Perko [200). Such an intersection of the stable and unstable manifolds leads to horseshoe structure associated with a fractal geometry and chaos (see Appx. I Wiggins (1988)).
Definition 4 (Heteroclinic orbit). Let p and q be two distinct saddle fixed points (or saddle cycles) of the map F. heteroclinic orhis Ohet = { JEEN connects two different fixed points p and q. - +00 and -
that is, 2x, + defined similarly to homoclinic intersections, and, like homoclinic points, inevitably lead to chaos (Wiggins 1988 Patra and Banerjee 2018).
- -c. Heteroclinic intersections (or points) can be
3.3 LOCATING UN-/STABLE MANIFOLDS
For our manifold detection algorithm, we exploit the PL structure of PLRNNs, which allows for an analytical construction of the manifolds within each linear subregion. We further assume that Fg is invertible (see above), such that the global un-/stable manifolds will be connected sets. We start by using a previously proposed algorithm, SCYFI (Eisenmann et al. (2023), to locate saddle (cyclic) points p of interest, for which we would like to construct the un-/stable manifolds. We then proceed as follows (Algo. T see Fig. for illustration of these steps in 1d):
also Pals al. (2024)),
4

---

Published as a conference paper at ICLR 2026
Algorithm 1 Manifold Construction
Hyperparameters: N₁: Number of support points to be sampled
Niter: Maximum number of iterations
Input: p: Periodic point
a E { 1, }}: flag indicating whether stable (+1) or unstable (-1) manifold is to be computed
0: PLRNN parameters
Output: {Qa) J: Support & eigenvectors parameterizing the manifold in each subregion
1: E0 -GEEEGGENVEC(O D(p))
D Compute eigenvectors p b Initialize set of manifold specifications Initialize set of seed & visisted subregions
2: Q t ((p,
3: K K Kiisited t e (GETSCCRECIONNDDG))
4: while n. S Niter A (isempty(Keee) do
5: for E Keeed do
LLoo over regions to expand from Sample N, points on the manifold Propagate forward/backward using Fg Connecting regions not yet visited > Loop over connecting regions A Compute eigenvalues
6:
S + SAMPLEPOINTS(Q2 .N.) S + PROPAGATETTONETTREGONSSS,, 0)
7:
8:
Koon(k) GETSURREGGIONID(S5)) \ K
9:
for j E Kam(()) do
10:
X + GETÉIGENVAL Dý)
11:
if [m())) + V (multiplicity geo (Ag) < multiplicity alg (A1)) then
12:
+ Q"QQu {
>> Determine manifold by kernel-PCA
else
14:
+ Q" U (
b Determine manifold by PCA
end if
16: end for
17:
Kvisited f- Kvisited U coo(())
b Update set of visited regions
18:
end for
19:
Kgeed + -Keon
Make all connecting regions new seed
20:
end while
21:
return {Q*)
D Manifold segments in all subregions
Locally, in the linear subregion harboring p. the unstable (resp. stable) manifold of p is simply given by the affine subspaces spanned by the eigenvectors with absolute eigenvalues within (stable) or outside of (unstable) the unit circle, and hence can easily be computed in closed form from the local Jacobian A + W (p)). Note this remains true for complex eigenvalues (associated with spiral points), as these always come conjugate pairs and hence always span a planar subspace.
2. To extend the un-/stable manifolds across the boundaries into neighboring linear regions, we randomly sample N₂ points on the manifold and propagate them forward by Fg (unstable) or backward by Fø (stable) into neighboring regions. The hyper-parameter N, is chosen such that with high likelihood all neighboring subregions into which the manifold extends are reached, which is mostly fulfilled if we choose N, kiniigh X d, the total number of neighboring subregions times the manifold dimension D For the inversion of Fg, specifically, we need to solve 2,-1 = (A + WD--1 - - h). solution to this eq. needs to be self-consistent, i.e. the signs of the z on the l.h.s. need to be consistent with the entries d;,r-1 in D.-1 on the r.h.s. To address this, we introduce a simple heuristic: 1) Perform a backward step using the current linear subregion D;; 2) perform forward step using the resulting candidate solution 2;-1; 3) if 2; = Fe (z 1), we are done; 4) if not, we take the candidate's 's linear subregion D( to attempt new inversion; 5) if this fails, we start checking the neighboring regions by iteratively flipping bits in D.-1 (see Appx. E. Algorithm 2 for details). The efficiency of this procedure comes from the fact that usually the candidate will lie in the correct linear subregion, and if this is not the case - it usually crosses only into neighboring regions.
3. The manifolds cannot change their dimensionality across subregions, but they may become curved along one or more directions as they cross into neighboring subregion (for instance, if the dynamics in the new subregion is determined by a spiral, see example in Fig. IT). In theory, we only need to compute as many points in each subregion required to uniquely anchor the manifold segment within that region. As we know the type of curvature qua Eq. ideally these
If unsure, one could increase N₂ stepwise until the number of subregions reached platcaus.
5

---

Published as a conference paper at ICLR 2026
would just be d + points for a d < M-dimensional manifold even in the curved case, see Appx. [ In practice, numerical issues (see below) and the type of approximation used may require sampling more points. If the manifold segment is non-curved, we use PCA as a quick and simple means to obtain the affine subspace from d + support points. If the manifold is curved, we use kernel-PCA (see Ap The 2 d eigenvectors with non-zero eigenvalues, together with least one support point, exactly define manifold segment Q7 in subregion k (see Appx. D for another type of approximation for the manifold segment in each subregion).
4.
We then iterate steps 1-3 for at most Nater times, recursively expanding the manifolds across the whole state space. The algorithm's 's final output will be for all K subregions the set of all eigenvectors with at least one support vector positioning each manifold segment.
Note that within each linear region the construction of the manifold is thus analytical from the 2 d support points, while across regions we recursively assemble the whole manifold following Eq. 3 In systems with highly disparate timescales, however, numerical issues may arise sampled support points mainly spread along the fast eigendirections with little variation along the slower eigendirections, leading to poor numerical approximation of a manifold. To balance the contributions from different directions we therefore adjust sampling density inversely to eigenvalue magnitude, with an additional hyperparameter e controlling this weighting (see Appx
Step Initial manifold
Step 2 Propagate to region Step 3
Determine Step
Sumple points 8
manifold
Equation (4)
Deternine
"Repeal
manifold
Process
Figure 1: One-dimensional illustration of the iterative procedure for computing stable manifolds with subregion boundaries indicated in purple-dashed. Step 1: The stable manifold (black) is initialized using the stable eigenvector E(+1) (blue) of the saddle point (green), and sample points (orange) are placed along it. Step These points are propagated until they enter another linear subregion. Step 3:
A new segment of the manifold is determined. Step 4: Repeating this process iteratively reconstructs the full global structure of the stable manifold. For trained model example see Fig. g
3.4 ENFORCING MAP INVERTIBILITY BY REGULARIZATION
In designing our algorithm, we relied on invertibility of the RNN map Fe in the sense that BID,-1 for which (F, (z,, D D--1)) = Zt. This is reasonable assumption as the flow map ((t, for most underlying ODE systems of interest, which we attempt approximate, is invert- ible due to Picard-Lindelöf (Prkko 2001}. However, empirically, invertibility of F' is not always guaranteed (depending on the quality of the approximation). Thus, we enforce this condition by regu- larization during RNN training (for training itself we used the previously proposed shPLRNN (Hess) ct al. 2023) or ALRNN (Brenner et al., 2024a). see sect. B.I , Fe is invertible, if the determinants of the Jacobian matrices of neighboring subregions have the same sign (sign condition) (Fujisawa jet al. 1972). This can be enforced (for any type of Rel. Reaeedd RNN) by adding the tollowing regularization term to the RNN training loss (commonly the MSE loss, see Appx.
Lag A
Σ VESup
max (0, dct(J:)),
(5)
|Siel
computed across small subset Sig of linear subregions, where J. A + WD; is the Jacobian in subregion i, and X regularization parameter. As shown in Fig. A (left), strategic sampling of only 1% of linear subregions, which are traversed by actual trajectories (cf. (Brenner al. 2024a)), is sufficient to ensure almost full invertibility of Fg across the whole state space. At the same time, it hardly affects runtime (Fig, A top-right) and reconstruction performance (Fig, Z bottom). Since invertible flows are an inherent property of many, if not most, DS of scientific interest, we would furthermore expect that this regularization does not hamper, or even improves, the reconstruction of DS from data, in particular if the systems carry intrinsic time reversibility (Huh et al. 2020). This is confirmed in Fig. BBC which shows that with the invertibility regularization in place, training of PLRNN on 10d damped nonlinear oscillator converges significantly faster good solution

---

Published as a conference paper at ICLR 2026
Damped linear Oscillator
Lerena-83
1.0
Time, EEocch
2.00
066
11
1006
1g
PE(20)
1d
cll
PF(20)
of
1.0
10F
0.5
10
10
20
300
40
50
20
30
40
10
M
Epoclas
Epochs
10"
50
1d
PB(20)
Epochs
150
p-0.023
Epuchs
40
p-0.122
10"
30
10-1
20
60
10
10
40
M
M
OUT On
Of
On
Figure 2: A) Top-left: Relative proportion of subregions with positive determinant of the Jacobian (J > 0), S+, as a function of latent space dimensionality M for different proportions (Syl//S) Medians = interquartile range are shown. Top-right: Runtime and reconstruction quality (bottom)
as a function of latent space dimensionality M for different proportions |See /|S of subregions for which invertibility was enforced by regularization different training runs on the Lorenz-63 63 system + SD are shown. Reconstruction quality was assessed through (dis- Jagreement in attractor geometry (bottom-left; Datsp. see Appx. F and 20-step-ahee ahead prediction error (bottom-right the number of training epochs when the invertibility regularization, Eq. 5 was turned off (blue) vs. on (orange), for damped nonlinear oscillator. Bottom: Convergence to predefined performance criterion (PE(20) S 0.5) was significantly faster with the regularization turned on vs. off. Median across 20 trained models, error bands interquartile range. C) Same as B for Lorenz-63 system.
0.1 exp(M) in Eq. 5. Means across 100
B) Top: 20-step-ahead prediction error, PE PE(20). as a function of
(p 0.023, Mann-W Whitney U-test) than without the regularization (Fig. 2B), while hardly affecting performance on other systems like the chaotic Lorenz-63 (Lorenz) 1963) system (Fig. 2C; p 0.122).
A
6.1
40h
5.0
2.5
0.9
2.5
0.6
6'
0.0
0.3
11
0.00
2.5
0.0
2.5
-0.3
5.01
40
5.00
0
-2.
(1.1
2.5
5.0
-0.5
01.31
01.5
3.0
40
221
220
40
51
2.5
0.3
2.5
5.0
Figure 3: Model validation. A) Two examples of saddle points (half-green) with stable (gray solid lines) and unstable (red solid lines) manifolds determined by our algorithm, and points (black/ red dots, respectively) sampled by the analytical resp. backward/ forward map, showing that these all fall onto the analytically determined manifolds. B) Basins of attraction (light green, confirmed by sampling initial conditions and tracing their trajectories) of stable fixed point (green dot) delineated by the stable manifold (black) of a 4-cycle (left) or 3-cycle (right).
4 DELINEATING BASINS OF ATTRACTION
Basic methods validation We first validate our algorithm on a simple toy example, a two- dimensional PL map for which we have analytical forms for the inverse and the fixed points (see Appx. H.1 for details). We chose parameters to produce simple test case with only a single saddle point. Fig. 3/ confirms that points produced by backward (resp. forward) iterating the map all lie on the stable (resp. unstable) manifold as determined by our algorithm. In this simple 2d example, the manifolds correspond to line segments. Changing the PL map's 's parameters slightly, we obtain stable fixed point coexisting with a saddle period-4 (Fig. BB. left) or period-3 (Fig. BB, right) cyclic point (Gardini al.) 2009). Tracing back the stable manifold of the period-4 period-3 saddle using our algorithm, we obtain the boundaries of the basins of attraction of the resp. fixed point (Fig. B). Fig. 33 further confirms this solution for the basin perfectly agrees with the 'classical' (not
7

---

Published as a conference paper at ICLR 2026
scalable) numerical approach of drawing initial conditions grid in state space, and observing their behavior the limit 0c.
In these simple 2d test cases the quality of the manifold reconstruction easily be verified visually. In higher dimensions, however, this will no longer be possible, and so we also introduce a simple metric to quantify the reconstruction quality. For point æt E Lf sampled randomly either from neighborhood U/ containing the relevant data range (see Appx. for details), or from the reconstructed manifold of saddle point p, 20 W nu, with (-1,+1} indicating whether considering the unstable or stable manifold, we compute
8 (mo) == minka20 || |F (mo) - pll
E [0, 1],
[D0 pll
For points randomly sampled from or from W (p) nu for the example above, Fig 66 in Appx G shows the distribution of 8G confirming & z 0 for points sampled on the manifolds, while off the manifold & naturally spreads across broader range. For the latter we assign an index Tu (
U) ]mmx] E {0, 1) with 8max max(8. (2o (P(p) N2)]]- In the following, report A, (Iu) 8 (xo (')p) nu) quality statistic, where (+) is the mean and 8, the median.
A
2
40
Z3
20
-2
3
-2
A2 A 1
72
8
Z1
B
Z1
2
2
Z
10
20
Z2
0.5
0.4
potential [mV]
23
-20
3.15
-30
22
-50
21
100
200
300
400
-1.0
-0.5
0.0
0.5
1.0
1.5
time
Z
Figure 4: A) Reconstruction of the Duffing system (ShPLRNN, M = 2, H- 10). Ground truth trajectories in blue, identified fixed points in green, and in gray the stable manifold of the saddle separating the two basins of attraction as determined by algorithm B) 3d subspace of the state space of an ALRNN (M 15, P = 6) trained on a 2-choice decision making task. Two point attractors (green) were identified, with the stable manifold (black/gray) of saddle (half-green) in the center separating the basins. The basin boundary consists of different segments, projected down from truly 15d space into 3d subspace (accounting for some of the 'folded' appearance) and approximated linearly for visualization. C) Reconstruction for a shPLRNN (M = 3, H 20) trained on the Lorenz63 system. The system has two saddle-spirals in the center of the lobes and saddle at the bottom. In black and red are the stable and unstable manifolds, respectively, of the saddle as identified by our algorithm (left), while on the right as computed by numerical continuation of the original system. The close agreement indicates the shPLRNN has correctly recovered the state space structure, although having been trained on data from the actual attractor only. D) ALRNN (M = 25, P= 6) trained on electrophysioloiica recordings. Left: Time series of membrane voltage (true: gray, simulated: black); right: 2d projection of the state space with stable manifold of a saddle (black) separating the basins of attraction of stable fixed point (green) and the 38-cycle (orange) corresponding to the spikes. Note that the true stable manifold is 24d curved object, which for visualization purposes is represented here by locally linear approximation in the shown 2d subspace.
Bistable Duffing system The Duffing 1918) system is a simple 2d nonlinear oscillator that can exhibit bistability between two spiral point attractors in certain parameter regimes (see Appx. H.3 for
8

---

Published as a conference paper at ICLR 2026
more details). We trained shPLRNN (M 2, H = =10; Hess et al. (2023)) this system using sparse teacher forcing (Mikhaeil et al. 2022), and used SCYFI (Eisenmann et al. 2023) to determine the fixed points. The basin boundary between the two spiral point attractors is the stable manifold of a saddle node in the center (Fig. AA), and - as computed by our algorithm - agrees with the trajectory flows of the true system in blue, as further confirmed by A. ~ 0.97 (Appx. Table 2 Fig Fig(
Multistable choice paradigm Simple models of decision making in the brain assume multistability between several speoccc attractor states, to which the system's 's state is driven as one or the other choice materializes (Wang) 2002). We trained an ALRNN (M = 15, P 6) (Brenner et al. 2024a) to perform a simple 2-choice decision making task taken from (Gerstner et al. 2014), and, as before, use SCYFI to find fixed points and Algorithm to determine the stable manifold of saddle separating the two basins of attraction corresponding to the two choices. The basin boundary consists of different planar and curved pieces and is approximately visualized in Fig. B in a 3d subspace (of the 1 5d system), together with the reconstructed system's 's fixed points (green) and some trajectories of the true system in blue (A, = 0.95 in this case, Appx. Table & Fig. FC). While visualization becomes tricky in these bigher-dimensional cases, computing distances to the manifold, using the construction in Appx. D provides further information: In this example, find that both stable fixed points are about equally far from the basin boundary (d ss 0.34 vs. d =s 0.32), giving the two possible choices about equal weight.
Lorenz-63 attractor The Lorenz-63 model of atmospheric convection is probably the most famous
example of a chaotic attractor. We reconstruct this system with shPLRNN (M = 3. H 20). Besides the chaotic attractor (blue) and two unstable spiral points (green), the system has saddle node for which we compute the stable and unstable manifolds (A, Rs 0.78, Appx. Table 2 & Fig 6D). Fig. C confirms that the manifolds computed by Algorithm for the trained shPLRNN (left), agrec even well with those determined by numerical integration of the original Lorenz ODE system (right). This example also proves that the shPLRNN faithfully captured the geometrical structure of the state space, beyond just reconstruction of the chaotic attractor itself.
Empirical example: Single cell recordings In Fig. 4D we trained an ALRNN (M = 25, P= 6) on membrane potential recordings from a cortical neuron (Hertäg et al. 2012). The trained ALRNN contains a 38-cycle which corresponds to the rhythmic spiking activity in the real cell. In addition has a stable fixed point and saddle whose stable manifold (determined by our algorithm, A, ~ 0.76) separates the stable cycle from the stable fixed point, as illustrated in Fig. 4D (right), where we visualized the manifold through locally linear approximation. Although we can compute the full global 24-dimensional manifold, its inherent curvature in the 25 dimensional state space makes impossible to visualize it directly. Computing distances to the manifold (Appx. D), however, we find that the oscillation at its lowest point is less sensitive to small perturbations than the cell's 's resting state, with minimum distance of d as 2.3 to the basin houndary, compared to d as 1.2 for the fixed point. Many types of cortical cells exhibit this type of bistability between spiking activity and a stable equilibrium near the resting or a more depolarized potential (Izhikevich) [0077] Durstewitz and Gabriel 2007), and this example demonstrates how our algorithm can be utilized to reveal the structure of the state space supporting this type of dynamics from real cells.
HOMO-/HETEROCIING ORBITS AND DETECTION OF CHAOS
The stable and unstable manifolds can be used to identify homoclinic and heteroclinic orbits, as defined in Sect. 3.2 (Def. 3 D. Intersections between the stable and unstable manifolds of a saddle p lead to homoclinic points, or to heteroclinic points of two saddles p + q. The existence of such points inevitably gives rise to a complex fractal geometry (a so-scalled horseshoe structure, see Appx. ii and thus chaos due to the mmale-Birkhofff Homoclinic Theorem (Wiggins) 1988, Simpson 2016). Finding such intersections is therefore highly illuminating for determining the dynamical behavior of a system. This is illustrated for a simple 2d PL map in Fig. 5/A which shows the homoclinic intersections identified by Algorithm [ For this 2d case, we can in fact analytically determine the presence of homoclinic points, as worked out in Appx. H.2 the results of which agree with Algorithm I Fig. [5B illustrates the whole resulting chaotic attractor, which happens lie on the unstable manifold in this case. Fig. 5C provides the bifurcation diagram as one varies the model's bias term as control parameter, and Fig. 5) the system's 's two Lyapunov exponents across the chaotic range
to

---

Published as a conference paper at ICLR 2026
of fı (confirming the presence of 'robust chaos', the Lyapunov exponents do not change across the chaotic regime (Banerjee et al., 1998)). Finally, in Appx. Fig. [10 we show the identification of heteroclinic intersection by our algorithm in tighddimessionl empirical example, human electrocardiograms (ECG), thus confirming the presence of chaos in this signal.
A
2
-1
-2
-2
1
-1
0
1
0.4
-02 0.2
0.4
-0.4
-0.2
0.2
h4
Figure 5: A) Stable (black) and unstable (red) manifolds of a saddle point (green dot) and their homoclinic intersections as identified by our algorithm. B) Structure of the chaotic attractor caused by these homoclinic intersections. C) Bifurcation diagram function bias parameter h1. D) Lyapunov exponents across the h1 -range for which the chaotic attractor exists.
6 CONCLUSIONS
Here we presented novel semi-analytical algorithm for determining stable and unstable manifolds of fixed and cyclic points in ReLU-basee RNNs. Within each linear subregion, the algorithm provides an analytical construction for the manifold segments, and then recursively assembles these segments into the global un-/stable manifolds by forward- or backward-iteratnng the RNN map Fg across subregions. These manifolds are profoundly important for studying an RNN's dynamical repertoire, illuminating dynamical mechanisms in underlying system reconstructed by the RNN, or the mechanisms by which an RNN solves a given task. Stable manifolds of saddle points segregate the state space into different basins of attraction, giving rise to the computationally important property of multistability (Durstewitz et al. 2000; |eudee et al. 2018). Intersections of stable and unstable manifolds, in heteroclinic points which produce fractal geometry and chaos (Wiggins) 1988). While here we developed these algorithms for the class of PLRNNS (Brenner et al. 2022
turn, lead to homo-
Hess et al.
2023
Brenner et al. (2024)), extensions to other PL systems, like SLDS (Linderman
et al. 2017 2016 or TLS (Parmelee et al. 2022. Morrison et al. 2016 Hahnloser and Seung, 2000) conceivable. SLDS, unlike PLRNNs, are not continuous across boundaries of linear regions, however, rendering them less suitable for DSR in general and making it more tricky to impose invertibility. TLNs, on the other hand, are defined in continuous rather than discrete time. Techniques for converting discrete into continuous time PLRNNs may thus be of help here (Monfared and Durstewitz] 2020). Sometimes restrictions are imposed on the weight matrices (like binary weights or mutual inhibition), leading nice theoretical results (Curto and Morrison [2233 Parmelec et al.) 2022 Morrison et al. (0166) but also curtailing the model's flexibility for approximating arbitrary DS.
Limitations arise in the context of chaotic dynamics, where invariant manifolds may fold into fractal structures. While points from these manifolds may still be sampled, the analytic construction through curved/planar segments spanned by support vectors breaks down, as the intricate, self-similar geometry of fractals cannot be captured this way. Nevertheless, we still may be able to retrieve some important dynamical characteristics by determining homo- or heteroc linic intersections, as discussed in Section [5] Another limitation is that in the worst case scenario the algorithm may scale as 2" with the number of linear subregions P. However, shown in (Brenner et al. 2024a), the number of subregions utilized by trained PLRNNs quickly saturates, suggesting an at most polynomial scaling if one restricts attention to the domain explored by the data (see also Fig. 99. Finally, although not strictly part of the present algorithm, whether one can recover all relevant manifolds also depends one whether all saddle points were detected in the first place, of course. In any case, we emphasize that this is the first algorithm for detecting un-/stable manifolds in ReL.Ubaaedd RNNs.
Limitations
Code /
folds
10

---

Published as a conference paper at ICLR 2026
7 ACKNOWLEDGEMENTS
This work was primarily supported by Samsung Advanced Institute of Technology, Samsung Elec- tronics Co., Ltd. Additional funding was provided through the German Research Foundation (DFG) within the FOR-5159 research cluster (Du 354/14-1) and the collaborative research center TRR 265, subproject A06. Z.M. also grateful the Bundesministerium für Forschung, Technologie und Raumfahrt (BMFTR, Federal Ministry of Research, Technology and Space) for funding through project OIDLITDSM, No. 011S24061.
REFERENCES
K. T. Alligood, T. D. Sauer, and J. A. Yorke. Chaos: An Introduction to Dynamical Systems. Springer, New York, NY, 1996.
G. V. Amaral, C. Letellier, and Luis A. Aguirre. Piecewise affine models of chaotic attractors: the Rossler and Lorenz systems, Chaos (Woodbury, N.Y.), 16(1), March 2006. ISSN 1054-1500. doi: 1.11063/.21149527
V. Avrutin, I.. Gardini, Sushko, and Tramontana. Continuous And Discontinuous Piecewise- smooth dineddimessonaal Maps: Invariant Sets And Bifurcation Structures. World Scientific, May 2019. ISBN 978-981-12-0471-5 Google-Boo Boall-BBook-IDD ltacDwAAQBAJ.
Banerjee, A. Yorke, and C. Grebogi. Robust chaos. Physical Review Letters, 80, 1998.
A. Bemporad, F. Borrelli, and M. Morari. Piecewise linear optimal controllers for hybrid systems. In Proceedings of the 2000 American Control Conference. ACC (IEEE Cat. No. 00CH36334), volume 2, pages 1190 1194. IEEE, 2000.
N. Bertschinger and T. Natschläger. Real-time computation at the edge of chaos in recurrent neural networks. Neural computation, 16(7):1413-1))) 1e111111114,,,,, 2004.
A. Bill, C. Guiver, H. Logemann, and S. Townley, Stability of nonnegative lur'e systems. SIAM Journal on Control and Optimization, 54(3):1176-122 2016.
M. Brenner, F. Hess, J. M. Mikhacil, L. F. Bereska, Z. Monfared, P. Kuo, and D. Durstewitz. Tractable Dendritic RNNs for Reconstructing Nonlinear Dynamical Systems. In Proceedings of the 39th International Conference Machine Learning, pages 2292 2320. PMLR, June 2022. ISSN: 2640-3498.
M. Brenner, J. Hemmer, Z. Monfared, and D. Durstewitz. Almost-linear rnns yield highly inter- pretable symbolic codes in dynamical systems reconstruction. In A. Globerson, L. Mackey, D. Belgrave, A. Fan, U. Paquet, J. Tomezak, and C. Zhang, editors, Advances in Neural Information Processing Systems, volume 37, pages 36829-36868. Curran Associates, Inc., newrips ccppaper.cpppppppeeeeee http://pprr_flllssppeeee//////// pddf
2024a. URL https::
M. Brenner, F. Hess, G. Koppe, and D. Durstewitz. Integrating multimodal data for joint generative modeling of complex dynamics. In Ruslan Salakhutdinov, Zico Kolter, Katherine Heller, Adrian Weller, Nuria Oliver, Jonathan Scarlett, and Felix Berkenkamp, editors, Proceedings of the 41 International Conference on Machine Learning, volume 235 of Proceedings of Machine Learning Research, pages 4482 4516. PMLR, 21-27 Jul 2024b. URL URRLtttts///poceeiiii .mrr html
press/v235/brene
V. Carmona, E. Freire, E. Ponce, and F. Torres. On simplifying and classifying piecewise-linear systems. IEEE Transactions on Circuits and Systems I: Fundamental Theory and Applications, 49 (5):609-62)) 2002.
M. Casdagli. Nonlinear prediction of chaotic time series. Physica D: Nonlinear Phenomena, 35(3): 335-356, 1989.
R. Cestnik and M. Abel. Inferring the dynamics of oscillatory systems using recurrent neural networks. Chaos: An Interdisciplinary Journal of Nonlinear Science, 29(6): 06 2019. ISSN 1054. 1500. doi: N0.1063/1.55 5096918. URL /ttt::/://aoi.orr/001103///////3

---

Published as a conference paper at ICLR 2026
S. Coombes, M. Şayli, R. Thul, R. Nicks, M. A. Porter, and Y. M. Lai. Oscillatory networks: Insights from piecewise-lincar modeling, SIAM Review, 66(4):619-679) 17-19999-9999 2024. doi: 10.1137/22M15544365 URL https:://d...og//.......
A. C. Costa, Ahamed, and G. Stephens. Adaptive, locally linear models of complex dynamics. Proceedings of the National Academy of Sciences, 1166))1501111100 2019.
C. Curto and K. Morrison. Graph rules for recurrent neural network dynamics. Notices of the American Mathematical Society, 70(04), 2023.
A. Dovzhenok and A. S. Kuznetsov. Exploring neuronal bistability at the depolarization block. PLoS One, 7(8):42211, I, August 2012.
G. Duffing. Erawungene Schwingungen bei veränderlicher Eigenfrequenz und ihre technische Bedentung. Number 41-42. Vieweg, 1918.
D. Durstewitz. state space approach for piecewiselinear recurrent neural networks for reconstructing nonlinear dynamics from neural measurements. PLoS Computational Biology. 13(6), 2017.
D. Durstewitz and T. Gabriel. Dynamical basis of irregular spiking in amda-driven prefrontal cortex neurons. Cerebral cortex, 17(4):994900, 2007.
D. Durstewitz, J. Seamans, and T. Sejnowski. Neurocompuuational models of working memory. Nature Neurosci, 3:1184-1191, 2000.
D. Durstewitz, G. Koppe, and M. Thurm. Reconstructing computational system dynamics from neural data with recurrent neural networks. Nature Reviews Neuroscience, 24(11):693-71; 2023.
L. Eisenmann, Z. Monfared, N. Göring, and D. Durstewitz. Bifurcations and loss jumps in RNN training. In A. Oh, T. Neumann, A. Globerson, K. Saenko, M. Hardt, and S. Levine, editors, Advances in Neural Information Processing Systems, volume 36, pages 70511-70547. Curran Associates, Inc., 2023.
M. Ester, H. Kriegel, J. Sander, X. Xu, et al. A density-based algorithm for discovering clusters in large spatial databases with noise. In kdd, volume 96, pages 226-231, 1996.
U. Feudel, AN. Pisarchik, and K. Showalter. Multistability and tipping: From mathematics and physics climate and train-minirevieww and preface to the focus issue. Chaos, 28, 2018.
Fujisawa, E. Kuh, and Ohtsuki. sparse matrix method for analysis of piecewise- lineer resistive networks. IEEE Transactions Circuit Theory, 19(6) -585844 1972. doi: 10. 109/TCT.1972. 1083550.
V. Ganatra and S. BBanerjee. Sketching 1d stable manifolds of 2d maps without the inverse. Interna- tional Journal of Bifurcation and Chaos, 32(08):225011 11, 2022.
L. Gardini, V. Avrutin, M. Schanz, et Connection between bifurcations on the Poincaré Equator and the dangerous bifurcations. 2009.
W. Gerstner, W. M. Kistler, R. Naud, and L. Paninski. Neuronal dynamics: From single neurons to networks and models of cognition. Cambridge University Press, 2014.
M. Golub and D. Sussillo. FixedPointFinder: A Tensorflow toolbox for identifying and characterizing fixed points in recurrent neural networks. Journal of Open Source Software, 3(30):0033, November 2018. ISSN 2475-9066.
A. Gu and T. Dao. Mamba: Linear-time sequence modeling with selective state spaces, arXiv preprint arXiv:2312.00752. 2023.
J. Guckenheimer and P. Holmes. Nonlinear oscillations, dynamical systems, and bifurcations of vector fields, volume 42. Springer Science Business Media, 2013.
R. Hahnloser and H. Seung. Permitted and forbidden sets in symmetric threshold- -incar networks. Advances in neural information processing systems, 13, 2000.
12

---

Published as a conference paper at ICLR 2026
B. Hasselblatt and A. Katok. Handbook of dynamical systems. Elsevier, 2002.
L. Hertäg, J. Hass, T. Golovko, and D. Durstewitz. An approximation the adaptive exponential integrate-and-fing neuron model allows fast and predictive fitting to physiological data. Frontiers in computational neuroscience, 6:62, 2012.
Hess, Monfared, M. Brenner, and D. Durstewitz. Generalized teacher forcing for learning chaotic dynamics. In Proceedings of the 40th International Conference on Machine Learning. volume 202 of Proceedings of Machine Learning Research, pages 13017-13049. PMLR, 23-29 Jul 2023. httphst:s//proeeeeel
M. W. Hirsch, S. Smale, and R. L. Devaney. Differential equations, dynamical systems, and introduction to chaos. Academic press, 2013.
T. Hofmann, B. Schölkopf, and A. J. Smola. Kerel methods in machine learning. 2008.
I. Huh, E. Yang, S. J. Hwang, and J. Shin. Time-reversal symmetric ode network. Advances in Neural Information Processing Systems, 1902199919907,, 2020.
A. R. Ives and V. Dakos. Detecting dynamical changes in nonlinear time series using locally linear state-space models. Ecosphere, 366):1-15, 2012.
E. M. Izhikevich. Dynamical systems neuroscience. MIT press, 2007.
H. Kantz and T. Schreiber. Nonlinear time series analysis, volume Cambridge university press, 2004.
G. E. Katz and A. Reggia. Using Directional Fibers Locate Fixed Points of Recurrent Neural Networks. IEEE Transactions on Neural Networks and Learning Systems, 29(6):363633646 August 2018. ISSN 2162-2388. doi: 10.1 Transactions on Neural Networks and Learning Systems.
Conference Name: IEEE
G. Koppe, S. Guloksuz, U. Reininghaus, and D. Durstewitz. Recurrent neural networks in mobile sampling and intervention. Schizophrenia bulletin, 45(2):272-276, 2019.
B. Krauskopf, H. M. Osinga, and Jorge Galán Viidne Numerical continuation methods for dynamical systems, volume 2. Springer, 2007.
Linderman, M. Johnson, A. Miller, R. Adams, Blei, and L. Paninski. Bayesian Learning and Inference in Recurrent Switching Linear Dynamical Systems. In Proceedings of the 20th International Conference on Artificial Inteliiencee and Statisties, pages 914 914-922. PMLR, April 2017. URL 2640-3498.
ngs mlr. www.mm............ wwwppprsss//// http:///////ineeeeee rman 7a. html ISSN:
W. Linderman, A. C. Miller, R. P. Adams, D. M. Blei, L. Paninski, and M. J. Johnson. Recurrent switching linear dynamical systems. arXiv:1610. 161.0066 [stat], October 2016. URL http:
//axxiv.org/bbe.. http://ab/161.....00 08466
arXiv: 1610.08466.
E. N. Lorenz. Deterministic nonperiodic flow. Journal of atmospheric sciences, 20(2):130-141, 1963.
L. McInnes, J. Healy, and Astels. hdbscan: Hierarchical density based clustering. The Journal of Open Source Software, 2(11), mar 2017. doi: 10.21105//ss.000255. URL https://doi. org/
20.2110582F joss. C0205
J. M. Mikhaeil, Z. Monfared, and D. Durstewitz. On the difficulty of learning chaotic dynamics with RNNs. In Advances Neural Information Processing Systems, 2022.
Monfared and Durstewitz. Transformation of ReL.U-based recurrent neural networks from discrete-time to continuous-time. In International Conference on Machine Learning, pages 6999- 7009. PMLR, 2020. ISSN: 2640-3498.
K. Morrison, A. Degeratu, Itskov, and C.a Curto. Diversity of emergent dynamics in competitive threshold-linear networks. arXiv preprint arXiv: 1605 1605.0446, 2016.
13

---

Published as a conference paper at ICLR 2026
J. Nassar, S. Linderman, M. Bugallo, and I. M. Park. Tree-structured recurrent switching linear dynam- ical systems for multi-scale modeling. In International Conference on Learning Representations, 2019.
A. Orvieto, S. L. Smith, A. Gu, A. Fernando, C. Gulcehre, R. Pascanu, and S. De. Resurrecting recurrent neural networks for long sequences. In International Conference on Machine Learning, pages
PMLR, 2023.
H. M. Osinga. Computing global invariant manifolds: Techniques and applications. In Proceedings of the International Congress of Mathematicians, volume 4, pages 1101-1123 2014.
M. Pals, A. E. Sagtekin, Pei, M. Gloeckler, and Macke. Inferring stochastic low-rank recurrent neural networks from neural data. Advances in Neural Information Processing Systems, 37: 18225-18264, 2024.
H. Park, B. Kim, C. M. Kang, C. C. Chung, and J. W. Choi. Sequence-to-sequene prediction of vehicle trajectory via Istm encoder-decoder architecture. In 2018 IEEE Intelligent Vehicles Symposium (IV), pages 1672-1678, 2018. doi: 10. 109//VS 20120001118000
C. Parmelee, S.a Moore, K. Morrison, and C. Curto. Core motifs predict dynamic attractors in combinatorial threshold-linear networks. PloS one, 17(3)ee0244456, 2022.
D. Patel and E. Ott. Using machine learning anticipate tipping points and extrapolate post-tipping dynamics of non-stationary dynamical systems. Chaos: An Interdisciplinary Journal of Nonlinear Science, 33(2). 2023.
M. Patra and S. Banerjee. Robust chaos in 3-d piecewise linear maps. Chaos: An Interdisciplinary Journal of Nonlinear Science, 28, 2018.
B. Peng, E. Alcaide, Q. Anthony, A. Albalak, S. Arcadinho, Biderman, H. Cao, X. Cheng, M. Chung, M. Grella, al. Rwkv: Reinventing rnns for the transformer era. arXiv preprini arXiv:2305. 13048, 2023.
U. Pereira-Obillioon J. Aljadeff, and N. Brunel. Forgetting leads to chaos in attractor networks. Phys. Rev. X, 13:01 1009, Jan 2023. doi: 1..100%Phy-Reex.133.110000 URL https:///iit onssoorg/oi//10 htt.d////11..1110////eeeeee 13. 011009
L. Perko. Differential Equations and Dynamical Systems. Springer New York, NY, 2001.
A. Platt, G. Penny, T. A. Smith, Tnn Chen, and H. D.I. Abarbanel. Constraining chaos: Enforcing dynamical invariants in the training of recurrent neural networks. arXiv preprint arXiv:2304. /2865, 2023.
M. Rabinovich, R. Huerta, and G. Laurent. Transient dynamics for neural processing. Science, 321 (5885):48-50 2008
I. Roy, M. Patra, and S. Banerjee. Shilnikov-type dynamics in three-dimensional piecewise smooth maps. Chaos Soliton. Fract., 133, 2020.
T. K. Rusch and Mishra. Coupled oscillatory recurrent neural network (cornn): An accurate and (gradient) stable architecture for learning long time dependencies. arXiv preprint arXiv:2010.00955 2020.
K. Rusch, S. Mishra, N. Erichson, and M. W. Mahoney. Long expressive memory for sequence modeling. In International Conference Learning Representations.
D. Schmidt, G. Koppe, Z. Monfared, M. Beutelspacher, and D. Durstewitz. Identifying nonlinear dynamical systems with multiple time scales and long-range dependencies. In proceedings of the Ninth International Conference on Learning Representations, ICLR, 2021.
B. Schölkopf, A. Smola, and K. Müller. Nonlinear component analysis kernel eigenvalue problem. Neural computation, (255)112991133 1998.
H. Siegelmann and E. D Sontag. On the computational power of neural nets. In Proceedings of the fifth annual workshop Computational learning theory, pages 440-449 1992.
14

---

Published as a conference paper at ICLR 2026
D. J. W. Simpson. Unfolding homoclinic connections formed by comer intersections in piecewise- smooth maps. Chaos: An Interdisciplinary Journal of Nonlinear Science. 26(7), 2016.
D. J. W. Simpson. How to compute multi-dimensional stable and unstable manifolds of piecewise- linear maps. arXiv preprint arXiv:2310.0900 09944, 2023.
M. Storace and De Feo. Piecewise-linear approximation of nonlinear dynamical systems. IEEE Transactions Circuits and Systems I: Regular Papers, 3144:830-842 2004.
L. Su, P. Seiler, J. Carrasco, and Sei Z. Khong. On the necessity and sufficiency of discrete-time -shea-zames-aee falb multipliers. Automatica, 150: 100:11000 10872, 2023.
S. Waitman, L. Bako, P. Massioni, G. Scorletti, and V. Fromion. Incremental stability of lur'e systems through piecewise-affine approximations. IFACC-PppesOnnlinne, 50(1):1673-1679, 2017.
X. Wang. Probabilistic decision making by slow reverberation in cortical circuits. Neuron, 36(5): 955-966 968, 2002.
X. Wang. Decision making in recurtent neuronal circuits. Neuron, (002):215-234
S. Wiggins, Global Bifurcation and Chaos. Springer-Verlage New York, 1988.
15

---

Published as a conference paper at ICLR 2026
APPENDIX:
A LLM USAGE AND CODE AVAILABILITY
Code for the algorithms developed here and results produced is available at /ltttt:://ggtttt com/Durstewit /aaaa//Deectt ingMani folds Large Language Models (LLMs) were used for implementation of standard routines and assistance with literature search.
B DERIVATION OF EQ.
The stable and unstable manifolds of saddles are invariant sets with all trajectories with initial conditions the manifold confined to the manifold for all (Alligood al. 1996 Perko 2001), Fo(W*) C Wa We recursively assemble these manifolds following Eq. B with points sampled from T following Eq. within each linear subregion. Here we derive this equation starting with linear system of the form 2++1 Lz;-
Case (1): If the geometric multiplicity equals the algebraic multiplicity, then L is diagonalizable and admits basis of linearly independent eigenvectors "1 Un. Thus, the orbits of the system can be expanded in terms of eigenvectors
= + 022//2 + +... c... C.
(7)
Case (2): If the geometric multiplicity of X is less than the algebraic multiplicity, then the eigenvalue said be defective and admits basis of generalized eigenvectors. this case, has Jordan block (A;); of size d > 1
N
Ja(h.)
0
with the 1-th power
A
() -
( X
( X a
0
X
(9 1
(x x
(N)
:
0 0
X
( x-1
0
X
Moreover, a generalized eigenvector wa * 0 of degree d, corresponding to the defective eigenvalue A₁ satisfies
- = but (L - w. O,
and L has d linearly independent generalized eigenvectors associated with A In fact, construct a chain generalized eigenvectors (wn, wa) such that
(L - Af) = wa-1, (L - AlJwa- = w
(L-AI)wy = wi.
(9)
where un is regular eigenvector. Given chain of length d, the orbit contribution corresponding to the Jordan block Ja();) is given by
12
dd-
defective 2
X (((c1u) C2 tw. + c3 ws
Cd (d - 1)! wa
(10)
2!
The full orbit of the system is linear combination of contributions from the non-defective eigneveaactiveeeeeealluuuus dy. < j # i n, and the defective eigenvalue A; according to
2,
Σ
Nyfo X
Σ
Cr (r - 1)! Wr.
(11)
16

---

Published as a conference paper at ICLR 2026
This formula describes exponential evolution for all non-defective eigenvalues and polynomially modified exponentials for the defective eigenvalue Aę.
PLRNN is not strictly lincar, but affine in each subregion, see Eq. 2 For an affine system of the form 2t+1 Lz, h, Eq. T4 consists of homogeneous part (determined by the eigenstructure of L) and a particular part due to the constant bias term h. In the presence of a defective eigenvalue X; with Jordan block of size d, the full orbit is given by
tr Cr (r 1)1
21 Σ GDS
X
E
Wr + Σtth.
L'h.
(12)
The final term Σt=1 h accounts for the cumulative effect of the bias, modifying the orbit away from purely exponential or polynomial-exponentaal behavior.
C ADDITIONAL DETAILS ON MANIFOLD CONSTRUCTION
Backtracking algorithm and other PLRNN variants Step 2 in sect. 3.3 outlines the steps for inverting the RNN map F& for PLRNNs of the form Eq. [ can easily be amended for other ReLU-based RNNs. For instance, for 1-hidden layer PLRNN, called the shallow PLRNN (shPLRNN) (Hess et al. 2023), given by
spelled out in detail in Algorithm 2 This
2141 Az, + W, D, (W2 h22 + h1
(13)
Algorithm 2 Backtracking Time Series in a ReLU based RNNs
1: <- an initial State
2: 0 <- Parameters
3: Initialize list: S [2T]
4: tort=:; do
5:
2 = S t]
6:
D,-1 + diag(2: 21 > 0)
Initialize as a diagonal matrix Perform backward step
7:
- F Dr
8:
& t --0 F(0,=111 =
D-Perform a forward step
9:
if t 2, then
10:
S t SU {2t-1}
>-If forward step is correct
else
12:
DI-1 + diag(z أا 0)
>Update D with new candidate Retry backward step Retry forward step
13:
4-1 (01(0,
2, x- F(0,2i11
15:
if ã = then
16:
+ SU {
b If forward step is correct
17:
else
18:
2 4- -rrryPrvviouRReiionn(. D D_pool, 24 18-11
> Try previous regions
19:
if 2 # then
20:
&, + Try Bitffiispps( 21, 21
b Hierarchically check neighbours
21:
if &, 2, then
22:
S - SU (
else
24:
return
end if
26:
27: end end if
28:
end if
29:
end for
30:
return trajectory
17

---

Published as a conference paper at ICLR 2026
where A E RMMM yields
W1 E R
W2 R h2 E R h1 E RM,
the inversion of this map
2 (A W:D: WiDi-1he h1)
(14)
Manifold modeling To generate initial conditions on the local manifold in Algorithm 1 we sample perturbations in the stable eigenspace around base point 2y (equal to p the initial linear subregion). To account for anisotropy in the vector field and potential undersampling of slow eigendirections, we first construct orthonormal basis ao Via the QR decomposition (E" = E Eo R) which preserves the span of the stable eigenspace while avoiding numerical degeneracy. Transforming the linear oper- ator into this basis yields the effective eigenvalues \; = (R diag(A),...... Ad, )) which govern contraction rates along the orthonormalized directions. Because systems with disparate timescales typically exhibit large variability along fast eigendirections but very limited variability along slow ones, naive sampling leads poorly conditioned point clouds. We therefore rescale perturbations inversely to the magnitude of the effective eigenvalues via = [A, where the hyperparameter c
controls how strongly slow directions are emphasized. We then draw N, perturbation coefficients from Gaussian mixture model (GMM), a ~ 3 E E N(, 1) o% {0.01, 0.1, 0.5), where the three scales σę ensure multi-scale coverage of the local manifold. Eigenvalue-basee rescaling is applied multiplicatively (ā; 8; a;)). Finally, the sampled points on the manifold are generated by (i.e., inconsistent with the active zx)) are rejected. This procedure yields a numerically well-conditioned set of initial points distributed across all stable eigendirections and
23 E Ex a. and points violating
multiple sampling scales.
In the case of curved manifolds, we used kernel-PCA (Schölkopf et al. 1998) for approximation in Algorithm Specifically, used RBF kernel, the squared exponential kernel K (ær, y) exp(-áll<> - )) for generally good approximation properties for smooth manifolds with curvature as described by 4 (Hofmann et al. 2008). Since K ((x, y) decays with Euclidean distance, kernel-PCA with an RBF kernel emphasizes local geometry, which is suitable for manifolds with curvature. Moreover, the RBF kernel is universal and capable of representing smooth curvature. While more specialized kernels (c.g., spectral periodic or polynomial-exponnttinn may sound like a theoretically better match, they tend to require careful tuning of system-specifie frequencies or decay rates. The RBF kernel, by contrast, offers a robust and broadly applicable alternative that performs reliably without region-dependent parameter adjustments.
spectral-periodic
function BACKWARDFORWARDCC.C D, 2: F F-(
3: F(O,2*)
4: return =²,²
5: end function
6: function TRYPPEVIOUSRRFIIONSS. D D_pooll 2,
7: for D E D D_pool do
8:
2¹, - BackwardForward((. D,2)
9: if 2 = 2 then return 2
10:
11: end
end for
13:
end function
function TRYBITHLIPS(O, 2,
15:
for k = 1: num mwm_relus do
D vereirisnnsss < generate_
17:
for D E D versions do
18:
x" + BackwardFoorward D,
19:
if 2 = then
20:
return ~
21:
end if end for
22:
23:
end for
24:
end function
18

---

Published as a conference paper at ICLR 2026
Fallback algorithm We are usually interested in the structure of the state space in the context of DS reconstruction where we train RNNs on time series observations from time- ume-continuous DS. In these cases, Algorithm[]] iis always applicable. More generally, however, maps may exhibit large jumps, with orbits "erratically" hopping around among subregions. In these cases the un-/stable manifolds may acquire more complicated, non-smooth geometry (as in Fig. 3 & Fig- 5 Algorithm formulates a 'fallback' procedure for such cases that works by perturbing seed points along the analytically defined local manifold, and then iterates I' generate a larger set of support vectors. As manifolds can re-enter the same subregion in multiple folds, we apply HDBSCAN (Mcinnes et al., 2017: Ester et al. 1996) to cluster support vectors into distinet segments. Although computationally more demanding. this fallback reliably captures manifolds with discontinuous or folding structures when sequential tracing as Algorithm| fails.
Algorithm Finding stable/unstable manifolds: fallback algorithm
1: (p, E)) <- SCYFI
>> p: Fixed Point, E: Eigenvectors
2: S(+1) 0
D-Stable Manifolds
3: S(-1) +- 0
>Unstable Manifolds
4: for i=1:N1 do
5:
20
=
For N₁ different initialisations
6:
for
0 111 F( do
7:
20 += nl-1) rand()
Perrurre into subspace
8:
end for
9:
[(-1) + Ge(ForvardTS(z))
10: S(-1) + -(-1) U {T(-1)}
end for
12: for
i=1:N2 do
20 =
For N₂ different initialisations
14:
for
+(+ E(+1) do
15:
20 += v1+1) rand()
b Perturb into subspace
16:
end for
17: T(+1) + Bet
18: S(+1) + s(+1) ([rl+1}}
19: end for
20: S(+1)
+ 0
Piecewise linear manifold fits
21: S(-1)
22: for each E Dn do
s ((+)] - s(+1) n Do
23:
D Go through all subregions
st-1l <- st-1l Do
24:
25:
(C(+1), ,,(--))) FIT((S(()) $Q -
+1)
Cluster points and fit
5[+1) + -51+ StHt
26:
5(-1) < 5(-1) (CC [h
28: end for
29: return (51+1), 51-11)
b Piecewise linear manifolds
D ALTERNATIVE CONSTRUCTION OF MANIFOLDS
Algorithm I constructs the local manifolds within each subregion using PCA or kernel-PCA. Here we provide an alternative parameterization that may be simpler use for some purposes, like calculating distances on and to the manifold. First note that Eq. 4 (cf. Appx. B) yields 1d curve per initial condition. invariant d-dimensional manifold (d > T) is smooth surface that contains infinitely many such curves. To describe the entire manifold, we need parameterization using d variables, rather than a family of I-dimensional curves.
Let the state space be partitioned into finitely many subregions J and
F(x) LLE mk,
Sk-
(15)
19

---

Published as a conference paper at ICLR 2026
Assume that there exists d-dimensional invariant manifold W in Ky S (1, ................................ K} subregions of F, and for every such subregion
W W +
E
(16)
We focus on a single fixed point p F, and assume that this fixed point lies in a unique subregion Ska. The invariant manifold W is the un-/stable manifold associated with this fixed point. Accordingly. the local manifold segment used initiate the construction centered at this fixed point. The manifold segments W k E IK, {ko),, in other subregions are obtained by propagating this local manifold under the dynamics of F. Thus, in all other subregions the local manifold segments are generally not centered at fixed point of the corresponding affine subsystem.
For each subregion Sg. choose reference point
क W W
(17)
which is assumed lie on the invariant manifold and serves local anchor of the manifold segment in subregion k. For instance, ø can be chosen as a point where single invariant 1d trajectory generated by intersets the subregion Sk. In the unique subregion that contains the fixed point, zi coincides with this fixed point.
Suppose that TĄ E GL R) is any invertible linear change of coordinates whose first d columns span the tangent space Tug W of the invariant manifold at x Denote by Ska, the unique subregion that contains the fixed point p and set x = p. p. the tangent space of the d-dimensional invariant manifold Wil1 111 is given by the corresponding invariant eigenspace of Lip (stable or unstable). Let
VED
(1) [kk
(d) Kg
E R.d
(18)
be basis matrix of this eigenspace. For any subregion Sk, let ot W be anchor point on the invariant manifold. Assume that there exists a finite sequence of points
Te = p. 31, 22
IN.. = o
(19)
such that
F+ F(#t);
1= 0, N
(20)
and x2 E W for all f. In particular, such a sequence can be obtained by iterating a point on an invariant 1-dimensional trajectory generated by Eq. A but the construction does not depend how the points on W+1 are generated, provided that they follow the dynamics. Then the tangent space of the invariant manifold at at is obtained by transporting the tangent space p along this trajectory
V (r(M--1) ... Lr( Vka-
(21)
The columns of V2 span Ta; W The matrix T₃ is obtained by taking the first d columns equal to the columns of V& and completing them with any additional (n - d) linearly independent vectors to form an invertible matrix.
d-dimensional un-/stable manifold segment W is itself a d-dimensional smooth surface embed- ded in _.. The Implicit Function Theorem guarantees that any d-dimensional smooth manifold that transversal to a complementary (n dimeensiona subspace, can be locally expressed as a graph
of smooth function h4 .dd Ra-d
We introduce local coordinates (u, v) by
x の + Th
u E HER, vER"-d
(22)
In these coordinates, the manifold segment in subregion & is represented by a graph
W = x + Th
u (hw(w))
u q Uk
(23)
with
(hx = 0. (0) 0,
(24)
20

---

Published as a conference paper at ICLR 2026
The condition (0) = indicates that the chosen coordinates are tangential to the invariant manifold at the anchor point z
For E Sy we have
Using the coordinate change associated with subregion k. mk
x = F(x) (-)x))a = 1 +
(25)
u
(26)
the image of (u, v) under the dynamics can be written as
f T
+ Th (L m+ xt/)
(27)
We denote the block decomposition by
T LTT
Ak
B& D2)
(L my xx/)
0.
(28)
Ck
bs
If both æ and F(x) belong to the same subregion Sk, then points on the manifold satisfy
hu),,,
ot hn(="))
(29)
Substituting (u) into Eq.
27 yields
Agu Bgh (u) + ag.
(30)
+ Cşu + D&hy (u) + b
(31)
Hence, the local invariance equation in subregion given by hx(Agu + Bahx(u) + ax) = Cgu + (xh by.
(32)
The affine terms ag and by appear only for at F p. that is, when the anchor at is not a fixed point of the local affine map.
Let a E W C Sy and
F(x)= E Sp:
i k.
(33)
In subregion j the manifold is represented by
n Tj
hz(
(34)
Invariance of the global manifold implies the existence of a local reparameterization map á (2)
(35)
such that
(u))((u) Lx ar + Ta ha (u())) )
2. + T3 h2 (R((
+ hy..
(36)
Eq. 36 couples the local graph parameterizations of neighboring subregions.
Differentiating Eq. 36 at u. 0 and using Dhg (0) 0, Dhy(0) = 0, yields
T, (DRs = L To
La
(37)
Since T& and T7 are invertible, DR (0) has full rank d. Therefore, the intrinsic dimension of the invariant manifold is preserved when passing from subregion & to subregion j. in agreement with the assumption that the manifold cannot change its dimensionality across subregions.
Global manifold. The invariant manifold of the piecewise- -ineer system is obtained as the union of the locally parameterized segments
K U 1=1
æ + Th ( (b)(w)
w+1
E Uk
(38)
with consistency between neighbouring subregions enforced by Eq. 36
21

---

Published as a conference paper at ICLR 2026
How to compute hú in practice We outline the standard series method which produces polynomial approximation of (u) for each subregion Sy. This method works for either stable unstable manifolds by choosing the corresponding tangent space basis T₃ and block matrices (Auj. Bkj. Ckj. Dxj), which link region & to its neighboring region j = r(k) following the manifold flow. The affine terms Ok; bR, appear when the anchor point oi is not fixed point of the local affine map.
Step 1: In each subregion Sg. develop he (u) as a multivariate Taylor series about u 0 (no linear term because tangency the first d columns TĄ at the anchor point ai.):
(u)
Σ 1a/22
1R,0 u",
(39)
where a E N is multi-index and u
m [1 Ui
Step 2: Substitute the series into the invariance equation linking region i to its image region j:
hy (Auzu+ u+ + (u)h+(u) + ax) = -CRzu Dajhu (u) by-
(40)
Expand both sides into multivariate Taylor series in u. For the left-hand side, first expand the argument U A, + Bxyly(u) + ag. then expand (U) using its own series expansion about U = O. Equate coefficients of monomials u" of homogeneous degrec |a to obtain equations for the unknown tensors a and hy.a.
Step 3: Solve order by order:
For order 1: The linear terms vanish automatically due Dh₂ 0, Dh2(0) and condition Ck₃ = 0 which must hold for the manifold to exist as a graph. This serves as consistency check.
For order 2: For |o:| = 2, the left-hand side contributes (A) u) 82 from the quadratic term of hy, plus terms involving ag and Bazhu(u). The right-hand side contributes Dxyhi u). This yields linear relation:
H (Auju) =Ds Dis H 82 + (terms from ax)
(41)
where H and H are the quadratic coefficient tensors. Here, (Akju) denotes the second tensor (Kronecker) power of the vector Akyu, and H ((2) (AgJu)
H (Arju, Atju) the evaluation of the quadratic coefficient tensor H(2) (i.e., the Hes-
sian of hy at u = 0) on the pair of vectors (Auy Ak u)): essence, this term 'measures" how the curvature of the manifold in region acts the direction obtained by mapping from region & into region j through the linear map Agi-
For orders é > The equation for hg involves:
Known lower-order coefficients of hy (from (u)hu(u) expanding the argument)
Known coefficients of hú up to order
(from expanding h; about U 0)
The affine offset as,, which shifts the expansion point and creates additional lower-order contributions
The unknown hi. appears linearly, multiplied by Dkj. yielding a linear system:
(hu, RHS, (known lower-order terms), where Lkj is a linear operator determined by Day and the eigenvalues of Auj-
(42)
Step 4: For each degree 2 2, solve the linear homological system for the coefficients he,a- Non- resonance conditions between the eigenvalues of Akj and Dx, ensure invertibility of the linear operator at each step. Resonances require special treatment (free parameters or solvability conditions), indicating that the manifold may not be unique may require higher-order analysis.
22

---

Published as a conference paper at ICLR 2026
Step 5: Repeat the process for all regions & visited by the manifold, ensuring consistency at region boundaries. The final result is a piecewise polynomial representation
W = { a0) + Ta (na((n)
U2 C R
(43)
for each subregion Sk, with the hy given by the computed Taylor coefficients.
E DETAILS ON RNN TRAINING
For DS reconstruction on the simulated and empirical benchmarks (Figs. 2 & 4. we used the shPLRNN (Hess ct al. 2023) and ALRNN (Brenner et al. (2024)) trained by sparse teacher forcing (Mikhaeil al. 2022),
-antrol-theoreticalyy motivated training method for ensuring that recon- structed systems exhibit the same long-term statistics (cf. Appx. as the underlying, data- anta-producing
DS. We use the standard MSE loss
(MSE X X)
M Σ dr. xe
(44)
to which the regularization Eq. 5 was added, where X = {x1TT) the training sequence of length T and X are the corresponding model predictions, obtained from the latent states through some observation function, d; = g(z,), simply taken be the identity in the settings considered here. Sparse teacher forcing replaces latent states zt by data-inferred states & = g (æx) every time steps via decoder model inversion. Training and experiments were performed on single CPU (Intel Xeon Gold 6132 with 512GB RAM or Intel Xeon Gold 6248 with 832GB RAM). with parameter configurations specified in Table The parameter settings are based on (Brenner al. 2024a). refined by grid search.
Parameter
Lorenz63 Fig. Z
Oscillator
Lorenz63
Duffing
Decision making
Lorenzñ3
Empirical
Fig. ALRNN
Fig. BC
Fig. 4
Fig.4b 4s ALRNN
Fig. 4
Fig. 4b
Model
ALRNN
ALRNN
shallow PLRNN
SPLERNN
ALRNN
M
10/20/30/50
40
30
2
15
3
25
#Rel.Us
M-3
15
6
6
Soquence length
200
200
Gaussian noise
0.0
0.0
0.01
0.05
0.02
Amern
0.0 I-epp-)
0.0/le15
0.0/le10
0.0
0.2
0.0
0.3
Batch Size
16
16
32
16
16
Epochs
10000
1000
1000
10000
20000
1000
2000
Learning rate
16
10
15
15
15
15
20
Table 1: Parameter configurations. teacher forcing interval.
F DS RECONSTRUCTION MEASURES
Measure for Geometrical Agreement (Dstp) Given probability distributions p(æ) across ground truth trajectories and ((x) across model-generated trajectories, Sssp defined as the Kullback-Leibh (KL) divergence
Dusop := (p(x) |l ((x))) = JXERN p(x) log p(xx) dx g(x)
(45)
For low- -ww-iiennioonl observation spaces, p(æ) and g(æ) can be estimated using binning (Koppe et al. 2019 Brenner et al. 2022), yielding the discrete approximation
(x(2:) 9x (x)
D = DKL (p(x)) ((xx) ~ Σ (x) log
(46)
23

---

Published as a conference paper at ICLR 2026
K m is the total number of bins, with bins per dimension. Px(x)) and (x)x)) are the normalized counts in bin & for ground truth and model-generated orbits, respectively. Here we used 30 bins per dimension
Prediction error (PE) The n-step prediction error is defined the mean squared error between ground truth data {a,} and n-step ahead predictions of the model { i.e.
PE(n)
1ettn Ez-n
(47)
N(T n)
G MANIFOLD QUALITY ASSESSMENT
Fig 6 iilustrates the sampling process used for assessing the manifold measure (left column), Eq. 6 and the distributions of this measure (right column) when sample points 2% from the manifold of a saddle p. Zo E Wo (p) nu, vs. generally from some neighborhood Zo E 2 as illustrated on the left. Table 2 furthermore reports for each of the examples used in the main paper the median 8, across values on and off the manifold (usually in tight agreement with the respective means), as well the average index (Iu) for points sampled from U, and the final quality statistic A, as defined in the main text.
We note that sometimes numerical issues can occur with applying Eq. 6 as the recursive iteration of F for large & will exponentially diverge along unstable directions (which we always have for saddle in both forward and backward time!) if 2t++ gets off track even so slightly. Eventually this will almost always the case for - 00, and hence the issue becomes particularly severe for stable eigendirections with absolute eigenvalues close to 1, since the slow movement in this case requires many iterations (during which errors accrue) to reach the saddle (as was the case for the empirical example in Fig. D). Another potential issue is that
points the unstable manifold may temporarily move close to the saddle under forward iteration the map; since we are looking for the least distance the fixed point p. this could further blur the differences between & measures of points starting on vs. off the manifold. A partial remedy might be to replace 6mas below Eq. 6 by 'softer' threshold, like the 99th percentile of the å distribution, which bit more lenient by allowing for some distributional outliers in the computation of A..
Table 2: Manifold Reconstruction Quality Statistics
On Manifold
Randomly
Example
Median 8,
Median 8,
(Iu)
A+1
Figure A (left): unstable:
1.6e
+ 1.6e
0.97 0.026
0.98
0.98
Figure A: Duffing
0.0066 + 0.002
0.46 + 0.25
0.97
0.97
Figure
Decision making task
0.00011 = 4.3e 52
0.36 # 0.18
0.95
0.95
Figure
Lorenz63
0.12 + 0.08
0.8 0.2
0.90
0.78
Figure D: Cortical neuron
6.5e 4.2e
0.59 0.16
0.76
0.76
24

---

Published as a conference paper at ICLR 2026
A
Nocmalized Counts
0.0
-25
-5.0
-2.
0.0
2.5
0.0
0.00
0.25
0,50
0.75
1.60
6.1
B
Nocnalized Counts
-2
6.
C
Norrallzeed Counts
8
6,1
D
Normalhed Counts
25
20
15
10
0.00
0.25
0.50
0.75
1.00
E
Normalized Counts
1.05
20
1.00
40
0.95
0.90
6,1
Figure 6: A) Left: Illustration the sampling for Fig. 3 (left), with 1000 samples (red dots) drawn from the unstable manifold (reddish) and 1000 points (orange dots) drawn randomly from box U as illustrated. Right: Histograms of & for 2o ((-1)(p) ni (red) and wo E U (orange). B) Left: Same as A for the bistable Duffing system from Fig. 4A with points (black) sampled from the stable manifold (gray) from environment 2/ (orange). Right: Histograms
in A. C) Left: Same as A for the choice model from Fig. 4B, with points (black) sampled from the stable manifold (gray) from environment l4 (orange). Right: Histograms
in A. D) Left: Same as A for the Loren63 system from Fig. 4 with points (black) sampled from the stable manifold (gray) vs. from an environment 2/ (orange). Right: Histograms potential recordings from Fig. 4D. with points (black) sampled from the stable manifold (gray) vs. from an environment 2 (orange). Note that this a 2d projection of 25-dimensional system, such that the orange box in this case only gives a rough idea of the sampling region surrounding the stable manifold. Right: Histograms as in A.
in A. E) Left: Same as for the membrane
25

---

Published as a conference paper at ICLR 2026
H BENCHMARK SYSTEMS FOR TESTING
H.1 ANALYTICAL TEST CASES
For our analytical test cases we chose example 2d PL maps from Gardini et al. (2009), which can reformulate precisely as a PLRNN, such that direct comparison with the ground truth is possible (i.e., without inference of the underlying DS as the numerical examples from Fig. 4). In Gardini et al. (209), the PL map considered has the 2d form
{ A X B, $ A, X + B, i
F(X)
(48)
with x the two-coordinate state vector, and the transformation matrices and bias vector
given by
A (-81 1). A. = (-8, a), ( B
h1 h2
This map can be reformulated as a PLRNN (Eq. 2 by defining the parameters A, ,h as:
A
Ti b1
2)
W
Tr T -8,
0). h
h11
+ 81
hz2
For the analytical test case Fig. 3A. we chose parameters: (left) T, = 1.2, 8, = 1.8. T2 = -0.3.
0.9, c = 11.5, 1.0 ./2 0.1, h1 -0.13, (right) T, = 1.26, 8, = 0.71, 7I -0.39,
5: = 0.91, c -0.44, d 0.56, h22 0.62, h1 = -0.28.
For the test case in Fig. 3
we adapt parameters from Figs. 3 & 4 in Gardini et al. (2009) as
follows: (left) T, - 1.85. 8, 0.9, TI -.3, 6, 0.9, 1, d 0, h12 0, h1 1. (right)
T 1.85, 8, 0.9, m = 0.9, 61 0.9, c= 1, 0, hı2 = 0, hi
For Fig. 5 we considered the following setting: Tr = 1.5, 8, = 0.75, -1.77, = 0.9, - 0.6.
0.15, h2 = 0.4, h1 = 0.7.
H.2 RECURSIVE ALGORITHM FOR DETERMINING HOMOCLINIC INTERSECTIONS
Here introduce another algorithm for finding homoclinic intersections in low-dimensional settings. We use this to validate our main algorithm in sect. 5 Specifically, to investigate the existence of homoclinic intersections for the map Eq. 488 we use an algorithm similar the one proposed in (Roy et al. 2020). As illustrated in Fig. 7 it is based on a recursive procedure as follows:
Let O2 = v2)) E L. an admissible saddle fixed point, and A. A, stable and unstable eigenvalues of Ac. Assume that e and e are the stable and unstable eigenlines generated by the associated stable and unstable eigenvectors, respectively. Analogous to
suppose that e,, hits the border at Po yo) E where yo is given by Eq. 60 2. Consider the image of Po and assume this point is on the right side of the border, i.e. Pı E S.
Sect. I.2
Since this point has the coordinate P1 (eyo + hi, yo + h2) follows that
cyo + h
(bih1 + (1- (Au d) - b2(ch2 + (1- ddh1)
(Au d)(1 - + De)
h1 0.
(49)
Moreover P1 - is the first fold point of the unstable manifold of Oz, and all its images will also be fold points.
3. Suppose that the orbit starting from Po will return to the border again. Let n be the the border return time defined as the minimum number of iterations needed for Po to cross the
26

---

Published as a conference paper at ICLR 2026
("
Σ
O2
Po
("(O))
Patl
P
Pa
Figure 7: Schematic diagram to illustrate the procedure for finding homoclinic intersections.
border and return the left side again, i.e., such that all iterations Pt, P2....,,PP P lie the right hand side while Pn L. Using recursive method similar to the one proposed in (Roy al. 2020), we can compute the n-th iteration, Pr, directly from Po. That is. there is no need to compute any other previous iterations P1,P,, P For this purpose, first we calculate Pn as
= = A Po + (AR 1) (A - h
(50)
Suppose that Ar has two distinct eigenvalues, then, according to Proposition the matrix A has the form Eq. 73 Hence, for Po (0, vo) we have
Anti dAn
CAn
(0) + PR((1) x
Pa
brAn
dAn Dr An-1
(10)
a h-d
/ - dAn -1
Dr Dr An-n-1 1 ) 111
cAn
br
1- arr
bnAn
dAn
(12)
(1-d)(Ant1-dA.. tob, An Nite
An, - Dr An
cAn yo
(x(1)
(51)
(dAn DrAa-1)un by (tt
h1 hr
PR(1)
where A, given by Eq. 74 and I R, DR. PR. are the trace, determinant and characteristic polynomial of Ar respectively; and
6" = (a, -a) DRAN-)) + (d- (dd-Dyy Di2)
(52)
Now, by Eq. 5T we can compute all iterations and thus find the border return time n for which the unstable manifold passes the border again.
4. Finally, we calculate P and check whether or not the points Pa+1 and P.. are on opposite sides of the stable eigenline Ex i.e. whether we obtain a homoclinic intersection. For this, let P, = (In. Un) TT Pat1 = (In+1, 3T+1) and consider L(x, y) given by Eq. 62 If L(z... L(I 1n41) manifold, while L L(x,,, Da) L(x Dn+1) implics at least one of the points P, or Pa+1 lies exactly the stable manifold. In both cases there exists homoclinie intersection.
0, then P., and Pa+1
opposite sides of the stable
27

---

if Published as a conference paper at ICLR 2026
Algorithm 4 Investigating Homoclinic Intersections in 2D
procedure LMYESTOAACHOOMOLLNNCEEEESSSSSSSS
2:
Oc = (2))
Admissible saddle fixed point
3:
(Ag,2a) t (Stable Eigenvalue, Unstable Eigenvalue)
4:
(L,, + (Stable Eigenline, Unstable Eigenline)
5:
Po 6 -0, E
> Unstable eigenline hits the border at Po > Image of Po on the right side of the border > First fold point of the unstable manifold Minimum iterations for Po return to the left side Compute P, directly from Po Check for homoclinic intersection
6:
Pi < (cyo + h1, dyo + l2)
7:
Pi E e=
8:
4 Border return time
9:
PR + o + (AR 1) (AA h
10:
Pa+1 Next iteration
11:
if ( Vn) -L(nn+1 Bn+1) < then
Return "Homoclinic intersection exists"
13: else
Return "No homoclinic intersection"
15:
end
end procedure
H.3
SIMULATED BENCHMARKS
Damped oscillator
The 2M -M--dimessional damped oscillator, with æ, y E RM,, defined by
f Vi
(53)
Ú: 0
B
7i: i= M,
(54)
For the illustration in Fig. 2 we used M = 10 with parameters
a [1.00, 0.80, 0.90, 0.30, 0.30, 0.25, 0.50, 0.89, 0.73, 0.21]
B | 0.02. -0.10, 0.01, -0.01, -0.011, 0.01, -0.06, -0.052, 0.03, -0.012),
7 [0.10, 0.05, 0.12, 0.07, 0.00, 0.03, 0.009, 0.04, 0.06, 0.0008).
Duffing system The Duffing system is a classical nonlinear DS that models a damped and peri- odically driven oscillator with a nonlinear restoring force. Originally introduced by Georg Duffing (Dufling 1918). the system is described by second-order differential equation of the form
+ 8á + ocr + Axx = 7 cos(wil)
where x represents the displacement, a damping coefficient, a and 3 determine the linear and nonlinear stiffness, respectively, and cos(wt) is external periodic forcing term. The Duffing system exhibits rich dynamical behavior, including periodic, quasi-periodic, and chaotic dynamics (depending the forcing), making it prototypical example in the study of nonlinear dynamics. Here we used for the parameters a = 11 = 0.1, 5 = 0.5, and y = 0.
Simple models of decision making in the brain assume multistability between several choice-specific attractor states, to which the system's 's state is driven as one or the other choice materializes (Wang 2002). Here we employ a simple multistable model from Gerstner et al. (004). defined by
Multistable decision making model
dhgi -hgi we (egg)()))) WEl 7h RI1, Tg dt
(55)
dhga
Tø
-hg weegg(he)) WEl 7h RI2.
(56)
dt
dhinh dl
Tinh
hiih wir(gn(hm)) g
(57)
28

---

Published as a conference paper at ICLR 2026
with
gm
(58)
1 e-h-8))
(59)
and parameters:
TÆ 10, 1,
Tinh 5,
WEE 16, 12 = 40, 7 =
WEI -15,
WIE 12,
1 = 40,
5.
Lorenz63 The Lorenz63 system (Lorenz. 1963) is continuous-time dynamical system originally developed to model atmospheric convection. It describes the evolution of three state variables governed by the nonlinear differential equations
dz1 a ( - -xi);
di dr2 zi(p
xz3) 22:
dt
dz3 2122 ,333,
dt
where 21, .2, and denote, respectively, the convection rate, horizontal temperature difference, and vertical temperature difference. The parameters a, P, and correspond to physical constants related to the Prandtl number, Rayleigh number, and system geometry.
For specific values, e.g. a = 10, = 28, and 3 = 3. the system exhibits chaotic dynamics. These settings give rise to the well-known "butterfly attractor," prime example of deterministic chaos in low-dimensional systems.
H.4
EMPIRICAL DATA (SINGLE CELL RECORDINGS)
For Fig. 4D. vitro membrane potential recordings from a single cortical neuron were used (Hertäg et al. 2012). Many types of cortical cells exhibit bistability between spiking activity and
stable equilibrium either near the resting potential or a more depolarized level (Izhikevich) [0077] Durstewitz and Gabriel 2007 Dovzhenok and Kuznetsov 2012), and this example demonstrates how our algorithm can be utilized to reveal the structure the state space supporting this type of dynamics from real cells.
29

---

Published as a conference paper at ICLR 2026
1 ADDITIONAL RESULTS
I.1
ADDITIONAL FIGURES
Propagate region
Repeat process
0.035
0.035
Zoom
0.06
Zoom out
2/
0.050
0.050
out
0.04
0.045
0.045
0.040
0.040
0.085
0.085
0.09
Z
Z
2
Figure 8: Illustration the iterative procedure for computing stable manifolds with subregion boundaries of the shPLRNN (M 2, H = 10) model in purple (dashed) Step 1: The stable manifold (black) is initialized using the stable eigenvector at the saddle point (half-green), and sample points (orange) placed along it. Step 2: These points subregion, where the flow field is evaluated. Step 3: The updated manifold is given by the first principal component of points and the flow. Step 4: Repeating this process iteratively reconstructs the full global structure of the stable manifold (black), overlaid with the underlying GT flow field (blue)
propagated until they reach new linear
A
B
1040
10
Runtime (seconds)
20
Our Nethod
Runtime (seconds)
Our Method
Runtime (secondis)
Method
10a
10
10¹
16
16
20
40
22
2
2
29
2'2
Dimension (M)
number of linear subregions
number of linear subregions
Figure 9: A) Algorithm runtime for determining the stable manifold of saddle point in the Duffing
system averaged across runs (error bars STDV). For a constant number of linear subregions 2", Algo. s runtime hardly increases a function of model size M for ALRNN (P 2), confirming it is not significantly affected by the number of model parameters per se (note logarithmic scaling on y-axis). Algorithm was run with constant N, 1000 and Niter 2 4. Parameter N, was chosen deliberately high (much higher than usually required) to simulate the worst-case scenario. Runtime was determined on Intel Core i5- 124P. B) For constant model size M, Algo. Its runtime was hardly affected by the number 2 of linear subregions for an ALRNN (M = 40) when the manifold construction is restricted to the of linear subregions explored by the data the Duffing system. Algorithm I was run with N, = 1000 and Nater = 2". Runtime was determined on an Intel Core 1240P. C) Algorithm runtime for determining the stable manifold of saddle point averaged across 3 runs (error bars = STDV) for the Lorenz63. For a consiant model size M, Algo. Its runtime was hardly affected by the number 2" of lincar subregions for an ALRNN (M 100) when the manifold construction is restricted the set of linear subregions explored by the data. The manifold was computed with N, = 5000 and Niter = 2""; as in A-B, N, was chosen deliberately high to cover the worst case. All models were trained on the Lorenz system. Runtime was determined on an Intel Xeon Gold 6248. We emphasize that in general the scaling may strongly depend on the system's 's actual dynamics and topological structure, such that general statements regarding scaling are therefore difficult.
30

---

Published as a conference paper at ICLR 2026
A
C
spike time mean((Si)
200
907
10-
1.25
2g
10
200
400
600
800
Time
Spike time
21
Figure 10: ALRNN (M = 80, P 4) trained on 5-dimensional embedding of a human ECG recording. A) Example of ALRNN--imulated ECG activity (black) and true human ECG (gray). B) Exponential growth in differences in spike timing between ALRNN-generater trajectories from two slightly different initial conditions as a function of spike time (ISI interspike interval). While the divergence is exponential (note the logarithmic scale of the y-axis), thus suggesting the presence of chaos, the divergence rate appears to be rather small. This impression is corroborated by the estimated slightly positive maximum Lyapunov exponent of 0.0074 + 0.00083 (mean SD), while the sum all Lyapunov exponents < 0, as estimated from the model Jacobians evolved for 10,000 time steps for 5 different initial conditions. It thus appears the system is weakly chaotic (Alligood et al. 1996 Guuckenheimer and Holmes 2013). but is difficult to conclusively establish this from these observations alone. C) Visualization of the presumably chaotic attractor for the real data (gray, via delay-embedding:: Kantz and Schreiber (2004)) and as generated by the ALRNN (black) in 3d subspace of the 80-dimensional state space. heteroclinic intersection (within the limits of numerical precision, = 10 10-10 10-10) of the saddle nodes in green was identified by our algorithm at the red dot, thus confirming the presence chaos in this system.
5.5
5.0
4.5
4.0
&
3.5
3.0
2.5
2.0
-4
-3
-2
X1
Figure 11: Illustration of how curvature can arise in the un-/stable manifolds as borders other linear subregions are crossed. The stable manifold (black) of the saddle Chalf-green is straight line in the first linear subregion (left hand side), but becomes curve as enters the second subregion (right from dashed purple border), where the dynamics follows a spiral. The parameters of the 2d PLRNN used for illustration were: = [0.93 0.0; 0.0 0.92]. W [0.26 0.080; -0.21 0.24),
h [0.43, -0.57].
31

---

Published as a conference paper at ICLR 2026
1.2 ANALYTICAL CONDITION FOR HOMOCLINIC INTERSECTION OF THE FIRST AND SECOND FOLD POINTS
Homo- or heteroclinic intersections between stable and unstable manifolds imply the existence of infinitely many such intersections, leading to a "horseshoe structure' and the existence of chaotic orbits (Wiggins 1988). In this subsection we derive the analytical results needed to confirm such intersections using Algorithm 1
Let 02 = (a2 y be an admissible saddle fixed point in the left sub-region Since O2 is saddle, Ag has one stable and one unstable eigenvalue Å and X respectively. Let us denote the line generated by the associated stable eigenvector by E and the line produced by the corresponding unstable eigenvector by e Since the unstable eigenvector va (2) (Ax," 1)" ey hits
the border z 0 at Po = (0. yo) E Σ where
yo = v2 22 1/
01
= by h1 + (1-a,) h2
(1 d) + cc2
bj
Au -d)
1 1ctDe
1- l'c De
(b. 11 + (1 (1-1)/2)) 41)12) (Au d) - -bi(ch2 + (1 dd/12)
(60)
(Au - d)(1 - +
The image of Po is the first fold point of the unstable manifold of O2, and all its images will also be fold points. Its coordinate is P1 (c3o h yo h22) The image of the first fold point is the second fold point P2 = (I2, 82) with coordinates
Z2 c(au d)yo (ar + 1)h1 + ch2
if cyo + h1 < 0
1/2 (bic + ) + byh1 + (d + 1)h2
12 = c(a, + d)yo + (a, + 1)h1 + ch2
if + h1 >0
(61)
1/2 = (b,c + )yo + b, + (d + 1)h2
Now we check whether or not the points P1 and P2 are on opposite sides of the stable eigenline E When P1 and P2 are opposite sides of E. then the unstable manifold must have intersected the stable manifold. Thus, we have a homoclinic intersection which implics the occurrence of chaotic dynamies. Since the stable eigenvector v, ( 1)TT
and the eigenline e passes through
O = VT) Ey can be computed as
X - d bi
(1 d) h1 + che
C-
r Dr
\₂ (~,-
-d)
d
by h1 + (1 - a1) h2
L(r,y) = 0.
(62)
b
1 L + De
Now there are two possibilities:
Case I: I:L( v1) L(x2, V2) < 0
If L(z1 v1) L(x2, w2) 0, then P1 and P2 are on opposite sides of the stable manifold, while L(I1, L(r2, 32) 0 implies that at least one of the points P1 and P2 lies exactly on the stable manifold. In both cases there exists homoclinic intersection, i.e., whenever we have
A~ (dyo + h2) - (cyo h1) + .M)
A d
^ d ((by),, dl d)
+ bi/rt1 + (d + 1)hz) - (o(aı, + d)yo + (ay, + 1)h1 + ch2)
+ +M) s 0,
(63)
32

---

Published as a conference paper at ICLR 2026
where M
(1-4) By bo 1 1-Tc+Dc
-
by I-T 1-rc+Dc
The homoclinic intersection point Phom =
(Thom: yhom) ) is the point of intersection between the line joining the two fold points and the stable eigenline Ey and hence given by
hhom = (222 - 11) + x1
(64)
Uhom
(y2 v1)33 81
where
A-d - M
21
8
(65)
(/22
1/)
(x2 r1)
bj
Finally, we must ensure that the intersection happens before the stable eigenline hits the border. For this, need have Thome x2 0, which implies the intersection point Phom and the fixed point O are on the same side of Σ.
Case II: (1) L(x2, 12)
L(x1, 11) L(x2, 12) > O, we need to check whether the unstable manifold intersects with the part of the stable manifold which ensues after folding along the y-axis. For this have to calculate more points of the global stable manifold. Since 7 is assumed to be invertible, the global stable manifold formed by the union of all preimages (inverses) of any rank of the local stable set (a segment of the local stable eigenline). Assume, under the action of T the line e maps to the y-axis and intersects it at Po = (0, Jo)T E E. Then Po is the first fold point of the stable manifold of O% and jo is given
by
jo = (be h + (1 - a1) (A. - d) -
bi 12)
+ (1 d)m1
(66)
(), d)(1 D D2)
The preimage of Po is the second fold point P_ = (#-1, 1-1-1) with coordinates
$_ D - cijo ch2 - -dh,)
4 (Po h)
if
50
Ú- = (ar00 + b111 a111)
De
Φ_1 DR cio + ch2 dhi)
lo
(Po
- h))
if
2 0
(67)
5- = Dr Da(a; (a, + b a,112)
DR
where
T (Po h) -dhn c
(Jo h2).
(68)
DC/RR DeiR DC/R
Now the line joining the two fold points Po and P_1 is given by L(x, v) 0 where
[(æ,v) = 30 7 Do
r.
(69)
L(x1, 91) L(x2, 32) < O, then P1 and P2 are on opposite sides of the stable manifold, and the unstable manifold intersects the stable manifold at Phom = hohm (hom) J with
Thom (x2 r2) 51
(70)
Thom (322 v1)33 81
where
3 (wo 1)) ( 3o0 20
(71)
(3/ 3/ v1)ã-1 - (i-1 - 30) x1)
We need to have Thome & 0, which means the intersection point Phom and the point P_, on the same side of Σ. Remark L.I. An analogous procedure can be performed to analytically obtain homoclinic intersections for the fixed point O E R.
33

---

Published as a conference paper at ICLR 2026
1.2.1 HOMOCLINIC INTERSECTIONS FOR FURTHER ITERATIONS OF FOLD POINTS
In order determine homoclinic intersections, may be necessary check further iterations of fold points for which we would need recursive procedure. This will involve calculating the n-th power of 2 matrix, and we therefore first prove the following Proposition:
Proposition 1.2. Let M
a
a) have two distinct eigenvalues
b
+ d Vaa -oa-dd+++ + 4bc
VT
4D
2L2
F
(72)
2
2
2
where I' and D are the trace and determinant of M. Then, for every E N
(Anti - An
cAn
M
(73)
bAn
dAp D DAn An
where
A, Ay N
(74)
t1 Av
Proof. In our case M is diagonalizable, so for b + 0
M = VAV
A1 A1-d
A2-d"
A1
0
b
d - A2 A1 d
(75)
d1 A2
1
1
0
A2
b
Therefore
M = V A" V-1
A--d
d=-d b
0
(
b
d- A2
\Y A2
1
0
xm
A -d d
x
A1 x2
(g - A2) (7 - A1)
(d= (2F X4)
b(1 >3)
Xi X2) 2z A1)
-(Ant3
AAat - d An
d d-rd+D An
bAn
d dA, - D DAn-1/
-(Art1
(An+1 - d An
cAn
(76)
bAn
d A - D A
If b = 0, then the eigenvalues of M are X and A2 d (a + d). Thus,
a 0
M VAV-1
(77)
0
0
and so (for n E N)
(
=
6
-
M = VA-V-L
0
d
c(a" -d ") d
(78)
d
which yields equation Eq. 73 for b 0.
34


Reviews:
[{"review_id": "JhFePs29Iu", "reviewer": "Reviewer_KGgx", "summary": "The paper proposes an algorithm for detecting stable and unstable manifolds of fixed and cyclic points in piecewise-linear RNNs with ReLU activation functions. The authors claim this is the first algorithm for detecting stable/unstable manifolds in ReLU-based RNNs. The paper then utilizes some dynamical system examples such as Lorenz-63 attractors and single-cell recordings to illustrate the proposed algorithm.", "questions": "See Weakness", "limitations": "", "rating": 8, "confidence": 4, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "- Clarity: detailed math background and derivations, and pseudocodes facilitate the understanding and reproducibility of the proposed approach.\n\n- Mathematically rigorous. \n\n- Validation on across synthetic and real-world problems. And the algorithm handles high-dimensional models with sublinear scaling in subregion sampling. Its efficiency enables previously infeasible analyses, such as identifying attractor basins in decision-making tasks or chaos signatures in electrophysiological data, bridging theoretical insights with empirical applications.", "weaknesses": "To begin with, I'd like to disclose that I also reviewed this paper in NeurIPS 2025 and I gave quite positive score. Below I just attach my weaknesses when I reviewed this paper in NeurIPS 2025, and the authors' rebuttal in NeurIPS 2025 have addressed nearly all my concerns.\n\n1. I would have liked a bit more concrete discussion on which parts of the algorithm are heuristic and which are exact. And when they are heuristic, how do potential hyper parameter choices affect the outcome, and how are they picked? I will give some concrete examples here:\n\n- Use of approximate algorithm for finding the fixed points: Is there any concern in case one doesn't find all fixed points?\n- Sampling seed points: how many are sampled? How are they sampled? It doesn’t seem trivial (or even always possible) to me to e.g., uniformly sample within one region of linear dynamics.\n- Propagating seed points. I think one is not always sure we reach all reachable manifolds from a given sub-region? (E.g., if we undersample seed points, or all seed points / most of the space ends up converging in one particular direction).\n- The fallback algorithm 3 still doesn’t have any guarantees of finding the right subspace, does it?\n\n2. From results shown it is unclear to me that the regularisation (Eq. 6) actually reaches the desired effect. Besides the effect on training, does it actually lead to networks that are more likely to be invertible? Is there some way to verify this, maybe by comparing the average condition numbers of the to be inverted matrices when running your manifold construction algorithm on a model trained with and without the regulariser?", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "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": "", "extra_scores": {}}, {"review_id": "KIjZ8e9z52", "reviewer": "Reviewer_rEg4", "summary": "The authors propose an algorithm to detect stable and unstable manifolds in piecewise linear RNNs. The main challenge is moving across regions, and the algorithm handles this by testing for self-consistency across the boundary, and rejecting inconsistent continuations. The authors use the method to locate boundaries between basins of attraction, and also use it to analyze a chaotic system.", "questions": "Line 133: question mark\n\nThe stable manifold is a high-D object in general. Figure 1 is a 1D case.\n\nLine 338: should this be figure 2?\n\nLine 365 – this is a relatively high-dimensional example. But the validation of the reconstruction is only done visually in a 3D projection. It would be useful to quantify the quality of the stable manifold. For instance, running the dynamics forward from nearby points.\n\nScalability. How does the algorithm scale with system dimensionality? In general, validation on high dimensional examples is not clear.\n\nFigure 5: How does the figure show that the intersections specifically are the chaotic attractor? Figure 5B shows points on all the unstable manifold – not just the intersections.", "limitations": "", "rating": 4, "confidence": 4, "soundness": 3, "presentation": 2, "contribution": 3, "strengths": "Identifying stable and unstable manifolds is an important problem. Most RNN analysis methods deal with fixed points or trajectories.\n\nLimiting the analysis to PLRNN provides a feasible algorithm.", "weaknesses": "Discrete vs. continuous time: The definitions in equation 4 (for instance) are of discrete time, where stability is determined by the norm of eigenvalues. But section 3.3 (and later) describe continuous time conditions – negative/positive real part.\n\nRNNs are usually high-dimensional, as mentioned by the authors. Most examples in the paper are low dimensional. And the high-dimensional examples are only visualized in low-D as a check of their validity.\n\nMany results are qualitative, and only visual inspection is provided as a measure of correctness.\n\nThe presentation of the paper is not very clear.  For instance, the argument about curved surfaces. If there is a complex eigenvalue pair, it is associated with two eigenvectors spanning a hyperplane. The description in Section 3.3 is of sampling points on this hyperplane. It does not describe following one specific trajectory (which is usually curved also for real eigenvalues).", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "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": "", "extra_scores": {}}, {"review_id": "87f6iXBs9g", "reviewer": "Reviewer_MbXp", "summary": "This paper discusses the problem of finding stable and unstable manifolds in a class of recurrent neural networks. The authors treat RNN's involving the ReLU activation function which allows them to treat the system as a piecewise linear system. From what I understand of the paper, the authors exploit the piecewise linear property of the system and use the eigendecomposition of the Jacobian as the basis for computing the manifolds. They have performed some tests on various systems, some simple, some more complex and their algorithm appears to work well.", "questions": "The examples look good, but to me it was not clear exactly what the algorithm returned when computing manifolds - was it just samples from the manifold and then one would deduce from these samples the likely shape of the manifold? This would of course be dependent on the density of samples used, particularly in critical regions of the state-space. Some more explanation and guidance would be useful.", "limitations": "", "rating": 4, "confidence": 3, "soundness": 2, "presentation": 3, "contribution": 2, "strengths": "I found the paper interesting and relevant. I have read numerous papers on computing stable/unstable manifolds in nonlinear systems, but these have typically been in continuous time and I think this paper is one of the few (or perhaps only) paper which treats this specific type of system. It should be observed that generally the computation of stable/unstable manifolds for nonlinear systems is a very hard problem and one which gets harder as the order of the system increases and the dynamics become more complex. \n\nThe writing in the paper was good in terms of well constructed prose and a well-presented manuscript. \n\nThe authors' algorithm appeared to work well on the examples tested.", "weaknesses": "My chief complaint about the paper is it meandered awkwardly between giving the reader quite basic details (and perhaps too much of these) in some places and then depriving the reader of significant details of the authors' contribution. In some places it was not clear what was obvious, what was the authors' work and what was new. The inclusion of a few more references at particular places would be useful. The paper is quite easy to understand up until Section 3.3 when I felt the authors needed to give significantly more detail. In particular I thought:-\n\n1. The authors should give more detail about Algorithm 1, which does not appear to have been described comprehensively. Indeed, it does not seem to referenced in the main text until Section 4. I would like to see a better, more comprehensive explanation of the ideas behind it, its construction, and a more detailed description of its steps and the various mathematical objects it contains. For example in Algorithm 1 $\\sigma$ is an index, but what values does it take (one can guess? One expects $Q$ to be the manifold, but in the algorithm it seems like $Q^{\\sigma}_k$ is a selection of points in that manifold? I am not sure. Also in Step 4 it seems like $n$ denotes iteration number, but is this actually $k$? These issues are probably resolvable, but without some sort of description, make the paper only possible to appreciate at a superficial level. \n\n2. I was not entirely sure how the SCYFI package was used. Does this provide the points $P$ for Algorithm 1? Please explain and perhaps also indicate how $N_{max}$ is chosen. \n\n3. In Section 3.4 the authors discuss regularization which certainly makes sense. However, the details of how this regularizing term is included is not clear: to which \"loss\" is the term (6) included? In the discussion below equation (6), the authors refer to Figure 3 (top), (bottom) and (right), but this does not seem to correspond with the figures included in the paper - there is only one row of figures in Figure 3. I was confused.\n\n4. I found the appendices quite long and strangely structured. It was not quite clear to me how central some of these were to the paper. For example the first few paragraphs of Appendix D seemed to be only weakly relevant to the paper. Section D.1 seemed relevant, but then Section D.2 seemed rather verbose and un-focused. The relevance of the appendices to the main text needs more focus. Finally, it seems that some sort of \"training\" is needed for the algorithm reported, according to the appendix, but this does not seem to be discussed in the main text. Some clarity here would be useful for the reader. \n\n5. Finally a comment about the literature. The system which the authors consider belong to a type of so-called \"Lur'e\" or \"Lurie\" or \"Lurye\" systems considered in the control theory literature. Furthermore the PLL RNN which they consider seems quite close to the piecewise linear systems again considered in the control theory literture. I am not asking the authors to cite any specific papers, but they may find it useful to consult some of this literature.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "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": "", "extra_scores": {}}, {"review_id": "M7u8yAdjv2", "reviewer": "Reviewer_KgUd", "summary": "This paper presents a new algorithm for identifying stable and unstable manifolds in piecewise-linear recurrent neural networks (PLRNNs). The authors leverage the linearity within each ReLU activation region to derive an efficient, semi-analytical method that constructs manifolds by iteratively propagating sampled points across subregions. This approach enables explicit detection of manifolds that delineate basins of attraction and reveal multistability or chaotic structures in RNN dynamics. The authors demonstrate the method on several systems including the Duffing oscillator, Lorenz-63 attractor, decision-making RNNs, and real neuronal recordings. The results show close agreement between the algorithmic manifolds and those computed analytically or numerically from the underlying dynamical systems, suggesting both accuracy and scalability beyond traditional continuation techniques.", "questions": "Can the authors please clarify the relationship of PLRNNs to TLNs, and situate their findings within this literature?", "limitations": "", "rating": 6, "confidence": 4, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "The paper is clearly written and well-organized. The proposed algorithm effectively leverages the piecewise-linear structure of PLRNNs, leading to an elegant and computationally efficient solution. The authors successfully ground their work in established theoretical principles from dynamical systems theory and neuroscience, making a strong case for its relevance to both machine learning and scientific modeling. The experimental section is broad and persuasive, spanning toy systems, classical chaotic dynamics, and real neural recordings, which together illustrate the versatility of the method and its capacity to illuminate the internal mechanisms of learned recurrent dynamics. The addition of invertibility regularization and supporting theoretical analysis further reflects the authors’ careful consideration of the method’s stability and practical applicability.", "weaknesses": "As far as I can tell, the paper does not engage with the existing literature on Threshold-Linear Networks (TLNs), particularly the extensive work by Carina Curto and collaborators. This omission is notable, as TLNs share strong conceptual and mathematical parallels with the proposed framework. The TLN literature provides rigorous theoretical results on how ReLU-based recurrent networks partition state space into distinct regions, each exhibiting its own dynamical behavior—precisely the type of structure this paper seeks to characterize.  \n\nFor instance, see [this paper](https://arxiv.org/pdf/2109.03198) and [this paper](https://arxiv.org/pdf/1605.04463).", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "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": "", "extra_scores": {}}]

Author rebuttal:

Decision: Accept (Poster)