[1]Peter J. Ehlers

[1]Daniel Soh

[1]Wyant College of Optical Sciences, University of Arizona, Tuscon, AZ, US

2]School of Electrical Engineering and Telecommunications, University of New
South Wales, Sydney, Australia

# Stochastic Reservoir Computers

<EMAIL_ADDRESS>Hendra I. Nurdin<EMAIL_ADDRESS>*  [

###### Abstract

Reservoir computing is a form of machine learning that utilizes nonlinear
dynamical systems to perform complex tasks in a cost-effective manner when
compared to typical neural networks. Many recent advancements in reservoir
computing, in particular quantum reservoir computing, make use of reservoirs
that are inherently stochastic. However, the theoretical justification for
using these systems has not yet been well established. In this paper, we
investigate the universality of stochastic reservoir computers, in which we
use a stochastic system for reservoir computing using the probabilities of
each reservoir state as the readout instead of the states themselves. In
stochastic reservoir computing, the number of distinct states of the entire
reservoir computer can potentially scale exponentially with the size of the
reservoir hardware, offering the advantage of compact device size. We prove
that classes of stochastic echo state networks, and therefore the class of all
stochastic reservoir computers, are universal approximating classes. We also
investigate the performance of two practical examples of stochastic reservoir
computers in classification and chaotic time series prediction. While shot
noise is a limiting factor in the performance of stochastic reservoir
computing, we show significantly improved performance compared to a
deterministic reservoir computer with similar hardware in cases where the
effects of noise are small.

Reservoir computing is a form of machine learning in which inputs are
sequentially fed into a nonlinear dynamical system, whose state is measured
and processed using an affine transformation that is trained to produce an
output that closely matches some target sequence of interest. This process of
feeding inputs into a system and extracting readouts is iterative, making
reservoir computers (RCs) particularly well-suited to time-series modeling
tasks [1, 2, 3, 4, 5, 6, 7], and has also been successful for classification
tasks [8, 9, 10, 11] among many possible uses [12]. The training performed on
the readouts is a simple linear regression that can be solved exactly, giving
RCs a relatively low computational cost in comparison to the neural networks
in use today. The reason that simple linear regression suffices for complex
time series modeling is that the nonlinear dynamical system provides readouts
that are informationally rich in the inputs, realizing a nonlinear
transformation of not only the current input but recent neighboring inputs as
well. Using a naturally occurring physical system for the reservoir is a
promising avenue toward fast and efficient computation of difficult tasks, as
the natural dynamics of the system would perform the nonlinear transformations
automatically, leaving only linear computations in the final output readout
for a standard computer to perform.

The primary limitation of the reservoir computing approach is that a
particular reservoir computing system is only suitable for a specific subset
of problems. This limitation stems from the fact that the transformations
applied by the reservoir are fixed. Consequently, a given reservoir will
consistently process inputs in the same manner, which restricts the variety of
potential outputs it can generate. Nevertheless, classes of RCs have been
shown to be universal approximating classes, such as state-affine systems [13]
and echo state networks (ESNs) [14], among others [15, 16, 17, 18, 19, 20, 21,
22, 23]. However, achieving an exactly accurate reservoir computing system may
be impractical as it may require a large reservoir along with an extreme level
of numerical precision. Consequently, one often must balance approximation
accuracy with the practical limitations on size.

The principal finding of this paper is that classes of stochastic echo state
networks — a type of stochastic RC — constitute universal approximating
classes. Universality is an important property for classes of RCs as it
guarantees that any task we want to solve can be approximated to arbitrary
precision using the architecture that defines the class. Without the
universality property, choosing a specific reservoir architecture in the hope
of approximating a certain target function could result in complete failure.
We provide a rigorously detailed mathematical proof that demonstrates
universality in a promising scalable architecture — the stochastic reservoir
computer. This is a critical milestone essential for nurturing and investing
in an important scalable class of RCs. As outlined in the next section, our
definition of a stochastic RC deviates from the conventional reservoir
computing model by utilizing the probabilities of each observed reservoir
state for readouts, rather than the states themselves. This modification
necessitates our detailed analysis of these probabilities to confirm that they
are complex, nonlinear functions of the inputs, requiring only an affine
transformation for universal approximation.

The use of stochastic RCs holds significant potential because the reservoir
can be designed such that the number of computational nodes and, therefore,
the possible outcomes scale exponentially with the hardware size. For example,
in a qubit-based quantum computer, one can perform a unitary rotation of a set
of $n_{q}$ qubits based on the input at a given time step, then measure the
qubits in the computational basis. By running this process many times, we can
obtain an estimate of the probability of measuring each of the $2^{n_{q}}$
distinct bit strings and use them as readouts for training to get an output.
The requirement that the stochastic RC must be run many times to get an
estimate of the probabilities is a potential setback, representing a tradeoff
between more computational nodes and machine runtime. However, this approach
becomes valuable in scenarios where a deterministic computer lacks sufficient
computational nodes to perform the necessary calculations. In such cases, a
stochastic machine can achieve an equivalent number of computational nodes
using a much smaller device. We examine two different designs for a stochastic
RC more closely in the Results section, where we define the exponential
scaling more concretely.

### Stochastic Reservoir Computers

In this section, we will give our precise definition of a stochastic RC. There
are two primary differences between this class of RCs and the deterministic
RCs that are well established in the literature. First, the stochastic
reservoir gives one of $M$ possible outcomes in
$\mathcal{X}=\\{x^{(1)},\dots,x^{(M)}\\}$ at each time step $k$. Second,
instead of using the outcomes directly as readouts for training to acquire an
output, we use the probabilities $P_{k,a}$ of observing a distinct outcome
$x^{(a)}\in\mathcal{X}$ labeled by an index $a\in\\{1,\dots,M\\}$ at each time
step $k$ as readouts instead. Since in a deterministic RC the measured state
of the system $x_{k}$ also serves as the readout used for training, we need to
introduce language to distinguish these two concepts. For stochastic RCs, we
will refer to the measured state $x_{k}$ as the outcome of the reservoir,
while the probabilities $P_{k,a}$ of observing an outcome $x^{(a)}$ are called
the readouts, as these are what we ultimately use for training.

The reason this construction is useful is because it is easy to achieve a
large number of outcomes $M$ by stringing together a small number of
stochastic systems with only a few outcomes individually and reading the
outcomes out as a vector $x_{k}$ at each time step. For example, a binary
stochastic process (a $p$-bit [24]) has only two possible outcomes, but $L$ of
these together have $M=2^{L}$ distinct outcomes, and in general a vector
outcome $x$ whose components each have $m$ possible outcomes will have as an
ensemble $M=m^{L}$ distinct outcomes. This makes it easy to generate a RC with
a large number of computational nodes using a fairly small device, as we will
illustrate in the Results section.

Mathematically, the action of the stochastic reservoir is described by

$\displaystyle x_{k+1}$ $\displaystyle=F(x_{k},u_{k})$ (1) $\displaystyle
F(x,u)$ $\displaystyle\in\mathcal{X}=\\{x^{(1)},\dots,x^{(M)}\\},$ (2)

Thus the outcome of a reservoir measurement at time step $k$ is restricted to
a specific set of $M$ possible vectors of length $L$, each with its own
probability of occurring. The probabilities involved in this process are
defined as

$\displaystyle P_{k,a}$ $\displaystyle\equiv P(x_{k}=x^{(a)})$ (3)
$\displaystyle p_{ab}(u)$ $\displaystyle\equiv P(F(x,u)=x^{(a)}|x=x^{(b)}).$
(4)

In other words, $P_{k,a}$ is the probability at time step $k$ of measuring
outcome $x^{(a)}$, while $p_{ab}(u)$ is the conditional probability of
measuring $F(x,u)=x^{(a)}$ with an input $u$ given that the starting outcome
was $x^{(b)}$. Note that the conditional probability has no dependence on $k$
as it describes the effect that the reservoir has on the probabilities of each
outcome, and the reservoir dynamics do not change over time.

Using these probabilities, we can then write a recursive definition of the
probabilities $P_{k,a}$ using the definition of conditional probability to get

$\displaystyle P_{k+1,a}$ $\displaystyle=P(F(x_{k},u_{k})=x^{(a)})$ (5)
$\displaystyle=\sum_{b=1}^{M}P(F(x_{k},u_{k})=x^{(a)}|x_{k}=x^{(b)})P(x_{k}=x^{(b)})$
(6) $\displaystyle=\sum_{b=1}^{M}p_{ab}(u)P_{k,b}$ (7)

This can be written more succinctly in matrix form as $P_{k+1}=p(u_{k})P_{k}$.
This means that the probabilities themselves form a Markov chain [25], in
particular a controlled discrete-time Markov chain [26] in which the input
sequence $(u_{k})_{k\in\mathbb{Z}}=\\{u_{k}\left|k\in\mathbb{Z}\right.\\}$
acts as the fixed control sequence for the Markov chain. The controlled
transition probability matrix $p(u)$ whose elements are defined by (4) is
determined by the dynamical reservoir. This means that, rather than using the
finite number of outcomes in the set $\mathcal{X}$ as readouts, we can instead
use the probability vectors $P_{k}$ as our readouts and train those to fit a
target series. Thus the controlled Markov chain forms the basis of stochastic
reservoir computing, and we can then write the RC equations as

$\displaystyle P_{k+1}$ $\displaystyle=p(u_{k})P_{k}$ (8)
$\displaystyle\hat{y}_{k}$ $\displaystyle=W^{\top}P_{k},$ (9)

where $\hat{y}_{k}$ is our estimate for some target value $y_{k}$ which is fit
by training the $M$-dimensional weight vector $W$.

In practice, a large stochastic RC that would be suitable for finding good
approximations to difficult problems will have an enormous state space that
will make an analytic calculation of the probabilities $P_{k}$ virtually
impossible. Instead, we can obtain estimates for $P_{k}$ by performing many
independent runs of the stochastic RC through the training and test data and
averaging over the results. In the Results section, we will examine the
effects of estimating $P_{k,a}$ using a finite number of measurements versus
using the exact probabilities.

### Stochastic Echo State Networks

Our proof of universality will focus on a specific type of stochastic
reservoir computing based off of the echo state network (ESN). The action of
the stochastic ESN is described by

$\displaystyle z_{k}$ $\displaystyle=Ax_{k}+Bu_{k}$ (10) $\displaystyle
x_{k+1,i}$ $\displaystyle=\Gamma(z_{k,i}),$ (11)

where $\Gamma(\zeta)\in\\{\xi^{(1)},\dots,\xi^{(m)}\\}$ is a scalar stochastic
activation function with $m$ outcomes that is a function of some scalar input
$\zeta$. The probabilities then evolve via a controlled Markov process
according to Eq. (8), and the outputs are derived from the probabilities as in
Eq. (9). The elements of the controlled transition matrix associated with the
ESN are defined by

$\displaystyle p_{a}(z)$ $\displaystyle=\prod_{i=1}^{L}\varrho_{a_{i}}(z_{i})$
(12) $\displaystyle p_{ab}(u)$ $\displaystyle=p_{a}(Ax^{(b)}+Bu),$ (13)

where $\varrho_{a_{i}}(z_{i})=P(\Gamma(z_{i})=\xi^{(a_{i})})$ is the
probability of measuring $\xi^{(a_{i})}\equiv x^{(a)}_{i}$ given the $i$th
component of the vector $z$, and $L$ is the dimension of each $x_{k}$ in the
ESN. Each $a_{i}$ is an integer between 1 and $m$ such that there is some
bijective mapping between the set
$\\{a_{i}\left|i\in\\{1,\dots,L\\}\right.\\}$ and the integers
$a\in\\{1,\dots,M=m^{L}\\}$. $\varrho(\zeta)$ is to be interpreted as the
probability distribution associated with the scalar stochastic function
$\Gamma(\zeta)$. We now have a controlled composite probability distribution
$p_{a}(z)$ that determines the probability of measuring outcome
$a\in\\{1,\dots,M=m^{L}\\}$ given a vector-valued input $z$, which decomposes
into a product of independent controlled probabilities
$\varrho_{a_{i}}(z_{i})$ associated with each component of both $z$ and the
measurement outcome $x^{(a)}$. Note that since we always work with bounded
inputs and because there are finitely many possible values of $x$, every
component of $z=Ax+Bu$ will also be bounded for any given $A$ and $B$.

## Results

### Universality of Stochastic ESNs

Our main result is a rigorous proof of the universality of stochastic echo
state networks. More precisely, we have proven that a class of stochastic echo
state networks that share a common activation function form a universal
approximating class, provided that the stochastic behavior of the activation
function has certain properties. This universality proof guarantees that,
given a specific activation function that meets the criteria for universality,
there will be a stochastic ESN that can be used to approximate any given task
to arbitrary precision. In practical terms, this means that if an efficient
design for implementing these stochastic ESNs were found and could be scaled
up relatively easily, then that design would be sufficient for building an ESN
that will successfully approximate any problem of interest. Two possible
candidates for such designs are presented and analysed in the following
sections.

The formal statement of our theorem establishing the universality of
stochastic ESNs is as follows:

######  Theorem 1 (Universality).

With the set of uniformly bounded sequences
$K_{R_{u}}\subset(\mathbb{R}^{n})^{\mathbb{Z}_{-}}$ and a weighted metric
$||\cdot||_{w}$ defined in Thm. 2, let $\mathcal{G}_{w}^{\varrho}$ be the
class of functionals generated by stochastic ESNs defined by a controlled
probability distribution $\varrho(\zeta)$ that is defined on the interval
$\zeta\in[-R_{\zeta},R_{\zeta}]$ for some $R_{\zeta}>0$, is continuous in
$\zeta$, and satisfies $\varrho^{\top}(\zeta_{1})\varrho(\zeta_{2})>0$ for all
$\zeta_{1},\zeta_{2}\in[-R_{\zeta},R_{\zeta}]$. If, for the invertible
subspace $\mathcal{V}_{R_{u}}^{p}$ with respect to the controlled transition
matrix $p(u)$ corresponding to this class of ESNs, there exists a vector
$V\in\mathcal{V}_{R_{u}}^{p}$ for which $V^{\top}\varrho(\zeta)$ is strictly
monotonic on the interval $[-R_{\mathrm{mono}},R_{\mathrm{mono}}]$ for some
$R_{\mathrm{mono}}>0$, then the class $\mathcal{G}_{w}^{\varrho}$ is dense in
the space of functionals over the compact metric space
$(K_{R_{u}},||\cdot||_{w})$.

A detailed proof of this theorem and a definition of an invertible subspace
with respect to a controlled transition matrix $p(u)$ is given in the Methods
section. We also give as a corollary to Theorem 1 the universality statement
for the specific case of $m=2$, which has fewer conditions and is relevant to
both the numerical examples in this work as well as for broader applications
such as in quantum computers, as well as another corollary for the
universality of a broader class of stochastic RCs.

### Numerical Examples

In the next couple of sections we will demonstrate the convergence and
performance of stochastic ESNs using a couple of examples. In both examples,
we use a stochastic activation function $\Gamma(\zeta)$ that outputs either
$0$ or $1$ with controlled probabilities $\varrho_{0}(\zeta)$ and
$\varrho_{1}(\zeta)=1-\varrho_{0}(\zeta)$ that have the required properties
which ensures that the class $\mathcal{G}_{w}^{\varrho}$ of stochastic ESNs
using $\varrho(\zeta)$ is a universal approximating class. The designs of the
two examples are also meant to emphasize the practicality of implementing a
stochastic RC using readily available hardware. In particular, they can be
scaled up in a relatively easy way to produce an exponentially scaling number
of computational nodes, with the drawback of potentially needing a larger
number of measurements to estimate the probabilities to the required level of
precision.

To test the performance of these stochastic ESNs, we will use two different
tasks, the first of which is the Sine-Square wave identification task [27].
This task requires the ESN to identify whether a given section of the input is
a sine wave or a square wave, yielding a output of $\hat{y}_{k}<0.5$ for a
sine wave and $\hat{y}_{k}>0.5$ for a square wave. The other task we will use
in this work is the Lorenz $X$ task. The Lorenz $X$ task requires the ESN to
approximate a chaotic dynamical system and predict what the sequence will be 1
time step into the future. In other words, using the input sequence
$(u_{k})_{k\in\mathbb{Z}}=(y_{k-1})_{k\in\mathbb{Z}}$, we want the ESN to
successfully predict $(y_{k})_{k\in\mathbb{Z}}$ in this task. Further details
about these tasks and the simulation parameters are given in Supplementary
Information A.3.

#### Qubit Reservoir Network

Figure 1: Diagram of the design for the qubit reservoir network.

The qubit reservoir network is an example of a stochastic echo state network
that utilizes simple, single-qubit gate operations applied to an array of
qubits to generate a nonlinear stochastic activation function. This is a new
hybrid quantum-classical approach that capitalizes on the ease of performing
linear operators on classical computers while also leveraging a quantum device
for rapid nonlinear transformations, similar to a variational quantum
algorithm. The routine for each time step of the ESN is illustrated in Fig. 1,
where the linear transformation is applied classically, and the $z$ values are
passed onto each qubit using a transformation operator of the form
$U(\theta)=e^{-i\sigma_{x}\theta}$. First, at each time step $k$ we use the
most recent array of qubit measurements $x_{k}$ and the current input $u_{k}$
to calculate $z_{k}=Ax_{k}+Bu_{k}$ on a classical computer. Then, we apply the
transformation operator $U(z_{k,i}+\frac{\pi}{2}x_{k,i}+\delta)$ to each
qubit, where $\delta$ is some constant shift parameter. The extra factor of
$\frac{\pi}{2}x_{k,i}$ resets the qubits to $|0\rangle$ from the states
$|x_{k,i}\rangle$ they collapsed to after the previous measurement. Finally,
each qubit is then measured in the computational basis, which is recorded as
$x_{k+1}$ and used immediately in the next time step. We then repeat this
cycle many times and record the number of times each distinct measurement
outcome $x^{(a)}$ appears at each time step $k$ to obtain estimates of their
probabilities $P_{k,a}$ of occurring. These probabilities are then used as the
readout vector for training and obtaining estimates
$\hat{y}_{k}=W^{\top}P_{k}$ of the target output $y_{k}$.

For this class of stochastic ESNs, the scalar activation function
$\Gamma(\zeta)$ is the measurement of a single qubit in the computational
basis after applying $U(z_{k,i}+\frac{\pi}{2}x_{k,i}+\delta)$. The conditional
probabilities associated with $\Gamma$ for a single node ESN are given by

$\displaystyle P(\Gamma=\xi^{\prime}|\xi,u)$
$\displaystyle=|\langle\xi^{\prime}|U(A\xi+Bu+\frac{\pi}{2}\xi+\delta)|\xi\rangle|^{2}$
(14)
$\displaystyle=|\langle\xi^{\prime}|U(A\xi+Bu+\delta)e^{-i\frac{\pi}{2}\sigma_{x}\xi}|\xi\rangle|^{2}$
(15)
$\displaystyle=|\langle\xi^{\prime}|U(Ax+Bu+\delta)(-i\sigma_{x})^{\xi}|\xi\rangle|^{2}$
(16) $\displaystyle=|\langle\xi^{\prime}|U(A\xi+Bu+\delta)|0\rangle|^{2},$
(17)

where in these equations $A$ and $B$ are scalars and we have used that
$\xi=\\{0,1\\}$ to simplify
$(\sigma_{x})^{0}|0\rangle=(\sigma_{x})^{1}|1\rangle=|0\rangle$. This has the
form $P(\Gamma=\xi^{(a)}|\xi,u)=\varrho_{a}(\zeta)$ with $\zeta=A\xi+Bu$, so
this is a stochastic echo state network with the controlled probability
distribution

$\displaystyle\varrho_{a}(\zeta)=|\langle a|U(\zeta+\delta)|0\rangle|^{2}$
(18)

with $a=0,1$. Since $\varrho_{1}(\zeta)=1-\varrho_{0}(\zeta)$ for any $\zeta$,
the system is fully specified by the controlled scalar probability

$\displaystyle\varrho_{1}(\zeta)=|\langle
1|U(\zeta+\delta)|0\rangle|^{2}=|\langle
1|e^{-i(\zeta+\delta)\sigma_{x}}|0\rangle|^{2}=\sin^{2}(\zeta+\delta).$ (19)

This class of stochastic ESNs can be shown to be universal when
$\delta=\frac{\pi}{4}$ and when the matrix $A$ and vector $B$ are chosen so
that all elements of $z_{k}=Ax_{k}+Bu_{k}$ satisfy $|z_{k,i}|<\frac{\pi}{4}$
at every time step $k$. We show this in Supplementary Information A.1.

Figure 2: The left plot shows average NMSE values for the Lorenz $X$ task,
while the right plot shows the average error percentage for the Sine-Square
wave identification task, both as a function of the number of detectors using
the qubit reservoir network. The blue dots correspond to the deterministic
ESN, the orange dots correspond to the stochastic ESN averaged over 100000
runs for the Lorenz $X$ task and 10000 runs for the Sine-Square wave task, and
the green dots correspond to the stochastic ESN results in the theoretical
limit of infinite statistics, using the exact probabilities. Each data point
represents the average over 100 random choices of $A$ and $B$ for the Lorenz
$X$ task and 1000 choices for the Sine-Square wave task, and the error bars
give the standard deviations of each data point over the samples of ESNs with
different choices of $A$ and $B$.

The left plot of Fig. 2 shows a comparison of NMSE values for the Lorenz $X$
task averaged over 100 different choices of $A$ matrices and $B$ vectors for
the ESN, while the right plot shows a comparison of the percentage of errors
recorded for the Sine-Square wave identification task averaged over 1000
choices. We see in the figure that the stochastic ESN using the exact
probabilities performs much better than the deterministic ESN using roughly
the same hardware. With finite statistics, the stochastic ESN still performs
well for a small number of detectors, but beyond 4 detectors the values of the
error measure seems to level off, and the deterministic ESN begins to
outperform it on average by 4 detectors in the Lorenz $X$ task and at 7
detectors in the Sine-Square wave task.

An explanation for this can be found by analysing how the linear regression
behaves under the noise introduced by the stochastic reservoir, which is
detailed in Supplementary Information A.2. The takeaway is that if
$N_{\mathrm{runs}}$ is not large enough to suppress the shot noise represented
by $P_{k}^{\mathrm{diag}}$ relative to all of the eigenvalues of $K$, then the
performance of the stochastic reservoir will be reduced compared to what it
could be with the exact probabilities.

#### Stochastic Optical Network

Figure 3: Diagram of the design for the stochastic optical network.

The second example of a stochastic echo state network is an optical setup
based off of the optical neural network proposed in Ref. [28]. The novel
approach in their design is the use of single photon detectors (SPDs) in the
hidden layer of their neural network. We use this concept in the design of a
stochastic ESN, where very low intensity lasers have their photon number
measured using SPDs, generating a stochastic activation function. In this
work, we focus exclusively on results using coherent laser light, but
incoherent beams may be used as well, as in the main result of Ref. [28]. In
contrast to the previous example, this design is almost entirely classical
aside from the quantized measurement of light and can be built today only
using well-established hardware that has existed for decades. This
demonstrates that while stochastic RCs have a deep connection to quantum
computing due to the stochastic nature of quantum measurement, they are not
inherently quantum devices.

The design of the device is shown in Fig. 3. Here, an array of very low
intensity coherent lasers which are either turned on or off at each time step
$k$ encodes the vector $x_{k}$, where $x_{k,i}\in\\{0,1\\}$ for every $k$ and
$i$. This laser array is passed through a matrix-vector multiplier and joined
with another laser array whose intensities are proportional to the input
$u_{k}$. This applies the linear part of the reservoir action to produce
$z_{k}=Ax_{k}+Bu_{k}$. Then, each laser in the array is measured with a single
photon detector, which enacts the stochastic activation function
$x_{k+1,i}=\Gamma(z_{k,i})$. In this setup, the SPD only measures whether or
not any light is detected. For coherent light, the probability that the vacuum
state is measured is given by $P(\Gamma=0|\alpha)=e^{-\alpha^{*}\alpha}$ for
the coherent state parameter $\alpha$, and therefore
$P(\Gamma=1|\alpha)=1-e^{-\alpha^{*}\alpha}$. If we include a constant offset
$d$ along with the input beams, then the controlled probability distribution
associated with $\Gamma$ will be

$\displaystyle\varrho_{1}(\zeta)=1-e^{-(\zeta+d)^{2}}$ (20)

This class of stochastic ESNs can be shown to be universal with $\delta=1$ and
when the matrix $A$ and vector $B$ are chosen so that all elements of
$z_{k}=Ax_{k}+Bu_{k}$ satisfy $|z_{k,i}|<1$ at every time step $k$. We show
this in Supplementary Information A.1.

Figure 4: The left plot shows average NMSE values for the Lorenz $X$ task,
while the right plot shows the average error percentage for the Sine-Square
wave identification task, both as a function of the number of detectors using
the stochastic optical network. The blue dots correspond to the deterministic
ESN, the orange dots correspond to the stochastic ESN averaged over 100000
runs for the Lorenz $X$ task and 10000 runs for the Sine-Square wave task, and
the green dots correspond to the stochastic ESN results in the theoretical
limit of infinite statistics, using the exact probabilities. Each data point
represents the average over 100 random choices of $A$ and $B$ for the Lorenz
$X$ task and 1000 choices for the Sine-Square wave task, and the error bars
give the standard deviations of each data point over the samples of ESNs with
different choices of $A$ and $B$.

The left plot of Fig. 4 shows a comparison of the NMSE values associated with
the Lorenz $X$ task, averaged over 100 different choices of $A$ matrices and
$B$ vectors for the ESN, while the right plot shows a comparison of the
percentage of errors recorded for the Sine-Square wave identification task
averaged over 1000 choices. We see in the figure that the stochastic ESN using
the exact probabilities once again performs much better than the deterministic
one using similar hardware. With finite statistics, we see the same initially
strong performance as in the previous example followed by a leveling off of
improvement, with the deterministic ESN outperforming it on average by 6
detectors in the Lorenz $X$ task and by 5 detectors in the Sine-Square wave
task.

## Discussion

In this work, we define stochastic RCs as RCs that produce stochastic
outcomes, where the probability of each outcome rather than the individual
instantiation of the outcome serves as the readout. We provided necessary and
sufficient conditions for the contraction of the conditional probabilities,
which directly lead to critical elements for class universality, specifically,
uniform convergence and fading memory in the corresponding stochastic RC. We
rigorously proved that stochastic echo state networks whose activation
functions obey certain criteria form a universal approximating class. Then, we
considered two example stochastic RC hardware platforms theoretically: a qubit
reservoir network and a stochastic optical network. For both hardware
platforms, we investigated two practical numerical task examples, namely the
Sine-Square wave identification task and the Lorenz $X$ prediction task. We
found in those examples that while the stochastic ESN using the exact
probabilities performed much better than a deterministic ESN using similar
hardware, having to get an estimate for the probabilities using a finite
number of runs limited the performance of these stochastic ESNs.

From the plots of the average error measures on both the Lorenz $X$ and the
Sine-Square wave tasks in Figs. 2 and 4, we directly compared the performance
of both of our examples. While the stochastic ESNs demonstrate superior
performance over the deterministic ESNs with the same number of detectors, the
finite number of trials for obtaining the probability distributions limited
the performance of the stochastic ESNs. Moreover, we saw that generally the
qubit reservoir network stochastic ESN was able to achieve equal to slightly
better fits in both tasks in the stochastic case than the optical network
stochastic ESN. On the other hand, the deterministic ESN seems to perform
better at the Sine-Square wave tasks using the optical network activation
function as opposed to the qubit reservoir activation function. Interestingly,
this causes the stochastic optical network to show a bit of a better
comparison to the deterministic case on the Lorenz $X$ task, while the qubit
reservoir network has a more favorable comparison on the Sine-Square wave
identification task. Since nearly everything was the same in this comparison
between the qubit network ESN and the optical network ESN except for the
activation function used, it suggests that different choices of activation
function lead to changes in performance that are task dependent.

In addition to the direct use of stochastic RCs to perform specific
computational tasks, these machines may also prove useful in exploring a
conceptual question in quantum machine learning. Both the exponential scaling
of the Hilbert space and entanglement have been claimed as potential sources
of quantum advantage [29, 30, 31, 32], and yet it is still unclear to what
extent these properties can provide any form of computational speedup compared
to a classical counterpart. Stochastic RCs as described here exhibit an
exponential scaling of the state space, but do not in general employ any
quantum entanglement, as can be seen with the examples in the Results section.
Comparing the performance of these stochastic RCs with a genuine quantum RC
that incorporates entanglement into its dynamics could clarify the advantages
of quantum entanglement. By controlling for the scaling of the Hilbert space,
this comparison would provide a clearer understanding of which aspects of
quantum mechanics are beneficial for computational purposes.

Looking more closely at Figs. 2 and 4, it appears that the error measures in
both cases for the stochastic ESNs using the exact probabilities also level
off after a certain point, much like the 10000 run ESNs as discussed in
Supplementary Information A.2. In fact, if we keep increasing the number of
detectors of the deterministic ESNs beyond the values shown in the figures, we
would see the same leveling off of performance as well. Considering that the
smallest eigenvalue of $K$ continues getting smaller and smaller as the number
of detectors rises, it seems as though machine precision errors may
substantially limit the performance of larger RCs. One way to potentially
avoid this would be to use an alternative to optimizing via matrix inversion,
such as gradient descent. However, gradient descent methods have their own
issues resolving small values [33, 34, 35] which would likely be present in
this context. Another option would be to find reservoirs that do not lead to
small eigenvalues in $K$. These RCs would be able to achieve a manageable
signal-to-noise ratio for all of its eigentasks (see Ref. [36]), but this
would require a much stronger theoretical handle on how reservoir dynamics
lead to strong fits, which is at present a topic in need of extensive
research.

## Methods

Our proof of universality for stochastic RCs will follow the strategy used in
Refs. [13, 16], proving that the class of stochastic RCs satisfies the
criterion of the Stone-Weierstrass Theorem [37, Theorem 7.3.1]. In the next
few sections, we will explicitly prove the universality of certain classes of
stochastic ESNs, which implicitly proves the universality of the class
containing all stochastic RCs with the contraction property because ESNs are a
subclass of RCs.

To use the Stone-Weierstrass theorem, we first need to provide conditions on
the controlled transition matrix $p(u)$ that guarantees that the stochastic RC
will have the uniform convergence and fading memory properties, which we will
describe in further detail in the next section. Then, we must establish that
our class of stochastic ESNs obeys a polynomial algebra on a compact metric
space and show that it separates points in that algebra. These properties are
what the next two sections will demonstrate. Finally, we can then provide the
proof of Theorem 1 that utilizes all of these properties in conjunction with
the Stone-Weierstrass theorem, and conclude with some useful corollaries.

### Convergence and Fading Memory

In this section, we will prove that stochastic reservoir computers are
uniformly convergent and have the fading memory property at the level of the
probability distributions. The strategy for proving this will utilize Theorem
3.1 of [14] at its core. This theorem requires us to define a uniformly
bounded family of sequences as well as a contracting reservoir map. A
uniformly bounded family of sequences $K_{R}$ is defined to be the subset of
sequences in $(\mathbb{R}^{n})^{\mathbb{Z}_{-}}$ that obey

$\displaystyle
K_{R}=\left\\{(u_{k})_{k\in\mathbb{Z}_{-}}\in(\mathbb{R}^{n})^{\mathbb{Z}_{-}}\right|\left.\leavevmode\nobreak\
||u_{k}||\leq R\leavevmode\nobreak\ \forall k\in\mathbb{Z}_{-}\right\\},$ (21)

where $R$ is a positive number. The vector norm $||\cdot||$ can be chosen
arbitrarily, but for what follows it will be convenient to choose the 1-norm
$||x||_{1}=\sum_{i}|x_{i}|$. With $B_{n}[R]$ defined to be the closed ball of
length $R$ with respect to the 1-norm in $\mathbb{R}^{n}$ centered on the
origin, the reservoir map $f:B_{N}[R_{x}]\times B_{n}[R_{u}]\rightarrow
B_{N}[R_{x}]$ is said to be contracting with respect to the 1-norm if for all
$v_{1},v_{2}\in B_{N}[R_{x}],\leavevmode\nobreak\ u\in B_{n}[R_{u}]$ we have

$\displaystyle||f(v_{1},u)-f(v_{2},u)||_{1}\leq\epsilon||v_{1}-v_{2}||_{1}$
(22)

for some $0<\epsilon<1$.

Now we are ready to state a version of Theorem 3.1 from [14] that will be of
later use.

######  Theorem 2 (Uniform Convergence and Fading Memory).

With uniformly bounded sequences
$K_{R_{x}}\subset(\mathbb{R}^{N})^{\mathbb{Z}_{-}}$ and
$K_{R_{u}}\subset(\mathbb{R}^{n})^{\mathbb{Z}_{-}}$, if a continuous reservoir
map $f:B_{N}[R_{x}]\times B_{n}[R_{u}]\rightarrow B_{N}[R_{x}]$ is contracting
w.r.t the 1-norm, then the reservoir system associated with $f$ has the
uniform convergence property, meaning that for any given input sequence
$(u_{k})_{k\in\mathbb{Z}_{-}}\in K_{R_{u}}$ there exists a unique solution
$x\in B_{N}[R_{x}]$ to the equation
$x_{k+1}=f(x_{k},u_{k}),\leavevmode\nobreak\ \forall k\in\mathbb{Z}$.
Furthermore, because $f$ is a continuous reservoir map, there exists a unique
filter $U_{f}:K_{R_{u}}\rightarrow K_{R_{x}}$ associated with $f$ that is
causal, time-invariant, and has the fading memory property, meaning that the
filter $U_{f}$ is a continuous map between the metric spaces
$(K_{R_{u}},||\cdot||_{w})$ and $(K_{R_{x}},||\cdot||_{w})$ using the weighted
norm $||u||_{w}=\sup_{k\in\mathbb{Z}_{-}}\\{||u_{k}w_{-k}||_{1}\\}$, where
$w:\mathbb{N}\rightarrow(0,1]$ is some decreasing sequence that asymptotically
approaches zero.

To use Theorem 2, we must first establish that we have a continuous reservoir
map between bounded spaces. With respect to the 1-norm, it is easy to see that
every probability distribution $P_{k}$ belongs to $B_{M}[1]$ since all $M$
elements of $P_{k}$ are non-negative and must sum up to 1 for all
$k\in\mathbb{Z}$. Eq. (8) maintains the 1-norm of the probabilities since the
conditional probabilities obey $\sum_{a=1}^{M}p_{ab}(u)=1$ for any $b$, so
multiplication by the controlled transition matrix $p(u)$ does indeed qualify
as a reservoir map $f:B_{M}[1]\times B_{n}[R_{u}]\rightarrow B_{M}[1]$. In
fact, this maps vectors on the unit sphere $S_{M}[1]\in B_{M}[1]$ onto
$S_{M}[1]$, with the unit sphere defined by
$S_{M}(1)=\left\\{x\in\mathbb{R}^{M}\left|||x||=\sum_{i=1}^{M}|x_{i}|=1\right.\right\\}$,
so we can more specifically define $f:S_{M}[1]\times B_{n}[R_{u}]\rightarrow
S_{M}[1]$ Assuming that $p(u)$ is continuous in $u\in B_{n}[R_{u}]$ for some
radius $R_{u}$, the linearity of Eq. (8) w.r.t $P_{k}$ then guarantees that
this reservoir map is continuous. We have now demonstrated that, given the
assumptions about the continuity of $p(u)$ and that all $u_{k}$’s are bounded,
the probabilities in a stochastic RC satisfy all criteria for Theorem 2 aside
from contraction.

Now we are ready to establish the contraction criteria for stochastic
reservoir computing. The following theorem is essentially a statement about
the strict contraction of finite-state Markov chains with respect to the
1-norm, which is a basic and fundamental result that is taught in courses on
Markov chains [38, 39]. However, lectures notes like Refs. [38, 39] generally
only provide sufficient conditions for contraction, while the conditions
provided here are slightly stronger and are also necessary for convergence in
the context of finite-state Markov chains w.r.t the 1-norm.

######  Theorem 3.

A stochastic reservoir map $f:S_{M}[1]\times B_{n}[R_{u}]\rightarrow S_{M}[1]$
defined by a controlled transition matrix $p(u)$ that is continuous as a
function of $u$ is contracting w.r.t. the 1-norm if and only if the matrix
$p^{\top}(u)p(u)$ has no zero elements for all $u\in B_{n}[R_{u}]$.

Since convergence criterion of the conditional probability matrix is a well-
known result, the contents of Thm. 3 are almost surely in previous literature,
and we do not provide a proof in this section. For completeness, we do provide
a proof in Supplementary Information A.4.

With this theorem, we may now establish the full criteria for the uniform
convergence and fading memory properties of the probability distributions in a
stochastic RC.

######  Lemma 1.

With the uniformly bounded sequence
$K_{R_{u}}\subset(\mathbb{R}^{n})^{\mathbb{Z}_{-}}$, if a stochastic reservoir
map $f:S_{M}(1)\times B_{n}[R_{u}]\rightarrow S_{M}(1)$ (using the 1-norm)
defined by the controlled transition matrix $p(u)$ is continuous in $u$ and
$p^{\top}(u)p(u)$ has strictly positive elements for all $u\in B_{n}[R_{u}]$,
then the reservoir system associated with $f$ has the uniform convergence and
fading memory properties with respect to the probabilities for obtaining each
outcome.

This is simply proven by combining Theorems 2 and 3.

### Polynomial Algebra

Our next goal is to show that the class of stochastic ESNs that obey the
uniform convergence property and the fading memory property with respect to a
null sequence $w$ defined as part of Thm. 2 forms a polynomial algebra. More
specifically, writing $\mathcal{X}=\\{x^{(1)},\dots,x^{(M)}\\}$, for a
stochastic reservoir $F:\mathcal{X}\times B_{n}[R_{u}]\rightarrow\mathcal{X}$
with an associated controlled transition matrix $p(u)$ where multiplication by
$p$ generates a map $p:S_{M}[1]\times B_{n}[R_{u}]\rightarrow S_{M}[1]$, we
define a functional $G_{p,W}:K_{R_{u}}\rightarrow\mathbb{R}$ associated with
$p(u)$ and the weight vector $W$. The recursive application of the probability
map $p$ using an input sequence $(u_{k})_{k\in\mathbb{Z}_{-}}\in K_{R_{u}}$
yields a unique sequence of probability distributions
$(P_{k})_{k\in\mathbb{Z}_{-}}\in(S_{M}[1])^{\mathbb{Z}_{-}}$. Multiplication
of one of these probabilities by the weight vector $W$ generates a map
$W:S_{M}[1]\rightarrow\mathbb{R}$. Thus the functional $G_{p,W}$ for a
stochastic RC is the map from a given input sequence
$(u_{k})_{k\in\mathbb{Z}_{-}}\in K_{R_{u}}$ to some output
$\hat{y}_{0}\in\mathbb{R}$, which is the inner product of $W$ with the
probability distribution generated by these inputs and the controlled
transition matrix $p(u)$ at some arbitrary time, taken to be $k=0$. In the
event that the target sequence for our RC is vector-valued, we can view the
output as a vector whose elements are functionals $G_{p,W_{i}}$, each using
the same controlled transition matrix $p(u)$ but different weight vectors
$W_{i}$. Thus without loss of generality we can focus on scalar outputs
knowing that the results will straightforwardly generalize to vector outputs
as well.

With the set $\mathcal{G}_{w}$ of these functionals, which are all continuous
under the weighted norm
$||u||_{w}=\sup_{k\in\mathbb{Z}_{-}}\\{||u_{k}w_{-k}||_{1}\\}$ for a null
sequence $w$, we say that $\mathcal{G}_{w}$ forms a polynomial algebra if for
any
$G_{p_{1},W_{1}}\left((u_{k})_{k\in\mathbb{Z}_{-}}\right),G_{p_{2},W_{2}}\left((u_{k})_{k\in\mathbb{Z}_{-}}\right)\in\mathcal{G}_{w}$
we have that $G_{p_{1},W_{1}}\left((u_{k})_{k\in\mathbb{Z}_{-}}\right)+\lambda
G_{p_{2},W_{2}}\left((u_{k})_{k\in\mathbb{Z}_{-}}\right)\in\mathcal{G}_{w}$
for any scalar $\lambda$ and
$G_{p_{1},W_{1}}\left((u_{k})_{k\in\mathbb{Z}_{-}}\right)*G_{p_{2},W_{2}}\left((u_{k})_{k\in\mathbb{Z}_{-}}\right)\in\mathcal{G}_{w}$
using standard scalar addition and multiplication. We specifically want to
focus on the set $\mathcal{G}_{w}^{\varrho}$ of functionals corresponding to
stochastic ESNs with a specific controlled probability distribution $\varrho$,
in which any member
$G_{A,B,W}^{\varrho}\left((u_{k})_{k\in\mathbb{Z}_{-}}\right)$ of the set
corresponds to a controlled probability matrix $p(u)$ defined by $\varrho,A,$
and $B$ given in Eq. (13).

We begin by noting that if we take the matrix $A$ of any stochastic ESN of
size $N$ to be block diagonal with blocks of size $N_{1}$ and $N_{2}=N-N_{1}$,
and also separate $B$ so that

$\displaystyle A=\left(\begin{matrix}A_{1}&0\\\
0&A_{2}\end{matrix}\right),\quad\quad B=\left(\begin{matrix}B_{1}\\\
B_{2}\end{matrix}\right)$ (23)

The resulting outcomes $x^{(b)}$ will decompose analogously into independent
vectors $x_{1}^{(b_{1})}$ and $x_{2}^{(b_{2})}$, with
$x^{(b)\top}=\left(\begin{matrix}x_{1}^{(b_{1})\top}&x_{2}^{(b_{2})\top}\end{matrix}\right)$.
These components separately obey their own reservoir equations, so that the
controlled composite probability distribution and transition matrix factors
such that

$\displaystyle p_{s,a_{s}}(z)$
$\displaystyle=\prod_{i=1}^{N_{s}}\varrho_{a_{s,i}}(z_{i})$ (24)
$\displaystyle p_{ab}(u)$
$\displaystyle=p_{1,a}(A_{1}x_{1}^{(b_{1})}+B_{1}u)p_{2,a}(A_{2}x_{2}^{(b_{2})}+B_{2}u),$
(25)

with $s\in\\{1,2\\}$ and where $a_{s,i}$ is defined so that
$\xi^{(a_{s,i})}\equiv x_{s,i}^{(a_{s})}$, meaning that, when equipped with
weight vectors $W_{1}$ and $W_{2}$, these subsystems are themselves stochastic
ESNs corresponding to functionals
$G_{A_{1},B_{1},W_{1}}^{\varrho}\left((u_{k})_{k\in\mathbb{Z}_{-}}\right)$ and
$G_{A_{2},B_{2},W_{2}}^{\varrho}\left((u_{k})_{k\in\mathbb{Z}_{-}}\right)$
belonging to the class $\mathcal{G}_{w}^{\varrho}$. This means that they
generate their own probability distributions $P_{1,k}$ and $P_{2,k}$, and it
is not difficult to see that the distributions $P_{k}$ for the original ESN
using $A$ and $B$ are related by $P_{k,a}=P_{1,k,a_{1}}P_{2,k,a_{2}}$ for all
$a_{1}\in\\{1,\dots,m^{N_{1}}\\}$ and $a_{2}\in\\{1,\dots,m^{N_{2}}\\}$, with
the outcome $a$ corresponding to measuring outcomes $a_{1}$ and $a_{2}$
simultaneously.

Form here, it is easy to show that $\mathcal{G}_{w}^{\varrho}$ forms a
polynomial algebra. Taking the product of the two functionals defined above
yields
$G_{A_{1},B_{1},W_{1}}^{\varrho}\left((u_{k})_{k\in\mathbb{Z}_{-}}\right)*G_{A_{2},B_{2},W_{2}}^{\varrho}\left((u_{k})_{k\in\mathbb{Z}_{-}}\right)=(W_{1}^{\top}P_{1,0})(W_{2}^{\top}P_{2,0})=(W_{1}\otimes
W_{2})^{\top}(P_{1,0}\otimes P_{2,0})$. However, $P_{1,0}\otimes P_{2,0}$ is
just $P_{0}$ for the larger ESN that obeys Eq. (23), so with $W=W_{1}\otimes
W_{2}$ it is clear that the product belongs to $\mathcal{G}_{w}^{\varrho}$.
For the additive property, we must consider
$G_{A_{1},B_{1},W_{1}}^{\varrho}\left((u_{k})_{k\in\mathbb{Z}_{-}}\right)+\lambda
G_{A_{2},B_{2},W_{2}}^{\varrho}\left((u_{k})_{k\in\mathbb{Z}_{-}}\right)=W_{1}^{\top}P_{1,0}+\lambda
W_{2}^{\top}P_{2,0}$. Since $P_{1,0}$ and $P_{2,0}$ are probability
distributions, we can define the vector $e_{N}$ of length $N$ such that
$e_{N,i}=1$ for all $i\in\\{1,\dots,N\\}$ so that
$e_{N_{1}}^{\top}P_{1,0}=e_{N_{2}}^{\top}P_{2,0}=1$. We can then expand
$W_{1}^{\top}P_{1,0}=(W_{1}\otimes e_{N_{2}})^{\top}(P_{1,0}\otimes P_{2,0})$
and $W_{2}^{\top}P_{2,0}=(e_{N_{1}}\otimes W_{2})^{\top}(P_{1,0}\otimes
P_{2,0})$ to get $W_{1}^{\top}P_{1,0}+\lambda
W_{2}^{\top}P_{2,0}=(W_{1}\otimes e_{N_{2}}+\lambda e_{N_{1}}\otimes
W_{2})^{\top}(P_{1,0}\otimes P_{2,0})$. Therefore, using the the larger ESN
that obeys Eq. (23) again but with $W=W_{1}\otimes e_{N_{2}}+\lambda
e_{N_{1}}\otimes W_{2}$, we see that the sum of the functionals is also in
$\mathcal{G}_{w}^{\varrho}$. Hence $\mathcal{G}_{w}^{\varrho}$ forms a
polynomial algebra for any controlled probability distribution function
$\varrho$.

### Separation

The next and most involved step is to prove that the class
$\mathcal{G}_{w}^{\varrho}$ of stochastic ESNs defined by a scalar probability
distribution function $\varrho$ separates points in its domain $K_{R_{u}}$ of
bounded input sequences. We say that $\mathcal{G}_{w}^{\varrho}$ has the
separation property if for any
$(u_{k})_{k\in\mathbb{Z}_{-}},(v_{k})_{k\in\mathbb{Z}_{-}}\in K_{R_{u}}$ with
$(u_{k})_{k\in\mathbb{Z}_{-}}\neq(v_{k})_{k\in\mathbb{Z}_{-}}$ we can always
find at least one functional $G_{A,B,W}^{\varrho}\in\mathcal{G}_{w}^{\varrho}$
such that $G_{A,B,W}^{\varrho}\left((u_{k})_{k\in\mathbb{Z}_{-}}\right)\neq
G_{A,B,W}^{\varrho}\left((v_{k})_{k\in\mathbb{Z}_{-}}\right)$. To prove this
for a class $\mathcal{G}_{w}^{\varrho}$, we will look at the simplest case of
single-node ESNs. This property will not hold for arbitrary probability
distributions $\varrho$, so we will have to place conditions on $\varrho$ that
are necessary for satisfying this property.

######  Definition 1.

The vector space $\mathcal{V}_{R}^{p}\subseteq\mathbb{R}^{n}$ for some natural
number $n$ is called the invertible subspace with respect to the controlled
transition matrix $p(u)$ that is a function of $u\in[-R,R]$ if it is the
largest vector space that satisfies
$p(u)\mathcal{V}_{R}^{p}=\mathcal{V}_{R}^{p}$ for all $u\in[-R,R]$.

######  Theorem 4 (Separation).

Let $\mathcal{G}_{w}^{\varrho}$ be the class of stochastic ESNs defined by a
controlled probability distribution function $\varrho(\zeta)$ that is defined
on an interval $\zeta\in[-R_{\zeta},R_{\zeta}]$ for some $R_{\zeta}>0$. Also
let $\mathcal{V}_{R_{u}}^{p}\subseteq\mathbb{R}^{m}$ be the invertible
subspace with respect to the controlled transition matrix $p(u)$ whose
elements are defined by $p_{ab}(u)=\varrho_{a}(A\xi^{(b)}+Bu)$ for $m$
distinct measurement outcomes $\xi^{(b)}$ labeled by $b\in\\{1,\dots,m\\}$ and
for $u\in[-R_{u},R_{u}]$. The class $\mathcal{G}_{w}^{\varrho}$ has the
separation property if for some $R_{\mathrm{mono}}>0$ there exists a vector
$V\in\mathcal{V}_{R_{u}}^{p}$ for which $V^{\top}\varrho(\zeta)$ is strictly
monotonic on $\zeta$ in the interval $[-R_{\mathrm{mono}},R_{\mathrm{mono}}]$.

###### Proof.

To prove that the class $\mathcal{G}_{w}^{\varrho}$ of stochastic ESNs
separates points in $K_{R_{u}}$, we only need to show that for any pair of
input sequences $(u_{k})_{k\in\mathbb{Z}_{-}}$ and
$(v_{k})_{k\in\mathbb{Z}_{-}}\neq(u_{k})_{k\in\mathbb{Z}_{-}}$ we can find at
least one member
$G_{A,B,W}^{\varrho}\left((u_{k})_{k\in\mathbb{Z}_{-}}\right)$ of the class
$\mathcal{G}_{w}^{\varrho}$ for which
$G_{A,B,W}^{\varrho}\left((u_{k})_{k\in\mathbb{Z}_{-}}\right)\neq
G_{A,B,W}^{\varrho}\left((v_{k})_{k\in\mathbb{Z}_{-}}\right)$. In what
follows, it will turn out that there is always a single-node functional that
separates any two points, so we will restrict our analysis to the case where
$A$ and $B$ are both scalars. Consider a single-node stochastic ESN that
generates a functional in $\mathcal{G}_{w}^{\varrho}$. The general probability
of obtaining a given outcome $\xi^{(a)}$ is given by $P_{k,a}$. For a
stochastic ESN, The evolution of these probabilities are given by

$\displaystyle P_{k+1,a}=\sum_{b=1}^{m}\varrho_{a}(A\xi^{(b)}+Bu_{k})P_{k,b},$
(26)

where $u_{k},A,$ and $B$ are all scalars and there are $m$ possible outcomes.
For vector inputs, $B$ would also become a vector, and the proof that follows
can be applied by taking all but one of the components of $B$ to be zero. We
write $P_{k}$ to refer to the vector whose components are the probabilities
$P_{k,a}$. The functional is defined by $G_{A,B,W}^{\varrho}(u)=W^{\top}P_{0}$
for some $m$-dimensional vector $W$.

We are given two sequences
$(u_{k})_{k\in\mathbb{Z}_{-}},(v_{k})_{k\in\mathbb{Z}_{-}}\in K_{R_{u}}$ for
which $(u_{k})_{k\in\mathbb{Z}_{-}}\neq(v_{k})_{k\in\mathbb{Z}_{-}}$, meaning
that there is at least one $k\in\mathbb{Z}_{-}$ for which $u_{k}\neq v_{k}$.
Let $k_{\delta}$ be the largest number in $\mathbb{Z}_{-}$ for which
$u_{k_{\delta}}\neq v_{k_{\delta}}$, so that $u_{k}=v_{k}$ for all
$k\in\\{-1,-2,\dots,k_{\delta}+1\\}$. Without loss of generality, we can take
$u$ to be the sequence with the greater value of $u_{k_{\delta}}$, so that
$u_{k_{\delta}}-v_{k_{\delta}}>0$. With
$p_{ab}(u_{k})=\varrho_{a}(A\xi^{(b)}+Bu_{k})$, we see that

$\displaystyle G_{A,B,W}^{\varrho}\left((u_{k})_{k\in\mathbb{Z}_{-}}\right)$
$\displaystyle=W^{\top}\left(\prod_{\kappa=1}^{-k_{\delta}-1}p(u_{-\kappa})\right)p(u_{k_{\delta}})P_{k_{\delta}}\left((u_{k})_{k\in\mathbb{Z}_{-}}\right)$
(27)
$\displaystyle=\widetilde{W}^{\top}p(u_{k_{\delta}})P_{k_{\delta}}\left((u_{k})_{k\in\mathbb{Z}_{-}}\right),$
(28)

where
$\widetilde{W}^{\top}=W^{\top}\left(\prod_{\kappa=1}^{-k_{\delta}-1}p(u_{-\kappa})\right)$.
In this expression we write the dependence of $P_{k_{\delta}}$ on the the
input sequence $(u_{k})_{k\in\mathbb{Z}_{-}}$ explicitly as
$P_{k_{\delta}}\left((u_{k})_{k\in\mathbb{Z}_{-}}\right)$, which is defined by

$\displaystyle
P_{k_{\delta}}\left((u_{k})_{k\in\mathbb{Z}_{-}}\right)=\mathop{\lim}_{l\rightarrow\infty}\left(\prod_{k=1}^{l}p(u_{k_{\delta}-k})\right)P_{k_{\delta}-l}\left((u_{k})_{k\in\mathbb{Z}_{-}}\right),$
(29)

where the limit is guaranteed to converge because of the contraction property
of $p(u_{k})$, so the dependence of
$P_{k_{\delta}}\left((u_{k})_{k\in\mathbb{Z}_{-}}\right)$ on
$P_{k_{\delta}-l}\left((u_{k})_{k\in\mathbb{Z}_{-}}\right)$ vanishes as
$l\rightarrow\infty$. Because of how we defined $k_{\delta}$, we also have
that
$G_{A,B,W}^{\varrho}\left((v_{k})_{k\in\mathbb{Z}_{-}}\right)=\widetilde{W}^{\top}p(v_{k_{\delta}})P_{k_{\delta}}\left((v_{k})_{k\in\mathbb{Z}_{-}}\right)$
using the same vector $\widetilde{W}$. Thus by defining

$\displaystyle\widetilde{G}_{A,B,W}^{\varrho}(u_{k},P)\equiv
W^{\top}p(u_{k})P$ (30)

as a function of the vector $W$, number $u_{k}$, and probability distribution
$P$, we see that the filters can be written as

$\displaystyle
G_{A,B,W}^{\varrho}\left((u_{k})_{k\in\mathbb{Z}_{-}}\right)=\widetilde{G}_{A,B,\widetilde{W}}^{\varrho}\left(u_{k_{\delta}},P_{k_{\delta}}\left((u_{k})_{k\in\mathbb{Z}_{-}}\right)\right),$
(31)

with a similar expression for $v$.

In general the controlled transition matrix $p$ may not be invertible, and so
the resulting weight vector $\widetilde{W}$ may not be completely arbitrary.
However, if we assume that $\widetilde{W}\in\mathcal{V}_{R_{u}}^{p}$ as
defined in Thm. 4, then we can always find another vector
$V\in\mathcal{V}_{R_{u}}^{p}$ such that
$V^{\top}p(u_{k})=\widetilde{W}^{\top}$ for any $u_{k}\in[-R_{u},R_{u}]$.
Since $V\in\mathcal{V}_{R_{u}}^{p}$ as well, we can repeat this process as
many times as necessary to find the vector $W\in\mathcal{V}_{R_{u}}^{p}$ that
satisfies
$W^{\top}\left(\prod_{\kappa=1}^{-k_{\delta}-1}p(u_{-\kappa})\right)=\widetilde{W}^{\top}$
for any product $\prod_{\kappa=1}^{-k_{\delta}-1}p(u_{-\kappa})$. Thus for any
sequence $u$ we can always find a vector $W$ that will generate an arbitrary
vector $\widetilde{W}\in\mathcal{V}_{R_{u}}^{p}$.

The dependence on the rest of the sequence through
$P_{k_{\delta}}\left((u_{k})_{k\in\mathbb{Z}_{-}}\right)$ in the expression
for $G_{A,B,W}^{\varrho}\left((u_{k})_{k\in\mathbb{Z}_{-}}\right)$ is very
difficult to analyze, but we can take advantage of the fact that it is a
probability distribution to eliminate the dependence on the rest of the
sequence. Because of Eq. (31), it will suffice to show that
$G_{A,B,W}^{\varrho}\left((u_{k})_{k\in\mathbb{Z}_{-}}\right)\neq
G_{A,B,W}^{\varrho}\left((v_{k})_{k\in\mathbb{Z}_{-}}\right)$ by demonstrating
that there is always at least one $\widetilde{W},A$, and $B$ such that
$\widetilde{G}_{A,B,\widetilde{W}}^{\varrho}(u_{k_{\delta}},P_{1})\neq\widetilde{G}_{A,B,\widetilde{W}}^{\varrho}(v_{k_{\delta}},P_{2})$
for all distributions $P_{1}$ and $P_{2}$. Taking
$R(\zeta)=\widetilde{W}^{\top}\varrho(\zeta)$ for any scalar
$\zeta\in[-R_{\zeta},R_{\zeta}]$, we have

$\displaystyle\widetilde{G}_{A,B,\widetilde{W}}^{\varrho}(u_{k},P)=\widetilde{W}^{\top}p(u_{k})P=\sum_{b=1}^{m}R(A\xi^{(b)}+Bu_{k})P_{b}.$
(32)

Since $R(\zeta)$ is a scalar function, taking an average over a set of values
$\mathcal{Z}=\\{\zeta_{1},\dots,\zeta_{m}\\}$ for the function will give a
value between the minimum and maximum values of $R(\zeta)$ for
$\zeta\in\mathcal{Z}$. Thus we have that

$\displaystyle\widetilde{G}_{A,B,\widetilde{W}}^{\varrho}(u_{k},P)\in[\min_{\xi}R(A\xi+Bu_{k}),\max_{\xi}R(A\xi+Bu_{k})],$
(33)

where $\min_{\xi}f(\xi)=\min\\{f(\xi^{(1)}),\dots,f(\xi^{(m)})\\}$ and
likewise for $\max_{\xi}f(\xi)$. Thus we need to show that the above intervals
corresponding to $u_{k_{\delta}}$ and $v_{k_{\delta}}$ do not overlap for some
choice of $A,B,$ and $\widetilde{W}$.

Let us assume that there is some $\widetilde{W}\in\mathcal{V}_{R_{u}}^{p}$
such that $R(\zeta)$ is strictly monotonic in an interval
$[-R_{\mathrm{mono}},R_{\mathrm{mono}}]$ for some $R_{\mathrm{mono}}>0$.
Choose $A>0$ and $B>0$ such that
$\zeta_{\mathrm{max}}=A|\xi|_{\mathrm{max}}+BR_{u}\leq R_{\mathrm{mono}}$,
where $|\xi|_{\mathrm{max}}=\max\\{|\xi^{(1)}|,\dots,|\xi^{(m)}|\\}$. Because
$R(\zeta)$ is strictly monotonic on this interval, it is also bijective, so
any interval of the form in Eq. (33) corresponds to an interval
$[A\xi_{\mathrm{min}}+Bu_{k},A\xi_{\mathrm{max}}+Bu_{k}]\subset[-R_{\mathrm{mono}},R_{\mathrm{mono}}]$.
Since we took $u_{k_{\delta}}>v_{k_{\delta}}$, we only need to show that

$\displaystyle
A\xi_{\mathrm{min}}+Bu_{k_{\delta}}>A\xi_{\mathrm{max}}+Bv_{k_{\delta}}\quad\implies\quad
u_{k_{\delta}}-v_{k_{\delta}}>\frac{A}{B}(\xi_{\mathrm{max}}-\xi_{\mathrm{min}}),$
(34)

which can always be achieved by choosing
$A<B\frac{u_{k_{\delta}}-v_{k_{\delta}}}{\xi_{\mathrm{max}}-\xi_{\mathrm{min}}}$.
Therefore, under the above monotonicity assumption, since we have only placed
finite upper bounds on the positive numbers $A$ and $B$, and because
$\widetilde{W}\in\mathcal{V}_{R_{u}}^{p}$, for any pair of sequences
$(u_{k})_{k\in\mathbb{Z}_{-}},(v_{k})_{k\in\mathbb{Z}_{-}}\in K_{R_{u}}$ we
can always find choices of $A,B,$ and $W$ that result in a functional
$G_{A,B,W}^{\varrho}\left((u_{k})_{k\in\mathbb{Z}_{-}}\right)$ that separate
$(u_{k})_{k\in\mathbb{Z}_{-}}$ and $(v_{k})_{k\in\mathbb{Z}_{-}}$. This means
that the class of stochastic ESNs $\mathcal{G}_{w}^{\varrho}$ contains members
that can distinguish between any two input sequences
$(u_{k})_{k\in\mathbb{Z}_{-}},(v_{k})_{k\in\mathbb{Z}_{-}}\in K_{R_{u}}$ and
therefore separates points in $K_{R_{u}}$. ∎

The case where there are $m=2$ outcomes only is the simplest case, which also
makes it the most well-studied case in terms of potential hardware and the
most practical case to physically implement. In particular, any qubit-based
quantum computing architecture can be used to realize this class of stochastic
ESNs. The conditions for separability also simplify in this case, so we find
it useful to state the $m=2$ version of Theorem 4 as a corollary.

######  Corollary 1.

The class of stochastic ESNs $\mathcal{G}_{w}^{\varrho}$ defined by a
controlled probability distribution
$\varrho(\zeta)=\left(\begin{matrix}\varrho_{1}(\zeta)&(1-\varrho_{1}(\zeta))\end{matrix}\right)^{\top}$
defined on the interval $\zeta\in[-R_{\zeta},R_{\zeta}]$ for some
$R_{\zeta}>0$ has the separation property if for some $R_{\mathrm{mono}}>0$ we
have that $\varrho_{1}(\zeta)$ is strictly monotonic on the interval
$[-R_{\mathrm{mono}},R_{\mathrm{mono}}]$.

###### Proof.

The controlled transition matrix $p(u)$ for this system has the form

$\displaystyle
p(\zeta_{1},\zeta_{2})=\left(\begin{matrix}\varrho_{1}(\zeta_{1})&\varrho_{1}(\zeta_{2})\\\
1-\varrho_{1}(\zeta_{1})&1-\varrho_{1}(\zeta_{2})\end{matrix}\right)$ (35)

The determinant of this matrix is given by

$\displaystyle\det(p(\zeta_{1},\zeta_{2}))=\varrho_{1}(\zeta_{1})(1-\varrho_{1}(\zeta_{2}))-\varrho_{1}(\zeta_{2})(1-\varrho_{1}(\zeta_{1}))=\varrho_{1}(\zeta_{1})-\varrho_{1}(\zeta_{2})$
(36)

With $\zeta_{1},\zeta_{2}\in[-R_{\mathrm{mono}},R_{\mathrm{mono}}]$ for the
value of $R_{\mathrm{mono}}>0$ for which $\varrho_{1}(\zeta)$ is assumed to be
strictly monotonic, the only way this determinant can be 0 is if
$\zeta_{1}=\zeta_{2}$. Thus with $\zeta_{1,2}=A\xi^{(1,2)}+Bu_{k}$ and $A>0$
for distinct measurement outcomes $\xi^{(1)}$ and $\xi^{(2)}$, there is no
value of $u$ for which $\zeta_{1}=\zeta_{2}$, so as long as we choose $A$ and
$B$ such that $\zeta_{1},\zeta_{2}\in[-R_{\mathrm{mono}},R_{\mathrm{mono}}]$
for all $u\in[-R_{u},R_{u}]$, we can guarantee that $p(u)$ is invertible for
all $u\in[-R_{u},R_{u}]$. This means that the invertible subspace
$\mathcal{V}_{R_{u}}^{p}$ is simply $\mathbb{R}^{2}$. So by choosing
$V=\left(\begin{matrix}1&0\end{matrix}\right)^{\top}\in\mathbb{R}^{2}$, we
have that $V^{\top}\varrho(\zeta)=\varrho_{1}(\zeta)$ is strictly monotonic on
the interval $[-R_{\mathrm{mono}},R_{\mathrm{mono}}]$, satisfying the
separability criterion established in Thm. 4. ∎

### Proof of Theorem 1

With polynomial algebra proven and sufficient criteria for separability given
for classes of stochastic ESNs, we are now finally ready to apply the Stone-
Weierstrass Theorem [37, Theorem 7.3.1] and prove the universality of these
classes.

###### Proof of Theorem 1.

The conditions that the distribution $\varrho(\zeta)$ is continuous in $\zeta$
and satisfies $\varrho^{\top}(\zeta_{1})\varrho(\zeta_{2})>0$ for all
$\zeta_{1},\zeta_{2}\in[-R_{\zeta},R_{\zeta}]$ ensures that the conditional
probabilities for the ESNs are contracting by Theorem 3, so that the ESNs have
the uniform convergence and fading memory properties by Theorem 2. To
elaborate, the controlled composite probabilities of a general stochastic ESN
with dimension $L$ are given in Eq. (13) and have the form
$p_{a}(z)=\prod_{i=1}^{L}\varrho_{a_{i}}(z_{i})$. The continuity of $p_{a}(z)$
follows directly from the continuity of $\varrho(\zeta)$, and the condition
for contraction is satisfied if $\sum_{a=1}^{m^{L}}p_{a}(z_{1})p_{a}(z_{2})>0$
for all $L$-dimensional vectors $z_{1},z_{2}\in[-R_{\zeta},R_{\zeta}]^{L}$.
The sum over $a$ can be decomposed into $L$ sums over the $a_{i}$’s, so we
have
$\sum_{a}p_{a}(z_{1})p_{a}(z_{2})=\prod_{i=1}^{L}\left(\sum_{a_{i}=1}^{m}\varrho_{a_{i}}(z_{1,i})\varrho_{a_{i}}(z_{2,i})\right)=\prod_{i=1}^{L}\left(\varrho^{\top}(z_{1,i})\varrho(z_{2,i})\right)$.
The positivity of this product immediately follows from the condition that
$\varrho^{\top}(\zeta_{1})\varrho(\zeta_{2})>0$ for all scalars
$\zeta_{1},\zeta_{2}\in[-R_{\zeta},R_{\zeta}]$, so every stochastic ESN that
uses the probability distribution $\varrho$ has uniform convergence and fading
memory. This is important because by Thm. 2 it guarantees that each ESN
corresponds to a unique, causal, and time-invariant filter $U_{\varrho}$ as
well as a unique functional
$G_{A,B,W}^{\varrho}\left((u_{k})_{k\in\mathbb{Z}_{-}}\right)\in\mathcal{G}_{w}^{\varrho}$
on the space of uniformly bounded sequences $K_{R_{u}}$. Also, we can ensure
that any $z_{k}=Ax^{(a)}+Bu_{k}$ is in the interval
$[-R_{\zeta},R_{\zeta}]^{L}$ for any outcome $x^{(a)}$ and any input
$u_{k}\in[-R_{u},R_{u}]$ with appropriate choices of the matrix $A$ and vector
$B$.

To apply the Stone-Weierstrass Theorem on the class of functionals
$\mathcal{G}_{w}^{\varrho}$, it must form a polynomial algebra, contain the
constant functionals, and separate points in $K_{R_{u}}$. We established in a
previous section that any $\mathcal{G}_{w}^{\varrho}$ forms a polynomial
algebra with no additional conditions. It is also easy to show that by taking
$W$ such that $W_{a}=C$ for all indices $a$, the functional
$G_{A,B,W}^{\varrho}\left((u_{k})_{k\in\mathbb{Z}_{-}}\right)=W^{\top}P_{0}=\sum_{a}CP_{0,a}=C$,
so any class $\mathcal{G}_{w}^{\varrho}$ does indeed contain the constant
functionals. Finally, by Theorem 4 the additional condition that for some
$R_{\mathrm{mono}}>0$ there exists a vector $V\in\mathcal{V}_{R_{u}}^{p}$ for
which $V^{\top}\varrho(\zeta)$ is strictly monotonic on the interval
$[-R_{\mathrm{mono}},R_{\mathrm{mono}}]$ guarantees that
$\mathcal{G}_{w}^{\varrho}$ separates points in $K_{R_{u}}$. These three
properties allow us to apply the Stone-Weierstrass theorem to show that the
class $\mathcal{G}_{w}^{\varrho}$ is dense on the space of functionals over
$(K_{R_{u}},||\cdot||_{w})$, and so one can always find a functional
$G_{A,B,W}^{\varrho}\left((u_{k})_{k\in\mathbb{Z}_{-}}\right)\in\mathcal{G}_{w}^{\varrho}$
that approximates some arbitrary functional
$G^{*}\left((u_{k})_{k\in\mathbb{Z}_{-}}\right)$ on $K_{R_{u}}$ to arbitrary
precision, making the class $\mathcal{G}_{w}^{\varrho}$ of stochastic ESNs a
universal approximating class. ∎

As noted in the previous section, it is particularly useful to consider the
case when $m=2$. Here we state Theorem 1 given this specific case as another
corollary, which results in some simplification of the conditions necessary
for universality.

######  Corollary 2.

With the set of uniformly bounded sequences
$K_{R_{u}}\subset(\mathbb{R}^{n})^{\mathbb{Z}_{-}}$ and a weighted metric
$||\cdot||_{w}$ defined in Thm. 2, let $\mathcal{G}_{w}^{\varrho}$ be the
class of functionals generated by stochastic ESNs defined by a controlled
probability distribution
$\varrho(\zeta)=\left(\begin{matrix}\varrho_{1}(\zeta)&(1-\varrho_{1}(\zeta))\end{matrix}\right)^{\top}$
that is defined on the interval $\zeta\in[-R_{\zeta},R_{\zeta}]$, where
$\varrho_{1}(\zeta)$ is continuous in $\zeta$ and either
$\varrho_{1}(\zeta)\neq 0$ or $\varrho_{1}(\zeta)\neq 1$ for all
$\zeta\in[-R_{\zeta},R_{\zeta}]$. If for some $R_{\mathrm{mono}}>0$ we have
that $\varrho_{1}(\zeta)$ is strictly monotonic on the interval
$[-R_{\mathrm{mono}},R_{\mathrm{mono}}]$, then the class
$\mathcal{G}_{w}^{\varrho}$ is dense in the space of functionals over the
compact metric space $(K_{R_{u}},||\cdot||_{w})$.

###### Proof.

This is proven using exactly the same reasoning as in the proof of Theorem 1,
but with the simpler condition for separation of points in Corollary 1 when
the probability distribution $\varrho$ has only two components. Additionally,
the contraction constraint that
$\varrho^{\top}(\zeta_{1})\varrho(\zeta_{2})>0$ for all
$\zeta_{1},\zeta_{2}\in[-R_{\zeta},R_{\zeta}]$ can be simplified by observing
that the only way that this product can be zero for 2-dimensional probability
distributions is if one of the vectors is
$\left(\begin{matrix}1&0\end{matrix}\right)^{\top}$ and the other is
$\left(\begin{matrix}0&1\end{matrix}\right)^{\top}$. Thus we can ensure that
$\varrho^{\top}(\zeta_{1})\varrho(\zeta_{2})>0$ for all
$\zeta_{1},\zeta_{2}\in[-R_{\zeta},R_{\zeta}]$ by requiring that
$\varrho_{1}(\zeta)$ is either never equal to $0$ or never equal to $1$ on
$[-R_{\zeta},R_{\zeta}]$. ∎

Finally, we conclude this section by stating the universality of the class of
all stochastic RCs as defined in Section Stochastic Reservoir Computers.

######  Corollary 3.

The class $\mathcal{G}_{w}$ of all stochastic RCs with a contracting
controlled transition matrix $p(u)$ is universal.

###### Proof.

Because every stochastic RC in $\mathcal{G}_{w}$ has contracting transition
matrix $p(u)$ it will have an associated functional in $\mathcal{G}_{w}$ that
is unique by Thm. 2, so the class is well-defined in terms of the individual
stochastic RCs. Any class $\mathcal{G}_{w}^{\varrho}$ of stochastic ESNs that
satisfies the conditions for universality in Theorem 1 is contained in
$\mathcal{G}_{w}$, so we can always find a member
$G_{A,B,W}^{\varrho}\left((u_{k})_{k\in\mathbb{Z}_{-}}\right)\in\mathcal{G}_{w}^{\varrho}\subset\mathcal{G}_{w}$
that will approximate any causal, time-invariant filter to an arbitrary error
tolerance. The fact that
$G_{A,B,W}^{\varrho}\left((u_{k})_{k\in\mathbb{Z}_{-}}\right)\in\mathcal{G}_{w}$
as well for all $G_{A,B,W}^{\varrho}\left((u_{k})_{k\in\mathbb{Z}_{-}}\right)$
shows that $\mathcal{G}_{w}$ is also a universal approximating class, so in
this way the universality of $\mathcal{G}_{w}$ is inherited from the
universality of $\mathcal{G}_{w}^{\varrho}$. ∎

##  Appendix A Supplementary Information

###  A.1 Universality of Proposed Designs

We can show that qubit reservoir networks are universal using Corollary 2. If
we choose $\delta=\frac{\pi}{4}$ and restrict ourselves to choices of matrices
$A$ and vectors $B$ of dimension $L$ for which $z_{k}=Ax_{k}+Bu_{k}$ are in an
interval $[-R_{\zeta},R_{\zeta}]^{L}$ at every time step $k$ where
$R_{\zeta}<\frac{\pi}{4}$, then we guarantee that
$0<\sin^{2}(z_{k,i}+\frac{\pi}{4})<1$ for any $k$ and $i$ since
$0<z_{k,i}+\frac{\pi}{4}<\frac{\pi}{2}$. Also, we note that
$\sin^{2}(\zeta+\frac{\pi}{4})$ is strictly increasing for all
$\zeta\in[-R_{\zeta},R_{\zeta}]$ if $R_{\zeta}<\frac{\pi}{4}$. Thus the
controlled probability $\varrho_{1}(\zeta)$ given in Eq. (19) for qubit
reservoir networks with $\delta=\frac{\pi}{4}$ meets the requirements for a
universal approximating class as long the linearly transformed outcomes
$z_{k}=Ax_{k}+Bu_{k}$ have components that are always in the interval
$[-R_{\zeta},R_{\zeta}]^{L}$ for some $R_{\zeta}<\frac{\pi}{4}$.

We can also show that our stochastic optical networks form a universal
approximating class. Comparing to the conditions for universality in Corollary
2, the controlled probability given in Eq. (20) satisfies that
$\varrho_{1}(\zeta)\neq 1$ for any finite $\zeta$. However, because it is an
even function of $\zeta$, it is only strictly monotonic on the intervals
$(-\infty,0)$ and $(0,\infty)$. Similar to the previous example, this means we
have to choose an offset $d>0$ and make sure that $A$ and $B$ are chosen so
that every component of $z_{k}$ is in some closed interval
$[-R_{\zeta},R_{\zeta}]$ where $R_{\zeta}<d$. Thus by choosing some $d>0$ and
$0<R_{\zeta}<d$ so that at every time step $k$ and vector component $i$ we
have $z_{k,i}\in[-R_{\zeta},R_{\zeta}]$, we obtain a universal approximating
class using this stochastic ESN design.

###  A.2 Noise Analysis

Figure 5: The left plot shows NMSE values for the Lorenz $X$ task, while the
right plot shows the error percentage for the Sine-Square wave identification
task, both as a function of the number of runs for the qubit reservoir
network. Each data point represents the relevant error measure of the
stochastic ESN with fixed $A$ and $B$ using 2 detectors.

Figure 6: The left plot shows NMSE values for the Lorenz $X$ task, while the
right plot shows the error percentage for the Sine-Square wave identification
task, both as a function of the number of runs for the stochastic optical
network. Each data point represents the relevant error measure of the
stochastic ESN with fixed $A$ and $B$ using 2 detectors.

For a realistic implementation of a stochastic RC, we will need to get an
estimate for the probabilities $P_{k,a}$ of obtaining an outcome labeled by
$a$ at time step $k$ by averaging over many runs of the RC. Our procedure for
doing this is as follows. Each run evolves the readouts $x_{k}$ according to
Eq. (1), without any averaging or error mitigation. After running the RC for
$N_{\mathrm{runs}}$ times, we then count the number of occurrences
$n_{k}^{(a)}$ of the outcome $x_{k}=x^{(a)}$ at time step $k$ across all of
the runs. For a stochastic reservoir with $M$ outcomes, we then estimate

$\displaystyle P_{k,a}\approx\frac{n_{k}^{(a)}+1}{N_{\mathrm{runs}}+M}.$ (37)

The additional factors of $+1$ and $+M$ come from a Bayesian analysis of the
problem assuming no prior information about the distribution, where the above
expression for $P_{k,a}$ represents the expected value of $P_{k,a}$ given
$N_{\mathrm{runs}}$ prior measurements of $x_{k}$ where $n_{k}^{(a)}$ of them
resulted in $x^{(a)}$. This result can be found in general texts on Bayesian
analysis such as Ref. [40]. The main benefit to using this estimator is that
the estimates will never be exactly equal to 0, making it possible to divide
by components of the probability vector without having to check for a null
probability.

In our numerical analysis, we observed that when estimating the values of
$P_{k,a}$ with this method, the performance of the stochastic RC seems to
level off after a certain number of nodes. We provide an explanation for this
phenomenon through an analysis of how the linear regression behaves under the
noise introduced by the stochastic reservoir. Obtaining the optimal solution
for the weight vector $W$ involves calculating the matrix
$K=\frac{1}{N}\sum_{k}P_{k}P_{k}^{\top}$ and its inverse $K^{-1}$, where the
sum over $k$ is taken over the $N$ training data steps. In general, $K$ has
eigenvalues that span a wide range of scales, from some values that are order
1 to some that are many orders of magnitude less. With a finite number of
shots, the matrix
$K_{\mathrm{est}}=\frac{1}{N}\sum_{k}P_{k,\mathrm{est}}P_{k,\mathrm{est}}^{\top}$
using the estimated probabilities $P_{k,\mathrm{est}}\approx P_{k}$ has an
expectation value given by

$\displaystyle\langle
K_{\mathrm{est}}\rangle=\left(1-\frac{1}{N_{\mathrm{runs}}}\right)K+\frac{1}{N_{\mathrm{runs}}}P_{k}^{\mathrm{diag}},$
(38)

where $(P_{k}^{\mathrm{diag}})_{ab}=\delta_{ab}P_{k,a}$ is the probability
vector $P_{k}$ broadcast into a diagonal matrix. This expression is closely
related to the eigentask analysis of the effects of noise in expressive
capacity introduced in Ref. [36], and for $N_{\mathrm{runs}}=1$ will lead to
some of the results of Ref. [41].

This performance hit is illustrated for a 2-detector example in Fig. 5, where
we plot the error measures for both tasks using the qubit reservoir network as
a function of $N_{\mathrm{runs}}$ for specific choices of $A$ and $B$. For the
qubit reservoir network on the Lorenz $X$ task, the smallest eigenvalue of $K$
is $\lambda_{\min}=3.73*10^{-5}$, while for the Sine-Square wave tasks the
smallest eigenvalue of $K$ is $\lambda_{\min}=5.24*10^{-4}$. We can then
expect heuristically that, with $N_{\mathrm{runs}}=100000$ for the Lorenz $X$
task and $N_{\mathrm{runs}}=10000$ in the Sine-Square wave task, we should be
able to resolve every eigenstate of $K$. This is exactly what we see in Fig.
5, where the error measures are only decreasing logarithmically with the
number of runs until about $N_{\mathrm{runs}}=8*10^{4}\approx
3.0/\lambda_{\min}$ for the Lorenz $X$ task and
$N_{\mathrm{runs}}=4*10^{3}\approx 2.1/\lambda_{\min}$, where they begin to
converge to the error values derived using the exact probabilities.

We can perform the same analysis on the stochastic optical network for the
2-detector examples shown in Fig. 6, where we plot the relevant error measures
as a function of $N_{\mathrm{runs}}$ for specific choices of $A$ and $B$ (the
same choices used in Fig. 5). For the Lorenz $X$ task, the smallest eigenvalue
of $K$ is $\lambda_{\min}=1.32*10^{-5}$, while for the Sine-Square wave tasks
the smallest eigenvalue of $K$ is $\lambda_{\min}=2.35*10^{-4}$. We can then
expect heuristically that, with $N_{\mathrm{runs}}>100000$ for the Lorenz $X$
task and $N_{\mathrm{runs}}=10000$ in the Sine-Square wave task, we should be
able to resolve every eigenstate of $K$. This is exactly what we see in Fig.
6, where the error measures decrease logarithmically with the number of runs
until roughly $N_{\mathrm{runs}}=1.5*10^{5}\approx 2.0/\lambda_{\min}$ for the
Lorenz $X$ task and $N_{\mathrm{runs}}=8*10^{3}\approx 1.9/\lambda_{\min}$,
where they begin to converge to the error values derived using the exact
probabilities. Note that for the Lorenz $X$ task the trend appears to have a
plateau from about $3*10^{2}$ to $2*10^{3}$ runs, which could have been taken
as a signal that the theoretical minimum NMSE for this ESN is about $0.24$ in
a realistic calculation where we have no prior knowledge of the true minimum.

We can compare the performance boost in these examples due to increased
statistics through Figs. 5 and 6. We see that in both cases the qubit
reservoir network reaches the infinite statistics values a bit earlier than
the stochastic optical network in both cases, which may explain the general
performance difference we see in the other figures. Otherwise, the two plots
are very similar, with the qubit reservoir perhaps showing a more steady
improvement as a function of $N_{\mathrm{runs}}$ than the more curved trends
for the optical network. Note that the values of $\lambda_{\min}$ we found for
the Lorenz $X$ task are more than an order of magnitude smaller than for the
Sine-Square wave task. This will influence where the stochastic ESN
performance stalls out in Figs. 2 and 4, and it would explain why we needed
$N_{\mathrm{runs}}=100000$ to get comparable results for the Lorenz $X$ task
while only $N_{\mathrm{runs}}=10000$ was needed for the Sine-Square wave task.

###  A.3 Simulation Details

####  A.3.1 Task Details

The Sine-Square wave identification task [27] requires the ESN to identify
whether a given section of the input is a sine wave or a square wave. The
target sequence for this task is a random string $\\{\tau(m)\\}$ of ones and
zeroes, with equal probability of choosing 1 or 0 for each $m$, where a zero
corresponds to a sine wave input while a one corresponds to a square wave
input. The period of both waves are tuned to be exactly 8 time steps, and at
the end of each period the input wave has a random chance to switch to the
other waveform. The formulas for theses waves in terms of $k$ are given by

$\displaystyle u_{k}$
$\displaystyle=\delta_{0,\tau\left(\mathrm{floor}\left(\frac{k}{8}\right)\right)}u_{k}^{(sin)}+\delta_{1,\tau\left(\mathrm{floor}\left(\frac{k}{8}\right)\right)}u_{k}^{(sqr)},$
(39) $\displaystyle u_{k}^{(sin)}$
$\displaystyle=\sin\left(\frac{\pi}{4}k\right),$ (40) $\displaystyle
u_{k}^{(sqr)}$ $\displaystyle=(-1)^{\mathrm{floor}\left(\frac{k}{4}\right)},$
(41)

where $\mathrm{floor}\left(s\right)$ is the largest integer $n_{s}$ such that
$n_{s}\leq s$ for any real number $s$. The task for the ESN is to distinguish
the sine and square wave parts of the input, yielding an output of
$\hat{y}_{k}<0.5$ for a sine wave and $\hat{y}_{k}>0.5$ for a square wave. In
this work, we generated one random sine-square wave input and used it for all
of the numerical work found in the following sections.

The other task we use in this work is the Lorenz $X$ task, which requires the
ESN to approximate a chaotic dynamical system described by $X(t)$ in the
Lorenz or equations:

$\displaystyle\frac{dX}{dt}(t)$ $\displaystyle=a(Y-X)$ (42)
$\displaystyle\frac{dY}{dt}(t)$ $\displaystyle=X(b-Z)-Y$ (43)
$\displaystyle\frac{dZ}{dt}(t)$ $\displaystyle=XY-cZ$ (44)

where we choose the standard values $a=10,b=28,$ and $c=8/3$. We numerically
approximate the solution to this equation using the fourth-order Runge-Kutta
method with $\delta t=1.0$ and initial conditions of $X(0)=0.5,Y(0)=0.1,$ and
$Z(0)=0.2$. We also run the solution for $1000$ steps before using it for the
task, or in other words the target sequence we use actually starts with
$y_{1000}$. The task for the ESN is to predict what the sequence will be 1
time step into the future. In other words, using the input sequence
$(u_{k})_{k\in\mathbb{Z}}=(y_{k-1})_{k\in\mathbb{Z}}$, we want the ESN to
successfully predict $(y_{k})_{k\in\mathbb{Z}}$.

####  A.3.2 Simulation Parameters

In the numerical analysis illustrated in Figs. 2 and 4, the left plot of both
figures shows a comparison of NMSE values for the Lorenz $X$ task averaged
over 100 different choices of $A$ matrices and $B$ vectors for the ESN, while
the right plot shows a comparison of the percentage of errors recorded for the
Sine-Square wave identification task averaged over 1000 choices. For the
deterministic ESNs, we use the exact same choices of $A$ and $B$, using the
probability function $\varrho_{1}(\zeta)$ as the deterministic activation
function, so that the reservoir equations are given by

$\displaystyle(x_{\mathrm{det}})_{k+1}$
$\displaystyle=\langle\Gamma(A(x_{\mathrm{det}})_{k}+Bu_{k})\rangle=\varrho_{1}(A(x_{\mathrm{det}})_{k}+Bu_{k})$
(45) $\displaystyle\hat{y}_{k}$ $\displaystyle=W(x_{\mathrm{det}})_{k}+C,$
(46)

where the scalar activation function $\Gamma(\zeta)$ and controlled
probability $\varrho_{1}(\zeta)$ are applied elementwise to their vector
arguments. In Fig. 2, we use $\varrho_{1}(\zeta)$ given in Eq. (19)
corresponding to the qubit reservoir network, while in Fig. 4 we use
$\varrho_{1}(\zeta)$ given in Eq. (20) for the stochastic optical network. For
the stochastic ESNs, we plot the cases where we average over a finite number
of runs as well as the theoretical limit of infinite runs using the exact
probabilities. We used 100000 runs in the finite case for the Lorenz $X$ task
and 10000 runs for the Sine-Square wave task.

For each run, we used 3000 training data steps performed after 100 steps of
initial startup time for the Lorenz $X$ task and 96 steps for the Sine-Square
wave task. The values for the error measures shown here are for a set of time
steps immediately after the training data set, totaling 500 steps for the
Lorenz $X$ task and 504 for the Sine-Square wave task. These same parameters
apply to the data in Figs. 5 and 6 from the previous section as well. However,
for the Sine-Square wave plots in these figures, each data point is an average
over 10 repetitions of the error percentage calculation, using the same number
of runs each time, to reduce the spread of the values in this plot.

###  A.4 Proof of Theorem 3

###### Proof.

Starting from the left-hand side of Eq. (22), we have for the stochastic
reservoir map defined by $p(u)$ acting on two probability vectors
$P_{1},P_{2}\in S_{M}[1]$ and an input $u\in B_{n}[R_{u}]$

$\displaystyle||f(P_{1},u)-f(P_{2},u)||_{1}=||p(u)(P_{1}-P_{2})||_{1}$ (47)

Using the triangle inequality, we can see that

$\displaystyle||p(u)(P_{1}-P_{2})||_{1}$
$\displaystyle=\sum_{a=1}^{M}\left|\sum_{b=1}^{M}p_{ab}(u)(P_{1}-P_{2})_{b}\right|$
(48)
$\displaystyle\leq\sum_{a=1}^{M}\sum_{b=1}^{M}\left|p_{ab}(u)(P_{1}-P_{2})_{b}\right|$
(49)
$\displaystyle=\sum_{b=1}^{M}\left(\sum_{a=1}^{M}p_{ab}(u)\right)\left|(P_{1}-P_{2})_{b}\right|$
(50)
$\displaystyle=\sum_{b=1}^{M}(1)\left|(P_{1}-P_{2})_{b}\right|=||P_{1}-P_{2}||_{1},$
(51)

where in the third line we used the non-negativity of probabilities.

For the absolute value of two scalars $c$ and $d$, the inequality
$|c+d|\leq|c|+|d|$ is saturated if and only if $c$ and $d$ are both non-
negative or both non-positive, or more succinctly if the product $cd$ is non-
negative. Thus for $n$ scalars the inequality
$|\sum_{i=1}^{n}c_{i}|\leq\sum_{i=1}^{n}|c_{i}|$ is saturated if and only if
all of the $c_{i}$’s are either non-negative or non-positive, or alternatively
if the products $c_{i_{1}}c_{i_{2}}$ are non-negative for all
$i_{1},i_{2}\in\\{1,\dots,n\\}$. This means that for our use of the triangle
inequality in Eq. (48), equality occurs if and only if for every
$a\in{1,\dots,M}$ the products
$p_{ab_{1}}(u)p_{ab_{2}}(u)(P_{1}-P_{2})_{b_{1}}(P_{1}-P_{2})_{b_{2}}\geq 0$
for all $b_{1},b_{2}\in\\{1,\dots,M\\}$. Therefore, the inequality is strict
if and only if for some $a\in\\{1,\dots,M\\}$ there is a pair of distinct
indices $b_{1},b_{2}\in\\{1,\dots,M\\}$ for which
$p_{ab_{1}}(u)p_{ab_{2}}(u)(P_{1}-P_{2})_{b_{1}}(P_{1}-P_{2})_{b_{2}}<0$. The
inequality in Eq. (48) will be strict if the triangle inequality is strict for
just one index $a$ because the sum over $a$ in Eq. (48) is over non-negative
quantities, so if any one of the terms in the sum fails to saturate the
triangle inequality, the entire sum will not saturate it, either.

For any two distinct probability distributions
$||P_{1}-P_{2}||_{1}=\sum_{a=1}^{M}|P_{1,a}-P_{2,a}|>0$ but
$\sum_{a=1}^{M}(P_{1,a}-P_{2,a})=1-1=0$. This implies that there is at least
one element of $P_{1}-P_{2}$ that is positive and at least one other that is
negative, so there will always be at least one pair of distinct indices
$b_{1},b_{2}$ for which $(P_{1}-P_{2})_{b_{1}}(P_{1}-P_{2})_{b_{2}}<0$.
However, if we consider the distributions $P_{1,a}=\delta_{ab_{1}}$ and
$P_{2,a}=\delta_{ab_{2}}$ with $b_{1}\neq b_{2}$, then we see that
$b_{1},b_{2}$ is the only pair of distinct indices for which
$(P_{1}-P_{2})_{b_{1}}(P_{1}-P_{2})_{b_{2}}<0$. Thus in order for
$p_{ab_{1}}(u)p_{ab_{2}}(u)(P_{1}-P_{2})_{b_{1}}(P_{1}-P_{2})_{b_{2}}<0$ to
hold for any distinct probability distributions $P_{1},P_{2}\in S_{M}(1)$, we
must have that $p_{ab_{1}}(u)p_{ab_{2}}(u)\neq 0$ for all pairs of distinct
indices $b_{1},b_{2}\in\\{1,\dots,M\\}$. Since $p_{ab}(u)$ is always non-
negative, these inequalities are all satisfied when $p_{ab}(u)>0$ for all
$b\in\\{1,\dots,M\\}$. Therefore, the inequality is strict for all distinct
distributions $P_{1},P_{2}\in S_{M}(1)$ if and only if for all pairs of
distinct indices $b_{1},b_{2}\in\\{1,\dots,M\\}$ we have for some
$a\in\\{1,\dots,M\\}$ that $p_{ab_{1}}(u)p_{ab_{2}}(u)>0$.

Finally, we note that because the controlled transition matrix $p(u)$ has non-
negative entries for all inputs, $\sum_{a=1}^{M}p_{ab_{1}}(u)p_{ab_{2}}(u)>0$
if and only if $p_{ab_{1}}(u)p_{ab_{2}}(u)>0$ for at least one
$a\in\\{1,\dots,M\\}$. But
$\sum_{a=1}^{M}p_{ab_{1}}(u)p_{ab_{2}}(u)=\left(p^{\top}(u)p(u)\right)_{b_{1}b_{2}}$,
so we can restate the condition for strict inequality using this matrix. Thus
we have that the triangle inequality is strict for all distinct probability
distributions $P_{1},P_{2}\in S_{M}(1)$ and all inputs $u\in B_{n}[R_{u}]$ if
and only if the matrix $p^{\top}(u)p(u)$ has strictly positive elements for
all $u\in B_{n}[R_{u}]$. Also, because each element of $p(u)$ is continuous
over all $u$ in the compact set $B_{n}[R_{u}]$, the image of $p(u)$ is also
compact. This means that $d=\inf_{u}p^{\top}(u)p(u)$ is contained in the image
of $p^{\top}(u)p(u)$ on $u\in B_{n}[R_{u}]$, and since $p^{\top}(u)p(u)$ has
strictly positive elements for all $u\in B_{n}[R_{u}]$ we have that $d>0$, and
so $p^{\top}(u)p(u)\geq d$ for some $d>0$. Therefore the contraction criterion

$\displaystyle||p(u)(P_{1}-P_{2})||_{1}\leq\epsilon||P_{1}-P_{2}||_{1}$ (52)

holds for a continuous probability matrix $p(u)$ and some $0<\epsilon<1$ if
and only if the matrix $p^{\top}(u)p(u)$ has strictly positive elements for
all inputs $u\in B_{n}[R_{u}]$ since
$||p(u)(P_{1}-P_{2})||_{1}<||P_{1}-P_{2}||_{1}$ for all distinct probability
distributions $P_{1},P_{2}\in S_{M}(1)$ and
$||p(u)(P_{1}-P_{2})||_{1}=||P_{1}-P_{2}||_{1}=0$ when $P_{1}=P_{2}$. ∎

In addition, we leave off with a remark that is intended to show that a fully
stochastic reservoir map is not difficult to achieve. This is equivalent to
the final, stronger convergence criterion given in Ref. [39].

######  Remark 1.

If for every $u\in B_{n}[R_{u}]$ there is at least one value of $a$ for which
$p_{ab}(u)\neq 0$ for all $b\in\\{1,\dots,M\\}$, then $p^{\top}(u)p(u)$ has
strictly positive elements for all inputs $u\in B_{n}[R_{u}]$.

Since
$\left(p^{\top}(u)p(u)\right)_{b_{1}b_{2}}=\sum_{a=1}^{M}p_{ab_{1}}(u)p_{ab_{2}}(u)$,
if there is at least one index $a$ for which $p_{ab}(u)\neq 0$ for all
$b\in\\{1,\dots,M\\}$, then $p_{ab_{1}}(u)p_{ab_{2}}(u)>0$ for all
$b_{1},b_{2}\in\\{1,\dots,M\\}$ as well. If this holds for every input $u$,
then this guarantees that $p^{\top}(u)p(u)$ always has strictly positive
elements for all $u\in B_{n}[R_{u}]$. Thus if we can guarantee that for every
$u$ there is at least one row of $p(u)$ for which every probability is
nonzero, then it will be contracting.

###  A.5 Note on a Recent Result

Stochastic RCs as we define them have been considered in a recent paper by
Polloreno [41]. In this work, it is shown that the information processing
capacity of a single run of a stochastic RC only scales polynomially with the
scale of the hardware. However, in our work we run the RC many times and get
estimates for the probabilities $P_{k,a}$ by counting the number of
occurrences of each outcome. The derivation of Theorem 2 in Ref. [41] is
predicated on the fact that the second moment of the outcome probabilities
$P_{k,a}$ is linear in these probabilities, but this is no longer true when
$P_{k,a}$ is estimated with an average over multiple shots. In Supplementary
Information A.2, we show that past a certain number of aggregated runs, the
estimated probabilities provide a very good estimate of the results using the
exact probabilities. Of course, it could still be argued that the number of
runs needed to achieve this level of accuracy may be exponentially scaling
with the scale of the hardware as well. In any case, we do not claim that
stochastic RCs have exponentially scaling computational complexity. We merely
state that the number of computational nodes of the reservoir is exponentially
scaling, and these nodes may or may not be used effectively for any given
computation.

## References

  *   * Antonik et al. [2018] Antonik, P., Gulina, M., Pauwels, J., Massar, S.: Using a reservoir computer to learn chaotic attractors, with applications to chaos synchronization and cryptography. Physical Review E 98(1), 012215 (2018)
  * Chen and Nurdin [2022] Chen, J., Nurdin, H.I.: Nonlinear autoregression with convergent dynamics on novel computational platforms. IEEE Transactions on Control Systems Technology 30(5), 2228–2234 (2022) https://doi.org/10.1109/TCST.2021.3136227
  * Lu et al. [2017] Lu, Z., et al.: Reservoir observers: Model-free inference of unmeasured variables in chaotic systems. Chaos 27, 041102 (2017)
  * Larger et al. [2017] Larger, L., et al.: High-speed photonic reservoir computing using a time-delay-based architecture: Million words per second classification. Physical Review X 7(1), 011015 (2017)
  * Canaday et al. [2018] Canaday, D., Griffith, A., Gauthier, D.J.: Rapid time series prediction with a hardware-based reservoir computer. Chaos: An Interdisciplinary Journal of Nonlinear Science 28(12), 123119 (2018)
  * Pathak et al. [2018] Pathak, J., et al.: Model-free prediction of large spatiotemporally chaotic systems from data: A reservoir computing approach. Physical review letters 120(2), 024102 (2018)
  * Rafayelyan et al. [2020] Rafayelyan, M., et al.: Large-scale optical reservoir computing for spatiotemporal chaotic systems prediction. Phys. Rev. X 10, 041037 (2020)
  * Antonik et al. [2020] Antonik, P., Marsal, N., Rontani, D.: Large-scale spatiotemporal photonic reservoir computer for image classification. IEEE Journal of Selected Topics in Quantum Electronics 26(1), 1–12 (2020) https://doi.org/10.1109/JSTQE.2019.2924138
  * Bianchi et al. [2021] Bianchi, F.M., Scardapane, S., Løkse, S., Jenssen, R.: Reservoir computing approaches for representation and classification of multivariate time series. IEEE Transactions on Neural Networks and Learning Systems 32(5), 2169–2179 (2021) https://doi.org/10.1109/TNNLS.2020.3001377
  * Carroll [2021] Carroll, T.L.: Optimizing reservoir computers for signal classification. Frontiers in Physiology 12, 685121 (2021)
  * Angelatos et al. [2017] Angelatos, G., Khan, S.A., Türeci, H.: Reservoir computing approach to quantum state measurement. Physical Review X 11, 041062 (2017)
  * Tanaka et al. [2019] Tanaka, G., Yamane, T., Héroux, J.B., Nakane, R., Kanazawa, N., Takeda, S., Numata, H., Nakano, D., Hirose, A.: Recent advances in physical reservoir computing: A review. Neural Networks 115, 100–123 (2019)
  * Grigoryeva and Ortega [2018a] Grigoryeva, L., Ortega, J.-P.: Universal discrete-time reservoir computers with stochastic inputs and linear readouts using non-homogeneous state-affine systems. Journal of Machine Learning Research 19(24), 1–40 (2018)
  * Grigoryeva and Ortega [2018b] Grigoryeva, L., Ortega, J.-P.: Echo state networks are universal. Neural Networks 108, 495–508 (2018) https://doi.org/10.1016/j.neunet.2018.08.025
  * Bishop [2022] Bishop, A.N.: Universal Time-Uniform Trajectory Approximation for Random Dynamical Systems with Recurrent Neural Networks. Preprint at arXiv:2211.08018 (2022)
  * Chen and Nurdin [2019] Chen, J., Nurdin, H.I.: Learning nonlinear input–output maps with dissipative quantum systems. Quantum Information Processing 18(7) (2019) https://doi.org/10.1007/s11128-019-2311-9
  * Gonon and Ortega [2020] Gonon, L., Ortega, J.-P.: Reservoir computing universality with stochastic inputs. IEEE Transactions on Neural Networks and Learning Systems 31(1), 100–112 (2020) https://doi.org/10.1109/TNNLS.2019.2899649
  * Hart et al. [2020] Hart, A., Hook, J., Dawes, J.: Embedding and approximation theorems for echo state networks. Neural Networks 128, 234–247 (2020) https://doi.org/10.1016/j.neunet.2020.05.013
  * Li et al. [2023] Li, B., Fong, R.S., Tiňo, P.: Simple Cycle Reservoirs are Universal. Preprint at arXiv:2308.10793 (2023)
  * Li and Yang [2023] Li, Z., Yang, Y.: Universality and approximation bounds for echo state networks with random weights. IEEE Transactions on Neural Networks and Learning Systems, 1–13 (2023) https://doi.org/10.1109/TNNLS.2023.3339512
  * Martínez-Peña and Ortega [2023] Martínez-Peña, R., Ortega, J.-P.: Quantum reservoir computing in finite dimensions. Phys. Rev. E 107, 035306 (2023) https://doi.org/10.1103/PhysRevE.107.035306
  * Nokkala et al. [2021] Nokkala, J., Martínez-Peña, R., Giorgi, G.L., Parigi, V., Soriano, M.C., Zambrini, R.: Gaussian states of continuous-variable quantum systems provide universal and versatile reservoir computing. Communications Physics 4(1) (2021) https://doi.org/10.1038/s42005-021-00556-w
  * Sugiura et al. [2024] Sugiura, S., Ariizumi, R., Asai, T., Azuma, S.-i.: Existence of reservoir with finite-dimensional output for universal reservoir computing. Scientific Reports 14, 8448 (2024) https://doi.org/10.1038/s41598-024-56742-7
  * Roques-Carmes et al. [2023] Roques-Carmes, C., et al.: Biasing the quantum vacuum to control macroscopic probability distributions. Science 381(6654), 205–209 (2023)
  * Ross [2023] Ross, S.M.: Introduction to Probability Models. 13e. Academic Press, Amsterdam, Boston (2023)
  * Meyn et al. [2009] Meyn, S., Tweedie, R.L., Glynn, P.W.: Markov Chains and Stochastic Stability. 2e. Cambridge University Press, Cambridge (2009)
  * Dudas et al. [2023] Dudas, J., Carles, B., Plouet, E., Mizrahi, F.A., Grollier, J., Marković, D.: Quantum reservoir computing implementation on coherently coupled quantum oscillators. npj Quantum Information 9(1), 64 (2023)
  * Ma et al. [2023] Ma, S.-Y., Wang, T., Laydevant, J., Wright, L.G., McMahon, P.L.: Quantum-noise-limited optical neural networks operating at a few quanta per activation. Preprint at arXiv:2307.15712 (2023)
  * Götting et al. [2023] Götting, N., Lohof, F., Gies, C.: Exploring quantumness in quantum reservoir computing. Phys. Rev. A 108, 052427 (2023) https://doi.org/10.1103/PhysRevA.108.052427
  * Kenigsberg et al. [2006] Kenigsberg, D., Mor, T., Ratsaby, G.: Quantum advantage without entanglement. Quantum Inf. Comput. 6(7), 606–615 (2006)
  * Kora et al. [2024] Kora, Y., Zadeh-Haghighi, H., Stewart, T.C., Heshami, K., Simon, C.: Frequency-dependent entanglement advantage in spin-network Quantum Reservoir Computing. Preprint at arXiv:2403.08998 (2024)
  * Mujal et al. [2023] Mujal, P., Martínez-Peña, R., Giorgi, G.L., Soriano, M.C., Zambrini, R.: Time-series quantum reservoir computing with weak and projective measurements. npj Quantum Information 9(1), 16 (2023)
  * McClean et al. [2018] McClean, J.R., Boixo, S., Smelyanskiy, V.N., Babbush, R., Neven, H.: Barren plateaus in quantum neural network training landscapes. Nature communications 9(1), 4812 (2018)
  * Kolen and Kremer [2001] Kolen, J.F., Kremer, S.C.: Gradient Flow in Recurrent Nets: The Difficulty of Learning Long Term Dependencies, pp. 237–243. Wiley-IEEE Press, New York (2001). https://doi.org/10.1109/9780470544037.ch14
  * Shalev-Shwartz et al. [2017] Shalev-Shwartz, S., Shamir, O., Shammah, S.: Failures of Gradient-Based Deep Learning. Preprint at arXiv:1703.07950 (2017)
  * Hu et al. [2023] Hu, F., Angelatos, G., Khan, S.A., Vives, M., Türeci, E., Bello, L., Rowlands, G.E., Ribeill, G.J., Türeci, H.E.: Tackling sampling noise in physical systems for machine learning applications: Fundamental limits and eigentasks. Phys. Rev. X 13, 041020 (2023) https://doi.org/10.1103/PhysRevX.13.041020
  * Dieudonné [1960] Dieudonné, J.: Foundations of Modern Analysis. Academic Press, New York (1960)
  * Lalley [2016] Lalley, S.: Markov Chains: Basic Theory. MarkovChains.pdf (2016)
  * Whitt [2013] Whitt, W.: Stochastic Models I: Markov Chains. lect1029.pdf (2013)
  * Heckerman [2008] Heckerman, D.: In: Holmes, D.E., Jain, L.C. (eds.) A Tutorial on Learning with Bayesian Networks, pp. 33–82. Springer, Berlin, Heidelberg (2008). https://doi.org/10.1007/978-3-540-85066-3_3
  * Polloreno [2023] Polloreno, A.M.: Limits to Reservoir Learning. Preprint at arXiv:2307.14474 (2023)