Weighted Recursive Least Squares: A Step-by-step Tutorial in Python¶
This notebook accompanies my blog post Derivation of a Weighted Recursive Linear Least Squares Estimator. In the post, we derived update equations that refine a linear least squares fit whenever new observations arrive, without solving the whole problem again. Here, we turn these equations into Python code.
We start with the simplest case, one new observation at a time, and implement the equations exactly as they appear in the post. Then we add one ingredient after the other and motivate each of them with a small example: the initial regularization, observation weights, forgetting of old observations, batches of observations and several outputs. After each step, we compare the recursive estimate with the direct least squares solution. This comparison is the best way to convince ourselves that the implementation does what it should. At the end, we collect everything in a small, documented class and use it for a few typical scenarios.
The notebook only needs NumPy and Matplotlib. It runs locally or in Google Colab (see the badge above).
Contents
- The goal: fitting a line to a stream of data
- The recursive update for a single observation
- The initialization and the role of $\rho$
- Weighting observations
- Forgetting old observations
- Rounding errors and the symmetry of $\mathbf{A}^{-1}$
- Processing batches of observations
- Several outputs at once
- Putting it all together: the
WeightedRLSclass - Usage examples
- Summary
An optional appendix shows a QR-based variant of the update, with examples in single precision (float32); it can also improve robustness for difficult double-precision problems.
Notation. The symbols follow the blog post. The table lists the names that we use for them in the code.
| Symbol | Python | Meaning |
|---|---|---|
| $\mathbf{x}_{n+1}$, $y_{n+1}$, $w_{n+1}$ | x, y, w | New input vector, its target and its weight |
| $\boldsymbol{\theta}_n$ | theta | Current coefficients (a matrix $\boldsymbol{\Theta}_n$ with one column per output from Section 8 on) |
| $\mathbf{A}_n^{-1}$ | A_inv | Inverse of the regularized, weighted normal matrix |
| $e_{n+1}$ | e | Prediction error of the new observation, before the update |
| $\boldsymbol{\delta}_{n+1}$ | delta | Gain vector |
| $\rho$ | rho | Initial regularization (ridge) coefficient, $\mathbf{A}_0=\rho\mathbf{I}$ |
| $\lambda$, $\lambda_n$ | lam, lam_n | Forgetting factor, and the factor used in the current update ($\lambda_0=1$) |
| $\mathbf{X}_\mu$, $\mathbf{y}_\mu$, $\mathbf{W}_\mu$ | X_mu, y_mu, w_mu | A batch of $\mu$ input vectors (as rows), their targets and weights (the diagonal of $\mathbf{W}_\mu$) |
| $\mathbf{e}_\mu$, $\boldsymbol{\Delta}_\mu$ | e_mu, Delta_mu | Prediction errors and gain matrix of a batch |
A few NumPy operations appear throughout: A @ B is the matrix product (for two vectors, their dot product), A.T is the transpose, np.outer(a, b) is the matrix $\mathbf{a}\mathbf{b}^T$, and X[:n] selects the first $n$ rows of X.
import time
import numpy as np
import matplotlib.pyplot as plt
from IPython.display import display
# Colorblind-safe colors and a calm plotting style for all figures
BLUE, ORANGE, AQUA = "#2a78d6", "#eb6834", "#1baf7a"
RAMP = ["#86b6ef", "#5598e7", "#2a78d6", "#1c5cab", "#0d366b"] # light to dark = early to late
INK, MUTED, POINTS = "#52514e", "#898781", "#bdbcb6"
plt.rcParams.update({
"figure.dpi": 120, "axes.prop_cycle": plt.cycler(color=[BLUE, ORANGE, AQUA]),
"axes.spines.top": False, "axes.spines.right": False, "axes.edgecolor": "#c3c2b7",
"axes.grid": True, "grid.color": "#e1e0d9", "grid.linewidth": 0.8,
"axes.titlesize": 11, "axes.titlelocation": "left", "axes.labelcolor": INK,
"xtick.color": MUTED, "ytick.color": MUTED, "xtick.labelcolor": INK, "ytick.labelcolor": INK,
"legend.frameon": False, "legend.fontsize": 9, "lines.linewidth": 1.8,
})
print("NumPy", np.__version__)
NumPy 2.2.6
To look at the actual numbers in the examples, we use a small helper, show, which displays vectors and matrices in the familiar mathematical notation, together with their dimensions. It is not needed for the algorithm itself and only uses IPython's display function, which is available in every Jupyter environment.
import re
def format_number(value, digits=3):
"""Formats a number for the formula output of `show`.
Args:
value: The number to format.
digits: Number of decimal places.
Returns:
The number as a string. Values that would be rounded to zero and very
large values are written with a power of ten.
"""
if value != 0 and (round(abs(value), digits) == 0 or abs(value) >= 10**5):
mantissa, exponent = f"{value:.{digits - 1}e}".split("e")
return rf"{mantissa}\cdot 10^{{{int(exponent)}}}"
return f"{value:.{digits}f}"
def show(*items, digits=3):
r"""Displays named numbers, vectors and matrices together with their dimensions.
In Jupyter and Colab, the output is rendered as formulas. Viewers that
cannot render formulas show a plain-text version instead.
Args:
*items: Pairs (name, value). The name is a LaTeX string such as
r"\boldsymbol{\theta}_n", and the value is a number or a NumPy
array. Vectors (1-D arrays) are shown as columns, as in the post.
digits: Number of decimal places.
"""
formulas, texts = [], []
for name, value in items:
a = np.asarray(value, dtype=float)
if a.ndim == 0: # a single number
body, space = format_number(a.item(), digits), r"\mathbb{R}"
else:
rows = a[:, None] if a.ndim == 1 else a # a vector is shown as a column
entries = r" \\ ".join(" & ".join(format_number(v, digits) for v in row) for row in rows)
body = r"\begin{bmatrix}" + entries + r"\end{bmatrix}"
if a.ndim == 1:
space = rf"\mathbb{{R}}^{{{len(a)}}}"
else:
space = rf"\mathbb{{R}}^{{{a.shape[0]}\times {a.shape[1]}}}"
formulas.append(rf"{name} = {body} \in {space}")
# plain-text version: remove the LaTeX markup from the name, e.g. \boldsymbol{\theta}_n -> theta_n
plain = re.sub(r"\\(?:boldsymbol|mathbf|text|hat)\{([^{}]*)\}", r"\1", name)
plain = plain.replace("\\ ", " ").replace("\\", "")
texts.append(f"{plain} (shape {a.shape}):\n{np.array2string(a, precision=digits, suppress_small=True)}")
display({"text/latex": r"$\displaystyle " + r",\qquad ".join(formulas) + "$",
"text/plain": "\n".join(texts)}, raw=True)
1. The Goal: Fitting a Line to a Stream of Data¶
Imagine that pairs $(x_i, y_i)$ arrive one after another, for example measurements from a sensor, and that we would like to know the best straight line $y \approx \theta_0 + \theta_1 x$ at any time. We write each input as a vector $\mathbf{x}_i = (1, x_i)^T$, where the constant $1$ belongs to the intercept $\theta_0$, and collect the first $n$ input vectors as rows of a matrix $\mathbf{X}_n$. The coefficients that fit the first $n$ observations best are given by the closed-form (weighted and regularized) least squares solution from the blog post,
$$\boldsymbol{\theta}_n = \big(\mathbf{X}_n^T \mathbf{W}_n \mathbf{X}_n + \rho \mathbf{I}\big)^{-1} \mathbf{X}_n^T \mathbf{W}_n \mathbf{y}_n .$$
For now, all weights are $1$ ($\mathbf{W}_n=\mathbf{I}$), and $\rho$ is a small number whose role we discuss in Section 3. We implement this direct solution first, because it will serve as our reference throughout the notebook.
def batch_solution(X, y, w=None, rho=0.0):
"""Computes the direct weighted least squares solution.
Solves the regularized normal equations (X^T W X + rho I) theta = X^T W y,
where W is the diagonal matrix of the observation weights.
Args:
X: Input vectors as rows, shape (n, k).
y: Targets, shape (n,) for one output or (n, m) for m outputs.
w: Observation weights, shape (n,). Defaults to all ones.
rho: Regularization coefficient rho >= 0.
Returns:
The coefficients, shape (k,) for one output or (k, m) for m outputs.
"""
w = np.ones(len(X)) if w is None else np.asarray(w, dtype=float)
WX = w[:, None] * X # W X: row i of X times w_i (w[:, None] turns w into a column)
A = X.T @ WX + rho * np.eye(X.shape[1]) # X^T W X + rho I, shape (k, k)
b = WX.T @ y # X^T W y, shape (k,) or (k, m)
return np.linalg.solve(A, b) # solves A theta = b; same as inv(A) @ b, but more accurate
Now we create a stream of $200$ noisy observations of the line $y = 1 + 0.6\,x$ and, for every $n$, compute the direct solution from the first $n$ observations. This is what the recursion will have to reproduce later. The first rows of the matrix $\mathbf{X}_n$ show the structure of the input vectors $(1, x_i)$:
rng = np.random.default_rng(1) # random number generator with a fixed seed
n_obs = 200
theta_true = np.array([1.0, 0.6]) # true intercept and slope
x_values = rng.uniform(-3, 3, n_obs) # 200 random x-values between -3 and 3
X = np.column_stack([np.ones(n_obs), x_values]) # input vectors (1, x_i) as rows, shape (200, 2)
y = X @ theta_true + rng.normal(0, 0.8, n_obs) # X @ theta_true is one value per row; plus noise
rho = 0.01
# Refit from scratch after every new observation (this is what we want to avoid)
refits = np.array([batch_solution(X[:n], y[:n], rho=rho) for n in range(1, n_obs + 1)])
show((r"\mathbf{X}_3", X[:3]), (r"\mathbf{y}_3", y[:3])) # the first three observations
show((r"\boldsymbol{\theta}_{200}", refits[-1])) # coefficients after all observations
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(11, 4.2))
# Left: the observations and the lines fitted after 2, 5, 20 and 200 of them
grid = np.array([-3.0, 3.0]) # a line is drawn through its two end points
ax1.scatter(x_values, y, s=14, color=POINTS, label="observations")
for color, n in zip(RAMP[1:], [2, 5, 20, 200]):
ax1.plot(grid, refits[n - 1, 0] + refits[n - 1, 1] * grid, color=color, label=f"fit after {n} observations")
ax1.plot(grid, theta_true[0] + theta_true[1] * grid, "--", color=INK, lw=1.2, label="true line")
ax1.set(xlabel="$x$", ylabel="$y$", title="The fitted line improves as data arrive")
ax1.legend(loc="upper left", fontsize=8)
# Right: both coefficients as functions of the number of observations
n_axis = np.arange(1, n_obs + 1)
ax2.plot(n_axis, refits[:, 0], label=r"intercept $\theta_0$")
ax2.plot(n_axis, refits[:, 1], label=r"slope $\theta_1$")
ax2.axhline(theta_true[0], ls="--", color=INK, lw=1.2, label="true values")
ax2.axhline(theta_true[1], ls="--", color=INK, lw=1.2)
ax2.set(xscale="log", xlabel="number of observations $n$", ylabel="coefficient",
title="Estimated coefficients")
ax2.legend()
fig.tight_layout()
plt.show()
Refitting after every observation works, but it has two drawbacks: we have to store all previous observations, and every refit builds and solves the complete problem again. The recursive least squares (RLS) update produces exactly the same coefficients from the previous coefficients, the matrix $\mathbf{A}_n^{-1}$ and the new observation alone.
2. The Recursive Update for a Single Observation¶
We begin with the special case at the end of the blog post: one new observation $(\mathbf{x}_{n+1}, y_{n+1})$ at a time. To keep things simple, all weights are $1$ and nothing is forgotten ($w_{n+1}=1$ and $\lambda_n=1$). The update then consists of four steps:
$$ \begin{aligned} e_{n+1} &= y_{n+1}-\mathbf{x}_{n+1}^T\boldsymbol{\theta}_n, \\ \boldsymbol{\delta}_{n+1} &= \frac{\mathbf{A}_n^{-1}\mathbf{x}_{n+1}}{1+\mathbf{x}_{n+1}^T\mathbf{A}_n^{-1}\mathbf{x}_{n+1}}, \\ \boldsymbol{\theta}_{n+1} &= \boldsymbol{\theta}_n+\boldsymbol{\delta}_{n+1}e_{n+1}, \\ \mathbf{A}_{n+1}^{-1} &= \mathbf{A}_n^{-1}-\boldsymbol{\delta}_{n+1}\mathbf{x}_{n+1}^T\mathbf{A}_n^{-1}. \end{aligned} $$
In words: we first predict the new target with the current coefficients and measure the prediction error $e_{n+1}$. The gain vector $\boldsymbol{\delta}_{n+1}$ then determines how strongly each coefficient reacts to this error, and the coefficients are corrected accordingly. Finally, $\mathbf{A}^{-1}$ is updated for the next observation. The recursion starts with $\boldsymbol{\theta}_0=\mathbf{0}$ and $\mathbf{A}_0^{-1}=\rho^{-1}\mathbf{I}$.
The Python function is a line-by-line translation:
def rls_update(theta, A_inv, x, y):
"""Performs one recursive least squares update (unit weight, no forgetting).
Args:
theta: Current coefficients theta_n, shape (k,).
A_inv: Current inverse A_n^{-1}, shape (k, k).
x: New input vector x_{n+1}, shape (k,).
y: New target y_{n+1} (a number).
Returns:
A tuple (theta, A_inv) with the new coefficients theta_{n+1}, shape
(k,), and the new inverse A_{n+1}^{-1}, shape (k, k).
"""
e = y - x @ theta # prediction error e_{n+1}; x @ theta is x^T theta
delta = A_inv @ x / (1.0 + x @ A_inv @ x) # gain vector delta_{n+1}, shape (k,)
theta = theta + delta * e # theta_{n+1} = theta_n + delta_{n+1} e_{n+1}
A_inv = A_inv - np.outer(delta, x @ A_inv) # A_{n+1}^{-1} = A_n^{-1} - delta_{n+1} (x_{n+1}^T A_n^{-1})
return theta, A_inv
One update in numbers. Before we run the recursion on the whole stream, we look at the very first update in detail. We start with $\boldsymbol{\theta}_0=\mathbf{0}$ and $\mathbf{A}_0^{-1}=\rho^{-1}\mathbf{I}=100\,\mathbf{I}$, write out the four steps of rls_update one by one, and display every intermediate result:
theta_0 = np.zeros(2) # theta_0 = 0
A_inv_0 = np.eye(2) / rho # A_0^{-1} = I / rho
x_1, y_1 = X[0], y[0] # the first observation
e_1 = y_1 - x_1 @ theta_0 # prediction error
delta_1 = A_inv_0 @ x_1 / (1.0 + x_1 @ A_inv_0 @ x_1) # gain vector
theta_1 = theta_0 + delta_1 * e_1 # new coefficients
A_inv_1 = A_inv_0 - np.outer(delta_1, x_1 @ A_inv_0) # new inverse
show((r"\mathbf{x}_1", x_1), (r"y_1", y_1), (r"e_1", e_1))
show((r"\boldsymbol{\delta}_1", delta_1), (r"\boldsymbol{\theta}_1", theta_1), (r"\mathbf{A}_1^{-1}", A_inv_1))
show((r"\mathbf{x}_1^T\boldsymbol{\theta}_1", x_1 @ theta_1)) # the prediction for x_1 after the update
Since $\boldsymbol{\theta}_0=\mathbf{0}$, the first prediction is $0$, and the prediction error equals the target. After the update, the model already reproduces the first target almost exactly. However, $\mathbf{A}_1^{-1}$ still contains large entries: a single observation only provides information along the direction of $\mathbf{x}_1$, and in the perpendicular direction, the coefficients are still only constrained by the regularization.
Now we feed all observations from Section 1 into the recursion, one at a time, and store the coefficients after every update:
theta = np.zeros(2) # theta_0 = 0
A_inv = np.eye(2) / rho # A_0^{-1} = I / rho
path = [] # the coefficients after each update
for x_i, y_i in zip(X, y): # zip walks through the rows of X and the entries of y in parallel
theta, A_inv = rls_update(theta, A_inv, x_i, y_i)
path.append(theta)
path = np.array(path) # shape (200, 2): one row per update
show((r"\boldsymbol{\theta}_{200}\ \text{(recursive)}", theta), (r"\boldsymbol{\theta}_{200}\ \text{(direct)}", refits[-1]),
digits=6)
show((r"\mathbf{A}_{200}^{-1}", A_inv))
What does a single update do? Each panel below shows the line before an update, the new observation with its prediction error $e_{n+1}$, and the line after the update.
fig, axes = plt.subplots(1, 3, figsize=(11, 3.7), sharey=True)
grid = np.array([-3.0, 3.0])
for ax, n in zip(axes, [2, 10, 100]):
before, after = path[n - 1], path[n] # theta_n and theta_{n+1} (path[0] is theta_1)
x_new, y_new = x_values[n], y[n] # the (n+1)-th observation (Python counts from 0)
y_pred = before[0] + before[1] * x_new # its prediction before the update
ax.scatter(x_values[:n], y[:n], s=14, color=POINTS, label="previous observations")
ax.plot(grid, before[0] + before[1] * grid, color=BLUE, label="line before the update")
ax.plot(grid, after[0] + after[1] * grid, color=AQUA, label="line after the update")
ax.plot([x_new, x_new], [y_pred, y_new], color=ORANGE, lw=1.4, label="prediction error $e_{n+1}$")
ax.scatter([x_new], [y_new], s=50, color=ORANGE, edgecolor="white", linewidth=1.5, zorder=3,
label="new observation")
ax.set(title=f"Update {n} → {n + 1}", xlabel="$x$", ylim=(-3.5, 5.5))
axes[0].set_ylabel("$y$")
axes[0].legend(loc="upper left", fontsize=8)
fig.tight_layout()
plt.show()
Early on, a single observation moves the line considerably. After a hundred observations, a similar prediction error hardly changes the line anymore, because the gain vector has become small: the entries of $\mathbf{A}^{-1}$ shrink as more data arrive, and the estimate becomes more certain.
Check: does the recursion reproduce the direct solution? The derivation promises that the recursion gives exactly the direct least squares solution after every single update, not only approximately at the end. We therefore compare the stored path with the refits from Section 1. The differences should be of the size of floating-point rounding errors, i.e. around $10^{-15}$ for coefficients of order one.
difference = np.abs(path - refits).max(axis=1) # largest difference of the two coefficients, per update
print(f"largest difference over all {n_obs} updates: {difference.max():.1e}")
fig, ax = plt.subplots(figsize=(6.5, 3.2))
# np.maximum(..., 1e-17) avoids log(0) for differences that are exactly zero
ax.semilogy(np.arange(1, n_obs + 1), np.maximum(difference, 1e-17), color=BLUE)
ax.set(xlabel="number of observations $n$", ylabel="largest absolute difference",
title="Recursive vs. direct solution: equal up to rounding errors")
fig.tight_layout()
plt.show()
largest difference over all 200 updates: 9.4e-16
3. The Initialization and the Role of $\rho$¶
The recursion needs a starting point. Before any data have arrived, the natural choice would be $\mathbf{A}_0=\mathbf{0}$ (there is no information yet), but a zero matrix has no inverse. The direct solution has the same problem as long as there are fewer observations than coefficients. A single point does not determine a line, since infinitely many lines pass through it, and mathematically this shows up as a singular matrix $\mathbf{X}^T\mathbf{X}$. For a single observation, its second row is a multiple of the first:
X_one, y_one = X[:1], y[:1] # a single observation, but two coefficients
show((r"\mathbf{X}_1^T\mathbf{X}_1", X_one.T @ X_one))
print("rank:", np.linalg.matrix_rank(X_one.T @ X_one), "(2 would be needed for an inverse)")
try:
batch_solution(X_one, y_one, rho=0.0) # without regularization, the solve fails
except np.linalg.LinAlgError as error:
print("without regularization:", error)
show((r"\mathbf{X}_1^T\mathbf{X}_1+\rho\mathbf{I}", X_one.T @ X_one + rho * np.eye(2)),
(r"\boldsymbol{\theta}_1", batch_solution(X_one, y_one, rho=rho)))
rank: 1 (2 would be needed for an inverse) without regularization: Singular matrix
The regularization term $\rho\mathbf{I}$ solves both problems: $\mathbf{A}_0=\rho\mathbf{I}$ is invertible with $\mathbf{A}_0^{-1}=\rho^{-1}\mathbf{I}$, and the direct solution exists from the very first observation. The price is a small bias, since $\rho$ acts like a mild preference for coefficients close to zero. Its influence fades as more observations arrive, because the data term $\mathbf{X}_n^T\mathbf{X}_n$ keeps growing while $\rho\mathbf{I}$ stays the same.
How should we choose $\rho$? The left plot below shows that a large $\rho$ holds the first estimates back towards zero. The right plot shows the opposite end: with an extremely small $\rho$, the initial matrix $\mathbf{A}_0^{-1}=\rho^{-1}\mathbf{I}$ contains huge numbers, the rounding errors of the first updates no longer vanish, and the recursion drifts away from the direct solution. For inputs of order one, values between about $10^{-3}$ and $1$ are a reasonable choice.
def run_rls(X, y, rho):
"""Processes all observations one at a time, starting from theta_0 = 0 and A_0^{-1} = I / rho.
Args:
X: Input vectors as rows, shape (n, k).
y: Targets, shape (n,).
rho: Initial regularization coefficient rho > 0.
Returns:
The coefficients after each update, shape (n, k).
"""
theta, A_inv = np.zeros(X.shape[1]), np.eye(X.shape[1]) / rho
path = []
for x_i, y_i in zip(X, y):
theta, A_inv = rls_update(theta, A_inv, x_i, y_i)
path.append(theta)
return np.array(path)
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(11, 4))
# Left: the estimated slope for three different values of rho
for rho_i in [0.01, 10.0, 100.0]:
ax1.plot(n_axis, run_rls(X, y, rho_i)[:, 1], label=rf"$\rho$ = {rho_i:g}")
ax1.axhline(theta_true[1], ls="--", color=INK, lw=1.2, label="true slope")
ax1.set(xscale="log", xlabel="number of observations $n$", ylabel=r"slope $\theta_1$",
title=r"A large $\rho$ holds the first estimates back")
ax1.legend()
# Right: the largest difference between the recursion and the direct solution, for many values of rho
rhos = np.logspace(-14, 2, 17) # 10^-14, 10^-13, ..., 10^2
deviation = []
for rho_i in rhos:
direct = np.array([batch_solution(X[:n], y[:n], rho=rho_i) for n in range(2, n_obs + 1)])
deviation.append(np.abs(run_rls(X, y, rho_i)[1:] - direct).max())
ax2.loglog(rhos, deviation, "o-", color=BLUE, markersize=4)
ax2.set(xlabel=r"regularization $\rho$", ylabel="largest difference to the direct solution",
title=r"A tiny $\rho$ amplifies rounding errors")
fig.tight_layout()
plt.show()
4. Weighting Observations¶
So far, every observation counts equally. In practice, some observations are more reliable than others. Suppose that our data come from two sensors: a precise one with noise standard deviation $\sigma=0.3$ and a cheap one with $\sigma=3$. Weighted least squares assigns each observation a weight $w_i$, and a common choice is the inverse noise variance $w_i=1/\sigma_i^2$, so that precise observations count much more.
According to the blog post, the weight $w_{n+1}$ of the new observation only changes the gain vector:
$$\boldsymbol{\delta}_{n+1} = \frac{w_{n+1}\mathbf{A}_n^{-1}\mathbf{x}_{n+1}}{1+w_{n+1}\mathbf{x}_{n+1}^T\mathbf{A}_n^{-1}\mathbf{x}_{n+1}} .$$
def rls_update(theta, A_inv, x, y, w=1.0):
"""Performs one recursive least squares update with an observation weight (no forgetting).
Args:
theta: Current coefficients theta_n, shape (k,).
A_inv: Current inverse A_n^{-1}, shape (k, k).
x: New input vector x_{n+1}, shape (k,).
y: New target y_{n+1} (a number).
w: Positive weight w_{n+1} of the new observation.
Returns:
A tuple (theta, A_inv) with the new coefficients theta_{n+1}, shape
(k,), and the new inverse A_{n+1}^{-1}, shape (k, k).
"""
e = y - x @ theta
delta = w * A_inv @ x / (1.0 + w * x @ A_inv @ x) # the only line that changes
theta = theta + delta * e
A_inv = A_inv - np.outer(delta, x @ A_inv)
return theta, A_inv
The next cell generates data from both sensors and runs the recursion with and without weights. Since a single data set can be lucky or unlucky, we also repeat the comparison on $50$ independent data sets and average the distance between the estimated and the true coefficients.
def two_sensor_data(n, rng):
"""Generates noisy observations of the line y = 1 + 0.6 x from two sensors.
Each observation comes from a precise sensor (noise standard deviation
0.3) or from a noisy one (standard deviation 3), each with probability 1/2.
Args:
n: Number of observations.
rng: NumPy random number generator.
Returns:
A tuple (X, y, sigma) with the input vectors (1, x_i) as rows, shape
(n, 2), the noisy targets, shape (n,), and the noise standard deviation
of each observation, shape (n,).
"""
x = rng.uniform(-3, 3, n)
sigma = np.where(rng.random(n) < 0.5, 0.3, 3.0) # which sensor measured each point
y = 1.0 + 0.6 * x + rng.normal(0, sigma) # each point gets the noise of its sensor
return np.column_stack([np.ones(n), x]), y, sigma
def error_path(X, y, w, rho=0.01):
"""Runs the weighted recursion and tracks the distance to the true coefficients.
Args:
X: Input vectors as rows, shape (n, 2).
y: Targets, shape (n,).
w: Observation weights, shape (n,).
rho: Initial regularization coefficient rho > 0.
Returns:
A tuple (errors, theta) with the distance between the estimate and the
true coefficients (1, 0.6) after each update, shape (n,), and the final
coefficients, shape (2,).
"""
theta, A_inv = np.zeros(2), np.eye(2) / rho
errors = []
for x_i, y_i, w_i in zip(X, y, w):
theta, A_inv = rls_update(theta, A_inv, x_i, y_i, w_i)
errors.append(np.linalg.norm(theta - theta_true)) # Euclidean distance to (1, 0.6)
return np.array(errors), theta
rng = np.random.default_rng(2)
X_s, y_s, sigma = two_sensor_data(300, rng)
show((r"\boldsymbol{\sigma}_{1:5}", sigma[:5]), (r"\mathbf{w}_{1:5}=1/\boldsymbol{\sigma}_{1:5}^2", 1 / sigma[:5]**2),
digits=2) # noise levels and weights of the first five points
n_show = 40 # the left plot shows the first 40 observations
_, theta_weighted = error_path(X_s[:n_show], y_s[:n_show], 1 / sigma[:n_show]**2)
_, theta_unweighted = error_path(X_s[:n_show], y_s[:n_show], np.ones(n_show))
show((r"\boldsymbol{\theta}_{40}\ \text{(weighted)}", theta_weighted),
(r"\boldsymbol{\theta}_{40}\ \text{(unweighted)}", theta_unweighted), (r"\boldsymbol{\theta}_{\text{true}}", theta_true))
# The same comparison on 50 independent data sets, averaged
data_sets = [two_sensor_data(300, rng) for _ in range(50)]
mean_error_weighted = np.mean([error_path(Xr, yr, 1 / sr**2)[0] for Xr, yr, sr in data_sets], axis=0)
mean_error_unweighted = np.mean([error_path(Xr, yr, np.ones(300))[0] for Xr, yr, sr in data_sets], axis=0)
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(11, 4.2))
# Left: the first 40 observations of both sensors and the two fits
precise = sigma[:n_show] < 1 # True for the points of the precise sensor
x_show, y_show = X_s[:n_show, 1], y_s[:n_show]
ax1.scatter(x_show[~precise], y_show[~precise], s=18, facecolors="none", edgecolors=MUTED, label="noisy sensor")
ax1.scatter(x_show[precise], y_show[precise], s=14, color=INK, label="precise sensor")
grid = np.array([-3.0, 3.0])
ax1.plot(grid, theta_weighted[0] + theta_weighted[1] * grid, color=BLUE, label="weighted RLS")
ax1.plot(grid, theta_unweighted[0] + theta_unweighted[1] * grid, color=ORANGE, label="unweighted RLS")
ax1.plot(grid, 1.0 + 0.6 * grid, "--", color=INK, lw=1.2, label="true line")
ax1.set(xlabel="$x$", ylabel="$y$", ylim=(-9, 11), title=f"Two sensors: fits after {n_show} observations")
ax1.legend(loc="upper left", fontsize=8, ncol=2)
# Right: the average distance to the true coefficients over 50 data sets
n_axis_s = np.arange(1, 301)
ax2.loglog(n_axis_s, mean_error_weighted, color=BLUE, label=r"weighted, $w_i = 1/\sigma_i^2$")
ax2.loglog(n_axis_s, mean_error_unweighted, color=ORANGE, label=r"unweighted, $w_i = 1$")
ax2.set(xlabel="number of observations $n$", ylabel="distance to the true coefficients",
title="Average over 50 data sets")
ax2.legend()
fig.tight_layout()
plt.show()
The weighted estimate is closer to the true line, and on average it reaches a given accuracy with far fewer observations. As before, we check that the weighted recursion reproduces the direct weighted solution after every update:
theta, A_inv, worst = np.zeros(2), np.eye(2) / rho, 0.0
w_s = 1 / sigma**2
for n, (x_i, y_i, w_i) in enumerate(zip(X_s, y_s, w_s), start=1): # n counts the observations so far
theta, A_inv = rls_update(theta, A_inv, x_i, y_i, w_i)
worst = max(worst, np.abs(theta - batch_solution(X_s[:n], y_s[:n], w_s[:n], rho)).max())
print(f"largest difference to the direct weighted solution: {worst:.1e}")
largest difference to the direct weighted solution: 3.4e-14
5. Forgetting Old Observations¶
Up to now, every observation keeps its weight forever. This is exactly right if the relationship between $x$ and $y$ never changes. In many applications, however, it does change (a sensor ages, or a process is adjusted), and old observations then describe a situation that no longer exists. Exponential forgetting multiplies the weights of all previous observations by a factor $0<\lambda\leq1$ in every update. An observation that is $a$ updates old therefore has the weight $\lambda^a$ (left plot below). These weights add up to about $1/(1-\lambda)$, which gives a useful rule of thumb for the memory of the estimator: with $\lambda=0.98$, it mainly relies on the last $50$ observations.
With forgetting, the blog post gives
$$ \begin{aligned} \boldsymbol{\delta}_{n+1} &= \frac{w_{n+1}\mathbf{A}_n^{-1}\mathbf{x}_{n+1}}{\lambda_n+w_{n+1}\mathbf{x}_{n+1}^T\mathbf{A}_n^{-1}\mathbf{x}_{n+1}}, \\ \mathbf{A}_{n+1}^{-1} &= \frac{1}{\lambda_n}\Big(\mathbf{A}_n^{-1}-\boldsymbol{\delta}_{n+1}\mathbf{x}_{n+1}^T\mathbf{A}_n^{-1}\Big), \end{aligned} $$
while $e_{n+1}$ and $\boldsymbol{\theta}_{n+1}$ are computed as before. Here, $\lambda_n=\lambda$ for all updates except the very first one, which uses $\lambda_0=1$: there is nothing to forget yet, and the initial regularization should start at its full value $\rho$.
def rls_update(theta, A_inv, x, y, w=1.0, lam_n=1.0):
"""Performs one weighted recursive least squares update with forgetting.
Args:
theta: Current coefficients theta_n, shape (k,).
A_inv: Current inverse A_n^{-1}, shape (k, k).
x: New input vector x_{n+1}, shape (k,).
y: New target y_{n+1} (a number).
w: Positive weight w_{n+1} of the new observation.
lam_n: Forgetting factor lambda_n of this update, in (0, 1].
Returns:
A tuple (theta, A_inv) with the new coefficients theta_{n+1}, shape
(k,), and the new inverse A_{n+1}^{-1}, shape (k, k).
"""
e = y - x @ theta
delta = w * A_inv @ x / (lam_n + w * x @ A_inv @ x) # lam_n replaces the 1 in the denominator
theta = theta + delta * e
A_inv = (A_inv - np.outer(delta, x @ A_inv)) / lam_n # dividing by lam_n forgets old information
return theta, A_inv
As an example, the true intercept and slope change abruptly after $400$ of $800$ observations. We run the recursion with three different forgetting factors:
rng = np.random.default_rng(3)
n_drift = 800
x_drift = rng.uniform(-3, 3, n_drift)
X_drift = np.column_stack([np.ones(n_drift), x_drift])
before_change = np.arange(n_drift)[:, None] < 400 # True for the first 400 observations
theta_drift = np.where(before_change, [1.0, 0.6], [-0.5, -0.4]) # true coefficients over time, shape (800, 2)
y_drift = np.sum(X_drift * theta_drift, axis=1) + rng.normal(0, 0.5, n_drift) # row-wise x_i^T theta_i
def run_rls_forgetting(X, y, lam, rho=0.01):
"""Processes all observations one at a time with a forgetting factor.
Args:
X: Input vectors as rows, shape (n, k).
y: Targets, shape (n,).
lam: Forgetting factor lambda in (0, 1].
rho: Initial regularization coefficient rho > 0.
Returns:
The coefficients after each update, shape (n, k).
"""
theta, A_inv = np.zeros(X.shape[1]), np.eye(X.shape[1]) / rho
path = []
for n, (x_i, y_i) in enumerate(zip(X, y)):
lam_n = 1.0 if n == 0 else lam # no forgetting in the very first update
theta, A_inv = rls_update(theta, A_inv, x_i, y_i, 1.0, lam_n)
path.append(theta)
return np.array(path)
paths = {lam: run_rls_forgetting(X_drift, y_drift, lam) for lam in [1.0, 0.98, 0.9]}
show(*[(rf"\boldsymbol{{\theta}}_{{800}}\ (\lambda={lam:g})", path_lam[-1]) for lam, path_lam in paths.items()],
(r"\boldsymbol{\theta}_{\text{true}}", theta_drift[-1])) # final estimates and the true values
fig, axes = plt.subplots(1, 3, figsize=(12, 3.9))
# Left: the weight lambda^a of an observation that is a updates old
ages = np.arange(150)
for lam in paths:
axes[0].plot(ages, lam ** ages, label=rf"$\lambda$ = {lam:g}")
axes[0].set(xlabel="age $a$ (updates since the observation)", ylabel=r"weight $\lambda^a$",
title="How much an old observation counts")
axes[0].legend()
# Middle and right: the estimated coefficients over time
for ax, j, name in [(axes[1], 0, r"Intercept $\theta_0$"), (axes[2], 1, r"Slope $\theta_1$")]:
for lam, path_lam in paths.items():
ax.plot(path_lam[:, j], lw=1.4, label=rf"$\lambda$ = {lam:g}")
ax.plot(theta_drift[:, j], "--", color=INK, lw=1.2, label="true value")
ax.set(xlabel="observation", title=name)
axes[2].legend()
fig.tight_layout()
plt.show()
Without forgetting ($\lambda=1$), the estimate after the change is an average over both situations and approaches the new values only very slowly. With $\lambda=0.98$, the estimate follows the change within about a hundred observations. With $\lambda=0.9$, it reacts even faster but is noticeably noisier, because it effectively uses only the last ten observations. Choosing $\lambda$ is therefore a trade-off between fast adaptation and low noise.
Check: which problem does the recursion solve now? With forgetting, the recursion should reproduce the direct weighted solution in which an observation of age $a$ has the weight $\lambda^a$. There is one subtlety, which the blog post discusses: the recursion forgets the initial regularization as well, because the old normal-matrix contribution, including the initial $\rho\mathbf{I}$, is multiplied by $\lambda$ in every update after the first, just like the old observation weights; the inverse is instead scaled by $1/\lambda$. After $n$ observations, the effective regularization is therefore $\rho\lambda^{n-1}$. The comparison below confirms this: the direct solution with the decayed $\rho$ agrees with the recursion up to rounding errors, whereas a fixed $\rho$ leads to small but clear differences.
lam = 0.98
path_lam = paths[lam]
diff_decayed = diff_fixed = 0.0
for n in range(1, n_drift + 1):
weights = lam ** np.arange(n - 1, -1, -1) # lambda^age for observations 1, ..., n (ages n-1, ..., 0)
decayed = batch_solution(X_drift[:n], y_drift[:n], weights, rho * lam ** (n - 1))
fixed = batch_solution(X_drift[:n], y_drift[:n], weights, rho)
diff_decayed = max(diff_decayed, np.abs(path_lam[n - 1] - decayed).max())
diff_fixed = max(diff_fixed, np.abs(path_lam[n - 1] - fixed).max())
print(f"effective regularization after {n_drift} observations: rho * lam^(n-1) = {rho * lam ** (n_drift - 1):.1e}")
print(f"direct solution with decayed rho: largest difference {diff_decayed:.1e}")
print(f"direct solution with fixed rho: largest difference {diff_fixed:.1e}")
effective regularization after 800 observations: rho * lam^(n-1) = 9.8e-10 direct solution with decayed rho: largest difference 2.4e-15 direct solution with fixed rho: largest difference 1.7e-03
6. Rounding Errors and the Symmetry of $\mathbf{A}^{-1}$¶
The matrix $\mathbf{A}_n$ is symmetric, and so is its inverse. In the update, we compute both $\mathbf{A}_n^{-1}\mathbf{x}_{n+1}$ (for the gain) and $\mathbf{x}_{n+1}^T\mathbf{A}_n^{-1}$ (for the new inverse). For a symmetric matrix, these two products contain the same numbers, so it is tempting to compute only one of them and reuse it. In exact arithmetic, this changes nothing. In floating-point arithmetic, however, the computed $\mathbf{A}^{-1}$ is never perfectly symmetric, and with forgetting this becomes a problem. Every update divides by $\lambda<1$, and in the variant with reuse, the tiny asymmetric part of the rounding errors is multiplied by $1/\lambda$ in every step, without anything that would damp it.
The following experiment runs three variants on a stream whose true coefficients never change, with $\lambda=0.95$: the literal implementation from Section 5, the variant that reuses $\mathbf{A}_n^{-1}\mathbf{x}_{n+1}$, and the same variant with a simple remedy, which replaces $\mathbf{A}^{-1}$ after every update by its symmetric part $\frac{1}{2}\big(\mathbf{A}^{-1}+(\mathbf{A}^{-1})^T\big)$.
def rls_update_reuse(theta, A_inv, x, y, w=1.0, lam_n=1.0, symmetrize=False):
"""Like rls_update, but computes A^{-1} x once and reuses it instead of x^T A^{-1}.
Args:
theta: Current coefficients theta_n, shape (k,).
A_inv: Current inverse A_n^{-1}, shape (k, k).
x: New input vector x_{n+1}, shape (k,).
y: New target y_{n+1} (a number).
w: Positive weight w_{n+1} of the new observation.
lam_n: Forgetting factor lambda_n of this update, in (0, 1].
symmetrize: If True, A_inv is replaced by its symmetric part after the update.
Returns:
A tuple (theta, A_inv) with the new coefficients and the new inverse.
"""
e = y - x @ theta
Ax = A_inv @ x # computed once ...
delta = w * Ax / (lam_n + w * x @ Ax)
theta = theta + delta * e
A_inv = (A_inv - np.outer(delta, Ax)) / lam_n # ... and reused here instead of x @ A_inv
if symmetrize:
A_inv = 0.5 * (A_inv + A_inv.T) # the remedy: keep A^{-1} exactly symmetric
return theta, A_inv
rng = np.random.default_rng(4)
n_sym, lam = 2000, 0.95
X_sym = np.column_stack([np.ones(n_sym), rng.uniform(-3, 3, n_sym)])
y_sym = X_sym @ theta_true + rng.normal(0, 0.5, n_sym)
def run_variant(variant):
"""Runs one of the three variants on the stream (X_sym, y_sym) with lambda = 0.95.
Args:
variant: "literal", "reuse" or "reuse + symmetrize".
Returns:
A tuple (asymmetry, slopes, A_inv) with the relative asymmetry of A^{-1}
after each update, shape (n,), the estimated slope after each update,
shape (n,), and the final A^{-1}, shape (2, 2).
"""
theta, A_inv = np.zeros(2), np.eye(2) / rho
asymmetry, slopes = [], []
for n, (x_i, y_i) in enumerate(zip(X_sym, y_sym)):
lam_n = 1.0 if n == 0 else lam
if variant == "literal":
theta, A_inv = rls_update(theta, A_inv, x_i, y_i, 1.0, lam_n)
else:
theta, A_inv = rls_update_reuse(theta, A_inv, x_i, y_i, 1.0, lam_n,
symmetrize=(variant == "reuse + symmetrize"))
# largest entry of A^{-1} - (A^{-1})^T, relative to the largest entry of A^{-1}
asymmetry.append(np.abs(A_inv - A_inv.T).max() / np.abs(A_inv).max())
slopes.append(theta[1])
return np.array(asymmetry), np.array(slopes), A_inv
results = {variant: run_variant(variant) for variant in ["literal", "reuse", "reuse + symmetrize"]}
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(11, 4))
for name, (asymmetry, slopes, _) in results.items():
ax1.semilogy(np.maximum(asymmetry, 1e-18), lw=1.4, label=name) # exact zeros are drawn at 1e-18
if name == "literal": # drawn wide and light, because the symmetrized variant lies on top of it
ax2.plot(slopes, lw=4, alpha=0.35, label=name)
else:
ax2.plot(slopes, lw=1.2, label=name)
ax1.set(xlabel="update", ylabel=r"relative asymmetry of $\mathbf{A}^{-1}$", ylim=(1e-19, 10),
title=r"Loss of symmetry ($\lambda$ = 0.95)")
ax1.legend()
ax2.axhline(theta_true[1], ls="--", color=INK, lw=1.2, label="true slope")
ax2.set(xlabel="update", ylabel=r"slope $\theta_1$", ylim=(-1, 2), title="Estimated slope")
ax2.legend(loc="lower left")
fig.tight_layout()
plt.show()
first_bad = np.argmax(results["reuse"][0] > 1e-3) # index of the first update above 0.1 %
print(f"reuse: asymmetry exceeds 0.1 % of A^-1 after {first_bad} updates")
show((r"\mathbf{A}^{-1}_{2000}\ \text{(literal)}", results["literal"][2]),
(r"\mathbf{A}^{-1}_{2000}\ \text{(reuse)}", results["reuse"][2]))
reuse: asymmetry exceeds 0.1 % of A^-1 after 634 updates
The literal implementation stays symmetric up to rounding errors (at about $10^{-16}$). The variant with reuse loses its symmetry exponentially fast and, after several hundred updates, produces useless estimates; its final $\mathbf{A}^{-1}$ above is far from symmetric, and its entries have exploded. With symmetrization, the asymmetry is exactly zero (it is drawn at the bottom of the left plot), and the estimates are as good as those of the literal implementation. The class in Section 9 reuses products for efficiency and therefore symmetrizes $\mathbf{A}^{-1}$ after every update; this costs almost nothing.
Symmetrization removes this asymmetry, but does not guarantee positive definiteness or prevent cancellation. For example, in double precision a one-feature update with $\rho=1$ and $x=10^8$ can round $A^{-1}$ to zero even though its exact value is about $10^{-16}$, preventing later learning. Scale the features to comparable, moderate magnitudes. The QR appendix shows an alternative that can also help difficult double-precision problems; for numerical stress tests, a direct QR- or SVD-based solution of the weighted, ridge-augmented least squares problem provides a useful reference alongside the normal-equation checks above.
7. Processing Batches of Observations¶
Observations often arrive in groups, for instance when a sensor delivers several measurements at once, or we may simply want to process a large data set faster. For a batch of $\mu$ observations, the blog post derives
$$ \begin{aligned} \mathbf{e}_\mu &= \mathbf{y}_\mu-\mathbf{X}_\mu\boldsymbol{\theta}_n, \\ \boldsymbol{\Delta}_\mu &= \mathbf{A}_n^{-1}\mathbf{X}_\mu^T\big(\lambda_n\mathbf{W}_\mu^{-1}+\mathbf{X}_\mu\mathbf{A}_n^{-1}\mathbf{X}_\mu^T\big)^{-1}, \\ \boldsymbol{\theta}_{n+\mu} &= \boldsymbol{\theta}_n+\boldsymbol{\Delta}_\mu\mathbf{e}_\mu, \\ \mathbf{A}_{n+\mu}^{-1} &= \frac{1}{\lambda_n}\Big(\mathbf{A}_n^{-1}-\boldsymbol{\Delta}_\mu\mathbf{X}_\mu\mathbf{A}_n^{-1}\Big). \end{aligned} $$
The structure is the same as before, but the prediction errors now form a vector $\mathbf{e}_\mu$, the gain becomes a $k\times\mu$ matrix $\boldsymbol{\Delta}_\mu$, and the scalar denominator turns into a $\mu\times\mu$ matrix that has to be inverted. $\mathbf{W}_\mu$ is the diagonal matrix with the $\mu$ weights of the batch. Again, the code is a direct translation:
def rls_batch_update(theta, A_inv, X_mu, y_mu, w_mu, lam_n=1.0):
"""Performs one recursive least squares update with a batch of observations.
Args:
theta: Current coefficients theta_n, shape (k,), or (k, m) for m outputs.
A_inv: Current inverse A_n^{-1}, shape (k, k).
X_mu: Input vectors of the batch as rows, shape (mu, k).
y_mu: Targets of the batch, shape (mu,), or (mu, m) for m outputs.
w_mu: Positive observation weights of the batch, shape (mu,).
lam_n: Forgetting factor lambda_n of this update, in (0, 1].
Returns:
A tuple (theta, A_inv) with the new coefficients theta_{n+mu} and the
new inverse A_{n+mu}^{-1}, with the same shapes as the inputs.
"""
e_mu = y_mu - X_mu @ theta # prediction errors, one per observation
S = lam_n * np.diag(1.0 / w_mu) + X_mu @ A_inv @ X_mu.T # the mu x mu matrix (np.diag builds W_mu^{-1})
Delta_mu = A_inv @ X_mu.T @ np.linalg.inv(S) # gain matrix, shape (k, mu)
theta = theta + Delta_mu @ e_mu
A_inv = (A_inv - Delta_mu @ X_mu @ A_inv) / lam_n
return theta, A_inv
One batch update in numbers. To see the dimensions of all quantities, we process the first three observations of Section 1 as one batch ($\mu=3$), starting from $\boldsymbol{\theta}_0$ and $\mathbf{A}_0^{-1}$ again. The result must agree with the direct solution after three observations:
theta_0, A_inv_0 = np.zeros(2), np.eye(2) / rho # start again from theta_0 and A_0^{-1}
X_mu, y_mu, w_mu = X[:3], y[:3], np.ones(3) # a batch of mu = 3 observations with unit weights
e_mu = y_mu - X_mu @ theta_0 # three prediction errors
S = np.diag(1.0 / w_mu) + X_mu @ A_inv_0 @ X_mu.T # the 3 x 3 matrix (lambda_0 = 1)
Delta_mu = A_inv_0 @ X_mu.T @ np.linalg.inv(S) # the 2 x 3 gain matrix
theta_3 = theta_0 + Delta_mu @ e_mu # coefficients after the batch
show((r"\mathbf{X}_\mu", X_mu), (r"\mathbf{y}_\mu", y_mu), (r"\mathbf{e}_\mu", e_mu))
show((r"\lambda_0\mathbf{W}_\mu^{-1}+\mathbf{X}_\mu\mathbf{A}_0^{-1}\mathbf{X}_\mu^T", S), (r"\boldsymbol{\Delta}_\mu", Delta_mu))
show((r"\boldsymbol{\theta}_3\ \text{(batch)}", theta_3), (r"\boldsymbol{\theta}_3\ \text{(direct)}", refits[2]), digits=6)
Check: batches give the same result as single observations. Without forgetting, it should not matter whether we feed the observations one at a time or in batches. We process the stream from Section 1 in batches of $\mu=10$ and compare the coefficients after each batch with the direct solution:
mu = 10
theta, A_inv, worst = np.zeros(2), np.eye(2) / rho, 0.0
for start in range(0, n_obs, mu): # start = 0, 10, 20, ...
X_mu, y_mu = X[start:start + mu], y[start:start + mu] # rows start, ..., start + 9
theta, A_inv = rls_batch_update(theta, A_inv, X_mu, y_mu, np.ones(len(X_mu)))
worst = max(worst, np.abs(theta - refits[start + len(X_mu) - 1]).max())
print(f"largest difference to the direct solution after each batch: {worst:.1e}")
largest difference to the direct solution after each batch: 5.1e-11
The differences are somewhat larger than for single observations, because the first batch starts from the large matrix $\mathbf{A}_0^{-1}=100\,\mathbf{I}$ and the $\mu\times\mu$ matrix is inverted explicitly. They are still far too small to matter in practice.
Forgetting in batch mode. With $\lambda<1$, the forgetting factor is applied once per batch. All observations of a batch therefore have the same age, and the memory is counted in batches rather than in observations. A batch recursion with $\lambda$ does not reproduce the single-observation recursion with the same $\lambda$. For a similar memory, we can use $\lambda^\mu$ per batch instead (for example $0.98^{10}\approx0.82$ for batches of ten).
Why batches? Speed. Every call of the update involves a few matrix operations, and in Python, the overhead of a call is often larger than the arithmetic itself. Processing many observations per call therefore saves time. On the other hand, the $\mu\times\mu$ matrix in the gain has to be inverted, and this cost grows like $\mu^3$. The next experiment measures how long it takes to process $20{,}000$ observations with $k=10$ features for different batch sizes. The exact times depend on the computer that runs the notebook.
Benchmark scope. The timing experiment uses the generic rls_batch_update for every batch size, including a $1\times1$ matrix inverse when $\mu=1$. Its speedup is relative to that implementation; the specialized single-observation update avoids this overhead. The final class also uses solve instead of the explicit inverse timed here.
rng = np.random.default_rng(5)
n_big, k_big = 20_000, 10
X_big = rng.normal(size=(n_big, k_big))
y_big = X_big @ rng.normal(size=k_big) + rng.normal(0, 0.1, n_big)
batch_sizes = [1, 2, 5, 10, 20, 50, 100, 200, 500, 1000]
seconds = []
for mu in batch_sizes:
theta, A_inv = np.zeros(k_big), np.eye(k_big)
start_time = time.perf_counter() # a precise clock for measuring durations
for start in range(0, n_big, mu):
X_mu, y_mu = X_big[start:start + mu], y_big[start:start + mu]
theta, A_inv = rls_batch_update(theta, A_inv, X_mu, y_mu, np.ones(len(X_mu)))
seconds.append(time.perf_counter() - start_time)
fig, ax = plt.subplots(figsize=(6.5, 3.6))
ax.loglog(batch_sizes, seconds, "o-", color=BLUE, markersize=5)
ax.set(xlabel=r"batch size $\mu$", ylabel="time in seconds", title="Processing 20,000 observations (k = 10)")
fig.tight_layout()
plt.show()
fastest = batch_sizes[int(np.argmin(seconds))]
print(f"fastest batch size in this run: {fastest} ({min(seconds):.3f} s, single observations: {seconds[0]:.3f} s)")
fastest batch size in this run: 50 (0.114 s, single observations: 2.901 s)
Moderate batch sizes are fastest. For very large batches, the $\mu\times\mu$ matrix dominates and the recursion becomes slower again.
8. Several Outputs at Once¶
Sometimes we want to predict several quantities from the same inputs, for example the horizontal and the vertical position of an object. In the blog post, this only requires replacing the targets by a matrix $\mathbf{Y}_\mu$ and the coefficient vector by a matrix $\boldsymbol{\Theta}_n$, each with one column per output. The gain and the update of $\mathbf{A}^{-1}$ do not depend on the targets at all. All outputs share them, so additional outputs reuse these computations. A single-observation update costs $O(k^2+km)$, and the coefficient matrix needs $O(km)$ storage; the additional cost is small when the number of outputs is small compared with the number of features.
For a single new observation, $\mathbf{Y}_\mu$ is a row vector with one target per output, and the update becomes
$$ \begin{aligned} \mathbf{E}_\mu &= \mathbf{Y}_\mu-\mathbf{x}_{n+1}^T\boldsymbol{\Theta}_n, \\ \boldsymbol{\Theta}_{n+1} &= \boldsymbol{\Theta}_n+\boldsymbol{\delta}_{n+1}\mathbf{E}_\mu, \end{aligned} $$
while $\boldsymbol{\delta}_{n+1}$ and $\mathbf{A}_{n+1}^{-1}$ are computed exactly as in Section 5. The row vector $\mathbf{E}_\mu$ contains the prediction errors of all outputs, and the coefficient update is an outer product: every column of $\boldsymbol{\Theta}$ is corrected in the same direction $\boldsymbol{\delta}_{n+1}$, scaled by the prediction error of its own output. In the code, only these two lines change:
def rls_update_multi(Theta, A_inv, x, y, w=1.0, lam_n=1.0):
"""Performs one recursive least squares update for m outputs.
Args:
Theta: Current coefficients Theta_n, one column per output, shape (k, m).
A_inv: Current inverse A_n^{-1}, shape (k, k).
x: New input vector x_{n+1}, shape (k,).
y: New targets, one per output, shape (m,).
w: Positive weight w_{n+1} of the new observation.
lam_n: Forgetting factor lambda_n of this update, in (0, 1].
Returns:
A tuple (Theta, A_inv) with the new coefficients Theta_{n+1}, shape
(k, m), and the new inverse A_{n+1}^{-1}, shape (k, k).
"""
e = y - x @ Theta # prediction errors of all outputs, shape (m,)
delta = w * A_inv @ x / (lam_n + w * x @ A_inv @ x) # the same gain as before, shape (k,)
Theta = Theta + np.outer(delta, e) # Theta_{n+1} = Theta_n + delta_{n+1} E_mu, shape (k, m)
A_inv = (A_inv - np.outer(delta, x @ A_inv)) / lam_n # the same update as before
return Theta, A_inv
For batches, the blog post uses $\mathbf{E}_\mu=\mathbf{Y}_\mu-\mathbf{X}_\mu\boldsymbol{\Theta}_n$ and $\boldsymbol{\Theta}_{n+\mu}=\boldsymbol{\Theta}_n+\boldsymbol{\Delta}_\mu\mathbf{E}_\mu$ with a $\mu\times m$ matrix $\mathbf{Y}_\mu$. Here, no new code is needed: rls_batch_update from Section 7 consists of matrix products only and works unchanged when theta and y_mu are matrices.
As an example, we learn a closed curve from noisy points. Each observation consists of a parameter $t$ (think of the time on a circular track) and the noisy position $(u, v)$ of the corresponding point on the curve. The inputs are $\mathbf{x}=(1,\cos t,\sin t,\cos 2t,\sin 2t,\ldots,\cos 4t,\sin 4t)^T$, so that $k=9$, and the two outputs are $u$ and $v$. Although the curve is far from a straight line, the model is linear in its coefficients, and that is all that least squares needs. The observations arrive one at a time, and we process them with rls_update_multi.
def fourier_features(t, order=4):
"""Computes the input vectors (1, cos t, sin t, ..., cos(order t), sin(order t)).
Args:
t: Curve parameters, shape (n,).
order: Highest multiple of t in the sines and cosines.
Returns:
The input vectors as rows, shape (n, 2 * order + 1).
"""
columns = [np.ones_like(t)]
for j in range(1, order + 1):
columns += [np.cos(j * t), np.sin(j * t)]
return np.column_stack(columns)
def heart(t):
"""Computes points of a heart-shaped closed curve.
Args:
t: Curve parameters in [0, 2 pi), shape (n,).
Returns:
The points (u, v) as rows, shape (n, 2).
"""
u = 12 * np.sin(t) - 4 * np.sin(3 * t)
v = 13 * np.cos(t) - 5 * np.cos(2 * t) - 2 * np.cos(3 * t) - np.cos(4 * t)
return np.column_stack([u, v])
rng = np.random.default_rng(6)
n_curve = 500
t_obs = rng.uniform(0, 2 * np.pi, n_curve)
X_curve = fourier_features(t_obs) # k = 9 inputs, shape (500, 9)
Y_curve = heart(t_obs) + rng.normal(0, 1.0, (n_curve, 2)) # m = 2 noisy outputs, shape (500, 2)
Theta, A_inv = np.zeros((9, 2)), np.eye(9) / rho # one column of coefficients per output
snapshots = {} # the coefficients after each update
for n, (x_i, y_i) in enumerate(zip(X_curve, Y_curve), start=1):
Theta, A_inv = rls_update_multi(Theta, A_inv, x_i, y_i) # x_i has shape (9,), y_i = (u, v) shape (2,)
snapshots[n] = Theta.copy() # copy, so that later updates do not change it
t_grid = np.linspace(0, 2 * np.pi, 400) # parameters for drawing smooth curves
fig, axes = plt.subplots(1, 4, figsize=(12, 3.6), sharex=True, sharey=True)
for ax, n in zip(axes, [10, 25, 60, 500]):
fitted = fourier_features(t_grid) @ snapshots[n] # (400, 9) @ (9, 2): 400 fitted points (u, v)
ax.scatter(*Y_curve[:n].T, s=10, color=POINTS, label="observations so far") # *...T unpacks u and v
ax.plot(*heart(t_grid).T, "--", color=INK, lw=1.1, label="true curve")
ax.plot(*fitted.T, color=BLUE, label="RLS estimate")
ax.set(title=f"after {n} observations", aspect="equal", xticks=[], yticks=[])
handles, labels = axes[0].get_legend_handles_labels()
fig.legend(handles, labels, loc="lower center", ncol=3)
fig.tight_layout(rect=(0, 0.08, 1, 1))
plt.show()
The learned coefficient matrix $\boldsymbol{\Theta}_{500}$ has one row per input ($1,\cos t,\sin t,\ldots,\sin 4t$) and one column per output ($u$ and $v$). Comparing it with the coefficients of the heart's formula shows that the recursion has found the curve's building blocks, apart from small deviations caused by the noise:
Theta_true = np.zeros((9, 2)) # the coefficients in the formula of heart()
Theta_true[[2, 6], 0] = [12, -4] # u = 12 sin t - 4 sin 3t
Theta_true[[1, 3, 5, 7], 1] = [13, -5, -2, -1] # v = 13 cos t - 5 cos 2t - 2 cos 3t - cos 4t
show((r"\boldsymbol{\Theta}_{500}", Theta), (r"\boldsymbol{\Theta}_{\text{true}}", Theta_true), digits=2)
Check: several outputs are just several single-output problems. Because the outputs share the inputs, the joint fit must agree with $m$ separate fits, one per output. It must also agree with the batch version from Section 7, which processes the same data in batches of five, and with the direct solution:
separate = []
for j in range(2): # fit each output on its own
theta_j, A_inv_j = np.zeros(9), np.eye(9) / rho
for x_i, y_i in zip(X_curve, Y_curve[:, j]):
theta_j, A_inv_j = rls_update(theta_j, A_inv_j, x_i, y_i)
separate.append(theta_j)
Theta_batch, A_inv_batch = np.zeros((9, 2)), np.eye(9) / rho # the same data in batches of 5
for start in range(0, n_curve, 5):
Theta_batch, A_inv_batch = rls_batch_update(Theta_batch, A_inv_batch, X_curve[start:start + 5],
Y_curve[start:start + 5], np.ones(5))
print(f"joint fit vs. separate fits: {np.abs(Theta - np.column_stack(separate)).max():.1e}")
print(f"single observations vs. batches: {np.abs(Theta - Theta_batch).max():.1e}")
print(f"joint fit vs. direct solution: {np.abs(Theta - batch_solution(X_curve, Y_curve, rho=rho)).max():.1e}")
joint fit vs. separate fits: 0.0e+00 single observations vs. batches: 4.2e-13 joint fit vs. direct solution: 5.3e-15
9. Putting It All Together: the WeightedRLS Class¶
We now have all ingredients, and the class below combines them:
- the initialization $\boldsymbol{\Theta}_0=\mathbf{0}$ and $\mathbf{A}_0^{-1}=\rho^{-1}\mathbf{I}$ (Section 3),
- positive observation weights (Section 4),
- exponential forgetting with $\lambda_0=1$ for the first update (Section 5),
- batches of any size, including single observations (Section 7), and
- one or several outputs (Section 8).
Compared with rls_batch_update, two implementation details differ. First, the product $\mathbf{X}_\mu\mathbf{A}_n^{-1}$ is computed once and reused, and the gain is obtained with np.linalg.solve instead of an explicit matrix inverse (the result is the same, but solving is cheaper and more accurate). Since this relies on the symmetry of $\mathbf{A}^{-1}$, the class symmetrizes it after every update, as discussed in Section 6. Second, the class checks its inputs, so that a wrongly shaped array or a non-positive weight raises an error instead of silently corrupting the estimate. The same code is available as the module weighted_rls.py.
Input-validation limits. Supply a positive integer batch_size, a finite positive ridge, nonempty batches for update, and exactly matching numbers of input rows, targets and weights. These conditions are not all enforced, and fit can fail after earlier batches have already changed the model. The module README documents the current edge cases.
from typing import Optional, Tuple
import numpy as np
from numpy.typing import ArrayLike, NDArray
class WeightedRLS:
"""Weighted recursive least squares with exponential forgetting.
Fits the linear model ``Y ≈ X @ theta`` and refines the coefficients
whenever new observations arrive, using the batch equations of the blog
post. Each call of :meth:`update` processes a batch of ``mu >= 1``
observations (the rows of ``X``) with optional positive weights, for one or
several outputs (the columns of ``Y``).
Attributes:
n_features: Number of entries k of an input vector.
n_outputs: Number of outputs m.
forgetting: Forgetting factor lambda in (0, 1].
ridge: Initial regularization coefficient rho > 0.
theta: Current coefficients, shape (n_features, n_outputs).
A_inv: Current inverse of the regularized, weighted normal matrix,
shape (n_features, n_features).
n_seen: Number of observations processed so far.
Example:
>>> rng = np.random.default_rng(0)
>>> X = np.column_stack([np.ones(100), rng.uniform(-1, 1, 100)])
>>> y = X @ [2.0, -1.0] + 0.01 * rng.normal(size=100)
>>> rls = WeightedRLS(n_features=2, ridge=1e-3)
>>> errors = rls.fit(X, y, batch_size=10)
>>> rls.theta.round(2).ravel()
array([ 2., -1.])
>>> rls.predict([1.0, 0.5]).round(2)
array([[1.5]])
"""
def __init__(self, n_features: int, n_outputs: int = 1,
forgetting: float = 1.0, ridge: float = 1.0) -> None:
"""Initializes the estimator with theta = 0 and A_inv = I / ridge.
Args:
n_features: Number of entries k of an input vector. If the model
needs an intercept, add a column of ones to the inputs.
n_outputs: Number of outputs m.
forgetting: Forgetting factor lambda in (0, 1]. With 1, all
observations keep their weight. With a smaller value, the
weights of all earlier batches are multiplied by lambda in
every update except the first one.
ridge: Initial regularization coefficient rho > 0, i.e.
A_0 = rho * I. With forgetting, it decays together with the
old observations.
Raises:
ValueError: If an argument is outside its valid range.
"""
if n_features < 1 or n_outputs < 1:
raise ValueError("n_features and n_outputs must be at least 1")
if not 0.0 < forgetting <= 1.0:
raise ValueError("forgetting must be in (0, 1]")
if not ridge > 0.0:
raise ValueError("ridge must be positive")
self.n_features = n_features
self.n_outputs = n_outputs
self.forgetting = forgetting
self.ridge = ridge
self.reset()
def reset(self) -> None:
"""Forgets all observations: sets theta = 0 and A_inv = I / ridge."""
self.theta: NDArray[np.float64] = np.zeros((self.n_features, self.n_outputs))
self.A_inv: NDArray[np.float64] = np.eye(self.n_features) / self.ridge
self.n_seen = 0
def update(self, X: ArrayLike, Y: ArrayLike,
weights: Optional[ArrayLike] = None) -> NDArray[np.float64]:
"""Learns from one batch of observations.
Args:
X: Input vectors as rows, shape (mu, n_features). A 1-D array is
interpreted as a single observation.
Y: Targets, shape (mu, n_outputs). For a single output, shape (mu,)
is accepted as well, and for a single observation, shape
(n_outputs,) or a scalar.
weights: Positive observation weights, shape (mu,), i.e. the
diagonal of W_mu. Defaults to all ones.
Returns:
The prediction errors ``Y - X @ theta``, shape (mu, n_outputs),
computed with the coefficients from *before* the update.
Raises:
ValueError: If the batch has the wrong shape, contains NaN or
infinite values, or has non-positive weights.
"""
X_mu, Y_mu, w_mu = self._check_batch(X, Y, weights)
lam_n = 1.0 if self.n_seen == 0 else self.forgetting # lambda_0 = 1: nothing to forget yet
E_mu = Y_mu - X_mu @ self.theta # prediction errors E_mu, shape (mu, m)
XA = X_mu @ self.A_inv # X_mu A_n^{-1}, shape (mu, k), used three times
S = lam_n * np.diag(1.0 / w_mu) + XA @ X_mu.T # lambda_n W_mu^{-1} + X_mu A_n^{-1} X_mu^T, (mu, mu)
# solve(S, XA) computes S^{-1} X_mu A_n^{-1} without inverting S explicitly. Because S
# and A^{-1} are symmetric, its transpose is the gain Delta_mu = A_n^{-1} X_mu^T S^{-1}.
Delta_mu = np.linalg.solve(S, XA).T # gain matrix, shape (k, mu)
self.theta = self.theta + Delta_mu @ E_mu # Theta_{n+mu} = Theta_n + Delta_mu E_mu
A_inv = (self.A_inv - Delta_mu @ XA) / lam_n # A_{n+mu}^{-1}
self.A_inv = 0.5 * (A_inv + A_inv.T) # keep A^{-1} exactly symmetric
self.n_seen += len(X_mu)
return E_mu
def fit(self, X: ArrayLike, Y: ArrayLike, batch_size: int = 1,
weights: Optional[ArrayLike] = None) -> NDArray[np.float64]:
"""Processes a whole data set in consecutive batches.
The fit continues from the current state; call :meth:`reset` first to
start from scratch.
Args:
X: Input vectors as rows, shape (n, n_features).
Y: Targets, shape (n, n_outputs), or (n,) for a single output.
batch_size: Number of rows per batch. The last batch may be
smaller.
weights: Positive observation weights, shape (n,). Defaults to
all ones.
Returns:
The prediction errors of all rows, shape (n, n_outputs), each
computed before its batch was learned.
"""
X = np.asarray(X, dtype=float)
Y = np.asarray(Y, dtype=float)
w = np.ones(len(X)) if weights is None else np.asarray(weights, dtype=float)
# rows i, ..., i + batch_size - 1 form one batch; slicing stops at the end of the data
errors = [self.update(X[i:i + batch_size], Y[i:i + batch_size], w[i:i + batch_size])
for i in range(0, len(X), batch_size)]
return np.vstack(errors) if errors else np.empty((0, self.n_outputs))
def predict(self, X: ArrayLike) -> NDArray[np.float64]:
"""Predicts the targets of the given inputs with the current coefficients.
Args:
X: Input vectors as rows, shape (n, n_features), or a single input
vector of shape (n_features,).
Returns:
The predictions ``X @ theta``, shape (n, n_outputs).
Raises:
ValueError: If X does not have n_features columns.
"""
X = np.atleast_2d(np.asarray(X, dtype=float)) # a single vector becomes one row
if X.shape[1] != self.n_features:
raise ValueError(f"X must have {self.n_features} columns")
return X @ self.theta
def _check_batch(self, X: ArrayLike, Y: ArrayLike, weights: Optional[ArrayLike]
) -> Tuple[NDArray[np.float64], NDArray[np.float64], NDArray[np.float64]]:
"""Brings a batch into standard shapes and validates it.
Args:
X: Input vectors, as described in :meth:`update`.
Y: Targets, as described in :meth:`update`.
weights: Observation weights or None, as described in :meth:`update`.
Returns:
A tuple ``(X_mu, Y_mu, w_mu)`` with the shapes (mu, k), (mu, m)
and (mu,).
Raises:
ValueError: If the batch has the wrong shape, contains NaN or
infinite values, or has non-positive weights.
"""
X_mu = np.atleast_2d(np.asarray(X, dtype=float)) # a single vector becomes one row
mu = X_mu.shape[0]
Y_mu = np.asarray(Y, dtype=float)
if Y_mu.ndim == 0 or (Y_mu.ndim == 1 and self.n_outputs == 1):
Y_mu = Y_mu.reshape(-1, 1) # scalar, or the mu targets of a single output
elif Y_mu.ndim == 1 and mu == 1:
Y_mu = Y_mu.reshape(1, -1) # the m targets of a single observation
w_mu = np.ones(mu) if weights is None else np.asarray(weights, dtype=float).reshape(-1)
if X_mu.ndim != 2 or X_mu.shape[1] != self.n_features:
raise ValueError(f"X must have shape (mu, {self.n_features})")
if Y_mu.shape != (mu, self.n_outputs):
raise ValueError(f"Y must have shape ({mu}, {self.n_outputs})")
if w_mu.shape != (mu,) or np.any(w_mu <= 0):
raise ValueError(f"weights must contain {mu} positive values")
if not (np.isfinite(X_mu).all() and np.isfinite(Y_mu).all() and np.isfinite(w_mu).all()):
raise ValueError("the batch must not contain NaN or infinite values")
return X_mu, Y_mu, w_mu
A short usage example, the same as in the docstring: we fit a line to $100$ observations in batches of ten and predict the value at $x=0.5$.
rng = np.random.default_rng(0)
X_demo = np.column_stack([np.ones(100), rng.uniform(-1, 1, 100)]) # column of ones for the intercept
y_demo = X_demo @ [2.0, -1.0] + 0.01 * rng.normal(size=100)
rls = WeightedRLS(n_features=2, ridge=1e-3)
errors = rls.fit(X_demo, y_demo, batch_size=10) # 10 batches of 10 observations
print("shape of the returned prediction errors:", errors.shape)
show((r"\boldsymbol{\theta}", rls.theta), (r"\mathbf{A}^{-1}", rls.A_inv))
show((r"\hat{y}(x=0.5)", rls.predict([1.0, 0.5]).item())) # predict expects the input vector (1, x)
shape of the returned prediction errors: (100, 1)
Final check. To make sure that all pieces fit together, we compare the class with the direct solution in several configurations: single observations and batches of ten (the last batch is smaller), one and three outputs, unit and random weights, with and without forgetting. With forgetting, the direct solution uses the weights $\lambda^a$, where the age $a$ is now counted in batches, and the decayed regularization $\rho\lambda^{M-1}$ after $M$ batches. The comparison is made after every update.
def direct_solution_with_forgetting(X, Y, w, rho, lam, batch_size):
"""Computes the direct solution of the problem that the recursion solves in batches.
An observation whose batch is a batches old has the weight w * lam^a, and
the regularization has decayed to rho * lam^(M - 1) after M batches.
Args:
X: Input vectors as rows, shape (n, k).
Y: Targets, shape (n, m).
w: Observation weights, shape (n,).
rho: Initial regularization coefficient rho > 0.
lam: Forgetting factor lambda in (0, 1].
batch_size: Number of rows per batch (the last batch may be smaller).
Returns:
The coefficients, shape (k, m).
"""
batch_index = np.arange(len(X)) // batch_size # 0, 0, ..., 1, 1, ...: the batch of each row
n_batches = batch_index[-1] + 1
ages = n_batches - 1 - batch_index # ages counted in batches
return batch_solution(X, Y, w * lam ** ages, rho * lam ** (n_batches - 1))
rng = np.random.default_rng(7)
X_test = np.column_stack([np.ones(103), rng.normal(size=(103, 3))]) # k = 4, 103 rows
Y_test = X_test @ rng.normal(size=(4, 3)) + 0.1 * rng.normal(size=(103, 3))
print("batch size | outputs | weights | lambda | largest difference")
for batch_size in [1, 10]:
for m in [1, 3]:
for weighted in [False, True]:
for lam in [1.0, 0.9]:
w = rng.uniform(0.2, 5, 103) if weighted else np.ones(103)
model = WeightedRLS(n_features=4, n_outputs=m, forgetting=lam, ridge=0.5)
worst = 0.0
for start in range(0, 103, batch_size):
stop = min(start + batch_size, 103)
model.update(X_test[start:stop], Y_test[start:stop, :m], w[start:stop])
direct = direct_solution_with_forgetting(X_test[:stop], Y_test[:stop, :m],
w[:stop], 0.5, lam, batch_size)
worst = max(worst, np.abs(model.theta - direct).max())
print(f"{batch_size:>10} | {m:>7} | {'random' if weighted else 'unit':>7} | {lam:>6} | {worst:.1e}")
batch size | outputs | weights | lambda | largest difference
1 | 1 | unit | 1.0 | 1.1e-15
1 | 1 | unit | 0.9 | 1.1e-15
1 | 1 | random | 1.0 | 8.9e-16
1 | 1 | random | 0.9 | 1.1e-15
1 | 3 | unit | 1.0 | 1.6e-15
1 | 3 | unit | 0.9 | 1.8e-15
1 | 3 | random | 1.0 | 2.7e-15
1 | 3 | random | 0.9 | 2.1e-15
10 | 1 | unit | 1.0 | 8.9e-16
10 | 1 | unit | 0.9 | 8.9e-16
10 | 1 | random | 1.0 | 1.0e-15
10 | 1 | random | 0.9 | 9.4e-16
10 | 3 | unit | 1.0 | 1.3e-15
10 | 3 | unit | 0.9 | 1.3e-15
10 | 3 | random | 1.0 | 1.3e-15
10 | 3 | random | 0.9 | 1.3e-15
10. Usage Examples¶
Finally, a few typical ways to use the class. The four examples cover single observations and batches, each with one and with several outputs. All of them use synthetic data, so the settings can easily be changed and the cells run again.
The Most Common Case: Single Observations, Slight Forgetting, No Weights¶
In many applications, the observations arrive one at a time, they are all equally reliable, and the underlying relationship changes slowly. Then we process one observation per update, keep the default weights and choose a forgetting factor slightly below $1$. A classic example is the prediction of a signal from its own past values: we predict the next value $y_t$ from the two previous values, i.e. with the input vector $\mathbf{x}_t=(y_{t-1},y_{t-2})^T$ and the model $y_t\approx\theta_1y_{t-1}+\theta_2y_{t-2}$. Our synthetic signal is a noisy oscillation whose frequency slowly increases, so the best coefficients change over time. In every step, the model first predicts the new value and then learns from it.
rng = np.random.default_rng(10)
n_signal = 3000
frequency = np.linspace(0.2, 1.0, n_signal) # slowly increasing frequency (radians per step)
theta1_true = 2 * 0.95 * np.cos(frequency) # the true coefficients of the signal
theta2_true = -0.95**2
signal = np.zeros(n_signal)
for t in range(2, n_signal):
signal[t] = theta1_true[t] * signal[t - 1] + theta2_true * signal[t - 2] + rng.normal()
results = {}
for lam in [1.0, 0.99]:
model = WeightedRLS(n_features=2, forgetting=lam, ridge=1.0)
predictions, coefficients = [], []
for t in range(2, n_signal):
x_t = [signal[t - 1], signal[t - 2]] # the two previous values
predictions.append(model.predict(x_t).item()) # predict first ...
model.update(x_t, signal[t]) # ... then learn from the new value
coefficients.append(model.theta[:, 0].copy())
results[lam] = (np.array(predictions), np.array(coefficients))
show((r"\boldsymbol{\theta}_{\text{end}}\ (\lambda=1)", results[1.0][1][-1]),
(r"\boldsymbol{\theta}_{\text{end}}\ (\lambda=0.99)", results[0.99][1][-1]),
(r"\boldsymbol{\theta}_{\text{true}}", [theta1_true[-1], theta2_true]))
steps = np.arange(2, n_signal)
fig, axes = plt.subplots(1, 3, figsize=(12, 3.9))
# Left: the last 80 values of the signal and their one-step predictions
last = slice(len(steps) - 80, len(steps))
axes[0].plot(steps[last], signal[2:][last], "o-", color=POINTS, markersize=3, lw=1, label="signal")
axes[0].plot(steps[last], results[0.99][0][last], color=ORANGE, lw=1.4, label=r"prediction, $\lambda$ = 0.99")
axes[0].set(xlabel="time step $t$", ylabel="$y_t$", ylim=(-7, 9), title="The last 80 steps")
axes[0].legend(loc="upper left", ncol=2, fontsize=8)
# Middle and right: the first coefficient and the moving average of the squared prediction errors
for lam, color in [(1.0, BLUE), (0.99, ORANGE)]:
predictions, coefficients = results[lam]
axes[1].plot(steps, coefficients[:, 0], color=color, lw=1.4, label=rf"$\lambda$ = {lam:g}")
squared = (signal[2:] - predictions) ** 2
moving_average = np.convolve(squared, np.ones(100) / 100, mode="valid") # mean of 100 consecutive values
axes[2].plot(steps[99:], moving_average, color=color, lw=1.4, label=rf"$\lambda$ = {lam:g}")
axes[1].plot(steps, theta1_true[2:], "--", color=INK, lw=1.2, label="true value")
axes[1].set(xlabel="time step $t$", title=r"Estimated coefficient $\theta_1$")
axes[1].legend()
axes[2].axhline(1.0, ls="--", color=INK, lw=1.2, label="noise variance")
axes[2].set(xlabel="time step $t$", ylabel="mean of the last 100 steps", title="Squared prediction error")
axes[2].legend()
fig.tight_layout()
plt.show()
for lam, (predictions, _) in results.items():
print(f"lambda = {lam:g}: mean squared prediction error over the last 500 steps "
f"{np.mean((signal[-500:] - predictions[-500:]) ** 2):.2f}")
lambda = 1: mean squared prediction error over the last 500 steps 2.90 lambda = 0.99: mean squared prediction error over the last 500 steps 1.02
Without forgetting, the estimated coefficient lags further and further behind, because the model tries to describe the whole history with a single set of coefficients, and its prediction error grows. With $\lambda=0.99$, which corresponds to a memory of about $100$ steps, the coefficient follows the slow change, and the prediction error stays close to the variance of the noise, which is the best that any predictor can achieve for this signal.
Single Observations and Several Outputs: Predicting the Next Position of a Moving Point¶
With several outputs, the typical setup is the same as in the previous example, except that each observation has several targets. Here, a point moves around the origin on a noisy circular path, and its angular speed increases slowly. Apart from noise, each new position is the previous position rotated by a small angle and pulled slightly towards the origin. We therefore predict both coordinates of the next position (two outputs) from the current position (two inputs). The learned $2\times2$ coefficient matrix describes the rotation, so we can read off the estimated rotation angle per step and compare it with the true one.
def rotation(angle):
"""Computes the matrix that rotates a point in the plane.
Args:
angle: Rotation angle in radians (counterclockwise).
Returns:
The rotation matrix, shape (2, 2).
"""
return np.array([[np.cos(angle), -np.sin(angle)], [np.sin(angle), np.cos(angle)]])
rng = np.random.default_rng(11)
n_steps = 3000
angle_true = np.linspace(0.1, 0.5, n_steps) # slowly increasing rotation per step
position = np.zeros((n_steps, 2))
for t in range(1, n_steps):
position[t] = 0.97 * rotation(angle_true[t]) @ position[t - 1] + rng.normal(0, 0.5, 2)
results_2d, final_theta = {}, {}
for lam in [1.0, 0.99]:
model = WeightedRLS(n_features=2, n_outputs=2, forgetting=lam, ridge=1.0)
predictions, angles = [], []
for t in range(1, n_steps):
predictions.append(model.predict(position[t - 1])[0]) # predict the next position ...
model.update(position[t - 1], position[t]) # ... then learn from it
Th = model.theta # 2 x 2 coefficient matrix
# For Th = r * rotation(a).T, the rotation angle a follows from the four entries:
angles.append(np.arctan2(Th[0, 1] - Th[1, 0], Th[0, 0] + Th[1, 1]))
results_2d[lam] = (np.array(predictions), np.array(angles))
final_theta[lam] = model.theta.copy()
# The model predicts the next position as a row vector: position[t] ≈ position[t-1] @ Theta,
# so the true coefficient matrix is the transpose of 0.97 * rotation(angle)
show((r"\boldsymbol{\Theta}_{\text{end}}\ (\lambda=0.99)", final_theta[0.99]),
(r"\boldsymbol{\Theta}_{\text{true}}", (0.97 * rotation(angle_true[-1])).T))
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(11, 4.4))
# Left: the last 30 positions, their predictions and a line from each prediction to its position
actual = position[1:][-30:] # the last 30 positions ...
predicted = results_2d[0.99][0][-30:] # ... and their one-step predictions
ax1.plot(*actual.T, color=POINTS, lw=1, zorder=1)
for a, p in zip(actual, predicted): # connect each prediction with its position
ax1.plot([p[0], a[0]], [p[1], a[1]], color=ORANGE, lw=0.8, alpha=0.6, zorder=2)
ax1.scatter(*actual.T, s=16, color=MUTED, zorder=3, label="positions")
ax1.scatter(*predicted.T, s=22, color=ORANGE, zorder=4, label=r"predictions, $\lambda$ = 0.99")
ax1.set(aspect="equal", xlabel="first coordinate", ylabel="second coordinate", title="The last 30 steps")
ax1.legend(fontsize=8, loc="upper left")
# Right: the estimated rotation angle per step
for lam, color in [(1.0, BLUE), (0.99, ORANGE)]:
ax2.plot(np.arange(1, n_steps), results_2d[lam][1], color=color, lw=1.4, label=rf"$\lambda$ = {lam:g}")
ax2.plot(angle_true, "--", color=INK, lw=1.2, label="true angle")
ax2.set(xlabel="time step $t$", ylabel="angle in radians", ylim=(0, 0.7), title="Estimated rotation per step")
ax2.legend()
fig.tight_layout()
plt.show()
for lam, (predictions, _) in results_2d.items():
squared_distance = np.sum((position[1:][-500:] - predictions[-500:]) ** 2, axis=1)
print(f"lambda = {lam:g}: mean squared prediction error over the last 500 steps {squared_distance.mean():.2f} "
f"(noise: {2 * 0.5**2:.2f})")
lambda = 1: mean squared prediction error over the last 500 steps 0.83 (noise: 0.50) lambda = 0.99: mean squared prediction error over the last 500 steps 0.54 (noise: 0.50)
As in the previous example, the model with $\lambda=0.99$ follows the slowly changing rotation, whereas the model without forgetting lags behind, and its predictions are worse. Both outputs share a single $2\times2$ matrix $\mathbf{A}^{-1}$.
A Whole Data Set in Batches: Fitting a Curve¶
RLS is not limited to straight lines. Any model that is linear in its coefficients works, for example a polynomial. If a data set is available at once (or arrives in chunks), we can process it in batches with fit. Here, we fit a cubic polynomial with the inputs $\mathbf{x}=(1,x,x^2,x^3)^T$ to noisy observations of $y=\sin(x)$, once with single observations and once in batches of $50$. Without forgetting, both give the same coefficients, but the batches are much faster.
def poly_features(x):
"""Computes the input vectors (1, x, x^2, x^3) of a cubic polynomial.
Args:
x: A number or an array of numbers, shape (n,).
Returns:
The input vectors as rows, shape (n, 4).
"""
x = np.atleast_1d(x) # a single number becomes an array of length 1
return np.column_stack([np.ones_like(x), x, x**2, x**3])
rng = np.random.default_rng(8)
x_sin = rng.uniform(-3, 3, 1000)
y_sin = np.sin(x_sin) + rng.normal(0, 0.3, 1000)
X_sin = poly_features(x_sin) # shape (1000, 4)
start_time = time.perf_counter()
rls_single = WeightedRLS(n_features=4, ridge=0.1)
rls_single.fit(X_sin, y_sin, batch_size=1) # 1000 updates with one observation each
seconds_single = time.perf_counter() - start_time
start_time = time.perf_counter()
rls_batch = WeightedRLS(n_features=4, ridge=0.1)
rls_batch.fit(X_sin, y_sin, batch_size=50) # 20 updates with 50 observations each
seconds_batch = time.perf_counter() - start_time
print(f"single observations: {seconds_single:.3f} s, batches of 50: {seconds_batch:.3f} s")
print(f"largest coefficient difference: {np.abs(rls_single.theta - rls_batch.theta).max():.1e}")
show((r"\boldsymbol{\theta}", rls_batch.theta)) # the four coefficients of 1, x, x^2 and x^3
fig, ax = plt.subplots(figsize=(6.5, 3.8))
x_grid = np.linspace(-3, 3, 200)
ax.scatter(x_sin, y_sin, s=6, color=POINTS, label="observations")
ax.plot(x_grid, rls_batch.predict(poly_features(x_grid))[:, 0], color=BLUE, label="cubic polynomial (RLS)")
ax.plot(x_grid, np.sin(x_grid), "--", color=INK, lw=1.2, label=r"true curve $\sin(x)$")
ax.set(xlabel="$x$", ylabel="$y$", ylim=(-2.5, 2.5), title="A cubic polynomial fitted in batches of 50")
ax.legend(loc="upper left", fontsize=8)
fig.tight_layout()
plt.show()
single observations: 0.146 s, batches of 50: 0.007 s largest coefficient difference: 1.2e-14
Batches, Several Outputs and Forgetting: Tracking a Changing Shape¶
In the last example, all features of the class come together. We observe noisy points of a closed curve in batches of $20$ and learn both coordinates at once, as in Section 8. After $1{,}500$ observations, the curve suddenly changes from the heart into an ellipse. We compare a model without forgetting with one that uses $\lambda=0.8$ per batch, which corresponds to a memory of about $1/(1-0.8)=5$ batches, or $100$ observations.
def ellipse(t):
"""Computes points of an ellipse.
Args:
t: Curve parameters in [0, 2 pi), shape (n,).
Returns:
The points (u, v) as rows, shape (n, 2).
"""
return np.column_stack([14 * np.sin(t), 10 * np.cos(t)])
rng = np.random.default_rng(9)
n_shape, mu, change = 3000, 20, 1500
t_shape = rng.uniform(0, 2 * np.pi, n_shape)
# the first 1500 points lie on the heart, the remaining ones on the ellipse
true_points = np.where(np.arange(n_shape)[:, None] < change, heart(t_shape), ellipse(t_shape))
Y_shape = true_points + rng.normal(0, 1.0, (n_shape, 2))
X_shape = fourier_features(t_shape)
models = {1.0: WeightedRLS(9, 2, forgetting=1.0, ridge=0.01),
0.8: WeightedRLS(9, 2, forgetting=0.8, ridge=0.01)}
snapshots = {lam: {} for lam in models}
for start in range(0, n_shape, mu): # batches of 20 observations
for lam, model in models.items():
model.update(X_shape[start:start + mu], Y_shape[start:start + mu])
if start + mu in (1500, 1700, 3000):
snapshots[lam][start + mu] = model.theta.copy()
fig, axes = plt.subplots(2, 3, figsize=(11, 7.6), sharex=True, sharey=True)
titles = {1500: "1,500 observations (before the change)", 1700: "1,700 observations", 3000: "3,000 observations"}
for row, (lam, color) in enumerate([(1.0, BLUE), (0.8, ORANGE)]): # one row per forgetting factor
for col, n in enumerate([1500, 1700, 3000]): # one column per point in time
ax = axes[row, col]
shape_now = heart if n <= change else ellipse
ax.scatter(*Y_shape[n - 100:n].T, s=8, color=POINTS, label="last 100 observations")
ax.plot(*shape_now(t_grid).T, "--", color=INK, lw=1.1, label="current true curve")
ax.plot(*(fourier_features(t_grid) @ snapshots[lam][n]).T, color=color, label=rf"estimate with $\lambda$ = {lam:g}")
ax.set(aspect="equal", xticks=[], yticks=[])
if row == 0:
ax.set_title(titles[n], fontsize=10)
if col == 0:
ax.set_ylabel(rf"$\lambda$ = {lam:g}")
handles, labels = axes[0, 0].get_legend_handles_labels()
handles.append(axes[1, 0].get_lines()[-1]); labels.append(axes[1, 0].get_lines()[-1].get_label())
fig.legend(handles, labels, loc="lower center", ncol=4)
fig.tight_layout(rect=(0, 0.05, 1, 1))
plt.show()
Without forgetting, the estimate after the change is a blend of both shapes, and even $1{,}500$ observations later it still carries the heart along. With $\lambda=0.8$ per batch, the estimate has adopted the ellipse after a few batches.
11. Summary¶
The core of the class consists of the batch equations of the blog post,
$$ \begin{aligned} \mathbf{E}_\mu &= \mathbf{Y}_\mu-\mathbf{X}_\mu\boldsymbol{\Theta}_n, \\ \boldsymbol{\Delta}_\mu &= \mathbf{A}_n^{-1}\mathbf{X}_\mu^T\big(\lambda_n\mathbf{W}_\mu^{-1}+\mathbf{X}_\mu\mathbf{A}_n^{-1}\mathbf{X}_\mu^T\big)^{-1}, \\ \boldsymbol{\Theta}_{n+\mu} &= \boldsymbol{\Theta}_n+\boldsymbol{\Delta}_\mu\mathbf{E}_\mu, \\ \mathbf{A}_{n+\mu}^{-1} &= \frac{1}{\lambda_n}\Big(\mathbf{A}_n^{-1}-\boldsymbol{\Delta}_\mu\mathbf{X}_\mu\mathbf{A}_n^{-1}\Big), \end{aligned} $$
with $\boldsymbol{\Theta}_0=\mathbf{0}$, $\mathbf{A}_0^{-1}=\rho^{-1}\mathbf{I}$ and $\lambda_0=1$. In WeightedRLS.update, they correspond to these lines:
lam_n = 1.0 if self.n_seen == 0 else self.forgetting # lambda_0 = 1: nothing to forget yet
E_mu = Y_mu - X_mu @ self.theta # prediction errors E_mu, shape (mu, m)
XA = X_mu @ self.A_inv # X_mu A_n^{-1}, shape (mu, k), used three times
S = lam_n * np.diag(1.0 / w_mu) + XA @ X_mu.T # lambda_n W_mu^{-1} + X_mu A_n^{-1} X_mu^T, (mu, mu)
# solve(S, XA) computes S^{-1} X_mu A_n^{-1} without inverting S explicitly. Because S
# and A^{-1} are symmetric, its transpose is the gain Delta_mu = A_n^{-1} X_mu^T S^{-1}.
Delta_mu = np.linalg.solve(S, XA).T # gain matrix, shape (k, mu)
self.theta = self.theta + Delta_mu @ E_mu # Theta_{n+mu} = Theta_n + Delta_mu E_mu
A_inv = (self.A_inv - Delta_mu @ XA) / lam_n # A_{n+mu}^{-1}
self.A_inv = 0.5 * (A_inv + A_inv.T) # keep A^{-1} exactly symmetric
A few practical recommendations follow from the examples:
- Add a column of ones to the inputs if the model needs an intercept.
- Choose $\rho$ small compared with the typical size of $\mathbf{x}^T\mathbf{x}$, but not extremely small (Section 3). For inputs of order one, values between $10^{-3}$ and $1$ work well.
- Use $\lambda=1$ if the relationship does not change. If it may change slowly, a forgetting factor slightly below $1$, such as $0.99$, is the usual choice; in general, choose $\lambda$ from the desired memory of about $1/(1-\lambda)$ updates (Sections 5 and 10). In batch mode, the memory is counted in batches.
- Observation weights are only needed if some observations are known to be more reliable than others (Section 4). Otherwise, the default weights of $1$ are the right choice.
- Keep $\mathbf{A}^{-1}$ symmetric, in particular with forgetting (Section 6).
- In double precision, the update of the class is accurate even for strongly correlated inputs. Only if you have to compute in single precision (
float32), use the QR-based update from the appendix. - Moderate batch sizes are usually fastest (Section 7).
- With forgetting, directions of the input space that receive no new information for a long time are forgotten as well, and the corresponding entries of $\mathbf{A}^{-1}$ grow. Make sure that the inputs keep varying, or choose $\lambda$ closer to $1$.
Appendix: A QR-based Update for Single Precision¶
This appendix is optional. Everything above uses double precision (float64, NumPy's default, with about $16$ significant digits), and in double precision the standard update of the class is accurate enough for the examples here, including the strongly correlated inputs below. On some hardware, however, for example on microcontrollers or GPUs, computations run in single precision (float32, about $7$ significant digits). There, the standard update can fail when inputs are strongly correlated, whereas a variant based on the QR decomposition keeps working. The QR variant can improve robustness in double precision as well; the difference is small in these examples, but can matter with more severe ill-conditioning, poor scaling or cancellation. Its benefit is therefore not limited to single precision.
What is a QR decomposition? Every matrix $\mathbf{S}$ with at least as many rows as columns can be written as the product $\mathbf{S}=\mathbf{Q}\mathbf{R}$ of two special matrices. The columns of $\mathbf{Q}$ have length $1$ and are perpendicular to each other, so that $\mathbf{Q}^T\mathbf{Q}=\mathbf{I}$. Multiplying a vector by such a matrix only rotates or reflects it, without changing its length, and therefore does not amplify rounding errors either. The matrix $\mathbf{R}$ is upper triangular, i.e. all its entries below the diagonal are zero, so that a system $\mathbf{R}\boldsymbol{\theta}=\mathbf{c}$ can be solved row by row, starting with the last row (back substitution). NumPy computes the decomposition with np.linalg.qr. A small example with a $3\times2$ matrix:
S_example = np.array([[1.0, 2.0], [1.0, 0.0], [1.0, -1.0]])
Q_example, R_example = np.linalg.qr(S_example) # S = Q R
show((r"\mathbf{S}", S_example), (r"\mathbf{Q}", Q_example), (r"\mathbf{R}", R_example))
show((r"\mathbf{Q}^T\mathbf{Q}", Q_example.T @ Q_example), (r"\mathbf{Q}\mathbf{R}", Q_example @ R_example))
Apart from rounding errors of about $10^{-16}$, $\mathbf{Q}^T\mathbf{Q}$ is the identity matrix and $\mathbf{Q}\mathbf{R}$ reproduces $\mathbf{S}$.
The idea. Instead of $\mathbf{A}_n^{-1}$, we keep an upper triangular matrix $\mathbf{R}_n$ with $\mathbf{R}_n^T\mathbf{R}_n=\mathbf{A}_n$, a kind of square root of $\mathbf{A}_n$, and the vector $\mathbf{c}_n=\mathbf{R}_n\boldsymbol{\theta}_n$. From the blog post, we know how $\mathbf{A}$ and $\mathbf{b}$ change with a new batch:
$$\mathbf{A}_{n+\mu}=\lambda_n\mathbf{A}_n+\mathbf{X}_\mu^T\mathbf{W}_\mu\mathbf{X}_\mu, \qquad \mathbf{b}_{n+\mu}=\lambda_n\mathbf{b}_n+\mathbf{X}_\mu^T\mathbf{W}_\mu\mathbf{y}_\mu .$$
Since $\mathbf{R}_n^T\mathbf{c}_n=\mathbf{R}_n^T\mathbf{R}_n\boldsymbol{\theta}_n=\mathbf{A}_n\boldsymbol{\theta}_n=\mathbf{b}_n$, both can be written with the stacked matrix $\mathbf{S}$ and the stacked vector $\mathbf{t}$,
$$\mathbf{A}_{n+\mu}=\mathbf{S}^T\mathbf{S}, \quad \mathbf{b}_{n+\mu}=\mathbf{S}^T\mathbf{t}, \qquad \mathbf{S}=\begin{bmatrix}\sqrt{\lambda_n}\,\mathbf{R}_n\\ \mathbf{W}_\mu^{1/2}\mathbf{X}_\mu\end{bmatrix}, \quad \mathbf{t}=\begin{bmatrix}\sqrt{\lambda_n}\,\mathbf{c}_n\\ \mathbf{W}_\mu^{1/2}\mathbf{y}_\mu\end{bmatrix}.$$
With the QR decomposition $\mathbf{S}=\mathbf{Q}\mathbf{R}_{n+\mu}$, the factor $\mathbf{Q}$ cancels because $\mathbf{Q}^T\mathbf{Q}=\mathbf{I}$: then $\mathbf{A}_{n+\mu}=\mathbf{R}_{n+\mu}^T\mathbf{R}_{n+\mu}$, and with $\mathbf{c}_{n+\mu}=\mathbf{Q}^T\mathbf{t}$, the normal equations $\mathbf{A}_{n+\mu}\boldsymbol{\theta}_{n+\mu}=\mathbf{b}_{n+\mu}$ simplify to the triangular system $\mathbf{R}_{n+\mu}\boldsymbol{\theta}_{n+\mu}=\mathbf{c}_{n+\mu}$. The recursion starts with $\mathbf{R}_0=\sqrt{\rho}\,\mathbf{I}$ and $\mathbf{c}_0=\mathbf{0}$.
This form is more robust because the condition number of $\mathbf{R}$ is only the square root of that of $\mathbf{A}$, and because orthogonal transformations do not amplify rounding errors. Neither $\mathbf{A}$ nor $\mathbf{A}^{-1}$ is ever formed. The price is a QR decomposition of a $(k+\mu)\times k$ matrix in every update, and the coefficients have to be computed from the triangular system whenever they are needed.
def rls_qr_update(R, c, X_mu, Y_mu, w_mu, lam_n=1.0):
"""Performs one QR-based recursive least squares update.
All computations use the floating-point precision of R (for example
float32), so the function can also run in single precision.
Args:
R: Current upper triangular factor R_n with R_n^T R_n = A_n, shape (k, k).
c: Current vector c_n = R_n theta_n, shape (k, m).
X_mu: Input vectors of the batch as rows, shape (mu, k).
Y_mu: Targets of the batch, shape (mu, m).
w_mu: Positive observation weights of the batch, shape (mu,).
lam_n: Forgetting factor lambda_n of this update, in (0, 1].
Returns:
A tuple (R, c) with the new factor R_{n+mu}, shape (k, k), and the new
vector c_{n+mu}, shape (k, m). The coefficients follow from R theta = c.
"""
X_mu, Y_mu, w_mu = (np.asarray(a, dtype=R.dtype) for a in (X_mu, Y_mu, w_mu)) # precision of R
s = np.sqrt(np.asarray(lam_n, dtype=R.dtype)) # sqrt(lambda_n)
sqrt_w = np.sqrt(w_mu)[:, None] # multiplying row i by sqrt(w_i) applies W_mu^{1/2}
S = np.vstack([s * R, sqrt_w * X_mu]) # the stacked matrix S, shape (k + mu, k)
t = np.vstack([s * c, sqrt_w * Y_mu]) # the stacked vector t, shape (k + mu, m)
Q, R_new = np.linalg.qr(S) # S = Q R_new, Q^T Q = I, R_new upper triangular
return R_new, Q.T @ t # R_{n+mu} and c_{n+mu} = Q^T t
Check: the QR recursion solves the same problem. In double precision, the QR recursion must also reproduce the direct solution after every update. We use the stream from Section 1 and treat each observation as a batch with $\mu=1$:
R, c = np.sqrt(rho) * np.eye(2), np.zeros((2, 1)) # R_0 = sqrt(rho) I and c_0 = 0
worst = 0.0
for n in range(n_obs):
R, c = rls_qr_update(R, c, X[n:n + 1], y[n:n + 1, None], np.ones(1)) # one observation, shapes (1, 2) and (1, 1)
theta_qr = np.linalg.solve(R, c)[:, 0] # solve R theta = c
worst = max(worst, np.abs(theta_qr - refits[n]).max())
print(f"largest difference to the direct solution: {worst:.1e}")
show((r"\mathbf{R}_{200}", R), (r"\mathbf{R}_{200}^T\mathbf{R}_{200}", R.T @ R),
(r"\mathbf{A}_{200}=\mathbf{X}_{200}^T\mathbf{X}_{200}+\rho\mathbf{I}", X.T @ X + rho * np.eye(2)), digits=2)
largest difference to the direct solution: 3.1e-15
Single versus double precision. Now we compare three variants on a stream of $20{,}000$ observations with the forgetting factor $\lambda=0.99$: the standard update of the class in double precision, the same update in single precision, and the QR update in single precision. The inputs are $(1, x, x+\varepsilon z)$, where $x$ and $z$ are random numbers. For small $\varepsilon$, the last two inputs are strongly correlated, since they differ only by the small term $\varepsilon z$. We measure the squared one-step prediction errors, which cannot fall below the variance of the noise, $10^{-4}$.
The class itself always computes in double precision, so we write its update for a single observation as a small function that keeps the precision of its inputs:
def standard_update(theta, A_inv, x, y, lam_n=1.0):
"""Performs the update of WeightedRLS for one observation with unit weight.
Unlike the class, the function keeps the floating-point precision of its
inputs, so it can also run in single precision.
Args:
theta: Current coefficients theta_n, shape (k,).
A_inv: Current inverse A_n^{-1}, shape (k, k).
x: New input vector x_{n+1}, shape (k,).
y: New target y_{n+1} (a number).
lam_n: Forgetting factor lambda_n of this update, in (0, 1].
Returns:
A tuple (theta, A_inv) with the new coefficients and the new inverse.
"""
e = y - x @ theta
Ax = A_inv @ x # reused, as in the class
delta = Ax / (lam_n + x @ Ax)
theta = theta + delta * e
A_inv = (A_inv - np.outer(delta, Ax)) / lam_n
return theta, 0.5 * (A_inv + A_inv.T) # symmetrized, as in the class
def prediction_errors(eps, n=20_000, lam=0.99, rho=0.01, seed=13):
"""Runs the three variants on the same stream and records their prediction errors.
Args:
eps: Difference between the third and the second input (small values
mean strongly correlated inputs).
n: Number of observations.
lam: Forgetting factor lambda.
rho: Initial regularization coefficient rho.
seed: Seed of the random number generator.
Returns:
A dict that maps the name of each variant to its squared one-step
prediction errors, shape (n,).
"""
rng = np.random.default_rng(seed)
x, z = rng.uniform(-1, 1, n), rng.normal(size=n)
X64 = np.column_stack([np.ones(n), x, x + eps * z])
y64 = X64 @ [0.5, 1.0, -1.0] + 0.01 * rng.normal(size=n)
X32, y32 = X64.astype(np.float32), y64.astype(np.float32) # the same data in single precision
theta64, A_inv64 = np.zeros(3), np.eye(3) / rho # standard update, float64
theta32 = np.zeros(3, dtype=np.float32) # standard update, float32
A_inv32 = np.eye(3, dtype=np.float32) / np.float32(rho)
R32 = np.sqrt(np.float32(rho)) * np.eye(3, dtype=np.float32) # QR update, float32
c32 = np.zeros((3, 1), dtype=np.float32)
one = np.ones(1, dtype=np.float32)
errors = {"standard, float64": [], "standard, float32": [], "QR, float32": []}
for t in range(n):
lam_n = 1.0 if t == 0 else lam
# predict first ...
errors["standard, float64"].append((y64[t] - X64[t] @ theta64) ** 2)
errors["standard, float32"].append((y64[t] - X32[t] @ theta32) ** 2)
errors["QR, float32"].append((y64[t] - X32[t] @ np.linalg.solve(R32, c32)[:, 0]) ** 2)
# ... then learn from the new observation
theta64, A_inv64 = standard_update(theta64, A_inv64, X64[t], y64[t], lam_n)
theta32, A_inv32 = standard_update(theta32, A_inv32, X32[t], y32[t], lam_n)
R32, c32 = rls_qr_update(R32, c32, X32[t:t + 1], y32[t:t + 1, None], one, lam_n)
return {name: np.array(values) for name, values in errors.items()}
comparison = {eps: prediction_errors(eps) for eps in [1e-2, 1e-4]}
fig, axes = plt.subplots(1, 2, figsize=(11, 4), sharey=True)
window = 500
# Lines that coincide stay visible: double precision is drawn wide and light, QR dashed on top
styles = {"standard, float64": dict(color=BLUE, lw=4, alpha=0.35),
"standard, float32": dict(color=ORANGE, lw=1.4),
"QR, float32": dict(color=AQUA, lw=1.4, ls=(0, (4, 2)))}
titles = {1e-2: r"Moderately correlated inputs ($\varepsilon$ = 0.01)", 1e-4: r"Strongly correlated inputs ($\varepsilon$ = 0.0001)"}
for ax, (eps, errors) in zip(axes, comparison.items()):
for name, values in errors.items():
moving_average = np.convolve(values, np.ones(window) / window, mode="valid") # mean of 500 consecutive values
ax.plot(np.arange(window, len(values) + 1), moving_average, label=name, **styles[name])
ax.axhline(1e-4, ls="--", color=INK, lw=1.2, label="noise variance")
ax.set(yscale="log", xlabel="time step", title=titles[eps])
axes[0].set_ylabel("squared prediction error (mean of 500)")
axes[0].legend(fontsize=8, loc="upper right")
fig.tight_layout()
plt.show()
for eps, errors in comparison.items():
summary = ", ".join(f"{name}: {np.mean(values[-5000:]):.1e}" for name, values in errors.items())
print(f"eps = {eps:g}, mean squared prediction error of the last 5000 steps: {summary}")
eps = 0.01, mean squared prediction error of the last 5000 steps: standard, float64: 1.0e-04, standard, float32: 1.0e-04, QR, float32: 1.0e-04 eps = 0.0001, mean squared prediction error of the last 5000 steps: standard, float64: 1.0e-04, standard, float32: 3.6e-02, QR, float32: 1.0e-04
With moderately correlated inputs (left), all three variants predict equally well. With strongly correlated inputs (right), the standard update in single precision breaks down: its prediction errors are often hundreds of times larger than the noise, and its $\mathbf{A}^{-1}$ even loses its positive definiteness at times. The QR update in single precision predicts as well as the standard update in double precision. In double precision, the standard update is accurate even for these inputs, which is why the class does not use the QR variant. If you have to compute in single precision, replace the update of the class by rls_qr_update and obtain the coefficients from $\mathbf{R}\boldsymbol{\theta}=\mathbf{c}$ whenever you need them.