\section{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~\cite{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~\cite{ShorFact}.
Therefore, considerable experimental effort has been devoted to realising full-fledged universal quantum computers~\cite{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~\cite{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~\cite{QuantumSimulator1,QuantumSimulator2,QuantumSimulator3}, and a restricted class of quantum dynamics, known as adiabatic dynamics, has also been applied to combinatorial optimisation problems~\cite{Nishimori,Farhi,Dwave,Dwave2}.
However, complex real-time quantum dynamics, which is one of the most difficult tasks for classical computers to simulate
~\cite{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 {\it reservoir}, generating transient dynamics that facilitates the separation of input states \cite{Transient}.
If the dynamics of the reservoir involve both adequate memory and nonlinearity \cite{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) \cite{Jaeger0} or models of neurons for liquid state machines \cite{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 \cite{Bucket}, analogue circuits and optoelectronic systems \cite{Laser0,Laser1,Laser1b,Laser2,Laser3,Laser4}, and neuromorphic chips \cite{Neuromorphic0}.
Recently, it has been reported that the mechanical bodies of soft and compliant robots have also been successfully used as a reservoir \cite{Helmut0,Kohei0,Kohei1,Kohei2,Kohei3,Ken1}.
In contrast to the refinements required by learning algorithms, such as in deep learning \cite{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~\cite{QRC}.
Since the idea of QRC has been proposed in Ref.~\cite{QRC},
its proof-of-principle experimental demonstration for non temporal tasks~\cite{QRC_ex} and performance analysis and improvement~\cite{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 \cite{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~\ref{sec:quantum_mechanics},
we will provide a pedagogical introduction to quantum mechanics
for those who are not familiar to it and fix our notation.
In Sec~\ref{sec:machine_learning}, we will briefly mention to
several machine learning techniques like,
linear and nonlinear regressions, temporal machine learning tasks
and reservoir computing.
In Sec~\ref{sec:quantum_machine_learning},
we will explain QRC and related approaches,
quantum extreme learning machine~\cite{QRC_ex} and quantum circuit learning~\cite{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.~\ref{sec:conclusion} is devoted to conclusion and discussion.

\section{Pedagogical introduction to quantum mechanics}
\label{sec: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.~\ref{sec:machine_learning}.

\subsection{Quantum state}
A state of a quantum system is described by a {\it state vector},
\begin{eqnarray}
|\psi \rangle =
\left(
\begin{array}{c}
c_1
\\
\vdots
\\
c_d
\end{array}
\right)
\end{eqnarray}
on a complex $d$-dimensional system $\mathbb{C}^d$,
where the symbol $| \cdot \rangle$ is called {\it ket}
and indicates a complex column vector.
Similarly, $\langle \cdot |$ is called {\it bra}
and indicates a complex row vector, and they are
related complex conjugate,
\begin{eqnarray}
\langle \psi | = |\psi \rangle ^{\dag}  =
\left(
\begin{array}{ccc}
c^*_1 & \cdots & c^*_d
\end{array}
\right).
\end{eqnarray}
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
\begin{eqnarray}
|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),
\end{eqnarray}
a quantum state in the $d$-dimensional system
can be described simply by
\begin{eqnarray}
|\psi \rangle = \sum _{i=1}^{d} c_i | i\rangle .
\end{eqnarray}
The state is said to be a {\it superposition state} of $|i\rangle$.
The coefficients $\{c_i\}$ are complex,
and called {\it complex probability amplitudes}.
If we measure the system in the basis $\{ |i \rangle \}$,
we obtain the measurement outcome $i$ with
a probability
\begin{eqnarray}
p_i = | \langle i | \psi \rangle |^2 =  | c_i | ^2 ,
\end{eqnarray}
and hence the complex probability amplitudes have
to be normalized as follows
\begin{eqnarray}
|\langle \psi | \psi \rangle |^2 = \sum _{i=1}^d |c_i|^2 =1.
\end{eqnarray}
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
\begin{eqnarray}
\sum _i a_i p_i = \langle \psi | A |\psi \rangle \equiv \langle A \rangle,
\end{eqnarray}
where we define an hermitian operator
\begin{eqnarray}
A = \sum _i a_i |i\rangle \langle i|,
\end{eqnarray}
which is called {\it 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.

\subsection{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{\"o}dinger equation, is
given by
\begin{eqnarray}
i \frac{\partial}{\partial t} |\psi (t)\rangle = H |\psi (t)\rangle .
\end{eqnarray}
This equation can be formally solved by
\begin{eqnarray}
|\psi (t) \rangle = e^{-i H t } |\psi (0)\rangle.
\end{eqnarray}
Therefore the time evolution is given by
an operator $e^{-i H t}$,
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$
\begin{eqnarray}
|\psi_{\rm out} \rangle  = U |\psi _{\rm in} \rangle .
\end{eqnarray}
In quantum computing, the time evolution $U$ is sometimes called
{\it quantum gate}.

\subsection{Qubits}
The smallest nontrivial quantum system is a two-dimensional quantum system $\mathbb{C}^2$, which is
called {\it quantum bit} or {\it qubit}:
\begin{eqnarray}
\alpha |0\rangle + \beta |1\rangle , \;\;\; (|\alpha|^2 + |\beta |^2 = 1).
\end{eqnarray}
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\}$,
\begin{eqnarray}
|x_1 \rangle \otimes | x_2 \rangle \otimes \cdots \otimes |x_n \rangle,
\end{eqnarray}
which is simply denoted by
\begin{eqnarray}
|x_1  x_2  \cdots x_n \rangle.
\end{eqnarray}
Then a state of the $n$-qubit system can be described as
\begin{eqnarray}
|\psi \rangle = \sum _{x_1 ,x_2 ,...,x_n} \alpha _{x_1 ,x_2, ...,x_n} |x_1  x_2  \cdots x_n \rangle.
\end{eqnarray}
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.

\subsection{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$:
\begin{eqnarray}
\rho = |\psi \rangle \langle \psi |.
\end{eqnarray}
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
\begin{eqnarray}
p_i = {\rm Tr}[ |i \rangle \langle i | \rho ],
\end{eqnarray}
where ${\rm Tr}$ is the matrix trace.
An expectation value of an observable $A$ is given by
\begin{eqnarray}
\langle A \rangle = {\rm Tr}[A\rho].
\end{eqnarray}
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
\begin{eqnarray}
p_i = \sum _k q_k {\rm Tr}[|i \rangle \langle i | \rho _k],
\end{eqnarray}
where $\rho _k = | \psi _k \rangle \langle \psi _k |$.
By using linearity of the trace function, this reads
\begin{eqnarray}
p_i = {\rm Tr}[|i \rangle \langle i | \sum _k q_k \rho _k].
\end{eqnarray}
Now we interpret that the density
operator is given by
\begin{eqnarray}
\rho = \sum _k q_k |\psi _k \rangle \langle \psi _k | .
\end{eqnarray}
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:
\begin{eqnarray}
\rho = \sum \lambda _i | \lambda _i \rangle \langle \lambda _i | ,
\end{eqnarray}
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
\begin{eqnarray}
\rho (t) = e^{-i H t} \rho (0) e^{iHt}
\end{eqnarray}
or
\begin{eqnarray}
\rho _{\rm out} = U \rho _{\rm in} U^{\dag}.
\end{eqnarray}
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
\begin{eqnarray}
\rho _{\rm out} = p U \rho U ^{\dag} + (1-p) V \rho V ^{\dag}.
\end{eqnarray}
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
\begin{eqnarray}
\sum _i {\rm Tr}[|i \rangle \langle i | \rho ] |i \rangle \langle i |
= \sum _i |i \rangle \langle i | \rho |i \rangle \langle i |.
\end{eqnarray}
Therefore if we define a map from a density operator to another,
which we call {\it superoperator},
\begin{eqnarray}
\mathcal{M} (\cdots) = \sum _i |i \rangle \langle i | (\cdots) |i \rangle \langle i |,
\end{eqnarray}
the above non-selective measurement (forgetting about the measurement outcomes)
is simply written by
\begin{eqnarray}
\mathcal{M} (\rho ).
\end{eqnarray}
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 ^{\dag} K_i = I$ with an identity operator $I$:
\begin{eqnarray}
\mathcal{K}(\rho) = \sum _i K_i \rho K^{\dag}_i.
\end{eqnarray}
The operators $\{K_i\}$ are called {\it Kraus operators}.

\subsection{Vector representation of density operators}
\label{subsec:vector_density}
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 {\it Hilbert-Schmidt inner product}, by
\begin{eqnarray}
{\rm Tr} [A^{\dag} B ] .
\end{eqnarray}
The operators on the $n$-qubit system can be
spanned by the tensor product of Pauli operators $\{ I, X, Y,Z\}^{\otimes n}$,
\begin{eqnarray}
P(\mb{i}) =\bigotimes_{k=1}^{n} \sigma_{i_{2k-1}i_{2k}}.
\end{eqnarray}
where
$\sigma_{ij}$ is the Pauli operators:
\begin{eqnarray*}
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),
\end{eqnarray*}
\begin{eqnarray}
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).
\end{eqnarray}
Since the Pauli operators
constitute a complete basis on the operator space,
any operator $A$ can be decomposed into a linear combination of $P(\mb{i})$,
\begin{eqnarray}
A = \sum _ {\mb{i}} a_{\mb{i}} P(\mb{i}).
\end{eqnarray}
The coefficient $a_{\mb{i}}$ can be calculated
by using the Hilbert-Schmidt inner product as follows:
\begin{eqnarray}
a_{\mb{i}} = {\rm Tr}[P(\mb{i})A]/2^n,
\end{eqnarray}
by virtue of the orthogonality
\begin{eqnarray}
{\rm Tr} [P(\mb{i})P(\mb{j})] /2^n = \delta _{\mb{i},\mb{j}}.
\end{eqnarray}
The number of the $n$-qubit Pauli operators $\{ P (\mb{i})\}$
is $4^n$, and hence
a density operator $\rho$ of the $n$-qubit system can be represented
as a $4^n$-dimensional vector
\begin{eqnarray}
\mb{r} =
\left(
\begin{array}{c}
r_{00...0}
\\
\vdots
\\
r_{11...1}
\end{array}
\right),
\end{eqnarray}
where $r_{00...0}=1/2^n$ because of ${\rm Tr}[\rho]=1$.
Moreover, because $P(\mb{i})$ is hermitian,
$\mb{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 $\mb{r}$:
\begin{eqnarray}
\rho' = \mathcal{K}(\rho) \Leftrightarrow \mb{r}' = K \mb{r},
\end{eqnarray}
where the matrix element is given by
\begin{eqnarray}
K_{\mb{i}\mb{j}} = {\rm Tr}[ P(\mb{i}) \mathcal{K} \left(P(\mb{j}) \right)]/2^n .
\end{eqnarray}
In this way,
a density operator $\rho$ and a quantum operation $\mathcal{K}$
on it can be represented by a vector $\mb{r}$ and a matrix $K$, respectively.

\section{Machine learning and reservoir approach}
In this section, we briefly introduce machine learning
and reservoir approaches.

\label{sec:machine_learning}
\subsection{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 {\it linear regression},
which models $f(x)$ as a linear function with respect to the input:
\begin{eqnarray}
f(x) = \sum_{i=1}^{d} w_i x_i + w_0.
\end{eqnarray}
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. {\it 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:
\begin{eqnarray}
 \min _{\{w_i\}} \sum _j (\sum_{i=0}^{d} w_i x^{(j)}_i  - y^{(j)})^2,
\end{eqnarray}
where we introduced a constant node $x_0 =1$.
This corresponds to solving a superimposing equations:
\begin{eqnarray}
\mathbf{y}  = \mathbf{X} \mathbf{w},
\end{eqnarray}
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} = U D V^{T}$ to be
\begin{eqnarray}
\mathbf{X}^{+} = V D U^{T}.
\end{eqnarray}

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}$:
\begin{eqnarray}
 W^{\rm in} x.
\end{eqnarray}
Then each element of the hidden nodes is now processed by
a nonlinear activation function $\sigma$ such as $\tanh$,
which is denoted by
\begin{eqnarray}
\sigma( W^{\rm in} x).
\end{eqnarray}
Finally the output is extracted by an output weight $W^{\rm out}$ ($1 \times N$ dimensional matrix):
\begin{eqnarray}
 W^{\rm out} \sigma( W^{\rm in} x).
\end{eqnarray}
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:
\begin{eqnarray}
 W^{\rm out} \sigma \left( W^{(l)} \cdots \sigma \left( W^{(1)} \sigma( W^{\rm in} x) \right) \right),
\end{eqnarray}
which is called multi-layer perceptron or deep neural network.

\subsection{Temporal task}
\label{subsec:temporal}
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,
\begin{eqnarray}
\bar y_k = f(\{ x_j\}_{j=1}^{k}).
\end{eqnarray}
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
\begin{eqnarray}
\mb{r}=\left(
\begin{array}{c}
r_1
\\
\vdots
\\
r_m
\end{array}
\right).
\end{eqnarray}
To process the input time series,
the nodes evolve by
\begin{eqnarray}
\mb{r}(k+1) = \sigma [W \mb{r}(k) +W^{\rm in} x_k],
\end{eqnarray}
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
\begin{eqnarray}
y_k = W^{\rm out} \mb{r}(k).
\end{eqnarray}
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.

\subsection{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., {\it reservoir}, and
by performing simple linear regression on it.

Let us first see a reservoir approach on a feedforward neural network,
which is called {\it extreme learning machine}~\cite{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:
\begin{eqnarray}
\{ \phi _1 (x), \phi_2 (x),..., \phi_N(x)\}.
\end{eqnarray}
Now the output is defined as a linear combination of these
\begin{eqnarray}
\sum _i w_i \phi _i (x) + w_0
\end{eqnarray}
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 {\it echo state network} is similar but employs the reservoir idea
for the recurrent neural network~\cite{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
\begin{eqnarray}
\sum _{k} (y_k - \bar y_k)^2.
\end{eqnarray}
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 {\it reservoir} for information processing,
namely, physical reservoir computing~\cite{Bucket,Laser0,Laser1,Laser1b,Laser2,Laser3,Laser4,Neuromorphic0,Helmut0,Kohei0,Kohei1,Kohei2,Kohei3,Ken1}.

\section{Quantum machine learning on near-term quantum devices}
\label{sec:quantum_machine_learning}
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~\cite{QRC_ex}.
While it is an opposite direction to reservoir computing,
we will also see  {\it quantum circuit learning} (QCL)~\cite{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.

\subsection{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
\begin{eqnarray}
|0\rangle ^{\otimes n}.
\end{eqnarray}
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:
\begin{eqnarray}
V(x)|0\rangle ^{\otimes n}.
\end{eqnarray}
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})$:
\begin{eqnarray}
e^{-i \theta Y} |0\rangle = \sqrt{x} |0\rangle + \sqrt{1-x} |1\rangle.
\end{eqnarray}
The expectation value of $Z$ with respect to $e^{-i \theta Y} |0\rangle$
becomes
\begin{eqnarray}
\langle  Z \rangle = 2x -1,
\end{eqnarray}
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$:
\begin{eqnarray}
\theta = \arccos(\sqrt{\phi(x)}).
\end{eqnarray}
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
\begin{eqnarray*}
V(x) |0\rangle ^{\otimes n} &=& (\sqrt{x} |0\rangle + \sqrt{1-x} |1\rangle )^{\otimes n}
\nonumber \\
&=& (1-x)^{n/2} \sum _{i_1,...,i_n}  \prod _k  \sqrt{ \frac{x}{1-x} }^{i_k}  |i_1,...,i_n\rangle .
\end{eqnarray*}
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 {\it correlated} operator like $Z_1 Z_2$,
we can obtain a second order nonlinear output
\begin{eqnarray}
\langle Z_1 Z_2 \rangle = (2x-1)^2
\end{eqnarray}
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$:
\begin{eqnarray}
\langle \psi | \Lambda_{1,2}(X) Z_1 \Lambda_{1,2}(X) |\psi \rangle  = \langle \psi | Z_1 Z_2 | \psi \rangle.
\end{eqnarray}
In general, an $n$-qubit unitary operation $U$ transforms the observable $Z$
under the conjugation into a linear combination of Pauli operators:
\begin{eqnarray}
U^{\dag} Z_1 U = \sum _{\mb{i}} \alpha _{\mb{i}} P(\mb{i}).
\end{eqnarray}
Thus if you measure the output of the quantum circuit after applying a unitary operation $U$,
\begin{eqnarray}
U V(x) |0\rangle ^{\otimes n},
\end{eqnarray}
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.

\begin{figure}
\centering
\includegraphics[width=80mm]{fig01.pdf}
\caption{The expectation value $\langle Z \rangle$
of the output of a quantum circuit as a function of the input $(x_0,x_1)$.}
\label{fig01}
\end{figure}
To see how the output behaves in a nonlinear way with respect to the input,
in Fig.~\ref{fig01},
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
\begin{eqnarray}
\theta _{2k} =  k \arccos(\sqrt{x_0})
\\
\theta _{2k+1} =  k \arccos(\sqrt{x_1})
\end{eqnarray}
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.

\begin{figure}
\centering
\includegraphics[width=80mm]{fig02.pdf}
\caption{(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.}
\label{fig02}
\end{figure}
Suppose the Pauli $Z$ operator is measured on each qubit as an observable.
Then we have
\begin{eqnarray}
z_i = \langle Z_i \rangle ,
\end{eqnarray}
for each qubit.
In quantum extreme learning machine,
the output is defined by taking linear combination of these $n$ output:
\begin{eqnarray}\label{LR}
y = \sum _{i=1}^{n} w_i z_i .
\end{eqnarray}
Now the linear readout weights $\{ w_i \}$ are tuned
so that the quadratic loss function
\begin{eqnarray}
L = \sum _j (y^{(j)}- \bar y^{(j)})^2
\end{eqnarray}
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. (\ref{LR}) 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 \cite{UAP_goto}.

Here we should note that a similar approach, quantum kernel estimation, has been taken in Ref.~\cite{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') \equiv \langle \Phi (x) | \Phi (x') \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.~\ref{fig02},
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.~\ref{fig02} (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.

\subsection{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.~\cite{QCL,FarhiQNN},
which we call {\it quantum circuit learning}.
Let us consider the same situation with quantum extreme learning machine.
The state before the measurement is given by
\begin{eqnarray}
UV(x) |0\rangle ^{\otimes n}.
\end{eqnarray}
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:
\begin{eqnarray}
U(\{ \phi_k \}) = \prod _{k} u( \phi_k ).
\end{eqnarray}
Thereby, the output from the quantum circuit
with respect to an observable $A$
\begin{eqnarray*}
\langle A (\{ \phi _k\},x) \rangle  =
\langle 0|^{\otimes n} V^{\dag}(x) U(\{ \phi _k \})^{\dag}
Z_i U(\{ \phi _k \})V(x) |0 \rangle ^{\otimes n}
\end{eqnarray*}
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
\begin{eqnarray}
L(\{ \phi _k\}) = \sum _j (\langle A (\{ \phi _k\},x^{(j)}) \rangle - y^{(j)})^2.
\end{eqnarray}
The gradient of the loss function can be obtained as follows:
\begin{eqnarray*}
\frac{\partial}{\partial \phi _l } L(\{ \phi _k\})
&=&\frac{\partial}{\partial \phi _l } \sum _j (\langle A (\{ \phi _k\},x^{(j)}) \rangle - y^{(j)})^2
\\
&=&\sum _j 2(\langle A (\{ \phi _k\},x^{(j)}) \rangle - y^{(j)}) \frac{\partial}{\partial \phi _l } \langle A (\{ \phi _k\},x^{(j)}) \rangle.
\end{eqnarray*}
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
\begin{eqnarray}
u(\phi_k) = W_k e^{-i (\phi _k/2) P_k},
\end{eqnarray}
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$  ~\cite{QCL,Mitarai19}:
\begin{eqnarray*}
&&\frac{\partial}{\partial \phi _l } \langle A (\{ \phi _k\},x^{(j)}) \rangle
\nonumber \\
&=& \frac{1}{2 \sin \epsilon} (
\langle A (\{ \phi _1 ,..., \phi _l + \epsilon, \phi _{l+1},... \},x^{(j)}) \rangle \\
&-& \langle A (\{ \phi _1 ,..., \phi _l - \epsilon, \phi _{l+1},... \},x^{(j)}) \rangle).
\nonumber \\
&&
\end{eqnarray*}
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:
\begin{eqnarray}
 \theta _l ^{(m+1)} = \theta _l ^{(m)}  -  \alpha \frac{\partial }{\partial \phi _l} L(\{ \phi _k\}).
\end{eqnarray}

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~\cite{QCL} and a similar idea~\cite{FarhiQNN},
several researches have been performed with
various types of parameterized quantum circuits~\cite{CircuitCentric,TensorNet,Qneural_Chen,generalized_TensorNet,QClassifier_Du} and various models and types of
machine learning including generative models~\cite{GeneQML_Benedetti,Gene_LiuWang}
and generative adversarial models~\cite{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~\cite{ExpressionPower}.
Experimentally feasible ways to measure an analytical gradient of
the parameterized quantum circuits have been investigated~\cite{Mitarai19,Xanadu_gradient,CalculusPQC}.
An advantage of using such a gradient for the parameter optimization
has been also argued in a simple setting~\cite{HarrowGrad},
while the parameter tuning becomes difficult because of the
vanishing gradient by an exponentially large Hilbert space~\cite{Plateaus}.
Software libraries for optimizing parameterized quantum circuits are now developing
~\cite{Xanadu,VQNet}.
Quantum machine learning on near-term devices,
especially for quantum optical systems,
is proposed in Refs~\cite{Opt_QML,CV_QML}.
Quantum circuit learning with parameterized quantum circuits
has been already experimentally demonstrated on
superconducting qubit systems~\cite{IBM_QML_ex,Rigetti_QML_ex} and
a trapped ion system~\cite{IonQ_QML_ex}.

\subsection{Quantum reservoir computing}
\begin{figure}[t]
\centering
\includegraphics[width=80mm]{fig04.pdf}
\caption{(a) Quantum reservoir computing. (b) Virtual nodes and temporal multiplexing.}
\label{fig04}
\end{figure}
Now we return to the reservoir approach and extend
quantum extreme learning machine from non temporal tasks to temporal ones,
namely, quantum reservoir computing~\cite{QRC}.
We consider a temporal task, which we explained in Sec.~\ref{subsec:temporal}.
The input is given by a time series
$\{ x_k\}_{k}^{L}$ and the purpose is to learn
a nonlinear temporal function:
\begin{eqnarray}
\bar y_k = f(\{ x_j \}_{j}^{k}).
\end{eqnarray}
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.~\ref{subsec:vector_density}.

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,
$\mb{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 {\it hidden} node of the network.
As we seen in Sec.~\ref{subsec:vector_density},
any physical operation can be written as a
linear transformation of the real vector by a $4^N \times 4^N$ matrix $W$:
\begin{eqnarray}
\mb{r}' = W \mb{r}.
\label{eq:time_evolution}
\end{eqnarray}

Now we see, from Eq.~(\ref{eq:time_evolution}), a time evolution
similar to the recurrent neural network,
$\mb{r}' = \tanh( W \mb{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 \mb{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
\begin{eqnarray}
\rho _{x_k} = \frac{I+(2 x_k - 1)Z}{2}.
\end{eqnarray}
For simplicity,
below we consider the case where
only one qubit is replaced for the input.
Corresponding matrix $S_{x_k}$
is given by
\begin{eqnarray*}
(S_{x_k})_{\mb{j}\mb{i}} = {\rm Tr}\left\{ P(\mb{j}) \frac{I+(2x_k-1)Z}{2}
\otimes {\rm Tr}_{\rm replace} [P(\mb{i})] \right\} /2^N,
\end{eqnarray*}
where ${\rm Tr}_{\rm replace}$ indicates
a partial trace with respect to the replaced qubit.
With this definition, we have
\begin{eqnarray}
\rho ' =  {\rm Tr}_{\rm replace}[\rho] \otimes \rho _{x_k}
\Leftrightarrow \mb{r}'  = S_{x_k} \mb{r}.
\end{eqnarray}

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^{-i H \tau}$ for a given time interval $\tau$.
Let us denote its representation on the vector space by $U_\tau$:
\begin{eqnarray}
\rho' = e^{-i H \tau} \rho e^{i H \tau}
\Leftrightarrow \mb{r}' = U_{\tau} \mb{r}.
\end{eqnarray}
Then, a unit time step is written as
an input-depending linear transformation:
\begin{eqnarray}
\mb{r}((k+1)\tau) =  U_\tau S_{x_k}  \mb{r}(k\tau).
\end{eqnarray}
where $\mb{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 {\it true nodes}, is
defined by a $M \times 4^N $ matrix $R$,
\begin{eqnarray}
\bar r_l (k\tau) = \sum _{\mb{i}} R_{l\mb{i}}  r_{\mb{i}}(k\tau).
\end{eqnarray}
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.~\ref{fig04} (a)):
\begin{eqnarray}
y_k = \sum _l W^{\rm out}_l \bar r_l (k\tau),
\end{eqnarray}
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.,
\begin{eqnarray}
\bar r_l = {\rm Tr}[Z_l \rho],
\end{eqnarray}
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~\cite{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.~\ref{fig04} (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.,
\begin{eqnarray}
\mb{r}(k\tau +(v/V)\tau) \equiv U_{(v/V)\tau} S_{x_{k}} \mb{r}(k\tau).
\end{eqnarray}
In total, now we have $N \times V$ nodes,
and the output is defined as their linear combination:
\begin{eqnarray}
y_k = \sum _{l=1}^{N} \sum _{v=1}^{V} W^{\rm out}_{j,v} \bar r_{l}(k\tau +(v/V)\tau).
\end{eqnarray}
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.~\cite{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~\cite{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.~\cite{QRC_OTOC,QRC_sig,QRC_aus}.
Specifically, in Ref.~\cite{QRC_sig},
the authors use quantum reserovir computing to
detect many-body entanglement by estimating nonlinear functions
of deinsity operators like entropy.

\subsection{Emulating chaotic attractors using quantum dynamics}
\begin{figure*}
\centering
\includegraphics[width=140mm]{fig03.pdf}
\caption{Demonstrations of chaotic attractor emulations. (a) Lorenz attractor. (b) Mackey-Glass system. (c) R{\"o}ssler attractor. (d) H{\'e}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.}
\label{fig03}
\end{figure*}
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
\begin{eqnarray}
\bar y_{k} = x_{k+1} = f(\{x_{j}\}_{j=1}^{k}).
\end{eqnarray}
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,
\begin{eqnarray}
\frac{dx}{dt} &=& a (y-x),
\\
\frac{dy}{dt} &=& x (b-z) - y,
\\
\frac{dz}{dt} &=& xy - cz,
\end{eqnarray}
with $(a,b,c)=(10,28,8/3)$,
(ii) the chaotic attractor of Mackey-Glass equation,
\begin{eqnarray}
\frac{d}{dt} x(t) = \beta \frac{x(t-\tau)}{1+x(t-\tau)^n} - \gamma x(t)
\end{eqnarray}
with $(\beta , \gamma , n) = (0.2, 0.1 , 10)$ and $\tau = 17$,
(iii) R{\"o}ssler attoractor,
\begin{eqnarray}
\frac{dx}{dt} &=& -y-z,
\\
\frac{dy}{dt} &=& x+ay,
\\
\frac{dz}{dt} &=& b+z(x-c),
\end{eqnarray}
with $(0.2,0.2,5.7)$,
and (iv)
H{\'e}non map,
\begin{eqnarray}
x_{t+1} = 1-1.4x_t +0.3x_{t-1}.
\end{eqnarray}
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
\begin{eqnarray}
H = \sum _{ij} J_{ij} X_i X_j + h Z_i,
\label{eq_Ising}
\end{eqnarray}
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.~\ref{fig03} for each of (a) Lorenz attractor,
(b) the chaotic attractor of Mackey-Glass system, (c) R{\"o}ssler attractor, and (d) H{\'e}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.

\section{Conclusion and Discussion}
\label{sec:conclusion}
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.

\section*{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.


\section{Introduction}

Prediction of traffic states and travel times evolution is a key component of any traffic monitoring system and decision support system.
Their accurate estimation is critical for freeway managers, especially when the network becomes congested.
This problem has been extensively investigated in the transportation literature using model-based, simulation-based and data-driven approaches \citep{VLAHOGIANNI2005211,mori_review_2015,Wang2018}.
For short-time prediction, model-based and simulation-based approaches use traffic flow models in conjunction with data assimilation techniques such as recursive Bayesian estimators to predict the traffic states and the resulting travel times \citep{vlahogianni_short-term_2014,mori_review_2015,wang_renaissance_2006,kumar_bus_2017}.
Most data-driven approaches use general purpose parameterized mathematical model such as linear regression \citep{rice_simple_2004}, Kalman filtering \citep{van_lint_online_2008,Nanthawichit2003}, particle filters \citep{Wang2007} support vector regression \citep{huang_deep_2014}, random forest, Bayesian network \citep{LI2019}, artificial neural networks \citep{adeli2001,van_lint_online_2008,VLAHOGIANNI2005211,XU2020,LI2017} and many other techniques to capture and learn from data the correlations between traffic variables (speed, travel-time) over space and time \citep{COIFMAN2002351,Polson2017,MA2020}.
As pointed out by \cite{yildirimoglu_experienced_2013}, who wrote a complete and useful state of the art of the estimation methods, these approaches suffer from various limitations.
To quote only a few, a common limitation is the spatio-temporal correlations that are mainly artificially selected \citep{XU2020}.
Beside, some of the existing methods resort to experienced travel times, i.e. travel times calculated by traveling a trajectory through the velocity field.
However, this information is rarely available in real time because experienced travel time is usually greater than the prediction horizon.

Consequently, the purpose of this paper is to the evolution of congestion, therefore, travel times with a simple, fully explainable, and practice-ready method.
The proposed method uses both historical and real time traffic information to calculate short-term congestion and travel time evolution forecast.
To this end, we used the concept of \textit{congestion map} to consider queue propagation rather than traffic states variables evolution, such as density or speed.
Figure \ref{fig:fig0} presents the mechanism of the global algorithm of the proposed method.

The fist step (\textcircled{\tiny{1}}) consists in partitioning historical information into $k$ clusters $C_k$ presenting similar characteristics based on the traffic patterns observed in the highway.
As in \cite{yildirimoglu_experienced_2013}, \cite{SHIVANAGENDRA2003285} and many other papers, we first resort to a Principal Component Analysis to reduce the number of variables.
Then, a Gaussian Mixture Model and a k-means algorithm are used to gather the historical data.
Note that both approaches are used to evaluate the sensitivity of the estimation method to the clustering process.
Then, the proposed method differs from the one of \cite{yildirimoglu_experienced_2013}: rather than considering the global behavior of the clusters, we try to identify which day $d_k$ within the cluster $C_k$ is the most representative of the group.
This so-called \textit{consensual day} is determined based on the congestion maps of the clustered days through a method derived from consensual learning technique \citep{FILKOV2004}.

Once a set of consensual days $D_k$ has been established, the second step of the method (\textcircled{\tiny{2}}) is devoted to real-time application.
Based on the first minutes / hours of a new observation day, measurements are processed in real-time to identify in this set which is the closest consensual day in $D_k$.
The recorded congestion map and the observed speed of this closest consensual day are then used to predict the congestion and travel times evolution of the next minutes / hours of the new observation day.
The main benefit of using past measurements to compute future traffic conditions is to ensure realistic values that are fully consistent with traffic dynamics.

\begin{figure}[H]
    \centering
        \includegraphics[width=\textwidth]{Fig0.pdf}
        \caption{Graphical representation of the proposed method}
        \label{fig:fig0}
\end{figure}

The remainder of the paper is organized as follows: Section 2 presents the case study and the dataset used in the paper, Section 3 introduces the prediction methods, Section 4 is devoted to the clustering process and Section 5 to the analysis of the results, while Section 6 includes a conclusion.

\section{Case study and dataset}

In this paper, we focus on the $M6$ highway near Lyon, France.
Figure \ref{fig:fig1} depicts a sketch of the site.
It is important to notice that this highway is used to access Lyon's city center through a tunnel, which is a recurrent, active bottleneck.
Moreover, this highway is one of the most important in France and favored by holidaymakers because it links Paris to the south of France (French Riviera) and the Alps.
This highway is thus also called \textit{Motorway of the Sun}.
Consequently, major traffic jams are always observed during holidays.

\begin{figure}[H]
    \centering
        \includegraphics[width=\textwidth]{Fig1.png}
        \caption{Sketch of the studied site}
        \label{fig:fig1}
\end{figure}

Regularly staggered loop detectors can be found on this highway section.
These detectors provide average flow, speed and occupancy rate per lane every $1$ minute.
In this study, we mainly focus on data from $9$ detectors of a $7$ km long section, see Figure \ref{fig:fig1}.
Accordingly, $9$ sections of length $\Delta l$ roughly centered on the detectors have been defined.
The maximal authorized speed varies from $90$ km/h to $70$ km/h within this section.
In the remainder of the paper, detectors are labeled by increasing positions from $0$ to $8$.
We consider that traffic conditions between detectors can be interpolated by using observations of the closest detector.
All data from January $2018$ to October $2018$ is available (except May).
We partitioned the data by day from $6$:$00$ to $22$:$00$ ($960$ observations per day).
Finally, data has been roughly cleaned up to remove unrealistic values or problem of acquisition.

\section{Methodology}

\subsection{Congestion maps}
To mainly focus on traffic dynamic rather than speed evolution, we used the concept of \textit{congestion map}.
Indeed, recorded values of speed can vary because of many local phenomenons (such as variations in driving behaviors, noises in measurements, etc.) but without correlation with the macroscopic dynamics of flow that rule traffic conditions.
One potential method, but far from perfect, to reduce this bias is to average the observations.
Here, we decide to use a more drastic method to only focus on two possible traffic states: free-flow and congestion.
For each loop detector $l \in [0,8]$, we consider therefore a variable $x_l$ that, at time $t$, is equal to $1$ if the observed speed $v_l(t)$ is lower than a congested speed threshold $v_{cong}$ (fixed here at $40$ km/h) and equal to $0$ otherwise.
It makes it possible to compute map $M_d$ of day $d$ as Boolean matrix of size "number of detectors" x "number of observations per day" composed of elements $x_l(t)$.
Figure \ref{fig:fig2}a shows the classical speed maps for $6$ randomly selected days of the case study.
Note that the darker the color the lower the speed and that the traffic flows from the top to the bottom of the graphs.
Different congestion patterns clearly appear with different length propagation and duration.
These observation are easier when focusing on the associated congestion maps $M_d$ in Figure \ref{fig:fig2}b.
Note that black color stands for $x_l(t) = 0$ (free-flow) and white color for $x_l(t) = 1$ (congestion).

\begin{figure}[H]
    \centering
        (a)
        \includegraphics[width=\textwidth]{Fig2a.png}
        (b)
        \includegraphics[width=\textwidth]{Fig2b.png}
        \caption{(a) Speed maps for $6$ different days of the historical dataset (traffic flows from 0 to 8), and (b) associated binary congestion maps $M_d$}
        \label{fig:fig2}
\end{figure}

\subsection{Clustering historical data}

The regularity of traffic events makes very useful the information that we can obtain from historical data and what can been learned from past situations.
Consequently, it is worth appealing to classify the different observation days.
Because the size of the initial dataset is very large ($8640$ variables for a single day), the first step, to speed up the clustering methods and obtain accurate results, is to perform a Principal Component Analysis (PCA) to reduce the dimensions of the observations \citep{SHIVANAGENDRA2003285}.
Notice that we consider here speed ${v_l}(t)$ for $l \in [0,8]$ and $t \in [0, 960]$ as the main vector of the PCA.
Since PCA is a usual method, we do not study here the results in details.

Then, the purpose of the second step of the method is to cluster similar days of the historical dataset.
To this end, we use two classical clustering algorithms: (i) the k-means with an Euclidean distance between observations to gather and (ii) a Gaussian Mixture Model (GMM) as proposed in \cite{yildirimoglu_experienced_2013}.

(i) The k-means algorithm is one of the most popular unsupervised learning methods, that aims to gather the data into groups of equal variance that minimize the inertia, i.e., within-cluster sum-of-squares.
With a fast computing time, results are easy to interpret, but k-means implicitly assumes that all clusters are spherical.
This shape could introduce a strong bias, especially for observations of highly non-linear phenomenons.

(ii) In the GMM method, the underlying idea is to consider that similar days constituting the different clusters follow normal distributions.
Consequently, the set of clusters, i.e., the partition, is ruled by a GMM.
This clustering method gives importance to the distribution of the data points and not only the distance between them.
Consequently, GMM is well adapted to our case because it provides clusters that may have different sizes and correlation within them.
Note that a full covariance matrix is assumed in GMM.

\subsection{Identifying consensual days}

Now, the idea is to determine which day of a cluster is the most representative of the group.
As already explained, the motivation is to use the representative days as the prediction for the coming days, see Figure \ref{fig:fig0}.
To this end, we have adapted the consensual learning method used in \cite{lopez_revealing_2017} to our specific case.
Thus, the representative day of a cluster is identified according to a distance based on the Rand index, i.e., the accuracy between congestion maps.
The Rand index between two maps $M_d$ and $M_p$ of days $d$ and $p$ is defined as the number of concomitant results among the total number of observations.
Here, we consider two observations as concomitant when, at a given detector $l$ and time $t$, both observations for day $d$ and $p$ are in the same state (free-flow / free-flow or congested / congested), i.e. $M_d(i,j) = M_p(i,j)$.
The distance that we used can be formulated as:

\begin{equation}
    R_M(d,p) = \frac{a}{a+b+c}
\end{equation}
where:
\begin{itemize}
    \item $a$ is the size of $\{x_l(t), M_d(l,t)-M_p(l,t) = 0 \}$ (free-flow / free-flow or congested / congested);
    \item $b = \vert \{x_l(t), M_d(l,t)-M_p(l,t) = 1\} \vert $ (free-flow / congested);
    \item $c = \vert \{x_l(t), M_d(l,t)-M_p(l,t) = -1\}\vert$ (congested / free-flow);
\end{itemize}

Then, we define the consensual day $d_k$ of a given cluster $C_k$ as the one that maximizes the sum of the Rand indices within a cluster:

\begin{equation}
    d_k = \text{arg}\,\max_{d \in C_k} \{\sum_{p \in C_k} R_M(d,p)\}
\end{equation}

Consequently, the set of consensual days $D_k$ can be determined for the whole historical dataset.

\subsection{Travel time prediction}

We are now going to take advantage of the consensual days to predict both congestion and travel times evolution in real time.
Consider a new day of observation $p$.
The time interval is discretized into periods of $\delta t$ minutes.
At time $t$, we try to determine which consensual day $d_k \in D_k$ is the closest to this new day $p$.
To this end, we only consider the last $\Delta t$ observations and build a partial congestion map $m_p(t-\Delta t,t)$ composed of $x_l(t)$ for $t \in [t-\Delta t,t]$ and $l \in [0,8]$.
This map is compared to the partial maps extracted from the congestion maps of the consensual days $m_{d_k}(t-\Delta t,t) = M_{d_k}(t'), \: \forall t' / t' \in [t-\Delta t,t] $ and $d_k \in D_k$.
The consensual day $d^*_p(t)$ that has the maximal Rand index, based on the partial maps, with the observation $p$ is selected:

\begin{equation}
    d^*_p(t) = \text{arg}\,\max_{d \in d_k} \{R_m(p,d_k)\}
\end{equation}

Now that $d^*_p(t)$ have been identified, we used the historical data $x_l(t)$ and $v_l(t)$ of $d^*_p$ to predict variables of day $p$ for the next time step $\delta t$.
This process is iterated at every $\delta t$, and congestion and speed maps can be built by gathering the data.
Notice that the optimal consensual day $d^*_p$ can (and surely will) change with time $t$.
Duration $\delta t$ (prediction horizon) and $\Delta t$ (learning period) belong to the parameters of the proposed method.
Thus, the prediction of the congestion map and the travel times (respectively) of day $p$ are (respectively):
\begin{equation}
    M^*_p(l,t) = M_{d^*_p(t)}(l,t), \: \forall l \in [0,8]
\end{equation}
and:
\begin{equation}
    v^*(l,t) = v_{d^*_p(t)}(l,t).
\end{equation}
where $M_{d^*_p(t)}(l,t)$ is the element of the observed congestion map of the consensual day $d^*_p(t)$ at time $t$ for detector $l$, and $v_{d^*_p(t)}(l,t)$ is the observed speed at time $t$ on the congestion day $d^*_p(t)$ for detector $l$.

The travel time $\tau$ at time $t$ is calculated as:
\begin{equation}
    \tau(t) = \sum_{l=0}^{9} \frac{\Delta x_{l}}{v_l(t)}
\end{equation}
where $\Delta x_{l}$ is the length of section $l$, see Figure \ref{fig:fig1}.
It is important to notice that it corresponds to the definition of the instantaneous travel times: the travel time at time $t$ is calculated based on the speed at the different detector locations at time $t$ \citep{YEON2008325}.
It may introduce some bias compared to the experience travel times, but, in our case, this bias is very limited because of the relatively short length of the case study (less than $10$ km).

In order to evaluate the global prediction method proposed in the paper, a cross-validation procedure can be used.
A simple holdout method is considered by randomly selecting $75 \%$ of the initial data as the learning set.
This training set is clustered into $K$ groups, and the associated $K$ consensual days are determined.
Then, travel times are predicted for the remaining $25 \%$ and compared with the observations to evaluate and validate the method.
Details are presented in the following section.

\subsection{Comparison with the method presented by \cite{yildirimoglu_experienced_2013}}

As already mentioned, our method share several common features with the work of \cite{yildirimoglu_experienced_2013}.
Consequently, this section proposes a short comparison of the two approaches.

To summarize, the method of \cite{yildirimoglu_experienced_2013} aims to identify groups of days presenting similar traffic conditions to build stochastic congestion maps, identifying the probability that a section of the test case is congested.
According to the newly available information, the stochastic congestion maps are used to revise the state of a priori knowledge.
Their work is organized as follows:
\begin{itemize}
    \item Historical data:
    \begin{itemize}
        \item Groups of similar days are identified with a PCAxGMM combination (speed is the observation vector);
        \item Stochastic congestion maps are produced for each group to highlight blocks of recurrent congestion;
    \end{itemize}
    \item Real-time: the goal is to identify the closest situation in the historical data to determine which bottlenecks are active. Then, potential active bottlenecks are integrated in the processing of real-time measurements to produce experienced travel times.
    \begin{itemize}
        \item Identification of the active blocks based on the comparison of the real-time congestion map with stochastic congestion maps;
        \item Travel times are produced by combining the real-time measurements of the speed at the different detectors and a correction based on the identification of potential bottleneck.
    \end{itemize}
\end{itemize}

In comparison, our method aims to identify groups of days presenting similar traffic conditions to determine a single representative day of each group.
The measurements of these representative days are then used in real-time to predict future traffic conditions.
The newly available information is only used to determine the closest representative days and is directly processed to product prevision of the travel times.
The method can be resumed as follows:
\begin{itemize}
    \item Historical data:
    \begin{itemize}
        \item Groups of similar days are identified with a PCAxGMM / k-means combination (speed in the observation vector);
        \item Identification of the most representative day of each group by consensual learning and congestion maps;
    \end{itemize}
    \item Real-time: the goal is to identify the closest consensual day and use past observations to predict future conditions.
    \begin{itemize}
        \item Identification of the closest consensual day by the comparison of real-time congestion map with the consensual congestion maps;
        \item Travel times are directly those observe for the identified consensual day.
    \end{itemize}
\end{itemize}

\section{Clustering of days with similar traffic conditions into the historical dataset}

\subsection{Determining the optimal number of clusters}

For both approaches, the number of groups $K$ has to be fixed exogenously, i.e., before performing the clustering of the dataset.
Even if there is no definitive answer, several methods exist to determine the optimal number of clusters.
These methods are either based on a criterion minimization/maximization (such as the elbow or averaged Silhouette methods \citep{ROUSSEEUW198753}) or on a statistical test (such as gap statistic method, Bayesian information criterion or Akaike information criterion).
However, these criteria are frequently not consistent between them, and it was the case for our study.
Consequently, we have decided to tailor our metrics to determine the optimal number of clusters $K$ that we need to predict the evolution of congestion and travel times.

The prerequisite is to determine clusters that gather a sufficient number of days (for example, $5\%$ of the test dataset) with similar traffic conditions (intra-cluster homogeneity) and that are different enough from one to the other to justify different groups (inter-cluster dissimilarity).
To this end, three metrics have been developed.
The intra-cluster homogeneity is appraised by calculating the average of the Rand index inside a cluster:
\begin{equation}
    \frac{1}{|C_k|} \sum_{k \in K} \frac{1}{|C_k|-1} \sum_{\forall (p,q) \in C_{k}, p \neq q} R_M(p,q),
\end{equation}
where $K$ is the number of clusters, $C_k$ is the cluster $k$, $d_k$ is the consensual day of cluster $C_k$, $p$ and $q$ are days belonging to $C_k$, and $|C_k|$ is the cardinal of cluster $C_k$.
The inter-cluster dissimilarity is evaluated by determining the average of the Rand index between each pairs of consensual days:
\begin{equation}
    \frac{1}{\binom{n}{2}} \sum_{\forall (p,p') in D_{k}} R_M(p,p').
\end{equation}
Finally, we also compared the number of clusters $n_{min}$ gathering more than $5$ days to the targeted number of cluster $K$.

In order to perform this analysis, the clustering process is iterated $20$ times to ensure the generality of the results.
Then, the averages of the three metrics presented above are calculated.

Figure \ref{fig:fig2b}a compares the intra-cluster homogeneity with inter-cluster dissimilarity in regards to various values of $K$ for both approaches.
It reveals that the inter-cluster dissimilarity is stable for $K > 10$.
Simultaneously, the intra-cluster homogeneity continues to increase but a very low rate.
In the meantime, Figure \ref{fig:fig2b}b depicts $n_{min}$ the number of clusters gathering more than $5$ days in function of $K$.
We also propose the evolution of the Silhouette score with the number of clusters (see Figure \ref{fig:fig2b}c).
It clearly shows that for $K$ bigger than $30$, both methods are unable to identify large clusters.
It is even worse for high values of $K$.
Consequently, it appears that $n_{opt} = 18$ constitutes a good balance to gather the data of our test site (black cross in Figure \ref{fig:fig2b}a, b, and c).
Especially, we can almost observe a local optimum for the average silhouette width at between $15$ and $20$ clusters.
$n_{opt} = 18$ ensures a sufficient number of clusters and a sufficient number of days per cluster (more than $13$ days on average) with an acceptable intra-cluster homogeneity and inter-cluster dissimilarity.
It also means that clusters can be modified or removed retrospectively without changing the metrics significantly.

\begin{figure}[H]
    \centering
        \includegraphics[width=\textwidth]{Fig3.png}
        \caption{(a) Comparison of intra-cluster homogeneity and  inter-cluster dissimilarity, (b) Number of large clusters vs $K$ for the k-means and the GMM methods, and (c) Evolution of the Silhouette score with the number of clusters}.
        \label{fig:fig2b}
\end{figure}

In addition to these averaged metrics, the stability of the clustering results is also studied.
To this end, the Rand index between successive partitions is calculated.
Consequently, if $(C_{n-1})$ is the $n-1$ partition and $(C_n)$ is the $n$ partition, the Rand index can be expressed as follows:
\begin{equation}
    R = \frac{a+b}{a+b+c+d}
\end{equation}
where:
\begin{itemize}
    \item a is the number of pairs of days that are in the same clusters in $(C_{n-1})$ and $(C_n)$;
    \item b is the number of pairs of days that are different clusters in $(C_{n-1})$ and $(C_n)$;
    \item c is the number of pairs of days that are in the same clusters in $(C_{n-1})$ and in different clusters in $(C_n)$;
    \item d is the number of pairs of days that are in different clusters in $(C_{n-1})$ and in the same clusters $(C_n)$.
\end{itemize}

Figure \ref{fig:fig2c} shows the evolution of the Rand index with the number of clusters for the two clustering methods.
From 10 to 50 clusters, less than $5\%$ of days change group when increasing the number of clusters.
It appears that the results are reasonably stable and confirm that $n_{opt} = 18$ is a satisfying trade-off between the different metrics.

\begin{figure}[H]
    \centering
        \includegraphics[width=\textwidth]{Fig_review_4.png}
        \caption{Evolution of the Rand Index between successive partitions}.
        \label{fig:fig2c}
\end{figure}
\subsection{Post-clustering analysis}

First of all, Figure \ref{fig:jour_con}a (b, respectively) shows the congestion maps of the consensual days $D_k$ that have been obtained for the k-means method (GMM, respectively).
Interestingly, the patterns are quite different from one day to the other, for both methods.
Note that $d_k$ have been sorted by decreasing size of $C_k$, i.e. $C_i$ gathers more days than $C_{i+1}$.
Notice that the consensual days of the k-means are denoted $d_i$ whereas those of GMM method are denoted $d'_i$.
Congestion maps of $d_1$ (\ref{fig:jour_con}a) for k-means clustering and of $d'_1$ for GMM method (\ref{fig:jour_con}b) clearly correspond to the cluster of free-flow days.
Situations with morning and evening peak hours can be identified (see, for example, $d_3$ and $d'_2$) or with only a morning peak hour (see, for example, $d_{4}$ and $d'_{15}$).
It is also important to notice that the two methods do not lead exactly to the same consensual days (and the same clusters).
In this case, they only have $7$ consensual days in common.

\begin{figure}[H]
    \centering
        (a) k-means method
        \includegraphics[width=\textwidth]{Fig4b.png}
        (b) GMM method
        \includegraphics[width=\textwidth]{Fig4a.png}
        \caption{Congestion map of consensual days $d_k$ identified with (a) k-means method and (b) GMM method}.
        \label{fig:jour_con}
\end{figure}

Figure \ref{fig:C3} shows the different congestion maps of the days gathered into cluster $C_3$.
Visually, the pattern is similar with always morning and evening peak hours.

\begin{figure}[H]
    \centering
        \includegraphics[width=\textwidth]{Fig5.png}
        \caption{Congestion maps cluster $C_3$ obtained with the k-means method}.
        \label{fig:C3}
\end{figure}

Interpreting the results of the clustering methods is always a difficult task because a relevant balance must be found between a trivial classification and useless and complex justifications of the hypothetical links between observations grouped into the same cluster.
In our case, the number of attributes available to analyze and explain the configuration of the computed clusters is very limited.
Consequently, the post-clustering analysis is mainly focused on the weekday and the month of the observations that are gathered into the same group.

Figure \ref{fig:JC} shows the distributions of both variables.
It turns out that distributions among days (Figure \ref{fig:JC}a and c) and months (Figure \ref{fig:JC}b and d) are very similar.
Note that several iterations have been performed for each method and that the histograms have been normalized.

\begin{figure}[H]
    \centering
        \includegraphics[width=\textwidth]{Fig6.png}
        \caption{Distribution of weekdays and months of identified consensual days for the (a-c) GMM method and (b-d) k-means method}.
        \label{fig:JC}
\end{figure}

To continue the post-clustering analysis, distributions of weekdays and months of the days gathered into the same clusters are now analyzed.
Because the results are very similar for both approaches, the study is now centered on the k-means method.
Let us focus on the six largest clusters of the partition, see Figure \ref{fig:distrib_kmeans} that shows the normalized distributions of the weekdays and the months.
Cluster $C_1$ mainly gathers Saturdays and Sundays, and the month of August, that corresponds to the summer holidays in France.
Without any surprise, the visual inspection of the congestion maps reveals that this group is composed of free-flow situations, see the congestion map of $d_1$ in Figure \ref{fig:jour_con}.
Cluster $C_2$ encompasses the weekdays, mainly during winter.
Furthermore, Cluster $C_3$ is very interesting because it almost only holds Saturdays of holidays: February/March for winter holidays and August for summer holidays.
Cluster $C_5$ gathers the beginning of the week (Mondays and Tuesdays), whereas Cluster $C_6$ the end of the week.
The comparison of the congestion maps of the different days, and $d_5$ and $d_6$ reveals different patterns: the midday hours are more congested for $C_5$.
Finally, cluster $C_4$ groups days with only evening congestion; it may explain the majority of Sundays in the set.

\begin{figure}[H]
    \centering
        \begin{tabular}{ccc}
            \includegraphics[width=0.3\textwidth]{Fig7a.png} & \includegraphics[width=0.3\textwidth]{Fig7b.png} &
            \includegraphics[width=0.3\textwidth]{Fig7c.png} \\
            \includegraphics[width=0.3\textwidth]{Fig7d.png} & \includegraphics[width=0.3\textwidth]{Fig7e.png} &
            \includegraphics[width=0.3\textwidth]{Fig7f.png}
        \end{tabular}
    \caption{Distribution of weekdays and months of the days of the $6$ largest clusters identified with k-means method}.
    \label{fig:distrib_kmeans}
\end{figure}

Because analyzing the results cluster by cluster is a tedious task, a convenient approach is to use an alluvial diagram.
Figure \ref{fig:sankey} shows the contribution of each weekday and month into the different clusters.
The careful analysis of the alluvial diagram completes the previous analysis.
For example, it appears that Sundays are only present in 5 clusters, such as the Saturdays.
On the contrary, each weekday can be found in $10$ different clusters.
A similar analysis can be performed with the contribution of the months to each cluster.

\begin{figure}[H]
    \centering
        \includegraphics[width=\textwidth]{Fig8.png}
        \caption{Alluvial diagram of contribution of days and months into each cluster}.
        \label{fig:sankey}
\end{figure}

These results are very encouraging: the GMM and the k-means methods capture physically similar situations.
Because the two approaches are fully explainable and rely on understandable metrics, all the identified clusters can be easily justified and interpreted.

\section{Results}
The global prediction method is now tested for our case study.
We mainly focus on two metrics: the predicted congestion maps and the travel times series.
As already explained, the data has been divided into a training set ($75\%$ of the initial data) and a validation set (the remaining $25 \%$).
Moreover, we decided to only focus on the k-means method for the sake of brevity in the presentations of the results.

\subsection{Comparison with existing method}

We first decide to compare our method with existing approaches to identify the optimal domain of application.
This task is tough for many reasons.
Notably, it necessitates selecting an existing method that is accurate enough, and which does not require tremendous work to be developed and calibrated.
Moreover, existing methods have different objectives and domains of applications. Consequently, we decided to compare our approach to a naive instantaneous method and a historical average method.
Besides, the comparison with two other versions of our method is performed to convince about the predictive power of our approach.

Therefore, we define the following methods:
\begin{itemize}
    \item $M0$ is a naive approach consisting in shifting the observation made at time $t$ to a horizon $\delta t$. Consequently, this method will be almost perfect for short $\delta t$;
    \item $M1$ is a historical average method consists in calculating average values of the historical data for each day of the week and each time period.
    \item $M2$ is the original method proposed in the paper based (Figure \ref{fig:fig0}) on the clustering of the historical data and the identification of consensual days;
    \item $M3$ is the same as $M2$ except that averaged congestion maps are calculated for each group instead of identifying a consensual day;
    \item $M4$ follows the same process as $M2$ except that historical data are not clustered.
    Each day of the historical dataset can be used as a prediction.
\end{itemize}

To compare these approaches, we use the following metrics:
\begin{itemize}
    \item congestion prediction: Rand indexes, i.e. accuracy, and F1-score between predicted congestion maps and ground truth;
    \item Travel times prediction: RMSE between prevision and observation.
\end{itemize}

The test set is composed of the 56 days that have not been used to perform the clustering (167 days in the historical data set).

\begin{figure}[H]
    \centering
        \includegraphics[width=1\textwidth]{Fig10b.png}
    \caption{Evolution of the boxplots of (a) Rand index, (b) F1-score, (c) MSE according to the horizon (in minutes) between predictions and observations}.
    \label{fig:comp1}
\end{figure}

Figure \ref{fig:comp1}a (respectively b) shows the evolution of the boxplots of the rand index (respectively F1-score) between predicted and observed congestion maps for the four different methods according to the horizon predictions.
It turns out that methods $M2$, $M3$, and $M4$ outperform the instantaneous method $M0$ for horizons bigger than $30$ minutes.
In the same time, these methods are better than the historical approach $M1$ for horizons smaller than $1$ hour.
Figure \ref{fig:comp1}c depicts similar results for the MSE between predicted and observed travel times, except that that $M1$ is almost the worst method between the five tested.
The comparison shed light on the optimal domain of application of our method.
Consequently, a horizon of $\delta t = 60$ minutes is used in the remaining of the paper.

It is also important to notice that methods $M1$, $M2$, and $M3$ have similar performance concerning the congestion propagation prediction.
This is not surprising because they are three variants of the same methods.
Besides, method $M2$ and $M1$ are slightly better than $M3$ for the travel times.
However, we prefer to retain method $M1$ because the clustering makes it more understandable (prediction can be related to qualitative transportation scenario), and the use of a consensual day produces more realistic travel times when focusing on the time series of the predicted travel times.

\subsection{Congestion propagation}
The ability of the proposed method to anticipate congestion propagation is now appraised by comparing the predicted congestion maps with the real observations.
Remember that the prediction method has three parameters: the duration $\Delta t$ that is used to make the prediction, the horizon $\delta t$ for which the prediction is made, and the congested speed threshold (fixed here to $40$ km/h).

Figure \ref{fig:fig5} shows the observed maps, the predicted congestion maps (for $\Delta t = 15$ min and $\delta t = 60$ min), and the difference between these maps for $6$ randomly selected days of the testing sample.
Note that the predicted congestion maps are the result from the repeated previsions at every time step (one minute in our case).
The visual inspection reveals that the differences are very low between prediction and ground truth, and this for significantly different shapes of congestion propagation (Figure \ref{fig:fig5}c, the blue color is when a congestion propagation is predicted, but the observation is free-flow, red for the opposite).
This qualitative analysis is very encouraging.

\begin{figure}[H]
    \centering
        (a) \hspace{2cm} (b) \hspace{2.5cm} (c) \hspace{2cm} (d)
        \includegraphics[width=\textwidth]{Fig9a.png}
        \includegraphics[width=\textwidth]{Fig9b.png}
        \includegraphics[width=\textwidth]{Fig9c.png}
        \includegraphics[width=\textwidth]{Fig9d.png}
        \includegraphics[width=\textwidth]{Fig9e.png}
        \includegraphics[width=\textwidth]{Fig9f.png}
        \caption{Comparison of (a) observed and (b)predicted congestion maps (black is for free-flow - FF conditions, beige for congestion - C Conditions), (c) difference between them (blue is for C / FF and red for FF /C , white and gray for correct predictions) and (d) operational metrics of prediction accuracy (red is for +1, blue -1, black - 2, beige +2 and gray for correct predictions)}
        \label{fig:fig5}
\end{figure}

To go further, we decided to calculate the accuracy of the prevision.
The values are particularly good since they are very close to $100\%$, see titles of Figure \ref{fig:fig5}.
This score might be imputed to the free-flow situations that are the most frequent and easy to predict in the congestion maps.
We also calculate the F1-score of the prevision, that are all satisfying even if this metrics is more sensitive to the wrong prediction.
However, it should be noted that the F1-score is calculated with the congested traffic conditions only.
It may introduces a bias when comparing days with significant different volume of congestion.
However, the forthcoming travel times evolution analysis will confirm that the results of the method are encouraging.

Keeping in mind that the method proposed here is tailored to be practice-ready and answer to operational needs, specific metrics are developed to evaluate the prediction accuracy.
The confidence that can be allowed by an operator to a decision support system, such as a traffic evolution predictor, is often based on simpler indicators than those used in a research paper.
A potential criterion is the decision support system's ability to predict in the same direction what can be observed in real-time by the operator.
For example, if an increase of the congestion is predicted for a horizon $\delta t$, will the operator really observe the congestion propagation in his.her monitoring system a time $\delta t$ later?

Consequently, only predictions of variations are now evaluated.
To this end, the difference $\Delta M_{d} = M_{d}(i+1,j)-M_{d}(i,j)$ can be calculated for all $i \in [1,n-1]$.
This variable is equal to $0$ if queue length remains stable, $1$ if the congestion propagates upstream, and $-1$ if the queue length reduces.
Then, $\Delta M_{d}$ can be compared to the prediction $\Delta M^*_{d}$.
It leads to the following table of the values of $\Delta M_{d}-M^*_{d}$.

\begin{table}[H]
\centering
\begin{tabular}{|c|c|c|c|c|}
\hline
\cellcolor[HTML]{EFEFEF} minus    & \multicolumn{4}{c|}{\textbf{$\Delta M*_{d}(i,j)$}}                                                                                \\ \hline
                             & \cellcolor[HTML]{C0C0C0}\textbf{-}  & \cellcolor[HTML]{C0C0C0}\textbf{0} & \cellcolor[HTML]{C0C0C0}\textbf{1} & \cellcolor[HTML]{C0C0C0}\textbf{-1} \\ \cline{2-5}
                             & \cellcolor[HTML]{C0C0C0}\textbf{0}  & 0                                  & -1                                 & 1                                   \\ \cline{2-5}
                             & \cellcolor[HTML]{C0C0C0}\textbf{1}  & 1                                  & 0                                  & 2                                   \\ \cline{2-5}
\multirow{-4}{*}{\textbf{$\Delta M_{d}(i,j)$}} & \cellcolor[HTML]{C0C0C0}\textbf{-1} & -1                                 & -2                                 & 0                                   \\ \hline
\end{tabular}
\caption{Comparison of the evolution of predicted and observed congestion lengths based on the congestion maps}
\label{tab:tab1}
\end{table}

Figure \ref{fig:fig5}d shows the colormap of $\Delta M_{d}- \Delta M^*_{d}$ for the $6$ random selected day.
It clearly reveals that only a few errors about the evolution direction are made by the method.
Especially if we compare to the rough differences of the prediction and the observation (Figure \ref{fig:fig5}d), we do not accumulate the errors.
In addition to this qualitative analysis, a global indicator can be computed: the ratio $\rho$ of accurate prediction over the total number of observations.

\begin{equation}
    \rho(d) = \frac{card \Big ( \{ \Delta M_{d}- \Delta M^*_{d}= 0 \} \Big ) }{n-1}
\end{equation}

Distribution of $\rho$ for the whole testing set is not shown because it is completely centered on the mean value ($\bar{\rho} = 98.2 \%$) with a minimal standard deviation ($0.4 \%$).
Therefore, the accuracy in the prediction of the evolution of the traffic conditions is excellent.
It is very encouraging for the potential integration of the method into an operational decision support system.
This will be the case for our test study.

\subsection{Travel times estimation}

As already mentioned, the proposed method can predict travel times by using the observed speeds of the consensual days.
The main benefit of this approach is to produce realistic travel times because predictions come from past observations.

Figure \ref{fig:fig6}a shows the travel times series for the $6$ previously studied days.
The orange curves correspond to the prediction, whereas the blue ones are the observations.
Note that a horizon $\delta t = 60$ min and a learning duration of $\Delta t = 15$ min have been used.

\begin{figure}[H]
    \centering
        \includegraphics[width=\textwidth]{Fig10.png}
        \caption{(a) Comparison of observed (in blue) and predicted travel times (in orange), (b) distribution of the absolute errors, (c) evolution of optimal consensual days along the prediction}
        \label{fig:fig6}
\end{figure}

Visually, the results are very good, and the trends of traffic evolution are well predicted.
The method can accurately capture the variations of travel times.
Besides, we have calculated the normalized RMSE between observations and predictions, and the criterions appear to be satisfying.
To complete the analysis, Figure \ref{fig:fig6}b shows the distributions of the errors between predictions and observations.
The proportion of the predictions that have a precision below $2$ minute stands from $40\%$ to $79\%$.
If we focus on the $3$ minutes window, this value increases to $68 \%$ to $87 \%$.
Finally, Figure \ref{fig:fig6}c highlights the variation of the optimal consensual day determined by the prevision method.
It shows out that a large subset of the consensual days is used to predict travel times throughout the day.

\section{Conclusion}

Prediction of congestion propagation and travel times evolution is still a topic actively studied in the literature.
This paper tries to make its contribution by proposing a simple method, which has the main benefit of being fully explainable compared to recent approaches based on machine learning or artificial intelligence methods.

The key component is the concept of \textit{congestion map}, which is a binary observation metric of traffic states on a highway.
Using this proxy reinforces the importance of focusing on traffic dynamics rather than on the numeric values of observation variables.
Thus, the historical dataset of a French $10$ km long freeway is classified into groups of days presenting similar traffic states.
The second innovation of the proposed method is to identify a \textit{consensual day} for each cluster.
According to a distance based on the Rand index, the objective is to determine the day that is the most representative of the cluster.
It makes it possible to find almost the expected traffic situations of this highway: morning and evening peak hours, only morning/evening peak hours, all day long congestion, free-flow day, holiday traffic, etc.

Once those consensual days have been determined, the method can be applied in real-time to predict congestion propagation and travel times evolution.
According to a real-time learning period, observations of a new day are compared to the consensual days' congestion maps.
The closest one is identified, and the congestion map and speeds that have been observed for this specific day are used to predict the behavior of the new day for a given horizon.
This very simple method gives encouraging results for both congestion and travel times evolution.
Especially, the comparison with naive methods (instantaneous and historical average) reveals that the proposed model is useful for longer prediction horizons.

Various future directions can be pursued.
The methodology can be improved by identifying the best duration to compare congestion maps, the accurate prediction horizon, and the congested speed threshold.
To perform the study, we naively tested different values, but a sensitivity analysis could be conducted.
A second improvement could be to use congestion maps that are no more dependent on the time of the day.
The idea is to focus only on the shape of the maps, i.e., shockwaves profiles.
These claims still need to be researched and validated.

\newpage
\section*{Acknowledgments}

This research work was carried out as part of the “Lyon Covoiturage Experimentation – LCE” project at the Technological Research Institute SystemX, and was supported with public funding within the scope of the French Program “Investissements d’Avenir”.

The authors thank the different members of the IRT-Sytem X for the fruitful discussion about the data clustering and the Metropole de Lyon, which provide the data for this study.
Plenty of other interesting datasets can be found on: https://data.grandlyon.com/

This paper was mainly written during the long night of the Covid-19 lockdown.

\newpage
\section*{List of notations}
\begin{table}[H]
\setcellgapes{2.5pt}
\renewcommand*{\arraystretch}{1}
\makegapedcells
\begin{tabular}{ | c  l  l | }
  \hhline{|===|}
  Variable			& Description 							&\\
  \hhline{|===|}
  $t$ 			& Time of the day					& \\
  $l$ & Loop detector index &\\
  $\Delta x_l$ & Length of section $l$ & \\
  $x_l(t)$ & Binary variable & \\
  $v_l(t)$ & Speed recorded at detector $l$ at time $t$ & \\
  $v_{cong}$ & Congestion speed threshold & \\
  $M_d$ & Boolean matrix of congestion map of day $d$ such as $M_d(l,t) = x_l(t)$ & \\
  $R_M(d,p)$ & Rand index between congestion maps of day $d$ and $p$ & \\
  $C_k$ & Cluster $k$ & \\
  $K$ & Number of clusters & \\
  $d_k$ & consensual day of cluster $C_k$ & \\
  $D_k$ & Set of consensual days & \\
  $\Delta t$ & Learning period & \\
  $\delta t$ & Horizon of prediction & \\
  $d^*_p(t)$ & Optimal consensual day for a new observation day $p$ at time $t$ & \\
  \hhline{|===|}
\end{tabular}
\label{notations}
\nomakegapedcells
\end{table}
\newpage