# Quantum reservoir computing: a reservoir approach toward quantum machine
learning on near-term quantum devices

Keisuke Fujii<EMAIL_ADDRESS>Graduate School of Engineering
Science, Osaka University, 1-3 Machikaneyama, Toyonaka, Osaka 560-8531, Japan.
Kohei Nakajima<EMAIL_ADDRESS>Graduate School of Information
Science and Technology, The University of Tokyo, Bunkyo-ku, 113-8656 Tokyo,
Japan

###### Abstract

Quantum systems have an exponentially large degree of freedom in the number of
particles and hence provide a rich dynamics that could not be simulated on
conventional computers. Quantum reservoir computing is an approach to use such
a complex and rich dynamics on the quantum systems as it is for temporal
machine learning. In this chapter, we explain quantum reservoir computing and
related approaches, quantum extreme learning machine and quantum circuit
learning, starting from a pedagogical introduction to quantum mechanics and
machine learning. All these quantum machine learning approaches are
experimentally feasible and effective on the state-of-the-art quantum devices.

quantum reservoir computing, quantum circuit learning, quantum extreme
learning machine, quantum machine learning, noisy intermediate-scale quantum
technology

##  I Introduction

Over the past several decades, we have enjoyed exponential growth of
computational power, namely, Moore’s law. Nowadays even smart phone or tablet
PC is much more powerful than super computers in 1980s. Even though, people
are still seeking more computational power, especially for artificial
intelligence (machine learning), chemical and material simulations, and
forecasting complex phenomena like economics, weather and climate. In addition
to improving computational power of conventional computers, i.e., more Moore’s
law, a new generation of computing paradigm has been started to be
investigated to go beyond Moore’s law. Among them, natural computing seeks to
exploit natural physical or biological systems as computational resource.
Quantum reservoir computing is an intersection of two different paradigms of
natural computing, namely, quantum computing and reservoir computing.

Regarding quantum computing, the recent rapid experimental progress in
controlling complex quantum systems motivates us to use quantum mechanical law
as a new principle of information processing, namely, quantum information
processing NielsenChuang ; FujiiText . For example, certain mathematical
problems, such as integer factorisation, which are believed to be intractable
on a classical computer, are known to be efficiently solvable by a
sophisticatedly synthesized quantum algorithm ShorFact . Therefore,
considerable experimental effort has been devoted to realising full-fledged
universal quantum computers UCSB1 ; UCSB2 . In the near feature, quantum
computers of size $>50$ qubits with fidelity $>99\%$ for each elementary gate
would appear to achieve quantum computational supreamcy beating simulation on
the-state-of-the-art classical supercomputers Preskill ; Boixo . While this
does not directly mean that a quantum computer outperforms classical computers
for a useful task like machine learning, now applications of such a near-term
quantum device for useful tasks including machine leanring has been widely
explored. On the other hand, quantum simulators are thought to be much easier
to implement than a full-fledged universal quantum computer. In this regard,
existing quantum simulators have already shed new light on the physics of
complex many-body quantum systems QuantumSimulator1 ; QuantumSimulator2 ;
QuantumSimulator3 , and a restricted class of quantum dynamics, known as
adiabatic dynamics, has also been applied to combinatorial optimisation
problems Nishimori ; Farhi ; Dwave ; Dwave2 . However, complex real-time
quantum dynamics, which is one of the most difficult tasks for classical
computers to simulate FujiiMorimae ; Fujiietal ; FujiiTamate  and has great
potential to perform nontrivial information processing, is now waiting to be
harnessed as a resource for more general purpose information processing.

Physical reservoir computing, which is the main subject throughout this book,
is another paradigm for exploiting complex physical systems for information
processing. In this framework, the low-dimensional input is projected to a
high-dimensional dynamical system, which is typically referred to as a
reservoir, generating transient dynamics that facilitates the separation of
input states Transient . If the dynamics of the reservoir involve both
adequate memory and nonlinearity Capacity , emulating nonlinear dynamical
systems only requires adding a linear and static readout from the high-
dimensional state space of the reservoir. A number of different
implementations of reservoirs have been proposed, such as abstract dynamical
systems for echo state networks (ESNs) Jaeger0  or models of neurons for
liquid state machines Maass0 . The implementations are not limited to programs
running on the PC but also include physical systems, such as the surface of
water in a laminar state Bucket , analogue circuits and optoelectronic systems
Laser0 ; Laser1 ; Laser1b ; Laser2 ; Laser3 ; Laser4 , and neuromorphic chips
Neuromorphic0 . Recently, it has been reported that the mechanical bodies of
soft and compliant robots have also been successfully used as a reservoir
Helmut0 ; Kohei0 ; Kohei1 ; Kohei2 ; Kohei3 ; Ken1 . In contrast to the
refinements required by learning algorithms, such as in deep learning DL , the
approach followed by reservoir computing, especially when applied to real
systems, is to find an appropriate form of physics that exhibits rich
dynamics, thereby allowing us to outsource a part of the computation.

Quantum reservoir computing (QRC) was born in the marriage of quantum
computing and physical reservoir computing above to harness complex quantum
dynamics as a reservoir for real-time machine learning tasks QRC . Since the
idea of QRC has been proposed in Ref. QRC , its proof-of-principle
experimental demonstration for non temporal tasks QRC_ex  and performance
analysis and improvement QRC_multiplex ; QRC_OTOC ; QRC_Higher_Order  has been
explored. The QRC approach to quantum tasks such as quantum tomography and
quantum state preparation has been recently garnering attention QRC_sig ;
QRC_sig1 ; QRC_sig2 . In this book chapter, we will provide a broad picture of
QRC and related approaches starting from a pedagogical introduction to quantum
mechanics and machine learning.

The rest of this paper is organized as follows. In Sec II, we will provide a
pedagogical introduction to quantum mechanics for those who are not familiar
to it and fix our notation. In Sec III, we will briefly mention to several
machine learning techniques like, linear and nonlinear regressions, temporal
machine learning tasks and reservoir computing. In Sec IV, we will explain QRC
and related approaches, quantum extreme learning machine QRC_ex  and quantum
circuit learning QCL . The former is a framework to use quantum reservoir for
non temporal tasks, that is, the input is fed into a quantum system, and
generalization or classification tasks are performed by a linear regression on
a quantum enhanced feature space. In the latter, the parameters of the quantum
system is further fine-tuned via the gradient descent by measuring an
analytically obtained gradient, just like the back propagation for feedforward
neural networks. Regarding QRC, we will also see chaotic time series
predictions as demonstrations. Sec. V is devoted to conclusion and discussion.

##  II Pedagogical introduction to quantum mechanics

In this section, we would like to provide a pedagogical introduction to how
quantum mechanical systems work for those who are not familiar to quantum
mechanics. If you already familiar to quantum mechanics and its notations,
please skip to Sec. III.

###  II.1 Quantum state

A state of a quantum system is described by a state vector,

$\displaystyle|\psi\rangle=\left(\begin{array}[]{c}c_{1}\\\ \vdots\\\
c_{d}\end{array}\right)$ (4)

on a complex $d$-dimensional system $\mathbb{C}^{d}$, where the symbol
$|\cdot\rangle$ is called ket and indicates a complex column vector.
Similarly, $\langle\cdot|$ is called bra and indicates a complex row vector,
and they are related complex conjugate,

$\displaystyle\langle\psi|=|\psi\rangle^{{\dagger}}=\left(\begin{array}[]{ccc}c^{*}_{1}&\cdots&c^{*}_{d}\end{array}\right).$
(6)

With this notation, we can writte an inner product of two quantum state
$|\psi\rangle$ and $|\phi\rangle$ by $\langle\psi|\phi\rangle$. Let us define
an orthogonal basis

$\displaystyle|1\rangle=\left(\begin{array}[]{c}1\\\ 0\\\ \vdots\\\ \vdots\\\
0\end{array}\right),...\;\;\;|k\rangle=\left(\begin{array}[]{c}0\\\ \vdots\\\
1\\\ 0\\\
\vdots\end{array}\right),...\;\;\;|d\rangle=\left(\begin{array}[]{c}0\\\
\vdots\\\ \vdots\\\ \vdots\\\ d\end{array}\right),$ (22)

a quantum state in the $d$-dimensional system can be described simply by

$\displaystyle|\psi\rangle=\sum_{i=1}^{d}c_{i}|i\rangle.$ (23)

The state is said to be a superposition state of $|i\rangle$. The coefficients
$\\{c_{i}\\}$ are complex, and called complex probability amplitudes. If we
measure the system in the basis $\\{|i\rangle\\}$, we obtain the measurement
outcome $i$ with a probability

$\displaystyle p_{i}=|\langle i|\psi\rangle|^{2}=|c_{i}|^{2},$ (24)

and hence the complex probability amplitudes have to be normalized as follows

$\displaystyle|\langle\psi|\psi\rangle|^{2}=\sum_{i=1}^{d}|c_{i}|^{2}=1.$ (25)

In other words, a quantum state is represented as a normalized vector on a
complex vector space.

Suppose the measurement outcome $i$ corresponds to a certain physical value
$a_{i}$, like energy, magnetization and so on, then the expectation value of
the physical valuable is given by

$\displaystyle\sum_{i}a_{i}p_{i}=\langle\psi|A|\psi\rangle\equiv\langle
A\rangle,$ (26)

where we define an hermitian operator

$\displaystyle A=\sum_{i}a_{i}|i\rangle\langle i|,$ (27)

which is called observable, and has the information of the measurement basis
and physical valuable.

The state vector in quantum mechanics is similar to a probability
distribution, but essentially different form it, since it is much more
primitive; it can take complex value and is more like a square root of a
probability. The unique features of the quantum systems come from this
property.

###  II.2 Time evolution

The time evolution of a quantum system is determined by a Hamiltonian $H$,
which is a hermitian operator acting on the system. Let us denote a quantum
state at time $t=0$ by $|\psi(0)\rangle$. The equation of motion for quantum
mechanics, so-called Schrödinger equation, is given by

$\displaystyle i\frac{\partial}{\partial t}|\psi(t)\rangle=H|\psi(t)\rangle.$
(28)

This equation can be formally solved by

$\displaystyle|\psi(t)\rangle=e^{-iHt}|\psi(0)\rangle.$ (29)

Therefore the time evolution is given by an operator $e^{-iHt}$, which is a
unitary operator and hence the norm of the state vector is preserved, meaning
the probability conservation. In general, the Hamiltonian can be time
dependent. Regarding the time evolution, if you are not interested in the
continuous time evolution, but in just its input and output relation, then the
time evolution is nothing but a unitary operator $U$

$\displaystyle|\psi_{\rm out}\rangle=U|\psi_{\rm in}\rangle.$ (30)

In quantum computing, the time evolution $U$ is sometimes called quantum gate.

###  II.3 Qubits

The smallest nontrivial quantum system is a two-dimensional quantum system
$\mathbb{C}^{2}$, which is called quantum bit or qubit:

$\displaystyle\alpha|0\rangle+\beta|1\rangle,\;\;\;(|\alpha|^{2}+|\beta|^{2}=1).$
(31)

Suppose we have $n$ qubits. The $n$-qubit system is defined by a tensor
product space $(\mathbb{C}^{2})^{\otimes n}$ of each two-dimensional system as
follows. A basis of the system is defined by a direct product of a binary
state $|x_{k}\rangle$ with $x_{k}\in\\{0,1\\}$,

$\displaystyle|x_{1}\rangle\otimes|x_{2}\rangle\otimes\cdots\otimes|x_{n}\rangle,$
(32)

which is simply denoted by

$\displaystyle|x_{1}x_{2}\cdots x_{n}\rangle.$ (33)

Then a state of the $n$-qubit system can be described as

$\displaystyle|\psi\rangle=\sum_{x_{1},x_{2},...,x_{n}}\alpha_{x_{1},x_{2},...,x_{n}}|x_{1}x_{2}\cdots
x_{n}\rangle.$ (34)

The dimension of the $n$-qubit system is $2^{n}$, and hence the tensor product
space is nothing but a $2^{n}$-dimensional complex vector space
$\mathbb{C}^{2^{n}}$. The dimension of the $n$-qubit system increases
exponentially in the number $n$ of the qubits.

###  II.4 Density operator

Next, I would like to introduce operator formalism of the above quantum
mechanics. This describes an exactly the same thing but sometimes the operator
formalism would be convenient. Let us consider an operator $\rho$ constructed
from the state vector $|\psi\rangle$:

$\displaystyle\rho=|\psi\rangle\langle\psi|.$ (35)

If you chose the basis of the system $\\{|i\rangle\\}$ for the matrix
representation, then the diagonal elements of $\rho$ corresponds the
probability distribution $p_{i}=|c_{i}|^{2}$ when the system is measured in
the basis $\\{|i\rangle\\}$. Therefore the operator $\rho$ is called a density
operator. The probability distribution can also be given in terms of $\rho$ by

$\displaystyle p_{i}={\rm Tr}[|i\rangle\langle i|\rho],$ (36)

where ${\rm Tr}$ is the matrix trace. An expectation value of an observable
$A$ is given by

$\displaystyle\langle A\rangle={\rm Tr}[A\rho].$ (37)

The density operator can handle a more general situation where a quantum state
is sampled form a set of quantum states $\\{|\psi_{k}\rangle\\}$ with a
probability distribution $\\{q_{k}\\}$. In this case, if we measure the system
in the basis $\\{|i\rangle\langle i|\\}$, the probability to obtain the
measurement outcome $i$ is given by

$\displaystyle p_{i}=\sum_{k}q_{k}{\rm Tr}[|i\rangle\langle i|\rho_{k}],$ (38)

where $\rho_{k}=|\psi_{k}\rangle\langle\psi_{k}|$. By using linearity of the
trace function, this reads

$\displaystyle p_{i}={\rm Tr}[|i\rangle\langle i|\sum_{k}q_{k}\rho_{k}].$ (39)

Now we interpret that the density operator is given by

$\displaystyle\rho=\sum_{k}q_{k}|\psi_{k}\rangle\langle\psi_{k}|.$ (40)

In this way, a density operator can represent classical mixture of quantum
states by a convex mixture of density operators, which is convenient in many
cases. In general, a positive and hermitian operator $\rho$ being subject to
${\rm Tr}[\rho]=1$ can be a density operator, since it can be interpreted as a
convex mixture of quantum states via spectral decomposition:

$\displaystyle\rho=\sum\lambda_{i}|\lambda_{i}\rangle\langle\lambda_{i}|,$
(41)

where $\\{|\lambda_{i}\rangle\\}$ and $\\{\lambda_{i}\\}$ are the eigenstates
and eigenvectors respectively. Because of ${\rm Tr}[\rho]=1$, we have
$\sum_{i}\lambda_{i}=1$.

From its definition, the time evolution of $\rho$ can be given by

$\displaystyle\rho(t)=e^{-iHt}\rho(0)e^{iHt}$ (42)

or

$\displaystyle\rho_{\rm out}=U\rho_{\rm in}U^{{\dagger}}.$ (43)

Moreover, we can define more general operations for the density operators. For
example, if we apply unitary operators $U$ and $V$ with probabilities $p$ and
$(1-p)$, respectively, then we have

$\displaystyle\rho_{\rm out}=pU\rho U^{{\dagger}}+(1-p)V\rho V^{{\dagger}}.$
(44)

As another example, if we perform the measurement of $\rho$ in the basis
$\\{|i\rangle\\}$, and we forget about the measurement outcome, then the state
is now given by a density operator

$\displaystyle\sum_{i}{\rm Tr}[|i\rangle\langle i|\rho]|i\rangle\langle
i|=\sum_{i}|i\rangle\langle i|\rho|i\rangle\langle i|.$ (45)

Therefore if we define a map from a density operator to another, which we call
superoperator,

$\displaystyle\mathcal{M}(\cdots)=\sum_{i}|i\rangle\langle
i|(\cdots)|i\rangle\langle i|,$ (46)

the above non-selective measurement (forgetting about the measurement
outcomes) is simply written by

$\displaystyle\mathcal{M}(\rho).$ (47)

In general, any physically allowed quantum operation $\mathcal{K}$ that maps a
density operator to another can be represented in terms of a set of operators
$\\{K_{i}\\}$ being subject to $K_{i}^{{\dagger}}K_{i}=I$ with an identity
operator $I$:

$\displaystyle\mathcal{K}(\rho)=\sum_{i}K_{i}\rho K^{{\dagger}}_{i}.$ (48)

The operators $\\{K_{i}\\}$ are called Kraus operators.

###  II.5 Vector representation of density operators

Finally, we would like to introduce a vector representation of the above
operator formalism. The operators themselves satisfy axioms of the linear
space. Moreover, we can also define an inner product for two operators, so-
called Hilbert-Schmidt inner product, by

$\displaystyle{\rm Tr}[A^{{\dagger}}B].$ (49)

The operators on the $n$-qubit system can be spanned by the tensor product of
Pauli operators $\\{I,X,Y,Z\\}^{\otimes n}$,

$\displaystyle
P(\mbox{\boldmath$i$})=\bigotimes_{k=1}^{n}\sigma_{i_{2k-1}i_{2k}}.$ (50)

where $\sigma_{ij}$ is the Pauli operators:

$\displaystyle I=\sigma_{00}=\left(\begin{array}[]{cc}1&0\\\
0&1\end{array}\right),X=\sigma_{10}=\left(\begin{array}[]{cc}0&1\\\
1&0\end{array}\right),$ $\displaystyle
Z=\sigma_{01}=\left(\begin{array}[]{cc}1&0\\\
0&-1\end{array}\right),Y=\sigma_{11}=\left(\begin{array}[]{cc}0&-i\\\
i&0\end{array}\right).$ (56)

Since the Pauli operators constitute a complete basis on the operator space,
any operator $A$ can be decomposed into a linear combination of
$P(\mbox{\boldmath$i$})$,

$\displaystyle
A=\sum_{\mbox{\boldmath$i$}}a_{\mbox{\boldmath$i$}}P(\mbox{\boldmath$i$}).$
(57)

The coefficient $a_{\mbox{\boldmath$i$}}$ can be calculated by using the
Hilbert-Schmidt inner product as follows:

$\displaystyle a_{\mbox{\boldmath$i$}}={\rm
Tr}[P(\mbox{\boldmath$i$})A]/2^{n},$ (58)

by virtue of the orthogonality

$\displaystyle{\rm
Tr}[P(\mbox{\boldmath$i$})P(\mbox{\boldmath$j$})]/2^{n}=\delta_{\mbox{\boldmath$i$},\mbox{\boldmath$j$}}.$
(59)

The number of the $n$-qubit Pauli operators $\\{P(\mbox{\boldmath$i$})\\}$ is
$4^{n}$, and hence a density operator $\rho$ of the $n$-qubit system can be
represented as a $4^{n}$-dimensional vector

$\displaystyle\mbox{\boldmath$r$}=\left(\begin{array}[]{c}r_{00...0}\\\
\vdots\\\ r_{11...1}\end{array}\right),$ (63)

where $r_{00...0}=1/2^{n}$ because of ${\rm Tr}[\rho]=1$. Moreover, because
$P(\mbox{\boldmath$i$})$ is hermitian, $r$ is a real vector. The superoperator
$\mathcal{K}$ is a linear map for the operator, and hence can be represented
as a matrix acting on the vector $r$:

$\displaystyle\rho^{\prime}=\mathcal{K}(\rho)\Leftrightarrow\mbox{\boldmath$r$}^{\prime}=K\mbox{\boldmath$r$},$
(64)

where the matrix element is given by

$\displaystyle K_{\mbox{\boldmath$i$}\mbox{\boldmath$j$}}={\rm
Tr}[P(\mbox{\boldmath$i$})\mathcal{K}\left(P(\mbox{\boldmath$j$})\right)]/2^{n}.$
(65)

In this way, a density operator $\rho$ and a quantum operation $\mathcal{K}$
on it can be represented by a vector $r$ and a matrix $K$, respectively.

##  III Machine learning and reservoir approach

In this section, we briefly introduce machine learning and reservoir
approaches.

###  III.1 Linear and nonlinear regression

A supervised machine learning is a task to construct a model $f(x)$ from a
given set of teacher data $\\{x^{(j)},y^{(j)}\\}$ and to predict the output of
an unknown input $x$. Suppose $x$ is a $d$-dimensional data, and $f(x)$ is one
dimensional, for simplicity. The simplest model is linear regression, which
models $f(x)$ as a linear function with respect to the input:

$\displaystyle f(x)=\sum_{i=1}^{d}w_{i}x_{i}+w_{0}.$ (66)

The weights $\\{w_{i}\\}$ and bias $w_{0}$ are chosen such that an error
between $f(x)$ and the output of the teacher data, i.e. loss, becomes minimum.
If we employ a quadratic loss function for given teacher data
$\\{\\{x^{(j)}_{i}\\},y^{(j)}\\}$, the problem we have to solve is as follows:

$\displaystyle\min_{\\{w_{i}\\}}\sum_{j}(\sum_{i=0}^{d}w_{i}x^{(j)}_{i}-y^{(j)})^{2},$
(67)

where we introduced a constant node $x_{0}=1$. This corresponds to solving a
superimposing equations:

$\displaystyle\mathbf{y}=\mathbf{X}\mathbf{w},$ (68)

where $\mathbf{y}_{j}=y^{(j)}$, $\mathbf{X}_{ji}=x^{(j)}_{i}$, and
$\mathbf{w}_{i}=w_{i}$. This can be solved by using the Moore-Penrose pseudo
inverse $\mathbf{X}^{+}$, which can be defined from the singular value
decomposition of $\mathbf{X}=UDV^{T}$ to be

$\displaystyle\mathbf{X}^{+}=VDU^{T}.$ (69)

Unfortunately, the linear regression results in a poor performance in
complicated machine learning tasks, and any kind of nonlinearity is
essentially required in the model. A neural network is a way to introduce
nonlinearity to the model, which is inspired by the human brain. In the neural
network, the $d$-dimensional input data $x$ is fed into $N$-dimensional hidden
nodes with an $N\times d$ input matrix $W^{\rm in}$:

$\displaystyle W^{\rm in}x.$ (70)

Then each element of the hidden nodes is now processed by a nonlinear
activation function $\sigma$ such as $\tanh$, which is denoted by

$\displaystyle\sigma(W^{\rm in}x).$ (71)

Finally the output is extracted by an output weight $W^{\rm out}$ ($1\times N$
dimensional matrix):

$\displaystyle W^{\rm out}\sigma(W^{\rm in}x).$ (72)

The parameters in $W^{\rm in}$ and $W^{\rm out}$ are trained such that the
error between the output and teacher data becomes minimum. While this
optimization problem is highly nonlinear, a gradient based optimization, so-
called back propagation, can be employed. To improve a representation power of
the model, we can concatenate the linear transformation and the activation
function as follows:

$\displaystyle W^{\rm
out}\sigma\left(W^{(l)}\cdots\sigma\left(W^{(1)}\sigma(W^{\rm
in}x)\right)\right),$ (73)

which is called multi-layer perceptron or deep neural network.

###  III.2 Temporal task

The above task is not a temporal task, meaning that the input data is not
sequential but given simultaneously like the recognition task of images for
hand written language, pictures and so on. However, for a recognition of
spoken language or prediction of time series like stock market, which are
called temporal tasks, the network has to handle the input data that is given
in a sequential way. To do so, the recurrent neural network feeds the previous
states of the nodes back into the states of the nodes at next step, which
allows the network to memorize the past input. In contrast, the neural network
without any recurrency is called a feedforward neural network.

Let us formalize a temporal machine learning task with the recurrent neural
network. For given input time series $\\{x_{k}\\}_{k=1}^{L}$ and target time
series $\\{\bar{y}_{k}\\}_{k=1}^{L}$, a temporal machine learning is a task to
generalize a nonlinear function,

$\displaystyle\bar{y}_{k}=f(\\{x_{j}\\}_{j=1}^{k}).$ (74)

For simplicity, we consider one-dimensional input and output time series, but
their generalization to a multi-dimensional case is straightforward. To learn
the nonlinear function $f(\\{x_{j}\\}_{j=1}^{k})$, the recurrent neural
network can be employed as a model. Suppose the recurrent neural network
consists of $m$ nodes and is denoted by $m$-dimensional vector

$\displaystyle\mbox{\boldmath$r$}=\left(\begin{array}[]{c}r_{1}\\\ \vdots\\\
r_{m}\end{array}\right).$ (78)

To process the input time series, the nodes evolve by

$\displaystyle\mbox{\boldmath$r$}(k+1)=\sigma[W\mbox{\boldmath$r$}(k)+W^{\rm
in}x_{k}],$ (79)

where $W$ is an $m\times m$ transition matrix and $W^{\rm in}$ is an $m\times
1$ input weight matrix. Nonlinearity comes from the nonlinear function
$\sigma$ applied on each element of the nodes. The output time series from the
network is defined in terms of a $1\times m$ readout weights by

$\displaystyle y_{k}=W^{\rm out}\mbox{\boldmath$r$}(k).$ (80)

Then the learning task is to determine the parameters in $W^{\rm in}$, $W$,
and $W^{\rm out}$ by using the teacher data
$\\{x_{k},\bar{y}_{k}\\}_{k=1}^{L}$ so as to minimize an error between the
teacher $\\{\bar{y}_{k}\\}$ and the output $\\{y_{k}\\}$ of the network.

###  III.3 Reservoir approach

While the representation power of the recurrent neural network can be improved
by increasing the number of the nodes, it makes the optimization process of
the weights hard and unstable. Specifically, the back propagation based
methods always suffer from the vanishing gradient problem. The idea of
reservoir computing is to resolve this problem by mapping an input into a
complex higher dimensional feature space, i.e., reservoir, and by performing
simple linear regression on it.

Let us first see a reservoir approach on a feedforward neural network, which
is called extreme learning machine ELM . The input data $x$ is fed into a
network like multi-layer perceptron, where all weights are chosen randomly.
The states of the hidden nodes at some layer is now regarded as basis
functions of the input $x$ in the feature space:

$\displaystyle\\{\phi_{1}(x),\phi_{2}(x),...,\phi_{N}(x)\\}.$ (81)

Now the output is defined as a linear combination of these

$\displaystyle\sum_{i}w_{i}\phi_{i}(x)+w_{0}$ (82)

and hence the coefficients are determined simply by the linear regression as
mentioned before. If the dimension and nonlinearity of the the basis functions
are high enough, we can model a complex task simply by the linear regression.

The echo state network is similar but employs the reservoir idea for the
recurrent neural network Jaeger0 ; Maass0 ; Reservoir , which has been
proposed before extreme learning machine appeared. To be specific, the input
weights $W^{\rm in}$ and weight matrix $W$ are both chosen randomly up to an
appropriate normalization. Then the learning task is done by finding the
readout weights $W^{\rm out}$ to minimize the mean square error

$\displaystyle\sum_{k}(y_{k}-\bar{y}_{k})^{2}.$ (83)

This problem can be solved stably by using the pseudo inverse as we mentioned
before.

For both feedforward and recurrent types, the reservoir approach does not need
to tune the internal parameters of the network depending on the tasks as long
as it posses sufficient complexity. Therefore, the system, to which the
machine learning tasks are outsourced, is not necessarily the neural network
anymore, but any nonlinear physical system of large degree of freedoms can be
employed as a reservoir for information processing, namely, physical reservoir
computing Bucket ; Laser0 ; Laser1 ; Laser1b ; Laser2 ; Laser3 ; Laser4 ;
Neuromorphic0 ; Helmut0 ; Kohei0 ; Kohei1 ; Kohei2 ; Kohei3 ; Ken1 .

##  IV Quantum machine learning on near-term quantum devices

In this section, we will see QRC and related frameworks for quantum machine
learning. Before going deep into the temporal tasks done on QRC, we first
explain how complicated quantum natural dynamics can be exploit as
generalization and classification tasks. This can be viewed as a quantum
version of extreme learning machine QRC_ex . While it is an opposite direction
to reservoir computing, we will also see quantum circuit learning (QCL) QCL ,
where the parameters in the complex dynamics is further tuned in addition to
the linear readout weights. QCL is a quantum version of a feedforward neural
network. Finally, we will explain quantum reservoir computing by extending
quantum extreme learning machine for temporal learning tasks.

###  IV.1 Quantum extreme learning machine

The idea of quantum extreme learning machine lies in using a Hilbert space,
where quantum states live, as an enhanced feature space of the input data. Let
us denote the set of input and teacher data by $\\{x^{(j)},\bar{y}^{(j)}\\}$.
Suppose we have an $n$-qubit system, which is initialized to

$\displaystyle|0\rangle^{\otimes n}.$ (84)

In order to feed the input data into quantum system, a unitary operation
parameterized by $x$, say $V(x)$, is applied on the initial state:

$\displaystyle V(x)|0\rangle^{\otimes n}.$ (85)

For example, if $x$ is one-dimensional data and normalized to be $0\leq x\leq
1$, then we may employ the $Y$-basis rotation $e^{-i\theta Y}$ with an angle
$\theta=\arccos(\sqrt{x})$:

$\displaystyle e^{-i\theta Y}|0\rangle=\sqrt{x}|0\rangle+\sqrt{1-x}|1\rangle.$
(86)

The expectation value of $Z$ with respect to $e^{-i\theta Y}|0\rangle$ becomes

$\displaystyle\langle Z\rangle=2x-1,$ (87)

and hence is linearly related to the input $x$. To enhance the power of
quantum enhanced feature space, the input could be transformed by using a
nonlinear function $\phi$:

$\displaystyle\theta=\arccos(\sqrt{\phi(x)}).$ (88)

The nonlinear function $\phi$ could be, for example, hyperbolic tangent,
Legendre polynomial, and so on. For simplicity, below we will use the simple
linear input $\theta=\arccos(\sqrt{x})$.

If we apply the same operation on each of the $n$ qubits, we have

$\displaystyle V(x)|0\rangle^{\otimes n}$ $\displaystyle=$
$\displaystyle(\sqrt{x}|0\rangle+\sqrt{1-x}|1\rangle)^{\otimes n}$
$\displaystyle=$
$\displaystyle(1-x)^{n/2}\sum_{i_{1},...,i_{n}}\prod_{k}\sqrt{\frac{x}{1-x}}^{i_{k}}|i_{1},...,i_{n}\rangle.$

Therefore, we have coefficients that are nonlinear with respect to the input
$x$ because of the tensor product structure. Still the expectation value of
the single qubit operator $Z_{k}$ on the $k$th qubit is $2x-1$. However, if we
measure a correlated operator like $Z_{1}Z_{2}$, we can obtain a second order
nonlinear output

$\displaystyle\langle Z_{1}Z_{2}\rangle=(2x-1)^{2}$ (89)

with respect to the input $x$. To measure a correlated operator, it is enough
to apply an entangling unitary operation like CNOT gate
$\Lambda(X)=|0\rangle\langle 0|\otimes I+|1\rangle\langle 1|\otimes X$:

$\displaystyle\langle\psi|\Lambda_{1,2}(X)Z_{1}\Lambda_{1,2}(X)|\psi\rangle=\langle\psi|Z_{1}Z_{2}|\psi\rangle.$
(90)

In general, an $n$-qubit unitary operation $U$ transforms the observable $Z$
under the conjugation into a linear combination of Pauli operators:

$\displaystyle
U^{{\dagger}}Z_{1}U=\sum_{\mbox{\boldmath$i$}}\alpha_{\mbox{\boldmath$i$}}P(\mbox{\boldmath$i$}).$
(91)

Thus if you measure the output of the quantum circuit after applying a unitary
operation $U$,

$\displaystyle UV(x)|0\rangle^{\otimes n},$ (92)

you can get a complex nonlinear output, which could be represented as a linear
combination of exponentially many nonlinear functions. $U$ should be chosen to
be appropriately complex with keeping experimental feasibility but not
necessarily fine-tuned.

Figure 1: The expectation value $\langle Z\rangle$ of the output of a quantum
circuit as a function of the input $(x_{0},x_{1})$.

To see how the output behaves in a nonlinear way with respect to the input, in
Fig. 1, we will plot the output $\langle Z\rangle$ for the input
$(x_{0},x_{1})$ and $n=8$, where the inputs are fed into the quantum state by
the $Y$-rotation with angles

$\displaystyle\theta_{2k}=k\arccos(\sqrt{x_{0}})$ (93)
$\displaystyle\theta_{2k+1}=k\arccos(\sqrt{x_{1}})$ (94)

on the $2k$th and $(2k+1)$th qubits, respectively. Regarding the unitary
operation $U$, random two-qubit gates are sequentially applied on any pairs of
two qubits on the 8-qubit system.

Figure 2: (a) The quantum circuit for quantum extreme learning machine. The
box with $theta_{k}$ indicates $Y$-rotations by angles $\theta_{k}$. The red
and blue boxes correspond to $X$ and $Z$ rotations by random angles, Each
dotted-line box represent a two-qubit gate consisting of two controlled-$Z$
gates and 8 $X$-rotations and 4 $Z$-rotations. As denoted by the dashed-line
box, the sequence of the 7 dotted boxes is repeated twice. The readout is
defined by a linear combination of $\langle Z_{i}\rangle$ with constant bias
term 1.0 and the input $(x_{0},x_{1})$. (b) (Left) The training data for a
two-class classification problem. (Middle) The readout after learning. (Right)
Prediction from the readout with threshold at 0.5.

Suppose the Pauli $Z$ operator is measured on each qubit as an observable.
Then we have

$\displaystyle z_{i}=\langle Z_{i}\rangle,$ (95)

for each qubit. In quantum extreme learning machine, the output is defined by
taking linear combination of these $n$ output:

$\displaystyle y=\sum_{i=1}^{n}w_{i}z_{i}.$ (96)

Now the linear readout weights $\\{w_{i}\\}$ are tuned so that the quadratic
loss function

$\displaystyle L=\sum_{j}(y^{(j)}-\bar{y}^{(j)})^{2}$ (97)

becomes minimum. As we mentioned previously, this can be solved by using the
pseudo inverse. In short, quantum extreme learning machine is a linear
regression on a randomly chosen nonlinear basis functions, which come from the
quantum state in a space of an exponentially large dimension, namely quantum
enhanced feature space. Furthermore, under some typical nonlinear function and
unitary operations settings to transform the observables, the output in Eq.
(96) can approximate any continuous function of the input. This property is
known as the universal approximation property (UAP), which implies that the
quantum extreme learning machine can handle a wide class of machine learning
tasks with at least the same power as the classical extreme learning machine
UAP_goto .

Here we should note that a similar approach, quantum kernel estimation, has
been taken in Ref. IBM_QML_ex . In quantum extreme learning machine, a
classical feature vector
$\phi_{i}(x)\equiv\langle\Phi(x)|Z_{i}|\Phi(x)\rangle$ is extracted from
observables on the quantum feature space $|\Phi(x)\rangle\equiv
V(x)|0\rangle^{\otimes n}$. Then linear regression is taken by using the
classical feature vector. On the other hand, in quantum kernel estimation,
quantum feature space is fully employed by using support vector machine with
the kernel functions
$K(x,x^{\prime})\equiv\langle\Phi(x)|\Phi(x^{\prime})\rangle$, which can be
estimated on a quantum computer. While classification power would be better
for quantum kernel estimation, it requires more quantum computational costs
both for learning and prediction in contrast to quantum extreme learning
machine.

In Fig. 2, we demonstrate quantum extreme learning machine for a two-class
classification task of a two-dimensional input $0\leq x_{0},x_{1}\leq 1$.
Class 0 and 1 are defined to be those being subject to
$(x_{0}-0.5)^{2}+(x_{1}-0.5)^{2}\leq 0.15$ and $>0.15$, respectively. The
linear readout weights $\\{w_{i}\\}$ are learned with 1000 randomly chosen
training data and prediction is performed with 1000 randomly chosen inputs.
The class 0 and 1 are determined whether or not the output $y$ is larger than
0.5. Quantum extreme learning machine with an 8-qubit quantum circuit shown in
Fig. 2 (a) succeeds to predict the class with $95\%$ accuracy. On the other
hand, a simple linear regression for $(x_{0},x_{1})$ results in $39\%$.
Moreover, quantum extreme learning machine with $U=I$, meaning no entangling
gate, also results in poor, $42\%$. In this way, the feature space enhanced by
quantum entangling operations is important to obtain a good performance in
quantum extreme learning machine.

###  IV.2 Quantum circuit learning

In the split of reservoir computing, dynamics of a physical system is not
fine-tuned but natural dynamics of the system is harnessed for machine
learning tasks. However, if we see the-state-of-the-art quantum computing
devices, the parameter of quantum operations can be finely tuned as done for
universal quantum computing. Therefore it is natural to extend quantum extreme
learning machine by tuning the parameters in the quantum circuit just like
feedfoward neural networks with back propagation.

Using parameterized quantum circuits for supervised machine leaning tasks such
as generalization of nonlinear functions and pattern recognitions have been
proposed in Refs. QCL ; FarhiQNN , which we call quantum circuit learning. Let
us consider the same situation with quantum extreme learning machine. The
state before the measurement is given by

$\displaystyle UV(x)|0\rangle^{\otimes n}.$ (98)

In the case of quantum extreme learning machine the unitary operation for a
nonlinear transformation with respect to the input parameter $x$ is randomly
chosen. However, the unitary operation $U$ may also be parameterized:

$\displaystyle U(\\{\phi_{k}\\})=\prod_{k}u(\phi_{k}).$ (99)

Thereby, the output from the quantum circuit with respect to an observable $A$

$\displaystyle\langle A(\\{\phi_{k}\\},x)\rangle=\langle 0|^{\otimes
n}V^{{\dagger}}(x)U(\\{\phi_{k}\\})^{{\dagger}}Z_{i}U(\\{\phi_{k}\\})V(x)|0\rangle^{\otimes
n}$

becomes a function of the circuit parameters $\\{\phi_{k}\\}$ in addition to
the input $x$. Then the parameters $\\{\phi_{k}\\}$ is tuned so as to minimize
the error between teacher data and the output, for example, by using the
gradient just like the output of the feedforward neural network.

Let us define a teacher dataset $\\{x^{(j)},y^{(j)}\\}$ and a quadratic loss
function

$\displaystyle L(\\{\phi_{k}\\})=\sum_{j}(\langle
A(\\{\phi_{k}\\},x^{(j)})\rangle-y^{(j)})^{2}.$ (100)

The gradient of the loss function can be obtained as follows:

$\displaystyle\frac{\partial}{\partial\phi_{l}}L(\\{\phi_{k}\\})$
$\displaystyle=$
$\displaystyle\frac{\partial}{\partial\phi_{l}}\sum_{j}(\langle
A(\\{\phi_{k}\\},x^{(j)})\rangle-y^{(j)})^{2}$ $\displaystyle=$
$\displaystyle\sum_{j}2(\langle
A(\\{\phi_{k}\\},x^{(j)})\rangle-y^{(j)})\frac{\partial}{\partial\phi_{l}}\langle
A(\\{\phi_{k}\\},x^{(j)})\rangle.$

Therefore if we can measure the gradient of the observable $\langle
A(\\{\phi_{k}\\},x^{(j)})\rangle$, the loss function can be minimized
according to the gradient descent.

If the unitary operation $u(\phi_{k})$ is given by

$\displaystyle u(\phi_{k})=W_{k}e^{-i(\phi_{k}/2)P_{k}},$ (101)

where $W_{k}$ is an arbitrary unitary, and $P_{k}$ is a Pauli operator. Then
the partial derivative with respect to the $l$th parameter can be analytically
calculated from the outputs $\langle A(\\{\phi_{k}\\},x^{(j)})\rangle$ with
shifting the $l$th parameter by $\pm\epsilon$ QCL ; Mitarai19 :

$\displaystyle\frac{\partial}{\partial\phi_{l}}\langle
A(\\{\phi_{k}\\},x^{(j)})\rangle$ $\displaystyle=$
$\displaystyle\frac{1}{2\sin\epsilon}(\langle
A(\\{\phi_{1},...,\phi_{l}+\epsilon,\phi_{l+1},...\\},x^{(j)})\rangle$
$\displaystyle-$ $\displaystyle\langle
A(\\{\phi_{1},...,\phi_{l}-\epsilon,\phi_{l+1},...\\},x^{(j)})\rangle).$

By considering the statistical error to measure the observable $\langle
A\rangle$, $\epsilon$ should be chosen to be $\epsilon=\pi/2$ so as to make
the denominator maximum. After measuring the partial derivatives for all
parameters $\phi_{k}$ and calculating the gradient of the loss function
$L(\\{\phi_{k}\\})$, the parameters are now updated by the gradient descent:

$\displaystyle\theta_{l}^{(m+1)}=\theta_{l}^{(m)}-\alpha\frac{\partial}{\partial\phi_{l}}L(\\{\phi_{k}\\}).$
(102)

The idea of using the parameterized quantum circuits for machine learning is
now widespread. After the proposal of quantum circuit learning based on the
analytical gradient estimation above QCL  and a similar idea FarhiQNN ,
several researches have been performed with various types of parameterized
quantum circuits CircuitCentric ; TensorNet ; Qneural_Chen ;
generalized_TensorNet ; QClassifier_Du  and various models and types of
machine learning including generative models GeneQML_Benedetti ; Gene_LiuWang
and generative adversarial models GAN_Situ ; GAN_Zeng ; GAN_Guzik . Moreover,
an expression power of the parameterized quantum circuits and its advantage
against classical probabilistic models have been investigated ExpressionPower
. Experimentally feasible ways to measure an analytical gradient of the
parameterized quantum circuits have been investigated Mitarai19 ;
Xanadu_gradient ; CalculusPQC . An advantage of using such a gradient for the
parameter optimization has been also argued in a simple setting HarrowGrad ,
while the parameter tuning becomes difficult because of the vanishing gradient
by an exponentially large Hilbert space Plateaus . Software libraries for
optimizing parameterized quantum circuits are now developing Xanadu ; VQNet .
Quantum machine learning on near-term devices, especially for quantum optical
systems, is proposed in Refs Opt_QML ; CV_QML . Quantum circuit learning with
parameterized quantum circuits has been already experimentally demonstrated on
superconducting qubit systems IBM_QML_ex ; Rigetti_QML_ex  and a trapped ion
system IonQ_QML_ex .

###  IV.3 Quantum reservoir computing

Figure 3: (a) Quantum reservoir computing. (b) Virtual nodes and temporal
multiplexing.

Now we return to the reservoir approach and extend quantum extreme learning
machine from non temporal tasks to temporal ones, namely, quantum reservoir
computing QRC . We consider a temporal task, which we explained in Sec. III.2.
The input is given by a time series $\\{x_{k}\\}_{k}^{L}$ and the purpose is
to learn a nonlinear temporal function:

$\displaystyle\bar{y}_{k}=f(\\{x_{j}\\}_{j}^{k}).$ (103)

To this end, the target time series $\\{\bar{y}_{k}\\}_{k=1}^{L}$ is also
provided as teacher.

Contrast to the previous setting with non temporal tasks, we have to fed input
into a quantum system sequentially. This requires us to perform an
initialization process during computation, and hence the quantum state of the
system becomes mixed state. Therefore, in the formulation of QRC, we will use
the vector representation of density operators, which was explained in Sec.
II.5.

In the vector representation of density operators, the quantum state of an
$N$-qubit system is given by a vector in a $4^{N}$-dimensional real vector
space, $\mbox{\boldmath$r$}\in\mathbb{R}^{4^{N}}$. In QRC, similarly to
recurrent neural networks, each element of the $4^{N}$-dimensional vector is
regarded as a hidden node of the network. As we seen in Sec. II.5, any
physical operation can be written as a linear transformation of the real
vector by a $4^{N}\times 4^{N}$ matrix $W$:

$\displaystyle\mbox{\boldmath$r$}^{\prime}=W\mbox{\boldmath$r$}.$ (104)

Now we see, from Eq. (104), a time evolution similar to the recurrent neural
network, $\mbox{\boldmath$r$}^{\prime}=\tanh(W\mbox{\boldmath$r$})$. However,
there is no nonlinearity such as $\tanh$ in each quantum operation $W$.
Instead, the time evolution $W$ can be changed according to the external input
$x_{k}$, namely $W_{x_{k}}$, which contrasts to the conventional recurrent
neural network where the input is fed additively $W\mbox{\boldmath$r$}+W^{\rm
in}x_{k}$. This allows the quantum reservoir to process the input information
$\\{x_{k}\\}$ nonlinearly, by repetitively feeding the input.

Suppose the input $\\{x_{k}\\}$ is normalized such that $0\leq x_{k}\leq 1$.
As an input, we replace a part of the qubits to the quantum state. The density
operator is given by

$\displaystyle\rho_{x_{k}}=\frac{I+(2x_{k}-1)Z}{2}.$ (105)

For simplicity, below we consider the case where only one qubit is replaced
for the input. Corresponding matrix $S_{x_{k}}$ is given by

$\displaystyle(S_{x_{k}})_{\mbox{\boldmath$j$}\mbox{\boldmath$i$}}={\rm
Tr}\left\\{P(\mbox{\boldmath$j$})\frac{I+(2x_{k}-1)Z}{2}\otimes{\rm Tr}_{\rm
replace}[P(\mbox{\boldmath$i$})]\right\\}/2^{N},$

where ${\rm Tr}_{\rm replace}$ indicates a partial trace with respect to the
replaced qubit. With this definition, we have

$\displaystyle\rho^{\prime}={\rm Tr}_{\rm
replace}[\rho]\otimes\rho_{x_{k}}\Leftrightarrow\mbox{\boldmath$r$}^{\prime}=S_{x_{k}}\mbox{\boldmath$r$}.$
(106)

The unitary time evolution, which is necessary to obtain a nonlinear behavior
with respect to the input valuable $x_{k}$, is taken as a Hamiltonian dynamics
$e^{-iH\tau}$ for a given time interval $\tau$. Let us denote its
representation on the vector space by $U_{\tau}$:

$\displaystyle\rho^{\prime}=e^{-iH\tau}\rho
e^{iH\tau}\Leftrightarrow\mbox{\boldmath$r$}^{\prime}=U_{\tau}\mbox{\boldmath$r$}.$
(107)

Then, a unit time step is written as an input-depending linear transformation:

$\displaystyle\mbox{\boldmath$r$}((k+1)\tau)=U_{\tau}S_{x_{k}}\mbox{\boldmath$r$}(k\tau).$
(108)

where $\mbox{\boldmath$r$}(k\tau)$ indicates the hidden nodes at time $k\tau$.

Since the number of the hidden nodes are exponentially large, it is not
feasible to observe all nodes from experiments. Instead, a set of observed
nodes $\\{\bar{r}_{l}\\}_{l=1}^{M}$, which we call true nodes, is defined by a
$M\times 4^{N}$ matrix $R$,

$\displaystyle\bar{r}_{l}(k\tau)=\sum_{\mbox{\boldmath$i$}}R_{l\mbox{\boldmath$i$}}r_{\mbox{\boldmath$i$}}(k\tau).$
(109)

The number of true nodes $M$ has to be a polynomial in the number of qubits
$N$. That is, from exponentially many hidden nodes, a polynomial number of
true nodes are obtained to define the output from QR (see Fig. 3 (a)):

$\displaystyle y_{k}=\sum_{l}W^{\rm out}_{l}\bar{r}_{l}(k\tau),$ (110)

where $W_{\rm out}$ is the readout weights, which is obtained by using the
training data. For simplicity, we take the single-qubit Pauli $Z$ operator on
each qubit as the true nodes, i.e.,

$\displaystyle\bar{r}_{l}={\rm Tr}[Z_{l}\rho],$ (111)

so that if there is no dynamics these nodes simply provide a linear output
$(2x_{k}-1)$ with respect to the input $x_{k}$.

Moreover, in order to improve the performance we also perform the temporal
multiplexing. The temporal multiplexing has been found to be useful to extract
complex dynamics on the exponentially large hidden nodes through the
restricted number of the true nodes QRC . In temporal multiplexing, not only
the true nodes just after the time evolution $U_{\tau}$, also at each of the
subdivided $V$ time intervals during the unitary evolution $U_{\tau}$ to
construct $V$ virtual nodes, as shown in Fig. 3 (b). After each input by
$S_{x_{k}}$, the signals from the hidden nodes (via the true nodes) are
measured for each subdevided intervals after the time evolution by
$U_{v\tau/V}$ ($v=1,2,...V$), i.e.,

$\displaystyle\mbox{\boldmath$r$}(k\tau+(v/V)\tau)\equiv
U_{(v/V)\tau}S_{x_{k}}\mbox{\boldmath$r$}(k\tau).$ (112)

In total, now we have $N\times V$ nodes, and the output is defined as their
linear combination:

$\displaystyle y_{k}=\sum_{l=1}^{N}\sum_{v=1}^{V}W^{\rm
out}_{j,v}\bar{r}_{l}(k\tau+(v/V)\tau).$ (113)

By using the teacher data $\\{\bar{y}_{k}\\}_{k}^{L}$, the linear readout
weights $W^{\rm out}_{j,v}$ can be determined by using the pseudo inverse. In
Ref. QRC , the performance of QRC has been investigated extensively for both
binary and continuous inputs. The result shows that even if the number of the
qubits are small like 5-7 qubits the performance as powerful as the echo state
network of the 100-500 nodes have been reported both in short term memory and
parity check capacities. Note that, although we do not go into detail in this
chapter, the technique called spatial multiplexing QRC_multiplex , which
exploits multiple quantum reservoirs with common input sequence injected, is
also introduced to harness quantum dynamics as a computational resource.
Recently, QRC has been further investigated in Refs. QRC_OTOC ; QRC_sig ;
QRC_aus . Specifically, in Ref. QRC_sig , the authors use quantum reserovir
computing to detect many-body entanglement by estimating nonlinear functions
of deinsity operators like entropy.

###  IV.4 Emulating chaotic attractors using quantum dynamics

Figure 4: Demonstrations of chaotic attractor emulations. (a) Lorenz
attractor. (b) Mackey-Glass system. (c) Rössler attractor. (d) Hénon map. The
dotted line shows the time step when the system is switched from teacher
forced state to autonomous state. In the right side, delayed phase diagrams of
learned dynamics are shown.

To see a performance of QRC, here we demonstrate an emulation of chaotic
attractors. Suppose $\\{x_{k}\\}_{k}^{L}$ is a discretized time sequence being
subject to a complex nonlinear equation, which might has a chaotic behavior.
In this task, the target, which the network is to output, is defined to be

$\displaystyle\bar{y}_{k}=x_{k+1}=f(\\{x_{j}\\}_{j=1}^{k}).$ (114)

That is, the system learns the input of the next step. Once the system
successfully learns $\bar{y}_{k}$, by feeding the output into the input of the
next step of the system, the system evolves autonomously.

Here we employ the following target time series from chaotic attractors: (i)
Lorenz attractor,

$\displaystyle\frac{dx}{dt}$ $\displaystyle=$ $\displaystyle a(y-x),$ (115)
$\displaystyle\frac{dy}{dt}$ $\displaystyle=$ $\displaystyle x(b-z)-y,$ (116)
$\displaystyle\frac{dz}{dt}$ $\displaystyle=$ $\displaystyle xy-cz,$ (117)

with $(a,b,c)=(10,28,8/3)$, (ii) the chaotic attractor of Mackey-Glass
equation,

$\displaystyle\frac{d}{dt}x(t)=\beta\frac{x(t-\tau)}{1+x(t-\tau)^{n}}-\gamma
x(t)$ (118)

with $(\beta,\gamma,n)=(0.2,0.1,10)$ and $\tau=17$, (iii) Rössler attoractor,

$\displaystyle\frac{dx}{dt}$ $\displaystyle=$ $\displaystyle-y-z,$ (119)
$\displaystyle\frac{dy}{dt}$ $\displaystyle=$ $\displaystyle x+ay,$ (120)
$\displaystyle\frac{dz}{dt}$ $\displaystyle=$ $\displaystyle b+z(x-c),$ (121)

with $(0.2,0.2,5.7)$, and (iv) Hénon map,

$\displaystyle x_{t+1}=1-1.4x_{t}+0.3x_{t-1}.$ (122)

Regarding (i)-(iii), the time series is obtained by using the fourth-order
Runge-Kutta method with step size 0.02, and only $x(t)$ is employed as a
target. For the time evolution of quantum reservoir, we employ a fully
connected transverse-field Ising model

$\displaystyle H=\sum_{ij}J_{ij}X_{i}X_{j}+hZ_{i},$ (123)

where the coupling strengths are randomly chosen such that $J_{ij}$ is
distributed randomly from $[-0.5,0.5]$ and $h=1.0$. The time interval and the
number of the virtual nodes are chosen to be $\tau=4.0$ and $v=10$ so as to
obtain the best performance. The first $10^{4}$ steps are used for training.
After the linear readout weights are determined, several $10^{3}$ steps are
predicted by autonomously evolving the quantum reservoir. The results are
shown in Fig. 4 for each of (a) Lorenz attractor, (b) the chaotic attractor of
Mackey-Glass system, (c) Rössler attractor, and (d) Hénon map. All these
results show that training is done well and the prediction is successful for
several hundreds steps. Moreover, the output from the quantum reservoir also
successfully reconstruct the structures of these chaotic attractors as you can
see from the delayed phase diagram.

##  V Conclusion and Discussion

Here we reviewed quantum reservoir computing and related approaches, quantum
extreme learning machine and quantum circuit learning. The idea of quantum
reservoir computing comes from the spirit of reservoir computing, i.e.,
outsourcing information processing to natural physical systems. This idea is
best suited to quantum machine learning on near-term quantum devices in NISQ
(noisy intermediate quantum) era. Since reservoir computing uses complex
physical systems as a feature space to construct a model by the simple linear
regression, this approach would be a good way to understand the power of a
quantum enhanced feature space.

## Acknowledgement

KF is supported by KAKENHI No.16H02211, JST PRESTO JPMJPR1668, JST ERATO
JPMJER1601, and JST CREST JPMJCR1673. KN is supported by JST PRESTO Grant
Number JPMJPR15E7, Japan, by JSPS KAKENHI Grant Numbers JP18H05472,
JP16KT0019, and JP15K16076. KN would like to acknowledge Dr. Quoc Hoan Tran
for his helpful comments. This work is supported by MEXT Quantum Leap Flagship
Program (MEXT Q-LEAP) Grant No. JPMXS0118067394.

## References

  * (1) R.P. Feynman, Simulating physics with computers, Int. J. Theor. Phys. 21, 467 (1982).
  * (2) M.A. Nielsen and I. L. Chuang, Quantum computation and quantum information, (Cambridge university press 2010).
  * (3) K. Fujii, Quantum Computation with Topological Codes -From Qubit to Topological Fault-Tolerance-, SpringerBriefs in Mathematical Physics (Springer-Verlag 2015).
  * (4) P. W. Shor, Algorithms for quantum computation: Discrete logarithms and factoring, In Proceedings of the 35th Annual Symposium on Foundations of Computer Science, 124 (1994).
  * (5) R. Barends et al., Superconducting quantum circuits at the surface code threshold for fault tolerance, Nature 508, 500 (2014).
  * (6) J. Kelly et al., State preservation by repetitive error detection in a superconducting quantum circuit, Nature 519, 66 (2015).
  * (7) J. Preskill, Quantum Computing in the NISQ era and beyond., Quantum 2, 79 (2018).
  * (8) S. Boixo et al., Characterizing quantum supremacy in near-term devices., Nature Physics 14, 595 (2018).
  * (9) J.I. Cirac and P. Zoller, Goals and opportunities in quantum simulation, Nat. Phys. 8, 264 (2012).
  * (10) I. Bloch, J. Dalibard, and S. Nascimbène, Quantum simulations with ultracold quantum gases, Nat. Phys. 8, 267 (2012).
  * (11) I. M. Georgescu, S. Ashhab, and F. Nori, Quantum simulation, Rev. of Mod. Phys. 86, 153 (2014).
  * (12) T. Kadowaki and H. Nishimori, Quantum annealing in the transverse Ising model, Phys. Rev. E 58, 5355 (1998).
  * (13) E. Farhi, J. Goldstone, S. Gutmann, J. Lapan, A. Lundgren, and D. Preda, A quantum adiabatic evolution algorithm applied to random instances of an NP-complete problem, Science 292, 472 (2001).
  * (14) T. F. Rønnow, Z. Wang, J. Job, S. Boixo, S. V. Isakov, D. Wecker, J. M. Martinis, D. A. Lidar, M. Troyer Defining and detecting quantum speedup, Science 345, 420 (2014).
  * (15) S. Boixo, T. F. Rønnow, S. V. Isakov, Z. Wang, D. Wecker, D. A. Lidar, J. M. Martinis, and M. Troyer, Evidence for quantum annealing with more than one hundred qubits, Nat. Phys. 10, 218 (2014).
  * (16) T. Morimae, K. Fujii, and J. F. Fitzsimons, Hardness of classically simulating the one-clean-qubit model, Phys. Rev. Lett. 112, 130502 (2014).
  * (17) K. Fujii, H. Kobayashi, T. Morimae, H. Nishimura, S. Tamate, and S. Tani, Power of Quantum Computation with Few Clean Qubits, Proceedings of 43rd International Colloquium on Automata, Languages, and Programming (ICALP 2016), pp.13:1-13:14.
  * (18) K. Fujii and S. Tamate, Computational quantum-classical boundary of complex and noisy quantum systems, Sci. Rep. 6, 25598 (2016).
  * (19) M. Rabinovich, R. Huerta, and G. Laurent, Transient dynamics for neural processing Science 321, 48 (2008).
  * (20) J. Dambre, D. Verstraeten, B. Schrauwen, and S. Massar, Information processing capacity of dynamical systems, Sci. Rep. 2, 514 (2012).
  * (21) H. Jaeger and H. Haas, Harnessing nonlinearity: predicting chaotic systems and saving energy in wireless communication, Science 304, 78 (2004).
  * (22) 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).
  * (23) C. Fernando and S. Sojakka, Pattern recognition in a bucket In Lecture Notes in Computer Science 2801, p. 588 (Springer, 2003).
  * (24) L. Appeltant, M. C. Soriano, G. Van der Sande, J. Danckaert, S. Massar, J. Dambre, B. Schrauwen, C. R. Mirasso, and I. Fischer, Information processing using a single dynamical node as complex system. Nat. Commun. 2, 468 (2011).
  * (25) D. Woods and T. J. Naughton, Photonic neural networks., Nat. Phys. 8, 257 (2012).
  * (26) L. Larger, M. C. Soriano, D. Brunner, L. Appeltant, J. M. Gutierrez, L. Pesquera, C. R. Mirasso, and I. Fischer, Photonic information processing beyond Turing: an optoelectronic implementation of reservoir computing, Optics Express 20, 3241 (2012).
  * (27) Y. Paquot, F. Duport, A. Smerieri, J. Dambre, B. Schrauwen, M. Haelterman, and S. Massar, Optoelectronic Reservoir Computing, Sci. Rep. 2, 287 (2012).
  * (28) D. Brunner, M. C. Soriano, C. R. Mirasso, and I. Fischer, Parallel photonic information processing at gigabyte per second data rates using transient states, Nat. Commun. 4, 1364 (2013).
  * (29) K. Vandoorne, P. Mechet, T. V. Vaerenbergh, M. Fiers, G. Morthier, D. Verstraeten, B. Schrauwen, J. Dambre, and P. Bienstman, Experimental demonstration of reservoir computing on a silicon photonics chip Nat. Commun. 5, 3541 (2014).
  * (30) A. Z. Stieg, A. V. Avizienis, H. O. Sillin, C. Martin-Olmos, M. Aono, and J. K. Gimzewski Emergent criticality in complex turing B-type atomic switch networks, Adv. Mater. 24, 286 (2012).
  * (31) H. Hauser, A. J. Ijspeert, R. M. Füchslin, R. Pfeifer, and W. Maass, Towards a theoretical foundation for morphological computation with compliant bodies Biol. Cybern. 105, 355 (2011).
  * (32) K. Nakajima, H. Hauser, R. Kang, E. Guglielmino, D. G. Caldwell, and R. Pfeifer, Computing with a Muscular-Hydrostat System, Proceedings of 2013 IEEE International Conference on Robotics and Automation (ICRA), 1496 (2013).
  * (33) K. Nakajima, H. Hauser, R. Kang, E. Guglielmino, D. G. Caldwell, and R. Pfeifer, A soft body as a reservoir: case studies in a dynamic model of octopus-inspired soft robotic arm Front. Comput. Neurosci. 7, 1 (2013).
  * (34) K. Nakajima, T. Li, H. Hauser, and R. Pfeifer, Exploiting short-term memory in soft body dynamics as a computational resource, J. R. Soc. Interface 11, 20140437 (2014).
  * (35) K. Nakajima, H. Hauser, T. Li, and R. Pfeifer, Information processing via physical soft body, Sci. Rep. 5, 10487 (2015).
  * (36) K. Caluwaerts, J. Despraz, A. Işçen, A. P. Sabelhaus, J. Bruce, B. Schrauwen, and V. SunSpiral, Design and control of compliant tensegrity robots through simulations and hardware validation, J. R. Soc. Interface 11, 20140520 (2014).
  * (37) Y. LeCun, Y. Bengio, and G. Hinton, Deep Learning Nature 521, 436 (2015).
  * (38) K. Fujii and K. Nakajima, Harnessing Disordered-Ensemble Quantum Dynamics for Machine Learning, Phys. Rev. Applied 8, 024030 (2017).
  * (39) M. Negoro et al., Machine learning with controllable quantum dynamics of a nuclear spin ensemble in a solid, arXiv:1806.10910 (2018).
  * (40) K Nakajima et al., Boosting computational power through spatial multiplexing in quantum reservoir computing, Phys. Rev. Applied 11, 034021 (2019).
  * (41) A Kutvonen, K Fujii, and T Sagawa, Optimizing a quantum reservoir computer for time series prediction, Sci Rep 10, 14687 (2020).
  * (42) Q. H. Tran and K. Nakajima, Higher-order quantum reservoir computing, arXiv:2006.08999 (2020).
  * (43) K. Mitarai et al., Quantum circuit learning, Phys. Rev. A 98, 032309 (2018).
  * (44) G.-B. Huang, Q.-Y. Zhu, and C.-K. Siew, Extreme learning machine: theory and applications, Neurocomputing 70, 489 (2006).
  * (45) D. Verstraeten, B. Schrauwen, M. D’Haene, and D. Stroobandt, An experimental unification of reservoir computing methods, Neural Netw. 20, 391 (2007).
  * (46) V. Havlicek et al, Supervised learning with quantum enhanced feature spaces, Nature 567, 209 (2019).
  * (47) E. Farhi and H. Neven, Classification with quantum neural networks on near term processors, arXiv:1802.06002 (2018).
  * (48) M. Schuld et al., Circuit-centric quantum classifiers, Phys. Rev. A 101, 032308 (2020).
  * (49) H. Chen et al., Universal discriminative quantum neural networks, arXiv:1805.08654 (2018).
  * (50) W. Huggins et al., Towards Quantum Machine Learning with Tensor Networks, Quantum Science and Technology 4, 024001 (2019).
  * (51) I. Glasser, N. Pancotti, and J. I. Cirac, From probabilistic graphical models to generalized tensor networks for supervised learning, arXiv:1806.05964 (2018).
  * (52) Y. Du et al., Implementable Quantum Classifier for Nonlinear Data, arXiv:1809.06056 (2018).
  * (53) M. Benedetti et al., Adversarial quantum circuit learning for pure state approximation, New J. Phys. 21, 043023 (2019).
  * (54) M. Benedetti et al., A generative modeling approach for benchmarking and training shallow quantum circuits, npj Quantum Information 5, 45 (2019).
  * (55) J.-G. Liu and L. Wang, Phys. Rev. A 98, 062324 (2018).
  * (56) H. Situ et al., Quantum generative adversarial network for generating discrete data, Information Sciences 538, 193 (2020).
  * (57) J. Zeng et al., Learning and Inference on Generative Adversarial Quantum Circuits, Phys. Rev. A 99, 052306 (2019).
  * (58) J. Romero and A. Aspuru-Guzik Variational quantum generators: Generative adversarial quantum machine learning for continuous distributions, arXiv:1901.00848 (2019).
  * (59) Y. Du et al., The expressive power of parameterized quantum circuits, Phys. Rev. Research 2, 033125 (2020).
  * (60) M. Schuld et al., Evaluating analytic gradients on quantum hardware, Phys. Rev. A 99, 032331 (2019).
  * (61) K. Mitarai and K. Fujii, Methodology for replacing indirect measurements with direct measurements, Phys. Rev. Research 1, 013006 (2019).
  * (62) J. G. Vidal, and D. O. Theis, Calculus on parameterized quantum circuits, arXiv:1812.06323 (2018).
  * (63) A. Harrow and N. John Low-depth gradient measurements can improve convergence in variational hybrid quantum-classical algorithms, arXiv:1901.05374 (2019).
  * (64) J. R. McClean et al. Barren plateaus in quantum neural network training landscapes, Nature Communications 9, 4812 (2018).
  * (65) V. Bergholm et al., PennyLane: Automatic differentiation of hybrid quantum-classical computations, arXiv:1811.04968 (2018).
  * (66) Z.-Y. Chen et al., VQNet: Library for a Quantum-Classical Hybrid Neural Network, arXiv:1901.09133 (2019).
  * (67) G. R. Steinbrecher et al., Quantum optical neural networks, npj Quantum Information 5, 60 (2019).
  * (68) N. Killoran et al., Continuous-variable quantum neural networks, Phys. Rev. Research 1, 033063 (2019).
  * (69) C. M. Wilson et al., Quantum Kitchen Sinks: An algorithm for machine learning on near-term quantum computers, arXiv:1806.08321 (2018).
  * (70) D. Zhu et al., Training of Quantum Circuits on a Hybrid Quantum Computer, Science Advances 5, 9918 (2019).
  * (71) S. Ghosh et al., Quantum reservoir processing, npj Quantum Information 5, 35 (2019).
  * (72) S. Ghosh, T. Paterek, and T. C. H. Liew, Quantum neuromorphic platform for quantum state preparation, Phys. Rev. Lett. 123, 260404 (2019).
  * (73) S. Ghosh et al., Reconstructing quantum states with quantum reservoir networks, IEEE Trans. Neural Netw. Learn. Syst., pp. 1–8 (2020).
  * (74) J. Chen and H. I. Nurdin, Learning Nonlinear Input-Output Maps with Dissipative Quantum Systems, Quantum Information Processing 18, 198 (2019).
  * (75) T. Goto, Q. H. Tran and K. Nakajima, Universal approximation property of quantum feature maps, arXiv: 2009.00298 (2020).