Skip to content

Readout Algorithms

The readout layer maps the reservoir state \(\mathbf{x}(t)\) to the output \(\mathbf{y}(t)\).

\[ \mathbf{y}(t) = \mathbf{W}_{out} \mathbf{x}(t) \]

Ridge Regression (Batch)

Ridge regression (also known as Tikhonov regularization) minimizes the squared error while penalizing large weights to prevent overfitting. rclib implements several strategies to solve this efficiently.

Primal Formulation (\(N \le T\))

When the number of neurons (\(N\)) is less than or equal to the number of samples (\(T\)), we solve the normal equations:

\[ \mathbf{W}_{out} = (\mathbf{X}^T \mathbf{X} + \alpha \mathbf{I})^{-1} \mathbf{X}^T \mathbf{Y} \]

Where \(\mathbf{X}\) is the \(T \times N\) state matrix, \(\mathbf{Y}\) is the \(T \times O\) target matrix, and \(\alpha\) is the regularization parameter. rclib uses optimized matrix-matrix multiplication (GEMM) to form the \(N \times N\) covariance matrix \(\mathbf{X}^T \mathbf{X}\) efficiently, leveraging multi-core parallelization.

Dual Formulation (\(N > T\))

When the reservoir is very large (\(N > T\)), the dual formulation is more efficient as it operates in the \(T \times T\) sample space:

\[ \mathbf{W}_{out} = \mathbf{X}^T (\mathbf{X} \mathbf{X}^T + \alpha \mathbf{I})^{-1} \mathbf{Y} \]

This approach significantly reduces computational complexity from \(O(N^3)\) to \(O(T^3)\) in underdetermined cases.

Adaptive Solver Selection

rclib automatically selects a solver based on the problem dimensions: - Primal Cholesky: Standard case (\(N \le T\)). - Dual Cholesky: High-dimensional case (\(N > T\)). - Implicit Conjugate Gradient: Very high-dimensional case (n_features >= 4,000) where explicit matrix formation is avoided.

When include_bias=True, rclib appends a constant feature to the state matrix. The current ridge implementation regularizes this appended bias weight in the same way as the other weights.

Recursive Least Squares (RLS) (Online)

RLS updates the weights recursively for each new data point. It maintains an inverse covariance matrix \(\mathbf{P}\).

  1. Gain Calculation: $$ \mathbf{k} = \frac{\mathbf{P} \mathbf{x}}{ \lambda + \mathbf{x}^T \mathbf{P} \mathbf{x} } $$
  2. Weight Update: $$ \mathbf{W} \leftarrow \mathbf{W} + \mathbf{k} \mathbf{e}^T $$ where \(\mathbf{e} = \mathbf{d} - \mathbf{W}^T \mathbf{x}\) is the prediction error.
  3. Covariance Matrix Update: $$ \mathbf{P} \leftarrow \lambda^{-1} (\mathbf{P} - \mathbf{k} \mathbf{x}^T \mathbf{P}) $$

rclib implements two strategies for the covariance update:

  1. Sequential Rank-1 Update: For single samples or small batches, it uses symmetric rank-1 updates to reduce computational cost: $$ \mathbf{P} \leftarrow \lambda^{-1} \left( \mathbf{P} - \frac{(\mathbf{P}\mathbf{x})(\mathbf{P}\mathbf{x})^T}{\lambda + \mathbf{x}^T \mathbf{P} \mathbf{x}} \right) $$

  2. Woodbury Rank-K Update (Mini-batch): For larger mini-batches (\(\lambda = 1.0\)), it leverages the Woodbury Matrix Identity to update the covariance matrix using matrix-matrix products (GEMM): $$ \mathbf{P}_{new} = \mathbf{P} - \mathbf{P} \mathbf{X}^T (\mathbf{I} + \mathbf{X} \mathbf{P} \mathbf{X}^T)^{-1} \mathbf{X} \mathbf{P} $$ This turns sequential \(O(B \cdot N^2)\) operations into high-density GEMM calls, significantly improving throughput on multi-core systems when the batch size \(B\) is sufficiently large.

For lambda < 1.0, rank_k_update falls back to sequential rank-1 updates so the result remains mathematically equivalent to ordinary RLS.

Least Mean Squares (LMS) (Online)

LMS is a stochastic gradient descent method.

\[ \mathbf{W} \leftarrow \mathbf{W} + \eta \mathbf{e} \mathbf{x}^T \]

Where \(\eta\) is the learning rate. It is computationally cheaper (\(O(N)\)) than RLS (\(O(N^2)\)) but typically converges slower.

When partial_fit receives more than one sample, rclib performs an averaged batch-gradient LMS update:

\[ \mathbf{W} \leftarrow \mathbf{W} + \frac{\eta}{B}\mathbf{X}^T(\mathbf{Y} - \mathbf{X}\mathbf{W}) \]

where \(B\) is the mini-batch size.