Recursive Least Squares

Recursive Least Squares (RLS) is a streaming, second-order training method that processes the training series one timestep at a time while converging dramatically faster than the online delta rule. It achieves this by maintaining a running estimate of the inverse correlation matrix of the reservoir state, using it to scale each weight update according to the local curvature of the loss surface.

The Training Problem

As with ridge regression, the goal is to find the readout weight matrix $W_{\text{out}} \in \mathbb{R}^{N_{\text{out}} \times N}$ such that $\hat{\mathbf{y}}_t = W_{\text{out}}\,\mathbf{x}_t$ approximates the target $\mathbf{y}_t$ over the training sequence. RLS solves this incrementally: at each timestep it updates $W_{\text{out}}$ using only the current reservoir state and target, without storing the full state matrix.

The Algorithm

RLS maintains an $N \times N$ matrix $P(t) \approx \bigl(\Phi^\top \Phi\bigr)^{-1}$, initialized as $P(0) = \delta^{-1} I$. At each timestep $t$:

1. Compute the Kalman gain:

2. Compute the prediction error:

3. Update the readout weights:

4. Update the inverse correlation matrix:

where $d(t) = \lambda + \mathbf{x}(t)^\top P(t-1)\,\mathbf{x}(t)$ is the gain denominator and $\lambda \in (0, 1]$ is the forgetting factor.

For the full derivation, see Training Theory.

Parameters

delta — Initial Prior

delta initializes $P(0) = \delta^{-1} I$. It encodes the prior uncertainty over the readout weights before any data is seen:

• Small delta (e.g., 0.001): Large initial $P$, wide prior. The algorithm adapts quickly at the start but may be noisy.
• Large delta (e.g., 1.0): Unit prior. A common default that works well for most tasks.
• Very large delta (e.g., 1000.0): Tight prior, slow initial adaptation.

delta = 1.0 is a sensible default for most reservoir tasks.

lambda — Forgetting Factor

lambda controls how much weight is given to older observations:

For most offline training scenarios where the training data is stationary, lambda = 1.0 is recommended. Set lambda < 1.0 only when the target distribution is known to drift over time.

SPIRES API

RLS training is performed with a single call to spires_train_rls():

spires_status spires_train_rls(
    spires_reservoir *r,
    const double *input_series,
    const double *target_series,
    size_t series_length,
    double delta,
    double lambda
);

Parameters:

Returns: SPIRES_OK on success, or an error status code.

Example Usage

#include <spires.h>
#include <stdio.h>
 
/* Assume reservoir r is already created */
/* Assume input_train[] and target_train[] are populated */
 
double delta  = 1.0;   /* unit prior */
double lambda = 1.0;   /* no forgetting — stationary task */
 
spires_status s = spires_train_rls(r, input_train, target_train, N_TRAIN, delta, lambda);
if (s != SPIRES_OK) {
    fprintf(stderr, "RLS training failed with status %d\n", s);
    spires_reservoir_destroy(r);
    return 1;
}
 
/* Reservoir is now trained -- run inference */
spires_reservoir_reset(r);
for (size_t t = 0; t < N_TEST; t++) {
    spires_step(r, &input_test[t * num_inputs]);
    double y[NUM_OUTPUTS];
    spires_compute_output(r, y);
}

Internal Implementation

When you call spires_train_rls(), SPIRES:

1. Resets the reservoir to its initial state and zeroes $W_{\text{out}}$.
2. Allocates the inverse correlation matrix $P \in \mathbb{R}^{N \times N}$ and workspace vectors.
3. Initializes $P = \delta^{-1} I$.
4. Processes each timestep in sequence: steps the reservoir, reads the state, and applies the four-step RLS update using BLAS dgemv, ddot, dscal, and dger.
5. Frees $P$ and workspace buffers. The trained $W_{\text{out}}$ is stored inside the reservoir.

The $P$ matrix is local to the training call and is not retained afterward. If you need to continue training from where you left off across multiple calls, use the online delta rule instead.

Memory Considerations

The dominant cost of RLS is the $P$ matrix:

• Memory: $N \times N \times 8$ bytes (double precision)
• For $N = 400$: ~1.3 MB
• For $N = 1000$: ~8 MB
• For $N = 2000$: ~32 MB

Unlike ridge regression, RLS does not store the full state matrix $\Phi$ (which costs $T \times N \times 8$ bytes). For long training sequences, RLS can therefore use substantially less memory than ridge while converging faster than the delta rule.

Comparison with Other Methods

← Ridge Regression | Online Delta Rule →