A Hopf physical reservoir computer

Physical reservoir computing utilizes a physical system as a computational resource. This nontraditional computing technique can be computationally powerful, without the need of costly training. Here, a Hopf oscillator is implemented as a reservoir computer by using a node-based architecture; however, this implementation does not use delayed feedback lines. This reservoir computer is still powerful, but it is considerably simpler and cheaper to implement as a physical Hopf oscillator. A non-periodic stochastic masking procedure is applied for this reservoir computer following the time multiplexing method. Due to the presence of noise, the Euler–Maruyama method is used to simulate the resulting stochastic differential equations that represent this reservoir computer. An analog electrical circuit is built to implement this Hopf oscillator reservoir computer experimentally. The information processing capability was tested numerically and experimentally by performing logical tasks, emulation tasks, and time series prediction tasks. This reservoir computer has several attractive features, including a simple design that is easy to implement, noise robustness, and a high computational ability for many different benchmark tasks. Since limit cycle oscillators model many physical systems, this architecture could be relatively easily applied in many contexts.

the nodes without the presence of any delay or feedback line 38 . This approach is studied less, though it makes the reservoir architecture much simpler.
Here, a Hopf oscillator is used as a physical reservoir. The Hopf oscillator can also be used as the building block for adaptive oscillators 39,40 , which can natively learn information without any training. The Hopf oscillator can exhibit limit cycle motion, which provides a source of memory by storing information in its dynamic states. Although a binary periodic masking function is popularly used for time-multiplexed reservoir 6,10,12 , noise can also be used as periodic mask 41 . In this paper, a Hopf oscillator PRC is constructed that uses a non-periodic stochastic mask. A Hopf oscillator physical reservoir computer is fabricated as an analog circuit, which is compared with Euler-Maruyama simulations 40,42,43 . This Hopf PRC can successfully complete benchmark machine learning tasks, including parity tasks, fundamental logic gate tasks 12 , nonlinear dynamic emulation tasks 4 , and various time series prediction tasks 44 . The information rate is used as the performance metric for logical tasks 11 , and the normalized mean square error (NMSE) is used for the emulation and time series tasks 4 .
The rest of the article is organized as follows. In "System equations for Hopf physical reservoir computer" section, the equations of motion for the stochastic Hopf oscillator PRC are presented. In "Mapping methodology" section, the methodology of mapping the oscillator's dynamics to an information processing scheme is discussed for an example task by using the Euler-Maruyama simulation. The effects of the pseudo-period and the noise on computational ability are discussed in "Pseudo-period and noise" section. In "Analog circuit experiment" section, the analog circuit experiment is described. In "Benchmark tasks for Hopf PRC" section, different benchmark tasks are performed with the numerical and experimental Hopf PRC, which includes logic tasks, emulation tasks of time series, and prediction tasks. The concluding remarks are stated in "Concluding remarks" section.

System equations for Hopf physical reservoir computer
The equations of motion for the Hopf oscillator are 45 : For this Hopf oscillator, x and y are the first and second states, respectively, and the sinusoidal forcing is given by A sin(�t + φ) . A list of the parameters is given in Table 1. The information is first encoded as an input, u(t), which will depend on the benchmark task being performed. The mask is defined by white Gaussian noise as: Here, σ is the noise amplitude, Ẇ is white Gaussian noise, and β is a positive bias. It should be noted that Ẇ does not exist, but its differential form, dW, does 46 .
To send information to the PRC to be processed, an external forcing function that contains the information signal, u(t), and the stochastic mask, m(t), is constructed as: This external forcing function is injected into both the amplitude of the sinusoidal forcing, A, and the parameter affecting the limit cycle radius, µ . Including this force, the equations for the Hopf PRC are written as:

Mapping methodology
To use the dynamics of the Hopf oscillator as a physical reservoir computer, the dynamics must first be mapped. To describe this mapping, an exclusive OR (XOR) logical task is used as an example. In this section, the Hopf PRC is simulated using an Euler-Maruyama scheme, since the mask is stochastic 42 . Shannon's information metric is used to quantify the performance of the reservoir when performing logical tasks, such as the XOR operation 11,43 . For this task, the binary "false" and "true" values are encoded as discrete negative ones and positive ones, respectively, in a discrete signal, r(z). r(z) is defined such that z ∈ Z + and r(z) ∈ {−1, +1} , which is depicted in Fig. 1a. To input this into a continuous dynamical system, these values are first mapped to a continuous input function, u(t), as follows: This function is depicted in Fig. 1b. T p is a constant pseudo-period, in which u(t) does not change its value. Thus, for the XOR logical task, the input function, u(t) ∈ {−1, +1} , is a random square wave with a pseudo-period, T p . This implies that each of the "true" (e.g., +1) or "false" (e.g., −1) values affect the system for an amount of time, T p . The mask function, m(t), is depicted in Fig. 1c. The Hopf PRC system described in Eq. (4) is numerically integrated using the Euler-Maruyama (EM) method, since the PRC is stochastic 42,47 . For these simulations, the integration time step, dt = 10 −5 seconds, the total simulation time in this case was 3000T p = 300 seconds, and T p = 0.1 sec. This example simulation is shown in Fig. 1.
The time history of the x state obtained from the simulation is depicted in Fig. 1d. Next, x(t) is re-scaled by subtracting the mean, µ x , and dividing by the standard deviation, σ x , using Eq. (6): is used for the subsequent steps. The time history of the X state is depicted in Fig. 1e. Next, equidistant nodes are created by dividing each pseudo-period, T p , equally into N(= 20) nodes, as shown in Fig. 1f. Over each pseudo-period, T p , the N node values are referred to as the nodal state, which is depicted in Fig. 1g.
The node matrix, S, is an N × K matrix; for this example, N = 20 is the number of nodes over a pseudoperiod, and K = 3000 is the total number of pseudo-periods. Truncating the final 20% of this S matrix ( 600T p ), a new matrix, L ( 480T p ) is formed, which will be used in the training process. Throughout this paper, next, the reservoir computer is trained using ridge regression, as in Eq. (7): A target signal (the M vector) is created from the encoded input based on a benchmark task, which in this case is XOR task. For each pseudo-period, there will be one target value that is found by performing the XOR operation between the inputs, r(z) and r(z − 1) . In this way, the target vector, M, is found for the XOR task. Linear regression based training is then applied to the nodal state matrix, L, to map it to the desired output using Eq. (7). In Eq. (7), w is the weight vector found after training, I is the identity matrix, = 10 −1 is the regularization parameter used to avoid over-fitting, and o(k) is the prediction of the reservoir computer at the kth pseudoperiod. The discrete input, r(z) and continuous input, u(t) are given in Fig. 2a and b respectively. Figure 2c shows this prediction along with the corresponding target signal. In the final step, the prediction is binarized since XOR is a binary task, which is depicted in Fig. 2d. It should be noted that a nonlinear dynamic emulation task would not require this final step of discretization.
For a logical task, the efficacy of the reservoir computer is quantified using Shannon's information rate 48 . The information rate, R, can be defined as follows: Here H(x) is the Shannon entropy, which denotes how much information is encoded in a signal. This can be defined as follows: In this equation, p i is the probability of getting a particular bit, i. H y (x) is the conditional entropy, which denotes the probability of getting an incorrect bit in the target signal: Here p i (j) = p(j|i) = p(i,j) j p(i,j) and p(i, j) is the joint probability distribution of the two variables, i and j, each of which can take a value of "1" or "− 1" for a logical task. i is a bit from the target, and j is a bit from the prediction. For the simulation depicted here, the parameters were set such that: www.nature.com/scientificreports/ The information rate, R, for this case was calculated to be 0.98 based on the prediction from the validation portion (not including in the training process). Due to the nature of this binary target signal, the Shannon entropy is 1.0, which marks the maximum value of the information rate for this task. It should be noted that the lower limit of R is zero, which would be achieved if every prediction was incorrect, while the upper limit of R depends on the task. For the parity tasks considered here, the upper limit of R is equal to one.

Pseudo-period and noise
In this section, the effects of pseudo-period and noise on the computational ability of the reservoir are explored. For this discussion, several parity tasks (defined in Eq. (12)) are used to understand the effects of the pseudoperiod and noise on the computational ability of the reservoir. The relationship between the pseudo-period, T p , and the natural frequency of the oscillator, ω 0 , is explored in Fig. 3 by using the 2nd and 4th order parity tasks. In Fig. 3a-c, the reservoir computer's performance is measured for three different values of T p while varying the natural period, 2π ω 0 . It is found that the reservoir has better performance when the pseudo-period is an integer multiple of the natural period of the oscillator. The reservoir's performance is studied using both resonance ( ω 0 = � ) and non-resonance ( ω 0 = � ) conditions. It is found that both cases can result in strong or weak computational ability depending on the fractional relationship between the natural period and the pseudo-period. However, maintaining this design can still fail to make a robust reservoir computer when T p is very low (e.g., T p = 0.05 seconds). For the remainder of the paper, combinations of T p and ω 0 are chosen such that the pseudo-period is an integer multiple of the natural period of the oscillator.
Noise is ubiquitous in physical systems. For this reason, noise is introduced into this system using a stochastic masking function. Figure 4 shows the relationship between the computational ability, as measured with R, the noise amplitude, σ , and the noise bias, β. The simulations presented in Fig. 4 are performed for the 4th order parity task (left) and 6th order parity task (right). The reservoir is found to be robust against a certain level of noise intensity, which demonstrates its potential to be implemented under the influence of environmental noise. However, increasing noise intensity does decrease the computational ability of the reservoir. This effect may be observed for a higher order task, which requires a longer memory (e.g., the 6th order parity task of Fig. 4). When β = 0 , the computational ability was the lowest. Since the non-periodic noise mask with increasing noise intensity deteriorates the computational ability, it should be noted that the Hopf reservoir computer can also be built by excluding the noise mask ( σ = 0).

Analog circuit experiment
To build a physical reservoir computer (PRC), an analog circuit implementation of Eq. (4) was designed, fabricated, and tested. The circuit's equations are given in Eq. (11): www.nature.com/scientificreports/ Here, V u is the input voltage, V m is the stochastic masking voltage, V µ is the limit cycle radius voltage, V ω 0 is the resonance constant voltage. V x and V y are the states, which correspond to states x and y in Eq. (4). The circuit implementation used TL082 operational amplifiers and AD633 multipliers in standard integrator network configurations. The error tolerance is 1% for the resistors and 2% for the capacitors. The continuous input function, V u , the stochastic masking function, V m , and the sinusoidal forcing, sin(�t + φ) , were created in MATLAB and sent to the circuit via a National Instrument (NI) cDAQ-9174. This cDAQ-9174 also collected the V x and V y states. A sampling frequency of 10 5 samples/s was used to collect data for all the experiments. The resistor values were chosen such that R 1 = 10 k and R 2 = 100k , and the capacitor values were chosen such that C = 0.1µ F. A simplified schematic is shown in Fig. 5.
The V x state will be treated in the same manner that the x state was treated in "Mapping methodology" section. That is, the V x state will be rescaled using Eq. (6), and then the rescaled state will be used to form the nodal state matrix, L. The target signal vector, M, will be created following the same process discussed in "Mapping methodology" section. Finally, Eq. (7) will be used to train the PRC to map input data to the desired output values. As an example, the analog circuit Hopf PRC was used to solve the XOR task as in the previous section, which is depicted in Fig. 6. The information rate, R, for this case was calculated to be 1.0 based on the prediction from the validation portion (not including in the training process).

Benchmark tasks for Hopf PRC
The Hopf PRC is numerically and experimentally tested with three benchmark tasks: (1) logic tasks, (2) emulation tasks of time series, and (3) prediction tasks. Logic tasks include the fundamental logic gate tasks and parity tasks of different orders. Emulation tasks of time series will test the PRC's ability to reproduce nonlinear auto (11)  Logic benchmark tasks. Parity tasks. The computing efficacy of the reservoir is first evaluated with parity benchmark tasks. Since it is a logical task, the input function, u(t), is generated with a random binary signal, r(z), as discussed in "Mapping methodology" section. The nth order parity function, P n , is defined by the following equation: As n increases, this task will require more memory and nonlinearity from the reservoir. As given in "Mapping methodology" section, Shannon's information metric is used to measure the performance of the PRC for logic tasks. For n = 1 , the first-order task does not require any memory from the input of the previous pseudo-period, so the task is linear. For n > 1 , the task is nonlinear, which demands that the reservoir computer must also possess memory and the nonlinear separation ability. In Fig. 7, the ability of the Hopf PRC to follow parity tasks of 2nd to 5th order, both experimentally and in simulations. The initial 4000T p = 400 seconds are used for training, and the final 1000T p = 100 seconds are used for testing. The performance difference between the PRC experiment and the simulation could be due to the presence of nonlinear circuit components in the analog circuit, which are not represented in Eq. (11). For instance, V u must jump between −1 and +1 , but this instantaneous change takes a finite amount of time in the circuit.
Fundamental logic gate tasks. The computing performance of the reservoir is also assessed with fundamental logic gates: NOT ( ¬ ), AND ( ∧ ), and OR ( ∨ ). The input function, u(t), is generated with a random binary signal as discussed in "Mapping methodology" section, and the Shannon's information metric is used again to measure the performance of this PRC. Figure 8 depicts the response of the Hopf PRC acting as fundamental logic gates, both experimentally and in simulations. In all cases, the Hopf PRC achieved an information rate that was maximal.
Emulation tasks. The reservoir is also evaluated with emulation tasks. The nonlinear auto-regressive moving average (NARMA) time series is used to test whether the reservoir possesses adequate nonlinearity and long time lags 4,24,26,28 . These tasks show the multi-tasking capability of the reservoir. NARMA tasks from the 2nd to 20th orders are used to test the reservoir. A NARMA task of order n is given in Eq. (13), where the initial target values are set to 0.19: In Eq. (13), M n is the target of the system. n is the order of NARMA task, (f 1 , f 2 , f 3 ) = ( 2.11 500 , 3.73 500 , 4.33 500 ) , and (α, ζ , γ , δ) = (0.3, 0.05, 1.5, 0.1) 26,28 . u(t) is the continuous input that is used to force the Hopf PRC, which is www.nature.com/scientificreports/ a function of a three sinusoidal functions. It should be noted that this formulation of the NARMA emulation task is non-standard. The u(t) given in Eq. (13) was used for other dynamic systems in which inertia played a large role 4,26 . Similarly, this non-standard NARMA task is used here to evaluate this analog circuit reservoir. The reservoir emulates this nonlinear function, but it should be noted that the correlation present in Eq. (13) does not allow a definitive evaluation of the long-term memory characteristics of this reservoir.
In the simulations and experiments, t = 0.1 seconds, and the sampling rate was 10 5 samples/second. Figure 9 shows several NARMA tasks. Instead of the information rate, the normalised mean square error (NMSE) is used to evaluate the performance of the reservoir computer for the NARMA tasks:  Prediction tasks. Santa Fe task. Time series forecasting is an important benchmark for a reservoir. The Santa Fe time series was first used in a time series forecasting competition as a benchmark test. The Santa Fe time series data set A is a univariate time series found from the recorded intensity of a chaotic far-infrared-laser 49 . The target signal is generated to predict the value at the next time step based on the values of the current and previous time steps. Figure 11a shows the Hopf PRC's performance on this laser time series, for both the experiment and the numerical simulations. NMSE is used as the performance metric.   52 . The heart rate, chest volume (respiration force), and blood oxygen concentration comprise the target.
For each of these time series, the target signal is again generated to predict the next step based on the values of the current and previous time steps. In each case, the original time series is normalized to use as the input. Figure 11b-d shows the reservoir computer's performance in predicting subsequent values of the heart rate, respiratory force, and blood oxygen concentration, respectively, through both experiments and numerical simulations. The NMSE is calculated in each case to evaluate the performance of the reservoir.
Sunspot prediction task. The prediction of the total number of sunspots ( S n ) is also a one-step time series prediction task similar to the Santa Fe time series 10 . Daily and monthly total sunspot numbers were used in one step forecasting purpose by the reservoir computer. The necessary data set is taken from WDC-SILSO, Royal Observatory of Belgium, Brussels 53 . Again, for each of the time series, the target signal is generated to predict the next value based on the value of the current and previous time steps, and the original time series is normalized to use as the input to the oscillator. Figure 12 (top) shows the reservoir's performance in predicting the next steps of the daily total counted sunspots, and Fig. 12 (bottom) shows the performance in predicting monthly counted sunspots. Again, the NMSE is used to evaluate the reservoir's efficacy for this task.