Gradient-free training of recurrent neural networks

## Abstract

Abstract Recurrent neural networks are a successful neural architecture for many time-dependent problems, including time series analysis, forecasting, and modeling of dynamical systems.
Training such networks with
backpropagation through time is a notoriously difficult problem because their loss gradients
tend to explode or vanish.
In this contribution, we introduce a computational approach to construct all weights and biases of a recurrent neural network without using gradient-based methods.
The approach is based on a combination of random feature networks and Koopman operator theory for dynamical systems.
The hidden parameters of a single recurrent block are sampled at random, while the outer weights are constructed using extended dynamic mode decomposition.
This approach alleviates all problems with backpropagation commonly related to recurrent networks.
The connection to Koopman operator theory also allows us to start using results in this area to analyze recurrent neural networks.
In computational experiments on time series, forecasting for chaotic dynamical systems, and control problems, as well as on weather data,
we observe that the training time and forecasting accuracy of the recurrent neural networks we construct are improved when compared to commonly used gradient-based methods.

## 1 Introduction

Recurrent neural networks (RNNs) are notoriously difficult to train because their loss gradients backpropagated in time tend to saturate or diverge during training, commonly referred to as the Exploding and Vanishing Gradient Problem (EVGP) (Pascanu et al., ; Schmidt et al., ).
In Equation 1, we denote a general non-linear system governed by maps $F$ and $g$, controlled by external inputs ${\bm{x}}_{t}$.
A typical RNN modeling approach would directly train all parameters of the networks ${\mathcal{F}},{\mathcal{G}}$ as well as all matrices $K,B,C,V$ on the right of Equation 1.

In our approach, we first construct a linear SSM based on Koopman operator theory (cf. Korda and Mezić ()). The matrix $K$ maps states from ${\bm{z}}_{t-1}$ to ${\bm{z}}_{t}$, with the input ${\bm{x}}_{t}$ affecting this dynamic through the matrix $B$. The observations are related to the states through the linear map $V$.

$$ $\underbrace{\begin{array}[]{rcl}{\bm{h}}_{t}&=&F({\bm{h}}_{t-1},{\bm{x}}_{t}) \\ {\bm{y}}_{t}&=&g({\bm{h}}_{t})\end{array}}_{\text{original non-linear system}} \leftrightarrow\underbrace{\begin{array}[]{rcl}{\bm{z}}_{t}&=&K{\bm{z}}_{t-1}+ B{\bm{x}}_{t}\\ {\bm{y}}_{t}&=&V{\bm{z}}_{t}\end{array}}_{\text{linear (Koopman) state-space model}}\leftrightarrow\underbrace{\begin{array}[]{rcl}{\bm{z}}_{t}&=&K{ \mathcal{F}}({\bm{h}}_{t-1})+B{\mathcal{G}}({\bm{x}}_{t})\\ {\bm{h}}_{t}&=&C{\bm{z}}_{t}\\ {\bm{y}}_{t}&=&V{\bm{z}}_{t}\end{array}}_{\text{non-linear RNN}}.$ (1) $$

We propose to use a non-linear neural network to map the low-dimensional state ${\bm{h}}$ into a higher dimension, through ${\mathcal{F}}({\bm{h}}):=\sigma(W{\bm{h}}+{\bm{b}})$, and, similarly, to map the input ${\bm{x}}$ to a high-dimensional input ${\mathcal{G}}({\bm{x}})=\sigma_{x}(W_{x}{\bm{x}}+{\bm{b}}_{x})$. This turns the linear SSM into a non-linear recurrent neural network. Generally, such a non-linear map is also the core idea for many Koopman operator-based models: start with the state ${\bm{h}}_{t-1}$ of a non-linear system, use a non-linear map to a higher-dimensional state ${\bm{z}}_{t-1}$, and then move forward in time in this high-dimensional space using a linear map $K$. This map is then an approximation of the linear Koopman operator of the original, non-linear system.

Figure: Figure 1: Illustration of the components of one recurrent block we construct in the paper. The state ${\bm{z}}_{t-1}$ enters on the left, and is processed through matrix $C$ and the neural network ${\mathcal{F}}=\sigma(W\cdot+{\bm{b}})$. We then advance in time to ${\bm{z}}_{t}$, using the Koopman matrix $K$ and the processed control inputs.
Refer to caption: x1.png

Crucially, in this contribution, we sample the weights $W$ and biases ${\bm{b}}$ from a specific, data-dependent probability distribution. The remaining free parameters in the matrices $K,B,C,V$ can then be found by solving a sequence of linear problems (cf. Figure 1).
This approach alleviates all issues related to the backpropagation of gradients and connects the idea of recurrent architectures to the Koopman operator of the underlying dynamical system. The latter allows us to prove convergence results about the approximation quality w.r.t. the number of neurons in ${\mathcal{F}}$.
We demonstrate the performance of this approach on several challenging examples with synthetic and real data.

## 2 Related work

Exploding and vanishing gradients. For all major types of RNNs, including LSTMs and GRUs, the dynamics and loss gradients of RNNs are closely linked. If the RNN dynamics converge to a stable fixed point or cycle, loss gradients will remain bounded, but they may vanish (Mikhaeil et al., ). Yet, established remedies (Hochreiter and Schmidhuber, ; Schmidt et al., ) can be used to effectively prevent their gradients from vanishing. However, in chaotic dynamics, gradients invariably explode, posing a challenge that cannot be mitigated through RNN architectural adjustments, regularization, or constraints; instead, it necessitates addressing the problem during the training process (Mikhaeil et al., ).
Bifurcations may also contribute to sudden jumps in loss observed during RNN training, potentially hindering the training process severely (Doya et al., ; Eisenmann et al., ). Eisenmann et al. () have demonstrated that specific bifurcations in ReLU-based RNNs are always associated with EVGP during training. Therefore, to harness the full potential of RNNs, the training algorithm needs careful design to tackle challenges posed by bifurcations and the possible emergence of EVGP.
Curse of memory. The existence of long-term memory adversely affects the learning process of RNNs (Bengio et al., ; Hochreiter et al., ; Li et al., ). This negative impact is captured by the concept of the “curse of memory”, which states that when long-term memory is present in the data, approximating relationships demands an exponentially large number of neurons, resulting in a significant slowdown in learning dynamics. Specifically, when the target relationship includes long-term memory, both the approximation and optimization of the learning process become very challenging (Li et al., ).
Loss function for chaotic and multistable dynamics reconstruction. To effectively train RNNs and evaluate reconstruction, it is crucial to carefully choose a proper loss function. The Mean Squared Error (MSE) is a commonly used loss function for reconstruction tasks.
MSE is derived under the assumption of Gaussian noise, and it may not be the most appropriate choice when dealing with chaotic systems or multistable dynamics, where the underlying noise characteristics may deviate from Gaussian distribution assumptions. It is not suitable as a test loss for chaotic dynamical systems due to their unpredictable behavior and abrupt changes (Wood, ).
In multistable systems, where there are multiple stable states, MSE may struggle to distinguish between these states. The loss function may not adequately penalize deviations between different attractors, leading to a less accurate reconstruction of the system’s multistable behavior.
Despite proposed alternatives, challenges persist, and an optimal loss function for reconstructing chaotic or multistable dynamics is still lacking (Ciampiconi et al., ).

Interpretability deficiency. Deep learning models are commonly regarded as “black boxes”,
and existing methods to comprehend the decision-making processes of RNNs offer restricted explanations or rely on local theories. The lack of theory behind analyzing the training algorithms of RNNs, as well as the training process itself, are outstanding open questions in the field (Redman et al., ) beyond models for linear systems (Datar et al., ).

Randomly choosing the internal network parameters.
The general idea of randomly choosing the internal parameters of neural networks is studied in random feature models (including deep architectures). Barron () and Rahimi and Recht () developed the basic theory, and Gallicchio and Scardapane () provide a review.
Reservoir computing (also called echo-state networks, cf. Jaeger and Haas ()) is potentially the closest idea to what we are proposing here. In a reservoir computer, the internal weights are also randomly sampled, and a recurrent, time-delayed model is constructed to approximate a given dynamical system.
This type of architecture has been used successfully to model chaotic systems (cf. Pathak et al. (); Gauthier et al. (); Bollt ()).
While there are similarities to the recurrent architecture (cf. Lukoševičius and Jaeger ()), the concept, as well as the architecture of a reservoir computer, is often treated separately from classic recurrent neural networks that are trained with backpropagation-in-time. In our work, we directly compute all parameters of classical recurrent neural networks without the time-delay component present in reservoir computers.

Koopman operator theory.
The Koopman operator is an object associated to every dynamical system. It evolves observables of the state of the system in time. This evolution is linear, which is the main reason the operator is employed and studied extensively for modeling dynamical systems (Mezić, , ; Korda and Mezić, ; Gallos et al., ). Many numerical approximation algorithms exist (Schmid, ; Williams et al., ; Li et al., ; Mezic, ; Schmid, ).
The dictionary for the approximation of the Koopman operator has been constructed with neural networks using gradient descent (Li et al., ) and using random features (Salam et al., ).
Reservoir computing has also been related to Koopman operator approximation by Bollt (); Gulina and Mauroy ().
To our knowledge, the relation of the Koopman operator to the weight matrices of recurrent neural networks has not been observed before. This is what we discuss in this work. We also provide a data-dependent probability distribution for the hidden parameters, which is the strongest deviation from the data-agnostic distributions (e.g., normal, uniform) typically used in reservoir computing.

## 3 Mathematical framework

We start by defining a general framework of RNNs and the underlying dynamical system before introducing sampling and its connection to the Koopman operator. Let $\mathcal{X}\subseteq\mathbb{R}^{d_{x}}$ be an input space, $\mathcal{Y}\subseteq\mathbb{R}^{d_{y}}$ an output space, and $\mathcal{H}\subseteq\mathbb{R}^{d_{h}}$ a state space. We assume that these spaces are associated with the measures $\mu_{x}$, $\mu_{y}$, and $\mu_{h}$, respectively. The underlying dynamical system is then defined through the evolution operator $F$, where we may be working in an uncontrolled system ${\bm{h}}_{t}=F({\bm{h}}_{t-1})$ or a controlled system ${\bm{h}}_{t}=F({\bm{h}}_{t-1},{\bm{x}}_{t})$. We also denote the input dataset $X=[{\bm{x}}_{1},{\bm{x}}_{2},\dots,{\bm{x}}_{N}]$ and the dataset of observations $Y=[{\bm{y}}_{1},{\bm{y}}_{2},\dots,{\bm{y}}_{N}]$. In this paper we are interested in RNNs modelling dynamical systems where the state is observable, and we therefore usually assume access to the dataset $H=[{\bm{h}}_{1},\dots,{\bm{h}}_{N}]$ as well as its copy after one time step $H^{\prime}=[{\bm{h}}^{\prime}_{1},\dots,{\bm{h}}^{\prime}_{N}]$, where ${\bm{h}}^{\prime}_{t}=F({\bm{h}}_{t})$, $t\in\left\{1,\dots,N\right\}$.

We denote activation functions as $\sigma\colon\mathbb{R}\rightarrow\mathbb{R}$, where we are mainly working with $\sigma=\tanh$ in this paper, as it is an analytic function and connects to SWIM (Bolager et al., ). Other functions such as ReLU are also a valid choice. The following definition outlines the models we consider.

###### Definition 1 .

Let $W_{h}\in\mathbb{R}^{M\times d_{h}}$, $W_{x}\in\mathbb{R}^{\hat{M}\times d_{x}}$, ${\bm{b}}_{h}\in\mathbb{R}^{M},{\bm{b}}_{x}\in\mathbb{R}^{\hat{M}}$, $C_{h}\in\mathbb{R}^{d_{h}\times M}$, $C_{x}\in\mathbb{R}^{d_{h}\times\hat{M}}$, and $V\in\mathbb{R}^{d_{y}\times d_{h}}$. For time step $t$ and ${\bm{h}}_{0}\in\mathcal{H}$, we define an RNN by

$$ $\displaystyle{\bm{h}}_{t}$ $\displaystyle=\sigma_{hx}(C_{h}\,\sigma(W_{h}\,{\bm{h}}_{t-1}+{\bm{b}}_{h})+C_ {x}\,\sigma(W_{x}\,{\bm{x}}_{t}+{\bm{b}}_{x})+{\bm{b}}_{hx}),$ (2) $\displaystyle{\bm{y}}_{t}$ $\displaystyle=V{\bm{h}}_{t}.$ (3) $$

###### Remark 1 .

For completeness we have added $\sigma_{hx}$ as an arbitrary activation function. We choose to set $\sigma_{hx}$ as the identity function to let us solve for the last linear layer in the procedure described below. Other activation functions such as the logit is possible as well.

The classical way to train this type of RNN is through iterative backpropagation, which suffers the aforementioned issues such as EVGP and high computational complexity. We instead start by sampling the hidden layer parameters to circumvent backpropagation, as explained next.

###### Definition 1 .

###### Remark 1 .

### 3.1 Sampling RNN

Sampled neural networks are neural networks where the parameters of the hidden layers are sampled from some distribution, and the last linear layer is either sampled or, more typically, solved through a linear solver. Following ideas from Galaris et al. () and Bolager et al. (), we sample the weights and biases of the hidden layers of both $F_{\mathcal{H}}$ and $F_{\mathcal{X}}$ by sampling pairs of points from the domain $\mathcal{H}$ and $\mathcal{X}$, and construct the weights and biases from these pairs of points. One then solves a linear (general) regression problem at the end for the last linear layer that maps to the next state. Concretely, let $\mathbb{P}_{\mathcal{H}}$ and $\mathbb{P}_{\mathcal{X}}$ be probability distributions over $\mathcal{H}^{2}$ and $\mathcal{X}^{2}$ respectively. For each neuron in the hidden layer of $F_{\mathcal{H}}$, sample $({\bm{h}}^{(1)},{\bm{h}}^{(2)})\sim\mathbb{P}_{\mathcal{H}}$ and set the weight $w$ and the bias $b$ of said neuron to

$$ $\displaystyle{\bm{w}}=s_{1}\frac{{\bm{h}}^{(2)}-{\bm{h}}^{(1)}}{\lVert{\bm{h}} ^{(2)}-{\bm{h}}^{(1)}\rVert^{2}},\quad b=-\langle{\bm{w}},{\bm{h}}^{(1)} \rangle+s_{2},$ (4) $$

where $\lVert\cdot\rVert$ and $\langle\cdot,\cdot\rangle$ are typically the Euclidean norm and inner product, and $s_{1},s_{2}\in\mathbb{R}$ are constants. Repeating the same procedure for all neurons in both $F_{\mathcal{H}}$ and $F_{\mathcal{X}}$. As we stick to networks with one hidden layer in this paper, we ignore the multilayer sampling here and direct the reader to Bolager et al. () for the full sample and construction procedure for an arbitrary number of hidden layers.

This sampling technique adapts the weights and biases to the underlying domain and constructs weights with direction along the data (see Appendix C for an example how this can be used to interpret the resulting network). Empirically, this has shown to be an improvement over using data-agnostic distributions such as the standard Gaussian one.

One can choose arbitrary probability distributions as $\mathbb{P}_{\mathcal{H}}$ and $\mathbb{P}_{\mathcal{X}}$, with uniform distribution being a common choice. For the supervised setting, Bolager et al. () also proposed a sampling distribution whose density captures the steepest gradients of the target function. For this paper, we sample with densities $p_{\mathcal{H}}$ and $p_{\mathcal{X}}$ proportional to

$$ $\displaystyle p_{\mathcal{H}}\propto\frac{\lVert F({\bm{h}}^{(2)})-F({\bm{h}}^ {(1)})\rVert}{\lVert{\bm{h}}^{(2)}-{\bm{h}}^{(1)}\rVert},\quad p_{\mathcal{X}} \propto 1,$ $$

respectively. Once the weights and biases are sampled, we must solve a general regression problem

$$ $\displaystyle[C_{h},C_{x}]=\operatorname*{arg\,min}_{\hat{C}_{h},\hat{C}_{x}} \sum_{n=1}^{N}\lVert(\hat{C}_{h}\,\sigma(W_{h}\,{\bm{h}}_{n}+{\bm{b}}_{h})+ \hat{C}_{x}\,\sigma(W_{x}\,{\bm{x}}_{n}+{\bm{b}}_{x})+{\bm{b}}_{hx})-{\bm{h}}^ {\prime}_{n}\rVert^{2}.$ (5) $$

To summarize, we define a sampled RNN as a model that is constructed by sampling weights of the hidden layer of the RNN and subsequently solving the regression problem in Equation 5.

### 3.2 Involving the Koopman operator

We already introduced the network ${\mathcal{F}_{M}\colon\mathbb{R}^{d_{h}}\rightarrow\mathbb{R}^{M}=\sigma(W_{h}
\cdot+{\bm{b}}_{h})}$, where $M$ is the number of neurons in its single hidden layer, and likewise with $\mathcal{G}_{\hat{M}}\colon\mathbb{R}^{d_{x}}\rightarrow\mathbb{R}^{\hat{M}}=
\sigma(W_{x}\cdot+{\bm{b}}_{x})$. We then project down to $\mathbb{R}^{d_{h}}$ by applying $C_{h}$ and $C_{x}$ respectively. To add more structure and interpretability to the matrices $C_{h}$ and $C_{x}$, we will set $C_{h}=CK$ and $C_{x}=CB$, where $K\in\mathbb{R}^{M\times M}$, $B\in\mathbb{R}^{M\times\hat{M}}$, and $C\in\mathbb{R}^{d_{h}\times M}$. We then end up with the function

$$ $\displaystyle{\bm{h}}_{t}=C{\bm{z}}_{t}=C(K\sigma(W_{h}\,{\bm{h}}_{t-1}+{\bm{b }}_{h})+B\sigma(W_{x}\,{\bm{x}}_{t}+{\bm{b}}_{x})).$ (6) $$

Note that if the output ${\bm{y}}$ differs from ${\bm{h}}$, we take advantage of the high dimensionality and rather set ${\bm{y}}_{t}=V{\bm{z}}_{t}$ than first projecting down by $C$. This completes the setting from Equation 1.

The reason for the splitting of matrices $C_{h}$ and $C_{x}$ is the following: the hidden layers $\mathcal{F}_{M}$ and $\mathcal{G}_{\hat{M}}$ map their respective input to a higher dimensional space. In a higher dimensional space, the possibly nonlinear evolution described by $F$ becomes more and more linear (Korda and Mezić, ). This evolution is then captured by the linear transformation of $K$ and $B$ before we map down to the state space through $C$. This also allows us to connect $K$ and $B$ to the Koopman theory applied to the dynamical system. Given a suitable functional space $\mathcal{F}$, the Koopman operator $\mathcal{K}\colon\mathcal{F}\rightarrow\mathbb{R}$ is defined as

$$ $\displaystyle[\mathcal{K}\phi]({\bm{h}})=(\phi\circ F)({\bm{h}}),\quad\phi\in \mathcal{F}.$ $$

The Koopman operator captures the evolution of the dynamical systems in the function space $\mathcal{F}$ and not in the state space itself. In most cases, this makes the operator infinite-dimensional, but in return, a linear operator. For an introduction to the Koopman operator and surrounding theory, see Appendix A. The matrices $K$ and $B$ can then be seen as an approximation of the Koopman operator, and the intuition is that choosing $M$ large enough means we can capture the linear evolution before we map down to the state space.

The matrices $W_{h}$ and $W_{x}$ are found by sampling each row i.i.d. from $\mathbb{P}_{h}$ and $\mathbb{P}_{x}$ respectively.
To estimate $C,K,$ and $B$, we use extended dynamic mode decomposition (EDMD) — a classical method to find a finite-dimensional approximation of the Koopman operator. We give a very brief description of EDMD here and give a more thorough introduction in Section A.1.
In the uncontrolled setting, this method picks out a dictionary $\tilde{\mathcal{F}}_{M}=\{\psi_{1},\dots,\psi_{M}\colon\psi_{i}\colon\mathcal{
H}\rightarrow\mathbb{R}\in\tilde{\mathcal{F}}\}$ and estimates the Koopman operator $\mathcal{K}$ using the data $H,H^{\prime}$ by minimizing

$$ $\displaystyle K=\operatorname*{arg\,min}_{\tilde{K}\in\mathbb{R}^{M\times M}} \sum_{n=1}^{N}\lVert\tilde{\mathcal{F}}_{M}({\bm{h}}^{\prime}_{n})-\tilde{K} \tilde{\mathcal{F}}_{M}({\bm{h}}_{n})\rVert,\quad{\bm{h}}^{\prime}_{n}\in H^{ \prime},{\bm{h}}_{n}\in H,$ $$

where $\tilde{\mathcal{F}}_{M}({\bm{h}})=[\phi_{1}({\bm{h}}),\dots,\phi_{M}({\bm{h}})
]^{\mkern-1.5mu\mathsf{T}}$. Letting $\tilde{\mathcal{F}}_{M}(H)=[\tilde{\mathcal{F}}_{M}({\bm{h}}_{1}),\dots,\tilde
{\mathcal{F}}_{M}({\bm{h}}_{N})]\in\mathbb{R}^{M\times N}$ and $\tilde{\mathcal{F}}_{M}(H^{\prime})=[\tilde{\mathcal{F}}_{M}({\bm{h}}^{\prime}
_{1}),\dots,\tilde{\mathcal{F}}_{M}({\bm{h}}^{\prime}_{N})]\in\mathbb{R}^{M
\times N}$, the approximation can then be written as

$$ $\displaystyle K=\tilde{\mathcal{F}}_{M}(H^{\prime})\,\tilde{\mathcal{F}}_{M}(H )^{+},$ (7) $$

where ^+ is the matrix pseudoinverse. Similarly, in the controlled setting, the approximation of $\mathcal{K}$ separated into matrices $K$ and $B$,

$$ $\displaystyle[K,B]=\tilde{\mathcal{F}}_{M}(H^{\prime})\,[\tilde{\mathcal{F}}_{ M}(H),\tilde{\mathcal{G}}_{\hat{M}}(X)]^{+},$ $$

where $\tilde{\mathcal{G}}_{\hat{M}}$ is the second dictionary mapping from $\mathbb{R}^{d_{x}}$ to $\mathbb{R}^{\hat{M}}$. For more on the Koopman operator in the controlled setting, as well as EDMD, see Section A.2. Regardless of whether uncontrolled or controlled, the mapping $C$ projects down to the state space from the high dimensional dictionary space and is approximated by minimizing $\lVert H-C\tilde{\mathcal{F}}_{M}(H)\rVert$, hence

$$ $\displaystyle C=H\tilde{\mathcal{F}}_{M}(H)^{+}.$ $$

We connect the EDMD algorithm to our recurrent network in Equation 6 by setting the functions $\tilde{\mathcal{F}}$ and $\tilde{\mathcal{G}}$ as hidden layers, i.e., $\tilde{\mathcal{F}}_{M}({\bm{h}})=\mathcal{F}_{M}({\bm{h}})=\sigma(W_{h}\,{\bm
{h}}+{\bm{b}}_{h})$ and $\tilde{\mathcal{G}}_{\hat{M}}({\bm{x}})=\mathcal{G}_{\hat{M}}({\bm{x}})=\sigma
(W_{x}\,{\bm{x}}+{\bm{b}}_{x})$. With this we see that approximating $K$ and $B$ in Equation 6 can be seen as a Koopman approximation with the dictionary being a hidden layer with $M$ and $\hat{M}$ neurons respectively.
This connection highlights the benefit of operating in a higher dimensional space. It also allows us to make use of Koopman theory in the next section as the EDMD approximation is known to converge to $\mathcal{K}$. Finally, it is important to notice that the resulting function in Equation 6 consists of two neural networks applied to the previous state and input. The hidden layers are sampled, and the outer matrices are constructed using linear solvers.
Hence, Equation 6 is a sampled
recurrent neural network. The methods above can be summarized by sampling weights in Algorithm 1, constructing the RNN in Algorithm 2, prediction for uncontrolled systems in Algorithm 3, and model predictive control in Algorithm 4.

Figure: Algorithm 1 Sampling weights and bias for a given dataset and probability distribution.

Figure: Algorithm 2 Sampling RNNS for the controlled setting with output ${\bm{y}}$.

Figure: Algorithm 3 Prediction of new trajectory uncontrolled inputs using sampled RNN.

Figure: Algorithm 4 Model predictive control (MPC) of system $F$ using LQR and sampled RNN.

### 3.3 Convergence of sampled RNNs

Under some conditions, the convergence of sampled RNNs for uncontrolled systems can be shown for arbitrary finite horizon predictions. The result shows convergence by estimating $K$ using Equation 7. This differs from the usual existence proofs for parameters of RNNs, and is possible due to the Koopman connection established in the previous section.

With $L^{2}\coloneqq L^{2}(\mathcal{H},\mu_{h})$ being the usual Lebesgue space, we can state the required assumptions.

###### Assumption 1 .

The assumptions on $\mu_{h}$, $\mathcal{F}_{M}$, and $F$ are the following.

- 1.
$\mu_{h}$ is regular and finite for compact subsets.
- 2.
Hidden layer $\mathcal{F}_{M}$ must fulfill $\mu_{h}\{{\bm{h}}\in\mathcal{H}\mid{\bm{c}}^{\mkern-1.5mu\mathsf{T}}\mathcal{F
}_{M}({\bm{h}})=0\}=0$, for all nonzero ${\bm{c}}\in\mathbb{R}^{M}$.
- 3.
The Koopman operator $\mathcal{K}\colon L^{2}\rightarrow L^{2}$ is a bounded operator.

The first two points are not very restrictive and hold for many measures and activation functions (such as $tanh$ activation function and the Lebesgue measure). The third assumption is common when showing convergence in Koopman approximation theory and holds for a broad set of dynamical systems (See Section B.1 for further discussion of all three points). Finally, we also require $\mathcal{H}$ to follow Definition 2 in Section B.1.

We now denote $L^{2}_{K}$ as the space of vector valued functions functions $f=[f_{1},f_{2},\dots,f_{K}]$, where $f_{i}\in L^{2}$ and $\lVert f\rVert=\sum_{k=1}^{K}\lVert f_{k}\rVert_{L^{2}}$. We let $F^{t}({\bm{h}}_{0})={\bm{h}}_{t}$ be the true state after time $t$, and $K_{N}$ be the solution of Equation 7, where $N$ data points have been used to solve the least square problem.

###### Theorem 1 .

Let $f\in L^{2}_{K}$, $H,H^{\prime}$ be the dataset with $N$ data points used in Equation 7, and Assumption 1 holds. For any $\epsilon>0$ and $T\in\mathbb{N}$, there exist an $M\in\mathbb{N}$ and hidden layers $\mathcal{F}_{M}$ and matrices $C$ such that

$$ $\displaystyle\lim_{N\rightarrow\infty}\int_{\mathcal{H}}\lVert CK_{N}^{t} \mathcal{F}_{M}-f\circ F^{t}\rVert_{2}^{2}d\mu_{h}<\epsilon,$ $$

for all $t\in[1,2,\dots,T]$.

For prediction of the system output, as the identity function $Id({\bm{h}})\mapsto{\bm{h}}$ is in $L^{2}_{d_{h}}$, the result above implies convergence of $\int_{\mathcal{H}}\lVert CK_{N}^{t}\mathcal{F}_{M}-F^{t}\rVert_{2}^{2}d\mu_{h}$. The proof can be found in Section B.1, and in Section B.2 we discuss the limitations of the result w.r.t. the controlled setting.

###### Assumption 1 .

###### Theorem 1 .

## 4 Computational experiments

We now discuss a series of experiments designed to illustrate the benefits and challenges of our construction approach. We compare our method to the state-of-the-art iterative gradient-based method called shPLRNN, which we explain in Appendix D. For the real-world weather dataset, we also compare our approach with a long short-term memory (LSTM) model. Furthermore, due to the similarities our method bears with reservoir models, we also compare with an established reservoir model, namely an echo state network (ESN), further explained in Appendix E.
Details on the datasets can be found in Appendix F, hyperparameters for all models are given in Appendix H, the evaluation metrics are explained in Appendix G and a further comparison discussion as well as the hardware details are provided in Appendix H. The code to reproduce the experiments from the paper, as well as an up-to-date code base, can be found at

[https://gitlab.com/felix.dietrich/swimrnn-paper](https://gitlab.com/felix.dietrich/swimrnn-paper)

[https://gitlab.com/felix.dietrich/swimrnn](https://gitlab.com/felix.dietrich/swimrnn)

In Table 1, we list the quantitative results of the experiments without control. Each entry stems from five different runs, where the random seed is changed in order to ensure a more robust result. We give the mean over these five runs, as well as the minimum and maximum among them.

**Table 1: Results from computational experiments. We report the training time and MSE (mean squared error) or EKL (empirical Kullback–Leibler divergence, see Appendix G) for a sampled RNN (our approach), a reservoir model ESN (see LABEL:{sec:reservoirs}), a state of the art backpropagation-based RNN called shPLRNN (see Appendix D), and a long short-term memory (LSTM) model.**
| Example | Model | Time [s]: avg (min, max) | MSE: avg (min, max) |
| --- | --- | --- | --- |
| Van der Pol | sampled RNN | $0.26$ ($0.18$, $0.35$) | 9.55e-4 (7.08e-4, 1.28e-3) |
|  | ESN | $3.76$ ($2.29$, $5.40$) | 1.58e-2 (1.15e-2, 2.07e-2) |
|  | shPLRNN | $217.93$ ($203.10$, $251.51$) | 1.39e-2 (5.66e-3, 3.00e-2) |
| 1D Van der Pol | sampled RNN | $0.29$ ($0.25$, $0.31$) | 5.06e-3 (1.57e-4, 1.58-2) |
| Weather (day) | sampled RNN | $4.87$ ($4.83$, $4.92$) | $2.239^{\circ}C$ ($2.088^{\circ}C$, $2.392^{\circ}C$) |
|  | LSTM | $378.10$ ($284.00$, $421.30$) | $2.531^{\circ}C$ ($2.183^{\circ}C$, $2.754^{\circ}C$) |
|  | shPLRNN | $321.76$ ($298.98$, $384.46$) | $2.296^{\circ}C$ ($1.803^{\circ}C$, $2.548^{\circ}C$) |
| Weather (week) | sampled RNN | $4.87$ ($4.81$, $4.90$) | $4.624^{\circ}C$ ($4.169^{\circ}C$, $4.867^{\circ}C$) |
|  | LSTM | $628.10$ ($580.70$, $648.50$) | $4.544^{\circ}C$ ($3.964^{\circ}C$, $4.893^{\circ}C$) |
|  | shPLRNN | $830.84$ ($796.55$, $847.32$) | $2.604^{\circ}C$ ($2.500^{\circ}C$, $2.801^{\circ}C$) |
| Example | Model | Time [s]: avg (min, max) | EKL: avg (min, max) |
| Lorenz-63 | sampled RNN | $1.67$ ($1.34$, $1.92$) | 4.36e-3 (3.66e-3, 5.36e-3) |
|  | ESN | $3.54$ ($2.87$, $4.47$) | 8.73e-3 (7.20e-3, 1.06e-2) |
|  | shPLRNN | $607.42$ ($581.39$,$650.56$) | 5.79e-3(4.41e-3,7.56e-3) |
| Rössler | sampled RNN | $5.36$ ($4.39$, $6.39$) | 1.57e-4 (5.86e-5, 3.82e-4) |
|  | ESN | $8.11$ ($7.94$, $8.31$) | 8.33e-5 (3.79e-5, 2.25e-4) |
|  | shPLRNN | $866.17$ ($848.56$, $939.06$) | 6.53e-4 (4.35e-4,1.09e-3) |

### 4.1 Simple ODEs: Van der Pol Oscillator

We consider the Van der Pol oscillator system for a simple illustration of our method. A sampled RNN with a $tanh$ activation and a single hidden layer of width 80 is used, and the prediction method follows Algorithm 3.
The model is evaluated on test data; the averaged error and training time are reported in Table 1. One trajectory from the test set is visualized in Figure 2.
It should be noted that predictions start with an initial condition from the test dataset which is used to make the first prediction, and afterwards continue using this prediction as an input to predict the next state, without information from the ground-truth dataset.
In the results we observe a very stable trajectory over a long prediction horizon, and furthermore all eigenvalues of the Koopman operator are inside the unit disk (see Appendix C), thus we are certain that the model is stable. This experiment is also significant due to the periodic nature of the system, which is captured with our model, although neural network architectures in general struggle to capture periodicity (Ziyin et al., ). Compared to the gradient-descent trained shPLRNN our method is much faster and achieves higher forecasting accuracy, with lower MSE as prediction error. Compared with an ESN, our model has a shorter fit time and a smaller error. In addition, the hyperparameter search is much simpler and less time consuming for our method since there are fewer hyperparameters to tune.

### 4.2 Example with time delay embedding: Van der Pol Oscillator

For many real-world examples, it is not possible to observe the full state of a system. Here, we use the same datasets as in the simple Van der Pol experiment (Section 4.1) but only consider the first coordinate $h_{1}$. We embed the data using a time-delay embedding of six followed by a principal component analysis (PCA) projection which reduces the dimensionality to two. A sampled RNN with $tanh$ activation and a single hidden layer of width 80 is trained. Predicted trajectories from the initial test dataset state are shown in the bottom row of the right column of Figure 2. The fit time and MSE error are provided in Table 1, they are fairly similar to the fit time and error for the example where the full state is observed indicating that this model also captures the true dynamics.

Figure: Figure 2: Comparison of true and predicted trajectories fror the Van der Pol experiments are shown for a test trajectory. Left: state space representation. Right: the top two rows show the full state system’s first and second coordinate from Section 4.1, and the bottom most row shows the partially observed system from Section 4.2.
Refer to caption: x2.png

### 4.3 Examples of chaotic dynamics: Lorenz and Rössler systems

Chaotic systems pose a challenging forecasting problem from the class of dynamical systems. As an example, we consider the well-known Lorenz system in the chaotic regime. A sampling network with a $tanh$ activation and a single hidden layer of width 200 is trained. Predictions for a test trajectory are visualized in Figure 3. To evaluate the model, we do not calculate MSE since it is not suitable for chaotic trajectories, but instead an empirical KL divergence (EKL), to compare the orbits. For details on this geometric measure see Appendix G. The averaged EKL and training time are reported in Table 1. Our sampled RNN achieves a comparable performance with the reservoir model. However, it should be noted that training a reservoir model on Lorenz data is a well-studied problem, and the choice of hyperparameters has been tuned carefully to achieve excellent performance. When choosing hyperparameters for our sampling RNN we found that the necessary effort is low, as there are not as many degrees of freedom as in a reservoir.
On the other hand, compared with the shPLRNN trained with gradient descent, we observe a much better performance both in terms of error and training time.

Figure: Figure 3: The results from the Lorenz experiment are shown for a test trajectory. Left: state space representation of true and predicted trajectories. Right: trajectories obtained from the Lorenz model described in Section 4.3.
Refer to caption: x3.png

Furthermore, we consider the Rössler system in the chaotic regime. We use a sampled RNN with 300 hidden layer nodes and $tanh$ activation. A predicted test trajectory is shown in Figure 4. Since the system is chaotic, we calculate an EKL for our model and report it along with the training time in Table 1.
Our model requires a slightly shorter fit time than the ESN and achieves comparable performance in terms of the EKL error.
As depicted in Table 1, our method for chaotic systems is significantly faster and more accurate at forecasting, with lower prediction error EKL, compared to the gradient-based model. Additionally, training with iterative methods requires many hyperparameters that need careful tuning to achieve optimal performance. Efficiently finding the best hyperparameters can be very challenging, especially when there is a small amount of training data for chaotic trajectories, making the training of such data more difficult. However, our training method performs effectively even with a very small amount of training data.

Figure: Figure 4: Trajectories from the Rössler experiment are shown for a test trajectory. Left: state space representation of true and predicted trajectories. Right: trajectories obtained from the Rössler model described in Section 4.3.
Refer to caption: x4.png

### 4.4 Example with control inputs: Forced Van der Pol Oscillator

We consider again the Van der Pol oscillator, where now the second coordinate, $h_{2}$, is controlled with an external input $x$, and using a sampled RNN model we perform model predictive control (MPC) as in Algorithm 4. We let $\mathcal{G}_{\hat{M}}$ be the identity function, and let $B$ map from $\mathcal{X}$ to $\mathbb{R}^{M}$ (adding nonlinearity for $x$ did not yield different results for this system, and is described in Section H.1.1). The sampled recurrent neural network is then the identity $\mathcal{G}_{\hat{M}}$ and $\mathcal{F}_{M}$ is a hidden layer of width 128 and $tanh$ activation. This network is then passed as a surrogate model to a linear-quadratic regulator (LQR). The network, in combination with the LQR, can successfully steer the state to the target state (see Figure 5). We consider five different runs, where only the random seed is varied, and obtain the mean controller cost to be 125.92 and the mean training time of 1.122 seconds. The $L^{2}$ norm of the state is also tracked over time, for five different runs we show the norms and the pointwise mean (over the runs) in Figure 5.
This experiment highlights a key advantage of our model, which allows for modelling a nonlinear system such as the Van der Pol oscillator using a linear controller such as LQR. This implies that the well-established tools from linear control theory can be applied to non-linear systems using our method.

Figure: Figure 5: Controlled (i.e. forced) Van der Pol experiment (Section 4.4) for initial condition ${\bm{h}}_{0}=[-1.5,-1]^{\mkern-1.5mu\mathsf{T}}$. Left: state space representation of controlled and uncontrolled trajectories. Right: $L^{2}$ norm of the controlled trajectory for five different runs and the $L^{2}$ norm of the target state.
Refer to caption: x5.png

### 4.5 Example with real-world data

Figure: Figure 6: The predictions of the best models on the test dataset for the horizon of one week: sampled RNN (left), LSTM (middle), and shPLRNN (right).
Refer to caption: x6.png

##### Weather data

We apply our approach to the climate data presented in TensorFlow (). The dataset contains a time series of 14 weather parameters recorded in Jena (Germany) between January 1st, 2009, and December 31st, 2016. We used the largest consecutive time range for the experiments and trained a sampled RNN to predict the temperature parameter (Figure 6). We also trained an autoregressive LSTM model described in TensorFlow () and a shPLRNN as iterative baselines.

The data offers freedom in choosing the sizes of the time delay and prediction horizons. We decided to fix the time delay to one week and set two separate experiments with a prediction horizon of one day and one week. In the sampled RNN and LSTM experiments, we performed a grid search for each model with hyperparameters specified in Section H.4.1.

Table 1 shows the averaged training time and error metrics for the two selected horizons (day and week). We observe that the sampled and iterative models perform similarly in the case of a shorter horizon, while sampling offers much faster training. When considering a one-week horizon, shPLRNN predictably outperforms LSTM, but sampled RNN is still orders of magnitude faster. We also note the prediction horizon does not influence the training time of the sampled model that agrees with the Algorithm 2. When comparing predictions for the longer horizon in Figure 6, we notice that the LSTM struggles to predict the high-frequency fluctuations of the measurement, but the sampled RNN and shPLRNN successfully capture them. Figure 6 also highlights the deficiency of the MSE metric because a low mean error does not always correspond to accurate predictions, as illustrated by all three models. Overall, we conclude that sampled RNNs can successfully capture chaotic real-world dynamics and produce results comparable to the iterative models while offering a significant speed-up in training.

## 5 Conclusion

We introduce an efficient and interpretable training method for recurrent neural networks by combining ideas from random feature networks and Koopman operator theory.

Benefits of the approach
In many examples, we demonstrate that we can train accurate recurrent networks orders of magnitude faster compared to networks trained by iterative methods.
We also observe that the training approach works with a very small amount of training data.
The direct connection to Koopman operator theory allows us to draw on existing theoretical results for dynamical systems, which we use to prove convergence in the limit of infinite width. In addition, the same connection allows us to apply linear control theory to non-linear systems.

Limitations compared to other methods
The training method we use involves the solution of a large, linear system. The complexity of solving this system depends cubically on the minimum number of neurons and the number of data points (respectively, time steps). This means if both the network and the number of data points grow together, the computational time and memory demands for training grow too quickly.
With backpropagation-in-time, the memory requirements are mostly because many gradients must be stored for one update pass.

Remaining challenges and future work
Recurrent networks are often used for tasks in computer vision and natural language processing. These tasks require network architectures beyond feed-forward networks, like convolutional neural networks or transformers. The sampling scheme we use to construct the hidden weights of the neurons is currently not useful in constructing parameters for such architectures, but it is an intriguing challenge to work towards sampling them.
Remaining challenges in the theoretical work include extending the theory shown in this paper to controlled systems, as well as bridging the Koopman theory for continuous dynamical systems and NeuralODEs.

## 6 Ethics statement

For comparatively small problems, our construction method is orders of magnitude faster to construct the parameters of a recurrent network. As neural networks are generally dual-use, our work potentially also allows faster training for misuse of this technology.
Still, we connect the construction of recurrent neural networks to Koopman operator theory and dynamical systems. This connection allows researchers in these fields to better understand the behavior, failure modes, and robustness of recurrent architectures. In addition, by sampling weights and bias from the input spaces it adds interpretability to the models which again lets users understand better the underlying networks used (see Appendix C for an example of this).
We believe that this far outweighs the potential downsides of misuse because recurrent architectures that are understood much better can also be regulated in a more straightforward way.

## 7 Reproducibility statement

The assumptions for the theoretical results are stated in the main paper and further discussed in the appendix. The complete, detailed proofs can be found in the Appendix B and an introduction to Koopman theory can be found in Appendix A. For the computational experiments, we provide all hyperparameters in Appendix H. All reported results that are based on five runs with different random seeds, to ensure robustness. The code is available open-source, and can be found at
[https://gitlab.com/felix.dietrich/swimrnn-paper](https://gitlab.com/felix.dietrich/swimrnn-paper).

## Acknowledgments

E.B., I.B., and F.D. are thankful to the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) - for funding through project no. 468830823, and also acknowledge the association to DFG-SPP-229. F.D. and A.C. are supported by the TUM Georg Nemetschek Institute - Artificial Intelligence for the Built World. Z.M. is grateful to the Bundesministerium für Bildung und Forschung (BMBF, Federal Ministry of Education and Research) for funding through project OIDLITDSM, No. 01IS24061. We also thank Lukas Eisenmann and Florian Hess for their helpful suggestions and guidance in training the shPLRNNs used in this paper. Furtheremore, we are grateful to Daniel Lehmberg for developing and maintaining the datafold library which was useful for developing the codebase and the sampled RNN experiments.

## Appendix

## Appendix A Koopman operator and extended dynamic mode decomposition

In this section we give a more thorough introduction to the Koopman operator and its use in dynamical system theory. In addition, we also walk the reader through extended dynamic mode decomposition (EDMD), which is the finite approximation of the Koopman operator we use in the main paper.

Consider a dynamical system $(\mathcal{H},F)$, where the state space $\mathcal{H}$ is a topological space and $F\colon\mathcal{H}\rightarrow\mathcal{H}$ is an evolution operator/map. We impose more structure on our system by requiring our state space $\mathcal{H}$ to be a measure space $(\mathcal{H},\mathcal{F},\mu)$, and $F$ to be $\mathcal{F}$-measurable. We start by considering $(\mathcal{H},F)$ to be a discrete-time system, which are the systems we mainly work with in this paper. We can then write the evolution as

$$ $\displaystyle{\bm{h}}_{t+1}\,=\,f({\bm{h}}_{t}),\quad{\bm{h}}_{t}\in\mathcal{H }\subseteq\mathbb{R}^{d_{h}},\quad t\geq 0.$ (8) $$

The analysis of the evolution in the original state space can be hard, especially when $F$ is non-linear. The Koopman theory takes a different approach, and looks at the evolution of observables (e.g. measurements of the states) instead of the states themselves. The observables one is working with are typically $\phi\colon\mathcal{H}\rightarrow\mathbb{C}$ in a suitable Hilbert space $\mathcal{H}$. The evolution of the observables is then captured by the Koopman operator $\mathcal{K}$,

$$ $\displaystyle[\mathcal{K}\,\phi]({\bm{h}}):=(\phi\circ F)({\bm{h}}).$ $$

As long as the space of observables is a vector space, the Koopman operator is linear and we may analyse the dynamical system with non-linear $F$ using spectral analysis, with the caveat that in the majority of cases $\mathcal{K}$ is infinite dimensional. The choice of the function space $\mathcal{F}$ is crucial, as $\phi\circ F$ must belong to $\mathcal{F}$ for all $\phi\in\mathcal{F}$. Assuming $F$ is measure-preserving — which is common in ergodic theory — one can address the issue by setting

$$ $\displaystyle\mathcal{F}=L^{2}(\mathcal{H},\mu_{h})\coloneq\left\{\phi\colon \mathcal{H}\rightarrow\mathbb{F}\,\,\middle|\,\,\left\lVert\phi\right\rVert_{L ^{2}(\mathcal{H},\mu_{h})}=\left(\int_{\mathcal{H}}\lvert\phi({\bm{h}})\rvert^ {2}\,\mu_{h}(d{\bm{h}})\right)^{\frac{1}{2}}<\infty\right\}.$ (9) $$

where $\mathbb{F}=\mathbb{R}$ or $\mathbb{C}$. As we consider spectral analysis in this section, we let $\mathbb{F}=\mathbb{C}$. Since $F$ is measure-preserving means $\mathcal{K}$ is an isometry and the issue is resolved. The map $\mathcal{K}$ might still not be well-defined, as $\phi_{1},\phi_{2}\in L^{2}(\mathcal{H},\mu_{h})$ may differ only on a null set, yet their images under $\mathcal{K}$, could differ over a set of positive measure. To exclude this possibility, $F$ must be $\mu_{h}$-nonsingular, meaning that for each $H\subseteq\mathcal{H}$, $\mu_{h}(F^{-1}(H))=0\,$ if $\mu_{h}(H)=0$.

Once the function space $\mathcal{H}$ is chosen, making sure that $\mathcal{K}$ is well-defined, we may apply spectral analysis. A Koopman eigenfunction $\varphi_{k}\in L^{2}(\mathcal{H},\mu_{h})$ corresponding to a Koopman eigenvalue $\lambda_{k}\in\sigma(\mathcal{K})$ satisfies

$$ $\displaystyle\varphi_{k}({\bm{h}}_{t+1})\,=\,\mathcal{K}\varphi_{\lambda}({\bm {h}}_{t})\,=\,\lambda_{k}\,\varphi_{\lambda}({\bm{h}}_{t}).$ $$

When the state space $\mathcal{H}\subseteq\mathbb{R}^{d}$, under certain assumptions of the space of eigenfunctions, we can evolve ${\bm{h}}$ using the spectrum of $\mathcal{K}$. More concretely, let $\Phi\colon\mathcal{H}\rightarrow\mathbb{C}^{d}$ be a vector of observables, where each observable $\phi_{i}({\bm{h}})={\bm{h}}_{i}$, where ${\bm{h}}_{i}$ is the $i$th component of ${\bm{h}}$. Assuming ${\phi_{i}\in\text{Span}\{\varphi_{k}\}\subset L^{2}(\mathcal{H},\mu_{h})\,}$, we have

$$ $\displaystyle\mathcal{K}\,\phi_{i}\,=\,\mathcal{K}\sum_{k}c_{k}\,\varphi_{k}\, =\,\sum_{k}c_{k}\,\mathcal{K}\,\varphi_{k}\,=\,\sum_{k}c_{k}\,\lambda_{k}\, \varphi_{k},$ $$

due to the linearity of $\mathcal{K}$. Iterative mapping of $\mathcal{K}$ yields

$$ $\displaystyle{\bm{h}}_{t}=\Phi({\bm{h}}_{t})\,=\,(\underbrace{[\mathcal{K}\phi ]\circ\cdots\circ[\mathcal{K}\phi]}_{t})({\bm{h}}_{0})=\sum_{k}\lambda^{t}_{k} \,\,\phi_{\lambda_{k}}({\bm{h}}_{0})\,\,{\bm{c}}^{\Phi}_{k}.$ (10) $$

This process is known as Koopman mode decomposition (KMD), and the vectors ${\bm{c}}^{\Phi}_{k}\in\mathbb{C}^{d}$ are referred to as Koopman modes associated with the observable $\Phi$. This reveals one of the true strengths of the Koopman theory.

Up until now we have considered a discrete dynamical system, but the Koopman theory can also be extended to continuous systems where ${\bm{h}}$ is a function of time $t$ and its evolution is given by

$$ $\displaystyle\dot{\bm{h}}(t)=F({\bm{h}}(t)),\quad{\bm{h}}(t)\in\mathcal{H},\,t \in\mathbb{R}_{\geq 0},$ $$

by using the flow map. For any $t$, the flow map operator denoted by $F^{t}\colon\mathcal{H}\rightarrow\mathcal{H}$ is defined as

$$ $\displaystyle{\bm{h}}(t)=F^{t}({\bm{h}})={\bm{h}}(0)+\int_{0}^{t}F({\bm{h}}( \tau))d\tau,$ $$

which maps from an initial condition ${\bm{h}}(0)$ to point on the trajectory $t$ time units away. We can then define the Koopman operator for each $t\in\mathbb{R}_{\geq 0}$ ,

$$ $\displaystyle[\mathcal{K}^{t}\,\phi]({\bm{h}})=(\phi\circ F^{t})({\bm{h}}),$ $$

where $\phi\in L^{2}(\mathcal{H},\mu_{h})$. The set of all these operators $\{\mathcal{K}^{t}\}_{t\in\mathbb{R}_{\geq 0}}$ forms a semigroup with an infinitesimal generator $\mathcal{L}$. With some assumption on the continuity of the semigroup, the generator is the Lie derivative of $\phi$ along the vector field $F({\bm{h}})$ and can be written as

$$ $\displaystyle[\mathcal{L}\,\phi]({\bm{h}})=\lim_{t\downarrow 0}\frac{[\mathcal {K}^{t}\,\phi]({\bm{h}})-\phi({\bm{h}})}{t}=\frac{d}{dt}\phi({\bm{h}}(t))\Big{ |}_{t=0}=\nabla\phi\cdot\dot{\bm{h}}(0)=\nabla\phi\cdot F({\bm{h}}(0)).$ $$

The eigenfunction and eigenvalue is in the continuous time case scalars and functions fulfilling

$$ $\displaystyle[\mathcal{K}^{t}\varphi_{k}]({\bm{h}})=e^{\lambda_{k}t}\,\varphi_ {k}({\bm{h}}),$ $$

where $\{e^{\lambda_{k}}\}$ are the eigenvalues of the semigroup. This allows us to use the Koopman theory for continuous dynamical systems as well, and possibly make the connection for NeuralODEs and Koopman theory in similar fashion we have done with discrete system Koopman operator and RNN.

As the Koopman operator is infinite dimensional makes it impossible to apply it directly, and raises the need for a method to create a finite approximation of $\mathcal{K}$ and its spectrum, namely the extended dynamic mode decomposition.

### A.1 Extended dynamic mode decomposition

As we are mostly working with discrete systems in this paper, we focus on approximating the Koopman operator $\mathcal{K}$ for discrete dynamical systems. The way to approximate $\mathcal{K}$ is by extended dynamic mode decomposition (EDMD), which is an algorithm that provides a data driven finite dimensional approximation of the Koopman operator $\mathcal{K}$ through a linear map $K$. The spectral properties of $K$ subsequently serve to approximate those of $\mathcal{K}$. Utilizing this approach enables us to derive the Koopman eigenvalues, eigenfunctions, and modes. Here, we provide a brief overview of EDMD. For further details, refer to Williams et al. [].

The core concept of EDMD involves approximating the operator’s action on $\mathcal{F}\,=\,L^{2}(\mathcal{H},\mu_{h})$ by selecting a finite dimensional subspace $\widetilde{\mathcal{F}}_{M}\subset\mathcal{F}$. To define this subspace, we start by choosing a dictionary $\mathcal{F}_{M}\,=\,\{\psi_{i}:\mathcal{H}\rightarrow\mathbb{R}\in\mathcal{H}
\mid i=1,\cdots,M\}$. We have

$$ $\displaystyle\mathcal{F}_{M}({\bm{h}})=\begin{bmatrix}\psi_{1}({\bm{h}})&\psi_ {2}({\bm{h}})&\cdots&\psi_{M}({\bm{h}})\end{bmatrix}^{\mkern-1.5mu\mathsf{T}},$ $$

and we let the finite dimensional subspace $\widetilde{\mathcal{F}}_{M}$ be

$$ $\displaystyle\widetilde{\mathcal{F}}_{M}\,=\,\text{Span}\{\psi_{1},\psi_{2}, \cdots,\psi_{M}\}\,=\,\{{\bm{a}}^{\mkern-1.5mu\mathsf{T}}\mathcal{F}_{M}\colon {\bm{a}}\in\mathbb{C}^{M}\}\subset\mathcal{F}.$ $$

The action of the Koopman operator on $\phi\in\widetilde{\mathcal{F}}_{M}$ due to linearity is

$$ $\displaystyle\mathcal{K}\phi={\bm{a}}^{\mkern-1.5mu\mathsf{T}}\mathcal{K} \mathcal{F}_{M}={\bm{a}}^{\mkern-1.5mu\mathsf{T}}\mathcal{F}_{M}\circ F.$ $$

Assuming that the subspace $\widetilde{\mathcal{F}}_{M}$ is invariant under $\mathcal{K}$, i.e., $\mathcal{K}(\widetilde{\mathcal{F}}_{M})\subseteq\widetilde{\mathcal{F}}$, we can write $\mathcal{K}\phi={\bm{b}}^{\mkern-1.5mu\mathsf{T}}\mathcal{F}_{M}$ for any $\phi\in\widetilde{\mathcal{F}}_{M}$. It follows that $\mathcal{K}|_{\widetilde{\mathcal{F}}_{M}}$ is finite dimensional and can be written as a matrix $K\in\mathbb{R}^{M\times M}$ such that ${\bm{b}}^{\mkern-1.5mu\mathsf{T}}={\bm{a}}^{\mkern-1.5mu\mathsf{T}}K$. The equality can be seen, where $\phi\in\widetilde{\mathcal{F}}_{M}$, through

$$ $\displaystyle K\phi={\bm{a}}^{\mkern-1.5mu\mathsf{T}}K\mathcal{F}_{M}={\bm{b}} ^{\mkern-1.5mu\mathsf{T}}\mathcal{F}_{M}=\mathcal{K}\phi$ $$

where one first applies linearity of $K$, then the invariant assumption. When $\widetilde{\mathcal{F}}_{M}$ is not an invariant subspace of the Koopman operator $\mathcal{K}$, $K$ becomes an approximation. Decomposing $\mathcal{K}\phi\,=\,{\bm{b}}^{\mkern-1.5mu\mathsf{T}}\mathcal{F}_{M}+\rho$, where $\rho\in L^{2}(\mathcal{H},\mu_{h})$, EDMD approximates $\mathcal{K}$ by using data

$$ $\displaystyle H=\begin{bmatrix}{\bm{h}}_{1}&{\bm{h}}_{2}&\cdots&{\bm{h}}_{N} \end{bmatrix},\quad H^{\prime}=\begin{bmatrix}{\bm{h}}_{1}^{\prime}&{\bm{h}}_{ 2}^{\prime}&\cdots&{\bm{h}}_{N}^{\prime}\end{bmatrix},$ $$

where ${\bm{h}}_{n}^{\prime}=F({\bm{h}}_{n})$. Then, to find $K$, the EDMD minimizes the following cost function

$$ $\displaystyle\mathcal{J}\,=\,\frac{1}{2}\sum_{n=1}^{N}\left\lVert\rho({\bm{h}} _{n})\right\rVert^{2}=\frac{1}{2}\sum_{n=1}^{N}\left\lVert{\bm{a}}^{\mkern-1.5 mu\mathsf{T}}\left(\mathcal{F}_{M}({\bm{h}}_{n}^{\prime})-K\mathcal{F}_{M}({ \bm{h}}_{n})\right)\right\rVert^{2}.$ (11) $$

Letting

$$ $\displaystyle\mathcal{F}_{M}(H)$ $\displaystyle=\begin{bmatrix}\mathcal{F}_{M}({\bm{h}}_{1}),\mathcal{F}_{M}({ \bm{h}}_{2})\dots,\mathcal{F}_{M}({\bm{h}}_{N})\end{bmatrix}\in\mathbb{C}^{K \times N}$ $$

and equivalently for $\mathcal{F}_{M}(H^{\prime})$, a solution to Equation 11 is given by

$$ $\displaystyle K=\mathcal{F}_{M}(H^{\prime})\,\mathcal{F}_{M}(H)^{+}$ (12) $$

where $\mathcal{F}_{M}(H)^{+}$ denotes the pseudo-inverse.

Upon obtaining $K$, we find approximations of eigenfunctions

$$ $\displaystyle\varphi_{k}=\xi_{k}\mathcal{F}_{M},$ $$

where $\lambda_{k}$ and $\xi_{k}$ are eigenvalue and left eigenvector of $K$ respectively.
Finally, denoting $\bm{\varphi}\colon\mathcal{H}\rightarrow\mathbb{C}^{K}$ as the function ${\bm{h}}\mapsto\begin{bmatrix}\varphi_{1}({\bm{h}}),\varphi_{2}({\bm{h}}),
\dots,\varphi_{K}({\bm{h}})\end{bmatrix}$, we approximate the Koopman modes by

$$ $\displaystyle C=\operatorname*{arg\,min}_{\tilde{C}\in\mathbb{C}^{d\times K}} \lVert\Phi(H)-\tilde{C}\bm{\varphi}(H)\rVert_{Fr}^{2}=\operatorname*{arg\,min} _{\tilde{C}\in\mathbb{C}^{d\times K}}\lVert H-\tilde{C}\bm{\varphi}(H)\rVert_{ Fr}^{2},$ $$

where $\lVert\cdot\rVert_{Fr}$ is the Frobenius norm and $\Phi({\bm{h}})={\bm{h}}$, as defined in previous section. We may now approximate Equation 10 as

$$ $\displaystyle{\bm{h}}_{t}=\Phi({\bm{h}}_{t})=C\,\Lambda^{t}\bm{\varphi}({\bm{h }}_{0}),$ $$

where $\Lambda$ is a diagonal matrix with the eigenvalues. Due to numerical issues, one typically predicts one time step at the time, and projects down to state space each time.

In many applications, one is not necessarily interested in the spectral analysis, but only prediction. One then typically sets $\mathbb{F}=\mathbb{R}$, solves for $K$ in exact way, but solves for $C$ as

$$ $\displaystyle C=\operatorname*{arg\,min}_{\tilde{C}\in\mathbb{R}^{d\times K}} \lVert\Phi(H)-\tilde{C}\mathcal{F}_{M}(H)\rVert^{2}_{Fr}=\operatorname*{arg\, min}_{\tilde{C}}\lVert H-\tilde{C}\mathcal{F}_{M}(H)\rVert^{2}_{Fr}.$ $$

For prediction one simply applies $K$ several times,

$$ $\displaystyle{\bm{h}}_{t}=CK^{t}\mathcal{F}_{M}({\bm{h}}_{0}).$ $$

In practice, one usually predicts step by step ${\bm{h}}_{t}$, that is, maps it down to the state space and maps back to the observable image space, before applying $K$ again. It is also worth noting that one may be more interested in mapping to an output of a function $f\in\mathcal{F}$ instead, and then one simply swaps $\Phi$ with $f$ when approximating $C$.

To conclude this introduction to EDMD, we do want to mention that the set of eigenfunctions we find through the EDMD algorithm comes with its own set of issues, such as spectral pollution, and efforts to mitigate certain issues has spawned extensions to EDMD, e.g., measure-preserving EDMD and residual DMD [Colbrook, , Colbrook et al., ].

### A.2 Controlled dynamical systems and the Koopman operator

Extending Koopman theory to controlled systems can be done in several ways, and we opt to follow Korda and Mezić [] and limit ourselves to linear controlled systems,

$$ $\displaystyle{\bm{h}}_{t}=F({\bm{h}}_{t-1},{\bm{x}}_{t})=A_{h}{\bm{h}}_{t-1}+A _{x}{\bm{x}}_{t},\quad{\bm{y}}_{t}=A_{y}{\bm{h}}_{t}.$ (13) $$

Rewriting the input and the evolution operator slightly allows us to apply the Koopman operator and its theory as described in previous sections. Let

$$ $\displaystyle\tilde{\bm{h}}=\begin{bmatrix}{\bm{h}}\\ \tilde{\bm{x}}\end{bmatrix},$ $$

where ${\bm{h}}\in\mathcal{H}$ and $\tilde{\bm{x}}\in\ell(\mathcal{X})$, with $\ell(\mathcal{X})$ is the space of all countable sequences $\{{\bm{x}}_{i}\}_{i=1}^{\infty}$ such that ${\bm{x}}_{i}\in\mathcal{X}$. Then we can rewrite the evolution operator $F\colon\mathcal{H}\times\mathcal{X}\rightarrow\mathcal{H}$ to $\tilde{F}\colon\mathcal{H}\times\ell(\mathcal{X})\rightarrow\mathcal{H}$, where

$$ $\displaystyle\tilde{\bm{h}}_{t}=\tilde{F}(\tilde{\bm{h}}_{t-1})=\begin{bmatrix }F({\bm{h}}_{t-1},\tilde{\bm{x}}(0))\\ \mathcal{S}\,\tilde{\bm{x}}\end{bmatrix},$ $$

where $\tilde{\bm{x}}(i)={\bm{x}}_{i}\in\tilde{\bm{x}}$ and $\mathcal{S}$ is the left-shift operator, i.e., $\mathcal{S}(\tilde{\bm{x}}(i))=\tilde{\bm{x}}(i+1)$. The Koopman operator $\mathcal{K}$ can be applied to $\tilde{F}$ with observables $\phi\colon\mathcal{H}\times\ell(\mathcal{X})\rightarrow\mathbb{C}$, and the rest of the Koopman theory follows.

When approximating the Koopman operator for controlled systems with EDMD, the dictionary we choose needs alteration due to the domain $\mathcal{H}\times\ell(\mathcal{X})$ being infinite dimensional. Korda and Mezić [] propose dictionaries that are both computable and enforce the linearity relationship assumed in Equation 13. The dictionaries to be considered are on the form

$$ $\displaystyle\psi_{i}({\bm{h}},\tilde{\bm{x}})=\psi_{i}^{(h)}({\bm{h}})+\psi_{ i}^{(x)}(\tilde{\bm{x}}),$ $$

where $\psi_{i}^{(x)}\colon\ell(\mathcal{X})\rightarrow\mathbb{R}$ is a linear functional and $\psi_{i}^{(h)}\in\mathcal{F}$. The new dictionary can be written as

$$ $\displaystyle\mathcal{F}_{M}=\{\psi_{1}^{(h)},\dots,\psi^{(h)}_{M}\},\quad \mathcal{G}_{\tilde{M}}=\{\psi_{1}^{(x)},\dots,\psi^{(x)}_{\tilde{M}}\}.$ $$

Note the number of observables $M$ and $\tilde{M}$ can differ, even though in the main paper we typically set $\tilde{M}=M$. If they differ, the matrix $B$ will map from $\mathbb{F}^{\tilde{M}}$ to $\mathbb{F}^{M}$. Once the dictionaries are set, we simply solve the optimization problem

$$ $\displaystyle\operatorname*{arg\,min}_{K\in\mathbb{F}^{M\times M},B\in\mathbb{ F}^{M\times\tilde{M}}}\frac{1}{2}\sum_{n=1}^{N}\lVert\mathcal{F}_{M}({\bm{h}}_ {n}^{\prime})-(K\mathcal{F}_{M}({\bm{h}}_{n})+B\mathcal{G}_{\tilde{M}}({\bm{h} }_{n}))\rVert^{2}$ $$

with the analytical solution being

$$ $\displaystyle[K,B]=\mathcal{F}_{M}(H^{\prime})[\mathcal{F}_{M}(H),\mathcal{G}_ {\tilde{M}}(H)]^{+}.$ $$

Approximating $C$ is done in the same manner as for uncontrolled systems, as it only needs to learn how to map from $\mathcal{F}_{M}(\mathcal{H})$ to $\mathcal{H}$. For further details, see Korda and Mezić [].

## Appendix B Theory

In this section we give the necessary assumptions and proofs for the theoretical results in the main paper. We start by defining the state space $\mathcal{H}$ and input space $\mathcal{X}$, following the setup from Bolager et al. []. Letting

$$ $\displaystyle d_{\mathbb{R}^{d_{z}}}({\bm{z}},A)=\inf\{d({\bm{z}},{\bm{a}}) \colon{\bm{a}}\in A\}.$ $$

where $d$ is the canonical Euclidean distance in the space $\mathbb{R}^{d_{z}}$. The medial axis is defined as

$$ $\displaystyle\text{Med}(A)=\{{\bm{h}}\in\mathbb{R}^{d_{z}}\colon\exists{\bm{p} }\neq{\bm{q}}\in A,\lVert{\bm{p}}-{\bm{z}}\rVert=\lVert{\bm{q}}-{\bm{z}}\rVert =d_{\mathbb{R}^{d_{z}}}({\bm{z}},A)\}$ $$

and the reach is the scalar

$$ $\displaystyle\tau_{A}=\inf_{{\bm{a}}\in A}d_{\mathbb{R}^{d_{z}}}({\bm{a}}, \text{Med}(A)),$ $$

i.e., the point in $A$ that is closest to the projection of points in $A^{c}$.

###### Definition 2 .

Let $\widetilde{\mathcal{H}}$ be a nonempty compact subset of $\mathbb{R}^{d_{h}}$ with reach $\tau_{\widetilde{\mathcal{H}}}>0$ and equivalently for $\widetilde{\mathcal{X}}\in\mathbb{R}^{d_{x}}$. The input space $\mathcal{H}$ is defined as

$$ $\displaystyle\mathcal{H}=\{{\bm{h}}\in\mathbb{R}^{d_{h}}\colon d_{\mathbb{R}^{ d_{h}}}({\bm{h}},\widetilde{\mathcal{H}})\leq\epsilon_{\mathcal{H}}\},$ $$

where $0<\epsilon_{\mathcal{H}}<\min\{\tau_{\widetilde{\mathcal{H}}},1\}$. Equivalently for $\mathcal{X}$,

$$ $\displaystyle\mathcal{X}=\{{\bm{x}}\in\mathbb{R}^{d_{x}}\colon d_{\mathbb{R}^{ d_{x}}}({\bm{x}},\widetilde{\mathcal{X}})\leq\epsilon_{\mathcal{X}}\},$ $$

where $0<\epsilon_{\mathcal{X}}<\min\{\tau_{\widetilde{\mathcal{X}}},1\}$.

###### Remark 2 .

This restriction to the type of state and input spaces we consider is sufficient to construct all neural networks of interest by choosing pair of points from the space in question, and constructing the weight and bias as in Equation 4. It is also argued in Bolager et al. [] that most interesting real-world application will contain some noise and make sure that the state and input spaces is approximately on the form given in the definition.

As we are not considering the eigenfunctions of the system, only prediction, and we are working with real valued neural networks, we are in the setting with $\mathbb{F}=\mathbb{R}$ in Equation 9.

###### Definition 2 .

###### Remark 2 .

### B.1 Uncontrolled systems

For uncontrolled systems we have that the evolution is described as ${\bm{h}}_{t}=F({\bm{h}}_{t-1})$. For the following theory, we have the following assumptions.

###### Assumption 2 .

The measure defining the space $\mathcal{F}=L^{2}(\mathcal{H},\mu_{h})$, we assume $\mu_{h}$ is regular and finite for compact subsets.

###### Assumption 3 .

The following assumptions are made for the dictionary and the underlying system $F$:

- 1.
Any dictionary $\mathcal{F}_{M}$ must fulfill $\mu_{h}\{{\bm{h}}\in\mathcal{H}\mid{\bm{c}}^{\mkern-1.5mu\mathsf{T}}\mathcal{F
}_{M}({\bm{h}})=0\}=0$, for all nonzero ${\bm{c}}\in\mathbb{R}^{M}$.
- 2.
$\mathcal{K}\colon\mathcal{F}\rightarrow\mathcal{F}$ is a bounded operator.

Assumption 2 is not very limited, as it holds for most measures we are interested in, such as measures absolutely continuous to the Lebesgue measure. For Item 3.1 we have the following result.

###### Lemma 1 .

Let $(\mathbb{R}^{d_{h}},\mathscr{B}(\mathbb{R}^{d_{h}}),\mu_{h})$ be a measurable space with $\text{supp}(\mu_{h})=\mathcal{H}$, and $\lambda$ be the Lebesgue measure for $\mathbb{R}^{d_{h}}$. If for all non-zero ${\bm{c}}\in\mathbb{R}^{d_{h}}$, the following holds:

- •
The set of functions $\{\psi_{1},\dots,\psi_{M}\}$ is linear independent.
- •
${\bm{c}}^{\mkern-1.5mu\mathsf{T}}\mathcal{F}_{M}=\sum_{m=1}^{M}c_{m}\psi_{m}$ is analytic on $\mathbb{R}^{d_{h}}$,
- •
$\mu_{h}\ll\lambda$,

then Item 3.1 holds. In particular, it holds when $\{\psi_{i}\}_{m=1}^{M}$ are independent tanh functions.

###### Proof.

Let ${\bm{c}}\in\mathbb{R}^{d_{h}}$ be any non-zero vector. As $\{\psi_{1},\dots,\psi_{M}\}$ are linear independent, we have ${\bm{c}}^{\mkern-1.5mu\mathsf{T}}\mathcal{F}_{M}\not\equiv 0$. As ${\bm{c}}^{\mkern-1.5mu\mathsf{T}}\mathcal{F}_{M}$ is analytic, we have the set

$$ $\displaystyle A=\{{\bm{h}}\in\mathbb{R}^{d_{h}}\mid{\bm{c}}^{\mkern-1.5mu \mathsf{T}}\mathcal{F}_{M}({\bm{h}})=0\},$ $$

to have measure zero, $\lambda(A)=0$, due to Proposition 1 by Mityagin []. We have $\lambda(A\cap\mathcal{H})\leq\lambda(A)=0$. Finally due to absolutely continuous measure $\mu_{h}$, we have

$$ $\displaystyle\mu_{h}(\{{\bm{h}}\in\mathcal{H}\colon{\bm{c}}^{\mkern-1.5mu \mathsf{T}}\mathcal{F}_{M}({\bm{h}})=0\})=\lambda(A\cap\mathcal{H})=0,$ $$

and hence Item 3.1 holds. As the linear projection and shift of bias is analytic on $\mathbb{R}^{d_{h}}$, tanh is analytic on $\mathbb{R}$, and analytic functions are closed under compositions, means Item 3.1 holds when $\mathcal{F}_{M}$ is a set of linearly independent neurons with tanh activation function.
∎

###### Remark 3 .

The requirement of the functions being analytic for the whole $\mathbb{R}^{d_{h}}$ can certainly be relaxed if necessary to an open and connected set $U$, s.t. $\mathcal{H}\subseteq U)$. Further relaxation can be made with some additional work. The result above also agrees with the claim made in Korda and Mezić [] about Item 3.1 holds for many measures and most basis functions such as polynomials and radial basis functions. Finally, we note that the independence requirement is easily true when we sample the neurons.

Item 3.2 is commonly enforced in Koopman theory when considering convergence of EDMD and for example holds when $F$ is Lipschitz, has Lipschitz invertible, and $\mu_{h}$ is the Lebesgue measure [Korda and Mezić, ].

We note

$$ $\displaystyle\mathcal{NN}_{[1,1:\infty]}=\bigcup_{M=1}^{\infty}\mathcal{NN}_{1 }(\mathcal{H},\mathbb{R}^{M})$ $$

is the space of all hidden layers with tanh activation function from $\mathcal{H}$. We continue with a result relating this space with $\mathcal{F}$.

###### Lemma 2 .

For tanh activation function and when Assumption 2 holds, then $\mathcal{NN}_{[1,1:\infty]}$ is dense in $\mathcal{F}$ and has a countable basis $\{\psi_{i}\in\mathcal{NN}_{[1,1:\infty]}\}_{i=1}^{\infty}$, both when the parameter space is the full Euclidean space and when constructed as in Equation 4.

###### Proof.

It is well known that such space is dense in $C(\mathcal{H},\mathbb{R}^{K})$ for any $K\in\mathbb{N}$. This holds both when the weight space is the full Euclidean space [Cybenko, , Pinkus, ] and when limited to the weight construction in Equation 4 [Bolager et al., ]. As $\mathcal{H}$ is compact, we have that $\mathcal{NN}_{[1,1:\infty]}$ is dense in $\mathcal{F}$. Furthermore, as $\mathcal{F}$ is a separable Hilbert space and metric space, there exists a countable subset $\{\psi_{i}\in\mathcal{NN}_{[1,1:\infty]}\}_{i=1}^{\infty}$ that is a basis for $\mathcal{F}$.
∎

The following lemma makes sure we can circumvent assumptions made in Korda and Mezić [], which requires on the dictionary in the EDMD algorithm to be an orthonormal basis (o.n.b.) of $\mathcal{F}$.

###### Lemma 3 .

Let $H,H^{\prime}$ be the dataset used in Equation 12. For every set of $M\in\mathbb{N}$ linearly independent functions $\mathcal{F}_{M}=\{\phi_{i}\}_{i=1}^{M}$ from a dense subset of $\mathcal{F}$ and any function $f=c^{\mkern-1.5mu\mathsf{T}}\mathcal{F}_{M}$, there exists a $\tilde{c}$ and matrix $V$ such that

$$ $\displaystyle\tilde{c}^{\mkern-1.5mu\mathsf{T}}\tilde{K}\tilde{\mathcal{F}}_{M }=c^{\mkern-1.5mu\mathsf{T}}K\mathcal{F}_{M}$ $$

and

$$ $\displaystyle\tilde{c}^{\mkern-1.5mu\mathsf{T}}\tilde{\mathcal{F}}_{M}=f=c^{ \mkern-1.5mu\mathsf{T}}\mathcal{F}_{M},$ $$

where $\tilde{\mathcal{F}}_{M}=[\tilde{\psi}_{1},\tilde{\psi}_{2},\dots,\tilde{\psi}_
{M}]$ are functions from an orthornormal basis $\{\tilde{\psi}_{i}\}_{i=1}^{\infty}$ of $\mathcal{F}$, and $K,\tilde{K}$ are the Koopman approximations for the dictionaries $\mathcal{F}_{M}$ and $\tilde{\mathcal{F}}_{M}$ respectively.

###### Proof.

As $\mathcal{F}$ is a separable Hilbert space and a metric space, there exists a countable basis $\{\phi_{i}\}_{i=1}^{M}\cup\{\phi_{i}\}_{i=M+1}^{\infty}$, and by applying the Gram-Schmidt process to the basis, we have an o.n.b. $\{\tilde{\psi}_{i}\}_{i=1}^{\infty}$. Any $M$ step Gram-Schmidt process applied to a finite set of linearly independent vectors, can be written as a sequence of invertible matrices $V=\prod_{j=1}^{M+1}V_{j}$. Each matrix $V_{j}$ for $j<M+1$ transforms the $j$th vector and the last matrix simply scales. Constructing such matrix $V$ applied to $\mathcal{F}_{M}$ yields $\tilde{\mathcal{F}}_{M}$. Setting $\tilde{c}^{\mkern-1.5mu\mathsf{T}}=c^{\mkern-1.5mu\mathsf{T}}V^{-1}$, which means

$$ $\displaystyle\tilde{c}^{\mkern-1.5mu\mathsf{T}}\tilde{\mathcal{F}}_{M}=c^{ \mkern-1.5mu\mathsf{T}}V^{-1}\tilde{\mathcal{F}}_{M}=c^{\mkern-1.5mu\mathsf{T} }\mathcal{F}_{M}=f.$ $$

Furthermore, we have

$$ $\displaystyle\tilde{c}^{\mkern-1.5mu\mathsf{T}}\tilde{K}\tilde{\mathcal{F}}_{M}$ $\displaystyle=c^{\mkern-1.5mu\mathsf{T}}V^{-1}[\tilde{\mathcal{F}}_{M}(H^{ \prime})\tilde{\mathcal{F}}_{M}(H)^{+}]V\mathcal{F}_{M}$ $\displaystyle=c^{\mkern-1.5mu\mathsf{T}}V^{-1}[V\mathcal{F}_{M}(H^{\prime}) \mathcal{F}_{M}(H)^{+}V^{-1}]V\mathcal{F}_{M}$ $\displaystyle=c^{\mkern-1.5mu\mathsf{T}}[\mathcal{F}_{M}(H^{\prime})\mathcal{F }_{M}(H)^{+}]\mathcal{F}_{M}=c^{\mkern-1.5mu\mathsf{T}}K\mathcal{F}_{M}.$ $$

∎

We are now ready to prove the Theorem 1 from the paper, namely the existence of networks for finite horizon predictions. We denote $\mathcal{F}^{K}$ as the space of vector valued functions functions $f=[f_{1},f_{2},\dots,f_{K}]$, where $f_{i}\in\mathcal{F}$ and $\lVert f\rVert=\sum_{k=1}^{K}\lVert f_{k}\rVert_{L^{2}}$. In addition, we let $K_{N}$ be the Koopman approximation for $\mathcal{F}_{M}$ where $N$ data points have been used to solve the least square problem.

###### Theorem 2 .

Let $f\in\mathcal{F}^{K}$, $H,H^{\prime}$ be the dataset with $N$ data points used in Equation 12, and Assumption 2 and Assumption 3 hold. For any $\epsilon>0$ and $T\in\mathbb{N}$, there exist an $M\in\mathbb{N}$ and hidden layers $\mathcal{F}_{M}$ with $M$ neurons and matrices $C$ such that

$$ $\displaystyle\lim_{N\rightarrow\infty}\int_{\mathcal{H}}\lVert CK_{N}^{t} \mathcal{F}_{M}-f\circ F^{t}\rVert_{2}^{2}d\mu_{h}<\epsilon,$ (14) $$

for all $t\in[1,2,\dots,T]$. In particular, there exist hidden layers and matrices $C$ such that

$$ $\displaystyle\lim_{N\rightarrow\infty}\int_{\mathcal{H}}\lVert CK_{N}^{t} \mathcal{F}_{M}-F^{t}\rVert_{2}^{2}d\mu_{h}<\epsilon.$ (15) $$

###### Proof.

W.l.o.g., we let $K=1$. Due to Lemma 2, we know there exist hidden layers $\mathcal{F}_{M}$ and vectors ${\bm{c}}$ such that $\lVert f_{m}-f\rVert_{L^{2}}^{2}<\epsilon_{2}$, where ${\bm{c}}^{\mkern-1.5mu\mathsf{T}}\mathcal{F}_{M}=f_{m}$ and

$$ $\displaystyle\epsilon_{2}<\frac{\epsilon}{2\cdot\max\{\lVert\mathcal{K}\rVert_ {op}^{2T},\lVert\mathcal{K}\rVert_{op}^{2}\}},$ $$

with $\lVert\cdot\rVert_{op}$ being the operator norm. This is possible due to Assumption 3 and Definition 2. We then have for any $t\in[1,2,\dots,T]$

$$ $\displaystyle\lim_{N\rightarrow\infty}\int_{\mathcal{H}}\lVert{\bm{c}}^{\mkern -1.5mu\mathsf{T}}K_{N}^{t}\mathcal{F}_{M}-\mathcal{K}^{t}\,f\rVert_{2}^{2}d\mu _{h}$ $\displaystyle\leq\lim_{N\rightarrow\infty}\int_{\mathcal{H}}\lVert{\bm{c}}^{ \mkern-1.5mu\mathsf{T}}K_{N}^{t}\mathcal{F}_{M}-\mathcal{K}^{t}\,f_{m}\rVert_{ 2}^{2}d\mu_{h}+\lVert\mathcal{K}^{t}\,f_{m}-\mathcal{K}^{t}\,f\rVert_{L^{2}}$ $\displaystyle\leq\lim_{N\rightarrow\infty}\int_{\mathcal{H}}\lVert{\bm{c}}^{ \mkern-1.5mu\mathsf{T}}K_{N}^{t}\mathcal{F}_{M}-\mathcal{K}^{t}\,f_{m}\rVert_{ 2}^{2}d\mu_{h}+\lVert f_{m}-f\rVert_{L^{2}}^{2}\max\{\lVert\mathcal{K}\rVert_{ op}^{2T},\lVert\mathcal{K}\rVert_{op}^{2}\}$ $\displaystyle<\frac{\epsilon}{2}+\frac{\epsilon}{2}=\epsilon,$ $$

where we use Theorem 5 in Korda and Mezić [] to bound $\lVert{\bm{c}}^{\mkern-1.5mu\mathsf{T}}K_{N}^{t}\mathcal{F}_{M}-\mathcal{K}^{t
}\,f_{m}\rVert_{2}^{2}d\mu_{h}$ (in theory we might need a larger $M$, which we simply set and the bound of $f_{m}-f$ still holds). From the convergence above, Equation 14 follows by definition of the Koopman operator. For Equation 15, simply note that $f({\bm{h}})={\bm{h}}$ is in $\mathcal{F}^{d_{h}}$ due to Definition 2, and the result holds.
∎

###### Assumption 2 .

###### Assumption 3 .

###### Lemma 1 .

###### Proof.

###### Remark 3 .

###### Lemma 2 .

###### Proof.

###### Lemma 3 .

###### Proof.

###### Theorem 2 .

###### Proof.

### B.2 Controlled systems

Proving convergence for controlled systems gives different challenges. The results above cannot easily be shown for controlled systems. The reason being that the dictionary space one uses is not a basis for the observables in the controlled setting, with the dynamical system extended by the left-shift operator, and the simplification made for EDMD in controlled systems. The results above may be extended, but the EDMD will not converge to the Koopman operator, but rather to $P^{\mu}_{\infty}\mathcal{K}_{\mathcal{F}_{\infty}}$, where $P^{\mu}_{\infty}$ is the $L^{2}(\mu)$ projection onto the closure of the dictionary space [Korda and Mezić, ]. However, results exists for the continuous controlled systems, with certain convergences for the generator. This is sadly not as strong as above, but an interesting path to connect RNNs/NeuralODEs to such theory [Nüske et al., ].

## Appendix C Interpretability using Koopman and sampling

Figure 7 shows two possible ways we can interpret our constructed RNN models, beyond what is usually possible for classical RNNs trained with gradient descent. On the left, we indicate which data pairs were chosen from the training set to construct neurons of the non-linear network ${\mathcal{F}}$. This can help to see if the “coverage” of the training set by neurons (resp. their associated data pairs) is reasonable, or if more neurons or data points are needed to cover highly non-linear regions.

On the right in Figure 7, we plot the locations of the eigenvalues of the Koopman matrix $K$ in our sampled RNN from Section 4.1. We can see that all eigenvalues are located on and inside the unit circle, indicating stable and oscillatory behavior.

Figure: Figure 7: Left: Interpretability of the hidden state network weights over the data domain in the Van der Pol example. The arrows (with initial $x_{1}$ and final points $x_{2}$) depict the pairs that were chosen during sampling. This means that they directly visualize where neurons are placed on the domain. Right: Eigenvalues of the Koopman matrix approximated with 80 neurons for the Van der Pol example.
Refer to caption: extracted/5966941/figures/VdP_interpretability_phase_space.png

## Appendix D ReLU-based RNNs for Dynamical Systems Modeling (DSM)

Most RNNs are parameterized discrete-time recursive maps of the given in Definition 1,

$$ $\displaystyle{\bm{h}}_{t}$ $\displaystyle=F_{hx}(F_{\mathcal{H}}({\bm{h}}_{t-1}),F_{\mathcal{X}}({\bm{x}}_ {t})),$ (16) $$

with latent states ${\bm{h}}_{t-1}$, optional external inputs ${\bm{x}}_{t}$. A piecewise linear RNN (PLRNN), introduced by Koppe et al. [], has the generic form

$$ $\displaystyle{\bm{h}}_{t}=W_{h}^{(1)}{\bm{h}}_{t-1}+W_{h}^{(2)}\sigma({\bm{h}} _{t-1})+{\bm{b}}_{0}+W_{x}{\bm{x}}_{t},$ (17) $$

where $\sigma({\bm{h}}_{t-1})=\max(0,{\bm{h}}_{t-1})$ is the element-wise rectified linear unit (ReLU) function, $W_{h}^{(1)}\in\mathbb{R}^{d_{h}\times d_{h}}$ is a diagonal matrix of auto-regression weights, $W_{h}^{(2)}\in\mathbb{R}^{d_{h}\times d_{h}}$ is a matrix of connection weights, the vector ${\bm{b}}_{0}\in\mathbb{R}^{d_{h}}$ represents the bias, and the external input is weighted by $W_{x}\in\mathbb{R}^{d_{h}\times d_{x}}$. Afterwards, Brenner et al. [] extended this basic structure by incorporating a linear spline basis expansion, referred to as the dendritic PLRNN (dendPLRNN)

$$ $\displaystyle{\bm{h}}_{t}=W_{h}^{(1)}{\bm{h}}_{t-1}+W_{h}^{(2)}\sum_{j=1}^{J} \alpha_{j}\,\sigma({\bm{h}}_{t-1}-{\bm{b}}_{j})+{\bm{b}}_{0}+W_{x}{\bm{x}}_{t},$ (18) $$

where $\{\alpha_{j},{\bm{b}}_{j}\}_{j=1}^{J}$ represents slope-threshold pairs, with $J$ denoting the number of bases. This expansion was introduced to increase the expressivity of each unit’s nonlinearity, thereby facilitating DSM in reduced dimensions. Moreover, Hess et al. [] proposed the following “1-hidden-layer” ReLU-based RNN, which they referred to as the shallow PLRNN (shPLRNN)

$$ $\displaystyle{\bm{h}}_{t}=W_{h}^{(1)}{\bm{h}}_{t-1}+W_{h}^{(2)}\sigma(W_{h}^{( 3)}{\bm{h}}_{t-1}+{\bm{b}}_{1})+{\bm{b}}_{0}+W_{x}{\bm{x}}_{t},$ (19) $$

where $W_{h}^{(1)}\in\mathbb{R}^{d_{h}\times d_{h}}$ is a diagonal matrix, $W_{h}^{(2)}\in\mathbb{R}^{d_{h}\times M}$ and $W_{h}^{(3)}\in\mathbb{R}^{M\times d_{h}}$ are rectangular connectivity matrices, and ${\bm{b}}_{1}\in\mathbb{R}^{M}$, ${\bm{b}}_{0}\in\mathbb{R}^{d_{h}}$ denote thresholds. The combination of Generalized Teacher Forcing (GTF) and shPLRNN results in a powerful DSM algorithm on challenging real-world data; for more information see Hess et al. []. When $M>d_{h}$, it is possible to rewrite any shPLRNN as a dendPLRNN by expanding the activation of each unit into a weighted sum of ReLU nonlinearities [Hess et al., ].

A clipping mechanism can be added to the shPLRNN to prevent states from diverging to infinity as a result of the unbounded ReLU nonlinearity

$$ $\displaystyle{\bm{h}}_{t}=W_{h}^{(1)}{\bm{h}}_{t-1}+W_{h}^{(2)}\big{[}\sigma(W _{h}^{(3)}{\bm{h}}_{t-1}+{\bm{b}}_{1})-\sigma(W_{h}^{(3)}{\bm{h}}_{t-1})\big{] }+{\bm{b}}_{0}+W_{x}{\bm{x}}_{t}.$ (20) $$

This guarantees bounded orbits under certain conditions on the matrix $W_{h}^{(1)}$ [Hess et al., ].

In our experiments with RNNs, the clipped shPLRNN is trained by GTF for DSM on the benchmark systems (see below Appendix F).

## Appendix E Reservoir Models: Echo State Networks

As our newly proposed method bears similarities to a reservoir computing architecture, we have trained reservoir models as part of our computational experiments. We used the simplest recurrent reservoir architecture, which is an echo state network (ESN) introduced by Jaeger and Haas []. Here we will briefly introduce the main ideas behind ESNs, and we refer the interested reader to the review by Lukoševičius and Jaeger [] for more details.

An ESN consists of a reservoir and a readout. The reservoir contains neurons which are randomly connected to inputs and are these are not trained. For a certain discrete time $t\in\mathbb{Z}_{\geq 0}$ denoting the inputs as ${\bm{h}}_{t}\in\mathbb{R}^{N_{h}}$ and output as ${\bm{y}}_{t}\in\mathbb{R}^{N_{y}}$, and the internal reservoir states as ${\bm{k}}_{t}\in\mathbb{R}^{N_{k}}$, the reservoir provides an update rule for the internal units as

$$ $\displaystyle{\bm{k}}_{t+1}=f\left(W^{in}{\bm{h}}_{t+1}+W{\bm{k}}_{t}+W^{back} {\bm{y}}_{t}\right),$ (21) $$

for an activation function $f$ and weight matrices $W^{in}\in\mathbb{R}^{N_{k}\times N_{h}}$, $W\in\mathbb{R}^{N_{k}\times N_{k}}$ and $W^{back}\in\mathbb{R}^{N_{k}\times N_{y}}$. After the reservoir comes the readout, which maps the inputs, reservoir states and outputs to a new output state

$$ $\displaystyle{\bm{y}}_{t+1}=f_{out}\left(W^{out}({\bm{h}}_{t+1},{\bm{k}}_{t+1} ,{\bm{y}}_{t})\right),$ (22) $$

where $f^{out}$ is the output activation, $W^{out}\in\mathbb{R}^{N_{y}\times N_{y}}$ are output weights and $({\bm{h}}_{t+1},{\bm{k}}_{t+1},{\bm{y}}_{t})$ denotes the concatenation of ${\bm{h}}_{t+1}$, ${\bm{k}}_{t+1}$ and ${\bm{y}}_{t}$.
In the readout the model learns the connections from the reservoir to the readout, for example via (regularized) regression. A so-called feedback connection allows for the readout values to be fed back into the reservoir, as shown in Equation 21, establishing a recurrent relation.

## Appendix F Benchmark systems and real-world data

##### Van der Pol

In 1927, Balthasar Van der Pol introduced a non-conservative oscillatory system with a nonlinear damping term to describe oscillations in a vacuum tube electrical circuit. The system is described as a two dimensional ODE

$$ $\displaystyle\begin{cases}\dot{h_{1}}&=h_{2}\\ \dot{h_{2}}&=\mu(1-h_{1}^{2})-h_{1}+h_{2},\end{cases}$ $$

where $\mu$ is a scalar parameter indicating the nonlinearity and the strength of the damping. In our experiment we set $\mu=1$, that is, in the limit cycle regime.
For the experiments in Section 4.1 and Section 4.2, training data are generated by solving an initial value problem for $t\in[0,20]$ with $\Delta t=0.1$ for 50 initial conditions, where each initial condition is random vector ${\bm{h}}_{0}\sim\text{Uniform}([-3,3]^{2})$. We used an explicit Runge-Kutta method of order 8 to solve the initial value problem. Validation and test data are generated similarly but for $t\in[0,50]$.

##### Lorenz-63

Devised by Edward Lorenz in 1963 [Lorenz, ] to model atmospheric convection, the Lorenz-63 system is defined as

$$ $\displaystyle\begin{cases}\dot{h_{1}}&=\sigma(h_{2}-h_{1})\\ \dot{h_{2}}&=h_{1}(\rho-h_{3})-h_{2}\\ \dot{h_{3}}&=h_{1}\,h_{2}-\beta h_{3},\end{cases},$ $$

where $\sigma,\rho,\beta$, are parameters that control the dynamics of the system. In our experiment, we set $\sigma=10$, $\beta=\frac{8}{3}$, and $\rho=28$, which means we are in the chaotic regime.
For the experiments in Section 4.3, training data are generated by solving an initial value problem for $t\in[0,5]$ with $\Delta t=0.01$ for 50 initial conditions, where each initial condition is random vector ${\bm{h}}_{0}\sim\text{Uniform}([-20,20]\times[-20,20]\times[0,50])$. The solver used explicit Runge-Kutta of order 8, likewise as the previous dataset. Validation and test data are generated similarly but for $t\in[0,50]$. We normalize datasets to scale the values to the range $[-3,3]$ to improve the training.

##### Rössler

Otto Rössler introduced the Rössler system in 1976 [Rössler, ] as a model that generates chaotic dynamics

$$ $\displaystyle\begin{cases}\dot{h_{1}}&=-h_{2}-h_{3}\\ \dot{h_{2}}&=h_{1}+\alpha h_{2}\\ \dot{h_{3}}&=\beta+h_{3}(h_{1}-\kappa),\end{cases}$ $$

where $\alpha$, $\beta$, $\kappa$, are parameters controlling the dynamics of the system. Here, we set $\alpha=0.15$, $\beta=0.2$, and $\kappa=10$, which puts the system in the chaotic regime. The setup for data generation is similar to the Lorenz example. For the experiments in Table 1, training data are generated by solving an initial value problem for $t\in[0,10]$ with $\Delta t=0.01$ for 50 initial conditions, where each initial condition is random vector ${\bm{h}}_{0}\sim\text{Uniform}([-20,20]\times[-20,20]\times[0,40])$. Again, an explicit Runge-Kutta method of order 8 was used. Validation and test data are generated similarly but for $t\in[0,200]$. We normalize datasets to scale the values to the range $[-3,3]$ to improve the training.

##### Forced Van der Pol Oscillator

As an example for a controlled system, we use the Van Der Pol oscillator with external input forcing $x$, and

$$ $\displaystyle\begin{cases}\dot{h_{1}}&=h_{2}\\ \dot{h_{2}}&=\mu(1-h_{1}^{2})h_{2}-h_{1}+x\\ \end{cases}$ $$

The data is obtained for $t\in[0,50\Delta t]$ with $\Delta t=0.05$, with 150 initial conditions, where ${\bm{h}}_{0}\sim\text{Uniform}([-3,3]^{2})$ and $x_{0}\sim\text{Uniform}([-3,3])$. The solver used here was using an explicit Runge-Kutta method of order 5(4). It is important to highlight that the control input data $x_{0}$ does not come from any controller with a particular target state, i.e. random control is applied to the trajectories in the training dataset.

##### Weather data

In this experiment, we follow TensorFlow [] and use the Jena Climate dataset [for Biogeochemistry, ]. The original data contains inconsistent date and time values, leading to gaps and overlaps between measurements. We extracted the longest consecutive time period and thus worked with the data between July 1st, 2010, and May 16th, 2013. We additionally downsampled the time series from the original 10-minute to 1-hour measurements. Then, the first 70% of records were used as the train set, the next 20% as the validation, and the remaining 10% as the test set. Identically to TensorFlow [], we pre-processed the features and added sin and cos time-embeddings of hour, day, and month. We plot the dataset, indicating the train-validation-test split with colors in Figure 8.

Figure: Figure 8: The weather dataset and its splits: train (blue), validation (orange), and test (green).
Refer to caption: x8.png

## Appendix G Evaluation measures

### G.1 Geometrical measure

The Kullback–Leibler divergence of two probability densities $p({\bm{x}})$ and $q({\bm{x}})$ is defined as

$$ $\displaystyle D_{KL}(p({\bm{x}})\lVert q({\bm{x}}))=\int_{{\bm{x}}\in\mathbb{R }}p({\bm{x}})\log\frac{p({\bm{x}})}{q({\bm{x}})}d{\bm{x}}.$ (23) $$

In order to be able to accurately evaluate also high-dimensional systems, we follow the approach used in Hess et al. [] and place Gaussian Mixture Models (GMM) on the along the true trajectory ${\bm{x}}$ and predicted $\hat{{\bm{x}}}$ trajectories, obtaining $\hat{p}({\bm{x}})=\frac{1}{T}\sum_{t=1}^{T}\mathcal{N}({\bm{x}},{\bm{x}}_{t},\Sigma)$ and $\hat{q}({\bm{x}})=\frac{1}{T}\sum_{t=1}^{T}\mathcal{N}({\bm{x}},\hat{{\bm{x}}}
_{t},\Sigma)$ for $T$ snapshots. Using the estimated densities, we consider a Monte Carlo approximation of Equation 23 by drawing $n$ random samples from the GMMs and obtain the density measure

$$ $\displaystyle D_{KL}(\hat{p}({\bm{x}})\lVert\hat{q}({\bm{x}}))\approx\frac{1}{ n}\sum_{i=1}^{n}\log\frac{\hat{p}({\bm{x}})}{\hat{q}({\bm{x}})}.$ (24) $$

We call this metric empirical KL divergence (EKL) in the manuscript.
To make our results comparable with Hess et al. [], we use $\sigma^{2}=1.0$ and $n=1000$.

## Appendix H Model details and comparison

### H.1 Sampled RNN

The implementation of our gradient-free training module for sampled RNNs was done in Python since the key tools for the algorithm already exist as Python libraries.
The algorithm requires the ability to sample weights and biases, thus we used the Python library swimnetworks by Bolager et al. []. Furthermore, we were able to alleviate the approximation of the Koopman operator, in the uncontrolled as well as controlled setting using some functionalities from the Python library datafold by Lehmberg et al. [].

Sampled RNNs have only a few hyperparameters: the number of nodes in the hidden layer, the activation function of the hidden layer, and a cutoff for small singular values in the least-squares solver. We often refer to the singular value cutoff hyperparameter as the regularization rate.
In the case of a sampled RNN with time delay, additional hyperparameters are the number of time delays and the number of PCA components, if used.

Sampled RNNs do not only have a low fit time but also a short hyperparameter tuning since there are only a few degrees of freedom. Thus, our newly proposed method offers efficiency beyond the short training time.

We empirically investigated the scaling of sampled RNNs as a function of the number of neurons. Similarly, as for the other results, we perform five runs and report the average as well as the minimum and maximum; the results can be seen in Figure 9. These are results for a dataset of fixed size, where we would expect that the only change in computation is the least-squares solver and additional samples of weights and biases.

Figure: Figure 9: Fit time over number of neurons for the Van der Pol experiment from Section 4.1.
Refer to caption: x9.png

In addition, we investigated the effect of increasing the number of neurons on training and validation performance. Results from five runs are shown in Figure 10. We notice that the errors decrease up to a certain number of neurons and stagnate afterwards.

Figure: Figure 10: Training and validation MSE error over number of neurons for the Van der Pol experiment from Section 4.1.
Refer to caption: x10.png

The Van der Pol experiment from Section 4.1 was used as one of the first prototypes of our method. Our initial experiments used a smaller dataset than the one described in Appendix F, namely with very short trajectories consisting of only three time steps. The sampled RNN was able to learn the full dynamics from this data. However, it was not possible to train the gradient-based shPLRNN with such a dataset; thus, we used longer trajectories to be able to compare our method with others. The final choice of hyperparameters for the hyperparameters is given in Table 2.

Training the 1D Van der Pol from Section 4.2 was done similarly, the number of hidden nodes and activation function remained unchanged. The hyperparameters are listed in Table 2. A time delay of six was chosen, however it was also possible to use only two time delays and get an excellent performance, but we opted against this choice because this was based on our knowledge that the true system is two dimensional. Thus, we opted for a higher number of time delays which then get reduced with the additional PCA transformation to mimic a real-world scenario where the true dimension of the problem is unknown.

The chaotic systems from Section 4.3 required more hidden nodes in order to obtain good prediction, as compared to the Van der Pol but nonetheless the process of hyperparameter tuning was simple. The final choice of hyperparameters for the Lorenz and Rössler problems is given in Table 2. We show a visual comparison of the predicted trajectories from the sampled RNNs, ESN and shPLRNNs for the Lorenz system in Figure 11.

Figure: (a) sampled RNN
Refer to caption: x11.png

With the forced Van der Pol example, presented in Section 4.4, we use the sampled RNN as a surrogate model to a controller, essentially making this a model-based control example. In order to be able to accurately control the state, the surrogate needs to capture the dynamics of the system with sufficient accuracy. The hyperparameters are listed in Table 2.
We opted for an LQR controller as our optimal controller, as a simple choice for linear systems. Although the Van der Pol oscillator is not a linear system, the surrogate sampled RNN captures the dynamics via a linear map with the Koopman operator, allowing minimizing the cost. For our LQR, $R=1$, and $Q=\text{diag}(10)$. The dataset we use contains randomly initialized trajectories evolved in time, as well as randomly chosen control inputs. We observed that a dataset with many shorter trajectories is more useful than a dataset with a few longer trajectories since for this problem diverse initial conditions are more informative of the vector field of the system, as opposed to long trajectories which become periodic. The necessary effort for tuning the sampled RNN hyperparameters and LQR controller parameters was low.

Although not reported in the manuscript as part of the experiments, in Section H.1.1 we discuss the forced Van der Pol model using a nonlinearity for the inputs.

The training of a sampling RNN on weather data was performed using time-delay embeddings to account for a possibly partial observation of the state. The objective is to train a model which can predict the temperature. To find the hyperparameters a grid search was performed, and this is documented in Table 7. The final choice of hyperparameters is given in Table 2. After training, predictions are done for fixed chunks of time, which are concatenated together afterwards. The number of time-delays dictates how many ground-truth datapoints are necessary to predict the next state, which gets concatenated with the previous predictions. We see this as a reasonable approach for weather prediction, as typically one would use the available information for a fairly long period of time, and predict for a fixed, likely shorter time horizon.

#### H.1.1 LQR with nonlinear lifting

In this section we describe shortly how we perform model predictive control described in Algorithm 4 with nonlinear lifting $\mathcal{G}_{\hat{M}}$, i.e. applying a sampling hidden layer on the input $x$. The experiment in the main paper was mainly done by setting $\mathcal{G}_{\hat{M}}$ to the identity, but adding nonlinearity with $\hat{M}>d_{x}$ is possible and may beneficial for future testing, even though for the particular system in the main paper, it yielded equivalent results as with $\mathcal{G}_{\hat{M}}\equiv Id$.

Following Algorithm 4, we lift ${\bm{x}}_{t}$ to $\mathbb{R}^{\hat{M}}$ through the hidden layer $\mathcal{G}_{\hat{M}}$. Then we solve for $B$ exactly as described in Algorithm 2 and fit the LQR. The LQR computes the control input, but now in the $\mathbb{R}^{\hat{M}}$ space. Before we can pass it to $F$, we approximate a linear map $P$, such that $P\mathcal{G}_{\hat{M}}({\bm{x}}_{t})\approx{\bm{x}}_{t}$. Once projected down, we can pass the projected value $\hat{\bm{x}}_{t}$ to $F$ and continue as usual.

#### H.1.2 Hyperparameters

**Table 2: Hyperparameters of sampled RNN models**
|  | Van der Pol | 1D Van der Pol | Lorenz | Rössler | forced Van der Pol | Weather |
| --- | --- | --- | --- | --- | --- | --- |
| Hidden layer width | $80$ | $80$ | $200$ | $300$ | $128$ | $256$ |
| Activation function | $tanh$ | $tanh$ | $tanh$ | $tanh$ | $tanh$ | $tanh$ |
| Regularization rate | 1e-8 | 0 | 1e-7 | 1e-4 | 1e-10 | 1e-6 |
| Time delays | - | $6$ | - | - | - | 168 |
| PCA components | - | $2$ | - | - | - | - |

#### H.1.3 Hardware

The machine used for training the sampled RNNs was 13th Gen Intel(R) Core(TM) i5-1335U @ 4.6 GHz (16GB RAM, 12 cores), no GPU hardware was used.

### H.2 ESN

For the implementation of ESNs we used the Python library reservoirpy by Trouvain et al. []. All models were trained using a single reservoir and a ridge regression readout. We consider only reservoir models without any warm-up phase in order to keep them comparable to our sampling RNN which does not have a warm-up.

For the chaotic systems, we were able to find guidelines in the literature for a suitable choice of hyperparameters specifically tailored to the Lorenz system (see [Viehweg et al., ]). With minor modifications and following the proposed guidelines, we found hyperparameters also for the Rössler system. The choice of hyperparameters is specified in Table 3, any hyperparameters not mentioned are set to the default value of reservoirpy.

For the Van der Pol problem many hyperparameter combinations were tried out until we were able to find a model with good performance and sufficient robustness to a change in the random seed. The hyperparameter combinations we considered are shown in Table 4. Due to the high expenses of a grid-search approach, a random search on a grid was employed and 1000 hyperparameter combinations were considered. The hyperparameter choice with the best performance on the validation dataset was selected and then evaluated on the test dataset. The final choice of hyperparameters is documented in Table 3, and any unmentioned hyperparameters are assumed to be set to the default value from the reservoirpy library.

#### H.2.1 Hyperparameters

The hyperparameters used for the ESN models are given in Table 3.

**Table 3: Hyperparameters of reservoir models**
|  | Van der Pol 4.1 | Lorenz 4.3 | Rössler 4.3 |
| --- | --- | --- | --- |
| Width/units | $500$ | $300$ | $500$ |
| Leak rate | $0.9$ | $0.3$ | $0.3$ |
| Spectral radius | $0.5$ | $1.25$ | $0.5$ |
| Input scaling | $0.05$ | $0.1$ | $0.1$ |
| Connectivity | $0.8$ | $0.1$ | $0.1$ |
| Inter connectivity | $0.2$ | $0.2$ | $0.2$ |
| Ridge regularization coeff. | 1e-10 | 1e-4 | 1e-8 |
| Warmup steps | 0 0 | 0 0 | 0 0 |

**Table 4: Hyperparameters used in random search on a grid for a Van der Pol reservoir model**
|  | Van der Pol 4.1 |
| --- | --- |
| Width/units | $100,200,500$ |
| Leak rate | $0.1,0.3,0.5,0.7,0.9$ |
| Spectral radius | $0.25,0.5,0.75,1,1.25,1.5,2,3,5$ |
| Input scaling | $0.05,0.1,0.5,1,1.5,2$ |
| Connectivity | $0.2,0.4,0.6,0.8,1$ |
| Inter connectivity | $0.2,0.4,0.6,0.8,1$ |
| Ridge regularization coeff. | 1e-4, 1e-6, 1e-8, 1e-10 |
| Warmup steps | 0 0 |

#### H.2.2 Hardware

The machine used for fitting the ESN models was 13th Gen Intel(R) Core(TM) i5-1335U @ 4.6 GHz (16GB RAM, 12 cores).

### H.3 shPLRNN

For an explanation of shPLRNN, see Appendix D.

#### H.3.1 Hyperparameters

We used the clipped shPLRNN trained by GTF. For the Van der Pol, Lorenz and the weather datasets we considered a fixed GTF parameter $\alpha$, while for the Rössler we considered an adaptive $\alpha$ (starting from an upper bound) as proposed by Hess et al. []. The code repository by Hess et al. [] was used to perform the computations.
The hyperparameters selected for all datasets are detailed in Table 5. Any hyperparameters not specified are set to their default values in the corresponding repository of Hess et al. []. For Rössler dataset, finding optimal hyperparameters was more challenging, and for training, we also utilized regularizations for the latent and observation models.

**Table 5: Hyperparameters of shPLRNN trained by GTF**
|  | Van der Pol 4.1 | Lorenz 4.3 | Rössler4.3 | Weather |
| --- | --- | --- | --- | --- |
| Hidden dimension | $35$ | $100$ | $50$ | $200$ |
| Batch Size | $32$ | $30$ | $50$ | $32$ |
| Sequence length | $37$ | $100$ | $150$ | $40$ |
| Teacher forcing interval | $25$ | $13$ | $25$ | $20$ |
| Epochs | $1300$ | $2000$ | $2000$ | $2500$ |
| GTF parameter $\alpha$ | $0.98$ | $0.3$ | $0.9$ (Upper bound) | $0.75$ |
| Latent model regularization rate | - | - | 1e-6 | - |
| Observation model regularization rate | - | - | 1e-4 | - |

#### H.3.2 Hardware

The hardware we used to iteratively train the clipped shPLRNN models includes an 11th Gen Intel(R) Core(TM) i7-11800H CPU @ 2.30GHz and 64.0 GB of RAM (63.7 GB usable).

### H.4 LSTM

We use Tensorflow [Abadi et al., ] to train a baseline LSTM for the weather data problem of predicting the temperature. An LSTM is trained similarly to the process by TensorFlow [].

After training the weather models, we compute the predictions with the specified horizon and concatenate them into a single time series. More specifically, if the horizon is set to one day with a time delay of one week, the prediction on a dataset split goes as follows. First, we use ground-truth values of days one to seven to predict the value for day eight. Then, we use the ground-truth values of days two to eight to predict the value for day nine. This process is repeated until we reach the end of the split. Then, the resulting predictions are concatenated and compared to the ground-truth measurements. This prediction process alligns with the one used for sampled RNNs to ensure a consistent comparison.

#### H.4.1 Hyperparameters

**Table 6: Hyperparameters used for LSTM for Section 4.5.**
|  | LSTM |
| --- | --- |
| Width/units | 64 |
| Learning rate | 5e-5 |
| Max epochs | 30 |
| Patience | 5 |

See Table 7 for details on the hyperparameter grid search for Section 4.5 performed for the sampled RNN and LSTM models.

**Table 7: Hyperparameters used in the grid search for training sampled RNN and LSTM for the temperature prediction in Section 4.5.**
|  | sampled RNN | LSTM |
| --- | --- | --- |
| Width/units | $32,64,128$ | $32,64,128,256,512$ |
| Regularization rate | 1e-10, 1e-8, 1e-6, 1e-4, 1e-2 | — |
| Learning rate | — | 1e-5, 5e-5, 1e-4, 5e-4, 1e-3, 5e-3 |

#### H.4.2 Hardware

For the experiments with the weather data in Section 4.5, we used a machine with AMD EPYC 7402 @ 2.80GHz (256GB RAM, 24 cores) and RTX 3080 Turbo (10GB VRAM, CUDA 12.0).

## Peer Reviews

confidence: 5
contribution: 2
initial_rating: 5
presentation: 2
questions: N/A
soundness: 4
strengths: - Interesting connection to Koopman operator.
- Interesting topic of trying to circumvent gradient based training.
summary: -  The paper proposes a training method for modeling recurrent neural networks without the use of gradient-based methods, such as backpropagation through time (BPTT), which suffers from exploding and vanishing gradients, mostly occurring in a system with chaotic dynamics. Building on concepts from random feature models, such as reservoir computing (echo-state networks), the paper proposes the random sampling of weights (W­­­ and b) in the RNN from a data-driven distribution. In addition the paper employs Koopman operator theory to find the outer weights of the RNN model, which map the current state to the next state. The Koopman operator theory maps the finite nonlinear transformation matrices (outer weights) to a linear infinite-dimensional space, where the extended dynamic mode decomposition (EDMD) method is used to find a finite-dimensional approximation of the Koopman operator. 
-  For model validation, they show some computational experiments comprising simple ODEs, such as the Van Der Pol Oscillator, chaotic dynamics (Lorenz and Rossler systems), and real-world examples involving weather data.
- Paper reports results from these computational experiments in the form of training time and error (MSE/KL Divergence). When compared to other models, such as an LSTM, ESN (echo-state network), and shPLRNN (state of the art backpropagation-based RNN), the proposed model (Sampled-RNN) achieves comparable performance, in terms of MSE and KL Divergence, and a faster training time.
weaknesses: - Sampling procedure not clearly explained. (E.g. 'As we stick to networks with one hidden layer in this paper, we ignore the multilayer sampling here and direct the reader to Bolager et al. (2023) for the full sample and construction procedure for an arbitrary number of hidden layers.')
- Paper presentation generally not clear. Hard to follow. Illustrative example: equation 1 is referred to before presented. Paper requires excessive ‘detective work’ to understand what they did and/or are talking about.
- How sampling approach differs from ESNs seems unclear, and a minor innovation at best.
- Although interesting, connection to Koopman operator theory does not seem novel.

---

confidence: 3
contribution: 3
initial_rating: 5
presentation: 3
questions: 1. “For completeness we have added sigma_hx as an arbitrary activation function. We choose to set sigma_hx as the identity function to let us solve for the last linear layer… other activation functions such as the logit are possible as well”. However this is not shown in the paper. Can the authors provide clarification on the effect of non-linear activations here?
2. In the weather task, the sampled RNN has a similar MSE to shPLRNNs and LSTMs for 1-day forecasting, but higher than both alternative models for week-long forecasting. However, MSE is higher than ESNs for the Rosseler system task. Can the authors provide some clarification on this?
4. Can the authors provide clarifications/results on how performance for chaotic systems might change based on how many samples are drawn for the sampled RNN (and the influence on dimensionality of the hidden layer)? 
5. The authors state: “The complexity of solving this system depends cubically on the minimum number of neurons and the number of data points (respectively, time steps). This means if both the network and the number of data points grow together, the computational time and memory demands for training grow too quickly. For BPTT, the memory requirements are mostly because many gradients must be stored for one update pass”. — can the authors provide a comparison of # of neurons vs # of time steps?
6. Can the authors provide details on how they infer the “data-dependent probability distribution” (especially for the weather data)? 
5. Additionally, can the authors expand on how weights and biases are sampled from the specific, data-dependent probability distribution (again for weather data) and how other parameters are computed using a sequence of linear equations?
soundness: 3
strengths: 1. Proof seems to be sound and theoretically correct given the assumptions stated by the authors, however i am not an very well versed with koopman operator theory and thus have provided a lower confidence score.
2. Their method (called sampled RNN) takes significantly less time to train than alternative state-of-the-art models such as ESNs, shPLRNNs, and LSTMs. Predictions from this model capture patterns in toy experiments as well as real-world data (like weather forecasting), and outperform current models (especially on problems with long horizons).
summary: This paper presents an alternative strategy to backpropagation through time (BPTT) by using a combination of Koopman operator theory and random feature networks instead of gradient-based techniques. This novel method avoids vanishing and exploding gradient problems, and outperforms BPTT in terms of training time and accuracy in a series of empirical comparisons, including time series, forecasting, control, and weather problems.
weaknesses: 1. There is little motivation for the problem in the introduction, and the paper jumps right into formalisms. 
2. The weights and biases need to be sampled from a data-dependent probability distribution, however it's unclear how feasible this is?
3. This method does not converge for controlled systems.

---

confidence: 4
contribution: 2
initial_rating: 6
presentation: 3
questions: The impact of exploding and vanishing gradients is not made clear in this manuscript. I understand their model outperforms LSTMs and other, gradient based models based on numerical values. Could the authors make the impact of EVGP more apparent?

Can the authors include more real-world datasets? As of now, all but 1 result is simulated, and additional examples would significantly strengthen the manuscript.
soundness: 3
strengths: Training recurrent neural networks is notoriously hard. I am partial to the authors approach of fitting an RNN to a a dynamical system by a smart sampling of hidden parameters and Koopman operator based modeling. Particularly, I am drawn to the use of Koopman operator in the context of RNNs. 

The paper is very well written and is a nice read, and the results comparing their gradient-free approach to trained models are impressive.
summary: The authors propose a novel way to construct recurrent neural networks. Their approach is two stepped, where they first generates hidden weights and biases according to a data-dependent distribution, and then construct read-out parameters by approximating the dimensional Koopman operator with dynamic mode decomposition.
weaknesses: It is unclear how much the Koopman operator and EDMD components contributed to model performance. Since the hidden weights initialisation schema is an application of previous work (Bolager et al. (2023)), I see the Koopman related work as the main conceptual innovation. However, based on the presented results it is hard to disentangle where improvements come from and I am not fully convinced of the added value of EDMD. I suggest the authors include two additional experiments: a. setting hidden weights randomly and learning read-out with EDMD and b. setting hidden weights based on Bolager et al. (2023) and learning read-out without EDMD. 

I am also not keen on the name and think it over-promises. Being able to train general recurrent neural networks without gradient descent or BPTT is an extremely ambitious goal, which this paper does not fulfils. As the authors explain in the (very much appreciated) limitation section, their approach does not immediately extend to RNN tasks relating to computer vision or NLP, thus the paper results are mostly regard dimensional dynamical system. I still believe their results are impressive, but the language, and specifically the title, needs to be toned down.

---

confidence: 4
contribution: 3
initial_rating: 8
presentation: 3
questions: The authors compare against ESNs on some problems and LSTM on others.  Is there any reason for these choices?

Can the authors provide more intuition about the weight initialization, i.e., Eq. 4 and the equation for the distribution p_H?
soundness: 4
strengths: The proposed approach fits simple nonlinear systems with high accuracy and very little computational time compared with using BPTT in RNNs.  This makes the approach a very appealing alternative for modeling and controlling such systems.  To my knowledge, the approach is novel (but see below), and open up a new perpective on random networks.
summary: The MS explores an alternative method for fitting recurrent neural networks, based on Koopman operator theory.  A single-layer neural network with *random weights* is used to map the state (and possibly a control input) to a higher-dimensional space, in which the dynamics are assumed to be linear.  Since the state is assumed fully observed, the low-D state at consecutive time steps can be mapped to the high-D state, and the state matrix fit with the normal equations (and similarly for the input matrix when there is a control input).  The map back to the low-dimensional state and the output matrix can be found likewise.  The intuition from Koopman operator theory is that there exists a set, possibly infinite-dimensional, of measurement functions that evolve linearly in time.  The authors shows that this method yields computationally cheap and highly accurate models of simple nonlinear systems.
weaknesses: From the point of view of implementation, the manuscript's random RNNs are a fairly minor variation on echo-state networks.  Indeed, there is a large literature on ESNs and reservoir computing, and I am surprised that this variation has not previously been explored.  Can the authors confirm this?

The introduction of Koopman operator theory does provide a nice intuition about why a linear dynamical system should exist in a higher-dimensional space.  But it seems that the trade-off it enables---linearity for higher dimension---limits the scope of application of this approach: As the authors note, the matrix inversion operation (in the normal equations) is expensive, and it is cubic precisely the dimension of the (large) measurement space.

My understanding of the theory is that, in general, not only is this space infinite-dimensional, but also there is no way to bound the number of required dimensions.  That is, M could be arbitrarily large.  Perhaps this is addressed in the proof in Appendix B, which I did not read closely.  Can the authors provide more insight here?

Finally, I found the paper somewhat hard to read.  This could be my fault, but there seem to be notational issues.  For example, near l. 139 the authors write h_t = F(h_{t-1}), and then later in the same paragraph, h_t' = F(h_t).  Is the idea that, in the first instance, h_t is the model state; where in the second instance, h_t is the *true* state (and therefore F(h_t) need not be equal to h_{t+1})?

Meta-review: In this work, the authors propose a gradient-free approach to training recurrent neural networks in low-dimensions. The proposed approach applies EDMD for output weight training to the sampling RNN framework of Bolager et al. (2023), bypassing backpropagation, and negating the exploding-vanishing gradient problem. Reviewers agreed that the manuscript is well-written, and that the approach is sound and yields good results. 

The more negative responses pointed out that the primary contribution was the incorporation of EDMD for the training procedure. The use of sampling strongly resembles ESNs, which can also be considered "gradient-free" in the same way as the proposed method. From this viewpoint, the new contribution pertains to a difference in architecture in the final layer, which requires the EDMD approach. This is not how the method is portrayed, however, as has been pointed out by reviewers. I disagree with this as a minor issue, or easily addressed by further experiments, as the exposition should be upfront about this, and the title not so broad to compensate. It is also strange that ESNs were not compared to in all settings. The authors commented on this in the discussion period, claiming that ESNs are best for chaotic systems, but if so, why not include them? Finally, I am surprised that reviewers did not critique the baselines, as I would not consider shPLRNNs, LSTMs or ESNs to be state of the art.

I am uncomfortable with recommending the paper for acceptance as-is, and would instead suggest that the authors think more deeply about the relationships between their proposed method and other existing approaches, including ESNs and the actual state-of-the-art, perhaps including other RNN approaches inspired by Koopman operators and DMD.

Decision: Reject