Graduate School of Information Science and Technology, Osaka University, 1-5 Yamadaoka, Suita-shi, Osaka, Japan e-mail<EMAIL_ADDRESS>
© The Author(s) 2024

H. Suzuki et al. (eds.), *Photonic Neural Networks with Spatiotemporal Dynamics*, https://doi.org/10.1007/978-981-99-5072-0\_10

In this chapter, we will focus on random sampling for this purpose. In statistics and machine learning, Monte Carlo (MC) integration is often used to compute expectations. The problem is to determine the expected value of a function *f* (*x*) over a probability distribution *p* denoted by

$$\mathbb{E}\_p[f(\mathbf{x})] = \int f(\mathbf{x}) p(\mathbf{x}) d\mathbf{x} \,. \tag{1}$$

In MC integration, a sequence of samples *x*(1) , *x*(2) ,..., *x*(*T* ) that follow the probability distribution *p* is generated, and then the expectation is approximated as the sample average:

$$\mathbb{E}\_p[f(\mathbf{x})] \simeq \frac{1}{T} \sum\_{t=1}^T f(\mathbf{x}^{(t)}) \,. \tag{2}$$

For example, MC integration is important for Bayesian inferences. When inferring the parameters θ of a probabilistic model *p*(*x*; θ ) from the data D, we often use the maximum likelihood method taking θ as an estimate that maximizes the likelihood *P*(D; θ ) of the data. Instead, Bayesian inference introduces a prior distribution π(θ ) for the parameters and applies Bayes' rule to obtain the posterior distribution *P*(θ ) ∝ π(θ )*P*(D; θ ) as the inference result. The advantage of assuming distributions in the parameters is that estimates with uncertainty can be obtained.

Because the variability in physical phenomena is often understood as a stochastic behavior, this sampling and MC integration is a promising application of physical phenomena in computation. Usually, a random number generator (RNG) is used for sampling. Typically, a pseudo-RNG is used to generate a sequence of numbers using deterministic computations; however, physical random numbers can also be used. In either case, the RNG must be precisely designed to guarantee the quality of the output. Although the use of stochastic behavior observed in physical phenomena is promising, it is also difficult to precisely control its probabilistic properties.

Herding, proposed by [2, 3], is an algorithm that can be used in the same manner as MC integration. However, it does not use RNGs, so we can expect it to be a possible method for avoiding such difficulties. In this chapter, we review the herding algorithm and its related studies, including the author's previous results. Then we also discuss the prospects of studying herding from perspectives that include the use of physical phenomena in computation.

The remainder of this chapter is organized as follows. In Sects. 1 and 2, we introduce herding, focusing on its aspects as a method ofMC integration. In particular, we introduce the herding algorithm in Sect. 1. In Sect. 2, the herded Gibbs algorithm is introduced as an application of herding.We discuss the improvement of the estimation error and its convergence and review the relevant literature. In Sect. 3, we consider herding from a different perspective, the maximum entropy principle. We explore the relationship between herding and entropy, including the role of high-dimensional nonlinear dynamics in herding. In Sect. 4, based on the discussion in the previous sections, we discuss the prospects of the study of herding, including the aspect of computation with the variability of physical phenomena, and finally, we provide concluding remarks in Sect. 5.

#### **1 Herding: Sample Generation Using Nonlinear Dynamics and Monte Carlo Integration**

In this section, we introduce the herding algorithm. We adopt a different approach to the introduction than that in the original work [2, 3].

#### *1.1 Deterministic Tossup*

We begin with a simple example. Consider a sequence of *T* random variables *x*(1) ,..., *x*(*T* ) ∈ R. LetE*p*[·] and Var*p*[·] denote the expectation and variance, respectively, of a function over a probability distribution *p*. Assume that each marginal distribution is equivalent to the distribution *p* whose expected value is E*p*[*x*] = μ and whose variance is Var*p*[*x*] = *V*. Then, the expected value of the sample mean

$$
\hat{\mu} = \frac{1}{T} \sum\_{t=1}^{T} \mathbf{x}^{(t)} \tag{3}
$$

satisfies E[ ˆμ] = μ and the variance is calculated as

$$\begin{split} \text{Var}\left[\hat{\mu}\right] &= \mathbb{E}\left[\left(\frac{1}{T}\sum\_{t=1}^{T}(\mathbf{x}^{(t)} - \mu)\right)^{2}\right] \\ &= \frac{1}{T^{2}}\sum\_{t=1}^{T}\sum\_{i'=1}^{T}\mathbb{E}[(\mathbf{x}^{(t)} - \mu)(\mathbf{x}^{(t)} - \mu)] \\ &= \frac{1}{T^{2}}\left(\sum\_{t=1}^{T}\mathbb{E}[(\mathbf{x}^{(t)} - \mu)^{2}] + 2\sum\_{k=1}^{T-1}\sum\_{t=1}^{T-k}\mathbb{E}[(\mathbf{x}^{(t)} - \mu)(\mathbf{x}^{(t+k)} - \mu)]\right). \end{split} (4)$$

Assume that *x*(1) ,..., *x*(*T* ) are not independently and identically distributed (i.i.d.) and have a stationary autocorrelation as

$$\mathbb{E}[(\mathbf{x}^{(t)} - \mu)(\mathbf{x}^{(t+k)} - \mu)] = V\rho\_k \,\, , \tag{5}$$

then the variance is calculated as

204 H. Yamashita

$$\text{Var}\left[\hat{\mu}\right] = \frac{V}{T}\left(1 + 2\sum\_{k=1}^{T-1} \left(1 - \frac{k}{T}\right)\rho\_k\right) \tag{6}$$

Considering these *x*(*t*) as the samples inMC integration, this suggests that the variance of the estimate is reduced when the number of samples *T* is large, and it can be further reduced when the autocorrelation ρ*k* is small or especially negative.

Consider a distribution of the binary random variable *x* ∼ *p*π with parameter π, representing a tossup that returns 1 with probability π and 0 with probability 1 − π. This distribution is the well-known Bernoulli distribution with the parameter π. Instead of generating a sequence of i.i.d. samples from *p*π , we consider generating a sequence that follows the deterministic process presented in Algorithm 1 and Fig. 1. According to the update rule (line 3) of the algorithm, w(*t*) is expected to be smaller after the step if *x*(*t*) = 1, and w(*t*) is expected to be larger after the step if *x*(*t*) = 0. Therefore, from the equation in line 2, we expect that this sequence will have negative autocorrelations.

Although the above arguments concern only ρ1, in fact the sample mean μˆ converges very quickly as *O*(1/*T* ). Let *C* = [−1 + π, π) and assume that the initial value w(0) is taken from *C*. Then, one can easily check that w(*t*) ∈ *C* holds for all *t*, for example, from that we can change the variable w(*t*) to w(*t*) to obtain an equivalent system including the rotation of the circle represented as

$$\mathbf{x}^{(t)} = \begin{cases} 1 & (w'^{(t-1)} \le \pi) \\ 0 & (w'^{(t-1)} > \pi) \end{cases},\tag{7}$$

$$w^{\prime(t)} = w^{\prime(t-1)} + \pi \bmod 1\,\,.\tag{8}$$

Adding the update formula (line 3) with respect to *t* = 1,..., *T* , we obtain

$$w^{(T)} - w^{(0)} = T\pi - \sum\_{t=1}^{T} x^{(t)}.\tag{9}$$

Then we evaluate the difference in average as

**Fig. 1** An example trajectory of the deterministic tossup (herding for Bernoulli distribution). In the left panel, off-diagonal dashed lines represent the maps used to update w(*t*) . The points on the diagonal dashed line indicate the w(*t*) for each *t*. The circles represent *x*(*t*+1) = 1 and crosses represent *x*(*t*+1) = 0. The gray area represents the condition for w(*t*) to output *x*(*t*+1) = 1. The right panel represents the trajectory as the function of *t*

$$|\pi - \hat{\mu}| = \frac{1}{T} |w^{(T)} - w^{(0)}|\,\,,\tag{10}$$

which converges as *O*(1/*T* ) because w(0) , w(*T* ) ∈ *C*.

#### *1.2 Herding*

We then generalize the deterministic tossup above to introduce herding. Let X be the sample space and ϕ*m* : X → R for *m* ∈ M be the feature functions defined therein, where M is the set of their indices. Suppose that we have parameters μ*m* ∈ R that specify the target moment values of ϕ*m* over output samples. The empirical moment values can be used as parameters if the dataset is available. We can also consider a situation in which only the aggregated observation is available, e.g., for privacy reasons or because the subject is microscopic so that individual observation is impossible. The problem considered here is the reconstruction of the distribution such that the expected value of the feature Eπ [ϕ*m*(*x*)] is equal to the given moment μ*m*. Herding is an algorithm proposed by Welling [2], which generates a sequence of samples *x*(1) , *x*(2) ,... that are diverse while satisfying the moment condition.

For ease of notation, we denote by*ϕ*(*x*) and*μ*the vectors of features and moments, respectively. We denote by ·, · the inner product of vectors and by *x*≡ √ *x*, *x* the square norm. The overall herding algorithm is presented in Algorithm 2. We refer to *w*(*t*) of the algorithm as the weight vector. As discussed below, herding is an extension of the above deterministic tossup where the sample mean of *ϕ* converges to *μ*. We can also consider the algorithm as a nonlinear dynamical system with the discrete-time variable *t*. As shown in Algorithm 2, the system includes optimization

#### **Algorithm 2** Herding algorithm

1: **for** *t* = 1,..., *T* **do** 2: *x*(*t*) = arg max *x*∈X *w*(*t*−1) , *ϕ*(*x*) 3: *w*(*t*) = *w*(*t*−1) + *μ* − *ϕ*(*x*(*t*) ) 4: **end for** 5: Output *x*(1) ,..., *x*(*T* )

problem as a component. If we obtain the unique solution to the optimization problem, then updating the weight vector is deterministic. However, the system behaves in a random or chaotic manner. This random-like behavior can generate a pseudo-random sequence of samples while satisfying the moment conditions. Thus, the generated sequence that reflects the given information can be used as samples drawn from the background distribution. We discuss the importance of the weight vector behavior later.

There are many studies related to herding; for example, it is extended to include hidden variables throughMarkov probability fields and is applied to data compression by learning the rate coefficients of the dynamical system of herding [4]. It is also combined with the kernel trick to sample from a continuous distribution [5]. However, these are not discussed in detail here.

#### *1.3 Convergence of Herding*

Let us consider the herding algorithm for the Bernoulli distribution. In this situation, the sample space is X = {0, 1} and the features are *ϕ*(*x*) = (φ0(*x*), φ1(*x*)) = (1 − *x*, *x*), whose moments are *μ* = (1 − π, π). The herding algorithm in this situation is equivalent to

$$\mathbf{x}^{(t)} = \begin{cases} 1 & (f^{(t-1)}(\mathbf{l}) \ge f^{(t-1)}(\mathbf{0})) \\ 0 & (f^{(t-1)}(\mathbf{l}) < f^{(t-1)}(\mathbf{0})) \end{cases},\tag{11}$$

$$w\_0^{(t)} = w\_0^{(t-1)} + (1 - \pi) - (1 - \pi)$$

$$\dots^{(t-1)} \quad \dots \quad \dots$$

$$
\hat{w}\_t = w\_0^{(t-1)} - \pi + x\_+ \tag{12}
$$

$$w\_1^{(t)} = w\_1^{(t-1)} + \pi - \ge \,\,,\tag{13}$$

where *f* (*t*−1) (*x*) ≡ *w*(*t*−1) , *φ*(*x*)= w(*t*−1) 0 + (w(*t*−1) 1 − w(*t*−1) 0 )*x*. This updated formula corresponds to the deterministic tossup (Algorithm 1), with the relation (w(*t*) 0 , w(*t*) 1 ) = (−w(*t*) , w(*t*) ); the two algorithms become equivalent if the initial value satisfies w(*t*) 0 = −w(*t*) 1 .

Generalizing and using this relation conversely, we can perform an analysis similar to that of the deterministic tossup in the general case. Suppose *x*(*t*) is chosen such that ϕ*m*(*x*(*t*) ) is small at step *t* and the target value is ϕ*m*(*x*(*t*) )<μ*m*. In this case, the corresponding weight w(*t*) *m* increases as

$$\mathbf{w}^{(t)} = \mathbf{w}^{(t-1)} + \boldsymbol{\mu} - \boldsymbol{\varphi}(\mathbf{x}^{(t)})\,,\tag{14}$$

which corresponds to line 3 of Algorithm 2. Therefore, ϕ*m* becomes more important in the optimization step and ϕ*m*(*x*(*t*+1) ) is expected to be larger in the next step. Similarly, ϕ*m*(*x*(*t*+1) )is expected to be smaller when ϕ*m*(*x*(*t*) )is large. In other words, the sequence ϕ*m*(*x*(*t*) ) is expected to have negative autocorrelation, which makes its average converge to μ*m* faster.

In addition, we obtain a convergence result for herding similar to that of the deterministic tossup. In particular, the average of the features for the generated sequence converges as follows:

$$\frac{1}{T} \sum\_{t=1}^{T} \varphi(\mathbf{x}^{(t)}) \to \boldsymbol{\mu} \,, \tag{15}$$

where the convergence rate is *O*(1/*T* ). This convergence is obtained by using the boundedness of *w*(*T* ) and equation

$$\left(\boldsymbol{\mathfrak{w}}^{(T)} - \boldsymbol{\mathfrak{w}}^{(0)} = T\boldsymbol{\mathfrak{w}} - \sum\_{t=1}^{T} \boldsymbol{\mathfrak{o}}\left(\boldsymbol{\chi}^{(t)}\right)\,. \tag{16}$$

This equation is obtained by summing both sides of (14) for *t* = 1,..., *T* .

Using the optimality for the optimization problem

$$\underset{\boldsymbol{x}\in\mathcal{X}}{\arg\max} \left< \boldsymbol{w}^{(t-1)}, \boldsymbol{\varphi}(\boldsymbol{x}) \right>\,,\tag{17}$$

we can guarantee the boundedness of *w*(*T* ) as the following summarized version of the proof in [2]: Let *R*(*t*) = *w*(*t*) . The change in the value of *R* is calculated as follows:

$$\begin{aligned} \| (\mathcal{R}^{(t)})^2 &= \| \mathfrak{w}^{(t)} \|^2 \\ &= \| \mathfrak{w}^{(t-1)} - \mathfrak{o}(\mathfrak{x}^{(t)}) + \mathfrak{h} \|^2 \\ &= (\mathcal{R}^{(t-1)})^2 - 2 \langle \mathfrak{w}^{(t-1)}, \mathfrak{o}(\mathfrak{x}^{(t)}) \rangle + 2 \langle \mathfrak{w}^{(t-1)}, \mathfrak{w} \rangle + \| \mathfrak{o}(\mathfrak{x}^{(t)}) - \mathfrak{w} \|^2 . \end{aligned} (18)$$

Note that *x*(*t*) is obtained by the optimization problem in (17); we define the following quantity

$$A = \min\_{\|\hat{\Psi}\|=1} \left( \langle \tilde{\boldsymbol{w}}, \boldsymbol{\varphi}(\hat{\boldsymbol{x}}) \rangle - \langle \tilde{\boldsymbol{w}}, \boldsymbol{\mu} \rangle \right) \,, \tag{19}$$

where *x*ˆ is defined as *x*ˆ = arg max *x w*˜ , *ϕ*(*x*). From the optimality of *x*ˆ, we obtain *A* > 0 under a mild assumption of *μ*. We also assume that the third term in (18) is bounded; that is, *ϕ*(*x*(*t*) ) − *μ* 2 ≤ *B* holds for some *B* ∈ R. This holds if the functions *ϕ* are bounded. Then, we obtain

$$(R^{(t)})^2 \le (R^{(t-1)})^2 - 2R^{(t-1)}A + B \,. \tag{20}$$

In other words, if *R*(*t*−1) > *R* holds for *R* ≡ *B*/2*A*, then *R*(*t*) is decreasing; *R*(*t*) < *R*(*t*−1) . On the other hand, if *R*(*t*−1) ≤ *R*, then *R*(*t*) is bounded. Thus, there are *R* ∈ R such that *R*(*t*) ≤ *R* always holds if *R*(0) ≤ *R* .

#### **2 Herded Gibbs: Model-Based Herding on Spin System**

Herding typically requires the feature moments*μ*as the inputs. However, in sampling applications, the probability model π is often available instead. In this section, we describe herded Gibbs (HG), an extension of herding that can be applied to such model-based situations.

#### *2.1 Markov Chain Monte Carlo Method*

Sampling independently from a general probability distribution π is difficult in general; thus the Markov chain Monte Carlo (MCMC) method is often used. In MCMC, we consider the state variable *x*(*t*) that is repeatedly updated to be the next state *x*(*t*+1) , where this transition is represented as a Markov chain. By appropriately designing this transition, the probability distribution of *x*(*t*) converges to π. In addition, the two samples *x*(*t*) and *x*(*t*+*T* ) become nearly independent for sufficiently large *T* ; the measure of such a time difference *T* is called the mixing time. If the transitions are designed such that the mixing time is short, then the MCMC becomes more effective for MC integration.

MCMC transitions are often designed to be local, i.e., the state *x*(*t*) moves only in its neighborhood at each step. In such a case, the mixing time strongly depends on the energy landscape of probability distribution π to be sampled. The energy landscape is represented as the graph of *E*(*x*), where the probability distribution is represented by the Gibbs-Boltzmann distribution, π(*x*) ∝ exp(−*E*(*x*)). In particular, if there are large energy barriers in the phase space, the mixing time increases because it takes time for *x*(*t*) to pass over them.

A commonly used method for designing MCMC is the Metropolis-Hastings method [6, 7], and it is a good example of a local transition. In the Metropolis-Hastings method, given the current state *x*, the state *y* is generated according to the proposal distribution *qx* (*y*). Typically, it is designed so that *y* is in the neighborhood of *x*. The proposed state *y* is accepted as the next state with an acceptance probability α ∈ [0, 1] and is rejected with probability 1 − α. If rejected, state *x* is kept as the next state. If the proposal distribution is symmetric, i.e., *qx* (*y*) = *qy* (*x*) ∀*x*, *y*, and the acceptance probability is defined as α = min(1, exp(−*E*(*y*))/ exp(−*E*(*x*))), the chain has π(*x*) ∝ exp(−*E*(*x*)) as its stationary distribution. Transition probabilities are defined such that a transition to *y* with higher energy *E*(*y*) has a lower acceptance probability.

Gibbs sampling [8] is another type of MCMC method. Let us denote by *p*(·|·)the conditional distribution for a distribution *p*. Suppose π is a probability distribution of *N* variables *x*1,..., *xN* , and we can compute the conditional distribution π(*xi* | *x*−*i*) for each *i*, where *x*−*i* is the vector of all variables except *xi* . Then, we can construct a Markov chain as presented in Algorithm 3.

In particular, in this chapter, we consider the following spin system: Let us consider *N* binary random variables, *x*1,..., *xN* ∈ {0, 1}. Suppose that the random variables V = {*x*1,..., *xN* } form a network *G* = (V, E), where E denotes the set of edges, and the joint distribution π is represented as

$$
\pi(\mathbf{x}) = \frac{1}{Z} \exp\left(-E(\mathbf{x})\right),
\tag{21}
$$

$$E(\mathbf{x}) = -\sum\_{(i,j)\in\mathcal{E}} W\_{ij}\mathbf{x}\_i\mathbf{x}\_j - \sum\_{i=1}^{N} b\_i\mathbf{x}\_i \; , \tag{22}$$

where *Z* is the normalization factor that makes the sum of the probabilities equal to one;

$$Z = \sum\_{\mathbf{x}} \exp(-E(\mathbf{x})) \,. \tag{23}$$

This probability model is widely known as the Boltzmann machine (BM), which is also discussed in Chap. 2.

For simplicity, let us suppose that there is no bias term, such that

210 H. Yamashita

$$E(\mathbf{x}) = -\sum\_{(i,j)\in\mathcal{E}} W\_{ij} x\_i x\_j \ . \tag{24}$$

Let N(*i*) be the index set of neighboring variables of *xi* on *G* and *x*N(*i*) be the vector of the corresponding variables. For the BM, the conditional distribution can be easily obtained as

$$\pi(\mathbf{x}\_i \mid \mathbf{x}\_{-i}) = \pi(\mathbf{x}\_i \mid \mathbf{x}\_{\mathcal{N}(i)}) = \frac{1}{Z\_i} \exp\left(\mathbf{x}\_i \left(\sum\_{j \in \mathcal{N}(i)} W\_{ij} \mathbf{x}\_j\right)\right),\tag{25}$$

where *Zi* is also the normalization factor. Therefore, sampling from this distribution can be easily performed using Gibbs sampling.

#### *2.2 Herded Gibbs*

Herding can be used as a deterministic sampling algorithm, but it cannot be applied directly to the BM because the input parameter*μ*is typically unavailable. By combining Gibbs sampling with herding, a deterministic sampling method for BMs called the herded Gibbs (HG) algorithm [9] is obtained. The structure of the algorithm is the same as that of Gibbs sampling, but the random update step for each variable *xi* is replaced by herding, as presented in Algorithm 4. The variable neighborhood *x*N(*i*) of *xi* has 2|N(*i*)| configurations, so let us give indices for such configurations. HG uses weight variables, denoted by w*i*,*j* for the *i*th variable *xi* and the *j*-th configuration of the variable neighborhood *x*N(*i*). We can compute the corresponding conditional probability given that *x*N(*i*) takes the *j*th state denoted by π*i*,*j* ≡ π(*xi* = 1 | *x*N(*i*) = *j*).

HG is not only a deterministic variant of Gibbs sampling but is reported to have better performance than Gibbs sampling [9]. For a function *f* of the spins *x* in the BM, let us consider estimating the expected value Eπ [ *f* (*x*)]. Theoretically, for the BM on the complete graph, the error of the estimate by HG decreases at a rate of *O*(1/*T* ), whereas it decreases at a rate of *O*(1/ √*T* ) for random sampling. Experimental results show that HG outperforms Gibbs sampling for image processing and natural language processing tasks.

When a BM is used as a probabilistic model, its parameters must be learned from data. The learning procedure is expressed as

$$W\_{ij} \leftarrow W\_{ij} + \mathbb{E}\_{\text{data}}[\mathbf{x}\_i \mathbf{x}\_j] - \mathbb{E}\_{\text{model}}[\mathbf{x}\_i \mathbf{x}\_j] \,, \tag{26}$$

where Edata[·] is the mean of the training data and Emodel[·] = Eπ [·] is the expected value in the model (21) with the current parameters *Wi j* . The exact calculation of the third term is typically difficult owing to the exponentially increasing number of terms, but MCMC can be used to estimate it [10–12]. HG can also be applied to the learning process via its estimation; the variance reduction by HG can have a positive effect on learning, as demonstrated for BM learning on handwritten digit images [13].

#### *2.3 Sharing Weight Variables*

The original literature on HG [9] also introduces the idea of "weight sharing." For a variable *xi* , let *j*, *j* be indices of different configurations with equal conditional probabilities π*i*,*j* = π*i*,*j* . For example, this can occur when the coupling coefficients are restricted to *Wi*,*j* ∈ {0, ±*c*} for *c* ∈ R; thus the conditional probability is determined only by counting neighboring variable values. In weight sharing, such configurations are classified into groups, and those in the same group, indexed by *y*, share a weight variable w*i*,*y* . We will also refer to these versions, which are based on general graphs and use weight sharing, as HG.

Eskelinen [14] proposed a variant of this algorithm called discretized herded Gibbs (DHG), which is presented in Algorithm 5. It uses *B* disjoint intervals of the conditional probability value dividing the unit interval [0, 1]. Let us call them bins and denote by *Cy* the *y*-th interval. At the time of the update, the conditional probability π(*xi* = 1 | *x*N(*i*)) is approximated by the representative value π˜*i*,*y* ∈ *Cy* , where *y* is the index of the interval to which the configuration *x*N(*i*) belongs. The weight variable used for the update is shared among several spin configurations with similar conditional probability values. The probability value used for the update is replaced by the representative value, but this replacement introduces an error in the distribution of the output samples. A trade-off exists between the computational complexity proportional to *B* and the magnitude of the error.

These ideas aim to reduce computational complexity by reducing the number of weights. However, Bornn et al. [9] reported a performance improvement for image

restoration using the BM, and Eskelinen [14] reported a reduction in the estimation error in the early iterations for the BM on a small complete graph.

#### *2.4 Monte Carlo Integration Using HG*

On complete graphs, HG has been shown to be consistent; the estimates converge to the true value associated with the target distribution π. On general graphs, however, this is not always the case. In particular, the error decay for HG in the general case has two characteristics: faster convergence in the early iterations and convergence bias in the later iterations. These are explained by the "weight sharing" and the temporal correlation of the weights.

Specifically, the accuracy of the sample approximation of the function *f* with HG was evaluated as follows [13]: For HG with a target BM π (21), let *P* be the empirical distribution of the obtained samples. Let the error be defined as

$$D \equiv \mathbb{E}\_P[f(\mathbf{x})] - \mathbb{E}\_\pi[f(\mathbf{x})] \,. \tag{27}$$

The magnitude of *D* is evaluated as

$$|D| \le \lambda\_{\text{cor}} D^{\text{cor}} + \lambda\_{\text{herding}} D^{\text{hering}} + \lambda\_{\text{approx}} D^{\text{approx}} + \lambda\_{\text{z}} D^{\text{z}}.\tag{28}$$

This is obtained by decomposing the estimation error with respect to a variable, where each term has the following meaning: Let *x* = *xi* with fixing *i*, and let *z* = *x*−*i* be all the variables except *xi* . Let *y* be the index of the weight used to generate *xi* (line 3 of Algorithm 4 or line 4 of Algorithm 5) at each step. *D*approx is the error term resulting from the replacement of the conditional probability π(*x* = 1 | *z*) by the representative value π˜*i*,*y* , which occurs in the DHG. *D*z is the error in the distribution of *z*, namely, the joint distribution of the variables other than *xi* , and expected to depend largely on the mixing time of the original Gibbs sampling.

The other two, *D*cor and *D*herding, are the terms most relevant to the herding dynamics. *D*herding is the term corresponding to the herding algorithm for the Bernoulli distribution or the deterministic tossup that decays by *O*(1/*T* ). *D*cor is a non-vanishing term. Thus, the error decay is characterized by the fast decay in early iterations dominated by *D*herding and its stagnation in later iterations dominated by *D*cor. *D*cor is the term owing to the temporal correlation of the weight variables in the HG and is represented by the following equation:

$$D^{\rm cor} = \sum\_{\mathbf{y}} P(\mathbf{y}) \sum\_{\mathbf{x}, \mathbf{z}} \left| P(\mathbf{x}, \mathbf{z} \mid \mathbf{y}) - P(\mathbf{x} \mid \mathbf{y}) P(\mathbf{z} \mid \mathbf{y}) \right| \,. \tag{29}$$

Let us consider the ideal case that the weight index *y* and the conditional distribution π(*x* | *z*) have a one-to-one correspondence. For DHG, this means that there are enough many small bins. Then, for the target distribution, *x* and *z* conditioned on *y* are independent; π(*x*, *z* | *y*) = π(*x* | *y*)π(*z* | *y*). This is because π(*x* | *y*) becomes constant on *z*. For the output distribution, however, this independence does not hold even under this assumption because the internal state w*i*,*y* determines the value of *x* and also has the temporal correlation with *z*. *D*cor evaluates the correlation between the conditional distributions of *x* and *z*.

This temporal correlation is mainly caused by the deterministic nature of the herding dynamics. To mitigate this, the algorithm can be modified by introducing some stochasticity in the transitions. This has been shown to solve the bias problem [13], but at the cost of a degradation in the estimation for small *T* .

#### *2.5 Rao-Blackwellization*

Let us further discuss this analysis with a more concrete example. Let us fix *i* as before. Let *f* (*x*) = *xi* and consider the accuracy of estimating the expected value Eπ [ *f* (*x*)] = Eπ [*xi*]. The estimation error is bounded as

$$\begin{split} |D| = |\mathbb{E}\_P[\mathbf{x}\_i] - \mathbb{E}\_\pi[\mathbf{x}\_i]| &= \left| \sum\_z P(\mathbf{z})P(\mathbf{x} = \mathbf{l} \mid \mathbf{z}) - \sum\_z \pi(\mathbf{z})\pi(\mathbf{x} = \mathbf{l} \mid \mathbf{z}) \right| \\ &\leq \sum\_z P(\mathbf{z}) \left| P(\mathbf{x} = \mathbf{l} \mid \mathbf{z}) - \pi(\mathbf{x} = \mathbf{l} \mid \mathbf{z}) \right| \\ &+ \sum\_z \left| P(\mathbf{z}) - \pi(\mathbf{z}) \right| \pi(\mathbf{x} = \mathbf{l} \mid \mathbf{z}) \ . \end{split} \tag{30}$$

Under some assumptions, the first term corresponds to *D*herding and the second term corresponds to *D*z . Therefore, we can evaluate the error decay as the first term decaying by *O*(1/*T* ) as in herding but is eventually dominated by the second term for a large *T* .

The results were compared with those of a similar algorithm as follows: First, an i.i.d. sample sequence *x*(*t*) is generated via random sampling. Subsequently, the expected value of *f* (*x*) = *xi* is approximated as the sample mean of the conditional probability π(*x* = 1 | *z*). In this case, the approximation error is bounded above as

$$|D| \le \sum\_{\mathbf{z}} |P(\mathbf{z}) - \pi(\mathbf{z})| \pi(\mathbf{x} = 1 \mid \mathbf{z}) \,, \tag{31}$$

which is equal to the term that is dominant in the HG case. The substitution of the conditional mean (or distribution) into an estimator is sometimes referred to as Rao-Blackwellization [15, 16], after the Rao-Blackwell theorem that guarantees the improvement in estimation accuracy from this substitution. In other words, herding can be viewed as a sample-based estimation method that achieves error reduction through Rao-Blackwellization, although it is implemented implicitly.

However, HG is not equivalent to this; the evaluation of the decaying error simultaneously holds for any other *x j* (*j* = *i*) without changing the algorithm. If the conditional probability π(*x* = 1 | *z*) is considered as a function of *z*, this function itself can have an approximation similar to that of (28). That is, we can decompose this function as

$$\pi(\mathbf{x} = 1 \mid \mathbf{z}) = \hat{f}(\mathbf{z}) + \sum\_{j \neq i} \beta\_j x\_j \; \; , \tag{32}$$

and evaluate each term in the summation by (28) for *f* (*x*) = *x j* . Then we find that the first term in (32) is dominant, which leads to the improvement of the error evaluation in HG.

We can also use this formula for i.i.d. random sampling to improve the estimation. However, in this case, it is necessary to determine the appropriate β*j* by estimating from the random samples obtained. However, in the case of herding, because this evaluation is valid for any β*j* , selecting β*j* is not necessary. Extending this discussion, we may be able to analyze the accuracy of estimation by herding in general; however, this is left open for future study.

#### **3 Entropic Herding: Regularizer for Dynamical Sampling**

In the previous section, we discussed the characteristics of the herding algorithm as a Monte Carlo numerical integration algorithm. In particular, there is a negative autocorrelation in the weight variables and samples, which is important for numerical integration. However, this does not cover all the characteristics of herding as a highdimensional nonlinear dynamics. This section discusses the connection between herding and the maximum entropy principle.

#### *3.1 Maximum Entropy Principle*

The maximum entropy principle [17] is a common approach to statistical inference. It states that the distribution with the greatest uncertainty among those consistent with the information at hand should be used. As the name implies, entropy is often used to measure uncertainty. Specifically, for a distribution with a probability mass (or density) function *p*, the (differential) entropy is defined as

$$H(p) = \mathbb{E}\_p[-\log p(\mathbf{x})]\,. \tag{33}$$

Furthermore, we assume that information is collected on features ϕ*m* : X → R indexed with *m* ∈ M, where X is the sample space and M is the set of feature indices. Assume further that the collected data is available as the mean μ*m* of each feature. According to the maximum entropy principle, the estimated distribution *p* should satisfy the condition that the expected value E*p*[ϕ*m*(*x*)] of the feature values is equal to the given value μ*m*. Therefore, the maximum entropy principle can be explicitly expressed as the following optimization problem:

$$\begin{aligned} \underset{p}{\text{maximize}} \quad & H(p) \\ \text{subject to} \quad \mathbb{E}\_p[\varphi\_m(\mathbf{x})] = \mu\_m \quad \forall m \in \mathcal{M} \end{aligned} \tag{34}$$

We obtain the condition that the solution should satisfy as follows: For simplicity, let X = {*x*1,..., *xN* } be a discrete set and, for each *i*, let *pi* denote the probability of the *i*-th state *xi* . The gradient of *H* is

$$\frac{\partial}{\partial p\_i} H(p) = -\log p\_i - 1\ . \tag{35}$$

Using the Lagrange multiplier method, we obtain the condition

$$p\_i \propto \exp\left(-\sum \theta\_m \varphi\_m(\mathbf{x}\_i)\right),\tag{36}$$

where θ*m* corresponds to the Lagrange multiplier for the moment condition E*p*[ϕ*m*(*x*)] = μ*m*. For continuous X, we can make a similar argument using the functional derivative.

In general, a family of distributions of the form (36) where the parameters θ*m* are unconstrained is called an exponential family, and is also called Gibbs-Boltzmann distribution. Optimizing the parameters θ*m* requires a learning algorithm similar to that of the Boltzmann machine (21), which is generally computationally difficult.

For example, if X = R and (ϕ1(*x*), ϕ2(*x*)) = (*x*, *x* 2), the maximum entropy principle gives the normal distribution as

216 H. Yamashita

$$p(\mathbf{x}) \propto \exp\left(-\theta\_1 \mathbf{x} - \theta\_2 \mathbf{x}^2\right)$$

$$\propto \exp\left(-\frac{(\mathbf{x} - m)^2}{2\sigma^2}\right),\tag{37}$$

where (*m*,σ) = (−θ1/2θ2, √1/2θ2).

If we use X ∈ {0, 1}*N* and ϕ*m*(*x*) = *xi x j* for all pairs *i* < *j* where *m* = (*i*, *j*), we obtain

$$p(\mathbf{x}) \propto \exp\left(-\sum\_{i$$

which is identical to the Boltzmann machine (21) without bias on the fully connected graph.

The output of the herding algorithm is also such that the moment condition E*p*[ϕ*m*(*x*)] = μ*m* ∀*m* ∈ M is satisfied in the limit of *T* → +∞, and is expected to be diversified due to the complexity of the herding dynamics. Therefore, we can expect the output sequence to follow, at least partially, the maximum entropy principle. The original literature on herding [2, 3] also describes its motive in the context of the maximum entropy principle.

#### *3.2 Entropic Herding*

We observed that herding is an algorithm closely related to the maximum entropy principle. Entropic herding, which is described in this section, is an algorithm that incorporates this principle in a more explicit manner.

Let us consider the same features ϕ*m* and target means μ*m* as above. Pseudocode for entropic herding is provided in Algorithm 6. Here, we introduce scale parameters *m* ≥ 0 for each condition *m* ∈ M, which are used to control the penalty for the condition E*p*[ϕ*m*] = μ*m* as described below. In addition, we introduce step-size parameters ε(*t*) ≥ 0 to the algorithm.

Similar to the original herding algorithm, entropic herding is an iterative process and the time-varying weight *am* for each feature is included in the system. Each

**Algorithm 6** Entropic herding

1: Choose an initial distribution *r*(0) and let η(0) *m* = E*r*(0) [ϕ*m*(*x*)] and *a*(0) *m* = *m*(η(0) *m* − μ*m*) 2: **for** *t* = 1,..., *T* **do** 3: Solve min*q*∈Q *m*∈M *a*(*t*−1) *m* E*q* [ϕ*m*(*x*)] − *H*(*q*) (Equation (39)) 4: Let *r*(*t*) be the minimizer

5: Update *a*(*t*) *m* = *a*(*t*−1) *m* + ε(*t*) *m*(E*r*(*t*) [ϕ*m*(*x*)] − μ*m*) − *a*(*t*−1) *m* (Equation (40))

6: **end for**

iteration of the algorithm, indexed by *t*, consists of two steps: the first step solves the optimization problem and the second step updates the parameters based on the solution. Unlike the original herding algorithm, the entropic herding algorithm outputs a sequence of distributions *r*(1) ,*r*(2) ,... instead of points. The two steps for each time step, derived later, are as follows:

$$r^{(t)} = \underset{q \in \mathcal{Q}}{\arg\min} \left( \left( \sum\_{m \in \mathcal{M}} a\_m^{(t-1)} \eta\_m(q) \right) - H(q) \right), \tag{39}$$

$$a\_m^{(t)} = a\_m^{(t-1)} + \varepsilon^{(t)} \left(\Lambda\_m(\eta\_m(r^{(t)}) - \mu\_m) - a\_m^{(t-1)}\right) \,, \tag{40}$$

where Q is the (sometimes parameterized) family of distributions and η*m*(*q*) ≡ E*q* [ϕ*m*(*x*)] is the feature mean for the distribution *q*.

Obtaining the exact solution to the optimization problem (39) is often computationally intractable. We can restrict the candidate distributions Q or allow suboptimal solutions of (39) to reduce the computational cost.

#### *3.3 Dynamical Entropic Maximization in Herding*

Entropic herding is derived from the minimization problem min*p* L(*p*) with the following objective function:

$$\mathcal{L}(p) = \frac{1}{2} \left( \sum\_{m \in \mathcal{M}} \Lambda\_m (\eta\_m(p) - \mu\_m)^2 \right) - H(p) \,. \tag{41}$$

The optimal solution of (41) is expressed with parameter *θ* ∗ as

$$p\_{\theta^\*} (\mathbf{x}) \propto \exp \left( - \sum\_{m \in \mathcal{M}} \theta\_m^\* \varphi\_m (\mathbf{x}) \right) \,. \tag{42}$$

This has the same form as the distribution obtained using the maximum entropy principle (36). However, parameter *θ* ∗ does not coincide with the optimal solution *θ* expressed in (36). Specifically, it satisfies the following equation:

$$
\boldsymbol{\theta}\_{m}^{\*} = \Lambda\_{m}(\eta\_{m}(\boldsymbol{p}\_{\boldsymbol{\theta}^{\*}}) - \boldsymbol{\mu}\_{m})\,. \tag{43}
$$

It roughly implies that the moment error becomes smaller when *m* is large.

Entropic herding can be interpreted as an approximate optimization algorithm for this problem with a restricted search space. Suppose it only considers *P* that has the following form:

$$P(\mathbf{x}) = \sum\_{t=0}^{T} \rho\_t r^{(t)}(\mathbf{x})\,,\tag{44}$$

where ρ*t* are the fixed component weights and each component *r*(*t*) is determined sequentially using the results of the previous steps *r*(0) ,...,*r*(*t*−1) .

Let L˜ be the function obtained by replacing the entropy term *H*(*p*) of L with the weighted average of the entropies of each component *H*˜ ≡ *T t*=0 ρ*t H*(*r*(*t*) ):

$$\tilde{\mathcal{L}} = \frac{1}{2} \left( \sum\_{m \in \mathcal{M}} \Lambda\_m \left( \eta\_m(P) - \mu\_m \right)^2 \right) - \sum\_{t=0}^T \rho\_t H(r^{(t)}) \,. \tag{45}$$

Since *H* ≥ *H*˜ holds due to the convexity of the entropy, L˜ is an upper bound on L. Equation (39) is obtained by optimizing *r*(*t*) with fixing the components *r*(0) ,...,*r*(*t*−1) to minimize L˜ (*P*(*t*) ), where *P*(*t*) ∝ *t t* =0 ρ*tr*(*t* ) is the tentative solution. The optimal condition is equivalent to (39), where the parameters are represented using the tentative solution *P*(*t*) as

$$a\_m^{(\iota -1)} = \Lambda\_m(\eta\_m(P^{(\iota -1)}) - \mu\_m) \,. \tag{46}$$

If we choose ρ*t* appropriately, the update formula (40) can be derived by (46) and a recursive relation between the tentative solutions for consecutive steps, *P*(*t*−1) and *P*(*t*) .

We expect that the formalization of entropic herding will help theoretical studies of the herding algorithm. Equation (41) contains the entropy term *H*(*p*), and we need to make it larger. Entropic herding does so through the following two mechanisms:

We rely on entropy maximization to reconstruct the target distribution, especially to reproduce the distributional characteristics not directly included in the inputs *ϕ*(*x*) and *μ*. In the case of the original non-entropic herding algorithm, this depends entirely on the complexity of the dynamics, which is difficult to analyze accurately. Entropic herding explicitly incorporates this concept.

The added entropy term regularizes the output sample distribution. From (41), the error of η*m* is small when *m* is large; however, the absolute value of *am* defined in (46) becomes large, and thus, the entropy values of the output components *r*(*t*) are expected to be small. Thus, a trade-off exists in the choice of *m*.

Additionally, setting *m* to +∞ implies ignoring the entropy terms in (39). This is equivalent to the original herding algorithm. That is, it is an extremum of entropic herding such that only (b) implicit diversification is used for entropy maximization.

#### *3.4 Other Notes on Entropic Herding*

Entropic herding yields the mixture distribution expressed in (44). The required number of components *T* of the mixture distribution is large when the target distribution is complex. In this case, the usual parameter fitting requires the simultaneous optimization of *O*(*T* ) parameters, which causes computational difficulties. On the other hand, entropic herding can determine *r*(*t*) sequentially. Therefore, the number of parameters to be optimized in each step can be kept small, which simplifies the implementation and improves its numerical stability.

In general, the optimization (39) is non-convex. Therefore, a local improvement strategy alone may lead to a local optimum. One possible solution to this problem is to prepare several candidate initial values. Unlike normal situations in optimization, herding solves optimization problems that differ only in their weights repeatedly. As a result, a good initial value can again become a good initial value again in later steps. Furthermore, the candidate initial values themselves can be improved during the algorithm. This introduces dynamics into them, which has been shown to exhibit interesting behavior [18].

#### **4 Discussion**

We reviewed the herding algorithm that exploits negative autocorrelation and the complexity of nonlinear dynamical systems. In this section, we present two perspectives on herding with a review of related studies.

#### *4.1 Application of Physical Combinatorial Optimization Using Herding*

To use physical phenomena in devices, we must precisely control them, particularly for computation. However, a certain amount of variability usually remains in both fabrication and operation, which hinders their practical use in computation. For example, quantum mechanics inherently involves probabilistic behavior, and current devices for quantum computing are also susceptible to thermal noise. Given that information processing in the human brain occurs through neurons that are energy efficient but noisy, perfect control may not be an essential element of information processing.

Sampling is a promising direction for the application of physical systems in computation, because such variability can be used as a source of complexity. However, completely disorganized behavior is not useful, and the system should be balanced with appropriate control; herding can be the basis for such control. Specifically, this incompletely controlled noisy device can be applied to the optimization step of the herding algorithm. Although the output may deviate from the optimal solution owing to the variability, this is compatible with the formulation of entropic herding, in which a point distribution is not expected. The gap of the output from the ideal solution may affect *R*, the upper bound of the weight norm *w* used in Sect. 1.3, but a larger *R* may be acceptable because it can be handled by increasing the number of samples with a trade-off in the overall cost of the algorithm as long as *R* does not diverge to infinity.

Conventional computer systems have important advantages, such as the versatility of the von Neumann architecture, established communication standards, and the efficiency from the long history of performance improvements. Therefore, even if a novel device that achieves a breakthrough in computation, it is expected to eventually be used as a hybrid system coupled to conventional computers, in which we expect that herding can serve as the interface. For example, herding and optimization with quantum computers such as quantum annealing (QA) [19, 20] or quantum approximate optimization algorithm (QAOA) [21] can be a good combination. As shown in Sect. 1, herding reduces the variance of the estimate through negative autocorrelation. However, this does not mean that the mixing time of MCMC is improved. For herding, when the objective function in the optimization step is multimodal and steep, finding the optimal solution should be difficult even for herding. When an optimization method that can overcome this difficulty, such as quantum annealing, is realized, its combination with herding is expected to lead to further improvements.

#### *4.2 Generative Models and Herding*

Among the latest machine learning technologies, deep generative models are one of the areas that have attracted the most attention in recent years. In machine learning, a generative model is a mathematical model that represents the process of generating points from the sample space, which represents the set of possible targets, such as images or sentences. For example, flow-based models [22] and diffusion models [23] are among the most famous deep generative models, and other popular methods such as variational autoencoders (VAE) [24] and generative adversarial networks (GAN) [25] use generative models as components. These models are known for their amazing performance in AI, but their behavior with a massive number of tunable parameters is generally a black-box.

One of the unique approaches used in machine learning and other new mathematical techniques that take full advantage of computer technology is the use of *implicit representation* of probability distributions through black-box generative processes. For example, the stochastic gradient method, which is the basis of deep learning, requires only its generative method as an algorithm, not a density function in analytic form. When there is a large amount of data, randomly selecting from it can be considered as generating a point from a true generative process. Therefore, it is possible to implement the minimization of the objective function defined using an inaccessible true probability distribution. This is also exemplified by the famous reparameterization trick for VAE [24], where the variational model to be optimized is represented as a generative model instead of an explicit density function. Deep reinforcement learning is another good example, where the probability distribution of the possible time series of agent and environment states is represented as records of past trials, and the agent is trained using randomly sampled records of past experience [26]. This way of using the models increases their affinity with evolving computing technologies and the big-data accumulated from the growing Internet.

On the other hand, one of the main purposes of using more traditional *explicit modeling* of the density functions using parameterized analytic forms is to correlate model parameters and knowledge. The simplest example of this is regression analysis, as expressed by *Y* ∼ β*X*. Here, the probability density function representing the model is fitted to the data, and whether or not the model parameter β is 0 is interpreted as an indication of whether the variables *X* and *Y* are related. We can also use a mathematical model to represent knowledge or hypotheses about a subject and use it to mathematically formalize the arguments based on that knowledge or hypotheses.

If we look at herding from the above perspective, we can find a unique feature of herding. Herding is the process of generating a sequence of points and can be thought of as generative models that implicitly represent a probability distribution. In addition, another important feature is that the algorithmic parameters *ϕ*(*x*) and *μ* can be directly associated with knowledge as parameters of a mathematical model. Herding is a combination of implicit representation and explicit modeling. Although random sampling methods such as MCMC can be used to generate a set of points in the same manner as herding, the process is derived from the explicit model of the distribution, as opposed to the implicit definition of probability distributions as in herding.

#### **5 Conclusion**

In this chapter, we reviewed the herding algorithm, which can be used in the same manner as random sampling, although it is implemented as a deterministic dynamical system. Herding is expected to play an important role in the application of physical phenomena in computation because of the two aspects discussed in this chapter: a Monte Carlo numerical integration algorithm in Sects. 1 and 2, which are related to the mathematical advantages of the herding algorithm, and the connection to the maximum entropy principle in Sect. 3, which is related to the application of physical phenomena with herding. In addition, as shown in Sect. 4, herding is an interesting algorithm not only from these perspectives. We expect that theoretical studies of herding and its applications will continue in both directions.

**Acknowledgements** The author would like to thank Hideyuki Suzuki and Naoki Watamura for their valuable discussions and comments on the manuscript.

#### **References**

**Open Access** This chapter is licensed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license and indicate if changes were made.

The images or other third party material in this chapter are included in the chapter's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the chapter's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.

**Photonic Reservoir Computing**

## **Reservoir Computing Based on Iterative Function Systems**

**Suguru Shimomura**

**Abstract** Various approaches have been proposed to construct reservoir computing systems. However, the network structure and information processing capacity of these systems are often tied to their individual implementations, which typically become difficult to modify after physical setup. This limitation can hinder performance when the system is required to handle a wide spectrum of prediction tasks. To address this limitation, it is crucial to develop tunable systems that can adapt to a wide range of problem domains. This chapter presents a tunable optical computing method based on the iterative function system (IFS). The tuning capability of IFS provides adjustment of the network structure and optimizes the performance of the optical system. Numerical and experimental results show the tuning capability of the IFS reservoir computing. The relationship between tuning parameters and reservoir properties is discussed. We further investigate the impact of optical feedback on the reservoir properties and present the prediction results.

#### **1 Introduction**

An artificial neural network (ANN) is a brain-inspired computing model and contributes to a wide field of information processing including image classification and speech recognition [1]. The ANN is represented by a network structure connected by weighted links. By optimizing the weight of the connections, ANNs have capabilities for desired information processing [2]. However, the optimization requires updating all connections in ANNs, and it is difficult to realize large-scale ANNs. Reservoir computing, which is a kind of recurrent neural network for processing time-series data emerged [3]. Typical models of reservoir computing are divided into an echo state network (ESN) and a liquid state machine [4, 5]. The idea of ESN has been employed for hardware implementation of recurrent neural network, and various architectures have been proposed (see Chap. 13). This chapter focuses on reservoir

S. Shimomura (B)

Graduate School of Information Science and Technology, Osaka University, 1–5, Yamadaoka, Suita, Osaka, Japan e-mail<EMAIL_ADDRESS>
© The Author(s) 2024

H. Suzuki et al. (eds.), *Photonic Neural Networks with Spatiotemporal Dynamics*, https://doi.org/10.1007/978-981-99-5072-0\_11

**Fig. 1** Model of reservoir computing

computing based on the ESN. A model of reservoir computing is shown in Fig. 1. An echo state machine consists of three layers: input, reservoir, and output layers. Nodes in the reservoir layer are connected, and the structure is a recurrent network. A state of the nodes in the reservoir layer at time *t*, *X*(*t*), is updated by

$$X(t+1) = f[\mathbf{W}\_{\text{res}}X(t) + \mathbf{W}\_{\text{in}}\boldsymbol{\mu}(t)],\tag{1}$$

where *u*(*t*) is the input signal at time *t*, and *f* is a nonlinear function such as a hyperbolic tangent and a sigmoid function. Each component of *X*(*t*) is transferred to the other nodes in the reservoir layer according to the connecting weight *W*res. After adding with the weighted input signals *W*in*u*(*t*), the nonlinear function is applied, and the next state is updated as *X*(*t* + 1). The connection weights *W*in between the input and reservoir layers and *W*res in the reservoir layer are fixed and not updated in learning process. The ESN is optimized by a linear regression of weights *W*out between the reservoir and the output layers. Owing to simple structure of the echo state network and low computational processing in the learning process, reservoir computing can be implemented as hardware.

Reservoir computing is an intriguing and dynamic research field, offering a wide range of possibilities for hardware implementation by leveraging diverse types of materials and phenomena [6]. To construct a reservoir computing system, it is essential to design and implement a reservoir as the hardware component, and its response must satisfy Eq. 1. Thus far, various types of reservoir computing systems with individual properties of the utilized phenomenon have been proposed. For instance, the dynamic motion of soft materials has been used to determine the response of a reservoir [7]. The interaction of spin-torque oscillators based on spin waves provides small-scale reservoir devices [8]. In the field of optics, recurrent network circuits implemented on silicon chips and time-delayed feedback loop systems utilizing optical fibers have been successfully employed as reservoirs [9]. Individual reservoirs exhibit a specific reservoir property, which is related to the prediction performance of time-series data [10]. For example, the coupling matrix *W*res shown in Eq. 1 is determined by the characteristics of utilized materials and phenomenon. Consequently, the performance may decrease depending on the prediction task. Tuning of the coupling matrix *W*res is crucial in optimizing the prediction performance of a reservoir computing system, allowing it to effectively address a wide range of problems. To achieve this, the integration of the tuning function within an optical system is imperative. However, once an optical reservoir computing system is deployed, its physical configuration becomes fixed. This fixed configuration poses a challenge for any subsequent tuning. To predict various types of time-series data, it is necessary to tune systems' parameters after the construction of the system.

Free-space optics, which expresses the coupling matrix as a light-transfer matrix, is a promising solution for optimizing the performance of RC after the construction of the system. The use of a spatial light modulator provides a flexible adjustment of the transfer matrix by controlling the wavefront. In optical reservoir computing, transmission through scattering media is used to multiply the signal by *W*res [11]. However, the controllability is limited by the SLM pixel size and pitch, which affect the representable coupling matrix.

In this chapter, we describe an optical RC approach using iterative function systems (IFSs) as a method to achieve optical tuning of the coupling matrix [12]. By employing optical affine transformation and video feedback, the coupling matrix can be flexibly tuned, allowing for the optimization of specific tasks.

#### **2 Iterative Function Systems**

For the adjustment of the coupling matrix *W*res, an optical fractal synthesizer (OFS) was employed as the tuning function. The OFS utilizes an optical computing system to generate a fractal pattern using a pseudorandom signal, which can be applied to various applications, including steam ciphers [13, 14]. Pseudo-random signals are generated based on an IFS using a collection of deterministic contraction mappings [15]. Figure 2 shows the generation of pseudorandom signals by IFS. IFS mapping comprises affine transformations of signals, including rotation, scaling, and shifting, as follows:

$$
\begin{bmatrix} \mathbf{x'}\\\mathbf{y'} \end{bmatrix} = \begin{bmatrix} \mathbf{s} \ \mathbf{0} \\\ \mathbf{0} \ \mathbf{s} \end{bmatrix} \begin{bmatrix} \cos \theta & -\sin \theta \\\ \sin \theta & \cos \theta \end{bmatrix} \begin{bmatrix} \mathbf{x} \\\ \mathbf{y} \end{bmatrix} + \begin{bmatrix} t\_x \\ t\_y \end{bmatrix}, \tag{2}
$$

where *x* and *y* are the coordinates after translation, *xandy* are those before translation, *s* is the scaling factor, θ is the rotation angle, and *tx* and *ty* are the translation parameters. The OFS can generate pseudorandom signals when *s* > 1, or fractal patterns when *s* < 1. Owing to the simultaneous processing of input images, IFS

**Fig. 2** Pattern generation by iterative function systems

provides spatio-parallel processing. The operations in Eq. 2 can be implemented by using two dove prisms and lenses, and the parameters can be tuned by adjusting the optical components. Moreover, the IFS allows for the duplication and transformation of signals in a recursive manner, enabling more intricate and complex pattern generation. This results in an increased number of degrees of freedom in the IFS operation. In the proposed system, the IFS is utilized for the operation of the coupling matrix *W*res, which can be tuned by controlling the optical setup including the rotation and tilt of the dove mirrors.

#### **3 Iterative Function System-Based Reservoir Computing**

Hardware can be implemented for tunable reservoir computing by utilizing IFS. We refer to reservoir computing using IFS as IFS reservoir computing.

Figure 3 shows the model of optical RC based on the IFS. The input signal *u*(*t*) at time *t* is multiplied by the coupling matrix *W*in and converted into a two-dimensional image. The reservoir state *X*(*t*), which is generated from the input signal, is assigned as the input of the IFS after undergoing an electric-optic conversion and multiplication quantization operation *B*. The signal is duplicated and combined after individual optical affine transformations. This signal processing is represented by the multiplication with the matrix *W*res, which can be adjusted using

**Fig. 3** Model of IFS reservoir computing [12]

the parameters of the optical affine transformation and the number of iterations. Following the addition of the input image *u*(*t*) to the transferred signal, the reservoir state *X*(*t*) is updated as

$$X(t+1) = \alpha f[B[W\_{\text{res}}B[X(t)]] + W\_{\text{in}}\mathfrak{u}(t)] + (1-\alpha)X(t), \quad 0 < \alpha < 1,\tag{3}$$

where α is a leaking rate, which determines the memory capacity of the signal in the reservoir layer. Thus, the sequence of reservoir states corresponding to a sequence of input signals is obtained. The output signal is generated by multiplying the reservoir state *X*(*t*) with a variable weight *W*out, as follows:

$$\mathbf{y}(t) = \mathbf{W}\_{\text{out}} \mathbf{X}'(t),\tag{4}$$

where *X*- (*t*) is a subset of pixels extracted from the reservoir state *X*(*t*). In reservoir computing, only the output connection weights *W*out is updated using a dataset with a sequence of input signals *u*(*t*) and the corresponding sequence of the reservoir state *X*(*t*). During the training phase, Ridge regression was adopted to optimize the output signal. The loss function *E* is expressed as

$$E = \frac{1}{n} \sum\_{t=1}^{n} \left( \mathbf{y}(t) - \hat{\mathbf{y}}(t) \right)^2 + \lambda \sum\_{i=1}^{N} \omega\_i^2,\tag{5}$$

where *n* is the number of training sets, *y***ˆ**(*t*) is the correct value, λ is the regulation parameter, and ω*i* is the *i*th element of *W*out. By minimizing the loss function *E*, *W*out is optimized.

#### **4 Prediction Performance of IFSRC**

Time-series data prediction is divided into two types: multistep and one-step-ahead predictions. The former is a task involving continuous predictions of the input signal by updating the input to the IFS reservoir computing based on the predicted value, and the latter involves predicting the input signal one-step-ahead based on the reservoir state at each time. To verify whether the proposed system can predict various types of tasks, the prediction performance was evaluated for both types of prediction tasks.

#### *4.1 Multi-step Ahead Prediction*

In the evaluation of the prediction performance for multistep ahead prediction, we employed the prediction of the Mackey-Glass equation which represents a chaotic signal and is used as a benchmark for time-series signal prediction [16]. The Mackey-Glass equation in this study is given by:

$$u(t+1) = au(t) + \frac{bu(t-\tau)}{c + u(t-\tau)^m} + 0.5,\tag{6}$$

where *a*, *b*, *c* and *m* are constants, and τ is the delay parameter. A dataset of 30,000 inputs and the next predicted values obtained from the equation were prepared, and *W*out was optimized by using Ridge regression which is the optimization method using Eq. 5. After optimization of *W*out, we assessed the system's ability to replicate a chaotic signal by inputting the predicted output into the system. The pixel size of reservoir state *X*(*t*) was set to 64 × 64 pixels.

Figure 4 shows the predicted results. The best parameters of the IFS reservoir to predict the Mackey-Glass equation were in Table 1. The IFS reservoir parameters are listed in Table 1. The inital time step in prediction phase was 300. The chaotic behavior of the Mackey-Glass equation was reproduced even for prediction phase.

**Fig. 4** Prediction of Mackey-Glass equation [12]

**Table 1** Simulation parameters in IFS reservoir computing to predict Mackey-Glass equation [12]

To evaluate the performance, the mean squared error (MSE) between the target signal and prediction output was estimated. The target signals were predicted for 261 time steps with a satisfactory MSE of <0.01. These results demonstrate the capability of IFS reservoir computing in predicting time-series data.

#### *4.2 Single-Step Ahead Prediction of Santa Fe Time-Series Data*

Single-step-ahead prediction is a task that predicts the next time signal from the input. To evaluate the performance of the system, we employed the Santa Fe time-series data, which requires memory to be predicted accurately. The Santa Fe time-series data, which models the behavior of a chaotic laser, is a widely recognized benchmark for evaluating reservoir computing systems [17]. The number of samples used for training and performing the test were 3,000, and 1,000, respectively. Figure 5a, b shows the targeted data and the prediction result. Table 2 presents the parameters of the IFS reservoir that exhibited the highest performance. Note that the best IFS parameter was different from that in case of the signal prediction of the Mackey-Glass equation. The system predicted signals similar to the label data. To evaluate the prediction performance, the normalized mean squared error was calculated between the predicted output and label. The definition of NMSE is described as follows:

$$\text{NMSE} = \frac{1}{n\sigma^2} \sum\_{t=1}^{n} (\mathbf{y}(t) - \hat{\mathbf{y}}(t))^2,\tag{7}$$

where *n* is the number of dataset, σ is the standard deviation of the inputs, *y*(*t*) is the prediction, and *y*ˆ(*t*) is the label value. As shown in Fig. 5c, the NMSE was 8.5 × 10−3. These results demonstrate that prediction performance can be improved by adjusting the IFS reservoir.

**Fig. 5** Prediction result of Santa Fe time-series **a** target data, **b** predicted data, and **c** difference between (**a**) and (**b**) [12]

**Table 2** Simulation parameters in IFS reservoir computing to predict Santa Fe time-series [12]

#### **5 Experimental Performance of IFS Reservoir Computing**

#### *5.1 Optical Setup*

To evaluate the hardware performance of the IFS reservoir computing system, an optical system featuring a video feedback system was constructed, as depicted in Fig. 6. First, the image representing the reservoir state *X*(*t*) is projected into the display (MIP3508, Prament, number of pixels: 480 × 320), and replicated by a beam splitter (BS). By using Dove prisms, optical affine transformations through rotation and tilt are processed to the individual images. The scaling factor is determined by the difference between the focal lengths of lenses L2 and L3. Individual images are

**Fig. 6** Optical setup of IFS reservoir [12]

combined by using a BS and pass through lens L4. Finally, the image is captured by the image sensor (S3–U3–123S6, FLIR, number of pixels: 4096 × 3000). The captured image is resized to perform a predetermined number of iterations and is fed back to the display. The finally-obtained image is updated to next reservoir state *X*(*t* + 1) by Eq. 3. After the processing, *X*(*t* + 1) is fed back to the display as the next IFS reservoir state. The same process is repeated, and learning is performed with pairs of reservoir states and label data. In the experiment, a region of 37 × 30 pixels in the display was sampled and used as the signals of the IFS reservoir to decrease the computational cost in Ridge regression. Moreover, a hyperbolic tangent function was used as the nonlinear function.

#### *5.2 Multi-step Ahead Prediction of Mackey-Glass Equation*

Figure 7 shows the predicted results for the Mackey-Glass equation by using the optical system. The IFS parameters used in the experiment are listed in Table 3, and each value was estimated from the obtained images. The initial status of the reservoir was set to zero, and the number of training data points was 30,000. The prediction output is the chaotic signal similar to the Mackey-Glass equation. This result shows that the IFS reservoir can perform the prediction of time-series data. However, the prediction point with MSE <0.01 was 85 steps, which is lower than the time step in numerical simulation. The reason is that the iteration parameter was fixed to 1 in the optical setup, and the captured image was resized for feedback to the display. Increasing the output signal from the reservoir layer improves the performance in the physical reservoir computing [18, 19]. However, too large size of *X*(*t*) takes computational cost in Ridge regression to optimize *W*out. Therefore it is important to

**Fig. 7** Prediction result of Mackey-Glass equation **a** target data, **b** predicted data, and **c** difference between (**a**) and (**b**) [12]

**Table 3** Experimental parameters in IFS reservoir computing [12]

adjust the resolution of an image sensor approximately depending on the time-series data to be predicted.

#### *5.3 Single-step Ahead Prediction of Santa Fe Time-Series Data*

Next, the one-step-ahead prediction of the Santa Fe time-series was evaluated. The number of data points for training and prediction was set to 3,000 and 1,000, respectively. The parameters used are listed in Table 3. Figures 8 show the label data, predicted signal, and their difference. Similar to the prediction of the Mackey-Glass equation shown in Fig. 7, the IFS reservoir computing generates a signal waveform similar to the target signal. From the difference, the NMSE is estimated as 0.033. Although the IFS parameters are not fine-tuned, the performance of the experimental IFS reservoir system is higher than that of existing physical reservoir computers [20,

**Fig. 8** Prediction result of Santa Fe time-series by optical setup. **a** Targeted data, **b** Predicted data, **c** difference between (**a**) and (**b**) [12]

21]. The results provide a promising perspective for IFS reservoir computing, which can tune the performance and flexibility in optical implementations.

#### **6 Relationship Between Performance and Spectral Radius**

To evaluate a property of reservoir computing, a spectral radius of the coupling matrix *W*res in the reservoir layer is often used [22]. The spectral radius is the largest absolute value of eigenvalues of a matrix and is defined as follows:

$$\rho(W) = \max(|\lambda\_i|, i = 1, 2, \dots, n),\tag{8}$$

where λ1, λ2,...,λ*n* are the eigenvalues of the matrix. The memory capacity increases as the spectral radius increases. In reservoir computing, a spectral radius less than one is preferred because the signal memory of the reservoir layer should be faded out [3].

In the IFS reservoir with leaking rate α, the coupling matrix was calculated as follows:

$$\mathbf{W} = \alpha \mathbf{W}\_{\text{res}} + (1 - \alpha)I,\tag{9}$$

where *I* denotes a unit matrix. To investigate the characteristics of the IFS reservoir, the spectral radius and the NMSE of one-step-ahead prediction for the Santa Fe time-series were calculated. The individual parameters were set to the values listed in

**Table 4** Combination of parameters in IFS reservoir computing [12]

**Fig. 9** Spectral radius and NMSE on one-step prediction for Santa Fe time-series [12]

Table 4, and the relationship between the spectral radius and NMSE were comprehensively verified. The size of the input image was 64 × 64, all the pixels were used for training, and the leaking rate was set to 1.0.

Figure 9 shows the relationship between the spectral radius and NMSE. Depending on the IFS parameters, the value of spectral radius is modulated, and the combination of scaling factors 0.8 and 1.0 generated a smaller spectral radius and improved the prediction performance. It was confirmed that the adjustment of IFS parameters provides modulation of coupling matrix *W*. In case of three and five iterations, the correlation coefficients were larger than 0.7, which indicates a relationship between the spectral radius and prediction performance. This result indicates that the one-stepahead prediction of the Santa Fe time-series does not require rich memory capability for the task.

#### **7 IFS Reservoir Computing with Optical Feedback Loop**

It was demonstrated that the spectral radius changed with the number of iterations, and the prediction performance changed accordingly. However, in the experiment, it was necessary to repeat the electronic feedback process to change the number of

**Fig. 10** Optical setup of IFS reservoir with an optical feedback loop

iterations. Therefore, optical feedback was introduced to realize optical control of the number of iterations. The experimental system is shown in Fig. 10. To facilitate optical feedback of the signals for multiple iterations, the combined image was transferred to BS1 through a relay lens. Consequently, a signal wherein the optical affine transformation is repeatedly executed can be generated. The combined matrix calculated in the experimental system is expressed as.

$$\mathbf{W}\_{\rm res} = (1 - \beta)\mathbf{A} + \beta(1 - \beta)\mathbf{A}^2 + \beta^2(1 - \beta)\mathbf{A}^3 + \dots \cdot \beta^{n-1}(1 - \beta)\mathbf{A}^n,\quad(10)$$

where β denotes the feedback rate of the light signal branched by the beam splitter, and *A* is a coupling matrix when the number of iterations is one. By realizing multistage iterations, the range of combined matrix values was expanded. The signal after passing through the IFS processing was detected by an image sensor. Subsequently, the reservoir state was updated based on Eq. 3, and fed back to the display. The same procedure is repeated to develop the status of the IFS reservoir.

For evaluation, one-step-ahead prediction of the Santa Fe time-series data was performed. The feedback rate β, the number of training data, and the number of prediction data were set to 0.5, 3,000, and 1,000, respectively. The optical parameters are listed in Table 3. Figure 11a–c show the prediction results for the time-series data, label data, and their differences. The NMSE value obtained under these conditions was 0.098. When the feedback loop signal was removed under the same conditions, the NMSE was 0.105, demonstrating the potential for improved prediction accuracy with feedback.

The values presented in Eq. 10 demonstrate a decrease with an increase in the order. This is owing to β < 1. The sensitivity of the image sensor was used to adjust

**Fig. 11 a** Label of Santa Fe time-series data. **b**, **d**, **f** One-step ahead prediction of the proposed system with an optical feedback loop and **c**, **e**, **g** the difference when the magnification of the gain in the image sensor was 1.0, 1.1, and 1.2, respectively

the signal acquisition of the calculated higher-order terms. Therefore, the change in the prediction ability due to the adjustment of the image sensor gain was verified.

Figure 11d, f show the prediction results when the gain was 1.1 and 1.2 times higher than that in case of Fig. 11b. In case of Fig. 11e The NMSE was 0.083 and the prediction accuracy improved. However, when the magnification of the gain was 1.2, the NMSE was 0.800, and the prediction accuracy decreased. When the gain is increased, a saturation of the light intensity occurs in the image sensor. . Consequently, a signal with effective prediction information cannot be obtained, and the prediction accuracy decreases. Therefore, it was demonstrated that it is necessary to adjust not only IFS parameters but also the sensitivity of the image sensor appropriately.

#### **8 Discussion**

The IFS reservoir computing allows for tuning of the parameters depending on the prediction task. Next step is to optimize individual parameters and maximize the performance. Various approaches have been suggested for optimizing hyperparameters in reservoir computing through computational processing [23–25]. In physical reservoir computing using FPGA, researchers have proposed methods for parameter tuning utilizing genetic algorithms [26]. Moreover, Bayesian estimation has been applied to realize more efficient parameters optimization compared to grid search methods [27]. Therefore, computational processing allows the optimization of hyperparameters efficiently, and the performance of reservoir computing can be optimized. In IFS reservoir computing, the number of IFS parameters is more than ten, which is twice as high as that in other studies. The number of parameters is corresponding to a degree of freedom in the tuning, and it is expected that higher prediction performance by the optimization can be realized. Furthermore, IFS reservoir computing provides the optimization of hyperparameters by adjusting the optical elements after construction of the system. By problem-specific parameter optimization, a reservoir computing predicting a wide range of time-series data can be built.

#### **References**

**Open Access** This chapter is licensed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license and indicate if changes were made.

The images or other third party material in this chapter are included in the chapter's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the chapter's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.

## **Bridging the Gap Between Reservoirs and Neural Networks**

**Masanori Hashimoto, Ángel López García-Arias, and Jaehoon Yu**

**Abstract** While reservoir computing is drawing attention, its applications are limited to small tasks. This chapter proposes a solution to this issue by introducing the Hidden-Fold Network, a recursive model with fixed random weights that resembles reservoir computing. The model is constructed based on recent discoveries in neural networks, namely the Strong Lottery Ticket Hypothesis and folding. By pruning an overparameterized model that is randomly initialized, it is possible to find accurate neural networks without the need for weight optimization. It is conjectured that residual networks may contain better subnetwork candidates for inference time when transformed into recurrent architectures, as they may be approximating unrolled shallow recurrent neural networks. This hypothesis is tested in image classification tasks, where subnetworks within the recurrent models are found to be more accurate and parameter-efficient than those within feedforward models, as well as the full models with learned weights.

#### **1 Introduction**

Due to their complex structure, deep neural networks are capable of various types of learning. However, when it comes to time-series problems, recurrent neural networks (RNNs) are commonly used, but as their size increases, the cost of training becomes an issue. Reservoir computing (RC) is a special type of RNNs that offer a solution to this problem. In RC, the weights of the intermediate layer are fixed, and only the output weights from the intermediate layer to the output layer are trained using low-cost learners such as linear learners, which reduces the overall training cost. Physical reservoir computing takes this a step further by implementing the intermediate layer nodes of reservoir computing as physical phenomena on hardware. This approach reduces the need for design and can achieve even greater energy efficiency

M. Hashimoto (B)

Department of Informatics, Kyoto University, Kyoto, Japan e-mail<EMAIL_ADDRESS>
Á. L. García-Arias · J. Yu Tokyo Institute of Technology, Yokohama, Japan

**Fig. 1** Approach to deep reservoir computing

by selecting appropriate physical phenomena for the intermediate layer. There are various approaches for implementing physical reservoir computing [1].

Currently, all RC models are single-layered and limited to processing small datasets like MNIST. However, although there is a desire to scale up RC to larger models capable of handling bigger datasets like CIFAR100, there are no methods to achieve this. Increasing the size of the reservoir is not feasible since it makes the output too complex for the simple classifier, and backpropagation cannot be utilized. Figure 1 illustrates our approach toward achieving deep reservoir computing (DRC). We begin by developing a method to make convolutional neural networks (CNNs) similar to RC and replacing CNN layers with reservoirs to construct a DRC.

This chapter explores the compatibility between CNN and RC using two recent concepts: the Strong Lottery Ticket Hypothesis [2, 3] and folding [4]. According to the Strong Lottery Ticket Hypothesis, high-performing neural networks can be obtained by pruning overparameterized dense models since they contain already available high-performing subnetworks. These subnetworks are sparse, random, and tiny, and can achieve competitive performance in vision tasks, making them suitable for efficient hardware implementation [5]. However, finding optimal connectivity patterns using current training methods can be challenging [6]. Folding, on the other hand, is based on the observation that a specific type of shallow RNN is equivalent to a deep residual neural network (ResNet) with weight sharing among layers [4], where the original ResNet is proposed in [7] and remains the backbone of many SOTA models. Implementing a folded RNN with significantly fewer parameters than the corresponding ResNet leads to similar performance.

Figure 2 shows our baseline idea. We fold a ResNet to then find a strong lottery ticket within the network with random fixed weights. We have found that the obtained recursive model, which is called Hidden-Fold Network (HFN), is similar to reservoir computing.

We have observed that ResNets have an inherent ability to learn ensembles of all possible unrollings of a shallow recurrent neural network. This restriction on

**Fig. 2** A sketch of hidden-fold network

the hypothesis space actually increases the number of potential subnetworks available at initialization time [8]. These resulting subnetworks are efficient in terms of parameters and demonstrate competitive performance on image classification tasks. Additionally, these subnetworks can be compressed to small memory sizes and have a high degree of parameter reusability, making them ideal candidates for inference acceleration.

#### **2 Background**

Driven by the increasingly powerful computational power promised by Moore's Law, artificial neural networks have grown in size, leading to the field of deep learning [9]. This trend was initially led by image classification models [10] and has continued as researchers enhance neural networks at the cost of larger models. Natural language processing and generative models have also joined the trend, offering impressive capabilities at the expense of immense size [11, 12]. However, while this is a convenient approach for cutting-edge research, the high computational cost of these models makes them impractical for real-world applications. Therefore, researchers are also exploring small and efficient models as an alternative trend of research.

Strong lottery tickets refer to a set of efficient neural networks obtained through a training process that combines learning, pruning, and weight quantization. The development of this approach, depicted in Fig. 3, is described in detail in Sect. 2.1. In this study, we introduce a method that converts a ResNet into a recurrent architecture to enhance the quality of the strong lottery tickets that can be extracted from it, where ResNet is introduced in Sect. 2.2.

**Fig. 3** Evolution of the training methods leading to the strong lottery ticket hypothesis

#### *2.1 Lottery Ticket Hypotheses*

Pruning is a widely used technique to compress trained neural networks into much smaller models by removing unnecessary weights [13–17]. By iteratively applying training and pruning, large portions of trained models can be removed without affecting accuracy, revealing that the original models are overparameterized. The sparsity of the resulting network can be leveraged for additional compression using entropy coding and for arithmetic optimization. This approach, combined with weight quantization, has resulted in highly efficient model compression schemes [18, 19] and specialized hardware neural accelerators [20].

In the past, it was observed that pruning did not separate the connectivity patterns from the pre-trained weights, making it impossible to reinitialize and train from scratch. However, a recent paper [21] introduced the concept of the Lottery Ticket Hypothesis (LTH), which states that overparameterized neural networks contain a subnetwork that can match the original model if trained in isolation. These subnetworks, known as winning tickets, are discovered by iteratively training, pruning, and resetting the remaining weights to their original value.

In an unexpected development during the analysis of the LTH, [2] discovered that learning weights is not essential: high-performing subnetworks exist within overparameterized neural networks at their randomly initialized state, which can be identified through pruning. Moreover, they proposed an algorithm for discovering these subnetworks by training a binary mask. Building upon this, [3] introduced a training algorithm and weight initialization method that produces sparse random subnetworks with competitive performance in image classification tasks.

Following a series of studies that explored the theoretical limits of necessary overparameterization to obtain these subnetworks [22–24], the notion of subnetworks obtained exclusively by pruning has been referred to as the *Strong Lottery Ticket Hypothesis* (SLTH). However, some researchers also refer to it as the multi-prize lottery ticket hypothesis [25], or simply as "hidden network."

Although the discovery of strong tickets is surprising, it has some related precedents. For example, extreme learning machines use fixed random weights in their hidden units, only learning the output layer. Reservoir computers also use random recurrent architectures in an analogous manner. Perturbative neural networks propose substituting convolutional layers with fixed additive random noise and a learned

linear combination. Training the batch normalization parameters of a fixed random network can achieve non-trivial accuracy. The binary neural network training method has been adapted to learn binary masks that, when applied to a trained model, extract subnetworks that perform well on untrained tasks.

#### *2.2 Residual Neural Networks*

The residual neural network (ResNet) [7], as illustrated in Fig. 4, is a popular architecture among the continuously expanding variety of neural network architectures. It serves as the backbone of many state-of-the-art (SOTA) models. ResNetis a deep convolutional neural network that follows a pyramidal feedforward structure. It comprises a convolutional pre-net, four stages of residual blocks, and a fully-connected post-net classifier. The first residual block of each stage downsamples the feature map and doubles the number of channels, adjusting the representation space size. The remaining blocks in a stage have the same size and shape and maintain the same representation space size.

Residual blocks in the ResNet architecture consist of batch normalization, ReLU activation, and convolutional layers concatenated together. In this chapter, we focus on the bottleneck block variant [7], which consists of three convolutional layers applied in sequence with kernel sizes of 1×1, 3×3, and 1×1, respectively. Each residual block has a skip-connection, which is an identity function running parallel to the block that adds the block's input to its output. In order to accommodate different representation space sizes, the downsampling blocks have a learnable layer in the skip-connection called a projection shortcut.

The original motivation for the ResNet architecture was to provide a clear path for backpropagation to reach a layer directly, solving the vanishing gradient problem. However, studies have shown that shuffling or removing residual blocks does not severely impact its performance, but rather there is a gradual degradation proportional to the amount of corruption introduced [26]. This phenomenon is not observed in other feedforward models where similar lesions result in critical performance loss. As a result, two alternative views of ResNet have been proposed: one suggests that it is an ensemble of shallow networks, while the other argues that it is an approximation of an unrolled shallow recurrent neural network. This work aims to reconcile these two views into a single coherent explanation.

Each residual block in ResNet can be thought of as a path divergence point, which leads to the interpretation of ResNet as an ensemble of all possible paths within it [26], resulting in 2*n* possible paths for a model with *n* residual blocks. Several improvements to ResNet have been proposed based on this perspective. For example, during training, removing random subsets of blocks makes the ensembled networks shallower and acts as regularization [27]. Furthermore, increasing the number of skipconnections per block increases the number of ensembled paths, thereby improving performance [28]. Another approach is to reduce network depth by increasing width, enabling the training of larger models for improved performance in less time [29].

An alternative explanation for the lesion and shuffling results is that all the blocks within a stage approximate the same function, as they have the same shape and partially receive the same inputs and gradients through the identity shortcuts. This implies that ResNet may naturally converge to the approximation of an unrolled shallow recurrent neural network, with each stage corresponding to a different hierarchical level of representation, including the downsampling block for feature map size adjustment. Meanwhile, the remaining blocks perform iterative refinement of features, according to proponents of this view [30].

The ensemble of unrollings in ResNet, which combines the two views mentioned above, the ensemble view and the unrolled iteration view, offers a vast search space to discover effective tickets within a model with limited parameters. Therefore, it is reasonable to assume that a ResNet with more recurrent approximations at its initial state may contain more potent tickets.

#### **3 Hidden-Fold Networks**

In this section, we present a technique for discovering a strong lottery ticket in a ResNet by first folding it. The approach, described in [8], yields Hidden-Fold Networks (HFNs), which outperform the strong tickets present in feedforward ResNets. Here, we provide details on the network's architecture and training. Refer to [8] for information on weight initialization and batch normalization.

#### *3.1 Folded ResNet Architecture*

In line with the unrolled iteration view, the chains of residual blocks with identical shapes in each stage approximate an iterative function. To achieve this, folding [4] is used, which converts these chains into recurrent blocks through weight sharing. In other words, *h* ≈ *g* ≈ *f* in Fig. 4 is explicitly transformed into *h* = *g* = *f* . Since applying the same functions in succession is equivalent to repeatedly applying one of them, the feedforward chain can be transformed into a single recurrent block, in a process opposite to time unrolling.

The downsampling block in ResNet has a different shape than the rest of the blocks, and therefore cannot be folded with them. Strategies to address this issue have been explored in previous works, such as removing the block to create an isotropic architecture or substituting it with a simpler block, as discussed in [4, 31]. However, this work does not modify the downsampling blocks based on the view that different stages correspond to different hierarchical levels of features, which are composed of downsampling blocks [30]. Instead, the rest of the blocks within a stage are folded into a single recurrent residual block that is iterated the same number of times as the original number of blocks, as illustrated in Fig. 5.

Folding has a dual effect on the search space. It restricts the hypothesis space of the model to iterative functions and reduces the number of parameters. Despite the exponential reduction in available subnetwork candidates, we argue that folded residual networks contain better performing strong tickets than their feedforward counterparts. If the weights of a ResNet naturally converge to approximations of unrolled iterative functions, then the strong tickets within it are likely approximations of recurrent networks. This restricts the number of relevant subnetworks to a small subset that includes consecutive blocks with similar random weights. By folding, all candidate subnetworks become recurrent, which increases the number of relevant subnetworks and their likelihood of containing stronger tickets.

Moreover, the parameter reduction occurs not only at inference time but also during training. By reducing the search space for strong tickets during training, folded tickets become easier to find. At inference time, the found subnetworks are even smaller and benefit from parameter reusability, making them ideal for efficient hardware implementation.

#### *3.2 Supermask Training*

Rather than optimizing the model's weights, the model is pruned to identify a highperforming subnetwork that is hidden within the randomly initialized folded model, which is referred to as an HFN. This connectivity pattern is learned by training a *supermask* [2], which is a pruning mask containing a binary element for each weight. During inference, the ticket is discovered by applying an element-wise product of the random weight tensor and the trained supermask.

This study adopts the edge-popup algorithm [3] for training the supermasks, as shown in Fig. 6a. The algorithm assigns a score to each weight, which is updated during backpropagation using straight-through estimation for the supermask (i.e., the supermask is not applied in the backward pass). The weights are sorted based on their scores, and the supermask is updated to include the weights with the highest top-k% scores and prune the rest. Although the value of top-k% is determined globally, the sparsity is enforced at the layer level. Folding does not impact this process; supermasks and scores are shared similarly to their corresponding weights, and backpropagation gradients are propagated through the unrolled model just like a feedforward model. Therefore, folded parameters receive distinct gradients from each iteration, as demonstrated in Fig. 6b.

#### **4 Experiments and Results**

In this section, we explore how to effectively integrate supermask training with a recurrent residual network and compare the outcomes with the baseline approaches outlined in Table 1.

**Table 1** Summary of the four methods compared on ResNet in this work

#### *4.1 Experimental Settings*

All experiments were implemented in PyTorch [32], using the original code of [3] available in their public repository [33]. The baseline model for all experiments was ResNet-50 [7] and its variations, as folding only applies to deep networks. The experiments were conducted on the CIFAR100 [34] dataset, which is relatively complex.

Unless stated otherwise, the experiments were conducted using the following methodology. We split the CIFAR100 dataset, consisting of 60*,* 000 images, into 45*,* 000 for training, 5*,* 000 for validation, and 10*,* 000 for the test set. Image preprocessing and augmentation were carried out as in [3]. We trained the models on CIFAR100 using stochastic gradient descent (SGD) with weight decay of 0*.*0005, momentum of 0*.*9, and batch size of 128 for 200 epochs. For models deeper than 100 layers or double width, we trained for an additional 100 epochs. The learning rate started at 0*.*1 and was reduced using cosine annealing with no warmup. We report the average of three runs of top-1 test accuracy scores measured at the highest scoring validation epoch. The standard deviation is shown with shaded areas on plots.

A vanilla ResNet-50 trained on CIFAR100 using an NVIDIA GeForce RTX 3090 requires 2*.*4 h. In comparison, the folded, HNN, and HFN versions of the same model require 2*.*2, 4*.*4, and 4*.*2 h, respectively.

#### *4.2 Results*

The comparison of accuracy and parameter count of various ResNet sizes trained with different methods on CIFAR100 is shown in Fig. 7a. HFN models are found to be the most parameter-efficient among the compared methods, with fewer parameters than equally performing models, and more accurate than models with similar parameter counts. Additionally, deeper and wider HFNs achieve the highest accuracies overall.

In addition, the superiority of HFNs is more pronounced when examining the model memory sizes under the compression scheme presented in [8], as shown in Fig. 7b. The memory size of HFNs is approximately half that of their feedforward counterparts, with ResNet-50 fitting into less than 2 MB. Moreover, the wider HFN models outperform the dense models that are more than 30× larger in size.

#### **5 Summary**

We explored the similarity between the Hidden-Fold Network, a recursive model with fixed random weights, and reservoir computing, as depicted in Fig. 8. In this chapter, we presented a method for folding and training ResNet with supermasks, as a first step toward deep reservoir computing. We also tested the conjecture that

**(a)** Accuracy vs. number of parameters.

**(b)** Accuracy vs. memory size.

**Fig. 7** Comparison of the four methods using different model sizes on CIFAR100. HFN is both the most parameter-efficient and the tiniest. RN and WRN are abbreviations of ResNet and Wide ResNet, respectively

**Fig. 8** The current status in this chapter and future work

recurrent residual networks have stronger tickets than their feedforward counterparts, which can be leveraged to significantly reduce the memory footprint while achieving comparable or superior accuracy to dense models. Since HFN's blocks are recurrent and unlearned, leaving all the learning load to a simple mask, the model bears strong resemblance to a deep reservoir computer. However, the second step of substituting convolutional layers with actual reservoir layers remains a future direction. Once this is achieved, we will bridge the gap between reservoirs and neural networks and approach deep reservoir computing that can solve complex tasks.

#### **References**

**Open Access** This chapter is licensed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license and indicate if changes were made.

The images or other third party material in this chapter are included in the chapter's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the chapter's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.

## **Brain-Inspired Reservoir Computing Models**

**Yuichi Katori**

**Abstract** This chapter presents an overview of brain-inspired reservoir computing models for sensory-motor information processing in the brain. These models are based on the idea that the brain processes information using a large population of interconnected neurons, where the dynamics of the system can amplify, transform, and integrate incoming signals. We discuss the reservoir predictive coding model, which uses predictive coding to explain how the brain generates expectations regarding sensory input and processes incoming signals. This model incorporates a reservoir of randomly connected neurons that can amplify and transform sensory inputs. Moreover, we describe the reservoir reinforcement learning model, which explains how the brain learns to make decisions based on rewards or punishments received after performing a certain action. This model uses a reservoir of randomly connected neurons to represent various possible actions and their associated rewards. The reservoir dynamics allow the brain to learn which actions lead to the highest reward. We then present an integrated model that combines these two reservoir computing models based on predictive coding and reinforcement learning. This model demonstrates how the brain integrates sensory information with reward signals to learn the most effective actions for a given situation. It also explains how the brain uses predictive coding to generate expectations about future sensory inputs and accordingly adjusts its actions. Overall, brain-inspired reservoir computing models provide a theoretical framework for understanding how the brain processes information and learns to make decisions. These models have the potential to revolutionize fields such as artificial intelligence and neuroscience, by advancing our understanding of the brain and inspiring new technologies.

Y. Katori (B)

The School of Systems Information Science, Future University Hakodate, 116-2 Kamedanakano-cho, Hakodate City 041-8655, Japan e-mail<EMAIL_ADDRESS>
#### **1 Introduction**

The brain's capacity to process sensory information and make decisions can be attributed to the intricate neural dynamics within a highly interconnected network of nonlinear elements known as neurons. However, the specific mechanisms underlying this framework are not yet fully understood. The primary objective of engineering and computer science researchers is to develop models that replicate the brain's information processing capabilities. A promising approach in this regard is the braininspired reservoir computing model, which has demonstrated effectiveness in diverse applications.

Reservoir computing (RC) is a framework for constructing recurrent neural networks that can model time-varying complex sensory signals [1, 2]. In the RC framework, recurrent connections are randomly and sparsely configured and do not require training. The readout connections from the reservoir are trained to reproduce a given target time series, reducing the network's computational cost. An important feature of RC is that it requires extremely low computational cost for learning because only the connections in the readout part are acquired through training. In addition, RC has many applications, including time-series generation and prediction, pattern recognition in time series, and robot control. The key requirement for RC is the presence of high-dimensional, that is, a large number of nodes or neurons that give rise to complex dynamics. Another important feature of the RC framework is the several possible physical implementations [3, 4], including electrical and optical systems, among other numerous possibilities. Provided that it has high dimensionality, nonlinearity, and echo state property, it can serve as a reservoir for computing. The framework of reservoir computing is being actively researched as an approach for modeling brain regions, such as the prefrontal cortex and cerebellum [5, 6].

Predictive coding is a theory of brain function, in which the brain processes sensory information by generating and updating internal models of the external world [7–9]. These models enable the brain to predict future sensory inputs. When the actual input deviates from the predicted input, prediction error signals are generated and sent back through the neural network to update the models. The iterative process of generating and updating predictions improves the accuracy of the model and reduces prediction error over time. Predictive coding is a widely accepted framework for understanding perception, attention, and learning in the brain and has been applied to various sensory modalities, including vision, audition, and touch. The predictive coding model is a key component of many neural network models of the brain and has been used to explain different neural phenomena, such as adaptation, attentional modulation, and perceptual illusions. Despite its success, however, the predictive coding model remains an active area of research with ongoing debates over the specific mechanisms and neural substrates underlying predictive coding in the brain [10].

The brain's reward system and the process of reinforcement learning are essential components of decision-making and learning. The reward system is a collection of neural circuits that processes information related to motivation, pleasure, and rewards. It plays a critical role in shaping behavior, such as learning and motivation, by providing feedback on the outcomes of an action. The neurotransmitter dopamine, which is released in response to a reward or the anticipation of a reward, is a key mediator of the activity of the reward system. In reinforcement learning, an agent learns to select actions based on rewards or punishments [12]. This involves learning to maximize long-term cumulative rewards by taking actions in an environment. Reinforcement learning algorithms often use trial and error to learn the optimal behavior, explore the environment, and observe the consequences of different actions. The brain's reward system and the process of reinforcement learning are closely related. The reward system provides feedback that drives the learning process, and reinforcement learning provides a framework for understanding how the brain learns to make decisions based on rewards. Computational models of reinforcement learning have been successful in explaining different behaviors, including goal-directed behavior, habit formation, and addiction [12].

In this chapter, we discuss brain-inspired reservoir computing models for sensorymotor information processing in the brain [13, 14]. First, we introduce the reservoir predictive coding model that corresponds to sensory information processing in the cerebral cortex [11]. Subsequently, we discuss the reservoir reinforcement learning model that corresponds to action learning based on rewards in the basal ganglia [14]. Finally, we present an integrated model that combines these two RC models based on predictive coding and reinforcement learning [14]. This integrated model has the potential to provide a more comprehensive understanding of the brain's information processing mechanisms.

#### **2 Reservoir-Based Predictive Coding Model**

Predictive coding is a neuroscience theory explaining how the brain processes sensory information by constantly making predictions about the world and updating them based on incoming data [7–9]. This concept is particularly relevant to the hierarchical organization of the visual system in the brain, which consists of multiple processing stages, each of which is responsible for detecting specific features of the visual input. For example, lower-level neurons may detect simple features such as edges, whereas higher-level neurons may identify more complex patterns or objects. The same principle applies to the architecture of CNNs, which have multiple layers that learn to extract increasingly complex features from input images.

Predictive coding posits that the brain actively generates predictions regarding sensory input at each level of the hierarchy. These predictions are based on information gathered from higher hierarchical levels and on previously learned internal models. These internal models, also known as generative models, make predictions and propagate them to lower levels via a top-down pathway. The difference between the actual input and the prediction, known as the prediction error, then propagates up the hierarchy in a bottom-up manner. This error signal helps the brain update its internal models and refine future predictions.

In the field of neural networks, the predictive coding with reservoir computing (PCRC) model proposed by Katori et al. [11] is a novel approach for processing time-varying sensory signals. The PCRC model employs a reservoir as the generative model for predictive coding, wherein the reservoir generates multidimensional, timevarying sensory signals. The prediction error is subsequently transmitted back to the reservoir, allowing for the rectification of the network's internal state. This model demonstrates the capability of reconstructing and predicting time-varying sensory signals.

The network architecture within each module comprises four key components: the prediction layer, input layer, prediction error layer, and reservoir (Fig. 1a). Within the module, the input signal located in the input layer is replicated in the prediction layer, which is facilitated by the complex motion of the reservoir. This prediction error is then fed back into the reservoir to minimize errors. During the training phase, the connection between the reservoir and the prediction layer is modulated using the first-order reduced and controlled error (FORCE) algorithm [15].

In the testing phase, the model operates in two distinct modes: error-driven and free-running. The error-driven mode involves feedback on the prediction error to the reservoir to further reduce the error. In contrast, the free-running mode does not involve the transmission of the prediction error to the reservoir, allowing for the autonomous operation of the reservoir. This dual-mode functionality highlights the versatility and adaptability of the PCRC model for processing time-varying sensory signals.

The PCRC module consists of a reservoir, prediction layer, input layer, and prediction error layer, which are mathematically described as follows: The membrane potential, or internal state, and the neuron activities within the reservoir are represented by *m* ∈ R*Nx* and *r* ∈ R*Nx* , respectively, where *Nx* denotes the size of the reservoir. The states of the reservoir are updated according to the following equations:

$$\begin{aligned} \mathfrak{m}(n+1) &= \mathfrak{m}(n) + \frac{1}{\tau} \{-\mathfrak{m}(n) + W\_{\text{rec}}\mathbf{r}(n) + W\_{\text{back}}\mathbf{y}(n) + \alpha\_{\varepsilon} W\_{\varepsilon} \mathbf{e}(n) - \mathfrak{b}(n)\} \,, \tag{1} \\\ r(n+1) &= \tanh(\beta\_{\mathfrak{m}} \mathfrak{m}(n+1)) \,, \tag{2} \end{aligned}$$

where *W*rec ∈ R*Nx*×*Nx* represents the matrix for recurrent connections, and τ is the time constant. The parameter β*m* scales the neuron activities. The reservoir receives inputs from the prediction layer *y*(*i*) ∈ R*Ny* through the feedback connection *W*back ∈ R*Nx*×*Ny* , the prediction error layer *e* ∈ R*Ny* with a coefficient α*e* that determines the error feedback strength and model operation mode, and the top-down input from the higher area network *b*(*i*) (*n*). The states of the prediction and prediction error layers are given by

$$\mathbf{y}(n) = \max(0, W\_{\text{out}} \mathbf{r}(n)) \; , \tag{3}$$

$$\mathfrak{e}(n) = \mathfrak{d}(n) - \mathfrak{y}(n) \,. \tag{4}$$

**Fig. 1** Network structure of the PCRC models. **a** Module of the predictive coding based on reservoir computing. **b** PCRC-based hierarchical model for the multimodal processing of the visual and auditory processing

In the error-driven mode (α*e* = 1), the reservoir is updated using the prediction error, and the state of the prediction layer follows the state of the input layer. In the free-running mode (α*e* = 0), the reservoir states are updated based on the internal dynamics, independent of the sensory input.

The network's configuration and learning process involve the following steps. The recurrent connections within the reservoir and the feedback connections from the prediction layer to the reservoir are configured in a random and sparse manner, with no need for training their connectivity. During the training phase, the network operates in error-driven mode, and the connections from the reservoir to the prediction layer are trained using the FORCE learning algorithm with a given time-series dataset. Recurrent connections *W*rec are set up using the following procedure. First, create a matrix *W*0 filled with zeros. Then, assign the non-zero values of either −1 or 1 to randomly chosen β*r* × *Nx* × *Nx* elements. Then, compute the spectral radius of *W*0: |ρ0|. Define *W*rec = α*r W*0/|ρ0|, where α*r* indicates the strength of recurrent connections. Feedback connections *W*back and *W*e are set up using the following procedure. Similar to *W*rec, generate a zero matrix *W*0 and assign non-zero values of −1 or 1 to the randomly selected β*b* × *Nx* × *Ny* elements, where β*b* specifies connectivity of the recurrent connection. Define *W*back = α*bW*0, with the strength of the feedback connections given by α*b*. Use the same procedure to generate *W*e with the coefficient α*e*.

The readout connections from the reservoir *Wout* are updated using the FORCE learning algorithm [15] as follows:

$$\mathbf{v}(n) = P(n)\mathbf{r}(n) \,, \tag{5}$$

$$P(n+1) = P(n) - \frac{\nu(n)\nu^T(n)}{1 + \nu^T(n)r(n)}\,,\tag{6}$$

$$W\_{\rm out}(n+1) = W\_{\rm out}(n) - \frac{\mathbf{e}\_{\rm }n)\mathbf{v}^{\rm T}(n)}{1 + \mathbf{v}^{\rm T}(n)r(n)}.\tag{7}$$

The initial value of *P*(*n*) ∈ R*Nx*×*Nx* is *P*(0) = *I* α *f* , where matrix *I* is an identity matrix and α *f* is a scaling parameter. Once the readout connection training is complete, the module can reconstruct the given input in error-driven mode and predict the input in free-running mode.

The PCRC-based hierarchical model for multimodal processing comprises three modules (Fig. 2b). Each module in the hierarchical model is distinguished by superscript (*i*), where *i* ∈ {*V*, *A*, *I*} denotes the visual, auditory, and integration modules, respectively.

The configuration and learning of this model were performed using the following steps: The recurrent and feedback connections were established in accordance with the previously described procedure. The connection matrices between the lower and higher levels, *U*(*A*) and *U*(*V*), are defined using the method below; their inverse matrices are used for dimensionality reduction.

Firstly, operate the lower area network (visual and auditory areas) in error-driven mode (α(*V*) *e* = 1 and α(*A*) *e* = 1) without top-down signals (α*td* = 0), and gather the time course of the reservoirs *r*(*V*) and *r*(*A*) in the state collecting matrices *R*(*V*) and *R*(*A*) , respectively. Next, compute the dimension reduction matrices *U*(*V*) and *U*(*A*). Assuming that *T* timesteps of reservoir states are collected in *R*(*i*) (*i* ∈ {*V*, *A*}), *R*(*i*) can be decomposed by principal component analysis (PCA) as *R*(*i*) = *S*(*i*) *UT* (*i*). Here, *S*(*i*) is a *T* × 20 matrix, and *U*(*i*) is an *N*(*i*) *x* × 20 matrix. The dimension reduction matrix *U* −1 (*i*) can be obtained as the pseudo-inverse matrix of *U*(*i*). Finally, connect the sensory modules (visual and auditory modules) and the integration using the obtained *U*(*i*), and operate the entire network in error-driven mode (α(*V*) *e* = α(*A*) *e* = 1) with α*td* > 0. The matrices *W*(*V*) out , *W*(*A*) out , and *W*(*I*) out are acquired using FORCE learning.

Within the hierarchical PCRC model for visual and auditory processing, the integration reservoir is responsible for reconstructing and predicting the compressed and concatenated states of the sensory reservoirs. Consequently, the integration reservoir is expected to reconstruct information from one modality using information from another modality. Both the visual and auditory reservoirs are driven by the prediction error on each sensory module and the integration reservoir.

The multimodal model is assessed using time-series data, consisting of pairs of hand-written digit images and their corresponding spoken number utterances. Three hand-written digit images ("2," "5," and "9") from the MNIST dataset are employed as visual signals [16]. Each image comprises 28 × 28 (784) grayscale pixels. These images undergo preprocessing via non-negative matrix factorization (NMF) and are converted into a 20-dimensional signal. Assuming *V* is an *L* × 784 matrix with each row representing an individual image, and *V* is a collection of *L* images, NMF decomposes *V* into two matrices: *V* = *H W*, where *H* is an *L* × 20 coefficient matrix and *W* is a 20 × 784 feature matrix. The transformed 20-dimensional vector serves as the input to the visual area network. The coefficient vector reconstructed by the PCRC module can be converted back into images using *W*.

In addition, linguistic data containing spoken number utterances from the Ti46 dataset are utilized as auditory signals [17]. This dataset comprises uncompressed audio data. Each dataset is preprocessed using a cochlear filter model [18], an auditory model that simulates sound propagation within the inner ear, and the conversion of acoustic energy into neural representations. The auditory signals are transformed into 55-dimensional signals. Figure 2 displays samples of the dataset. In the auditory signal, the initiation of spoken number utterances exhibits jitter, starting anywhere from 60 to 90 timesteps. The corresponding visual signals are presented from 80 to 160 timesteps without jitter.

After training, the network is expected to reconstruct sensory information from one modality based on input signals originating from the other modality. In the subsequent analysis, the focus is on reconstructing visual information in the presence of corresponding auditory signals. In this case, the auditory and integration reservoirs operate in error-driven mode, whereas the visual reservoir functions in free-running mode.

In the association process, a given auditory signal is initially presented to the input layer of the auditory area. At this time, the reservoir maintains a silent state; as no signal is formed in the prediction layer, a significant prediction error occurs. This prediction error serves as a trigger, inducing the motion of the auditory reservoir and generating the auditory signal in the prediction layer. Subsequently, the prediction error gradually decreases. A spatial pattern reflecting the temporal pattern in the auditory signal is represented within the reservoir. This information is spatially compressed and conveyed to the integration reservoir, where only the auditory information is input. The prediction layer in the integration area initially remains

**Fig. 2** Evaluation of the PCRC-based multimodal model with auditory and visual signals. **a** Training phase: All modules are operated in error-driven mode. The model is trained with the datasets comprising time-series pairs of visual signals (hand-written digits) and auditory signals (spoken utterances of corresponding digits). Each signal is displayed for 100 timesteps. Visual signals undergo preprocessing using NMF, resulting in a 20-dimensional signal that serves as sensory input for the visual area network. Auditory signals are preprocessed with a cochlear filter and converted into 55-dimensional signals. These signals are then provided as sensory input to the auditory area network. **b** Cross-modal association from the auditory signal to the visual signal after the training phase: Auditory and the integration modules are operated in error-driven mode, whereas the auditory module is operated in free-running mode. The model is driven by an auditory signal, and the corresponding visual image appears in the prediction layer of the visual area

silent, and the prediction error triggers the activity of the integration reservoir. As the integration reservoir begins to move, predictions are generated to compensate for the prediction error. At this time, both the auditory and corresponding visual signals are generated. Because there is no signal coming from the lower visual layer, the prediction error regarding the visual information is larger. This prediction error is then transmitted to the visual area, inducing activity in the visual reservoir. Based on the fluctuations in the visual reservoir, visual signal prediction is performed. In summary, a visual signal corresponding to the input auditory signal is generated in the prediction layer of the visual area.

During the processing of multidimensional complex time courses, the proposed hierarchical model combines the mechanisms of temporal structure accumulation and spatial pattern compression. The input signal is reconstructed using a reservoir that

**(a)**

captures the temporal structure of the signal within its high-dimensional nonlinear dynamics. Subsequently, the high-dimensional state vector in the reservoir, which encompasses a short history of the signal, undergoes spatial compression and is transferred to the integration area network. This combination of accumulation and compression results in a higher-order abstraction of the intricate time course. In cross-modal association, the processes of compression and abstraction are reversed, allowing the generation of sensory information through expansion and instantiation.

#### **3 Reservoir-Based Reinforcement Learning Model**

Reinforcement learning is a type of learning in which an agent learns to choose actions based on rewards or punishments, with the goal of maximizing long-term cumulative rewards by taking actions in an environment. Reinforcement learning algorithms often employ trial and error to learn optimal behavior, explore the environment, and observe the consequences of various actions. The learning process is fueled by feedback from the environment, providing information on the outcomes of an agent's actions. Among the various reinforcement learning approaches, TD-learning is a model-free technique that merges ideas from dynamic programming and Monte Carlo methods [12, 19]. It estimates the value function (expected future reward) by learning from the difference between consecutive predictions, which is known as the temporal difference error. This method enables agents to learn online and update their value estimates incrementally as new experiences are acquired, making it particularly suitable for learning in dynamic environments.

In recent years, RL has been combined with deep learning to create deep reinforcement learning (DRL), which has achieved remarkable success in solving complex control tasks with high-dimensional sensory inputs, such as images and sounds [20]. DRL algorithms, such as Deep Q-Networks (DQN) [20], proximal policy optimization (PPO) [21], and actor-critic methods [22, 23], have been successfully applied to various applications, such as video game playing, robotics, and autonomous driving.

The two important frameworks within RL are Markov decision processes (MDPs) [24] and partially observableMarkov decision processes (POMDPs) [25].MDPs are a mathematical framework used to model decision-making problems in reinforcement learning, where the environment's state transitions and rewards are assumed to be Markovian; that is, the future state depends only on the current state and action taken and not on previous states or actions. Although MDPs have been successfully applied to various problems, they exhibit certain limitations, particularly in partially observable environments.

In real-world situations, an agent may not have full access to the environment's state owing to noisy sensors, occlusions, or other factors. This lack of complete information about the environment's state can lead to suboptimal decision-making, as the agent cannot accurately estimate the value of different actions. This is where POMDPs play a significant role, extending the MDP framework to handle environments with partial observability.

POMDPs are a generalization of MDPs that consider uncertainty in perceiving the environment's state. Instead of using the environment's true state, the agent maintains a belief state, which is a probability distribution over the possible environment states. The belief state is updated as the agent takes actions and receives observations, allowing it to make better-informed decisions even with incomplete information. However, solving POMDPs is generally more computationally demanding than MDPs owing to the increased complexity associated with maintaining and updating belief states.

One approach to address POMDP is the use of reservoirs. Reservoirs, which store not only the current state but also the history of sensory inputs reflecting the environmental state in high-dimensional state vectors, can be expected to function effectively in POMDP environments by compensating for information that cannot be directly observed. In the following, we introduce how the reservoir reinforcement learning model, which is a model that reads action values from the reservoir states where the history of sensory information is accumulated, can effectively function in POMDP environments.

The proposed model consists of a sensory layer, a dynamical reservoir, and an output layer (Fig. 3a). The reservoir receives sensory input from the environment through the sensory layer and generates action values in the output layer, which are then converted to action commands. The reservoir state, comprising *Nx* neurons, is denoted by *x*(*t*) ∈ R*Nx* . The dynamical reservoir state *x*(*t*) evolves as follows:

$$\mathbf{x}(t+\Delta t) = \mathbf{x}(t) + \frac{\Delta t}{\tau\_x} \left( -\mathbf{x}(t) + W^{\text{in}} \mathbf{u}(t+\Delta t) + W^{\text{rec}} \mathbf{r}(t) + W^{\text{back}} \mathbf{q}(t) \right), \tag{8}$$

where τ*x* is the time constant; *W*in is the *Nx* × *Nu* sensory matrix from the sensory layer to the dynamical reservoir, *W*rec is the *Nx* × *Nx* recurrent weight matrix in the dynamical reservoir, and *W*back is the *Nx* × *Ny* feedback matrix from the output layer to the dynamical reservoir. These weight matrices *W*in, *W*rec, and *W*back are randomly and sparsely generated and remain fixed. The neuron firing rate in the dynamical reservoir,*r*(*t*), is defined as*r*(*t*) = *fr*(β*x*(*t*)), where β specifies the firing rate responses and *fr*(*x*) = tanh(*x*). The output layer state, denoted by *q*(*t*) ∈ R*Ny* , represents the *Ny* action values and is specified according to:

$$\boldsymbol{q}(t) = \boldsymbol{W}^{\text{out}} \boldsymbol{r}(t) \,. \tag{9}$$