Bifurcations and loss jumps in RNN training
Authors: Lukas Eisenmann
Zahra Monfared
Niclas Alexander Göring
Daniel Durstewitz
Venue: NeurIPS 2023 spotlight (2023)
Keywords: dynamical systems
bifurcations
Recurrent Neural Networks
attractors
training algorithm
BPTT
exploding and vanishing gradient problem
nonlinear dynamics
time series
TL;DR: We show that bifurcations cause loss jumps in RNN training and introduce an efficient algorithm for locating fixed points, cycles, and bifurcation manifolds in RNNs

Abstract: Recurrent neural networks (RNNs) are popular machine learning tools for modeling and forecasting sequential data and for inferring dynamical systems (DS) from observed time series. Concepts from DS theory (DST) have variously been used to further our understanding of both, how trained RNNs solve complex tasks, and the training process itself. Bifurcations are particularly important phenomena in DS, including RNNs, that refer to topological (qualitative) changes in a system's dynamical behavior as one or more of its parameters are varied. Knowing the bifurcation structure of an RNN will thus allow to deduce many of its computational and dynamical properties, like its sensitivity to parameter variations or its behavior during training. In particular, bifurcations may account for sudden loss jumps observed in RNN training that could severely impede the training process. Here we first mathematically prove for a particular class of ReLU-based RNNs that certain bifurcations are indeed associated with loss gradients tending toward infinity or zero. We then introduce a novel heuristic algorithm for detecting all fixed points and $k$-cycles in ReLU-based RNNs and their existence and stability regions, hence bifurcation manifolds in parameter space. In contrast to previous numerical algorithms for finding fixed points and common continuation methods, our algorithm provides $\textit{exact}$ results and returns fixed points and cycles up to high orders with surprisingly good scaling behavior. We exemplify the algorithm on the analysis of the training process of RNNs, and find that the recently introduced technique of generalized teacher forcing completely avoids certain types of bifurcations in training. Thus, besides facilitating the DST analysis of trained RNNs, our algorithm provides a powerful instrument for analyzing the training process itself.

Bifurcations and loss jumps in RNN training
Lukas Eisenmann" 12", Zahra Monfared ',", Niclas Göring'. and Daniel Durstewitz 1,2,3
@zi-mannheim.de Department of Theoretical Neuroscience, Central Institute of Mental Health, Medical Faculty Mannheim, Heidelberg University, Mannheim, Germany Faculty of Physics and Astronomy, Heidelberg University, Heidelberg, Germany Interdisciplinnry Center for Scientific Computing, Heidelberg University "These authors contributed equally
[lukas.eisenn
Abstract
Recurrent neural networks (RNNs) are popular machine learning tools for modeling and forecasting sequential data and for inferring dynamical systems (DS) from observed time series. Concepts from DS theory (DST) have variously been used further our understanding of both, how trained RNNs solve complex tasks, and the training process itself. Bifurcations are particularly important phenomena in DS, including RNNs, that refer topological (qualitative) changes in a system's 's dynamical behavior as one or more of its parameters are varied. Knowing the bifurcation structure of an RNN will thus allow to deduce many of its computa- tional and dynamical properties, like its sensitivity to parameter variations or its behavior during training. In particular, bifurcations may account for sudden loss jumps observed in RNN training that could severely impede the training process. Here we first mathematically prove for a particular class of ReLU-based RNNs that certain bifurcations are indeed associated with loss gradients tending toward infinity or zero. We then introduce a novel heuristic algorithm for detecting all fixed points and k-cycles in ReLU-based RNNs and their existence and stability regions, hence bifurcation manifolds in parameter space. In contrast to previous numerical algorithms for finding fixed points and common continuation methods, our algorithm provides exact results and returns fixed points and cycles up to high orders with surprisingly good scaling behavior. We exemplify the algorithm the analysis of the training process of RNNs, and find that the recently introduced technique of generalized teacher forcing completely avoids certain types of bifurca- tions in training. Thus, besides facilitating the DST analysis of trained RNNs, our algorithm provides powerful instrument for analyzing the training process itself.
Introduction
Recurrent neural networks (RNNs) are common and powerful tools for learning sequential tasks or modeling and forecasting time series data BI 65 9 26 Typically, RNNs are "black boxes," whose inner workings are hard to dissect. Techniques from dynamical system theory (DST) can significantly aid in this effort, as RNNs are formally discrete-tima dynamical systems (DS) [37]. A better understanding of how RNNs solve their tasks is important for detecting failure modes and designing better architectures and training algorithms. In scientific machine learning, hand, RNNs are often employed for reconstructing unknown DS from sets of time series observations
the other
63 23 55 e.g. in climate or disease modeling. In this case, the RNN is supposed to provide an approximation to the flow of the observed system that reproduces all its dynamical properties, e.g., cyclic behavior in climate epidemiological systems. In such scientific or medical
10 30 64,
37th Conference on Neural Information Processing Systems (NcurIPS 2023).

---

settings, detailed understanding of the RNN's dynamics and its sensitivity to parameter variations is in fact often crucial.
Beyond understanding the dynamical properties of a once trained RNN, it may also be of interest to know how its dynamical repertoire changes with changes in its parameters. The parameter space of any DS is partitioned into regions (or sets) with topologically different dynamical behaviors by bifurcation curves (or, more generally, manifolds). Such bifurcations, qualitative changes in the system dynamics due to parameter variations, may not only be crucial in scientific applications where we use RNNs, for instance, to predict tipping points in climate systems or medical scenarios like sepsis detection [45]. They also pose severe challenges for the training process itself as qualitative changes in RNN dynamics may go hand with sudden jumps in the loss landscape [16] [44. Although methods from DST have significantly advanced the field in recent years, especially with regards to
algorithms for reconstructing nonlinear DS from data [55] 46 67, progress still hampered by the lack of efficient tools for analysing the dynamics of higher-dimensional RNNs and their bifurcations. In partieular, methods are needed for exactly locating geometrical objects like fixed points or cycles in an RNN's state space, but current numerical techniques do not scale well to higher-dimensional
scenarios and provide only approximate solutions 229 21 61].
The contributions of this work are threefold: After first providing an introduction into bifurcations of serrewise-linear (ReLU-based) RNNs (PLRNNs), which have been extensively used for reconstruct- ing DS from empirical data [10 30 we mathematically prove that certain bifurcations during the training process will indeed cause loss gradients diverge to infimity, resulting in abrupt jumps in the loss, while others will cause them to vanish. As RNNs are likely to undergo several bifurcations during training on their way from some initial parameter configuration toward a dynamics that successfully implements any given task, this poses severe challenges for RNN training and may be one of the reasons for exploding and vanishing gradients [16 44 25 We then create a novel, efficient heuristic algorithm for exactly locating all fixed points and k-cycles in PLRNNs, which can be used to delineate bifurcation manifolds in higher-dimensional systems. Our algorithm finds these dynamical objects in many orders of magnitude less time than an exhaustive search would take. Using this algorithm, we demonstrate empirically that steep cliffs in loss landscapes and bifurcation curves indeed tightly overlap, and that bifurcations in the system dynamics accompanied by sudden loss jumps. Finally, we prove and demonstrate that the recently introduced technique of generalized teacher forcing (GTF) [24] completely eliminates certain types of bifurcation in training, providing an explanation for its efficiency.
2 Related Work
DS analysis of trained RNNs In many areas of science one is interested in identifying the nonlinear
DS that explains a set of observed time series [12 23 A variety of purely data-driven machine
learning approaches have been developed for this purpose 11 [5], but mostly RNNs are used for
the goal of reconstructing DS from measured time series 110
55
30 64 63
46
45
In scientific
settings in particular, but also often in engineering applications, we seek a detailed understanding of the system dynamics captured- or the dynamical repertoire produced - by a trained RNN [62] 33 61 [4] 13). To analyze an RNN's dynamics, typically its fixed points are determined numerically, for instance by optimizing some cost function [21 34 61], by numerically traversing directional fibers [29). or by co-training a switching linear dynamical system [56). These techniques, however, often scale poorly with dimensionality, provide only approximate solutions, and are not designed for detecting other dynamical objects like k-cycles. This is, however, crucial for understanding the complete dynamical repertoire of trained RNN. PLRNNS are of particular interest in this context because their piecewise linear nature makes some of their properties analytically accessible [55 [0) which is tremendous advantage from the perspective of DST. In this work, we develop an efficient heuristic algorithm for locating PLRNN's fixed points exactly, as well as its k-cycles.
10
55
30 38
Bifurcations and loss jumps in RNN training The idea that bifurcations in RNN dynamics could impede the training process is not new [4]. Doya [16], to our knowledge, was the first to point out that even in simple single-unit RNNS with sigmoid activation function (saddle-node) bifurcations may occur as an RNN parameter is adapted during training. This may not only cause an abrupt jump in training loss, but could lead to situations where it is impossible, even in principle, to reach the training objective (the desired target output), as across the bifurcation point there is discrete change

---

Pascanu et al. [44] discussed similar associations between steep cliffs in RNN loss functions and bifurcations. Although profound for training success, this topic received surprisingly little attention over the years. Haschke and Steil [22] extended previous work by a more formal treatment of bifurcation boundaries in RNNs, and Marichal et al. [35] examined fold bifurcations in RNNs. The effects of bifurcations and their relation to exploding gradients in
in network behavior [6, 68 4
gated recurrent units (GRUs) was investigated in Kanai et al. [27]. Ribeiro et al. [50] looked at the
connection between dynamics and smoothness of the cost function, but failed to find link between bifurcations and jumps in performance. In contrast, Rehmer and Kroll [99 observed large gradients at bifurcation boundaries and concluded that bifurcations optimization. To the best of our knowledge, we are, however, the first to formally prove a direct link between bifurcations and the behavior of loss gradients, and to derive systematic and efficient algorithmic procedure for identifying bifurcation manifolds for a class of ReI.U-based RNNs.
indeed cause problems in gradient-baseds
3 Theoretical analysis
In this paper we will focus on PLRNNS as one representative of the wider class of ReLU-bbsed RNNs, but similar derivations and algorithmic procedures could be devised for any type of ReLU-based RNN (in fact, many other types of ReLU-basee RNNs could be brought into the same functional form as PLRNNs; e.g. [TO). We will first, in sect. 3.1 provide some theoretical background on bifurcations in PLRNNs, and illustrate how existence and stability regions of fixed points and cycles could be analytically computed for low dimensional (2d) PLRNNS. In sect. 3.2 we will then state two theorems regarding the association of bifurcations and loss gradients in training. It turns out that for certain types of bifurcations exploding or vanishing gradients are inevitable in gradient-based training procedures like Back-Propagation Through Time (BPTT).
3.1 Bifurcation curves in PLRNN parameter space
The PLRNN, originally introduced as a kind of discrete time neural population model [8] has the general form
2, = Fe(z:-1, = A2-1+ W ((z=1)) + Css: + h
(1)
where 2, E R. is the latent state vector and are system parameters consisting of diagonal matrix A E R MxM cuutorregression weights), off-diagonal matrix W E R MxM (coupling weights), ( max(2;-1, 0) is the element-wise rectified linear unit (ReLU) function, E RM constant bias term, and E K represents external inputs weighted by C' E Rx The original formulation of the PLRNN stochastic 18 30), with Gaussian noise term added to eq. (I. but
here we will consider the deterministic variant.
Formally, like other ReLU-based RNNs, PLRNNS constitute piecewise linear (PWL) maps, sub- class of piecewise smooth (PWS) discrete-time DS. Define Date) = diag diag(dor())) as a diagonal matrix with indicator vector dex()) == (d,, d2, ... ,drr such that dm (2m,t) =: dm = 1 whenever Zm,t 0, and zero otherwise. Then (T can be rewritten
21 = Fo(z = (A + W DD + h Wa(-1) 21-1 + h,
(2)
where we have ignored external inputs 8, for simplicity. There are in general 2 different con- figurations for matrix D and hence for matrix W Wo(e-1)- dividing the phase space into 2M sub-regions separated by switching manifolds (see Appx. A.I.I for more details).
Recall that tixed points of map z; 1)) are defined as the set of points for which we have 2 = Fo(2*), and that the type (node, saddle, spiral) and stability of a fixed point can be read off from the eigenvalues of the Jacobian Jr (=--1)
0z, evaluated at 2"
3 47 Similarly,
dz.-1
dz.-1
a k-cycle of map Fg is a periodic orbit {2] 22 ....... 2t} such that each of the periodic points
2 1.... kt is distinct, and is a solution to the equation 2 = (=()) i.c. the & times iterated
map Fg. Type and stability of a k-cycle are then determined via the Jacobian I[ J++k-r II 0zits-r_ @zctz-1, Solving these equations and computing the corresponding Jacobians 0zt-1
Ōžz=k-r-1
thus allows to determine all existence and stability regions of fixed points and cycles, where the latter a subset of the former, bounded by bifurcation curves (see Appx. A.I for more formal details).
3

---

To provide specific example, assume M = 2 and fix - for the purpose of this exposition parameters
W12 W22 = O, such that we have
(111
(111 + Wi
W = Wn
Wn = Wa =
(3)
G22
W21 0:22
L.e. only one border which divides the phase space into two distinct sub-regions (see Appx. A.I.I.
For this setup, Fig. IA provides examples of analytically determined stability regions for two low order cycles in the (d11, "11 + W11) pparameere plane (see Appx. A.Ip. Note that there are regions in parameter space where two or more stability regions overlap: In these regions we have multi- mulli--tabillty. the co-existence different attractor states in the PLRNN's state space.
As noted above, bifurcation curves delimit the different stability regions in parameter space and hence associated with abrupt changes in the topological structure of a system's 's state space. In general, there are many different types of bifurcations through which dynamical objects can come into existence, disappear, change stability (see, e.g., [3] 40 47), the most common ones being saddle node, transcritical, pitchfork, homoclinic, and Hopf bifurcations. In comparison with smooth systems, bifurcation theory of PWS (or PWL.) maps includes additional dynamical phenomena related the existence of borders in the phase space [6) Border-collision bifurcatiion (BCBs) arise when for a PWS map specific point of an invariant set collides with a border and this collision leads to qualitative change of dynamics [5 6 42). More specifically, a BCB occurs, if for PWS map ž₁ = Fe (z a fixed point or k-cycle either crosses the switching manifold Σ, := f2 E R" : c, 2 = 0} transversely at A = 6* and its qualitative behavior changes in the event, or if it collides the border with another fixed point or k-cycle and both objects disappear 41] Degenerate transcritical bifurcations (DTBs) occur when fixed point or periodic point of cycle tends to infinity and one of its eigenvalues tends to 1 by variation of a parameter. Specifically, let k, k 2 1, be a fixed point or a k-cycle with the periodic points {zt 22..... ................................ 26., and assume A denotes an eigenvalue of the Jacobian matrix at the periodic point 2;" i E (1,2, ,k). Then I' undergoes DTB at O 8", if (²*) - +1 and || + c0. F& undergoes a degenerate flip bifurcation (DFB), iff X² = -1 and the map F* has locally, in some neighborhood of 2;. infinitely many 2-cycles at 9 = 0". A center bifurcation (CB) occurs, if I& has pair of complex conjugate eigenvalues \1 and locally becomes center at the bifurcation value f = 0", i.e. if its eigenvalues are complex and lie on the unit circle ((A1,2(0*) Another important class of bifurcations are multiple attractor bifurcations (MABs), discussed in more detail in Appx. A.1 (see Fig. (S3). In addition to existence and stability regions of fixed points and cycles, in Appx. A.1.2 A.1.4] we also illustrate how to analytically determine the types of bifurcation curves bounding the existence and stability regions of k-cycles (DTB, DFB, CB and BCB curves; see also Fig. S66 for a 1d example).
1).
3.2 Bifurcations and loss jumps in training
Here we will prove that major types of bifurcations discussed above are always associated with exploding or vanishing gradients in PLRNNs during training, and hence often with abrupt jumps in the loss. For this we may assume any generic loss function [(O), like a negative log- lilliielioood or a mean-squared-error (MSE) loss, and a gradient-based training technique like BPTT 152] or -eweureent-Laaaeeeeeee
(RTRL) that involves a recursion (via chain rule) through loss terms across time. The first theorem establishes that DTB inevitably causes exploding gradients.
Theorem 1. Consider a PLRNN of the form (2) with parameters O = (A, W, h). Assume that it has stable fixed point or k-cycle & (k: 2 1) with Br, as its basin of attraction. If I undergoes a degenerate transcritical bifurcation (DTB) for some parameter value = Gg E 0, then the norm of the PLRNN loss grudient, 2C. tends to infinity at 0 = 0o for every 21 E Br i.e.
00
lime ||
8C, CC.
Proof. See Appx. A.2.I
However, bifurcations may also cause gradients to suddenly vanish, as it is the case for BCB as established by our second theorem:
Theorem 2. Consider PLRNN of the form
(2) with parameters 0 (A, W. h). Assume that has a stable fixed point or k-cycle F (k 2 1) with Br, as its basin of attraction. If undergoes a border collision bifurcation (BCB) for some parameter value 0 = 0, E 9, then the gradient of the loss function, ac, vanishes at 0 0g for every 21 E Br i.e. limg 02, 0.
50
4

---

A
0
DFB+BCD
-1
SRC
-2
+-00
-2
3CB
33
01
-3
Sxc
4
-1.5
-0.5
0
0.5
1
1.5
-1.5
-0.5
Q11
1.5
C
DTB
BCB
DFB
DFB
D
2.0
1.5
1.0
0.5
-2
0.0
-0.5
-1.5 -1.0
0.0
0.5
1.0
Q1
21
Figure A) Analytically calculated stability regions for 2-cycle (SRe, red), a 3-cycle (Sx C2, blue), and their intersection (yellow) in the (011, 411 + wi1 -parameter plane for the system eq. (3) with 0222 0.2, W21 0.5. B) Same determined by SCYFI, with bifurcation curves bordering the stability regions labeled by the type of bifurcation (DTB Degenerate Transcritical Bifurcation, BCB= Border Collision Bifurcation, DFB Degenerate Flip Bifurcation). C) Bifurcation graph (showing the stable cyclic points in the 21 coordinate) along the cross-section in A indicated by the gray line, illustrating the different types of bifurcation encountered when moving in and out of the various stability regions in A. D) State space at the point denoted 'D' in A (for
0.253, an + W11 -2833), where the 2-cycle (red) and 3-cycle (blue) co-exist for the same parameter settings, with their corresponding basins of attraction indicated by lighter colors.
Proof. See Appx. A.2.2
Corollary 1. Assume that the PLRNN (2) has stable fixed point F1 with Br, as its basin of attraction. I undergoes a degenerate flip bifurcation (DFB) for some parameter value Og E then this will always coincide with a BCB of 2-cycle, and result limg 21 E Br,
8C2 for every
de
Proof. See Appx. A.2.3
Hence, certain bifurcations will inevitably cause gradients suddenly explode or vanish, and often induce abrupt jumps in the loss (see Appx. AA.3.2 for when this will happen for a BCB). We emphasize that these results imputs. Since inputs do not affect the Jacobians in eqn. (70). (7I) and (76). they do not change the theorems (even if they would affect the Jacobians, Theorem 2 would be unaltered, and Theorem coulld be amended in straightforward way). Furthermore, since we are addressing bifurcations that occur during model training, from this angle inputs may simply treated as either additional parameters (if piecewise constant) or states of the system (without changing any of the mathematical derivations). fact. mathematically, any non- non-autonomous dynamical system (RNN with inputs) can always and strictly be reformulated as an autonomous system (RNN without inputs), see B 47 66].
general and hold for systems of any dimension,
well in the presence of
4 Heuristic algorithm for finding PLRNN bifurcation manifolds
4.1 Searcher for fixed points and cycles (SCYFI): motivation and validation
sect. 3.1 and Appx. A.1.2 A.1.4 we derived existence and stability regions for fixed points and low order (I: 3) cycles in 2d PLRNNS with specific parameter constraints analytically. For higher-order cycles and higher-dimensional PLRNNS (or any other ReLU-type RNN) this is no longer feasible due to the combinatorial explosion in the number of subregions that need to be considered M and increase. Here we therefore introduce an efficient search algorithm for finding all k--cccles of given PLRNN, which we call Searcher for Cycles and Fixed points: SCYFI (Algorithm []. Once all k-cycles (k > 1) have been detected on some parameter grid, the stability regions of these objects and thereby the bifurcation manifolds can be determined. k--cycles were defined in sect. 3.1
5

---

Algorithm SCYFI
The algorithm is iteratively run with k = 1. 1... Kmaws with Kmar the max. order of cycles tested Input: PLRNN parameters A, W, h;
L = (Ln ) collection of all sets Ln = {(2 (m)) } M= M of all lower order n-cycles
discovered so far, where Mn is the number of found cycles {z( (m))
l==1 of order n, with
corresponding ReLU derivative matrices {D( (D( in 11=1
Parameters:
Nout: max. number of random initialisations;
Nen: max. number of iterations
Output: LU Lx: Lx: set of all discovered k-cycles
1: Lx =
2: +
3: while i < Nout..... do
4: Select subregions Dinit = (D1, D2,..., D at random with replacement
5: c 0
6:
while c < Nin do
7:
Solve eq. for cycle candidate {z/ [/ with Wi(k-rr defined in eq. 2 based
on D Determine
8:
(Di )k based on the signs of the corresponding components of { Jim
9:
if D
[D (self-consistency) and v1 S S k, ((z (1m) (- E f
li=1
{ (m) k- g { then
1=1
10:
LA + La U ((
11:
12:
else
13:
Dinit (Di 1=
14:
end
15:
16:
end while
17:
8
18:
end while
and for the PLRNN, eq. are given by the set of k-periodic points (2i z where
2 Wss(k-))
ΣI Work-.) + h
(4)
if (1-II-- = Wo(k-r)) is invertible (if not, we are dealing with a bifurcation or a continuous set 1,... Ak- - 1, with corresponding matrices Wn(n A+WD; Now, if the diagonal entries D; are consistent with the signs of the corresponding states =,m.. i.e. if d 1 if a > 0 and d is a true cycle of eq. (2). otherwise we call it virtual. To find a k-cycle, since an M M-dimensional PLRNN harbors M different linear sub-regions, there approximately 2 different combinations of configurations of the matrices DI, I = 1 k. to consider (strictly, mathematical constraints rule out some of these possibilities. e.g. not all periodic points can lie within the same sub-region/orthant)
of fixed points). The other periodic points are = F 1
otherwise for all {2), ..... =1]
Clearly, for higher-dimensional PLRNNs and higher cycle orders exhaustively searching this space becomes unfeasible. Instead, found that the following heuristic works surprisingly well: First, for some order and a random initialization of the matrices D,, I 1......... 1...k, generate a cycle candidate by solving eq. 4.. If each of the points 2/ I 1. k. is consistent with the diagonal entries in the corresponding matrix 1 1... k and none of them is already the current library of cyclic points, then true k-cycle has been identified, otherwise the cycle virtual (or super-set of lower-order cycles). We discovered that the search becomes extremely efficient, without the need to exhaustively consider all configurations, if new search loop is re-initialized at the last visited
virtual cyclic point (i.e., all inconsistent entries ( in the matrices D,, I = 1... k are flipped, (1) 1 (0) to bring them into agreement with the signs of the solution points, z; [2 of eq.

---

4 thus yielding the next initial configuration). It is straightforward to see that this procedure almost surely converges if ouu is chosen large enough, see Appx. [A.244] The whole procedure is formalized in Algorithm | and the code is available at https::/githhu. com Cetomm/Dustwwwtee.. ttiar//SC SCYFT
To validate the algorithm, can compare analytical solutions as derived in sect. 3.1 to the output of the algorithm. To delineate all existence and stability regions, the algorithm searches for all k-cycles up to some maximum order K along a fine grid across the (a11. a1 - w1 1 parameeter plane. A bifurcation happens whenever between two grid points cycle appears, disappears, or changes stability (as determined from the eigenvalue spectrum of the respective k**-order Jacobian). The results of this procedure are shown in Fig. [B, illustrating that the analytical solutions for existence and stability regions precisely overlap with those identified by Algorithm [
(see also Fig. S7).
4.2 Numerical and theoretical results on SCYFI's scaling behavior
Because of the combinatorial nature of the problem, it is generally not feasible to obtain ground truth settings in higher dimensions for SCYFI to compare to. To nevertheless assess its scaling behavior, we therefore studied two specific scenarios. For an exhaustive search, the expected and median numbers of linear subregions n until object of interest (fixed point or cycle) is found, i.e. the number of { constellations that need to be inspected until the first hit, are given by
2 +1 N ++1 m + m +11 πT =imin
(2Ms n
1
22 m
E|n]
min n E N
(5)
m
with m being the number of existing k-cycles and N the total number of combinations, as shown in Ahlgren 1 (assuming no prior knowledge about the mathematical limitations when drawing regions). The median Ti as well as the actual median number of to-be-searched combinations required by SCYFI find at least one k-cycle is given for low-dimensional systems in Fig. 2A as a function of cycle order k. and can be seen to be surprisingly linear as confirmed by linear regression fits to the data (see Fig. legend for details). To assess scaling as function of dimensionality M, we explicitly constructed systems with one known fixed point (see Appx. A for details) and determined the number n of subregions required by SCYFI to detect this embedded fixed point (Fig. 3B. ). general, the scaling depended on the system's 's eigenspectrum, but for reasonable scenarios was polynomial or
even sublinear (Fig. 2B, see also Fig. (S8. In either case, the number required SCYFI iterations scaled much more favorably than would be expected from an exhaustive search.
B 10'
lop-cg
linear subregions
15
linear aubeegoom
10a
mapton
quediati:
Ejs]
10
10a
Oowl -i 1.25
Dowl 2.33
5
10a
00
10
20 30
40
50
20
40
60
80
100
120
cycle order
Figure 2: A) Number of linear subregions searched until at least one cycle of order k was found by SCYFI (blue) vs. the median number ग7 an exhaustive search would take by randomly drawing combinations without replacement (black) as function of cycle order (M Each data point represents the median of 50 different initializations across 5 different PLRNN models. Error bars median absolute deviation. Linear regression fit using weighted least-squares (R² ~ 0.998. p < 10- 10-30). B) Number of linear subregions searched until specific fixed point was found as function of dimensionality M for different eigenvalue spectra (see Appx. A.3. for details).
fixed).
How could this surprisingly good scaling behavior be explained? As shown numerically in Fig. S9 when we initiate SCYFI in different randomly selected linear subregions. it converges to the sub- regions including the dynamical objects of interest exponentially fast, offsetting the combinatorial explosion. A more specific and stronger theoretical result about SCYFI's convergence speed can be obtained under certain conditions on the parameters (which agrees nicely with the numerical results in Fig. 2). It rests on the observation that SCYFI is designed to move only among subregions containing virtual or actual fixed points or cycles, based on the fact that it always reinitialized with the next virtual fixed (cyclic) point in case the consistency check fails. The result can be stated as follows:
7

---

Theorem 3. Consider PLRNN of the form (2) with parameters 0 = {A. W.h). h). Under certain
conditions on (for which |||| [W] < 1). SCYFI will converge in at most linear time.
Proof. See Appx. A.2.5
Loss landscapes and bifurcation curves
5.1
Bifurcations and loss jumps in training
1.0
4.0
4.5
C 60
80
0.5
3.5
3.0 g
100
60
60
0.0
2.0
2.5
1050
44
1.5
-0.5
1
20
20
0.5
-1.0
-1.0
-.5
Arke
0.0
0.5
1.0
%
500 1000
Epoch
1500
%
500
Epoch
1000
1500
Figure 3: A) Logarithm of gradient norm of the loss in PLRNN parameter space, with ground truth parameters centered at (0, Superimposed in yellow are bifurcation curves computed by SCYFI, and in green two examples of training trajectories from different parameter initial conditions (indicated by the stars). Red dots indicate the bifurcation crossing time points shown in C & D. B) Relief plot of the loss landscape from A to highlight the differences in loss altitude associated with the bifurcations. C) Loss during the training run represented by the green trajectory labeled C in A and B. Red dot indicates the time point of bifurcation crossing corresponding to the red dot in A and B. D) Same for trajectory labeled D in A and B.
Fig. provides 2d toy example illustrating the tight association between the loss landscape and bifurcation curves, as determined through SCYFI, for a PLRNN trained by BPTT on reproducing specific 16-cycle. Fig. 3AA depicts contour plot of the gradient norms with overlaid bifurcation curves in yellow, while Fig. 3 shows the MSE loss landscape as a relief for better appreciation of the sharp changes in loss height associated with the bifurcation curves. Shown in green are two trajectories from two different parameter initializations traced out during PLRNN training in parameter space, where training was confined to only those two parameters given in the graphs (i.e.. all other PLRNN parameters were kept fixed during training for the purpose of this illustration). As confirmed in Fig. 3C & D, as soon as the training trajectory crosses the bifurcation curves in parameter space, huge jump in the loss associated with sudden increase in the gradient norm occurs. This illustrates empirically and graphically the theoretical results derived sect.
Next we illustrate the application of SCYFI on a real-world example, learning the behavior of a rodent spiking cortical neuron observed through time series measurements of its membrane potential (note that spiking is highly nonlinear behavior involving fast within-spike and much slower between-spike time scales). For this, we constructed a 6-dimensional delay embedding of the membrane voltage [53 (28). and trained a PLRNN with one hidden layer (cf. eq. 6 using BPTT with sparse teacher forcing (STF) [37] to approximate the dynamics of the spiking neuron (see Appx. A.3.2 for a similar analysis on biophysical neuron model). With M = 6 latent states and H = 20 hidden dimensions, the trained PLRNN comprises 20 different linear subregions and |0] 272 parameters, much than the toy example considered above. Fig. 4A gives the MAE loss as function of training epoch (i.e., single SGD updates), while Figs. 4 & Cillustrate the well-trained behavior in time (Fig. 4B) and in 2-dimensional projection of the model's state space obtained by PCA (Fig. 4C). The loss curve exhibits several steep jumps. Zooming into one of these regions (Fig. A indicated by the red box) and examining the transitions in parameter space using SCYFI, we find they are indeed produced by bifurcations, with an example given in Fig. 4D. Note that we are now dealing with high-dimensional state and parameter spaces, such that visualization of results becomes tricky. For the bifurcation diagram in Fig. 4D we therefore projected all extracted k-cycles (A 2 1) onto line given by the PCA-derived maximum eigenvalue component, and plotted this function of training epoch!" [ Since SCYFI extracts all k-cycles and their cigenvalue spectrum, we can also determine the type of bifurcation that caused the jump. While before the loss jump the PLRNN already produced time series quite similar to those of the physiologically recorded cell (Fig. 90,
'f course, very many of the PLRNN parameters may change from epoch to the next.

---

DTB (cf. Theorem 1 produced catastrophic forgetting of the learned behavior with the PLRNN's states suddenly diverging to minus infinity (Fig. 4F; Fig. STO provides an example of BCB during PLRNN training, and Fig. [STT an example of a DFB). This illustrates how SCYFI can be used analyze the training process with respect to bifurcation events also for high-dimensional real-world examples, well as the behavior of the trained model (Fig CC).
A
order
patential [mV]
-10
Cota
1.5
10sc
1.0
-30
4000
0.3
B/C
-50
-0.5
1.0x10
3.0x10'
5.0x10'
200
epach
PC2
D
E
potential [mv]
-10
10000111 [mv]
5.000 10
-20
-6
-30
-1.50x10
337.5
-00
22.00x10*
generatedd Sall
-42.5
-2.50x1n
-45.0
4800
200
50
100
time
Figure 4: A) Loss across training epochs for a PLRNN with one hidden layer trained on electro- physiological recordings from cortical neuron. Red box zooms in one of the training phases with huge loss jump, caused by a DTB. Letters refer to selected training epochs in other subpanels. B) Time series of true (gray) and PLRNN--immulated (black) membrane potential in the well trained regime (see A). C) All fixed points and cycles discovered by SCYFI for the well-trained model in state space projected onto the first two principle components using PCA. Filled circles represent stable and open circles unstable objects. The stable 39-cy 39-cycle corresponds to the spiking behavior. D) Bifurcation diagram of the PLRNN as a function of training epoch around the loss peak A. Locations of stable (filled circles) and unstable (open circles) objects projected onto the first principle component. E) Model behavior as in B shortly before the DTB and associated loss jump (from the epoch indicated in A, D). F) Model behavior as in B right around the DTB (diverging to -00).
5.2 Implications for designing training algorithms
What are potential take-homes of the results in sects. 3.2 & 5.1 for designing RNN training algorithms? One possibility is to design smart initialization or training procedures that aim to place or push an RNN into the right topological regime by taking big leaps in parameter space whenever the current regime is not fit for the data, rather than dwelling within a wrong regime for too long. These ideas discussed in bit more depth in Appx Appx A.3.3 with proof of concept in Fig. S12
1.00
Mww.daaama Mnmm My
Wh
0.25
0.00
150 200
epoch
Figure 5: Example loss curves during training a PLRNN (M 10) on a 2d cycle using gradient descent, once without GTF (a = 0, blue curve) but gradient clipping, and once with GTF (a 0.1). Note that without GTF there are several sharp loss jumps associated with bifurcations in the PLRNN parameters, while activating GTF leads to smooth loss curve avoiding bifurcations. Note: For direct comparability both loss curves were eut off at and then scaled [0. Th.The absolute loss is much lower for GTR

---

More importantly, however, we discovered that the recently proposed technique of "generalized teacher forcing (GTF)' 24 tends to circumvent bifurcations in RNN training altogether, leading to much faster convergence illustrated in Fig- 5 The way this works is that GTF, by trading off forward-iterated RNN latent states with data-inferred states according to specific annealing schedule during training (see Appx. A.2.6 tends to pull the RNN directly into the right dynamical regime. In fact, for DTBs we can strictly prove these will never occur in PLRNN training with the right adjustment of the GTF parameter:
Theorem 4. Consider a PLRNN of the form (2) with parameters O = {A, W. h). Assume that it
has stable fixed point k-cycle y (k 2 1) that undergoes a degenerate transcritical bifurcation
(DTB), for some parameter value 0- = 0g E 0.
(i) If |||]]| +|W||| 1, then for any GTF parameter 0 < < 1, GTF controls the system, avoiding a DTB and, hence, gradient divergence at Oo-
(i) If ||A|| + | = r > 1, then for any 1 < 1, GTF prevents a DTB and, hence, gradient divergence at 0g.
Proof. See Appx. A.2.6
As this example illustrates, we may be able to amend training procedures such as to avoid specific types of bifurcations.
6 Discussion
DS theory 3 47 58 is increasingly appreciated in the ML/AI community as a powerful mathematical
framework for understanding both the training process of ML models [49] [54 44 [6] well as the behavior of trained models [62] 33 [1]. While the latter is generally useful for understanding how a trained RNN performs given ML task, with prospects of improving found solutions, it is in fact imperative in areas like science or medicine where excavating the dynamical behavior and repertoire of trained models yields direct insight into the underlying physical, biological, or medical processes the model is supposed to capture. However, application of DS theory is often not straightforward, especially when dealing with higher-dimensional systems, and commonly requires numerical routines that may only find some of the dynamical objects of interest, and also only approximate solutions. One central contribution of the present work therefore was the design of novel algorithm, SCYFI, that can exactly locate fixed points and cycles of wide class of ReLU-based RNNs. This provides an efficient instrument for the DS analysis of trained models, supporting their interpretability and explainability.
surprising observation was that SCYFI often finds cycles in only linear time, despite the com- binatorial nature of the problem, a feature shared with the famous Simplex algorithm for solving linear programming tasks [36] 32 57]. While we discovered numerically that SCYFI for empirically relevant scenarios converges surprisingly fast, deriving strict theoretical guarantees is hard, and so far we could establish stronger theoretical results on its convergence properties only under specific assumptions on the RNN parameters. Further theoretical work is therefore necessary to precisely understand why the algorithm works so effectively.
In this work we applied SCYFI to illuminate the training process itself. Since RNNs are themselves DS, they are subject different forms of bifurcations during training as their parameters are varied under the action of a training algorithm (similar considerations may apply to very deep NNs). It has been recognized for some time that bifurcations in RNN training may give rise sudden jumps in the loss [44 16], but the phenomenon has rarely been treated more systematically and mathematically. Another major contribution of this work thus was to formally prove a strict connection between three types of bifurcations and abrupt changes in the gradient norms, and to use SCYFI to further reveal such events during PLRNN training on various example systems. There are numerous other types of bifurcations (e.g., center bifurcations, Hopf bifurcations etc.) that are likely to impact gradients during training, only for subset of which we could provide formal proofs here. As we have demonstrated, understanding the topological and bifurcation landscape of RNNs could help improve training algorithms and provide insights into their working. Hence, a more general understanding of how various types of bifurcation affect the training process in diverse range of RNN architectures is a promising future avenue not only for our theoretical understandings of RNNs, but also for guiding future algorithm design.
10

---

Acknowledgements This work was supported by the German Research Foundation (DFG) through individual grants Du 354/10-1 & Du 354/15-1 to DD, within research cluster FOR-5159 ("Resolving prefrontal flexibility"; Du 354/14-1), and through the Excellence Strategy EXC 2181/1 -390900948 (STRUCTURES). We also thank Mahasheta Patra for lending us code for graphing analytically derived bifurcation diagrams and state spaces.
References
[1] J. Ahlgren. The probability distribution for draws until first success without replacement. arXiv preprint arXiv: 1404. 1161, 2014.
[2] A. Alahi, K. Goel, V. Ramanathan, A. Robicquet, L. Fei-Fei, and S. Savarese. Social Istm: Human trajectory prediction in crowded spaces. In Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition (CVPR), June 2016.
[3] Kathleen T. Alligood, Tim D. Sauer, and James A. Yorke. Chaos: An Introduction Dynamical Systems. Springer, New York, NY, 1996.
[4] S. Arora, Z. Li, and A. Panigrahi. Understanding gradient descent on the edge of stability in deep learning. In International Conference on Machine Learning, pages 948-1024. PMLR, 2022.
[5] V. Avrutin, M. Schanz, and Banerjee. Occurrence of multiple attractor bifurcations in the two-dimensional piecewise linear normal form map. Nonlinear Dynamics, 67:293-307, 2012.
[6] V. Avrutin. L. Gardini, I. Sushko, and F. Tramontana. Continuous and Discontinuous Piecewise Smooth One-Dimensional Maps: Invariant Sets and Bifurcation structures, volume 1. World Scientific, Singapore, 2019.
[7] Banerjee, MS. Karthik, Yuan, and A. Yorke. Bifurcations in one-dimensional piecewise smooth maps-theory and applications in switching circuits. IEEE Transactions Circuits and Systems I: Fundamental Theory and Applications, 47(3):389-39) 2000.
[8] Y. Bengio, P. Simard, and P. Frasconi. Learning long- term dependencies with gradient descent is difficult. IEEE Transactions on Neural Networks, 5(2):157-166, 1994. doi: 10.1109/72.29918I.
[9] S.h Bouktif, A. Fiaz, A. Ouni, and M. A. Serhani. Multi-sequence Istm-rnn deep learning and metaheuristics for electric load forecasting. Energies, 13(2), 2020. ISSN 1996-1073. doi: 10.3390/en 3020391.
[10] M. Brenner, F. Hess, J. M. Mikhaeil, L. F. Bereska, Z. Monfared, P. Kuo, and D. Durstewitz. Tractable Dendritic RNNs for Reconstructing Nonlinear Dynamical Systems. Proceedings of the 39th International Conference Machine Learning, pages 2292-233222 2320. PMLR, June 2022. ISSN: 2640-3498.
[11] S. L. Brunton, J. L. Proctor, and J. N. Kutz. Discovering governing equations from data by sparse identification of nonlinear dynamical systems. Proceedings of the National Academy of Sciences, 1 11(15))3932-3937 2016.
[12] K. Champion, B. Lusch, J. N. Kutz, and S. L. Brunton. Data-driven discovery of coordinates and governing equations. Proceedings of the National Academy of Sciences, 116645):22445-2251,, 2019. doi: 10. 1073/aa// maaaas..99999 19069951 16.
[13] V. Chatziafratis, S. G. Nagarajan, and Panageas. Better width trade-offs for neural networks through the lens of dynamical systems. In Hal Daumé III and Aarti Singh, editors, Proceedings of the 37th International Conference on Machine Learning, volume 119 of Pro- ceedings of Machine Learning Research, pages 1469-1475 PMLR, 13-18 Jul 2020. URL .20m.hmml
[14] V. Chatziafratis, I. Panageas, and S. Sanford, C.and Stavroulakis. On scrambling phenomena for randomly initialized recurrent networks. Advances in Neural Information Processing Systems, 35:18501-18513, 2022.
11

---

[15] B. M. de Silva, K. Champion, M. Quade, J. Loiseau, J. N. Kutz, and L. Brunton. Pysindy: python package for the sparse identification of nonlinear dynamics from data, 2020. URL
https: http//aww..rgg.o.oo...oo abtgg///s/// 2004. 08424
[16] K. Doya. Bifurcations in the learning of recurrent neural networks. Proceedings 1992 IEEE International Symposium on Circuits and Systems, 6:2777-2780, 1992.
[17] D. Durstewitz. Implications of synaptic biophysics for recurrent network dynamics and active memory. Neural Nerworks, 22(8):1 189-1800. 2009. ISSN 0893-6080. doi: https://do..or
n0.1016/..eeunee2200...0... 07 6. URL https: / //www/ sci encedi rwwt.ccom/cccceee.............ot sci science/ 62.2 Cortical Microcircuits.
art att.ccllepp11/003350000000
[18] D. Durstewitz. state space approach for piecewiselinear recurrent neural networks for reconstructing nonlinear dynamics from neural measurements. PLOS Computational Biology, 13(6), 2017.
[19] A. Ganguli and S. Banerjee. Dangerous bifurcation at border collision: When does it occur? Physical Review E, 71, 2005.
[20] L. Gardini, D. Fournier-Prunaret, and P. Chargé. Border collision bifurcations in a two- dimensional piecewise smooth map from a simple switching circuit. Chaos: An Interdisciplinary Journal of Nonlinear Science, 21, 2011.
[21] M. Golub and Sussillo. FixedPointFinder: A Tensorflow toolbox for identifying and characterizing fixed points in recurrent neural networks. Journal of Open Source Software, 3 (31):1003, November 2018. ISSN 2475-9066.
[22] R. Haschke and J. J. Steil. Input space bifurcation manifolds of recurrent neural networks. Neurocomputing, 64:25-38, March 2005. ISSN 0925-2312. doi: (.001//..neuomm2200
030.
[23] D.I Hernandez, A. K. Moretti, Z.g Wei, S. Saxena, Cunningham, and Paninski. Nonlinear Evolution via Spatially-Dependedl Linear Dynamics for Electrophysiology and Calcium Data, June 2020, 1 [es, q-bio, stat].
[24] F. Hess, Z. 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. URL https: hhap::://proceeingge
mmr...r...........mm.........m.m mlr Satesss/v02//hesee httml
[25] S. Hochreiter and J. Schmidhuber. Long short-term memory. Neural Computation, 9(8): 1735-1780, 1997.
[26] T. Hossen, A. S. Nair, R. A. Chinnathambi, and P. Ranganathan. Residential load forecasting using deep neural networks (dnn). In 2018 North American Power Symposium (NAPS), pages 1-5, 2018. doi: 10. 109/NAPS.2001 2018.8600549
[27] S. Kanai, Y. Fujiwara, and Iwamura. Preventing Gradient Explosions in Gated Recurrent Units. In Advances in Neural Information Processing Systems, volume 30. Curran Associates, Inc., 2017.
[28] H. Kantz and T. Schreiber. Nonlinear time series analysis, volume 7. Cambridge university press, 2004.
[29] G. E. Katz and J. A. Reggia. Using Directional Fibers to Locate Fixed Points of Recurrent Neural Networks. JEEE Transactions on Neural Networks and Learning Systems, 29):3366 3646, August 2018. ISSN 2162-2388. doi: 10. 1099TNNI..2001..2770. IEEE Transactions on Neural Networks and Learning Systems.
Conference Name:
[30] G. Koppe, H. Toutounji, P. Kirsch, S. Lis, and D. Durstewitz. Identifying nonlinear dynamical systems via generative recurrent neural networks with applications to fmri. PLOS Computational Biology, 15(8) 2019.
12

---

[31] T. Lang and M. Rettenmeier. Understanding consumer behavior with recurrent neural networks. 2017.
[32] D. G. Luenberger and Y. Ye. Linear and Nonlinear Programming. Springer International Publishing, 2016. doi: 100079978--3311-18882222
[33] N. Maheswaranathan and D. Sussillo. How recurrent networks implement contextual processing in sentiment analysis. In Proceedings of the 37th International Conference Machine Learning. volume 1 119 of Proceedings of Machine Learning Research, pages 6608 6619. PMLR, 13-18 Jul 2020.
[34] N. Maheswaranathan, AH. Williams, MD. Golub, S. Ganguli, and D. Sussillo. Reverse engineering recurrent networks for sentiment classification reveals line attractor dynamics. Advances in neural information processing systemav, 32: 15755,, 2019.
[35] R Marichal, JD Piñeiro, and E González. Study of fold bifurcation in discrete recurrent neural network. In Proceedings of the World Congress on Engineering and Computer Science, volume 2, 2009.
[36] N. Megiddo. On the complexity of linear programminng page 225-268. Econometric Society Monographs, Cambridge University Press, 1987. doi: 110017/CCOLL05511444446
[37] J. M. Mikhaeil, Z. Monfared, and D. Durstewitz. On the difficulty of learning chaotic dynamics with RNNs. In Advances in Neural Information Processing Systems, 2022.
[38] Z. Monfared and D. Durstewitz. Existence of n-cycles and border-collision bifurcations in piecewise- linear continuous maps with applications to recurrent neural networks. Nonlinear Dynamics, 101(2): 2020.
[39] Z. Monfared and D. Durstewitz. Transformation ReLU-based recurrent neural networks from discrete-time to continuous-time. In International Conference on Machine Learning, pages 6999-7009. PMLR, 2020. ISSN: 2640-3498.
[40] Z. Monfared, Z. Dadi, N. Miladi Lari, and Z. Afsharnezhad. Existence and nonexistence of periodic solution and hopf bifurcation of a tourism- tourism-based social-ecological system. Optik, 127 (22): 199110998--188,,, 2016.
[41] H. E. Nusse and J. A. Yorke. Border-collision bifurcations including "period two to period three" for piecewise smooth systems. Physica D: Nonlinear Phenomena, 57(1-2):39-57, 1992.
[42] H.E. Nusse and J.A. Yorke. Border-collision bifurcations for piecewise smooth one dimensional maps. Int. Bifure. Chaos, 5(1):189-20) 207, 1995.
[43] S. H. Park, B. Kim, C. M. Kang, C. C. Chung, and J. W. Choi. Sequence-to-sequenee prediction of vehicle trajectory via Istm encoder-decoder architecture. In 2018 IEEE Intelligent Vehicles Symposium (IV), pages 1672-1678, 2018. doi:
2018 8500658.
[44] R. Pascanu, T. Mikolov, and Y. Bengio. On the difficulty of training recurrent neural networks. In Proceedings of the 30th International Conference on Machine Learning, volume 28 of Proceedings Machine Learning Research, pages 1310-1318, Atlanta, Georgia, USA, 17-19 Jun 2013. PMLR.
[45] D. Patel and E. Ott. Using machine learning to anticipate tipping points and extrapolate to post-tipping dynamics of non-stationary dynamical systems, 2022. URL httphhtps///aaatt org/abs 2207 222772002711
[46] Pathak, B. Hunt, M. Girvan, Z.n Lu, and E. Ott. Model-free prediction of large spatiotempo- rally chaotic systems from data: A reservoir computing approach. Phys. Rev. Lett., 120:024102, Jan 2018. doi: 10.1 Reeehyyyeeeee 120.024102.
[47] L. Perko. Differential Equations and Dynamical Systems. Springer New York, NY, 2001.
[48] P. Razvan. On Recurrent and Deep Neural Networks. PhD thesis, 2014. Université de Montréal.
13

---

[49] A. Rehmer and A. Kroll. The effect of the forget gate on bifurcation boundaries and dynamics in Recurrent Neural Networks and its implications for gradient- optimization. In 2022 International Joint Conference on Neural Networks (JJCNN), pages 01-08, July 2022. doi: 10.1 (009//CCNN55064.22222229999 28922299 ISSN: 2161-4407
[50] A. H. Ribeiro, K. Tiels, Lu. A. Aguirre, and T. Schön. Beyond exploding and vanishing gradients: analysing RNN training using attractors and smoothness. In Proceedings of the Twenty Third International Conference on Artificial Intelligence and Statistics,, pages 2370-2380. PMLR, June 2020. URL httpp::/pprcceedaiiit ISSN: 2640-3498.
mlr /tt.
[51] S. H. Rudy, Brunton, L.. Proctor, and J. N. Kutz. Data-driven discovery of partial differential equations. Science Advances, 3(4))e1602614, 2017. doi: 10.1 (0.1126cciadv.160226
[52] D. Rumelhart, G. E. Hinton, and R. J. Williams. Learning representations by back propapropgggting errors. nature, 323(6088)5333-536 1986.
[53] Sauer, J. Yorke, and M. Casdagli. Embedology. Journal of statistical Physics, 65:579-616 1991.
[54] A. M. Saxe, J. L. McClelland, and S. Ganguli. Exact solutions to the nonlinear dynamics of learning in deep linear neural networks, 2013. htt/:/tppp///xx...o/o//////
6120
[55] 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.
[56] Smith, S. Linderman, and D. Sussillo. Reverse engineering recurrent neural networks with Jacobian switching linear dynamical systems. In Advances in Neural Information Processing Systems, volume 34, pages 16700-16713. 13. Curran Associates, Inc., 2021.
[57] D. A. Spielman and S. Teng. Smoothed analysis of algorithms: Why the simplex algorithm usually takes polynomial time, 2001. URL httth:://aaxiiii org/ htng//bs/ccss 0c1101050
[58] S. H. Strogatz. Nonlinear Dynamics and Chaos: Applications to Physics, Biology, Chemistry, and Engineering: With Applications to Physics, Biology, Chemistry and Engineering. CRC Press, 2015.
[59] Sushko and Gardini. Center bifurcation for a two-dimensional piecewise linear map, in business cycles dynamics. Models and Tools, eds. Puu, T. & Sushko, I. (Springer-Veriee NY), pages 49-78, 2006.
[60] I. Sushko and L. Gardini. Center bifurcation for two-dimensional border-collision normal form. International Journal of Bifurcation and Chaos, 18:9--0105 2008
[61] D. Sussillo and O. Barak. Opening the Black Box: Now-Dimensional Dynamics in High- Dimensional Recurrent Neural Networks. Neural Computation, 25(3):626 649, March 2013. ISSN 0899-7667, 1530-888X. doi: 10. 10.116//ECC_a_MM
[62] E. Turner, K. V. Dabholkar, and O. Barak. Charting and navigating the space of solutions for recurrent neural networks. In Advances in Neural Information Processing Systems, volume 34, pages 35320-25333. Curran Associates, Inc., 2021.
[63] P. R. Vlachas, W. Byeon, Z. Y. Wan, T. P. Sapsis, and P. Koumoutsakos. Data-driven forecasting of high-dimensional chaotic systems with long short-term memory networks. Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences, 474(2213):201 70844, 2018. doi: 10. 00.1099/ssaa.201..0444
[64] P.R. Vlachas, J. Pathak, B.R. Hunt, T.P. Sapsis, M. Girvan, E. Ott, and P. Koumoutsakos. Backpropagation algorithms and reservoir computing recurrent neural networks for the forecasting of complex spatiotemporal dynamics. Neural Networks. 126:191-217, 2020.
14

---

[65] C. Wang, D. Han, Q. Liu, and S. Luo. A deep learning approach for credit scoring of peer- to-peer lending using attention mechanism Istm. IEEE Access, 7:2161-2168, 2019. doi: 10. 1109/CCESS.20 20......281288.8811
[66] D. Zhang, H.and Liu and Z. Wang. Controlling chaos: suppression, synchronization and chaotification. Springer Science & Business Media, 2009.
[67] Xun Zheng, Manzil Zaheer, Amr Ahmed, Yuan Wang, Eric P. Xing, and Alexander J. Smola. State space LSTM models with particle MCMC inference. CoRR, abs/1711.11179, 2017. URL
attp://arxis: org abs 1tt.//111..1177111
[68] X. Zhu, Z. Wang, X. Wang, M. Zhou, and R. Ge. Understanding edge-of-stability training dynamics with a minimalist example. In The Eleventh International Conference on Learning Representations, 2023. URL
15

---

A
Appendix
A.1
Analysis of bifurcation curves
A.1.1 PLRNNS
The standard PLRNN
18 given in eq. in sect. 3.1 was defined by
2 F0(
- + W Wo(zz-1) + Cs:
where o 1
max 0). There are various extensions of this basic architecture like the dendPLRNN [T] or the 'shallow PLRNN' (shPLRNN) (24, as used in sect. [5.1 for training on single cell membrane potentials. The latter is essentially 1-hidden-layer version of the form
2, with W E R&-x and W2 RLXM L 2 M, connectivity matrices, h1 E RM,, h2 E R bias terms, and all other parameters and variables as in eq. ① While this formulation is beneficial for training, the shPLRNN can essentially be rewritten in standard PLRNN form (see (24)).
= 1 + Wio 2 + h2) + Cs. + h1;
(6)
Assume that Datt) = diag(da()) is a diagonal matrix with an indicator vector do())) = (d1, d2. dx)) such that dm(2m,, =: dm = whenever 2m, 0, and zero otherwise. Then eq.
(I) can be rewritten
(A W (D(t-1)) + Css + h =: Wn(-1) #1 + Cst + h.
Let us ignore the inputs for simplicity. There are 2 different configurations for matrix Dn(t11) and so 2 different forms for matrix Wn(t-1) in the system
Fo Wn(--)) (t-t-1)2-11 2t-1
+ h.
Thus, the phase space of the system is divided into 2 sub-regions corresponding to the indexed matrices
WR + W Dn.
== 1,2,
(8)
see
38
39
for more details. For M = 2, assuming
W11
9)
0
W
0
(9)
1/21
in
we have
a11
0
Woi
Wns =
A,
022
(a11 + W11
0
Wa Was
(10)
W21
022
Hence, for this parameter constellation, the map simplifies there exists only border which divides the phase space into two distinct sub-regions, such that (7) can be rewritten as a map of the form
=1.
((21,-1122,-11)
F2
a1
c
21,t-1
111
+
; 21,-1 S
b
22,1-1
h22
Ac
(11)
2
h1
T, 22,-111 =
+
; -1,111 2 20
b
22,0-2
h22
Ar
16

---

with = d11, a, d11 +11+ b, W21, d = 222. b: c 0. The map (11) is PWL dynamical system whose phase space is split into left and right half-planes (sub-regions) by the borderline Σ (z2-axis). Note that bifurcation curves of the 2d PLRNN 7 in the (411, @11 + wn >-parmeeter space can be determined analogous to those of the PWL map 11) in the (ar, a, J-arammteeee space.
Another way to simplify the PLRNN to a 2d (M = 2) PWL map with just single border is to (Ø1 (31,--1) B x2,--1" where R
remove one of the ReLU nonlinearities and define o( and Φ1 is some variant of the ReLU function such as the leaky or parametrie ReLU given by
2;
>
(1
(a E R).
(12)
a
2 00
Then Dst)) diag(d1, B) such that
1;
21, > 0
di (1(211)) =: d1
(13)
a;
21,1 50
and so
W = Wa = (111 + oww
Bw12
d)
= br
O2221
022 + BB2222
d11 + W1
Bw12
ar
Wn = Wn
(14)
W/1
022 + Bu022
by d)
This gives another example of a representative of 2d PWL maps with only one border defined in eq. (TT) We are pointing this out because eq. [11) is generic system considered more widely in the
discrete dynamical systems literature | and also was the basis for the analyses below and in
Fig. 0
A.1.2 Fixed points of the map
11 and their bifurcations
For at ar. b,, br. C, d, h1, h2 E R, the map has the following two fixed points
OCTR- = (2) C/R L/R (1-
(1
d) hr ch2 byrr h1 + (1 - a1/r) h2
a/rr - b1re (1- d)(1 - a1/r) byjr ) (15)
T
(1- d)
The fixed points Og and OR exist iff of < 0 and xF Hence, the existence regions of admissible fixed points are
respectively; otherwise they are virtual.
Eos
{ her, cc,,)
(1 - d) h1 + cc2 (1 d)(1 a1) bi
=0}.
Eor
{ (1 {
(h, h2, ar br, e, d)
- d) h1 + ccl22
0
(16)
(1 d)(1 - ar) - br C
Let DL/R be the determinant and Te/R the trace of A C/TR, and
Pa/R(A)= x² (a,ar + d)A + a1/r b1/r c X2-TT./7 + DejR.
(17)
its characteristic polynomial. The corresponding eigenvalues are given by
\1.2(Oc/1)) a1/r a1/r + + Vlayp - d)₂ + +byjcc Te/R + V VēIR -
3-RR-DD0m 4D2/R
(18)
2
which are always real for b1,, c 2 O, while for b1, < 0 they are real provided that lays d >
//-bb For complex conjugate eigenvalues of A A/RR obviously |X|2 D2/RR- Thus computing the real eigenvalues, the stability condition for the fixed points determined as
-(1 + DE/R) < TE/R < 1 + DE/R.
(19)
Accordingly, the stability region of the fixed points Oc and OR can be obtained by PE/R (+1) 1 F (at/r + d) + 01/r d byr > 0 and DC/RR < as
SE/R = { (1i h2, ai/rr at,
bysrs c, d) E EOcra layr d b1sr c << 1,
17

---

1 + (ar/+ + atfr d bifr <>0p of.
(20)
Note that when De/RR < O, all the eigenvalues real and there cannot be any spiralling orbit.
Remark 1. Consider the PLRNN
with M 2. For the parameter setting i.e.
011
(011 + W11
W Wn
Ac
Wi WΩ
AR.
022
W21 a22
(21)
the two fixed points OC/R
(=
.22 122 2 are given by
h1
h22
h1
U21 h (1 a1 wn1) h2
T
O =
OR
01
1 - (22
01
W11
(1 a12) a1 wi1)
(22)
Hence, the existence regions of admissible fixed points are
Eo = ={ (hi, a22)l
h1 h-< -=<00
of
EOx (th1,ass
h1
a00 022; a22,11))
0
i-an ai
1-a11 a11 W11
and their stability regions can be obtained as
Sc = (+=(thn,en, (111))))022)) a22) E Eoz en a22 < 1, 1 + (an + azz) + am a22 > o0.
(23)
ST. {{ (h,, (11, wi) E Eor (an + wn )a2 1, + (an + wi a22) + (an + wi Jaza >
0
Remark 2. Ifb c= O, then 20,2(021R)) are real and the stability regions Sc/R become
al/ (1h ai/r bi/rs d) E EOeER by,, -1 0,-1, --150/1 1, -1 - (24)
SC/R
The fixed points are regular saddles for all parameters that belong to
e, d) E Eoeya airr > 1, aifr d airr + 1 b1j < air a.rd). d)
atr byr byr
(25)
and in this case (0 > 1, (2(OR/)) one negative eigenvalue) if parameters are in
Furthermore, they are flip saddles (i.e., with
(bi, (bi, h2, ater byrr c, d) E Eo |a/r + d > 1, air d < bisr c < air d + air + d + 1 U
h11, al/r b1r, Eocyr ad-a d + by airr d + aifr + d + 1,
0 0<ayr + +d 11, 01/ )
(26)
for which (1(O)/) 1, -1 < A2(Oc/R) < 0. as well as in
(hi S0r, byr, E Bogya uijr + d 5 -1, 01/1 d < bur < ai/r d Aifr d + 1)U
h2,, a1r, byrr C, d) Fozyr aijr + ai/r + < byre e < ayrd ai/r + 1
-1< + <<0 <0)
(27)
18

---

such that 0 < <1(00/R)) <1, 12(O2/R) < 11
When byrr c < O and lay/ d) < / byy the eigenvalues are complex conjugates and both Og and OR are spirally attracting (attracting focus) if at/y d d-bddbb--. b1/r c. < 1. In this case, if 01/r + d 0 then they are clockwise spiral, while for air + d < 0 the spiralling motion will be counterelockwise. Moreover, for a - byre c > 1 they are repelling foci. Finally, for atir d bt/r 1, the fixed points are locally centers and they undergo a CB at the following boundaries:
Cc = (ar by, c, d)) byc < O, |a - dl < 2V bi as b bc==1.
CR { ( c, by 0, aa, d < 2V bb, e, U2 d by c=
(28)
At these boundaries, the fixed points lose their stability with a pair of complex conjugate eigenvalues crossing the unit circle. For the parameters belonging to CL/R. the Jacobian JC/R is rotation matrix whose determinant is equal to 1. In this case, J//R can be determined by a rotation number which is either rational (2) or irrational (p). Therefore, in some neighborhood of Oz/R. there is region filled with invariant ellipses such that they are periodic with period p (if the rotation number a rational number 2 quasiperiodic (if the rotation number is an irrational number p); for more information see 660 59
For (1 d) h1 + ch2 o O, at the boundary
Tc (hi, h2, at, bi, c, d)| 1 at - d a aid b.c
(29)
the fixed point Og undergoes a DTB, since, if the parameters tend to Tc, then O2 + to0 and x(Ox) - 1. Similarly, for (1 d) h4 + ch2 # 0, DTB occurs for the fixed point OR at the boundary
TR = -
( (1, h2, ar, by, c, d)) a, d or d do.dd-bee by b,
(30)
A DTB of fixed point results in its disappearance, as in this case the fixed point becomes virtual which may lead to changes in the global dynamics [6]. Furthermore, the BCB curves are given by
&c ( {bi ar, bi, e, (1 d)(1 a1) b 0, (1 d)hi hi + ch2 = =o]. (31)
and
ER = = (1,,,,,,ar, ar, br, c, d) (1 d)(1 a,) b, + 0, (1 d) hi + ch2 = 0).
(32)
In addition, the DFB curves for the fixed points Oc and OR are
Fc (h1, at bi, d)) at d + at by bic=
FR (hi,, ar, by d) 1 + d a. by c= 0f.
(33)
Remark 3. The existence regions Eo, and For are bounded by the BCB curves &c and &x- Remark 4. The stability regions SE/R of fixed points (eq. (20)) are bounded by the DTB curves
TE/R (eqn. (29) and (30)) the DFB curves F FE/R (eq. 33 and the CB curves 0ijr d buy c 1.
For instance, for d I, are illustrated in Fig. SP a- In this case, the stability regions only exist for byjr c < 0 and aifr - d -2 Moreover, shown in Fig.
(1b)., Sl ()), for d = 0.01,
these stability regions can exist for both cases bIjo < ai/r d < -22 bIfy (in blue), and byr < 0, 01/r d 2b/- (in green), but not byfr e 2 O. Furthermore, ife 1, there are stability regions Sc/R for the two cases bitr < 0, aifr-tar - d < -
(in blue), and byj, c > O (in purple); see Fig. S2 a). Finally, when O. as explained in Remark the stability regions have the form (24). i.e.
(hr, h2, ai/r, bifes e, E Foorr c= b1rr ER, -1 S 0t/r: d d
SL/R
h2,
(34)
which are displayed in Fig. S2 b).
19

---

2
0
-2
-
2
2
0
-2
-2
-1
0
-2
-2
bije
bitt
Figure S1: Stability regions SC/R- Left: for d = 1; right: for d = 0.01. The case by/r < 0,
airr -2 /-i/r is plotted in blue, and the case b1 0, 01/r d > 2/ 2 is drawn
in green.
2
1
2
0
0
-2
10
0
-1
0
0
2
-2
-2
-10
-2
byt
01/r
bi/r
01/r
Figure S2: Stability regions Sc/R- Left: for e = 1; the case by, C < 0, 01/r - d < 22/b11/ is
plotted in yellow, and the case by/, >> 0 is drawn in purple. Right: for = 0.
Since the system is linear map in each sub-region L and R. there cannot be any n-cycle, n 2 2, with all periodie points on only one linear side. So, all period-n orbits have both letters L and R in their symbolic sequence.
A.1.3 2-cycles of the map
(11) and their bifurcations
The 2-cycle ORe of the map where
is determined by solving the equation T T2oTR(22,22) = (21, (21,22)
al ar + by ar b + b, d
a cc
21
cch + h1 (as +)
T2)
(35)
d bbb
22
by hi h22 + 1)
In this case if (I - Jc JR) is invertible, then the solution (21, 22) = (z1 (1) 2/ 2 is given by
((1 d)h; h22 (a: + b c
(2f (1) 2/ (1))== (a, - (a by c(b d2 al ar +
h2 + d - as ar by 4: a, d a
(b by + b, d + at by d by b,
(ap d - ca by -c0
d² aq a, +
(36)
Also Tx (2 (2(1), 20))) (z) 2222 yields
(2) (2) T
(=(2) 2/ 12))
(a, ((1 ta
+
by c c+11
b, c) (a: by c(b + b,) d2 a1 a, +
h2 + d - at a, - by c 41 a, d + b, h1 (br + a b, + b d + ar by by by c)
(a, d b, c) by c) c(b: b,) d² a a, +
(37)
Hence, the existence region of the 2-cycle ORe is
Eorc = { (hii ar br, c, d)
((1 d)h, + c h2) (at d ai by c
> 0,
(a, d - b, c) d - by c) c(br + b.) d2 at 0r +
20

---

(a, d d b, d ())) - cbi + b,) - d2 - 4 +1 < of
((1 d)hı + + h2) (a, (a,
d ar by c
(38)
aq + b
act cd)
The characteristic polynomial of JRC Jc JR.
is given by
a, b,
de +bc
Porc (A) - aa bic by c)) (a, c) c).
(39)
and
DRC (a, d by (a) bi c),
PoRe (1) = (a, d b, c) - c(bi b,) - - ar +
PoRe( (-1) (a, b, () by c(bi b,) d ar +1,
(1,((RE)) ai + +ch +
2
a + b/ c) + (b, d212 2( a, by (b, + by + a, bi)
2
(40)
Thus, the stability region of ORc is
(1) 122 01/r, (/1,,/ by/r: e, d) E Eoxc - < (a, d, c) -aadd----/// 11,
SRC
(a, -a, by by - < <(bi + b,) + a <
d da b
(41)
In addition, for ((1- - d)h1 + c h2) (as + d + ai d by c 1) + 0, the set
TRC = ={ (h2, h2, at ar, bi, br, e, d) as or + t c b c d as or dl - by by 6² + a by ed
t bi c -1- 11
(42)
is the DTB curve for the 2-cycle ORC. As in this case, for the parameter values belonging to TRC- the points of the 2-cycle ORC tend too, and the corresponding eigenvalue tends one. Moreover, for (1 - d'hn + c l22 # O, the BCB curves ORe can be computed as
SRc han bi, e, d) az + a
by e + -1
Se (hn by c, or + d + a, d - b, c+ - 1+
(43)
Note that here the condition (1 - d)h1 + h2 * 0 guarantees a regular BCB in the sense that only periodic point of ORe collides with the switching boundary; for more details see [6]. Besides,
FRC =f (h1, 122 as, ar, br, br, d) aia, + b.c + b,c + of + as ad + by b,c²
=10 0).
as by ed a, bied + 1
(44)
is the DFB curve the 2-cycle ORE-
Remark 5. One can see that for (1 - dyh + b2 # 0 the DFB curves of the fixed points Oc and OR (Fe and FR) and the BCB boundaries of the 2-cycle ORC (ete and &2c) are the same. this case, the DFB of the fixed points can lead to the (attracting) 2-cycle ORc-
21

---

A.1.4 3-cycles of the map and their bifurcations
Here, we investigate the existence, stability and bifurcation structure of maximal or basic 3-cycles. Note that for the continuous map (II). basic n-cycles ORC--= (n 2 3) exist in pairs with their complementary cycles R c2-22_)). and they appear via BCBs such that one of them may be attracting and the other repelling 20 38). this case, BCB of basic cycles demonstrates a non-smooth fold
bifurcation which includes stable basic orbit and an unstable nonbasic orbit [6] [ [9]. Furthermore,
the complementary orbits can have nonempty stability regions such that, similar to the basic orbits,
they are bounded by curves of BCBs, DTBs and DFBs 6 51.
Basic 3-cycles ORce and their complementary cycles Orac. The basic 3-cycle ORea can be obtained from the equation Te Te T2( (z1, 22) where
Te o Te TR(21,22) - (rr(a + c) + by (arc
dat + by c) + d( c cd) d(d + bi b) c(a bi bi b(d)
21
(r + by + a 6 d)
hi (hi(brc+ ar + 1) + 1) h2 (are od 1))
+
(45)
hi (b + br + 1)) + h2 (1 c d(d + 1) +
If - JI-Jz JR) is invertible, then the solution 22) = (zi 222 ) is
+ ch2) G1 G2 T
22 =
(46)
G
G
where
G1 = af + and aî 2a bled - abbic + and + and + a + 2/ bied + bic + dl + + 1,
G of a, d² c(aa b a b, a bi + d(2 b.)) + (a, d b, c)( d b c)² + 1,
G2 = h2 + b1h + dh + d h2 + a + bich2 + bi - aj + b,d - b, chi
aia,dh2 + abb hh + byb, 0² h2 a, b/ h2 az + afb d bi b, d²h,
+ a a by h a, ch2 - a, b ch2 + a, by h + b, chi b, cdh2 + ai a bich2
+ + by dhi - a: b
by cdh,
+ 2as a, cd2 2a, by by cdh1.
(47)
Further
Tx( ,2 (1)) = (=) (2) (() = (11 (11- djdn + ch2) K2
G
(2) 22)) = -=13) ( 10) ((1 + ch2) H1 H2 T
Te(=
(48)
G
G
where
K₁ = a, + d + aa, + bic + a,d + az + + aja,d abbcc b,cd + aja, + bib,
a,b,ed - abbicd + 1,
K₂ h2 + b, + dhv + d + abbeli + b,cha + bidh, ara + abb + byd
+ bbb,ch2 + a, h1 + byb, h2 -a af h2 + 16 - abbic - bccc + abiddhy
+ by by chi - bi cd h2 bia, - a bib,ch1 - a, by c dh hbb,cd, + a bid h
ot b,b,,,, a, bdc + cdhz + a by cd by cdli - aha,
H. = at + + ana + and + b,c + and + de + arand - a, bic - bred + d + bbb,
adb hed + 1,
H2 he + bihi + dh2 + d + h2 + bibii + biche + b,dhi a2a,t2 + ch1
22

---

+ b1 afa, 2h2 - biedh + ajbid + bf - bfcch1 - ad h2 + b2
+ aaabii a1bch a1bch2 + dbb + qua,ich22 - abbib,c abbed h2
a arbib, b eddh, + bac + aja, abbıb,cdh1.
(49)
Therefore, the existence region of the 3-cycle OR is given by
EORca {( h2 at, b,, c, ((1 d)h1 ch2) Hh
((1 d)h, + G1 0,
((1 d)hj + +ch2)
0,
G
G
o.
(50)
G
where G, G1, K1 and H1 are defined in and polynomial of
(49
On the other hand, the characteristic
( arf + be
c) cdd)
dia + hc) d(ar c + ed)
JRC = J2 JR
6, bb a bjd)
d(d + b,c) + c(a b + by b(d))
PoRER (A) x2 -afafr + d clabir arb + a, + d(2b: + by)) (a,d brc)(ard byc)?
(51)
According to
DRC = (a: d - bi b)) da,b,,
PORCa (1) aja, d c(a:b: abb, + a, d(2b; + b,))
PoRca aj d cla:bi ar,, + a, d(2b: + b,)) + (a,d - bb,))aad bio)² + 1,
(52)
the stability region of the 3-cycle ORCa is given by
Spe² = { (he hi, airr byr d) E Foxcol - << (a bi)
c) <1,
- (aad byc) (a, b,c) 1 afa, + d + c(abbrr ab, + a, + d(2b: + b,))
( d b) c) b1
(53)
Furthermore, for (1- d1-d + h2 and G1, G2, KI, K2, H1, Ha * the DTB curve for the 3-cycle ORE² is
TRC² = =( (h1 h2 as, ar by br c, d) by of or of d3 + a. 2b, aibre + 20, a by cd
+ b. a b, E2 a, bie²d 2 + a,bi b, cd d - = 0
(54)
For (1 d1-1)h1 + 122 #0
Skaa= hi, K1 a, + d + af + d + aaa
- bab &f by d -bb.d-abbdd-aaa---aa----
= { (n = -ott:aoooooooo--a--oo----------- dl + a a, d
23

---

or bi by cidd at a, de + by c² a bd ar bied byycd+1- =1= (55)
are (regular) BCB curves of ORCs. Furthermore, the set
FRe² = { h2 at ar, by, br, c, d)) - b, ot ed? + + + a. + 2b, 2a, arbie
0)
+ abbic + brarc b, + ur bi ed 2bred + a + by ed + d +1= +1=
(56)
the DFB curve of OR.a. As noted, the basic 3-cycle OR exists in pair with its complementary
cycle ORac, c. Moreover, the existence region of ORa can easily be found by interchanging the letters
C and R in all notations of the equations (45)- (49) and considering
< O,
> 0,
21 O.
(57)
Further, the stability region of the 3-cycle ORac for the parameter values satisfying 57 is given by
SRAC { (11z. h2 afrr byyr c. d) - < (a, d b, c) by c) 1,
(a,d b,c) ( - blc) - < az d c(arbr a,d + aib, d(2b, bi))
(a, d by c) c)
(58)
Notice that whenever the stable 3-cycle ORCa exists, its complementary orbit ORa also exists, but it is unstable. Furthermore, for (1 - d)h1 + e h2 o 0 both the 3-cycles Orce and ORag appear at the same BCB curves (55). On the other hand, the DTB and DFB curves of the 3-cycle ORac are given by
= (hi, =( (bi,ha,au h2, at ar b, by, d. bra ajazd aia; b, cdd b,
TREE
a by + biac by b. a 19 à cb + arb e + d - = 0
(59)
and
FRac = ={ (hi, h2 ar ar br, by, c,d) bia,ed? + araeed + aia + 2bra,be b,
+ a.bbc + a - by 82,0 + az D,.2² d + 2b + a by c cddd + d3 +1 = 0
(60)
respectively.
A.1.5 Multiple attractor bifurcations (MABs) of the map
To detect multiple attractor bifurcations for the map II a straightforward way is to determine the overlapping stability regions of different periodic orbits. For instance, as shown in Fig. [TA for e = 0.8, d = 0.2, b -0.4, b, = 0.5, two stability regions SRÆ and SREa overlap in the (ar, a, p-paammeerr plane (or in the (a11, an wi wwi )-parameter space for the 2d PLRNN (7). Their overlapping region, displayed in yellow, reveals the structure of the (at a )-) -parameter plane. This helps find various MABs. Assuming h2 0 and varying h from negative value a positive one, an MAB of the form
O
O + OR
(61)
occurs in the overlapping region. An example of this kind of bifurcation is illustrated in Fig. A3A.
24

---

Moreover, Fig. S4 indicates, for c 0.9, d 0.3, by = -0.6, b, = -1.44 there are nonempty overlapping regions SRC SR.z and SR.² n nSR. This leads to the occurrence of two different MABs given by (6T) and
02 OR + ORCZ,
(62)
for h2 and 11 changing from negative to positive values. Both of these bifurcations are shown in
Fig. [S3/ and Fig. 3S3B, associated with the points Pnand P2 in Fig. [S4 Note that in Fig. [S4] all the points P. P3 and P4 belong to the overlapping region SR.² n SR (in sky blue). These points are 0.6, b, 1.54, a, 11.8, h2 0, and
related to the parameter values 0.9, d = 0.3, they only differ in the parameter 0₁. In this case, one can see that fixing all parameters and changing only the parameter a,, from P22 Pa. the basins of attraction change. The corresponding basins of attraction for these three points are demonstrated in Fig. S3B (right) and Fig. SS for h1 0.5 (after the bifurcation).
A 2
-2
-0.5
0
0.5
.4
-3
h
B
0.4
-0.2
-1
-0.4
-0.3
-0.2
-0.1
0
0.1
0.2
0.3
0.4
0.5
-2 22
-1.5
-1
-0.5
0.5
h,
0. A) Left: Bifurcation dia- gram for 0₁ -0.44 and 1.8 corresponding the point Pi in Fig. S4 Right: Multistabilitty of the fixed point O and the 3-cycle O after the bifurcation and their basins attraction at h = 0.5. B) Left: Bifurcation diagram for di -0.35 and a, 2.2 corresponding to the point P2 in Fig. S4 Right: Multistability of the 2-cycle O and the 3-cycle O after the bifurcation
Figure S3: MAB at c= 0.9, d 0.3, b = -0.6, = h2
and their basins of attraction at h 0.7.
P.
-1
Gx
Sec
-2
Fre
-3
-3
780
P1
Tech
4
0.5 Figure S4: Analytically calculated stability regions for different parameter setting than used in Fig. I Left: Analytically calculated stability regions SR, SRE and SRE², shown green, red and blue, respectively, in the (011, a11 + wi])-parameter plane for az2 0.3, #22 -1.4. The overlapping regions SRC O SR.2² and SR.² OSR. representing multi-stable regimes, are given in yellow and sky blue. Right: Bifurcation curves for the same parameter settings as determined by SCYFL. Note that SCYFI identifies additional structure (regions demarcated by gray curves) not included in analytical derivations.
-1.5 -1
-0.5
011
0
0.5
-1.5 -1 -0.5 an
25

---

0
-1.5
-0.5
2
Figure S5: Multistability of the fixed point O and the 3-cycle O at c 0.9, d = 0.3, b -0.6, b, = -.544 a, 1.8, h2 0 after the bifurcation and their basins of attraction at hı 0.5. Left: a = 0.06 (point P3 in Fig. S4h: right: Ot 0.26 (point P4 in Fig. SS4.
A.2 Proofs theorems
A.2.1 Proof of theorem
Proof. Let L()) be some loss function employed for PLRNN training that decomposes in time as
Σi=1 Le. Then
dc
T
OC,
Σ
00
00
0C.
0h. dz,
(63)
00
Oz: 00
Denoting the Jacobian of system
2 at time t by
J. 0Fo(z1-1)
8z:
(64)
02,-1
O2--1
we have
Oz 0z.
0z
I+ zt
θz,-1
0+
+
(65)
00 0z,-1
00
00
09
00
where O+ denotes the immediate partial derivative (see 48 for more details), Assume that I' is
k-cycle (A: 2 1) of (2). Thus, a3 is a set of temporally successive periodic points
Ps := (Z -1 Ze-k - (k-1)) = (Z F(z..))....................).))) (Z+.x)}
(66)
such that all of them are fixed points of
Zi+b = (2.) = (o(
(67)
and & is the smallest such positive integer (for k = 1, Pi is fixed point of Fe). Similar to tangent vector Oz be computed as
665 the
Ōzikk Ōz. 0 21th
k-1
Ōz
0 Z1th
II
Jetk-r
(68)
00
dz,
00
00
-=0
Thus, for 2,-x Fd (Z(.x) we have
azu. de
0+ 2p12
II
J.
+
(69)
08
08
Accordingly
adj I Jenk-)
k-1
-1 0 21zm
g+
gtz 2
Ject_
(70)
00
00
PITt-d Jponee (1)
09
r=0
26

---

where PIT:=a - (1) det
I-
8-1 J
Moreover, from
(63) and (70) we have
aL1 00
OLi adj Ozyen
k-1
0+ 2x=x
II J.ierr
(71)
PIr Je (1)
Now, suppose that I' undergoes DTB, such that the fixed cyclic points given by (66 tend to infinity and one of their eigenvalues tends to for some parameter value = 0o- This implies
fir..! J. (1) becomes zero at = 0o and so, due to (70).
0z,.k
goes to infinity. Therefore
00
the norm of the loss gradient, loss function.
&C, ae || tends to infinity at A 8o which results in abrupt jump in the
Let (Zty Ztzg Pta be orbit which converges Kk i.e.
lim d(zt,. [x) n-+00
(72)
Then there exists neighborhood U of T& and &o sub-sequences {z 1m=1, (tam+) mm= (1tem4((-)) 1 of the sequence {zt, Ja such that all these sub-sequences belong to and
a) Zlem+s
F (Z 11+ 0,1,2, ... . 1.
b) lim Zlem+s = 2. 8 == 0,,,,,... ... 1,
c) for every Zt,, E U there some 8 (0,1,2, ... , 1-13 such that 2t, E (2
100
This implies for every zt, E U with zt,, E { Nem m= there exists some ñ E such that Zt,, Zlents and lim Ziento 2, -s Consequently, there exists some Ñ E N such that for every ñ 2 both ZiaRio and 2q.t_a belong to the same sub-region and so the matrices Wi(H24+)) and Wn(t*=sa (8 E 1,, ... k 1}) are identical. Without loss of generality, let s = 0. Since
21x(111) = F* (z(1))),
dzic(at)
azta
II
J J 00
(73)
00
00
0z
Ozica
On the other hand, 1om
lima +00
which results in
azen
k-1
g+ Z(e(n41)
lima-so
I
II Je.n-r)
lim
00
00
adj
I - k- Jern-r
g+ Ztz(at1)
lima-sx
(74)
PTT+ Jacker (1)
00
This means as ñ - oc, for any orbit converging to T₃ the norm of the loss gradient tends to infinity 0 Oo which completes the proof.
A.2.2 Proof of theorem
Proof. Let T& be k--yycle (k: 2 1) of 2 defined by periodic points (66). Suppose further that Ig undergoes BCB for parameter value 0 0g. Hence, one of its periodic points, e.g.
27

---

-,- collides with one border. Therefore, 1mt = 0 for some 1 < m S M by the definition of discontinuity boundaries in 38
39
Similar to the proof of Theorem
for 0 A, W we have
adj
1--1-=0 Jook-
0+ 2p
(75)
00
P (1)
00
in which
0+ 2pnk1 = (p.m) Discokk 2
8wmm
2+ Zpck-1
m), Ren
(76)
00mm
where 1 (m,,)) is an M x M indicator matrix with a 1 for the (n, m)) 'h entry and 0 everywhere else.
a+ z1.k-1 becomes the vector at 0 fp- Consequently,
Since 2m grt
0 0.. due (76) 00
0z and so 0w vanishes at 0 Og. Now it can be shown that at 0 0g the loss gradient
80
goes to zero for every E Br, (the proof is similar to the last part of the proof of Theorem D D
A.2.3 Proof of corollary 1
Proof. For M = 2, let hn # O. Then the DFB curves of the fixed point l31 coincide with the BCB curves of the 2-cycle ORc of the form
Fi Eke {(h1, a11, a22)
a11 + 022 d11022 0),
(77)
or
F2 = 2x = ((h1, fu, a11 W11, W21, 422) + 01l + Wi1 + ag2 + (an + wi1)a22 = 0}. (78)
For > 2, assume that I = { } is fixed point of the system, i.e.
z = (I- Wnup)) -1 h adj(I adj(I-Waurp) h
(79)
PI-Wnep (1)
where PI Woop (1) is the characteristic polynomial of I T-Ways at 1. Let us denote the first Wate)
row of the adjoint matrix of I Wa(r) by adj(I Wa(g)) If adj( Wx())) h j O. then we can analogously demonstrate that the DFB curves of the fixed point align with the BCB curves of the 2-cyeles. This implies that, in accordance with Theorem DFBs of fixed points will also lead to vanishing gradients in the loss function.
A.2.4 Convergence of SCYFI
To ensure that SCYFI almost surely converges, we can simply choose the number of random initializations (i.e., Nout in algorithm [} large enough such that every linear subregion will have been sampled with probability almost T.More precisely, drawing uniformly from the 2 different Da-matrices (linear subregions) for initialization, the probability that particular subregion has not been drawn after 7' repetitions is p (1 2 )". Hence, in order to ensure that all 2 subregions
have been visited with probability
we need r 2 In(1 gdz iterations. Choosing Nout = r,
In(e)
we can thus ensure that SCYFI was initialized in each subregion with probability almost 1, and thus, in the limit, will have probed all subregions for dynamical objects. This argument extends to k-cycles by replacing M by 2 above (strictly, a more precise bound for 2 2 is given by (k-1) (2 1) = 2 (k-1), due to the fact that the PLRNN 2 is linear map each subregion and, hence, cannot have any k-cyclls with all periodic points in only one subregion).
28

---

A.2.5 Proof of theorem 3
Proof. We examine the convergence and scaling behavior of SCYFI for fixed points. A similar argument applies k-cycles where k 1.
Let z be a fixed point of the system, i.e.
zi (I Wa(et)) h.
(80)
zi is true fixed point iff
(d a) 2m,t > O,
Vm E {1,2,-. , [M]
(81)
where Dote (i)) ... , ,dr and < a < is a positive real constant.
For examining SCYFI's efficiency, here we focus on two scenarios that impose specific constraints on parameters 0; other cases remain to be investigated.
Case (I) Let R be randomly generated matrix with uniformly distributed entries in the interval [0, 1), and h * O be random vector with all its components being non-negative. For an arbitrary
C 0, we set
A
diag R),
2 + R
+
- diagia()))
(82)
2 + R2 +
Then
diag R
|||||
2+ R| +
2 | +
| diag( (l))
|R] + diag (R)
1+ |R|
W|| 2 + R +
(83)
2 R +
2 +
[ +
and so ||A|+ + ||| < 1. Therefore
Vt W moll = ||A+ W $ |||||| + ||]|| Dacoll s ||||+ || < 1,
(84)
and so
Vt ((aa)) S Wao <
(85)
In this case, for any N, we also have
II W
s II
Woue ll s (14I++IWI)
+ )* < 1.
(86)
==1
This ensures the stability of all fixed points and k-cycles of the system.
According to 85 we have
- (-Wa(
Σ 0
W 90g)
I Waon)) + (i())
(87)
Hence, all the elements of (I W are positive, and so zm,r; > 0 for every i This implies that all true and virtual fixed points exist within singular sub-region. Additionally, only one fixed point is true, while all the other fixed points are virtual.
Case (II) Let h = (h1, h2, ... h be random vector with all hm uniformly distributed in (0, and
Bmin min (hm hm E 1 5 m < M} 0,
0 < Bmin S 1,
Bmar max hm An E h 1 < m < M > 0,
0 < Smax 1 1.
(88)
29

---

Assume further that R1 is randomly generated matrix with uniformly distributed entries in the interval (-1, 0], and for M
Bmnn
W
RR diag diag(Ri))
(89)
M + R2 +E
Consider
Bmen
Omma = max [|aal] wij E W).
0 S Omma
(90)
M + [ +
and S C. {1, ... Mj = such that K 2--Ccrr((S)
<
2 Suppose that R2
diag(ri, .... TM) is randomly chosen diagonal matrix with Tm uniformly distributed in (-1, 1)
for m E I 1 S. and the other elements (m E S) uniformly distributed in (r* 1,0) where r" (M-1)atmas Amaz
Since
Dmin
05 (M- m0x Binax
(M - -1) Bmar
(M -1)
(M MM-1 1)
S
<1, (91)
Amin
M + |R| +
M + R +E
M
80 -1 15 - 1< O.
A
R2.
(92)
2+ |R +
then
||
|||||
2 R2 +
2 [ +
[|] Bmin
R diag( Ri)
|R diag R
+ R
(93)
M + R2 +
M +
R +
2 + R
+
which implies ||A|| + W| < 1. We set large enough to satisfy the condition
Wn(e))
Waet I+ Wa(z))
(94)
On the other hand, for any / we have
011 W21d (t) Ward (t)
W12d2(t)
Wisds (t)
wiMdM (t)
a22
W23da (t)
w2arda (t)
Wom A M (D()
(t)d2(t)
033
waMdda(t)
(95)
wMid((1)
wM2d2(1)
wM)
aMas
Hence
M Wmy dz(ti)hg = (1 - ama) hm
M
2m.ij (1 + amm) hm +
[mmll dj(ti)h..
(96)
1=1 jmm
jmm
Since for every 1
M | (d((t)hg <
Wmsl hys
(97)
#=1 j7m
1=1 jmm
80
M
Em, 2 (1+ 4mm) hm
Wmy hi
Vm I.
(98)
j=1 4mm
30

---

Moreover ass E (r* - 1, 0), for every E S, and thus
M
Omar Omann Bman
/ s hi
j= [was
j=
(M 1) Pmaz
j
3#3
age 111
(99)
Smin
9min
ha
Therefore, due to (98) and (99p, 2,,17 > O for every L and E S. This means that all true and virtual fixed points only exist within a relatively small number of sub-regions, denoted as K 2M- < 2M Given our specific initalization of 0, in both cases (I) and (II) there is set of different sub-regions, each associated with unique D2(1) matrix. We refer to the entire set of these matrices as
Dx (D1, Dk).
(100)
SCYFI, by its definition, only moves within the sub-regions that have virtual and true fixed points, continuing until it discovers a true fixed point (or gets stuck a virtual cycle). Thus, it can iterate between J < K sub-regions
D- (D.. ... D₂} S DR:
(101)
or within the set of virtual fixed points
2, = (21, (21,...,z.)
(102)
In case (I), there is only one true fixed point. Since all virtual fixed points are located within the same single sub-region, SCYFI's initialization will naturally position it within the correct linear region, requiring no more than iteration. Hence, it needs at most 2 iterations to find the true fixed point. Consequently, SCYFI's scaling is constant.
For case (II), if we suppose that SCYFI follows the virtual/true fixed point structure of the underlying system in these K sub-regions, the necessity for the probability of discovering the fixed point to be close to 1, specifically
is to have
In(e)
In(c)
N 2
In(1
=
in(1
1
(103)
2 [Mrcaad
iterations. Since
S (S) I M so K 2 2 and In(1 k) = K
For e* 2 6, let
N In(e")
Infe)
1. then
[n
In(1 x
N = [ In(E*) In(1 R
In(E*) In(1 )
In _ == cK,
(104)
e
which implies the number of iterations is bounded from above. If, for every M, we choose K small enough, then the upper bound will stay within a linear growth.
A.2.6 Proof of theorem
In GTF (24), during training RNN latent states are replaced by weighted sum of forward propagated states 2, F=( and fata-inferred states e; G (a:) (obtained by inversion of the decoder
model Gz)):
1 a)2, + az
(105)
where 0 a 1 is the GTF parameter (usually adaptively regulated in training, (44). This leads to the following factorization of Jacobians in PLRNN training:
0z
0z, 0E,-1
FF 0E.1 = (1 (1- = (1 - aaaa Way- (106)
JGFF
0z:-1
0ž,-1 0z:-1
02-11
dz--1
Proof. (1) Since ||||| + W < 1, we have
Vt
Waw ||
A + W W < ||||| ||]]| Powll $ [|4|| + |W] S 1,
(107)
31

---

and so
4t p W mo)) S Wo 1
(108)
where p denotes the spectral radius of matrix. this case, for any n E N, we also have
( $
II W Wo(t)) ((4))
S III | Woolll $ [[Al+IWI]) s
(109)
Now, for any 0 the product of Jacobians under GTF is
II 1=1
GTF = (1 - a)" II
Jr
a1 II
Wne
(110)
1=1
1=1
and
GTF 4
= (1-1)) a)) II W Wo(s) = (1 a0-"a)ppp
Way))
(1 a)" 11
(111)
Hence o [1 JCTF has no eigenvalue equal to 1 and so no DTB can occur (see definition of DTB in sect.
which implies for any and < < 1, the product II.
J..
GTP
3
(ii) Let || + W = > 1. then for any n E N we have
Jo
S rn
(112)
and thus
II J
GTFF
(1 a)" II W v
[11 a)) r)
(113)
Since 0 <
1, inserting
1 < a == a into the r.h.s. of 113 again gives
[[==
GTF <
for any n E N, implying that no DTB can occur.
32

---

A.3
Additional results
A
2
OTB
DFB+BCB
DFB
SR
DTB
ar
B
DFB
nn
2
ar
Figure S6: Analytically calculated stability regions for 2-cycle SRC (red), and a fixed point SR. (green) for the 1d skew tent map, defined in the figure, in the parameter plane given by (a,, qj). B) Bifurcation diagram along the cross section indicated by the gray line A, showing BCB and DFB occurring simultancously at -1 and DTB occurring at a, 1.
1.0
F 5.0
0.5
0.0
0.0
-2.5
0.2
0,0
0.4
Q11
011
Figure S7: Results from SCYFI (blue) versus analytical results (orange) for the 2-cycle in Fig. (A11 + W11 = 2). Eigenvalues (left) and location in state space (right) for one of the cyclic points. This confirms that fixed point locations and eigenvalues computed in closed-form and via SCYFI exactly agree, they should.
A.3.1 Scaling analysis
Although the results presented in Fig. suggest that SCYFI's scaling behavior much better than theoretically expected, the fact that is hard to obtain ground truth comparisons for high-dimensional systems (because of the combinatorial explosion) generally makes an extensive empirical analysis difficult. For Fig. / we therefore focused on scenarios for which we can also provide analytical curves for an exhaustive search strategy (eq. (5)) and where we then either examined scaling with cycle-order & for rather low-dimensional systems, or where we explicitly embedded fixed points to search for which allowed us to move to very high dimensionality M. In general we observed that the scaling behavior also depended on the PLRNN's matrix norms and the eigenspectrum of the embedded fixed points, we constructed different scenarios where we varied these factors as well.
To construct fairly well behaved case with low matrix norms, we randomly generated matrices R with uniformly distributed entries the interval [-1, 1] and then normalized by its maximal eigenvalue: We PLRNN parameters A 1 diag R and W
= = diag(R)) R)), and Xupes
Amaz
chose h uniformly the interval [-50, 50). For each of 10 different such systems, we fixed the number of outer loops and inner loops (Nouts N., in Algorithm I such that fixed point would be detected in at least 50/75 independent runs the algorithm, and then determined the total number of linear regions (i.c., across all N different initializations) the algorithm needed to cycle through
33

---

to detect stable fixed point. We also ensured that across all different runs this stable fixed point would be the same. in accordance with our assumptions. The resulting scaling behavior was well fitted by doubly-logarihmmi curve of the form ([In(M)) + (R² R 0.913, < 10-1). This low-matrix norm scenario with a stable fixed point may be seen as kind of lower bound on the scaling.
To embed specific fixed point z", we again start with a matrix R described above and take A diag( Rt and W = (R- diag R)). We then minimize
min A.W.h
x= - ((AA+ W Dn(-)) . x" +h)
(114)
subject to A staying diagonal and W off-diagonal (we observed that adding small Gaussian noise term to the right appearance x eq. [14 which decayed proportionally to the learning rate improved numerical stability in the optimization process). The such constructed PLRNNs generally have several fixed points, but to compute n we only search for the inserted fixed point x" (making eq. (5) directly applicable). This way produced 5 10 systems. initializing R with values in 0.2,, 0.2] (orange curve in Fig.23) 33 or [-1.0, 1.0] (blue curve in Fig. 23), thus effectively restricting certain range.
the eigenspectrum of the fixed point as well as the matrix norms of the PLRNN However, since matrix norms may change during optimization, eq. (114), our procedure is not strictly guaranteed to produce eigenspectra and matrix norms within desired range, which is crucial especially for the first scenario where we wanted to keep norms within a 'typical range' (see below). So here, to ensure consistency among drawn systems and with our assumptions, the mean absolute eigenvalue of the embedded fixed points was kept close to 0.31 = 0.05 and the mean maximum absolute eigenvalue close to 1.25 + 0.13. For 75% of the resulting systems spectral matrix norms were within the range [1.0. 3.0). While this produced matrix spectra typical for trained PLRNNs (> 95% out of 361 PLRNNs trained on various benchmarks and data had spectral matrix norms within [1.0, 3.0)), the second initialization range resulted in unnaturally large matrix norms and hence may be seen as providing a kind of upper bound on SCYFI's scaling behavior. Fig. 58 shows the best case (left: purple curve in Fig. 2 and typical (right; orange curve in Fig. 2 scaling scenarios on linear scale to better expose the scaling behavior and function fits.
6
300
linear subregions
linear subregions
5
200
4
100
-3
0
25
50
75
100
125
50
100
Figure $8: Zoom-ins linear axes of the scenarios with doubly-logarithmic (left: R² ~ 0.913. p
10 10-4 and quadratic (right; R² x 0.925, p < 10 10-5) scaling behaviors from Fig. 2B.
34

---

A
1000
Mean wlinear regions Linear region
750
500
250
0.0
2.5
5.0
7.5
10.0
SCYFI iterations
1000
750
500
250
0.0
2.5
5.0
7.5
10.0
SCYFI iterations
Figure 89: A) Initializing SCYFI in wide array of different subregions (different colors), it quickly converges within just a few iterations to the same set of linear subregions which contain the dynamical objects of interest (fixed points in this case). B) The number of different subregions explored by SCYFI when started from different initializations shrinks exponentially fast with the number of iterations. Shown are means (+ stdv) from 10 different systems with M 10.
A.3.2 Loss jumps & bifurcations in PLRNN training on biophysical model simulations
Here we provide an additional illustration of how SCYFI can be used to dissect bifurcations in model training. For this, we produced time series of membrane voltage and gating variable from biophysical neuron model [17. on which we trained dendPLRNN [IO] using BPTT 52 with
2 bases) has 28 different linear subregions and | 124 parameters. Fig. SIOA gives the MSE loss as function of training epoch (i.e., single SGD updates). The loss curve exhibits several steep jumps. Zooming into these points and examining the transitions in parameter space using SCYFI, we find they are indeed produced by bifurcations, with an example given in Fig. [SOBB. As we had done for Fig. f the main text, since the state and parameter spaces are very high dimensional, for the bifurcation diagram in Fig.SSIOp SSOB all extracted k-cycles (k 2 1), including fixed points, were projected onto a line given by the PCA-derived maximum eigenvalue component, and plotted as a function of training epoch. For the example Fig. STOB, we found that BCB (Theorem) underlies the transition the qualitative dynamics of the PLRNN as training progresses. Fig. STOC illustrates the dendPLRNN dynamics just before (left) and right after (right) the bifurcation point highlighted in Fig. S1OB, together with time series from the true system.
sparse teacher forcing (STF) 37 The dendPLRNN used (M latent states,
More generally, whether a bifurcation associated with vanishing gradients produces a loss jump depends on the system's dynamics before and after the bifurcation point. In the case of BCBs, one possible scenario involves change in stability, fixed point (or cycle) can loose stability as it passes through the bifurcation point. The maximum Lyapunov exponent of an unstable fixed point (or cycle) is positive, resulting in exploding gradients right after the bifurcation point 37 and consequently to very steep slope in the loss function near the bifurcation point as in Fig.
illustrated in Fig .STO) During BCB, a stable
S10
35

---

A
B
0.08
0.2
0.06
0.1
D
Loss
PC1
0.0
0.02
-0.1
-0.2
0
2000
4000
2400
2450
2500
2550 2600
epoch
epoch
C
0.0
0.0
: 0.5
0.5
1.0
1.0
0
50
100 time
150
200
50
100 time
150
200
Figure $10: A) Loss across training epochs for dendPLRNN (M - states, B 2 bases) trained on a biophysical neuron model in limit cycle (spiking) regime. Red dots indicate training epochs just before and after loss jump for which time graphs are given in C and D. B) Bifurcation diagram of the dendPLRNN as function of training epoch, with all state space locations of stable (filled circles) and unstable (open eircles) objects projected onto the first principle component. The loss jump in A is accompanied by bifurcation from fixed point to eyclic behavior. C) Time series of the voltage variable (x1) of the biophysical model (gray) and that predicted by the dendPLRNN (black) before the bifurcation event indicated in B. D) Same directly after the bifurcation event.
A
0.175 S
0.150
0.125
0.100
0.075
170
180
190
200
210
220
epoch
B
DFB+BCB
0.2
0.1
0.0
0.
180
190
200
220
epoch
Figure S11: Loss jump induced by a degenerate flip bifurcation (DFB). A) Loss during a training run of PLRNN (M 5) 2-cycle. The gray line indicates a loss jump corresponding to DFB and a simultancously occurring border collision bifurcation (BCB). B) Bifurcation diagram of the PLRNN, with the DFB and BCB leading the destruction of the fixed point and the emergence of 2-cycle as indicated by the gray line.
A.3.3 Dealing with bifurcations in RNN training
Here are some additional thoughts on how RNN training algorithms could possibly be modified to deal with bifurcations. If the algorithm finds itself during training in a parameter regime which does exhibit the right topological structure, does make sense further dwell within that regime, or possibly anywhere within the vicinity of the current parameter estimate. Unlike standard SGD, the algorithm should therefore perhaps take large leaps parameter space as soon it gets stuck in non-suitable dynamical regime. One possibility to implement this is through 'loak-ahead' mechanism that probes for topological properties of regions not visited so far. While fully fleshing out this idea is beyond this paper, proof of concept that this may speed up convergence is provided Fig. S12 Along similar lines, if we knew the model's full bifurcation structure in parameter space ahead of time, could simply pick parameter set which corresponds to the right dynamics
36

---

describing patterns in the data best. While of course it will in general not be feasible to chart the whole bifurcation structure before training (this is in sense the whole point of training algorithm), it may be possible to design smart initialization procedures based on this insight, e.g. probing topological regimes randomly selected points parameter space before starting training and initializing with parameters that produce a desired type of dynamics (e.g., cyclic behavior) to begin with.
A 1.0
B 1.0
look-ehesa look ahezd
look-ahead look- ahoad
0.9
0.9
Loss
Loss
0.8
0.8
0.7
0.7
1000
1000
epoch
epoch
Figure S12: A) Example loss curves for RNNs trained electrophysiological recordings by BPTT without (blue) vs. with (black) 'look-ahead" (the look-ahead function checks whether there would be bifurcation away from stable fixed point when taking 10x the current gradient step). Dashed yellow line indicates the epoch at which the look-ahead step was executed. B) Average across 6 loss of RNNs trained without (blue) vs. with (black) look-ahead. Error bands SEM.
37


Reviews:
[{"review_id": "8w4LfFI271", "reviewer": "Reviewer_qxJb", "summary": "In this paper, the authors provide a thorough analysis of the bifurcation structure of RNNs and its relationship with discontinuities of the loss function during the training process. The authors also present a novel algorithm for detecting fixed points and k-cycles in ReLU-based RNNs, providing a new tool for understanding the computational and dynamical properties of these networks. ", "questions": "Can the authors say something about other kinds of bifurcations and their connection to loss function behavior during training?\n\nOther questions are listed in the \"weaknesses\" sections", "limitations": "The authors should address the limitations more openly, both in terms of the presented algorithm and in terms of the scope of their theoretical results.", "rating": 6, "confidence": 3, "soundness": 3, "presentation": 4, "contribution": 3, "strengths": "This paper nicely connects dynamical systems theory (DST) and machine learning by studying the problem of RNN training. This subject is known to be a challenging one, and the authors apply DST techniques to prove a rigorous result about the connection between bifurcations and jumps in the loss function during training, a question I would see as quite natural but rather unexplored in this context.\n\nThe paper is nicely written and accessible to the broad audience. Results are clearly stated and novel.", "weaknesses": "From the theoretical standpoint, the authors limit their discussions to two types of bifurcations. A more general and comprehensive discussion would be a nice addition to the paper.\n\nThe proposed algorithm appears to be purely heuristic and clashes a bit with the more rigorous first part of the paper.\n\nThe numerical experiments are a bit limited. The 2d example only qualitatively fits with the stated results: the bifurcation points do not seem to be exactly the points where the phase transition occurs. Is this due to numerical approximation?", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the referee for the very supportive and overall positive assessment of our work!\n\n**Weaknesses**\n\n*“Limited discussion of bifurcation types”*: Yes, true, originally there were two types of bifurcation for which we could offer rigorous proofs, so we focused on them in our presentation. To obtain general results for all types of bifurcations is virtually impossible within the scope of such a paper.  \nHowever, we now managed to derive another proof for a degenerate flip bifurcation (DFB), which we will add to the revised manuscript. In the case of a DFB, gradients also vanish at the bifurcation point (Theorem stated below, proof will be added to ms.). We can also provide another experimental example of this during training, see Fig. R5 in the uploaded PDF, as well as new results on center bifurcations (CB). Based on these new results and observations, we are happy to expand our discussion of other bifurcations in sects. 3.1 and 6 and provide a more general perspective on this topic.\n\n*Proposition*. Assume that the PLRNN (2) has a stable fixed point $\\Gamma_1$ with $B_{\\Gamma_1}$ as its basin of attraction. If $\\Gamma_1$ undergoes a degenerate flip bifurcation (DFB) for some parameter value $\\theta=\\theta_0 \\in \\boldsymbol{\\theta}$, then  $\\lim_{\\theta\\to\\theta_0} || \\frac{\\partial L_t}{\\partial \\theta} || \\to 0$ for every $Z_1 \\in B_{\\Gamma_1}.$\n\n*”Algorithm heuristic”*: Yes, this is true, at the time of submission the algorithmic approach rested on a heuristic, yet a powerful one. We also would like to note that historically there were many powerful algorithms invented (like Simplex or various training tweaks in deep learning) where a more theoretical understanding sometimes lagged behind by (many) years. \nHowever, in the meantime we made some progress here: We can now strictly prove that the algorithm will converge with probability almost 1 under some general conditions on the PLRNN parameters, and will add this to our revision. We also can empirically demonstrate that the algorithm converges exponentially fast to those particular subregions in state space which harbor the dynamical objects of interest (fixed points or cycles), see Fig. R1 in the provided PDF document.\n\n*”Only qualitative fits/ examples limited”*: We will add further examples on bifurcations (DFB and CB) to the Appx., see above and Fig. R5 in the provided PDF.\nAnd yes, in the 2d example in Fig. 3 (which we assume the referee is referring to) the apparent mismatch is a purely numerical (and visualization) issue, partly due to the fact that we cannot apply SCYFI for each and every BPTT updating step (note that here this was done mainly for purposes of illustration, while usually one would use SCYFI only for once trained models such that this issue would not occur).\n\n**Question**\n\nYes, see above, we now have a proof also for Degenerate Flip Bifurcations (DFBs), and will extend our discussion and examples on other types of bifurcation.\n\n**Limitations**\n\nSome of the limitations we think we had discussed, but we agree there are other ones that were perhaps less explicit. So we will add a more explicit pg. on current limitations of the algorithm and approach to the ‘Conclusions’ sect. (incl. the current gaps in theoretical understanding, and the limitation to certain classes of bifurcations)."}, {"review_id": "QUmT06ECk3", "reviewer": "Reviewer_1SGy", "summary": "The main focus of the paper is piecewise-linear RNNs, that is ReLU-based RNNs called PLRNNs, and how certain bifurcations that take place during their training phase will have an adverse effect on their gradients and hence will cause sudden jumps in the loss functions. This helps shed light onto the difficulty of training RNNs and yields new perspectives for the problem of exploding/vanishing gradients.\n\nMoreover, the paper offers a novel way for locating fixed points and k-cycles in PLRNNs which can be helpful in visualizing and exploring the properties of these networks.\n\nThe main technical contributions are Theorem 1 and Theorem 2 which provide sufficient conditions under which the loss gradient will go to infinity or will vanish. These sufficient conditions are directly linked to the bifurcation properties of the PLRNN.", "questions": "Q1: Is there a simple 1-d example that is even simpler that the 2d example you provide? And that helps illustrate your two main theorems?\nQ2: In more applied settings, RNNs often operate in large dimensions. Could you comment on how your results can be helpful in these setting as well?", "limitations": "No significant limitations.", "rating": 7, "confidence": 4, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "+the paper presents much-needed theory on a topic that is interesting both for theoreticians and practitioners\n+conceptually, the work is very interesting and helps clarify an important problem in the literature of RNNs\n+the distinction between the two cases of exploding and vanishing gradient through the lens of bifurcation diagrams is interesting\n+their algorithm for discovering fixed points and k-cycles, though heuristic in nature, seems to be performing well.", "weaknesses": "-the algorithm for finding fixed points doesn't come with guarantees but that's perhaps not the most important weakness.\n\nThe main weakness is in presentation and in the literature review that make the paper hard to read in certain parts. The topic on training RNNs and the difficulty that arises has been studied a lot, but more recently there have been efforts to explain certain phenomena through the lens of dynamical systems.\n\nIn terms of presentation, I believe the paper would benefit by instantiating the theorems 1 and 2, on specific simple examples. You do sth related in lines 142-184 but I would create a small subsection \"Simple Examples\" to help illustrate the points. I believe Figure 1 is hard to understand and would perhaps break it into 2 different figures.\n\nIn terms of comparison to related works, recently we have seen several works that analyze behavior of the networks at the edge of stability or using the notion of k-cycles and chaos to illustrate why RNNs are hard to train. See papers below:\n\n[1] On Scrambling Phenomena for Randomly Initialized Recurrent Networks by Chatziafratis, Panageas, Sanford, Stavroulakis\n\n[2] \"Better Depth-Width Trade-offs for Neural Networks through the lens of Dynamical Systems\" by Chatziafratis, Nagarajan, Panageas,\n\n[3] Understanding Edge-of-Stability Training Dynamics with a Minimalist Example by Xingyu Zhu, Zixuan Wang, Xiang Wang, Mo Zhou, Rong Ge\n\n[4] Understanding Gradient Descent on Edge of Stability in Deep Learning by Sanjeev Arora, Zhiyuan Li, Abhishek Panigrahi\n\nYour paper studies low-dimensional settings, especially in Fig.4 you provided a 2d toy example that helps illustrating the tight association between the loss landscape and bifurcation curves, whereas [1] and [2] studies the 1-dimensional setting giving characterizations from Li-Yorke chaotic systems, a standard notion from the dynamical systems literature. Papers [3] and [4] study the dynamics of GD and establish new connections with the notion of \"sharpness\"; while the emphasis is more on deep learning, the papers are very related as they analyze the dynamics through eigenvalues analyses.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the referee for the very supportive and constructive feedback!\n\n**Weaknesses**\n\n*”Theoretical guarantees for algorithm”*: Obtaining rigorous theoretical guarantees is challenging in the current case (and sometimes, historically, theoretical understanding of convergence and other properties of an algorithm lagged behind by many years, e.g. Simplex or various training tweaks in DL). \nHowever, in the meantime we made some progress here: We can now strictly prove that the algorithm will converge with probability almost 1 under some general conditions on the PLRNN parameters, and will add this to our revision. Furthermore, we can add some numerical results on this topic: As shown in Fig. R1 in the included PDF file, when we initiate SCYFI in different randomly selected linear subregions, it converges to the subregions including the dynamical objects of interest *exponentially fast*. This numerical observation further helps to explain its efficiency. \n\n*”Presentation”*: Yes, this is a fair point. We have now substantially reworked Fig. 1 to exactly demonstrate the kind of situations our theorems are dealing with, see Fig. R2 in included PDF file. We will furthermore rewrite sect. 3.1, moving some of the current material (like that on multistability) to the Appx. and putting a stronger focus on the bifurcations and situations our paper specifically deals with.\n\n*”Literature covered”*: We thank the Referee for providing additional relevant literature (some of which we are indeed familiar with) which we will incorporate into our Related Work section, apologies for this omission! \n\n**Questions**\n\n*Q1*: In principle yes, one can construct 1d cases which essentially boil down to the skew tent map. See Fig. R3 in the included PDF where we provide one such example for different bifurcations in this 1d case. However, strictly, the smallest type of PLRNN which really makes sense is 2d (e.g., because of the off-diagonal nature of the connectivity matrix $\\textbf{W}$ in the original formulation). Also, a 2d representation allows us to bring across certain concepts, like multi-stability in state space or the association with loss landscapes in Fig. 3, a bit better we think. While we are happy to include 1d examples as in the attached PDF, our preference therefore would be to stick mainly with 2d examples.\n\n*Q2*: Our theorems in fact apply to systems of any dimension. Low-dimensional examples in Figs. 1 \\& 3 were chosen mainly for purposes of illustration and visualization, but the insights gained from them remain valid in higher-dimensional settings as well. We will make this more explicit in our revision. Also note that Fig. 4 indeed provides one example of a much higher-dimensional setting (likewise Fig. S6 in the Appx.). Its state space and bifurcation graph were only projected into 2d or 1d, respectively, for purposes of visualization. \n"}, {"review_id": "zzKoBXhtEk", "reviewer": "Reviewer_T1NJ", "summary": "The paper analyses a class of ReLU-RNNs (PLRNNs) from the perspective of bifurcations in dynamical systems. This gives rise to new insights into the training process of RNNs (as specific bifurcations are associated with exploding or vanishing gradients -- a major issue in training RNNs on long-term dependency data), but also into the behavior of the trained RNN. The paper further proposes a novel heuristic algorithm to detect bifurcation manifolds in PLRNNs. Finally, the paper provides empirical evidence demonstrating the theoretical findings.\n", "questions": "As already stated in the weaknesses: is there a scenario, where the gained insights can positively change/influence the training procedure of the PLRNN?\n\nSince deep feedforward neural networks are also dynamical systems (simply without forcing), can the provided theory also be applied to (residual) feedforward neural networks with ReLU activations?", "limitations": "As already explained above, the only limitation/weakness is the limited significance of this work. If the authors can convince me of its significance, I'm happy to further increase my score.", "rating": 7, "confidence": 4, "soundness": 4, "presentation": 3, "contribution": 2, "strengths": "The paper is very well written and easy to follow. The idea of applying bifurcation analysis in such a rigorous and broad manner is novel up to my knowledge. Moreover, the authors manage to leverage this perspective to gain new and interesting insights into the training behavior of a broad class of RNNs. The reviewer further appreciates that the paper does not over-claim but rather truthfully demonstrates and discusses its contributions.", "weaknesses": "There are not many weaknesses. The paper is novel and interesting. However, one might argue that besides a novel theoretical perspective into the training process (restricted to the exploding and vanishing gradient problem, which has been studied in RNNs for several decades now), there is very little significance of this work. Although one can find some theoretical justification using the proposed SCYFI algorithm on why the training exhibits loss jumps, this has no influence on the training itself. More concretely, even though one might know a-priori that the RNN will run into a bifurcation associated with exploding or vanishing gradients, this does not change the training procedure. \nTherefore, in order to argue that the proposed method is not only novel, interesting and well explained, the reviewer would highly appreciate if the authors can provide an experiment or application, where the newly gained insights into the training process has an influence on how the RNN is trained in the end. ", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the referee for the very constructive and fair feedback, and the supportive response!\n\n**Weaknesses**\n\n*”Limited significance for design of training algorithms”*: Although the exploding or vanishing gradient problem has been studied for a while now, viewing the training process from the perspective of bifurcation theory adds to it and can, in our minds, indeed have significant implications for designing novel training algorithms as it highlights the importance of topological features for guiding training. There are at least three lines along which improvements are conceivable, along two of them we can already provide some new results:\n\n1) If the algorithm finds itself during training in a parameter regime which does not exhibit the right topological structure, it does not make sense to dwell within that regime or anywhere within the vicinity of the current parameter estimate. An implication of this observation is that, unlike standard SGD procedures, we should take large leaps in parameter space right away as soon as we find ourselves within a non-proper dynamical regime, possibly by building in a ‘look-ahead’ mechanism that probes for topological properties of regions not visited so far. While fully fleshing out these ideas is beyond a single paper, for this latter strategy we now provide a proof of concept in Fig. R4 of the uploaded PDF file, which shows that a topological look-ahead mechanism can lead to much faster convergence to the desired dynamical regime.\n\n2) Another potentially important insight is that if we knew the model’s full bifurcation structure in parameter space ahead of time, we could simply pick, right away, a parameter set which corresponds to the right dynamics describing patterns in the data best. While of course it will in general not be feasible to chart the whole bifurcation structure before training (this is in a sense the whole point of a training algorithm), we may be able to design smart initialization procedures based on this insight, e.g. probing topological regimes at randomly selected points in parameter space before even starting training, or testing for bifurcations along maximum eigendirections of the system. In fact, the very recently proposed technique of generalized teacher forcing (GTF) often appears to quickly pull the RNN into the right topological regime:\n\n3) We discovered that a recently proposed technique for training RNNs on chaotic dynamical systems, dubbed ‘generalized teacher forcing (GTF)’ by the authors (Hess et al. 2023 ICML), can be adjusted such as to avoid degenerate transcritical (and possibly other) bifurcations (DTB) during training altogether. Specifically, we can strictly prove (Theorem stated below, proof will be included in revised paper) that a DTB will never occur in training if we adjust the GTF parameter $0< \\alpha <1$. This is of course just one example, but it demonstrates that we can actually harvest or amend training procedures such as to avoid specific types of bifurcation. More generally, we observed that GTF tends to circumvent bifurcations, as illustrated in new Fig. R6 in the PDF, leading to much faster convergence to suitable parameter regimes. The way this works is that GTF, by trading off RNN latent states with actual observations according to a specific annealing schedule during training, tends to pull the RNN states directly into the right dynamical regime. Our techniques can explain how this works, and possibly improve on it.\n\nEqually important, however, in our minds, although we showcase SCYFI on bifurcations during training problems, its applications are much wider and general: In most scientific and medical applications of RNN algorithms we are not only interested in obtaining good predictions, but in understanding exactly how and why the RNN model works. By finding all fixed points and cycles (up to some order), as well as bifurcation manifolds, SCYFI adds a level of *interpretability* and *explainability* to trained RNNs that is difficult to obtain otherwise (as, e.g., reflected in Fig. 4). It allows for in depth analysis of trained models and for dissecting their inner workings that helps to understand mechanisms and limitations in trained models. We believe this is a very important application, in particular in medical and scientific contexts. We will clarify and work this out much more in our revision.\n\n*Proposition 1*. Consider a PLRNN of the form (2) with parameters $\\boldsymbol{\\theta}=\\\\{A, W, h\\\\}$. Assume that it has a stable fixed point or $k$-cycle $\\Gamma_k$ $(k \\geq 1)$ that undergoes a degenerate transcritical bifurcation (DTB) for some parameter value $\\theta=\\theta_0 \\in \\boldsymbol{\\theta}$. \\\n(i) If $||A + W || \\leq   1$, then for any GTF parameter $0 < \\alpha < 1$, GTF controls the system, avoiding DTB and so gradient divergence at $\\theta_0$. \\\n (ii) If $||A + W||  =   r > 1$, then for any $ 1 - \\frac{1}{r} < \\alpha < 1 $, GTF prevents DTB and so gradient divergence at $\\theta_0$.\n\n**Questions**\n\n*Q1*: Yes, see above.\n\n*Q2*: This is an interesting idea, and the answer should be yes: As parameters of a DNN smoothly vary, there may be bifurcations in the forward dynamics across multiple layers that are significant for the training process. We will comment on this in our Discussion.\n\n**Limitations**\n\nSee above: We now provide several lines, and provide preliminary results in the included PDF, along which bifurcation theory can help to improve training algorithms.\n"}, {"review_id": "ybanQ4DFS9", "reviewer": "Reviewer_v6FA", "summary": "The paper studies the dynamical systems properties of piecewise-linear recurrent neural networks (PLRNNs). Specifically:\n- Section 3.1 (and its accompanying Appendix A.1-A.5) analytically categorizes the fixed points of 2-dimensional PLRNNs, and analyzes its bifurcations by explaining how the fixed points (and cycles) change as a function of the network parameters. The analysis involves computing the Jacobian of the recurrent unit and classifying fixed points by using its spectrum. The paper specifically identifies instances of multi-stability, where where multiple attractor regions exist.\n- Theorems 1 and 2 in Section 3.2 show that a PLRNN whose unit undergoes a degenerate transcritical bifurcation (DTB) or a border collision bifurcation (BCB) has gradients trending towards infinite and zero respectively as the bifurcation occurs in parameter space. \n- Section 4 introduces the SCYFI algorithm for identifying fixed points and cycles in higher dimensional PLRNNs. SCYFI is a heuristic algorithm that finds k-cycles by using exploring the exponentially large space of linear regions from a randomly initialized space of regions. They find empirically that the algorithm finds fixed points much faster than brute-force algorithms that search naively search each subregion.\n- Section 5 applies SCYFI and their bifurcation analysis to the prediction of a neuron spiking with a PLRNN and the resulting detection of a stable 39-cycle that predicts the spike. They observe a bifurcation in the training procedure where the stable 39-cycle is eradicated due to a DTB that causes the loss to spike (as predicted by Theorem 1's analysis of exploding gradients.) ", "questions": "The paper chooses to ignore external inputs $s_t$ to its RNNs, which seems like a very significant omission, given that the primary applications of RNNs in the past have been to NLP and other sequence-learning tasks. Since [some results](http://proceedings.mlr.press/v97/hanin19a/hanin19a.pdf) suggest that adding random noise RNNs can dramatically change the structure of fixed points, are the results expected to hold in regimes with noisy or structured inputs $s_t$ included? \n\nAs written, the SCYFI algorithm is somewhat ambiguous to the reviewer:\n- $z_n$ and $D_n$ terms are ambiguous. $z_n $ appears to be a known cycle ahead of time (line 1), component obtained by solving equation (7) (line 8), and a cycle already stored in $\\mathcal{L}_{k - m}$ (line 9). Which of these sets are being compared in the RHS of line 9? \n- Why does the algorithm start out by presuming to know some $z_n$ and $D_n$, or is that just setting up notation for future lines? (line 1) \n- When solving equation (7), what precisely is being solved for?  The equation relates $z_k^*$, $W_{\\Omega(k-r)}$, and $h$. What is known ahead of time, and what must be solved for? Also, how do we solve the equation? If we are solving for $D_n$ by solving for $W_{\\Omega(n)}$, then it appears that the equation is not linear.\n- Does the algorithm only work if every $\\mathcal{L}_n$ for $n < k$ is exactly correct? While it's difficult to prove convergence results due to the searching of an exponential space, some kind of theoretical guarantee about what is necessary for the algorithm to be successful would be useful, even if just in an appendix.\n", "limitations": "Besides the omission of $s_t$ from the analysis, limitations of the paper appear to be well-documented by the authors.", "rating": 6, "confidence": 3, "soundness": 2, "presentation": 2, "contribution": 3, "strengths": "As far as the reviewer is aware, the paper provides novel analysis of bifurcations on PLRNNs and identifies innovative applications of their techniques through their neuroscience application. The analysis appears to be correct. The theoretical results effectively draw connections between bifurcation theory and gradient-based training and are validated experimentally with the loss jump coinciding the the DTB bifurcation in Section 5.", "weaknesses": "As of now, the SCYFI algorithm as written is highly ambiguous, and I cannot recommend acceptance as long as that holds. (These ambiguities are noted in the next part of the review.) I am happy to reconsider if these concerns and concerns about presentation are satisfactorily addressed, given the innovative nature of the results.\n\n**Note: the score of the paper was increased from 4 to 6, because the authors rewrote the algorithm more precisely.**\n\nThe results of the paper appear to be novel from a dynamical systems perspective, but their current presentation of the results may be inaccessible for the ML community of NeurIPS. Given the relative novelty of work on the intersection of dynamical systems and ML, clear explanation of dynamical systems concepts and their relevance to ML is key.\n- In particular, the paper could more explicitly trace the importance of being able to identify k-cycles and bifurcations when training neural networks. The application in Section 5 is interesting and it could be discussed more prominently in the intro to motivate the study of k-cycles for ML researchers looking for more justification.\n- The paper would benefit from more rigorously defining the bifurcations under consideration, like BCBs and DTBs, as well as having a visualization in the parameter space of a toy example. \n- The plots in Figure 1 may be unintuitive to NeurIPS readers with a limited dynamical systems background, especially as Figure 1B is not discussed in the text. This applies more broadly to figures in the main body: including more discussion of what concretely is represented in each plot (while potentially pushing some plots to the appendix) would make for a clearer read.\n\nMinor typos and notes: \n- $\\Gamma$ is notationally overloaded: In 3.1, it refers the spectrum of a fixed point. In 3.2, it is a $k$-cycle.\n- On l238, the code's location is referred to by \"[PLACEHOLDER].\" \n- The first paragraph of Section 5 refers to Figure 4, but likely means Figure 3.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the referee for her/his constructive comments and for bringing up a couple of important points we agree need to be addressed.\n\n**Weaknesses**\n\n*”Improving presentation”*: We will 1) take up the referee’s suggestions and add a pg. to the Introduction which primes more clearly for the importance of detecting bifurcations and cycles, 2) revise sect. 3.1 to introduce these concepts more succinctly while deferring some of the mathematical details to the Appx. \n\n*”Toy example \\& better discussion of BCB/DTB, Fig. 1 unintuitive”*: Yes, fair point. We addressed this now by thoroughly reworking Fig. 1 (see Fig. R2 in included PDF file) that now specifically features BCBs and DTBs to bring it into closer agreement with the main results. We will also rewrite sect. 3.1 (see above) and add math. precise def. of BCBs and DTBs. Some of the other material in sect. 3.1, which may have been more tangential to the main messages of our paper, will instead be deferred to the Appx. Figs. 2-4 are all important though, we believe, but we will add further clarification to the legends and the discussion in the text.\n\n**Questions**\n\n*Q1 (“Ignoring inputs”)*: The presence of inputs does not change any of our results. Since inputs do not affect the Jacobians in eqn. (71), (72) and (77), they do not change our two theorems (even if they would affect the Jacobians, Theorem 2 would be unaltered, and Theorem 1 could be amended in a straightforward way). Furthermore, since we are addressing bifurcations that occur during model *training*, from this angle inputs may simply be treated as either additional parameters (if piecewise constant) or states of the system (without changing any of the mathematical derivations). In fact, mathematically, any non-autonomous dynamical system (RNN with inputs) can always and strictly be reformulated as an autonomous system (RNN without inputs), see (Alligood et al., 1996; Perko, 2001; Zhang et al., 2009, Controlling Chaos, Springer).\nWe will make this very explicit now in sects. 3-5 of the revised paper. We also thank the referee for the interesting reference!\n\n*Q2 (“SCYFI algorithm ambiguous”)*: Thank you for pointing out, we fully agree the pseudo-code could have been presented in a much clearer way, and include a thoroughly updated version in the general rebuttal (this will be formatted more nicely in the paper, the openreview boxes were a bit limited here).\nIn brief (lines refer to original version): Line 1 was simply meant to introduce the notation, the cycle is not known a priori. In eq. 7 (line 8) we obtain a *candidate cycle*, not a true cycle yet. And in line 9, the set may not include that specific cycle yet, but there may be other cycles of order $k$ already detected in previous iterations (there are rare instances where the algorithm may stumble across the same cycle via a different route, and we found it’s more cheaply just to check for it again rather than avoiding this to begin with). Furthermore, in line 9 then a consistency check is performed: The candidate cycle obtained by solving eg. 7 is only a *true cycle* 1) if its corresponding set of \\\\(D^\\*\\_l\\\\) matrices agrees with \\\\(D_\\{init}\\\\), and 2) if it is not a superset of an already detected lower-order cycle (or otherwise already in the list).\nWe hope this becomes all much clearer from our revised pseudo-code as provided in the general response to all referees further above.\n\n*Q3*: As noted above, in line 1 we just set up the notation and do not presume to know any specific $z_n$ or $D_n$ ahead of time.\n\n*Q4*: In eq. 7 we solve for a *candidate cycle* \\\\( z^\\*\\_k \\\\) on the l.h.s. (now indicated also with a \\* in the pseudo-code), given all the terms on the r.h.s.: The PLRNN parameters $A$, $W$ and $h$ are provided as inputs, and the set of matrices $D_n$ needed to solve this eq. is the initial guess $D_{init}$, which determines $W_{\\Omega(n)}$. The next line then checks whether the initial guess $D_{init}$ is indeed consistent with the signs of $z^*_k$ in the respective linear subregions (so that the solution in eq. 7 is indeed self-consistent and the cycle is a true one).\nThis should (hopefully) all be much clearer now from our revised pseudo-code that can be found at the top in the general rebuttal.\n\n*Q5*: Because the cycle candidates are determined analytically and the next step just checks consistency, the solutions will in fact always be exact (down to computer precision).\n\n*Q6 (convergence results)*: Obtaining rigorous theoretical guarantees is challenging in the current case (and sometimes, historically, theoretical understanding of convergence and other properties of an algorithm lagged behind by many years, eg. Simplex or various training tweaks in DL like drop-out). \nHowever, in the meantime we made some progress here: We can now strictly prove that the algorithm will converge with probability almost 1 under some general conditions on the PLRNN parameters, and we will add the theorem and proof to our revision. Furthermore, we can add some numerical results on this topic: As shown in Fig. R1 in the included PDF file, when we initiate SCYFI in different randomly selected linear subregions, it converges to the subregions including the dynamical objects of interest *exponentially fast*. This numerical observation further helps to explain its efficiency. \nThe scaling behavior demonstrated in Fig. 2B further alludes to that: Except for quite artificial and extreme situations that we never encountered in a large database of PLRNNs trained on empirical data (see Appx. sect. A.7), the scaling is always well behaved. We will make this more explicit in sect. 4 where Fig. 2 is discussed.\n\n*Limitations*: See our response to this point above.\n\n**Minor points \\& typos**:\n\n1) Will be changed for clarity.\n2) Was meant as a placeholder for our github site where we will make the code publicly available, not yet included to prevent deanonymization.\n3) Thanks for spotting this, will be corrected!\n"}]

Author rebuttal:
We thank all four referees for their generally positive evaluation and supportive response! We appreciate the detailed reading of our manuscript, and the many helpful comments pointing out where our paper needs clarification or extension of results. Besides the detailed reply to all issues raised in our responses below, we have included a PDF with further material and figures.

In brief, here we

- provide a new theorem and numerical results on SCYFI’s convergence properties
- extend the bifurcation results by providing a new theorem and examples for two other important bifurcations (degenerate flip bifurcation (DFB) and center bifurcation (CB))
- updated the pseudo-code for SCYFI for clarification (see below)
- completely reworked Fig. 1 for clarification and to bring it into closer relation with the rest of the text
- provide an example of bifurcations in a related 1-d system
- provide a new theorem and numerical illustrations to demonstrate the implications of bifurcation theory for designing novel training algorithms

We hope these new results together with the point-by-point replies below sufficiently cover all the remaining concerns.

Here is a clarified version of the SCYFI pseudo-code:

**************************
## Algorithm 2 SCYFI
**************************
The algorithm is iteratively run with \\( k = 1 \ldots K_{max} \\), where \\( K_{max} \\) is the max. order of cycles tested \
**Input:** PLRNN parameters \\( A \\), \\( W \\), \\( h \\);
\
&emsp;&emsp; \\( \mathcal{L} = \\{ \mathcal{L}\_n \\}\_{n=1}\^{k-1} \\): collection of all sets \\( \mathcal{L}\_n = \\{ \\{ z\^{(m) }\_{l} \\}\_{l=1}\^{n} \\}\^{M_n}\_{m=1} \\) of all lower order \\( n \\)-cycles discovered so far, where \\( M_n \\) is the number of found cycles \\( \\{ z\^{(m)}\_{l} \\}\_{l=1}\^{n} \\) of order \\( n \\), with corresponding ReLU derivative matrices \\( \\{ D\^{(m)}\_{l} \\}\_{l=1}\^{n} \\).

**Parameters:**
\
&emsp;&emsp;  \\( N_{out} \\): max. number of random initializations; \
&emsp;&emsp;  \\( N_{in} \\): max. number of iterations

**Output:** \\( \mathcal{L} \cup \mathcal{L}\_k \\); \\( \mathcal{L}_k \\): set of all discovered \\( k \\)-cycles

1: \\( \mathcal{L}\_k = \\{ \\} \\) \
2: \\( i \rightarrow 0 \\) \
3: **while** \\( i < N_{out} \ldots \\)  do \
4: &emsp; Select \\( k \\) subregions \\( D_{init} = \\{ D_1, D_2, \ldots, D_{k} \\} \\) at random with replacement \
5:  &emsp;  \\( c \rightarrow 0 \\) \
6:  &emsp;  **while** \\( c < N_{in} \ldots \\)  do \
7: &emsp; &emsp; Solve eq. (7) for a cycle candidate \\( \\{ z^*_{l} \\}\_{l=1}^{k} \\) with \\( W_{\Omega(k-r)} \\) as defined in equation (2) based on \\( D_{init} \\) \
8: &emsp; &emsp; Determine \\( \\{ D^*_{l} \\}\_{l=1}^{k} \\) based on the signs of the corresponding components of \\( \\{ z^*_{l} \\}\_{l=1}\^{k} \\) \
9: &emsp; &emsp; **if**  \\( D_{init} = \\{ D^*_{l} \\}\_{l=1}^{k} \\) (*self-consistency*) &  \\( \forall 1 \leq s \leq k, \forall \\{ z^{(m)}\_{l} \\}\_{l=1}\^{k-s+1} \in \mathcal{L}\_{k-s+1}: \\{ z^{(m)}_{l} \\}\_{l=1}\^{k-s+1} \not\subseteq  \\{ z^*_l \\}\_{l=1}\^{k} \\) **then** 

10: &emsp; &emsp; &emsp; \\( \mathcal{L}\_k \rightarrow \mathcal{L}\_k  \cup \\{ \\{ z^*_{l} \\}\_{l=1}\^{k} \\} \\)  \
11:  &emsp; &emsp; &emsp;  \\( i \rightarrow c \rightarrow 0 \\) \
12:   &emsp; &emsp;  **else** \
13:  &emsp; &emsp; &emsp; \\( D_{init} \rightarrow \\{ D^*_{l} \\}\_{l=1}\^{k} \\) \
14:   &emsp; &emsp;  **end if** \
15:  &emsp; &emsp;  \\( c \rightarrow c+1 \\) \
16:   &emsp; **end while** \
17: &emsp;  \\( i \rightarrow i+1 \\) \
18:  **end while** 
**************************************


Meta-review: All reviewers agree that the paper discusses an interesting and challenging problem and brings new insights and a significant contribution to the problem of learning RNNs based on dynamical systems theory. Also during the discussion phase, the authors were able to provide critical feedback to the comments given by the reviewers, which was helpful in addressing their technical concerns. As a whole, I recommend acceptance (splotlight) of the paper for this submission.
Decision: Accept (spotlight)