Identifying nonlinear dynamical systems with multiple time scales and long-range dependencies
Authors: Dominik Schmidt
Georgia Koppe
Zahra Monfared
Max Beutelspacher
Daniel Durstewitz
Venue: ICLR 2021 Spotlight (2021)
Keywords: nonlinear dynamical systems
recurrent neural networks
attractors
computational neuroscience
vanishing gradient problem
LSTM

Abstract: A main theoretical interest in biology and physics is to identify the nonlinear dynamical system (DS) that generated observed time series. Recurrent Neural Networks (RNN) are, in principle, powerful enough to approximate any underlying DS, but in their vanilla form suffer from the exploding vs. vanishing gradients problem. Previous attempts to alleviate this problem resulted either in more complicated, mathematically less tractable RNN architectures, or strongly limited the dynamical expressiveness of the RNN. 
Here we address this issue by suggesting a simple regularization scheme for vanilla RNN with ReLU activation which enables them to solve long-range dependency problems and express slow time scales, while retaining a simple mathematical structure which makes their DS properties partly analytically accessible. We prove two theorems that establish a tight connection between the regularized RNN dynamics and their gradients, illustrate on DS benchmarks that our regularization approach strongly eases the reconstruction of DS which harbor widely differing time scales, and show that our method is also en par with other long-range architectures like LSTMs on several tasks.

Published as a conference paper at ICLR 2021
IDENTIFYING NONLINEAR DYNAMICAL SYSTEMS WITH MULTIPLE TIME SCALES AND LONG-RANGE DEPENDENCIES
Dominik Schmidt**, Georgia Koppe 2* Zahra Monfared', Max Beutelspacher', Daniel Durstewitz!.3
ABSTRACT
A main theoretical interest in biology and physics is to identify the nonlinear dy- namical system (DS) that generated observed time series. Recurrent Neural Net- works (RNNs) are, in principle, powerful enough to approximate any underlying DS, but in their vanilla form suffer from the exploding vs, vanishing gradients problem. Previous attempts to alleviate this problem resulted either in more com- plicated, mathematically less tractable RNN architectures, or strongly limited the dynamical expressiveness of the RNN. Here we address this issue by suggesting a simple regularization scheme for vanilla RNNs with ReLU activation which en- ables them to solve long-range dependency problems and express slow time scales, while retaining a simple mathematical structure which makes their DS properties partly analytically accessible. We prove two theorems that establish tight con- nection between the regularized RNN dynamics and its gradients, illustrate on DS benchmarks that regularization approach strongly eases the reconstruction of DS which harbor widely differing time scales, and show that our method is also en par with other long-range architectures like LSTMs on several tasks.
INTRODUCTION
Theories in the natural sciences often formulated in terms of sets of stochastic differential or dif- ference equations, i.e. as stochastic dynamical systems (DS). Such systems exhibit a range of com- mon phenomena, like (limit) cycles, chaotic attractors, or specific bifurcations, which are the subject of nonlinear dynamical systems theory (DST; Strogatz (2015); Ott (2002)). A long-standing desire to retrieve the generating dynamical equations directly from observed time series data (Kantz & Schreiber, 2004), and thus to 'automatize' the laborious process of scientific theory building to some degree. A variety of machine and deep learning methodologies toward this goal have been introduced in recent years (Chen et al., 2017: Champion et al., 2019; Ayed et al., 2019; Koppe et al., 2019; Hamilton et al., 2017: Razaghi Paninski, 2019; Hernandez et al., 2020). Often these are based on sufficiently expressive series expansions for approximating the unknown system of genera- tive equations, such as polynomial basis expansions (Brunton et al., 2016: Champion et al., 2019) or recurrent neural networks (RNNs) (Vlachas et al., 2018; Hernandez et al., 2020; Durstewitz, 2017; Koppe et al., 2019). Formally, RNNs are (usually discrete-time) nonlinear DS that are dynamically universal in the sense that they can approximate to arbitrary precision the flow field of any other DS on compact sets of the real space (Funahashi & Nakamura, 1993; Kimura & Nakano, 1998; Hanson & Raginsky, 2020). Hence, RNNs seem like good choice for reconstructing in this sense of dynamically equivalent behavior - the set of governing equations underlying real time series data.
However, RNNs in their vanilla form suffer from the 'vanishing or exploding gradients' problem (Hochreiter & Schmidhuber, 1997; Bengio et al., 1994): During training, error gradients tend to either exponentially explode or decay away across successive time steps, and hence vanilla RNNs face severe problems in capturing long time scales or long-range dependencies in the data. Specially designed RNN architectures equipped with gating mechanisms and linear memory cells have been proposed for mitigating this issue (Hochreiter & Schmidhuber, 1997; Cho et al., 2014). However, from a DST perspective, simpler models that can be more casily analyzed and interpreted in DS
Department of Theoretical Neuroscience, "Clinic Psychiatry and Psychotherapy. Central Institute of Mental Health, Medical Faculty Mannheim, Heidelberg University FFacuity of Physics and Astronomy, Heidelberg University & Bernstein Center Computational Neuroscience
"These authors contributed equally 'Corresponding author: danies e durstewi tzØzi -mannhe de

---

Published as a conference paper at ICLR 2021
terms (Monfared & Durstewitz, 2020a:b), and for which more efficient inference algorithms exist that emphasize approximation of the true underlying DS (Koppe et al., 2019; Hernandez et al., 2020; Zhao & Park, 2020), would be preferable. More recent solutions to the vanishing vs. exploding gradient problem attempt to retain the simplicity of vanilla RNNs by initializing or constraining the recurrent weight matrix to be the identity (Le et al., 2015), orthogonal (Henaff et al., 2016; Helfrich et al., 2018) or unitary (Arjovsky et al., 2016). While merely initialization-bart based solutions, however, may be unstable and quickly dissolve during training, orthogonal or unitary constraints, on the other hand, are too restrictive for reconstructing DS, and more generally from computational perspective as well (Kerg et al., 2019): For instance, neither chaotic behavior (that requires diverging directions) nor multi-stabilitty, that is the coexistence of several distinct attractors, are possible.
Here we therefore suggest a different solution to the problem which takes inspiration from computa- tional neuroscience: Supported by experimental evidence (Daie et al., 2015; Brody et al., 2003), line or plane attractors have been suggested as a dynamical mechanism for maintaining arbitrary infor- mation in working memory (Seung, 1996; Machens et al., 2005), goal-related active form of short- term memory. A line or plane attractor is a continuous set of marginally stable fixed points to which the system's state converges from some neighborhood, while along the line itself there is neither con- nor divergence (Fig. IA). Hence, a line attractor will perform a perfect integration of inputs and re- tain updated states indefinitely, while a slightly detuned line attractor will equip the system with arbi- trarily slow time constants (Fig. IB). This latter configuration has been suggested as a dynamical ba- sis for neural interval timing (Durstewitz, 2003; 2004). The present idea is to exploit this dynamical setup for long short-term memory and arbitrary slow time scales by forcing part of the RNN's sub- space toward plane (line) attractor configuration through specifically designed regularization terms.
Specifically, our goal here is not so much to beat the state of the art on long short-term memory tasks, but rather to address the exploding vs. vanishing gradient problem within a simple, dynamically tractable RNN, optimized for DS reconstruction and interpretation. For this we build on piecewise- linear RNNs (PIRNNs) (Koppe et al., 2019; Monfared & Durstewitz, 2020b) which employ ReI.U activation functions. PLRNNS have simple mathematical structure (see eq. 1) which makes them dynamically interpretable in the sense that many geometric properties of the system's state space can in principle be computed analytically, including fixed points, cycles, and their stability (Suppl. 6.1.2; Koppe et al. (2019); Monfared & Durstewitz (2020a)), i.e. do not require numerical techniques (Sussillo & Barak, 2013). Moreover, PLRNNS constitute type of piecewise linear (PWL) map for which many important bifurcations have been comparatively well characterized (Monfared & Durstewitz, 2020a; Avrutin al., 2019). PLRNNS can furthermore be translated into equivalent continuous time ordinary differential equation (ODE) systems (Monfared & Durstewitz, 2020b) which comes with further advantages for analysis, e.g. continuous flow fields (Fig. IA.B).
We retain the PLRNN's structural simplicity and analytical tractability while mitigating the explod- ing vs. vanishing gradient problem by adding special regularization terms for subset of PLRNN units to the loss function. These terms are designed to push the system toward line attractor config- urations, without strictly enforeing them, along some - but not all directions in state space. We further establish tight mathematical relationship between the PLRNN dynamics and the behay- ior its gradients during training. Finally, we demonstrate that our approach outperforms LSTM and other, initialization-based. methods on number of 'classical' machine learning benchmarks (Hochreiter & Schmidhuber, 1997). Much more importantly in the present DST context, we demon- strate that our new regularization-suppore inference efficiently captures all relevant time scales when reconstructing challenging nonlinear DS with multiple short- and long-range phenomena.
RELATED WORK
Dynamical systems reconstruction. From a natural science perspective, the goal of reconstructing or identifying the underlying DS is substantially more ambitious than (and different from) building system that 'merely' yields good ahead predictions: In DS identification we require that the inferred model can freely reproduce (when no longer guided by the data) the underlying attractor geometries and state space properties (see section 3.5, Fig. S2; Kantz& & Schreiber (2004)).
Earlier work using RNNs for DS reconstruction (Roweis & Ghahramani, 2002; Yu et al., 2005) mainly focused on inferring the posterior over latent trajectories Z = { } given time se- ries data X {x1,-..,, BT), p( Z x), and on ahead predictions (Iu et al., 2017), as does much of
2

---

Published as a conference paper at ICLR 2021
A
B
Bu=l A
inputs
67
22
add +
Ca=1
W1y11
idx
.,E(0.3)
G
hy=-1
21
Figure 1: A)-B): Illustration of the state space of a 2-unit RNN with flow field (grey) and nullcfines (set of points at which the flow of one of the variables vanishes, in blue and red). Insets: Time graphs of 21 for T = 30 000. A) Perfect line attractor. The flow converges to the line attractor, thus retaining states indefinitely in the absence of perturbations, as illustrated for 3 example trajectories (green). B) Slightly detuned line attractor. The system's state still converges toward the "attractor ghost", but then very slowly crawls up within the 'attractor tunnel' (green trajectory) until it hits the stable fixed point at the intersection of nullclines. Within the tunnel, flow velocity is smoothly regulated by the gap between nullclines, thus enabling arbitrary time constants. C) Simple 2-unit solution to the addition problem exploiting the line attractor properties of ReLUUs. The output unit serves as a perfect integrator (see Suppl. 6.1.1 for complete parameters).
the recent work on variational inference of DS (Duncker et al., 2019; Zhao & Park, 2020; Hernan- dez et al., 2020). Although this enables insight into the dynamics along the empirically observed trajectories, both - posterior inference and good ahead predictions - do not per se guarantee that the inferred models can generate the underlying attractor geometries on their own (see Fig, S2, Koppe et al. (2019)). In contrast, if fully generative reconstruction of the underlying DS in this latter sense were achieved, formal analysis simulation of the resulting RNN equations could provide much deeper understanding of the dynamical mechanisms underlying empirical observations (Fig. 1 C).
Some approaches geared toward this latter goal of full DS reconstruction make specifie structural assumptions about the form of the DS equations ( ('white box approach'; Meeds et al. (2019); Raissi (2018); Gorbach et al. (2017)), e.g. based on physical or biological domain knowledge, and focus on estimating the system's 's latent states and parameters, rather than approximating an unknown DS based on the observed time series information alone ('black box approach'). Others (Trischler & D' D'Eleuterio, 2016; Brunton et al., 2016; Champion et al., 2019) attempt to approximate the flow field, obtained e.g. by numerical differentiation, directly through basis expansions or neural net- works. However, numerical derivatives are problematic for their high variance and other numerical issues (Raissi, 2018; Baydin et al., 2018; Chen et al., 2017). Another factor to consider is that in many biological systems like the brain the intrinsic dynamics are highly stochastic with many noise sources, like probabilistic synaptic release (Stevens, 2003). Models that do not explicitly account for dynamical process noise (Ayed et al., 2019; Champion et al., 2019; Rudy et al., 2019) are therefore less suited and more vulnerable to model misspecificcatinn. Finally, some fully probabilistic models for DS reconstruction based on GRU (Fraccaro et al., 2016), LSTM (Zheng al., 2017; Vlachas et al., 2018), or radial basis function (Zhao & Park, 2020) networks, are not easily interpretable and amenable to DS analysis in the sense defined in sect. 3.3. Most importantly, none of these previous approaches consider the long-range dependency problem within more easily tractable RNNs for DS.
Long-range dependency problems in RNNs. Error gradients in vanilla RNNs tend to either explode or vanish due to the large product of derivative terms that results from recursive application of the chain rule over time steps (Hochreiter, 1991; Bengio et al., 1994: Hochreiter Schmidhuber, 1997). To address this issue, RNNs with gated memory cells (Hochreiter & Schmidhuber, 1997; Cho et al., 2014) have been specifically designed, but their more complicated mathematical structure makes them less amenable to systematic DS analysis. Even simple objects like fixed points of these sys- tems have to be found by numerical techniques (Sussillo & Barak, 2013; Jordan et al., 2019). Thus, approaches which retain the simplicity of vanilla RNNs while solving the exploding vs. vanishing gradients problem would be desirable. Recently, Le et al. (2015) observed that initialization of the recurrent weight matrix W to the identity in ReLU-based RNNs may yield performance par with LSTMs on standard machine learning benchmarks. Talathi & Vartak (2016) expanded on this idea by initializing the recurrence matrix such that its largest absolute eigenvalue is 1. Later work en-
3

---

Published as a conference paper at ICLR 2021
forced orthogonal (Henaff al., 2016; Helfrich et al., 2018; Jing et al., 2019) unitary (Arjovsky et al., 2016) constraints on the recurrent weight matrix during training. While this appears to yield long-term memory performance sometimes superior to that of LSTMs (but see (Henaff et al., 2016)), these networks are limited in their computational power (Kerg et al., 2019). This may be a conse- quence of the fact that RNNs with orthogonal recurrence matrix are quite restricted in the range of dynamical phenomena they can produce, c.g. chaotic attractors are not possible since (locally) diverging -igen-directions are disabled.
Our approach therefore is to establish line/plane attractors only along some but not all directions in state space, and to only push the RNN toward these configurations but not strictly enforce them, such that convergence or (local) divergence of RNN dynamics is still possible. We furthermore im- plement these concepts through regularization terms in the loss functions, rather than through mere initialization. This way plane attractors are encouraged throughout training without fading away.
MODEL FORMULATION AND THEORETICAL ANALYSIS
3.1 BASIC MODEL FORMULATION
Assume we are given two multivariate time series S = { and = one we will denote as 'imputs' (S) and the other as 'outputs' (X). the 'classical' (supervised) machine learning setting, we usually wish to map S on X through a RNN with latent state equation 2₁ Fe (2,-1, 8t) and outputs I ~ Px (z)) as for instance in the 'addition problem' (Hochreiter & Schmidhuber, 1997). In DS reconstruction, in contrast, usually have dense time series X from which we wish to infer (unsupervised) the underlying DS, where S may provide an additional forcing function or sparse experimental inputs or perturbations. While our focus in this paper is on this latter task, DS reconstruction, we will demonstrate that our approach brings benefits in both these settings.
Here we consider for the latent model PLRNN (Koppe et al.. 2019) which takes the form
2, Az,1 W (Zz2-1) Cs, h + Et: e NO, E),
(1)
where 2, E Mx1 is the hidden state (column) vector of dimension M, A E M a diagonal and W R MMM an off-diagonal matrix, 8, E Rx1 the external input of dimension K, C E R Mxx the input mapping, h E x1 bias, and E, a Gaussian noise term with diagonal covariance matrix diag (E) E R. The nonlinearity Φ0z) is ReLU, Φ(z); = max(0, 2), E {l..... M). This specific formulation represents a discrete-time time version of firing rate (population) models as used in computational neuroscience (Song et al., 2016; Durstewitz, 2017: Engelken et al., 2020).
We will assume that the latent RNN states are coupled to the actual observations 2, through a simple observation model of the form
Dt Bg(z)) + nu: n. N)
in the case of observations F; Nx1 where B E RNM is factor loading matrix, 9 some (usu- ally monotonic) nonlinear transfer function (e.g., ReLU),, and diag(T) R. the diagonal covari- ance matrix of the Gaussian observation noise, or through a softmax function in case of categorical observations Ti,t E {0, 1} (see Suppl. 6.7.7 for details).
3.2 REGULARIZATION APPROACII
First note that by letting A I, W = 0, and h 0 in eq. 1, every point in space will be marginally stable fixed point of the system, leading it to perform a perfect integration of external inputs in parametric working memory (Machens et al., 2005; Brody et al., 2003).' This is similar in spirit to Le et al. (2015) who initialized RNN parameters such that it performs an identity mapping neuroscientifially motivated network architecture (eq. 1) that enables the identity mapping across the variables' entire support, 2i,t E [-00, +oo]. which we conjecture will be of advantage for establishing long short-term memory properties, 2) we encourage
for zi,t > 0. However, here 1)
Note that this very property of marginal stability required for input integration also makes the system sensitive to noise perturhations directly the manifold attractor. Interestingly, this property has indeed been observed experimentally for real neural integrator systems (Major al., 2004; Mizumori & Williams, 1993).
4

---

Published as a conference paper at ICLR 2021
this mapping only for subset Mreg M of units (Fig. S1), leaving others free to perform arbitrary computations, and 3) we stabilize this configuration throughout training by introducing a specific L2 regularization for parameters A, W, and h in eq. 1. When embedded into larger, (locally) convergent system, we will call this configuration more generally a manifold attractor.
That way, divide the units into two types, where the regularized units serve as memory that tends to decay very slowly (depending on the size of the regularization term), while the remaining units maintain the flexibility approximate any underlying DS, yet retaining the simplicity of the original PLRNN (eq. 1). Specifically, the following penalty is added to the loss function (Fig. S1);
Meng
Myng
Myng
Lreg TA Σ (Ar,- 1)2 TW 1=1
W
h
J
(Recall from sect. 3.1 that A is diagonal and W is an off-diagonal matrix.) While this formulation allows us to trade off, for instance, the tendency toward manifold attractor (A - I, h. - 0) vs. the sensitivity to other units' inputs (W - 0), for all experiments performed here a common value, TA TW Th was assumed for the three regularization factors. We will refer to (21 ... .M.a)) as the regularized ('memory') subsystem, and to (1 .M. as the non-rege non-regularized ( ('compu- tational') subsystem. Note that in the limit oo exact manifold attractors would be enforced.
3.3 THEORETICAL ANALYSIS
We will now establish a tight connection between the PLRNN dynamics and its error gradients. Similar ideas appeared in Chang et (2019), but these authors focused only on fixed point dynam- ies, while here we will consider the more general case including cycles of any order. First, note that by interpretabillty of model eq. 1 we mean that it easily amenable a rigorous DS analysis: As shown in Suppl. 6.2.2, we can explicitly determine all the system's fixed points and cycles and their stability. Moreover, as shown in Monfared & Durstewitz (2020b). we can under certain conditions transform the PLRNN into an equivalent continuous-tima (ODE) piecewisc-lings system, which brings further advantages for DS analysis.
Let us rewrite eq. in the form
21 = F(x-1) = (A W DoX--1))
WR(--1) 21-1 + h.
(4)
where D 1) is the diagonal matrix of outer derivatives of the ReLU function evaluated at # (see Suppl. 6. 6.1) and we ignore external inputs and noise terms for now. Starting from some initial condition 21, we can recursively develop 2r as (see Suppl. 6.1.2 for more details):
1-1
ZT FT = II W Wor-- 21 +
Σ
II
W a(T--) + 11]n h.
(5)
Likewise, for some common loss function L(A, W. h) - -t Lr, we can recursively develop the derivatives w.r.t. weights Wink (and similar for components of A and h) as
ac
T OL. Ōz.
Oz. 1 (m. DD(1-1) 2t-1
with Ow
(6)
Owin
Oz
Owik
Oz2
Σ II E=1
Wa(co) ((m,k) Douc-))
II W Wo(t-)) Dwink
where (m,,) an M M indicator matrix with 1 for the (m, k)'th entry and 0 everywhere else. Observing that eqs. 5 and 6 contain similar product terms which determine the system's 's long-term behavior, our first theorem links the PLRNN dynamics to its total error gradients:
Theorem 1. Consider a PLRNN given by eq. 4, and assume that it converges to a stable fixed point, say 2 := =1, or a k-cycle (k 1) with the periodic points [Zy.k, Zirk + oc. Suppose that, 2 and E {0, 10.1..... where Wi) denotes the Jacobian of the system at z and maaz indicates the largest
for 1), o omanl W W = Ware*=
5

---

Published as a conference paper at ICLR 2021
singular value of a matrix. Then, the 2-norms of the tensors collecting all derivatives, Azm
der @aW 1
Bxr
12- will be bounded above, i.e. will not diverge for - 00.
BA
2'
ah
Proof.
See Suppl. sect. 6.1 (subsection 6.1.3).
While Theorem is a general statement about PLRNN dynamics and total gradients, our next the- orem more specifically provides conditions under which Jacobians linking temporally distant states 2T and 2t. 7 > I, will neither vanish nor explode in the regularized PLRNN:
Theorem 2. Assume a PLRNN with matrix A + W partitioned as in Fig. S1, i.e. with the first Mreg rows corresponding to those of an M x M identity matrix. Suppose that the son-regularized subsystem (2Mreg-1....2 .MM, if considered in isolation, satisfies Theorem 1, i.e. converges to a k- cycle with k 2 1. Then, for the full system (21 (2M).. the 2-norm of the Jucobians connecting temporally distal states and will be bounded from above and below for all T > t, L.e. oc >
Pup 2
8zy = HIEAST Wa(s] 2 2 Plow > 0. In particular, for state variables Zer and 0=
2
i.e. that connect states from the 'memory' to those of the 'computational subsystem, one also has do >> Aup 2 as T - oc, i.e. these derivatives will never vanish nor explode.
Zzt such that i E {Mreg + +,...,,M , and j E
Basg Bays 2 Noow >
Proof. See Suppl. sect. 6.1 (subsection 6.1.4).
The bounds Pup. Plow. Aups Ntow given in Suppl. sect. 6.1.4. We remark that when the regular- ization conditions are not exactly met, i.c. when parameters A and W slightly deviate from those in Fig. S1, memory (and gradients) may ultimately dissipate, but only very slowly, actually required for temporal processes with very slow yet not infinite time constants (Fig. IB).
3.4 TRAINING PROCEDURES
For the (supervised) machine learning problems, all networks were trained by stochastic gradient descent (SGD) to minimize the squared-error loss between estimated and actual outputs for the addition and multiplication problems, and the cross entropy loss for sequential MNIST (see Suppl. 6.1.7). Adam (Kingma & Ba, 2014) from PyTorch package (Paszke et al., 2017) was used as the optimizer, with learning rate of 0.001, gradient clip parameter of 10, and batch size of 500. SGD was stopped after 100 epochs and the fit with the lowest loss across all epochs was taken, except for LSTM which was allowed run for up to 200 epochs as it took longer converge (Fig. S10). For comparability, the PLRNN latent state dynamics eq. was assumed to be deterministic in this IN in eq. For the regularized PLRNN (rPLRNN),
setting (i.e., Σ 0), g(z)) 2, and penalty eq. 3 was added to the loss function. For the (unsupervised) DS reconstruction problems, the fully probabilistic, generative RNN eq. I was considered. Together with eq. 2 (where we take g(z;) o(2,)) this gives the typical form of nonlinear state space model (Durbin & Koopman, 2012) with observation and process noise, and an Dnpectation-Maxiiiaane efficiently exploits the model's piecewise linear structure (Durstewitz, 2017; Koppe et al., 2019) was used to solve for the parameters by maximum likelihood. Details are given in Suppl. 6.5.5. All code used here will be made openly available at https:/ggihub.com/DDutteteee wtwp://wwww..a/gg........
(EM) algorithm that
3.5 PERFORMANCE MEASURES
For the machine learning benchmarks we employed the same criteria as used for optimization (MSE or cross-entropy, Suppl. 6.1.7) performance metrics, evaluated across left-out test sets. In ad- dition, we report the relative frequency Peorrect of correctly predicted trials across the test set (see Suppl. 6.1.7 for details). For DS reconstruction problems, is not sufficient or even sensible to judge a method's 's ability to infer the underlying DS purely based on some form of (ahead- spredicciiniiiii error like the MSE defined on the time series itself (Ch.12 12 in Kantz. & Schreiber (2004)). Rather, we require that the inferred model can freely reproduce (when no longer guided by the data) the un- derlying attractor geometries and state space properties. This is not automatically guaranteed for a model that yields agrecable ahead predictions Wood (2010)). We therefore followed Koppe et al. (2019) and used the between true and reproduced probability distributions across states in state space to quantify how well an inferred PLRNN captured the underlying dynamics, thus assessing the agreement in attractor geometries (cf. Takens (1981); Sauer et al. (1991)) (see Suppl. 6.1.6 for more details).
time series (Fig. S2A; cf. Koppe et al. (2019);
divergence

---

Published as a conference paper at ICLR 2021
A
Addition problem
Multiplication problem
Sequential MNIST
10
10
CDUS endrupy luss
10
10
10
10
10
10-3
10-1
10-7
10
20
50.70
150 250 100
20
30
50
70
150
250 400
900
RANN
IRMW
upRNV
oRVV
LSTM
PLRAN
IPLRAV PPLRN LaPPLRNN
KNN
upRNNN (Tolothi Vastak, 2016) LSTM (Rochreiter Schmichhhher 1997)
IPLKNN
IRNN (Le nl., 2015)
12pPLRNN
(Norontaov
PLRNN (Kogge 2019)
PLLRNN
0.8
0.8
0.8
Porrett
0..6
0.66
0.1
0.1
0.4
0.2
RAW IRVV
apRNN aRVA
LSTM
PLRAN IPLRAV *PLRNN
L2pPLRN
Comparison of FPLRNN (r Mass 0.5, cf. Fig. S3) to other methods for A) addition problem, B) multiplication problem and C) sequential MNIST. Top row gives loss as function of time series length T (error bars SEM, n 2 5), bottom row shows relative frequency of correct trials. Note that better performance (lower values in top row, higher values in bottom row) reflected more rightward shift of curves. Dashed lines indicate chance level, black dots in indicate individual repetitions.
Figure 2:
B
C
1212
16
normalliod MSR
50-500
14
I
0.015
0.015
MSE
12.
0.01
0.01
0.005
0.006
10
1el
1d2
1c0
164
165
1cl
1x2
160
1ok
168
1el
162
163
1c5
D
E
Lrue
rogulzsteed
generated
latent states lm
0
22
0.15
0.5
0.45
0.15
0.45
0.8
0.15
0.3
0.65
0.6
time (6)
time
time
Figure Reconstruction of 2-time scale DS in limit cycle regime. A) KL divergence (DKL) between true and generated state space distributions. Globally diverging system estimates were re- moved. Average MSE between power spectra true and reconstructed DS and C) split according to low (S 50 Hz) and high (> 50 Hz) frequency components. Error bars SEM (n = 33). D) Ex- ample (best) generated time series (red=reconstruction with t = 3) See Fig. S5A for variable n. E) Dynamics of regularized and non-regularized latent states for the example in D.
NUMERICAL. EXPERIMENTS
4.1 MACHINE LEARNING BENCHMARKS
Although not our prime interest here, we first examined how the rPLRNN would fare on supervised machine learning benchmarks where inputs (S) are to be mapped onto target outputs (X) across long

---

Published as a conference paper at ICLR 2021
time spans (i.e., requiring long short-term maintenance of information), namely the addition and multiplication problems (Talathi & Vartak, 2016; Hochreiter & Schmidhuber, 1997), and sequential MNIST (LeCun et al., 2010). Details of these experimental setups are in Suppl. 6. Performance of the FPLRNN (eq. 1, eq. 3) on all 3 benchmarks was compared to several other models summarized in Suppl. Table 1. To achieve meaningful comparison, all models have the same number M = 40 (based on Fig. S3) of hidden states (which gives LSTMs overall about 4 times as many trainable parameters). On all three problems the rPLRNN outperforms all other tested methods, including LSTM, iRNN (RNN initialized by the identity matrix as in Le et al. (2015)), and a version of the orthogonal RNN (oRNN; Vorontsov et al. (2017)) (similar results were obtained for other settings of M and batch size). LSTM performs even worse than iRNN and iPLRNN (PLRNN initialized with the identity as the iRNN), although had 4 times as many parameters and was given twice as many epochs (and thus opportunities) for training, as it also took longer to converge (Fig. S10). In addition, the iPLRNN tends to perform slightly better than the IRNN on all three problems, suggesting that the specific structure eq. 1 of the PLRNN that allows for manifold attractor across the variables' full range may be advantageous to begin with, while the regularization further improves performance.
4.2 NUMERICAL EXPERIMENTS ON DYNAMICAL SYSTEMS WITH DIFFERENT TIME SCALES
While it is encouraging that the rPLRNN may perform even better than several previous approaches to the vanishing vs. exploding gradients problem, our major goal here to examine whether our regularization scheme would help with the (unsupervised) identification of DS that harbor widely different time scales. To test this, we used biophysical, bursting cortical neuron model with one voltage (V) and two conductance recovery variables (see Durstewitz (2009),, one slow (h) and one fast (n; Suppl. 6.1.8) Reproduction of this DS is challenging since it produces very fast spikes on top of a slow nonlinear oscillation (Fig. 3D). Only short time series (as in scientific data) of length T = 1500 from this model were provided for training. rPLRNNs with M = {8. .... 18} states were trained, with the regularization factor varied within T E {0. 10', 10?, 10², 10 100110"//T. Note that for r = 0 (no regularization), the approach reduces to the standard PLRNN (Koppe et al., 2019).
Fig. 3A confirms our intuition that stronger regularization leads to better DS reconstruction as assessed by the KL. divergence between true and generated state distributions (similar results were obtained with ahead-prediction
errors as metric, Fig. S4A), accompanied by a likewise decrease in the MSE between the power spectra of true (suppl. eq. 55) and generated (rPLRNN) voltage traces (Fig. 3B). Fig. 3D gives example of voltage traces (V) and the slower of the two gating variables (h; see Fig. S5A for variable n) freely simulated (i.e., sampled) from the autonomously running rPLRNN. This illustrates that our model is in principle capable of capturing both the stiff spike dynamics and the slower oscillations in the second gating variable at the same time. Fig. 3C provides more insight into how the regularization worked: While the high frequency components (> 50 Hz) related to the repetitive spiking activity hardly benefited from increasing there was a strong reduction in the MSE computed on the power spectrum for the lower frequency range (< 50 Hz), suggesting that increased regularization helps to map slowly evolving components of the dynamics. This result more general as shown in Fig. S6 for another DS example. In contrast, an orthogonality (Vorontsov et al., 2017) or plain L2 constraint on weight matrices did not help at all on this problem (Fig. S4B).
Further insight into the dynamical mechanisms by which the rPLRNN solves the problem can be obtained by examining the latent dynamies: As shown in Fig. 3E (see also Fig. S5), regularized states indeed help map the slow components of the dynamics, while non-regularized states focus on the fast spikes. These observations further corroborate the findings in Fig. 3C and Fig. S6C.
4.3 REGULARIZATION PROPERTIES AND MANIFOLD ATTRACTORS
In Figs. 2 and 3 we demonstrated that the rPLRNN able to solve problems and reconstruct dy- namics that involve long-range dependencies. Figs. 3A.B furthermore directly confirm that solutions improve with stronger regularization, while Figs. 3C,E give insight into the mechanism by which the regularization works. To further verify empirically that our specific form of regularization, eq. 3, is important, Fig. 2 also shows results for a PLRNN with standard L2 norm on a fraction of Mreg/M 0.5 states (L2pPLRNN). Fig. S7 provides additional results for PLRNNs with L2 norm on all weights and for vanilla L2-regularized RNNs. All these systems fell far behind the perfor- mance of the rPLRNN all tasks tested. Moreover, Fig. 4 reveals that the specific regularization
8

---

Published as a conference paper at ICLR 2021
54
B 0A
indetiee Exequnes
0.3
TayeTa00 1a0
estabox
mand]]Meee 11
rlative Enconacy
02
0.2
9
84
1.2
1:5p
3el 160
165
1.5
1.25 1T
Dloe around F1s
1500
Die assudd ETs
Figure 4: A) Distribution of maximum absolute eigenvalues À of Jacobians around fixed points for iPLRNN for different T and L2PLRNN trained on bursting neuron DS. B) Absolute deviations of max. [^] from 1 (using for each system the one eigenvalue with smallest deviation). C) Same as for addition problem for rPLRNN (T = 5) vs. standard, fully 1.2- (1.2f), and partially 1.2 (1.2p)- regularized PLRNN. D) Same as B for the models from C. Error bars stdv. See also Fig. S8.
proposed indeed encourages manifold attractors, and that this is not achieved by a standard L2 reg- ularization: In contrast to L2PLRNN, the regularization factor is increased, more and more of the maximum absolute eigenvalues around the system's fixed points (computed according to eq. 8, sect. 6.1.2) cluster on near 1, indicating directions of marginal stability in state space. Also, the deviations from I become smaller for strongly regularized PLRNNs (Fig. 4B,D), indicating a higher precision in attractor tuning. Fig. S9 in addition confirms that FPLRNN parameters are increasingly driven toward values that would support manifold attractors with stronger regularization. Fig. 3E furthermore suggests that both regularized and non-regularized states are utilized to map the full dynamics. But how should the ratio Mreg/M M be chosen in practice? While for the problems here this meta-parann pera-parrrettr was determined through 'classical' grid-search and cross-validation, Figs. S3 C- E suggest that the precise setting of Mreg/M is actually not overly important: Nearly optimal performance achieved for a broader range Mreg/M /M e [0.3, 0.6] on all problems tested. Hence, in practice, setting M /M 0.5 should mostly work fine.
CONCLUSIONS
In this work we introduced a simple solution to the long short-term memory problem in RNNs that retains the simplicity and tractability of PLRNNs, yet does not curtail their universal computational capabilities (Koiran et al., 1994; Siegelmann & Sontag, 1995) and their ability to approximate arbi- trary DS (Funahashi & Nakamura, 1993; Kimura & Nakano, 1998; Trischler & D D'Eleuterio, 2016). We achieved this by adding regularization terms to the loss function that encourage the system to form a 'memory subspace' (Seung, 1996; Durstewitz, 2003) which would store arbitrary values for, if unperturbed, arbitrarily long periods. At the same time we did not rigorously enforce this con- straint, which allowed the system to capture slow time scales by slightly departing from a perfect manifold attractor. In neuroscience, this has been diseussed
dynamical mechanism for regu- lating the speed of flow in DS and learning of arbitrary time constants not naturally included qua RNN design (Durstewitz, 2003; 2004) (Fig. IB). While other RNN architectures, including vanilla RNNs, can, in principle, also develop line attractors to solve specific tasks (Maheswaranathan et al., 2019), they are generally much harder to train to achieve this and may exhibit less precise attractor tuning (cf. Fig. 4), which is needed to bridge long time scales (Durstewitz, 2003). Moreover, part of PLRNN's latent space was not regularized at all, leaving the system enough degrees of freedom for realizing arbitrary computations or dynamics (see also Fig. S11 for chaotic example). We showed that the rPLRNN is en par with or outperforms initialization-base approaches, orthogonal RNNs, and LSTMs on number of classical benchmarks. More importantly, however, the regular- ization strongly facilitates the identification of challenging DS with widely different time scales in PLRNN- based algorithms for DS reconstruction. Similar regularization schemes as proposed here (eq. 3) may, in principle, also be designed for other architectures, but the convenient mathematical form of the PLRNN makes their implementation particularly powerful and straightforward.
ACKNOWLEDGEMENTS
This work was funded by grants from the German Research Foundation (DFG) to DD (Du 354/10-1, Du 354/8-2 within SPP 1665) and to GK (TRR265: A06 & B08), and under Germany's 's Excellence Strategy EXC-2181 390900948 ("Structures').

---

Published as a conference paper at ICLR 2021
REFERENCES
Martin Arjovsky, Amar Shah, and Yoshua Bengio. Unitary evolution recurrent neural networks. In Proceedings of the 33rd International Conference on Machine Learning, 2016.
Viktor Avrutin, Laura Gardini, Iryna Sushko, and Fabio Tramontana. Continuous and Discontimuous Srecewise-Smooth One- Ohe-Dimensional Maps. WORLD SCIENTIFIC, 2019.
Ibrahim Ayed, Emmanuel de Bézenac, Arthur Pajot, Julien Brajard, and Patrick Gallinari. Learning Dynamical Systems from Partial Observations. arXiv preprint, 2019.
Atlım Güneş Baydin, Barak Pearlmutter, Alexey Andreyevich Radul, and Jeffrey Mark Siskind. Automatic differentiation in machine learning: A survey. Journal of Machine Learning Research, 18:1-43, 2018
Yoshua Bengio, Patrice Simard, and Paolo Frasconi. Learning long-term dependencies with gradient descent is difficult. IEEE Transactions on Neural Networks, 5(2):157-166, 1994.
Carlos D Brody, Ranulfo Romo, and Adam Kepecs. Basic mechanisms for graded persistent activ- ity: discrete attractors, continuous attractors, and dynamic representations. Current Opinion in Neurobiology, 13(2):204-211 211, 2003.
Steven L. Brunton, Joshua L. Proctor, and J. Nathan Kutz. Discovering governing equations from data by sparse identification of nonlinear dynamical systems. Proceedings of the National Academy of Sciences of the United States of America, 1 113(55)333223937, 2016.
Kathleen Champion, Bethany Lusch, Nathan Kutz, and Steven L. Brunton. Data-driven discovery of coordinates and governing equations. arXiv preprint, 2019.
Bo Chang, Minmin Chen, Eldad Haber, and Ed H. Chi. AntisymmetricRNN: A Dynamical System View on Recurrent Neural Networks. arXiv preprint, 2019.
Shizhe Chen, Ali Shojaie, and Daniela M. Witten. Network Reconstruction From Nigh-Dimensional Ordinary Differential Equations. Journal of the American Statistical Association, 1 112(520): 1697- 1707, 2017.
Kyunghyun Cho, Bart Van Merriënboer, Caglar Gulcehre, Dzmitry Bahdanau, Fethi Bougares, Hol- ger Schwenk, and Yoshua Bengio. Learning phrase representations using RNN encoder-decoder for statistical machine translation. Proceedings of the Conference on Empirical Methods in Natural Language Processing, pp. 1724-1734. Association for Computational Linguistics (ACL), 2014.
Kayvon Daie, Mark S. Goldman, and Emre R.F. Aksay. Spatial patterns of persistent neural activity vary with the behavioral context of short-term memory. Neuron, 8544)8847-860 2015.
Lea Duncker, Gergo Bohner, Julien Boussard, and Maneesh Sahani. Learning interpretable continuous-time models of latent stochastic dynamical systems. arXiv preprint, 2019.
James Durbin and Siem Jan Koopman. Time Series Analysis by State Space Methods. Oxford
University Press, 2012.
Daniel Durstewitz. Self-Organd Organizinn Neural Integrator Predicts Interval Times. Journal of Neuro- science, 23(12) 2003.
Daniel Durstewitz. Neural representation of interval time. NewroReport, 15(5):745-749, 2004.
Daniel Durstewitz. Implications of synaptic biophysics for recurrent network dynamics and active memory. Neural Networks, 2288):1189-1200, 2009.
Daniel Durstewitz. A State Space Approach for Piecewise-Lincar Recurrent Neural Networks for Reconstructing Nonlinear Dynamics from Neural Measurements. PLoS Computational Biology, 130554 2017.
Rainer Engelken. Fred Wolf, and L. F. Abbott. Lyapumov spectra of chaotic recurrent neural net- works. arXiv preprint, 2020.
10

---

Published as a conference paper at ICLR 2021
Ludwig Fahrmeir and Gerhard Tutz. Multivariate Statistical Modelling Based Generalized Linear Models. Springer New York, 2001.
Jianqing Fan and Qiwei Yao. Nonlinear Time Series. Springer New York, 2003.
Marco Fraccaro, Søren Kaae Sønderby, Ulrich Paquet, and Ole Winther. Sequential neural models with stochastic layers. In Advances in Neural Information Processing Systems, 2016.
Ken-ichi Funahashi and Yuichi Nakamura. Approximation of dynamical systems by continuous time recurrent neural networks. Neural Networks, 6(6):801-806, 1993.
Nico S. Gorbach, Stefan Bauer, and Joachim M. Buhmann. Scalable variational inference for dy- namical systems. In Advances Neural Information Processing Systems, 2017.
Franz Hamilton, Alun L. Lloyd, and Kevin B. Flores. Hybrid modeling and prediction of dynamical systems. PLOS Computational Biology, 13(7):e 2017.
Joshua Hanson and Maxim Raginsky. Universal simulation of stable dynamical systems by recurrent neural nets. In Proceedings of Machine Learning Research, volume 120, pp. 384-392, 2020.
Kyle E. Helfrich, Devin Whimott, and Qiang Ye. Orthogonal recurrent neural networks with scaled Cayley transform. Proceedings of the 35th International Conference on Machine Learning, 2018.
Mikael Henaff, Arthur Szlam, and Yann LeCun. Recurrent orthogonal networks and long- memoor tasks. In Proceedings of the 33rd International Conference on Machine Learning, 2016.
Daniel Hernandez, Antonio Khalil Moretti, Ziqiang Wei, Shreya Saxena, John Cunningham, and Liam Paninski. Nonlinear evolution via spatially-dependent linear dynamics for electrophysiol- ogy and calcium data. NBDT, 3(3), 2020.
John R. Hershey and Peder A. Olsen. Approximating the Kullback Leibler divergence between gaussian mixture models. In 2007 IEEE International Conference on Acoustics, Speech and Signal Processing, volume 4, pp. IV-317-1V-320, 2007.
Sepp Hochreiter. Untersuchungen zu dynamischen neuronalen Netzen. Diploma thesis, Institut für Informatik, Lehrstuhl Prof. Brauer, Technische Universität München, 1991.
Sepp Hochreiter and Jürgen Schmidhuber. Long short-term memory. Neural computation, 9(8): 1735 80, 1997.
Li Jing, Caglar Gulcehre, John Peurifoy, Yichen Shen, Max Tegmark, Marin Soljacic, and Yoshua Bengio. Gated Orthogonal Recurrent Units: On Learning to Forget. Neural Computation, 31()): 765-783, 2019.
lan D. Jordan, Piotr A. Sokol, and II Memming Park. Gated recurrent units viewed through the lens of continuous time dynamical systems. arXiv preprint, 2019.
Holger Kantz and Thomas Schreiber. Nonlinear Time Series Analysis. Cambridge University Press, 2. edition, 2004.
Giancario Kerg, Kyle Goyette, Maximilian P. Touzel, Gauthier Gidel, Eugene Vorontsov, Yoshua Bengio, and Guillaume Lajoie. Non-normal Recurrent Neural Network (nnRNN): learning long time dependencies while improving expressivity with transient dynamics. arXiv preprint, 2019.
M. Kimura and R. Nakano. Learning dynamical systems by recurrent neural networks from orbits. Newral Networks, 1199):15891159) 1998.
Diederik P. Kingma and Jimmy Ba. Adam: Method for Stochastic Optimization. arXiv preprint, 2014.
Pascal Koiran, Michel Cosnard, and Max H. Garzon. Computability with low-dimensional dynami- cal systems. Theoretical Computer Science, 132:1 13-128, 1994.

---

Published as a conference paper at ICLR 2021
Georgia Koppe, Hazem Toutounji, Peter Kirsch, Stefanie Lis, and Daniel Durstewitz. Identifying nonlinear dynamical systems via generative recurrent neural networks with applications to fMRI. PLOS Computational Biology, 15(8)=e 2019.
Quoc V. Le, Navdeep Jaitly, and Geoffrey E. Hinton. A Simple Way to Initialize Recurrent Networks of Rectified Linear Units. arXiv preprint, 2015.
Yann LeCun, Corinna Cortes, and CJ Burges. Mnist handwritten digit database. ATT Labs [Online). 2010.
Available:
Zhixin Lu, Jaideep Pathak. Brian Hunt, Michelle Girvan, Roger Brockett, and Edward Ott. Reservoir observers: Model-free inference of unmeasured variables in chaotic systems. Chaos, 27)4), 2017.
Christian K. Machens, Ranulfo Romo, and Carlos D. Brody. Flexible control of mutual inhibition: A neural model of two-interval discrimination. Science, 307(5712): 2005. doi: 10. 1166/ccience. 104171.
Niru Maheswaranathan, Alex Williams, Matthew Golub, Surya Ganguli, and David Sussillo. Re- verse engineering recurrent networks for sentiment classification reveals line attractor dynamics. In H. Wallach, H. Larochelle, A. Beygelzimer, F. d'Alché-Buc, E. Fox, and R. Garnett (eds.), Advances in Neural Information Processing Systems, volume 32, pp. 15696-15705. Curran As- sociates, Inc., 2019.
G. Major, R. Baker, E. Aksay, B. Mensh, H. Seung, and D. W. Tank. Plasticity and tuning by visual feedback of the stability of neural integrator. Proceedings of the National Academy of Sciences, 101())):7739-771 77444, 2004.
P.E. McSharry, G.D. Clifford, L. Tarassenko, and L.A. Smith. A dynamical model for generating synthetic electrocardiogram signals. IEEE Transactions on Biomedical Engineering, 50(3):289- 294. March 2003.
Ted Meeds, Geoffrey Roeder, Paul Grant, Andrew Phillips, and Neil Dalchau. Efficient amortised Bayesian inference for hierarchical and nonlinear dynamical systems. In Proceedings of the 36th International Conference Machine Learning, 2019.
SJ Mizumori and JD Williams. Directionally selective mnemonic properties of neurons in the lateral dorsal nucleus of the thalamus of rats. The Journal of Neuroscience, 13(9):4015-402 1993.
Zahra Monfared and Daniel Durstewitz. Existence of n-cycles and border-collision bifurcations in -ineer continuous maps with applications to recurrent neural networks. Nonlinear Dynamics, 103 2020a.
Zarah Monfared and Daniel Durstewitz. Transformation of ReLU-based recurrent neural networks from discrete-time to continuous-time, In Proceedings of the 37th International Conference Machine Learning. 2020b.
Edward Ott. Chaos in Dynamical Systems. Cambridge University Press, 2002.
Adam Paszke, Sam Gross, Soumith Chintala, Gregory Chanan, Edward Yang, Zachary DeVito, Zeming Lin, Alban Desmaison, Luca Antiga, and Adam Lerer. Automatic differentiation in PyTorch. In NIPS Autodiff Workshop. 2017.
Maziar Raissi. Deep hidden physics models: Deep learning nonlinear partial differential equa- tions. Journal of Machine Learning Research, 19:1-24, 2018.
Hooshmand Shokri Razaghi and Liam Paninski. Filtering normalizing flows. 4th workshop on Bayesian Deep Learning (NeurIPS 2019), 2019.
Sam Roweis and Zoubin Ghahramani. Learning Nonlinear Dynamical Systems Using the Exxpetation-Maxiiiaatti Algorithm. In Kalman Filtering and Neural Networks, chapter 6, pp- 175-220. -iley-Blackwell, 2002.
Samuel H. Rudy, Steven L. Brunton, and J. Nathan Kutz. Smoothing and parameter estimation by soft-adherence to governing equations. Journal of Computational Physics, 398: 108860, Decem- ber 2019.
12

---

Published as a conference paper at ICLR 2021
Tim Sauer, James A. Yorke, and Martin Casdagli. Embedology. Journal of Statistical Physies, 65 (3-4):579-60 616, 1991.
H. Sebastian Seung. How the brain keeps the eyes still. Proceedings of the National Academy of Sciences, (3(23):13339-133
1996.
Hava T. Siegelmann and Eduardo Sontag. On the computational power of neural nets. Journal of Computer and System Sciences, 50(1):132-150, 1995.
H. Francis Song, Guangyu R. Yang, and Xiao Jing Wang, Training Excitatory-Inhiiitory Recurrent Neural Networks for Cognitive Tasks: A Simple and Flexible Framework. PLoS Computational Biology, 12(2):1-30, 2016.
Charles F Stevens. Neurotransmitter release central synapses. Neuron, 40(2):381-38) 2003.
Steven Strogatz. Nonlinear Dynamics and Chaos: Applications Physics, Biology, Chemistry, and Engineering: With Application to Physics, Biology, Chemistry and Engineering. CRC Press, 2015.
David Sussillo and Omri Barak. Opening the black box: low-dimensional dynamics in high- dimensional recurrent neural networks. Neural computation, 25(3):626-44 2013.
Floris Takens. Detecting strange attractors in turbulence. In Lecture Notes in Mathematics, pp- 366-381. Springer Berlin Heidelberg, 1981.
Sachin S. Talathi and Aniket Vartak. Improving performance of recurrent neural network with ReLU nonlinearity. ICLR Workshop submission, 2016.
Adam Trischler and Gabriele M.T. D'Eleuterio. Synthesis of recurrent neural networks for dy- namical system simulation. Neural Networks, 80:67-78, 2016.
Pantelis R. Vlachas, Wonmin Byeon, Zhong Y. Wan, Themistoklis P. Sapsis, and Petros Koumout- sakos. 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), 2018.
Eugene Vorontsov, Chiheb Trabelsi, Samuel Kadoury, and Chris Pal. On orthogonality and learning recurrent networks with long term dependencies. Proceedings of the 34th International Confer- ence on Machine Learning Reseurch, 70:3570-3578, 2017.
Simon N Wood. Statistical inference for noisy nonlinear ecological dynamic systems. Nature, 466: 1102-1104, 2010.
Byron M. Yu, Afsheen Afshar, Gopal Santhanam, Stephen I. Ryu, Krishna V. Shenoy, and Maneesh Sahani. Extracting dynamical structure embedded neural activity. In Advances in Neural Information Processing Systems 18, pp.
MIT Press, 2005.
Rose Yu, Stephan Zheng, Anima Anandkumar, and Yisong Yue. Long-Term Forecasting using Higher-Order Tensor RNNs. arXiv preprint, 2019.
Yuan Zhao and II Memming Park. Variational online learning of neural dynamics. arXiv preprint, 2020.
Xun Zheng, Manzil Zaheer, Amr Ahmed, Yuan Wang, Eric P. Xing, and Alexander J. Smola. State Space LSTM Models with Particle MCMC Inference. arXiv preprint, 2017.
13

---

Published as a conference paper at ICLR 2021
APPENDIX
6.1 SUPPLEMENTARY TEXT
6.1.1 Simple exact PLRNN solution for addition problem
The exact PLRNN parameter settings (cf. eq. 1, eq. 2) for solving the addition problem with units (ef. Fig. IC) are as follows:
A 0) w - :) ! h-
( ( -1 0)
(7)
6.1.2 Computation of fixed points and eycles in PLRNN
Consider the PLRNN in the form of eq. 4. For clarity, let define dn()) (d1, d2, ,dv) as an indicator vector with dm (2m.r)
dm 1 for all states 2m,t > 0 and zeros otherwise, and Doc) = diag(ddag (da()) as the diagonal matrix formed from this vector. Note that there are at most 2 distinct matrices WIX(e) as defined in eq. 4, depending on the sign of the components of z,.
If h = 0 and Wa(e) is the identity matrix, then the map F becomes the identity map and so every point will fixed point of F. Otherwise, the fixed points of can be found solving the equation F(z*1) z=1 as
= (I Wa(rer)) h2 H h,
(8)
2t-1 21-1-1, if dct(I Wa((r))) = PWarer, (1) 7 0, i.e. Wa(ee) has no eigenvalue equal to 1. Stability and type of fixed points (node, saddle, spiral) can then be determined from the eigenvalues of the Jacobian A WD = Water) (Strogatz (2015)).
where ==1
For 1, solving f* (2*h) = one can obtain a k-cycle of the map F with the periodic points fz* F(z*h F2 th) ).... , 1x For this, we first compute F as follows:
2: FZz-11) Wa(c-1) Zc-1 + h,
*+1 = (2 = F(z.) Wa(t) Wa(--)) zt-1 + (Wa(1) + 1)h.
2142 = (2)-1) = F( = Wo(+H) WoO Wn(t-1) #1-1
(Wag4))) Wao Wo(t-1) + I)h,
k+1
# k k-j+2
214(k-1) = Fh (z(-1) =
II
W () 2:-1 + 1=2
Σ
II
WaoHa-0
+
(9)
12
in which I1 k+1 W (o((+(k-1))) W Woo+(k22) 21) Wn((+(k-3)) Wn(-1)-
1-2
Assuming t + (k- 1) = tk then the k-cycle is given by the fixed point of the k-times iterated map Fk as
k
k-j+1
k
zk
W Warek 1)
Σ
II
W Wa + +1
h H h,
(10)
j=2
i=1
Zpsk = Zpek provided that I - - Wa is invertible. That is where tk
det I I- Wo( = PIL Ware (1) # 0 and I Wo == W has no
eigenvalue equal 1. As for the fixed points, we can determine stability of the k-cycle from the cigenvalues the Jacobians [ Wa
14

---

Published as a conference paper at ICLR 2021
It may also be helpful spell out the recursions in eq. and eq. 6 in section 3.3 in bit more detail. Analogously to the derivations above, for t = 1.2..... T we can recursively compute 22, 23:...; 2T
(T E N)
22 = F(z1) = Wa(1)= 21 + h:
23 (21) F(2)) Wa()) Wa(1) (Wa(2) + h)h,
27 = (21) = F(2r-1) = 1 W2(-22) .. Wa()
+ (WarT-) Wn(T-2)- Wn(2)
+ Wa(--) Wn(r-2) .. Wa(3) ... W1(--1) + 1)h
T-1
T-2 T T-j-1
W W + 1 h
= II Wo(T- ž1
1=1
j-1
II W Wo(T-s) 21
Σ
II W W I h.
(11)
i=1
Likewise, write out the derivatives eq. 6 more explicitly as
dz.
(F(z-1) = 1 (m, D zT-1 + (A W WDntt-1)) dwin
Uzi-1
dw
Owmn
= (m,,) Da(e-1) 2 (A+ W WDatt 1)) 1 (m,,)) 2)((t2) Zt-2
0z1-2
(A + WDD((-11) W (D(((22)) Owmk
(m, D 2 (A W WD 20
+ +AA + WD B( ) ( + WD (m,,)) (21(t-3) 2t-
+A W WD (2 1)) A WD (A (A+ W ()())))) Ozi-3
Owmn
1-2
j-1
Σ [II ) 1(m, Dast-1)) Zt-j
(m,,)) Da 2e-1 +
2
im1
1-2
Oz2
II
W ((t-s)) Dwin
(12)
where 300m
02.2 144 * m and 0zm.2 Swme dx The derivatives w.r.t. the elements of A and h can be expanded in similar way, only that the terms Ds(e) 2; on the
O32
021.2
@ZM.2 with Twon
Puin
last line of eq. 12 need to be replaced by just 2; for au.. in these cases, the indicator matrix will be the diagonal matrix (m),)))
,2,, and by just a vector I's for 8sc. (also, Sh.m
6.1.3 Proof of Theorem 1
To state the proof, let us rewrite the derivatives of the loss function L(W, A,h) ΣT La in the following tensor form:
ac
iC.
0L. 0L. 0z:
where
(13)
aw
Ow'
aw
Oz OW
15

---

Published as a conference paper at ICLR 2021
for which the 3D tensor
UzL aW
02g.t
Ōz. aw
JW
(14)
8EM.
JW
of dimension M. M x M, consists of all the gradient matrices
0=:
02s. BoriM Decr Zurg
2or11
Buri- asLa Twrg-
izit OW
Parg1
Zhega
i i= 1,2, .. ,M,
(15)
821,
Bas, Bu
Bust Bwaa
Ow
Baraz
where wi, E RM is a rewtor
Now, suppose that (21, #2, .............................. .. is an orbit of the system which converges stable fixed point,
L.. lim 2T xx Then
T-200
lim 27
lim (WaTT-1) zT-1 h) = x² Wattery 222 + h,
(16)
T
T000
and so
lim (Warr-1) x²1 Wa 211 T-
(17)
Assume that lim ( (Wa(r-)) = L. Since eq. 17 holds for every ="1, then substituting z=1 T-200 ef (1,0, 0) in eq. 17, can prove that the first column of L equals the first column of W Waee-.. Performing the same procedure for z=1 = e!_ 2,3, ..., M,, yields
lim Wae--1) = Wn(er). T-000
(18)
Also, for every E1 <i<
lim Wor-i- Wa(el), T-200
(19)
i.c.
Ve > 0 ENEN s... T T-i2N N Wa(T-) W Se.
(20)
Thus,
Wo(T-a)
| < Wa( Wa gives
Vc <>0 0 EN EN s.t. T 1 N W s Wox +c.
(21)
Since
- 1> 2 T > N,
VE 0
Worr s W +c,
N.
(22)
Hence
T-N
T-N
T-N W < |||| Wo [(]
Ve > 0
II
Werr
<
II
(23)
et W < 1, it is concluded that
W Wageey/ ||
< 1, then for any < 1, considering ē E
T-N
([Warr,))) +6) T-N 0. (24)
lim 7-400
II
W Worr-s
= lim
W Wa
lim T-x
T-x
=1
1=1
16

---

Published as a conference paper at ICLR 2021
Therefore
T-1
lim T-000
II
W ()(-() = 0.
(25)
the orbit {21, 22, #3....; tends to stable k-cycle (k: > 1) with the periodic points
(=*h), pk-1l(- (2 (z* F 2 (2=*h) ) F(=*))) (Z 2 " = ((k-))))))
then, denoting the stable k k-cycle by
F& = { 8z-t
(k- 1)r R -1 - ...
(26)
we have
lim d(zr, T-000
0.
(27)
Hence, there exists a neighborhood U of Ia and sub-sequences {zT.. Jm=1: {zT... 1
(2Tania (ik 1) 100 of the sequence { 1 such that these sub- sequences belong to and
(i) TTan+
F (2 -0,1,2
(ii) lim Tient ss 0,,,2,
1,
T-
(iii) for every 2T E there is some $ E {0. 1,2, .... 13 such that 2y {2Tent Imml-
In this case, for every 2T with zT E {*Ten* Jaz1 we have lim 2T Zt-k- for some
T-2
8 = 0, 1,2, ... h Therefore, continuity of F implies that lim F(zr) = F(z and so
T
lim (Waer) 2T + h) = Wag* -s) 2 + h. T-x
(28)
Thus, similarly, we can prove that
3s E {0, ... , 1---1p s.t.
lim Wa(r) = Wo T
(29)
Analogously, for every E (1 <i <i< < 00)
3s, E {0,1,2, ... , 1} s.t. lim Wa(r-i- = Wa(re-s;))
(30)
TH00
On the other hand, generality, assuming
W Wor+ ا|له < 1 for all 8; E 1,,1,2,... 1}. So, without loss of
max 05
| Wooe +1)l
Wa(ex)) < 1,
(31)
again obtain some relations similar to eq. 23-eq. 25 for ts k 1.
Since {zr-1 _ is convergent sequence, it is bounded, i.e. there exists real number for all T E N. Furthermore, || | x 1 for all T. Therefore, by eq. 12
such that |2 and eq. 23 (for k- 1)
0=T
Owink
T-1 J-1 Wacr-s) ((m,k) DocT- ZT-j
-(mm)) DD((-1)) +1+
Σ j=2
II
T-1
+
Wi(T-i) D1(1) 21
(32)
im1
T-1
S 1|
W n((T
|| +
II
W R(T-
||
$
q 1
7-1
Σ
(/ +2) +z)(-))
+ (| -c)=-1
T-1 ||
(33)
j=2
17

---

Published as a conference paper at ICLR 2021
Thus, by Wa +ē 1, we have
dzr
Wa +ē
lim
<
(34)
W Watetyl
Le., by eq. 14 and eq. 15, the 2-norm of total gradient matrices and hence (explode) under the assumptions Theorem
lla will not diverge
Analogously, we can prove that a ll and
Ozr 2h
ll will not diverge either. Since, similar in
the derivations above, can be shown that relation eq. 34 is true for Jamm
with q. where the upper bound of || as {#rIT=1 is convergent. Furthermore, relation eq. 34 also holds for
dzr
8zr 3hm
with q q= 1.
Remark 2.1. By eq. 24 the Jacobian parts
8zr 0z.
connecting any two states 2r and Zb T > t,
will not diverge either.
Corollary 2.1. The results of Theorem I are also true if Waseo a normal matrix with no eigen- value equal to one.
Proof. If W Wo is normal, then |W = p(Wa(e))) < 1 which satisfies the conditions of Theorem I.
6.1.4 Proof of Theorem 2
Let A, W and Dack) < k T, be partitioned follows
I Legg
O Anreg )
Oreg
O
D reg
OT
A
W =
D =
O
S
W W nreeg
D Düreg
(35)
O Mroy x Mrug = Oreg E Rreg x X O,S E where T Mrry = I Treg E Mragg XMang
R(M ((M-) M A AiMong 4: Mreg
A Anreg E (M) (M- is diagonal sub-
matrix, W I.M +1MM) = W Wnreg= E RR(M M (M- an off-diagonal sub-matrix (cf. Fig. S1). Moreover, DisrggMMrn X.reg == Di E Mreo X M and D D nreg E (M X(M X(M Meegl are diagonal sub-matrices. Then, we have
{ [M.y I:M, +1:M) ==
II Woag)) 1CkST
Ireg
о Diveg
I Ireg
OT
II tckST
D reg
Anreg
Wnreg D
S D reg
W preg
tckST
Ireg
OT
(36)
SD 11 reg
+ EJ
(Tlesksts-1 W nreg)
D reg
I
W nreg
Therefore, considering the 2-norm, we obtain
II WR( ICRST
0z:
I Ireg
<00.
S SD4+ (+1 + + (Theasstt-- Wk S D4+1
Whreg)
IIcckT
Wh
nreg
(37)
18

---

Published as a conference paper at ICLR 2021
Moreover
8zr
max (( p
II
Wooo)
II
Wa())
(38)
Ōzi
t<kET
tCkET
where Wr-t Neskk
W nrey' Therefore, eq. 37 and eq. 38 yield
Ozr
1 Plow
S
S Pup < co.
Oz
Furthermore, we assumed that the non-regularized subsystem (#M =M), if considered in isolation, satisfies Theorem 1. Hence, similar to the proof of Theorem is concluded that
limm
W nreg
O Onreg-
(39)
T-+00
On the other hand, by definition of Di(k), for every t < & s T, we have
D S and so
|| D &
||||| D
S |||
(40)
reg
which, in accordance with the the assumptions of Theorem nregl implies
by convergence of
E I W
b+j-1
++j-1
lim T T-000
S D
++1
Σ
II
W S D
S||s|| 1+ Jim T
II
[
(reall)
j=2
k=t+1
j=2 k=(+1
< ||S|| M Moreg-
(41)
Thus, denoting Q = DD
EJ= (Thessets-1 Whreg S D nreg
from eq. we deduce
rng
reg
that
Amaz lim TQT ()) = lim p(QTQ) Q) S lim T-000
= lim T
|| 2 S (IS) Mareg) ?
T
T-000
(42)
Now, T tends to cc, then eq. 37, eq. 39 and eq. 42 result
Ozr
Ireg
O nnegg (
O
))
Plow
Oman
=
Amar (Treg
lim (Q Q)) T
Q
Pup < oc.
(43)
Remark 2.2. |f||S|| 0, then
Bzz Oz
1 as T
xo.
6.1.5 Details on EM algorithm and DS reconstruction
For DS reconstruction we request that the latent RNN approximates the true generating system of equations, which is taller order than learning the mapping S + X or predicting future values in time series (cf. sect. 3.5). This point has important implications for the design of models, inference algorithms and performance metrics the primary goal is DS reconstruction rather than "mere' time series forecasting. 3 In this context consider the fully probabiliitic, generative RNN eq. Together with eq. 2 (where we take g(2,) = o(2,)) this gives the typical form of a nonlinear
By reconstructing the governing equations we mean their approximation in the sense of the universal approximation theorems for DS (Funahashi & Nakamura. 1993; Kimura Nakano, 1998), i.c. such that the hehavior of the reconstructed system hecomes dynamically equivalent that the true underlying system. "I this context we also remark that models which include longer histories of hidden activations (Yu et al., 2019), as in many statistical time series models (Fan & Yao, 2003). not formally valid DS models anymore since they violate the uniqueness of flow in state space (Strogatz, 2015).
19

---

Published as a conference paper at ICLR 2021
state space model (Durbin & Koopman, 2012) with observation and process noise. We solve for the parameters 0 = {A, W. C, h. po. Σ, B. by maximum likelihood, for which an efficient Expectation-Me DooctaationMMaximizztiioo (EM) algorithm has recently been suggested (Durstewitz, 2017; Koppe et al., 2019), which we will summarize here. Since the involved integrals are not tractable, we start off from the evidence-lower bound (ELBO) to the log-likelihood which can be rewritten in various useful ways:
log (( 2 -E log X, + (
llogp( DEL ( po (o((
=: (0,
(44)
In the E-step, given current estimate 6" for the parameters, we seek to determine the posterior Po () which we approximate by a global Gaussian g ( instantiated by the maximizer (mode) Z" of pa (2 as an estimator of the mean, and the negative inverse Hessian around this maximizer as an estimator of the state covariance, i.e.
E X ~ Z = arg max log X)
arg max log pe (X + log pe (Z) log pe (X)]
z
arg max llog pe X Z + log pe pe Z)]
(45)
z
since Z integrates out in po (x) (equivalently, this result can be derived from a Laplace approxima- tion to the log-likelihood, log p( X (0) ~ log x Po (Z*) 2 log L' const, where L" is the Hessian evaluated at the maximizer). We solve this optimization problem by a fixed-point iteration scheme that efficiently exploits the model's piecewise linear structure, detailed below.
Using this approximate posterior for po X based on the model's piecewise-linear structure most of the expectation values Ezzg [Ø(=)], Exnug [0(z)=], and Ex~g [o(z))(z)T] could be solved for (semi-)ane Jemll-yanalytically (where is the concatenated vector form of Z. see below). In the M-step, we seek 0 == arg maxg L(O, ))),, assuming proposal density 9 to be given from the E-step, which for a Gaussian observation model amounts to a simple linear regression problem (see Suppl. eq. 49). To force the PLRNN to really capture the underlying DS in its governing equations, we use previously suggested (Koppe et al., 2019) stepwise annealing protocol that gradually shifts the burden of fitting the observations X from the observation model eq. 2 to the latent RNN model eq. during training, the idea of which is to establish a mapping from latent states Z to observations X first, fixing this, and then enforcing the temporal consistency constraints implied by eq. while accounting for the actual observations.
Now we briefly outline the fixedpooint-iteration algorithm for solving the maximization problem in eq. 45 (for more details Durstewitz (2017); Koppe et al. (2019)). Given Gaussian latent PLRNN and a Gaussian observation model, the joint density p(X, Z) will be piecewise Gaussian, hence eq. 45 piecewise quadratic in Z. Let us concatenate state variables across m and into
one long column vector = (21,1, 2M.1.-... 21,7: 2M,T) arrange matrices A, W into
T
large MT MT block tri-diagonal matrices, define d, == (I=,,>0: =2,10 ... 1 1 >0) an indicator vector with 1 for all states 2m,t 0 and zeros otherwise, and Dn diag(dn) as the diagonal matrix formed from this vector. Collecting all terms quadratic, linear, or constant in z, we can then write down the optimization criterion in the form
1 12T (Uo DRU1 U.Do + Dn)
@h(z) = - 2
IT (vo + Davi) - (v + DAvi) T x]
const.
(46)
In essence, the algorithm now iterates between the two steps:
Given fixed D3. solve
= =U DDU, + +U Da + DDU,Da) -1. -(0 +Dvu))
(47)
2. Given fixed z-, recompute Dn
20

---

Published as a conference paper at ICLR 2021
until either convergence one of several stopping criteria (partly likelihood- lkeeiihodd-based, partly to avoid loops) is reached. The solution may afterwards be refined by one quadratic programming step. Numerical experiments showed this algorithm to be very fast and efficient (Durstewitz, 2017; Koppe et al., 2019). At z", an estimate of the state covariance is then obtained as the inverse negative Hessian,
V= (Uo + DRU1 + UJ Dn + DR)
(48)
In the M-step, using the proposal density q" from the E-step, the solution to the maximization problem 0* = arg max L(0, q²), can generally be expressed in the form
0
(>
[s,33]] [[
0" = Σ [anB!]]
(49)
where, for the latent model, eq. 1, ax = 2e and Br = | ol (z1-1) 1) sT.] RRM and
for the observation model, eq. 2, a; = ie and Br = g (2t).
6.1.6 More details on DS performance measure
As argued before (Koppe et al., 2019; Wood, 2010), in DS reconstruction we require that the RNN captures the underlying attractor geometries and state space properties. This does not necessarily entail that the reconstructed system could predict future time series observations more than a few time steps ahead, and vice versa. For instance, if the underlying attractor is chaotic, even if we had the exact true system available, with a tiny bit of noise trajectories starting from the same initial condition will quickly diverge and ahead-prediction errors become essentially meaningless as DS performance metric (Fig. S2B).
To quantify how well an inferred PLRNN captured the underlying dynamics we therefore followed Koppe et al. (2019) and used the Kullback-Libler divergence between the true and reproduced prob- ability distributions across states in state space, thus assessing the agreement in attractor geometries (cf. Takens (1981); Sauer et al. (1991)) rather than in precise matching of time series,
K
(R)
Pirue (ar)
DKL (Porue (x))lge (x =)) R
Perue (x) log
(50)
Poen (x|x)
where Ptrue(2)) is the true distribution of observations across state space (not time!), Pgen (z|²) is the distribution of observations generated by running the inferred PLRNN, and the sum indicates spatial discretization (binning) of the observed state space. We emphasize that pltd (x|²) is obtained from freely simulated trajectories, i.e. drawn from the prior p(z) specified by eq. 1, not from the inferred posteriors p(2 [rrain). In addition, to assess reproduction of time scales by the inferred PLRNN, the average MSE between the power spectra of the true and generated time series was computed, as displayed in Fig. 3B-C
The measure DKL. introduced above only works for situations where the ground truth Perue (X) is known. Following Koppe et al. (2019), we next briefly indicate how a proxy for DKL may be obtained empirical situations where no ground truth is available. Reasoning that for well recon- structed DS the inferred posterior Pine(2 given the observations should be good representative of the prior generative dynamics Pgen(2) (2), one may use the Kullback-Leibler divergence between the distribution over latent states, obtained by sampling from the prior density P (2), and the (data-
constrained) posterior distribution Pint(z|a (z|2z (where E IRM x1 system's state space:
and x E R I), taken across the
-f IZERMxX
Pint(2/22) log Pinr(2/1z). dz Pgea (2)
DKL (pin (2]
(51)
As evaluating this integral is difficult, one could further approximate Pint (2(x)) and Pgen by Gaussian mixtures across trajectories, i.c. Pan((2/z) ss $ -T= p(zz (1:T) and Pgen(z) se 1. EE= p(2) 1)),,,, where the mean and covariance of p(z, (1rT) and p(zi =-)) are obtained by marginalizing over the multivariate distributions p(Z and Pgen (Z), respectively, yielding [| [( and covariance matrices Var(z, m.r) and Var(zi|2)-) Supplementary
21

---

Published as a conference paper at ICLR 2021
eq. 51 may then be numerically approximated through Monte Carlo sampling (Hershey & Olsen, 2007) by
Pint (znt(z1)) (z))))) |x))
DKL. (pint (z)) ~ log
Pgen (z()))
~ Pinf (z|
(52)
Alternatively, there is also variational approximation of eq. 51 available (Hershey & Olsen, 2007):
T
Σ j=1 e eeDruulwtaaeeerrrllleeeeeeee
I
1
DKL
Dyariational (2tnt(//2)) (z)) ~
log
(53)
T
Ek=1 e - t))
where the KI. divergences in the exponentials are among Gaussians for which we have an analytical expression.
6.1.7 More details on benchmark tasks and model comparisons
We compared the performance of our rPLRNN to the other models summarized in Suppl. Table on the following three benchmarks requiring long thort-term maintenance of information (Talathi & Vartak (2016); Hochreiter & Schmidhuber (1997)): 1) The addition problem of time length T consists of 100000 training and 000 test samples of 2 x input series S ( ST): where entries 81, E [0, 1] are drawn from uniform random distribution and $2, E {0, 1) contains zeros except for two indicator bits placed randomly at times f < 10 and t2 < T/2. Constraints on
t and f2 are chosen such that every trial requires a long memory of at least T/2 time steps. At the last time step T, the target output of the network is the sum of the two inputs in 81,: indicated by the l-entries in 82,;, FT target
81,41 + 81.12- 2) The multiplication problem is the same as the addition problem, only that the product instead of the sum has to be produced by the RNN as an output at time T, IT 60 000 training and 000 28 x 28 test images of hand written digits. To make this a time series problem, in sequential MNIST the images are presented sequentially, pixel-by-pixel, scanning lines from upper left to bottom- right,, resulting in time series of fixed length T = 784.
Earget
81,41 81,f2 3) The MNIST dataset (LeCun al., 2010) consists of
For training on the addition and multiplication problems, the mean squared-error loss across R 2
samples, L = i ER= (# - -xy [(,)))
between estimated and actual outputs was used, while
(n) the enoss-y loss L = ER= (- I1 2 log(p!7)) was employed for sequential MNIST, where
Pi, == Do (xi.e = 1|zx)= Bi.
Bj Bj.
(54)
with Fut 10. Σg Z4, 1. We remark that long the observation model takes the form of a generalized linear model (Fahrmeir & Tutz, 2001), as assumed here, meaning may be assigned to the latent states 2m by virtue of their association with specific sets of observations z, through the factor loading matrix B. This adds another layer of model interpretability (besides its accessibility in DS terms).
The large error bars in Fig. at the transition from good to bad performance result from the fact that the networks mostly learn these tasks in an all-or-none fashion. While the rPLRNN in gen- eral outperformed the pure initialization-base models (iRNN, npRNN, iPLRNN), confirming that a manifold attractor subspace present initialization may be lost throughout training. we conjec- ture that this difference in performance will become even more pronounced as noise levels or task complexity increase.
22

---

Published as a conference paper at ICLR 2021
Table 1: Overview over the different models used for comparison
NAME
DESCRIPTION
RNN
Vanilla ReLU based RNN
iRNN
RNN with initialization Wo = and ho (Le et al., 2015)
npRNN
RNN with weights initialized to a normalized positive definite matrix with largest eigenvalue of and biases initialized zero (Talathi & Vartak, 2016)
PLRNN
PLRNN as given in eq. (Koppe et al., 2019)
iPLRNN
PLRNN with initialization Au I, Wo 0 and ho 0
rPLRNN
PLRNN initialized as illustrated in Fig. S1, with additional regular-
ization term (eq. 3)
LSTM
Long Short-TTerm Memory (Hochreiter & Schmidhuber, 1997)
ORNN
ReLU RNN with W regularized toward orthogonality (W WT I)
(Vorontsov et al., 2017)
L2RNN
Vanilla RNN with standard L2 regularization all weights
L2pPLRNN
PLRNN with (partial) standard L2 regularization for proportion Mre/ M 0.5 of units (i.e., pushing all terms in A and W for these units to 0)
L2APLRNN
PLRNN with (full) standard 1.2 regularization on all weights
6.1.8 More details on single neuron model
The neuron model used in section 4.2 is described by
-Cm V (V FL)) + mo (V)(V EENa)
vVnn(V - Fxx) + gMh(V - Fxx)
+ (V) ENMDA)
(55)
h 1100 (V) -
(56)
Th
à noc (V)
(57)
(V) [1 + .33e 06251 -00250
(58)
where C'm refers to the neuron's 's membrane capacitance, the ga to different membrane conductances, F. to the respective reversal potentials, and m, h, and n are gating variables with limiting values given by
{Me. Roo, h2o = [1 + el ((Vi. VV) (x
(59)
Different parameter settings in this model lead to different dynamical phenomena, including regular
spiking, slow bursting or chaos (see Durstewitz. (2009) for details). Parameter settings used here were: Cm = 6FF, 9r 8 mS. F,, = -80 mV INo 20 mS, F. = 60 VnNa = -20 V.
kNo 15, 9K 10, ÉK -90 mV, VyK -25 m. kK 5, Tn 1 ms. gM = 25 mS,
Vhas = -55m V, KM = 5, Tr = 200 ms, 9NMDA = 10.2 mS.
23

---

Published as a conference paper at ICLR 2021
6.2
SUPPLEMENTARY FIGURES
0
0
0
0
0
X X
1OOXXX
00000
0
0
0
0000
0
0
0
X
0
AX
X
0
X
X
0
X
A
W
h
Illustration of the attp:iddddattaattteeeaaaaraaaeae matrix A, coupling matrix W, and bias terms h. Regularized values are indicated in red, crosses mark arbitrary values (all other values set to 0 as indicated).
Figure S1:
for the PLRNN's auto-regression
A
-20
-20
400
-60
0.1
0.3
Lime (6)
Lime (6)
0.06
0.06
0.04
0.04
0.5
D
-50
V
Figure S2: MSE evaluated between time series is not good measure for DS reconstruction. A) Time graph (top) and state space (bottom) for the single neuron model (see section 4.2 and Suppl. 6.1.8) with parameters in the chaotic regime (blue curves) and with simple fixed point dynamics in the limit (red line). Although the system has vastly different limiting behaviors (attractor geometries) in these two cases, as visualized in the state space, the agreement in time series initially seems to indicate a perfect fit. B) Same as A) for two trajectories drawn from exactly the same DS (i.c., same parameters) with slightly different initial conditions. Despite identical dynamics, the trajectories immediately diverge, resulting in a high MSE. Dash-dotted grey lines top graphs indicate the point from which onward the state space trajectories were depicted.
24

---

Published as a conference paper at ICLR 2021
A
10
10-1
10
10
MSE
0.1
MSE
MSE
10-2
10
10"4
10
10.0
10
30
50
0.3
1
20
50
100
0.1
0.5
0.7
0.9
Meag/ /M
10
1.5
cross entropy
MSE
0.6
10-4
0.3
0.1
0.3
0.5
0.7
0.1
0.3
0.5
0.7 0.9
Men /M
Ma
Figure S3: Performance of the rPLRNN for different A) numbers of latent states M, B) values of т, and C-E) proportions Mreg/Au /M of regularized states. A-C are for the addition problem, for the multiplication problem, and E for sequential MNIST. Dashed lines denote the values used for the results reported section 4.1.
A
B
PE
OPLRNN
10
:PLRNN
L2PLRNN
1XPPLRXN
363
1n4
16h
<- 1500
T. 1500
Figure S4: A) 20-step-ahead prediction error between true and generated observations for rPLRNN function of regularization B) KL divergence (DRL. between true and generated state space distributions for orthogonal PLRNN (oPLRNN; i.e., the PLRNN with the 'manifold attractor regu- larization' replaced by an orthogonality regularization, (A (A + W)" - I), well for the partially (L2p) and fully (L2f) standard L2-regularized PLRNNs (i.e., with all weight parameters for all (L.2f) or only fraction Mrog/M /M of states (L2p) driven to 0). Note that the quality of the DS reconstruction does not significantly depend on the strength of regularization т, or becomes even slightly worse, for the oPLNNN, L2pPLRNN and L21PLRNN. Globally diverging estimates were removed.
25

---

Published as a conference paper at ICLR 2021
A
gensrraedd
-3
-2
0.6
time (s)
time (5)
cime (4)
B
lascutt states
-3-
simo (6)
timse (6)
sime
2-
Insent staze
-10
0.3
0.45
0.75
0.15
0.65
0.6
0.78
0.3
0.45
Simo (8)
tiroe (6)
simo
Figure S5: A) Reconstruction of fast gating variable n (rightmost) not shown in Fig. 3D. For com- pleteness and comparison, other variables have been re-plotted from Fig. 3D as well. B) Example of reconstruction voltage (V, left) and slow gating (h, center) observations, and underlying latent state dynamics (right) for oPLRNN (with orthogonality regularization on A + W, see Fig. S4 leg- end). C) Example V (left) and h (center) observations for standard PLRNN, and underlying latent state dynamics (right). In general, both the standard and the oPLRNN tended to produce many fixed point solutions.
those cases where this was not the case, the standard PLRNN tended to reproduce only the fast components of the dynamics as in the example in C (in agreement with the results in Figs. 3C & 3E), while the oPLRNN tended to capture only the slow components as in the example in B (as expected from the fact that the orthogonality constraint tends produce solutions similar to those obtained for the regularized states only, of. Fig. 3E).
26

---

Published as a conference paper at ICLR 2021
A
B
0.08
0.5
bortawlized MSE
0.03
0.3
0.02
0.2
0.01
1e5
1ol
2000
2600
3600
D3
Kencroted
17
-1
-2
(5)
Lizur
time
Figure $6: Reconstruction of DS with multiple time scales like fast spikes and slow T-waves (simu- lated ECG signal, see McSharry al. (2003)). A) KL divergence (DKL) between true and generated state space distributions as function of r. Unstable (globally diverging) system estimates were re- moved. B) Average MSE between power spectra (slightly smoothed) of and reconstructed DS.
C) Average normalized MSE between power spectra of true and reconstructed DS split according to low (S 2.5 Hz) and high (> 2.5 Hz) frequency components. Error bars SEM in all graphs. D) Example of (best) generated time series (standardized, red=reconstruction with 1000/3600),
A
Addition problem
B
Multiplication problem
Sequential MNIST
2.5
10a2
10-1
CIOSS entrrpy loss
10
MSE
10
10
1.5
10
10-1
10-5
10-4
10
0.5
10
10-7
20
30
70
150 250 400
20
50
150
250
400
900
RNN
LSTM
PLINN
Lam LRW
1op RNNN
12RUN
RNN 12PPLRNN
LSTM (Hochhetter & Schusichaber, 1997)
iPLRNN
12pPLRNN
L2RNN
0.8
0.8
Comec
Porret
0.6
0.2
0.2
0
20
30
50
70
150 250
400
20
30
550
150 250
400
900
RNW
LSTM
PPLRNN
LaPLONS
1SPPLRNN
LORNN
T
Figure S7: Same Fig. 2, illustrating performance for L2RNN (vanilla RNN with L.2 regulariza- tion on all weights) and L.2PPLRNN (PLRNN with 1.2 regularization on all weights) on the three problems shown in Fig. 2. Note that the L2IPLRNN is essentially not able to learn any of the tasks, likely because conventional L2 norm drives the PLRNN parameters away from a manifold attractor configuration (as supported by Fig. 4 and Fig. S8). Results for rPLRNN, vanilla RNN, L2pPLRNN, and LSTM have been re-plotted from Fig. for comparison.
27

---

Published as a conference paper at ICLR 2021
A
B
PLRNN
IPLRNN
relative frequency
10
12pPLRNN
0.6
min A
0.2
0.5
1,20
L.p
PL
rPL
X maa around FPs
D
PLRNN
relative frequency
#LLRNN
0.8
0.8
L2PLRNN
12PPLRNN
0.6
0.6
0.4
0.4
0.2
0.2
0
-0.2
0
0.5
1.5
L2f
L2p
PL
rPL
Albaa around FPs
Figure S8: Same as Fig. 4 for A, B) multiplication problem, and C, D) sequential MNIST, Error bars stdv.
0.1
2
0.9
0.05
1
0.7
0.05
1
0ox
-0.1
0
1el
1e2
163
let
le5
Tel
162
168
led
165
lcl
lc2
1e3
le4
le5
7-
1500
T-
1500
T 1500
Figure S9: Effect of regularization strength 7 on rPLRNN network parameters (cf. eq. 1) (regular- ized parameters for states < Mreg. eq. in red). Note that some of the non-regularized network parameters (in blue) appear to systematically change as well as 7 is varied.
28

---

Published as a conference paper at ICLR 2021
10a
RNN
L2RNN
RNNN
ppRNN
PLRNN
10
iPLRNN
rPLRNN
LSTM
LOss
100
10
50
100
150
200
epoch
Figure S10: Cross-entropy loss as function of training epochs for the best model fits the sequential MNIST task. Note that LSTM takes longer to converge than the other models. LSTM training was therefore allowed to proceed for 200 epochs, after which convergence was usually reached, while training for all other models was stopped after 100 epochs. Also note that although for the best test performance on seq. MNIST shown here LSTM slightly supersedes rPLRNN, on average rPLRNN performed better than LSTT (as shown in Fig. 2C), despite having much fewer trainable parameters (when LSTM was given about the same number of parameters as IPLRNN, i.c. M M/4, its performance fell behind even more).
true
reconstructed
N
-2
-2
-2
2
-2
0
2
Figure S11: Example reconstruction of a chaotic system, the famous 3d Lorenz equations, by the
rPLRNN, Left: True state space trajectory of Lorenz system; right: trajectory simulated by rPLRNN 100/T, M 14) after training on time series of length T 1000 from the Lorenz system.
29


Reviews:
[{"review_id": "zL_CGqs57HE", "reviewer": "AnonReviewer1", "summary": "", "questions": "", "limitations": "", "rating": 7, "confidence": 5, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "This paper proposes a type of regularization for piecewise linear RNNs (PLRNNs) that encourages the network to learn line or plane attractors. The paper argues, through mathematical analysis of the regularized network as well as numerical experiments, that this regularization alleviates the vanishing and exploding gradient problem and allows PLRNNs to reconstruct nonlinear dynamical systems with multiple timescales from noisy observations.\n\nMajor concern:\nI found the paper interesting. My main concern has to do with the motivation of learning line or plane attractors. The paper argues that the proposed regularizer will improve performance *through a specific mechanism: encouraging plane attractors in the dynamics*. While the results in the experiments section are impressive, the paper as far as I can tell does not establish that the regularized PLRNNs have learned plane attractors.\n\nFor example, it could be that simply adding l2 regularization on all of the weights would also lead to better performance (rather than the specific regularization proposed), or perhaps the better performance arises from some other (undiscovered) mechanism. In fact, to show that the specific form of l2 regularization is what is useful, I think having a PLRNN with standard l2 regularization (applied to all of the weights) is a critical missing baseline.\n\nTo convince me that the benefit is really due to the given motivation (encouraging plane attractors), I want to see evidence that regularized PLRNNs have learned plane (or line) attractors, compared to unregularized PLRNNs (or better yet, PLRNNs with l2 regularization applied to all of the weights). This can be demonstrated in a number of ways, for example by showing eigenvalues of the Jacobian around fixed points of trained PLRNNs with and without the proposed regularization, on the three tasks in Fig. 2. If this was conclusively demonstrated, I would happily increase my rating.\n\nFor comparison, recent work showed that RNNs of multiple types (including vanilla RNNs and LSTMs) learned line attractors when solving an NLP task [1]. These line attractors were found in the networks after training, and did not have any special regularization encouraging them. I would appreciate if the authors would add some discussion comparing their regularization to these discovered line attractors.\n\n[1] Maheswaranathan et al, NeurIPS 2019 (http://papers.nips.cc/paper/9700-reverse-engineering-recurrent-networks-for-sentiment-classification-reveals-line-attractor-dynamics)\n\nOther concerns:\n- How is Mreg (the number of regularized dimensions) chosen? How does the performance vary as a function of Mreg?\n- What's g in eq(2)?\n- I was a little surprised that the vanilla RNN worked so well in Fig 2A. I would have expected an LSTM to work as well, if not better. Do the authors have intuition for why this is?\n\nMinor typos:\n- Quotes around automatize in first and last paragraphs of section 1 are incorrect. should be `automatize' and `classical'. in latex, use the backtick character (`) for the first quote.\n- The abbreviation \"RNN\" is used as if it is plural, but it is more commonly singular. (e.g. RNNs seem like instead of RNN seem like; or RNNs in their vanilla form instead of RNN in their vanilla form)", "extra_scores": {}, "extra_text": {}}, {"review_id": "Gddl40eQTz", "reviewer": "AnonReviewer4", "summary": "", "questions": "", "limitations": "", "rating": 8, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "The paper explores a very important question in dynamical system identification of how to make recurrent neural networks (RNNs) learn both long-term and short-term dependencies without the gradient vanishing or exploding limitation. They suggest using piece-wise linear RNNs (PLRNNs) with a novel regularization technique.\n\nThe paper is well written and is very thorough with the necessary theoretical foundation, numerical experiments and analysis. \n\nI think the theory and results of this paper are significant and will be relevant to further our understanding of RNNs and system identification.\n\nMajor points:\n\n1) L2 weight regularization can be easily applied to any of the RNN models used in the experiments. While other weight initialization schemes were compared to the paper's proposed model (rPLRNN), none of the other RNN models had similar regularization. This will shed some light on whether it is indeed the proposed regularization that matters or the full proposed model with PLRNN and a mix of regularized and non-regularized units.\n\n2) It is not clear to me how one can choose the correct ratio of regularized vs unregularized units in the model. While the amount of regularization clearly helps in reducing training error as shown in Figure 3, increasing the ratio of regularized units in Figure S3C did not help the error past 0.1 and then larger values resulted in large increases such that the error at ratio 1 is equivalent to the error at ratio 0. Perhaps this observation is specific to the addition problem, but I feel that a discussion of the effect of this ratio on performance should be included for clarity. Additionally, the ratio of regularized units with best performance could potentially be different for different regularization amounts.\n\nMinor Point:\ng is not defined in equation 2", "extra_scores": {}, "extra_text": {}}, {"review_id": "F6thXffhEj", "reviewer": "AnonReviewer3", "summary": "", "questions": "", "limitations": "", "rating": 7, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "This paper proposes a regularization scheme for training vanilla Relu RNN to tackle the exploding and vanishing gradients issue. The work eases the analysis of RNN in the dynamical system point of view and connects the RNN dynamics and gradient theoretically. The experiments show the competitive performance comparing to LSTM.\n\nThe vanilla RNNs simplify the analysis as dynamical systems without those gates in LSTMs and GRUs while suffering exploding or vanishing gradients. The idea of tackling the issue  while keeping the simplicity sounds very interesting and useful.\n\nI am leaning on accepting this manuscript if the authors could address my concerns.\n\n- Sec 3.2. The authors mentioned the connection between the particular setting of RNN and working memory. This setting leads to a system without any autonomous dynamics that is not a typical model for working memory (e.g. attracting fixed-points, line attrators and etc.). I disagree that the space is neurally stable since all the states are sensitive to perturbation and cannot persist a stable memory.\n\n- Sec 3.2. Following the so-called \"neurally stable\" setting, why the term on A in Eq.3 only regularizes the diagonal rather than penalizes the deviation from the identity matrix? This regularization does not lead to the tendency of A -> I mentioned in the text.\nFurthurmore, the regularization does not guarantee either A -> I or W -> 0 so that the resulting Mreg subspace does not have the \"memory\" property.\n\n- Eq.3. How is Mreg determined?\n\nMinor:\n\n- iPLRNN is used before defined.\n\n", "extra_scores": {}, "extra_text": {}}, {"review_id": "BVQgVSSJN5d", "reviewer": "AnonReviewer2", "summary": "", "questions": "", "limitations": "", "rating": 6, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "The paper proposes a novel regularization term to PLRNN. PLRNN has nice numerical properties given its simple mathematical structure, but is able to capture complicate dynamics. It's also easy to establish a theoretical connection between PLRNN dynamics and the behavior of its gradients, which is nice. Given such a dynamic model, they design novel L2 terms pushing partial parameters towards zero, which leads to a line/plane attractor that allows slow time constants for long short-term memory. I think the idea is quite novel and interesting.\n\nA few concerns:\n1. It might be a bit jumpy for readers outside neuroscience to establish the connection between working memory and short-term memory. I would say making this point more explicit in the intro would be helpful.\n2. Lack of definition of G at the beginning of 3.1.\n3. It would be helpful to point out that eq 1 is PLRNN (with full name spell spelled out).\n4. Section 4.1 is a bit confusing to me. It's unclear that why the higher the mse/cross entropy is, the better is the model should be; also the lower the Pcorrect value is, the better the model is. Moreover, I think it's still necessary to give the definition of all the model names included in the main text. Some are missing, e.g. iPLRNN, oRNN, maybe just very briefly. Just expect that not all readers would read the appendix.\n5. The paper seems to care more about \"interpretable\", which is not clearly reflected in the paper. Figure 3 only shows the reconstruction. But it would be more interesting to visualize the latents of this neuron model's dynamic. It's mentioned that M={8,...,18} states were trained. How do they look like? Some are more line-attractors and some are more related to the fast spiking dynamics? Which state number is finally picked? Figure 3 only shows how influential \\tao is. but there should be other insightful ablation studies to be done in order to understand the regularization term as well, i.e. how many latent states, how to split the two types, etc.\n6. It would also be helpful to visualize the reconstruction and the latents for other RNN/LSTM models for the neuron model as well. That helps to show the advantage of the proposed model more clearly. \n\nIn sum, I think the idea is quite interesting and practically useful. I also appreciate the theoretical analysis. But the experiment section is confusing by missing some explanations. And the presentation of the neuron model is not sufficient enough to prove that rPLRNN find interesting and interpretable dynamics. That's why I think it's a bit below a good paper. \n\n\n\n\n", "extra_scores": {}, "extra_text": {}}]

Author rebuttal:

Meta-review: This paper describes a clever new class of piecewise-linear RNNs that contains a long-time scale memory subsystem. The reviewers found the paper interesting and valuable, and I agree. The four submitted reviews were unanimous in their vote to accept. The theoretical insights and empirical results are impactful and would be suitable for spotlight presentation.
Decision: Accept (Spotlight)