Training Theory

In reservoir computing, the recurrent weight matrix $W$ and input weight matrix $W_{\text{in}}$ are fixed at initialization. The only learned component is the readout weight matrix $W^{\text{out}}$, which maps the high-dimensional reservoir state to the desired output. Because the readout is linear, training reduces to solving a system of linear equations — a convex optimization problem with a unique global minimum and no gradient pathology.

This page derives the three training algorithms implemented in SPIRES: ridge regression (batch, offline), recursive least squares (streaming, second-order), and the online delta rule (streaming, first-order).

Problem Setup

Suppose we drive the reservoir with a training input sequence $\{u(t)\}_{t=1}^{T}$ and collect the corresponding reservoir states $\{\mathbf{x}(t)\}_{t=1}^{T}$, discarding an initial transient of $T_0$ steps to allow the reservoir to wash out its initial conditions.

We arrange the states and targets into matrices:

where $\Phi$ is the state matrix (or design matrix) and $Y$ is the target matrix. Each row of $\Phi$ is a snapshot of the reservoir state at one time step.

Ridge Regression (Batch Training)

Ordinary Least Squares

The simplest approach is to find $W^{\text{out}}$ that minimizes the sum of squared errors:

where $\|\cdot\|_F$ is the Frobenius norm. Taking the gradient and setting it to zero:

Solving:

This is the ordinary least squares (OLS) solution. It is the unique minimum of a convex quadratic, so there are no local minima or saddle points.

The Overfitting Problem

In reservoir computing, the state dimensionality $N$ is typically large (hundreds to thousands of neurons), while the number of training samples $T - T_0$ may be comparable in magnitude. When $N$ is large relative to the number of samples, OLS exhibits severe overfitting: the readout learns noise in the training data and generalizes poorly.

Mathematically, the problem is that $\Phi^\top \Phi$ becomes ill-conditioned (has near-zero eigenvalues), causing the inverse $(\Phi^\top \Phi)^{-1}$ to amplify noise. The resulting weight vector has large magnitude, indicating that the readout is fitting to tiny fluctuations in the reservoir state.

Tikhonov Regularization (Ridge Regression)

The solution is to add a penalty on the magnitude of the weights. Ridge regression adds an $L_2$ regularization term:

where $\lambda > 0$ is the regularization parameter. The penalty $\lambda \|W\|_F^2$ discourages large weights, preferring solutions that are smooth combinations of reservoir states.

Taking the gradient and setting it to zero:

Solving:

This is the ridge regression solution, also known as Tikhonov regularization in the numerical analysis literature. The term $\lambda I$ shifts all eigenvalues of $\Phi^\top \Phi$ away from zero, ensuring the matrix is well-conditioned and the inverse is stable.

Interpretation of $\lambda$

The regularization parameter $\lambda$ controls the bias-variance trade-off:

• $\lambda \to 0$: Recovers OLS. Low bias, high variance. The readout fits the training data closely but may not generalize.
• $\lambda \to \infty$: The readout approaches $W^{\text{out}} \to 0$. High bias (always predicting near zero), zero variance.
• Optimal $\lambda$: Balances training error and weight magnitude. Typically found by cross-validation or the AGILE optimizer.

The effect on the eigenspectrum of $\Phi^\top \Phi$ is transparent. If $\Phi^\top \Phi$ has eigenvalues $\sigma_1^2 \geq \sigma_2^2 \geq \cdots \geq \sigma_N^2$, then the ridge solution shrinks each component by a factor of:

Large eigenvalues (strong signal directions) are barely affected. Small eigenvalues (noise directions) are suppressed toward zero. This is precisely the behavior needed in reservoir computing, where a few principal components of the reservoir state carry signal and the rest carry noise.

Computational Cost

Computing $W^{\text{out}}_{\text{ridge}}$ requires:

1. Forming $\Phi^\top \Phi \in \mathbb{R}^{N \times N}$: $O(T N^2)$ operations.
2. Forming $\Phi^\top Y \in \mathbb{R}^{N \times N_{\text{out}}}$: $O(T N N_{\text{out}})$ operations.
3. Solving $(\Phi^\top \Phi + \lambda I) W = \Phi^\top Y$: $O(N^3)$ operations (or $O(N^2)$ per output with Cholesky factorization).

The dominant cost is $O(T N^2 + N^3)$. For typical reservoir sizes ($N \sim 100\text{--}1000$), this is fast — seconds on modern hardware. SPIRES uses LAPACKE for the linear solve, leveraging optimized BLAS routines.

SVD Alternative

An equivalent formulation uses the singular value decomposition (SVD) of $\Phi = U \Sigma V^\top$:

This form is more numerically stable when $\Phi^\top \Phi$ is extremely ill-conditioned, though the Cholesky approach is faster when applicable. SPIRES uses the normal equations (Equation 4) by default.

Online Delta Rule (Streaming Training)

Motivation

Ridge regression requires collecting all training data before computing $W^{\text{out}}$. For applications where data arrives continuously or where the target distribution changes over time, an online (incremental) training algorithm is needed.

The Delta Rule

The online delta rule updates $W^{\text{out}}$ after each time step based on the prediction error:

where:

• $\eta > 0$ is the learning rate
• $\mathbf{y}(t) = W^{\text{out}} \mathbf{x}(t)$ is the current prediction
• $\mathbf{y}_{\text{target}}(t)$ is the target at time $t$
• $\mathbf{x}(t)$ is the current reservoir state

This is a stochastic gradient descent step on the squared error loss $\|\mathbf{y}_{\text{target}} - W^{\text{out}}\mathbf{x}\|^2$. The update is:

Properties

Convergence. For a fixed reservoir driven by a stationary input, the delta rule converges to the least-squares solution as $t \to \infty$, provided $\eta$ satisfies:

In practice, $\eta$ is set to a small constant (e.g., $10^{-4}$ to $10^{-2}$).

Computational cost per step. Each update requires $O(N \cdot N_{\text{out}})$ operations — a single outer product. This is negligible compared to the $O(N^2)$ cost of the reservoir state update.

Tracking non-stationary targets. Because the delta rule continuously adapts, it can track slow changes in the target distribution. This makes it suitable for online learning and adaptive control applications.

No regularization. The basic delta rule does not include explicit regularization. Weight decay can be added:

This is the online equivalent of ridge regression, with $\lambda$ controlling weight decay.

Recursive Least Squares (Online, Second-Order)

Motivation

The delta rule converges slowly because it uses only first-order (gradient) information. Each step moves in the direction of steepest descent, but the step size is the same in every direction of the weight space regardless of how curved the loss surface is. In directions where the loss changes slowly, the delta rule makes unnecessarily cautious updates; in directions where it changes sharply, the same step size can overshoot.

Recursive Least Squares (RLS) corrects this by maintaining a running estimate of the inverse correlation matrix of the reservoir state, $P(t) \approx \bigl(\Phi^\top \Phi\bigr)^{-1}$. This allows each step to be scaled according to the local curvature of the loss surface — a second-order method — achieving the same asymptotic solution as ridge regression in far fewer steps.

Derivation

RLS can be derived as a recursive implementation of the normal equations. Suppose at time $t-1$ we have accumulated the inverse correlation matrix:

where $\lambda \in (0, 1]$ is a forgetting factor that exponentially discounts older observations, and $\delta > 0$ initializes $P(0) = \delta^{-1} I$.

When a new sample $\bigl(\mathbf{x}(t),\, \mathbf{y}_{\text{target}}(t)\bigr)$ arrives, the Sherman–Morrison–Woodbury identity gives a rank-1 update for $P$ without re-inverting the full matrix:

where $\mathbf{e}(t) = \mathbf{y}_{\text{target}}(t) - W^{\text{out}}(t-1)\,\mathbf{x}(t)$ is the prediction error and $\mathbf{k}(t) \in \mathbb{R}^N$ is the Kalman gain vector.

Symmetric Shortcut for the $P$ Update

Since $P$ is symmetric, $\mathbf{x}(t)^\top P(t-1) = \bigl(P(t-1)\,\mathbf{x}(t)\bigr)^\top$. Before normalization, $P(t-1)\,\mathbf{x}(t)$ was already computed as $\lambda\, d(t)\, \mathbf{k}(t)$, where $d(t) = \lambda + \mathbf{x}(t)^\top P(t-1)\,\mathbf{x}(t)$ is the denominator. Substituting:

so Equation 10 simplifies to a single outer product:

This avoids a second matrix-vector product and is the form used in the SPIRES implementation.

Properties

Convergence. For $\lambda = 1$ and stationary inputs, RLS converges to the ordinary least-squares solution exactly after $N$ steps (in exact arithmetic), compared to the $O(1/\eta)$ steps required by the delta rule. In practice, convergence is reached in far fewer than $T$ steps for most reservoir tasks.

Computational cost per step. Each update requires:

• One matrix-vector product $P\mathbf{x}$: $O(N^2)$
• One scalar and two rank-1 updates: $O(N^2)$

Total: $O(N^2)$ per step, dominated by the $P$ update.

Memory. Storing $P \in \mathbb{R}^{N \times N}$ requires $O(N^2)$ memory — $8N^2$ bytes in double precision. For $N = 500$, this is 2 MB; for $N = 2000$, 32 MB. This is the principal cost relative to the delta rule, which requires only $O(N \cdot N_{\text{out}})$.

Forgetting factor. When $\lambda < 1$, older samples are exponentially discounted with effective memory horizon $\approx 1/(1-\lambda)$ steps. This allows RLS to track slowly non-stationary targets, at the cost of increased variance in the weight estimates. For stationary problems, $\lambda = 1$ is recommended.

Initialization. $P(0) = \delta^{-1} I$ encodes the prior that all weight directions are equally uncertain with variance $\delta^{-1}$. A larger $\delta^{-1}$ (smaller $\delta$) produces a more diffuse prior and faster initial adaptation. A common choice is $\delta = 1$ (unit prior).

Batch vs. Streaming: Comparison

When to Use Which

Use ridge regression when:

• The training data is fixed and available in advance.
• An exact global optimum is desired.
• The dataset fits comfortably in memory ($O(TN)$ bytes).

Use RLS when:

• Data arrives as a stream but fast convergence is required.
• $O(N^2)$ memory is affordable (typical for $N \lesssim 2000$).
• The target distribution is slowly non-stationary ($\lambda < 1$).

Use the online delta rule when:

• Memory is severely constrained and $O(N^2)$ is not affordable.
• The reservoir is very large and the per-step $O(N^2)$ cost of RLS is prohibitive.
• The application is embedded or resource-constrained and convergence speed is secondary.
• The target distribution changes rapidly and continuous tracking is more important than fast convergence.

In practice, RLS is the preferred online method when memory permits. The delta rule is reserved for deployments where the $O(N^2)$ memory of RLS is not available.

Connection to Reservoir Computing Theory

The fact that only the readout layer is trained has profound theoretical implications:

Universality. A reservoir with the echo state property and a sufficiently rich state space can approximate any fading-memory filter to arbitrary accuracy, provided the readout is powerful enough. Since linear readouts are universal approximators when the reservoir provides a sufficiently rich feature expansion, the combination of a random reservoir + linear readout is a universal function approximator for temporal signals.

No vanishing/exploding gradients. Because no gradients are propagated through the recurrent dynamics, the training is immune to the vanishing and exploding gradient problems that plague backpropagation-through-time in conventional RNNs.

Convexity. The ridge regression objective (Equation 3) is a strictly convex function of $W^{\text{out}}$. There is a unique global minimum, no local minima, and no saddle points. Training is deterministic and reproducible.

Separation of dynamics and learning. The reservoir dynamics (determined by $W$, $W_{\text{in}}$, neuron model, and topology) and the learning (determined by $W^{\text{out}}$ and $\lambda$) are completely decoupled. This allows them to be designed and optimized independently, which is the architectural advantage exploited by the AGILE optimizer.

References

1. Jaeger, H. (2001). The “echo state” approach to analysing and training recurrent neural networks. GMD Report 148.
2. Lukoševičius, M., & Jaeger, H. (2009). Reservoir computing approaches to recurrent neural network training. Computer Science Review, 3(3), 127—149.
3. Tikhonov, A. N., & Arsenin, V. Y. (1977). Solutions of Ill-Posed Problems. Winston & Sons.
4. Bishop, C. M. (2006). Pattern Recognition and Machine Learning, Chapter 3. Springer.
5. Widrow, B., & Hoff, M. E. (1960). Adaptive switching circuits. IRE WESCON Convention Record, Part 4, 96—104.
6. Haykin, S. (2002). Adaptive Filter Theory (4th ed.), Chapter 13: Recursive Least-Squares Estimation. Prentice Hall.
7. Sussillo, D., & Abbott, L. F. (2009). Generating coherent sequences by learning in continuous-time recurrent networks. Neuron, 63(4), 544—557. (FORCE learning — RLS applied to reservoir readout training.)

← Spectral Radius | Case Studies: Spoken Digit Recognition →