Gradient-free training of recurrent neural networks

We would like to thank the reviewer for the review, and we believe certain points in the review have greatly strengthened the current manuscript. Below, we make a few comments on the weaknesses pointed out, as well as answer the reviewer's questions.

[W1] Indeed, the connection between Koopman operator approximation methods and reservoir computing (RC) / echo-state networks (ESNs) has not been explored extensively. We are only aware of the two papers we cite in the section on related work (Bollt (2021); Gulina & Mauroy (2020)). Bollt (2021) only briefly mentions the connection of RC to DMD (i.e. without nonlinear observables), while Gulina & Mauroy (2020) use the hidden state of the reservoir as a dictionary to approximate the Koopman operator - but do not use the Koopman operator for prediction in combination with the reservoir, i.e., as part of the (linear) flow map, as we propose it here.

[W2] In terms of the inverse operation, the complexity of the operation is indeed the number of neurons $M$ cubed. Which is the same complexity for computing the $C$ matrix. The difference between the two are $O(2M^3)$, instead of $O(M^3)$. This certainly adds some additional computational time, but in our experiments, we observe that these differences are negligible. In particular, they are not noticeable when compared to backpropagation.

[W3] This is correct; it still does not provide the required number of neurons and is a hyperparameter that needs tuning. In the experiments, we observe that the required neurons to get good results are in the hundreds (which also means the addition of $K$ is not too computationally heavy, cf. previous weakness). In terms of theory, one has results with certain rates and probabilistic bounds for feedforward neural networks, but the issue is in adding this to the Koopman theory. We are still working on this, and we hope to understand better in the future.

[W4] We thank the reviewer for pointing this out. There is indeed a mistake in the manuscript. When discussing one arbitrary time step, we consistently write $h_t$ to $F(h_t)=h\_{t+1}$. However, for the dataset, we use $h_n$ and $F(h_n) = h_n'$ to denote the nth data point and the resulting output after applying the evolution operator. This is because denoting the data as $$ h_1, ..., h_T $, h\_{t+1} = F(h_t)$, implies that the data is from one trajectory, which is not the case. The data consists of $N$ points $H$ and $N$ points $H'$ after applying the evolution operator once to $H$.

The mistake was on line 143 and has been updated. Whenever we discuss the specific points in our dataset, we refer to $h_n$ and $h_n’$; for arbitrary points, we denote $F(h_t)=h\_{t+1}$.

[Q1] These decisions were driven mainly by the modeling trends for certain tasks. For example, we consider ESNs to be prominent when it comes to modeling chaotic systems, while we consider LSTMs to be more typical for real-world multivariate time series forecasting with partial observation. In addition, LSTM is a known benchmark for this particular data.

[Q2] Following equation 4, we construct the weight as the direction between the two points we sampled and then put the bias on the line between the original points. To set the bias between two points means that after the weight is applied, part of the points on the line between the two points are less than the bias, and the rest are greater than. This means we can capture the direction of the underlying manifold/input space and set the hyperplane (if using ReLU) or the zero line (tanh) in interesting areas of the manifold. And if the manifold/input space changes drastically, the weight space also adapts.

For the distribution, we construct it by considering the input space and the output space. It emphasizes a pair of points close in the input space that is far away in the output space, such as where the gradient of the true function is high. This means the weights and biases are put in interesting areas of the manifold, where the function changes rapidly.

If the data is not in a supervised setting, as with the controlled input $x$, we sample uniformly among the possible pair of points. This still considers the data through equation 4, but not to the same extent.

We would like to thank the reviewer for a great review, that has indeed improved the updated paper in our opinion. Below you will find a few remarks to the weaknesses you pointed out, as well as answer to your questions.

[W1] We agree with the reviewer that such results should be added. We are still running experiments for this, and hope to return very soon with some updated results.

[W2] This is a fair criticism, and we certainly are not intending on claiming more than we can show for; hence the limitation section as well. We have reviewed the title, and decided upon adding "for low dimensional data", which is what we demonstrate at the current moment. We still find the contribution very valuable, and look forward to tackle high dimensional tasks such as vision and NLP. We also added an additional nod to the limitation in the end of the introduction.

[Q1] Training RNNs with gradient descent algorithms is notoriously hard, especially if the data are governed by long and complicated dependencies. In particular, classical RNNs trained to model chaotic systems, regardless of their architecture, will inevitably suffer from the exploding gradient problem $1$ . Consequently, gradients grow exponentially during backpropagation, leading to extremely large weight updates. Thus, the model's parameters can diverge, causing training to become unstable, or the model may fail to converge, or it may be very difficult to find suitable hyperparameters, and so on. Moreover, if gradients vanishes during backpropagation, the weight updates become negligible, making it difficult for the model to learn from the data, or the model struggles to learn long-term dependencies, as the influence of earlier time steps diminishes quickly, or the model may fail to capture complex patterns in the data, leading to underfitting. Therefore, EVGP results in RNNs becoming unable to learn the underlying temporal patterns.

$1$ Chaos paper: J. Mikhaeil, Z. Monfared, D. Durstewitz, On the difficulty of learning chaotic dynamics with RNNs, NeurIPS, (2022).

[Q2] We have trained a sampled RNN with an additional real-world dataset, and added this to the appendix. Specifically we use a dataset containing household electric power consumption by UC Irvine (https://archive.ics.uci.edu/dataset/235/individual+household+electric+power+consumption), which was used in other publications working with real-world timeseries data (https://arxiv.org/pdf/1509.05142, https://proceedings.neurips.cc/paper_files/paper/2017/file/1b5230e3ea6d7123847ad55a1e06fffd-Paper.pdf).
Further details and results are given in appendix I. We observe good performance for our method on this dataset as well.

Dear authors,

As I have stated in my original review and first response, in my view the main innovation of the work is using EDMD to train RNN output weights, as the initialization for the input and hidden weights is a direct application of Bolager et al. (2023). In my initial review, I requested evidence demonstrating EDMD's specific contribution to RNN training effectiveness compared to standard approaches. The authors' response did not provide experimental results addressing this core concern.

Therefore, I recommend rejection and am lowering my score to a 3 from 5.

Thank you for the insightful review, which we believe have strengthened the updated paper. Below we make some thoughts and remarks on the weaknesses you point out, as well as answers to your questions.

[W1] Our thinking behind the introduction is to quickly point out some known problems with training RNNs before giving a quick intuition of our final model. We then enter a more extensive related work section, contextualizing and motivating the problems we quickly mention in the introduction. However, we agree with the reviewer that the change of pace is perhaps a bit too rapid, and we have rewritten part of the introduction to give some more motivation and a less formal introduction to the method we will describe in the paper.

[W2] We use two different data-dependent probability distributions, briefly described around line 182. They are chosen based on whether we work in a supervised setting. As the latter distribution is simply a uniform distribution over the input data, this is always possible to derive. This is still a data-dependent probability, as it still captures the underlying properties of the input space. That is, two different input spaces create two different weight/bias spaces, meaning the distribution we sample from differs.

[W3] It is true it is harder to show convergence for controlled systems, which we briefly discuss in Appendix B.2. This is also known as a problem in the Koopman community. However, this is when we show convergence given the specific way we estimate K, which means we are also saying something about the learning algorithm. This is even stronger as a result, as it holds for arbitrary finite horizon prediction, without saying we need arbitrary long trajectories as data. If one merely talks about existence and does not consider the algorithm nor the data, then known results of universal approximation of RNNs can be applied. That being said, we want to understand better convergence for controlled systems, also theoretically, and not just demonstrate it experimentally as we do in the paper.

[Q1] If the state space we map to differs from the Euclidean space, e.g., $[0,1]^d$, then we may choose $\sigma_hx$ to map from the Euclidean space and down to $[0,1]^d$. For the specific example, the final approximation is changed from linear to logistic regression. This means we need an iterative method such as gradient descent, but we only need to compute gradients for the final layer; hence, there is no backpropagation, which is still faster than BPTT. Of course, if one can formulate the approximation of the outer layer as a linear regression problem, as we do in the rest of the paper, it is preferred because it is much faster and gradient-free.

[Q2] There is often a trade-off between accuracy and training time, and our goal is to minimize this and even be better when possible.

More concretely, for chaotic systems, we would like to clarify that we did not use MSE as a measure of error for chaotic systems like the Rössler system; instead, we utilized EKL (empirical Kullback-Leibler divergence). Note that no single evaluation metric may offer a complete picture for assessing model performance for complex or chaotic data. As noted in $1$ , one may consider multiple evaluation metrics to evaluate model performance, particularly in chaotic or complex data. This paper provides several metrics, such as D_stsp, D_H, and PE, to evaluate model performance. In some cases, it is probable that while one or two metrics favor a particular method, others may indicate higher errors (for instance, see Table 1 in $1$ ). With this in mind, in the case of the Rössler example, we acknowledge that the EKL for our method is slightly higher than that of ESNs. However, we would argue that the EKL between the models is fairly similar, and we can conclude that the performance is comparable, as well as our method being faster than the other methods.

One reason the LSTM and shPLRNN produce better MSE for the weather data could lie in their predictions’ (in)stability. As shown in Figure 6 of the paper, the LSTM model doesn't capture fluctuations in the data, while shPLRNN and sampled models do. We can also see that fluctuations of the sampled model have a higher frequency than those of shPLRNN. Predicting higher frequencies is a good strategy when forecasting a short horizon, but it can lead to overfitting on a longer horizon. Thus, we can possibly see these results as underfitting from the LSTM and overfitting from the sampled model in the case of week-long forecasting.

$1$ GTF paper: F. Hess, Z. Monfared, M. Brenner, D. Durstewitz, Generalized Teacher Forcing for Learning Chaotic Dynamics, ICML, (2023).

To answer all your questions adequately, we continue below and hope the reviewer finds this acceptable. Thanks!

We thank the reviewer for the review which we believe have improved our new version of the manuscript. Below we add a few comments on their weakness points:

[W1] We chose to refer to the original paper for multilayer as we are not using more than one hidden layer for each network we sample. We still provide the details for how a pair of points form the weights and bias and how we sample them by describing two probability distributions on line 183. We still recognize that it certainly is not easy to briefly understand the complete algorithm.

[W2] Thank you for pointing this out. We agree that, in particular, the introduction can be hard to follow. To make it clearer and avoid jumping straight into formalism, we have rewritten part of the introduction to alleviate the issue.

[W3] The sampling approach approximates the Koopman map, which would be a random matrix for an ESN. This comes with three benefits: 1) increased interpretability and better tools to analyze the network --- both theoretically but also through spectral properties of the approximated Koopman matrix, and 2) significantly decrease in the hyperparameter search space, 3) the possibility of including linear control tools such as LQR (see forced Van der Pol in paper).

Furthermore, we consider the data-dependent sampling distribution to be a key element in our method, which can reduce the number of necessary hidden units, as well as less hyperparameter tuning, and thus a shorter training time compared to ESNs. We confirmed this in numerical experiments.

[W4] There is a lot of research on numerical approximation of the Koopman operator, including papers on random Fourier features used as dictionary in EDMD (as cited in the related work).

However, there is very little research that connects recurrent neural networks with Koopman operator theory. Our paper provides one of the first theoretical results in this direction - and we are not aware of any paper that directly ties recurrent network weights to Koopman operator approximations.

We also refer to our response to reviewer 4 here, because it provides related information:
The connection between Koopman operator approximation methods and reservoir computing (RC) / echo-state networks (ESNs) has not been explored extensively. We are only aware of the two papers we cite in the section on related work (Bollt (2021); Gulina & Mauroy (2020)). Bollt (2021) only briefly mentions the connection of RC to DMD (i.e. without nonlinear observables), while Gulina & Mauroy (2020) use the hidden state of the reservoir as a dictionary to approximate the Koopman operator - but do not use the Koopman operator for prediction in combination with the reservoir, i.e., as part of the (linear) flow map, as we propose it here.