{'fraction_non_alphanumeric': 0.02753013737033922, 'fraction_numerical': 0.008662741799831791, 'mean_word_length': 4.585812715314751, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Model-free inference of unseen attractors: Reconstructing phase space features from a single noisy trajectory using reservoir computing

Reservoir computers are powerful tools for chaotic time series prediction. They can be trained to approximate phase space flows and can thus both predict future values to a high accuracy, as well as reconstruct the general properties of a chaotic attractor without requiring a model. In this work, we show that the ability to learn the dynamics of a complex system can be extended to systems with co-existing attractors, here a 4-dimensional extension of the well-known Lorenz chaotic system. We demonstrate that a reservoir computer can infer entirely unexplored parts of the phase space: a properly trained reservoir computer can predict the existence of attractors that were never approached during training and therefore are labelled as unseen. We provide examples where attractor inference is achieved after training solely on a single noisy trajectory.

readouts can also provide universal computational properties. 3,4 It was soon discovered that reservoir computers are highly capable at predicting chaotic time series, improving accuracy by orders of magnitude compared to previous methods. 5 By feeding the reservoir a chaotic time series and choosing the target to be the next point of the chaotic trajectory, the reservoir computer learns to perform so-called one-step-ahead prediction. Intriguingly, after training for one-step-ahead prediction, such a reservoir can be put into a 'closed loop configuration,' where its own prediction can be used as the next input. In this way, the reservoir computer becomes an autonomous dynamical system that can accurately continue the original data. Generally, in a chaotic attractor, the predicted time series does not agree with the specific original training data forever, as even tiny differences will eventually lead to divergence of trajectories.
However, the autonomously created time series can reproduce the general structure of the chaotic attractor it was trained on even for long time scales. This is sometimes referred to as the 'climate prediction' property, as opposed to the short-term prediction error, which analogously is referred to as 'weather prediction' property. It is worth emphasizing that the reservoir is reproducing the source system in a model-free manner, i.e., without any knowledge about the origin of the training data.
Chaotic time series prediction using reservoir computing continues to be an active field of research with several noteworthy results. Pathak et al. 6 showed that properly trained reservoir computers can be used to reconstruct the largest Lyapunov eigenvalues inside the chaotic attractor. Chen et al. 7 showed that such reservoir computers also reproduce several typical geometric metrics. Estebanez et al. 8 demonstrated how noise can be used to improve the long-term attractor reconstruction and Zhu et al. 9 showed how delayed-feedback control can be combined with reservoir computing to find the unstable periodic orbits embedded inside the chaotic attractor.
Furthermore, unstable fixed points can be found using reservoir computing, even when the trajectory never visits them during training. 4, 10 Kim et al. 11 recently showed that a reservoir computer can correctly predict bifurcations if the training is done with the bifurcation parameter as an input explicitly fed to the reservoir. Remarkably, they only used training trajectories taken from time series data generated below the bifurcation point to teach the reservoir the influence of the bifurcation parameter. This all hints at the ability of a properly trained reservoir computer to reproduce the phase space structure of the original system accurately and even its control parameter dependence when in autonomous mode.
We demonstrate that the ability of phase space reconstruction can be exploited even further. We show that a properly trained reservoir computer cannot only reconstruct the attractor it is trained on -but it can even infer other unseen co-existing attractors in phase space and reconstruct their structure. By unseen co-existing attractors we specifically mean those attractors in a multistable system, whose basin of attraction was never reached during training so that no direct traces of them can be found in the original time series data. To this end, we test the ability of a reservoir computer to infer unseen attractors for a four-dimensional chaotic extension of the Lorenz system with coexisting attractors in two scenarios. The first scenario comprises a torus solution co-existing with a pair of symmetric limit cycles. In the second scenario a pair of symmetric chaotic attractors coexists with a torus.
Attractor reconstruction is arguably related to the field of nonlinear system identification, which aims to derive the governing equations from a sample set of time series data. In particular, techniques such as Sparse Identification of Nonlinear Dynamics (SINDy), 12 ResNet-based approaches, 13 auto-regressive moving average (ARMA) models and variants (ARMAx and NAR-MAx) are capable of producing accurate model descriptions in some cases. 14 However, all of these either require previous knowledge about the general structure of the system or multiple sets of training trajectories. Moreover, those that only estimate the vector field would have to be combined with an integrator to facilitate one-step-ahead prediction.
In contrast, we show that reservoir computers, operated in autonomous mode, can serve as model-free surrogates of target systems even outside of their training region with minimal input, i.e., they can learn to reconstruct unseen attractors after learning from a single, noisy trajectory.
Thus, reservoir computers can satisfactorily mimic target systems in cases where the training data is noisy, not a lot of training data is available, and nothing concrete about the shape of the underlying differential equation is known. Here, the reservoir computers do not learn the governing equations of the original system: Instead, they learn how to integrate and propagate along trajectories. Thus, the reservoir computers are learning the phase space flow without formulating any intermediate model. From recent studies it is known that attractor reconstruction can sometimes fail even for parameter-sets optimized for one-step-ahead prediction. 15,16 This also applies to the cases presented here, where a reservoir computer infers the existence of other attractors. Lacking a way to verify the prediction on the original system, this failure mode is harder to detect because the original data does not contain any sampling from the unseen attractors. Obtaining a quantitative estimate of the quality of the attractor reconstruction without knowing the target system remains an open problem.
In the following, we first introduce the model used and the target system chosen. We then analyze the reservoir computers' ability to reproduce unseen attractors. We find that for a set of parameters the reservoir computer succeeds in predicting the existence of a torus and symmetric limit cycle with low errors. In a second scenario, involving two chaotic attractors and a torus, it is more difficult to succeed. We analyze cases of partial success and discuss what the current limitations are.

A. Echo State Network
We use a continuous time version of an echo state network based on ordinary differential equations similar to those used by Lu et al. 17 The reservoir computer consists of a network of N coupled real-valued nodes with the total state X ∈ R N describing all nodes evolving according tȯ where the random matrices W res ∈ R N×N and W in ∈ R N×U describe the internal and the input weights, respectively. The input gain is given by scalar G. The time-varying input u(t) ∈ R U is a step-function with period θ , with the amplitude of each plateau corresponding to one point of data of the discrete source series u k ∈ R U×K . B ∈ R N is a random bias vector. The network and input matrices are sparse and the details of the random initialization and the list of all parameters are given in Appendix A.
The state X of the network is periodically recorded at intervals θ and used to construct the state matrix S ∈ R K×N+1 . We add a bias column into the state matrix such that a row S k is finally given by where t k = T w + kθ with T w being the 'washout' or 'warmup' time used to remove the influence of the starting state of the reservoir. The prediction Y ∈ R K×U of the reservoir is then given by A training series of K training elements is used to drive the system to generate S training , which, in turn, is used to determine the output weights W out ∈ R N+1×U . The error between the prediction Y and the known true targetsŶ are minimized using an L 2 -norm. To reduce the amount of over-fitting, we use Tikhonov regularization as described in Appendix C. The regularization strength η is an important parameter as it controls the level of detail the training data is approximated to. High η will fit the training data too crudely, rendering attractor inference impossible. Conversely, low η will lead to over-fitting and poor generalization.
Once the reservoir computer is successfully trained, we use it to probe other parts of the phase space of the target system. To this end, we feed the reservoir computer the beginning of a different ground-truth transient, which for now is required. We then observe how it evolves in autonomous mode. See Appendix C for details.

B. Target system by Li and Sprott
We study the properties of unseen attractors and the ability of a reservoir computer to infer their existence. As our target dynamical system, we use a 4-dimensional extension of the Lorenz attractor as proposed by Li and Sprott, 18 based on earlier work by Gao and Zhang, 19 whose dynamics is given byẋ which has only two parameters a and b. Here, we add additive independent identically distributed Gaussian noise terms ξ i with mean 0 and unit variance, E[ξ i (t)ξ j (t )] = δ i j δ tt . For sake of simplicity, the different variables share the same noise strength σ . This simplification is acceptable because all variables share the same order of magnitude. For σ = 0, the system is noise-free and deterministic. We chose this system because it shows coexisting attractors of different complexities depending on parameters. This allows us to test the ability of reservoir computers to infer the existence of unseen attractors in various scenarios involving periodic, quasiperiodic and chaotic attractors.
By initializing the system (4)- (7) in different initial conditions, we can reach different attractors, as each trajectory will eventually reach one of the stable attractors. The volume of phase space from which trajectories lead to a certain attractor is called its basin of attraction. Notably, the system (4)-(7) never features any fixed points, and hence the attractors are called hidden, and the basins of all attractors studied in this work are fractal. 18

C. Error estimates
We use a quantitative measure for the quality of attractor reconstruction. For the reconstructed and target attractors, we take the time average for each variable x , y , z , u , and the time average of the absolute values |x| , |y| , |z| , |u| . The absolute values help differentiate when averages become 0, such as for periodic states centred around the origin.
We determine the differences ∆ between prediction x and the ground-truthx, normalizing by the average ground-truth absolute averages and similarly for y, z, and u. From this, we calculate an error estimate ∆ att for the target attractor as the root of the sum of all squares While the quantitative error ∆ att does not capture the full picture, it allows us to easily discriminate between (partially) successful and failed attractor inference. To judge the power of a particular reservoir computer for a given scenario, we sum ∆ att for all existing target attractors to the total error ∆ tot via The applicability of other geometric error measures can also be considered, in particular the symmetric Hausdorff distance and the average Euclidean distances per point between the real and inferred attractor. However, the former is unsuited because it only measures the worst point. The latter suffers from requiring too many data points. Our sets sample the chaotic attractors too thinly, poorly covering rare trajectories, whose points in turn dominated the average Euclidean distance.
In contrast, the simple error measure ∆ att of Eq. (10) uses averages and is quite robust, even for small data set sizes.
We optimize the meta-parameters of the reservoir computers with the help of this error measure ∆ att . In this context, a meta-parameter is any parameter of the reservoir, i.e., a parameter that does not refer to the Li-Sprott oscillatory system. Some meta-parameters differ between simulations for the chaotic versus the limit cycle scenarios. We obtain optimal values via a rough grid search for The starting values for the transients used for training and details of numerical integration of the source system are described in Appendix B. We always train on trajectories with 11, 000 data points that stay within one basin of attraction. Of these, 1, 000 points are used as initial warm-up time T w .
We then train a continuous time echo state network of Eq. (1) with N = 300 on the Li-Sprott system with the phase space and attractors as shown in Fig. 1. We designate one of the two limit cycles as our training attractor and initialize a transient that converges to this limit cycle. Using this transient, we drive the reservoir computer, recording the activation of all the nodes. As per standard approach for attractor reconstruction, our training target is a one-step ahead prediction.
The optimal output matrix W out is found and used subsequently.
For our goal of having a system that can find unseen attractors, we use the autonomous operation of the reservoir computer. For this, we feed the reservoir its own prediction. In our case, additional care must be taken to avoid a training failure. First, it is typical that the beginning of the training data does not directly translate into the state matrix S training . In fact, we want the state X of our reservoir to be only dependent on the drive signal u(t). Thus, to avoid any influence of the initial state X(0) of the reservoir on the training, the first T w of u(t) is used to 'wash out' or 'warm up' the reservoir. We use 1, 000 of the 11, 000 training data points for this. But the training limit cycle is highly attracting. Therefore, any reasonably long 'wash-out' period T w leads to the loss of all information contained in the initial transient of the Li-Sprott system. The reservoir will effectively be trained only on the target limit cycle itself.
Second, when the training data is too low dimensional, the reservoir does not get to 'see' the full shape of the dynamical flow. Even for relatively short T w , the remaining training data only consists of the 1-dimensional limit cycle. Even when training succeeds in achieving a low fitting error on the training set, the system is susceptible to being unstable in the autonomous mode. Because certain phase space directions might never be seen during training, the output of the reservoir computer in autonomous mode tends to diverge in those directions.
Fortunately, these problems can be overcome with a modification that makes the training both more stable and makes the whole procedure more appropriate for realistic problems: by the inclusion of noise in the training data. When σ > 0, the training data produced by the Li-Sprott system (4)- (7) is in a sense permanently transient. Even small noise is enough to allow the trajectory to at least partially diverge from the pure limit cycle, effectively sampling all directions in the 4-dimensional phase space. While not explored in detail, there is likely a trade off between training length, noise strength and regularization that influences the performance of the reservoir.
We ensure that the noise does not lead to any attractor hopping in the training data set by visual inspection. We also tried using additive noise combined with a noise-free time series as has been studied previously 8 , but did not achieve success, likely due the fact that we are not only targeting chaotic attractors.
In comparison to Kim et al., 11 we do not use multiple parallel training trajectories. Instead we only train on a single, noisy trajectory. Also, compared to their work, we use a much coarser sampling of the original time-series as they used every point of the source data without any sampling and even included the Runge-Kutta auxiliary terms.
We always feed the reservoir all 4 variables of the Li-Sprott system both during training and the start of the autonomous mode, i.e., in this work there are no unseen degrees of freedom. Therefore, the reservoir obtains the full information about the target system in each time step. In principle, this should imply that no memory is needed and that removing the recurrent connections by setting W res = 0 in an approach similar to an extreme learning machine could suffice. 20 We performed a few preliminary investigations with W res = 0 and find that performance worsened. However, a detailed study of such an approach is outside the scope of this work. In general, we believe that the results presented in this manuscript can be extended to partially observed target systems in which case memory would be necessary. Furthermore, our chosen architecture is compatible with experimental reservoir computing, where the memory of the reservoir is an important aspect.
Inferring the torus and symmetric limit cycle Figure 2 shows a case of successful training and attractor inference: The system is trained on a single noisy trajectory of with K training = 10 4 points and noise strength σ 6 × 10 −3 . For testing, we feed the system the first 1, 000 points of a ground-truth reference trajectory (see Appendix C for details) for each of the three attractors (top, middle and bottom row in Fig. 2) and then put it into autonomous mode generating 10, 000 points on its own. We always use noise free (σ = 0) transients for initializing the reservoir in the autonomous mode and as targets.
The red dots in Fig. 2 represent the reconstructed attractor by the reservoir, while the black dots show the true attractors of the Li-Sprott system, the latter are almost completely covered by the former. The system learns the attractor it is trained on (Fig. 2 a, b) with a low error of ∆ att = 8.4×10 −3 . Moreover, it can also predict with high fidelity the existence and full shape of the second, symmetric limit cycle (∆ att = 7.6 × 10 −3 ) and of the torus (Fig. 2 e, f, ∆ att = 5.6 × 10 −2 ).
The total error in this case was ∆ tot ≈ 0.06.
All attractors have complex shapes. It had been shown that a reservoir computer trained on a chaotic system can implicitly learn the position of surrounding and defining unstable fixed points and the implicit unstable periodic orbits within. 4,9,10 Figure 2 clearly shows a case where we can go beyond this by demonstrating how an autonomous reservoir computer can infer entirely separate attractors in a target system. In particular, because the attractors of an autonomous system are always disjoint, this means that the system has to infer the properties of parts of the phase space that were neither part of the training data nor are connected to it via any trajectories, unlike unstable structures which might be connected via heteroclines. Here, the reservoir computer must make a plausible inference of the global phase space flow of the dynamical system from local knowledge.
It is important to note that there is quite a variation in quality of the attractor reconstruction depending on both the target as well as the randomly generated reservoir. In the following, we address a more difficult parameter region of the Li-Sprott system and demonstrate how partial successes and inference failure manifest themselves in such cases.

Training in chaos and quasiperiodic dynamics
Here, we change the target of the training. Using the Li-Sprott system with a = 6.0 and b = 0.1, we enter a regime where two symmetric chaotic attractors coexist with a torus. Once again, we train the system on one of the attractors using a single noisy trajectory. We find that attractor reconstruction and inference in this regime are significantly harder. In particular, we never achieve a total error ∆ tot < 2, which is almost two orders-of-magnitude larger than in the previous case.
When training on one of the chaotic attractors, the system has to infer the existence of its symmetric counterpart and a torus. Despite extensive numerical investigation, we do not find a case where the reservoir computer fully succeeds. Figure 3 shows one of the most successful results with an error of ∆ tot = 2.4: the system is able to learn the chaotic region it is trained on (Fig. 3 a, b, ∆ att = 0.14) and it correctly infers the existence and type of the two other attractors, as depicted in Fig. 3 c-f. These other reconstructed attractors show obvious deviations. For the reconstructed symmetric chaotic region, the dynamics of the x, y, and z variables are, in fact, almost accurately predicted with errors ∆ i at most 0.18. However, the u variable shows a persistent positive off-set from the ground truth, as can be seen in Fig. 3 c, e, and this is reflected by higher errors ∆ u ≈ −∆ |u| ≈ 0.43. Similarly, the torus (Fig. 3 e, f) is mismatched in the u-dimension with ∆ u ≈ 0.7. In addition, the shape is distorted increasing the error of the average absolute values ∆ |x| ≈ 0.6, ∆ |y| ≈ 1.0, ∆ |z| ≈ 1.9 and ∆ |u| ≈ −0.3.
The deviation in the u variable for the unseen attractors is likely due to the large separation between the training and target attractor: the two chaotic regions are clearly separated from the torus and from each other in the u-variable. We assume that this large separation makes an accurate inference harder as compared to the case of the intertwined limit cycles and torus studied earlier. The reservoir is required to correctly infer the phase space flow over longer distances, which accumulates the visible deviations between ground truth and prediction in Fig. 3 e, f. Furthermore, the different time scales cannot be ignored in this case. While the limit cycle and torus dynamics of the previous section allowed for a sampling that was reasonably adequate for all oscillations involved, this is not the case here. In particular, the chaotic dynamics show both slow oscillations on a time scale of 100 and fast oscillations on the order of 3. Similarly, the torus contains slow oscillations around 6 and fast oscillations of order 120. This is a separation of almost two orders-of-magnitude which makes choosing a suitable sampling rate a difficult compromise.
We chose a sampling interval of 0.2, mostly adapted to the fast oscillations. Nevertheless, even in this more difficult case, the system correctly predicts the type, shape and rough position of the unseen attractors. If knowing the exact u position and torus shape is not required, the reservoir successfully predicts the unseen attractors.
We find many cases of partial success and failure. Two types of obvious failures in particular are detectable even in a model-free setting. First, if the reservoir fails to reproduce the training data set, it is safe to assume that its predictions cannot be considered dependable. We see this type of failure remarkably often, and the fact that some randomly generated reservoir topologies work worse than others is a known issue. 15  We show the subset of simulations with the lowest total error ∆ tot in Fig. 4. The histogram shows a clear bimodal distribution. While the different predictions do not cleanly fall into classes, after investigating the source of this bimodality, we conjecture that the first peak is related to predictions that approximate the targets as good as possible. The second peak corresponds to cases in which inference fails and only the training region is predicted to be stable, i.e., the multistability can not be inferred.
Ultimately, the error measures ∆ always compress information. When there is a significant deviation between ground truth and prediction, many cases of partial success arise with different qualities. We show an additional example to highlight how partial success can manifest itself. Figure 5 shows an example where inference deviates from the ground truth by one attractor exhibiting a different type and another not being detected. Figure 5 and Fig. 3 only differ in their randomly generated topology. As in the successful trial, the reservoir succeeds in reproducing the training attractor (Fig. 5 a, b) with a low error of ∆ att = 0.09 and therefore, barring additional information, both the attractor inference of Fig. 5 and Fig. 3 appear equally valid. Nevertheless, the prediction of the two unseen attractors is less successful for Fig. 5: The total error is slightly larger with ∆ tot ≈ 2.5. Inference of the torus is not fully achieved and instead a limit cycle is predicted Similarly, we can sometimes observe cases where a torus or limit cycle is predicted with slightly distorted shape. Furthermore, this limit cycle roughly follows the outline of the torus, indicating that it correctly reflects some truth about the underlying phase-space. It might also be possible, that this limit cycle exists in the original system but is only slightly unstable.
In contrast, the inference of the chaotic region fails for Fig. 5 c, d with ∆ att = 1.5. The reservoir does not even detect the existence of an attracting structure and converges to the same attractor as  For now, we have no method for preventing this kind of failure. If a method could be devised that can judge the quality of inferred attractors, one could use this to discriminate between successful and failing reservoirs. One potential avenue is a scenario for which one does know two attractors of the target system: In this case, one can use one of them as a training data set, and another as an independent test. For systems with only a single training data set, independent error estimation remains an open question for now.
A final open question regards the role of noise in the training data: minute details of the phase space flow that might be important for long range inference are washed out by the noise. We therefore tried to reduce the noise as much as possible. However, a certain level of noise is necessary for our scheme to succeed. Systematically investigating the relationship between necessary noise levels (as compared to computational noise levels), training series length, inference distance and regularization strength are interesting challenges for future work to give insights into these important aspects.

IV. DISCUSSION
We report the successful use of an autonomously operated reservoir computer for attractor reconstruction. This includes not only the reconstruction of the training attractor, but also of attractors that the system never sees during training. Moreover, we perform training using only a single noisy trajectory. This demonstrates that reservoir computers are powerful tools for modelfree system analysis and attractor prediction, as there are no assumptions put into the reservoir about the target system.
We show how the system is able to correctly predict the existence and shape of a pair of symmetric intertwined complex limit cycles and a torus in a 4-dimensional extension of the Lorenz system. We then explore more difficult scenarios involving two chaotic regions and a torus. For this problem, we obtain partial success. And even for cases when we cannot find reservoirs that fully succeeded in attractor inference, some useful partial information is still extracted such as the approximate location of an attractor, or just the existence of 'some kind' of attractor.
Several open questions remain to further develop this approach: First, and foremost, the fraction of randomly generated reservoirs that fail is quite high. For non-optimized parameters, it can be more than 95%. Even for optimized networks, the percentage remains high: For the parameters used in Fig. 2 about 1500  Over 1300 even had a deviation of ∆ i ≥ 100 in at least one variable indicating that the trajectory failed to converge. We speculate, that this due to larger separation in phase space of the attractors, as well as the larger difference in time scales of the involved oscillations. The widely separated time scales might be better tackled and the prediction errors alleviated by using different reservoirs adapted to the different time scales. This remains a question for future studies.
Another remaining question concerns the full understanding of the role of noise. 14 In our scheme, noise in the source system helps to explore all directions in phase space and enables the system training to succeed. However, noise also results in a destruction of fine details. Finding the best trade-off between noise and accuracy remains a future challenge, which will give important hints as to the kind of information the reservoir is exploiting when inferring the existence of unseen attractors. Conversely, it will be fruitful to study how to achieve improved training on noise-free data or data with minimal noise.
We speculate that the Ridge regression parameter likely plays a crucial role. It controls how much importance the reservoir assigns to small differences in the training data, and to what extent it tries to reproduce those. In fact, the Ridge regression parameter might control the 'model complexity': training with strong regularization should lead to simpler models, while training with weak regularization allows for more complexity, but is also more susceptible to overfitting. This also relates to the influence of finite precision and its role in attractor inference.
Taking a broader view, the reservoir has to perform its task based on a finite amount of noisy data. Thus, the amount of information that can be extracted about 'unseen parts' of the phase space should also have some fundamental limit. We expect that a more complex target dynamical system will also require a longer training sequence for a reservoir to be able to emulate it. Questions about the relationship of data length, target system complexity and related fundamental limits will give further insights about how close to an 'optimal' prediction can be made by a particular reservoir computer.
Furthermore, in a real-world scenario, a way to fully explore the phase space without the need of initial transients from other regions is vital. We suspect that a rough initialization with artificial data such as constant values is possible but have not tested this proposition.
With the examples shown in our work, we provide a proof-of-principle what can be achieved.
A reservoir computer can definitely infer the existence of unseen attractors with varying degrees of success. As such, this further proves their suitability as model-free substitutes for a target system.

DATA AVAILABILITY STATEMENT
The data that support the findings of this study are available from the corresponding author upon reasonable request.
Appendix A: Construction of the random echo state network The echo state networks used in this work always consists of N = 300 real-valued nodes, whose vectorial state X evolves according tȯ The real-valued reservoir weight matrix W res , input weight matrix W in and bias vector B are random but fixed and initialized as follows: The bias vector B consists of N = 300 uniformly distributed random numbers in [−b, b]. W in is a 300 × 4 dimensional matrix for our system. Each of the N = 300 rows has only a single non-zero entry, whose column is randomly chosen with equal probabilities. The real-valued entry itself is drawn from a uniform random distribution in [−1, 1] for each row. The strength of the input weights is globally controlled via the parameter G.
Finally, the reservoir weight matrix W res is a sparse matrix with sparsity ρ. We used ρ = 0.1 in all presented simulations. Larger ρ were tried but did not yield any immediately noticeable differences, and as sparse matrices are faster to simulate, we chose the smallest acceptable ρ. The non-zero values are initially distributed in [−1, 1]. After creation, we calculate the largest real part of the eigenvalues Re(Λ) max of the resulting matrix W sparse . We then divide the entire matrix and multiply by a margin factor, such that the resulting reservoir weight matrix W res possesses a maximal real part of any eigenvalue with Re(λ ) max , that are related through In our simulations, Re(λ ) max is either 0.95 or 0.99, depending on what is found to be better. To reach the two limit cycles for a = 2.0, b = 0.8, we start the time series at (±5, ±1, 1 ± 1).

Appendix C: Reservoir computing technical details
We integrate (1) with a fourth order Runge-Kutta integrator. The state of all nodes X n ∈ R is initialized as 0 in each element. The integration time step is 0.1. We let the echo state network evolve without input for 300 time units or 3000 steps, and then switch to the "warm up" or "wash out" procedure with input. The input sequence u k is fed as a piece-wise constant function u(t) ∈ R 4 with interval lengths θ . We use an interval length of θ = 2.5 found through meta-parameter scans, albeit the exact length does not critically influence the performance. Each component of u(t) corresponds to one of the four variables of the Li-Sprott-variant of the Lorenz system (4)-(7).
We use the first 1, 000 points encoded in u(t) to remove any influence of the starting state, i.e., T w = 2, 500. After that, we start recording the system state for 10 4 inputs, i.e., 25, 000 time units, in accordance to (2), i.e., every θ time units the state of all 300 nodes is recorded. We always use the last integration point in each piece-wise constant interval of u(t) to record the maximally large reaction to the input. We obtain a state matrix S with 9, 999 rows and N + 1 = 301 columns.
We create a target vector Y with 9, 999 rows and 4 columns. Each column contains one variable of the input u k shifted one step into the future. Because both input and output are 4-dimensional, the output weight matrix W out is of size (N + 1) × 4. We then solve the following equation for where η is the regularization factor of the Tikhonov regularization and I N+1 the identity matrix of appropriate size N + 1. S T is the transpose of the state matrix.

Autonomous operation
For the autonomous operation, we use the output weights W out learned during training. First, we reset the system state X to the point it was after the first 300 time units of input-free evolution.
We then feed the first 1, 000 ground truth points of a noise-free transient u k leading to one of the three target attractors. We observe the system state X at the end of this washout-period, and construct state vector S i (X) ∈ R 301 that looks like recording the row of a state matrix (2). But instead of saving its entries, we directly multiply with W out , which yields 4 values corresponding to the 4 dimensions of the predicted output stream. We treat these 4 values as the next element of a self-generating input sequenceũ k and keep track ofũ k . With the new response X generated, we repeat the process for 10, 000 steps. The seriesũ k then is the reconstructed attractor as inferred or learned by the echo state network.