Seeing double with a multifunctional reservoir computer

Flynn, Andrew; Tsachouridis, Vassilios A.; Amann, Andreas

doi:10.1063/5.0157648

Multifunctional biological neural networks exploit multistability in order to perform multiple tasks without changing any network properties. Enabling artificial neural networks (ANNs) to obtain certain multistabilities in order to perform several tasks, where each task is related to a particular attractor in the network’s state space, naturally has many benefits from a machine learning perspective. Given the association to multistability, in this paper, we explore how the relationship between different attractors influences the ability of a reservoir computer (RC), which is a dynamical system in the form of an ANN, to achieve multifunctionality. We construct the “seeing double” problem in order to systematically study how a RC reconstructs a coexistence of attractors when there is an overlap between them. As the amount of overlap increases, we discover that for multifunctionality to occur, there is a critical dependence on a suitable choice of the spectral radius for the RC’s internal network connections. A bifurcation analysis reveals how multifunctionality emerges and is destroyed as the RC enters a chaotic regime that can lead to chaotic itinerancy.

Multifunctionality is an emerging area of interest in reservoir computing that opens the door to many new machine learning application areas. Multifunctionality describes the ability of a neural network to perform multiple tasks without changing any network properties. However, there is still relatively little known about how multifunctional reservoir computers (RCs) perform certain tasks, such as reconstructing coexistences of attractors. Questions regarding the relationship between the attractors and how many attractors can coexist become important. In this paper, we shed some light on these issues by exploring the dynamics of a RC when trained to solve the “seeing double” problem. To be more specific, this problem involves training a RC to reconstruct a coexistence of a pair of two-dimensional circular orbits that rotate in opposite directions and can be moved closer together and overlap in state space. We find that the range of RC training parameters where multifunctionality is achieved greatly diminishes as the amount of overlap between the orbits increases. In the extreme case when the orbits are completely overlapping, we examine the nature of the basins of attraction and explore the bifurcation structure of the trained RC. By virtue of its simplicity, the seeing double problem enables us to study the fundamentals and improve the formalism behind how a RC achieves multifunctionality.

I. INTRODUCTION

A “reservoir computer” (RC) is a dynamical system that can be realized as an artificial neural network (ANN) and trained to solve certain machine learning (ML) problems. As outlined in Ref. 1, Verstraeten et al.² coined the term RC to unify the following shared philosophies of Jaeger³ and Maass et al.⁴: instead of training all the weights in a network it is sufficient to optimize only the weights of a readout layer to solve a given task. This ideological shift in training ANNs stems from the design of a RC’s internal layer, known as the “reservoir,” whose role is to enable the RC to respond to a sequence of driving inputs that are related to a particular problem. Depending on the nature of the problem, before the RC is trained, it is described as an open-loop (driven) system and its response to a sequence of driving inputs is used to train the RC by producing a readout layer that can, for instance, transform the open-loop RC to a closed-loop (autonomous) system or enable the open-loop RC to generate a desired response to a different input. Using this approach, RCs have been trained to perform tasks such as, decoding neuronal information in brain–machine interfaces,⁵ model-free control of nonlinear systems,^6,7 predicting critical transitions,^8,9 and attractor reconstruction.^10,11

Within the context of attractor reconstruction, which involves training an open-loop RC to become a closed-loop RC that reconstructs the long-term dynamics of a given attractor, the following question naturally emerges, what is the “reconstruction capacity” of a given RC design? In other words, can a RC be trained to have more than one function and become multifunctional? If so, how many attractors can a given RC reconstruct simultaneously and what are the limiting factors? It is questions like these that motivated the translation of “multifunctionality” from biological neural networks (BNNs) to ANNs using a RC.¹² Here it was shown that, in a similar fashion to the operational principles of multifunctional neural networks in the brains of humans and other animals, a RC can also be trained to harness its inherent ability to exhibit multistable dynamics and achieve multifunctionality by performing multiple tasks without needing to change any network properties to do so. It is important to make clear that, in the present paper, multifunctionality is regarded as a feature of only the trained closed-loop multistable RC.

While it has been demonstrated that multifunctional RCs can possess exotic multistabilities like, for example, reconstructing coexistences of chaotic attractors from different dynamical systems,^12–15 several issues arise that do not ordinarily occur in the case of reconstructing only a single attractor. For instance, as shown in Flynn et al.,¹⁴ when training different types of RCs to reconstruct a coexistence of the Lorenz and Halvorsen chaotic attractors, the RCs fail to achieve multifunctionality when the training data describing these attractors are too close together and begin to overlap, i.e., share common regions of state space. However, there is a limit to the level of insight about the relationship between multifunctionality and overlapping training data that can be gained from the above results given this arbitrary choice of chaotic attractors. This subsequently motivates the need for a reductionist approach to study this issue of overlap in a systematic manner.

In this paper, we introduce the “seeing double” problem as a means to explore the dynamics of a RC when trained to achieve multifunctionality in a paradigmatic case of overlapping training data. The seeing double problem involves training a RC to reconstruct a coexistence of two circular orbits that rotate in opposite directions. While these are relatively simple dynamical objects for a RC to reconstruct, we find that the closer together these orbits are the more problematic it becomes for the closed-loop RC to achieve multifunctionality.

In Sec. II, we further discuss the basics of multifunctionality and outline the procedure we use to train a RC to become multifunctional. In Sec. III, we describe how the training data are generated. The results of our numerical experiments are detailed in Secs. IV and V. Here we show that as the orbits are moved closer together, the range of spectral radii of the RC’s internal layer, $ρ$ ⁠, where multifunctionality is achieved greatly diminishes. Despite these difficulties, we find that for a small range of $ρ$ values the closed-loop RC achieves multifunctionality when the orbits are completely overlapping. We examine the nature of the trained closed-loop RC’s basins of attraction in this extreme case and unearth the existence of several “untrained attractors,” attractors which exist in the trained closed-loop RC’s state space but were not present during the training. We explore the interplay between these untrained attractors and the reconstructed orbits and identify the particular bifurcations responsible for the rise and fall of multifunctionality. The relationship between symmetry and the seeing double problem is discussed in Appendix B. Several other dynamical phenomena are also encountered during our analysis, such as routes to chaos and evidence of chaotic itinerancy. We provide some concluding remarks in Sec. VI.

II. MULTIFUNCTIONAL RESERVOIR COMPUTERS

A. Introduction to multifunctionality

Multifunctionality is the term used to describe neural networks that are capable of changing their dynamics on demand of a given duty without needing to alter their synaptic properties. This is a fundamental feature of many neural architectures and is considered to be key to the survival of certain species over time. Multifunctionality is typically found in BNNs with a small number of neurons which are used to perform a set of mutually exclusive tasks. For example, in switching between swimming and crawling motions¹⁶ or regular breathing, sighing, and gasping.^17,18 Multifunctionality has been an active area of research in neuroscience since the mid-1980s with seminal work published by Mpitsos and Cohan¹⁹ and Getting,²⁰ followed by review papers by Dickinson²¹ and Marder and Calabrese²² and more recently reviewed by Briggman and Kristan.²³

At the core of the above examples is a BNN that is capable of changing its activity patterns based on a given input without altering any synapses. From the perspective of dynamical systems, we can say that a multifunctional neural network possesses a multistability. Each task that the network performs there is, in this sense, an attractor associated with it. This attractor is in coexistence with several other attractors in the network’s state space, each of which are distinctly related to the tasks that the network performs.

This multistability interpretation of events dates back to the seminal work of Mpitsos and Cohan,¹⁹ where the researchers use the term multifunctionality to describe networks in which neurons have multiple dynamical regimes but without altering the strength of the synapses. Through Definition 1 in Sec. II C we provide a more precise definition of multifunctionality in the context of reservoir computing.

It is natural to consider that many more network systems may have the ability to achieve multifunctionality. Where this becomes immediately relevant is in the domain of ML. By harnessing multifunctionality from a ML perspective, it can be used to unlock additional computational capabilities of certain ANNs which would otherwise have remained dormant.

Moreover, training ANNs to achieve multifunctionality is highly advantageous from a practicality point of view as it has the potential to broaden the network’s functional capacity to solve multiple tasks using the same set of trained weights. As multifunctionality brings multistability into the world of ML, this, in turn, opens the door for many new application areas.

The development of multifunctional RCs came as a result of recognizing that since a multifunctional neural network in principle resembles a system with a coexistence of attractors, and that a RC can be trained to reconstruct the dynamics of an attractor, then by extending this formalism a RC was trained to exhibit multifunctionality by reconstructing a coexistence of attractors.¹²

B. Reservoir computing

The RC that is studied throughout this paper was presented by Lu et al.¹¹ To outline the steps in how this RC is trained to reconstruct a coexistence of attractors and achieve multifunctionality, let us first consider the case of training this RC to reconstruct the dynamics of a single attractor, $A \subset R^{D}$ ⁠, given access to a trajectory on $A$ described by a vector $u (t) \in A$ ⁠.

There are two stages involved in training this RC, a listening stage and a training stage and we outline the specifics of these stages after introducing the RC. In both stages, the open-loop RC, defined as the nonautonomous dynamical system in Eq. (referred to as the open-loop system in Sec. I), is driven by

u (t)

for

0 \leq t \leq t_{train}

and its response to this driving input is found by generating solutions of the following:

\begin{aligned} \dot{r} (t) & = γ [- r (t) + \tanh (M r (t) + σ W_{i n} u (t))], \end{aligned}

(1)

\begin{aligned} r (0) & = 0^{T} . \end{aligned}

(2)

In Eq. ,

r (t) \in R^{N}

describes the state of the open-loop RC at a given time

t

and

N

is the number of artificial neurons in the network. Solutions of Eq. are computed using the fourth order Runge–Kutta method with time step

τ

⁠.

γ

is a decay-rate parameter arising from the derivation of this RC from the initial discrete-time design proposed by Jaeger.³ The

\tanh

“activation function” is a pointwise operation and is defined as

\tanh (\cdot) : R^{N} \to R^{N}

⁠. The topology of the network is described by the adjacency matrix,

M \in R^{N \times N}

⁠. The input strength parameter,

σ

⁠, and the input matrix,

W_{i n} \in R^{N \times D}

⁠, when multiplied together represent the weight given to

u (t)

⁠, as it is projected into the open-loop RC. The steps involved in constructing

M

and

W_{i n}

are outlined in Appendix A.

In our results, the spectral radius of $M$ ⁠, which is denoted by $ρ$ ⁠, is shown to play a crucial role in training this RC design to solve the seeing double problem and achieve multifunctionality. $ρ$ is tuned by rescaling the elements of $M$ such that the maximum of the absolute values of its eigenvalues is $ρ$ ⁠. $ρ$ has also been a key parameter in our previous results on training a RC to achieve multifunctionality.^12–14

We now describe the specifics of the listening and training stages we mentioned earlier. The listening stage is described as the duration of time where Eq. (1) is driven by $u (t)$ for $0 \leq t \leq t_{listen}$ ⁠, where $t_{listen} < t_{train}$ is chosen such that $r (t)$ is determined by a history of driving inputs and is no longer dependent on its initial condition.

The training stage is the time where Eq. is driven by

u (t)

for

t_{listen} \leq t \leq t_{train}

⁠. In this paper, training consists of finding a “readout function/layer” defined as,

\hat{ψ} (\cdot) : R^{2 N} \to R^{D}

⁠, that replaces

u (t)

in Eq. and is written as

\hat{ψ} (r (t)) = W_{o u t} q (r (t)),

(3)

where

W_{o u t} \in R^{D \times 2 N}

is the “readout matrix” and

q (r (t)) \in R^{2 N}

is given by

q (r (t)) = (\begin{matrix} r (t) \\ r^{2} (t) \end{matrix}),

(4)

where

r^{2} (t) = {(r_{1}^{2} (t), r_{2}^{2} (t), \dots, r_{N}^{2} (t))}^{T}

⁠. From Eq. it then follows that Eq. can be rewritten as,

\hat{ψ} (r (t)) = W_{o u t}^{(1)} r (t) + W_{o u t}^{(2)} r^{2} (t),

(5)

where

W_{o u t}^{(1)}

is the “linear readout matrix,” and

W_{o u t}^{(2)}

is the “square readout matrix” which breaks the symmetry in Eq. when replacing

u (t)

with

\hat{ψ} (r (t))

after the training and prevents the occurrence of “mirror-attractors” which can impede the ability of the RC to reconstruct attractors.^13,24

W_{o u t}

is determined by the ridge regression technique which consists of minimizing the following expression:

\frac{1}{t^{*} - l^{*}} [\sum_{i = l^{*}}^{t^{*}} | | W_{o u t} q (r [i]) - u [i] | |_{2}^{2} + β | | W_{o u t} | |_{2}^{2}],

(6)

with respect to

W_{o u t}

⁠, where

l^{*} = t_{listen} / τ

and

t^{*} = t_{train} / τ

⁠,

r [i]

and

u [i]

are discrete-time samples of the continuous-time variables

r (t)

and

u (t)

at discrete time

i = t / τ

where, in this paper,

τ

is chosen as the time step of the integration. The corresponding time series,

{[r [i]]}_{i = l^{*}}^{t^{*}}

and

{[u [i]]}_{i = l^{*}}^{t^{*}}

⁠, are constructed by sampling

r (t)

and

u (t)

at intervals of length

τ

(discrete-time intervals of length

1

⁠).

β

is the regularization parameter and the purpose of the

β | | W_{o u t} | |_{2}

term in Eq. is to modify the linear least-squares regression to reduce the magnitudes of elements in

W_{o u t}

in order to discourage overfitting.

Minimizing Eq. involves computing its partial derivative with respect to

W_{o u t}

and setting the resulting expression equal to

0

⁠.

W_{o u t}

that minimizes Eq. is given by

W_{o u t} = Y X^{T} {(X X^{T} + β I)}^{- 1},

(7)

where

X = [\begin{array}{cccc} q (r [l^{*}]) & q (r [l^{*} + 1]) & \dots & q (r [t^{*}]) \end{array}],

(8)

is the “response data matrix” and,

Y = [\begin{array}{cccc} u [l^{*}] & u [l^{*} + 1] & \dots & u [t^{*}] \end{array}],

(9)

is the “input data matrix” and

I

is the identity matrix.

After the training,

u (t)

in Eq. is replaced by

\hat{ψ} (r (t))

and in Eq. we now define the closed-loop RC (referred to as the trained closed-loop system in Sec. I) as the following autonomous dynamical system:

\begin{aligned} \dot{\hat{r}} (t) & = γ [- \hat{r} (t) + \tanh (M \hat{r} (t) + σ W_{i n} W_{o u t} q (\hat{r} (t)))], \end{aligned}

(10)

\begin{aligned} \hat{r} (0) & = r (t_{t r a i n}), \end{aligned}

(11)

where

\hat{r} (t)

denotes the state of the closed-loop RC at a given time

t

⁠. While

\hat{r} (t)

and

r (t)

are both

N

-dimensional vectors, to distinguish the dynamics of

\hat{r} (t)

from

r (t)

⁠,

\hat{r} (t)

is defined as

\hat{r} (t) \in S

where

S

is referred to as the “RC’s state space” and is used henceforth when discussing the dynamics of the closed-loop RC. By computing solutions of Eq. , predictions of

u (t)

for

t > t_{train}

⁠, denoted as

\hat{u} (t)

⁠, are given by

\hat{u} (t) = \hat{ψ} (\hat{r} (t)) .

(12)

Again, while both

u (t)

and

\hat{u} (t)

are

D

-dimensional vectors, to distinguish the dynamics of

\hat{u} (t)

from

u (t)

⁠,

\hat{u} (t)

is defined as

\hat{u} (t) \in P

where

P

is referred to as the “projected state space” and is used henceforth when discussing the dynamics of the closed-loop RC when projected via

\hat{ψ}

⁠.

We say the closed-loop RC has achieved attractor reconstruction when the long-term dynamical characteristics of $\hat{u} (t)$ are indistinguishable from $u (t)$ ⁠. In this scenario, there exists an attractor $S \subset S$ such that when the state of the closed-loop RC approaches $S$ and is projected from $S$ to $P$ via $\hat{ψ} (\cdot)$ ⁠, the dynamics of the “reconstructed attractor,” $\hat{A} \subset P$ ⁠, will resemble the dynamics of $A$ in the long-term. By resembling the long-term dynamics it is meant that, for instance, $A$ and $\hat{A}$ will have nearly identical Poincaré sections when computed for the same region of $R^{D}$ and $P$ as $t \to \infty$ for $t > t_{train}$ ⁠.

C. Training a RC to become multifunctional

Using the terminology we have established, we now define multifunctionality in the context of having successfully trained the closed-loop RC in Eq. (10) to reconstruct a coexistence of $n$ attractors, $A_{1} \subset R^{D}$ ⁠, $A_{2} \subset R^{D}$ ⁠, $\dots$ ⁠, $A_{n} \subset R^{D}$ ⁠, given input training data from each of, $u_{(A_{1})} (t) \in A_{1}$ ⁠, $u_{(A_{2})} (t) \in A_{2}$ ⁠, $\dots$ ⁠, $u_{(A_{n})} (t) \in A_{n}$ ⁠, for $0 < t \leq t_{train}$ ⁠.

Definition 1

The closed-loop RC, defined in Eq. (10), achieves multifunctionality if for every $A_{j}$ ⁠, for $j = 1, \dots, n$ ⁠, there exists a corresponding attractor $S_{j} \subset S$ such that the projection of each $S_{j}$ from $S$ to $P$ via $\hat{ψ}$ resembles the long-term dynamics of the respective $A_{j}$ ⁠.

The projection of a given $S_{j}$ from $S$ to $P$ is also referred to as the reconstructed attractor ${\hat{A}}_{j} \subset P$ ⁠. What Definition 1 says it that in order for a RC to achieve multifunctionality it is necessary for it to inherit the ability to perform multiple tasks using the same readout matrix, $W_{o u t}$ ⁠, and without changing any other structural properties of the RC.

The steps involved in computing $W_{o u t}$ from Sec. II B are adapted for the case of multifunctionality. For simplicity, let us first consider the case of reconstructing a coexistence of two attractors, $A_{1} \subset R^{D}$ and $A_{2} \subset R^{D}$ ⁠.

We first drive the open-loop RC in Eq. with input

u_{(A_{1})} (t) \in A_{1}

for

0 < t \leq t_{train}

and then repeat for

u_{(A_{2})} (t) \in A_{2}

⁠. The responses of the open-loop RC to these driving inputs are denoted by

r_{(A_{1})} (t)

and

r_{(A_{2})} (t)

⁠. It is important to highlight that

M

⁠,

W_{i n}

⁠, and all training parameters remain identical when generating

r_{(A_{1})} (t)

and

r_{(A_{2})} (t)

in order to be consistent with the descriptions of multifunctionality from the biological perspective and the reservoir computing analogue in Sec. II A. The same readout function design as in Eq. is used and the ridge regression approach for computing

W_{o u t} \in R^{D \times 2 N}

in Eq. is modified to account for additional attractors. This modified ridge regression approach consists of minimizing the following expression with respect to

W_{o u t}

⁠:

\begin{aligned} \frac{1}{t^{*} - l^{*}} [\sum_{i = l^{*}}^{t^{*}} | | W_{o u t} q (r_{(A_{1})} [i]) - u_{(A_{1})} [i] | |_{2}^{2} \\ + \sum_{i = l^{*}}^{t^{*}} | | W_{o u t} q (r_{(A_{2})} [i]) - u_{(A_{2})} [i] | |_{2}^{2} + β | | W_{o u t} | |_{2}^{2}], \end{aligned}

(13)

where

r_{(A_{1})} [i]

⁠,

u_{(A_{1})} [i]

⁠,

r_{(A_{2})} [i]

⁠, and

u_{(A_{2})} [i]

are discrete-time samples of

r_{(A_{1})} (t)

⁠,

u_{(A_{1})} (t)

⁠,

r_{(A_{2})} (t)

⁠, and

u_{(A_{2})} (t)

and are constructed using the same convention outlined in Sec. II B. The regularization parameter

β

plays the same role as discussed in Sec. II B.

W_{o u t}

that minimizes Eq. is

W_{o u t} = Y_{C} X_{C}^{T} {(X_{C} X_{C}^{T} + β I)}^{- 1},

(14)

where

X_{C} = [X_{(A_{1})}, X_{(A_{2})}]

⁠,

Y_{C} = [Y_{(A_{1})}, Y_{(A_{2})}]

with

X_{(A_{1})}

and

X_{(A_{2})}

constructed as in Eq. for the corresponding

r_{(A_{1})} [i]

and

r_{(A_{2})} [i]

and similarly for

Y_{(A_{1})}

as in Eq. and

Y_{(A_{2})}

for the associated

u_{(A_{1})} [i]

and

u_{(A_{2})} [i]

⁠. In Eq. ,

I

is the identity matrix of the appropriate size.

The closed-loop RC is then described by the same autonomous continuous-time dynamical system as in Eq. (10) using $W_{o u t}$ that is obtained from Eq. (14). If multifunctionality is achieved, then we refer to the resulting multistable closed-loop RC as the “multifunctional RC.”

As per Definition 1, for each $A_{1}$ and $A_{2}$ ⁠, there exist two coexisting attractors $S_{1}, S_{2} \subset S$ such that when the state of the multifunctional RC approaches either $S_{1}$ or $S_{2}$ and is projected from $S$ to $P$ using $\hat{ψ} (\cdot)$ as in Eq. (3), then the long-term dynamics of the reconstructed attractors, ${\hat{A}}_{1}$ and ${\hat{A}}_{2}$ ⁠, will resemble the long-term dynamics of $A_{1}$ and $A_{2}$ ⁠.

Therefore, to reconstruct the dynamics of either $A_{1}$ or $A_{2}$ ⁠, the multifunctional RC needs to be initialized with the corresponding $\hat{r} (0) = r_{(A_{1})} (t_{t r a i n})$ or $r_{(A_{2})} (t_{t r a i n})$ ⁠, or some other initial condition that is in the basin of attraction of the corresponding $S_{1}$ or $S_{2}$ ⁠.

The same steps as outlined in Sec. II C can be extended in order to produce a multifunctional RC that reconstructs a coexistence of $n$ attractors, $A_{1}, A_{2}, \dots, A_{n}$ ⁠. All that is required is to produce the corresponding $X_{C} = [X_{(A_{1})}, X_{(A_{2})}, \dots, X_{(A_{n})}]$ and $Y_{C} = [Y_{(A_{1})}, Y_{(A_{2})}, \dots, Y_{(A_{n})}]$ and solve for $W_{o u t}$ as in Eq. (14).

D. Summary of previous work

As mentioned previously, multifunctionality was first translated from BNNs to ANNs in Flynn et al.¹² Here, it was shown that the same RC design as in Eq. (1) can be trained to reconstruct a coexistence of chaotic attractors from a multistable system, different parameter settings of the same system, and two systems with different nonlinearities. The issues of choosing appropriate training parameters, differences in the timescales of the chaotic attractor’s dynamics, and the existence of untrained attractors are also investigated. Independent of these results, Herteux and Ráth²⁴ also demonstrated how to train a RC to reconstruct a coexistence of chaotic attractors using the same training technique described in Sec. II C. The phenomena which arise when training Eq. (1) to reconstruct a coexistence of attractors that are related by certain symmetries are explored in Flynn et al.,¹³ the analytical and numerical results presented here show how a RC can adapt to and take on the symmetry of the trained system spontaneously. In Flynn et al.¹⁴, the ability of different types of RCs to achieve multifunctionality in the presence of overlapping training data is assessed. The present paper extends on the above results in many ways. Moreover, it is thanks to the seeing double problem’s simplicity that we are able to more rigorously examine the dynamics and improve the formalism behind how a RC achieves multifunctionality.

Since the first results regarding multifunctionality in a RC,¹² it has been discovered that a given $W_{o u t}$ can enable a RC to perform more than one task using different principles to those involving multifunctionality. In this context, the so-called “parameter-aware” RCs have been trained to interpolate and extrapolate the dynamics of a system such as those studied in Kim et al.,²⁵ Kong et al.,^26,27 and Fan et al.²⁸ To do this, an additional parameter input channel is incorporated into the RC’s architecture during training so that, in this case, the closed-loop RC performs a task associated with a given parameter value and $W_{o u t}$ remains fixed after the training. In these examples of parameter-aware RCs, it is shown that the RC can go through similar bifurcations to the attractor from the dynamical system used to generate the training data without necessarily needing to be trained on data after bifurcation takes place. In future work, it would be interesting to investigate the dynamical principles behind how this is achieved and how this relates to the phenomenon of untrained attractors as presented in Flynn et al.¹² These attractors form a key aspect of our analysis in Sec. V. Furthermore, while multifunctionality has so far not been studied using these parameter-aware RCs, we comment that a “multifunctional parameter-aware” RC may further enhance the reconstruction capacity of a given RC.

III. SEEING DOUBLE

The “seeing double” problem is introduced in this section. This numerical experiment is designed as a means to systematically induce the issues related to multifunctionality and overlapping training data.

This task involves training a RC to reconstruct a coexistence of attractors that follow trajectories on two circular orbits of equal radius, rotate in opposite directions, and the centers of these orbits can be moved either closer together or further apart. When these orbits are overlapping, the RC is therefore required to distinguish between points in $R^{D}$ that are common to both cycles in order to exhibit multifunctionality.

By virtue of the seeing double problem’s simplicity, a greater emphasis can be placed on examining how multifunctionality arises in a RC. At the same time, despite this simplicity, an enormity of elaborate dynamics are encountered and studied in Secs. IV and V.

A. Numerical experiment setup

The training data are computed via

u (t) = (\begin{matrix} x (t) \\ y (t) \end{matrix}) = (\begin{matrix} b_{x} \cos (t) + x_{c e n} \\ b_{y} \sin (t) + y_{c e n} \end{matrix}),

(15)

for

t = 0, τ, 2 τ, \dots,

using the time step

τ = 0.01

⁠. The resultant

u (t)

generated from Eq. resembles a trajectory around a circle of radius

b = | b_{x} | = | b_{y} |

and centered at

(x_{c e n}, y_{c e n})

⁠. The sets,

C_{A}

and

C_{B}

⁠, are produced using Eq. for specific values of

b_{x}, b_{y}, x_{c e n},

and

y_{c e n}

⁠. The corresponding input signals that are used to drive the open-loop RC in Eq. are denoted by

u_{(C_{A})} (t)

and

u_{(C_{B})} (t)

⁠. The open-loop RC’s response to

u_{(C_{A})} (t)

and

u_{(C_{B})} (t)

are denoted as

r_{(C_{A})} (t)

and

r_{(C_{B})} (t)

for

0 \leq t \leq t_{t r a i n}

⁠. The training procedure outlined in Sec. II C is used to produce the corresponding

X_{(C_{A})}, X_{(C_{B})}, Y_{(C_{A})},

and

Y_{(C_{B})}

that are used to compute

W_{o u t}

in Eq. . This

W_{o u t}

is then used in the closed-loop RC setup given by Eq. .

We say that this closed-loop RC achieves multifunctionality and solves the seeing double problem once it reconstructs a coexistence of two attractors, $S_{A}$ and $S_{B}$ ⁠, that exist in $S$ and resemble $C_{A}$ and $C_{B}$ when projected to $P$ using $\hat{ψ} (\cdot)$ as in Eq. (3) with $W_{o u t}$ as mentioned above. As per the same convention used earlier, these projected dynamics of $S_{A}$ and $S_{B}$ are referred to as the reconstructed attractors and denoted by ${\hat{C}}_{A}$ and ${\hat{C}}_{B}$ ⁠. To reconstruct the dynamics of $C_{A}$ or $C_{B}$ using this multifunctional RC, we initialize Eq. (10) with $\hat{r} (0) = {\hat{r}}_{(C_{A})} (0) = r_{(C_{A})} (t_{listen})$ or $\hat{r} (0) = {\hat{r}}_{(C_{B})} (0) = r_{(C_{B})} (t_{listen})$ or some known point in the basin of attraction of $S_{A}$ or $S_{B}$ ⁠. The subsequent states of Eq. (10) when approaching $S_{A}, S_{B} \subset S$ (⁠ ${\hat{C}}_{A}, {\hat{C}}_{B} \subset P$ ⁠) are written as ${\hat{r}}_{(C_{A})} (t)$ and ${\hat{r}}_{(C_{B})} (t)$ ⁠.

In this paper, $b_{x}$ and $b_{y}$ are set to $b_{x} = b_{y} = 5$ to create $C_{A}$ and for $C_{B}$ ⁠, $b_{x} = - 5$ and $b_{y} = 5$ ⁠. The centers of $C_{A}$ and $C_{B}$ are brought closer together by changing the values of $(x_{c e n}, y_{c e n})$ ⁠. In order to simplify the experiment further, $C_{A}$ and $C_{B}$ are designed to be centered at $(x_{c e n}, 0)$ and $(- x_{c e n}, 0)$ ⁠, respectively. Therefore, by changing $x_{c e n}$ the centers of these cycles are moved equidistantly along the line $y = 0$ ⁠. In this case, an overlapping region between $C_{A}$ and $C_{B}$ exists whenever $| x_{c e n} | \leq b = 5$ ⁠, i.e., $C_{A} \cap C_{B} \neq \emptyset$ $\forall | x_{c e n} | \leq 5$ ⁠. Furthermore, $C_{A}$ and $C_{B}$ are said to be “entirely/completely overlapping” when $x_{c e n} = 0$ ⁠. In this extreme case, the only difference between $C_{A}$ and $C_{B}$ is the direction of rotation on both cycles.

B. Illustrating issues of overlap

The orbits $C_{A}$ and $C_{B}$ produced from Eq. (15) are depicted in Fig. 1 with $u_{(C_{A})}$ and $u_{(C_{B})}$ plotted in the $(x, y)$ -plane for a given $x_{c e n}$ ⁠. The green arrows indicate the direction of rotation in which an orbit along these circles evolves upon in time.

FIG. 1.

Illustrating how the centers of C A and C B move in R D for changes in x c e n. In (a) C A ∩ C B = ∅ for x c e n = − 5.5, and in (b) C A ∩ C B ≠ ∅ when x c e n = − 3. Green arrows indicate a trajectory on these cycles. Red arrows indicate a potential path to failure.

View large Download slide

Illustrating how the centers of $C_{A}$ and $C_{B}$ move in $R^{D}$ for changes in $x_{c e n}$ ⁠. In (a) $C_{A} \cap C_{B} = \emptyset$ for $x_{c e n} = - 5.5$ ⁠, and in (b) $C_{A} \cap C_{B} \neq \emptyset$ when $x_{c e n} = - 3$ ⁠. Green arrows indicate a trajectory on these cycles. Red arrows indicate a potential path to failure.

Figure 1(a) shows that the cycles do not overlap when $x_{c e n} = - 5.5$ ⁠. However, when $| x_{c e n} | \leq 5$ ⁠, like in Fig. 1(b) for $x_{c e n} = - 3$ ⁠, the orbits share similar regions of the $(x, y)$ -space. This presents a challenge for the resulting closed-loop RC to remain on the correct path when its state approaches the point where these cycles intersect as indicated by the green arrows at the crossing point in Fig. 1(b). If the RC fails to distinguish between these overlapping cycles then, for example, the state of the RC may move onto a different path [as indicated by the red arrow in Fig. 1(b)] or potentially result in an attractor merging crisis scenario.

Therefore, when there is an overlap between $C_{A}$ and $C_{B}$ ⁠, the state of the open-loop RC must contain sufficient knowledge of its previous states during the training so that after the training the closed-loop RC remains on the correct path when approaching a crossing between $C_{A}$ and $C_{B}$ ⁠. In Sec. IV it is found that by changing the spectral radius, $ρ$ ⁠, of the RC’s internal connections, which as a result tunes the weight given to the RC’s previous states, it is possible to overcome the issue of overlap.

IV. SEEING DOUBLE WITH A MULTIFUNCTIONAL RC

The focus of this section is to explore the dynamics of the closed-loop RC in Eq. (10) after it is has trained to solve the seeing double problem. Here we show that as the difficulty of the seeing double problem is varied according to changes in $x_{c e n}$ (the amount of overlap between $C_{A}$ and $C_{B}$ and the nearer they are to one another), the more “memory” is required of the RC and consequently the spectral radius, $ρ$ ⁠, plays a more significant role. In this instance, we do not use the word memory in a quantitative sense but in how it relates to a property of this RC design that, during the training, by increasing $ρ$ up to a certain point the state of the RC depends not only on the current input but also previous inputs. In other words, the RC has a greater ability to recall information from the past.

The results in this section are generated by training the RC using the method outlined in Sec. II C using the input data, $u_{(C_{A})} (t)$ and $u_{(C_{B})} (t)$ constructed as described in Sec. III A. The RC training and design parameters are specified in Table II and the same $M$ and $W_{i n}$ matrices are used to produce the results presented throughout this paper.

While the results that follow are computed using just one random realization of $M$ and $W_{i n}$ ⁠, we do not claim that our results are universal across all possible $M$ and $W_{i n}$ ⁠. At the same time, the analysis tools we use throughout this paper provide significant insight toward how a RC solves the seeing double problem and are applicable to all random realizations of $M$ and $W_{i n}$ ⁠. Where appropriate, we highlight certain features of our results that appear across several random realizations of $M$ and $W_{i n}$ and other features that may be heavily dependent on the particular $M$ and $W_{i n}$ that is used to generate our results.

A. Exploring multifunctionality in (x_cen, ρ)-plane

We now explore the long-term dynamics of the closed-loop RC [described in Eq. (10)] when initialized with $\hat{r} (0) = r_{(C_{A})} (t_{listen})$ and $\hat{r} (0) = r_{(C_{B})} (t_{listen})$ ⁠, as viewed from $P$ ⁠, for a given $x_{c e n} \in [- 10, 10]$ and $ρ \in [0.1, 2.5]$ ⁠.

1. Prediction analysis

After the training, a number of outcomes are possible. If the training is successful, then ${\hat{C}}_{A} \approx C_{A}$ and ${\hat{C}}_{B} \approx C_{B}$ and the closed-loop RC’s predicted trajectory on either $C_{A}$ or $C_{B}$ will remain on ${\hat{C}}_{A}$ or ${\hat{C}}_{B}$ in the long-term (⁠ $t \to \infty$ ⁠). However, given the nature of the seeing double problem, it is also possible for the predicted trajectory to switch from ${\hat{C}}_{A}$ to ${\hat{C}}_{B}$ or vice versa if the RC fails to lift the corresponding driving input data, $u_{(C_{A})} (t)$ and $u_{(C_{B})} (t)$ ⁠, into separate regions of $S$ ⁠. Furthermore, if the RC fails to reconstruct $C_{A}$ or $C_{B}$ then the dynamics of Eq. (10) can, for example, decay to a fixed point, some other limit cycle, or never settle down to periodic behavior and display evidence of the RC entering a chaotic regime.

An error metric known as the “roundness” is used to assess the accuracy of the RC’s reconstruction of a given orbit. The roundness of a given cycle, $C$ ⁠, is denoted by $δ (C)$ and is found by calculating the difference between the radii of the maximum and minimum circles needed to enclose and inscribe the cycle. For simplicity, the centers of these maximum and minimum circles are taken to be located at the center of the given target orbit. It was found empirically that if the relative roundness of a given cycle, $δ_{r e l} (C) = δ (C) / b$ ⁠, was $< δ_{r e l}^{+} = 0.25$ ⁠, then the RC is in general able to successfully reconstruct the desired orbit with a high level of accuracy. This error metric is typically used in the production of circular shaped materials, like plastic tubes or copper pipe, and the formula to calculate the roundness is found in most machining manuals in line with the International Organisation for Standardisation (ISO).

The RC’s reconstruction of a given orbit is characterized as the particular attractor that the state of the closed-loop RC is approaching for $t_{predict} - t^{*} \leq t \leq t_{predict}$ where $t_{predict} = 600$ and $t^{*} = 40$ as viewed from $P$ ⁠. The result of this is displayed in Figs. 2(a)–2(b) where each point in the $(x_{c e n}, ρ)$ -plane is colored according to the specifications in Table I.

FIG. 2.

View large Download slide

Characterizing the prediction of (a) $C_{A}$ and (b) $C_{B}$ in the $(x_{c e n}, ρ)$ -plane. Coloring scheme is as outlined in Table I.

TABLE I.

Color scheme for prediction analysis in Fig. 2.

Blue

Correct cycle is reconstructed and

δ_{r e l} (C) < δ_{r e l}^{+}

⁠.

Yellow

Prediction switches to the other cycle.

Magenta

Prediction does not settle down to periodic motion

and label this as some attractor,

S

⁠.

Black

Prediction tends to some other limit cycle

or

δ_{r e l} (C) \geq δ_{r e l}^{+}

⁠, this attractor is labeled as, LC.

Green

Prediction decays toward some fixed point labeled as FP.

In order to deduce the regions of the $(x_{c e n}, ρ)$ -plane in which multifunctionality was achieved, an additional picture is provided in Fig. 3 where the common blue regions of Figs. 2(a) and 2(b) are identified and the larger of the $δ_{r e l} (C)$ values when assessing the accuracy of the closed-loop RC’s reconstruction of $C_{A}$ and $C_{B}$ ⁠, defined as $δ_{r e l}^{m a x} = max (δ_{r e l} ({\hat{C}}_{A}), δ_{r e l} ({\hat{C}}_{B}))$ ⁠, is plotted for each point in the $(x_{c e n}, ρ)$ -plane.

FIG. 3.

View large Download slide

Regions of multifunctionality in the $(x_{c e n}, ρ)$ -plane with error $δ_{r e l}^{m a x}$ for $C_{A}$ and $C_{B}$ rotating in opposite directions. Green and cyan curves highlight where $λ_{MAX} > λ_{C}$ for the first time for ${\hat{C}}_{A}$ and ${\hat{C}}_{B}$ ⁠, respectively.

2. Regions of multifunctionality

The relationship between overlapping training data, multifunctionality, and $ρ$ becomes evident in Figs. 2 and 3 where a “Goldilocks effect” is found to occur. By this, it is meant that if $ρ$ is too small or too large then the closed-loop RC fails to achieve multifunctionality. When the orbits are furthest apart, we see that the closed-loop RC is multifunctional even at small values of $ρ$ ⁠. However, as the orbits are moved closer together and begin to overlap (i.e., when $C_{A} \cap C_{B} \neq \emptyset$ ⁠), then the closed-loop RC requires a much larger $ρ$ in order to distinguish between the cycles and prevent the predicted trajectories from decaying to a fixed point or some other limit cycle. At the same time, if $ρ$ is too large, the closed-loop RC fails to reconstruct either orbit regardless the investigated values of $x_{c e n}$ and displays evidence of chaos.

We comment that from further inspection there are some subtle differences in the results for different $M$ and $W_{i n}$ ⁠; however, this Goldilocks effect is a consistent feature.

Furthermore, it is important to note that a near symmetric set of results occur about $x_{c e n} = 0$ in all Figs. 2 and 3, some insight to this is provided in Appendix B 1.

3. $C_{A}$ and $C_{B}$ rotating in the same direction

To provide a comparison, the same RC (constructed via the parameters in Table II) is trained to reconstruct a coexistence of $C_{A}$ and $C_{B}$ except, in this case, both cycles are rotating in the same direction. This is done by setting $b_{x} = b_{y} = 5$ for both $C_{A}$ and $C_{B}$ in Eq. (15). We examine the regions of multifunctionality in the $(x_{c e n}, ρ)$ -plane and these are computed using the same procedure that produced Fig. 3.

TABLE II.

RC design and training parameters used to generate the results displayed in Sec. IV.

N	P	ρ	σ	γ	β	t_listen	t_train
1000	0.04	[0.1, 2.5]	0.2	5	10⁻²	200	400

N	P	ρ	σ	γ	β	t_listen	t_train
1000	0.04	[0.1, 2.5]	0.2	5	10⁻²	200	400

The result of the above is presented in Fig. 4 where a similar phenomenon to Fig. 3 is shown to occur, once the overlap is introduced then multifunctionality is achieved only when $ρ$ is sufficiently large. However, the main difference between Figs. 4 and 3 is the trivial case for $x_{c e n} = 0$ as $C_{A}$ and $C_{B}$ are the exact same and the closed-loop RC has no issue reconstructing the dynamics for $ρ \in [0.1, 1.7]$ ⁠. For $ρ > 2.1$ ⁠, the RC fails to become multifunctional $\forall x_{c e n} \in [- 10, 10]$ ⁠.

FIG. 4.

View large Download slide

Regions of multifunctionality in the $(x_{c e n}, ρ)$ -plane with error $δ_{r e l}^{m a x}$ for $C_{A}$ and $C_{B}$ rotating in the same direction. Green and cyan curves highlight where $λ_{MAX} > λ_{C}$ for the first time for ${\hat{C}}_{A}$ and ${\hat{C}}_{B}$ ⁠, respectively.

B. Transition to chaotic dynamics

1. Dynamics before chaos

In Figs. 3 and 4, we find that for a given value of $x_{c e n}$ and $ρ$ ⁠, the smallest error, i.e., the smallest $δ_{r e l}^{m a x}$ value, occurs in many instances before the closed-loop RC’s behavior becomes chaotic for large $ρ$ ⁠. This transition to chaos is quantified in Figs. 3 and 4 by the green and cyan curves which show the location where the largest Lyapunov exponent (LLE), denoted by $λ_{MAX}$ ⁠, first becomes larger than $λ_{C} = 0.01$ when increasing $ρ$ for a given $x_{c e n}$ ⁠. We do this in order to indicate, within some measurement error, that $λ_{MAX}$ is positive.

Figure 4 shows that this transition to chaotic dynamics appears in the same vicinity of $ρ$ values as it does in Fig. 3 for $| x_{c e n} | ≳ 1.5$ ⁠. As the same RC design is used in both cases, this suggests that the location of this transition may also arise as an intrinsic property of the RC rather than solely being related to the nature of the task itself.

While the above comment is only speculation, some progress has been made in answering questions of this nature. For instance, Jiang and Lai²⁹ conducted many extensive experiments which provide significant insight toward assessing the role of $ρ$ in the performance of a RC when training 100 random realizations of this RC in attractor reconstruction problems. The central result from Jiang and Lai’s experiments is that, for a given RC, there exists an interval of $ρ$ values where optimal or near-optimal performance is achieved. A similar phenomenon is seen to occur in Figs. 3 and 4 which compliments the work of Jiang and Lai and also offers additional information as Figs. 3 and 4 both show that as the difficulty of the task is varied (by changing $x_{c e n}$ ⁠) then the width of these intervals change accordingly.

We remark that in studies involving the training of input-driven RCs [after the training the readout function does not replace the driving input like in Eq. (10)], a similar phenomenon known as the “edge of chaos” or the “edge of stability” and has been the focus of many investigations which suggest that RCs achieve optimal computational capacity at the edge of chaos.^30–32 However, Carroll³³ presents several examples of RCs which perform better on certain problems just prior to crossing the edge of chaos.

2. Dynamics after chaos

In some cases, after the transition to chaotic dynamics, an interesting sequence of events is found to occur when increasing $ρ$ past the point of where multifunctionality is no longer achieved. As an example, Fig. 5 illustrates the changes in the dynamics of Eq. (10) as viewed from $P$ as $ρ$ is increased from $1.864$ to $1.944$ in the case where $x_{c e n} = 6.5$ ⁠. In relation to Figs. 2(a) and 2(b), the long-term dynamics of the closed-loop RC when trained at these $ρ$ values and initialized from ${\hat{r}}_{(C_{A})}$ or ${\hat{r}}_{(C_{B})}$ are color-coded as magenta, indicating that the state of Eq. (10) does not settle down to periodic motion.

FIG. 5.

View large Download slide

Closed-loop RC’s reconstruction of $C_{A}$ and $C_{B}$ when $x_{c e n} = 6.5$ for $ρ = 1.864$ (top left), $1.904$ (top right), $1.944$ (middle left). Middle right: Reconstructed $x$ variable time traces highlight the nature of the RC’s itinerant dynamics for $ρ = 1.944$ ⁠. Distribution of residence times the RC spends in either ${\tilde{C}}_{A}$ (bottom left) or ${\tilde{C}}_{B}$ (bottom right) when computing a trajectory of the closed-loop RC up to $t = 1 000 000$ ⁠.

FIG. 6.

View large Download slide

RC’s reconstruction of $C_{A}$ and $C_{B}$ when $x_{c e n} = 0$ for $ρ = 0.80$ ⁠, $1.10$ ⁠, $1.25$ ⁠, $1.40$ ⁠, and $1.70$ ⁠.

Figure 5 shows that for $ρ = 1.864$ and $1.904$ ⁠, the dynamics of both ${\hat{C}}_{A}$ and ${\hat{C}}_{B}$ begin to occupy larger and larger regions of $P$ ⁠. This behavior continues up to a point where the closed-loop RC loses multistability and the state of Eq. (10) begins to switch from one orbit to the other and back again, shown here for $ρ = 1.944$ ⁠.

While further analysis is required in order to determine the route of these switching dynamics, it is possible that the closed-loop RC has entered into the dynamical regime known as “chaotic itinerancy.” Chaotic itinerancy^34–36 describes the dynamical phenomenon of an autonomous switching process whereby the state of a given autonomous dynamical system switches between several “attractor ruins” or “quasi-attractors.” These were previously coexisting attractors that are now connected to form one global attractor. The quasi-attractors retain much of their original features except trajectories on these quasi-attractors leak into each other. In the case of Fig. 5 for $ρ = 1.944$ ⁠, ${\hat{C}}_{A}$ and ${\hat{C}}_{B}$ exhibit these quasi-attractor properties. The state of the closed-loop RC wanders between visiting different regions of $P$ associated with attractor ruins of the previously coexisting ${\hat{C}}_{A}$ and ${\hat{C}}_{B}$ ⁠. The corresponding time traces of the reconstructed $x$ variable (as shown in the middle right in Fig. 5) provide a clearer picture of the switching dynamics between attractor ruins as emphasised by the arrows in the associated $P$ (middle left panel in Fig. 5). After simulating the dynamics of Eq. (10) up to $t = 1 000 000$ in this metastable regime from $\hat{r} (0) = r_{(C_{A})} (t_{t r a i n})$ ⁠, we find that there are over $200$ switchings between the two quasi-attractors. In the bottom panels, we also plot the distribution of time spent in each quasi-attractor, labeled here as ${\tilde{C}}_{A}$ and ${\tilde{C}}_{B}$ ⁠. These plots show that the state of the closed-loop RC spends a far greater amount of time in ${\tilde{C}}_{A}$ as opposed to ${\tilde{C}}_{B}$ ⁠. In future work, we aim to explore the relationship between $ρ$ and the resulting distribution of residence times in each quasi-attractor.

It is important to stress here that there is no external input used to generate the switching dynamics shown in Fig. 5 as it is an intrinsic property of the closed-loop RC after it loses multifunctionality and subsequently multistability. This differs from previous cases where RCs were trained to switch between learning different chaotic and periodic attractors based on an external input like, for example, in the work of Inoue et al. in Ref. 37 and Lu and Bassett in Ref. 38. In these scenarios, $W_{o u t}$ is updated online in response to a given driving input from an individual attractor and only one attractor is shown to exist at each point in time.

C. Achieving multifunctionality when x_cen = 0

1. Dynamics in $P P$

From a dynamical systems perspective, it is remarkable that the closed-loop RC reconstructs a coexistence of $C_{A}$ and $C_{B}$ when these share similar regions of state space. However, the most striking result from Fig. 3 is that the closed-loop RC reconstructs a coexistence of these orbits even when they are completely overlapping, i.e., when $x_{c e n} = 0$ ⁠. In this extreme scenario, multifunctionality is achieved for a small range of $ρ$ values, specifically for $1.1 \leq ρ \leq 1.4$ ⁠. The most accurate reconstruction of both orbits is achieved when $ρ = 1.25$ ⁠.

Figure 6 highlights how values of $δ_{r e l}^{m a x}$ differ in regions where multifunctionality was achieved, for instance, $δ_{r e l}^{m a x} \approx 0.062, 0.019, 0.125$ for $ρ = 1.1, 1.25, 1.4$ ⁠, respectively. These account for the slightly off-circular reconstructions of $C_{A}$ for $ρ = 1.1$ and $C_{B}$ for $ρ = 1.4$ as shown in Fig. 6.

Examples where the closed-loop RC fails to achieve multifunctionality just outside the small window of $ρ$ values identified in Fig. 3 are also illustrated in Fig. 6 for $ρ = 1.7$ and $0.8$ ⁠. Here the chaotic-like behavior of the closed-loop RC is seen for $ρ = 1.7$ ⁠. What is important to note in this example is that shape swept out by the closed-loop RC’s trajectories when initialized from both $\hat{r} (0) = r_{(C_{A})} (t_{listen})$ and $\hat{r} (0) = r_{(C_{B})} (t_{listen})$ is bounded to regions nearby the original circular orbits, however the direction of rotation is not preserved.

Figure 6 also shows that when the reconstruction of both $C_{A}$ and $C_{B}$ fails for $ρ = 0.8$ then the state of Eq. (10) approaches the same limit cycle which bears no resemblance to either $C_{A}$ and $C_{B}$ ⁠. The role this limit cycle plays in how the closed-loop RC solves the seeing double problem is assessed in Sec. V.

Much of the following results in this paper focus on exploring the dynamics of the closed-loop RC in this “simplest hardest” case, “simplest” in terms of the dynamics being limit cycles and “hardest” as the training data are entirely overlapping. We comment that, from analysis not shown in the present paper, the small range of $ρ$ values where multifunctionality is achieved is relatively consistent across several random realizations of $M$ and $W_{i n}$ ⁠.

2. Dynamics in $S S$

Before discussing these further results, we wish to provide some insight toward the internal dynamics of the multifunctional RC and how the behavior of the individual neurons differ when reconstructing either $C_{A}$ or $C_{B}$ in the case where $x_{c e n} = 0$ ⁠. Figure 7 depicts a representative example of how the dynamics of the same set of neurons can significantly differ in $S$ as the state of the multifunctional RC approaches either $S_{A}$ or $S_{B}$ while the only difference between the projected dynamics (the dynamics of the reconstructed attractors, ${\hat{C}}_{A}$ and ${\hat{C}}_{B}$ ⁠), is the direction of rotation. To place a stronger emphasis on the fact that we are now examining the dynamics of $S_{A}$ and $S_{B}$ ⁠, in this section we refer to the corresponding states of Eq. (10) as ${\hat{r}}_{(S_{A})} (t)$ and ${\hat{r}}_{(S_{B})} (t)$ ⁠. Figure 7 highlights how the state of the multifunctional RC’s neurons $9, 25, 50, 88, 103, 250, 400$ ⁠, and $650$ evolve over time [starting from when Eq. (10) is initialized with $\hat{r} (0) = {\hat{r}}_{(S_{A})} (0) = r_{(C_{A})} (t_{t r a i n})$ or $\hat{r} (0) = {\hat{r}}_{(S_{B})} (0) = r_{(C_{B})} (t_{t r a i n})$ ] for the case of $ρ = 1.25$ whose dynamics in $P$ are shown in Fig. 6. It is found that many neurons behave similarly to each other when the state of the multifunctional RC approaches either $S_{A}$ or $S_{B}$ ⁠, see the dynamics of ${\hat{r}}_{650}$ in Fig. 7 as an example. At the same time, Fig. 7 shows that the behavior of many other neurons may bare no resemblance to one another. While all neurons exhibit periodic dynamics, some neurons oscillate with multiple local maxima, however, what is consistent across all neurons is that they each have the same period of rotation.

FIG. 7.

View large Download slide

Illustrating how the dynamics of individual neurons can differ depending on whether the multifunctional RC approaches either $S_{A}$ or $S_{B}$ in order to reconstruct $C_{A}$ or $C_{B}$ when $x_{c e n} = 0$ for $ρ = 1.25$ ⁠.

In order to provide further insight to the dynamics of the individual neurons, we now continue our study of the multifunctional RC in the above case where $ρ = 1.25$ and $x_{c e n} = 0$ ⁠. We now examine how the dynamics of these neurons differs when the state of the multifunctional RC approaches the same point on each reconstructed attractor in $P$ ⁠, i.e., when $x_{{\hat{C}}_{A}} (t) - x_{{\hat{C}}_{B}} (t) = 0$ where $x_{{\hat{C}}_{A}} (t)$ and $x_{{\hat{C}}_{A}} (t)$ denote the location of the $x$ -variable in $P$ at a given time $t$ [starting from when Eq. (10) is initialized] when reconstructing $C_{A}$ and similarly for $C_{B}$ ⁠. Note that the associated $y$ -variables are always equal for all values of $t$ ⁠, this is why we need only look at the $x$ -variable. Using the same notation convention as for $x_{{\hat{C}}_{A}}$ and $x_{{\hat{C}}_{B}}$ ⁠, in the top panel of Fig. 8 we plot how the difference between ${\hat{r}}_{S_{A}}$ and ${\hat{r}}_{S_{B}}$ evolves over time. To convey how these $N = 1000$ neurons behave, we do this by generating a histogram at each time step of the simulation which describes how many components of ${\hat{r}}_{S_{A}} (t) - {\hat{r}}_{S_{B}} (t)$ have a particular value at a given time $t$ ⁠. These values are confined to be within $- 2$ and $2$ (which is the maximum and minimum allowed values due to the $\tanh$ activation function in the RC design) in intervals of $0.02$ (which is the bin width of the histograms and is chosen to show sufficient details of the dynamics at an appropriate resolution). In the lower panel, we plot the evolution of the corresponding difference between $x_{{\hat{C}}_{A}}$ and $x_{{\hat{C}}_{B}}$ ⁠.

FIG. 8.

View large Download slide

For $x_{c e n} = 0$ and $ρ = 1.25$ ⁠, the top panel illustrates the difference between the state of the multifunctional RC when approaching $S_{A}$ or $S_{B}$ ⁠, the lower panel depicts the difference between the respective $x$ -variables when reconstructing $C_{A}$ or $C_{B}$ ⁠.

From Fig. 8, it is evident that at times when $x_{{\hat{C}}_{A}} (t) - x_{{\hat{C}}_{B}} (t) \approx 0$ ⁠, even though many neurons fire identically the crucial aspect for multifunctionality is that not all neurons fire identically as evidenced by the relatively wide distribution of the values of the components of ${\hat{r}}_{S_{A}} (t) - {\hat{r}}_{S_{B}} (t)$ at these particular times. We also see here that the greater the distance between points on each reconstructed attractor, i.e., when $| x_{{\hat{C}}_{A}} - x_{{\hat{C}}_{B}} |$ grows larger, the greater the difference between ${\hat{r}}_{S_{A}}$ and ${\hat{r}}_{S_{B}}$ (the wider the distribution of the firing patterns).

D. Projected basins of attraction for x_cen = 0

The focus of this section is to provide some insight toward the nature of the basins of attraction created by the closed-loop RC in order to solve the seeing double problem in the extreme case of $x_{c e n} = 0$ ⁠.

1. Method of generating the basins

Considering the dimension of the closed-loop RC (⁠ $N = 1000$ ⁠), to determine the precise location of the basins of attraction in $S$ for both $S_{A}$ and $S_{B}$ is highly expensive. Instead, to get some insight about the structure of these basins, we examine how the state of the closed-loop RC evolves from the open-loop RC’s representation of different points in $R^{D}$ ⁠. From this, it is possible to see how $P$ is split to accommodate a coexistence of ${\hat{C}}_{A}$ and ${\hat{C}}_{B}$ when $C_{A}$ and $C_{B}$ are completely overlapping.

To be more specific, for a given $ρ$ and $x_{c e n} = 0$ ⁠, we drive Eq. (1) with a point in $R^{D}$ from $t = 0$ to $t = t_{listen}$ ⁠. After this time, the state of Eq. (1) converges to a representation of the point in $R^{N}$ ⁠. This response vector is labeled as $r_{D} = r (t_{listen})$ ⁠. The corresponding closed-loop RC is then initialized with $\hat{r} (0) = r_{D}$ and the long-term behavior of Eq. (10) is characterized from this initial condition.

2. Basins for ρ = 1.25

The projected basins of attraction for $ρ = 1.25$ ⁠, which is the case where multifunctionality was achieved with the best accuracy, is shown in Fig. 9. Each point in Fig. 9 describes which attractor Eq. (10) approaches in $P$ for a given point in $R^{D}$ ⁠. If the state of Eq. (10) approaches ${\hat{C}}_{A}$ ⁠, then the point is colored blue and similarly for ${\hat{C}}_{B}$ ⁠, the point is colored orange. There is one exception, when Eq. (1) is driven with $(0, 0)$ ⁠, its state remains at the origin and so too does Eq. (10), this is colored gray. The reconstructed attractors, ${\hat{C}}_{A}$ and ${\hat{C}}_{B}$ ⁠, are both plotted in black.

FIG. 9.

Basins of attraction for C ^ A and C ^ B as viewed from P for ρ = 1.25. The respective attractors are plotted in black.

View large Download slide

Basins of attraction for ${\hat{C}}_{A}$ and ${\hat{C}}_{B}$ as viewed from $P$ for $ρ = 1.25$ ⁠. The respective attractors are plotted in black.

Figure 9 shows how ${\hat{C}}_{A}$ and ${\hat{C}}_{B}$ share $P$ and that the number of points which converge to either attractor is not evenly distributed. In particular, most of the points closest to the location of the reconstructed attractors in $P$ tend toward ${\hat{C}}_{A}$ ⁠. Furthermore, there is no single curve which entirely divides the boundary between the basins of ${\hat{C}}_{A}$ and ${\hat{C}}_{B}$ in Fig. 9, at the same time, a more regular or complex pattern may occur in $S$ ⁠. Nonetheless, the method used to produce Fig. 9 provides an indication that the basin of attraction for the attractor $S_{A}$ in $S$ ⁠, which when projected to $P$ describes ${\hat{C}}_{A}$ ⁠, is of a larger volume than the basin for $S_{B}$ in $S$ ⁠. We comment that this result may not be the case across all random realizations of $M$ and $W_{i n}$ ⁠.

3. Basins for ρ = 0.9, 0.8, 0.4, 0.3

In order to gain a broader understanding of the closed-loop RC’s dynamics when it fails to distinguish between $C_{A}$ and $C_{B}$ ⁠, we now study the different basins of attraction which exist before the closed-loop RC exhibits multifunctionality. The same procedure which generated the basins for $ρ = 1.25$ is now used to find the projected basins of attraction in $P$ when setting $ρ = 0.9, 0.8, 0.4$ ⁠, and $0.3$ in the case of $x_{c e n} = 0$ ⁠.

Prior to achieving multifunctionality, Fig. 2 shows that the closed-loop RC is able to reconstruct $C_{B}$ before $C_{A}$ ⁠. For $ρ = 0.9$ ⁠, Fig. 10(a) illustrates that in addition to ${\hat{C}}_{B}$ there are other basins of attraction for a limit cycle labeled as LC $_{1}$ and a quasi-periodic torus labeled as T $_{1}$ ⁠. So while the aim is to reconstruct a coexistence of $C_{A}$ and $C_{B}$ ⁠, there are other attractors that can be approached from an abundance of initial conditions. These untrained attractors, are attractors that exist in $P$ but were not present during the training, these attractors were also found in Flynn et al.¹²

FIG. 10.

Basins of attraction for C ^ B, LC 1, T 1, LC 2, FP 1, FP 2, FP 3, and FP 4 as viewed from P when setting x c e n = 0 and (a) ρ = 0.9, (b) ρ = 0.8, (c) ρ = 0.4, and (d) ρ = 0.3. The respective attractors are plotted in black.

View large Download slide

Basins of attraction for ${\hat{C}}_{B}$ ⁠, ${LC}_{1}$ ⁠, $T_{1}$ ⁠, ${LC}_{2}$ ⁠, ${FP}_{1}$ ⁠, ${FP}_{2}$ ⁠, ${FP}_{3}$ ⁠, and ${FP}_{4}$ as viewed from $P$ when setting $x_{c e n} = 0$ and (a) $ρ = 0.9$ ⁠, (b) $ρ = 0.8$ ⁠, (c) $ρ = 0.4$ ⁠, and (d) $ρ = 0.3$ ⁠. The respective attractors are plotted in black.

In the case where neither $C_{A}$ or $C_{B}$ can be reconstructed, for example, when $ρ = 0.8$ as seen in Figs. 2 and 6, Fig. 10(b) shows that there are basins of attraction for two different limit cycles, LC $_{1}$ from Fig. 10(a) and another limit cycle denoted as LC $_{2}$ ⁠. So at one point before achieving multifunctionality, there is a coexistence of two different limit cycles in $P$ ⁠.

When setting $ρ = 0.4$ ⁠, Fig. 10(c) reveals that four fixed points ${FP}_{1}$ ⁠, ${FP}_{2}$ ⁠, ${FP}_{3}$ ⁠, and ${FP}_{4}$ are the only stable attractors which populate $P$ ⁠. Far away from the origin it can be seen that $P$ is split into nearly equal quadrants. However, closer to the origin there is a greater interaction between the basins resulting in the emergence of a fractal pattern. A similar result is found when setting $ρ = 0.3$ in Fig. 10(d) with fractal basins boundaries appearing closer to the origin. For further reading on fractal basin boundaries and the difficulties these cause in predicting the long-term behavior of a dynamical system see: Grebogi et al.^39,40 and McDonald et al.⁴¹

The symmetrical properties of the basins discussed in the present section are analyzed in Appendix B.

V. SEEING DOUBLE DYNAMICS

Section IV outlines where multifunctionality is achieved in the $(x_{c e n}, ρ)$ -plane by assessing the long-term dynamics of the closed-loop RC given by Eq. (10). However, from Figs. 2(a), 2(b), and 10, it appears that the closed-loop RC must overcome the influence of several untrained attractors before either $C_{A}$ or $C_{B}$ can be reconstructed. From our previous work in Flynn et al.,¹² it is clear that there are many advantages to treating Eq. (10) like any other dynamical system as by exploring the closed-loop RC’s dynamics through a bifurcation analysis further advancements on, for instance, how the RC learns to solve a given problem can be made. In this section, the interplay between the untrained attractors and the reconstructed attractors is investigated in greater detail and the particular bifurcations through which the RC learns how to solve the seeing double problem are identified.

A. Bifurcation analysis for x_cen = 0

In Sec. IV D, multiple fixed points, limit cycles, and a torus were found to coexist in $P$ at different $ρ$ values. However, according to the results shown in Fig. 9, these attractors seemingly vanish from $P$ as multifunctionality is achieved. Moreover, for a small change in $ρ$ there is also a small change in the location of these fixed points and the nature of these limit cycles. This intriguing behavior warrants further investigation in order to establish the suggested link between these untrained attractors and the reconstructed attractors.

Taking the above statements into account, the changes in the dynamics of the fixed points, ${FP}_{1}$ ⁠, ${FP}_{2}$ ⁠, ${FP}_{3}$ ⁠, and ${FP}_{4}$ ⁠, shown in Figs. 10(c)–10(d) are now tracked with respect to changes in $ρ$ ⁠. The same tracking procedure in Flynn et al.¹² is used here which consists of the following: the changes in the dynamics of a given attractor with respect to, for instance, $ρ$ is tracked by repeating the process of incrementally changing $ρ$ ⁠, retraining, and initializing the state of the closed-loop RC with the attractor corresponding to the previous $ρ$ ⁠. If there is a change in the dynamics of the closed-loop RC, for instance, a bifurcation from a fixed point to limit cycle or a period-doubling (P-D) bifurcation, then the dynamics of the resulting stable attractor is tracked. Figure 11 shows the results of this.

FIG. 11.

Stable attractors in P for ρ ∈ [ 0.01 , 1.47 ].

View large Download slide

Stable attractors in $P$ for $ρ \in [0.01, 1.47]$ ⁠.

While the particular bifurcation structure in Fig. 11 may be heavily dependent on the choice of $M$ and $W_{i n}$ ⁠, overall, Fig. 11 provides a clear picture on the challenges the closed-loop RC overcomes in order to achieve multifunctionality. For instance, Fig. 11 shows that for small $ρ$ the closed-loop RC is unable to create any time-varying dynamics and instead produces two pairs of antisymmetric fixed points. As $ρ$ is increased, Fig. 11 establishes that these branches of mirrored fixed points are drawn closer together and two limit cycles are born in quick succession of one another. These are the same limit cycles, ${LC}_{1}$ and ${LC}_{2}$ ⁠, that were found previously in Fig. 10(b).

The role played by these fixed points and limit cycles is akin to the formation of a memory in the brain. As the closed-loop RC learns the correct memory, in this case the correct attractors, along the way its dynamics begins to bear a stronger resemblance to ${\hat{C}}_{A}$ and ${\hat{C}}_{B}$ as the influence of the training data grows until attractor reconstruction is achieved. Once $ρ$ is sufficiently large, Fig. 11 shows that the closed-loop RC is first able to reconstruct $C_{B}$ and subsequently for larger $ρ$ ⁠, the closed-loop RC achieves multifunctionality. A more comprehensive breakdown on each of the major events described above is outlined in the remainder of Sec. V.

1. Bifurcations of FP₁ and FP₂ to LC₁, and of FP₃ and FP₄ to LC₂

Figure 12 provides a more detailed picture on the bifurcations from ${FP}_{1}$ and ${FP}_{2}$ to ${LC}_{2}$ at $ρ = 0.4935$ in the left hand panel and from ${FP}_{3}$ and ${FP}_{4}$ to ${LC}_{1}$ at $ρ = 0.4931$ in the middle panel. What is plotted here are the location of the fixed points and the local maxima and minima of the limit cycles in the $(ρ, x)$ -plane. When tracking the evolution of these limit cycles for decreasing $ρ$ ⁠, it is found that ${LC}_{2}$ can no longer be tracked for $ρ < 0.4905$ and similarly for ${LC}_{1}$ for $ρ < 0.4847$ ⁠.

FIG. 12.

Left and middle panels illustrate the bifurcations from FP 1 and FP 2 to LC 1, and from FP 3 and FP 4 to LC 2, for ρ ∈ [ 0.2 , 0.8 ] in the ( ρ , x )-plane. Right panel shows the location of the fixed points and limit cycles in P just before they can no longer be tracked.

View large Download slide

Left and middle panels illustrate the bifurcations from ${FP}_{1}$ and ${FP}_{2}$ to ${LC}_{1}$ ⁠, and from ${FP}_{3}$ and ${FP}_{4}$ to ${LC}_{2}$ ⁠, for $ρ \in [0.2, 0.8]$ in the $(ρ, x)$ -plane. Right panel shows the location of the fixed points and limit cycles in $P$ just before they can no longer be tracked.

The panel on the right of Fig. 12 displays the location of the fixed points and limit cycles just before and after the respective bifurcations take place on each fixed point and limit cycle. It also appears here that the location of the fixed points just prior to the bifurcation is related to the amplitude of oscillation for both limits cycles after bifurcation.

Figure 13 shows that as $ρ$ decreases and approaches the point at which ${LC}_{1}$ and ${LC}_{2}$ can no longer be tracked, the period, $T (\cdot)$ ⁠, of both limit cycles tends to infinity. $T (\cdot)$ is plotted here in units of $t / 2 π$ to provide a comparison to the period of $C_{A}$ and $C_{B}$ which correspond to $T (\cdot) = 1$ ⁠.

FIG. 13.

Period, T ( ⋅ ), of LC 1 and LC 2 in units of t / 2 π vs ρ.

View large Download slide

Period, T $(\cdot)$ ⁠, of ${LC}_{1}$ and ${LC}_{2}$ in units of $t / 2 π$ vs $ρ$ ⁠.

Figure 13 also shows that as $ρ$ increases, T(⁠ ${LC}_{1}$ ⁠) and T(⁠ ${LC}_{2}$ ⁠) converge to $1$ ⁠. This provides further indication that while there does not appear to be a smooth transition from either ${LC}_{1}$ or ${LC}_{2}$ to ${\hat{C}}_{A}$ or ${\hat{C}}_{B}$ ⁠, these untrained attractors bear a stronger resemblance to the desired dynamics as $ρ$ increases.

Taking all of the above information on ${LC}_{1}$ and ${LC}_{2}$ into account, there are two potential bifurcations which describe how these transitions from two fixed points to one limit cycle takes place. On one hand, since there appears to be a small range of $ρ$ values where ${FP}_{1}$ ⁠, ${FP}_{2}$ ⁠, and ${LC}_{2}$ coexist for $ρ \in (0.4905, 0.4935)$ and similarly where ${FP}_{3}$ ⁠, ${FP}_{4}$ ⁠, and ${LC}_{1}$ coexist for $ρ \in (0.4847, 0.4931)$ ⁠, then a homoclinic bifurcation may be responsible for the birth of these stable limit cycles. On the other hand, since these windows of multistability appear for only a small range of $ρ$ values, their existence may be an artefact of the tracking procedure used to compute the evolution of these fixed points and limit cycles with regard to changes in $ρ$ ⁠. If there is no window of multistability, then a saddle-node infinite period (SNIPER) bifurcation may instead be responsible for the transition from antisymmetric pairs of fixed points to the corresponding limit cycle.

The above quandary identifies a pitfall of this tracking procedure as it is unable to determine when, for instance, an attractor becomes unstable. Instead, this approach relies on a suitably chosen step size of the bifurcation parameter, in this case $ρ$ ⁠, to account for when a bifurcation takes place. If the step size is sufficiently small, then this method can, for example, continue to track unstable solutions of the given dynamical system for a certain range of the bifurcation parameter values. Despite these issues, the results presented in this section provide significant insight toward the role of untrained attractors, how this closed-loop RC solves the seeing double problem and achieves multifunctionality.

2. Bifurcations of LC₁ and LC₂ and the birth of ${\hat{C}}_{B}$

Figure 11 also illustrates how the limit cycles, ${LC}_{1}$ and ${LC}_{2}$ ⁠, evolve in $P$ with regard to changes in $ρ$ ⁠. An extensive account of this behavior and the events which result in the death of these stable limit cycles and the creation of the first reconstructed cycle, ${\hat{C}}_{B}$ ⁠, is provided in Fig. 14.

FIG. 14.

Top left panel shows how the local x maxima of LC 1, LC 2, and C ^ B change for changes in ρ. Top right panel shows the attractors in P for ρ = 0.9. Second, third, and fourth rows show snapshots of C ^ B, LC 2, and LC 1 before they can no longer be tracked.

View large Download slide

Top left panel shows how the local $x$ maxima of ${LC}_{1}$ ⁠, ${LC}_{2}$ ⁠, and ${\hat{C}}_{B}$ change for changes in $ρ$ ⁠. Top right panel shows the attractors in $P$ for $ρ = 0.9$ ⁠. Second, third, and fourth rows show snapshots of ${\hat{C}}_{B}$ ⁠, ${LC}_{2}$ ⁠, and ${LC}_{1}$ before they can no longer be tracked.

The evolution of the local maxima of the closed-loop RC’s stable solutions in the projected $x$ coordinate (denoted as $x_{m}$ ⁠) for $ρ \in (0.848, 0.948)$ is plotted in the top left panel of Fig. 14. Here, it can be seen that the quasi-periodic torus, $T_{1}$ ⁠, shown earlier in Fig. 10(a) comes into existence due to ${LC}_{2}$ undergoing a torus bifurcation at $ρ = 0.8746$ ⁠. In light of this, the top right panel provides an updated picture on the stable attractors which, in Fig. 10(a), were shown to coexist in $P$ for $ρ = 0.9$ ⁠. The color scheme in the top right panel is consistent with the bifurcation diagram shown in the top left panel instead of assigning the light gray color to the torus as in Fig. 10(a).

The second row of images in Fig. 14 provides additional insight toward the sequence of events which result in the creation/destruction of ${\hat{C}}_{B}$ as $ρ$ increases/decreases. Moving from left to right, these images show the result of tracking the evolution of ${\hat{C}}_{B}$ as $ρ$ decreases. The dashed trajectory indicates that ${\hat{C}}_{B}$ is found to either become unstable or no longer exist between $ρ = 0.8937$ and $0.8936$ and the state of the closed-loop RC is on a transient to ${LC}_{1}$ ⁠. The closed-loop RC then remains on ${LC}_{1}$ when $ρ$ is decreased further as illustrated by the image on the right for $ρ = 0.8935$ ⁠.

The third row of images in Fig. 14 depicts the previously mentioned torus bifurcation and the subsequent death of this torus as $ρ$ increases. By increasing $ρ$ from $0.9073$ to $0.9074$ ⁠, the torus can no longer be tracked. Following this bifurcation, the state of the closed-loop RC departs from the unstable torus (as illustrated by the dashed line) and approaches ${\hat{C}}_{B}$ and, according to the tracking algorithm, remains on ${\hat{C}}_{B}$ when $ρ$ is increased as shown in the right panel for $ρ = 0.9075$ ⁠.

The bottom row of images in Fig. 14 shows that by increasing $ρ$ from $0.9419$ to $0.9420$ ${LC}_{1}$ can no longer be tracked. As ${\hat{C}}_{B}$ is the only remaining stable solution of the closed-loop RC for $ρ = 0.9420$ ⁠, the state of the closed-loop RC has no choice but to embark on a transient from the previously stable LC $_{1}$ to ${\hat{C}}_{B}$ (plotted as the dashed line) and remains on ${\hat{C}}_{B}$ as $ρ$ is increased to $0.9421$ ⁠.

FIG. 15.

Top panel shows how the local x maxima of C ^ B and C ^ A change for changes in ρ. Second and third rows show snapshots of C ^ B and C ^ A before they can no longer be tracked.

View large Download slide

Top panel shows how the local $x$ maxima of ${\hat{C}}_{B}$ and ${\hat{C}}_{A}$ change for changes in $ρ$ ⁠. Second and third rows show snapshots of ${\hat{C}}_{B}$ and ${\hat{C}}_{A}$ before they can no longer be tracked.

3. The rise and fall of multifunctionality

While Fig. 11 highlights the range of $ρ$ values where multifunctionality is achieved, Fig. 15 expands on this and shows how the rise and fall of multifunctionality begins with the arrival of ${\hat{C}}_{A}$ at $ρ = 1.0831$ and death of ${\hat{C}}_{B}$ at $ρ = 1.4297$ ⁠.

The top panel of Fig. 15 picks up where Fig. 14 leaves off by continuing to track the changes in the dynamics of ${\hat{C}}_{B}$ with regard to changes in $ρ$ ⁠. The second row of images in Fig. 15 highlights the circumstances in which ${\hat{C}}_{B}$ can no longer be tracked. By increasing $ρ$ from $1.4296$ to $1.4297$ ⁠, it is shown here that ${\hat{C}}_{B}$ no longer has a fixed amplitude of oscillation. Following this, by increasing $ρ$ to $1.4298$ ⁠, ${\hat{C}}_{B}$ is no longer a stable limit cycle and the state of the closed-loop RC is instead on a transient to ${\hat{C}}_{A}$ and remains on ${\hat{C}}_{A}$ for increasing $ρ$ ⁠.

The third row of images in Fig. 15 depicts how ${\hat{C}}_{A}$ is born. Shown from right to left is the result of tracking the evolution of ${\hat{C}}_{A}$ for decreasing $ρ$ ⁠. Here it can be seen that between $ρ = 1.0831$ and $1.0830$ ⁠, ${\hat{C}}_{A}$ can no longer be tracked and the state of the closed-loop RC subsequently tends toward ${\hat{C}}_{B}$ and remains on ${\hat{C}}_{B}$ for smaller $ρ$ values.

While Fig. 6 shows an instance in which the closed-loop RC fails to achieve multifunctionality for $ρ = 1.7$ ⁠, the top panel of Fig. 15 reveals the means in which this comes about. Following the death of ${\hat{C}}_{B}$ at $ρ = 1.4298$ ⁠, ${\hat{C}}_{A}$ is the only stable solution present in $P$ ⁠. However, ${\hat{C}}_{A}$ is found to undergo a torus bifurcation at $ρ = 1.4711$ ⁠. Following this, several windows of periodic and chaotic behavior appear and the closed-loop RC remains chaotic for $ρ > 1.6$ ⁠. The first instance of chaos for increasing $ρ$ is in agreement with the corresponding Lyapunov analysis in Fig. 3. As $ρ$ is increased beyond this point, the range of local maxima in the $x$ variable, as plotted here, continues to grow. As a result, the dynamics of the closed-loop RC loses its resemblance to the original driving input.

Furthermore, the bifurcation diagram presented in Fig. 15 provides additional information on the results shown previously in Fig. 6. Figure 15 shows that the trajectories taken by the closed-loop RC when initialized with either $r_{(C_{A})} (t_{listen})$ or $r_{(C_{B})} (t_{listen})$ for $ρ = 1.7$ are trajectories on the same chaotic attractor. However, we do not claim that the mechanisms responsible for the loss of multifunctionality are universal across all possible $M$ and $W_{i n}$ as we find, from further analysis not presented in the present paper, examples where ${\hat{C}}_{A}$ is first to disappear as $ρ$ increases.

4. Floquet analysis

The results presented thus far in Sec. V illustrate some of the hurdles that the closed-loop RC overcomes in order to achieve multifunctionality. However, given the nature of the seeing double problem, it is also possible to provide some quantitative insight into how multifunctionality is achieved by conducting a Floquet analysis. This involves computing the eigenvalues of the “monodromy matrix,”

Q

⁠, to classify the stability of a particular attractor using a trajectory of one period on the attractor.

Q

is the solution of

\dot{Q} (t) = J (t) Q (t), Q (0) = I_{N},

(16)

after one period,

T

⁠.

J (t)

is the Jacobian matrix, which in this case, is the Jacobian of the closed-loop RC in Eq. and

J (t)

is computed at each point along a trajectory of one period on the attractor. The eigenvalues of

Q

⁠, denoted by

μ_{i}

⁠, are called the Floquet multipliers. In the case of the problem at hand, by solving Eq. with input of one period from, for instance,

{\hat{C}}_{A}

then the largest multiplier,

μ_{1}

⁠, will have magnitude

= 1

and all other

| μ_{i} | < 1

⁠. While this provides information on when

{\hat{C}}_{A}

and

{\hat{C}}_{B}

are stable at certain

ρ

values, the particular bifurcation through which the attractors become stable/unstable can also be determined. Furthermore, dynamical features which the tracking algorithm (used to compute the bifurcation diagrams shown previously) fails to notice may be unearthed by conducting this Floquet analysis. For more information on Floquet theory, see Nayfeh and Balachandran.⁴²

Taking this into account, Fig. 16 shows how the five largest Floquet multipliers evolve for changes in $ρ$ when solving Eq. (16) with drive from ${\hat{C}}_{A}$ and ${\hat{C}}_{B}$ in the top and bottom panels, respectively. In each case, $| μ_{1} | = 1$ at all $ρ$ values where ${\hat{C}}_{A}$ and ${\hat{C}}_{B}$ were found to exist as stable attractors.

FIG. 16.

Magnitude, real part, and imaginary part of the first i = 5 largest Floquet multipliers, μ i, for changes in ρ for C ^ A (top row) and C ^ B (bottom row). Dashed vertical lines highlight the ρ values where C ^ A and C ^ B were found to exist between.

View large Download slide

Magnitude, real part, and imaginary part of the first $i = 5$ largest Floquet multipliers, $μ_{i}$ ⁠, for changes in $ρ$ for ${\hat{C}}_{A}$ (top row) and ${\hat{C}}_{B}$ (bottom row). Dashed vertical lines highlight the $ρ$ values where ${\hat{C}}_{A}$ and ${\hat{C}}_{B}$ were found to exist between.

As $ρ$ decreases, Fig. 16 shows that Re $(μ_{2}) \to 1$ and since all Im $(μ_{i}) = 0$ it can be deduced that ${\hat{C}}_{A}$ and ${\hat{C}}_{B}$ are born through a saddle-node bifurcation. For ${\hat{C}}_{B}$ ⁠, Fig. 16 demonstrates that Re $(μ_{2})$ and Re $(μ_{3})$ become equal at $ρ = 1.3690$ and beyond this point the corresponding imaginary components of $μ_{2}$ and $μ_{3}$ grow according to Im $(μ_{2}) = - Im (μ_{3})$ until $| μ_{2} |$ and $| μ_{3} |$ tend to $1$ resulting in ${\hat{C}}_{B}$ becoming unstable at $ρ = 1.4927$ ⁠. Based on this information and the results shown in Figs. 15 and 16 indicates that ${\hat{C}}_{B}$ becomes unstable through a subcritical torus bifurcation. The behavior of the Floquet multipliers computed for ${\hat{C}}_{A}$ in Fig. 16 establishes that a supercritical torus bifurcation is responsible for the ${\hat{C}}_{A}$ ’s change in dynamics from a stable limit cycle to stable torus at $ρ = 1.4711$ as illustrated in Fig. 15. Figure 16 also reveals further information which Fig. 15 fails to capture. For instance, a symmetry breaking bifurcation on ${\hat{C}}_{A}$ is found to occur at $ρ = 1.4450$ as Re $(μ_{2})$ sharply increases to $1$ and quickly decreases from $1$ ⁠.

FIG. 17.

Result of tracking the evolution of stable attractors in P for ρ ∈ [ 0.48 , 0.52 ] and x c e n ∈ [ 0.575 , − 0.575 ]. Curves represent the end point of tracking stable solutions. Filled regions describe the nature of the closed-loop RC’s stable solutions.

View large Download slide

Result of tracking the evolution of stable attractors in $P$ for $ρ \in [0.48, 0.52]$ and $x_{c e n} \in [0.575, - 0.575]$ ⁠. Curves represent the end point of tracking stable solutions. Filled regions describe the nature of the closed-loop RC’s stable solutions.

B. Bifurcation analysis in the (x_cen, ρ)-plane

Section V A provides an extensive account of the particular bifurcations the closed-loop RC goes through in order to achieve multifunctionality when $x_{c e n} = 0$ ⁠. These results suggest that for increasing $ρ$ ⁠, on one hand, the dynamics of ${FP}_{1}$ ⁠, ${FP}_{2}$ ⁠, ${FP}_{3}$ ⁠, and ${FP}_{4}$ which lead to the creation of ${LC}_{1}$ and ${LC}_{2}$ ⁠, symbolizes a set of challenges for the closed-loop RC to overcome in order to achieve multifunctionality. On the other hand, as $ρ$ increases it can also be argued that the role these fixed points and limit cycles play is to sculpt $P$ for the arrival of ${\hat{C}}_{A}$ and ${\hat{C}}_{B}$ ⁠. Furthermore, given the ubiquity of the Goldilocks effect across multiple random realizations of $M$ and $W_{i n}$ as mentioned in Sec. IV A, it is worthwhile to study the behavior of these untrained attractors at different $x_{c e n}$ values.

In this section, many of the tools applied in Sec. V A are used to explore the dynamics of the closed-loop RC in the $(x_{c e n}, ρ)$ -plane. The results presented in this section reveal how the relationships between the fixed points and limit cycles mentioned above evolve with respect to changes in $x_{c e n}$ nearby where these six attractors coexist for a small range of $ρ$ values when $x_{c e n} = 0$ in Figs. 11 and 12. These results also provide insight on the different roads to $C_{A}$ and $C_{B}$ ’s reconstruction as conceivably carved out by these fixed points and limit cycles.

1. Tracking stable solutions in the (x_cen, ρ)-plane

The result of tracking the changes in the dynamics of the closed-loop RC’s stable solutions for $ρ \in [0.48, 0.52]$ and $x_{c e n} \in [0.575, - 0.575]$ is presented in Fig. 17. These results are generated by tracking the evolution of ${LC}_{1}$ and ${LC}_{2}$ for changes in $x_{c e n}$ beginning nearby the bifurcations at $ρ = 0.4847$ and $0.4905$ ⁠, which spawn ${LC}_{1}$ and ${LC}_{2}$ ⁠, respectively when $x_{c e n} = 0$ ⁠, as illustrated in Figs. 11 and 12. Similarly, the evolution of ${FP}_{1}$ ⁠, ${FP}_{2}$ ⁠, ${FP}_{3}$ ⁠, and ${FP}_{4}$ is also tracked with respect to changes in $x_{c e n}$ nearby the bifurcations at $ρ = 0.4935$ and $0.4931$ that result in the death of these fixed points. Once these attractors can no longer be tracked for changes in $x_{c e n}$ ⁠, this location in the $(x_{c e n}, ρ)$ -plane is recorded and added to the corresponding curve in Fig. 17. The same process is repeated for incremental changes in $ρ$ for $ρ \in [0.48, 0.52]$ ⁠. We did not find any additional attractors in the portion of the $(x_{c e n}, ρ)$ -plane shown in Fig. 17 that were unrelated to those studied previously in Sec. V A.

Each curve and filled region in Fig. 17 has a specific meaning which we outline below.

The labels and color scheme of each curve in Fig. 17 are consistent with illustrations of the attractors studied in the previous section. While, for example, the stable period-1 limit cycle, ${LC}_{1}$ ⁠, is found to undergo P-D bifurcations at several locations in the $(x_{c e n}, ρ)$ -plane, the resulting stable period-2 limit cycle is also referred to as ${LC}_{1}$ ⁠. This convention, which was also taken in relation to ${LC}_{2}$ and $T_{1}$ in Fig. 14, is applied even when, for instance, an infinite sequence of P-D bifurcations have taken place when tracking the dynamics of the original period-1 limit cycle. What is not shown nor highlighted along these curves are the continuation of the corresponding codim-1 bifurcation curves and location of codim-2 points associated with the dynamics of certain attractors in the $(x_{c e n}, ρ)$ -plane. This decision was taken in order to avoid overcrowding Fig. 17 with information that may distract from its main message of tracking how the attractors mentioned above behave in the $(x_{c e n}, ρ)$ -plane. Instead, a more detailed breakdown on some of the bifurcations shown in Fig. 17 is provided later in the present section.

The color scheme of the filled regions in the portion of the $(x_{c e n}, ρ)$ -plane shown in Fig. 17, as defined by the area the above curves enclose, describes the nature of the closed-loop RC’s stable solutions in $P$ for a given $x_{c e n}$ and $ρ$ ⁠.

2. General comments on RC’s dynamics in the (x_cen, ρ)-plane

The primary message of Fig. 17 continues the central narrative of this paper, the greater the overlap between $C_{A}$ and $C_{B}$ ⁠, the closed-loop RC requires larger $ρ$ values to recall more information from the past in order to achieve multifunctionality and the greater the influence of these fixed points and limit cycles until reconstruction occurs. Figure 17 further substantiates this statement through examining the behavior of the stable solutions in $P$ for the selected $x_{c e n}$ and $ρ$ values. What this reveals is the wider role these attractors play in terms of how the closed-loop RC solves the seeing double problem.

One of the main features of Fig. 17 is the symmetry in the closed-loop RC’s dynamics about $x_{c e n} = 0$ and this is studied in much greater detail in Appendix B.

3. Behavior of FP₃ and FP₄

Concentrating first on ${FP}_{3}$ and ${FP}_{4}$ ⁠, Fig. 17 reveals that as the magnitude of $x_{c e n}$ increases these fixed points play a diminishing role in $P$ ⁠. At $ρ = 0.48$ ⁠, these fixed points cease to coexist for $| x_{c e n} | > 0.055$ and both no longer exist for $| x_{c e n} | > 0.13$ ⁠. The inset plot on the right hand side of Fig. 17 highlights that for increasing $ρ$ along the ${FP}_{3}$ and ${FP}_{4}$ curves from $x_{c e n} = 0.13$ and $- 0.13$ ⁠, respectively, the ${FP}_{3}$ and ${FP}_{4}$ curves become point-like at $x_{c e n} \approx - 0.035$ and $0.035$ ⁠, respectively, and can no longer be tracked for $ρ > 0.495$ ⁠.

4. Behavior of FP₁ and FP₂

Figure 17 reveals that as the magnitude of $x_{c e n}$ increases, ${FP}_{1}$ and ${FP}_{2}$ begin to play a more prominent role in the closed-loop RC’s dynamics in comparison to ${FP}_{3}$ and ${FP}_{4}$ for the particular $ρ$ values explored here in the $(x_{c e n}, ρ)$ -plane.

While Fig. 17 shows that the coexistence of ${FP}_{1}$ and ${FP}_{2}$ comes to an end at $ρ = 0.48$ for $| x_{c e n} | > 0.07$ ⁠, these fixed points exist far beyond the limits of the $x_{c e n}$ axis displayed here. These fixed points occupy the largest combined area in the portion of the $(x_{c e n}, ρ)$ -plane shown here. For increasing $ρ$ along the curves which enclose the regions of the $(x_{c e n}, ρ)$ -plane where ${FP}_{1}$ and ${FP}_{2}$ are both stable, the same fate as ${FP}_{3}$ and ${FP}_{4}$ awaits ${FP}_{1}$ and ${FP}_{2}$ as these curves become point-like at $x_{c e n} \approx 0.34$ and $- 0.34$ ⁠, respectively.

As a further comment, the inset plot on the right hand side of Fig. 17 provides additional insight to the bifurcation diagrams in Figs. 11 and 12 and draws attention to how the coexistence between ${FP}_{1}$ and ${FP}_{2}$ ⁠, and ${FP}_{3}$ and ${FP}_{4}$ at $x_{c e n} = 0$ unfolds with respect to changes in $x_{c e n}$ ⁠.

5. Behavior of LC₁

Focusing next on ${LC}_{1}$ ⁠, Fig. 17 illustrates the limited role this limit cycle plays in the closed-loop RC’s dynamics as the magnitude of $x_{c e n}$ increases from $0$ ⁠. For decreasing $ρ$ ⁠, the range of $x_{c e n}$ values where ${LC}_{1}$ is stable also decreases. By increasing $| x_{c e n} |$ from $0$ for $ρ = 0.48$ ⁠, a region in the $(x_{c e n}, ρ)$ -plane where all four fixed points coexist, LC $_{1}$ is first encountered as a stable solution at $x_{c e n} \approx \pm 0.05$ ⁠. ${LC}_{1}$ briefly coexists with these four fixed points until ${FP}_{3}$ and ${FP}_{4}$ can no longer be tracked as previously mentioned at $x_{c e n} = \mp 0.055$ and likewise for ${FP}_{1}$ and ${FP}_{2}$ at $x_{c e n} = \pm 0.07$ ⁠. This leaves ${LC}_{1}$ in coexistence with ${FP}_{2}$ and ${FP}_{3}$ for $x_{c e n} \in [0.071, 0.097]$ and in coexistence with ${FP}_{1}$ and ${FP}_{4}$ for $x_{c e n} \in [- 0.071, - 0.097]$ ⁠. ${LC}_{1}$ can no longer be tracked for $| x_{c e n} | > 0.097$ at $ρ = 0.48$ ⁠.

As a further comment, as $ρ$ decreases, the brown curves on either side of $x_{c e n} = 0$ ⁠, which outline where ${LC}_{1}$ is stable, appear to be converging to a point along the ${FP}_{3}$ and ${FP}_{4}$ curves. This suggests that a codim-2 bifurcation takes place in the $(x_{c e n}, ρ)$ -plane outside of the regions explored in Fig. 17.

Prior to no longer being able to track ${LC}_{1}$ for increasing the magnitude of $x_{c e n}$ ⁠, Fig. 17 also illustrates that ${LC}_{1}$ begins to exhibit chaos along the dotted brown curve $\forall ρ$ values shown here. The extent of this chaos in the $(x_{c e n}, ρ)$ -plane is described by the area enclosed by this dotted curve and the ${LC}_{1}$ curve, and is further emphasized by the hatched pattern. It is shown here that, on one hand, these windows of chaos widen as $ρ$ increases, on the other hand, the inset plot on the left hand side of Fig. 17 provides a more detailed picture on how, this small region of chaos narrows as $ρ$ decreases.

Figure 18 sheds some light on ${LC}_{1}$ ’s change from periodic to chaotic dynamics. For $ρ = 0.515$ ⁠, it is shown here that as $x_{c e n}$ increases, ${LC}_{1}$ undergoes a series of P-D bifurcations. An infinite sequence of these bifurcations occurs between $x_{c e n} = 0.1301$ and $0.1302$ which gives rise to the chaotic dynamics shown in Fig. 18 for $x_{c e n} = 0.1303$ ⁠.

FIG. 18.

View large Download slide

Evolution of ${LC}_{1}$ in $P$ for changes in $x_{c e n}$ with $ρ = 0.515$ ⁠.

6. Behavior of LC₂

Figure 17 lays bare the intricate relationship between ${LC}_{2}$ ⁠, $ρ$ ⁠, and $x_{c e n}$ ⁠. Several events occur along the ${LC}_{2}$ curve as the magnitude of $x_{c e n}$ increases. The inset plot on the right hand side highlights the first of these incidents with the loss of its coexistence with the four fixed points and ${LC}_{1}$ ⁠. This continues with larger $| x_{c e n} |$ values and at $ρ = 0.502$ ⁠, ${LC}_{2}$ is in coexistence with only ${FP}_{1}$ for $x_{c e n} < - 0.115$ and FP $_{2}$ for $x_{c e n} > 0.115$ until $C_{B}$ is reconstructed at larger $x_{c e n}$ values.

It is also shown in Fig. 17 that nearby these points in the $(x_{c e n}, ρ)$ -plane where ${LC}_{2}$ loses its coexistence with ${LC}_{1}$ ⁠, the ${LC}_{2}$ curve itself reaches its first of several local maxima. In order to continue tracking ${LC}_{2}$ beyond this point, $ρ$ is required to decrease as $| x_{c e n} |$ increases until the curve reaches its first of several local minima at $(x_{c e n}, ρ) = (\pm 0.16, 0.49)$ where a P-D bifurcation occurs. Subsequent P-D bifurcations are encountered along the ${LC}_{2}$ curve for slightly increasing $| x_{c e n} |$ and $ρ$ ⁠. The hatched region of the $(x, ρ)$ -plane enclosed by the ${LC}_{2}$ curve and the nearby dotted curve characterize the small window of chaotic solutions produced by an infinite sequence of these P-D bifurcations.

Figure 19(a) provides an example of the closed-loop RC’s dynamics nearby one of these P-D cascades for $ρ = 0.501$ ⁠. It is shown here that when tracking the evolution of the original stable period-1 limit cycle, successive P-D bifurcations are encountered for increasing $x_{c e n}$ ⁠. Figure 19(a) shows examples of the corresponding period-1 dynamics at $x_{c e n} = 0.1500$ ⁠, period-2 dynamics at $x_{c e n} = 0.1660$ ⁠, and period-4 dynamics at $x_{c e n} = 0.1665$ ⁠. In this scenario the period-8 limit cycle, whose dynamics are shown here for $x_{c e n} = 0.1667$ ⁠, undergoes an infinite sequence of P-D bifurcations for a small increase in $x_{c e n}$ ⁠. This produces the chaotic dynamics shown here for $x_{c e n} = 0.1669$ ⁠. This region of chaos in Fig. 17 ends once the ${LC}_{2}$ curve and dotted curve meet at $(x_{c e n}, ρ) = (0.178, 0.507)$ ⁠. For increasing $| x_{c e n} |$ beyond this point, it is found that ${LC}_{2}$ is a period-2 limit cycle before it can no longer be tracked. However, subsequent P-D bifurcations occur along the ${LC}_{2}$ curve for increasing $| x_{c e n} |$ ⁠.

FIG. 19.

View large Download slide

Evolution of ${LC}_{2}$ in $P$ for changes in $x_{c e n}$ with $ρ = 0.501$ in (a) and $0.511$ in (b).

Following an infinite sequence of these P-D bifurcations, what Fig. 17 also shows is that along the LC $_{2}$ curve, chaos reappears at $(x_{c e n}, ρ) = (\pm 0.35, 0.508)$ ⁠. The dotted magenta curve which stems from this point distinguishes between when LC $_{2}$ is periodic or chaotic as further emphasized by the hatched pattern. An example of this particular route to chaos is provided in Fig. 19(b). Furthermore, although it is not shown here, these P-D bifurcation do not only occur along the LC $_{2}$ curve, there is a locus which connects each given set of P-D bifurcation points in the $(x_{c e n}, ρ)$ -plane in a similar way to the dotted magenta curve mentioned above.

What Fig. 19(b) depicts is the result of tracking the evolution of ${LC}_{2}$ ’s dynamics for changes in $x_{c e n}$ when $ρ = 0.511$ ⁠. The remnants of the previous small window of chaos, which occurs after the first local minima along the ${LC}_{2}$ curve, are visible here with a period-2 limit cycle found for $x_{c e n} = 0.2$ ⁠. The major difference between this period-2 limit cycle and the example illustrated in Fig. 19(a) for $ρ = 0.501$ is the prominence of the sharp turn around the max $y$ value on the period-2 limit cycle in Fig. 19(b). This remains a distinctive feature of ${LC}_{2}$ when tracking its dynamics as $x_{c e n}$ is further increased and persists after each of the subsequent P-D bifurcations that result in chaos. It is also worth noting that this kink along ${LC}_{2}$ ’s dynamics occurs nearby the location of ${FP}_{2}$ and ${FP}_{1}$ for the corresponding positive and negative $x_{c e n}$ values, respectively.

${LC}_{2}$ remains chaotic as $| x_{c e n} |$ increases from $0.35$ and, depending on the sign of $x_{c e n}$ ⁠, continues its coexistence with either ${FP}_{1}$ or ${FP}_{2}$ for the range of $x_{c e n}$ values shown in Fig. 17. Much like ${LC}_{1}$ with ${FP}_{3}$ and ${FP}_{4}$ ⁠, Fig. 17 also suggests that the ${LC}_{1}$ curve may connect with the ${FP}_{1}$ and ${FP}_{2}$ curves beyond the range of $x_{c e n}$ values shown here.

As $ρ$ increases along this dotted curve of chaos, which begins at $(x_{c e n}, ρ) = (\pm 0.35, 0.508)$ ⁠, there is a point where these sequences of P-D bifurcations are no longer directly responsible for the onset of chaotic dynamics. For instance, Fig. 20(a) shows that when tracking the evolution of ${LC}_{2}$ for $ρ = 0.515$ ⁠, by increasing $x_{c e n}$ from $0.3623$ to $0.3624$ ⁠, there is a sudden change in dynamics from a period-2 limit cycle to a chaotic attractor. The same is shown Fig. 20(b) for $ρ = 0.517$ and $x_{c e n}$ increases from $0.3558$ to $0.3559$ ⁠. The only difference between these examples of ${LC}_{2}$ ’s dynamics prior to the dawning of chaos is the lack of a period-4 limit cycle for $ρ = 0.517$ ⁠.

FIG. 20.

View large Download slide

Evolution of ${LC}_{2}$ in $P$ for changes in $x_{c e n}$ with $ρ = 0.515$ in (a) and $ρ = 0.517$ in (b).

7. Behavior of ${\hat{C}}_{B}$

Figure 17 also illustrates regions in the $(x_{c e n}, ρ)$ -plane where ${\hat{C}}_{B}$ comes into existence. For much of its presence in Fig. 17, ${\hat{C}}_{B}$ is in coexistence with either, LC $_{2}$ ⁠, or one of ${FP}_{2}$ or ${FP}_{1}$ depending on the sign of $x_{c e n}$ ⁠, or ${LC}_{2}$ and ${FP}_{2}$ or ${FP}_{1}$ ⁠.

At the limits of the $x_{c e n}$ -axis in Fig. 17 (when $| x_{c e n} | = 0.575$ ⁠), the nature of ${\hat{C}}_{B}$ differs greatly depending on $ρ$ ⁠. In particular, the dotted orange curve describes the boundary between when ${\hat{C}}_{B}$ exhibits periodic or chaotic behavior, and it is within the hatched region that ${\hat{C}}_{B}$ is chaotic. For decreasing $| x_{c e n} |$ ⁠, these two curves approach one another and the window of chaos becomes arbitrarily small at $| x_{c e n} | = 0.45$ ⁠. By slightly decreasing $| x_{c e n} |$ from this point in the $(x_{c e n}, ρ)$ -plane, a sudden increase in $ρ$ is required in order to reconstruct ${\hat{C}}_{B}$ ⁠.

Figure 21 illustrates the changes in ${\hat{C}}_{B}$ ’s dynamics when tracking its evolution for decreasing $x_{c e n}$ when $ρ = 0.504$ ⁠. Here it can be seen that ${\hat{C}}_{B}$ is a period-1 limit cycle for $x_{c e n} = 0.7000$ and is also not perfectly circular. ${\hat{C}}_{B}$ undergoes several P-D bifurcations with examples of period-2, period-4, and period-8 dynamics shown here for $x_{c e n} = 0.6400, 0.5200,$ and $0.5205$ ⁠. An infinite sequence of P-D bifurcations occurs as $x_{c e n}$ is decreased further resulting in the chaotic dynamics like those shown here for $x_{c e n} = 0.4965$ ⁠.

FIG. 21.

Evolution of C ^ B in P for changes in x c e n with ρ = 0.504.

View large Download slide

Evolution of ${\hat{C}}_{B}$ in $P$ for changes in $x_{c e n}$ with $ρ = 0.504$ ⁠.

Above all, what Figs. 21 and 17 show is that while ${\hat{C}}_{B}$ is stable its dynamics do not always resemble $C_{B}$ even to the point when ${\hat{C}}_{B}$ is not in any way circular, has a much larger period, or even exhibits chaos. At the same time, it is also worth noting that after a number of P-D bifurcations, for instance the period-4 limit cycle shown in Fig. 21, portions of the trajectory on ${\hat{C}}_{B}$ begin to resemble a more circular orbit with similar characteristics to $C_{B}$ ⁠. This suggests that these P-D bifurcations may play an important role in how the closed-loop RC improves its reconstruction of $C_{B}$ ⁠.

VI. CONCLUSION

The results presented here show that for the closed-loop RC described in Eq. (10) to achieve multifunctionality there is a crucial dependence on $ρ$ ⁠, the spectral radius of the RC’s internal connections. The intricacies of this relationship are revealed in Sec. IV A when the cycles $C_{A}$ and $C_{B}$ are moved closer together. As $C_{A}$ and $C_{B}$ begin to overlap, Fig. 3 shows that the closed-loop RC requires a greater amount of memory in order to successfully distinguish between the orbits. A “Goldilocks” effect is also found whereby if $ρ$ is too small or too large then multifunctionality will not occur.

From our results, we have obtained a greater understanding of multifunctionality in a RC. In the context of the formalism developed in Sec. II C, we have shown that while, for instance, $A_{1} \cap A_{2} \neq \emptyset$ (there can be an overlap in the training data), for a closed-loop RC to achieve multifunctionality, we require that the open-loop RC responds uniquely to each of the driving inputs $u_{(A_{1})} (t)$ and $u_{(A_{2})} (t)$ so that at no time $t$ during the training $r_{(A_{1})} (t) = r_{(A_{2})} (t)$ even if $u_{(A_{1})} (t) = u_{(A_{2})} (t)$ ⁠.

While Figs. 3 and 4 show that the closed-loop RC achieves its best performance just prior to exhibiting chaotic dynamics, Fig. 5 provides some evidence that beyond this transition to chaos the closed-loop RC enters into the dynamical regime of chaotic itinerancy. From a different point of view, multifunctionality can be used as a precursor to generate a specifically designed chaotically itinerant orbit between several attractor ruins thereby providing a suitable means to realize Tsuda’s thoughts on the role of “chaotic itinerancy” in cognitive neurodynamics in Tsuda.⁴³ Tsuda puts forward chaotic itinerancy as a potential dynamical mechanism to describe how the brain transitions between performing different tasks or recounting different memories characterized by each attractor ruin. We intend to study this connection between multifunctionality and chaotic itinerancy and how can multifunctionality be exploited as a data-driven modelling tool in future work.

Section IV D focuses on how the closed-loop RC constructs basins of attraction from the perspective of $P$ for ${\hat{C}}_{A}$ and ${\hat{C}}_{B}$ when they are entirely overlapping, i.e., when $x_{c e n} = 0$ ⁠. The nature of $P$ prior to the closed-loop RC achieving multifunctionality is also explored. Figure 10 unveils the coexistence of several fixed points with fractal basin boundaries, limit cycles and tori in $P$ ⁠, these are untrained attractors which were not involved in the training. This discovery provides the basis for conducting the bifurcation analyses carried out in Sec. V in order to determine the role these attractors play in how the closed-loop RC solves the seeing double problem.

This bifurcation analysis is divided into two studies, the results of the first are outlined in Sec. V A which explores how changes in $ρ$ effects the closed-loop RC’s dynamics in the case of $x_{c e n} = 0$ ⁠. The results of the second study are presented in Sec. V B where the closed-loop RC’s dynamics in a portion of the $(x_{c e n}, ρ)$ -plane are assessed. Here the role of the untrained attractors in the learning dynamics of the closed-loop RC becomes apparent. The bifurcations which result in the rise and fall of multifunctionality are identified through the results of the Floquet analysis shown in Fig. 16. The main take home message from both of these bifurcation analyses is the same, the more $C_{A}$ and $C_{B}$ overlap, the more memory the the closed-loop RC requires in order to exhibit multifunctionality and the greater the influence of these untrained attractors until reconstruction occurs.

Even though the bifurcation analysis tools used in this paper have provided significant insight into how the closed-loop RC learns, performs certain tasks, and interacts with untrained attractors, it is likely that only half the story is being told as the unstable solutions are not accounted for. As mentioned in Flynn et al.,¹² these issues can be alleviated by adapting well-known continuation software like AUTO⁴⁴ to the case of reservoir computing.

While the bifurcation analysis in Sec. V is conducted for only one random realization of $M$ and $W_{i n}$ ⁠, it is also possible that multifunctionality may be achieved in various other scenarios for different choices of $M$ and $W_{i n}$ ⁠. For instance, Morra et al.⁴⁵ investigate the dynamics of the same RC design as in Eq. (1), when trained to solve the seeing double problem in the case of $x_{c e n} = 0$ ⁠, using the same technique that produced Fig. 11. It is shown that for a much smaller number of neurons (⁠ $N = 500$ as opposed to $1000$ ⁠) there is a greater interaction and coexistence between the trained attractors and the untrained attractors in comparison to the results presented in Fig. 11. This may indicate that as $N$ is increased, the less the untrained attractors interact and coexist with the trained attractors. We remark that even though these bifurcation diagrams are produced for smaller $N$ and different random realizations of $M$ and $W_{i n}$ ⁠, they share similar characteristics features in terms of there being a coexistence of four fixed points at small $ρ$ ⁠, and that there is a bifurcation between these fixed points and two limit cycles. The results presented in Morra et al.⁴⁵ also show that by constructing the weights of $M$ based on a region of a fruit-fly’s connectome (where $N = 426$ ⁠) and repeating the same experiment, there is a slightly different bifurcation structure. What is interesting here, in comparison to the results presented in this paper, is that for different designs of $M$ and for $N = 426$ as opposed to $N = 1000$ ⁠, multifunctionality is achieved for a much greater range of $ρ$ values.

We comment that bifurcation analysis on ANNs similar to the closed-loop RC studied throughout this paper have so far mainly concentrated low dimensional ANNs such as the work of Doya,⁴⁶ Beer,^47,48 and Ceni et al.⁴⁹ where the results from a bifurcation analysis of a RC provides the insight needed to solve the “flip-flop problem.” By conducting a bifurcation analysis of low-dimensional RCs and studying how the dynamics of the RC changes as the dimension increases may provide further insight into the behavior of higher dimensional RCs and establish how much multifunctionality a given RC can have.

We do not claim that the training method outlined in Sec. II C is the optimal method for reconstructing a coexistence of any number of attractors. It is clear from the results presented in this paper that even for the closed-loop RC to reconstruct a coexistence of two relatively simple dynamical objects, $C_{A}$ and $C_{B}$ ⁠, it is the relationship between them that is the more crucial aspect. By making further progress on establishing other limiting factors which impact the reconstruction capacity of a given closed-loop RC, it may be possible to extend the current record of the seven coexisting reconstructed chaotic attractors presented in Herteux.¹⁵

In summary, this paper provides significant insight toward how a closed-loop RC achieves multifunctionality in a paradigmatic setting with a particular focus on exploring the dynamics of the closed-loop RC as it solves the seeing double problem. However, there are still many unanswered questions concerning the reconstruction capacity of a given closed-loop RC and the limiting factors to the interesting dynamics a multifunctional RC can produce. While this paper improves the language surrounding the basic mechanisms needed for multifunctionality to occur from previous studies, in future work we aim to extend the description of how RCs achieve attractor reconstruction using the formalism of generalized synchronization, like in Lu et al.,¹¹ to the case of multifunctionality. Furthermore, we also aim to investigate the role of untrained attractors in how RCs achieve attractor reconstruction in more general cases. From an industry and applications perspective, the study of untrained attractors is an important issue for critical safety applications when, for instance, constructing digital twins of real-world systems as the operator needs to ensure that no unintentional dynamics can be triggered when noise is introduced. Moreover, the existence of untrained attractors raises questions about what behavioral assurances can be guaranteed for closed-loop RCs of this form.

ACKNOWLEDGMENTS

This work was funded by the Irish Research Council Enterprise Partnership Scheme (Grant No. EPSPG/2017/301). We thank the reviewers of this paper for their very helpful comments. We would also like to thank Andrea Ceni, Lou Pecora, Tom Carroll, Ling-Wei Kong, Joschka Herteux, Daniel Köglmayr, Oliver Heilmann, Christoph Räth, Jacob Morra, and many fellow researchers from UCC, for their influential conversations and input when discussing the contents of this paper as it has evolved over the past few years.

AUTHOR DECLARATIONS

Conflict of Interest

The authors have no conflicts to disclose.

Author Contributions

Andrew Flynn: Conceptualization (equal); Data curation (equal); Formal analysis (equal); Funding acquisition (equal); Investigation (equal); Methodology (equal); Validation (equal); Writing – original draft (equal); Writing – review & editing (equal). Vassilios A. Tsachouridis: Funding acquisition (equal); Supervision (supporting); Writing – review & editing (supporting). Andreas Amann: Conceptualization (supporting); Formal analysis (supporting); Funding acquisition (equal); Investigation (supporting); Methodology (supporting); Supervision (lead); Validation (equal); Writing – original draft (supporting); Writing – review & editing (equal).

DATA AVAILABILITY

The data that support the findings of this study are available from the corresponding author upon reasonable request.

APPENDIX A: RC DESIGN AND TRAINING PARAMETERS

In this paper, $M$ is constructed with an Erdös–Renyi topology where each of the non-zero elements are then replaced with a random number between $- 1$ and $1$ ⁠, the matrix is subsequently scaled to a specific spectral radius, $ρ$ ⁠. To elaborate, $M$ is designed such that each element is chosen independently to be nonzero with probability $P$ (i.e., sparsity $= P$ or degree $= N / P$ ⁠) and these nonzero elements are chosen uniformly from $(- 1, 1)$ ⁠. This random sparse matrix is rescaled such that the spectral radius, which is the maximum of the absolute values of its eigenvalues, is $ρ$ ⁠. This topology is chosen in order to provide the RC with enough dynamical flexibility to solicit multistable dynamics in tasks requiring multifunctionality. Similarly, the input matrix, $W_{i n}$ ⁠, is designed such that each row has only one nonzero randomly assigned element, chosen uniformly from $(- 1, 1)$ ⁠.

The numerical results presented throughout this paper are generated using the same $M$ and $W_{i n}$ ⁠. From further analysis and results which are not shown here, we remark that there are relatively small quantitative changes to our results when using different initializations of $M$ and $W_{i n}$ however the main characteristics of our results remain similar.

APPENDIX B: SYMMETRY AND SEEING DOUBLE

It is important to note that in Figs. 9 and 10 there are certain symmetries present in the projected basins of attraction.

In Figs. 9, 10(a), and 10(b), the closed-loop RC approaches the same attractor whether it is initialized with its representation of either the point $(x^{'}, y^{'})$ or $(- x^{'}, - y^{'})$ ⁠. However, there is a slightly different situation in Figs. 10(c) and 10(d). Using the terminology developed by Herteux and Räth,²⁴ when the closed-loop RC is initialized with its representation of the point $(x^{'}, y^{'})$ then $\hat{r} (t)$ settles down to a particular fixed point, for example, ${FP}_{1}$ or ${FP}_{3}$ ⁠, if instead the closed-loop RC was initialized with its representation of $(- x^{'}, - y^{'})$ then $\hat{r} (t)$ settles down to the “mirror-attractor” of that fixed point, which in this case is ${FP}_{2}$ or ${FP}_{4}$ ⁠.

Despite breaking the symmetry in this closed-loop RC setup with the square terms in Eq. (4) in the design of $W_{o u t}$ ⁠, there are certain properties of the seeing double training data that allow for the emergence of mirror-attractors in a similar way to some of the results presented in Flynn et al.¹³ regarding the “symmetry kills the square” phenomenon. The manner in which these scenarios arise in the case of the seeing double problem are outlined in the remainder of this Appendix.

1. Seeing double and $W_{o u t}^{(2)} \to 0$

In order for mirror-attractors to exist, the symmetry breaking terms in Eq. (10), i.e., the $r^{2} (t)$ term in Eq. (4) and the corresponding $W_{o u t}^{(2)}$ term in Eq. (5), cannot exist.²⁴

In the case of the seeing double problem, there are two conditions which contribute to forcing $W_{o u t}^{(2)} \to 0$ that need to be satisfied simultaneously. These are the anti-symmetric nature of the individual circular orbits, as described by Eqs. (28) and (29) in Flynn et al.,¹³ and the subsequent anti-symmetric nature of the combined training data set described in Eqs. (16) and (17) in Flynn et al.¹³ which determine $W_{o u t}$ ⁠.

The first of the conditions identified above, which relates to Eqs. (28) and (29) in Flynn et al.,¹³ holds when the RC is trained to reconstruct only

C_{A}

or

C_{B}

for

x_{c e n} = 0

as

u_{(C_{A})} (t) = - u_{(C_{A})} (t + \frac{T}{2}) & u_{(C_{B})} (t) = - u_{(C_{B})} (t + \frac{T}{2}),

(B1)

is satisfied

\forall t

so long as the open-loop RC responds to this driving input such that

r_{(C_{A})} (t) = - r_{(C_{A})} (t + \frac{T}{2}) & r_{(C_{B})} (t) = - r_{(C_{B})} (t + \frac{T}{2}),

(B2)

where

T

is the period of oscillation on each circular orbit.

However, complications arise with the second condition that relates to Eqs. (16) and (17) in Flynn et al.¹³ and concerns the structure of the resulting

X_{C}

and

Y_{C}

matrices used in the case of multifunctionality. On one hand, this condition holds

\forall x_{c e n}

when

C_{A}

and

C_{B}

are rotating in the same direction as the RC is effectively being trained to reconstruct mirror-attractors. To be more specific, the input data,

u (t)

⁠, from, for instance,

C_{A}

when centered at

(x_{c e n}, 0)

⁠, is equal to

- u (t)

when

C_{A}

is centered at

(- x_{c e n}, 0)

⁠, denoted as

C_{A}^{-}

⁠. As the RC is trained to reconstruct a coexistence of

C_{A}

and

C_{B}

for a given

x_{c e n}

⁠, and

C_{B} = C_{A}^{-} = - C_{A}

then

u_{(C_{B})} (t) = - u_{(C_{A})} (t)

and

\begin{aligned} (\begin{matrix} r_{(C_{B})} (t) \\ r_{(C_{B})}^{2} (t) \end{matrix}) = (\begin{matrix} r_{(C_{A}^{-})} (t) \\ r_{(C_{A}^{-})}^{2} (t) \end{matrix}) = (\begin{matrix} - r_{(C_{A})} (t) \\ {(- r_{(C_{A})} (t))}^{2} \end{matrix}) = (\begin{matrix} - r_{(C_{A})} (t) \\ r_{(C_{A})}^{2} (t) \end{matrix}) . \end{aligned}

(B3)

Therefore, the resulting

X_{C}

and

Y_{C}

training data matrices are of the form described in Eqs. (16) and (17) in Flynn et al.¹³ and as a consequence,

W_{o u t}^{(2)} = 0

⁠.

On the other hand, when $C_{A}$ and $C_{B}$ are rotating in opposite directions there is no mirror symmetry present in the combined training data sets, it is only when $x_{c e n} = 0$ that “symmetry kills the square” as the combined training data matrices from both $C_{A}$ and $C_{B}$ satisfy necessary conditions. As a further comment, if instead the RC was trained to reconstruct a coexistence of $C_{A}$ and $C_{B}$ when these are of different radii and $x_{c e n} = 0$ ⁠, so long as Eqs. (28) and (29) in Flynn et al.¹³ are satisfied for both $u_{(C_{A})} (t)$ and $u_{(C_{B})} (t)$ then $W_{o u t}^{(2)} = 0$ as the RC is trained to provide a coexistence of attractors which are all anti-symmetric to themselves. In return, this places no constraint on the behavior of the corresponding response attractors, $S_{A}$ and $S_{B}$ ⁠, in $S$ as these do not necessarily need to be mirror-attractors of one another as already shown in Fig. 9.

2. ρ vs $W_{o u t}^{(2)}$

Given the insight gained about the role of $ρ$ in the behavior of $W_{o u t}^{(2)}$ from Fig. 7 in Flynn et al.,¹³ the same procedure is applied here for the present case of the seeing double problem when $x_{c e n} = 0$ and the result of this is illustrated in Fig. 22(a). There is a threshold value of $ρ$ associated with the symmetry-breaking bifurcation discussed in Sec. V A 4 which corresponds with $W_{o u t}^{(2)}$ coming into existence in Fig. 22(a). Since $W_{o u t}^{(2)} = 0$ for all $ρ$ values less than this critical value, each solution of the closed-loop RC is required to either have a mirror-attractor or is anti-symmetric to itself, whether or not multifunctionality is achieved. While the behavior of $W_{o u t}^{(2)}$ in Fig. 22(a) is qualitatively similar to Fig. 7 in Flynn et al.¹³ there are some differences, mainly that $W_{o u t}^{(2)}$ reappears at smaller $ρ$ and its elements grow much faster and larger in magnitude over a similar interval of $ρ$ ⁠.

FIG. 22.

Histograms of the W o u t ( 2 ) elements for ρ ∈ [ 0.1 , 2.5 ] when in (a) x c e n = 0, in (b) x c e n = − 1, and in (c) x c e n = − 5.

View large Download slide

Histograms of the $W_{o u t}^{(2)}$ elements for $ρ \in [0.1, 2.5]$ when in (a) $x_{c e n} = 0$ ⁠, in (b) $x_{c e n} = - 1$ ⁠, and in (c) $x_{c e n} = - 5$ ⁠.

Nevertheless, the existence of mirror-attractors in Figs. 10(c) and 10(d), which led toward this discussion in the first place, and that $W_{o u t}^{(2)} = 0$ for the corresponding $ρ$ values shown in Fig. 22(a), together show that while both driving inputs possess the symmetry described in Eq. (B1) and the open-loop RC may respond to these driving inputs according to Eq. (B2), it is still possible that $W_{o u t}^{(2)} = 0$ even if the resulting attractors in $P$ are not periodic but are instead mirrored fixed points.

When $x_{c e n} \neq 0$ the symmetry in the training data is broken, which is much like shifting the location of the mirror-attractors as in Flynn et al.¹³ The effect this has on the elements of $W_{o u t}^{(2)}$ for $ρ \in [0.1, 2.5]$ is shown in Figs. 22(b)–22(c). Similar to the results depicted in Fig. 2 in Flynn et al.,¹³ the central message of Figs. 22(b)–22(c) is that as the magnitude of $x_{c e n}$ increases so too does the magnitude of $W_{o u t}^{(2)}$ ’s elements.

3. ${(W_{o u t}^{(1)})}_{x_{c e n}} = {(W_{o u t}^{(1)})}_{- x_{c e n}}$ and ${(W_{o u t}^{(2)})}_{x_{c e n}} = - {(W_{o u t}^{(2)})}_{- x_{c e n}}$

While it is shown in Fig. 22(a) that $W_{o u t}^{(2)} = 0$ for $x_{c e n} = 0$ ⁠, there is another “symmetry kills the square” phenomenon and related to the structure of the seeing double training data that is generated as follows.

Building on the statements made earlier which led toward the relationship established in Eq. ,

u_{(C_{A})} (t) = - u_{(C_{A}^{-})} (t)

and the resulting input training data matrix,

Y

⁠, for a given

x_{c e n}

and

- x_{c e n}

are related in the following way:

Y_{x_{c e n}} = - Y_{- x_{c e n}} .

(B4)

Therefore, in order for the training to be successful and

\hat{C_{A}} = - {\hat{C_{A}}}^{-}

⁠, then the following is also required:

{(W_{o u t})}_{x_{c e n}} (\begin{matrix} r_{(C_{A})} (t) \\ r_{(C_{A})}^{2} (t) \end{matrix}) = - {(W_{o u t})}_{- x_{c e n}} (\begin{matrix} r_{(C_{A}^{-})} (t) \\ r_{(C_{A}^{-})}^{2} (t) \end{matrix}),

(B5)

where

{(W_{o u t})}_{x_{c e n}}

is the resulting

W_{o u t}

matrix trained using data from

C_{A}

and

C_{B}^{-}

when centered at

(x_{c e n}, 0)

and

(- x_{c e n}, 0)

⁠, respectively, for a given

x_{c e n}

⁠, and similarly for

{(W_{o u t})}_{- x_{c e n}}

with

C_{A}^{-}

and

C_{B}

for the corresponding

- x_{c e n}

⁠.

Using the information from Eq. , this allows Eq. to be rewritten as

{(W_{o u t})}_{x_{c e n}} (\begin{matrix} - r_{(C_{A}^{-})} (t) \\ r_{(C_{A}^{-})}^{2} (t) \end{matrix}) = - {(W_{o u t})}_{- x_{c e n}} (\begin{matrix} r_{(C_{A}^{-})} (t) \\ r_{(C_{A}^{-})}^{2} (t) \end{matrix}),

(B6)

from which it follows that

\begin{aligned} - {(W_{o u t}^{(1)})}_{x_{c e n}} r_{(C_{A}^{-})} (t) + {(W_{o u t}^{(2)})}_{x_{c e n}} r_{(C_{A}^{-})}^{2} (t) \\ = - {(W_{o u t}^{(1)})}_{- x_{c e n}} r_{(C_{A}^{-})} (t) - {(W_{o u t}^{(2)})}_{- x_{c e n}} r_{(C_{A}^{-})}^{2} (t), \end{aligned}

(B7)

\begin{aligned} ⟹ ({(W_{o u t}^{(2)})}_{x_{c e n}} + {(W_{o u t}^{(2)})}_{- x_{c e n}}) r_{(C_{A}^{-})}^{2} (t) \\ = ({(W_{o u t}^{(1)})}_{x_{c e n}} - {(W_{o u t}^{(1)})}_{- x_{c e n}}) r_{(C_{A}^{-})} (t), \end{aligned}

(B8)

\begin{aligned} ∴ {(W_{o u t}^{(1)})}_{x_{c e n}} = {(W_{o u t}^{(1)})}_{- x_{c e n}}, \\ and {(W_{o u t}^{(2)})}_{x_{c e n}} = - {(W_{o u t}^{(2)})}_{- x_{c e n}} . \end{aligned}

(B9)

The above analytical results in Eq. (B8) are confirmed numerically in Fig. 23 by examining the behavior of the resulting ${(W_{o u t})}_{x_{c e n}}$ matrices after training the RC with data from $C_{A}$ and $C_{B}$ for $x_{c e n} = 3$ in Fig. 23(a) and $x_{c e n} = - 3$ in Fig. 23(b) where $ρ = 1.45$ in both cases. The elements of the resulting matrix from the element-wise addition of ${(W_{o u t})}_{x_{c e n} = 3}$ and ${(W_{o u t})}_{x_{c e n} = - 3}$ is depicted in Fig. 23(c) which shows that ${(W_{o u t}^{(2)})}_{x_{c e n}} + {(W_{o u t}^{(2)})}_{- x_{c e n}} = 0$ ⁠. Similarly, in Fig. 23(d), the element-wise subtraction of ${(W_{o u t})}_{x_{c e n} = - 3}$ from ${(W_{o u t})}_{x_{c e n} = 3}$ results in ${(W_{o u t}^{(1)})}_{x_{c e n}} - {(W_{o u t}^{(1)})}_{- x_{c e n}} = 0$ ⁠.

FIG. 23.

View large Download slide

${(W_{o u t})}_{x_{c e n}}$ elements with $ρ = 1.45$ for $x_{c e n} = 3$ in (a) and $x_{c e n} = - 3$ in (b). Result of element-wise addition of ${(W_{o u t})}_{x_{c e n} = 3} + {(W_{o u t})}_{x_{c e n} = - 3}$ in (c) and element-wise subtraction of ${(W_{o u t})}_{x_{c e n} = 3} - {(W_{o u t})}_{x_{c e n} = - 3}$ in (d). Each column represents the weight, $w_{i, j}$ ⁠, given to the $j^{t h}$ component of $q (r (t)) = {(r (t), r^{2} (t))}^{T}$ ⁠. The distinction is made between which columns belong to $W_{o u t}^{(1)}$ and $W_{o u t}^{(2)}$ ⁠. Note, the matrix plots may vary depending on the PDF viewer used, it is recommended to use “Adobe Reader” or “Chrome PDF Viewer” to see these images as intended.

While these results are interesting from an analytical and numerical perspective, their influence has been a common factor throughout many of the results presented in this paper and arise due to the constraint in Eq. (B8) as imposed by the structure of the training data.

REFERENCES

1.

K.

Nakajima

and

I.

Fischer

,

Reservoir Computing

(

Springer

,

2021

).

Google Scholar

Crossref

2.

D.

Verstraeten

,

B.

Schrauwen

, and

D.

Stroobandt

, in Proceedings of the 16th Annual ProRISC Workshop (Veldhoven, The Netherlands, 2005), pp. 454–459.

3.

H.

Jaeger

, “Bonn, Germany: German national research center for information technology,” GMD Technical Report 148 (2001).

4.

W.

Maass

,

T.

Natschläger

, and

H.

Markram

,

Neural Comput.

14

,

2531

(

2002

).

https://doi.org/10.1162/089976602760407955

Crossref

PubMed

5.

D.

Sussillo

,

P.

Nuyujukian

,

J. M.

Fan

,

J. C.

Kao

,

S. D.

Stavisky

,

S.

Ryu

, and

K.

Shenoy

,

J. Neural Eng.

9

,

026027

(

2012

).

https://doi.org/10.1088/1741-2560/9/2/026027

Crossref

PubMed

6.

T.

Waegeman

,

F.

wyffels

, and

B.

Schrauwen

,

IEEE Trans. Neural Networks Learn. Syst.

23

,

1637

(

2012

).

https://doi.org/10.1109/TNNLS.2012.2208655

Crossref

7.

D.

Canaday

,

A.

Pomerance

, and

D. J.

Gauthier

,

J. Phys.: Complexity

2

,

035025

(

2021

).

https://doi.org/10.1088/2632-072X/ac24f3

Crossref

8.

S. H.

Lim

,

L.

Theo Giorgini

,

W.

Moon

, and

J. S.

Wettlaufer

,

Chaos

30

,

123126

(

2020

).

https://doi.org/10.1063/5.0023764

Crossref

PubMed

9.

L.-W.

Kong

,

H.-W.

Fan

,

C.

Grebogi

, and

Y.-C.

Lai

,

Phys. Rev. Res.

3

,

013090

(

2021

).

https://doi.org/10.1103/PhysRevResearch.3.013090

Crossref

10.

H.

Jaeger

and

H.

Haas

,

Science

304

,

78

(

2004

).

https://doi.org/10.1126/science.1091277

Crossref

PubMed

11.

Z.

Lu

,

B. R.

Hunt

, and

E.

Ott

,

Chaos

28

,

061104

(

2018

).

https://doi.org/10.1063/1.5039508

Crossref

PubMed

12.

A.

Flynn

,

V. A.

Tsachouridis

, and

A.

Amann

,

Chaos

31

,

013125

(

2021

).

https://doi.org/10.1063/5.0019974

Crossref

PubMed

13.

A.

Flynn

,

J.

Herteux

,

V. A.

Tsachouridis

,

C.

Räth

, and

A.

Amann

,

Chaos

31

,

073122

(

2021

).

https://doi.org/10.1063/5.0055699

Crossref

PubMed

14.

A.

Flynn

,

O.

Heilmann

,

D.

Köglmayr

,

V. A.

Tsachouridis

,

C.

Räth

, and

A.

Amann

, in 2022 International Joint Conference on Neural Networks (IJCNN) (IEEE, Piscataway, NJ, 2022), pp. 1–8.

15.

J.

Herteux

, “The influence of the activation function on reservoir computers,” Master’s thesis (Ludwig-Maximilians-Universität München, 2021).

16.

K. L.

Briggman

and

W. B.

Kristan

,

J. Neurosci.

26

,

10925

(

2006

).

https://doi.org/10.1523/JNEUROSCI.3265-06.2006

Crossref

PubMed

17.

S.

Lieske

,

M.

Thoby-Brisson

,

P.

Telgkamp

, and

J.

Ramirez

,

Nat. Neurosci.

3

,

600

(

2000

).

https://doi.org/10.1038/75776

Crossref

PubMed

18.

I. R.

Popescu

and

W. N.

Frost

,

J. Neurosci.

22

,

1985

(

2002

).

https://doi.org/10.1523/JNEUROSCI.22-05-01985.2002

Crossref

PubMed

19.

G. J.

Mpitsos

and

C. S.

Cohan

,

J. Neurobiol.

17

,

517

(

1986

).

https://doi.org/10.1002/neu.480170512

Crossref

PubMed

20.

P. A.

Getting

,

Annu. Rev. Neurosci.

12

,

185

(

1989

).

https://doi.org/10.1146/annurev.ne.12.030189.001153

Crossref

PubMed

21.

P. S.

Dickinson

,

Curr. Opin. Neurobiol.

5

,

792

(

1995

).

https://doi.org/10.1016/0959-4388(95)80108-1

Crossref

PubMed

22.

E.

Marder

and

R. L.

Calabrese

,

Physiol. Rev.

76

,

687

(

1996

).

https://doi.org/10.1152/physrev.1996.76.3.687

Crossref

PubMed

23.

K. L.

Briggman

and

W.

Kristan

, Jr.,

Annu. Rev. Neurosci.

31

,

271

(

2008

).

https://doi.org/10.1146/annurev.neuro.31.060407.125552

Crossref

PubMed

24.

J.

Herteux

and

C.

Räth

,

Chaos

30

,

123142

(

2020

).

https://doi.org/10.1063/5.0028993

Crossref

PubMed

25.

J. Z.

Kim

,

Z.

Lu

,

E.

Nozari

,

G. J.

Pappas

, and

D. S.

Bassett

,

Nat. Mach. Intell.

3

,

316

(

2021

).

https://doi.org/10.1038/s42256-021-00321-2

Crossref

26.

L.-W.

Kong

,

H.

Fan

,

C.

Grebogi

, and

Y.-C.

Lai

,

J. Phys.: Complexity

2

,

035014

(

2021

).

https://doi.org/10.1088/2632-072X/ac0b00

Crossref

27.

L.-W.

Kong

,

Y.

Weng

,

B.

Glaz

,

M.

Haile

, and

Y.-C.

Lai

,

Chaos

33

,

033111

(

2023

).

https://doi.org/10.1063/5.0138661

Crossref

PubMed

28.

H.

Fan

,

L.-W.

Kong

,

Y.-C.

Lai

, and

X.

Wang

,

Phys. Rev. Res.

3

,

023237

(

2021

).

https://doi.org/10.1103/PhysRevResearch.3.023237

Crossref

29.

J.

Jiang

and

Y.-C.

Lai

,

Phys. Rev. Res.

1

,

033056

(

2019

).

https://doi.org/10.1103/PhysRevResearch.1.033056

Crossref

30.

N.

Bertschinger

and

T.

Natschläger

,

Neural Comput.

16

,

1413

(

2004

).

https://doi.org/10.1162/089976604323057443

Crossref

PubMed

31.

J.

Boedecker

,

O.

Obst

,

J. T.

Lizier

,

N. M.

Mayer

, and

M.

Asada

,

Theory Biosci.

131

,

205

(

2012

).

https://doi.org/10.1007/s12064-011-0146-8

Crossref

PubMed

32.

C.

Teuscher

,

BioSystems

00

,

104693

(

2022

).

https://doi.org/10.1016/j.biosystems.2022.104693

Crossref

33.

T. L.

Carroll

,

Chaos

30

,

121109

(

2020

).

https://doi.org/10.1063/5.0038163

Crossref

PubMed

34.

I.

Tsuda

,

E.

Koerner

, and

H.

Shimizu

,

Prog. Theor. Phys.

78

,

51

(

1987

).

https://doi.org/10.1143/PTP.78.51

Crossref

35.

K.

Ikeda

,

K.

Otsuka

, and

K.

Matsumoto

,

Prog. Theor. Phys. Suppl.

99

,

295

(

1989

).

https://doi.org/10.1143/PTPS.99.295

Crossref

36.

K.

Kaneko

,

Physica D

41

,

137

(

1990

).

https://doi.org/10.1016/0167-2789(90)90119-A

Crossref

37.

K.

Inoue

,

K.

Nakajima

, and

Y.

Kuniyoshi

,

Sci. Adv.

6

,

0000

(

2020

).

https://doi.org/10.1126/sciadv.abb3989

Crossref

38.

Z.

Lu

and

D. S.

Bassett

,

Chaos

30

,

063133

(

2020

).

https://doi.org/10.1063/5.0004344

Crossref

PubMed

39.

C.

Grebogi

,

E.

Ott

, and

J. A.

Yorke

,

Phys. Rev. Lett.

50

,

935

(

1983

).

https://doi.org/10.1103/PhysRevLett.50.935

Crossref

40.

C.

Grebogi

,

S. W.

McDonald

,

E.

Ott

, and

J. A.

Yorke

,

Phys. Lett. A

99

,

415

(

1983

).

https://doi.org/10.1016/0375-9601(83)90945-3

Crossref

41.

S. W.

McDonald

,

C.

Grebogi

,

E.

Ott

, and

J. A.

Yorke

,

Physica D

17

,

125

(

1985

).

https://doi.org/10.1016/0167-2789(85)90001-6

Crossref

42.

A. H.

Nayfeh

and

B.

Balachandran

,

Applied Nonlinear Dynamics: Analytical, Computational, and Experimental Methods

(

John Wiley & Sons

,

2008

).

Google Scholar

43.

I.

Tsuda

,

Curr. Opin. Neurobiol.

31

,

67

(

2015

).

https://doi.org/10.1016/j.conb.2014.08.011

Crossref

PubMed

44.

E. J.

Doedel

,

T. F.

Fairgrieve

,

B.

Sandstede

,

A. R.

Champneys

,

Y. A.

Kuznetsov

, and

X.

Wang

, “Auto-07p: Continuation and bifurcation software for ordinary differential equations” Technical Report. (Concordia University, 2007); available at http://cmvl.cs.concordia.ca/auto/.

45.

J.

Morra

,

A.

Flynn

,

A.

Amann

, and

M.

Daley

, arXiv:2306.01885 (2023).

46.

K.

Doya

, in Proceedings 1992 IEEE International Symposium on Circuits and Systems (IEEE, 1992), Vol. 6, pp. 2777–2780.

47.

R. D.

Beer

,

Adapt. Behav.

3

,

469

(

1995

).

https://doi.org/10.1177/105971239500300405

Crossref

48.

R. D.

Beer

,

Biol. Cybern.

116

,

501

(

2022

).

https://doi.org/10.1007/s00422-022-00938-5

Crossref

PubMed

49.

A.

Ceni

,

P.

Ashwin

, and

L.

Livi

,

Cognit. Comput.

12

,

330

(

2020

).

https://doi.org/10.1007/s12559-019-09634-2

Crossref

Blue	Correct cycle is reconstructed and $δ_{r e l} (C) < δ_{r e l}^{+}$ ⁠.
Yellow	Prediction switches to the other cycle.
Magenta	Prediction does not settle down to periodic motion
	and label this as some attractor, $S$ ⁠.
Black	Prediction tends to some other limit cycle
	or $δ_{r e l} (C) \geq δ_{r e l}^{+}$ ⁠, this attractor is labeled as, LC.
Green	Prediction decays toward some fixed point labeled as FP.

Seeing double with a multifunctional reservoir computer

I. INTRODUCTION

II. MULTIFUNCTIONAL RESERVOIR COMPUTERS

A. Introduction to multifunctionality

B. Reservoir computing

C. Training a RC to become multifunctional

D. Summary of previous work

III. SEEING DOUBLE

A. Numerical experiment setup

B. Illustrating issues of overlap

IV. SEEING DOUBLE WITH A MULTIFUNCTIONAL RC

A. Exploring multifunctionality in (xcen, ρ)-plane

1. Prediction analysis

2. Regions of multifunctionality

3. C A and C B rotating in the same direction

B. Transition to chaotic dynamics

1. Dynamics before chaos

2. Dynamics after chaos

C. Achieving multifunctionality when xcen = 0

1. Dynamics in P P

2. Dynamics in S S

D. Projected basins of attraction for xcen = 0

1. Method of generating the basins

2. Basins for ρ = 1.25

3. Basins for ρ = 0.9, 0.8, 0.4, 0.3

V. SEEING DOUBLE DYNAMICS

A. Bifurcation analysis for xcen = 0

1. Bifurcations of FP1 and FP2 to LC1, and of FP3 and FP4 to LC2

2. Bifurcations of LC1 and LC2 and the birth of C ^ B

3. The rise and fall of multifunctionality

4. Floquet analysis

B. Bifurcation analysis in the (xcen, ρ)-plane

1. Tracking stable solutions in the (xcen, ρ)-plane

2. General comments on RC’s dynamics in the (xcen, ρ)-plane

3. Behavior of FP3 and FP4

4. Behavior of FP1 and FP2

5. Behavior of LC1

6. Behavior of LC2

7. Behavior of C ^ B

VI. CONCLUSION

ACKNOWLEDGMENTS

AUTHOR DECLARATIONS

Conflict of Interest

Author Contributions

DATA AVAILABILITY

APPENDIX A: RC DESIGN AND TRAINING PARAMETERS

APPENDIX B: SYMMETRY AND SEEING DOUBLE

1. Seeing double and W o u t ( 2 ) → 0

2. ρ vs W o u t (2)

3. ( W o u t ( 1 ) ) x c e n = ( W o u t ( 1 ) ) − x c e n and ( W o u t ( 2 ) ) x c e n = − ( W o u t ( 2 ) ) − x c e n

REFERENCES

Citing articles via

Resources

Explore

pubs.aip.org

Connect with AIP Publishing

This Feature Is Available To Subscribers Only

A. Exploring multifunctionality in (x_cen, ρ)-plane

3. $C_{A}$ and $C_{B}$ rotating in the same direction

C. Achieving multifunctionality when x_cen = 0

1. Dynamics in $P P$

2. Dynamics in $S S$

D. Projected basins of attraction for x_cen = 0

A. Bifurcation analysis for x_cen = 0

1. Bifurcations of FP₁ and FP₂ to LC₁, and of FP₃ and FP₄ to LC₂

2. Bifurcations of LC₁ and LC₂ and the birth of ${\hat{C}}_{B}$

B. Bifurcation analysis in the (x_cen, ρ)-plane

1. Tracking stable solutions in the (x_cen, ρ)-plane

2. General comments on RC’s dynamics in the (x_cen, ρ)-plane

3. Behavior of FP₃ and FP₄

4. Behavior of FP₁ and FP₂

5. Behavior of LC₁

6. Behavior of LC₂

7. Behavior of ${\hat{C}}_{B}$

1. Seeing double and $W_{o u t}^{(2)} \to 0$

2. ρ vs $W_{o u t}^{(2)}$

3. ${(W_{o u t}^{(1)})}_{x_{c e n}} = {(W_{o u t}^{(1)})}_{- x_{c e n}}$ and ${(W_{o u t}^{(2)})}_{x_{c e n}} = - {(W_{o u t}^{(2)})}_{- x_{c e n}}$