* (29) M. Luscher, Symmetry Breaking Aspects of the Roughening Transition in Gauge Theories, Nucl. Phys. B 180 (1981) 317–329.
  * (30) R. Sundrum, Hadronic string from confinement, hep-ph/9702306.
  * (31) S. S. Gubser, I. R. Klebanov, and A. M. Polyakov, Gauge theory correlators from noncritical string theory, Phys. Lett. B428 (1998) 105–114, [hep-th/9802109].
  * (32) J. Polchinski and M. J. Strassler, Deep inelastic scattering and gauge / string duality, JHEP 05 (2003) 012, [hep-th/0209211].
  * (33) A. Dymarsky, K. Farnsworth, Z. Komargodski, M. A. Luty, and V. Prilepina, Scale Invariance, Conformality, and Generalized Free Fields, JHEP 02 (2016) 099, [arXiv:1402.6322].
  * (34) J. Kaplan, “Lectures on AdS/CFT from the Bottom Up.”
  * (35) W. Thirring, Lagrangian formulation of the zachariasen model, Phys. Rev. 126 (May, 1962) 1209–1212.
  * (36) N. G. Deshpande and X.-G. He, Unparticle Realization Through Continuous Mass Scale Invariant Theories, Phys. Rev. D 78 (2008) 055006, [arXiv:0806.2009].
  * (37) G. Dvali and C. Gomez, Quantum Information and Gravity Cutoff in Theories with Species, Phys. Lett. B 674 (2009) 303–307, [arXiv:0812.1940].
  * (38) O. Aharony, A Brief review of ’little string theories’, Class. Quant. Grav. 17 (2000) 929–938, [hep-th/9911147].
  * (39) L. Randall and R. Sundrum, An Alternative to compactification, Phys. Rev. Lett. 83 (1999) 4690–4693, [hep-th/9906064].
  * (40) K. R. Dienes and B. Thomas, Dynamical Dark Matter: I. Theoretical Overview, Phys. Rev. D 85 (2012) 083523, [arXiv:1106.4546].
  * (41) K. R. Dienes and B. Thomas, Dynamical Dark Matter: II. An Explicit Model, Phys. Rev. D 85 (2012) 083524, [arXiv:1107.0721].
  * (42) S. S. Gubser, AdS / CFT and gravity, Phys. Rev. D63 (2001) 084017, [hep-th/9912001].
  * (43) T. Shiromizu, K.-i. Maeda, and M. Sasaki, The Einstein equation on the 3-brane world, Phys. Rev. D 62 (2000) 024012, [gr-qc/9910076].
  * (44) P. Binetruy, C. Deffayet, U. Ellwanger, and D. Langlois, Brane cosmological evolution in a bulk with cosmological constant, Phys. Lett. B 477 (2000) 285–291, [hep-th/9910219].
  * (45) A. Hebecker and J. March-Russell, Randall-Sundrum II cosmology, AdS / CFT, and the bulk black hole, Nucl. Phys. B 608 (2001) 375–393, [hep-ph/0103214].
  * (46) D. Langlois, L. Sorbo, and M. Rodriguez-Martinez, Cosmology of a brane radiating gravitons into the extra dimension, Phys. Rev. Lett. 89 (2002) 171301, [hep-th/0206146].
  * (47) D. Langlois and L. Sorbo, Bulk gravitons from a cosmological brane, Phys. Rev. D 68 (2003) 084006, [hep-th/0306281].
  * (48) A. Karch, E. Katz, D. T. Son, and M. A. Stephanov, Linear confinement and AdS/QCD, Phys. Rev. D74 (2006) 015005, [hep-ph/0602229].
  * (49) U. Gursoy and E. Kiritsis, Exploring improved holographic theories for QCD: Part I, JHEP 02 (2008) 032, [arXiv:0707.1324].
  * (50) U. Gursoy, E. Kiritsis, and F. Nitti, Exploring improved holographic theories for QCD: Part II, JHEP 02 (2008) 019, [arXiv:0707.1349].
  * (51) S. S. Gubser and A. Nellore, Mimicking the QCD equation of state with a dual black hole, Phys. Rev. D 78 (2008) 086007, [arXiv:0804.0434].
  * (52) A. Falkowski and M. Perez-Victoria, Electroweak Breaking on a Soft Wall, JHEP 12 (2008) 107, [arXiv:0806.1737].
  * (53) B. Batell and T. Gherghetta, Dynamical Soft-Wall AdS/QCD, Phys. Rev. D78 (2008) 026002, [arXiv:0801.4383].
  * (54) B. Batell, T. Gherghetta, and D. Sword, The Soft-Wall Standard Model, Phys. Rev. D78 (2008) 116011, [arXiv:0808.3977].
  * (55) J. A. Cabrer, G. von Gersdorff, and M. Quirós, Soft-Wall Stabilization, New J. Phys. 12 (2010) 075012, [arXiv:0907.5361].
  * (56) G. von Gersdorff, From Soft Walls to Infrared Branes, Phys. Rev. D82 (2010) 086010, [arXiv:1005.5134].
  * (57) J. A. Cabrer, G. von Gersdorff, and M. Quirós, Suppressing Electroweak Precision Observables in 5D Warped Models, JHEP 05 (2011) 083, [arXiv:1103.1388].
  * (58) E. Megías and M. Quirós, Gapped Continuum Kaluza-Klein spectrum, JHEP 08 (2019) 166, [arXiv:1905.07364].
  * (59) P. Brax and C. van de Bruck, Cosmology and brane worlds: A Review, Class. Quant. Grav. 20 (2003) R201–R232, [hep-th/0303095].
  * (60) J. W. York, Jr., Role of conformal three geometry in the dynamics of gravitation, Phys. Rev. Lett. 28 (1972) 1082–1085.
  * (61) G. W. Gibbons and S. W. Hawking, Action Integrals and Partition Functions in Quantum Gravity, Phys. Rev. D 15 (1977) 2752–2756.
  * (62) E. Megías, G. Nardini, and M. Quirós, Cosmological Phase Transitions in Warped Space: Gravitational Waves and Collider Signatures, JHEP 09 (2018) 095, [arXiv:1806.04877].
  * (63) T. Tanaka and Y. Himemoto, Generation of dark radiation in bulk inflaton model, Phys. Rev. D 67 (2003) 104007, [gr-qc/0301010].
  * (64) S. Fichet, On holography in general background and the boundary effective action from AdS to dS, JHEP 07 (2022) 113, [arXiv:2112.00746].
  * (65) B. Batell and T. Gherghetta, Holographic mixing quantified, Phys. Rev. D76 (2007) 045017, [arXiv:0706.0890].
  * (66) S. L. Dubovsky, V. A. Rubakov, and P. G. Tinyakov, Brane world: Disappearing massive matter, Phys. Rev. D 62 (2000) 105011, [hep-th/0006046].
  * (67) L. Di Pietro, V. Gorbenko, and S. Komatsu, Analyticity and Unitarity for Cosmological Correlators, arXiv:2108.01695.
  * (68) S. Fichet, Dressing in AdS and a Conformal Bethe-Salpeter Equation, arXiv:2106.04604.
  * (69) S. S. Gubser and I. Mitra, Double trace operators and one loop vacuum energy in AdS / CFT, Phys. Rev. D 67 (2003) 064018, [hep-th/0210093].
  * (70) P. Brax, S. Fichet, and P. Tanedo, The Warped Dark Sector, Phys. Lett. B 798 (2019) 135012, [arXiv:1906.02199].
  * (71) R. Zwicky, A Brief Introduction to Dispersion Relations and Analyticity, in Quantum Field Theory at the Limits: From Strong Fields to Heavy Quarks, pp. 93–120, 2017.  arXiv:1610.06090.
  * (72) K. Hinterbichler and J. Khoury, Symmetron Fields: Screening Long-Range Forces Through Local Symmetry Restoration, Phys. Rev. Lett. 104 (2010) 231301, [arXiv:1001.4525].
  * (73) P. Callin and F. Ravndal, Higher order corrections to the Newtonian potential in the Randall-Sundrum model, Phys. Rev. D 70 (2004) 104009, [hep-ph/0403302].
  * (74) S. Fichet, Braneworld effective field theories — holography, consistency and conformal effects, JHEP 04 (2020) 016, [arXiv:1912.12316].
  * (75) E. Megías and M. Quirós, The Continuum Linear Dilaton, Acta Phys. Polon. B 52 (2021) 711, [arXiv:2104.10260].
  * (76) U. Gursoy, E. Kiritsis, L. Mazzanti, and F. Nitti, Holography and Thermodynamics of 5D Dilaton-gravity, JHEP 05 (2009) 033, [arXiv:0812.0792].
  * (77) L. Randall and M. D. Schwartz, Quantum field theory and unification in AdS5, JHEP 11 (2001) 003, [hep-th/0108114].
  * (78) A. Costantino and S. Fichet, Opacity from Loops in AdS, arXiv:2011.06603.
  * (79) E. Megías and M. Quirós, Analytical Green’s functions for continuum spectra, JHEP 09 (2021) 157, [arXiv:2106.09598].
  * (80) S. J. Smullin, A. A. Geraci, D. M. Weld, J. Chiaverini, S. P. Holmes, and A. Kapitulnik, New constraints on Yukawa-type deviations from Newtonian gravity at 20 microns, Phys. Rev. D72 (2005) 122001, [hep-ph/0508204]. [Erratum: Phys. Rev.D72,129901(2005)].
  * (81) W. Hu, Structure Formation with Generalized Dark Matter, APJ 506 (Oct., 1998) 485–494, [astro-ph/9801234].
  * (82) M. Kopp, C. Skordis, D. B. Thomas, and S. Ilić, Dark Matter Equation of State through Cosmic History, Phys. Rev. Lett. 120 (2018), no. 22 221102, [arXiv:1802.09541].
  * (83) D. B. Thomas, M. Kopp, and K. Markovič, Using large-scale structure data and a halo model to constrain generalized dark matter, Mon. Not. Roy. Astron. Soc. 490 (2019), no. 1 813–831, [arXiv:1905.02739].
  * (84) T. Lin, Dark matter models and direct detection, PoS 333 (2019) 009, [arXiv:1904.07915].
  * (85) P. Brax, S. Fichet, and G. Pignol, Bounding Quantum Dark Forces, Phys. Rev. D97 (2018), no. 11 115034, [arXiv:1710.00850].
  * (86) A. Costantino, S. Fichet, and P. Tanedo, Effective Field Theory in AdS: Continuum Regime, Soft Bombs, and IR Emergence, Phys. Rev. D 102 (2020), no. 11 115038, [arXiv:2002.12335].
  * (87) S. Fichet, Opacity and effective field theory in anti–de Sitter backgrounds, Phys. Rev. D100 (2019), no. 9 095002, [arXiv:1905.05779].
  * (88) I. Antoniadis, A. Arvanitaki, S. Dimopoulos, and A. Giveon, Phenomenology of TeV Little String Theory from Holography, Phys. Rev. Lett. 108 (2012) 081602, [arXiv:1102.4043].
  * (89) P. Cox and T. Gherghetta, Radion Dynamics and Phenomenology in the Linear Dilaton Model, JHEP 05 (2012) 149, [arXiv:1203.5870].


# Adapting reservoir computing to solve the Schrödinger equation

L.Domingo  Instituto de Ciencias Matemáticas (ICMAT); Campus de Cantoblanco
UAM; Nicolás Cabrera, 13-15; 28049 Madrid (Spain)  Departamento de Química;
Universidad Autónoma de Madrid; CANTOBLANCO - 28049 Madrid (Spain)  Grupo de
Sistemas Complejos; Universidad Politécnica de Madrid; 28035 Madrid (Spain)
<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>J. Borondo
Departamento de Gestión Empresarial; Universidad Pontificia de Comillas ICADE;
Alberto Aguilera 23; 28015 Madrid (Spain)  AgrowingData; Navarro Rodrigo 2 AT;
04001 Almería (Spain)  F. Borondo  Instituto de Ciencias Matemáticas (ICMAT);
Campus de Cantoblanco UAM; Nicolás Cabrera, 13-15; 28049 Madrid (Spain)
Departamento de Química; Universidad Autónoma de Madrid; CANTOBLANCO - 28049
Madrid (Spain)

###### Abstract

Reservoir computing is a machine learning algorithm that excels at predicting
the evolution of time series, in particular, dynamical systems. Moreover, it
has also shown superb performance at solving partial differential equations.
In this work, we adapt this methodology to integrate the time-dependent
Schrödinger equation, propagating an initial wavefunction in time. Since such
wavefunctions are complex-valued high-dimensional arrays the reservoir
computing formalism needs to be extended to cope with complex-valued data.
Furthermore, we propose a multi-step learning strategy that avoids overfitting
the training data. We illustrate the performance of our adapted reservoir
computing method by application to four standard problems in molecular
vibrational dynamics.

Reservoir computing is a machine learning method that has proven very useful
in predicting the evolution of time series, such as weather forecast, stock
index or the dynamics of chaotic systems. In this work, we adapt this machine
learning algorithm to predict the evolution of quantum systems, whose dynamics
is determined by the Schrödinger equation. In this context, the quantum states
are described using large complex-valued arrays, which makes the prediction
difficult with classical reservoir computing mainly due to overfitting. For
this reason, we propose a new reservoir computing-based methodology that
allows to efficiently propagate quantum states in time. We illustrate the
method and its performance with several quantum systems.

##  I Introduction

Recurrent neural networks (RNN) Rumelhart, Hinton, and Williams (1986);
Hochreiter and Schmidhuber (1997) are a deep learning tool used to deal with
sequential data, such as text or time series. Within this context, RNNs have
proven to give optimal results, due to their ability to keep the knowledge
learnt from past data. However, RNNs are known to be computationally expensive
during training, and such training may be complicated due to the problem of
exploding gradients during backpropagation Doya (2000). To solve these
problems a new trend of RNN, called Reservoir Computing (RC), was designed.

Indeed, Jaeger et al Jaeger (2001, 2007) proved that as long as the reservoir
fullfils certain properties, only training the readout layer is enough to
obtain excellent performance in many tasks. The concept of RC was proposed
simultaneously as Echo State Networks (ESNs) Jaeger (2001) in the field of
machine learning, and as Liquid State Machines (LSMs) Maass, Natschläger, and
Markram (2002) in computational neuroscience. Actually, reservoir computing
has proven useful for chaotic time series prediction Pathak  _et al._ (2017,
2018). The underlying idea is that a well-trained reservoir is able to
reproduce the attractor of the originating chaotic dynamics.

Recently, relevant advances have been made in developing RC-based algorithms.
Multiple stacked reservoirs have been used to increase the reservoir
performance Akiyama and Tanaka (2019); Qiao  _et al._ (2017). Also, RC has
been adapted to predict time series with short training data sets Chen  _et
al._ (2020). For predicting low dimensional time series, a RC algorithm
requiring no random matrices and fewer hyperparameters has shown to be
equivalent to standard RC Gauthier  _et al._ (2021).

In the field of quantum chemistry, machine learning methods are proving to
have advantages over the standard computational chemistry approaches Ferguson
_et al._ (2020), especially in the case of high-dimensional systems. In
particular, machine learning has been used to solve the time-independent
Schrödinger equation for the electronic Pfau  _et al._ (2020); Hermann,
Schätzle, and Noé (2020) and vibrational (ground state Manzhos, Yamashita, and
Carrington (2009); Pavlov, Serdyuk, and Ustinov (2019); Mills, Spanner, and
Tamblyn (2017) and high-energy states Domingo and Borondo (2021)) parts of the
total wavefunction within the Born-Oppenheimer approximation.

In this work, we aim to use RC to propagate wave packets with time, that is,
to solve the time-dependent Schrödinger equation in continuous quantum
systems. To do so, the RC model needs to be adapted to work with
wavefunctions, which are complex-valued high-dimensional matrices. We also
propose here a new learning strategy that allows propagating wavefunctions
while reducing the overfitting of the training data. Such learning strategy
consists of a two-step training of the readout layer, where the reservoir is
shown how predicting unseen data affects the evolution of the internal states.
The learning algorithm is adapted to prevent fast error propagation during the
test phase. A more detailed explanation of the algorithm is given in Sect.
IV.2.

The organization of this paper is as follows. In Sect. II we review the
original formulation of RC introduced in Ref. Jaeger, 2001. In Sect. III we
present the quantum setting and the challenges that it presents to standard
RC. Section IV describes the method we propose to overcome these challenges.
Section V describes the four quantum vibrational systems that are used for
illustration in this work. The corresponding results are presented in Sect.
VI. Finally, Sect. VII ends the paper by summarizing the main conclusions of
the present work, and presenting an outlook for future work..

##  II Reservoir computing

In the RC framework, the learning complexity of the algorithm is reduced to
performing a linear regression. The key point is to design a dynamical system
(the internal network) that learns the input-output dynamics. For this reason,
the internal states of the network are called echo states, since they can be
thought of as an echo of their past Jaeger (2001). Such internal states are
unambiguously determined by the input and output of the network. The structure
of the network is shown in Fig. 1,

Figure 1: Architecture of a reservoir computing model. Only the readout layer
$W^{\text{out}}$ is learnt during training.

where

  * •

$W^{\text{in}}\in M(\mathbb{R})_{N\times K}$ gives the weights from the input
units to the internal states,

  * •

$W\in M(\mathbb{R})_{N\times N}$ gives the weights between the different
internal states,

  * •

$W^{\text{out}}\in M(\mathbb{R})_{L\times N}$ gives the wights from internal
states to output units, and

  * •

$W^{\text{back}}\in M(\mathbb{R})_{N\times L}$ gives the weights from the
output units to the internal states,

and where

  * •

$u(t)=(u_{1}(t),u_{2}(t),\ldots,u_{K}(t))$ is a $K$-dimensional vector giving
the input units at time $t$,

  * •

$x(t)=(x_{1}(t),x_{2}(t),\ldots,x_{N}(t)))$ is an $N$-dimensional vector
giving the internal states at time $t$, and

  * •

$y(t)=(y_{1}(t),y_{2}(t),\ldots,y_{L}(t)))$ is an $L$-dimensional vector
giving the output units at time $t$.

Notice that the matrices $W^{\text{in}}$, $W$ and $W^{\text{back}}$ are fixed,
so that they do not change during training. The only learnable parameters are
the weights $W^{\text{out}}$. The steps to train the ESN are the following:

  1. 1.

Generate the reservoir matrices $(W,W^{\text{in}},W^{\text{back}})$ randomly.

  2. 2.

Update the internal states by teacher forcing:

$\displaystyle\tilde{x}(t)$
$\displaystyle=f((W^{\text{in}}u(t)+Wx(t-1)+W^{\text{back}}y_{\text{teach}}(t-1))$
$\displaystyle x(t)$ $\displaystyle=(1-\alpha)x(t-1)+\alpha\tilde{x}(t)$ (1)

where $y_{\text{teach}}$ is the output that we want out network to predict,
$\alpha\in(0,1]$ is the leaking rate and $f$ is the activation function.
Usually, $f(\cdot)=\tanh(\cdot)$, which is applied component-wise.

  3. 3.

Discard a transient of $t_{\text{min}}$ states to guarantee the convergence of
the reservoir dynamics.

  4. 4.

Find the readout matrix $W^{\text{out}}$ by minimizing the mean square error:

$\displaystyle MSE(y,y_{\text{teach}})=$
$\displaystyle\frac{1}{T-t_{\text{min}}}\sum_{t=t_{\min}}^{T}\Big{[}f^{{\text{out}}^{-1}}(y_{\text{teach}}(t))$
(2) $\displaystyle-W^{\text{out}}x(t)\Big{]}^{2},$

where $f^{\text{out}}$ is the output activation function (usually
$f^{\text{out}}=Id$). Notice that this step only requires performing a linear
regression.

The steps to make the predictions after training the network are the
following:

  1. 1.

Given an input $u(t)$, update the state of the reservoir:

$\displaystyle\tilde{x}(t)$
$\displaystyle=f[(W^{\text{in}}u(t)+Wx(t-1)+W{{}^{\text{back}}}y(t-1)]$ (3)
$\displaystyle x(t)$ $\displaystyle=(1-\alpha)x(t-1)+\alpha\tilde{x}(t)$

where $y(t-1)$ is the prediction of the output at time $t-1$.

  2. 2.

Compute the prediction $y(t)$:

$y(t)=f^{\text{out}}\Big{[}W^{\text{out}}x(t)\Big{]}.$ (4)

Let us explain now these steps in more detail. To ensure that the RC model
works, the internal network must fulfil the echo state properties. Roughly
speaking, this means that the internal states must have the state and input
forgetting property. That is, for $t$ large enough, $x(t)$ does not depend on
$x(0),u(0)$ or $y_{\text{teach}}(0)$. For this reason, in step 2 of the
training phase, we dismiss the initial values of the internal states, which
could be influenced by the initial parameters of the reservoir. The spectral
radius $\rho(W)$ of the internal matrix $W$ also influences the performance of
the method. Larger $\rho(W)$ provide longer memory. However, for $\rho(W)>>1$,
the reservoir will likely have no echo states. In fact, it is a sufficient
condition that $\rho(W)<1$ to have echo states, even though reservoirs with
$\rho(W)$ slightly larger than 1 may also give optimal results Jaeger (2001).
Also, the matrix $W$ should be sparse and in-homogeneous so that the internal
states contain a diverse set of trajectories. On the other hand, the leaking
rate $\alpha$ influences the velocity of the dynamics of the output. A fast-
changing output should be trained with greater values of $\alpha$. However,
the RC model is also frequently used without the leaky integration, which is a
special case obtained by setting $\alpha=1$ and thus $\tilde{x}(t)=x(t)$.

Once the internal states have been computed, the learning phase only consists
of training a linear model to find the mapping $W^{\text{out}}$ from the
internal states to the output. This linear model is usually a ridge
regression, which minimizes the following expression

$\displaystyle MSE_{r}(y,y_{\text{teach}})=$
$\displaystyle\frac{1}{T-t_{\text{min}}}\sum_{t=t_{\min}}^{T}\Big{(}f^{{\text{out}}^{-1}}(y_{\text{teach}}(t))$
(5)
$\displaystyle-W^{\text{out}}x(t)\Big{)}^{2}+\gamma||W^{\text{out}}||^{2}.$

This ridge regression is used to reduce overfitting during training. Other
methods, such as adding noise to the input during the teacher forcing phase,
or adding a random constant input can also be used to reduce overfitting.
Finally, notice that the RC model can also be used without the input layer
when we aim to predict a time series without explanatory variables. In this
case, the term $W^{\text{in}}u(t)$ is removed from Eqs. (1) and (1).

##  III Reservoir computing challenges in quantum problems

###  III.1 The quantum setting

In this work, we aim to solve the time-dependent Schrödinger equation given by

$i\hbar\frac{\partial\psi(\vec{x},t)}{\partial t}=\hat{H}\psi(\vec{x},t),$ (6)

where $\hat{H}$ is the Hamiltonian which describes the quantum system. The
solution of Eq. (6), given an initial wavefunction $\psi(\vec{x},0)$, provides
the time-evolution of the associated quantum state. The propagation of a
quantum state with certain (mean) energy allows the calculation, by Fourier
transformation, of the eigenenergies and eigenfunctions of the quantum system
around the same energy. Full details details on how to obtain the
eigenenergies and eigenfunctions from $\psi(\vec{x},t)$ are provided in Sect.
V. Notice that this method present great advantages over the usual variational
method in the case of excited states, since it does not require the
calculation of the low lying states Domingo and Borondo (2021).

In this work, we aim to develop a RC model that can be trained with the short-
term evolution of a wave packet and then use RC to predict its longer-term
evolution. The wavefunction $\psi(\vec{x},t)$ is a complex-valued function of
the spatial coordinates $\vec{x}$, which are usually multidimensional. In the
RC framework, $\psi(\vec{x},t)$ is represented as a set of matrices
$\\{\psi(\vec{x},t)\\}_{t}$, where each matrix contains the values of
$\psi(\vec{x})$ at time $t$ in a grid of points spanning $\vec{x}$. The use of
$\psi(\vec{x},t)$ as the target of the reservoir presents two challenges in
the RC original formulation, which are described in the next subsections.

###  III.2 Complex numbers

The usual RC framework is done with real numbers. That is, all the inputs
$u(t)$, outputs $y(t)$, internal states $x(t)$, and weight matrices
$W,W^{\text{in}},W^{\text{back}}$ and $W^{\text{out}}$ are real-valued.
However, the target time series for this quantum problem is a wavefunction
$\psi(\vec{x},t)$ which, in general, takes complex values. Therefore, we need
to extend update of the internal state and the learning algorithm to complex-
valued data.

###  III.3 High dimensional data

The target data is represented as a matrix, where each entry is the value of
the wavefunction in a discretized spatial grid. The size of this matrix
increases exponentially with the dimension of the physical system, that is,
the dimension of $\vec{x}$. In the usual RC setting, we use a large reservoir
$W$, compared to the dimension of the input. However, in the quantum problems,
the target can have a very high dimension, and thus it is not feasible to
design a reservoir much larger than the target size. Moreover, the bigger the
reservoir, the more data are needed to train the linear model. If the
complexity of the method is is too large, the linear model may overfit the
training data. Therefore, we need to adapt the RC framework to avoid this
overfitting when dealing with high-dimensional data.

##  IV Reservoir computing advances for quantum data

This section presents the modifications done in the standard RC framework in
order to adapt to the quantum setting.

###  IV.1 Complex-valued ridge regression

We begin by presenting the extension of the RC framework to complex-valued
arrays. The update of the internal state consists of matrix multiplications
and an application of a non-linear function, in the following way

$\displaystyle\tilde{x}(t)$ $\displaystyle=$ $\displaystyle
f[(W^{\text{in}}u(t)+Wx(t-1)+W^{\text{back}}y_{\text{teach}}(t-1)]$
$\displaystyle x(t)$ $\displaystyle=$
$\displaystyle(1-\alpha)\;x(t-1)+\alpha\tilde{x}(t)$ (7)

As long as the function $f$ is defined for complex numbers, the previous
equation holds for complex-valued arrays. In this work, we propose the use of
two activation functions:

  * •

$f_{1}(x)=\tanh(x)$

  * •

$f_{2}(x)=\tanh(\Re(x))+i\;\tanh(\Im(x))$

Also, we set $f^{\text{out}}=Id$, where $Id$ is the identity function. Once
the internal states have been calculated, the matrix $W^{\text{out}}$ is
calculated by minimizing the MSE with $L^{2}$ regularization. This corresponds
to performing a complex-valued ridge regression

$\displaystyle MSE_{r}(y,y_{\text{teach}})=$
$\displaystyle\frac{1}{T-t_{\text{min}}}\displaystyle\sum_{t=t_{\min}}^{T}\Big{|}f_{\text{out}}^{-1}(y_{\text{teach}}(t))$
(8) $\displaystyle-W^{\text{out}}x(t)\Big{|}^{2}+\alpha|W^{\text{out}}|^{2},$

where $|\cdot|$ denotes the $L^{2}$ norm in complex values. In matrix form
this equation becomes

$\displaystyle MSE_{r}(y,y_{\text{teach}})=$
$\displaystyle(f_{\text{out}}^{-1}(\vec{y}_{\text{teach}})-W^{\text{out}}X)^{*}\cdot(f_{\text{out}}^{-1}(\vec{y}_{\text{teach}})-W^{\text{out}}X)$
(9) $\displaystyle+\alpha||W^{\text{out}}||^{2},$

where ∗ denotes the conjugate transpose, and $X$ is the matrix containing the
internal states. In the case of real values, the conjugate transpose is just
the transpose. This linear model has a closed solution

$W^{\text{out}}=\Big{(}X^{*}X+\alpha\mathbb{I}\Big{)}^{-1}\cdot\Big{(}X^{*}f_{\text{out}}^{-1}(\vec{y}_{\text{teach}})\Big{)},$
(10)

and therefore, an analytical solution for the linear model with complex-valued
data can be computed.

###  IV.2 Multiple-step training

Figure 2: Architecture of multi-step training of the reservoir computing
model.

The second and main challenge of dealing with quantum wavefunctions is the
high-dimensionality of the target matrices. In order to capture all the
underlying information of the dynamical system, the size of the reservoir $W$
should be at least of the order of magnitude of the target data. According to
Hoeffding’s theorem Mohri, Rostamizadeh, and Talwalkar (2018), models with
high complexity require a large amount of data to reduce the variance of the
model predictions, following the (Hoeffding’s) inequality

$E_{out}\leq E_{in}+\mathcal{O}\Big{(}\sqrt{\frac{C}{N}}\Big{)},$ (11)

where $E_{out}$ is the error in the test set, $E_{in}$ is the error in the
training set, $C$ is a notion of complexity and $N$ is the number of data
samples. In our case, the size of the reservoir is similar to the number of
data samples, and thus the method is likely to produce overfitting. To reduce
this effect, we propose the use of a multi-step learning training. The
training steps in this method, depicted schematically in Fig. 2, are the
following:

  1. 1.

Split the training input and output into two parts.

$\displaystyle\vec{y}=$ $\displaystyle\begin{pmatrix}y(1)\\\ \vdots\\\
y(t_{1})\\\ y(t_{1}+1)\\\ \vdots\\\
y(T)\end{pmatrix}=\begin{pmatrix}y_{1}(1)\\\ \vdots\\\ y_{1}(t_{1})\\\
y_{2}(t_{1}+1)\\\ \vdots\\\ y_{2}(T)\\\
\end{pmatrix}=\begin{pmatrix}\vec{y_{1}}\\\ \vec{y_{2}}\\\ \end{pmatrix},$
$\displaystyle\vec{u}=$ $\displaystyle\begin{pmatrix}u(1)\\\ \vdots\\\
u(t_{1})\\\ u(t_{1}+1)\\\ \vdots\\\
u(T)\end{pmatrix}=\begin{pmatrix}u_{1}(1)\\\ \vdots\\\ u_{1}(t_{1})\\\
u_{2}(t_{1}+1)\\\ \vdots\\\ u_{2}(T)\\\
\end{pmatrix}=\begin{pmatrix}\vec{u_{1}}\\\ \vec{u_{2}}\\\ \end{pmatrix}$ (12)

  2. 2.

Find the internal states of the reservoir $x(1),\cdots x(t_{1})$ with
$\vec{y_{1}}$ and $\vec{u_{1}}$, and perform a ridge regression to find
$W^{\text{out}}$.

  3. 3.

Evolve the internal states of the reservoir $x(t_{1}+1),\ldots,x(T)$ by
letting the reservoir predict the values of $\vec{y_{2}}$.

  4. 4.

With all the internal states $x(1),\cdots,x(t_{1}),x(t_{1}+1),\ldots x(T)$
retrain the linear model to correct $W^{\text{out}}$.

When the linear model overfeeds the training data, small changes in the
internal states due to prediction errors produce high errors in the prediction
of the wavefunctions. Therefore, the errors tend to propagate fast until the
reservoir is not able to correctly propagate the quantum state. In our
proposed method, we allow the reservoir evolve after training by predicting a
part of the training data. This simulates the real situation found during
testing. During this predicting phase, the internal states are modified using
the predicted output instead of the exact ones (see Eq. 1). The predicted
outputs have some errors, which are transmitted to the internal states, which
are the used to retrain the linear model. This learning strategy shows the
reservoir how small predicting errors modify the internal states during the
test phase. This reduces overfitting, decreasing the test error significantly.

##  V Studied quantum systems

In this section we describe the quantum mechanical calculations, and we also
present the quantum systems chosen to be studied in this work to illustrate
the performance of our multi-step RC method.

###  V.1 Quantum mechanical calculations

In all cases, the training procedure is the same, and it consists of the
following steps:

  1. 1.

Generate training and test data by numerically integrating the Schrödinger
equation. For this purpose, we used the widespread Fast Fourier Transform
method proposed by Kosloff and Kosloff Kosloff and Kosloff (1983).

  2. 2.

Train the RC and multi-step RC models. Then, compute the predicted
wavefunctions $\psi(\vec{x},t)$ and evaluate the errors in both cases.

From the time evolution $\psi(\vec{x},t)$ of a suitable initial quantum state
$\psi_{0}(\vec{x})$ the eigenenergies and eigenfunctions of the quantum system
with similar energies to that of $\psi(\vec{x},0)$ can be calculated, by
subsequent Fourier transform at the quantized energies. The eigenenergies and
eigenstates completely describe the stationary states of the quantum system.
Therefore, the error of these eigenenergies gives another measure of the
performance of our RC models. To obtain these eigenstates we have to perform
two additional steps Zdánská and Moiseyev (2004).

  1. 3.

Calculate the energy spectrum $I(E)$. First we calculate the time-correlation
function

$A(t)=<\psi(t)|\psi(0)>=\int_{\vec{x}}\;\psi(t)^{*}\psi(0)\;d\vec{x}$ (13)

Then, we obtain the energy spectrum $I(E)$ as the Fourier transform of the
time-correlation function

$I(E)=\int_{t}\;A(t)\;e^{iEt/\hbar}\;dt.$ (14)

The eigenenergies $E_{n}$ correspond to the peaks of the energy spectrum
function.

  2. 4.

Calculate the eigenfunctions $\phi_{n}(\vec{x})$ by performing the Fourier
transform of $\psi(\vec{x},t)$ at the eigenenergies

$\phi_{n}(\vec{x})=\int_{t}\;\psi(\vec{x},t)\;e^{-iE_{n}t/\hbar}\;dt.$ (15)

  3. 5.

Evaluate the errors produced by the reservoir computing and multi-step
reservoir computing models in the eigenstates and eigenenergies.

Now we present the different quantum systems used for this study, which
include three one-dimensional quantum systems and one two-dimensional system.

###  V.2 1D systems

We aim to solve the time-dependent Schrödinger equation (6) for 1D systems

$i\hbar\frac{\partial\psi(x,t)}{\partial t}=\hat{H}\psi(x,t)$ (16)

In all cases, we set $\hbar=1$ and $m=1$ without loss of generality. We use
5000 time steps both for training and for test. The spatial domains are chosen
so that they contain the wavefunctions at all times $t$. The systems under
study are the following:

  1. 1.

Harmonic oscillator. In this case, the Hamiltonian $\hat{H}$ is

$\hat{H}(x)=\frac{1}{2}\frac{\partial^{2}}{\partial
x^{2}}+\frac{\omega^{2}}{2}(x-x_{0})^{2}.$ (17)

We choose $x_{0}=0$ and $\omega=1$. The (computationally effective) spatial
domain is $[-10,10]$ covered with a grid of 200 equidistant points. That is,
$\psi(x,t)$ for a fixed $t$ is a vector of size 200. The time interval between
time steps is $dt=0.002$ The initial wavefunction $\psi_{0}(x)$ is the minimum
uncertainty Gaussian wave packet

$\psi(x,0)=\left(\frac{1}{\pi}\right)^{1/4}e^{-(x-x_{0})^{2}/2}\;e^{ip_{0}x},$
(18)

which is centered at $x_{0}=0$, and it has a width of $\displaystyle\Delta
x=1/(2\sqrt{2})$, a phase/momentum $p_{0}=3.35$, and a (mean) energy $E=6.0$.

  2. 2.

Morse potential. The Morse Hamiltonian adequately represent the potential
interaction of a diatomic molecule. We write the corresponding Hamiltonian
$\hat{H}$ as

$\hat{H}(x)=\frac{1}{2}\frac{\partial^{2}}{\partial
x^{2}}+D_{e}\left(e^{-2a(x-x_{0})}-2e^{-a(x-x_{0})}\right),$ (19)

where $x$ is the distance between atoms, $x_{0}$ is the corresponding
equilibrium bond distance, $D_{e}$ is the well depth (defined relative to the
dissociated atoms energy value), and $a$ is a parameter controlling the
curvature of the potential at its minimum. We choose $D_{e}=7$, $a=0.09$ and
$x_{0}=0$. The spatial domain is $[-10,25]$ with 140 points. The time interval
between time steps id $dt=0.007$. The initial wavefunction $\psi(x,0)$ is
again the minimum uncertainty packet (Eq. 18) centered at $x_{0}=0$, with a
phase $p_{0}=2$, and an energy $E=-4.8$.

  3. 3.

Polynomial potential The last quantum system considered corresponds to a
quartic polynomial potential

$\hat{H}(x)=\frac{1}{2}\frac{\partial^{2}}{\partial
x^{2}}+\alpha_{0}+\alpha_{1}x+\alpha_{2}x^{2}+\alpha_{3}x^{3}+\alpha_{4}x^{4},$
(20)

where $\alpha_{0}=-0.5,\alpha_{1}=0.14,\alpha_{2}=0.09,\alpha_{3}=-0.01,$ and
$\alpha_{4}=0.001$. The spatial domain here is $[-15,25]$ with 150 points. The
time interval between time steps is $dt=0.007$. The initial wavefunction
$\psi(x,0)$ is again the minimum uncertainty packet (Eq. 18) centered at
$x_{0}=0$, with a phase $p_{0}=2$, and an energy $E=1.8$

###  V.3 2D system

Apart from the 1D systems we also study one 2D system, consisting of a 2D
harmonic oscillator potential, whose Hamiltonian is given by

$\hat{H}(x,y)=\frac{1}{2}\left(\frac{\partial^{2}}{\partial
x^{2}}+\frac{\partial^{2}}{\partial
y^{2}}\right)+\frac{\omega_{x}^{2}}{2}(x-x_{0})^{2}+\frac{\omega_{y}^{2}}{2}(y-y_{0})^{2}.$
(21)

We choose $\omega_{x}^{2}=1,$ $\omega_{y}^{2}=\sqrt{2}$ and $x_{0}=y_{0}=0$.
Notice that this choice of frequencies prevents the appearance of resonances
in the dynamics and degeneracies in the eigenenergies. The initial
wavefunction $\psi(x,y,0)$ is the minimum uncertainty Gaussian wave packet
with energy $E=3.2$

$\psi_{0}(x,y,0)=\left(\frac{1}{\pi}\right)^{1/8}e^{-(x-x_{0})^{2}/(4\Delta
x)^{2}-(y-y_{0})^{2}/(4\Delta y)^{2}}\;e^{i\,p_{0}(x-y)},$ (22)

which is centered at $x_{0}=y_{0}=0$, has widths of $\Delta x=\Delta y=0.9$,
and a phase $p_{0}=1.75$. The spatial domain is $[-5,5]\times[-5,5]$ with 50
points in each dimension. Therefore, $\psi(x,y,t)$ for each $t$ is a matrix of
size $50\times 50$. Notice that the dimension of the data is of the order of
the training data. Therefore, the RC model is very likely to produce
overfitting, which motivates the multi-step learning we propose.

Quantum system | $\rho(W)$ | $\alpha$ | $t_{\text{min}}$ | $\gamma$ | $N$ |  $W$ density
---|---|---|---|---|---|---
Harmonic oscillator 1D | 1.10 | 0.015 | 300 | 0.1 | 2500 | 0.015
Morse potential | 0.75 | 0.015 | 500 | 0.5 | 1500 | 0.015
Polynomial potential | 0.75 | 0.017 | 50 | 0.05 | 1500 | 0.008
Harmonic oscillator 2D | 1.10 | 0.06 | 300 | 0.5 | 2000 | 0.015
Table 1: Reservoir computing training parameters, as defined in Sect. II, for
the four studied quantum systems studied in this work (see Sect. V).

Table 1 shows the training parameters, as defined in Sect. II, used for each
of the quantum systems considered in this work. The same parameter are used
both for the standard RC and multi-step RC models. For all systems we used
$f(x)=\tanh(\Re(x))+i\;\tanh(\Im(x))$, and $f^{\text{out}}=\mathbb{I}$. For
the multi-step RC we used 85% of the training data for the first training step
and 15% for the second step.

##  VI Results and discussion

In this section, we show present and discuss the results obtained with our
multi-step RC method for the four models described in Sect. V, which are
compared with those rendered by the standard RC model in order to obtain an
estimation of its performance.

Quantum system | Model | MSE $\psi(\vec{x},t)$ | MSE energies
---|---|---|---
Harmonic oscillator 1D | Standard RC | $6\cdot 10^{-5}$ | $3\cdot 10^{-4}$
| Multi-step RC | $2\cdot 10^{-5}$ | $1\cdot 10^{-5}$
Morse | Standard RC | $7\cdot 10^{-4}$ | $4\cdot 10^{-5}$
| Multi-step RC | $1\cdot 10^{-4}$ | $2\cdot 10^{-5}$
Polynomial | Standard RC | $2\cdot 10^{-4}$ | $3\cdot 10^{-5}$
| Multi-step RC | $8\cdot 10^{-5}$ | $8\cdot 10^{-6}$
Harmonic oscillator 2D | Standard RC | $3\cdot 10^{-4}$ | $4\cdot 10^{-3}$
| Multi-step RC | $2\cdot 10^{-5}$ | $8\cdot 10^{-5}$
Table 2: Mean squared error (MSE) for the wavefunctions $\psi(\vec{x},t)$ and
eigenenergies for the four quantum systems considered in this work (see Sect.
V) using the two studied RC models described in Sects. II and IV.

In the first place, we show in Table 2 the mean square error (MSE) of the
wavefunctions $\psi(\vec{x},t)$ and (mean) energies, for the different systems
chosen to study. Let us remark that these results constitute a very strict way
of gauging the performance of our method, since the quantum properties are
evaluated in all space using $\psi^{*}\psi$ as the probability density, which
means that only the regions where the wavefunction is high are meaningful. We
see that for all cases the MSE of both the wavefunctions and energies is
smaller when using the multi-step RC model. The largest difference in
performance appears in the 2D harmonic oscillator, where the MSE is one order
of magnitude smaller for the wavefunctions and two orders of magnitude smaller
for the energies. Therefore, only the multi-step learning RC algorithm can
correctly recover the eigenfunctions and eigenenergies of such quantum system.
Notice that the training data for the 2D harmonic oscillator consist of
matrices of size $50\times 50$ (2500 entries), while the 1D data are vectors
of size 100 to 200. Therefore, it is harder to train the RC model on 2D data,
since the reservoir size is of the same order of magnitude as the amount of
training data. Thus, the standard RC method overfits the training data and
does not generalize well to the test data. On the contrary, the multi-step
learning strategy refits the linear readout $W^{\text{out}}$ by showing the
reservoir how to adapt to unseen test data. This prevents having a fast error
propagation during the test phase.

Figure 3: Time evolution of the mean square error (MSE) for the harmonic, Morse, and polynomial oscillator time-dependent wavefunctions computed with the standard and multi-step reservoir computing methods. Figure 4: Time evolution of the mean square error (MSE) of the 2D harmonic oscillator wavefunctions for the standard RC and the multi-step RC. Quantum system | Predicted eigenenergies
---|---
Harmonic oscillator 1D | 2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 8.5, 9.5
Morse | -6.8326, -6.5040, -6.1834, -5.8710, -5.5666,
| -5.2704, -4.9822, -4.7022, -4.4302, -4.1664,
| -3.9106, -3.6630, -3.4234, -3.1920
Polynomial potential | 0.1556, 0.6187, 1.0800, 1.5426, 2.0087,
| 2.4801, 2.9582, 3.4453
Harmonic oscillator 2D | 3.0946, 3.2838, 3.4730
Table 3: Eigenenergies obtained from the spectra calculated from the predicted
wavefunctions $\psi(\vec{x},t)$ computed using the multi-step RC, for all four
studied quantum systems (see Sect. V). Figure 5: Energy spectra obtained using
the predicted wavefunctions from the multi-step RC, for all four studied
quantum systems. The orange peaks indicate the eigenenergies of each for each
of quantum system. Moreover, the numerical values of the eigenenergies are
reported given in Table 3.

In Figs. 3 and 4 we present the results of Table 2 in a graphic way, by
plotting the MSE of the predicted wavefunctions $\psi(\vec{x},t)$ as a
function of time. We see that due to error propagation the MSE tends to
increase with time, except for some random fluctuations. For all the studied
systems the MSE increases slower with time when using the multi-step RC.
Moreover, we see that for the 1D and 2D harmonic oscillators and the Morse
Hamiltonian, the MSE from the multi-step RC seems to stabilise, while the MSE
of the standard RC keeps increasing. Therefore, for a fixed error tolerance,
the multi-step RC method allows to predict a longer time evolution than the
standard RC.

After studying the performance of the RC models when propagating the
wavefunctions in time, we can recover the eigenenergies and eigenfunctions
around a certain energy. Figure 5 shows the spectra for the different studied
systems obtained using the predicted $\psi(\vec{x},t)$ with the multi-step
reservoir computing model. The peaks of the spectrum appear at the
eigenenergies of each quantum system. These eigenenergies are close to the
energy of the initial state $\psi(\vec{x},0)$ of the system. For example, the
initial state of the 1D harmonic oscillator is a minimum uncertainty Gaussian
wavepacket [see Eq. (18)] with mean energy $E=6$. Accordingly, the obtained
eigenenergies go from $E=2.5$ to $E=9.5$, so they are centered around the
initial energy $E=6$. Therefore, integrating the time-dependent Schrödinger
equation allows recovering the eigenstates around certain energy, without
needing to compute all the eigenstates with lower energy, as it happens in the
usual variational method. Using the eigenenergies, we recover their associated
eigenfunctions by computing the Fourier transform of the wavefunction
$\psi(\vec{x},t)$. Figures 6, 7, and 8 show the predicted and exact
eigenfunctions for the 1D harmonic oscillator, the Morse Hamiltonian and the
polynomial potential systems, and Fig. 9 shows the predicted and exact
eigenfunctions for the 2D harmonic oscillator. The exact eigenstates are
calculated analytically for the harmonic oscillator (both 1D and 2D) and the
Morse potential, and numerically (using the variational method Domingo and
Borondo (2021)) for the polynomial potential. We see that in all cases the
predicted eigenfunctions are in very good agreement with the exact ones. This
fact confirms that the multi-step RC method can correctly propagate a quantum
wavepacket with time. All in all, our results show that our adaptation of the
RC method allows to accurately integrate the time-dependent Schrödinger
equation, and then obtain, also with a high accuracy, the associated
eigenenegies and eigenfunctions. This is possible thanks to the fact that our
method helps prevent overfitting when working with high-dimensional data thus
allowing the use of the adapted RC method for time propagation in quantum
systems.

Figure 6: Exact and predicted eigenfunctions for the harmonic oscillator
system (17). Figure 7: Same as Fig. 6 for the Morse Hamiltonian system (19.)
Figure 8: Same as Fig. 6 for the polynomial Hamiltonian system (20).

Figure 9: Predicted (top row) and exact (bottom row) wavefunctions for the 2D
harmonic oscillator system of Eq. (22).

##  VII Conclusions and outlook

Despite being based on a quite simple training framework, RC has shown
remarkable performance in various benchmark tasks, such as time series
forecasting Gao  _et al._ (2021) and image recognition Jalalvand  _et al._
(2016). Moreover, Pathak et al. Pathak  _et al._ (2017, 2018) demonstrated
that ESN approaches are useful for chaotic time series prediction.

In this work, we have extended the RC framework to solve quantum problems. For
this purpose, we propagate wavepackets with a certain energy, and then
performed Fourier transform to obtain the eigenenergies around the initial
energy and their corresponding eigenfunctions. The use of RC computing to
integrate the time-dependent Schrödinger equation presented two main
challenges. The first one is that the data (the quantum states
$\psi(\vec{x},t)$) is complex-valued. To overcome this problem, we proposed to
use a complex-valued activation function. Also, we extended the regularized
linear model to complex data by extending the ridge regression to the complex
domain. The second challenge is derived from the high-dimensionality of the
input-output data. These data are usually represented as a matrix containing
the values of the wavefunctions in a spatial grid. The size of the matrix
increases exponentially with the dimension of the quantum system. Therefore,
we need large reservoirs in order to propagate the quantum states, and it can
easily happen that the size of the reservoir is similar or even larger than
the number of training samples. In this case, the model is likely to produce
overfitting, even when using a regularized model. When the model overfits the
training data, it is unable to generalize to unseen data, and errors in the
predictions propagate fast. To reduce the effect of this problem, we propose
the use of a multi-step learning algorithm to train the RC model. This
algorithm consists of splitting the training data into two sets. We first
train the reservoir with the first $t_{1}$ time steps. Then, we make
predictions, updating the internal states of the reservoir, with the second
part of the training data. In this way, we teach the reservoir how the
predictions affect the evolution of its internal states. Last, we retrain the
linear readout $W^{\text{out}}$. Since the reservoir has seen how predicting
modifies its internal states, small changes to the internal states due to
prediction errors will not propagate leading to large errors in the predicted
wavefunctions. This multi-step RC model is designed to work when the dimension
of the input/output is larger or comparable to the size of the reservoir.
Therefore, it could be useful to train RC models with any high-dimensional
data, other than quantum wavefunctions.

As an illustration, we have applied our method to four quantum systems: three
1D systems and one 2D system. It is observed that the MSE of the propagated
wavefunctions increases slower with time when using multi-step learning. This
fact is critical in the 2D system, which has higher-dimensional data, this
leading to more overfitting. In this case, the standard RC model was not able
to correctly reproduce the eigenenergies of the system, while the multi-step
learning could. Moreover, the multi-step learning RC also allowed to recover
the eigenfunctions of all the quantum systems, proving that the method can
correctly predict the time evolution of the wavepackets and the corresponding
eigenstates.

Once the efficiency of our multi-step RC has been proved, the present work can
be extended by application to other interesting problems. This future work
could include, among others, the application to more complex and realistic
quantum systems, the computation of eigenstates in a high lying energy window
Revuelta  _et al._ (2020), or the calculation of the so-called scarred
functions Revuelta  _et al._ (2013) that play a very important role in the
field of quantum chaos.

###### Acknowledgements.

The project that gave rise to these results received the support of a
fellowship from "la Caixa" Foundation (ID 100010434). The fellowship code is
LCF/BQ/DR20/11790028. This work has also been partially supported by the
Spanish Ministry of Science, Innovation and Universities, Gobierno de España,
under Contracts No. PGC2018-093854-BI00, ICMAT Severo Ochoa CEX2019-000904-S;
and by the People Programme (Marie Curie Actions) of the European Union’s
Horizon 2020 Research and Innovation Program under Grant No. 734557.

## Data Availability Statement

The data that support the findings of this study are openly available in
https://github.com/laiadc/RC_quantum.

## References

  * Rumelhart, Hinton, and Williams (1986) D. E. Rumelhart, G. E. Hinton, and R. J. Williams, “Learning internal representations by error propagation,” in  _Parallel Distributed Processing: Explorations in the Microstructure of Cognition, Volume 1: Foundations_ (MIT Press, Cambridge, MA, 1986) pp. 318–362.
  * Hochreiter and Schmidhuber (1997) S. Hochreiter and J. Schmidhuber, “Long short-term memory,” Neural Comput. 9, 1735 (1997).
  * Doya (2000) K. Doya, “Bifurcations in the learning of recurrent neural networks,” Proc. IEEE Int. Symp. Circuits and Systems 6, 2777 (2000).
  * Jaeger (2001) H. Jaeger, “The "echo state" approach to analysing and training recurrent neural networks-with an erratum note,” German National Research Center for Information Technology GMD Technical Report 148 (2001).
  * Jaeger (2007) H. Jaeger, “Echo state network,” Scholarpedia 2, 2330 (2007).
  * Maass, Natschläger, and Markram (2002) W. Maass, T. Natschläger, and H. Markram, “Real-time computing without stable states: a new framework for neural computation based on perturbations.” Neural Comput. 14, 2531 (2002).
  * Pathak  _et al._ (2017) J. Pathak, Z. Lu, B. Hunt, M. Girvan, and E. Ott, “Using machine learning to replicate chaotic attractors and calculate Lyapunov exponents from data,” Chaos 27, 121102 (2017).
  * Pathak  _et al._ (2018) J. Pathak, B. Hunt, M. Girvan, Z. Lu, and E. Ott, “Model-free prediction of large spatiotemporally chaotic systems from data: A reservoir computing approach,” Phys. Rev. Lett. 120, 024102 (2018).
  * Akiyama and Tanaka (2019) T. Akiyama and G. Tanaka, “Analysis on characteristics of multi-step learning echo state networks for nonlinear time series prediction,” in  _2019 International Joint Conference on Neural Networks (IJCNN)_ (2019) pp. 1–8.
  * Qiao  _et al._ (2017) J. Qiao, F. Li, H. Han, and W. Li, “Growing echo-state network with multiple subreservoirs,” IEEE Trans. Neural Netw. Learn. Syst. 28, 391 (2017).
  * Chen  _et al._ (2020) P. Chen, R. Liu, K. Aihara, and L. Chen, “Autoreservoir computing for multistep ahead prediction based on the spatiotemporal information transformation,” Nat. Commun. 11, 4568 (2020).
  * Gauthier  _et al._ (2021) D. Gauthier, E. Bollt, A. Griffith, and W. Barbosa, “Next generation reservoir computing,” Nat. Commun. 12, 5564 (2021).
  * Ferguson  _et al._ (2020) A. Ferguson, J. Hachmann, T. Miller, and J. Pfaendtner, “Virtual special issue on machine learning in physical chemistry,” Journal of Physical Chemistry B 124, 9767 (2020).
  * Pfau  _et al._ (2020) D. Pfau, J. Spencer, A. Matthews, and W. Foulkes, “Ab initio solution of the many-electron Schrödinger equation with deep neural networks,” Phys. Rev. Res. 2, 033429 (2020).
  * Hermann, Schätzle, and Noé (2020) J. Hermann, Z. Schätzle, and F. Noé, “Deep-neural-network solution of the electronic Schrödinger equation,” Nat. Chem. 12, 891 (2020).
  * Manzhos, Yamashita, and Carrington (2009) S. Manzhos, K. Yamashita, and T. Carrington, “Using a neural network based method to solve the vibrational Schrödinger equation for H2O,” Chem. Phys. Lett. 474, 217 – 221 (2009).
  * Pavlov, Serdyuk, and Ustinov (2019) A. Pavlov, J. Serdyuk, and A. Ustinov, “Machine learning and the Schrödinger equation,” J. Phys.: Conf. Series 1236, 012050 (2019).
  * Mills, Spanner, and Tamblyn (2017) K. Mills, M. Spanner, and I. Tamblyn, “Deep learning and the Schrödinger equation,” Phys. Rev. A 96, 042113 (2017).
  * Domingo and Borondo (2021) L. Domingo and F. Borondo, “Deep learning methods for the computation of vibrational wavefunctions,” Commun. Nonlinear Sci. Numer. Simul. 103, 105989 (2021).
  * Mohri, Rostamizadeh, and Talwalkar (2018) M. Mohri, A. Rostamizadeh, and A. Talwalkar,  _Foundations of Machine Learning_ , 2nd ed. (The MIT Press, 2018).
  * Kosloff and Kosloff (1983) D. Kosloff and R. Kosloff, “A Fourier method solution for the time dependent Schrödinger equation as a tool in molecular dynamics,” J. Comput. Phys. 52, 35 (1983).
  * Zdánská and Moiseyev (2004) P. Zdánská and N. Moiseyev, “Complex autocorrelation function and energy spectrum by classical trajectory calculations,” J. Chem. Phys. 121, 6175 (2004).
  * Gao  _et al._ (2021) R. Gao, L. Du, O. Duru, and K. F. Yuen, “Time series forecasting based on echo state network and empirical wavelet transformation,” Appl. Soft Comput. 102, 107111 (2021).
  * Jalalvand  _et al._ (2016) A. Jalalvand, W. De Neve, R. Van de Walle, and J.-P. Martens, “Towards using reservoir computing networks for noise-robust image recognition,” in  _2016 International Joint Conference on Neural Networks (IJCNN)_ (2016) pp. 1666–1672.
  * Revuelta  _et al._ (2020) F. Revuelta, E. Vergini, R. M. Benito, and F. Borondo, “Short-periodic-orbit method for excited chaotic eigenfunctions,” Phys. Rev. E 102 (2020).
  * Revuelta  _et al._ (2013) F. Revuelta, R. Benito, F. Borondo, and E. Vergini, “Using basis sets of scar functions,” Phys. Rev. E 87, 042921 (2013).