# A Meta-learning Approach to Reservoir Computing: Time Series Prediction with
Limited Data

Daniel Canaday  Potomac Research, LLC, Alexandria, VA 22314, USA  Andrew
Pomerance  Potomac Research, LLC, Alexandria, VA 22314, USA  Michelle Girvan
Department of Physics, University of Maryland, College Park, Maryland, 20742,
USA

###### Abstract

Recent research has established the effectiveness of machine learning for
data-driven prediction of the future evolution of unknown dynamical systems,
including chaotic systems. However, these approaches require large amounts of
measured time series data from the process to be predicted. When only limited
data is available, forecasters are forced to impose significant model
structure that may or may not accurately represent the process of interest. In
this work, we present a Meta-learning Approach to Reservoir Computing (MARC),
a data-driven approach to automatically extract an appropriate model structure
from experimentally observed “related” processes that can be used to vastly
reduce the amount of data required to successfully train a predictive model.
We demonstrate our approach on a simple benchmark problem, where it beats the
state of the art meta-learning techniques, as well as a challenging chaotic
problem.

##  1 Introduction

Time series prediction is fundamentally and practically important across a
wide variety of domains, and approaches to this task vary wildly in complexity
and sophistication. Some examples include parameter estimation based on
physical models [1], filtering techniques [2, 3], autoregressive-moving
average models and their variants [4, 5, 6], and other machine learning models
such as support vector machines [7], deep learning [8, 9], and recurrent
neural networks [10, 11, 12].

Generally speaking, there is a trade-off inherent to these techniques: either
significant structure needs to be imposed by a physical or mathematical model,
or large amounts of training data is required for machine learning-based
methods. We are interested in this work in predicting time series with
extremely limited data and poor theoretical models. “Extremely limited” can be
defined as the case where not enough information is available to reduce our
uncertainty about the true dynamics that generate the time series; this,
naturally, forces us to impose some structure. Since our theoretical
understanding is poor, we would like to automate an inductive process that
learns the appropriate structure from other, related time series. This is an
example of knowledge transference, which is concerned with developing
techniques to efficiently adapt solutions to similar problems and has received
considerable interest in a variety of contexts [13, 14, 15, 16], including
time series prediction [17, 18, 19]. This problem and potential solutions are
often referred to as meta-learning or “learning about learning” [20, 21, 22,
23].

In this work, we introduce Meta-learning Approach to Reservoir Computing
(MARC), a novel meta-learning technique whose goal is to enable accurate time
series forecasting with vastly reduced training data. Instead of imposing
domain-specific model constraints, MARC learns appropriate model structure
from a library of related examples in an automated fashion, illustrated
schematically in Fig. 1. There are two distinct learning phases. In the first
learning phase, a recurrent neural network model known as a reservoir computer
(RC) is trained separately on each available time series in the library,
yielding an associated predictive model. While the RC model is an extremely
flexible time series modelling framework, a relatively large number
($10^{2}-10^{4}$) of parameters must be learned to get satisfactory results,
which increases the amount of training data required. Although the number of
RC parameters is large, we hypothesize that the trained RC models actually
reside on a manifold in the RC parameter space of much lower dimensionality by
virtue of the relatedness of the library processes. The second learning phase,
therefore, uses a feedforward deep autoencoder neural network [24, 25, 26] to
project the low-dimensional manifold of learned RC feature vectors onto a
relatively small set of latent representation parameters that encode the
predictive models. After the latent representation is learned, we can search
this space to identify candidate RC feature vectors that explain new, incoming
data. By constraining our attention to this reduced manifold of relevant RC
features, we drastically reduce data requirements while avoiding overfitting
because we are searching within the space of reasonable RC models rather than
the space of all RC models.

An illustrative example that will be revisited in Sec. 4.1 is that of modeling
a sine wave from only a few points on the curve. While this is a simple
problem, it has been previously used by others for algorithm validation and
thus allows quantitative comparison of our algorithm to state-of-the-art few-
shot meta-learning algorithms [27, 15]. In this toy regression problem, there
are two parameters that define the time series–the amplitude and phase. If we
know that the signal to be predicted is sinusoidal, very few noiseless
observations are required. However, if we were to use a standard machine
learning approach without any prior knowledge of the signal’s character, we
would need significantly more observations to generate a usable predictive
model. This is, in essence, the idea behind MARC: rather than searching for
the best fit among all RC models, we search only on the manifold of RC models
that approximate sinusoidal dynamics. Instead of explicitly fitting to a
sinusoid, MARC accomplishes this in a model-free manner by learning (via the
autoencoder) the space of sinusoidal RCs from other examples of sinusoidal
dynamics.

Figure 1: A schematic illustration of the MARC algorithm. At left, a library
of related time series (possibly from a combination of measurements and
approximate simulated processes) is used to train an associated library of RC
models. These RC features are used to train an autoencoder with a low-
dimensional hidden layer. At right, the decoder portion of the autoencoder is
used to map latent variables onto RC features, which then generates a time
series. The latent variables are optimized such that this generated time
series best agrees with a signal consisting of a limited amount of time series
data to be predicted.

Our proposed approach stands in contrast to previous meta-learning approaches
to time-series prediction problems. Many state-of-the-art meta-learning
algorithms, such as MAML ([15]), HSML ([28]), LEO ([29]) and others, use meta-
knowledge to constrain initial parameters (or sets of initial parameters) from
which a regressive model can be fine-tuned to a new task. By contrast, MARC
uses meta-knowledge (the library) to constrain the complete set of parameters
for a predictive model. This is possible because, as we discuss in the
following sections (see, in particular, Fig. 2), MARC creates a latent
representation of the unknown degrees of freedom that distinguish library
members from one another. Effectively, this allows us to do curve-fitting
without ever accessing the true unknown parameters–rather, we fit within the
latent space, from which the entire distribution of time series is
represented.

The rest of this paper is organized as follows. In Sec. 2, we describe the
scope of the problem space, including defining “related time series,” and
quantifying the data requirements. In Sec. 3, we describe the three parts of
the MARC framework: the RC model for time series prediction, compressing the
resulting RC models with autoencoders, and then searching the latent variable
space for usable predictive models with limited data. In Sec. 4, we
demonstrate the practical application of the proposed algorithm with a number
of simulated systems. Finally, we conclude with observations and
considerations for further research in Sec. 5.

##  2 Definitions and Scope

We assume that knowledge of a class of processes is available in the form of a
library of time series. We define a library L as a set of members
$\textbf{l}^{i}$, where each $\textbf{l}^{i}$ is a set of time-labeled
observations of the state of a dynamical system. That is,
$\textbf{L}=\\{\textbf{l}^{i}\\}$, and, as an example, each $\textbf{l}^{i}$
might be a discrete, evenly-spaced time sampling of a finite duration,
continuous-time, state observation $\textbf{y}^{i}(t)$ from a library process
$i$, where the library process $i$ and the process generating the signal to be
predicted are assumed to be related. Although it is not necessary for our MARC
technique, it is conceptually useful to consider the relatedness of the
library members as modeled as an ensemble of systems where each
$\textbf{y}^{i}(t)$ obeys

$\displaystyle\dot{\textbf{x}}^{i}$ $\displaystyle=$ $\displaystyle
f(\textbf{x}^{i};\textbf{v}^{i})$ $\displaystyle\textbf{y}^{i}$
$\displaystyle=$ $\displaystyle g(\textbf{x}^{i};\textbf{v}^{i})$ (1)

subject to the initial condition $\textbf{x}^{i}(0)=\textbf{x}^{i}_{0}$ and an
unobserved parameter vector $\textbf{v}^{i}$. We call $f$ and $g$ the
dynamical and observation functions, respectively, and $\textbf{x}^{i}$ and
$\textbf{y}^{i}$ the state-space and observation variables, respectively. With
these definitions, each $\textbf{l}^{i}$ is a function of a particular initial
condition $\textbf{x}_{0}^{i}$ and parameter vector $\textbf{v}^{i}$.
Considering $f$ and $g$ to be fixed but unknown, the library members can be
thought of as instances from an ensemble of realizations of a random process
characterized by a probability distribution
$p(\textbf{x}_{0}^{i},\textbf{v}_{0}^{i})$. In what follows, our discussion
will be in the context of relatedness as given by Eq. 2.

We call the signal to be predicted the test signal and denote it
$\textbf{s}(t)$. We assume that $\textbf{s}(t)$ originates from a dynamical
system given by Eq. 2 with initial condition $\textbf{x}_{0}^{s}$ and
parameter vector $\textbf{v}^{s}$. A test signal is then only distinguished
from a library member in that it is much shorter, and its parameter vector
and/or initial condition may be outside of the support of
$p(\textbf{x}_{0}^{i},\textbf{v}_{0}^{i})$, thus requiring any algorithm that
learns about L to extrapolate.

There are several cardinalities that are relevant to these definitions, namely
$\dim(\textbf{x)}$, $\dim(\textbf{y})$, $\dim(\textbf{v})$, $|\textbf{L}|$,
$|\textbf{l}^{i}|$, and $|\textbf{s}|$, which are the dimension of the state
space, the dimension of the observation space, the number of parameters, the
number of library members, the number of observation time-points per library
member, and the number of observations in the test signal, respectively. We
consider an interesting regime to be where $|\textbf{s}|\ll|\textbf{l}^{i}|$.
More exactly, we consider a regime where each $\textbf{l}^{i}$ is sufficiently
long that a useful predictive model can be obtained, but the same is not true
for s, hence the need to incorporate the knowledge from L. On the other hand,
we must also have
$|\textbf{s}|\dim(\textbf{y})>\dim(\textbf{x})+\dim(\textbf{v})$, otherwise
the problem is underdetermined even if L were to provide complete information
about $f$ and $g$.

Although we aim to be concrete and concise in our above definitions of
“relatedness”, we believe the applicability of MARC is quite broad. Reservoir
computing methods are capable of representing a wide range of dynamical
phenomena beyond what can be described by Eq. 1, including noisy dynamics
([30]), nonautonomous systems ([31, 32]), control systems ([33, 34]),
discrete-time maps ([32]), partial differential equations ([12]), delay
differential equations ([35]), and more. Generalizations of the version of
MARC presented here are straightforward and only involve modifying the form of
RC used to represent library members. Moreover, MARC can handle multi-modal
situations in which the library members have different functional forms
(Appendix E).

##  3 The MARC Framework

In this section, we describe the MARC algorithm (summarized in Algorithm 1)
and framework in more complete detail. As discussed above, the key idea is to
1) obtain RC features corresponding to an accurate predictive model for each
library member, 2) train a deep autoencoder to learn a latent representation
of these RC features, and then 3) fit, within latent variable space, a set of
decoded RC features that best explains a small sample of new data.

Algorithm 1 MARC algorithm

1:Library $\textbf{L}=\\{\textbf{l}^{i}\\}$, set of time series related to the
process of interest

2:each $\textbf{l}^{i}\in\textbf{L}$ is long enough for successful RC training

3:$T$, RC training procedure

4:$R$, RC prediction procedure

5:Short test signal s

6:Choose RC training hyperparameters and time scale $c$

7:Randomly assign W, $\textbf{W}_{in}$, and b

8:for all $\textbf{l}^{i}\in\textbf{L}$ do

9: $(\textbf{W}_{out}^{i},\textbf{r}_{0}^{i})\leftarrow
T(\textbf{l}^{i};\textbf{W}$, $\textbf{W}_{in}$, $\textbf{b},c)$

10:end for

11:Let $A=D\circ E$ be the autoencoder

12:initialize autoencoder parameters $\Theta$ and reconstruction error

13:while reconstruction error too large do

14: update parameters $\Theta$

15: error
$\leftarrow|(\textbf{W}_{out}^{i},\textbf{r}_{0}^{i})-A(\textbf{W}_{out}^{i},\textbf{r}_{0}^{i})|^{2}$

16:end while

17:Initialize latent variable e and prediction error

18:while prediction error too large do

19: update e

20: error $\leftarrow|R(D(\textbf{e});\textbf{W}$, $\textbf{W}_{in}$,
$\textbf{b},c)-\textbf{s}|^{2}$

21:end while

22:$(\textbf{W}_{out}^{s},\textbf{r}_{0}^{s})\leftarrow D(\textbf{e})$

23:Generate predictions using $\textbf{W}_{out}^{s},\textbf{r}_{0}^{s}$

###  3.1 Learning Time Series Representations with Reservoir Computing

The first step in the MARC algorithm is to learn independent representations
of each of the available time series in the library. For this purpose, we use
a type of recurrent neural network known as a reservoir computer (RC), a class
of machine learning techniques often applied to processing temporal data [35,
36, 37]. We refer to the appendix for details on the RC model used in this
work and [38] for a comparison to deep learning techniques for time series
processing.

The general RC procedure is as follows. Given an input time series
$\textbf{u}(t)$, a recurrent neural network known as the reservoir is driven
with $\textbf{u}(t)$ and evolves according to

$c\dot{\textbf{r}}(t)=-\textbf{r}(t)+\textbf{tanh}\left(\textbf{W}\textbf{r}(t)+\textbf{W}_{in}\textbf{u}(t)+\textbf{b}\right)$
(2)

where $\textbf{r}(t)$ is the state of the reservoir at time $t$, $c$ is a time
constant, W is a fixed matrix of recurrent connections, $\textbf{W}_{in}$ is a
fixed input matrix, and b is a fixed bias vector. After an appropriate warm-up
time, the state $\textbf{r}(t)$ of the reservoir is recorded and ridge
regression is used to identify a time-independent, linear map such that
$\textbf{W}_{out}\textbf{r}(t)\approx\textbf{o}^{d}(t)$, where
$\textbf{o}^{d}(t)$ is a desired output time series. This procedure, that maps
an input time series u to an initial state $\textbf{r}_{0}$ and output weight
matrix $\textbf{W}_{out}$ is the procedure $T(\textbf{u};\textbf{W}$,
$\textbf{W}_{in}$, $\textbf{b},c)$ in Algorithm 1. The training procedure has
several hyperparameters that are discussed in Appendix A.

For forecasting problems, the desired output is the same as the input signal,
$\textbf{o}^{d}(t)=\textbf{u}(t)$. After the training period, future values of
the input signal are predicted by replacing the input in Eq. 2 with the
predicted output $\textbf{W}_{out}\textbf{r}$. This forms an autonomous
predictive model with dynamics governed by

$\displaystyle\dot{\textbf{r}}=f_{RC}(\textbf{r};\textbf{W}_{out})$
$\displaystyle=$
$\displaystyle-\textbf{r}/c+\tanh\left(\textbf{W}\textbf{r}+\textbf{W}_{in}\textbf{W}_{out}\textbf{r}+\textbf{b}\right)/c$
$\displaystyle\textbf{o}=g_{RC}(\textbf{r};\textbf{W}_{out})$ $\displaystyle=$
$\displaystyle\textbf{W}_{out}\textbf{r}.$ (3)

Integrating these equations yields a prediction for a given $\textbf{r}_{0}$
and $\textbf{W}_{out}$, and is the function
$R(\textbf{W}_{out},\textbf{r}_{0};\textbf{W}$, $\textbf{W}_{in}$,
$\textbf{b})$ used in Algorithm 1. We note that Eq. 3 has the same form as Eq.
2 from Sec. 2, with the state x replaced with the reservoir state r, the
observation y replaced with the output o, the parameter vector v replaced with
the parameter matrix $\textbf{W}_{out}$, and the initial state
$\textbf{x}_{0}$ replaced with the initial reservoir state $\textbf{r}_{0}$.
Thus, it is here that we see that a single, fixed RC with a distribution
$p(\textbf{r}_{0},\textbf{W}_{out})$ of initial conditions and parameters is
capable of representing a distribution of time-series generated from
$p(\textbf{x}_{0},\textbf{v})$. We will exploit this connection in the next
section to learn a latent representation of the underlying library parameters.

###  3.2 Compressing RC Representations with Autoencoders

As described in Sec. 2, a library member $\textbf{l}^{i}$ is determined by a
particular parameter vector $\textbf{v}^{i}$ and initial condition
$\textbf{x}_{0}^{i}$, so that we have
$\textbf{l}^{i}=L(\textbf{x}_{0}^{i},\textbf{v}^{i})$. Similarly, the training
process $T$ identifies a map
$(\textbf{r}_{0}^{i},\textbf{W}_{out}^{i})=T(\textbf{l}^{i})$ from library
members to the associated reservoir initial condition $\textbf{r}_{0}^{i}$ and
output weight matrix $\textbf{W}_{out}^{i}$; this tuple of vectors is the RC
features that define a predictive RC model. In particular, we see that the
composition of dynamics and training yield a relationship between the
parameters and initial conditions of the underlying system to the
corresponding RC features. That is, we have

$(\textbf{r}_{0},\textbf{W}_{out})=T\circ L(\textbf{x}_{0},\textbf{v}).$ (4)

For a typical application, the number of reservoir nodes $N$ can be on the
order of thousands or more. If $m$ is the dimension of the desired outputs,
then the dimension of the RC feature space is $(m+1)N$ and can be extremely
large. However, Eq. 4 reveals that if a set of RCs is trained on a library,
then the RC features that encode dynamics approximating time series in the
class generated by Eq. 2 actually lie on a manifold of much smaller dimension
$\dim(\textbf{x}_{0})+\dim(\textbf{v})$ within the larger space. If this
manifold can be identified, then the problem of training an RC to predict a
new time series that is similar to the library becomes a much lower-
dimensional problem, and therefore can be accomplished with much less data.

Of course, the function $T\circ L$ cannot be known without direct measurements
of $(\textbf{x}_{0},\textbf{v)}$ which we assume to be unknown. However, an
approximate, latent representation of this function can be learned with deep
autoencoders [39, 40]. Mathematically, the autoencoder is learning an identity
function that is the composition of an encoding function
$\textbf{e}=E(\textbf{r}_{0},\textbf{W}_{out})$ and a decoding function

$D(\textbf{e})=(\hat{\textbf{r}}_{0},\hat{\textbf{W}}_{out})\approx(\textbf{r}_{0},\textbf{W}_{out})$
(5)

which maps latent variable vectors e to RC feature space. We refer to the
appendix for details on the autoencoder architecture, training, and
hyperparameter considerations for the numerical experiments presented in this
work.

###  3.3 Optimizing Latent Variable Space

Even though a reservoir in feedback mode defined by Eq. 3 with fixed input and
recurrent matrices is capable of representing a wide range of dynamical
systems depending on $\textbf{W}_{out}$, latent variable vectors decoded by
Eq. 5 yield RC features that correspond to RCs that generate signals similar
to those in the library. To use this restricted space of RCs to predict the
future evolution of the test signal s, we need a way to identify which e
decodes to $(\hat{\textbf{r}}_{0},\hat{\textbf{W}}_{out})$ that best matches
the short observed signal s.

Recall that, as each $(\textbf{x}_{0}^{i},\textbf{v}^{i})$ generates a time-
series $\textbf{l}^{i}$, each $(\textbf{r}_{0},\textbf{W}_{out})$ can be used
to generate a reservoir time series over the same time points on which s is
defined, which we denote by
$\hat{\textbf{s}}=R(\textbf{r}_{0},\textbf{W}_{out})$. If we restrict our
attention to RC features that are decoded latent variable vectors, then we
have $\hat{\textbf{s}}=R\circ D(\textbf{e)}$. The final step for forming our
predictive model of s is to find the RC feature vectors by

$\displaystyle\hat{\textbf{e}}$ $\displaystyle=$
$\displaystyle\text{argmin}_{\textbf{e}}|R\circ D(\textbf{e})-\textbf{s}|^{2}$
(6) $\displaystyle(\hat{\textbf{r}}_{0},\hat{\textbf{W}}_{out})$
$\displaystyle=$ $\displaystyle D(\hat{\textbf{e}})$

Because $\dim(\textbf{e})$ is relatively small, Eq. 6 can be efficiently
minimized with numerical algorithms such as Nelder-Mead [41] or dual annealing
[42].

Finally, we note that the composition of the decoding process (Eq. 5) and the
RC prediction process (Eq. 3) forms a map from latent variable space to a
time-series that is parallel to the process by which library members are
generated from Eq. 1. Thus, if the meta-learning process succeeds, the map
$T\circ D:\textbf{e}\rightarrow\textbf{l}^{i}$ is a latent representation of
the map $L:(\textbf{x}_{0},\textbf{v})\rightarrow\textbf{l}$. We will see
explicitly that this is the case for a toy example in the next section.

##  4 Numerical Results

In this section, we apply MARC to a couple of example problems, both to
demonstrate the efficacy of our approach as well as illustrate some concepts
discussed in the preceding sections. The parameters of the library and the
test signals as well as the RCs, autoencoder, and latent space optimizer can
be found in the Appendix.

###  4.1 Sinusoidal Regression

We first return to the sinusoidal regression problem discussed in the
introduction, examined before in various forms in the context of few-shot
meta-learning. As in [15] and [27], we consider a library of sine waves with
unit frequency and amplitudes and initial phases uniformly sampled from
$U[0.1,5.0]$ and $U[0,\pi]$, respectively. As in the other works, we consider
the test signal to be defined on $-5\leq t\leq 5$. However, we consider the
library signals to be defined on the longer interval $-6\leq t\leq 5$ to
ensure that the RC models are valid on the interval $-5\leq t\leq 5$–see
Appendix for more details on this issue. The library L contains 100 examples,
sampled from the distribution of initial conditions. Each member
$\textbf{l}^{i}$ contains a set of 1,000 equally-spaced observations. Each
test signal s contains 10 randomly-sampled observations of a different sine
wave. Details on the RC and autoencoder training are summarized in the
appendix.

To evaluate the performance of the MARC algorithm, we report the mean squared
error (MSE) over the complete sine wave in Table 1, and visualize typical
examples in Fig. 2. In panels (a) and (b), we show how the true parameters
(panel (a)) map to latent variable space (panel (b)). In panels (c) and (d) we
give two examples: from just 10 irregularly-spaced observations in s, the
results of the MARC algorithm is clearly a sinusoidal wave that is a close
approximation to the true data.

In Table 1, we compare MARC to MAML [15] and MeLA [27], two generic
metalearning algorithm proposals that report results for this task. We find
that MARC outperforms these algorithms by two orders of magnitude, but we note
some caveats to this comparison. In MARC, we use all available data (i.e., the
entire sine wave sample for each library member) to train one RC per library
member one time. On the other hand, MAML and MeLA repeatedly sample only 10
observations from library members in the inner loop of their metalearning
approach. Because of the different usages of the library data, it is hard to
directly compare the total data usage of MARC to MAML/MeLA. Moreover,
MAML/MeLA both use feedforward networks to perform nonlinear regression in
this task, whereas we are using RCs for time series processing and interpret
the sine wave as a time series. As a result, these two approaches were applied
to other deep learning architectures for other tasks, but the scope of this
work is limited to time series prediction. Finally, we note that MARC takes
tens of minutes on a workstation-class laptop without GPU acceleration.

| MSE
---|---
MAML [15] | 0.208
MeLA [27] | 0.129
MARC (this work) | 0.0016
Table 1: A comparison of the MSE in the sinusoidal regression problem for our
algorithm and two state-of-the-art meta-learning approaches, as reported in
[27].  Figure 2: Examination of the MARC algorithm applied to a toy regression
problem. Upper left: 10,000 randomly sampled phase/amplitude points, plotted
in the true 2-d parameter space, colored by location. Upper right: The first
two principal components of latent variable space, where the color of each
point corresponds to its original phase and amplitude. As we can see, the RC
features map onto a smooth, compact manifold within latent variable space.
Bottom: Two examples of the MARC algorithm applied to the sinusoidal
regression problem. The blue curve represents the ground truth, while the
dashed red line is the optimal reservoir time series $\hat{\textbf{s}}$. The
10 points, randomly sampled from the ground truth to form s, are indicated by
black triangles.

###  4.2 Lorenz Prediction

We now consider a more challenging example in the form of the Lorenz-63 [43]
system. The Lorenz-63 system is a three dimensional system of equations with
three parameters; the long time behavior yields the famous “butterfly
attractor” and it is a paradigmatic example of chaos used to benchmark time
series prediction algorithms.

Chaotic systems, due to their fundamental property of exponentially divergent
trajectories, are inherently difficult to predict, and any amount of error in
the predictive model will eventually lead to uncorrelated predictions.
Therefore, to quantify the performance of the MARC algorithm in this section,
we report a valid prediction length $T_{valid}$, defined as the amount of time
before the squared error first exceeds the variance of the target signal
(details in the Appendix). Intuitively, $T_{valid}$ defines the amount of time
that the predictive model remains correlated with the true signal. It is best
understood in units of the Lyapunov time $\Lambda_{L}$: small perturbations
grow as $e^{t/\Lambda_{L}}$. The Lyapunov time of the Lorenz-63 system is
$\Lambda_{L}\approx 1.104$s [44].

For this Lorenz-63 example, we construct a library of 100 example Lorenz
attractors, where the parameters are drawn from a normal distribution around
their canonical values, and each of the library members $\textbf{l}^{i}$ is a
set of 5,000 observations separated by $\Delta t=0.01$, or, in terms of the
Lyapunov time, about $50\Lambda_{L}$. When trained with that much data, the
baseline RCs that go into the autoencoder yield predictions with mean valid
times approximately 2.87 $\Lambda_{L}$. Rather than randomly sampling the
temporal domain, the test signals s consist of 10 data points separated by
$\Delta t=0.01$, totalling approximately $0.1\Lambda_{L}$. When trained with
such limited data data, RCs do a poor job at both short-term forecasting
(remaining valid for only 0.03 $\Lambda_{L}$ and at replicating the long-term
behavior of the system. However, as seen in Fig. 3, the MARC algorithm yields
a predictive model that is valid for more than a Lyapunov time and maintains
the butterfly attractor of the Lorenz system. Over 100 test signals, the MARC
prediction remains valid for an average of $1.06\Lambda_{L}$, more than 10
times longer than the observed portion of the attractor.

Figure 3: An example of the MARC algorithm applied to the Lorenz system. Left:
An example Lorenz attractor (blue) generated by the result of the MARC
approach using the 10 observed data points (red). As can be seen, the decoded
RC yields a textbook example, like decoded RCs in Sec. 4.1 all yield sines.
Right: The three components of the Lorenz system plotted vs. time illustrating
the valid time calculation, both the true system (blue) and the MARC-based
prediction (dashed red). The 10 points, sequentially sampled from the ground
truth to form s, are in green. The prediction remains valid for approximately
1.5 Lyapunov times, or until the dashed, vertical line, after which the
systems become decorrelated due to the presence of chaos.

##  5 Conclusions and Future Work

In this work, we have introduced a novel meta-learning approach to time series
prediction from limited training data from the process of interest. The core
idea behind our MARC method is that by using a flexible, general machine-
learning model with a well-behaved training algorithm, the model features
resulting from training a class of related processes lies on a compact, lower-
dimensional manifold. The manifold can be effectively learned with standard
autoencoders, and the latent variables can be searched to find a model that
matches a small amount of test data. The resulting model can be used for
prediction and avoids overfitting, because there are only a few free
parameters to adjust to new data. We have demonstrated our approach by
applying it to a simple regression example and obtaining results that greatly
exceed state-of-the-art gradient-based techniques. We then apply it to a
nontrivial example involving chaos, demonstrating the wide applicability of
the technique. We have explicitly explored some of the properties of latent
variable space and given a precise notion to the scope of situations in which
we expect our algorithm to apply.

There are several avenues to improve the results we have obtained. First, the
most apparent difficulty with our approach is the need for ultra-high-
precision autoencoders. For the results in this work, thousands of epochs were
trained to obtain an autoencoder with MSE on the order of $10^{-4}$ or
smaller. It is likely that better predictive models can be obtained if this
error were to be even further reduced, due to the high sensitivity of the
predicted signals $\hat{\textbf{s}}$ on errors in $\textbf{W}_{out}$ and
$\textbf{r}_{0}$. On the other hand, the RC algorithm itself might be refined
to be more robust to these errors, such that more easily obtained autoencoder
performance would be satisfactory. This avenue is suggested by some
preliminary analysis on the dependence of the MARC algorithm on the RC
regularization parameter (see Appendix), which suggests that values which are
much higher than optimal for conventional RC actually perform best in this
scenario. This may be because increasing the penalty for variance in
$\textbf{W}_{out}$ makes more error tolerable.

Second, the generic nature of the optimization routine in Sec. 3.2 can be
improved by further study. Nelder-Mead is a Jacobian-free method chosen mostly
for expedience, but given that the error $\hat{\textbf{s}}-\textbf{s}$ is the
function of a recurrent neural network (the RC) and a feedforward neural
network (the autoencoder), we could explicitly compute the Jacobian and use
other optimization schemes to more efficiently identify the correct latent
variable e.

Thirdly, the issue of choosing an appropriate latent variable dimension is an
issue. In the presented model system examples, we knew the proper
dimensionality of the manifold; in practice with real data, this dimension is
unknown. In our numerical experiments, we used a slightly larger latent
variable dimension than strictly necessary to accomodate loss due to the
autoencoder (see Appendix), but a principled approach to select this crucial
hyperparameter is needed before widespread practical application. Choosing a
dimension that is too large greatly increases the burden on the optimizer so
that the probability of getting stuck in a local optimium increases, while
choosing a dimension that is too small leads to a lack of expressiveness in
the RC model that may miss important dynamics.

Notwithstanding these issues, we are optimistic that the ideas behind MARC can
have far-reaching application beyond the focus of time series prediction
presented in this work. Firstly, it may be the case that RCs trained for other
tasks (such as control [45, 46], signal classification [47, 48], and anomaly
detection [49, 50]) or through reinforcement learning [51] can be similarly
compressed with autoencoders so that new, related systems can be controlled or
classified with greatly reduced training data requirements. Secondly, we have
focused here on RCs–both because of our initial focus on time series
prediction and due to attractive numerical properties–but the general scheme
of Sec. 3 may be applied to other to machine learning models, such as deep
network architectures, so that other tasks such as computer vision may benefit
from this approach. Finally, the ideas behind MARC may be a tool in scientific
discovery: systems with different dynamics cluster in latent variable space, a
behavior which may be useful in exploratory data analysis. Furthermore, in
experimental cases where data collection is hard and/or expensive, identifying
regions of latent variable space that are sparsely populated may be useful to
guide future data collection. Finally, RCs have been proposed as a
computationally-efficient method to replace integration of stiff ODEs [52,
53]; combining these ideas with a metalearning parameter interpolation scheme
could greatly reduce the computational requirements of scientific computing.

## Acknowledgements

This work was supported by DARPA under contract W31P4Q-20-C-0077. We thank
Edward Ott, Brian Hunt, and Jiangying Zhou for useful comments and discussion.
The views, opinions, and/or findings expressed are those of the author(s) and
should not be interpreted as representing the official views or policies of
the Department of Defense or the U.S. Government.

## References

  * [1] Ottar N Bjørnstad, Bärbel F Finkenstädt, and Bryan T Grenfell.  Dynamics of measles epidemics: estimating scaling of transmission rates using a time series sir model.  Ecological monographs, 72(2):169–184, 2002.
  * [2] Howard Musoff and Paul Zarchan.  Fundamentals of Kalman filtering: a practical approach.  American Institute of Aeronautics and Astronautics, 2009.
  * [3] Tae Woo Joo and Seoung Bum Kim.  Time series forecasting based on wavelet filtering.  Expert Systems with Applications, 42(8):3868–3874, 2015.
  * [4] S Sp Pappas, L Ekonomou, D Ch Karamousantas, GE Chatzarakis, SK Katsikas, and P Liatsis.  Electricity demand loads modeling using autoregressive moving average (arma) models.  Energy, 33(9):1353–1360, 2008.
  * [5] Scott H Holan, Robert Lund, Ginger Davis, et al.  The arma alphabet soup: A tour of arma model variants.  Statistics Surveys, 4:232–274, 2010.
  * [6] Durdu Ömer Faruk.  A hybrid neural network and arima model for water quality time series prediction.  Engineering applications of artificial intelligence, 23(4):586–594, 2010.
  * [7] Kyoung-jae Kim.  Financial time series forecasting using support vector machines.  Neurocomputing, 55(1-2):307–319, 2003.
  * [8] Hongju Yan and Hongbing Ouyang.  Financial time series prediction based on deep learning.  Wireless Personal Communications, 102(2):683–700, 2018.
  * [9] Zhongyang Han, Jun Zhao, Henry Leung, King Fai Ma, and Wei Wang.  A review of deep learning models for time series prediction.  IEEE Sensors Journal, 2019.
  * [10] Min Han, Jianhui Xi, Shiguo Xu, and Fu-Liang Yin.  Prediction of chaotic time series based on the recurrent predictor neural network.  IEEE transactions on signal processing, 52(12):3409–3416, 2004\.
  * [11] Felix A Gers, Douglas Eck, and Jürgen Schmidhuber.  Applying lstm to time series predictable through time-window approaches.  In Neural Nets WIRN Vietri-01, pages 193–200. Springer, 2002.
  * [12] Jaideep Pathak, Brian Hunt, Michelle Girvan, Zhixin Lu, and Edward Ott.  Model-free prediction of large spatiotemporally chaotic systems from data: A reservoir computing approach.  Physical review letters, 120(2):024102, 2018.
  * [13] Christiane Lemke, Marcin Budka, and Bogdan Gabrys.  Metalearning: a survey of trends and technologies.  Artificial intelligence review, 44(1):117–130, 2015.
  * [14] Matthias Feurer, Jost Springenberg, and Frank Hutter.  Initializing bayesian hyperparameter optimization via meta-learning.  In Proceedings of the AAAI Conference on Artificial Intelligence, volume 29, 2015.
  * [15] Chelsea Finn, Pieter Abbeel, and Sergey Levine.  Model-agnostic meta-learning for fast adaptation of deep networks.  In International Conference on Machine Learning, pages 1126–1135. PMLR, 2017.
  * [16] Chelsea Finn, Aravind Rajeswaran, Sham Kakade, and Sergey Levine.  Online meta-learning.  In International Conference on Machine Learning, pages 1920–1930. PMLR, 2019.
  * [17] Boris N Oreshkin, Dmitri Carpov, Nicolas Chapados, and Yoshua Bengio.  Meta-learning framework with applications to zero-shot time-series forecasting.  arXiv preprint arXiv:2002.02887, 2020.
  * [18] Thiyanga S Talagala, Rob J Hyndman, George Athanasopoulos, et al.  Meta-learning how to forecast time series.  Monash Econometrics and Business Statistics Working Papers, 6:18, 2018.
  * [19] Christiane Lemke and Bogdan Gabrys.  Meta-learning for time series forecasting and forecast combination.  Neurocomputing, 73(10-12):2006–2016, 2010.
  * [20] Ivan Olier, Noureddin Sadawi, G Richard Bickerton, Joaquin Vanschoren, Crina Grosan, Larisa Soldatova, and Ross D King.  Meta-qsar: a large-scale application of meta-learning to drug design and discovery.  Machine Learning, 107(1):285–311, 2018.
  * [21] Brian Gaudet, Richard Linares, and Roberto Furfaro.  Terminal adaptive guidance via reinforcement meta-learning: Applications to autonomous asteroid close-proximity operations.  Acta Astronautica, 171:1–13, 2020.
  * [22] Ahmed Abbasi, Conan Albrecht, Anthony Vance, and James Hansen.  Metafraud: a meta-learning framework for detecting financial fraud.  Mis Quarterly, pages 1293–1327, 2012.
  * [23] Jorge Kanda, Carlos Soares, Eduardo Hruschka, and Andre De Carvalho.  A meta-learning approach to select meta-heuristics for the traveling salesman problem using mlp-based label ranking.  In International Conference on Neural Information Processing, pages 488–495. Springer, 2012.
  * [24] Yann LeCun, Yoshua Bengio, and Geoffrey Hinton.  Deep learning.  nature, 521(7553):436–444, 2015.
  * [25] Diederik P Kingma and Max Welling.  Auto-encoding variational bayes.  arXiv preprint arXiv:1312.6114, 2013.
  * [26] Michael Tschannen, Olivier Bachem, and Mario Lucic.  Recent advances in autoencoder-based representation learning.  arXiv preprint arXiv:1812.05069, 2018.
  * [27] Tailin Wu, John Peurifoy, Isaac L Chuang, and Max Tegmark.  Meta-learning autoencoders for few-shot prediction.  arXiv preprint arXiv:1807.09912, 2018.
  * [28] Huaxiu Yao, Ying Wei, Junzhou Huang, and Zhenhui Li.  Hierarchically structured meta-learning.  In International Conference on Machine Learning, pages 7045–7054. PMLR, 2019.
  * [29] Andrei A Rusu, Dushyant Rao, Jakub Sygnowski, Oriol Vinyals, Razvan Pascanu, Simon Osindero, and Raia Hadsell.  Meta-learning with latent embedding optimization.  arXiv preprint arXiv:1807.05960, 2018.
  * [30] Roland S Zimmermann and Ulrich Parlitz.  Observing spatio-temporal dynamics of excitable media using reservoir computing.  Chaos: An Interdisciplinary Journal of Nonlinear Science, 28(4):043118, 2018.
  * [31] H Jaeger.  Advances in neural information processing systems.  In Advances in Neural Information Processing Systems, volume 15, pages 593–600. MIT Press, 2003.
  * [32] Dhruvit Patel, Daniel Canaday, Michelle Girvan, Andrew Pomerance, and Edward Ott.  Using machine learning to predict statistical properties of non-stationary dynamical processes: System climate, regime transitions, and the effect of stochasticity.  Chaos: An Interdisciplinary Journal of Nonlinear Science, 31(3):033149, 2021.
  * [33] Matthias Salmen and Paul G Ploger.  Echo state networks used for motor control.  In Proceedings of the 2005 IEEE international conference on robotics and automation, pages 1953–1958. IEEE, 2005.
  * [34] Jaemann Park, Bongju Lee, Seonhyeok Kang, Pan Young Kim, and H Jin Kim.  Online learning control of hydraulic excavators based on echo-state networks.  IEEE Transactions on Automation Science and Engineering, 14(1):249–259, 2016.
  * [35] Herbert Jaeger.  The “echo state” approach to analysing and training recurrent neural networks-with an erratum note.  Bonn, Germany: German National Research Center for Information Technology GMD Technical Report, 148(34):13, 2001.
  * [36] Wolfgang Maass, Thomas Natschläger, and Henry Markram.  Real-time computing without stable states: A new framework for neural computation based on perturbations.  Neural computation, 14(11):2531–2560, 2002.
  * [37] David Verstraeten, Benjamin Schrauwen, Michiel d’Haene, and Dirk Stroobandt.  An experimental unification of reservoir computing methods.  Neural networks, 20(3):391–403, 2007.
  * [38] Pantelis R Vlachas, Jaideep Pathak, Brian R Hunt, Themistoklis P Sapsis, Michelle Girvan, Edward Ott, and Petros Koumoutsakos.  Forecasting of spatio-temporal chaotic dynamics with recurrent neural networks: A comparative study of reservoir computing and backpropagation algorithms.  arXiv preprint arXiv:1910.05266, 2019.
  * [39] Pierre Baldi.  Autoencoders, unsupervised learning, and deep architectures.  In Proceedings of ICML workshop on unsupervised and transfer learning, pages 37–49. JMLR Workshop and Conference Proceedings, 2012.
  * [40] Geoffrey E Hinton and Ruslan R Salakhutdinov.  Reducing the dimensionality of data with neural networks.  science, 313(5786):504–507, 2006.
  * [41] John A Nelder and Roger Mead.  A simplex method for function minimization.  The computer journal, 7(4):308–313, 1965.
  * [42] Y Xiang, DY Sun, W Fan, and XG Gong.  Generalized simulated annealing algorithm and its application to the thomson model.  Physics Letters A, 233(3):216–220, 1997.
  * [43] Edward N Lorenz.  Deterministic nonperiodic flow.  Journal of atmospheric sciences, 20(2):130–141, 1963.
  * [44] Divakar Viswanath.  Lyapunov exponents from random fibonacci sequences to the lorenzequations.  Technical report, Cornell University, 1998.
  * [45] Daniel M Canaday, Andrew Pomerance, and Daniel J Gauthier.  Model-free control of dynamical systems with deep reservoir computing.  Journal of Physics: Complexity, 2021.
  * [46] Eric Aislan Antonelo and Benjamin Schrauwen.  On learning navigation behaviors for small mobile robots with reservoir computing architectures.  IEEE transactions on neural networks and learning systems, 26(4):763–780, 2014.
  * [47] Miguel Angel Escalona-Morán, Miguel C Soriano, Ingo Fischer, and Claudio R Mirasso.  Electrocardiogram classification using reservoir computing with logistic regression.  IEEE Journal of Biomedical and health Informatics, 19(3):892–898, 2014.
  * [48] Gouhei Tanaka, Ryosho Nakane, Toshiyuki Yamane, Seiji Takeda, Daiju Nakano, Shigeru Nakagawa, and Akira Hirose.  Waveform classification by memristive reservoir computing.  In International Conference on Neural Information Processing, pages 457–465. Springer, 2017.
  * [49] Oliver Obst, X Rosalind Wang, and Mikhail Prokopenko.  Using echo state networks for anomaly detection in underground coal mines.  In 2008 International Conference on Information Processing in Sensor Networks (ipsn 2008), pages 219–229. IEEE, 2008.
  * [50] Qing Chen, Anguo Zhang, Tingwen Huang, Qianping He, and Yongduan Song.  Imbalanced dataset-based echo state networks for anomaly detection.  Neural Computing and Applications, 32(8):3685–3694, 2020.
  * [51] Hanten Chang and Katsuya Futagami.  Reinforcement learning with convolutional reservoir computing.  Applied Intelligence, 50(8):2400–2410, Aug 2020.
  * [52] Marios Mattheakis, Hayden Joy, and Pavlos Protopapas.  Unsupervised reservoir computing for solving ordinary differential equations.  arXiv preprint arXiv:2108.11417, 2021.
  * [53] Ranjan Anantharaman, Yingbo Ma, Shashi Gowda, Chris Laughman, Viral Shah, Alan Edelman, and Chris Rackauckas.  Accelerating simulation of stiff nonlinear systems using continuous-time echo state networks.  arXiv preprint arXiv:2010.04004, 2020.
  * [54] Mantas Lukoševičius, Herbert Jaeger, and Benjamin Schrauwen.  Reservoir computing trends.  KI-Künstliche Intelligenz, 26(4):365–371, 2012.
  * [55] Andrei Nikolaevich Tikhonov, AV Goncharsky, VV Stepanov, and Anatoly G Yagola.  Numerical methods for the solution of ill-posed problems, volume 328.  Springer Science & Business Media, 2013.
  * [56] Trishul Chilimbi, Yutaka Suzue, Johnson Apacible, and Karthik Kalyanaraman.  Project adam: Building an efficient and scalable deep learning training system.  In 11th $\\{$USENIX$\\}$ Symposium on Operating Systems Design and Implementation ($\\{$OSDI$\\}$ 14), pages 571–582, 2014.
  * [57] Rich Caruana, Steve Lawrence, and Lee Giles.  Overfitting in neural nets: Backpropagation, conjugate gradient, and early stopping.  Advances in neural information processing systems, pages 402–408, 2001.

##  Appendix A Reservoir Computing

In this appendix, we provide additional details for our implementation of RC
in the numerical examples in this work, including hyperparameter
considerations for the specific experiments. Recall from Sec. 3.1 that the
first step in the MARC algorithm is to learn independent predictive models for
each available time series. Thus, we require a flexible machine learning
framework that is capable of representing a range of dynamical systems.
Further, we require the training algorithm to produce model features that are
smooth, deterministic functions of the training data such that Eq. 4 holds. A
commonly-used type of RC known as an echo-state network (ESN) ([35]) satisfies
these requirements and is used in this work, though we note that other machine
learning frameworks may be suitable as well.

###  A.1 Echo State Networks

In the ESN formulation of RC, the reservoir is a recurrent neural network
model whose state $\textbf{r}(t)$ is a column vector with
$\textbf{r}^{T}(t)=[r_{1}(t),r_{2}(t),...,r_{N}(t)]$ where the scalar state of
each network node $l$ is $r_{l}(t)$ and the reservoir dynamics are defined by
Eq. 2. In the equations, $c$ is a time constant, W is a matrix of recurrent
connections, $\textbf{W}_{in}$ is an input-connectivity matrix, b is a bias
vector, and $\textbf{W}_{out}$ is the trained output matrix. Importantly, we
emphasize that the network dynamics, i.e., $c$, W, $\textbf{W}_{in}$, and b
are held constant during training. Note also that the trainable weights
$\textbf{W}_{out}$ do not impact the internal reservoir dynamics and,
therefore, can be trained after observing the reservoir response
$\textbf{r}(t)$ to the input signal $\textbf{u}(t)$.

There are a number of different possibilities for constructing the the
network-defining matrices W, $\textbf{W}_{in}$, and b. In this work, we do so
according to the following procedures. We choose elements of the connectivity
matrix W to be nonzero with probability $k/N$, where $N$ is the number of
network nodes and $k$ is the mean in-degree of the network. Nonzero elements
are chosen from a normal distribution with zero mean and an arbitrary
variance. The elements of W are then rescaled to fix the spectral radius to
$\rho_{\textbf{W}}$. This scaling step is important to ensure that the
dynamics of $\textbf{r}(t)$ “forget” its initial condition and previous inputs
and to help match the time scale of the RC dynamics to the input signal [35].
The matrix $\textbf{W}_{in}$ and the vector b and have elements drawn from a
symmetric uniform distribution with maximum value $\sigma_{in}$ and
$\sigma_{\textbf{b}}$, respectively. The former hyperparameter determines the
coupling strength of the input to the reservoir, while the latter promotes a
diversity of individual node dynamics.

Generally speaking, $c$ describes the dominant time-scale of the reservoir
dynamics. Since the RC is used in this work to model a time series, $c$ should
be chosen to roughly match the dominant time-scale $\textbf{u}(t)$. Because
$\sigma_{in}$ scales the input signal, it is generally taken to be
approximately the inverse of the scale (standard deviation or maximum absolute
value) of the input signal. Finally, $\sigma_{\textbf{b}}$ is chosen between 0
and 1.

###  A.2 Training

For a review of RC methods, we refer the reader to [54]. As outlined in the
reference, RC training is most often achieved by linear regression with
Tikhonov regularization [55]. Let R be the matrix of discrete observations of
the reservoir state space. That is, the $n^{th}$ column of R is
$\textbf{r}(T_{init}+n\Delta t)$, where $T_{init}$ is the warm-up time,
$\Delta t$ is the sampling period, and $n$ is an integer. Similarly, let
$\textbf{O}^{d}$ be a matrix whose $n^{th}$ column is the desired output at
time $(T_{init}+n\Delta t)$. Then the output matrix $\textbf{W}_{out}$ is
chosen to minimize the loss given by

$|\textbf{W}_{out}\textbf{R}-\textbf{O}^{d}|^{2}+\alpha|\textbf{W}_{out}|^{2},$
(7)

where $|\textbf{M}|^{2}$ for a rectangular matrix M is defined as the trace of
$\textbf{MM}^{T}$ and $\alpha$ is a small, positive regularization parameter
that prevents overfitting. The regularization $\alpha$ is commonly selected
with cross-validation techniques. However, we find that the MARC algorithm
benefits from a significantly larger value of $\alpha$ than what is optimal
for using RC to predict a single time series. (This may be because less
sensitive $\textbf{W}_{out}$ values produce a manifold that is smoother and
easier for the autoencoder to learn.)

Conveniently, the minimum to Eq. 7 can be expressed in the closed form given
by

$\textbf{W}_{out}=(\textbf{R}^{T}\textbf{R}+\alpha\textbf{I})^{-1}\textbf{R}^{T}\textbf{O}^{d}.$
(8)

As the right-hand-side of Eq. 8 is the composition of continuous functions,
this makes explicit the claim in the main text that the trained RC features
are continuous functions of the training data.

The RC hyperparameters for the sine-regression and Lorenz-63 experiments are
presented in Tab. 2 and 3, respectively.

##  Appendix B RC Compression with Autoencoders

To obtain a mapping from a low-dimensional latent space to the continuous set
of RC features, we compress the features through an autoencoder with a low-
dimensional intermediate layer. Thus, the autoencoder $A$ is the composition
of an encoding function $E$ and a decoding function $D$, only the latter of
which is used for fitting RC features to a new test signal. A schematic of
this architecture is presented in Fig. 6.

For the examples in this paper, we construct the autoencoder with five hidden
layers, the middle one being the smallest size and identified as e. More
specifically, the autoencoder contains an input layer of size $(m+1)N$,
followed by two hidden layer, followed by the latent layer of size
$dim(\textbf{e})$, followed by two more hidden layer and finally an output
layer of size $(m+1)N$. The first six layers have sigmoid activation
functions, while the last layer is linear. All of the layers are dense and
have no intra-layer connections. The network is trained using the ADAM
algorithm [56] with learning rate fixed at 0.0001. We employ an early-stopping
procedure [57] to avoid over-fitting, where 10% of the training samples are
reserved for validation, and training ends when the validation error worsens
or does not improve for a number of epochs. For the sine-regression
experiment, the successive hidden layers are of dimension 200, 200, 4, 200,
200, respectively. For the Lorenz-63 experiment, these dimensions are 600,
200, 7, 200, 600, respectively.

Figure 4: The autoencoder architecture. An autoencoder takes
$(\textbf{r}_{0},\textbf{W}_{out})$ as input and produces an approximation to
$(\textbf{r}_{0},\textbf{W}_{out})$ as its output. The decoder layer
approximates $(\textbf{r}_{0},\textbf{W}_{out})$ from a small-dimensional
latent space variable e, mimicking the process of generating
$(\textbf{r}_{0},\textbf{W}_{out})$ from $(\textbf{x}_{0},\textbf{v})$.

##  Appendix C Optimizing in Latent Variable Space

Recall that, as each $(\textbf{x}_{0}^{i},\textbf{v}^{i})$ generates a time-
series $l^{i}$, each $(\textbf{r}_{0},\textbf{W}_{out})$ can be used to
generate a reservoir time series over the same time points on which $s$ is
defined, which we denote by
$\hat{\textbf{s}}=R(\textbf{r}_{0},\textbf{W}_{out})$. If we restrict our
attention to RC features that are decoded latent variable vectors, then we
have $\hat{\textbf{s}}=R\circ D(\textbf{e)}$. The final step for forming our
predictive model of s is to find the RC feature vectors by minimizing the loss
in Eq. 6.

Because $dim(\textbf{e})$ is small, Eq. 6 can be efficiently minimized with a
variety of numerical algorithms. For this purpose, we use an efficient global
search algorithm based on simulated dual annealing [42]. For alle xperiments
in this paper, we implement this algorithm with a maximum iteration limit of
100 and a local search based on the Nelder-Mead method [41].

library, test parameter | value |  RC parameter | value |  autoencoder parameter | value
---|---|---|---|---|---
$|\textbf{L}|$ | 100 | $N$ | 1,000 | $dim(\textbf{h}_{1})$ | 200
$|\textbf{l}^{i}|$ | 1,000 | $k$ | 100 | $dim(\textbf{h}_{2})$ | 200
$|\textbf{s}|$ | 10 | $\rho_{\textbf{W}}$ | 1.0 | $dim(\textbf{e})$ | 4
|  | $\sigma_{in}$ | 1.0 |  |
|  | $\sigma_{\textbf{b}}$ | 1.5 |  |
|  | $c$ | $1.0$ |  |
|  | $\alpha$ | $5e^{-6}$ |  |
|  | $T_{init}$ | 1.0 |  |
Table 2: The hyperparameters for the MARC algorithm used for the sinusoidal regression experiments described in Sec. 4.1. The same hyperparameters (except for $dim(\textbf{e})$) are used for the multi-modal experiments described in Appendix E. library, test parameter | value |  RC parameter | value |  autoencoder parameter | value
---|---|---|---|---|---
$|\textbf{L}|$ | 1,000 | $N$ | 1,000 | $dim(\textbf{h}_{1})$ | 600
$|\textbf{l}^{i}|$ | 5,000 | $k$ | 100 | $dim(\textbf{h}_{2})$ | 200
$|\textbf{s}|$ | 10 | $\rho_{\textbf{W}}$ | 0.8 | $dim(\textbf{e})$ | 7
|  | $\sigma_{in}$ | 0.05 |  |
|  | $\sigma_{\textbf{b}}$ | 0.5 |  |
|  | $c$ | $1.0$ |  |
|  | $\alpha$ | $5e^{-4}$ |  |
|  | $T_{init}$ | 10.0 |  |
Table 3: The hyperparameters for the MARC algorithm used for the Lorenz
prediction experiments described in Sec. 4.2.

##  Appendix D Details on the Lorenz Experiment

In Sec. 4.2, we discussed the challenging problem of forecasting the Lorenz-63
system from limited data. The Lorenz-63 system is described by the following
dynamical equations:

$\displaystyle\dot{x}_{1}$ $\displaystyle=$ $\displaystyle v_{1}(x_{2}-x_{1})$
$\displaystyle\dot{x}_{2}$ $\displaystyle=$ $\displaystyle
x_{1}(v_{2}-x_{3})-x_{2}$ (9) $\displaystyle\dot{x}_{3}$ $\displaystyle=$
$\displaystyle x_{2}x_{2}-v_{3}x_{3}$ $\displaystyle\textbf{y}=\textbf{x}.$
(10)

The parameters $v_{1}$, $v_{2}$, and $v_{3}$ are sampled from the uniform
distributions $U[10,15]$, $U[28,42]$, and $U[8/3,12/3]$, respectively. Initial
conditions are randomly sampled from the attractor by allowing the system to
propagate for $100\Lambda_{L}$ before recording the time series. The
hyperparameters used for the RCs for the numerical results in Sec. 4.2 are
presented in Table 3.

##  Appendix E Applicability to Multi-modal Tasks

In this section, we discuss the applicability of the MARC algorithm to multi-
modal tasks where the library L consists of different types of functions. In
particular, we consider the experiment from [28] where each $\textbf{l}^{i}$
and s are generated, with equal probability from the following four functions:

$\displaystyle y(x)=A\text{sin}(\omega x+b)$ (11) $\displaystyle
y(x)=A_{l}x+B_{l}$ (12) $\displaystyle y(x)=A_{q}x^{2}+B_{q}x+C_{q}$ (13)
$\displaystyle y(x)=A_{c}x^{3}+B_{c}x^{2}+C_{c}x+D_{c}$ (14)

, where the coefficients are drawn from the distributions $A\leftarrow
U[0.1,0.5]$, $\omega\leftarrow U[0.8,1.2]$, $b\leftarrow U[0,2\pi]$,
$A_{l}\leftarrow U[0.1,0.5]$, $B_{l}\leftarrow U[0.1,0.5]$, $A_{q}\leftarrow
U[0.1,0.5]$, $B_{q}\leftarrow U[0.1,0.5]$, $C_{q}\leftarrow U[0.1,0.5]$,
$A_{c}\leftarrow U[0.1,0.5]$, $B_{c}\leftarrow U[0.1,0.5]$, $C_{c}\leftarrow
U[0.1,0.5]$, and $D_{c}\leftarrow U[0.1,0.5]$.

We train the MARC algorithm on a library of size $\textbf{L}=1,000$. The
latent dimension $dim(\textbf{e})$ is increased to 7 to account for the
increased degrees of freedom in the library, but hyperparameters are otherwise
as in the sine regression experiments. Across 800 test signals, we find the
mean MSE (0.394) to be comparable to results reported in [28], being exceeded
only by methods that are designed with multi-modal distributions in mind.

Upon closer inspection of MARC results on this toy problem, we see that the
mean error is dominated by a few outliers among the many test signals. In
fact, the median MSE is several orders of magnitude better at 0.006, and that
errors appear to be normally distributed in logarithmic space, as illustrated
in Fig. 5. Further depicted in Fig. 5 are the median and worst-case
performances on the test signals. We see here that the median performance is
quite good, while worst-case scenarios are driven by predictions that
extrapolate poorly and appear to lose the the sinusoidal or polynomial
character of the library members. We also see that, generally, the cubic test
signals are most difficult for MARC to learn simultaneously with the
sinusoidal, linear, and quadratic examples.

We further investigate the source of poor predictions in Fig. xx where we plot
the relationship between the optimizer loss (the last step of MARC) and the
final MSE. It is clear that poor performance is predominantly a result of not
finding an adequate global minimum solution to Eq. 6.

The above observations, together with the fact that the first two steps (RC
training and autoencoding) of MARC perform exceptionally well on this problem
(as measured by autonomous RC prediction error and autoencoding error,
respectively), suggest that the difficulty of applying MARC to the multi-modal
problem arises from the optimization step. This is perhaps due to an increase
in local minima to Eq. 6. This example further motivates a more careful
consideration of latent-space optimization that is tuned to behave well with
the MARC algorithm. It also suggests possible refinements to MARC, such as
adding a term to the loss in Eq. that penalizes deviation from the library-
created manifold, or a hierarchical structure to the latent space in the
spirit of HSML. These ideas are the subject of ongoing research.

Figure 5: Illustrative examples of the MARC algorithm applied to the multi-
modal example. The first, second, third, and fourth rows correspond to the
sine, line, quadratic, and cubic test signals, respectively. The first
(second) column depicts the test with the median (worst) MSE out of 800 total
experiments. Black triangles indicate the observed points in s. Blue lines
indicate the ground truth, while dashed red lines indicate the MARC
prediction. The third column contains the distribution of errors in
logarithmic space. Figure 6: Optimizer loss versus final MSE for the multi-
modal experiments. Black, blue, green, and red dots correspond to sines,
lines, quadratic functions, and cubic functions, respectively.