---
abstract: 'The ability of discrete-time nonlinear recurrent neural networks to store time-varying small input signals is investigated by mean-field theory. The combination of a small input strength and mean-field assumptions makes it possible to derive an approximate expression for the conditional probability density of the state of a neuron given a past input signal. From this conditional probability density, we can analytically calculate short-term memory measures, such as memory capacity, mutual information, and Fisher information, and determine the relationships among these measures, which have not been clarified to date to the best of our knowledge. We show that the network contribution of these short-term memory measures peaks before the edge of chaos, where the dynamics of input-driven networks is stable but corresponding systems without input signals are unstable.'
author:
- Taichi Haruna
- Kohei Nakajima
title: 'Optimal short-term memory before the edge of chaos in driven random recurrent networks'
---

Introduction {#sec:intro}
============

Natural and artificial high-dimensional nonlinear dynamical systems can be used as resources for real-time computing. By nonlinearly mapping time-varying input signals into a high-dimensional space, the signals can be learned in a supervised manner if the dynamical systems have enough ability to store the signals in their present state and separate different signals [@Jaeger2001; @Maass2002]. A high computational performance can be achieved by tuning only the weights of linear connections to the output layer while keeping the parameters of the dynamical systems fixed [@Jaeger2004; @Verstraeten2007; @Lukosevicius2009; @Pathak2018]. Such dynamical systems called reservoirs can be artificial recurrent neural networks (RNNs) or physical systems, such as optical media [@Larger2012; @Larger2017], nanoscale magnetization dynamics [@Torrejon2017; @Tsunegi2019], soft materials [@Nakajima2018], and quantum systems [@Fujii2017].

As mentioned above, a requirement for real-time computing is the ability to memorize past input signals. Such short-term memory of dynamical systems has been studied extensively by assessing a quantity called memory capacity [@Jaeger2002; @Dambre2012]. For input-driven RNNs, it has been suggested that the part of memory capacity representing indirect memory through network takes a maximum value near the edge of chaos, namely, near the critical boundary between the stable and unstable dynamical regimes [@Bertschinger2004; @Boedecker2012]. Near criticality, different inputs are expected to lead to different states while suppressing the influence of the initial conditions. Hence, it seems reasonable for a dynamical system to be near the critical point for optimal memory capacity. However, it has also been pointed out that the dependence on network parameters is not straightforward based on a systematic numerical simulation [@Farkas2016].

For linear RNNs, detailed analytic studies of memory capacity can be performed for both discrete-time [@White2004; @Rodan2011] and continuous-time systems [@Hermans2010]. The ability to predict future inputs, which is complementary to memory capacity, has also been studied in linear systems with correlated input signals [@Marzen2017]. However, the memory capacity of nonlinear RNNs is difficult to study by analytical methods [@Ganguli2008]. Recently, Schuecker et al. [@Schuecker2018] successfully derived an analytical expression for memory capacity for continuous-time nonlinear RNNs [@Sompolinsky1988] in which each neuron is driven by independent input signals following a white-noise Gaussian process. Toyoizumi and Abbott [@Toyoizumi2011] analytically calculated the signal-to-noise ratio, which is equivalent to the inverse of memory capacity at the limit of zero input strength, for discrete-time nonlinear RNNs driven by a common time-varying input signal.

In this paper, we analytically investigate the memory capacity of discrete-time nonlinear RNNs called echo state networks (ESNs) [@Jaeger2001] by a mean-field theory when the strength of input signals is small but non-zero. The main idea of our approach is that the conditional probability density of the present state of a neuron given a past input signal can be approximately calculated from a functional derivative with respect to past input signals under the assumption of a small input strength. Once we obtain this conditional probability density, it is straightforward to derive the memory capacity and other alternative memory measures, such as mutual information and Fisher information [@Ganguli2008]. We show that all three measures of short-term memory through network behave similarly and take a maximum value before the edge of chaos, where the dynamics is stable in the presence of input signals but unstable in the absence of input signals. We also discuss the breakdown of the mean-field theory for calculating memory measures in the ordered regime and show that the linear approximation provides good predictions.

Results {#sec:results}
=======

Echo State Networks {#sec:esn}
-------------------

We consider ESNs consisting of $N$ artificial neurons. The state of neuron $i$ at discrete time step $t$ is denoted $x_i(t)$. We assume that the time evolution of state $x_i(t)$ is governed by $$x_i(t+1)=f(a_i(t)), 
\label{eq:1}$$ where $f$ is an activation function. $a_i(t)$ is the activation potential of neuron $i$ at time step $t$ given by $$a_i(t)=\sum_{j=1}^N w_{ij} x_j(t) + u_i s(t), 
\label{eq:2}$$ where $s(t)$ is a time-dependent input signal, $w_{ij}$ is a time-independent weight of the connection from neuron $j$ to neuron $i$, and $u_i$ is a time-independent weight representing the strength of the coupling from the input signal to neuron $i$. We use the matrix and vector notations $W:=\left( w_{ij} \right)_{1 \leq i,j \leq N}$, $\bm{u}:=(u_1,u_2,\dots,u_N)^\textrm{T}$, and $\bm{x}(t):=(x_1(t),x_2(t),\dots,x_N(t))^\textrm{T}$.

In the following analytical calculations and numerical simulations, the activation function is assumed to be a sigmoid function satisfying $\lim_{a \to \pm \infty} f(a)=\pm 1$, $f(-a)=-a$ and $f'(0)=1$. In particular, we adopt $f(a)=\textrm{erf}(\frac{\sqrt{\pi}}{2}a)$ owing to its analytical tractability. $w_{ij}$ are chosen independently at random from an identical Gaussian distribution with mean zero and variance $g^2/N$, where $g^2>0$ is a control parameter. For simplicity, $u_i$ are assumed to be independent variables taking $\pm 1$ with probability $\frac{1}{2}$. Since our primary concern is the memory capacity of ESNs, we consider an independent and identically distributed Gaussian input signal $s(t)$ with mean zero and variance $s^2$. All numerical results in this paper were obtained in the following way unless otherwise stated. We simulated ESNs with $N=1000$ artificial neurons over $40$ trials. For a single trial, each quantity (for example, stationary variance of $x_i(t)$) was calculated from its values over $10^5$ time steps after discarding the initial $10^4$ time steps. Then, averages were obtained over all artificial neurons and all trials. All infinite sums appearing in the following sections were evaluated by truncation at the $500$-th term.

![ Numerical results (marks) and mean-field predictions (solid lines) for the stationary variance of $x_i(t)$ are shown as functions of $g^2$ for $s^2=0.01$ (red), $s^2=0.02$ (green), and $s^2=0.04$ (blue). Vertical broken lines indicate the mean-field predictions of the critical value of $g^2$ obtained from Eq. (\[eq:4\]) for $s^2=0$ ($g^2=1$, black), $s^2=0.01$ ($g^2 \approx 1.39$, red), $s^2=0.02$ ($g^2 \approx 1.50$, green), and $s^2=0.04$ ($g^2 \approx 1.64$, blue). []{data-label="fig:1"}](fig1.eps){width="8cm"}

The mean-field theory of ESNs [@Cessac2007; @Massar2013; @Molgedey1992; @Toyoizumi2011] makes it possible to calculate the stationary variance of $x_i(t)$ and the largest Lyapunov exponent in the limit $N \to \infty$. It assumes that $x_i(t)$ are independent and identically distributed random variables. They are also assumed to be independent of $w_{ij}$ and $u_i$. This assumption can be justified in the limit $N \to \infty$ when there is no input signal. Since $f$ is odd, we can self-consistently assume that the mean of $x_i(t)$ is equal to zero. Let $\sigma^2(t)=\langle x_i(t)^2 \rangle$ be the variance of $x_i(t)$, where $\langle \cdots \rangle$ indicates the average over trials with the same $w_{ij}$ and $u_i$ but possibly different realizations of the input signal $s(t)$ and initial conditions. By the central limit theorem, $a_i(t)$ follows a Gaussian distribution with mean zero and variance $g^2 \sigma^2(t)+s^2+ O(N^{-\frac{1}{2}})$, where we use $u_i^2=1$. Neglecting the $O(N^{-\frac{1}{2}})$ term, the variance of $a_i(t)$ does not depend on specific realizations of $w_{ij}$ [@Toyoizumi2011]. In the following, we omit quantities that approach $0$ as $N \to \infty$ unless otherwise stated. Consequently, $\sigma^2(t)$ follows the following recurrence equation [@Massar2013]: $$\begin{aligned}
\sigma^2(t+1)
&= \int_{-\infty}^\infty da \, f^2(a) \frac{\exp\left( - \frac{a^2}{2\Sigma^2(t)} \right)}{\sqrt{2\pi \Sigma^2(t)}} \nonumber\\
&= -1 + \frac{4}{\pi} \arctan \left( \sqrt{1+\pi \Sigma^2(t)} \right),
\label{eq:3}\end{aligned}$$ where $\Sigma^2(t)=g^2 \sigma^2(t)+s^2$. By numerically solving Eq. (\[eq:3\]) with $\sigma^2(t+1)=\sigma^2(t)=\sigma^2$, we can obtain the mean-field prediction of the stationary variance of $x_i(t)$. We write $\Sigma^2=g^2 \sigma^2+s^2$. These values of $\sigma^2$ and $\Sigma^2$ have been used for plotting theoretical results later in the paper.

The largest Lyapunov exponent derived from the mean-field theory is [@Massar2013; @Molgedey1992; @Toyoizumi2011] $$\begin{aligned}
\lambda
&= \frac{1}{2} \log \left( g^2 \int_{-\infty}^\infty da \, f'^2(a) \frac{\exp\left( - \frac{a^2}{2\Sigma^2} \right)}{\sqrt{2\pi \Sigma^2}} \right) \nonumber\\
&= \frac{1}{2} \log \frac{g^2}{\sqrt{1+\pi \Sigma^2}}. 
\label{eq:4}\end{aligned}$$ Here, $\lambda>0$ indicates that a small perturbation to a state of the system leads to exponential growth, while $\lambda<0$ implies that the perturbation eventually becomes undetectable. The dynamics is called chaotic or unstable in the former case and called ordered or stable in the latter. When an input signal is absent, the boundary between chaos and stability $\lambda=0$ corresponds to $g^2=1$. The presence of input signals shifts the boundary towards the chaotic side [@Massar2013; @Molgedey1992].

In Fig. \[fig:1\], we confirm that the mean-field prediction of the stationary variance of $x_i(t)$ agrees well with the result obtained by numerical simulations. We note that the difference is negligible, even for the input-driven regime $g^2 \ll 1$ where the mean-field assumption is expected to be violated. This will be explained when we discuss the breakdown of the mean-field theory in the calculation of memory capacity.

Memory Capacity {#sec:mc}
---------------

The memory capacity of an ESN is defined as the quality of the optimal linear estimator of the past input $s(t-n)$ using the present state of neurons $x_i(t)$. Following previous work [@Schuecker2018; @Toyoizumi2011], we assume a sparse readout, namely, there are $K=O(1)$ readout neurons $1 \leq i \leq K$ and consider linear readout $\hat{s}(t)=\sum_{i=1}^K v_i x_i(t)$. Given time-delay $n$ ($n=1,2,\dots$), the weights $v_i$ are determined by minimizing the mean squared error between $s(t-n)$ and $\hat{s}(t)$ over a sufficiently long time period $T$. The optimal mean squared error as a function of time-delay $n$ is called the memory function and is given by [@Dambre2012; @Jaeger2002] $$M_n = \frac{\sum_{i,j=1}^K \langle s(t-n)x_i(t) \rangle_T d_{ij} \langle x_j(t)s(t-n) \rangle_T }{\langle s(t-n)^2 \rangle_T},
\label{eq:5}$$ where $d_{ij}$ is the $(i,j)$-th element of the matrix $C^{-1}$, which is the inverse of the matrix $C=\left( \langle x_i(t)x_j(t) \rangle_T \right)_{1 \leq i,j \leq K}$, and $\langle \cdots \rangle_T$ indicates the time average over the period of length $T$. The memory capacity is defined as the sum of $M_n$ over all time-delays $n$: $$M=\sum_{n=1}^\infty M_n. 
\label{eq:6}$$ It is known that $M \leq K$ holds [@Jaeger2002; @Dambre2012].

To compute $M_n$ by the mean-field theory, we replace the time average $\langle \cdots \rangle_T$ for $T \to \infty$ by the average over trials $\langle \cdots \rangle$ in stationary states. Since $\langle x_i(t)x_j(t) \rangle$ for $i \neq j$ vanishes as $N \to \infty$, the contribution of the off-diagonal terms of $C^{-1}$ to Eq. (\[eq:5\]) can be neglected for $N \to \infty$ by the sparse readout assumption $K=O(1)$. Thus, in the mean-field calculation, Eq. (\[eq:5\]) is just $K$ times $M_n$ for $K=1$. Since $$M_n = \frac{\langle x_i(t)s(t-n) \rangle^2}{\sigma^2 s^2} 
\label{eq:7}$$ when $K=1$, the task to obtain $M$ reduces to calculating $\langle x_i(t)s(t-n) \rangle$ for $n=1,2,\dots$. In the following, we perform the calculation by assuming that the strength of the input signal is small, namely, $s^2 \ll 1$.

![image](fig2.eps){width="16cm"}

The main idea to calculate $\langle x_i(t)s(t-n) \rangle$ is that conditioning of $a_i(t)$ on $s(t-n)$ ($n=0,1,2,\dots$) can be regarded as a small perturbation to dependence of $a_i(t)$ on $s(t-n)$, when $s^2 \ll 1$. We assume that $x_i(t)$ are independent and identically distributed and are also independent of $w_{ij}$ and $u_i$ even after conditioning. By the central limit theorem, $a_i(t)$ conditioned on $s(t-n)=c$ follows a Gaussian distribution in the limit of large $N$. Hence, it is sufficient to calculate its mean $\mu_{n,c,i}$ and variance $\Sigma^2_{n,i}$ to determine the conditional probability density $P_{a_i(t)|s(t-n)}(a|c)$.

First, we calculate the mean of $a_i(t)$ given $s(t-n)=c$. We regard $a_i(t)$ as a functional of stochastic variables $s(t),s(t-1),\dots,s(t-n),\bm{x}(t-n)$. We write $a_i(t)=\mathcal{F}_i \left[\bm{s}_{0:n}(t), \bm{x}(t-n) \right]$, where $\bm{s}_{0:n}(t)=(s(t),s(t-1),\dots,s(t-n))$. We consider a norm defined by the average over trials $\langle \cdots \rangle$ ($\|X\|:=\sqrt{\langle X^2 \rangle}$ for a stochastic variable $X$). Conditioning of $a_i(t)$ on $s(t-n)=c$ corresponds to replacing argument $s(t-n)$ with the constant stochastic variable $c$. If $c^2,s^2 \ll 1$, then $\| c - s(t-n)\|^2=\langle c - s(t-n) \rangle^2=c^2+s^2 \ll 1$. Thus, $a_i(t)$ given $s(t-n)=c$ can be approximated by the following:

$$\begin{aligned}
\mathcal{F}_i \left[\bm{s}_{0:n-1}(t),c,\bm{x}(t-n) \right]
&\simeq \mathcal{F}_i \left[\bm{s}_{0:n}(t),\bm{x}(t-n) \right] + \frac{\delta \mathcal{F}_i \left[\bm{s}_{0:n}(t),\bm{x}(t-n) \right] }{\delta s(t-n)} \left( c - s(t-n) \right) \nonumber\\
&= a_i(t) +  \frac{\delta a_i(t)}{\delta s(t-n)} \left( c - s(t-n) \right)
\label{eq:8}\end{aligned}$$

By taking the average over trials, we have $$\mu_{n,c,i}= \left\langle \frac{\delta a_i(t)}{\delta s(t-n)} \right\rangle c - \left\langle \frac{\delta a_i(t)}{\delta s(t-n)}s(t-n) \right\rangle. 
\label{eq:9}$$ By applying the mean-field assumption, we find (Appendix \[sec:a\]) $$\left\langle \frac{\delta a_i(t)}{\delta s(t-n)} \right\rangle = \left( V^n \bm{u} \right)_i, 
\label{eq:10}$$ where $$V=\frac{1}{\sqrt{1+\frac{\pi}{2}\Sigma^2}}W.
\label{eq:11}$$ Since $f$ is an odd function, we have (Appendix \[sec:b\]) $$\left\langle \frac{\delta a_i(t)}{\delta s(t-n)}s(t-n) \right\rangle = 0. 
\label{eq:12}$$ From Eqs. (\[eq:10\]) and (\[eq:12\]), the mean-field theory predicts $$\mu_{n,c,i}=\left( V^n \bm{u} \right)_i c 
\label{eq:13}$$ when $c^2, s^2 \ll 1$.

Second, the variance of $a_i(t)$ given $s(t-n)=c$ can be obtained as follows. The variance of $a_i(t)$ can be expressed as $$\begin{aligned}
\Sigma^2 
%&= \int_{-\infty}^{\infty} da P_{a_i(t)}(a) a^2 \nonumber\\
&= \int_{-\infty}^{\infty} da \, \int_{-\infty}^{\infty} dc \, P_{a_i(t) | s(t-n)}(a|c) P_{s(t-n)}(c) a^2 \nonumber\\
&= \Sigma_{n,i}^2 + \left( V^n \bm{u} \right)_i^2 s^2. 
\label{eq:14}\end{aligned}$$ Thus, we have $$\Sigma_{n,i}^2 = \Sigma^2 - \left( V^n \bm{u} \right)_i^2 s^2. 
\label{eq:15}$$ Note that the population variance of $\left( V^n \bm{u} \right)_i^2$ takes a nonzero finite value even in the limit of large $N$, as we will see in the linear case (Appendix \[sec:e\]) when we discuss the breakdown of the mean-field theory in the ordered regime. This implies that the value of $\Sigma_{n,i}^2$ depends on $i$ or, equivalently, realizations of $W$ and $\bm{u}$, even after discarding the $O(N^{-\frac{1}{2}})$ term. Another related remark is that Eq. (\[eq:15\]) holds only in the limit of small $s^2$. Otherwise, the right-hand side may become negative even in the limit of large $N$, since $\left( V^n \bm{u} \right)_i$ follows a Gaussian distribution with a variance of $O(1)$ and thus can take an arbitrarily large value.

In Fig. \[fig:2\], the mean of $a_1(t)$ given $s(t-n)=0.1$ and the conditional probability density $P_{a_1(t)|s(t-9)}(a|c)$ are shown for single specific realizations of $W$ and $\bm{u}$. Here, we set $s^2=0.01$. The numerical results are obtained by first generating a single orbit of length $10^6$ time steps after discarding the initial $10^4$ time steps and then sampling the value of $a_1(t)$ with $0.1 \leq s(t-n) <0.11$ for each $n$. The theoretical values for $\mu_{n,0.1,1}$ and $\Sigma_{n,1}^2$ are calculated from Eqs. (\[eq:11\]), (\[eq:13\]), and (\[eq:15\]), where $W$ and $\bm{u}$ are the same as those used in the numerical simulation. We can see that the numerical results and the theoretical predictions agree well.

![image](fig3.eps){width="16cm"}

![image](fig4.eps){width="16cm"}

Using Eqs. (\[eq:13\]) and (\[eq:15\]), we obtain (Appendix \[sec:c\]) $$\begin{aligned}
%&\langle x_i(t)s(t-n) \rangle \nonumber\\
%&= \int_{-\infty}^{\infty} da \, \int_{-\infty}^{\infty} dc \, P_{a_i(t-1) | s(t-n)}(a|c) P_{s(t-n)}(c) f(a)c \nonumber\\
%&= \frac{\left( V^{n-1} \bm{u} \right)_i s^2}{\sqrt{1+\frac{\pi}{2}\Sigma^2}} 
\langle x_i(t)s(t-n) \rangle = \frac{\left( V^{n-1} \bm{u} \right)_i s^2}{\sqrt{1+\frac{\pi}{2}\Sigma^2}}
\label{eq:16}\end{aligned}$$ for $n=1,2,\dots$. From Eqs. (\[eq:7\]) and (\[eq:16\]), we have $$M_n = \frac{\left( V^{n-1} \bm{u} \right)_i^2 s^2}{\sigma^2 \left( 1+\frac{\pi}{2}\Sigma^2 \right)}
\label{eq:17}$$ for $n=1,2,\dots$. The population average of $M_n$, which is equivalent to the average over realizations of $W$ and $\bm{u}$, is (Appendix \[sec:d\]) $$\mathbf{E}\left[ M_n \right] = \frac{r^{n-1} s^2}{\sigma^2 \left( 1+\frac{\pi}{2}\Sigma^2 \right)} = \frac{r^n s^2}{g^2 \sigma^2}, 
\label{eq:18}$$ where $$r=\frac{g^2}{1+\frac{\pi}{2}\Sigma^2}. 
\label{eq:19}$$ The population average of $M$ is $$\mathbf{E}\left[ M \right]=\sum_{n=1}^\infty \mathbf{E}\left[ M_n \right] = \frac{r s^2}{g^2 \sigma^2 (1 - r)}. 
\label{eq:20}$$ $M$ can be decomposed into two parts [@Farkas2016; @Schuecker2018]: the direct memory $M_1$ and the indirect memory through network $M_\textrm{net}:=M-M_1$. We call the latter the *network memory capacity*. The population average of the latter is $$\mathbf{E}\left[ M_\mathrm{net} \right] = \frac{r^2 s^2}{g^2 \sigma^2 (1 - r)} = r \mathbf{E}\left[ M \right]. 
\label{eq:21}$$

Figure \[fig:3\] (a) shows the population average of $M_n$ for different values of $g^2$ with $s^2=0.01$. The population averages of $M$ and $M_\mathrm{net}$ are shown in Fig. \[fig:3\] (b) and (c), respectively. $\mathbf{E}\left[ M_\mathrm{net} \right]$ peaks in the range $1 < g^2 < g_*^2$, where $g_*^2$ is the value of $g^2$ such that $\lambda=0$. The exact location of the maximum point depends on the value of $s^2$ and shifts to a larger value as $s^2$ increases. In the mean-field theory, $\mathbf{E}\left[ M_\mathrm{net} \right]$ is given as the product between $\mathbf{E}\left[ M \right]$ and $r$. Hence, its qualitative behavior can be understood from those of $\mathbf{E}\left[ M \right]$ and $r$ (Fig. \[fig:3\] (b) and (d), respectively). Since $\mathbf{E}\left[ M \right]$ is a measure of the linear short-term memory, it is expected to decrease as the nonlinearity of the system increases. On the other hand, $r$ can be interpreted as an effective measure of the nonlinear response of the system, which reaches saturation for sufficiently large $g^2$, since the activation function $f$ is a sigmoid function. Indeed, $r \to \frac{2}{\pi}$ as $g^2 \to \infty$, since $\sigma^2 \to 1$ as $g^2 \to \infty$ in Eq. (\[eq:19\]) (However, this cannot be seen from Fig. \[fig:3\] (d) because the range of $g^2$ shown is restricted upto $2$).

The mean-field predictions and the numerical results agree well over the whole range of $g^2$ for $\mathbf{E}\left[ M \right]$. However, there is a clear discrepancy for $\mathbf{E}\left[ M_\mathrm{net} \right]$ in the ordered regime (Fig. \[fig:3\](c)). This is due to the breakdown of the mean-field theory. That is, the assumption that $x_i(t)$ are independent and identically distributed and are also independent of $w_{ij}$ and $u_i$ is violated when the ESN dynamics is driven by input signals. Indeed, in a certain range of $g^2$ in the ordered regime ($0.2 < g^2 < 0.7$ in Fig. \[fig:4\] (a) where $s^2=0.01$), the numerically obtained values of $\mathbf{E}\left[ M_\mathrm{net} \right]$ do not approach the mean-field value as the system size $N$ increases. To understand the quantitative influence of the violation of the mean-field assumption on $\mathbf{E}\left[ M_\mathrm{net} \right]$, we consider the regime $g^2 \ll 1$, where the activation function $f$ can be approximated by the identity function $f(x)=x$. When $g^2 \ll 1$, we can approximately calculate $\mathbf{E}\left[ M_n \right]$ without the mean-field theory. Since both the mean-field theory and the linear approximation lead to $\mathbf{E}\left[ M \right]=1$ for $g^2 \ll 1$, the difference in $\mathbf{E}\left[ M_\mathrm{net} \right]$ is reduced to that in $\mathbf{E}\left[ M_1 \right]$. The linear approximation predicts (Appendix \[sec:e\]) $$\mathbf{E}\left[ M_1 \right] 
= s^2 \mathbf{E}\left[ \frac{1}{\sigma_i^2} \right] 
\simeq 1-g^2 + 2\frac{\left(1-g^2\right)^2}{1+g^2}g^4, 
\label{eq:22}$$ where $\sigma_i^2$ is the stationary variance of $x_i(t)$. Note that we have $\mathbf{E}\left[ \sigma_i^2 \right]=\frac{s^2}{1-g^2}$ in both the mean-field theory and the linear approximation. Indeed, in the mean-field theory, Eq. (\[eq:3\]) reduces to $\sigma^2=\Sigma^2=g^2\sigma^2+s^2$ when $f(x)=x$. The equation for the linear approximation is derived in Appendix \[sec:e\] (Eq. (\[eq:a14\])). However, the former predicts $\mathbf{E}\left[ \frac{1}{\sigma_i^2} \right]=\frac{1}{\sigma^2}+O(N^{-\frac{1}{2}})$. When calculating $\mathbf{E}\left[ M_1 \right]$ for $g^2 \ll 1$, the mean-field theory fails to capture the variance of $\sigma_i^2$.

We compare the values of $\mathbf{E}\left[ M_1 \right]$ obtained from the numerical simulation, the linear approximation, and the mean-field theory in Fig. \[fig:4\] (b). Although the mean-field line does not fit the numerical result, the linear approximation can explain it well.

Mutual Information and Fisher Memory {#sec:mi}
------------------------------------

![image](fig5.eps){width="16cm"}

![ Population average of $J_\mathrm{ot}$ (Eq. (\[eq:28\])) calculated based on the mean-field theory under the same assumption as for Eq. (\[eq:24\]). To adjust the scale of the vertical axis, $J_\mathrm{ot}$ is multiplied by $s^2$. []{data-label="fig:6"}](fig6.eps){width="8cm"}

Once we obtain the conditional probability density $P_{a_i(t)|s(t-n)}(a|c)$ for $n=0,1,2,\dots$, we can immediately calculate the mutual information between $x_i(t+1)$ and $s(t-n)$ as $$I(x_i(t+1);s(t-n))=I(a_i(t);s(t-n)) \simeq \frac{1}{2} \log \frac{\Sigma^2}{\Sigma_{n,i}^2} 
\label{eq:23}$$ when the mean-field assumption is valid. We would like to take the population average of Eq. (\[eq:23\]). Recall that we assumed $s^2 \ll 1$. Let us suppose $\Sigma^2=O(1)$ in the limit $s^2 \to 0$. In particular, this holds when $g^2>1$. Then, we can approximate the population average of $I(x_i(t+1);s(t-n))$ as $$\begin{aligned}
&\mathbf{E}\left[ I(x_i(t+1);s(t-n)) \right] \nonumber\\
&\simeq \frac{1}{2} \log \frac{\Sigma^2}{\mathbf{E}\left[ \Sigma_{n,i}^2 \right]} = \frac{1}{2} \log \frac{\Sigma^2}{\Sigma^2 - r^n s^2}. 
\label{eq:24}\end{aligned}$$

Let us consider the summation of Eq. (\[eq:23\]) over $n=1,2,\dots$ defined by $$I_\mathrm{ot} = \sum_{n=1}^\infty I(x_i(t+1);s(t-n)), 
\label{eq:25}$$ where the subscript $\mathrm{ot}$ indicates the mutual information between the future state and a one-time past input. Note that the $n=0$ term is not included in the summation. Thus, $I_\mathrm{ot}$ is a measure of network short-term memory analogous to $M_\mathrm{net}$ based on the mutual information. The population average $\mathbf{E}\left[ I_\mathrm{ot} \right]$ calculated based on the mean-field theory is shown in Fig. \[fig:5\](a) and is compared with the numerical results. Since the direct numerical estimate of the mutual information between $x_i(t+1)$ and $s(t-n)$ for all $n=1,2,\dots,500$ is computationally hard, we estimated the mutual information from the correlation coefficient between $a_i(t)$ and $s(t-n)$ assuming that $P_{a_i(t)|s(t-n)}(a|c)$ is Gaussian, which is valid both in the linear regime and the mean-field regime. As in the case of $\mathbf{E}\left[ M_\mathrm{net} \right]$, $\mathbf{E}\left[ I_\mathrm{ot} \right]$ also takes a maximum value in the range $1 < g^2 < g_*^2$.

As we have seen in the calculation of memory capacity, the mean-field theory is not applicable to the linear regime $g^2 \ll 1$. Indeed, although the discrepancy between the numerical results and the mean-field predictions appears to be small on the scale of Fig. \[fig:5\](a), the calculation of $I_\mathrm{ot}$ based on the linear approximation (Appendix \[sec:e\]) provides much better fits to the numerical results than the mean-field theory for $g^2 \ll 1$, as shown in Fig. \[fig:5\](b).

![image](fig7.eps){width="16cm"}

![The mean-field values of $g^2$ ($g_\chi^2$, $g_\iota^2$, $g_\phi^2$, and $g_\omega^2$) that attain the maxima of the five network memory measures ($\mathbf{E}\left[ M_\mathrm{net} \right]$, $\mathbf{E}\left[ I_\mathrm{ot} \right]$, $\mathbf{E}\left[ J_\mathrm{ot} \right]$, and $\mathbf{E}\left[ I \right]$ ($\mathbf{E}\left[ J \right]$), respectively) are shown as functions of $s^2$. Note that $\mathbf{E}\left[ I \right]$ and $\mathbf{E}\left[ J \right]$ peak at the same value of $g^2$. The critical line $g^2=g_*^2$ is also shown. []{data-label="fig:8"}](fig8.eps){width="8cm"}

![$I=\lim_{n \to \infty}I(x_i(t+1);\bm{s}_{1:n}(t))$ can be represented as the difference between $I_1=I(x_i(t+1);\bm{x}(t))$ and $I_2=\lim_{n \to \infty}I(x_i(t+1);\bm{x}(t) | \bm{s}_{1:n}(t))$ (Eq. (\[eq:36\])). The population-averaged values of these quantities calculated from the mean-field theory under the same assumption as for Eq. (\[eq:24\]) are shown for $s^2=0.01$. []{data-label="fig:9"}](fig9.eps){width="8cm"}

Another familiar information-theoretic memory measure is the Fisher information [@Cover2006; @Ganguli2008]. Here, we regard the past input $s(t-n)$ as a parameter and consider the Fisher information with respect to the conditional probability density $P_{x_i(t+1)|s(t-n)}(x|c)$, namely, information about $s(t-n)$ contained in $x_i(t+1)$. Since the activation function $f$ is invertible and the Fisher information is invariant under an invertible transformation of stochastic variables, we can use $P_{a_i(t)|s(t-n)}(a|c)$ to calculate the Fisher information. When $s^2 \ll 1$ and the mean-field assumption is valid, the Fisher information for $s(t-n)$ contained in $x_i(t+1)$ is calculated as $$\begin{aligned}
&J(n) \nonumber\\
&= \int_{-\infty}^\infty da \, \left( \frac{\partial}{\partial c} \log P_{a_i(t)|s(t-n)}(a|c) \right)^2 P_{a_i(t)|s(t-n)}(a|c) \nonumber\\
&\simeq \frac{\left( \frac{\partial \mu_{n,c,i}}{\partial c} \right)^2}{\Sigma_{n,i}^2} = \frac{\left( V^n \bm{u} \right)_i^2}{\Sigma_{n,i}^2}. 
\label{eq:26}\end{aligned}$$ The population average of Eq. (\[eq:26\]) can be approximately obtained under the same assumption as for Eq. (\[eq:24\]) as follows: $$\mathbf{E}\left[ J(n) \right] \simeq \frac{\mathbf{E}\left[ \left( V^n \bm{u} \right)_i^2 \right] }{\mathbf{E}\left[ \Sigma_{n,i}^2 \right]} = \frac{r^n}{\Sigma^2-r^n s^2}. 
\label{eq:27}$$ We define the network Fisher memory with respect to a one-time past input by $$J_\mathrm{ot} = \sum_{n=1}^\infty J(n). 
\label{eq:28}$$ As in the case of $I_\mathrm{ot}$, we exclude the $n=0$ term representing direct memory from the sum in the right-hand side of Eq. (\[eq:28\]). The population average $\mathbf{E}\left[ J_\mathrm{ot} \right]$ calculated based on the mean-field theory under the same assumption as for Eq. (\[eq:24\]) is shown in Fig. \[fig:6\]. $\mathbf{E}\left[ J_\mathrm{ot} \right]$ behaves qualitatively similarly to $\mathbf{E}\left[ M_\mathrm{net} \right]$ and $\mathbf{E}\left[ I_\mathrm{ot} \right]$ at least in the range where the mean-field theory is valid. Note that there is a close relationship between the mean-field predictions of memory function, mutual information, and Fisher information through Eqs. (\[eq:18\]), (\[eq:24\]) and (\[eq:27\]): $\mathbf{E}\left[ I(x_i(t+1);s(t-n)) \right] = \frac{1}{2}\log \left( 1 - \left( 1-\frac{s^2}{\Sigma^2} \right) \mathbf{E}\left[ M_n \right] \right) = \frac{1}{2} \log \left( 1+\mathbf{E}\left[ J(n) \right]s^2 \right)$.

The derivation of the conditional probability density $P_{a_i(t)|s(t-n)}(a|c)$ under the mean-field assumption can be extended in a straightforward manner to the conditioning on a set of past inputs at multiple time steps. In particular, the conditional probability density $P_{a_i(t)|\bm{s}_{1:n}(t)}(a|\bm{c})$ can be approximated as a Gaussian distribution with mean $\mu_{1:n,\bm{c},i}$ and variance $\Sigma_{1:n,i}^2$ where $\bm{s}_{1:n}(t)=(s(t-1),s(t-2),\dots,s(t-n))$, $\bm{c}=(c_1,c_2,\dots,c_n)$, $$\mu_{1:n,\bm{c},i} = \sum_{k=1}^n \left(V^k \bm{u} \right)_i c_k
\label{eq:29}$$ and $$\Sigma_{1:n,i}^2 = \Sigma^2 - \sum_{k=1}^n \left(V^k \bm{u} \right)_i^2 s^2. 
\label{eq:30}$$

An alternative network memory measure to $I_\mathrm{ot}$ is the limit $n \to \infty$ of the mutual information between $x_i(t+1)$ and $\bm{s}_{1:n}(t)$. When the mean-field assumption is valid, it is given by $$\begin{aligned}
I =& \lim_{n \to \infty}I(x_i(t+1);\bm{s}_{1:n}(t)) \nonumber\\
=& \lim_{n \to \infty}I(a_i(t);\bm{s}_{1:n}(t)) \simeq \frac{1}{2} \log \frac{\Sigma^2}{\Sigma_{1:\infty,i}^2}, 
\label{eq:31}\end{aligned}$$ where $\Sigma_{1:\infty,i}^2=\Sigma^2 - \sum_{k=1}^\infty \left(V^k \bm{u} \right)_i^2 s^2$. Under the same assumption as for Eq. (\[eq:24\]), its population average is approximately given by $$\mathbf{E}\left[ I \right] \simeq \frac{1}{2} \log \frac{(1-r)\Sigma^2}{(1-r)\Sigma^2 - r s^2}. 
\label{eq:32}$$

Similarly, we can consider an alternative network memory measure based on the Fisher information matrix with respect to $P_{a_i(t)|\bm{s}_{1:n}(t)}(a|\bm{c})$. When the mean-field assumption is valid, the $(k,l)$-th element of the Fisher information matrix is given by $$J_{kl} \simeq \frac{\frac{\partial \mu_{1:n,\bm{c},i}}{\partial c_k}\frac{\partial \mu_{1:n,\bm{c},i}}{\partial c_l}}{\Sigma_{1:n,i}^2} = \frac{\left(V^k \bm{u} \right)_i \left(V^l \bm{u} \right)_i}{\Sigma_{1:n,i}^2}. 
\label{eq:33}$$ The Fisher memory with respect to the whole past input history is defined by $$J = \lim_{n \to \infty}\sum_{k=1}^n J_{kk}
\label{eq:34}$$ and its approximate population average is found to be $$\mathbf{E}\left[ J \right] \simeq \frac{r}{(1-r)\Sigma^2 - r s^2} 
\label{eq:35}$$ under the same assumption as for Eq. (\[eq:24\]). Note that $\mathbf{E}\left[ I \right]$ and $\mathbf{E}\left[ J \right]$ are related by $\mathbf{E}\left[ I \right] \simeq \frac{1}{2}\log \left( 1 + \mathbf{E}\left[ J \right] s^2 \right)$. Figure \[fig:7\] shows Eqs. (\[eq:32\]) and (\[eq:35\]) for $s^2=0.01, 0.02$, and $0.04$. Both $\mathbf{E}\left[ I \right]$ and $\mathbf{E}\left[ J \right]$ take maximum values at points close to those for $\mathbf{E}\left[ M_\mathrm{net} \right]$, $\mathbf{E}\left[ I_\mathrm{ot} \right]$, and $\mathbf{E}\left[ J_\mathrm{ot} \right]$ as long as the mean-field theory is valid. Figure \[fig:8\] summarizes the maximum points of these network memory measures in the range $0<s^2\leq 0.04$ obtained from the mean-field theory.

Finally, we remark on the behavior of $I$. Because $x_i(t+1)$ and $\bm{x}(t)$ are conditionally independent given $\bm{s}_{1:n}(t)$, we obtain $$\begin{aligned}
&I(x_i(t+1);\bm{x}(t)) \nonumber\\
&= I(x_i(t+1);\bm{s}_{1:n}(t)) + I(x_i(t+1);\bm{x}(t) | \bm{s}_{1:n}(t)). 
\label{eq:36}\end{aligned}$$ The first term on the right-hand side of Eq. (\[eq:36\]) becomes $I$ by taking the limit $n \to \infty$. The second term $I(x_i(t+1);\bm{x}(t) | \bm{s}_{1:n}(t))$ will be negligible when the system is driven by input signals ($g^2 \ll 1$). On the other hand, the chaotic dynamics dominates as $g^2 \to \infty$ and $I(x_i(t+1);\bm{x}(t) | \bm{s}_{1:n}(t))$ will approach $I(x_i(t+1);\bm{x}(t))$. Thus, as $g^2$ varies from $0$ to $\infty$, $I$ will increase together with $I(x_i(t+1);\bm{x}(t))$ in the ordered regime but will decrease in the sufficiently chaotic regime. Hence, $I$ is expected to take a maximum value between the two extremes. In Fig. \[fig:9\], the population-averaged values of the three terms in Eq. (\[eq:36\]) in the limit $n \to \infty$ calculated from the mean-field theory under the same assumption as for Eq. (\[eq:24\]) are shown, where $I(x_i(t+1);\bm{x}(t))=\frac{1}{2} \log \frac{\Sigma^2}{s^2}$ due to the mean-field assumption.

Discussion {#sec:discussion}
==========

The three network memory measures studied in this paper take maximum values in the ordered regime for ESNs with small input signals. The value of $g^2$ that attains the maximum is always greater than $1$, which is the boundary between the ordered and chaotic regimes in the corresponding autonomous system. However, it is far from the critical $g_*^2$ (Fig. \[fig:8\]). In previous work, it was argued that the maximal Fisher information can be used to detect the edge of chaos [@Prokopenko2011; @Livi2017]. Our results suggest that such an approach is not necessarily effective for driven dynamical systems.

In the context of physical reservoir computing [@Larger2012; @Larger2017; @Torrejon2017; @Tsunegi2019; @Nakajima2018; @Fujii2017], it is generally difficult to tune the parameters of a given physical system for optimal computational performance. An alternative method is to choose an optimal input strength. Although the analysis presented in this paper assumes that $s^2$ is small, our results theoretically suggest that tuning the input strength is meaningful. For example, the value of $\mathbf{E}\left[ M_\mathrm{net} \right]$ for $s^2=0.02$ is greater than those for $s^2=0.01$ and $s^2=0.04$ around $g^2 \approx 1.3$ in Fig. \[fig:3\] (c).

Toyoizumi and Abbott [@Toyoizumi2011] analytically showed that the signal-to-noise ratio of ESNs decreases rapidly on the left side of the criticality $g^2=1$ when inputs are absent, but decreases much slowly on the right side. They suggested that high computational performance can be achieved without fine tuning in the latter. Our results confirm this expectation because all of the short-term memory measures in the presence of inputs peak near $g^2=1$ in the region $1 < g^2 < g_*^2$.

In general, dynamical regimes of autonomous systems beyond stable fixed points are candidates for computational resources. For example, RNNs with sinusoidal activation functions achieve a high computational performance in the non-chaotic window regions after transition to chaos occurs in their autonomous dynamics [@Marquez2018]. An online supervised learning algorithm for RNNs proposed by Sussillo and Abbott [@Sussillo2009] exhibits its best performance when their autonomous dynamics is adjusted to the chaotic region not far from the critical point where chaotic dynamics can be suppressed by input signals. Schuecker et al. [@Schuecker2018] showed that the network memory capacity for continuous-time nonlinear RNNs peaks in the ordered regime with $g^2>1$, which is consistent to our result. They argued that the dynamic suppression of chaos (DSC), which results from the fact that the onset of local instability precedes that of asymptotic instability, contributes to optimal information processing. However, DSC cannot occur in discrete-time ESNs where the two onsets coincide. In ESNs, the shift of the critical $g_*^2$ toward the chaotic regime induced by input signals is solely due to a mechanism called the static suppression of chaos (SSC), which increases the frequency with which an orbit visits the contracting region of the phase space. Unlike SSC, DSC is conjectured to occur based on fast switching among different unstable directions caused by input signals [@Schuecker2018]. However, ESNs with leaky neurons [@Jaeger2007] are expected to exhibit DSC [@Schuecker2018]. Future analyses of network memory measures for leaky ESNs based on the presented theory could deepen the understanding of the relationship between DSC and the information processing ability of dynamical systems.

It has been suggested that there exists a trade-off between nonlinearity of dynamical systems and their memory capacity [@Dambre2012]. Inubushi and Yoshimura [@Inubushi2017] theoretically investigated the trade-off in terms of how nonlinearity degrades small initial differences of input signals. The mean-field theory presented in this paper makes it possible to study the trade-off when input strength is small by directly calculating the nonlinear memory capacity proposed by Damble et al. [@Dambre2012]. Performing the detailed calculation is also left as future work.

Derivation of Eq. (\[eq:10\]) {#sec:a}
=============================

We have $$\begin{aligned}
&\left\langle \frac{\delta a_i(t)}{\delta s(t-n)} \right\rangle \nonumber\\
&= \sum_{j_1,\dots,j_n} w_{i j_1} \dots w_{j_{n-1}j_n}u_{j_n}\left\langle \prod_{k=1}^n f'(a_{j_k}(t-k)) \right\rangle. 
\label{eq:a0}\end{aligned}$$ Since $a_{j_k}(t-k)$ $(k=1,2,\dots,n)$ are independent and $a_{j_k}(t-k) \sim N(0,\Sigma^2)$ by the mean-field assumption, we have $$\begin{aligned}
\left\langle \prod_{k=1}^n f'(a_{j_k}(t-k)) \right\rangle 
&= \prod_{k=1}^n  \left\langle f'(a_{j_k}(t-k)) \right\rangle \nonumber\\
&= \left( \int_{-\infty}^\infty da \, f'(a) \frac{\exp\left( - \frac{a^2}{2\Sigma^2} \right)}{\sqrt{2\pi \Sigma^2}} \right)^n \nonumber\\
&= \left( \frac{1}{\sqrt{1+\frac{\pi}{2}\Sigma^2}} \right)^n, 
\label{eq:a1}\end{aligned}$$ where $f'(a)=e^{-\frac{\pi}{4}a^2}$ is used for $f(a)=\textrm{erf}(\frac{\sqrt{\pi}}{2}a)$. Equation (\[eq:10\]) follows from Eqs. (\[eq:11\]), (\[eq:a0\]) and (\[eq:a1\]).

Derivation of Eq. (\[eq:12\]) {#sec:b}
=============================

We have $$\begin{aligned}
&\left\langle \frac{\delta a_i(t)}{\delta s(t-n)} s(t-n) \right\rangle \nonumber\\
&= \sum_{j_1,\dots,j_n} w_{i j_1} \dots w_{j_{n-1}j_n}u_{j_n} \left\langle \left( \prod_{k=1}^n f'(a_{j_k}(t-k)) \right) s(t-n) \right\rangle. 
\label{eq:a2}\end{aligned}$$ It is sufficient to show $$\left\langle \left( \prod_{k=1}^n f'(a_{j_k}(t-k)) \right) s(t-n) \right\rangle = 0.
\label{eq:a3}$$ We set $$g_n(y,z_n)=y+u_{j_n}z_n
\label{eq:a4}$$ and

$$\begin{aligned}
g_{n-k}(y,z_{n-k},\dots,z_n) 
= \sum_{j_{n-k+1}} w_{j_{n-k}j_{n-k+1}} f(g_{n-k+1}(y,z_{n-k+1},\dots,z_n)) + u_{j_{n-k}}z_{n-k}
\label{eq:a5}\end{aligned}$$

for $k=1,\dots,n-1$. Let us introduce $$G(z_1,\dots,z_n)=z_n \int_{-\infty}^\infty Dy \, \prod_{k=1}^n f'(g_k(y,z_k,\dots,z_n)),
\label{eq:a6}$$ where $Dy=\frac{dy}{\sqrt{2\pi g^2 \sigma^2}}e^{-\frac{y^2}{2 g^2 \sigma^2}}$. The left-hand side of Eq. (\[eq:a3\]) can be written as $$\begin{aligned}
&\left\langle \left( \prod_{k=1}^n f'(a_{j_k}(t-k)) \right) s(t-n) \right\rangle \nonumber\\
&= \int_{-\infty}^\infty Dz_1 \dots \int_{-\infty}^\infty Dz_n \, G(z_1,\dots,z_n),
\label{eq:a7}\end{aligned}$$ where $Dz_k=\frac{dz_k}{\sqrt{2\pi s^2}}e^{-\frac{z_k^2}{2 s^2}}$. We shall show that $G$ is an odd function with respect to $(z_1,\dots,z_n)$, namely, $$G(-z_1,\dots,-z_n)=-G(z_1,\dots,z_n) 
\label{eq:a8}$$ holds. This yields the desired result. First, note that $g_{n-k}(y,z_{n-k},\dots,z_n)$ is odd with respect to $(y,z_{n-k},\dots,z_n)$ for $k=0,1,\dots,n-1$. Namely, we have $$g_{n-k}(-y,-z_{n-k},\dots,-z_n) = -g_{n-k}(y,z_{n-k},\dots,z_n). 
\label{eq:a9}$$ Indeed, Eq. (\[eq:a9\]) can be proved by mathematical induction. First, $g_n(-y,-z_n)=-y-u_{j_n}z_n=-g_n(y,z_n)$ for $k=0$. Assume that $g_{n-k}(-y,-z_{n-k},\dots,-z_n)=-g_{n-k}(y,z_{n-k},\dots,z_n)$ for $k=0,1,\dots,n-2$. Then, we obtain

$$\begin{aligned}
g_{n-(k+1)}(-y,-z_{n-(k+1)},\dots,-z_n) &= \sum_{j_{n-k}} w_{j_{n-(k+1)}j_{n-k}} f(g_{n-k}(-y,-z_{n-k},\dots,-z_n)) - u_{j_{n-(k+1)}}z_{n-(k+1)} \nonumber\\
&= \sum_{j_{n-k}} w_{j_{n-(k+1)}j_{n-k}} f(-g_{n-k}(y,z_{n-k},\dots,z_n)) - u_{j_{n-(k+1)}}z_{n-(k+1)} \nonumber\\
&= - \sum_{j_{n-k}} w_{j_{n-(k+1)}j_{n-k}} f(g_{n-k}(y,z_{n-k},\dots,z_n)) - u_{j_{n-(k+1)}}z_{n-(k+1)} \nonumber\\
&= -g_{n-(k+1)}(y,z_{n-(k+1)},\dots,z_n), 
\label{eq:a10}\end{aligned}$$

where we applied the induction hypothesis for the third equality and we used the fact that $f$ is an odd function for the fourth equality. Now, Eq. (\[eq:a8\]) is obtained by $$\begin{aligned}
&G(-z_1,\dots,-z_n) \nonumber\\
&= -z_n \int_{-\infty}^\infty Dy \, \prod_{k=1}^n f'(g_k(y,-z_k,\dots,-z_n)) \nonumber\\
&= -z_n \int_{-\infty}^\infty Dy \, \prod_{k=1}^n f'(g_k(-y,-z_k,\dots,-z_n)) \nonumber\\
&= -z_n \int_{-\infty}^\infty Dy \, \prod_{k=1}^n f'(-g_k(y,z_k,\dots,z_n)) \nonumber\\
&= -z_n \int_{-\infty}^\infty Dy \, \prod_{k=1}^n f'(g_k(y,z_k,\dots,z_n)) \nonumber\\
&= -G(z_1,\dots,z_n),
\label{eq:a11}\end{aligned}$$ where we used Eq. (\[eq:a9\]) for the third equality and the fact that $f'$ is an even function for the fourth equality.

Derivation of Eq. (\[eq:16\]) {#sec:c}
=============================

We have $$\begin{aligned}
&\langle x_i(t)s(t-n) \rangle \nonumber\\
&= \int_{-\infty}^{\infty} da \, \int_{-\infty}^{\infty} dc \, P_{a_i(t-1) | s(t-n)}(a|c) P_{s(t-n)}(c) f(a)c \nonumber\\
&= \frac{\left( V^{n-1} \bm{u} \right)_i s^2}{\sqrt{2\pi}\Sigma^3} \int_{-\infty}^\infty da \, e^{-\frac{a^2}{2\Sigma^2}} f(a) a. 
\label{eq:a20}\end{aligned}$$ By the change of variable $y=\frac{a}{\Sigma}$, the integral on the right-hand side of Eq. (\[eq:a20\]) can be calculated as $$\begin{aligned}
\int_{-\infty}^\infty da \, e^{-\frac{a^2}{2\Sigma^2}} f(a) a 
&= \Sigma^2 \int_{-\infty}^\infty dy \, e^{-\frac{y^2}{2}} f(\Sigma y) y \nonumber\\
&= \Sigma^2 \int_{-\infty}^\infty dy \, \frac{d}{dy}\left( - e^{-\frac{y^2}{2}}\right)f(\Sigma y) \nonumber\\
&= \Sigma^3 \int_{-\infty}^\infty dy \, e^{-\frac{\left( 1+\frac{\pi}{2}\Sigma^2 \right)y^2}{2}} \nonumber\\
&= \frac{\sqrt{2\pi}\Sigma^3}{\sqrt{1+\frac{\pi}{2}\Sigma^2}}. 
\label{eq:a21}\end{aligned}$$ From Eqs. (\[eq:a20\]) and (\[eq:a21\]), we obtain Eq. (\[eq:16\]).

Derivation of Eq. (\[eq:18\]) {#sec:d}
=============================

Since $\left( V^n \bm{u} \right)_i^2 = \left( \frac{1}{1+\frac{\pi}{2}\Sigma^2} \right)^n \left( W^n \bm{u} \right)_i^2$, it is sufficient to show $\mathbf{E} \left[ \left( W^n \bm{u} \right)_i^2 \right] = g^{2n} + O(\frac{1}{N})$ for $n=0,1,\dots$.

Let $w_{ij}^{(n)}$ be the $(i,j)$-th element of $W^n$. We have $$\begin{aligned}
\mathbf{E} \left[ \left( W^n \bm{u} \right)_i^2 \right]
&= \sum_{j,k} \mathbf{E} \left[ w_{ij}^{(n)} w_{ik}^{(n)} \right] \mathbf{E} \left[ u_j u_k \right] \nonumber\\
&= \sum_j \mathbf{E} \left[ \left( w_{ij}^{(n)} \right)^2 \right], 
\label{eq:a22}\end{aligned}$$ where we used $\mathbf{E} \left[ u_j u_k \right]=\delta_{jk}$ and $\delta_{jk}$ is the Kronecker delta. The population average of $\left( w_{ij}^{(n)} \right)^2$ is given by $$\mathbf{E} \left[ \left( w_{ij}^{(n)} \right)^2 \right] 
= \frac{g^{2n}}{N} + O\left(\frac{1}{N^2}\right), 
\label{eq:a23}$$ because $w_{ij} \sim N(0,\frac{g^2}{N})$ are independent. Thus, we obtain $$\mathbf{E} \left[ \left( W^n \bm{u} \right)_i^2 \right]=g^{2n} + O\left(\frac{1}{N}\right). 
\label{eq:a24}$$

Linear Approximation {#sec:e}
====================

When $g^2 \ll 1$, we can approximate $f(a)=a$ and obtain $$\bm{x}(t+1) = \sum_{n=0}^\infty s(t-n)W^n \bm{u}. 
\label{eq:a12}$$ Thus, $$\sigma_i^2 = \langle x_i(t+1)^2 \rangle = s^2 \sum_{n=0}^\infty \left( W^n \bm{u} \right)_i^2
\label{eq:a13}$$ holds. Ignoring the $O(\frac{1}{N})$ terms, the mean and variance of $\left( W^n \bm{u} \right)_i^2$ for $n \geq 1$ are given by $g^{2n}$ and $2g^{4n}$, respectively, because $w_{jk}$ and $u_j$ are independent and $w_{jk} \sim N(0,\frac{g^2}{N})$ and $u_j=\pm 1$. The population average of $\sigma_i^2$ is given by $$\mathbf{E} \left[ \sigma_i^2 \right] = s^2 \sum_{n=0}^\infty g^{2n} = \frac{s^2}{1-g^2}. 
\label{eq:a14}$$ Its variance is $$\begin{aligned}
\mathrm{Var} \left[ \sigma_i^2 \right] 
&= s^4 \left( \mathbf{E} \left[ \left( \sum_{n=0}^\infty \left( W^n \bm{u} \right)_i^2 \right)^2 \right] - \mathbf{E} \left[ \sum_{n=0}^\infty \left( W^n \bm{u} \right)_i^2 \right]^2 \right) \nonumber\\
&= s^4 \sum_{n=0}^\infty \mathrm{Var} \left[ \left( W^n \bm{u} \right)_i^2 \right] + O\left(\frac{1}{N}\right) \nonumber\\
&= 2s^4 \sum_{n=1}^\infty g^{4n} = \frac{2 g^4 s^4}{1-g^4}, 
\label{eq:a15}\end{aligned}$$ where we used $\mathbf{E}\left[ \left( W^m \bm{u} \right)_i^2 \left( W^n \bm{u} \right)_i^2 \right] = \mathbf{E}\left[ \left( W^m \bm{u} \right)_i^2 \right] \mathbf{E}\left[ \left( W^n \bm{u} \right)_i^2 \right] + O(\frac{1}{N})$ for $m \neq n$. From Eqs. (\[eq:a14\]) and (\[eq:a15\]), we obtain $$\begin{aligned}
\mathbf{E} \left[ \frac{1}{\sigma_i^2} \right] 
&\simeq \frac{1}{\mathbf{E} \left[ \sigma_i^2 \right]} \left( 1 + \frac{\mathrm{Var} \left[ \sigma_i^2 \right]}{\mathbf{E} \left[ \sigma_i^2 \right]^2} \right) \nonumber\\
&= \frac{1-g^2}{s^2} \left( 1 + 2 \frac{1-g^2}{1+g^2}g^4 \right). 
\label{eq:a16}\end{aligned}$$ Equation (\[eq:22\]) follows from Eq. (\[eq:a16\]).

Similarly, we can compute the mutual information between $x_i(t+1)$ and $s(t-n)$ in the linear regime $g^2 \ll 1$ as follows. Let $X:= \sum_{m=0}^\infty \left( W^m \bm{u} \right)_i^2$ and $X_n:=\sum_{m=0,m \neq n}^\infty \left( W^m \bm{u} \right)_i^2$. We have $$\mathbf{E} \left[ I(x_i(t+1);s(t-n))\right] = \frac{1}{2} \mathbf{E} \left[ \log \frac{X}{X_n} \right]. 
\label{eq:a17}$$ Since $\mathbf{E} \left[ \log X \right]$ can be approximated as $$\mathbf{E} \left[ \log X \right] \simeq \log \mathbf{E} \left[ X \right] - \frac{1}{2}\frac{\mathrm{Var}\left[ X \right]}{\mathbf{E} \left[ X \right]^2} 
\label{eq:a18}$$ and a similar approximation can be obtained for $\mathbf{E} \left[ \log X_n \right]$, we obtain $$\begin{aligned}
&\mathbf{E} \left[ I(x_i(t+1);s(t-n))\right] \nonumber\\
&\simeq \frac{1}{2} \log\frac{\mathbf{E} \left[ X \right]}{\mathbf{E} \left[ X_n \right]} + \frac{1}{4} \left( \frac{\mathrm{Var}\left[ X_n \right]}{\mathbf{E} \left[ X_n \right]^2} - \frac{\mathrm{Var}\left[ X \right]}{\mathbf{E} \left[ X \right]^2} \right), 
\label{eq:a19}\end{aligned}$$ where $\mathbf{E} \left[ X \right] = \frac{1}{1-g^2}$, $\mathbf{E} \left[ X_n \right] = \frac{1}{1-g^2} - g^{2n}$, $\mathrm{Var}\left[ X \right] = \frac{2g^4}{1-g^4}$ and $\mathrm{Var}\left[ X_n \right] = \frac{2g^4}{1-g^4} - 2g^{4n}$. $\mathbf{E}\left[ I_\mathrm{ot} \right]$ can be computed by summing Eq. (\[eq:a19\]) over $n \geq 1$.

The authors are grateful to the anonymous reviewers for their comments and suggestions that improved the manuscript. TH was supported by JSPS KAKENHI Grant Number JP18K03423. KN was supported by JSPS KAKENHI Grant Number JP18H05472 and by MEXT Quantum Leap Flagship Program (MEXT Q-LEAP) Grant Number JPMXS0118067394. This work is partially based on results obtained from a project commissioned by the New Energy and Industrial Technology Development Organization (NEDO).

[35]{}ifxundefined \[1\][ ifx[\#1]{} ]{}ifnum \[1\][ \#1firstoftwo secondoftwo ]{}ifx \[1\][ \#1firstoftwo secondoftwo ]{}““\#1””@noop \[0\][secondoftwo]{}sanitize@url \[0\][‘\
12‘\$12 ‘&12‘\#12‘12‘\_12‘%12]{}@startlink\[1\]@endlink\[0\]@bib@innerbibempty @noop [“,” ]{} (),  @noop [****,  ()]{} @noop [****,  ()]{} @noop [****,  ()]{} @noop [****,  ()]{} @noop [****,  ()]{} @noop [****,  ()]{} @noop [****,  ()]{} @noop [****,  ()]{} @noop [****,  ()]{} @noop [****,  ()]{} @noop [****,  ()]{} @noop [“,” ]{} (),  @noop [****,  ()]{} @noop [****,  ()]{} @noop [****,  ()]{} @noop [****,  ()]{} @noop [****,  ()]{} @noop [****,  ()]{} @noop [****,  ()]{} @noop [****,  ()]{} @noop [****,  ()]{} @noop [****,  ()]{} @noop [****,  ()]{} @noop [****,  ()]{} @noop [****,  ()]{} @noop [****,  ()]{} @noop [****,  ()]{} @noop [**]{} (, , ) @noop [****,  ()]{} @noop [****,  ()]{} @noop [****,  ()]{} @noop [****,  ()]{} @noop [****,  ()]{} @noop [****,  ()]{}
