High Dimensional Hybrid Reservoir Computing

Predicting two-dimensional spatiotemporal chaotic patterns with

optimized high-dimensional hybrid reservoir computing
T. Nakano,1 S. Baur,2, 3 and C. Räth2, 3
1) Deutsches Zentrum für Luft- und Raumfahrt (DLR), Institut für KI Sicherheit, Rathausallee 12, 53757 Sankt Augustin,
Germany
2) Deutsches Zentrum für Luft- und Raumfahrt (DLR), Institut für KI Sicherheit, Wilhelm-Runge-Str. 10, 89081 Ulm,
Germany
3) Department of Physics, Ludwig-Maximilians-Universität, Schellingstraße 4, 80799 Munich,
Germany
(*Electronic mail: [Link]@[Link])
(Dated: 7 January 2025)
As an alternative approach for predicting complex dynamical systems where physics-based models are no longer re-
liable, reservoir computing (RC) has gained popularity. The hybrid approach is considered an interesting option for
arXiv:2501.02369v1 [[Link]] 4 Jan 2025

improving the prediction performance of RC. The idea is to combine a knowledge-based model (KBM) to support the
fully data-driven RC prediction. There are three types of hybridization for RC, namely full hybrid (FH), input hybrid
(IH) and output hybrid (OH), where it was shown that the latter one is superior in terms of the accuracy and the ro-
bustness for the prediction of low-dimensional chaotic systems. Here, we extend the formalism to the prediction of
spatiotemporal patterns in two dimensions. To overcome the curse of dimensionality for this very high-dimensional
case we employ the local states ansatz, where only a few locally adjacent time series are utilized for the RC-based
prediction. Using simulation data from the Barkley model describing chaotic electrical wave propagation in cardiac
tissue, we outline the formalism of high-dimensional hybrid RC and assess the performance of the different hybridiza-
tion schemes. We find that all three methods (FH, IH and OH) perform better than reservoir only, where improvements
are small when the model is very inaccurate. For small model errors and small reservoirs FH and OH perform nearly
equally well and better than IH. Given the smaller CPU needs for OH and especially the better interpretability of it,
OH is to be favored. For large reservoirs the performance of OH drops below that of FH and IH. Generally, it may
be advisable to test the three setups for a given application and select the best suited one that optimizes between the
counteracting factors of prediction performance and CPU needs.

The prediction of high-dimensional spatiotemporal pat- tional needs4–6 . RC is by default a fully data-driven approach
terns emerging from nonlinear complex systems like and is expected to learn the dynamics of the underlying
turbulent flows, excitable media and earth systems, is an system from the data. The hybrid approach7 is considered as
essential task in vatious fields of science. Especially in an interesting option to improve the prediction. The idea is to
physics and engineering approximate or reduced order combine a knowledge-based model (KBM, e.g. an imperfect
models of the underlying dynamical phenomenon are governing equation) as a support to the fully data-driven
available in many cases. While these models are still too prediction by the reservoir. The hybrid approach can even be
inaccurate to be used for precise predictions, they may pursued if no governing equations are available. In this case
be combined with fully data-driven, AI-based methods one can derive proxies of the governing equations from the
to allow for precise predictions. In this paper we extend data with e.g. SINDy8 or a causality analysis9 and continue
these hybridization techniques to high-dimensional, with hybrid RC. Also this approach leads to improvements
spatially extended dynamical systems using reservoir in the prediction performance especially when the reservoir
computing (RC) as AI-method. We thoroughly test parameters are not optimized10 .
different setups for hybrid RC and discuss the respective The combination of data-driven and model-based elements
performance gain. Our results suggest that the prediction in hybrid reservoir computing can be done at the input or
of spatiotemporal patterns is significantly improved by output layer of the RC or both of them. The three setups
hybrid RC allowing for a whole new range of conceivable are respectively called, input hybrid (IH), output hybrid
applications. (OH) and full-hybrid (FH). Some studies have already been
performed, for example, on the input-hybrid by Shahi et
al.11 , full-hybrid by Pathak et al.7 and on the output-hybrid
by Doan et al.12,13 . Duncan et al.14 analyzed and compared
the different setups in a systematic manner and showed more
I. INTRODUCTION recently the superiority of the output-hybrid setup in terms
of in terms of accuracy, robustness and interpretabiliy of the
results for a set of three-dimensional chaotic model systems.
Reservoir computing (RC)1–3 has gained popularity as an
alternative solution to predict complex dynamical systems
where physics-based models are no more reliable, because The capability of RC to predict extended, high-dimensional
it combines superior forecasting results with little computa- spatiotemporally chaotic systems is also of great interest in
High Dimensional Hybrid Reservoir Computing 2

complex systems research. The predicted systems in the above It describes the dynamics of two coupled variables U(t) and
mentioned studies of hybrid RC have low dimensionalities V (t) that depend on the four coefficients ε, a, b, D. The
ranging form 3 to 6. To treat high-dimensional systems one numerical solution of the two partial differential equations
possible solution is to use a large reservoir (i.e. a large num- has been performed with an 1st order Euler method for the
ber of reservoir nodes), but this could make the training un- time integration of ∆t = 0.01. A grid of 80 × 80 points with
feasible. This is called “curse of dimensionality”, as one can ∆x = 0.1 constantly all over the grid was used. Thus the two-
also observe in many other machine learning methods. Par- dimensional data sets consists of 6400 simultaneously evolv-
litz et al. 15 suggested a parallel prediction approach based on ing time series. A second order spatial discretization was used
local states (LS), where only a few locally adjacent time se- for the approximation of the Laplace operator. The boundary
ries are utilized for the prediction. This procedure is repeated condition was set to the no-flux condition. The coefficients
for all input time series to be predicted. Pathak et al. ap- were set to be D = 0.02, a = 0.75, b = 0.06 and ε = 0.08.
plied this LS-approach to RC16 and showed its efficiency for The source code for the simulation is available19 by courtesy
a one-dimensional system (Kuramoto-Sivashinsky) consisting of Zimmermann and Parlitz. Fig.2 shows the evolved U(t) at
of 512 time series. Zimmermann et al. 17 applied RC with a randomly captured moment in the simulation.
LS to (two-dimensional) spatially extended systems, namely It is known from previous studies20 that adding noise to in-
models describing excitable media (Barkley and the Bueno- put time series data improves the long-term prediction at the
Orovio-Cherry-Fenton model) and performed a prediction be- cost of the short-term prediction quality. Thus a normally dis-
tween the variables of the system. The two variables U(t) and tributed noise with standard deviation of σSD,noise was added
V (t) of the system were cross-predicted, i.e. predicted from for the training and the synchronization of the training. Given
each other. Their prediction can thus rather be considered as a that σSD,input is the standard deviation of the input time series,
mapping or a reconstruction at each time step than a true tem- we set to σSD,noise = ασSD,input , when the noise-ratio α is a
poral prediction in the usual sense. hyper-parameter.
Wikner et al. 18 proposed an approach that combines hy- In the hybrid approach, the ε-model was employed. This
brid RC with local states and called it Combined Hybrid- was originally suggested by Pathak et al.7 . The idea is to
Parallel Prediction (CHyPP). In the study, CHyPP was com- vary a parameter of the governing equations of the system and
pared to the non-hybrid RC and the non-parallel hybrid then use it as the KBM. Here the same numerical integration
RC. CHyPP showed a better prediction quality and required method as for the ground truth (first order Euler method in our
smaller dataset for the training compared to the other two ap- case) was employed. By doing so, an imperfect predictor is
proaches. However, the studied system was one-dimensional artificially created. In our case, we multiplied the coefficient
(again Kuramoto-Sivashinsky) with 1024 dimensions at high- ε in the Barkely model by the factor of (1 + e) which gives
est and only output hybrid was considered. εe = ε(1 + e). e is called a model error. ε in the equation will
There is a high expectation for RC and hybrid RC to be be replaced by εe for the KBM. Although the model error is
capable of predicting spatially extended, much higher dimen- often denoted by ε as the this approach named "ε-model", we
sional systems for which only proxies of a reliable physics- denote here the model error e to avoid the confusion with ε in
based model are available. Atmospheric models and many the equations of the Barkley model.
phenomena in fluid dynamics are examples.
In this paper we extend the formalism of hybrid Reservoir
Computing to the prediction of high-dimensional, spatially III. RESERVOIR COMPUTING
extended two-dimensional data sets by using Local States. We
calculate and validate the predictions with IH, OH and FH and Let u denote the u-dimensional time series used for the
systematically study the results for different sizes of the reser- training and the synchronization. During the prediction, the
voir and varying model mismatches and compare them with previous prediction yr (t − ∆t) is used as the input for the next
the reservoir only results. prediction. Let ũ(t) denote all these inputs time series. In the
case of a classical, non-hybrid RC, a xdim -dimensional input
vector x(t) for the reservoir is given as x(t) = ũ(t). In the
II. BARKLEY MODEL case of hybrid RC the input vector becomes a function of the
input time series, x(t) = finp (ũ(t)). finp depends on the type
As outlined in Zimmermann et al. 17 we use the (cubic) of hybrid RC and is further specified below.
Barkley model describing chaotic electrical wave propagation RC is a type of recurrent neural networks (RNNs). The
in cardiac tissue as two-dimensional model system. It is a non- uniqueness of the RC is that the weights in the reservoir and
linear chaotic system which is described by partial differential in the input layer are fixed during the training. The reservoir
equations with a 2nd order diffusion term, is a sparse random network consisting of rdim nodes with the
average node degree κ. The random connections between the
  nodes and their strengths (weights) are described by the rdim ×
∂U 1 V +b rdim -dimensional adjacency matrix A. The weights are scaled
= D · ∇2U + U(1 −U) U −
∂t ε a to have a spectral radius ρ. The input time series x(t) are
(1)
∂V connected to the nodes of the network by the input matrix Win
= U 3 −V. which is a sparse rdim × xdim dimensional matrix. Concretely,
∂t
High Dimensional Hybrid Reservoir Computing 3

the input data x(t) are used to define the reservoir state r(t +
1) with the one in the previous time step r(t),
r(t + ∆t) = tanh [Ar(t) + Win x(t)] .
A and Win are initialized once before the training and then
kept fixed. In order to break detrimental symmetries in the
reservoir equations21 we introduce a non-linear transforma-
tion of the reservoir state r̃(t),
T  T
r̃ = r, r 2 = r1 , r2 , ..., rdim , r12 , r22 , ..., rdim
2

.
Then the reservoir state r̃(t) is again transformed to hdim - FIG. 2. U(t) given by the simu-
FIG. 1. Schematic illustration of lation with the Euler method at a
dimensional vector h(t) that can depend on the transformed the different hybrid methods con- randomly chosen instant
reservoir response r̃(t) and the input data ũ(t), h(t) = sidered in this study
fout (r̃(t), ũ(t)). Again in the classical RC h(t) is simply
given by h(t) = r̃(t), whereas in the hybrid case the contri-
butions of the KBM are added via a proper definition of fout . The output matrix Wout now has the shape of
The output matrix Wout connects fout to the target output udim × (rdim + Kdim ). Conversely, the input layer is left
yt (t). Wout is determined in the training by a ridge regression, unchanged in the OH method, i.e. x(t) = ũ(t).

The full hybrid (FH) method is simply a method where both

Wout = min [∥Wout h(t) − yt (t)∥ + β ∥Wout ∥] .
Wout the input- and output-hybrid are simultaneously employed.
Fig.1 schematically illustrates the three hybrid RC methods.
β is the regularization term to avoid an over-fitting. Once
In our example the input time series have the dimension of
trained, the output yr can be given by yr = Wout h(t). The
ũ(t) = 80×80×2 per time step, which is the multiplication of
prediction is done progressively using the prediction of the
the number of the grid points and the number of the variables.
previous time step.
Kdim has also 80 × 80 × 2 dimensions.
 
r(t + ∆t) = tanh Ar(t) + Win xpred (t) .
V. LOCAL STATES

IV. HYBRID APPROACH

Performing the training and the prediction on a 80 × 80 × 2
data-set at each time step can be computationally expensive
We employ three different hybrid approaches as outlined in and one can get trapped in the “curse of dimensionality”, as
the study by Duncan et al.14 . The KBM can be represented as previously mentioned. In a preliminary study, we attempted to
a function being applied on an input time series u and produc- perform the training and prediction with one single RC which
ing an imperfect next-step prediction, covers the whole 80 × 80 × 2 domain. This approach requires
a significantly long dataset for the training. The required com-
K(u) ≈ u(t + ∆t). putational memory is also huge. With limited computational
resource and limited data length, we got very poor prediction
In the input hybrid (IH) method, the reservoir input x(t) performance. These results further suggest the use of Local
is given as the concatenation of the input ũ(t) and the KBM States (LS) for the prediction of high-dimensional extended
output K(ũ(t)), systems. The single RC approach won’t be further considered
  in this paper.
ũ(t) LS is a parallel prediction approach in which only the local
x(t) = finp (ũ(t)) = .
K(ũ(t)) behavior of the system is considered by taking into account
This extends the input dimension of a classical (non-hybrid) only neighboring grid points. To predict the variables U(t)
RC, xdim = udim , to xdim = udim + Kdim . Therefore the input and V (t) at the grid point of (i, j), we use the information of
matrix Win is now a rdim × (udim + Kdim ) matrix. At the the neighbouring points in a σ × σ sized square, including
output layer only the reservoir response and its square is used, (i, j) itself. Fig.3 shows the concept. Physically this approach
i.e. h(t) = r̃(t) is justified, when the interaction lengths in the dynamical sys-
tem under study is limited and covered by the local neighbor-
In the output hybrid (OH) method, the KBM prediction hood. σ = 3 has been chosen from the results of a hyper-
K(ũ(t)) is only fed into the output layer, i.e. h(t) = r(t) in parameter study. This gives a σ 2 × 2 dimensional vector as
the classical (non-hybrid) RC is transformed to input, taking both variables U(t) and V (t). In our study, a
  reservoir has been initialized, trained and performed the pre-
r(t) diction at each point in the domain (i, j). All these procedures
h(t) = fout (r(t), ũ(t)) = .
K(ũ(t)) have been done separately for each point on the grid. This
High Dimensional Hybrid Reservoir Computing 4

FIG. 4. Training and prediction sections in the ensemble experiment.

nT = 2, nP = 3 in this study.

FIG. 3. Local States: the training and the prediction are done for the
red point taking only the information of the red point itself and the steps. The prediction quality was measured by the valid time.
orange points.
Table.I summarizes the examined values and the obtained op-
timal values. When a parameter was examined (varied), the
other parameters were fixed to the initial value. In the fol-
means that 80 × 80 independent reservoirs have been used. lowing section, the optimal values were used unless otherwise
described. The adjacency matrix A and the input matrix Win
were randomly initialized only once at the beginning and fixed
VI. VALID TIME
during the entire study.

We evaluated the quality of the prediction by the valid time.

The valid time is the time during which one can have a predic- TABLE I. Hyper-parameter study
tion with an error smaller than a certain acceptable level. The Parameters Examined values Initial value Optimal value
valid time tv is given as the duration before the normalized rdim ∈ {200, 400, 500, 600} 500 400
ρ ∈ {0.1, 0.3, 0.5, 0.6, 0.7, 0.8, 1.0, 1.2, 1.5} 1.0 0.5
time-dependent error e(t), σ ∈ {3, 5, 7} 5 3
α ∈ 10−x , x = {4, 5, 6} 10−4 10−6
∥y(t) − yr (t)∥ β ∈ 10−x , x = {2, 3, 4, 5, 6, 7, 8} 10−6 10−6
e(t) =
⟨∥y(t)∥2 ⟩1/2

exceeds a threshold value emax . In this study emax was set to

0.2. IX. PREDICTION RESULTS

The Barkley model has been simulated for 40,400 steps.

VII. ENSEMBLE EXPERIMENT METHOD
The first Ttransi = 2, 000 steps were discarded to avoid the tran-
sition regime in the simulation. The next Ttrain,sync = 200 steps
The quality of the prediction by a RC is, as many other ma- were used for the synchronization for the training, and the
chine learning methods, dependent on the specific time series next Ttrain = 30, 000 steps were used for the training. Again
sections used for the training and the prediction. Therefore Tpred,sync = 200 steps were used for the synchronization for
it is important to asses the performance of the prediction in a the prediction, and finally Tpred = 8, 000 steps are employed
statistically meaningful way, such as an ensemble experiment. for the prediction. The computation time for the training plus
This method was employed in Figs.7, 8 and 10 in Sec. IX. the prediction is typically around 24 hours with a CPU of the
In total nT trainings have been performed. Each training- 12th Gen Intel(R) Core(TM) i9-12900K. The reservoir dimen-
section iT is composed of NTD + NTS + NT time steps; NTS sion rdim and the training Ttraining and prediction steps Tpred are
steps of train-synchronization, NT of training. The first NTD the major factors for the computation time.
steps of each section are discarded to avoid the influence from Fig.5 shows the reservoir-only prediction (non-hybrid ap-
the last section. For a each training iT , nP times prediction has proach), and the comparison to the ground truth at a randomly
been done. Each prediction section consists of NPD +NPS +NP chosen moment, the instant A (Lyapunov Time tλmax = 1.4)
time steps, where the first NPD time steps are again discarded. and the instant B (tλmax = 2.4). At the instance A, at least
The example in Fig.4 shows the case of nT = 2, nP = 3. qualitatively the prediction captures well the dynamics of the
system. The prediction error is obviously larger at the instant
B than A.
VIII. HYPER-PARAMETER STUDY Fig.6 shows the prediction by the OH and the comparison
to the ground truth at the same randomly chosen moment as
A hyper-parameter study has been conducted to improve before, i.e. at instant A with tλmax = 1.4 and at instant B
the prediction quality for the Barkley model introduced in with tλmax = 2.4. The model error was set to e = 0.1 and
[Link]. The evaluation was done by a single training and a rdim = 400. In both instants, at least qualitatively the predic-
single prediction. The training has been done with 10.000 tion captures again well the dynamics of the system. As ex-
High Dimensional Hybrid Reservoir Computing 5

instant A, Lyapunov Time tλmax = 1.4 instant A, Lyapunov Time tλmax = 1.4

instant B, Lyapunov Time tλmax = 2.4 instant B, Lyapunov Time tλmax = 2.4

FIG. 5. Prediction by the reservoir-only, U(t): (a)-(c), V (t): (d)-(f), FIG. 6. Same as Fig. 5 but prediction by the OH.
ground truth, prediction and error e(t) = |yt (t) − yr (t)|, from left to
right, where yt (t) is the ground truth, yr (t) is the prediction.
cates that the increase of the reservoir dimension rdim makes
the contribution of the KBM smaller in the OH. This is why
pected the prediction error is larger at the instant B. At both the increase of the reservoir dimension rdim in the OH actu-
instances, the prediction error is alleviated compared to the ally lowers the prediction performance. In case of the IH, the
reservoir-only case in Fig.5, which implies that the prediction input matrix Win has rdim × (udim + Kdim ) = rdim × 36 dimen-
was significantly improved with the help of the OH. sions in total which can be split up inot rdim × udim = rdim × 18
In Fig.7-10, the prediction performance in terms of the for the reservoir part rdim and × Kdim = rdim × 18 for the KBM
model error e and the reservoir dimension rdim was investi- part. This means that the increase of the reservoir dimension
gated. The ensemble experiment has been done with nT = rdim does not make the contribution of the KBM smaller. This
3, nP = 6, 18 experiments in total (see [Link] for the defini- is why the IH shows rather better performance with a larger
tion). As in Fig.5 and Fig.6, the data lengths were set to be dimensions. The FH seems to take the advantages of both the
Ttransi = 2, 000, Ttrain,sync = 200, Ttrain = 30, 000, Tpred,sync = IH and the OH, where it is noteworthy that neither at low di-
200, Tpred = 8, 000. The adjacency matrix A and the input mensions nor at high dimensions of the reservoir FH strongly
matrix Win has been initialized only once at the first training- exceeds the prediction performance of OH or IH respectively.
section iT = 1 and have been kept for all the rest of the sec- In Fig.8-10, the valid time of the different hybrid methods
tions. as a function of the model error e is compared, with the reser-
Fig.7 shows the valid time tv of the different hybrid meth- voir dimension rdim = 100, 300 and 500 respectively. All hy-
ods versus the reservoir dimension rdim . The reservoir-only brid methods outperform the reservoir-only at smaller model
case actually performs worse than the KBM-fitted prediction errors e. This is especially true for the case of a small reservoir
at lower dimensions but as the dimension increases its perfor- rdim = 100, where reservoir-only completely fails to make ac-
mance is improved. The IH improves the valid time around 1 curate predictions. As expected, the hybrid methods get more
compared to the reservoir-only case at the tested range of the inaccurate as the model error e increases. For the OH, when
dimension rdim . The OH and the FH perform already well in the model error e is large, the correction by the reservoir part
the case of smaller reservoirs. The prediction performances at the output matrix Wout is not enough and it simply tries
decreases for OH with increasing the dimension rdim of the to get rid of the contribution from the KBM. No contribution
reservoir. As explained in [Link], the output matrix Wout has from the KBM means actually nothing else than a non-hybrid
udim × (rdim + Kdim ) dimensions in total, udim × rdim for the reservoir. This is why at large errors such as e = 10 or 100, we
reservoir part and udim × Kdim for the KBM part. This indi- see that the OH is much closer to the reservoir-only, compared
High Dimensional Hybrid Reservoir Computing 6

to the two other methods. In case of the IH, the large error by
the KBM will be corrected by the entire model, namely the
input matrix Win , the reservoir and the output matrix Wout .
Obviously this effect is larger than Wout . This seems to mit-
igate more the influence of the large error KBM, compare to
the OH. The FH inherits this advantage of the IH and is thus
more robust than the OH.
The elapsed time for a training plus a prediction for each
case in Fig.7 is shown in Fig.11. The computation has been
done by the same CPU previously mentioned (12th Gen In-
tel(R) Core(TM) i9-12900K). Each case has been executed by
a core. No multi-core computing has been used. The worksta-
tion containing the above CPUs were shared with other users.
For this reason, one should take into consideration that the
following elapsed time measurement has not been done in an
ideal environment to measure the computation speed. One
can see a tendency such that the elapsed time is in the order of
Reservoir − only < OH ≈ IH < FH. Without having the ex-
tra computation related to the hybrid method, the Reservoir- FIG. 7. Reservoir dimension rdim vs valid time with e = 0.1. The
only is the fastest, as expected. the OH and the IH do not valid time is estimated by calculating the median of all the nT × nP
show a significant difference. The FH is the slowest due to predictions and displayed with the corresponding standard deviations
(lower/higher quartile) across each dimension. The KBM-fitted pre-
its heavy computation for its inner IH as well as the OH. Also
diction is also displayed for the comparison.
as expected, the reservoir dimension rdim is clearly the major
factor to increase the computation time. These observations
indicate that from the point of view of CPU needs it is a bet-
ter choice to use the OH or the IH with a smaller reservoir
dimension rdim .
From the above observations, one can say that the OH or the
FH of a small reservoir dimension such as rdim = 50, 300, 100
is the best choice for this use-case, when the model error e is
small. On the other hand, the IH or FH with higher dimensions
such as rdim = 400 or more tend to perform better than the
OH. It should also be noticed that a higher reservoir dimension
takes more time for the training and the prediction. Further,
the FH is more time-consuming compared to the IH and OH
due to its more complex and larger structure.
Fig.12 and 13 visualize the ratio of the contribution (Reser-
voir/KBM) in the output matrix Wout in the OH. On the one
hand, we want to check the correct functioning of the OH
setup. On the other hand, we want to investigate how the
contributions from the (data-driven) RC part and the KBM
part are distributed - especially when the model error is in-
creased. The model error was set to e = 0, 0.1, 5, 100, respec- FIG. 8. Model error e vs valid time for the reservoir dimension rdim =
tively. For e = 0, the noise ratio α and the regularization term 100. The valid time is estimated by calculating the median of all the
β in the training data set were set to 0 to check if the KBM nT × nP predictions and displayed with the corresponding standard
part works correctly as designed. The ensemble experiment deviations (lower/higher quartile) across each dimension. The KBM-
has been performed with 20 experiments in total. We followed fitted prediction is also displayed for the comparison.
the method described in [Link] with 5 training-sections and
no prediction-sections, (nT = 5 and nP = 0). The adjacency
matrix A was initialized 4 times randomly (nA = 4) and 5
training sections were performed in each A (nT = 5). During not have the error ε in the KBM. The increase of the reservoir
the 5 consecutive training sections, A was kept constant. With contribution with the increase of the error in the model con-
e = 0, the trained model refers completely to the side of the firms the expected functioning of the OH. It is remarkable that
KBM, as expected. In the prediction of U the contribution of for the Barkley model the contributions from the KBM remain
the reservoir part increases compared to the KBM as the error larger than the one from RC even for very large model errors.
increases. In the prediction of V , this tendency is much less This means that even very inaccurate models do significantly
pronounced and the KBM contribution is still larger than the contribute to the prediction and are thus worth being included
one from the reservoir even at e = 100. This is because V does in a hybrid RC approach.
High Dimensional Hybrid Reservoir Computing 7

(a) e = 0
FIG. 9. Same as Fig. 8 but for the reservoir dimension rdim = 300.

(b) e = 0.1

FIG. 12. Contribution of Reservoir/KBM in Wout for no model errors

FIG. 10. Same as Fig. 8 but for the reservoir dimension rdim = 500.
e = 0 (a) and a model error of e = 0.1 (b). The contribution is esti-
mated by calculating the median of all the nA × nT predictions and
displayed with the corresponding standard deviations (lower/higher
quartile) across each dimension.

X. CONCLUSIONS AND OUTLOOK

We developed a framework for hybrid Reservoir Com-

puting supported by Local States for high dimensional sys-
tems. The predictions by different methods (non-hybrid,
input-hybrid, output-hybrid and full-hybrid) have been con-
ducted with the Barkley model simulating non-linear chaotic
excitable medium 80 × 80 × 2 = 12, 800 dimensions. The
performance of each method were investigated. RC success-
fully predicted the Barkley model with and without the hybrid
approach. The subsequent analyses showed that all hybrid
methods significantly improve the prediction performance.
FIG. 11. Elapsed time for one training and prediction for the ensem- The OH and the FH of a small reservoir dimension such as
ble experiment of Fig.7
rdim = 50, 100 is the best choice for this use case, when the
High Dimensional Hybrid Reservoir Computing 8

simulations in fluid dynamics and atmospheric dynamics as

well as excitable media.

ACKNOWLEDGEMENT

We thank R. Zimmermann and U. Parlitz for providing us

access to the simulation code of the Barkeley model.

DATA AVAILABILITY STATEMENT

The data that support the findings of this study are available
from the corresponding author upon reasonable request.

(a) e = 5
REFERENCES

1 H. 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,
13 (2001).
2 W. Maass, T. Natschläger, and H. Markram, “Real-time computing without

stable states: A new framework for neural computation based on perturba-

tions,” Neural computation 14, 2531–2560 (2002).
3 H. Jaeger and H. Haas, “Harnessing nonlinearity: Predicting chaotic sys-

tems and saving energy in wireless communication,” Science 304, 78–80

(2004).
4 A. Chattopadhyay, P. Hassanzadeh, and D. Subramanian, “Data-driven pre-

dictions of a multiscale lorenz 96 chaotic system using machine-learning

methods: reservoir computing, artificial neural network, and long short-
term memory network,” Nonlinear Processes in Geophysics 27, 373–389
(2020).
5 S. Bompas, B. Georgeot, and D. Guéry-Odelin, “Accuracy of
neural networks for the simulation of chaotic dynamics: Pre-
cision of training data vs precision of the algorithm,” Chaos:
(b) e = 100 An Interdisciplinary Journal of Nonlinear Science 30 (2020),
10.1063/5.0021264, 113118, [Link]
pdf/doi/10.1063/5.0021264/14629142/113118_1_online.pdf.
FIG. 13. Same as Fig. 12 but for large model errors, i.e. e = 5 (a) 6 S. Shahi, F. H. Fenton, and E. M. Cherry, “Prediction of chaotic time se-
and e = 100 (b). ries using recurrent neural networks and reservoir computing techniques: A
comparative study,” Machine Learning with Applications 8, 100300 (2022).
7 J. Pathak, A. Wikner, R. Fussell, S. Chandra, B. R. Hunt, M. Girvan, and

E. Ott, “Hybrid forecasting of chaotic processes: Using machine learning in

model error e is small. Given the smaller CPU needs for OH conjunction with a knowledge-based model,” Chaos: An Interdisciplinary
and especially the better interpretability of it, OH is to be fa- Journal of Nonlinear Science 28 (2018).
8 S. L. Brunton, J. L. Proctor, and J. N. Kutz, “Discovering governing equa-
vored. On the other hand, the IH or FH with higher dimen-
sions such as rdim = 400 or more tend to perform better than tions from data by sparse identification of nonlinear dynamical systems,”
Proceedings of the National Academy of Science 113, 3932–3937 (2016),
the OH. It was also shown that a higher reservoir dimension arXiv:1509.03580 [[Link]].
takes more time for the training and the prediction. The FH 9 H. Ma, A. Haluszczynski, D. Prosperino, and C. Räth, “Identifying causal-
setup is also more time-consuming as compared to the IH and ity drivers and deriving governing equations of nonlinear complex sys-
OH. Generally, it may thus be advisable to test the three se- tems,” Chaos 32, 103128 (2022).
10 F. Köster, D. Patel, A. Wikner, L. Jaurigue, and K. Lüdge, “Data-informed
tups for a given application and select the best suited one that
reservoir computing for efficient time-series prediction,” Chaos 33, 073109
optimizes between counteracting factors of prediction perfor- (2023).
mance and CPU requirements. 11 S. Shahi, C. D. Marcotte, C. J. Herndon, F. H. Fenton, Y. Shiferaw, and
In summary, the combination of locals sates and hybrid RC E. M. Cherry, “Long-time prediction of arrhythmic cardiac action poten-
expands the application potential of hybrid RC to the class of tials using recurrent neural networks and reservoir computing,” Frontiers in
high-dimensional, spatially extended complex systems. Con- physiology 12, 734178 (2021).
12 N. A. K. Doan, W. Polifke, and L. Magri, “Physics-informed echo state
ceivable application that we will consider in the future are networks for chaotic systems forecasting,” in Computational Science–ICCS
time evolutionary non-linear systems that are represented in 2019: 19th International Conference, Faro, Portugal, June 12–14, 2019,
the form of two- or three-dimensional image, such as flow Proceedings, Part IV 19 (Springer, 2019) pp. 192–198.
High Dimensional Hybrid Reservoir Computing 9

13 N. A. K. Doan, W. Polifke, and L. Magri, “Short-and long-term predic- 17 R. S. Zimmermann and U. Parlitz, “Observing spatio-temporal dynamics
tions of chaotic flows and extreme events: a physics-constrained reservoir of excitable media using reservoir computing,” Chaos: An Interdisciplinary
computing approach,” Proceedings of the Royal Society A 477, 20210135 Journal of Nonlinear Science 28 (2018).
(2021). 18 A. Wikner, J. Pathak, B. Hunt, M. Girvan, T. Arcomano, I. Szunyogh,
14 D. Duncan and C. Räth, “Optimizing the combination of data-driven and A. Pomerance, and E. Ott, “Combining machine learning with knowledge-
model-based elements in hybrid reservoir computing,” Chaos: An Interdis- based modeling for scalable forecasting and subgrid-scale closure of large,
ciplinary Journal of Nonlinear Science 33 (2023). complex, spatiotemporal systems,” Chaos: An Interdisciplinary Journal of
15 U. Parlitz and C. Merkwirth, “Prediction of spatiotemporal time series Nonlinear Science 30 (2020).
based on reconstructed local states,” Physical review letters 84, 1890 19 “Zimmermann’s source code of the barkley model,” [Link]

(2000). com/zimmerrol/rcp_spatio_temporal, accessed: 2024-03-25.

16 J. Pathak, B. Hunt, M. Girvan, Z. Lu, and E. Ott, “Model-free prediction of 20 S. Baur and C. Räth, “Predicting high-dimensional heterogeneous time se-

large spatiotemporally chaotic systems from data: A reservoir computing ries employing generalized local states,” Phys. Rev. Res. 3, 023215 (2021).
approach,” Physical review letters 120, 024102 (2018). 21 J. Herteux and C. Räth, “Breaking symmetries of the reservoir equa-

tions in echo state networks,” Chaos: An Interdisciplinary Journal of

Nonlinear Science 30, 123142 (2020), [Link]
pdf/doi/10.1063/5.0028993/14110522/123142_1_online.pdf.