{'fraction_non_alphanumeric': 0.0351840359001289, 'fraction_numerical': 0.007829283429608059, 'mean_word_length': 4.537404176579434, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 6, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Autonomous Learning by Dynamical Systems with Inertial or Delayed Feedbacks

Dynamical systems can autonomously adapt their organization so that the required target dynamics is reproduced. In the previous Rapid Communication [Phys. Rev. E 90,030901(R) (2014)], it was shown how such systems can be designed using delayed feedbacks. Here, the proposed method is further analyzed and improved. Its extension to adaptable systems, where delays are absent and inertial feedbacks are instead employed, is suggested. Numerical tests for three different models, including networks of phase and amplitude oscillators, are performed.


I. INTRODUCTION
In contrast to modern computers, information processing in animals and humans does not involve digital encoding and numerical computation. Instead, the analog mode is employedthe animal brain is essentially a dynamical system that can emulate various other dynamical systems in its environment and thus predict their behavior [1]. The question is whether one can construct dynamical systems that would be much more simple then neural networks, but nonetheless able to perform the tasks of emulation and prediction.
The current study is a continuation of previous work on engineering of dynamical systems where evolutionary methods have been employed for systems design. In 1988, it was proposed that, by simulated annealing, oscillatory dynamical systems with prescribed frequency spectra can be built [2]; subsequently this was confirmed in numerical simulations [3]. For flow distribution networks that represent an idealization of the signal transduction system of a biological cell, it was shown that, by using stochastic Metropolis algorithms, networks with prescribed static response patterns can be designed and made furthermore robust against various structural perturbations [4] or noise [5] (see also review [6]). These results could be further extended to flow network models with dynamic responses [7]. For coupled phase oscillators, it was shown that, by running evolutionary optimization, networks that are better entrained by pacemakers could be designed [8,9].
Model genetic oscillatory networks with prescribed oscillation periods could be designed by evolutionary optimization [10] and made robust against knock-outs of genes, removal of regulatory interactions, and introduction of noise [11]. Moreover, genetic regulatory networks with predefined adaptive dynamical responses [12,13] or giving rise to definite stationary expression patterns [14] could be constructed in a similar way. Stochastic Monte Carlo optimization methods with replica exchange were employed to design phase oscillator networks with improved synchronization tolerance against heterogeneities [15] or noise [16].
Application of evolutionary optimization methods in engineering of dynamical systems can be viewed as representing supervised learning: an external agent (actually, a digital computer) generates structural perturbations (i.e., the mutations), monitors their effects on system's performance and decides whether a particular mutation should be accepted or not.
In contrast to this, evolution leading to the designed dynamics should be performed entirely within the evolving system itself in the case of autonomous (or non-supervised) learning.
Through an adaptation process, the system learns to adapt its internal structure in such a way that the desired performance is reached.
In our earlier Rapid Communication [17], a general scheme of autonomous learning by dynamical systems with delayed feedbacks has been proposed and illustrated using networks of coupled phase oscillators. Now, we provide a detailed analysis of the previously proposed method, consider its possible improvements and extensions and give further application examples.

II. METHODS
We consider a system with variables x = (x 1 , x 2 , . . . , x N ) whose dynamics is governed by which depend on parameters w = (w 1 , w 2 , . . . , w K ). Suppose that this system is responsible in a selected application for execution of a function which needs to be optimized through the evolution of parameters w. This function can be quantified by a performance vector R = (R 1 , R 2 , . . . , R M ) and, as we assume, this vector can be expressed as a function (or, generally, a functional) of system's variables, i.e. as R = F (x). The ideal target performance of the system (1) is known and corresponds to a certain vector R 0 . Thus, the deviation ǫ of the actual performance from the target performance of the system, i.e. the optimization error, can be defined as ǫ = |R − R 0 | or, explicitly, as The error vanishes, ǫ = 0, when the system reaches the target.
Our aim is to introduce slow intrinsic dynamics of the parameters w which would lead to the minimization of error ǫ. To do this, we assume that parameters w are themselves dynamical variables and their temporal evolution is governed by stochastic differential delay equations [17] where ξ i (t) are independent random white noises with ξ(t) = 0 and ξ α (t)ξ β (t ′ ) = 2δ αβ δ(t − t ′ ). Hence, the system consists of the guided subsystem (1) which is controlled by the slow steering subsystem (3). The steering subsystem is persistently sensing the current and the previous performances of the guided subsystem and its deviations ǫ from the ideal target performance.
The evolution equation ( On the other hand, if the performance has become worse (δǫ > 0), the steering subsystem moves the parameters in the opposite direction. As the target performance is gradually approached, the effects of noise should be however reduced. Therefore, we choose S(ǫ) = S 0 ǫ, so that the noise strength is proportional to the deviation from the target performance.
In the previous publication [17], we assumed that the parameter τ remains constant during the evolution, i.e. we have taken τ (ǫ) = τ 0 . Such choice is not however optimal because the error variation δǫ = ǫ(t) − ǫ(t − ∆) within the delay time can be much smaller than the error itself and, as a result, diffusion can dominate over the drift. In the present study, we choose τ (ǫ) = τ 1 (ǫ + µ) where µ ≪ 1. Hence, the drift becomes accelerated when the error ǫ is decreased. The small parameter µ is introduced to prevent the divergence when the target performance is reached and ǫ = 0.
In the considered learning scheme, there should be a separation of time scales between the guided and the steering subsystems. According to equation (3), the characteristic time for the evolution of the parameters w is τ (ǫ)/δǫ = τ 1 (ǫ + µ)/δǫ. Because, during the evolution, we typically have δǫ ≪ ǫ, this implies that τ (ǫ)/δǫ ≫ τ 1 . When the characteristic time scale for the dynamics of the guided subsystem is chosen as unity, the separation of time scale is, on the average, ensured if we choose the parameter τ 1 to be of order one.
As noted already in the Rapid Communication [17], the scheme employing fixed time delays may be sometimes difficult to implement at the hardware level, because, when it is employed, the system needs to store in its memory the previous performances (i.e., the errors) and the previous parameter values. In such situations, learning schemes with inertial feedbacks can be used instead.
The inertial feedback learning scheme is formulated by further extending the system, so that it includes an additional set of parameters w´and an auxilliary error variable ǫ´. The evolution equations for parameters w are and the dynamics of the additional variables is determined by equations and These equations involve a characteristic time for the additionally introduced variables that is taken to be the same as the previous delay time ∆. The equation (1) for dynamical variables x is not changed.
To clarify the relationship between the inertial and delayed feedback versions, we note that equations (5) and (6) can be formally integrated, yielding and Therefore, for example, the auxilliary error ǫ ′ at time t is obtained by summing contributions from different delayed moments t − s taken with the exponentially decreasing weights exp(−s/∆). Hence, in contrast to the previous scheme with a fixed delay ∆, the inertial version effectively involves distributed time delays with the average equal to ∆.
To compare the operation and test the efficiency of different learning methods, three model problems will be used. Below we formulate these problems and specify the error functions and the parameter evolution equations in each case.
A. Learning to maintain a given steady state. Under strong environmental variations, the maintainance of a definite steady state may involve adaptation at the structural level, so that the parameters of a system become appropriately modified. As an illustration, we consider a simple problem of a system with gradient dynamics that learns to relax to a given stationary state.
The system is characterized by dynamical variables x = (x 1 , x 2 , . . . , x N ) and parameters w = (w 1 , w 2 , . . . , w K ). The dynamical variables change with time according to equations For simplicity, we assume that the potential function U is given by the parabolic form Hence, the dynamics represents linear relaxation to a stationary state with we have We assume that variations of x and the values of the parameters w are restricted to the interval (−R max , R max ).
Suppose that the task of learning is to change the parameters so that a given target steady sate R 0 is asymptotically reached. For this problem, the error function can be chosen as Thus defined, the error is always less than one.
Then, the evolution of parameters w is described by equations (3) in the scheme with the delayed feedbacks and by equations (4), (5) and (6) under the inertial learning scheme.
Although the solution of this simple optimization problem can be easily constructed, it is still interesting to see how an autonomous learning system would deal with it.
B. Learning to synchronize. As our second example, we choose the problem that has already been taken in the previous Rapid Communication [17]. The system represents a network of coupled phase oscillators and it has to learn under what set of connection weights the state with a predefined synchronization level is achieved.
The dynamical variables are the phases φ = (φ 1 , φ 2 , ..., φ N ) and the parameters are the connection weights w ij . The phase dynamics is described by the classical Kuramoto model where ω i is the natural frequency of an oscillator i. We assume that the interactions are symmetric, w ij = w ji and their weights can be positive or negative.
The synchronization is quantified by the order parameter Because it fluctuates with time, the degree of synchronization at time t is better characterized by an average over a time interval T , Note that R(t) can vary from zero to one and the value R(t) = 1 corresponds to the state with complete phase synchronization, when the phases of all oscillators are identical.
The task of learning in this case is to reach a structural state where synchronization at an arbitrary predefined level P with 0 < P < 1 is achieved. This has to be done by only rearranging the weights (and possibly changing their signs), whereas the total interaction strength N i,j=1 |w ij (t)| remains fixed. The error ǫ can be defined as When fixed delays are employed, the evolution equations for the connection weights are [17] dw ij (t) Here, a term is added to the right hand side of the equations to control the total strength of the interactions. The average absolute weight v(t) is defined as and this average weight should stay close to a given value W . The magnitude of the coefficient λ determines how strictly this condition should hold. In contrast to the previous study, the dependences τ (ǫ) = τ 1 (ǫ + µ) and S(ǫ) = S 0 ǫ will now be used.
When learning with inertial feedbacks is employed, evolution of connections weights is governed by equations and the dynamics of the auxilliary error variable ǫ ′ is given by The same definition of v(t) and the same dependences τ (ǫ) and S(ǫ) are then used. Note that the scheme still involves a short delay T needed for the collection of data to determine the running time-averaged synchronization order parameter R(t).
C. Learning to avoid the amplitude death. As the third example, a Kuramoto-

Suzuki network of interacting amplitude oscillators is chosen. Previous investigations have
shown [19] that such networks possess a rich dynamics. Particularly, amplitude death for a fraction of oscillations in networks with random connections could be observed. The task of learning may be to modify the structure of a network so that the fraction of oscillators undergoing the ampltude death is minimized.
In this model, the state of an oscillator j is described by its complex amplitude W j and the oscillator dynamics is governed by equations where the coefficients c 0 , c 1 and c 2 determine properties of individual oscillators and the parameter κ specifies the strength of coupling between them.
The coupling is characterized by the Laplacian matrix L with the elements Here, T is the adjacency matrix of the network whose elements are T ij = 1 if there is a link from j to i, and T ij = 0 otherwise. If a link exists, it is further characterized by its positive weight w ij . We assume that the adjacency matrix represents a random Erdös-Rényi symmetric network with connectivity p. It will be fixed during the evolution and only the connection weights w ij will be changed.
When the weights are all equal to unity, this system was previously investigated in Ref. [19]. It was found that, under certain conditions, some of the oscillators in the system undergo the amplitude death, i.e their oscillation amplitude drops down to (almost) zero.
Our aim is to distribute the connection weights over a given random network in such a way that the number of oscillators with the amplitude death is reduced.
We define the fraction f of oscillators with the oscillator death as where h is a fixed threshold and Θ(x) is the step function, i.e. Θ(x) = 1 for x > 0 and Θ(x) = 0 otherwise. This fraction fluctuates with time and, to define the error ǫ(t), averaging over a time interval T is performed, To ensure that the connection weights remain positive during the evolution, the following procedure is applied: We introduce additional variables g ij , such that w ij = g 2 ij . Thus, the weights stay positive whereas the variables g can change their signs. With this convention, the previous schemes for the autonomous learning can be applied.
We require that the average connection weight, is approximately conserved during the evolution and remains close to the unity. Here, L is the total number of connections in the considered fixed random network.
In the inertial learning scheme, the evolution of the variables g is governed by stochastic Note that, similar to our second model, we have included into these equations a term that imposes the approximate conservation of the average connection weight v(t). The parameter λ determines how strictly this condition is satisfied.
The evolution for the additional network variables g ′ ij and for the auxilliary error variable

IV. RESULTS
In this section, results of numerical simulations of autonomous learning for three problems, that were formulated above, will be presented and analyzed.
A. Learning to maintain a given steady state. This example will be used to discuss the effects of the dependence τ (ǫ) and to compare the efficiency of the delay and inertial schemes. In the simulations, we take a system with N = 5 variables and the target R 0 i = 10 for i = 1, 2, ..., 5. The initial conditions for variables x i are randomly chosen within the interval [−R max : R max ], with R max = 30. In the delay scheme, the initial conditions for the parameters an their memory are w i = x i . In the inertial scheme, we have w i = w ′ i = x i and ǫ ′ = ǫ at t = 0. In both cases the initial errors are computed with the initial configuration of the system. The delay parameter is ∆ = 10 and S 0 = 1. The system is integrated numerically with a stochastic Heun algorithm with the time step of 0.001.
Typical simulation results are displayed in Fig. 1. Here, the dependence of error ǫ on time is shown for the delayed and inertial feedback schemes with τ (ǫ) = τ 0 (Fig. 1 a and b) or τ (ǫ) = τ 1 (ǫ + µ) (Fig. 1 c and d). To facilitate the comparison, the same pseudo-random number series are used in all these simulations, so that the noise realizations are identical in all of them. Moreover, the same initial conditions are employed. We have τ 0 = 0.18, τ 1 = 1, and µ = 0.0001.
When the parameter τ remains constant during the evolution, in both schemes the system finds a relatively good solution at the time about 8 × 10 4 . However, after that the errors grow again and the system moves away from the target performance. This behavior is repeated again and again, so that the target performance is never stabilized. Obviously, such bouncing is caused by the effects of noise that become dominant when the errors are small. In contrast to this, the errors gradually decrease with time when the dependence τ (ǫ) = τ 1 (ǫ + µ) is employed and the target performance is asymptotically reached (Fig.   1c,d). The convergence to the target performance is somewhat better for the scheme with the delayed feedbacks.
These differences are further illustrated by Fig. 2 where the dependence of the weights on the time for the same four evolutions is shown. When τ (ǫ) = τ 0 , the weights approach the target values, but they cannot keep such values and the system soon escapes from the optimal state. On the other hand, when τ (ǫ) = τ 1 (ǫ + µ), the target values are found faster step dt = 0.01. As initial conditions, we set all the weights and memory as w ij = w ′ ij = W ; and the initial errors ǫ = ǫ ′ is computed with the initial system configuration. Finally, initial oscillator phases are randomly chosen between zero and 2π.
As seen in Fig. 3 where the dependence of the synchronization order parameter R on time is plotted, the system can learn to reproduce the required synchronization level of w ′ ij = w ij = 1 for the elements with T ij = 1. Finally, the errors ǫ ′ = ǫ at t = 0 are computed using the initial values. from the center of the plane, avoiding the amplitude death.

V. DISCUSSION
Dynamical systems can adapt their internal organization, e.g. the connection weights in a network, in such a way that a desired target dynamics is reproduced. The adaptation is performed through a built-in feedback loop. In this loop, the actual and the target performances of a system are compared and the resulting error signal is used to modify the structure of the system. Thus, through an autonomous evolution process, the system learns to generate a definite output.
In this study, an improved version of the previously proposed scheme with delayed feedbacks [17] was considered and a novel autonomous learning scheme with inertial feedbacks was formulated. Numerical tests for selected optimization problems have shown that both methods produce satisfactory results. Since the proposed methods are empirical, it would be important to check them for a broader range of applications too.
Although only their computer simulations have been so far performed, the proposed methods may also be implemented, partially or completely, at the "hardware" level. The guided dynamical system (1) can be of physical, chemical or biological origin. Its equation of motion and even a complete set of its dynamical variables may be unknown.
Essentially, such a real-world system can represent a "black box" whose output signals are controlled by varying a parameter set.
The steering system (2) should monitor the outputs and control the parameters of the guided system. It is not necessary that this is done by a computer in a digital way. Instead, this system can represent, for instance, an electronic circuit where the dynamics corresponding to the equations (2) is physically reproduced. Note that, in old analog computers, even the systems with much more complex dynamics could be implemented by electronic means.
Moreover, the two components, steering and guided, can also be incorporated into one physical device. Through an internal adaptation process, such a device would be able to generate the required outputs and to maintain these outputs despite large variations of environmental conditions or the occurrence of faults.
Finally, a question can be asked whether similar mechanisms are already used by biological systems, even so small as a single living cell. In the latter case, the steering slow component can be, for example, the genetic subsystem, whereas the guided fast component can represent the metabolic subsystem of a cell. This is an interesting topic for further research.
One of the authors (P.K.) acknowledges financial support from SeCTyP-UNCuyo (project