Insight into delay based reservoir computing via eigenvalue analysis

In this paper we give a profound insight into the computation capability of delay based reservoir computing via an eigenvalue analysis. We concentrate on the task-independent memory capacity to quantify the reservoir performance and compare these with the eigenvalue spectrum of the dynamical system. We show that these two quantities are deeply connected, and thus the reservoir computing performance is predictable by analysing the small signal response of the reservoir. Our results suggest that any dynamical system used as a reservoir can be analysed in this way. We apply our method exemplarily to a photonic laser system with feedback and compare the numerically computed recall capabilities with the eigenvalue spectrum. Optimal performance is found for a system with the eigenvalues having real parts close to zero and off-resonant imaginary parts.


Introduction
Reservoir computing is a novel approach for time-dependent tasks in machine learning. First introduced by Jaeger [1] and inspired by the human brain [2], it utilizes the inherent computational capabilities of dynamical systems. Very recently the universal approximation property has also been shown for a wide range of reservoir computers, which solidifies the concept as a broad applicable scheme [3].
A new and sophisticated approach to the reservoir computing scheme was introduced by Appeltant et al in [30], where a single dynamical node under the influence of external feedback utilizes a time-multiplexed reservoir. The spatially extended network structure of classical reservoirs is no longer needed with this scheme, which reduces the complexity in reservoir hardware in exchange for processing speed. A schematic sketch is shown in figure 1. Realizations with a single delayed reservoir [31][32][33][34][35][36] give a first glimpse over the potential of this idea for, e.g. time-series predictions [28,37], equalization tasks on nonlinearly distorted signals [38], and fast word recognition [39]. A general analysis, introduced by Dambre et al [40], was also used to quantify the task-independent computational capabilities of semiconductor lasers [41]. For an overview, we refer to [42][43][44].
A lot of research was already invested in order to develop a deeper understanding of reservoir computing systems, however, effective measures that allow to predict the performance are still missing. In this paper we want to fill this gap by providing a scheme that allows to predict general trends of the performance using the eigenvalue spectra of the dynamical system (the reservoir) without input. As an example reservoir, we chose a laser that is subjected to optical self-feedback. We use the Lang-Kobayashi system, which is an established model for a semiconductor laser with delayed external feedback. We calculate the total memory capacity as well as the linear and nonlinear contributions using the method derived in [40] and compare the results with the computed eigenvalue spectrum of the system, where we discover a clear connection. In particular, a high linear memory capacity is found for systems, where a large number of eigenvalues are close to criticality (with small negative real parts) and non-resonant (with imaginary parts not-resonant to the input timescale).
The paper is structured as follows. First, we give an overview of the methods used for calculating the memory capacity and the eigenvalue spectrum in section 2. After that, we present our results and discuss the impact of the eigenvalues on the performance and different nonlinear recall contributions first for a reservoir formed by a solitary laser and then by a laser with external cavity.

Methods
The reservoir computing scheme employs the idea of a dynamical reservoir, which projects input information into a high dimensional phase space. The nonlinear response of the reservoir is then used by a linear readout to approximate a specific task depending on the input. Often the reservoir consists of many nodes with relatively simple dynamics (for example, tanh-function [1]) in which the input enters via a weighted matrix. Afterward, the response is read out and linearly combined to generate an output. The idea is to minimize the Euclidean distance between the generated output and the target. This approach is particularly resourceful for time-dependent tasks, because the dynamical system which forms the reservoir has memory and thus acts as a memory kernel.
The modified approach introduced by [30] uses a single node with delay as a reservoir, in which the output dimensions are distributed over time. A mask g is used to feed the input into the system, which produces a high dimensional response. These responses are saved over time and used for the linear readout approximation. A sketch of the setup is shown in figure 1. In the following, we will give a short overview of the quantities and notations used in this paper. We also refer to our previous works [45,46], where a detailed explanation of how the reservoir setup is operated and task-independent memory capacities are computed is given.

Time-multiplexed reservoir computing
Let us briefly remind the main ingredients of the time-multiplexed reservoir computing scheme [30,45,46]. An input vector u ∈ R L enters the system componentwise at times t l = lT, l = 1, . . . , L, L being the number of sample points. The time between two inputs t l+1 − t l is called the clock cycle T and describes the period length in which one input u l is applied to the system. Inside each interval of one clock cycle T, a T-periodic mask function g is applied on the inputs (see figure 1). The mask g is piecewise-constant on N V intervals, each of length θ = T/N V corresponding to N V virtual nodes. The values of the mask function g play the same role as the input weights in spatially extended reservoirs, with the difference that the input weights are now distributed over time.
The system responses are collected in the state matrix S ∈ R L × R NV , where N V is the dimension of the measured system's state. More specifically, the elements of the state matrix are S ln = s(lT + nθ) with n = 1, . . . , N V , and l = 1, . . . , L, where s(t) ∈ R is the state of the dynamical element of the reservoir at time t, e.g. a variable of the delay system in simulations, or laser intensity in an experimental realization.
A linear combination of the state matrix is given by Sw, where w ∈ R M is a vector of weights. Such a combination is trained to find a least square approximation to some target vectorŷ: where ∥ · ∥ 2 is the Euclidean norm, and λ T is a Tikhonov regularization parameter. A solution to this problem is known to satisfy: when S T S + λ T I is invertible. In the case of our Lang-Kobayashi model, since the physical system is intrinsically noisy, we used the state noise regularization [1,47] and set λ T = 0. This is reasonable, as noise dominates very small dependencies in the given training data set, which the linear readouts would otherwise try to fit. It also gives a more realistic threshold for the precision of the state readouts compared to numerical precision. Comparisons of simulations without noise and with Tikhonov regularization to a noisy system without Tikhonov regularization yielded similar results. To quantify the system's performance, we use the normalized root mean square error (NRMSE) between the approximation y = Sw and the targetŷ: where var(ŷ) is the variance of the target valuesŷ = (ŷ 1 , . . . ,ŷ L ).

Memory capacity
Dambre et al have shown in [40] that the computational capability of a reservoir system can be quantified via an orthonormal set of basis functions on a sequence of inputs. Here we give a recap of the used quantities introduced in [45]. In particular, the capacity to fulfil a certain task is given by: The capacity equals 1 if y =ŷ and the reservoir computer computes perfectly the task; C = 0 if it can not compute it at all, and in between 0 and 1 if it is partially capable to fulfil the task. In section appendix A, we explain how equation (3) follows from the corresponding expression in [40]. Further, following Dambre et al [40], we use finite products of normalized Legendre polynomials P d as a basis of the Hilbert space of all possible transformations (thus tasks with targetsŷ) on an input sequence {u} = {u −L , . . . , u −3 , u −2 , u −1 }. As inputs into the system, we use uniformly distributed random numbers u l , which are independent and identically drawn in [−1, 1]. This yields uncorrelated inputs and thus uncorrelated memory capacities. After feeding the input sequence {u} of random numbers into the system, it yields a reservoir response S. Formally, the memory capacity (equation (3)) is defined for an infinitely long sequence L → ∞. To approximate it numerically, we use L = 250 000. In order to describe a task, the target vectorŷ is defined as: where {d} = {d 1 , ..., d I } is a sequence of degrees such that the Legendre polynomial P d i (u −i ) of degree d i is applied to the input u −i . The product of all such polynomials is used to generate the task (target vectorŷ). The collection of all tasks (4) for any possible degree sequence {d} is the Hilbert space of all possible transformations [40]. Further, to define the linear and nonlinear memory capacities, one uses special tasks, for which the sum of the degrees i d i is constant:ŷ Clearly, there are many such possible tasks for all sequences {d} with d = i d i . The memory capacity MC d of degree d is defined as the sum of the capacities Cŷ computed using equation (3) for all tasks (5) of degree d: The well known linear memory capacity corresponds to d = 1. The total memory capacity is then given by the memory capacities MC d of all degrees d.
It was shown in [40] that MC is limited by the readout-dimension N V , which equals the number of virtual nodes N V . An intuitive explanation is the following. The linear readout Sw of the reservoir computing scheme can be considered a linear combination of the columns of the state matrix S. Thus the amount of dimensions this basis can approximate is given by the number of linearly independent readouts. If the systems states are linearly independent, it can at most approximate N V different dimensions, which is in our case N V different tasks constructed from equation (4). A more rigorous explanation is given by Dambre et al in [40].

NARMA10
In addition to memory capacities, we evaluate the NRMSE of the NARMA10 task. NARMA10 is an often used benchmark test that combines linear and nonlinear memory transformations. It is given by the following iterative formula: A n−i + 1.5u n−9 u n + 0.1.
Here, A n is an iteratively given number and u n is an independent and identically drawn uniformly distributed random number in [0, 0.5]. The reservoir is fed with the random numbers u n and has to predict the value of A n + 1 .

Lang-Kobayashi model
We use the Lang-Kobayashi laser as an example reservoir. This is a model applicable for semiconductor lasers with external feedback operating with low feedback strength. The Lang-Kobayashi equations have been studied widely, modelling successfully semiconductor lasers [48,49] exhibiting complex dynamics and bifurcation scenarios [49][50][51]. The dimensionless equations of motion are given by [52]: The parameters scaling was chosen as in [53] with a modification to allow for the information input. The system time is normalized to the photon lifetime. Here, E is the complex electric field, N is the charge carrier inversion, and ξ describes spontaneous emission modelled by Gaussian white noise, g is the masking function, I is the input, α the amplitude-phase coupling, κ is the feedback strength, ϕ the feedback phase, τ the delay time, D noise the noise amplitude and P + ηI(t)g(t) is the pump current, composed of P a constant pump level and η the input strength of the information fed into the system via electric injection, which is small with respect to P. I(t)g(t) is the piecewise constant input function, which contains the data set values multiplied with a mask function. T LK is the time scale ratio, modelling class B laser behaviour for sufficiently large T LK seen in conventional semiconductor lasers. T LK ≪ 1 models class A behaviour observed in quantum dot lasers [54][55][56] or quantum cascade lasers [57,58]. Note that the threshold pump current for the solitary laser is at P th = 0, while P th changes with κ according to P th = −κ.

Calculating eigenvalue spectrum
The eigenvalues of any dynamical system describe the dynamics for small perturbations around a linearized point. Because reservoir computers are often operated close to a stable fixpoint the eigenvalue spectrum of its linearized system can be analysed. The goal of this paper is to find a relation between the nonlinear memory recallability and the eigenvalue spectrum. The latter can be computed with much less numerical effort and could then be used to predict good parameter ranges for reservoir setups. It also gives insight into the timescales of the eigendirections of the system, which contain information on the memory kernel of the reservoir. To compute the eigenvalue spectrum, we used two methods: the first method is an analytical approximation in the long delay limit [59], while the second relies on numerical computation with the DDE-biftool software package [60][61][62].
To begin with, we give a short overview of the first method from [59], which provides an approximation of the spectrum of long-delay systems. In [53], it was applied to the Lang-Kobayashi system. As delay-based reservoir computing is mostly used with a long delay compared to the local dynamics, this is a valid approximation that gives a general tool to analyse the reservoirs of such type.
The characteristic equation for the eigenvalues is obtained through the linearization around a steady state x * , and it reads as: with some constant matrices A and B, and I is the identity matrix. For large τ , its solutions can be decomposed into two parts, in which one scales as ℜ(λ) ∼ 1/τ , also called the pseudocontinuous spectrum, and a strongly unstable spectrum with the scaling ℜ(λ) ∼ 1 with ℜ(λ) > 0. The strongly unstable spectrum is absent for reservoir computing applications, since, otherwise, the reservoir's state is strongly unstable, and the echo state (or fading memory) property [1] is lost to a large extend. Formally, the condition for the absence of the strong unstable spectrum is the stability (all eigenvalues have negative real parts) of the linearization matrix A of instantaneous terms (see equation 11). That means that either the operating point of the reservoir shifts to another attractor or the reservoir is unstable. Hence, we focus on the pseudocontinuous spectrum, which can be obtained by introducing the Ansatz: where γ and µ are two new real variables. Substituting equation (12) into (11) one gets in the leading order: Equation (13) is a polynomial with respect to e −γ−iµτ . If Y j (µ) are solutions of this polynomial, then: and γ j (µ) = − ln |Y j (µ)| are the rescaled real parts of the eigenvalues from the pseudocontinuous spectrum. More exactly, the curves γ j (µ)/τ + iµ are approximated with the eigenvalues for large τ .
In the case of the Lang-Kobayashi system, there are two solutions for the real parts of the pseudocontinuous eigenvalue spectrum, see the derivation in [53], where ϵ = T −1 LK , and A 2 l = P−N * 2N * +1 is the constant intensity, with the corresponding inversion N * at the external cavity mode (ECM). ECMs are the solutions of the Lang-Kobayashi system of the form E = A l e iωt , N = N * with constant A l and N * , which play the role of the equilibria. Due to the S 1 symmetry of the system, each of these solutions can be transformed into an equilibrium with the corresponding characteristic equation (11). The Lang-Kobayahi system possesses many ECMs, however, for the case α = 0, we consider the ECM with N * = −κ, which is the most stable [53]. In general our method of linearization applies to all fixpoint solutions of any system, thus α ̸ = 0 can be analysed as long as the system is in fixpoint. For α ̸ = 0 the system often starts to jump between different solutions, where some are not of fixpoint nature. Thus comparisons between the different operating points is not as simple and the focus on the essential new method does get lost, why we choose α = 0.
In order to approximate the imaginary parts µ of the pseudocontinuous spectrum, we consider the argument of (14) and obtain: where µ j,k is the imaginary part of a kth eigenvalue on the jth branch. For the purpose of this paper, we need an approximation of the eigenvalues around the origin. As one can simply show using (17), for large τ , these eigenvalues (their imaginary parts) can be approximated as: as soon as k/τ ≪ 1, see also [63,64] for more detailed estimations. In the case of the Lang-Kobayashi system, the roots Y j (0) are real, hence, we have either argY j (0) = 0 or π. This leads to: with ν = 0 for argY j (0) = 0, and ν = 1 for argY j (0) = π. Hence, all imaginary parts µ j,k are integer multiples of π/τ . In particular, for any T ≈ jτ , j ∈ N (T ≈ jτ /2 for ν = 0, respectively), the product µ j,k T is proportional to an integer number of π. This kind of resonance occurs for all considered eigenvalues (independent of k), and it plays an important role in the total memory loss of the reservoir, which is discussed in section 3.2 below. The second method for computing the eigenvalue spectrum is based on DDE-biftool [60,61], which is a path-continuation package for Matlab capable of numerically computing the eigenvalues. In our case, we compute the first 100 eigenvalues with the highest real parts to compare these with results from the memory capacities.
We also consider the case of no feedback κ = 0. This yields a solitary semiconductor laser system that can be tuned from being an effectively 1-dimensional problem (Class A-like for T LK ≪ 1) to a 2D problem (Class B-like for T LK ≫ 1) [65]. We will use a linearization and numerical evaluation of the eigenvalue problem. Even though the laser system is 3-dimensional, it possesses the S 1 symmetry E → Ee iϕ allowing to reduce the dimension by one.
We would like to emphasize that the eigenvalue method would also apply to a reservoir computer with a different form of information input, e.g. optical injection in the case of the Lang-Kobayashi system, since the analysis is performed without any reservoir computer input. As long as the input is a small perturbation to the system, the responses are fully described by the linearized system.

Simulation description
Simulations have been performed in C++ with standard libraries, except for linear algebra calculations, which were done via the linear algebra library 'Armadillo' [66]. A Runge-Kutta 4th order method was applied to integrate the delay-differential equation given by equations (9) and (10) numerically, with an integration step ∆t = 0.01 in time units of the photon lifetime. The noise strength is D noise = 10 −7 in all simulations. After simulating the system without reservoir inputs to let transients decay, a buffer time of 100 000 inputs was applied (this is excluded from the training process). In the training process, 250 000 inputs were used to have sufficient statistics (N = 250 000). Afterwards, the memory capacities were calculated, whereby a testing phase is not necessary, because the independent and identically drawn uniformly distributed inputs u are statistically equal if drawn for training or testing phases. All possible combinations of the Legendre polynomials up to degree D = 5 and 500 input steps into the past were considered (i = −500). Capacities C d y {u} below 0.001 were excluded because of finite statistics. For calculating the matrix inverse, the Moore-Penrose pseudoinverse from the C++ linear algebra library 'Armadillo' was used. In the case of the NARMA10 task, 25 000 inputs for training and testing were used. For the piecewise-constant T-periodic masking function g independent and identically distributed random numbers between [0, 1] were used.
For all simulations, the input strength η was fixed to 0.01. The small input strength was used to guarantee the linear answers of the reservoir and, hence, the relevance of the eigenvalue analysis.

Geometrical intuition
In this paper we will use two quantities Φ and Λ to approximate the memory capacity properties of the reservoir computer. For these two values we would like to give a geometrical intuition, shown figure 2.
The first value Φ = ℑ(λ)T we call the relative angular distance between two inputs, where ℑ denotes the imaginary part of the eigenvalue λ. Here λ is a critical eigenvalue, i.e. one having its real part close to 0. Φ geometrically describes the angular distance between two distance vectors δs 1 and δs 2 of the system's state s 1 and s 2 at two instances in time separated by one clock cycle interval T. If this relative angular distance is a multiple of π the responses tend to overlap, reducing the separatibility of the inputs, thus degrading the reservoir computer performance.
The second quantity Λ = e ℜ(λ)T = |δs 2 |/|δs 1 | describes the distance reduction between two perturbed states, where ℜ denotes the real part of the eigenvalue λ. Λ describes the contraction of the system's state towards a new fixpoint due to a new reservoir input. To distinguish two responses s 1 and s 2 for two different inputs u 1 and u 2 , the distance |δs 1,2 | (see figure 2) between two responses should be large enough. On the other hand, if the reaction of the system is very fast, i.e. very negative eigenvalues, the system has a high echo state property and thus low memory capacities for any inputs longer than a few (in the worst case even longer than one) steps into the past. If the remaining information of the input nth steps back is degrading very fast (very negative eigenvalues), the systems capability to recall is lowered, and at some point reaches the level of the system noise. The distance reduction Λ = e ℜ(λ)T gives a good estimation for both of these properties.
In this paper we will show that both quantities together pinpoint to well performing reservoir computer setups.

Results
The following section is structured as follows. First, we will discuss a Lang-Kobayashi laser with κ = 0, i.e. a solitary laser system, as a reservoir to simplify and depict general concepts. Afterward, we will activate the delay and look at the full Lang-Kobayashi system as a reservoir computer.

Laser without feedback
We first consider a solitary semiconductor laser system as a reservoir. One has to think of the virtual nodes not to be located on the delay line, but rather as time separated readouts of the system state, that are used in a linear combination. We set κ = 0 in equations (9), (10) and use 10 virtual nodes (N V = 10). For the considered parameter values P = 0.05 and α = 0 and without input and noise (η = 0, D noise = 0), the system's solution converges to a single stable ECM, for which we compute the two eigenvalues. The two eigenvalues are plotted in figure 3(b) as a function of T LK , which gives from left to right the transition from Class A to Class B laser.
To compare the two eigenvalues with the recall capability of the laser, we plot the computed linear, quadratic, cubic, and total memory capacities in figure 3(a). The memory capacities do not change significantly for T LK ≲ 3 where the system corresponds to a class A laser with an adiabatic approximation of the charge carriers. For these parameter values, as one can see from the real parts of the eigenvalues, one eigendirection is considerably faster than the other, and thus can be ignored. At T LK ≈ 2, the transition from a Class A laser to a class B laser appears, whose steady state solution is a focus. The additional degree of freedom of the charge carrier dynamics leads to an increase of the total memory capacity from about 5 to the theoretical maximum of 10. Figure 3(c) shows the angular distance Φ = ℑ(λ)T taken modulo π, which is based on the rotation ℑ(λ)T of a small perturbation vector in the 2D phase space during the evolution over the time-interval T (see figure 2). The discontinuities of Φ in figure 3(c) (indicated with vertical dashed purple lines in figure 3(a)) correspond to resonances, i.e. integer numbers of half-a-circle rotations. Comparing the memory capacity at these points in the class B regime, one observes dips in the linear memory capacity and slight changes in the higher-order memories. This effect is pronounced if at the same time real parts of the eigenvalues are close to 0. Since the degradation of the linear memory coincides with the discontinuities in Φ mod π, it can be linked to an overlapping of the systems responses, and to a decreasing linear separatibility of the output. We would like to emphasize that even though the system has no optical feedback (κ = 0) the dynamical system still can act as a reservoir with very short-memory. This comes from the fact that the reaction of the system is not instantaneous yielding a memory kernel of a few inputs into the past. The memory of the reservoir is limited by the real part of the largest critical eigenvalue. This is due to the fact, that the real part of the largest eigenvalue yields the timescale on which small perturbations to the fixpoint exponentially decay.
For a larger picture of the resonance effects at Φ ≈ kπ (k ∈ N), a 2D parameter scan was done as a function of the timescale ratio T LK and the clock cycle T (shown in figure 4). The linear, quadratic, cubic, and total memory capacities are color-coded in panel (a)-(d). Bright regions in (a)-(d) correspond to high memory capacities, while dark blue to low memory capacities. The black dashed line shows the scan from figure 3. Purple solid lines show the parameter values where Φ ≈ kπ. The influence of the angular distance Φ is most prominent in figure 4(b), where dips are visible in the linear memory capacity. Its influence on the higher-order memory capacity is also detectable, but harder to describe, as both quadratic and cubic memories either decrease or increase depending on the resonance line.
The solid red lines denote parameter values, where the distance reduction Λ = e ℜ(λ)T = D noise = 10 −7 for the two eigenvalues λ 1,2 of the solitary laser system. The two red arrows indicate the direction in the parameter space for decreasing ℜ(λ) and thus decreasing Λ. The distance reduction shows a decrease in the memory capacities for a decrease in Λ. This rises from the fact that lower Λ correspond to faster eigendirections and thus faster echo state properties.
Combining the information about the two quantities Φ and Λ and comparing it with the memory capacity, we can pinpoint to well performing reservoir computers for the class B and class A laser system. Namely, the linear memory capacity has larger values in the absence of resonances Φ ≈ kπ and for values of Λ closer to 1. We now want to expand this knowledge to the case of a laser with external feedback.

Laser with feedback
We now expand our results to the infinite-dimensional phase space of a semiconductor system with delay, i.e. the Lang-Kobayashi system. In [9,22,26,45,46] it was shown that resonances between τ and T often decrease memory capacity and thus reservoir computing performance. Here we look at this phenomenon from another point of view, namely, as a resonance between T and the imaginary parts of the eigenvalues. We use the resonance property described in section 2.5: for certain resonant values of T, the product ℑ(λ)T is proportional to an integer number of π for all critical eigenvalues simultaneously. We computed the first 100 eigenvalues using DDE-biftool for the Lang-Kobayashi system. By superimposing all Φ i , where i ∈ 0, 1, 2..., N  is the index of the ith eigenvalue, we evaluate the resonance effects of the strongest N eigendirections by computing the average angular distance:Φ and compare the results with the linear, quadratic, and total memory capacities. The comparison of the memory capacities andΦ is shown in figure 5, where a 2D parameter scan is plotted in the parameter plane of the delay time τ and the clock cycle T. Bright regions in (a-c) correspond to high memory capacities, A match with lower total memory capacities, especially for the linear memory is clear. For a reservoir to be applicable to many tasks, a higher total memory capacity is desirable, thus a decrease in the total memory capacity decreases the applicability of the reservoir to many different tasks. Our results support the fact that the clock cycle T should be chosen to be off resonant of the delay time τ . The eigenvalue analysis gives an additional explanation and intuition for why this is the case. Taking into account the resonance effect and our results from [45], we set the delay time τ ≈ √ 2T for all following simulations. As we have seen in section 3.1 for a 2D reservoir, the reservoir performance decreases when the real part of the eigenvalues becomes strongly negative. In such a case, the reservoir 'forgets' the input too fast. Here we extend this idea to the case of the infinite-dimensional reservoir.
As long as the perturbation from the information fed into the system is small enough, one can think of all eigenvalues and their corresponding eigendirections as the available phase space of the reservoir computer. Thus, a higher phase space volume can lead to a more promising reservoir computer. We introduce the average of the distance reduction Λ by:Λ It describes the average distance reduction of the N slowest eigendirections. Since only a finite number of complex eigenvalues lie to the right of a line parallel to the imaginary axis, all eigendirections except a finite number are strongly contracting, i.e. possess strongly negative real parts [67]. This implies the possibility of considering only a finite number N of eigenvalues in equation (21). Figure 6 depicts the memory capacities and the distance reductionΛ as a function of the timescale ratio T LK , or in other words, the evolution of the memory capacities along the transition from a Class A to a Class B laser system with delayed feedback. Similarly to the case without feedback in figure 3, the memory capacity stays about constant for T LK ≲ 2, and increases when the additional dimensions become available by the reservoir for T LK ≳ 2. The increase ofΛ coincides with the increase of the linear memory. The higher orders show a similar trend, but are, in general, more involved and should be investigated more deeply. Thus, the knowledge of the eigenvalues provides a qualitative prediction of the linear memory capacity.
To give a broader overview, we perform a 2D parameter scan along the feedback strength κ and pump P (shown in figure 7) and plot the linear, quadratic, and total memories as a colour-code. Bright regions correspond to high memory capacities, and dark to low memory capacities. Additionally, in figure 7(d), the sum of the average distance reductionΛ is color-coded within the same 2D parameter plane spanned by κ and P. Comparing the three memory capacity scans withΛ, we can see a close relationship between them. Thus,Λ is a very good indicator for choosing well-performing reservoir computers. This saves a lot of computational efforts, as the eigenvalues can be computed in a fraction of the time needed to compute the memory capacities.   figure 7, which correspond to parameters close to and well above threshold. The first parameter set (figure 8(i)) corresponds to an eigenvalue spectrum for an operation point close above the threshold with a low power output. Here the laser system possesses more eigenvalues with real parts close to 0, thus it has many slowly contracting eigendirections which meansΛ is closer to 1. CalculatingΛ for the first parameter set yieldsΛ = 0.75. The second parameter set (figure 8(ii)) corresponds to a laser that is operated high above threshold. This laser has less slowly contracting eigendirections, i.e.Λ is closer to 0. CalculatingΛ for the second parameter set yieldsΛ = 0.6. Now we use the insights gained from the distance reduction Λ and from the angular distance Φ and test the reservoir computer performances for the two parameter sets from figures 8(i) and (ii), marked as black crosses in figure 7. The performance is quantified by evaluating both the memory capacity and the prediction error (NRMSE) for the NARMA10 task shown in figures 9(b) and (a). On the horizontal axis, we change the number of virtual nodes N V , i.e. we increase the number of readout dimensions, which should, naively thinking, increase the performance of the reservoir. While we do this, we keep the distance between the virtual nodes θ the same for 5 different cases of θ from θ = 1 up to θ = 5, shown in black and red with decreasing brightness for the optimized (figure 8(i)) and not-optimized (figure 8(ii)) point respectively. Increasing the virtual node distance θ should reduce the linear dependency of the nodes, as the time between two responses is increased. This is obviously dependent on the reaction time of the system, which is also given by the eigenvalues of the system. Thus the influence of increasing θ on the slowly contracting eigendirections (figure 8(i)) is pronounced compared to the one with fast eigendirection ( figure 8(ii)). The increase of θ also effectively increases the clock cycle T = N V θ and thus the delay time τ = 1.41T. We want to emphasize, that this does not alter the general trend of P = −0.095 (figure 8(i)) having many slowly contracting eigendirections compared to P = 0.095 (figure 8(ii)).
The results indicate that even though the number of virtual nodes increases, the NARMA10 error for the case with the low distance reductionΛ (ii) does not go below 0.45. On the other hand the case with the higĥ Λ, our optimal case (i), reaches very small errors below 0.15, a factor 3 better than the low distance reduction Λ (ii) case. We also want to emphasize that the simulation was done for a high noise value of D noise = 10 −7 .
Simulating the system without any noise D noise = 0 NARMA10 errors (NRMSE) of down to 0.05 were reached. We conclude that a high distance reduction is very beneficial for the performance.
As a dashed black line in figure 9(a) we additionally show the minimum reached by a linear regression without a reservoir. Every reservoir setup with results below this line has higher performance and thus can be considered an improvement on the NARMA10 task. We included this here to emphasize the reduction of the NARMA10 error (NRMSE) by the eigenvalue analysis, for which an improvement of about a factor 4 is reached to the linear regression without reservoir.
The total memory capacities for the two cases (i) and (ii) are shown in figure 9(b). We can see the same trend: the memory capacity reaches a limit for the case with lowΛ, whereas the improved case with the highestΛ increases in its memory capacity further for higher N V . The results suggest that operation points with high distance reductionΛ (solid lines in figure 9) pinpoint to well-performing reservoir computers.

Conclusion
We have shown that the eigenvalue spectrum analysis of a dynamical system used as a reservoir (e.g. a laser described by the Lang-Kobayashi system) is capable of predicting good reservoir computing operation points. Because of the available analytical and numerical tools for the description of the eigenvalue spectrum, such analysis can be readily applied for different dynamical systems, which are used as reservoirs with operating points close to an equilibrium. The eigenvalue method is of magnitudes faster to compute and could help in numerically predicting good reservoir computers for experimental setups.
Due to the relation between the eigenvalue spectrum and the performance of the delay-based reservoir computing, the central message of this paper is twofold: first, the eigenvalues must be off-resonant, where the resonance condition is given in terms of the imaginary parts of the eigenvalues. Namely, the product of the latter with the input clock-cycle should be away from values of kπ. Importantly, such resonances appear for all critical eigenvalues at almost the same parameter values, due to general properties of the spectrum of delay systems with large delays [59]. Therefore, such an off-resonant condition plays an important role even when the reservoir's effective dimensionality is high.
The second conclusion is that, for optimal performance, the spectrum must be close to criticality. This closeness is measured by the real part of the eigenvalue spectrum, which should be close to zero and negative. In this paper, we propose the average distance reduction as a measureΛ for such closeness, which is given by equation (21).
For future work we would like to analyse clearer properties of the reservoir via the eigenvalue spectrum. We believe that the timescale analysis, i.e. eigenvalue analysis, has more advantages. The eigenvalues do not only yield the timescale on which the system forgets, but also the timescale on which the system learns (for accordingly small inputs). Thus, through the eigenvalues we can define a system which is either fast and has short term memory or we can define a system like we did in this paper, which is slow and has more memory capacity further into the past.

Appendix A. Memory capacity expression (3)
We show how the expression for the memory capacity from [40] can be rewritten in the form of equation (3). From Dambre et al [40], the capacity to approximate a target dataŷ is given by: Here s i is the ith readout of the M system responses for the lth input-target pair, ⟨ν⟩ L = 1 L L l=1 ν l is the average over all input-output paris L and ⟨s i s j ⟩ −1 L is the inverse of ⟨s i s j ⟩ L . One can insert the average over all input-output pairs yielding: The first termŷ l s i,l is the ith system response to the lth input-output pair (ith column and lth row) multiplied with the lth target. Summing over all input-output pairs L, this is the same as the ith entry of the matrix product: