The RLS Algorithm in Python: Recursive Least Squares, Its Equivalence to the Kalman Filter, and the Forgetting Factor

Introduction

https://yuhi-sa.github.io/en/posts/20260310_adaptive_filter/1/ and https://yuhi-sa.github.io/en/posts/20260501_lms_nlms/1/ covered LMS and NLMS adaptive filters in depth. Both articles ended by noting that “RLS converges faster than LMS but at higher computational cost” and that “RLS is mathematically equivalent to a Kalman filter,” without deriving either claim. This article follows up: we derive RLS (Recursive Least Squares) from the matrix inversion lemma, implement it from scratch in Python, numerically confirm its exact equivalence to a Kalman filter under specific conditions, and quantify how the forgetting factor \(\lambda\) trades off steady-state accuracy against tracking speed.

Problem Setup: Exponentially Weighted Least Squares

Unlike LMS, which stochastically approximates the instantaneous gradient, RLS exactly minimizes the weighted sum of squared errors over all past observations. For an FIR filter of length \(M\) with coefficients \(\mathbf{w}\) , define the cost function at time \(n\) as

where \(\lambda \in (0, 1]\) is the forgetting factor, exponentially discounting the contribution of older samples. \(\lambda = 1\) weights all samples equally (ordinary least squares); \(\lambda < 1\) emphasizes recent samples, giving an adaptive estimate.

Setting the derivative of \(J(\mathbf{w}, n)\) with respect to \(\mathbf{w}\) to zero yields the normal equations

Inverting \(\mathbf{R}(n)\) from scratch every step costs \(O(M^3)\) . The matrix inversion lemma below reduces this to an \(O(M^2)\) recursive update.

O(M²) Recursion via the Matrix Inversion Lemma

\(\mathbf{R}(n)\) updates recursively as

Letting \(P(n) = \mathbf{R}(n)^{-1}\) and applying the Sherman-Morrison formula (matrix inversion lemma)

with \(A = \lambda \mathbf{R}(n-1)\) and \(\mathbf{u} = \mathbf{v} = \mathbf{x}(n)\) gives

The key point: \(P(n)\) is computed from \(P(n-1)\) using only matrix products and a scalar division, never a direct matrix inversion. Defining a Kalman-gain-like quantity

lets us write equation (5) compactly as \(P(n) = \frac{1}{\lambda}\bigl[P(n-1) - \mathbf{k}(n)\mathbf{x}(n)^\top P(n-1)\bigr]\) . The weight update is

Per-step cost is dominated by computing \(\mathbf{k}(n)\) , giving \(O(M^2)\) — heavier than LMS/NLMS’s \(O(M)\) .

Python Implementation

import numpy as np

def rls_filter(xs, d, M, lam=1.0, delta=1.0):
    """Adaptive filtering via the RLS algorithm

    Parameters
    ----------
    xs : ndarray, shape (N, M)
        Input vector at each time step
    d : ndarray, shape (N,)
        Desired signal
    M : int
        Filter length
    lam : float
        Forgetting factor (0 < lam <= 1)
    delta : float
        Regularization for the initial covariance (P(0) = I / delta)

    Returns
    -------
    W : ndarray, shape (N, M)
        Filter coefficient trajectory over time
    """
    w = np.zeros(M)
    P = np.eye(M) / delta
    W = np.zeros((len(d), M))
    for n in range(len(d)):
        x = xs[n]
        Px = P @ x
        k = Px / (lam + x @ Px)
        e = d[n] - x @ w
        w = w + k * e
        P = (P - np.outer(k, x @ P)) / lam
        W[n] = w
    return W

np.outer(k, x @ P) implements the rank-1 update term in equation (5); Px / (lam + x @ Px) implements the gain in equation (6).

Experiment 1: Convergence Speed vs. LMS

For an \(M=8\) -tap FIR system identification problem (SNR 20 dB), we compared LMS (step size \(\mu=0.05\) ) and RLS (\(\lambda=1\) ) learning curves averaged over 200 Monte Carlo trials.

Samples needed to reach twice the noise floor (MSE ≈ 0.0207): LMS took 57 samples, RLS took 27 — roughly half as many. This numerically confirms the textbook claim that LMS converges linearly (rate depends on the input’s eigenvalue spread) while RLS converges quadratically (reaches the normal-equation solution in the shortest possible time, independent of input statistics). The trade-off is RLS’s \(O(M^2)\) per-step cost versus LMS’s \(O(M)\) , which drives the choice between them in practice.

RLS’s Equivalence to the Kalman Filter

Recall the Kalman filter update equations (prediction, innovation, gain, update) derived in https://yuhi-sa.github.io/en/posts/20260224_kalman_filter/1/. RLS can be formulated as a special case of the Kalman filter:

RLS as a state-space model:

Writing out the Kalman filter update equations (\(K_k = P_{k|k-1}H^\top S_k^{-1}\) , etc.) under this mapping reproduces the RLS update equations (6)-(7) exactly. We implemented both and applied them to the same data to compare the coefficient trajectories directly.

def kalman_param_filter(xs, d, M, R=1.0, delta=1.0):
    w = np.zeros(M)
    P = np.eye(M) / delta
    W = np.zeros((len(d), M))
    for n in range(len(d)):
        x = xs[n]
        Px = P @ x
        S = x @ Px + R          # innovation covariance
        K = Px / S               # Kalman gain
        e = d[n] - x @ w
        w = w + K * e
        P = P - np.outer(K, x @ P)
        W[n] = w
    return W

Running both on the same \(M=4\) -tap, 500-sample system identification problem with matched \(\lambda = 1\) , \(R = 1\) , and \(\delta\) :

max |RLS - KF| over all n, all taps (lambda=1, R=1): 0.0
final RLS weights: [ 0.99811908 -0.49835497  0.29976788  0.19937108]
final KF  weights: [ 0.99811908 -0.49835497  0.29976788  0.19937108]
true weights     : [ 1.  -0.5  0.3  0.2]

The two trajectories are bit-for-bit identical at every time step and every tap, numerically confirming that RLS is a special case of the Kalman filter. This equivalence is practically useful: instead of tuning the forgetting factor, one can switch to the more general Kalman-filter framework and model non-stationarity explicitly through the process noise \(Q\) .

The Forgetting Factor λ: Steady State vs. Tracking Speed

The intuition behind \(\lambda < 1\) is that it approximates a Kalman filter with \(Q > 0\) . To test this, we constructed a non-stationary system whose true coefficients jump abruptly partway through (at sample 1000) and ran RLS with several values of \(\lambda\) .

\(\lambda = 1\) gives the smallest steady-state error (statistically most efficient, since it uses all past data equally) but tracks the change most slowly — it has barely converged to the new value even 100 samples after the jump. Conversely, \(\lambda = 0.9\) has a slightly larger steady-state error (effectively using only the last ~10 samples) but tracks changes very quickly, nearly converging within 100 samples. This experiment quantitatively confirms that \(\lambda\) directly trades off steady-state accuracy against the ability to track a non-stationary environment.

Computational Cost vs. LMS/NLMS

Reproducing the comparison table from https://yuhi-sa.github.io/en/posts/20260310_adaptive_filter/1/:

Because RLS iteratively updates the matrix \(P\) , finite-precision arithmetic can accumulate rounding errors that make \(P\) asymmetric or non-positive-definite, eventually causing divergence. In practice this is mitigated by symmetrizing \(P \leftarrow (P + P^\top)/2\) each step, or by using QR-decomposition-based RLS variants. The implementation in this article is a straightforward, pedagogical one; embedded applications running continuously for long periods need these numerical-stability safeguards.

Which to Use

Related Articles

• Adaptive Filters (LMS/RLS): Theory and Python Implementation - Broader treatment of RLS’s role, including noise-cancellation applications. This article is a focused follow-up on RLS’s derivation and its Kalman equivalence.
• LMS / NLMS Algorithms: Theory and Python Implementation - A companion article with a detailed derivation of LMS/NLMS convergence conditions and steady-state error.
• Kalman Filter: Theory and Python Implementation - The state-space filter shown numerically identical to RLS in this article.
• Nonlinear Kalman Smoothers in Python: Extended RTS Smoother and Unscented RTS Smoother - A further extension of the Kalman family; the RLS equivalence deepens the intuition for these nonlinear variants.
• Wiener Filter: Optimal Linear Filtering Theory and Python Implementation - The optimal solution that both RLS and LMS iteratively approximate.
• Autocorrelation Function: Theory and Python Implementation - \(\mathbf{R}\) in the normal equations \(\mathbf{R}(n)\mathbf{w}(n) = \mathbf{r}(n)\) is exactly the input autocorrelation matrix.
• SGD and Adam: Theory and Comparison - The deep-learning counterpart to LMS’s stochastic gradient descent framework; RLS’s second-order approach (closer to natural gradient) is a useful contrast.
• Digital Filter Design Guide: Three Axes for Selection and Comparison - Useful as a selection guide against fixed-design filters.

References

• Haykin, S. (2014). Adaptive Filter Theory (5th ed.). Pearson. Chapters 9-13.
• Sayed, A. H. (2008). Adaptive Filters. Wiley-IEEE Press. Chapter 5 (RLS).
• Kalman, R. E. (1960). A new approach to linear filtering and prediction problems. Journal of Basic Engineering, 82(1), 35-45.
• Haykin, S. (2001). Kalman Filtering and Neural Networks. Wiley. Chapter 1 (RLS-Kalman equivalence).