cs.DCFeb 23, 2026
SaveFast Linear Reservoirs via Diagonalization
Organizations: Mnemosyne, DI-ENS · ENS PSL, Inria Center of Bordeaux University, LaBRI, IMN · Mnemosyne · Inria Center of Bordeaux University, LaBRI, IMN ·
Abstract
We introduce a diagonalization-based optimization for Linear Echo State Networks (ESNs) that reduces the per-step computational complexity of reservoir state updates from quadratic to linear. By reformulating reservoir dynamics in the eigenbasis of the recurrent matrix, the recurrent update becomes a set of independent element-wise operations, eliminating the matrix multiplication. We further propose three methods to use our optimization depending on the situation: (i) Eigenbasis Weight Transformation (EWT), which preserves the dynamics of standard and trained Linear ESNs, (ii) End-to-End Eigenbasis Training (EET), which directly optimizes readout weights in the transformed space and (iii) Direct Parameter Generation (DPG), that bypasses matrix diagonalization by directly sampling eigenvalues and eigenvectors, achieving comparable performance to standard Linear ESNs. Across all experiments, both our methods preserve predictive accuracy while offering significant computational speedups, making them a replacement for standard Linear ESNs computations and training, and suggesting a shift of paradigm in linear ESN towards the direct selection of eigenvalues.
Figures & tables
Figure 1: Standard architecture of an Echo State Network (ESN). In the Reservoir Computing paradigm, the input weights , the internal recurrent weights , and the optional feedback weights are randomly initialized and kept fixed. The input is projected into a high-dimensional space within the reservoir, which utilizes its recurrent connections to maintain a deterministic state representation over time. Only the readout weights are trained to decode this internal state and produce the final output.
Figure 2: Comparison of standard computation and proposed optimizations for different reservoir sizes. (i) Generation Step: Evaluates three initialization methods. Normal generates a standard linear reservoir with a weight matrix . Diagonalization generates a standard and then diagonalizes it (applicable to EWT/EET in section 4 ). DPG directly generates a diagonal matrix (see section 4 ). (ii) Reservoir Step: Compare the standard computation with the Diagonal one. Because EWT, EET, and DPG share an identical diagonal structure after the generation step, they are represented together as Diagonal . (iii) Readout Step: Displays only a single curve because the computational cost is identical across all methods. With the implementation proposed in Appendix A , the training of the readout can be performed with real matrices, equating the readout cost of the standard method. (Note) For reservoir and readout steps, the duration displayed is for a single time step. To evaluate the total duration over the whole sequence, these values must be multiplied by the total number of time steps.
Figure 3: Comparison of eigenvalue distributions in the complex plane. On the first column (left) the spectrum derived from a standard random reservoir matrix . On the second column the spectrum generated via Uniform Distribution (Algorithm 1 ). On the third column the spectrum generated via the deterministic Golden Distribution method (Algorithm 3 ) without noise, and on the fourth column the specturm generated via the Noisy Golden Distribution . We observe that the Noisy Golden distribution (right) achieves a significantly more homogeneous coverage of the unit disk than the Uniform Distribution or its non-noisy counterpart, effectively matching the spectral density of the standard reservoir (left).
Figure 4: Illustration of the Multiple Superimposed Oscillators (MSO) time series for the task (MSO5). The target signal is generated by summing five distinct sinusoidal components. The complete sequence of 1000 time steps is partitioned into a training part of 400 steps (blue) which includes an initial 100-step washout used to discard transient reservoir dynamics (dashed blue), a validation of 300 steps (green) for hyperparameter tuning, and a final test of 300 steps (pink) to evaluate predictive performance.
| Hyper-parameter | Values considered |
|---|---|
| Reservoir size ( ) | 100 |
| Input scaling ( ) | 0.01, 0.1, 1 |
| Leaking rate ( lr ) | 0.1, 0.3, 0.5, 0.7, 0.9, 1.0 |
| Spectral radius ( ) | 0.1, 0.3, 0.5, 0.7, 0.9, 1.0 |
| Ridge regression regularization ( ) |
Table 1: Hyper-parameters values considered for model selection on the MSO tasks reported from Gallicchio et al. (2017)
| Task | Normal | Diagonalized | Uniform Dist. | Golden Dist. | Noisy Golden Dist. | Sim Dist. |
|---|---|---|---|---|---|---|
| ( ) | ( ) | |||||
| MSO1 | 1.65e-14 | 1.58e-14 | 5.85e-14 | 2.49e-14 | 4.77e-14 | 3.56e-14 |
| MSO2 | 2.55e-13 | 2.78e-13 | 2.28e-13 | 1.45e-13 | 2.39e-13 | 2.44e-13 |
| MSO3 | 5.42e-12 | 9.14e-12 | 4.49e-12 | 9.07e-12 | 6.14e-12 | 8.37e-12 |
| MSO4 | 1.39e-10 | 5.77e-10 | 3.64e-10 | 7.22e-11 | 6.93e-11 | 2.28e-10 |
| MSO5 | 2.75e-09 | 4.03e-08 | 2.95e-08 | 5.24e-10 | 1.63e-09 | 1.87e-08 |
Table 2: Performance comparison on the Multiple Superimposed Oscillators (MSO) across 10 independent seeds. Metric is the Root Mean Square Error (RMSE). Normal corresponds to a standard , Diagonalized to the EET method, others corresponds to the DPG method with different initialization strategies.
Figure 5: Visualization of spectral importance via readout weights on the MSO task. The reservoir eigenvalues are plotted in the complex plane, where the marker size is proportional to the absolute value of the corresponding weight in (normalized between 0 and 1). Larger points identify the specific eigenvalues that contribute most significantly to the model’s prediction for this task. Points are the same on the left and on the right, just the size of points is changing: on the right, points close to the circle are non longer visible which may create an illusion that the point distribution was contracted, but that’s not the case.
Figure 6: Memory capacity as a function of delay across different reservoir sizes. The four subplots illustrate the memory capacity evaluated against the delay for reservoir sizes N=100, 300, 600, and 1000. As expected, increasing the reservoir size brings an overall improvement in memory capacity. The performance of five different methods is compared: the standard Echo State Network with a classic weight matrix ( Normal ), Direct Parameter Generation (DPG) utilizing a Uniform Distribution , DPG utilizing the deterministic Golden Distribution , and DPG utilizing the Sim Distibution which uses eigenvalues extracted from a standard random matrix combined with randomly generated eigenvectors . Observing the results, the Diagonal Uniform method exhibits a more balanced degradation profile than the Normal baseline, though they remain roughly equivalent and intersect near a memory capacity of 0.5. Notably, the Golden Distribution initialization systematically outperforms the standard Normal method across all reservoir sizes. Conversely, the Sim Distribution curve closely follows the Normal baseline, suggesting that the eigenvectors play a secondary role compared to eigenvalues, even if it demonstrates a small but consistent underperformance.
Figure 7: Memory capacity as a function of reservoir connectivity. The four subplots illustrate the memory capacity evaluated across varying levels of matrix connectivity for reservoir sizes . To ensure a task of moderate and proportional difficulty, the requested delay for each configuration is chosen so that the memory capacity yields 0.5 when the connectivity is one (to chose the corresponding value we looked at Figure 6 ). The graphs compare the performance of the standard weight matrix approach ( Normal ) against its diagonalized counterpart ( Diagonalization , corresponding to the EET and EWT methods), with the secondary axis tracking their absolute performance gap ( Difference ). As expected, both methods experience a severe drop in memory capacity at extremely low connectivity levels; as the matrix becomes very sparse, the network loses the dimensionality required to maintain temporal representations. Moreover, the results highlight a specific minimum connectivity threshold which depends of the reservoir size. Below this connectivity, the Diagonalization method begins to systematically underperform the Normal baseline. At such extreme sparsity, the standard matrix lacks sufficient structural complexity for an effective eigendecomposition, causing the eigenspectrum to collapse and lose essential dimensional richness. Above this minimal threshold, the diagonalized approach successfully matches the standard baseline’s performance.
More in this topic
Reservoir Computing
Reservoir computing uses recurrent dynamical systems to transform temporal inputs into features that support efficient learning with simple readouts.
See all 107 papers