# Linear Reservoir: A Diagonalization-Based Optimization

Romain de Coudenhove, Yannis Bendi-Ouis, Anthony Strock, Xavier Hinaut

► To cite this version:

Romain de Coudenhove, Yannis Bendi-Ouis, Anthony Strock, Xavier Hinaut. Linear Reservoir: A Diagonalization-Based Optimization. 2026. hal-05521416

HAL Id: hal-05521416

https://hal.science/hal-05521416v1

Preprint submitted on 20 Feb 2026

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# Linear Reservoir: A Diagonalization-Based Optimization

**Romain de Coudenhove**

ENS PSL, Inria Center of Bordeaux University, LaBRI, IMN

ROMAIN.DE.COUDENHOVE@ENS.PSL.EU

**Yannis Bendi-Ouis**

Inria Center of Bordeaux University, LaBRI, IMN

YANNIS.BENDI-OUIS@INRIA.FR

**Anthony Strock**

Department of Psychiatry & Behavioral Sciences, Stanford University School of Medicine Stanford

ASTROCK@STANFORD.EDU

**Xavier Hinaut**

Inria Center of Bordeaux University, LaBRI, IMN

XAVIER.HINAUT@INRIA.FR

**Editor:** Preprint JMLR

## 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 $O(N^2)$ to $O(N)$. 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 than standard Linear ESNs. Across all experiments, both our methods preserve predictive accuracy while offering significant computational speedups, making them a replacement of standard Linear ESNs computations and training, and suggesting a shift of paradigm in linear ESN towards the direct selection of eigenvalues.

**Keywords:** Reservoir Computing, Linear Echo State Networks, Diagonalization, Computational Complexity, Recurrent Neural Networks, Parallelization, Efficient Training, Eigenvalue Decomposition, Spectral Initialization, Linear Recurrent Units

## 1 Introduction

Reservoir Computing (RC), most notably implemented as Echo State Networks (ESN), offers a compelling alternative to traditional Recurrent Neural Networks (RNNs) by keeping internal reservoir weights fixed and training only the output layer. While this approach simplifies the training procedure to a linear regression problem, the reservoir state update remains computationally expensive, scaling with quadratic complexity ($O(N^2)$) where N is the number of neurons). Consequently, there is a strong motivation to reduce this complexity without sacrificing performance. While early RC research emphasized the necessity of non-linear reservoir dynamics, recent studies have demonstrated the surprising efficacy of linear reservoirs. Theoretically, it has been shown that a linear reservoir coupled with a non-linear readout (such as a polynomial function or neural network) is a universal approximator
(Gonon and Ortega, 2019). Furthermore, linear dynamics are particularly well-suited for maximizing memory capacity, as analyzed in Dambre et al. (2012). Several architectures have leveraged these properties. For instance, Legendre Memory Units (LMUs) utilize linear dynamics derived from Legendre polynomials to project input history onto an orthogonal basis, thereby extending memory capacity (Voelker et al., 2019). Similarly, recent works by Peter Tino (Simon Fong et al., 2025; Verzelli et al., 2021; Tino, 2020) have explored the dynamical properties of linear models, clarifying the precise conditions under which linear ESNs can match or outperform their non-linear counterparts. Despite these advantages, the computational cost of matrix multiplication remains a bottleneck for large reservoirs. To address this, approaches such as Next Generation Reservoir Computing (Gauthier et al., 2021) have proposed replacing the explicit simulation of reservoir dynamics with static non-linear feature expansions, offering significant speedups. However, abandoning the recurrent dynamics entirely may not always be desirable as it limits the model to a fixed size window.

In this work, we present a method to reduce computational complexity while maintaining the explicit recurrence of the reservoir. We focus on optimizing Linear ESNs through matrix diagonalization, a transformation that converts costly matrix multiplications into efficient element-wise operations, reducing the complexity from $O(N^2)$ to $O(N)$. Beyond immediate computational efficiency, our approach is motivated by a fundamental advantage of linear systems: the ability to parallelize computation across the entire input sequence (see Appendix B). This principle of parallelizable state propagation aligns with recent advances in modern state-space models such as Mamba (Gu and Dao, 2024) and parallelized LMU implementations (Chilkuri and Eliasmith, 2021). First, Section 2 formalizes the Linear ESN framework and highlights the computational cost of standard dynamics. Section 3 details our diagonalization-based reformulation which reduces the update complexity to $O(N)$. Section 4 presents three practical implementations of this optimization, including the direct generation of spectral parameters. Finally, Section 5 validates these methods on standard benchmarks, followed by a discussion on the importance of spectral properties in Section 6.

## 2 Linear Echo State Networks: Definition

An Echo State Network (ESN) is a type of recurrent neural network particularly effective in learning complex dynamical systems. In its linear form, the network is fully characterized by a set of matrices:

* $W \in \mathbb{R}^{N \times N}$, the reservoir transition matrix

* $W_{\text{in}} \in \mathbb{R}^{D_{\text{in}} \times N}$, the input matrix

* $W_{\text{fb}} \in \mathbb{R}^{D_{\text{out}} \times N}$, the feedback matrix (optional)

* $W_{\text{out}} \in \mathbb{R}^{N' \times D_{\text{out}}}$, the readout matrix

Here, $N$ denotes the number of neurons in the reservoir, while $D_{\text{in}}$ and $D_{\text{out}}$ represent the dimensionality of the input and output signals, respectively. The dimension $N'$ corresponds to the size of the extended state vector used by the readout. Depending on the chosen configuration, this vector is formed by concatenating the reservoir state $r(t)$ with optional
<figure><img src="image_2.png" /><figcaption>Figure 1: Standard architecture of an Echo State Network (ESN). In the Reservoir Computing paradigm, the input weights $W_{in}$, the internal recurrent weights $W$, and the optional feedback weights $W_{fb}$ 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 $W_{out}$ are trained to decode this internal state and produce the final output.</figcaption></figure>

components: a constant bias, or the previous output $y(t-1)$. Consequently, $N'$ varies based on the inclusion of these additional terms alongside the $N$ reservoir units.

This structure leverages a fixed, randomly initialized reservoir $W$, where only the output weights $W_{out}$ are optimized. Let the time series input be $(u(t))_{t=1}^T \in \mathbb{R}^{T \times D_{in}}$, where $u(t) \in \mathbb{R}^{D_{in}}$ represents the input vector at time step $t$. For every time step $t$, the reservoir produces an internal state $r(t) \in \mathbb{R}^N$, and the readout produces an output $y(t) \in \mathbb{R}^{D_{out}}$.

## 2.1 Reservoir step

The internal state $r(t)$ evolves according to the following dynamics. This recursive update propagates information through time via a weighted combination of past states, current inputs, and past outputs.

$$r(t) = r(t-1)W + u(t)W_{in} + y(t-1)W_{fb} \quad (1)$$

**Note on optional input bias:** An input bias term can be added to the reservoir. However, since its effect is mathematically equivalent to adding a constant dimension of 1 to the input signal $u(t)$, we omit the input bias for the remainder of this work.

## 2.2 Readout step

The output $y(t)$ at time $t$ is a linear combination of an extended state vector (which includes a bias, the previous output, and the current reservoir state) and the readout weight matrix $W_{out}$.

$$y(t) = [1 \ y(t-1) \ r(t)] W_{out}, \quad \text{for } t = 1, \dots, T \quad (2)$$
$$W_{\text{out}} = \begin{bmatrix} W_{\text{out,bias}} \\ W_{\text{out,out}} \\ W_{\text{out,res}} \end{bmatrix} \qquad (3)$$

**Note on optional features:** Depending on the implementation, both bias and previous output are optional, only the current reservoir state is mandatory. Consequently $W_{\text{out,bias}}$ and $W_{\text{out,out}}$ are optional too.

**Note on initial states:** According to these equations, the model requires an initial state for the reservoir and for the output if there is a feedback. By specifying these initial conditions, the entire reservoir and output trajectories are uniquely determined for a given input sequence. Numerous methods exist to initialize those vectors, but in most of the case they are initialized either by setting them to the zero vector, or by using a preliminary « warmup » phase where the network runs for a short period to stabilize the state dynamics before training starts. In this context, we set the initial conditions to zero.

## 2.3 Leaking Rate

The leaking rate $lr \in [0, 1]$ is a hyperparameter in reservoir computing that controls the temporal decay of internal states. It introduces a form of memory retention by blending the previous state and the potential new state at each time step. A smaller value of $lr$ results in a slower update and consequently a bigger memory capacity, while $lr = 1$ yields the fastest update and consequently the smallest memory capacity. In practice in Linear ESNs, the leaking rate modifies the transition matrix $W$ by mixing it with the identity:

$$W^{(lr)} = lr \cdot W + (1 - lr) \cdot I \qquad (4)$$

The input and feedback weights are similarly rescaled:

$$W_{\text{in}}^{(lr)} = lr \cdot W_{\text{in}}$$

$$W_{\text{fb}}^{(lr)} = lr \cdot W_{\text{fb}}$$

As a result, the full reservoir dynamics under leaking rate become:

$$r(0) = 0, \qquad (5)$$

$$r(t) = r(t - 1) W^{(lr)} + u(t) W_{\text{in}}^{(lr)} + y(t) W_{\text{fb}}^{(lr)} \qquad (6)$$

Leaking rate only reparametrize the original matrices $W$, $W_{\text{in}}$ and $W_{\text{fb}}$, which does not alter the structure of the network. Thus, the exact same optimization method (section 3) can be used with $W^{(lr)}$, $W_{\text{in}}^{(lr)}$ and $W_{\text{fb}}^{(lr)}$.

## 2.4 Learning $W_{\text{out}}$

Given input-output pairs $\{(u(t), y(t))\}_{t=1,...,T}$, where $T$ is the total number of timesteps used for training, we construct:

$$X(t) = [1 \ y(t-1) \ r(t)] \in \mathbb{R}^{N'} \qquad (7)$$

$$Y(t) = [y(t)] \in \mathbb{R}^{D_{\text{out}}} \qquad (8)$$
$W_{out}$ is computed using standard ridge regression with $\alpha$ the regularization parameter:

$$W_{out} = (X^T X + \alpha I)^{-1} X^T Y \qquad (9)$$

## 2.5 Computational Complexity

Here, we analyze each major component of the ESN architecture in terms of time complexity, explicitly accounting for sparsity (1 − connectivity) where relevant.

* **$W_{in}$ and $W_{fb}$ generation:** One of the input weight matrix $W_{in} \in \mathbb{R}^{D_{in} \times N}$ or the feedback weight matrix is optional. These matrices are typically generated with a connectivity $c_{in}/c_{fb}$ (where each weight represents a non-zero connection with probability $c_{in}/c_{fb}$). The non-zero weights are sampled from a specific distribution (e.g., normal or uniform). The complexity of generating these matrices are $O(D_{in}N)$ and $O(D_{out}N)$.

* **$W$ generation and Spectral Radius scaling:** The reservoir weight matrix $W \in \mathbb{R}^{N \times N}$ is similarly generated with a connectivity $c_r$. A critical step is scaling $W$ to achieve the desired spectral radius $\rho_0$. The computational cost depends heavily on the algorithm used to compute the leading eigenvalues: For dense matrices, standard eigenvalue decomposition algorithms are $O(N^3)$. For sparse matrices, iterative methods such as the Implicitly Restarted Arnoldi Method (IRAM) are commonly use (Lehoucq et al., 1998). The complexity depends on the sparsity and convergence rate, scaling between $O(N^2)$ and $O(N^3)$.

* **Reservoir step:** This recursive update involves matrix-vector multiplications. By exploiting sparse matrix operations, the complexity per time step is reduced to the number of non-zero elements: $O(c_r N^2 + c_{in} D_{in} N + c_{fb} D_{out} N)$ per time step. Therefore, the total cost over $T$ time steps is $O(T(c_r N^2 + c_{in} D_{in} N + c_{fb} D_{out} N))$.

* **Readout Layer:** The output computation $y(t) = X(t)W_{out}$ involves a dense matrix multiplication. With $W_{out} \in \mathbb{R}^{N' \times D_{out}}$, the complexity is $O(N' D_{out})$. This cost is negligible compared to reservoir dynamics in long sequences, since $D_{out} \ll N$.

The computational complexity of ESN primarily depends on the reservoir step. This cost dominates due to the recurrent nature of the update, which must be computed sequentially $T$ times. However, it is important to note that the use of sparse matrix can dramatically reduce the effective time complexity compared to dense ones.

## 3 Diagonalization-Based Optimization

In this section, we introduce an optimization leveraging the diagonalization of the reservoir matrix $W$. By decomposing $W = PDP^{-1}$, we reformulate the reservoir step into a pointwise form, where each state component evolves independently according to its corresponding eigenvalue. This simplification drastically reduce the amount of computation needed to compute the activity and output of a linear ESN.
## 3.1 Core Transformation

**Theorem 1.** Change-of-basis of the ESN dynamics

(i) Let $W$ be the reservoir matrix and $P \in GL_n(\mathbb{C})$ a basis. Let us denote as $[\cdot]_P$ the transformation into the basis $P$, and express both weights and state of the ESN into this new basis.

Transformed reservoir matrix $[W]_P = P^{-1}WP$

Transformed reservoir state $[r(t)]_P = r(t)P$

Transformed input weights $[W_{in}]_P = W_{in}P$

Transformed feedback weights $[W_{fb}]_P = W_{fb}P$

$$ \text{Transformed output weights} \quad [W_{\text{out}}]_P = \begin{bmatrix} W_{\text{out,bias}} \\ W_{\text{out,out}} \\ [W_{\text{out,res}}]_P \end{bmatrix} \stackrel{\text{def}}{=} P^{-1} W_{\text{out,res}} = \begin{bmatrix} I & 0 \\ 0 & P^{-1} \end{bmatrix} W_{\text{out}} $$

(ii) Then the reservoir step becomes:

$$ [r(0)]_P = 0 \tag{10} $$

$$ [r(t)]_P = [r(t-1)]_P [W]_P + u(t) [W_{in}]_P + y(t-1) [W_{fb}]_P \tag{11} $$

(iii) And the readout step becomes:

$$ y(t) = W_{\text{out,bias}} + y(t-1) W_{\text{out,out}} + [r(t)]_P [W_{\text{out,res}}]_P \tag{12} $$

(iv) Using the notation of $X(t) = [1 \ y(t-1) \ r(t)]$, we define:

$$ [X(t)]_P = \begin{bmatrix} 1 & 0 \\ 0 & P \end{bmatrix} \quad (13) $$

It is then possible to learn the readout weights directly in $P$:

$$ [W_{\text{out}}]_P = \left( [X]_P^T [X]_P + \alpha \begin{bmatrix} I & 0 \\ 0 & P^T P \end{bmatrix} \right)^{-1} [X]_P^T Y \tag{14} $$

## Proof

$$ \begin{align*} [r(t)]_P &= r(t) P \\ &= r(t-1) W P + u(t) W_{\text{in}} P + y(t-1) W_{\text{fb}} P \\ &= r(t-1) P [W]_P P^{-1} P + u(t) [W_{\text{in}}]_P + y(t-1) [W_{\text{fb}}]_P \\ &= [r(t-1)]_P [W]_P + u(t) [W_{\text{in}}]_P + y(t-1) [W_{\text{fb}}]_P \end{align*} $$

$$ \begin{align*} y(t) &= W_{\text{out,bias}} + y(t-1) W_{\text{out,out}} + r(t) W_{\text{out,res}} \\ &= W_{\text{out,bias}} + y(t-1) W_{\text{out,out}} + r(t) P P^{-1} W_{\text{out,res}} \\ &= W_{\text{out,bias}} + y(t-1) W_{\text{out,out}} + [r(t)]_P [W_{\text{out,res}}]_P \end{align*} $$
$$
\begin{align*}
[W_{\text{out}}]_P &= \begin{bmatrix} I & 0 \\ 0 & P^{-1} \end{bmatrix} W_{\text{out}} \\
&= \begin{bmatrix} I & 0 \\ 0 & P^{-1} \end{bmatrix} (X^T X + \alpha I)^{-1} X^T Y \\
&= \begin{bmatrix} I & 0 \\ 0 & P^{-1} \end{bmatrix} \left( \left( [X]_P \begin{bmatrix} I & 0 \\ 0 & P^{-1} \end{bmatrix} \right)^T \left( [X]_P \begin{bmatrix} I & 0 \\ 0 & P^{-1} \end{bmatrix} \right) + \alpha I \right)^{-1} \left( [X]_P \begin{bmatrix} I & 0 \\ 0 & P^{-1} \end{bmatrix} \right)^T Y \\
&= \begin{bmatrix} I & 0 \\ 0 & P^{-1} \end{bmatrix} \left( \begin{bmatrix} I & 0 \\ 0 & (P^T)^{-1} \end{bmatrix} [X]_P^T [X]_P \begin{bmatrix} I & 0 \\ 0 & P^{-1} \end{bmatrix} + \alpha I \right)^{-1} \begin{bmatrix} I & 0 \\ 0 & (P^T)^{-1} \end{bmatrix} [X]_P^T Y \\
&= \begin{bmatrix} I & 0 \\ 0 & P \end{bmatrix}^{-1} \left( \begin{bmatrix} I & 0 \\ 0 & P^T \end{bmatrix}^{-1} [X]_P^T [X]_P \begin{bmatrix} I & 0 \\ 0 & P \end{bmatrix}^{-1} + \alpha I \right)^{-1} \begin{bmatrix} I & 0 \\ 0 & P^T \end{bmatrix}^{-1} [X]_P^T Y \\
&= \left( \begin{bmatrix} I & 0 \\ 0 & P^T \end{bmatrix} \left( \begin{bmatrix} I & 0 \\ 0 & P^T \end{bmatrix}^{-1} [X]_P^T [X]_P \begin{bmatrix} I & 0 \\ 0 & P \end{bmatrix}^{-1} + \alpha I \right) \begin{bmatrix} I & 0 \\ 0 & P \end{bmatrix} \right)^{-1} [X]_P^T Y \\
&= \left( \begin{bmatrix} I & 0 \\ 0 & P^T \end{bmatrix} \begin{bmatrix} I & 0 \\ 0 & P^T \end{bmatrix}^{-1} [X]_P^T [X]_P \begin{bmatrix} I & 0 \\ 0 & P \end{bmatrix}^{-1} \begin{bmatrix} I & 0 \\ 0 & P \end{bmatrix} + \alpha \begin{bmatrix} I & 0 \\ 0 & P^T \end{bmatrix} \begin{bmatrix} I & 0 \\ 0 & P \end{bmatrix} \right)^{-1} [X]_P^T Y \\
&= \left( [X]_P^T [X]_P + \alpha \begin{bmatrix} I & 0 \\ 0 & P^T P \end{bmatrix} \right)^{-1} [X]_P^T Y
\end{align*}
$$

<figure><img src="image_3.png" /></figure>

## 3.2 A special case of a diagonalizable matrix

**Corollary 2. Diagonalization Optimization**

Assume the reservoir matrix $W$ is diagonalizable. Let $W = P [W]_P P^{-1}$ where:

* $P \in GL_n(\mathbb{C})$: Matrix of eigenvectors

* $[W]_P = \text{diag}(\Lambda \stackrel{\text{def}}{=} (\lambda_1, \dots, \lambda_N))$: Diagonal matrix of eigenvalues

Then the reservoir step become pointwise:

$$ [r(t)]_P = [r(t-1)]_P \odot \Lambda + u(t) [W_{in}]_P + y(t-1) [W_{fb}]_P \quad (15) $$

where $\odot$ denotes element-wise multiplication.

**Proof** Applying the change of basis to the standard reservoir update equation yields

$$ [r(t)]_P = [r(t-1)]_P [W]_P + u(t) [W_{in}]_P + y(t-1) [W_{fb}]_P $$

Since $[r(t-1)]_P$ is a row vector and $[W]_P = \text{diag}(\Lambda)$, their matrix product strictly simplifies to the element-wise multiplication $[r(t-1)]_P \odot \Lambda$. ■

<figure><img src="image_4.png" /></figure>
### 3.3 Apply $W_{in}$ after Reservoir Steps

This subsection demonstrates that with our method the input matrix $W_{in}$ can actually be applied after the temporal update. This proves that the network depends merely on $\Lambda$, further justifying our focus on optimizing the spectral properties. For what follows, we suppose a diagonal reservoir without feedback weights ($W_{fb} = 0$) and with an initial state $r(0) = 0$.

**Lemma 3.**

$$r(t) = \sum_{i=1}^{t} (u(i)W_{in}) \odot \Lambda^{t-i}$$

where the power of $\Lambda$ denotes element-wise exponentiation, not a matrix product.

**Proof** By recurrence on $t$. For $t = 0$: $r(0) = 0$ by definition. The property holds. Assume the property is true for a given $t \ge 0$. For $t + 1$:

$$ 
\begin{aligned}
r(t+1) &= r(t) \odot \Lambda + u(t+1)W_{in} \\
&= \left( \sum_{i=1}^{t} (u(i)W_{in}) \odot \Lambda^{t-i} \right) \odot \Lambda + u(t+1)W_{in} \odot \Lambda^0 \\
&= \sum_{i=1}^{t} (u(i)W_{in}) \odot \Lambda^{t+1-i} + u(t+1)W_{in} \odot \Lambda^{t+1-(t+1)} \\
&= \sum_{i=1}^{t+1} (u(i)W_{in}) \odot \Lambda^{t+1-i}
\end{aligned}
$$

The property holds for $t + 1$, concluding the proof. ■

According to this formula, to compute the reservoir state, we first apply $W_{in}$ to the input signal and then apply $\Lambda$. Actually, this can be done in the reverse order. Let $R(t) \in \mathbb{C}^{D_{in} \times N}$ be a state matrix defined by:

$$ 
\begin{aligned}
R(0) &= 0 \\
R(t) &= R(t-1) \odot \begin{bmatrix} \Lambda \\ \vdots \\ \Lambda \end{bmatrix} + [u(t)^T \quad \dots \quad u(t)^T]
\end{aligned}
$$

Here, the first term broadcasts the row vector $\Lambda$ across all $D_{in}$ rows, and the second term repeats the column vector $u(t)^T$ across all $N$ columns to form a $D_{in} \times N$ matrix.

**Lemma 4.** Let $row_d(\cdot)$ denote the $d$-th row of a matrix.

$$row_d(R(t)) = \sum_{i=1}^{t} u(i)_d \Lambda^{t-i}$$
**Proof** By recurrence on $t$ for a fixed row $d$. For $t = 0$: $row_d(R(0)) = 0$ by definition. The property holds. Assume true for $t \ge 0$. For $t + 1$:

$$ 
\begin{aligned}
\mathrm{row}_d(R(t+1)) &= \mathrm{row}_d(R(t)) \odot \Lambda + u(t+1)_d \Lambda^0 \\
&= \left( \sum_{i=1}^{t} u(i)_d \Lambda^{t-i} \right) \odot \Lambda + u(t+1)_d \Lambda^0 \\
&= \sum_{i=1}^{t+1} u(i)_d \Lambda^{t+1-i}
\end{aligned}
$$

<figure><img src="image_5.png" /></figure>

**Theorem 5.** The standard reservoir state $r(t)$ can be recovered by element-wise weighting and summation:

$$ r(t) = \sum_{d=1}^{D_{in}} \mathrm{row}_d(W_{in}) \odot \mathrm{row}_d(R(t)) $$

Or more compactly, letting **1** be a column vector of ones of size $D_{in} \times 1$:

$$ r(t) = \mathbf{1}^T(W_{in} \odot R(t)) $$

**Proof** Using the linearity of the summation and detailing the intermediate associative steps:

$$ 
\begin{aligned}
r(t) &= \sum_{i=1}^{t} (u(i)W_{in}) \odot \Lambda^{t-i} \\
&= \sum_{i=1}^{t} \left( \sum_{d=1}^{D_{in}} u(i)_d \mathrm{row}_d(W_{in}) \right) \odot \Lambda^{t-i} \\
&= \sum_{i=1}^{t} \sum_{d=1}^{D_{in}} (u(i)_d \mathrm{row}_d(W_{in})) \odot \Lambda^{t-i} \\
&= \sum_{d=1}^{D_{in}} \sum_{i=1}^{t} \mathrm{row}_d(W_{in}) \odot (u(i)_d \Lambda^{t-i}) \\
&= \sum_{d=1}^{D_{in}} \mathrm{row}_d(W_{in}) \odot \left( \sum_{i=1}^{t} u(i)_d \Lambda^{t-i} \right) \\
&= \sum_{d=1}^{D_{in}} \mathrm{row}_d(W_{in}) \odot \mathrm{row}_d(R(t)) \\
&= \mathbf{1}^T(W_{in} \odot R(t))
\end{aligned}
$$

<figure><img src="image_6.png" /></figure>
### 3.3.1 APPLICATION: DIAGONAL LINEAR ESN

Let $(\Lambda \in \mathbb{C}^N, W_{in} \in \mathbb{C}^{D_{in} \times N}, W_{out} \in \mathbb{C}^{N \times D_{out}})$ be a diagonal linear ESN, fed with an input signal $(u(t))_{1 \le t \le T}$. According to the previous theorem:

$$y(t) = r(t)W_{out} = [\mathbf{1}^T(W_{in} \odot R(t))] W_{out}$$

While this formulation is mathematically exact for any dimensions, optimizing $W_{out}$ directly from $R(t)$ without fixing $W_{in}$ first remains complex when $D_{in} > 1$ or $D_{out} > 1$. However, this establishes a fundamental property: the general temporal dynamics are entirely captured by $R(t)$ and $\Lambda$, completely independently of the input weight $W_{in}$. This property becomes highly actionable when $D_{in} = 1$ and $D_{out} = 1$ (see Appendix C).

## 3.4 Computational Advantages

This approach reduces computational costs after initial preprocessing:

* **Eigendecomposition:** $O(N^3)$ (one-time preprocessing cost), compared to the spectral radius computation, which can be between $O(N^2)$ and $O(N^3)$ as explained above.

* $[W_{in}]_P$ and $[W_{fb}]_P$ computation: $O(N^2 D_{in})$ and $O(N^2 D_{fb})$

* **Transformed reservoir step:** $O(N(D_{in} + D_{out}))$ (reduced from $O(N(c_r N + c_{in} D_{in} + c_{fb} D_{out}))$ with $D_{in}, D_{out} \ll N$ (reservoir looking to project the inputs into a high dimensional space).

* **Transformed readout step:** $O(N' D_{out})$ (unchanged)

The key advantage is the reduction of the reservoir step complexity from $O(N^2)$ to $O(N(D_{in} + D_{out})$ per time step, making this approach highly efficient despite the initial $O(N^3)$ preprocessing cost. Figure 2 shows the comparison between standard linear reservoir computation and our optimisation in function of the number of neurons in the reservoir for each processing step. We can also note that even with the use of complex numbers (which should doubles the number of parameters with a real and an imaginary part), the training of the readout can be performed with real matrices (see Appendix A) as in the standard computations, leading to no additional cost.
<figure><img src="image_7.png" /><figcaption>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 W. Diagonalization generates a standard W 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.</figcaption></figure>

# 4 Methods

This section presents three methods for efficient linear reservoir design and training. The first, Eigenbasis Weight Transformation (EWT), enables efficient inference by transforming pre-trained weights $W_{out}$ into the eigenbasis without retraining. The second, End-to-End Eigenbasis Training (EET), enables full end-to-end training by directly computing $[W_{out}]_P$ in the transformed space, which drastically reduces computational cost during training. Finally, Direct Parameter Generation (DPG) eliminates the need for initial reservoir matrix W by directly constructing its spectral components $\Lambda$ and eigenvector basis $P$ and then computing $[W_{out}]_P$ in the transformed space.

## 4.1 Training Data Preparation

Given input-output pairs $\{(u(t), y(t))\}_{t=1,...,T_k}$, we construct:

$$X(t) = [1 \ y(t-1) \ r(t)] \in \mathbb{R}^{N'} \quad (\text{extended state}) \qquad (16)$$

$$Y(t) = [y(t)] \in \mathbb{R}^{D_{out}} \quad (\text{target output}) \qquad (17)$$

where T is the length of the sequence.
## 4.2 Method 1: Eigenbasis Weight Transformation (EWT)

If $W_{out}$ has already been computed using standard ridge regression:

$$W_{out} = (X^T X + \alpha I)^{-1} X^T Y \quad (18)$$

Then the transformed output weights are (see Theorem 1 in 2):

$$[W_{out}]_P = \begin{bmatrix} I & 0 \\ 0 & P^{-1} \end{bmatrix} W_{out} \quad (19)$$

## 4.3 Method 2: End-to-End Eigenbasis Training (EET)

If $W_{out}$ has not already been computed, it is possible to train directly $[W_{out}]_P$.

### Proposition: Direct Output Weight Computation

Let $[X]_P = [1 \quad y(t-1) \quad [r(t)]_P] \in \mathbb{C}^{N'}$ be the concatenated transformed reservoir states. It is then possible to obtain $W_{out}$ by training $[W_{out}]_P$ directly.

$$[W_{out}]_P = ([X]_P^T [X]_P + \alpha \begin{bmatrix} I & 0 \\ 0 & P^T P \end{bmatrix})^{-1} [X]_P^T Y \quad (20)$$

## 4.4 Method 3: Direct Parameter Generation (DPG)

Rather than generating $(W, W_{in})$ and then diagonalizing like the EET method, this section explores directly generating $(\Lambda, [W_{in}]_P, P)$ to avoid the initial diagonalization cost. We called that method Direct Parameter Generation (DPG).
<figure><img src="image_8.png" /><figcaption>Figure 3: Comparison of eigenvalue distributions in the complex plane. On the first column (left) the spectrum derived from a standard random reservoir matrix $W$. 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 spectrum 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).</figcaption></figure>

### 4.4.1 UNIFORM DISTRIBUTION

The spectral structure of the reservoir weight matrix $W$ plays a fundamental role in determining both the stability and the dynamics of its activity. A well-designed spectrum enables rich temporal dynamics, essential for tasks involving sequence modeling and prediction. Drawing upon insights from random matrix theory (particularly the classical work of Edelman and Kostlan (1995)), we observe that for an $N \times N$ random matrix with i.i.d. Gaussian entries, the expected number of real eigenvalues scales as:

$$ \mathbb{E}[N_{\text{real}}] \sim \sqrt{\frac{2N}{\pi}} \qquad (21) $$

We propose an approach to directly construct $P$ and $\Lambda$ without the need of generating $W$ by explicitly engineering the eigenvector matrix $P$ such that its associated eigenvalue spectrum $\Lambda$ matches a desired distribution with a controlled number of real and complex conjugate pairs. To this end, we introduce two complementary algorithms that respectively generate a valid $\Lambda$ and $P$.

Algorithm 1 generates a structured, random eigenvalue spectrum $\Lambda = (\lambda_1, \dots, \lambda_N)$ by balancing the number of real and complex conjugate pairs. Specifically, it ensures that
<caption>**Algorithm 1 Random generation of the eigenvalues**</caption>

1: **function** RANDOMEIGENVALUES($N$, sr)
2: $N_{real} \leftarrow \lfloor \sqrt{\frac{2N}{\pi}} \rfloor$
3: **if** $N_{real} \neq n \bmod 2$ **then**
4: $N_{real} \leftarrow N_{real} + 1$
5: **end if**
6: $\Lambda_{real}$: array of size $N_{real}$
7: $R2$: array of size $(N - N_{real})/2$
8: $\Theta$: array of size $(N - N_{real})/2$
9: $\Lambda_{cpx}$: array of size $(N - N_{real})/2$
10: Sample $\Lambda_{real} \sim \text{Uniform}(-\text{sr}, \text{sr})$
11: Sample $U \sim \text{Uniform}(0, 1)$
12: Sample $\Theta \sim \text{Uniform}(0, \pi)$
13: $\Lambda_{cpx} \leftarrow \text{sr} \cdot \sqrt{U} \cdot \exp(i \cdot \Theta)$
14: **return** $\Lambda_{real}, \Lambda_{cpx}, \text{conj}(\Lambda_{cpx})$
15: **end function**

exactly $N_{real}$ eigenvalues are real and since $W$ is real, the remaining $(N - N_{real})$ eigenvalues are $(N - N_{real})/2$ complex conjugate pairs. To guarantee this property, the algorithm enforces that $N_{real}$ and $N$ have the same parity (by incrementing $N_{real}$ by one if necessary), so that $N - N_{real}$ is always even and the complex eigenvalues can consistently be arranged into conjugate pairs. The radial component of each complex pair is sampled from $\sim \sqrt{U}$ (square root of a uniform distribution), while the angular component is uniformly drawn from $[0, \pi)$. An example of this distribution is provided in Figure 3.

<caption>**Algorithm 2 Random generation of the eigenvectors**</caption>

1: **function** RANDOMEIGENVECTORS($N, N_{real}$)
2: $P$: array of shape $(N, N)$
3: **for** $i \leftarrow 0$ to $N_{real} - 1$ **do**
4: &nbsp;&nbsp;&nbsp;&nbsp;Sample $v \sim \text{Normal}(0, 1)$
5: &nbsp;&nbsp;&nbsp;&nbsp;$v \leftarrow v / \|v\|_2$
6: &nbsp;&nbsp;&nbsp;&nbsp;$P[i, i] \leftarrow v$
7: **end for**
8: **for** $k \leftarrow 0$ to $(N - N_{real})/2 - 1$ **do**
9: &nbsp;&nbsp;&nbsp;&nbsp;Sample $v_R \sim \text{Normal}(0, 1)$
10: &nbsp;&nbsp;&nbsp;&nbsp;Sample $v_I \sim \text{Normal}(0, 1)$
11: &nbsp;&nbsp;&nbsp;&nbsp;$v \leftarrow (v_R + iv_I) / \|v_R + iv_I\|_2$
12: &nbsp;&nbsp;&nbsp;&nbsp;$P[i, N_{real} + k] \leftarrow v$
13: &nbsp;&nbsp;&nbsp;&nbsp;$P[i, (N + N_{real})/2 + k] \leftarrow \text{conj}(v)$
14: **end for**
15: **return** $P$
16: **end function**
Algorithm 2 constructs a random matrix of eigenvectors $P \in GL_n(C)$. Since $W$ is real, we sample eigenvectors associated with real eigenvalues and with complex-conjugate eigenvalue pairs in different ways. For real eigenvalues, we sample real-valued eigenvectors with entries drawn from a standard Gaussian distribution. For complex-conjugate eigenvalue pairs, we sample the first eigenvector as a complex-valued vector whose real and imaginary parts have entries drawn from a standard Gaussian distribution, and the second one as its conjugate. This sampling ensures that the resulting set of vectors forms a basis so that $P$ is invertible.

### 4.4.2 GOLDEN DISTRIBUTION

<caption>**Algorithm 3 Random generation of the golden eigenvalues**</caption>

1: **function GOLDENEIGENVALUES(N, sr)**
2: $N_{\text{real}} \leftarrow \lfloor \sqrt{\frac{2N}{\pi}} \rfloor$
3: $N_{\text{real}} \leftarrow (N - N_{\text{real}}) \bmod 2$
4: $N_{\text{cpx}} \leftarrow \frac{n-n_{\text{real}}}{2}$
5:
6: $\Lambda_{\text{real}}$: array of size $N_{\text{real}}$
7: Sample $\Lambda_{\text{real}} \sim \text{Uniform}(-1, 1)$
8:
9: $\Lambda_{\text{cpx}}$: array of size $N_{\text{cpx}}$
10: $i \leftarrow 0$
11: $k \leftarrow 0$
12: $v \leftarrow \text{Uniform}(0, 2)$
13: **while** $i < N_{\text{cpx}}$ **do**
14: $k \leftarrow k + 1$
15: $v \leftarrow (v + 3 - \sqrt{5}) \bmod 2$
16: **if** $v < 1$ **then**
17: $\Lambda_{\text{cpx}}(i) \leftarrow \sqrt{\frac{k}{2n_{\text{cpx}}}} \exp(i\pi v)$
18: $i \leftarrow i + 1$
19: **end if**
20: **end while**
21:
22: $\lambda \leftarrow \frac{\text{sr}}{\max(|\Lambda_{\text{real}}|, |\Lambda_{\text{cpx}}|)}$
23: $\Lambda_{\text{real}} \leftarrow \lambda \cdot \Lambda_{\text{real}}$
24: $\Lambda_{\text{cpx}} \leftarrow \lambda \cdot \Lambda_{\text{cpx}}$
25:
26: Noise: array of size $N_{\text{cpx}}$
27: Sample Noise $\sim \text{Normal}(0, \sigma) + i \cdot \text{Normal}(0, \sigma)$
28: $\Lambda_{\text{cpx}} \leftarrow \Lambda_{\text{cpx}} + \text{Noise}$
29:
30: **return** $\Lambda_{\text{real}}, \Lambda_{\text{cpx}}, \text{conj}(\Lambda_{\text{cpx}})$
31: **end function**
Algorithm 3 introduces a deterministic variation for generating the eigenvalue spectrum, aiming to improve the coverage of the complex plane compared to *Uniform Distribution*. While the partition between real and complex eigenvalues follows the same statistical scaling as Method 3 (where $N_{real} \approx \sqrt{2N/\pi}$), the placement of complex eigenvalues utilizes a deterministic sequence based on the Golden Ratio. Specifically, the complex eigenvalues are arranged in a phyllotaxis spiral pattern. The angular position is updated using a step derived from the Golden Angle (via the term $3-\sqrt{5}$), while the magnitude scales with the square root of the index to ensure a constant density distribution across the unit disk. This approach mitigates the clustering and gaps often observed with pseudo-random generators, providing a more homogeneous spectral structure. An example of this distribution is provided in Figure 3.

# 5 Experiments and Results

In this section, we evaluate our proposed methods on two standard benchmarks: the Multiple Superimposed Oscillators (MSO) and Memory Capacity (MC). Our experiments demonstrate that our approaches (notably Direct Parameter Generation (DPG) combined with the deterministic Golden distribution) achieve equivalent or superior performance compared to standard baselines, while also exposing the structural limitations of matrix diagonalization at extremely low connectivity levels.

## 5.1 Multiple Superimposed Oscillators

<figure><img src="image_9.png" /><figcaption>Figure 4: Illustration of the Multiple Superimposed Oscillators (MSO) time series for the $K = 5$ task (MSO5). The target signal $U_5(t)$ 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.</figcaption></figure>
The Multiple Superimposed Oscillators (MSO) task (Gallicchio et al., 2017; Koryakin et al., 2012) is a well-established benchmark used to evaluate the ability of recurrent networks to capture and reproduce complex, multi-frequency temporal dynamics. It consists of a time series generated by summing K discrete-time sinusoidal components:

$$U_K(t) = \sum_{k=1}^{K} \sin(\alpha_k t) \qquad (22)$$

where K denotes the total number of superimposed sinusoids and $\alpha_k$ represents their respective angular frequencies. Following the setup introduced by Gallicchio et al. (2017), we utilize a defined set of 12 distinct frequencies for the $\alpha_k$ parameters: $\alpha_1 = 0.2, \alpha_2 = 0.331, \alpha_3 = 0.42, \alpha_4 = 0.51, \alpha_5 = 0.63, \alpha_6 = 0.74, \alpha_7 = 0.85, \alpha_8 = 0.97, \alpha_9 = 1.08, \alpha_{10} = 1.19, \alpha_{11} = 1.27, \alpha_{12} = 1.32$. To systematically evaluate the performance of our models across varying levels of signal complexity, we define 12 distinct sequential tasks based on the number of components $K \in \{1, 2, \dots, 12\}$. For a given task K, the target signal $U_K$ is the sum of the first K sinusoids. The network is trained to predict the next time step $U_K(t+1)$ given the current ground truth input $U_K(t)$. For each of the 12 tasks, the generated time series is partitioned into three distinct period: a training set of $T_{train} = 400$ time steps, a validation set of $T_{valid} = 300$ time steps, and a test set of $T_{test} = 300$ time steps. To mitigate the impact of the arbitrary initial state of the reservoir, the first 100 steps of the training phase are discarded as a washout (or warmup) and are not used to compute the readout weights. To ensure the statistical robustness of our comparative analysis, the evaluation for each task $U_K$ is averaged over 10 seeds. For each test, an exhaustive grid search is conducted to determine the optimal configuration. The complete set of hyperparameter considered is detailed in Table 1. Furthermore, by leveraging Theorem 5, we significantly accelerate the hyperparameter search for all diagonal reservoir methods. Because the reservoir states can be computed independently of the input scaling factor, we only need to collect these states once per seed rather than recomputing them for each of the three input scaling values evaluated in our grid search. This mathematical property effectively divides the state computation time by a factor of three compared to the standard Normal baseline.

<table><caption>Table 1: Hyper-parameters values considered for model selection on the MSO tasks reported from Gallicchio et al. (2017)</caption><thead><tr><td>Hyper-parameter</td><td>Values considered</td></tr></thead><tbody><tr><td>Reservoir size (N)</td><td>100</td></tr><tr><td>Input scaling (scale<sub>in</sub>)</td><td>0.01, 0.1, 1</td></tr><tr><td>Leaking rate (lr)</td><td>0.1, 0.3, 0.5, 0.7, 0.9, 1.0</td></tr><tr><td>Spectral radius (ρ)</td><td>0.1, 0.3, 0.5, 0.7, 0.9, 1.0</td></tr><tr><td>Ridge regression regularization (α)</td><td>10<sup>-11</sup>, 10<sup>-10</sup>, ..., 10<sup>0</sup></td></tr></tbody></table>
<table><caption>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 W, Diagonalized to the EET method, others corresponds to the DPG method with different initialization strategies.</caption><thead><tr><th>Task</th><th>Normal</th><th>Diagonalized</th><th>Uniform Dist.</th><th>Golden Dist.<br>(σ = 0)</th><th>Noisy Golden Dist.<br>(σ = 0.2)</th><th>Sim Dist.</th></tr></thead><tbody><tr><td>MSO1</td><td>1.65e-14</td><td><b>1.58e-14</b></td><td>5.85e-14</td><td>2.49e-14</td><td>4.77e-14</td><td>3.56e-14</td></tr><tr><td>MSO2</td><td>2.55e-13</td><td>2.78e-13</td><td>2.28e-13</td><td><b>1.45e-13</b></td><td>2.39e-13</td><td>2.44e-13</td></tr><tr><td>MSO3</td><td>5.42e-12</td><td>9.14e-12</td><td><b>4.49e-12</b></td><td>9.07e-12</td><td>6.14e-12</td><td>8.37e-12</td></tr><tr><td>MSO4</td><td>1.39e-10</td><td>5.77e-10</td><td>3.64e-10</td><td>7.22e-11</td><td><b>6.93e-11</b></td><td>2.28e-10</td></tr><tr><td>MSO5</td><td>2.75e-09</td><td>4.03e-08</td><td>2.95e-08</td><td><b>5.24e-10</b></td><td>1.63e-09</td><td>1.87e-08</td></tr><tr><td>MSO6</td><td><b>7.38e-09</b></td><td>2.54e-08</td><td>2.07e-08</td><td>1.07e-08</td><td>1.18e-08</td><td>6.16e-08</td></tr><tr><td>MSO7</td><td><b>2.96e-08</b></td><td>9.48e-08</td><td>7.16e-08</td><td>5.98e-08</td><td>5.36e-08</td><td>8.55e-08</td></tr><tr><td>MSO8</td><td><b>2.75e-08</b></td><td>9.68e-08</td><td>3.57e-07</td><td>1.15e-07</td><td>6.44e-08</td><td>1.41e-07</td></tr><tr><td>MSO9</td><td><b>4.98e-08</b></td><td>2.69e-07</td><td>4.33e-07</td><td>1.69e-07</td><td>1.03e-07</td><td>1.63e-07</td></tr><tr><td>MSO10</td><td>4.65e-07</td><td>3.32e-07</td><td>4.15e-07</td><td>2.31e-07</td><td><b>1.61e-07</b></td><td>2.73e-07</td></tr><tr><td>MSO11</td><td>5.62e-07</td><td>7.38e-07</td><td>1.85e-06</td><td>7.49e-07</td><td><b>3.16e-07</b></td><td>3.71e-07</td></tr><tr><td>MSO12</td><td>9.71e-07</td><td>2.98e-06</td><td>1.34e-06</td><td>1.01e-06</td><td><b>8.44e-07</b></td><td>2.63e-06</td></tr></tbody></table>

The results presented in Table 2 show an interesting distribution of the best scores (in bold) across the different evaluated methods. This diversity suggests that our optimization approaches, particularly those based on Direct Parameter Generation (DPG) with Golden Distribution, perform at least as well as the Normal method (standard linear ESN). The slight RMSE variations from one method to another are naturally explained by the structural differences in the weights: since the DPG approach bypasses the generation of an explicit matrix W, the eigenvalue spectrum and the eigenvector configuration inevitably differ from those of a standard random matrix. Across the 12 MSO tasks, *Diagonal Golden Noised* ($\sigma$ = 0.2) achieves the best score 4 times (MSO 4, 10, 11, 12), in par with the Normal baseline, which also dominates 4 tasks (MSO 6, 7, 8, 9). If we strictly isolate these two methods for a direct comparison, they are in a perfect tie (6 to 6). This strong competitiveness indicates that our strategy of generating a noisy "Golden" spectrum constitutes a highly relevant alternative to classical initialization. Conversely, it can be noted that the *Diagonal Sim* method is the only one to never attain the top rank, although it closely follows the general trend with scores near those of the leaders. Furthermore, the comparison among DPG methods reveals that the addition of noise to the *Golden Distribution* consistently outperforms the strict *Golden Distribution* version. This phenomenon can be understood intuitively. Visually and statistically, the injection of noise allows the deterministic Golden distribution to smooth its structure and much more faithfully mimic the natural empirical eigenvalue distribution of a true random matrix. These results validate the DPG approach paired with a *Noisy Golden Distribution*. They empirically confirm that it is possible to directly generate a structured eigenvalue distribution combined with random eigenvectors to guarantee a performance level equivalent, and occasionally superior, to the standard method.
<figure><img src="image_10.png" /><figcaption>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 $W_{out}$ (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.</figcaption></figure>

## 5.2 Memory Capacity

The Memory Capacity (MC), introduced by Jaeger Jaeger (2001), is a fundamental metric used to quantify the short-term memory of a recurrent neural network. It measures the reservoir’s ability to reconstruct past continuous inputs from its current internal state. For an independent and identically distributed (i.i.d.) input sequence $u(t)$ and a trained output unit $y_k(t)$, the determination coefficient for a specific delay $k$ is given by the squared correlation coefficient:

$$d(u(t-k), y_k(t)) = \frac{\operatorname{cov}^2(u(t-k), y_k(t))}{\sigma^2(u(t))\sigma^2(y_k(t))} \qquad (23)$$

The $k$-delay short-term memory capacity of the network is then formally defined as the maximum of this determination coefficient optimized over the output weights $w_k^{out}$:

$$MC_k = \max_{w_k^{out}} d[w_k^{out}](u(t-k), y_k(t)) \qquad (24)$$

In our experiments, we utilize this established task to test the practical memory capacity of our methods. We use reservoirs without a leaking rate and with the spectral radius strictly set to 1. We then evaluate both the $MC_k$ scores and the computational duration of the different initialization methods across various values of reservoir units $N$ and requested delay $k$.
<figure><img src="image_11.png" /><figcaption>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 $W$ (Normal), Direct Parameter Generation (DPG) utilizing a *Uniform Distribution*, DPG utilizing the deterministic *Golden Distribution*, and DPG utilizing the *Sim Distribution* which uses eigenvalues extracted from a standard random matrix $W$ combined with randomly generated eigenvectors $P$. 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.</figcaption></figure>
<figure><img src="image_12.png" /><figcaption>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 $N = 100, 300, 600$, and $1000$. 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 choose 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 *W* becomes very sparse, the network loses the dimensionality required to maintain temporal representations. Moreover, the results highlight a specific minimum connectivity threshold which depends on the reservoir size. Below this connectivity, the *Diagonalization* method begins to systematically underperform the *Normal* baseline. At such extreme sparsity, the standard matrix *W* 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.</figcaption></figure>
# 6 Conclusion

Recent work in reservoir computing has highlighted two complementary trends that motivate the present study. On the one hand, several contributions have stressed the need to reduce the computational cost of reservoir-based models, either by simplifying the dynamics or by replacing recurrent computations with alternative feature constructions (Gauthier et al., 2021). On the other hand, a number of theoretical and empirical studies have shown that linear reservoirs already possess strong dynamical and memory-related properties, and can constitute competitive building blocks for temporal processing (Dambre et al., 2012; Gonon and Ortega, 2019; Voelker et al., 2019; Tino, 2020; Verzelli et al., 2021). In this context, our objective is not to revisit the modeling capabilities of linear reservoirs, but to revisit their implementation, and to investigate how their recurrent dynamics can be preserved while being expressed in a computationally more suitable form.

In this study, we showed that the dynamics of Linear Echo State Networks can be reformulated through matrix diagonalization, turning costly recurrent updates into efficient pointwise operations. This diagonalization-based optimization preserves the predictive power of standard linear ESNs while reducing computational complexity from $O(N^2)$ to $O(N)$. Beyond the theoretical speedup, our empirical results demonstrate that both the Eigenbasis Weight Transformation (EWT) and End-to-End Eigenbasis Training (EET) retain numerical stability and accuracy, with negligible differences from standard linear ESN outputs. We further introduced the Direct Parameter Generation (DPG) approach, which bypass diagonalization entirely by directly sampling eigenvalues and eigenvectors, allowing a better control of the eigenvalues that controls the dynamics of the reservoir, while replicating or even outperforming the behavior of classical Linear ESNs.

Our investigation into the diagonalization of Linear ESNs reveals that the behavior of these models is almost entirely governed by their spectral properties. The success of the Direct Parameter Generation (DPG) approach, which bypasses the generation of an explicit weight matrix $W$ to directly sample eigenvalues, suggests that the explicit matrix structure of the reservoir is secondary to the distribution of its eigenvalues. This finding aligns with the observation that structured spectral coverages, such as the deterministic "Golden" distribution, can yield performance superior to purely stochastic initializations, particularly for tasks requiring long-term memory, though its superiority for tasks with different dynamic requirements remains an open question.

This spectral perspective is further reinforced by analyzing the trained readout weights in the eigenbasis. As illustrated in Figure 5, we observe a strong heterogeneity in the contribution of spectral components: only a subset of eigenvalues is associated with large output weights, indicating that the readout mechanism naturally selects the specific frequency components or decay rates relevant to the task. This fundamental reliance on the spectrum is mathematically corroborated by our demonstration that the input weights $W_{in}$ can be applied entirely after the recurrent updates. Since the temporal states depend strictly on the eigenvalues $\Lambda$, it proves that the memory and dynamical richness of the network are governed by its spectral configuration. Consequently, DPG positions itself as a rational approach to reservoir initialization rather than a mere computational convenience. It allows for the explicit selection of spectral components, shifting the paradigm from arbitrary matrix generation to controlled spectral initialization.
While our method offers significant advantages, it is important to note that the diagonalization of a pre-existing matrix (EWT) relies on the assumption that $W$ is diagonalizable and incurs a one-time $O(n^3)$ preprocessing cost. Although this cost is amortized over long sequences, the DPG framework elegantly circumvents this bottleneck by avoiding matrix diagonalization entirely, making it highly scalable. Finally, the independent evolution of state variables in the eigenbasis has implications beyond sequential efficiency. Since each variable evolves independently, our formulation allows for the parallelization of state updates across the temporal dimension (see Appendix B). This property mirrors recent advances in state-space models like Mamba (Gu and Dao, 2024) or parallelized LMUs (Chilkuri and Eliasmith, 2021), suggesting that diagonalized Linear ESNs could leverage similar hardware optimizations. For future application, it could be interesting to adapt the methods with non linear readout (Gonon and Ortega, 2019).

These contributions suggest a more efficient use of linear reservoirs avoiding unnecessary computations, and a shift of paradigm in linear reservoir initialization towards the direct selection of spectral distributions. In future works, it would be interesting to explore the various possible distributions for the eigenvalues, in order to systematically investigate their impact on reservoir dynamics and model performance.
# Appendix A. Practical Implementation: Memory Reinterpretation

To maximize computational efficiency, we implement a hybrid approach: the reservoir update is performed by reinterpreting part of memory as complex for element-wise operations. Conversely, the training and readout steps remain strictly in the real domain to avoid supplementary complexity.

## A.1 Basis Construction

Assume the reservoir matrix $W \in \mathbb{R}^{n \times n}$ is diagonalizable. Let $P \in \text{GL}_n(\mathbb{C})$ be the matrix of eigenvectors and $\Lambda$ the diagonal matrix of eigenvalues, partitioned as follows:

* $\Lambda = \text{diag}(\lambda_1, \dots, \lambda_{n_r}, \mu_1, \bar{\mu}_1, \dots, \mu_{n_i}, \bar{\mu}_{n_i})$, where $\lambda_j \in \mathbb{R}$ and $\mu_k \in \mathbb{C} \setminus \mathbb{R}$.

* $P = [u_1, \dots, u_{n_r}, v_1, \bar{v}_1, \dots, v_{n_i}, \bar{v}_{n_i}]$ is the associated matrix of eigenvectors.

* $\Lambda_{\text{real}} = (\lambda_1, \dots, \lambda_{n_r})$ and $\Lambda_{\text{cpx}} = (\mu_1, \dots, \mu_{n_i})$.

Let define $Q$ as:

$$Q = [u_1, \dots, u_{n_r}, \text{Re}\{v_1\}, \text{Im}\{v_1\}, \dots, \text{Re}\{v_{n_i}\}, \text{Im}\{v_{n_i}\}]$$

Applying Theorem 1 with the change of basis matrix $Q$ yields the transformed system. However, in the $Q$-basis, the reservoir matrix $[W]_Q = Q^{-1}WQ$ is no longer diagonal; it is block-diagonal. Performing the update $[r(t-1)]_Q[W]_Q$ would typically require $2 \times 2$ matrix multiplications for the complex components, losing the efficiency of the diagonal representation.

## A.2 The Memory View Trick

To bypass the computational bottleneck of block-diagonal matrix multiplications, we propose an efficient implementation based on memory reinterpretation. The core idea is to process the state update in the complex domain where $[W]_P$ is purely diagonal, but to project it back into the real $Q$-basis for the readout phase. Since modern numerical libraries store complex numbers as pairs of contiguous real floating-point values, we treat the exact same physical memory block as either real or complex depending on the operation context. This allows us to perform the complex part of the reservoir update in the "complex world" and the real part of the reservoir update in the "real world" with fast element-wise operations (⊙), while keeping the readout training only in the "real world." While processing the state $r$ entirely in complex arithmetic would normally double the number of parameters required for the readout (due to separate real and imaginary weights), our approach guarantees that the readout cost remains perfectly equivalent to that of a standard real-valued ESN.

We define $[r]_Q^{\text{real}}$ and $[r]_Q^{\text{cpx}}$ as two different slices of the memory block $[r]_Q$, with $[r]_Q^{\text{cpx}}$ interpreted either as a list of couples reals number, or a list of complex number (containing real and imaginary part):

* $[r]_Q^{\text{real}} = [r]_Q[1 : n_r]$

* $[r]_Q^{\text{cpx}} = [r]_Q[n_r + 1 : n] \cdot \text{view}(\mathbb{C})$
Using this dual representation, the reservoir state update can be computed efficiently through separate element-wise multiplications for the real and complex partitions, before accumulating the input and feedback terms in the real domain.

**Reservoir Update Step:**

$$[r]_Q^{\text{real}} \leftarrow [r]_Q^{\text{real}} \odot \Lambda_{\text{real}}$$

$$[r]_Q^{\text{cpx}} \leftarrow [r]_Q^{\text{cpx}} \odot \Lambda_{\text{cpx}}$$

$$[r]_Q \leftarrow [r]_Q + u(t)[W_{\text{in}}]_Q + y(t-1)[W_{\text{fb}}]_Q$$

**Proof**

We define the transformation matrix $\mathcal{Z} = \text{diag}(I_{n_r}, Z, \dots, Z)$, where $Z$ is the change-of-basis matrix mapping a complex conjugate pair to its real and imaginary parts:

$$Z = \frac{1}{2} \begin{bmatrix} 1 & 1 \\ -i & i \end{bmatrix}$$

Then, $Q = P\mathcal{Z}$. By definition of the transformation bases, the relation $[r]_Q = rQ = rP\mathcal{Z} = [r]_P\mathcal{Z}$ always holds. We can expand this explicitly:

$$[r]_Q = [[r]_P^{(\lambda_1)} \quad \dots \quad [r]_P^{(\lambda_{n_r})} \quad \text{Re}\{[r]_P^{(\mu_1)}\} \quad \text{Im}\{[r]_P^{(\mu_1)}\} \quad \dots \quad \text{Re}\{[r]_P^{(\mu_{n_i})}\} \quad \text{Im}\{[r]_P^{(\mu_{n_i})}\}]$$

When applying the memory views, the data partitions correspond exactly to the real eigenvalues and the complex eigenvalues in the original P-basis:

$$[r]_Q^{\text{real}} \stackrel{\text{def}}{=} ([r]_Q)[1:n_r] = [[r]_P^{(\lambda_1)} \quad \dots \quad [r]_P^{(\lambda_{n_r})}]$$

$$[r]_Q^{\text{cpx}} \stackrel{\text{def}}{=} ([r]_Q)[n_r + 1:n_r].\text{view}(\mathbb{C}) = [[r]_P^{(\mu_1)} \quad \dots \quad [r]_P^{(\mu_{n_i})}]$$

Because the memory perfectly aligns with the real and imaginary components of the P-basis projection, performing the element-wise multiplications $[r]_Q^{\text{real}} \leftarrow [r]_Q^{\text{real}} \odot \Lambda_{\text{real}}$ and $[r]_Q^{\text{cpx}} \leftarrow [r]_Q^{\text{cpx}} \odot \Lambda_{\text{cpx}}$ on the specific views is strictly equivalent to performing the full complex element-wise multiplication in the original P-basis:

$$[r]_P \leftarrow [r]_P \odot \Lambda \qquad (25)$$

By reverting to the matrix multiplication notation, this operation directly corresponds to the state transition in the Q-basis:

$$[r]_Q \leftarrow [r]_Q[W]_Q \qquad (26)$$

Finally, by accumulating the input real value and feedback projections directly onto the real memory block $[r]_Q$, the complete reservoir step is correctly and efficiently evaluated as:

$$[r]_Q \leftarrow [r]_Q[W]_Q + u(t)[W_{\text{in}}]_Q + y(t-1)[W_{\text{fb}}]_Q \qquad (27)$$

<figure><img src="image_13.png" /></figure>
## A.3 Readout and Training

While the reservoir step partly exploits complex arithmetic to achieve optimal parallel efficiency, the readout and training are performed entirely in the Q-basis. Following Theorem 1 (iii) and (iv), the readout weights are learned via:

$$y(t) = W_{\text{out,bias}} + y(t-1) W_{\text{out,out}} + [r(t)]_Q [W_{\text{out,res}}]_Q \quad (28)$$

$$[W_{\text{out}}]_Q = ([X]_Q^T [X]_Q + \alpha \begin{bmatrix} I & 0 \\ 0 & Q^T Q \end{bmatrix})^{-1} [X]_Q^T Y \quad (29)$$

## Appendix B. Parallelization

We focus on Linear ESN instead of Non-Linear ESN because the sequential computation in Non-Linear ESNs is difficult to parallelize: each state $x(t)$ is obtained by applying a nonlinear transformation to the previous state $x(t-1)$. One can formally unroll the dynamics to express $x(t)$ directly in terms of past inputs, but this does not by itself enable parallel computation. To parallelize across time, the unrolled expression would need to be algebraically simplified to eliminate intermediate states. For generic nonlinear activations, however, such a simplification is not generally available. With polynomial activations, the size of the resulting closed-form expression can grow exponentially with sequence length, making this approach impractical. Even for a simple activation function like ReLU, Montúfar et al. (2014) demonstrated that it can lead to an exponential number of linear regions, thereby constraining the potential for algebraic simplification in non-linear models. However, with linear reservoirs, each state can be expressed as:

$$x(t) = \sum_{i=1}^{t} u(i) W_{\text{in}} W^{t-i} \quad (30)$$

where the contribution of each input $u_i$ to state $x_t$ can be computed independently. This independence allows each input's influence (or echo) to be evaluated separately at every time step, facilitating extensive parallelization on modern GPU architectures and leading to significant speedups, particularly for long sequences.

## Appendix C. Special case of $D_{\text{in}} = D_{\text{out}} = 1$

This section follows the reasoning present in section 3.3.1. Suppose now that the dimension of the input signal and the dimension of the output signal are $D_{\text{in}} = D_{\text{out}} = 1$. Since $D_{\text{in}} = 1$, $R(t)$ is effectively a row vector $1 \times N$.
$$
\begin{align*}
y(t) &= r(t)w_{\text{out}} \\
&= (W_{\text{in}} \odot R(t))w_{\text{out}} \\
&= \sum_{p=1}^{N} R(t)_p (w_{\text{in}})_p (w_{\text{out}})_p \\
&= R(t)(w_{\text{in}}^T \odot w_{\text{out}})
\end{align*}
$$

Let $\mathcal{L}$ be the objective loss function (e.g., Mean Squared Error) and $Y(t)$ the truth signal.
The standard readout $w_{\text{out}}$ is typically trained such that:

$$
w_{\text{out}} = \underset{\beta \in \mathbb{C}^{N \times 1}}{\operatorname{argmin}} \mathcal{L} \left( (r(t)\beta, Y(t))_{1 \le t \le T} \right)
$$

**Note:** If an L2 penalty (Ridge regression) is applied to regularize the readout, changing the parameterization alters the problem. Penalizing the transformed composite weights instead of $\|w_{\text{out}}\|^2$ will yield a different effective regularization, meaning standard Ridge regression is not perfectly equivalent under this transformation. However, for unregularized optimization, the following holds:

**Theorem 6.** If no weight in $w_{in}$ is 0, then the optimization can be performed directly on the unweighted states $R(t)$:

$$
\gamma^* = \underset{\gamma \in \mathbb{C}^{N \times 1}}{\operatorname{argmin}} \mathcal{L} \left( (R(t)\gamma, Y(t))_{1 \le t \le T} \right)
$$

Once the optimal $\gamma^*$ is found, the original output weights can be uniquely recovered via
$w_{\text{out}} = \gamma^* \oslash w_{\text{in}}^T$, where $\oslash$ denotes element-wise division.

**Proof** Let $\gamma = w_{in}^T \odot w_{out}$. Since $w_{in}$ contains no zeros, there is a bijection between $w_{out}$ and $\gamma$. The output equation becomes:

$$
\begin{align*}
y(t) &= R(t)(w_{\text{in}}^{T} \odot w_{\text{out}}) \\
&= R(t)\gamma
\end{align*}
$$

Thus, we can optimize for $\gamma$ directly using the states $R(t)$, which do not depend on $w_{in}$. $\blacksquare$

Under these hypotheses, the effective readout weights can be learned without explicitly
instantiating $w_{in}$ during the state collection phase.

**Acknowledgments and Disclosure of Funding**

Acknowledgments
References

Narsimha Reddy Chilkuri and Chris Eliasmith. Parallelizing legendre memory unit training. In *International conference on machine learning*, pages 1898–1907. PMLR, 2021.

Joni Dambre, David Verstraeten, Benjamin Schrauwen, and Serge Massar. Information processing capacity of dynamical systems. *Scientific reports*, 2(1):514, 2012.

Alan Edelman and Eric Kostlan. How many zeros of a random polynomial are real? *Bulletin of the American Mathematical Society*, 32(1):1–37, 1995.

Claudio Gallicchio, Alessio Micheli, and Luca Pedrelli. Hierarchical temporal representation in linear reservoir computing. In *Italian Workshop on Neural Nets*, pages 119–129. Springer, 2017.

Daniel J Gauthier, Erik Bollt, Aaron Griffith, and Wendson AS Barbosa. Next generation reservoir computing. *Nature communications*, 12(1):5564, 2021.

Lukas Gonon and Juan-Pablo Ortega. Reservoir computing universality with stochastic inputs. *IEEE transactions on neural networks and learning systems*, 31(1):100–112, 2019.

Albert Gu and Tri Dao. Mamba: Linear-time sequence modeling with selective state spaces. In *First conference on language modeling*, 2024.

Herbert Jaeger. Short term memory in echo state networks. 2001.

Danil Koryakin, Johannes Lohmann, and Martin V Butz. Balanced echo state networks. *Neural Networks*, 36:35–45, 2012.

Richard B Lehoucq, Danny C Sorensen, and Chao Yang. ARPACK users' guide: solution of large-scale eigenvalue problems with implicitly restarted Arnoldi methods. SIAM, 1998.

Guido Montúfar, Razvan Pascanu, Kyunghyun Cho, and Yoshua Bengio. On the number of linear regions of deep neural networks. In Z. Ghahramani, M. Welling, C. Cortes, N. Lawrence, and K.Q. Weinberger, editors, *Advances in Neural Information Processing Systems*, volume 27. Curran Associates, Inc., 2014. URL https://proceedings.neurips.cc/paper_files/paper/2014/file/fa6f2a469cc4d61a92d96e74617c3d2a-Paper.pdf.

Robert Simon Fong, Boyu Li, and Peter Tiño. Linear simple cycle reservoirs at the edge of stability perform fourier decomposition of the input driving signals. *Chaos: An Interdisciplinary Journal of Nonlinear Science*, 35(4), 2025.

Peter Tino. Dynamical systems as temporal feature spaces. *Journal of Machine Learning Research*, 21(44):1–42, 2020.

Pietro Verzelli, Cesare Alippi, Lorenzo Livi, and Peter Tiño. Input-to-state representation in linear reservoirs dynamics. *IEEE Transactions on Neural Networks and Learning Systems*, 33(9):4598–4609, 2021.
Aaron Voelker, Ivana Kajić, and Chris Eliasmith. Legendre memory units: Continuous-time representation in recurrent neural networks. *Advances in neural information processing systems*, 32, 2019.