Inferring Dynamical Systems with Long-Range Dependencies through Line Attractor Regularization
Authors: Dominik Schmidt
Georgia Koppe
Max Beutelspacher
Daniel Durstewitz
Venue:  (2020)
Keywords: Recurrent Neural Networks
Nonlinear State Space Models
Generative Models
Long short-term memory
vanishing/exploding gradient problem
Nonlinear dynamics
Interpretable machine learning
Time series analysis

Abstract: Vanilla RNN with ReLU activation have a simple structure that is amenable to systematic dynamical systems analysis and interpretation, but they suffer from the exploding vs. vanishing gradients problem. Recent attempts to retain this simplicity while alleviating the gradient problem are based on proper initialization schemes or orthogonality/unitary constraints on the RNN’s recurrency matrix, which, however, comes with limitations to its expressive power with regards to dynamical systems phenomena like chaos or multi-stability. Here, we instead suggest a regularization scheme that pushes part of the RNN’s latent subspace toward a line attractor configuration that enables long short-term memory and arbitrarily slow time scales. We show that our approach excels on a number of benchmarks like the sequential MNIST or multiplication problems, and enables reconstruction of dynamical systems which harbor widely different time scales.

Under review as a conference paper at ICLR 2020
INFERRING DYNAMICAL SYSTEMS WITH LONG-RANGE DEPENDENCIES THROUGH LINE ATTRACTOR REGU- LARIZATION
Anonymous authors Paper under double-blind review
ABSTRACT
Vanilla RNN with ReLU activation have simple structure that is amenable to sys- tematic dynamical systems analysis and interpretation, but they suffer from the ex- ploding vanishing gradients problem. Recent attempts retain this simplicity while alleviating the gradient problem are based on proper initialization schemes orthooonaliyy/mitary constraints the RNN's recurrence matrix, which, how- ever, comes with limitations to its expressive power with regards to dynamical systems phenomena like chaos or multi-stabiliit. Here, we instead suggest reg- ularization scheme that pushes part of the RNN's latent subspace toward a line attractor configuration that enables long short-term memory and arbitrarily slow time scales, We show that our approach excels on a number of benchmarks like the sequential MNIST or multiplication problems, and enables reconstruction of dynamical systems which harbor widely different time scales.
INTRODUCTION
Theories of complex systems in biology and physics often formulated in terms of sets of stochas- tie differential or difference equations, i.e. as stochastic dynamical systems (DS). A long-standing desire is retrieve these generating dynamical equations directly from observed time series data (Kantz& & Schreiber, 2004). A variety of machine and deep learning methodologies toward this goal have been introduced in recent years (Chen al., 2017: Champion et al., 2019: Jordan et al., 2019; Duncker et al., 2019; Ayed et al., 2019; Durstewitz, 2017: Koppe et al., 2019), many of them based on recurrent neural networks (RNN) which can universally approximate any DS (i.e., its flow field) under some mild conditions (Funahashi & Nakamura, 1993; Kimura & Nakano, 1998). However, vanilla RNN as often used this context are well known for their problems capturing long-term dependencies and slow time scales in the data (Hochreiter & Schmidhuber, 1997; Bengio et al., 1994). In DS terms, this is generally due to the fact that flexible information maintenance over long periods requires precise fine-tuning of model parameters toward 'line attractor' configurations (Fig. 1), a concept first propagated in computational neuroscience for addressing animal perfor- mance in parametric working memory tasks (Seung, 1996; Seung al., 2000; Durstewitz, 2003). Line attractors introduce directions of zero-flow into the model's 's state space that enable long-term maintenance of arbitrary values (Fig. 1). Specially designed RNN architectures equipped with gat- ing mechanisms and (linear) memory cells have been suggested for solving this issue (Hochreiter & Schmidhuber, 1997; Cho al., 2014). However, from DS perspective, simpler models that can more easily be analyzed and interpreted in DS terms, and for which more efficient inference algorithms exist that emphasize approximation of the true underlying DS would be preferable.
Recent solutions to the vanishing vs. exploding gradient problem attempt to retain the simplicity of vanilla RNN by initializing or constraining the recurrent weight matrix to be the identity (Le al., 2015), orthogonal (Henaff et al., 2016; Helfrich et al., 2018) 1) unitary (Arjovsky et al., 2016). In this way, in a system including piecewise linear (PL) components like rectified-linear units (ReLU), line attractor dimensions are established from the start by construction or ensured throughout training by a specifically parameterized matrix decomposition. However, for many DS problems, line attractors instantiated by mere initialization procedures may be unstable and quickly dissolve during training. On the other hand, orthogonal or unitary constraints are too restrictive for reconstructing DS, and more generally from computational perspective well (Kerg et al., 2019): For instance, neither

---

Under review as a conference paper at ICLR 2020
A
B
C
=
Au=1
inputs
Wiz=1
Car1 1
ide e
===
C2.11
hy=1
21
Fig. 1: Line attractors for solving long-time-scale problems. A)-B): Illustration of the state space of a 2-unit RNN (converted into a continuous time ODE, see Suppl. 7.1.2) with flow field (grey) and nullclines (set of points 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 from all directions and is exactly zero the line, thus retaining states indefinitely in the absence of perturbations, as illustrated for 3 example trajectories (green) started from different initial conditions. B) Slightly detuned line attractor (cf. Durstewitz (2003)). The system's state still converges toward the 'line attractor ghost' (Strogatz, 2015), but then very slowly crawls up within the "attractor tunnel' (green trajectory) until 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. Note that along other, not illustrated dimensions of the system's state space the flow may still evolve freely in all directions. C) Simple 2-unit solution to the addition problem exploiting the line attractor properties of ReLUs in the positive quadrant. The output unit serves perfect integrator, while the input unit will only convey those input values to the output unit that are accompanied by a II in the second input stream (see 7.1.1 for complete parameters).
chaotic behavior (that requires diverging directions) nor settings with multiple isolated fixed point or limit cycle attractors are possible.
Here therefore suggest different solution to the problem, by pushing (but not strictly enforeing) ReLU-based, piecewise-lined RNN (PLRNN) toward line attractor configurations along some (but not all) directions in state space. We achieve this by adding special regularization terms for subset of RNN units to the loss function that promote such a configuration. We demonstrate that our approach outperforms,
en par with, LSTM and other, initialization-based methods on a number of ccaassccll' machine learning benchmarks (Hochreiter & Schmidhuber, 1997). More importantly, we demonstrate that while with previous methods was difficult to capture slow behavior in DS that exhibits widely different time scales, our new regularizatioonssupprre inference efficiently captures all relevant time scales.
RELATED WORK
Long-range dependency problems in RNN. Error gradients in vanilla RNN trained by some form of gradient descent, like back-propagation through time (BPTT, Rumelhart et al. (1986)), 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 al., 1994; Hochreiter & Schmidhuber, 1997). Formally, RNN 2, = Fe (z,-1, 8t) are discrete time dynamical systems that tend either converge, c.g. to fixed point or limit cycle attractors, or diverge (to infinity or as in chaotic systems) over time, unless parameters of the system are precisely tuned to create directions of zero-flow in the system's 's state space (Fig. 1), called line attractors (Seung, 1996; Seung et al., 2000; Durstewitz, 2003). Convergence of the RNN in general implies vanishing and global divergence exploding gradients. To address this issue, RNN with gated memory cells have been specifically designed (Hochreiter & Schmidhuber, 1997: Cho al., 2014), but these are complicated and tedious to analyze from a DS perspective. Recently, Le et al. (2015) observed that initialization of the recurrent weight matrix W to the identity in ReLU-based RNN may yield performance en par with LSTMs on standard machine learning benchmarks. For a ReLU with activity 2, 2 0. zero bias and unit slope, this results in the identity mapping, hence line attractor configuration. Talathi
2

---

Under review as a conference paper at ICLR 2020
& Vartak (2016) expanded this idea by initializing the recurrence matrix such that its largest absolute eigenvalue is 1, arguing that this would leave other directions in the system's 's state space free for computations other than memory maintenance. Later work enforced orthogonal (Henaff et al., 2016; Helfrich et al., 2018; Jing et al., 2019) or unitary (Arjovsky et al., 2016) constraints on the recurrent weight matrix during training. While this appears to yield long-term memory performance superior to that of LSTMs, these networks are limited in their computational power (Kerg et al., 2019). This may be a consequence of the fact that RNN with orthogonal recurrence matrix are quite restricted in the range of dynamical phenomena they can produce, e.g. chaotic attractors are not possible since diverging eigen-directionn direetiressonns are disabled. Our approach therefore to establish line 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 divergence of RNN dynamics is still possible. We furthermore implement these concepts through regularization terms in the loss functions, such that they are encouraged throughout training unlike when only established through initialization.
Dynamical systems reconstruction. From a natural science perspective, the goal of reconstructing the underlying DS fundamentally differs from building system that 'merely' yields good ahead pre- dictions, as in DS reconstruction 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 RNN for DS identification (Roweis & Ghahramani, 2002; Yu al., 2006) mainly focused on inferring the posterior over latent trajectories Z == ( 2r} given time series data X = {xI,..., 2T J. pl Z X). and on ahead predictions (Lu al., 2017), hence did not show that inferred models can generate the underlying attractor geometries on their own. Others (Trischler & D D'Eleuterio, 2016; Brunton et al., 2016) attempt to approximate the flow field, obtained e.g. by numerical differentiation, directly through basis expan- sions or neural networks, but numerical derivatives are problematic for their high variance and other numerical issues (Raissi, 2018; Baydin et al., 2018; Chen et al., 2017). Some approaches assume the form of the DS equations basically to be given (Raissi, 2018; Gorbach et al., 2017) and focus on estimating the system's 's latent states and parameters, rather than approximating an unknown DS based the observed time series information alone. In many biological systems like the brain the intrinsic dynamics are highly stochastic with many noise sources, like probabilistic synaptic release (Stevens, 2003), such that models that do not explicitly account for dynamical process noise (Cham- pion et al., 2019; Rudy et al., 2019) may be less suitable. Finally, some fully probabilistic models for DS reconstruction based on GRU (Fraccaro al. (2016), cf. Jordan et al. (2019)), LSTM (Zheng et al., 2017), or radial basis function (Zhao & Park. 2017) networks are not easily interpretable and amenable to DS analysis. Most importantly, none of these previous approaches considers the long-range dependency problem within more easily tractable RNN for DS reconstruction.
3 MODEL FORMULATION AND OPTIMIZATION APPROACHES
3.1 MODEL AND PRELIMINARIES
Assume we are given two multivariate time series S = { and X = {x one we will denote as 'inputs' (S) and the other as 'outputs' (X). We will first consider the 'classical' (supervised) machine learning setting where we wish to map S on X through RNN with latent state equation
2, = Fg (z.-188t) as for instance in the "addition problem' (Hochreiter & Schmidhuber, 1997). In DS reconstruction, in contrast, we usually have a 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.
The latent RNN we consider here takes the specific form
2 = A + Wo(zz-1) + Css, + h h Et, E, ~ NO E),
(1)
x1 is the hidden state (column) vector of dimension M, A E RM M a diagonal and W E MMX an off-diagonal matrix, E IR 1xx1 the external input of dimension K, C M the imput mapping, h E Mx a bias, and E a Gaussian noise term with diagonal covariance ma- trix diag(E) Σ) E R.M The nonlinearity σ(z) is ReLU, Φ(z); max(0, 2,) (1....,M] M). This specific formulation is originally motivated by firing rate (population) models in computa- tional neuroscience (Song al., 2016; Durstewitz, 2017), where latent states 2, may represent
where 2, E RM
3

---

Under review as a conference paper at ICLR 2020
membrane voltages currents, A the neurons' passive time constants, W the synaptic coupling among neurons, and ((-) the voltage-to- tottte transfer function. However, for a RNN in the form 2, W (z,-)) + h. note that the simple change of variables y/ W 10-1(a [zz h) will yield the more familiar form yt = á (Wy y h) (Beer, 2006).
Besides its neuroscience motivation, note that by letting A = I, W 0. h 0, get strict line attractor system across the variables' whole support which we conjecture will be of advantage for establishing long short-term memory properties. Also can solve for all of the system's fixed points analytically by solving the equations xt (I A Dn) h, with Dg as defined in Suppl. 7.1.2, and can determine their stability from the cigenvalues of matrix A+ + W WD Dz. We could do the same for limit cycles, in principle, which are fixed points of the r-times iterated map Fjj al- though practically the number of configurations to consider increases exponentially as 2 2 Finally, we remark that a discrete piecewise- linear system can, under certain conditions, be transformed into an equivalent continuous-time (ODE) piecewise-linear system Ç = (C()()) s(t) (Suppl. 7.1.2, Ozaki (2012)); in the sense that if ((t) = 2t. then Ç(t + At) = 2t+1 after a defined time step At. These are among the properties that make PLRNNs more amenable to rigorous DS analysis than other RNN formulations.
We will assume that the latent RNN states are coupled to the actual observations æ, through simple observation model of the form
t = Bz( + ne T ~ I)
in the case of real-valued observations x, E Nx1, where B E NR a factor loading matrix and dTag(T) R. the diagonal covariance matrix of the Gaussian observation noise,
i = fi (5., 1) = (()))((BBBB)))) (c:+.))))
Σ
Br
in the case of multi-categorical observations Fi,t E (0,1), E =1.
3.2 REGULARIZATION APPROACH
We start from a similar idea as Le et al. (2015), who initialized RNN parameters such that it performs an identity mapping for 2,,t 2 O. However, 1) we use a neuroscientifically motivated network architecture (eq. 1) that enables the identity mapping across the variables whole support, 2,, [-o0, +o0), 2) we encourage this mapping only for a subset M S 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.
That way, we 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 RNN model (eq. 1). Specifically, the following penalty is added to the loss function (Fig. S1):
Moex
MERK
Mesk
I Trog TA Σ (Ari-17 1)2 1) TW
W2
+ Th bg
12
While this formulation allows us to trade off, for instance, the tendency toward line attractor (A + I, h + 0) vs. the sensitivity to other units" inputs (W - 0), for all experiments performed here common value, TA TW Th J, was assumed for the three regularization factors.
3.3 OPTIMIZATION PROCEDURE FOR MACHINE LEARNING BENCHMARRS
For comparability with other approaches like LSTMs (Hochreiter & Schmidhuber, 1997) or iRNN (Le et al., 2015), will assume that the latent state dynamics eq. are deterministic (i.e., Σ = 0), will take g(x,) = 2, and T = In in eq. 2 (leading to an implicit Gaussian assump- tion with identity covariance matrix), and will stochastic gradient descent (SGD) for training
4

---

Under review as a conference paper at ICLR 2020
- ((r)) between esti- to minimize the squared-error loss across R samples, L = En 21 mated and actual outputs for the addition and multiplication problems, and the cross entropy loss C= -> (- I-E Fir log( log(P,,T) for sequential MNIST, to which penalty eq. was added for the regularized PLRNN (fPLRNN). We used the Adam algorithm (Kingma & Ba, 2014) from the PyTorch package (Paszke et al., 2017) with learning rate of 0.001, a gradient clip parameter of 10, and batch size of 16. In all cases, SGD is stopped after 100 epochs and the fit with the lowest loss across all epochs is chosen.
3.4 OPTIMIZATION PROCEDURE FOR DYNAMICAL SYSTEMS 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 predicting future values in time series (ef. sect. 3.5). This point has important implications for the design of models, inference algorithms and performance metrics if the primary goal is DS reconstruction rather than 'mere' time series forecasting. In this context we consider the fully probabilistic, generative RNN eq. 1.
Together with eq. 2 (where take g(z)) of z,)) this gives the typical form of nonlinear state space model (Durbin & Koopman, 2012) with observation and process noise. We solve for the pa- rameters 0 {A, W. T.C, h. B. I) by maximum likelihood, for which efficient Expectation- Maximization (EM) algorithm has recently been suggested (Durstewitz, 2017; Koppe et al., 2019), which will briefly summarize here. Since the involved integrals are not tractable, start off from the evidence-lower bound (ELBO) to the log- ligllikei which can be rewritten in various useful ways:
log p( X x00) 2 Ez [log po ( X, Z) + H (1( (|))) =logpt X 0) DKL (7( X I|pp (Z)))) = ((.gg)
(5)
In the E-step, given a current estimate 0" for the parameters, we seek to determine the posterior Pe ( which we approximate by a global Gaussian g x)) instantiated by the maximizer (mode) Z of pe (Z as an estimator of the mean, and the negative inverse Hessian around this maximizer as an estimator of the state covariance, i.c.
Z X R Z"
arg max log Z X) arg max [log X 2) + log ()) log (X))
z
z
= arg max [log Po X Z log Po Z)],
(6)
since Z integrates out in Po Po (X) (equivalently, this result can be derived from a Laplace approxima- tion to the liglliiel,hood, log p( X (0)= ~ log [X +
(Z*) 2 log L' I +onstt 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 (sce Suppl. 7.1.3,
Pog po
Durstewitz (2017); Koppe et al. (2019)).
Using this approximate posterior for (a( X). based on the model's piecewise-lineat structure most of the expectation values Ezny [0(2)] Ezny [=) and Ez~g [o(2) (6(z)) )))())))))) T). could be solved for isemi-)analltiin
(where the concatenated vector form of Z, in Suppl. 7.1.3). In the M- step, we seek l == arg max L(O, (0044) ) assuming proposal density 9º to be given from the E-step,
which for Gaussian observation model amounts simple linear regression problem (see Suppl. eq. 23). To force the PLRNN to really capture the underlying DS in its governing equations, we use previously suggested (Koppe 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 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.
3.5 PERFORMANCE MEASURES
Measures of prediction error. For the machine learning benchmarks we employed the same criteria as used for optimization (MSE or cross-entropy, sect. 3.3) as performance metrics, evaluated across left-out test sets. In addition, report the relative frequency Prorrect of correctly predicted trials
5

---

Under review as a conference paper at ICLR 2020
A
Addition probiem
B
Multiplication problem
C
Sequential MINIS
2.5
10-
10
10-3
10"7
catzopy
MSE:
10
10
1.5
10
10
10"9
10
20
30
4050
100
400
4050
70
400
-RNN
IRNV
mpRMN
PLANV
iPLRNN
#PLRN
LSTM
LaBNN
EPRNN
PLRNN
iPLANN
IPLRNN
1STM
12HNN
Perret
0.8
0.6
Porret
0.6
0.4
0.2
0.2
0.2
20
30
4050
70
100
150 250
400
20
30
40.0
70
100
150 250
400
#RNN IRNV
spRiVN
PLRNN
iPLRNN
PLERN
LSTM
LORNN
Fig. 2: Comparison of rPLRNN (T = 5, Mrer 0.5, cf. Fig. S3) to other methods for A) addition problem, B) multiplication problem and C) sequential MNIST. Top row gives loss a function of time series length T (error bars SEM), bottom row shows relative frequency of correct trials. Dashed lines indicate chance level, black dots in C indicate individual repetitions the experiment.
M
across the test set. A correct trial in the addition and multiplication task is defined as an absolute prediction error smaller than 0.04 (analogous to Talathi & Vartak (2016)), while correct trial in the sequential MNIST data set is defined as one for which the largest probability p:- = max Prr
indicated the correct class 2,- .
Agreement in attractor geometries. From DS perspective, it is not sufficient judge a method's ability to infer the underlying DS purely based on some form (ahead- ppreddctioni error like the MSE defined on the time series itself (Ch.12 in Kantz & Schreiber (2004)). Rather, we require that the inferred model can freely reproduce (when no longer guided by the data) the underlying attractor geometries and state space properties. This is not automatically guaranteed for a model that yields agreeable ahead predictions on time series. Vice versa, if the underlying at- tractor is chaotic, with tiny bit of noise even trajectories starting from the will quickly diverge and thead-predicctinn eπтоrs are not even meaningful as performance metric (Fig. S2A). To quantify how well an inferred PLRNN captured the underlying dynamics we there- fore followed Koppe et al. (2019) and used the Dullback-Leible divergence between the true and reproduced probability distributions across states in state space, thus assessing the agreement in tractor geometries (cf. Takens (1981); Sauer et al. (1991)) rather than in precise matching of time series,
sensible to
initial condition
K Plad (x) log
Perue (x)
DEL. (Porue (PPol (2)2)) R
D (x)zz)
where Ptrue (x) is the true distribution of observations across state space (not time!), (x|=)) is the distribution of observations generated by rumning the inferred PLRNN, and the sum indicates spatial discretization (binning) of the observed state space (see Suppl. 7.1.4 for more details). We
emphasize that Pgen (x) is obtained from freely simulated trajectories, i.e. drawn from the prior p(z), not from the inferred posteriors p(z Ptraia). (The form of p(z) is given by the dynamical Gaussians' structure, see Koppe al. (2019).) In
(k)
model eq. and has 'mixture of piecewise-Gauussians addition, to assess reproduction of time scales by the inferred PLRNN, we computed the average correlation between the power spectra of the true and generated time series.

---

Under review as a conference paper at ICLR 2020
Table 1: Overview over the different models used for comparison
NAME
DESCRIPTION
RNN
Vanilla ReLU based RNN
L2RNN
Vanilla ReLU RNN with standard L2 regularization on all weights
iRNN
RNN with initialization W, = and ho = (Le et al., 2015)
npRNN
RNN with weights initialized to a normalized positive definite matrix with largest eigenvalue of and biases initialized to zero (Talathi &
Vartak, 2016)
PLRNN
PLRNN as given in eq. (Koppe et al., 2019)
iPLRNN
PLRNN with initialization Ao I, Wo- and hg 0
rPLRNN
PLRNN initialized as illustrated in Fig. S1, with additional regulariza- tion term (eq. 4) during training
LSTM
Long Short-Term Memory (Hochreiter & Schmidhuber, 1997)
4 NUMERICAL EXPERIMENTS
4.1 MACHINE LEARNING BENCHMARKS
We compared the performance of our rPLRNN other models on the following three benchmarks requiring long short-term maintenance of information (as in Talathi & Vartak (2016) and Hochreiter & Schmidhuber (1997)): 1) The addition problem of time length T consists of 100 000 training and 10 000 test samples of 2 x T input series S = {8),..., ST)] where entries 81 E [0, 1] are drawn from uniform random distribution and 82,: E 1) contains zeros except for two indicator bits placed randomly at times 1 < 10 and tz < T/2. Constraints on t1 and f/ are chosen such that every trial requires a long memory of 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 $1. indicated by the 1-entries in $2.;, 27 target $1,6, $1,t 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, ZTT dataset (LeCun & Cortes, 2010) consists of 60 000 training and 10 000 28 28 test images of hand written digits. To make this 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.
target 81,t, 181,t 3) The MNIST
On all three benchmarks we compare the performance of the rPLRNN (eq. 1) to several other models summarized in Table 1. To achieve meaningful comparison, all models have the same number of hidden states M, except for the LSTM, which requires three additional parameters for each hidden state and hence has only M/ hidden states, yielding the overall same number of trainable parame- ters as for the other models. In all cases, M = 40, which initial numerical exploration suggested to be a good compromise between model complexity (bias) and data fit (variance) (Fig. S3).
Fig. 2 summarizes the results for the machine learning benchmarks. As can be seen, on the addition and multiplication tasks, and in terms of either the MSE percentage correct, our rPLRNN outper- forms all other tested methods, including LSTMs. Indeed, the LSTM performs even significantly worse than the iRNN and the iPLRNN. The large error bars in Fig. result from the fact that the networks mostly learn these tasks in an all-or-none fashion, i.c. either learn the task and succeed in almost 100 percent of the cases or fail completely. The results for the sequential MNIST problem are summarized in Fig. 2C. While in this case the LSTM outperforms all other methods, the rPLRNN is almost en par with it. In addition, the iPLRNN outperforms the iRNN Similar results were obtained for M 100 units (M 25, respectively, for LSTM; Fig. S6). While the rPLRNN in general outperformed the pure initialization-banl based models (iRNN, npRNN, iPLRNN), confirming that a line attractor subspace present at initialization may be lost throughout training, conjecture that this difference in performance will become even more pronounced as noise levels or task complexity increase.
7

---

Under review as a conference paper at ICLR 2020
A
B
0.02
0.02
50-000 3t
0.015
normoklioed MSE
0.015
1460 Be
MSE
0.01
0.005
0.000
200
D
gemeraced
0.15
0.3
01:45
0.6
0.15
0.3
0.45
0.6
0.15
0.3
0.45
0.5
(4)
(+)
Fig. 3: Reconstruction of 2-time scale DS (biophysical bursting neuron model) in limit cycle regime. A) KL divergence (DKL) between true and generated state space distributions as a func- tion of т. Unstable (globally diverging) system estimates were removed. B) Average MSE be- tween power spectra of true and reconstructed DS. C) Average normalized MSE between power spectra of true and reconstructed DS split according low (S 50 Hz) and high (> 50 Hz) fre- quency components. Error bars = SEM in all graphs. D) Example of (best) generated time series (red=reconstructiin with r 2
4.2 NUMERICAL. EXPERIMENTS ON A DYNAMICAL SYSTEM WITH DIFFERENT TIME SCALES
Here our goal was to examine whether our regularization approach would also help with the iden- tification of DS that harbor widely different time scales. By tuning systems in the vicinity of line attractors, multiple arbitrary time scales can be realized in theory (Durstewitz, 2003). To test this, we used a biophysically motivated (highly nonlinear) bursting cortical neuron model with one voltage and two conductance recovery variables (see Durstewitz (2009)), one slow and one fast (Suppl. 7.1.5). Reproduction of this DS challenging since it produces very fast spikes on top of slow nonlinear oscillation (Fig. 3D). Time series of standardized variables of length T 1500 were generated from this model and provided as observations to the rPLRNN inference algorithm. PPLRNN with M {8... ... 118 states were estimated, with the regularization factor varied within T E (0, 10 102, 103, 10a, //5000.
Fig. 3A confirms our intuition that stronger regularization leads to better DS reconstruction sessed by the KL divergence between true and generated state distributions. This decrease in DRL. accompanied by likewise decrease in the MSE between the power spectra of true (Suppl. eq. 27) and generated (rPLRNN) voltage traces as T increased (Fig. 3B). Fig. 3D gives an example of volt- age traces and gating variables freely simulated (i.e., sampled) from the generative rPLRNN trained with r 3. illustrating that our model is in principle capable of capturing both the stiff spike dy- namies and the slower oscillations in the second gating variable 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 benefitted from increasing т, there was a strong re- duction in the MSE computed on the power spectrum for the lower frequency range (S 50 Hz), suggesting that increased regularization helps to map slowly evolving components the dynamics.
CONCLUSIONS
In this work we have introduced a simple solution to the long short-term memory problem RNN that on the one hand retains the simplicity and tractability of vanilla RNN, yet on the other hand does not curtail the universal computational capabilities of RNN (Koiran et al., 1994; Siegelmann & Sontag, 1995) and their ability approximate arbitrary DS (Funahashi & Nakamura, 1993; Kimura & Nakano, 1998; Trischler D' Eleeuerrio, 2016). We achieved this by adding regularization terms to the loss function that encourage the system to form 'memory subspace", that is, line attractor
8

---

Under review as a conference paper at ICLR 2020
dimensions (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 constraint which has important implications for capturing slow time scales in the data: It allows the RNN to slightly depart from a perfect line attractor, which has been shown to constitute a general dynamical mech- anism for regulating the speed of flow and thus the learning of arbitrary time constants that not naturally included qua RNN design (Durstewitz, 2003; 2004). This is because as we come infinites- imally close to a line attractor and thus bifurcation in the system' parameter space, the flow along this direction becomes arbitrarily slow until it vanishes completely in the line attractor configuration (Fig. 1). Moreover, part of the RNN's latent space was not regularized all. leaving the system enough degrees of freedom for realizing arbitrary computations or dynamics. We showed that the rPLRNN is en par with or outperforms initialization-base of classical benchmarks, and, more importantly, that the regularization strongly facilitates the iden- tification of challenging DS with widely different time scales in ELRNN-basee algorithms for DS reconstruction. Future work will explore wider range of DS models and empirical data with diverse temporal and dynamical phenomena. Another future direction may be to replace the EM algorithm by black-box variational inference, using the paaammtttrriztionn trick for gradient descent (Kingma & Welling, 2013; Rezende et al., 2014; Chung et al., 2015). While this would come with better scaling in M, the number of latent states (the scaling in T is linear for EM well, see Paninski et al. (2010)), the EM used here efficiently exploits the model's 's piecewise linear structure in finding the posterior over latent states and computing the parameters (see Suppl. 7.1.3). It may thus be more accurate and suitable for smaller-scale problems where high precision is required, as often encountered in neuroscience or physics.
approaches and LSTMs on a number
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 GK (TRR265: A06 & B08). We would like to cordially thank Dr. Zahra Monfared for her careful reading of the manuscript and her thoughtful suggestions.
REFERENCES
Martin Arjovsky, Amar Shah, and Yoshua Bengio. Unitary evolution recurrent neural networks. In Proceedings of The 33rd International Conference on Machine Learning, 2016. URL http:
/ //proceedings , mlr. www.mm...eees...... press/v43 arjovsky16 html.
Ibrahim Ayed, Emmanuel de Bézenac, Arthur Pajot, Julien Brajard, and Patrick Gallinari. Learning Dynamical Systems from Partial Observations. arXiv preprint, 2019. URL http: 1902. 11136.
Atilmm Güneş Baydin, Barak A. Pearlmutter, Alexey Andreyevich Radul, and Jeffrey Mark Siskind. Automatic differentiation in machine learning: A survey, Journal of Machine Learning Research, 18:1 43,, feb 2018 ISSN 15337928. URL /tretp://axxvv..oot org abs 11020 05767.
Randall D. Beer. Parameter space structure of continuous-time recurrent neural networks. Neural computation, 18(12):3009 51, 2006. doi: 10.1 22.1166//00.22006188..1 ncbi. nlm. mih. http://uummmd/d//eee/e//
URL http: www.
Yoshua Bengio, Patrice Simard, and Paolo Frasconi. Learning long-term dependencies with gradient descent is difficult. IEEE Transactions Neural Networks, 5(2):157-166, 1994. doi: 10. 109/ 72.279181. URL https: /oii.ooooooo//oooo//// www.ol.g//.o.oooooooooo
Steven L. Brunton, Joshua L. Proctor, J. Nathan Kutz, and William Bialek. Discovering governing equations from data by sparse identification nonlinear dynamical systems. Proceedings of the National Academy of Sciences of the United States of America, 133155)3932-3337, 2016. doi: 10.10733pnas.15111 15133841111177813 URL pattps:/ /www//wwwww. pnas /ttp///ottt//////////////
Kathleen Champion, Bethany Lusch, J. Nathan Kutz, and Steven L. Brunton. Data-driven discovery of coordinates and governing equations. arXiv preprint, 2019. URL http: /ttp/viivvvo/// /rg// abs/1904 02107.

---

Under review as a conference paper at ICLR 2020
Shizhe Chen, Ali Shojaie, and Daniela M. Witten. Network Reconstruction From High- High-Dimensiioal Ordinary Differential Equations. Journal of the American Statistical Association, 1126520):1697- 1707, 2017. doi: 001108001621159922066122299 01 6261499.. 2 01 6. 12291 97.
108C/ URL http://
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. In Proceedings of the Conference on Empirical Methods in Natural Language Processing, pp. 1724-1734. Association for Computational Linguisties (ACL), 2014.
Junyoung Chung. Kyle Kastner, Laurent Dinh, Kratarth Goel, Aaron C Courville, and Yoshua Ben- gio. recurrent latent variable model for sequential data. In C. Cortes, N. D. Lawrence, D. D. Lee. M. Sugiyama, and R. Garnett (eds.). Advances in Neural Information Processing Systems 28, pp. 2980-2988. Curran Associates, Inc., 2015. URL http: pppers nips cc cc/cccc/pa//oco /ptpp//pe//////////// 5653- a- recurrent latent- venttt-aribblee model - for sequential iatieaattaalaaaaaaaaaaaaaaaaaaaa pdf.
Lea Duncker, Gergo Bohner, Julien Boussard, and Maneesh Sahani. Learning interpretable continuous-time models of latent stochastic dynamical systems. arXiv preprint, 2019. URL. http://arxiv..rrt abs
James Durbin and Siem Jan Koopman. Time Series Analysis by State Space Methods. Oxford University Press, May 2012. doi: 10. 10.1099/4cmcdas000/790999944 001.17....000000 1,001.. httpp::: http:////.....//////rooooo:oooo. oso/ /tt.::::////77109111111117
wr..org///lo 10. 1093/ acprof
001. 0001.
Daniel Durstewitz. Orga-Ongg Neural Integrator Predicts Interval Times. Journal of Neuro- science, 23(12):5342-5353, 2003. URL https:tt//////t http::://wwww.. http:p://www...nossccchccccc
12/5342.
Daniel Durstewitz. Neural representation of interval time. NeuroReport, 15(5)7445-749, April 2004. doi: 044197/00001156-20440990--000000 40001755-2004040900000
URL httpss///ac...og//......
Daniel Durstewitz. Implications of synaptic biophysics for recurrent network dynamics and active
memory. Neural Newwork,, (28)):189-12000 oct 2009, doi: 10.10166/j www.fj.eeee..0000..000.0 URL https: wwc.lo..orrg/ooo 1016/115510
www.nentt........ 2009. 10.09000700 07. 016.
Daniel Durstewitz. A State Space Approach for Piecewise-Linear Recurrent Neural Networks for Reconstructing Nonlinear Dynamies from Neural Measurements. PLoS Computational Biology, 1306))e0005542, 2017. doi: 10. Iiii 1371/journal pcbi..cii....... 1005542.
URL https: /oo. /ttp::::////ii.o/////////
Marco Fraccaro, Søren Kaae Sønderby, Uirich Paquet, and Ole Winther. Sequential neural models with stochastic layers. Advances in Neural Information Processing Systems, 2016.
Ken-ichi Funahashi and Yuichi Nakamura. Approximation of dynamical systems by continuous time doi: 10. (11///0083-60000000 80125-x. URL htt.s:tpt//// /ttpas:///o... org/10 1 01 6/ /0/809933 933 6080 (05) 80125-x.
recurrent neural networks. Neural Nenworks, 666)8801-806, 1993.
Nico S. Gorbach, Stefan Bauer, and Joachim M. Buhmann. Scalable variational inference for dynamical systems. In Guyon, U. V. Luxburg, S. Bengio, H. Wallach, R. Fergus, S. Vish- wanathan, and R. Garnett (eds.), Advances in Neural Information Processing Systems 30. pp. 4806-4815. Curran Associates, Inc.. 2017. URL http::/ 7066-sca www.aaaaaeeleaaaaiieaaoaaaaaaaae
nips cc/paper/ wore--nineeenee---oo-deeoeeeeee- -yssssttsssssssssssss.m.. .ddf.
Kyle E. Helfrich, Devin Whimott, and Qiang Ye. Orthogonal recurrent neural networks with scaled Cayley transform. 35th International Conference Muchine Learning, ICMI. 2018, 5, 2018.
Mikael Henaff, Arthur Szlam, and Yann LeCun. Recurrent orthogonal networks and long- tasks. Maria Florina Balcan and Kilian Q. Weinberger (eds.), Proceedings of The 33rd In- ternational Conference on Machine Learning, volume 48 of Proceedings of Machine Learn- ing Research, pp. 2034-2042, New York, New York, USA, 20-22 Jun 2016. PMLR. URL
http:/ ppwpppooceddiiee ngs ml wal.r.prssss.... v48 00/4///eaaa fffefff1. h.mll.
10

---

Under review as a conference paper at ICLR 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 Sig- nal Processing ICASSP '07, volume 4, pp. IV-317-IV-320, April 2007. doi: 10.1109/CCSSP. 2007.366913. URL https:: /oo//ddooooooooooooo oog.ioorg..... 1109/ TCASSP 2007 2007.366 3666913
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): doi: 10.1 1e2////co.199997 1997.9.8.77355. URL http://www.
1735-80, 1997. ISSN 0899-7667. ncbi. nlm. nih. gov gov/pu/mee////// 2oibmed/9377271 9377276.
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, apr 2019. ISSN 1530888X. doi: 10.1162/neco.a.0111
lan D. Jordan, Piotr A. Sokol, and Memming Park. Gated recurrent units viewed through the lens of continuous time dynamical systems. arXiv preprint, abs/1906.01005, 2019. URL http:
/ / /arxivvv org/ /tapg//bb///ee /066. .0....1000 01005.
Holger Kantz and Thomas Schreiber. Nonlinearr Time Series Analysis. Cambridge University Press, 2. edition, 2004, doi: 10. 10.1017CBB0978011155777 URL https:://ci..orgg
CBO978051175579B.
Giancarlo 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.
URL http: /tap:://aaxxii org/ 1955 12080.
M. Kimura and R. Nakano. Learning dynamical systems by recurrent neural networks from orbits. Neural Networks, 1199:1589-1599 1998. doi: (000////8833-000988800000000 //doi org/1/.r111..1 11w1/1011110010 6080 ( (98) (99800000 00098-7.
URL https:
Diederik P. Kingma and Jimmy Ba. Adam: Method for Stochastic Optimization. arXiv preprint, 11412. 6980.
2014. URL http::tt//// http:///xxvvv..g///////
Diederik P Kingma and Max Welling. Auto-Encoding Variational Bayes. (MI):1-14, 2013. URL
http: /ttt:://arxiv..oo ht/p/aarvvvvv... abs 11122 6114.
Pascal Koiran, Michel Cosnard, and Max H. Garzon. Computability with low-dimensional dynami- cal systems. Theoretical Computer Science, 122:113-122 1994.
Georgia Koppe, Hazem Toutounji, Peter Kirsch, Stefanie Lis, and Daniel Durstewitz. Identifying nonlinear dynamical systems via generative recurrent neural networks with applications IMRI. PLOS Computational Biology, 1508):e000726,, 2019, doi: 10. 1...33//jornal.pp/ 1ct.....c........... ht https: //t::://ooioe// orgg 1371/ 1tt.://.1n/////r/aa////aa journal.....o pcbi.. 1007263.
1007263. URL
Quoe V. Le, Navdeep Jaitly, and Geoffrey E. Hinton. A Simple Way Initialize Recurrent Networks of Rectified Linear Units. arXiv preprint, 2015. URL http: http://aaaxi...oooob///b//// 00941.
Yann LeCun and Corinna Cortes. MNIST handwritten digit database. 2010. URL ARLhttp://yyn.
lecun. com/ weeeee/mniiiiii
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. doi: 10.1063/1 4979665.
Tohru Ozaki. Time Series Modeling of Neuroscience Data. CRC Press, 2012. doi: 1201/// 1527. URL https: http:///:ooooooo/././oo wwt.://///........... 0b/1152
Liam Paninski, Yashar Ahmadian, Daniel Gil Ferreira, Shinsuke Koyama, Kamiar Rahnama Rad, Michael Vidne, Joshua Vogelstein, and Wei Wu. A new look at state-space models for neural data. Journal of Computational Neuroscience, 29(1-2):107-126, 2010. ISSN 03295313. doi:
10. 10.1007741187770000017991

---

Under review as a conference paper at ICLR 2020
Adam Paszke, Sam Gross, Soumith Chintala, Gregory Chanan, Edward Yang, Zachary Vitito, 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 of nonlinear partial differential equa- tions. Journal of Machine Learning Research, 19:1-24, 2018. URL http::/jmlz... papers/v19/15 1 18tm4..t.m.....
.tmll.
Danilo Jimenez. Rezende, Shakir Mohamed, and Daan Wierstra. Stochastic Backpropagation and Approximate Inference in Deep Generative Models. 2014. URL http://arxiii orgg / 1401. 4082.
Sam Roweis and Zoubin Ghahramani. Learning Nonlinear Dynamical Systems Using the Expeetation-Maxiiizatiioe Algorithm. In Kalman Filtering and Neural Networks, chapter 6, pp. 175-220. Wiley- Blackwel, 2002. ISBN 9780471221544 doi: 10. 10.100220471221146..hh https: /t////nnnnniraaeeeeeeeeeeee wiley- com/ doii http///a/111.100000/21110
URL
/0t......00000002225555555 /ttpa//......//////////////////0 ch..
Samuel Rudy, Steven Brunton, and Nathan Kutz. Smoothing and parameter estimation by
soft-adherence to governing equations. Journal of Computational Physics, 398: 108660, Decem- ber 2019. doi: j0.1016//..pp.2019.. 108860. URL https: /ttp::://doiooorg//////
jap.
2019.108860
David E. Rumelhart, Geoffrey E. Hinton, and Ronald J. Williams. Learning representations by back-propagatiing errors. Nature, 323:533-536, 1986. doi: 10.1038/323533a0. URL https:
Tim Sauer, James A. Yorke, and Martin Casdagli. Embedology. Journal of Statistical Physics, 65
(3-4):579 616, November 1991. doi: 10.1007/6601053745.. URL https:://oi.orrtt
10C7/bf01053745.
H. Sebastian Scung. How the brain keeps the eyes still. Proceedings of the National Academy of Sciences, 13333):13339-13344 1996. doi: 10.1073/pnas.93.22 93. 23. 13339. URL https: //ts::///ol.
11tg//1.103//////// 1073/0077asss 93.23.0 13339.
H. Sebastian Seung, Daniel D. Lee, Ben Y. Reis, and David W. Tank. Stability of the memory of eye position in recurrent network of conductance-based model neurons. Neuron, 2000. doi: 10. URL https: /ttp:://doii.
81155-1.
Hava T. Siegelmann and Eduardo D. Sontag. On the computational power of neural nets. Journal of Computer and System Sciences, 5301)113 February 1995. doi: 10. 19906//ss.1995.101.. 1913. URL https: /ttp:://doaa www.oooogoo....oooo 1tt.oo...ooosssseeee////oos
1013.
H. Francis Song, Guangyu R. Yang, and Xiao Jing Wang. Training Excitatory-Inhibice Recurrent Neural Networks for Cognitive Tasks: A Simple and Flexible Framework. PLoS Computational Biology, 12(2):1-30, 2016. doi: 10.137 //o.////lllll...i....oo 1004999 URL https: /doi.org////oooooooooooooo.
10. 18.1137//00 /037////oraa....m wwtppnaa.ppp.b.....m ptb.: 1004792.
Charles Stevens. Neurotransmitter release at central synapses. Neuron, 40(2):881-388, Oc- tober 2003. doi: 10.101640896-6233000064444 62796-6773 (03) 00643-3.
URL https::/aci.org/////
Steven H. Strogatz. Nonlinear Dynamics and Chaos: Applications to Physics, Biology, Chem- istry, and Engineering: With Applications Physics, Biology, Chemistry and Engineering.
CRC Press, 2015. doi: 10. 10.1201/978042999955..
URL fttps:://ci... http::///cc..org11..1000
9780429492563.
Floris Takens. Detecting strange attractors in turbulence. In Lecture Notes in Mathematies, pp- 366-381. Springer Berlin Heidelberg, 1981. doi: 10.1007/bth00011924 URL https:://o/o. 10.1100000//00 1tt.11007/bf0001992000
Sachin S. Talathi and Aniket Vartak. Improving performance of recurrent neural network with ReLU nonlinearity. ICLR Workshop submission, 2016. URL https: / arxiv. wrw.or///aao... aarg//bss/5111 03771.
12

---

Under review as a conference paper at ICLR 2020
Adam P. Trischler and Gabriele M.T. D'Eleuterio. Synthesis of recurrent neural networks for dynam- ical system simulation. Neural Networks, 80:67-78, 2016. doi: n0.1016fj.ecumet2216.. 201.......................
URL http: / doi. 10 1ttp://...11/..... neunet. wtt.uneet2205......... 04 .011.
Byron M. Yu, Afsheen Afshar, Gopal Santhanam, Stephen I. Ryu, Krishna V. Shenoy, and Maneesh Sahani. Extracting dynamical structure embedded in neural activity. In Y. Weiss, B. Schölkopf. and C. Platt (eds.). Advances in Neural Information Processing Systems 18, pp. 1545-1552. MIT Press, 2006. URL http://papent nips. cc cc/paper/ 2 ingeetiainnddnmmmmmmmmmmmm cal wwwerettreeemmeeddd------e- pdf.
Yuan Zhao and II Memming Park. Variational Joint Filtering. arXiv preprint, 2017. URL http:
//aa///aaa/ //arxiv..rog wwg.wvo.r./aooooo
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. URL http://
arxiv. www. abs 11t.11199111
13

---

Under review as a conference paper at ICLR 2020
7
SUPPLEMENTARY MATERIAL
7.1 SUPPLEMENTARY TEXT
7.1.1 Simple exact PLRNN solution for addition problem
The exact PLRNN parameter settings (ef. eq. 1) for solving the addition problem with 2 units (cf.
Fig. IC) are as follows:
A - 0 w ( ( ( -1
B 0)
(8)
7.1.2 Conversion from discrete to continuous time PLRNN
Under some conditions we can translate the discrete into an equivalent continuous time PLRNN. Using Dace) defined below (7.1.3) for single time step L, we can rewrite (ignoring the noise term and inputs) PLRNN eq. I in the form
1++1 F(z)) = W + h, with Wa(e) = A + W Da())
where (2(t) = [m/ Zm, > 0} is the set of all unit indices with activation larger at time t. To convert this into equivalent (in the sense defined in eq. 11) system of (piecewise) ordinary differential equations (ODE), we need find parameters Wa and h,
GGC WnC(() + h,
(10)
such that
20 (()) - = = S(A(),
(11)
where At is the time step with which the empirically observed time series X was sampled. From these conditions it follows that for each of the E (1....,,M) subregions (orthants) defined by fixed index sets S2" C (I....,M. M} we must have
1 Ac CatC((t h dt,
(A + WD =) 20 h
(12)
where we assume that Do is constant for one time step, i.e. between 0 and At. We approach this by first solving the homogeneous system using the general ansatz for systems of linear ODEs,
(A + WD
1) 1
W Σ
oy dt
(13)
b
St
=
Σ
Ck U
cdal dl
(14)
Eondue (Sude-1)
(15)
d
Wo _ Σ Eougetsse
(16)
cdi
V
(17)
Am
W = VAV-1
(18)
where we have used 2g Ex UR lines 15 and 16. Hence infer matrix Wn. from the cigendecomposition of matrix Wa by letting dn = + log Ae. where d are the eigenvalues of Wa', and reassembling
W
V
log log(A)V-
(19)
Al
14

---

Under review as a conference paper at ICLR 2020
We obtain the general solution for the inhomogeneous by requiring that for all fixed points x² = F(*") of the map eq. we have G(z*) = 0. Using this we obtain
h n- W
(20)
Assuming inputs 8, be constant across time step At, can apply the same transformation to input matrix C. Fig. S5 illustrates the discrete to continuous PLRNN conversion for nonlinear oscillator.
Note that in the above derivations we have assumed that matrix Wn- can be diagonalized, and that all its eigenvalues are nonzero (in fact, Wh should not have any negative real eigenvalues). In general, not every discrete-time PLRNN can be converted into continuous-time ODE system in the sense defined above. For instance, we can have chaos in a nonlinear map, while we need at least a 3d ODE system to create chaos (Strogatz, 2015).
7.1.3 More details on EM algorithm
Here briefly outline the fixed- point-iteratioon algorithm for solving the maximization problem in eq. 6 (for more details see Durstewitz. (2017); Koppe et al. (2019)). Given a Gaussian latent PLRNN and Gaussian observation model, the joint density p( X Z) will be piecewise Gaussian, hence eq. 6 piecewise quadratic in Z. Let us concatenate all state variables across m and & into one long column vector
.2........ 2M,1, : 21,7: ZM.T), arrange matrices A, W into large MT MT block tri-diagonal matrices, define da == (1 =1=,1>0 22,100 ... 1 = as an indicator vector with a for all states 2m,t and zeros otherwise, and Dn diag(dg, as the diagonal matrix formed from this vector. Collecting all terms quadratic, linear, or constant in 2, we can then write down the optimization criterion in the form
(x(2) =-11=² (Uo DDU1 + U.Do + DoU,Da) 2 20000 + Davi)
(wo 1 2] + const.
(21)
In essence, the algorithm now iterates between the two steps:
1. Given fixed Dg, solve x" = -Up + DRU1 + +U + DDU,Da) (wo + DQvi)
2. Given fixed 2", recompute Da
until either convergence or one of several stopping criteria (partly skelihooddbassed partly to avoid loops) is reached. The solution may afterwards be refined by one quadratic programming step. Numerical experiments showed this algorithm be very fast and efficient (Durstewitz, 2017; Koppe et al., 2019). At x², estimate of the state covariance is then obtained as the inverse negative Hessian,
== + DRU, + +U + DnU,Da)
(22)
In the M-step, using the proposal density g² from the E-step, the solution to the maximization problem arg max L(O, q²), can generally be expressed the form
0"
Σ E (an
Σ E (3,,77
(23)
where, for the latent model, eq. 1, 04 2, and B, [z] 0( 2, 2+* and for the observation model, eq. 2, a; = and Br = (zt).
7.1.4 More details on DS performance measure
The measure DKI. introduced in the main text for assessing the agreement in attractor geometries only works for situations where the ground truth Purue (X) is known. Following Koppe et al. (2019), here we would like to briefly indicate how proxy for DKL may be obtained in empirical situa- tions where no ground truth is available. Reasoning that for a well reconstructed DS the inferred posterior (1 m)) given the observations should be a good representative of the prior generative
15

---

Under review as a conference paper at ICLR 2020
dynamics Pgen Pgen(2), one may use the Kullback-Leible divergence between the distribution over la- tent states, obtained by sampling from the prior density (z)n(2), and the (data-come constraind)))))) posterior distribution (where & E MM and a E IR Nx1x)), taken across the system's 's state space:
DKL (Pine ( [[) = /
=ERMX1 Pint (a((2/)) log Pint(z|ar) dz
(24)
Pgea (2)
As evaluating this integral is difficult, one could further approximate (=)[(z and Pgen (z) by Gaussian mixtures across trajectories, i.c. Pinr(z x)) x + 1L (T::)) and Pgen (z) R
T, (0-1((21|zz-1) where the mean and covariance of p(z; [(rr) and ((zz)zı11 are obtained by marginalizing
1
the multivariate distributions p Z X) and Pgen (Z), respectively, yielding E[z, | 2| i). and covariance matrices Var(z, .)r) and Var (1 Supplementary eq. 24 may then be numerically approximated through Monte Carlo sampling (Hershey & Olsen, 2007) by
DKL (Pint (z)=(( (Pgen(()))))
2 n2
1 "
Σ
log
Pint (2() xx)
26) (tnr(z)
(25)
Pgen (z(1)
For digh-dimensional state spaces, for which MC sampling becomes challenging, there is luckily variational approximation of eq. 24 available (Hershey & Olsen, 2007):
T DRu (
Dyaciational (pint (z)r((/())111>) (1) (=))
KL
1 T
log j=1
(26)
k= e -
where the KL divergences in the exponentials are among Gaussians for which we have an analytical expression.
7.1.5
More details on single neuron model
The neuron model used in section 4.2 is described by
-Cm V gr FI.)) + (V) - FNa)) gKn(V Fx)
+ gv Ex) + (V) - ENMDA)
(27)
Nox (V) -h
h
(28)
Th Noo (V)
(29)
Tn
(V)
[1 .3ee - -06-666677
(30)
where Cm refers to the neuron's membrane capacitance, the g. to different membrane conductances, E. the respective reversal potentials, and m, h, and n are gating variables with limiting values given by
(Mias Teoy 1100) = [1 V)
(31)
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: C.. 6FF, gr. 8 mS, F.z = -80 mV, INa = 20 mS, Ex = 60 VANa = --20 V. kNa = 15, 9K 10 Ex = 990 mV, VAK = -25 m kx = 5, 1 ms, gM = 25 mS,
VAM = mm KM = 5, Tr = 200 ms, 9NMDA = mS.
16

---

Under review as a conference paper at ICLR 2020
7.2
SUPPLEMENTARY FIGURES
0
0
0
0
0
XXXBOO
BOOXX
0
0
0
0
000
0
0
0
X
0
XXXX
X
0
0
0
0
X
X
A
W
h
Fig. S1: Illustration of the L.2-regularization for the PLRNN's auto-regression matrix A, coupling matrix W, and bias terms h. Regularized values are indicated in red, crosses mark arbitrary values (all other values set to as indicated).
A
B
0
-20
-40
-60
-80
0.1
0.2
0.3
0.1
0.2
0.3
time (6)
time (6)
0.06
0.06
0.04
0.04
0..
0.5
0 50
50
V
V
Fig. S2: MSE evaluated between time series graph (top) and state space (bottom) for the single neuron model (see section 4.2 and Suppl. 7.1.5) 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 time series initially seems to indicate perfect fit. B) Same in A) for two trajectories drawn from exactly the same DS (i.e., same parameters) with slightly different initial conditions. Despite identical dynamics, the trajectories immediately diverge, resulting in high MSE. Dash-dotted grey lines in top graphs indicate the
good measure for DS reconstruction. A) Time
point from which onward the state space trajectories were depicted.
17

---

Under review as a conference paper at ICLR 2020
A
0.16
10"
10"
MSE
0.1
MSE
MSE
10
10-3
10-4
10
10
0.05
10
30
50
0.3
10 20
50
100
0.1
0.3
0.5
0.7
0.9
M
Moas/ MM
Fig. S3: Performance of the rPLRNN the addition problem for different A) numbers of latent states M. B) values of r and C) proportions Mrag/M /M. Dashed lines denote the values used for the results reported in section 4.1
0.8
0.0
non-2og
-0.1
.006 .067 .667 6.67 66.7
Fig. S4: Effect of regularization strength on rPLRNN network parameters (cf. 1) (regularized parameters for states m Mreg. eq. 1, in red). Note that some of the non-regularized network parameters (in blue) appear to systematically change well as is varied.
k 2-
10
20
25
30
35
40
45
50
55
60
2
10
15
20
25
30
35
40
45
50
55
60
t.
Fig. S5: Illustaaioon of conversion of discrete into continuous time PLRNN for a PLRNN emula- tion of the nonlinear van-der-Pos Pol oscillator. Shown are the first two latent dimensions. Red lines: continuous solution; blue circles: discrete solution; black bars: perturbations (external inputs).
18

---

Under review as a conference paper at ICLR 2020
A
Addition problem
B
Multiplication peoblem
C
Sequential MNIST
2.5
outropy
MS1
10-2
10
10
10-
10
20
30
40300
70
100
150
250
400
20
30
40.50
70
100
150
250
400
IRNN
*PLRNN (100)
xPLRNN (40)
LSTM
IRNN
(PLRNN (300)
iPERNN (40)
18T3
0.8
0.8
0.8
Pomen
0.6
Pames
0.8
0.4
0.2
0.2
0.2
0
0
0
20
30
1050
70
100
150
250
400
20
30
4030
70
100
150
250
400
IRNN
IPERAN (100)
iPERNN (10)
LSTM
T
T
Fig. S6: Comparison of xPLRNN (r Ming 0.5, cf. Fig. S3) for M 40 and M 100 to iRNN (M 100) and LSTM (M = 100 (4) for A) addition problem, B) multiplication problem and function of time series length T (error bars SEM),
C) sequential MNIST. Top row gives loss bottom row shows relative frequency of correct trials. Dashed lines indicate chance level, black dots in C indicate individual repetitions of the experiment. Note that the rPL.NNN does not improve for
M 100 vs. M- 40.
19


Reviews:
[{"review_id": "rJxOVbICFS", "reviewer": "AnonReviewer3", "summary": "", "questions": "", "limitations": "", "rating": 6, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "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 regularization scheme to improve the ability of RNNs in capturing long-range dependency in the latent space. The proposed model then uses EM for inference and achieves superior performance in sequence modeling tasks over LSTMs and iRNNs. \n\n+ The motivation of the line attractor is novel and effective. The special RNN model studied in the paper has strong connections with neuro-dynamics models. \n+ The paper is well motivated and clearly written. The illustration about the line attractor is particularly interesting. \n+ Good experimental performance on multiple sequence modeling tasks including addition, multiplication and sequence MNIST.\n\n- The paper is building a generative model for sequences. It’s not clear to me why VAE or variational RNN type of approaches cannot be used in this setting. One might think of replacing the Gaussian prior with more complex distributions. The inference procedure can also be significantly simplified with variational inference.\n- The step-wise annealing together with EM inference is not scalable, which prohibits the model from applying to large-scale sequence modeling tasks. \n\nMinor comments\n- \" All code used in this work is freely available on the github site ... . \" Remove this sentence\n", "extra_scores": {"experience_assessment": "I have published one or two papers in this area.", "review_assessment:_thoroughness_in_paper_reading": "I read the paper thoroughly.", "review_assessment:_checking_correctness_of_experiments": "I carefully checked the experiments.", "review_assessment:_checking_correctness_of_derivations_and_theory": "I carefully checked the derivations and theory."}, "extra_text": {}}, {"review_id": "SklBfmHutH", "reviewer": "AnonReviewer2", "summary": "", "questions": "", "limitations": "", "rating": 3, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "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 presents a regulariser that encourages the formation of line attractors in RNNs. This regulariser works on a neuroscience-motivated formulation of RNN, bringing the Jacobian of the dynamic system close to identity. This paper is well-written, with a good coverage of background material from machine learning to computational neuroscience.  While the alleviating the exploding and vanishing gradient problem for simple RNNs is an interesting direction, I think the empirical results are not sufficient to support claims in the paper.\n\nThe paper starts by criticising initialisation and reparametrisation-based techniques. However, I am not convinced why such methods limit the expressive power of RNNs. In fact, one may argue that initialisation is a milder constraint compared with an explicit regulariser, since regularisation affects the entire learning process. It seems that such initialisation requires less tuning (i.e., just identity) compared with the regulariser weights (a rather wide range of choices). Either a theoretical justification or strong empirical results are required to support this claim. However, both are missing in the current paper. \n\nExperiments:\n\nFirst, the baseline results for even toy problems (e.g., addition) are unclear. Despite the results in the paper show the advantage of the proposed method, direct comparison with results from other papers are missing. For example,  when the T > 150 the results from Le et al, 2015 (A Simple Way to Initialize Recurrent Networks of\nRectified Linear Units) were much better compared with baseline results in the paper. This could be due to the smaller size of the models (40 vs 100 hidden units in Le et al.). For clear and direct comparison in this case, models with comparable size should be used in these experiments. Despite this, I wonder why the performance of addition and multiplication seem even worse than the much smaller model reported in Hochreiter and Schmidhuber 1997? (see table 7 and 8)? \n\nActually, it would be helpful to test the proposed method on more practical tasks such as language (at least synthetic ones) and speech modelling, which would bring more impact on the wider community.\n\nA few technical questions:\n\n- Is the form of eq. 1 necessary, or can the method be adapted for the more standard formulation of RNN used in machine learning? It seems that A can be interpreted as a skip connection\n\n- Eq.3 can simply be referred to as “softmax”\n\n- Please comment on the algorithm in section 3.4 in comparison with more standard variational approaches, such as stochastic variational inference with reparametrisation as in variational RNNs (e.g., Chung, et al., 2015, A Recurrent Latent Variable Model for Sequential Data). Is “stepwise annealling” always necessary?\n\n- Eq. 7 as a measurement of the match between trajectories still depends on z. Is there an additional expectation over p(z)? Is “freely simulated trajectories” the prior over trajectories? If so, what’s the form of this prior?", "extra_scores": {"experience_assessment": "I have published one or two papers in this area.", "review_assessment:_checking_correctness_of_derivations_and_theory": "I assessed the sensibility of the derivations and theory.", "review_assessment:_checking_correctness_of_experiments": "I carefully checked the experiments.", "review_assessment:_thoroughness_in_paper_reading": "I read the paper at least twice and used my best judgement in assessing the paper."}, "extra_text": {}}, {"review_id": "ByxJFNRRYB", "reviewer": "AnonReviewer1", "summary": "", "questions": "", "limitations": "", "rating": 1, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "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": "Overview\nThis paper proposes a type of regularization for recurrent networks, with the goal of encouraging particular dynamical structures (in this case, line attractors) in the dynamics of the networks. The proposed regularization penalty is only applied to a subset of the recurrent units; the motivation for this is to allow neurons not contained in the subset to learn different structures. The paper applies this regularization method on three example machine learning sequence tasks: an addition task, a multiplication task, and sequential MNIST classification, as well as on learning a 2-D dynamical system model of a bursting neuron with two different timescales.\n\nMajor comments\nI have a number of serious concerns about the paper's motivation, logic, and experiments:\n\n- First, the paper motivates the proposed regularization as a way to encourage the network to have line attractor dynamics. In particular, the paper dismisses gated architectures as not being interpretable, stating that LSTMs and GRUs \"are complicated and tedious to analyze from a DS perspective.\" (pg 2). However, there is recent work both theoretical (https://arxiv.org/abs/1906.01005) and empirical (https://arxiv.org/abs/1907.08549) that analyzes these gated architectures from a DS perspective. In particular, these papers demonstrate that LSTMs and GRUs are perfectly capable of learning line attractors. Given that it is possible to analyze gated architectures as dynamical systems, the overall motivation of the paper is much weaker.\n\n- The paper proposes a squared penalty on subsets of the weights in the recurrent network as a way to encourage line attractor dynamics. However, it is not clear to me that this is sufficient. In particular, unless the subset of the network that implements the line/plane attractor is completely disconnected from the rest of the network, then the overall dynamics may not contain a line attractor (the units will interact with the unregularized units). Also, the proposed regularization penalty only penalizes the diagonal elements of the A matrix to be close to 1--but shouldn't the off diagonal elements also be penalized to be close to zero?\n\n- Moreover, the paper makes no mention of the Jacobian of these recurrent networks. The eigenvalues of the Jacobian of the recurrent networks determine the behavior of the linearized system around fixed points--specifically, eigenvalues with real part close to 1 will exhibit slow dynamics (approximate line attractors along those dimensions). It seems to me that a much more natural way of encouraging line attractor dynamics is to place a regularization penalty on the Jacobian itself (which is analytically tricky, but numerically more plausible with modern autodifferentiation software). Regardless, the authors should compare the eigenvalues of the recurrent networks' Jacobian when using their regularization method vs without it. Does the proposed regularization encourage the Jacobian of the resulting networks to have eigenvalues close to 1?\n\n- The paper compares the proposed method with a number of vanilla RNNs with different initializations, and an LSTM. However, a critical missing baseline is simply an RNN with l2 regularization on the weights (standard regularization in the literature). This baseline is important to determine if the proposed regularization simply helps because it is an l2 penalty on the weights (note that none of the other baselines have regularization).\n\n- The paper motivates the method as trying to study line attractor dynamics, but then does not apply them to tasks where line attractors are required. For example, the addition and multiplication tasks require discrete memories, not line attractors. The bursting neuron approximation (2D dynamical system) also does not involve a line attractor. However, there definitely exist tasks both in neuroscience (e.g. sensory integration in decision making, path integration in navigation, etc.) and in machine learning (c.f. https://arxiv.org/abs/1906.10720) that use or require line attractors. The motivation of the paper would be much better tested on these tasks.\n\nMinor comments\n- The authors comment at the beginning of page 4 that by setting A=I, W=0, and h=0, the network contains a line attractor, but the more precise language would be to state that the network contains an N-dimensional plane attractor, where N is the number of units. Typically, 'line attractor' refers to a 1-dimensional manifold of fixed points along which the system can integrate inputs, but perturbations off of the line attractor are not remembered (decay back to the line attractor).", "extra_scores": {"experience_assessment": "I have published one or two papers in this area.", "review_assessment:_thoroughness_in_paper_reading": "I read the paper thoroughly.", "review_assessment:_checking_correctness_of_experiments": "I carefully checked the experiments.", "review_assessment:_checking_correctness_of_derivations_and_theory": "I carefully checked the derivations and theory."}, "extra_text": {}}]

Author rebuttal:

Meta-review: The paper proposes  an interesting idea to leave a very simple form for piecewise-linear RNN, but separate units in to two types, one of which acts as memory. The "memory" units are penalized towards the linear attractor parameters, i.e. making elements of $A$ close to 1 and off-diagonal of $W$ close to $1$. 
The benchmarks are presented that confirm the efficiency of the model.
The reviewer opinion were mixed; one "1", one "3" and one "6"; the Reviewer1 is far too negative and some of his claims are not very constructive, the "positive" reviewer is very short. Finally, the last reviewer raised a question about the actual quality on the results. This is not addressed. Although there is a motivation for such partial regularization, the main practical question is how many "memory neurons" are needed. I looked through the paper - this addressed only in the supplementary, where the value of $M_{reg}$ is mentioned (=0.5 M). For $M_{reg} = M$ it is the L2 penalty; what happens if the fraction is 0.1, 0.2, ... and more? A very crucial hyperparameter (and of course, smart selection of it can not be worse than L2RNN). This study is lacking. In my opinion, one can also introduce weights and sparsity constraints on them (in order to detect the number of "memory" neurons more-or less automatically). Although I feel this paper has a potential, it is not still ready for publication and could be significantly improved.
Decision: Reject