Stability Analysis of Reservoir Computers Dynamics via Lyapunov Functions

TL;DR

This paper applies Lyapunov functions to analyze the nonlinear stability of reservoir computers, identifying regions of stability and linking stability to the polynomial structure of the reservoir dynamics.

Contribution

It introduces a Lyapunov-based method for stability analysis of reservoir computers and analytically determines stability regions for both continuous and discrete systems.

Findings

Training error is lower within the stability region.

Stability is influenced by the polynomial structure of the reservoir dynamics.

Nonzero coefficients for odd and even powers are important for polynomial dynamics.

Abstract

A Lyapunov design method is used to analyze the nonlinear stability of a generic reservoir computer for both the cases of continuous-time and discrete-time dynamics. Using this method, for a given nonlinear reservoir computer, a radial region of stability around a fixed point is analytically determined. We see that the training error of the reservoir computer is lower in the region where the analysis predicts global stability but is also affected by the particular choice of the individual dynamics for the reservoir systems. For the case that the dynamics is polynomial, it appears to be important for the polynomial to have nonzero coefficients corresponding to at least one odd power (e.g., linear term) and one even power (e.g., quadratic term).

Equations

$${\dot{r}}_i(t)=f(r_i(t),\mathbf{\theta })+\sum _{j=1}^M{A}_{ij}{r}_j(t)+w_{i}s(t),i=1,2,\cdots ,M$$

$${\dot{r}}_i(t)=f(r_i(t),\mathbf{\theta })+\sum _{j=1}^M{A}_{ij}{r}_j(t),i=1,2,\cdots ,M$$

$$V(\text{r})=\frac{1}{2}{\text{r}}^T\text{r},$$

$$\frac{\partial V}{\partial \text{r}}\text{f}(\text{r},\mathbf{\theta })=$$

$$=$$

$$\le$$

$$\sum _i{r}_{i}f(r_i,\mathbf{\theta })+{\alpha }_{\max }\Vert \text{r}{\Vert }^2$$

$$\le$$

$$=$$

$$\left(K(c,\mathbf{\theta })+{\alpha }_{\max }\right)\Vert \text{r}{\Vert }^2\le 0$$

$$K(c,\mathbf{\theta })\le -{\alpha }_{\max }$$

$$\underset{K}{\min }$$

$$r_{i}f(r_i)-K{r}_i^2\le 0,\forall r_i\in [-c,c]$$

$$K^{\ast }(c,\mathbf{\theta }){r_i^{\ast }}^2-r_i^{\ast }f(r_i^{\ast },\mathbf{\theta })$$

$$K^{\ast }(c,\mathbf{\theta })$$

$$K^{\ast }(c,\mathbf{\theta })=\max \{\begin{array}{l}\frac{f(c,\mathbf{\theta })}{c},\\ \frac{f(-c,\mathbf{\theta })}{-c},\\ f^{\prime }(0,\mathbf{\theta }),\\ \frac{f(r_i^{\ast },\mathbf{\theta })}{r_i^{\ast }},\text{ where }{r}_i^{\ast }\text{ is the solution of }\\ r_i^{\ast }{f}^{\prime }(r_i^{\ast },\mathbf{\theta })-f(r_i^{\ast },\mathbf{\theta })=0,r_i^{\ast }\in [-c,c],r_i^{\ast }\ne 0\end{array}$$

$$K^{\ast }(c,\mathbf{\theta })\le -{\alpha }_{\max }$$

$$f(r_i,\mathbf{\theta })=p_1{r}_i+p_2{r}_i^2+p_3{r}_i^3,$$

$${\dot{r}}_i(t)=p_1{r}_i+p_2{r}_i^2+p_3{r}_i^3+\sum _{j=1}^M{A}_{ij}{r}_j+w_{i}s(t).$$

$$K^{\ast }(c,\mathbf{\theta })$$

$$r_i(n+1)=f(r_i(n),\mathbf{\theta })+\sum _{j=1}^M{A}_{ij}{r}_j(n)+w_{i}s(n),i=1,2,\cdots ,M,$$

$$r_i(n+1)=f(r_i(n),\mathbf{\theta })+\sum _{j=1}^M{A}_{ij}{r}_j(n),i=1,2,\cdots ,M,$$

$$V(\text{r})=\Vert \text{r}\Vert$$

$$V(f(\text{r},\mathbf{\theta })+A\text{r})-V(\text{r})$$

$$=\Vert f(\text{r},\mathbf{\theta })+A\text{r}\Vert -\Vert \text{r}\Vert$$

$$\Vert f(\text{r},\mathbf{\theta })+A\text{r}\Vert -\Vert \text{r}\Vert \le \Vert K(c,\mathbf{\theta })\text{r}+A\text{r}\Vert -\Vert \text{r}\Vert$$

$$|K(c,\mathbf{\theta })I+A|\Vert \text{r}\Vert -\Vert \text{r}\Vert \le 0$$

$$|K(c,\mathbf{\theta })+{\gamma }_i|\le 1,$$

$$|K+{\gamma }_i^{re}|\le \sqrt{1-({\gamma }_i^{im}{)}^2}$$

$$-\sqrt{1-({\gamma }_i^{im}{)}^2}\le K+{\gamma }_i^{re}\le \sqrt{1-({\gamma }_i^{im}{)}^2}$$

$${\rho }_{c-}^{-}\le K^{-}(c,\mathbf{\theta })\le K^{+}(c,\mathbf{\theta })\le {\rho }_{c+}^{+}$$

$$0\le \left|\frac{f(r_i,\mathbf{\theta })}{r_i}+{\gamma }_i^{re}\right|$$

$${\rho }_{c-}^{-}\le K^{-}(c,\mathbf{\theta })\le$$

$${\rho }_{c-}^{-}\le K^{-}(c,\mathbf{\theta })\le \frac{f(r_i,\mathbf{\theta })}{r_i}$$

Full text

While nonlinearity appears to be a fundamental component of reservoir computers, not much research has been performed to analyze stability of the nonlinear dynamics of these systems. In this paper, we use a Lyapunov design method to estimate the basin of attraction of a fixed point for the dynamics of a generic reservoir computer. Our nonlinear stability analysis unveils a trade-off between the need for global stability, which is achievable by linear dynamics alone, and the need for higher-order terms of the dynamics, which could in turn compromise stability.

I Introduction

A reservoir computer (RC) is a complex nonlinear dynamical system that is used for processing and analyzing empirical data, see e.g. jaeger2001echo ; schrauwen2007overview ; natschlager2002liquid ; maass2002real ; martinenghi2012photonic ; brunner2013parallel ; nakajima2015information ; hermans2015photonic ; vinckier2015high ; duport2016fully ; larger2017high , modeling of complex dynamical systems suykens2012artificial , speech recognition crutchfield2010introduction , learning of context free and context sensitive languages rodriguez2001simple ; gers2001lstm , the reconstruction and prediction of chaotic attractors lu2018attractor ; zimmermann2018observing ; antonik2018using ; jaeger2004harnessing ; pathak2017using ; pathak2018model , image recognition jalalvand2018application , and control of robotic systems graves2004biologically ; robinson1994application ; lukovsevivcius2012reservoir .

A typical RC consists of a set of nodes coupled together to form a network. Each node of the RC evolves in time in response to an input signal that is externally provided to the reservoir. An output signal is then generated from the time evolutions of the RC nodes. In a RC, the output connections (those that connect the RC nodes to the output) are trained to produce a best fit between the output signal and a training signal related to the original input signal. On the other hand, the connections between the nodes of the reservoir are constant parameters of the system. As a result, RCs are easier to analyze than other machine learning tools for which all the connections are typically trained.

The functions of RCs mainly depend on two factors; (i) nonlinearity of the nodal dynamics which is needed to process the information in the input signal and (ii) linear memory to boost the excitability of the RC dynamics dambre2012information . Though earlier works have shown that maximizing linear memory is important in information processing jaeger2002short ; ganguli2008memory ; buonomano2009state ; white2004short , more recent works have shown that the performance of a reservoir computer is related to consideration of both factors (i) and (ii) verstraeten2007experimental ; inubushi2017reservoir ; marzen2017difference . In addition, the performance of a reservoir computer is also affected by a number of other factors, including the reservoir adjacency matrix, i.e., the strengths of the connections between the RC nodes, and the dynamic range of the input signals verstraeten2009quantification .

From linear systems theory, a dynamical system is reliable and safe when it is stabilizable around some operating point sontag2013mathematical . Previous research has used linear stability to assess the stability of RCs around this operating point marcus1989stability ; tanaka1995stability ; belair1996frustration . However, as the RC dynamics requires nonlinearity for its proper operation, a linear stability analysis of the RC dynamics around a specific point is not sufficient. This motivates us to develop a nonlinear stability analysis, based on Lyapunov functions, which we use to characterize the basin of attraction of the desired operating point. As assessing stability of the nonlinear system when forced by an external stimulus is contingent on the particular stimulus provided, we characterize nonlinear stability of the unforced RC dynamics. If the desired operating point is found to be globally asymptotically stable, then stability is independent of the particular stimulus, as long as it is bounded.

We compute a constant $c$-radius region around the operating point such that the dynamics of the system remains bounded inside this region. This can be done by choosing the parameters of the system and the type of nonlinearity, with the goal of possibly enhancing the performance of the reservoir computer. We consider different types of nonlinear dynamics at the network nodes, e.g., polynomial, in either continuous time or discrete time. This differs from the common approach in the literature, where the nodal dynamics is chosen to have a squashing type nonlinearity (in most cases a sigmoid function, e.g., $\tanh()$) so that the states of the nonlinear system always remain bounded inside some region liu2001stability ; saul2000attractor ; yang2014exponential . (A squashing function is defined as a function that is monotonic and bounded within a small range. For example, the function $\tanh()$ squashes the argument to the interval $[-1,1]$.)

Few theoretical works have investigated how the underlying stability of a nonlinear RC affects its performance. Reference ganguli2008memory showed that the total memory capacity, a measure of performance, is related to the size of the network. This analysis did not consider how the underlying stability of the system is related to the total memory capacity. In Ref. verstraeten2010memory , an optimal parameter setting for the reservoir was studied, but no direct relation between the dynamic properties in terms of nonlinear stability and the reservoir performance was provided. Previous work verstraeten2009quantification found that the performance of the RC was improved when the condition number of the Jacobian of the reservoir dynamics was small. The references listed herein, and others, motivated us to perform a rigorous dynamical investigation of the nonlinear stability of RCs.

In this paper, for a given nonlinear RC, we determine the $c(\bm{\theta})$-region of stability, where $c$ is the radius around a fixed point and $\bm{\theta}$ is the set of parameters of the reservoir. We use Lyapunov design methods haddad2011nonlinear to find the $c(\bm{\theta})$-region where a nonlinear reservoir computer is stable, safe, reliable and its performance is predictable. Lyapunov design methods are used widely in controls engineering to design controllers that achieve qualitative objectives, such as stabilizing a system or maintaining a system’s state in a desired operating range sontag2013mathematical . To the best of our knowledge, the use of a Lyapunov-based design approach with respect to the performance of a reservoir computer is novel. In Ref. perkins2002lyapunov , a Lyapunov function has been used to design the controller of a reinforcement learning system, but the paper does not show how stability of the RC affects its performance. In this article, first we assume that the input signal is normalized and scaled properly, the eigenvalues of the adjacency matrix satisfy certain constraints, and second we design the $c(\bm{\theta})$-region of stability by using a Lyapunov design method. Our nonlinear stability analysis provides insight into the effects of the RC parameters on its performance.

In Sec. II we lay out the general theory to assess the basin of attraction within which the RC dynamics is stable for both continuous-time dynamics and discrete-time dynamics. In Sec. III we investigate the relation between our stability predictions and the RC performance in terms of the computed training error. Finally, conclusions are presented in Sec. IV.

II Methods

II.1 Reservoir Computer with Continuous Time Dynamics

We consider the dynamics of a reservoir computer as continuous-time carroll2019network ,

[TABLE]

and the unforced (without input) reservoir dynamics,

[TABLE]

where $r_{i}(t)\in\mathbb{R},i=1,2,\cdots,M$ denotes the state of node $i$ at time $t\in\mathbb{R}$, $f(r_{i},\bm{\theta})$ is the nonlinear nodal dynamics at node $i$, the adjacency matrix $A=\{A_{ij}\}$ indicates the coupling from node $j$ to node $i$, $s(t)\in\mathbb{R}$ is an input signal and $\textbf{w}=\left[w_{1},w_{2},\cdots,w_{M}\right]$ is a vector describing the coupling of input signal $s(t)$ with node $i$. We assume that $r_{i}^{*}(t)=0$ is a (linearly) stable fixed point of the system in Eq. (2) and the input signal $s(t)$ in Eq. (1) is normalized to have zero mean and standard deviation equal to one.

II.1.1 Lyapunov Function and $c(\bm{\theta})$-region design

We define a Lyapunov function $V:\mathcal{D}(c)\subseteq\mathbb{R}^{M}\rightarrow\mathbb{R}$ for the unforced dynamical equation in Eq. (2),

[TABLE]

where $\mathcal{D}(c)=\left\{\textbf{r}=(r_{1},r_{2},\cdots,r_{M})\in\mathbb{R}^{M}:\|\textbf{r}\|\leq c\right\}$ is the phase space region that is included in a hypersphere of radius $c$, centered at the origin. Here $V(\textbf{r})>0$ for $\textbf{r}\neq 0$ and $V(\textbf{0})=0$. Then for all $\textbf{r}\in\mathcal{D}(c)\backslash\textbf{0}$,

[TABLE]

where $A_{s}$ is the symmetric part of the matrix $A$, $\alpha_{\max}$ is the largest real eigenvalue of the matrix $A_{s}$, and $\bm{\theta}$ is the set of parameters that completely characterize the nonlinear function $f$. Now we introduce a quadratic upper bound to the term $r_{i}f(r_{i},\bm{\theta})$. Let us consider that the term $K(c,\bm{\theta})r_{i}^{2}$ is such a quadratic upper bound, that is $r_{i}f(r_{i},\bm{\theta})\leq r_{i}^{2}K(c,\bm{\theta})$, where $K(c,\bm{\theta})$ is a scalar function of $c$ and $\bm{\theta}$. The inequality in Eq. (4e) can now be written as,

[TABLE]

According to the second Lyapunov stability theorem haddad2011nonlinear , the system in Eq. (2) is stable within $\mathcal{D}$ if the following inequality holds,

[TABLE]

Note that for each $c$, the term on the left hand side of Eq. (6b) depends on the dynamics and parameters of the individual nodes and the term on the right hand side of Eq. (6b) depends on the network topology. Thus Eq. (6b) effectively decouples the stability problem into two terms that can be adjusted independently of each other: the nodal dynamics and the network topology.

We now provide a definition of $c(\bm{\theta})$-region stability of a reservoir computer.

Definition 1

A nonlinear reservoir computer is $c(\bm{\theta})$-region stable if $K(c,\bm{\theta})\leq-\alpha_{\max}$.

Note that this is a sufficient (but not necessary) condition for a reservoir to be stable inside the region $\mathcal{D}(c)$. Also, if $c\rightarrow\infty$, the system is globally asymptotically stable.

We note that for any constant $K\leq K(c,\bm{\theta})$, the system is $c(\bm{\theta})$-region stable. We will thus attempt to find an upper bound $K(c,\bm{\theta})r_{i}^{2}$ to $r_{i}f(r_{i},\bm{\theta})$ that is as tight as possible. To find the minimal $K\leq K(c,\bm{\theta})$, we define an optimization problem,

[TABLE]

As $r_{i}$ lies within the closed interval $[-c,c]$, the constraint in Eq. (7) must be satisfied along a continuum. However, we know that the constraint in Eq. (7b) achieves equality for some $r_{i}^{*}\in[-c,c]$ at the optimal solution $K^{*}(c,\bm{\theta})$,

[TABLE]

or equivalently,

[TABLE]

The optimal coefficient $K^{*}(c,\bm{\theta})$ is chosen as the term of maximum value in the set,

[TABLE]

where the four cases in Eq. (10) are the possible maxima of the ratio $f(r_{i},c)/r_{i}$ over a closed interval. The stability condition for the reservoir computer is

[TABLE]

From the inequality in Eq. (11) and the definition of $K^{*}(c,\bm{\theta})$ in Eq. (10), we can find $c_{\max}$ for which $K^{*}(c_{\max},\bm{\theta})=-\alpha_{\max}$ and the $c(\bm{\theta})$-region is determined by $\mathcal{D}(c_{\max})=\{\textbf{r}=(r_{1},r_{2},\cdots,r_{M})\in\mathbb{R}^{M}:\|\textbf{r}\|\leq c_{\max}\}$. We call $c_{\max}$ the radius of the region $\mathcal{D}(c_{\max})$. If $\lim\limits_{c\rightarrow\infty}K^{*}(c,\bm{\theta})<-\alpha_{\max}$ then we say the system is globally stable (less formally we say $c_{\max}=\infty$), while if $K^{*}(0,\bm{\theta})>-\alpha_{\max}$ then the system is unstable. In the next subsection, we will find the $c(\bm{\theta})$-region for the case that the nonlinear nodal dynamics is described by a polynomial.

II.1.2 Polynomial Type Nonlinearity

We now consider an example for which the reservoir computer consists of $M$ homogeneous nodes and the nodal dynamics of each node $i,i=1,2,\cdots,M$ is defined by the following third-order polynomial function carroll2019network ; carroll2019mutual ,

[TABLE]

Polynomial functions are a very general way to express nonlinearity.

The dynamical equation that governs the evolution of each node $i$ is,

[TABLE]

Here, $p_{1}$, $p_{2}$ and $p_{3}$ are the coefficients of the polynomial. In this case, the set of parameters $\bm{\theta}=\{p_{1},p_{2},p_{3}\}$. The origin $r_{i}=0$ is a fixed point for the dynamics and is linearly stable if the largest real part of the eigenvalues of the matrix $(A-p_{1}I)$ is negative. Therefore, in what follows, we fix $p_{1}$ so as to ensure that the matrix $(A-p_{1}I)$ is Hurwitz and then we characterize the basin of attraction as a function of the remaining parameters $p_{2}$ and $p_{3}$.

According to Eq. (10), the scalar function $K^{*}(c,\bm{\theta})$ can be obtained as,

[TABLE]

We note that the condition $K^{*}(c,\bm{\theta})=-\alpha_{\max}(A_{s})$ determines the radius $c_{\max}$ of the sphere $\mathcal{D}$ for which the reservoir computer is $c$-region stable. Global stability is achieved when $K^{*}(c,\bm{\theta})$ remains upper bounded by $-\alpha_{\max}(A_{s})$ as $c\rightarrow\infty$.

Here, we provide an example to explain how to find the $c(\bm{\theta})$-region for a simple reservoir computer with $M=2$ nodes. We set $p_{1}=-3$, $p_{3}=-1$, $A=\begin{bmatrix}0&1\\ -1&0\end{bmatrix}$ and let the parameter $p_{2}$ vary. In Fig. 1(A), we plot $K^{*}(c,\bm{\theta})$ versus $c$ for different values of $p_{2}$. The solid black line is the constant-ordinate line at $-\alpha_{\max}=0$. For this example, we observe that $K^{*}(c,\bm{\theta})$ is symmetric about the parameter $p_{2}=0$. For $p_{2}=\pm 4$ , $c_{max}=1$ which is represented by a black dot where the curves for $p_{2}=\pm 4$ cross $-\alpha_{\max}$. For $p_{2}=\pm 1,\pm 1$, $K^{*}(c,\bm{\theta})$ reaches a constant below $-\alpha_{\max}$ as $c$ grows, indicating that the basin of attraction has infinite radius. Figure 1(B) considers the case that $p_{2}=\pm 1$. We see two different regions in the $r_{1}(0),r_{2}(0)$-plane distinguished by two colors: the red region indicates the initial conditions from which the system’s time evolution approaches the origin as time grows and the yellow region indicates the initial conditions from which the system’s time evolution does not converge to the origin, in which case the dynamics converges to either another fixed point, or a limit cycle, or any other attractor other than the origin. The black circle is the solution of $\|\textbf{r}\|=c_{\max}$, for $p_{2}=\pm 4$. In Figs. 1(C) and 1(D) we plot the trajectory $r_{2}(t)$ versus $r_{1}(t)$ when the system is evolved from a typical initial condition from within the red and the yellow region, respectively.

II.2 Reservoir Computer with Discrete Time Dynamics

We now turn to the dynamics of a reservoir computer with discrete time dynamics,

[TABLE]

which in the unforced case becomes,

[TABLE]

where $r_{i}(n)\in\mathbb{R}$, $i=1,2,\cdots,M$, denotes the state of the node $i$ of the reservoir at time step $n$, $f(r_{i},\bm{\theta})$ is the nonlinear nodal dynamics of node $i$, the adjacency matrix $A=\{A_{ij}\}$ indicates the pattern of connectivity between the network nodes, $s(n)\in\mathbb{R}$ is the input signal at time step $n$ and $\textbf{w}=\left[w_{1},w_{2},\cdots,w_{M}\right]$ is a vector that describes the coupling of the input signal $s(n)$ to each one of the nodes. The input signal $s(n)$ in Eq. (15) is normalized to have mean [math] and standard deviation equal to $1$ lu2017reservoir ; carroll2019network .

Hereafter we assume that the operating fixed point for the dynamics Eq. (16) coincides with the origin (however, this assumption can be removed; see the example that follows for the case of a sigmoid nonlinearity.)

II.2.1 Lyapunov Function and $c(\bm{\theta})$-region Design

We define a Lyapunov function $V:\mathcal{G}(c)\subseteq\mathbb{R}^{M}\rightarrow\mathbb{R}$

[TABLE]

where $\mathcal{G}(c)=\left\{\textbf{r}=(r_{1},r_{2},\cdots,r_{M})\in\mathbb{R}^{M}:\|\textbf{r}\|\leq c\right\}$ and $\|\textbf{r}\|=\sqrt{r_{1}^{2}+\cdots+r_{M}^{2}}$. Here $V(\textbf{r})>0$ for $\textbf{r}\neq 0$ and $V(\textbf{0})=0$. Then $\text{ for all }\textbf{r}\in\mathcal{G}(c)\backslash\textbf{0}$,

[TABLE]

We seek to find a scalar function $K(c,\bm{\theta})$ such that $f(r_{i},\bm{\theta})\leq K(c,\bm{\theta})r_{i}$ which also satisfies the inequality in Eq. (18b),

[TABLE]

According to the Lyapunov stability theorem for discrete time dynamics haddad2011nonlinear , the system is stable only if,

[TABLE]

where $\gamma_{i}=\gamma_{i}^{re}+j\gamma_{i}^{im}$ is an eigenvalue of the matrix $A$ and $j=\sqrt{-1}$. The above inequality can be written as,

[TABLE]

We see that there is both an upper bound and a lower bound for $K(c,\bm{\theta})$. Hence, there are two critical eigenvalues: $\gamma_{c+}$ and $\gamma_{c-}$ which are the eigenvalues closest to the positive side of the unit circle and closest to the negative side of the unit circle when moving only along the real axis, respectively. This concept is displayed graphically in Fig. 2 where $\gamma_{2}=\gamma_{c+}$ and $\gamma_{5}=\gamma_{c-}$. The maximum distance the eigenvalue $\gamma_{c+}$ can shift to the right is $\rho_{c+}^{+}=\sqrt{1-(\gamma_{c+}^{im})^{2}}-\gamma_{c+}^{re}$ and the maximum distance the eigenvalue $\gamma_{c-}$ can shift to the left is $\rho_{c-}^{-}=-(\sqrt{1-(\gamma_{c-}^{im})^{2}}+\gamma_{c-}^{re})$. Thus there exist two scalar functions denoted by $K(c,\bm{\theta})=K^{-}(c,\bm{\theta})\leq 0$ and $K(c,\bm{\theta})=K^{+}(c,\bm{\theta})\geq 0$ such that,

[TABLE]

An illustration is presented in Fig. 2 which shows how to find the critical eigenvalues $\gamma_{c+}$ and $\gamma_{c-}$. In Fig. 2, several eigenvalues of some hypothetical adjacency matrix $A$ are shown inside the unit circle. For each eigenvalue, we compute $\rho_{i}^{+}$ and $\rho_{i}^{-}$. From the table we see that $\gamma_{2}$ is the critical eigenvalue $\gamma_{c+}$ and $\gamma_{5}$ is the critical eigenvalue $\gamma_{c-}$.

Now using the fact that all the nodes are homogeneous, from inequality Eq. (19), we can write,

[TABLE]

From the inequalities in Eqs. (22) and (23), it follows that,

[TABLE]

Thus, we find $K^{-}(c,\bm{\theta})$ and $K^{+}(c,\bm{\theta})$ such that

[TABLE]

As we want tight upper and lower bounds on $\frac{f(r_{i},\bm{\theta})}{r_{i}}$, we seek to find $K^{+}(c,\bm{\theta})$ and $K^{-}(c,\bm{\theta})$ that solve the following two optimization problems,

[TABLE]

and

[TABLE]

A solution ${K^{-}}^{*}(c,\bm{\theta})$ to the problem in Eq. (26) and ${K^{+}}^{*}(c,\bm{\theta})$ to the problem in Eq. (27) must satisfy their respective constraints exactly for some $r_{i}^{+*}$ and $r_{i}^{-*}$,

[TABLE]

or equivalently,

[TABLE]

We can find ${K^{-}}^{*}(c,\bm{\theta})$ as

[TABLE]

and we can find ${K^{+}}^{*}(c,\bm{\theta})$ as

[TABLE]

Once we obtain ${K^{-}}^{*}(c,\bm{\theta})$ and ${K^{+}}^{*}(c,\bm{\theta})$, we can find $c_{\max}^{+}$ and $c_{\max}^{-}$ such that ${K^{-}}^{*}(c_{\max}^{-},\bm{\theta})=\rho_{c-}^{-}$ and ${K^{+}}^{*}(c_{\max}^{+},\bm{\theta})=\rho_{c+}^{+}$, respectively. The $c(\bm{\theta})$-region for the discrete time RC can be determined as $\mathcal{G}(c_{\max})=\{\textbf{r}=(r_{1},r_{2},\cdots,r_{M})\in\mathbb{R}^{M}:\|\textbf{r}\|\leq c_{\max}\}$, where $c_{\max}=\min\{c_{\max}^{-},c_{\max}^{+}\}$.

Example: $\tanh()$ type nonlinearity

We choose the nodal dynamics to be

[TABLE]

The dynamics of node $i$ is described by,

[TABLE]

Here $\bm{\theta}=\{p_{1},p_{2}\}$. We see that for $s(n)=0$, the origin is a fixed point, which is stable if all the eigenvalues of the matrix $(A+p_{1}p_{2}I)$ are inside the unit circle. The constant functions ${K^{-}}^{*}(c,\bm{\theta})$ and ${K^{+}}^{*}(c,\bm{\theta})$ can be found as,

[TABLE]

and

[TABLE]

Now if $p_{1}<0$, then $p_{1}p_{2}={K^{-}}^{*}(c,\bm{\theta})\leq 0$ and if $p_{1}>0$, then $0\leq{K^{+}}^{*}(c,\bm{\theta})=p_{1}p_{2}$ for any choice of $c$.

II.2.2 Lyapunov Function and $c(\bm{\theta})$-region Design for Non-homogeneous Nodal Dynamics

One generalization of Eq. (16) is to the case of non-homogeneous nodal dynamics,

[TABLE]

Without loss of generality we retain the assumption that the above set of equations has a fixed point at the origin. In the case a fixed point exists that is different from the origin, this assumption can be removed by applying a coordinate transformation that moves the fixed point to the origin (see the example that follows for the case of a sigmoid nonlinearity). Now according to the Lyapunov function analysis described in section B.1, for each node $i$ we find scalar functions $K_{i}^{-*}(c,\bm{\theta}_{i})$ and $K_{i}^{+*}(c,\bm{\theta}_{i})$ which satisfy,

[TABLE]

We can find ${K_{i}^{-}}^{*}(c,\bm{\theta}_{i})$ as

[TABLE]

and we can find ${K_{i}^{+}}^{*}(c,\bm{\theta}_{i})$ as

[TABLE]

Now we define the scalar function ${K^{+}}^{*}(c,\bm{\theta})$ as,

[TABLE]

and the scalar function ${K^{+}}^{*}(c,\bm{\theta})$ as,

[TABLE]

Once we obtain ${K^{-}}^{*}(c,\bm{\theta})$ and ${K^{+}}^{*}(c,\bm{\theta})$, we can find $c_{\max}^{+}$ and $c_{\max}^{-}$ such that ${K^{-}}^{*}(c_{\max}^{-},\bm{\theta})=\rho_{c-}^{-}$ and ${K^{+}}^{*}(c_{\max}^{+},\bm{\theta})=\rho_{c+}^{+}$, respectively. Then the $c$-region for the reservoir computer can be determined as $\mathcal{G}(c_{\max})=\{\textbf{r}-\textbf{q}^{*}:\|\textbf{r}-\textbf{q}^{*}\|\leq c_{\max}\}$, where $c_{\max}=\min\{c_{\max}^{+},c_{\max}^{-}\}$.

Example: Sigmoid Nonlinearity

Here we choose the nodal dynamics to be

[TABLE]

The unforced dynamics of each node $i$ is,

[TABLE]

Here the parameters are $\bm{\theta}=\{p_{1},p_{2}\}$. For this example, we see that the origin is not a fixed point. Instead a nonzero fixed point exists: $q_{i}^{*}$, $i=1,...,M$. Let $\bar{r}_{i}=r_{i}-q_{i}^{*}$ and the system in Eq. (43) can be transformed to the form of Eqs. (36), where $\bar{f}_{i}(\bar{r}_{i}(n),\bm{\theta}_{i})=f((\bar{r}_{i}(n)+q_{i}^{*},\bm{\theta})+\sum_{j}A_{ij}q_{j}^{*}-q_{i}^{*}$. According to Eq. (38) and (39), we find the scalar functions ${K_{i}^{-}}^{*}(c,\bm{\theta}_{i})$ and ${K_{i}^{+}}^{*}(c,\bm{\theta}_{i})$ as follows,

[TABLE]

and

[TABLE]

Now we define the scalar function ${K^{-}}^{*}(c,\bm{\theta})$ as follows,

[TABLE]

and the scalar function ${K^{+}}^{*}(c,\bm{\theta})$ as follows,

[TABLE]

Once we obtain ${K^{-}}^{*}(c,\bm{\theta})$ and ${K^{+}}^{*}(c,\bm{\theta})$, we can find $c_{\max}^{+}$ and $c_{\max}^{-}$ such that ${K^{-}}^{*}(c_{\max}^{-},\bm{\theta})=\rho_{c-}^{-}$ and ${K^{+}}^{*}(c_{\max}^{+},\bm{\theta})=\rho_{c+}^{+}$, respectively. Then the $c$-region for the reservoir computer can be determined as $\mathcal{G}(c_{\max})=\{\textbf{r}-\textbf{q}^{*}:\|\textbf{r}-\textbf{q}^{*}\|\leq c_{\max}\}$, where $c_{\max}=\min\{c_{\max}^{+},c_{\max}^{-}\}$.

III Results

For our numerical simulations, in both continuous-time and discrete-time, we construct the adjacency matrix $A$ as follows: (i) We set the entries of the initial matrix $A$ to be equal to $A_{ij}=1-\delta_{ij}$, where $\delta_{ij}$ is the Kronecker delta, $i,j=1,...,M$. (ii) $50\%$ of the off-diagonal entries of the matrix $A$ are chosen randomly and set to zero. (iii) $50\%$ of the remaining nonzero entries of the matrix $A$ are chosen randomly and are flipped from $+1$ to $-1$. (iv) Finally, the adjacency matrix $A$ is normalized so that the absolute value of the largest real part of its eigenvalues is equal to $0.5$.

III.1 Training Error of a Reservoir Computer

The training error, $\Delta_{RC}$, is used to quantify how well the training signal $g(t)$ ($g(n)$, in the case of discrete dynamics) can be reconstructed from the input signal $s(t)$ ($s(n)$). Lower values of $\Delta_{RC}$ indicate a better performance of the reservoir computer. In the continuous-time case, the training signal and the input signal are discretized in time and are thus treated as sequences. Before driving the reservoir computer by the input signal $s(t)$ ($s(n)$), both $s(t)$ and $g(t)$ ($s(n)$ and $g(n)$) are normalized to have mean equal to zero and standard deviation equal to one lu2017reservoir . However, such a choice is completely arbitrary; the amplitude of the input signal can be conveniently reduced in case one finds the (stable) reservoir dynamics to become unstable when driven by the input signal. Next, we present the procedure to compute the training error in the case of discrete-time dynamics (the procedure for the case of continuous-time dynamics is analogous.) We set the number of nodes of the reservoir to $M=100$. When the reservoir is driven by the input signal $s(n)$, the first 2000 time steps are discarded as a transient. The next $N=10,000$ time steps from each node are recorded and combined into the $N\times(M+1)$ matrix,

[TABLE]

The last column of the matrix $\Omega$ is set to 1 to account for any constant offset in the fitting. We then introduce

[TABLE]

where ${\textbf{h}}=\left[h(1),h(2),\cdots,h(N)\right]$ is the fit to the training signal ${\textbf{g}}=\left[g(1),g(2),\cdots,g(N)\right]$ and $\textbf{k}=\left[k_{1},k_{2},\cdots,k_{M+1}\right]^{T}$ is the vector of coefficients. The vector k is obtained from the minimum-norm solution to the linear least squares problem,

[TABLE]

The training error is computed as,

[TABLE]

where the symbol $\langle\cdot\rangle$ is computed by using the following formula,

[TABLE]

where $\textbf{X}=\left[X(1),X(2),\cdots,X(N)\right]$ and $\mu=\frac{1}{N}\sum_{i=1}^{N}X(i)$.

III.2 Results for Continuous-Time Dynamics

We now consider continuous-time and use the general polynomial function for the individual nodal dynamics,

[TABLE]

of which Eq.  (13) is an example. We keep $p_{1}$ constant ($p_{1}=-3$) as we are mainly interested in the effects of the nonlinear terms. In what follows we will study the effects of varying different pairs of $p_{i}$ coefficients, $i>1$, of the polynomial (53) (in addition to setting $p_{1}=-3$.) For different cases we will indicate the pairs of coefficients we are focusing on, with the understanding that all the remaining coefficients are set to zero. We see that carefully choosing the form of the polynomial (53) is important and that setting certain coefficients $p_{i}$ to zero can dramatically worsen the performance of the reservoir, even when this is not directly predicted by the stability analysis.

Figure 3 provides a visual assessment of the training error of the reservoir computer in terms of the parameters $p_{2}$ and $p_{3}$. First we construct the matrix $A$ as described at the beginning of this section. For the input signal $s(t)$ we choose the $x$ signal from a Lorenz chaotic attractor, while the training signal $g(t)$ is the Lorenz $z$ signal pathak2017using . The input signal $s(t)$ and the training signal $g(t)$ are normalized to have mean equal to zero and standard deviation equal to one. In Fig. 3 we plot the training error $\Delta_{RC}$ as a function of the parameters $p_{2}$ and $p_{3}$. The color represents variations in log scale of the training error from dark blue (small) to dark red (large). The solid black curve represents the boundary between sets of parameters such that the RC is globally stable (below the black line) versus sets of parameters such that $c_{\max}$ is finite (above the black curve). In other words, the region under the black curve determines the part of the parameter space $(p_{2},p_{3})$ for which the reservoir computer is globally stable and the RC is successful in performing the computation (though the training error may vary by different orders of magnitude as the parameters change.) On the other hand, above the black curve, it is difficult to assess how the system behaves. For example, the system could be globally stable, or locally stable, in which case depending on the particular choice of the input, the system dynamics might approach a different attractor or might be driven to $\pm\infty$. Another interesting observation is the presence of a tiny triangular region above the black curve where the reservoir computer performs well. The white region is the area of the parameter space for which the system goes to $\pm\infty$ when it is driven by the input signal and the reservoir computer fails in performing the computation. For $p_{3}>2$ the RC dynamics diverges independent of the choice of $p_{2}$. We also notice that the training error is symmetric about the $p_{2}=0$-axis and the training error is very high (almost equals to 1) at $p_{2}=0$, which indicates that the reservoir computer fails to capture and transfer the information from input to output if the quadratic term is absent from the nodal dynamics. Note that this seems to indicate there are two distinct requirements for the reservoir to work properly: (i) it needs to operate in the area of global stability and (ii) the $p_{2}$ coefficient needs to be nonzero. We observe a similar behavior in Fig. 4 where the input signal $s(t)$ is the $x$- component and the training signal is the $y$-component of the Duffing chaotic attractor chang2017chaotic . Figure 7 shows in further detail the results in Fig. 3 for $p_{3}=-4$ while preserving the nodal dynamics, the adjacency matrix $A$, etc. Figure 7(A) is a plot of $K^{*}(c\rightarrow\infty,\bm{\theta})$ versus the parameter $p_{2}$ when $p_{3}=-4$ and the black line is the constant-ordinate line at $-\alpha_{\max}$ of the symmetric part of the matrix $A$. In Fig. 7(B), we plot the inverse of the radius of the $c(\bm{\theta})$-region ($\frac{1}{c_{\max}}$) versus the parameter $p_{2}$ for the particular case of $p_{3}=-4$. Here $\frac{1}{c_{\max}}=0$ indicates that the system is globally stable and $\frac{1}{c_{\max}}=\infty$ indicates that the system is unstable. Intermediate values of $\frac{1}{c_{\max}}$ indicate that the system is $c$-radius stable. In Fig. 7(C), we present a box plot of the training error ($\Delta_{RC}$) versus the parameter $p_{2}$ ($p_{3}=-4$). In this simulation, we consider 100 different realizations of the matrix $A$. The training signal $s(t)$ is the $x$-component and the input signal $g(t)$ is the $z$-component of the Lorenz attractor.

We run some additional simulations to investigate the importance of carefully choosing the nonzero coefficients $p_{i}$ of the polynomial (53) on the reservoir performance. These are shown in Figs. 5 and 6. The case shown in Fig. 5 is that of a polynomial with linear, cubic and quartic order terms (but no quadratic term.) Similarly to the case when linear, quadratic and cubic terms were present, we do not see the distinct boundary between sets of parameters such that the RC is globally stable (below the black line in . 3) versus sets of parameters such that $c_{\max}$ is finite (above the black curve in Fig. 3). In Fig. 5 we plot the training error $\Delta_{RC}$ as a function of the parameters $p_{3}$ and $p_{4}$. The color represents variations in log scale of the training error from dark blue (small) to dark red (large). The solid black curves represent the level curves for different values of $c_{\max}$ computed from $K^{*}(c_{\max},\bm{\theta})=-\alpha_{\max}(A_{s})$ in the parameter space $(p_{3},p_{4})$. We observe that the reservoir computer performs well for those values of $p_{3}$ and $p_{4}$ for which $c_{\max}=1$, except in the case in which $p_{4}=0$, which is analogous to what previously shown in Fig. 3 for $p_{2}=0$. Moreover, we see for a case in which the polynomial had first order, third order, and fifth order terms but no even power terms that the training error was always very high (shown in Fig. 6). We wish to emphasize that while global stability could be achieved by a proper choice of the parameter $p_{1}$ alone (with all the other $p_{i}$, $i>1$, equal zero), our simulations show the importance of setting certain coefficients $p_{i}$, $i>1$, in the polynomial (53) to be non-zero and in particular it appears to be important to have nonzero coefficients corresponding to at least one odd power and one even power.

III.3 Results for Discrete-Time Dynamics with Sigmoid Nonlinearity

For the numerical simulations of a reservoir computer with discrete-time dynamics, we set the nodal dynamics as in Eq.  (15) with $f(r_{i},\bm{\theta})=\frac{p_{1}}{1+e^{-p_{2}r_{i}}}$. We choose the matrix $A$ as we described at the beginning of the section. We set the input signal and training signal from the Lorenz chaotic attractor and compute the training error as described in Sec. III.A. We consider 100 realizations of the matrix $A$ but keep the input and the training signal unchanged. We compute the fixed point $q_{i}^{*}$ of Eq. (43) and the scalar functions $K^{-*}(c,\theta)$ and $K^{+*}(c,\theta)$ by following the Eqs. (44)-(47). In Fig. 8, we plot the training error $\Delta_{RC}$ in log scale as a function of the parameters $p_{1}$ and $p_{2}$. The solid black curve represents ${K^{+}}^{*}(c\rightarrow\infty,\bm{\theta})=\rho^{+}_{c}$ as a function of $p_{1}$ and $p_{2}$. The dashed black curve represents ${K^{-}}^{*}(c\rightarrow\infty,\bm{\theta})=\rho^{-}_{c}$ as a function of $p_{1}$ and $p_{2}$. The region between the two black curves determines the part of the parameter space $(p_{1},p_{2})$ for which the reservoir computer is globally stable and the reservoir computer is successful in performing the computation, while the training error may vary by different orders of magnitude as the parameter changes. On the other hand outside of this region, it is hard to assess how the system behaves. For example, the system could still be globally stable, or locally stable, in which case depending on the particular choice of input, the system dynamics might approach a different attractor or might be driven to $\pm\infty$. Figure 9 displays the results in Fig. 8 in further detail for the cross-section $p_{2}=0.5$. Figure 9(A) shows ${K^{+}}^{*}(c\rightarrow\infty,\bm{\theta})$ (magenta) and ${K^{-}}^{*}(c\rightarrow\infty,\bm{\theta})$ (green) as functions of $p_{1}$ for constant $p_{2}=0.5$. The black solid line is the constant-ordinate line at $\rho_{c}^{+}$ and the dashed black line is the constant-ordinate line at $\rho_{c}^{-}$. For $-4\leq p_{1}\leq 4$, the system is globally stable and the reservoir computer successfully performs the computation. In Fig. 9(B), we plot the training error ($\Delta_{RC}$) versus the parameter $p1$ for the particular case of $p_{2}=-0.5$. We notice that when $-4<p_{1}<0$, the training error is a bit high but the performance of the RC is consistence, and when $0<p_{1}<4$ the RC performs very well.

IV Conclusion and Discussion

In this paper we have used the Lyapunov design method to analyze the nonlinear stability of a generic reservoir computer for both continuous-time and discrete-time dynamics. Our analysis presented in the paper is simple, input independent and yet is able to predict with a certain efficacy the actual performance of the reservoir in terms of computed training error, see e.g. Fig. 3 and 4. Making the analysis input dependent presents the major drawback that the input may not be known a priori and one usually wants the reservoir to be able to process different input signals. Our approach allows computation of the $c({\bm{\theta}})$-radial region of stability about a desired fixed point, where $c$ is the radius of the stability region and ${\bm{\theta}}$ is a set of parameters for the nodes’ individual dynamics. For each $c$ our approach decouples the effects of the individual nodal dynamics from those of the network topology.

For the case of continuous-time dynamics, we have considered a general form of polynomial nonlinearity. We have derived a scalar function $K^{*}(c,{\bm{\theta}})$ which determines the region within the parameter space for which the system is globally stable and the reservoir performance is typically enhanced. Moreover, we have found that the particular type of nonlinearity matters. It is known from the literature that a RC requires nonlinearity to perform well, see e.g., jaeger2004harnessing . Here we have found additional evidence that (i) the performance of a reservoir typically worsens when one of the $p_{i}$ coefficients in the polynomial is set to zero and (ii) it is usually important to have nonzero coefficients corresponding to at least one odd power (e.g., linear) and one even power (e.g., quadratic or quartic order.) These observations hold for both cases that the input and training signals are generated by the Lorenz and the Duffing attractors.

For the case of discrete-time dynamics, a sigmoid function is used for the nodal dynamics. In this case, two scalar functions $K^{-*}(c,{\bm{\theta}})$ and $K^{+*}(c,{\bm{\theta}})$ determine the region on the parameter space for which the system is globally stable and the reservoir performance is enhanced.

Our plots in Figs. 3, 4, and 8 show a remarkable connection between the area of global stability in the parameter space predicted by the analysis and the performance of the RC as measured by the training error. In particular it appears that the training error worsens considerably when a change of the parameters causes loss of global stability. Our nonlinear stability analysis unveils a trade-off between the need for global stability, which is achievable by linear dynamics alone and the need for higher-order terms in the dynamics, which could in turn compromise stability. While fundamental insight into the exact role the nodal dynamics and the adjacency matrix have on the performance of the reservoir is an open area of research, the manual adjustment of the parameters of the reservoir computer is important to determine its dynamic regime. Our nonlinear stability analysis allows us to find the region within the parameter space for which satisfactory sets of parameters may be selected. One is able to then perform a brute search from within this region for adequate sets of parameters. Moreover, by reducing the parameter space to only a finite region, this analysis empowers us to design an optimization problem to find the optimal set of parameters to maximize the performance of a reservoir computer.

Acknowledgments

This work was supported by the National Science Foundation through Grant No. 1727948, the Office of Naval Research through Grant No. N00014-16-1-2637, and the Defense Threat Reduction Agency through Grant No. HDTRA1-12-1-0020.