{'fraction_non_alphanumeric': 0.04684932325435478, 'fraction_numerical': 0.021512444351489437, 'mean_word_length': 4.431678026614735, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 8, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Transformation of ReLU-based recurrent neural networks from discrete-time to continuous-time

Recurrent neural networks (RNN) as used in machine learning are commonly formulated in discrete time, i.e. as recursive maps. This brings a lot of advantages for training models on data, e.g. for the purpose of time series prediction or dynamical systems identification, as powerful and efficient inference algorithms exist for discrete time systems and numerical integration of differential equations is not necessary. On the other hand, mathematical analysis of dynamical systems inferred from data is often more convenient and enables additional insights if these are formulated in continuous time, i.e. as systems of ordinary (or partial) differential equations (ODE). Here we show how to perform such a translation from discrete to continuous time for a particular class of ReLU-based RNN. We prove three theorems on the mathematical equivalence between the discrete and continuous time formulations under a variety of conditions, and illustrate how to use our mathematical results on different machine learning and nonlinear dynamical systems examples.


Introduction
Recurrent neural networks (RNN) are popular devices in machine learning and AI for tasks that require processing and prediction of temporal sequences, like machine translation (Sutskever et al., 2014), natural language processing (Kumar et al., 2016;Zaheer et al., 2017), or tracking of moving objects in videos (Milan et al., 2017).
More recently, in the natural sciences, biology and physics in particular, RNNs were also introduced as powerful tools for approximating the unknown nonlinear dynamical system (DS) that produced a set of empirically observed time series, i.e. for identifying the data-generating nonlinear DS in a completely data-driven, bottom-up way (Durstewitz, 2017;Koppe et al., 2019;Razaghi & Paninski, 2019;Vlachas et al., 2018;Zhao & Park, 2017). Theoretically, it has been proven that (continuous) RNN can approximate the flow field of any other nonlinear DS to arbitrary precision on compact sets of the real space under some mild conditions (Funahashi & Nakamura, 1993;Hanson & Raginsky, 2020;Kimura & Nakano, 1998;Trischler & D'Eleuterio, 2016).
RNNs, in the form most widely used in machine learning, constitute discrete-time DS defined by a recursive transition rule (difference equation) which maps the network's activation states among consecutive time steps, z t = F θ (z t−1 , s t ), where θ are parameters of the system and {s t } is a sequence of external inputs. This formulation is highly advantageous for training RNN on observed data sequences since efficient variational inference and Expectation-Maximization algorithms, which do not require numerical integration of nonlinear ODE, exist for discrete-time systems (Durstewitz, 2017;Koppe et al., 2019;Razaghi & Paninski, 2019;Zhao & Park, 2017). In many scientific contexts, like neuroscience (Koch & Segev, 2003), ODE systems are often highly nonlinear and stiff and thus would require more involved implicit numerical integration schemes to achieve accurate solutions (Koch & Segev, 2003;Ozaki, 2012;Press et al., 2007). Furthermore, empirical data are always sampled at discrete time points, such that in most cases this assumption may not be too limiting for the purpose of model inference.
On the other hand, natural systems evolve in continuous time, and hence most mathematical theories in physics and biology are formulated in continuous time. Thus, for the purpose of theoretical analysis of models inferred from experimental data a continuous-time formulation of the system with biologically or physically meaningful time constants would be preferable. Moreover, a continuous time ODE system enables to analyze properties of the DS under study that are much more difficult or impossible to arXiv:2007.00321v1 [math.DS] 1 Jul 2020 assess in discrete time. For instance, a continuous time DS comes with a flow field that enables to visualize the systems dynamics more easily (e.g. Fig.3), and it enables to smoothly interpolate between observed data points and thus to assess the solution at arbitrary time points. This is of advantage in particular when observations come at irregular event times (Chen et al., 2018). Also, separating DS into slow and fast subsystems (separation of time scales) is a powerful analysis tool (Durstewitz & Gabriel, 2007;Rossetto et al., 1998) that is not readily available for discrete-time systems. More generally, the fact that an ODE system is usually smooth almost everywhere eases the mathematical study of many phenomena, like those of periodic and non-periodic solutions, stable and unstable manifolds of fixed points, bifurcations, or stability of solutions more generally (Absil, 2006). In fact, finding explicit solutions is often impossible even for 1-dimensional recursive maps, such that one often has to revert to graphical methods like cobwebs (Absil, 2006;Haschke, 2004).
It would therefore be highly desirable to find ways in which discrete time RNN models could be embedded into continuous time systems without changing their phase space. Such an embedding requires that for a given discrete-time Haschke, 2004). In general, finding such an embedding is not possible for nonlinear discrete-time systems, while the reverse problem, although not trivial (Ozaki, 2012), is much easier and there are different ways for obtaining a discrete time system from a continuous one 1 . The present paper will address this issue using a specific formulation of RNNs, namely piecewise-linear RNNs (PLRNN) that employ rectified-linear units (ReLU) as their activation function. PLRNN are universal in terms of their dynamical repertoire (Koiran et al., 1994;Siegelmann & Sontag, 1995;Lu et al., 2017), can extract long-term dependencies in sequential data just like LSTMs can (Schmidt et al., 2020), and in particular have been used previously to infer nonlinear DS from time series data (Durstewitz, 2017;Koppe et al., 2019). We show that for this specific class of RNN models mathematically equivalent ODE systems, in the sense defined above, can be derived under almost all conditions. We will exemplify these results on a couple of machine learning and DS model systems, including an ODE solution to the well-known addition problem (Hochreiter & Schmidhuber, 1997), limit cycle and chaotic dynamics, and on a PLRNN inferred from empirical time series (human functional magnetic resonance imaging [fMRI] 1 The most important of such methods is the Poincaré map, where a continuous time will be reduced to a discrete time system by successive intersections of the flow of the continuous system with the Poincaré section, such that most dynamical properties of the original system will be preserved. data; (Koppe et al., 2019)).

Related work
In some recent papers (de Brouwer et al., 2019;Jordan et al., 2019) continuous-time ODE approximations of discrete RNN were sought based on inverting the forward Euler rule for numerically solving continuous ODE systems. A related idea is that of Neural ODE (Chen et al., 2018), where the flow is given by a (deep) neural network (cf. (Pearlmutter, 1989)) to yield a system continuous across space (layers) or time (see also (Abarbanel et al., 2018) for closely related ideas). In (Chang et al., 2019) the reverse approach is taken for obtaining a discrete time formulation that preserves certain properties of the ODE system. None of this work, however, explicitly considered the problem of finding an equivalent continuous-time description of a discrete-time RNN. In fact, apart from the fact that a simple forward Euler rule is known to be rather inaccurate and unstable for integrating stiff ODE systems (Press et al., 2007), a naïve 'inversion' will generally not result in a mathematically equivalent ODE system in the sense defined further above, i.e. the resulting system usually will not have the same phase space and temporal behavior. Consequently, much of the previous work was less aimed at finding a mapping between discrete and continuous time networks, but rather to formulate the problem of neural network training in continuous time and/or space to begin with to exploit advantages of an ODE formulation in one way or the other.
Rather, a true translation of a discrete into a continuous RNN may be achieved by taking the continuous time limit , which, however, can be highly nontrivial or impossible for nonlinear systems (Ozaki, 2012). (Ozaki, 2012), while focusing on the continuous-to-discrete case, also briefly discusses some ideas on this reverse direction of a discrete-to-continuous mapping for nonlinear DS. Only approximations are considered, however, that may work well only for certain cases (RNN, in particular, were not discussed), while here we seek exact equivalence according to the definition further above. In (Nykamp, 2019) a discrete logistic equation is transformed to a continuous-time system by taking the temporal limit. It is, however, not possible to apply such a method for all discrete DS, and even if a transformation is possible, the discrete DS may have different dimensions than the continuous time equivalent; for instance, chaotic behavior is possible in a 1-dimensional recursive map but requires 3 dimensions in an ODE system (Strogatz, 2015). Specifically for linear systems, the necessary and sufficient conditions for embedding a discrete-time homogeneous linear system in a continuous-time system have been studied in (Reitmann, 1996). Here, we are interested more generally in discrete-time non-homogeneous systems which are piecewise linear, i.e. piecewise linear (ReLU-based) recurrent neural networks (PLRNN). We will show how to embed a PLRNN into an equivalent continuous-time ODE system, and how the dynamics of the PLRNN is directly connected to the dynamics of the corresponding ODE system. To the best of our knowledge this is the first time a method is introduced for converting discrete into continuous time RNNs in a mathematically exact way, i.e., such that the systems are mathematically and dynamically equivalent without using any approximations or numerical techniques.

Preliminaries
In the following we will collect some results which we will need for our derivations. Theorem 1. Consider the non-homogeneous systeṁ (1) Then φ(t) = e At (with φ(0) = I) is the fundamental matrix solution for linear systemẋ = Ax, and the solution of system (1) has the form (2) Proof. See (Perko, 1991).
Proof. The Taylor expansion of matrix B has the form: and such that B is invertible iff it does not have any zero eigenvalue. So, by (4)-(5), B is invertible iff λT / ∈ 2iπZ * .
Logarithm of real matrices. For a complex matrix a logarithm (not necessarily unique) will exist iff it is invertible (Higham, 2008). Real matrices do not always have a real logarithm. However, the following theorem guarantees the existence of a real logarithm for a real matrix.
Theorem 2. A real matrix A ∈ R n×n has a real logarithm if and only if (I) A is invertible, and (II) every k × k Jordan block associated with a negative eigenvalue occurs an even number of times in the Jordan form of A.
Corollary 2.1. Due to theorem 2 a real nonsingular 2 × 2 matrix A with two eigenvalues λ 1 and λ 2 will have a real logarithm in the following cases: (1) λ 1 and λ 2 are complex conjugate, i.e. λ 1,2 = a ± ib, b = 0. In this case A and log(A) have the Jordan where e α±iβ = a ± ib.
Remark 2.1. Suppose that A is a real nonsingular 2 × 2 matrix which has two equal and negative eigenvalues λ 1 = λ 2 < 0. If A has the Jordan form λ 1 1 0 λ 2 , then it will not have a real logarithm.
where φ(Z t ) = max(Z t , 0) is the element-wise rectified linear unit (ReLU) transfer function, Z t = (z 1t , · · · , z M t ) T ∈ R M denotes the neural state vector at time t = 1 · · · T , the diagonal entries of A = diag(a 11 , · · · , a M M ) ∈ R M ×M represent (linear) auto-regression weights, W ∈ R M ×M is a matrix of connection weights (sometimes assumed to be off-diagonal, i.e. with diagonal elements equal to zero, e.g. (Koppe et al., 2019)), and h is a bias term (Durstewitz, 2017;Koppe et al., 2019).
The discrete-time PLRNN (6) can be represented in the form where , depending on the sign of the components of Z t . That is, the phase space of system (7) is separated such that in each sub-region the dynamics are governed by a different linear map (cf. Fig. S1), i.e.
All sub-regions S Ω i corresponding to (9) together with all switching boundaries Σ ij =S Ω i ∩S Ω j between every pair of successive sub-regions S Ω i and S Ω j , with i, j ∈ {1, 2, · · · , 2 M }, are formally defined in Suppl. sect. 8. Note that map (7) is continuous, but has many discontinuities in the Jacobian across the switching boundaries Σ ij (for more details please see Suppl. sect. 8). It is easy to see that for every pair of matrices D Ω k which differ only in one diagonal entry, their corresponding matrices W Ω k will differ in only one column.

Transformation from discrete to continuous-time
As noted in the Introduction, measurements of physical or biological systems are always carried out at discrete times separated by finite time steps ∆t (often with a constant sampling rate), and efficient algorithms are available for inferring discrete RNN models from such data (Durstewitz, 2017;Koppe et al., 2019;Razaghi & Paninski, 2019;Zhao & Park, 2017). However, the state of natural dynamical systems is usually better described in continuous time, andfurthermore -continuous-time formulations enjoy a number of key advantages (Chen et al., 2018;Haschke, 2004): First, since for ODE systems trajectories are continuous curves rather than collections of single points, some dynamical properties can be determined more easily. For instance, in discrete-time systems it is not possible to distinguish quasiperiodic from periodic orbits with a large period. Likewise, for ODE systems we have a continuous phase portrait (almost) everywhere, and solutions are defined for any arbitrary time point. Finally, some types of analysis are much easier to do in continuous rather than discrete time. For example, performing a change of variables to compute probability distributions within normalizing flows can be more convenient in continuous-rather than discrete-time systems (Chen et al., 2018). 2 In the following we will state a set of theorems which show how to convert a discrete into a continuous PLRNN by transforming discrete-time system (9) on every sub-region into an equivalent continuous-time system. In this way, a system of piecewise ordinary differential equations will be assigned to the discrete-time system (9) on R M . Here we will just state our major results, while all details of the proofs will be given in Suppl. sect. 8. Note that the following theorems settle the problem for one time step ∆t taken by the PLRNN, from which, however, results for more time steps immediately follow. with and time step ∆t. Suppose that W Ω k is invertible and has no eigenvalue equal to one, i.e. P W Ω k (1) = 0, where P W Ω k denotes the characteristic polynomials of W Ω k .
(1) There exists a continuous-time systeṁ which is equivalent to (10) on [t 0 , t 0 +∆t] in the sense that Moreover, in this caseW Ω k is also an invertible matrix and has no eigenvalue equal to one, i.e. PW Ω k (1) = 0. Also, Furthermore, if for W Ω k each of its Jordan blocks associated with a negative eigenvalue occurs an even number of times, thenW Ω k will be a real matrix.
(2) If W Ω k is both invertible and diagonalizable, thenW Ω k will be invertible and diagonalizable too.
Corollary 3.1. The results of theorem 3 are also true if W Ω k is a positive-definite matrix with P W Ω k (1) = 0.
Proof. Let W Ω k be a positive-definite matrix. Then its determinant is positive, which implies that it is invertible, thus satisfying the conditions of theorem 3. Note that if W Ω k is also Hermitian (or symmetric for real matrices), all eigenvalues of W Ω k are real and it is diagonalizable as well (Bhatia, 2007).
In theorem 3 it is assumed that W Ω k has no eigenvalue equal to one. But there are some PLRNNs with interesting computational properties in the form of system (10), for which W Ω k has at least one eigenvalue equal to one (see Example 2). Hence, more generally we are also interested in converting such neural networks from discrete-to continuous-time. The next two theorems are stated to address this problem.
Theorem 4. Consider system (10) and assume that W Ω k is invertible, diagonalizable, and has at least one eigenvalue equal to 1, i.e. P W Ω k (1) = det(I − W Ω k ) = 0. Then, there exists an equivalent, in the sense defined in equation (13), continuous-time system (12) for (10) on [t 0 , t 0 + ∆t] such thatW Ω k is diagonalizable, but not invertible, and where n represents the number of eigenvalues of W Ω k which are equal to 1 (or the number of eigenvalues ofW Ω k equal to zero). Also, if each Jordan block of W Ω k associated with a negative eigenvalue occurs an even number of times, theñ W Ω k will be real.
W Ω k in theorem 4 must be diagonalizable, but for some computationally interesting PLRNNs W Ω k is not diagonalizable.
The following theorem is stated and proved to address this issue.
Theorem 5. Let W Ω k in system (10) be invertible and Then, there exists an equivalent, in the sense of equation (13), continuous-time system (12) for (10) Further suppose that each Jordan block of W Ω k associated with a negative eigenvalue occurs an even number of times, thenW Ω k is a real matrix.

Application examples
In the following we will illustrate how to use our mathematical results for four specific PLRNN systems of relevance in dynamical systems theory and machine learning. These include examples for a nonlinear oscillator (limit cycle), the 'addition problem' introduced by (Hochreiter & Schmidhuber, 1997) to probe long short-term-memory capacities of RNN, an example of a chaotic system (Lorenz attractor), and a PLRNN inferred from empirical (human fMRI) time series. Matlab code for all these examples is available at github.com/DurstewitzLab/contPLRNN. , provided online at github.com/DurstewitzLab). The Jacobian matrix of this system is always invertible and has no eigenvalue equal to one in any of the sub-regions S Ω i . Hence we can use theorem 3 to convert this system from discrete-to continuous-time. Fig. 1A illustrates time graphs overlaid for the discrete (blue circles) and continuous (red curves) PLRNN, while Fig. 1B depicts a 2d section of the system's continuous phase space with corresponding flow field. Note that there is a perfect agreement between the discrete and continuous solutions for the set of times at which discrete-PLRNN outputs are defined, while at the same time the continuous PLRNN smoothly interpolates between the discrete-time values. Also note that this agreement continues across different subregions S Ω k induced by the ReLU function.
As an example of a specific DS analysis that is much easier in the continuous than in the discrete time system we consider a special type of bifurcation (i.e., a point in the system's parameter space where the dynamic abruptly changes), the so-called grazing bifurcation of periodic orbits. It occurs in piecewise smooth continuous-time systems when a periodic orbit tangentially intersects ('grazes') with a switching boundary. A related bifurcation, the so-called bordercollision bifurcation, also occurs in discrete-time systems when a k-cycle collides with one border. However, in the discrete case, finding the specific bifurcation point can be very challenging, especially in high dimensions, since it amounts to solving highly nested nonlinear equations of the general form F k (Z, b)−Z = 0 for large k (where F k is the k-times iterated map and b a bifurcation parameter), and determining among the solutions that particular point that agrees with the conditions of the bifurcation. In the continuous case, in contrast, one can relatively straightforwardly solve the implied system of equations (see Suppl. section 8 (subsection 8.6) for more details), and hence converting the discrete into the continuous PLRNN offers a big advantage. An example for a grazing bifurcation in the continuous-time PLRNN emulation of the van-der-Pol system is shown in Fig.  2. Bifurcation phenomena like these are of great practical importance and may also have fundamental implications for training RNN (Doya, 1992), since they imply a sudden switch in the temporal structure of the system's behavior as the bifurcation point is crossed.
A simple discrete-time 2-unit PLRNN which (approximately) solves this task is the one with parameters Applying definition (8), W Ω k := A + W D Ω k , we have Hence, every W Ω k is invertible and diagonalizable, but has one eigenvalue equal to one, satisfying the conditions of theorem 4. Translating this system into continuous time using theorem 4, Fig. 3A displays the two system variables in both continuous and discrete time. Note that while z 2 directly responds to the random inputs, its activity will be integrated (summed) by z 1 whenever the second (indicator) inputs are on, i.e. s 2 = 1, which is the case here for the two temporal intervals [100, 105] and [400, 405], as can be seen by z 2 crossing the black dashed line (i.e., when z 2 > 0). Visualizing the continuous system's phase space gives some insight into how the RNN solved the addition problem ( Fig. 3B): The z 1 -line forms a line attractor for z 2 = −0.995/(1 − 0.01), with the flow converging toward this line from all directions, and zero flow right on this line, thus yielding a dimension in state space along which arbitrary values could be stored (Durstewitz, 2003;Seung, 1996). A series of sufficiently strong inputs to the RNN (i.e., whenever s 2 = 1) will push the system away from its current to a new position on the line attractor where it will remain until the second supra-threshold series of inputs arrives, in this manner integrating the s 1 inputs accompanied by 1's in s 2 (see also (Schmidt et al., 2020)). In this specific example, the first series of s 1 inputs sums up to ≈ 2.7 (as marked by the left green circle on the z 1 -line at z 2 ≈ −1) and the second to ≈ 3.4, and the PLRNN correctly reports the total sum of inputs (right green circle) in its final position on the line attractor.
Example 3. As an example for a system with chaotic dynamics we chose a PLRNN emulation (M = 10) of the 3d Lorenz system, i.e. a PLRNN trained to reproduce the dynamics of the Lorenz equations within the chaotic regime (taken from (Koppe et al., 2019)). In this case, we could apply theorem 3 to accomplish the transformation to continuous time, as all matrices W Ω k as defined in (8) were invertible with no eigenvalue equal to 1. Fig. 4 confirms that the continuous-time PLRNN agrees with the discrete-time solution also in this case of chaotic behavior.
Example 4. As a final example we applied the continuoustime transform developed here to a PLRNN inferred from empirical time series, namely human fMRI data. A discrete PLRNN with M = 10 latent states was used for this purpose that had been trained on multivariate (N = 20) time series of Blood Oxygenation Level Dependent (BOLD) signals recorded from human subjects performing a cognitive Figure 1. Transformation of a discrete PLRNN emulating the nonlinear van-der-Pol oscillator into a continuous-time ODE system. A) Time graphs for two of the system's 10 variables (unit activations). A zoom-in is provided for z2 to better highlight how the continuous solution interpolates between the discrete time points. Blue = discrete PLRNN, red = continuous PLRNN. B) Continuous 2d-subspace of the 10-dimensional state space corresponding to the two variables shown in A, with flow fields (black arrows) and the system's trajectory on the limit cycle (red = continuous, blue circles = discrete); note that since this is only a 2d section of a 10-variable system, convergence to the limit cycle cannot be fully assessed from the (z1, z2) vector field. (2) bif , where the system's trajectory tangentially touches the border Σ (center panel) and the behavior changes from simple-periodic (left panel) to a period-2 limit cycle (right panel). task while lying in a fMRI scanner. Details on the experimental procedure and task and on PLRNN training can be found in (Koppe et al., 2019) (briefly, an Expectation-Maximization algorithm was used for PLRNN training to maximize the evidence-lower-bound [ELBO] of the data log-likelihood). Theorem 3 could be applied in this case to translate the discrete-time system to continuous time, as all system's matrices met the conditions of this theorem. As Fig.  5A demonstrates, the continuous-time PLRNN (red curves) smoothly interpolates between the predictions produced by the discrete PLRNN (blue circles), and thus smoothly extrapolates among the observed data points when started from an initial condition inferred from the experimental data (Fig.  5B, green circles).

Conclusions
The aim of the present article was to show how to convert discrete-time into mathematically equivalent continuoustime RNN. As pointed out in the Introduction, sect. 1, and sect. 4.2, while RNN are usually trained in discrete time, a continuous-time description enjoys multiple advantages when it comes to analyzing the inferred systems and linking them to scientific theories. Here we examined such a transformation for a particular class of RNN based on rectified-linear unit (ReLU) activation functions. ReLU transfer functions are by now the most common choice in the deep learning community due to their piecewise constant gradients (Goodfellow et al., 2016;Lin & Jegelka, 2018;Montfar et al., 2014), easing gradient-descent based algorithms. ReLU-based RNN are computationally and  . Note that since this system is chaotic, time graphs will never precisely overlap since any small numerical difference will lead to exponential divergence of trajectories (here, the recurrence matrices were only close to real for the continuous system and imaginary parts were set to 0). B) 3d-projection of the 10-dimensional PLRNN state space corresponding to the three observation variables, exposing the butterfly-wing-type structure of the chaotic Lorenz attractor in discrete (blue circles) and continuous (red lines) time.
dynamically universal (Kimura & Nakano, 1998;Koiran et al., 1994;Lu et al., 2017), can outperform LSTMs on long short-term memory problems (Le et al., 2015;Schmidt et al., 2020), and can be efficiently inferred from time series data for approximating the underlying (nonlinear) dynamical system (Koppe et al., 2019). Hence our focus on the class of piecewise-linear RNN is not very restrictive, but instead encompasses a whole powerful family of RNN architectures and algorithms. Specifically, we proved that such a conversion from discrete to continuous time is possible under a variety of conditions. These include situations where one or more of the eigenvalues of the systems Jacobian are equal to 1, as required for long short-term maintenance (Example 2), or when we have cycles in the discrete case, leading to complex eigenvalues in the continuous case, as we illustrated in our specific applicational examples. Future work may address the problem of discrete-to-continuous-time conversion more generally, for arbitrary nonlinear activation functions, or may attempt to find useful solutions for those cases where a direct translation is not easily possible, e.g. when the recurrence matrix W Ω k is not invertible, when the ODE system corresponding to the recursive map would need to have different dimensionality (as may occur, e.g., for low-dimensional chaotic maps like the the logistic map), or when the time step ∆t is not constant. Another valuable extension, especially in the context of latent variable models, would be to stochastic DS.

Acknowledgements
This work was funded by the German Science Foundation (DFG) through individual grant Du 354/10-1 to DD, and via the Excellence Cluster 'Structures' at Heidelberg University (EXC-2181 390900948). We thank Dr. Georgia Koppe for kindly lending us the fMRI data and corresponding PLRNN parameters used for the empirical example in Fig.5.

Discontinuity boundaries
Consider map (7) and two sub-regions S Ω i and S Ω j (i, j ∈ {1, 2, · · · , 2 M }) as defined in Section 4 (subsection 4.2). Suppose that subindices i − 1 and j − 1 ofŜ i−1 andŜ j−1 are associated with i − 1 = i 1 i 2 · · · i M and j − 1 = j 1 j 2 · · · j M . Then S Ω i and S Ω j are two successive subregions with the switching boundary Moreover, Σ ij is a closed set (Σ ij = Σ ij ) and Σ ij = and Furthermore, it can be proven that 8.3. Proof of theorem 3 (1) Without loss of generality let t 0 = 0. Assume that there exists an equivalent continuous-time system for (10) Transformation of ReLU-based recurrent neural networks from discrete-time to continuous-time Figure S1. Example of subregions S Ω k and related matrices D Ω k for M = 2.
on [0, ∆t], in the form of equation (12). By equivalency in the sense of equation (13), we must have According to theorem 1, the solution of system (10) on [0, ∆t] is IfW Ω k is invertible, then and thus Furthermore, since we have Putting conditions (31) in yields Equation (38) has to hold for all Z 0 including Z 0 = 0. Hence, it is deduced that According to (39), matrix W Ω k should be invertible and cannot have any zero eigenvalue. Also, sinceW Ω k is invertible, it does not have any zero eigenvalue, which implies that W Ω k has no eigenvalue equal to one. Then, P W Ω k (1) = 0, which means I − W Ω k is invertible and (39) becomes equivalent to (14). Now, consideringW Ω k andh as in (14), we can obtain the desired equivalent continuous-time system (12) for (10) on [0, ∆t]. It is just required to prove that every fixed point Z * of map (10) is also an equilibrium point of system (12), and (14) is a solution of (38) for all Z * . For this purpose, let Z * be a fixed point of (10), then Z * must be an equilibrium of (12), i.e.
From (40) and (41) it is concluded that which shows that (39) or, equivalently, (14) is a solution of (38) for all Z * satisfying both relations (40) and (41). Finally, let each Jordan block of W Ω k associated with a negative eigenvalue occur an even number of times. Then, by theorem (2), the logarithm of real matrix W Ω k , i.e. the matrixW Ω k defined in (14), will be real.
(2) Let W Ω k be diagonalizable, then where E k = diag(λ k 1 , λ k 2 , · · · , λ k M ) and V is the matrix of eigenvectors of W Ω k . Since W Ω k is also invertible, by (14) W such that which completes the proof.

Proof of theorem 4
Again we prove the theorem for t 0 = 0 without loss of generality. Suppose that there is the equivalent continuous-time system (12) for (10) with non-invertible and diagonalizable matrixW Ω k . Similar to the proof of the previous theorem, relations (31) and (32) must hold for (10) and (12). On the other hand, non-invertibility and diagonalizability ofW Ω k demand that it has at least one eigenvalue equal to zero and where O n×n is a zero matrix corresponding to zero eigenvalues (n denotes the number of zero eigenvalues) and C is an invertible matrix corresponding to nonzero eigenvalues ofW Ω k . Therefore, for relation (33) we obtain In this case, relation (37) becomes Inserting conditions (31) into (47) gives Denoting and considering equality (48) for all Z 0 , particularly for Z 0 = 0, yields Since we can simplify h in (50) and rewrite it as or equivalently In addition, we can writẽ Thereforẽ which is equivalent to (15).
Finally, from W Ω k = eW Ω k ∆t it is deduced that W Ω k is invertible. It is only necessary to prove that for every point Z * satisfying both equations (40) and (41), i.e. equations relation (15) which demonstrates that (56) meets solution (52).
If every Jordan block of W Ω k associated with a negative eigenvalue occurs an even number of times, then theorem (2) guarantees thatW Ω k will be real. Also, similar to the proof of theorem 3, it is easy to see thatW Ω k will be diagonalizable when W Ω k has no negative real eigenvalues.

Proof of theorem 5
Let t 0 = 0 without loss of generality and assume there exists the equivalent continuous-time system (12) for (10), for which matrix W Ω k is non-invertible. Then, relations (31) and (32) must hold for (10) and (12), analogously to the proofs of the previous theorems. Also, by similar reasoning we have ζ(∆t) = eW Ω k ∆t ζ(0) + eW Ω k ∆t ∆t 0 e −W Ω k τ dτ h .
Then, by proposition 1, ∆t 0 e −W Ω k τ dτ is invertible and so   which is equal to equation (17). The last point which still has to be proven is that equation (56) meets solution (17) or, identically, (59), for every Z * . SinceW Ω k is non-invertible, it can be written in the following Jordan form: where B is a strictly upper triangular matrix and C is an invertible matrix. Then which completes the proof.
Finally, due to theorem (2),W Ω k will be real, provided that each Jordan block of W Ω k related to a negative eigenvalue occurs an even number of times.
Remark. In theorem 5, by (64) On the other hand, since C is invertible, det I − e −C∆t = 0. Therefore, ∆t 0 e −W Ω k τ dτ is invertible if and only if ∆t 0 e −Bτ dτ is invertible. Thus, for invertibility of ∆t 0 e −W Ω k τ dτ , it is required that relation (60) holds only for any pair of eigenvalues of B.