---
title: "The Variance of Least Squares"
subtitle: "Math 123: Mathematical Aspects of Data Analysis — Fall 2026"
date: 2026-10-01
number-sections: false
jupyter: python3
format:
html:
toc: true
toc-depth: 2
code-line-numbers: true
code-tools: true
out-width: "100%"
column: page
execute:
echo: true
draft: false
---
In [Lecture 8](../../lectures/lecture-08.qmd#thm-least-squares-moments) we showed that, under the [standard assumptions](../../lectures/lecture-08.qmd#def-standard-assumptions), least squares is unbiased, $\Ebb[\widehat{\bftheta}_{\mathrm{LS}}]=\bftheta$.
We then used the SVD $\bfX=\sum_{i=1}^d\sigma_i\bfu_i\bfv_i^\top$ to see [where the least squares error lives](../../lectures/lecture-08.qmd#cor-least-squares-error-svd):
$$
\Var(\bfv_i^\top\widehat{\bftheta}_{\mathrm{LS}})=\frac{\sigma^2}{\sigma_i^2},
\qquad
\Ebb\big[\|\widehat{\bftheta}_{\mathrm{LS}}-\bftheta\|_2^2\big]=\sigma^2\sum_{i=1}^d\frac{1}{\sigma_i^2}.
$$
Here we check both facts by simulation, on a problem where most of the singular values are tiny: sharpening an out-of-focus photograph.
Least squares will be right *on average*, but a single estimate will not even be close.
## Blurring is a matrix
A greyscale photograph is an array of pixel brightnesses between $0$ (black) and $1$ (white), as in the [greyscale image example](../../lectures/lecture-03.qmd#exm-greyscale-image).
An out-of-focus camera does not record each pixel's own brightness.
The light meant for one pixel spreads onto its neighbours, so each recorded pixel is a weighted average of the true pixels around it.
Start with a single row of $d$ pixels, $\bftheta\in\Rbb^d$.
A *Gaussian blur* of width $w$ replaces pixel $j$ by
$$
(\bfX\bftheta)_j=\sum_{k=1}^d g(j-k)\,\theta_k,
\qquad
g(m)=\frac{1}{\sqrt{2\pi}\,w}e^{-m^2/(2w^2)}.
$$
This is a weighted average of the pixels near $j$, with bell-shaped weights that sum to about $1$ and are concentrated within $2w$ or so pixels of $j$.
It is linear in $\bftheta$, so it is multiplication by the $d\times d$ matrix $\bfX$ with entries $X_{jk}=g(j-k)$, whose row $j$ is a bump centred on the diagonal.
(Pixels beyond the ends of the row count as black.)
The camera's sensor then adds a little noise, so we observe
$$
\bfy=\bfX\bftheta+\bfvarepsilon,
\qquad
\varepsilon_j\ \text{independent},\ \mathcal N(0,\sigma^2).
$$
This is the linear model under the standard assumptions, with the blur as the design matrix and the sharp row as the unknown parameter.
*Deblurring* is the problem of estimating $\bftheta$ from $\bfy$.
Here $n=d$, and $\bfX$ will turn out to be invertible, so least squares fits the data exactly and simply undoes the blur: $\widehat{\bftheta}_{\mathrm{LS}}=\bfX^{-1}\bfy$.
```{python}
import base64
import io
import numpy as np
import plotly.graph_objects as go
from matplotlib import colormaps
from PIL import Image
from plotly.subplots import make_subplots
from skimage import data
rng = np.random.default_rng(123)
# The blur: row j of X is a Gaussian bump of width w centred on pixel j.
d = 128
w = 2.0
pixels = np.arange(d)
X = np.exp(-(pixels[:, None] - pixels[None, :])**2 / (2 * w**2)) / (np.sqrt(2 * np.pi) * w)
sigma = 0.01 # the sensor noise: 1% of the range from black to white
print("row 64 of X, near the diagonal =", np.round(X[64, 58:71], 3))
print(f"sum of row 64 = {X[64].sum():.6f}")
```
```{python}
#| code-fold: true
#| code-summary: "Figure style and image helpers"
# The same muted palette as the other examples, which reads correctly against
# both the light and the dark version of this site.
INK = "#7c7c78"
GRID = "rgba(128, 128, 128, 0.28)"
BLUE, ORANGE = "#2a78d6", "#eb6834"
GAP = 6 # transparent pixels between images placed side by side
def style(fig, height=420, title=None, showlegend=True):
"""Apply the shared layout to a plot with axes."""
fig.update_layout(
template="simple_white",
paper_bgcolor="rgba(0,0,0,0)",
plot_bgcolor="rgba(0,0,0,0)",
font=dict(color=INK, size=13),
title=dict(text=title, x=0, xanchor="left", font=dict(size=15)),
margin=dict(l=60, r=30, t=60 if title else 40, b=50),
height=height,
showlegend=showlegend,
hovermode="closest",
legend=dict(orientation="h", y=-0.18, x=0),
)
fig.update_xaxes(gridcolor=GRID, linecolor=GRID, zerolinecolor=GRID, showgrid=True)
fig.update_yaxes(gridcolor=GRID, linecolor=GRID, zerolinecolor=GRID, showgrid=True)
return fig
def png_uri(array):
"""Encode a greyscale array as a PNG data URI, clipping brightness to [0, 1].
NaN pixels become transparent. A PNG is far smaller in the page than a heatmap.
"""
rgba = (colormaps["gray"](np.clip(np.nan_to_num(array), 0, 1)) * 255).astype(np.uint8)
rgba[np.isnan(array), 3] = 0
buffer = io.BytesIO()
Image.fromarray(rgba, "RGBA").save(buffer, "PNG")
return "data:image/png;base64," + base64.b64encode(buffer.getvalue()).decode()
def side_by_side(images):
"""Place d x d images in a row, separated by transparent pixels."""
out = np.full((d, len(images) * (d + GAP) - GAP), np.nan)
for j, image in enumerate(images):
out[:, j * (d + GAP):j * (d + GAP) + d] = image
return out
def label_images(fig, labels):
"""Put one label under each image of a side_by_side picture and hide the axes."""
fig.update_xaxes(tickvals=[j * (d + GAP) + (d - 1) / 2 for j in range(len(labels))],
ticktext=labels, showgrid=False, zeroline=False, showline=False, ticks="")
fig.update_yaxes(showticklabels=False, showgrid=False, zeroline=False, showline=False,
ticks="")
return fig
```
### Blurring a photograph
To blur a whole photograph $\bfTheta\in\Rbb^{d\times d}$, blur every column and then every row: $\bfTheta\mapsto\bfX\bfTheta\bfX^\top$.
This is again linear in the $d^2$ pixels, so it fits the same model, and least squares undoes both blurs: $\widehat{\bfTheta}=\bfX^{-1}\bfY\bfX^{-\top}$.
```{python}
# scikit-image's "cameraman" photograph, averaged down from 512 x 512 to 128 x 128.
photo = data.camera().reshape(d, 4, d, 4).mean(axis=(1, 3)) / 255
# Blur every column, then every row, then add sensor noise to every pixel.
blurred = X @ photo @ X.T + sigma * rng.standard_normal((d, d))
# Least squares: solve X Z = Y to undo the column blur, then the same for the rows.
naive = np.linalg.solve(X, np.linalg.solve(X, blurred).T).T
print(f"brightness of the photo: {photo.min():.2f} to {photo.max():.2f}")
print(f"brightness of the blurred photo: {blurred.min():.2f} to {blurred.max():.2f}")
print(f"brightness of the least squares fit: {naive.min():.1e} to {naive.max():.1e}")
```
```{python}
#| label: fig-photos
#| fig-cap: "The sharp photograph, the blurred photograph with sensor noise of standard deviation 0.01, and the least squares estimate of the sharp photograph, with brightness clipped to the range from black to white."
#| fig-alt: "Three greyscale images side by side. The first is a sharp photograph of a man with a camera on a tripod. The second is the same photograph, blurred, with a dark frame around the edges. The third is black and white static with no trace of the photograph."
#| code-fold: true
#| code-summary: "Plot the photographs"
fig = go.Figure(go.Image(source=png_uri(side_by_side([photo, blurred, naive])),
hoverinfo="skip"))
style(fig, height=330, showlegend=False)
label_images(fig, ["sharp", "blurred, with noise", "least squares estimate"])
fig.update_layout(margin=dict(l=20, r=20, t=20, b=40))
fig.show()
```
The noise in the middle picture is invisible.
A standard deviation of $0.01$ is about $2.5$ of the $256$ grey levels in a typical image file, and just rounding to those $256$ levels would already add noise with standard deviation about $0.001$.
The blur also looks mild: the cameraman is still recognisable, and the dark frame is the black that the blur pulls in from beyond the edges.
Yet the least squares estimate is static, with "brightnesses" in the tens of trillions.
To see why, we go back to a single row, where everything is a vector in $\Rbb^{128}$.
## One row
Take row $50$ of the photograph, which runs across the cameraman's coat, his hands, and his camera.
```{python}
row = 50
theta = photo[row]
y = X @ theta + sigma * rng.standard_normal(d)
```
```{python}
#| label: fig-row
#| fig-cap: "Row 50 of the photograph: the true brightness theta (orange, dashed), the blurred row X theta (blue), and the observed data y, which adds the sensor noise (grey)."
#| fig-alt: "A line chart of brightness against pixel position from 0 to 127. A jagged orange dashed line with sharp jumps is the true row. A smoother blue line follows it but rounds off every jump. Grey dots lie almost exactly on the blue line."
#| code-fold: true
#| code-summary: "Plot the row"
fig = go.Figure()
fig.add_trace(go.Scatter(
x=pixels, y=theta, mode="lines", name="true row θ",
line=dict(color=ORANGE, width=2, dash="dash"),
hovertemplate="pixel %{x}<br>θ = %{y:.3f}<extra></extra>",
))
fig.add_trace(go.Scatter(
x=pixels, y=X @ theta, mode="lines", name="blurred row <b>X</b>θ",
line=dict(color=BLUE, width=2),
hovertemplate="pixel %{x}<br>Xθ = %{y:.3f}<extra></extra>",
))
fig.add_trace(go.Scatter(
x=pixels, y=y, mode="markers", name="data <b>y</b>",
marker=dict(color=INK, size=6, opacity=0.8),
hovertemplate="pixel %{x}<br>y = %{y:.3f}<extra></extra>",
))
style(fig)
fig.update_xaxes(title="pixel")
fig.update_yaxes(title="brightness")
fig.show()
```
The blur rounds off every sharp edge, and the noise is too small to see.
Now undo the blur.
```{python}
theta_hat = np.linalg.solve(X, y)
theta_noiseless = np.linalg.solve(X, X @ theta) # what we would get with no noise
print(f"||theta|| = {np.linalg.norm(theta):.2f}")
print(f"||theta_hat - theta|| = {np.linalg.norm(theta_hat - theta):.2e}")
print(f"||X theta_hat - y|| = {np.linalg.norm(X @ theta_hat - y):.0e}")
print(f"with no noise, ||error|| = {np.linalg.norm(theta_noiseless - theta):.0e}")
```
```{python}
#| label: fig-row-estimate
#| fig-cap: "The least squares estimate of row 50 (blue) and the true row (orange, dashed). On this scale, the true row is a flat line at zero."
#| fig-alt: "A line chart against pixel position from 0 to 127. A blue line oscillates up and down from one pixel to the next, with amplitude in the millions. An orange dashed line, the true row, is flat along zero."
#| code-fold: true
#| code-summary: "Plot the estimate"
fig = go.Figure()
fig.add_trace(go.Scatter(
x=pixels, y=theta_hat, mode="lines", name="least squares estimate",
line=dict(color=BLUE, width=1.5),
hovertemplate="pixel %{x}<br>estimate = %{y:,.0f}<extra></extra>",
))
fig.add_trace(go.Scatter(
x=pixels, y=theta, mode="lines", name="true row θ",
line=dict(color=ORANGE, width=2.5, dash="dash"),
hovertemplate="pixel %{x}<br>θ = %{y:.3f}<extra></extra>",
))
style(fig)
fig.update_xaxes(title="pixel")
fig.update_yaxes(title="brightness")
fig.show()
```
The estimate reproduces the data, $\bfX\widehat{\bftheta}=\bfy$ up to rounding error, and with no noise it would recover $\bftheta$ almost exactly.
So the trouble is the noise: one percent of noise in $\bfy$ became an error hundreds of thousands of times the size of $\bftheta$.
## The singular vectors of a blur
Since $X_{jk}=g(j-k)$ depends only on the distance $|j-k|$, the matrix $\bfX$ is symmetric.
Its eigenvalues also turn out to be positive, so its SVD is the same as its eigendecomposition: $\bfu_i=\bfv_i$, and
$$
\bfX\bfv_i=\sigma_i\bfv_i.
$$
Blurring $\bfv_i$ gives back $\bfv_i$, shrunk by the factor $\sigma_i$.
```{python}
U, s, Vt = np.linalg.svd(X)
V = Vt.T
print("largest singular values =", np.array2string(s[:3], precision=4))
print("smallest singular values =", np.array2string(s[-3:], precision=2))
print(f"condition number = {s[0] / s[-1]:.1e}")
print(f"max |U - V| = {np.abs(U - V).max():.0e}")
```
```{python}
#| label: fig-svd
#| fig-cap: "Left: the singular values of the blur, on a logarithmic scale. Right: four of its singular vectors. The later the singular vector, the faster it oscillates and the smaller its singular value."
#| fig-alt: "Left, a curve falling from 1 at index 1 to below one billionth at index 128 on a logarithmic axis, slowly at first and then ever faster. Right, four stacked line charts of singular vectors 1, 3, 10, and 128 against pixel position. Vector 1 is a single smooth hump, vector 3 has two sign changes, vector 10 oscillates nine times, and vector 128 flips sign at every pixel under a hump-shaped envelope."
#| code-fold: true
#| code-summary: "Plot the singular values and vectors"
shown = [1, 3, 10, d]
fig = make_subplots(
rows=4, cols=2, column_widths=[0.45, 0.55], horizontal_spacing=0.1, vertical_spacing=0.05,
specs=[[{"rowspan": 4}, {}], [None, {}], [None, {}], [None, {}]],
subplot_titles=["singular values σ<sub>i</sub>"] + [f"<b>v</b><sub>{i}</sub>" for i in shown],
)
fig.add_trace(go.Scatter(
x=np.arange(1, d + 1), y=s, mode="lines+markers", line=dict(color=BLUE, width=2),
marker=dict(size=4), hovertemplate="i = %{x}<br>σ<sub>i</sub> = %{y:.2e}<extra></extra>",
), row=1, col=1)
for r, i in enumerate(shown, start=1):
fig.add_trace(go.Scatter(
x=pixels, y=V[:, i - 1], mode="lines", line=dict(color=BLUE, width=1.5),
hovertemplate=f"pixel %{{x}}<br>v<sub>{i}</sub> = %{{y:.3f}}<extra></extra>",
), row=r, col=2)
fig.update_yaxes(range=[-0.16, 0.16], showticklabels=False, row=r, col=2)
fig.update_xaxes(showticklabels=r == 4, row=r, col=2)
style(fig, height=520, showlegend=False)
fig.update_yaxes(type="log", exponentformat="power", row=1, col=1)
fig.update_xaxes(title="i", row=1, col=1)
fig.update_xaxes(title="pixel", row=4, col=2)
fig.update_annotations(font=dict(color=INK, size=14))
fig.show()
```
The first singular vector is a single smooth hump, and the later ones oscillate faster and faster: $\bfv_i$ changes sign $i-1$ times, until $\bfv_{128}$ flips sign at every pixel.
Blurring averages neighbouring pixels.
Averaging barely changes a slowly varying vector, so $\sigma_1\approx1$.
A rapidly oscillating vector averages to almost nothing, because its positive and negative entries cancel, so $\sigma_{128}$ is below $10^{-8}$.
This is what it means for a picture to be blurry: the fine detail is almost gone.
Now write least squares in the SVD, using the [pseudoinverse](../../lectures/lecture-08.qmd#def-pseudoinverse):
$$
\widehat{\bftheta}_{\mathrm{LS}}=\sum_{i=1}^d\frac{\bfu_i^\top\bfy}{\sigma_i}\bfv_i,
\qquad
\bfu_i^\top\bfy=\underbrace{\sigma_i\,\bfv_i^\top\bftheta}_{\text{signal}}+\underbrace{\bfu_i^\top\bfvarepsilon}_{\text{noise}}.
$$
The blur shrinks the signal in each coefficient of the data by $\sigma_i$.
It does nothing to the noise: since $\bfu_i$ is a unit vector, $\bfu_i^\top\bfvarepsilon$ has variance $\sigma^2$ for every $i$ by [linear maps of random vectors](../../lectures/lecture-08.qmd#prp-linear-maps-moments).
@fig-picard compares the two.
```{python}
#| label: fig-picard
#| fig-cap: "For row 50, the size of each coefficient of the data |u_i^T y| (grey), the size of its signal part sigma_i |v_i^T theta| (blue), and the noise level sigma (orange, dashed), on a logarithmic scale. Past about i = 35 the signal falls below the noise, and the data coefficients are pure noise."
#| fig-alt: "A chart against the index i from 1 to 128 on a logarithmic vertical axis. A blue line falls steadily from about 1 to below one billionth. Grey dots follow the blue line until about i = 35, and then level off in a band around one hundredth, along an orange dashed horizontal line at 0.01."
#| code-fold: true
#| code-summary: "Plot the coefficients"
index = np.arange(1, d + 1)
fig = go.Figure()
fig.add_trace(go.Scatter(
x=index, y=s * np.abs(V.T @ theta), mode="lines", name="signal σ<sub>i</sub>|<b>v</b><sub>i</sub><sup>T</sup>θ|",
line=dict(color=BLUE, width=2),
hovertemplate="i = %{x}<br>signal = %{y:.1e}<extra></extra>",
))
fig.add_trace(go.Scatter(
x=[1, d], y=[sigma, sigma], mode="lines", name="noise level σ",
line=dict(color=ORANGE, width=2.5, dash="dash"), hoverinfo="skip",
))
fig.add_trace(go.Scatter(
x=index, y=np.abs(U.T @ y), mode="markers", name="data |<b>u</b><sub>i</sub><sup>T</sup><b>y</b>|",
marker=dict(color=INK, size=7, opacity=0.8, line=dict(color="white", width=1)),
hovertemplate="i = %{x}<br>data = %{y:.1e}<extra></extra>",
))
style(fig, height=440)
fig.update_xaxes(title="i")
fig.update_yaxes(title="size of coefficient", type="log", exponentformat="power")
fig.show()
```
For the first $35$ or so directions, the data coefficients follow the signal.
After that, the signal is buried under the noise: the blur has erased what $\bftheta$ looks like along those $\bfv_i$, and the data contain only noise there.
Least squares does not know this.
It divides each of those noise coefficients by $\sigma_i$, which is as small as $6\times10^{-9}$, and the result is the oscillating mess of @fig-row-estimate.
The corollary makes this quantitative before we simulate anything.
```{python}
var_theory = sigma**2 / s**2 # Var(v_i^T theta_hat), from the corollary
mse_theory = var_theory.sum() # E||theta_hat - theta||^2
print(f"sigma^2 / sigma_i^2 runs from {var_theory[0]:.1e} to {var_theory[-1]:.1e}")
print(f"E||theta_hat - theta||^2 = {mse_theory:.1e}")
print(f"||theta||^2, for comparison = {theta @ theta:.1f}")
```
## Many experiments
The blur and the sharp row stay fixed, so all of the randomness is in $\bfvarepsilon$.
We repeat the experiment $M=2000$ times: draw fresh noise $\bfvarepsilon^{(k)}$, set $\bfy^{(k)}=\bfX\bftheta+\bfvarepsilon^{(k)}$, and compute $\widehat{\bftheta}^{(k)}=\sum_i(\bfu_i^\top\bfy^{(k)}/\sigma_i)\bfv_i$.
Stacking the $\bfy^{(k)}$ as columns computes all $M$ estimates at once.
```{python}
M = 2000
Y = (X @ theta)[:, None] + sigma * rng.standard_normal((d, M)) # column k is y^(k)
thetas = (V @ ((U.T @ Y) / s[:, None])).T # row k is theta_hat^(k)
errors = thetas - theta
coords = errors @ V # column i holds v_i^T (theta_hat - theta)
```
@fig-variance compares the sample variance of each $\bfv_i^\top\widehat{\bftheta}^{(k)}$ over the $M$ experiments with the prediction $\sigma^2/\sigma_i^2$.
```{python}
#| label: fig-variance
#| fig-cap: "The sample variance of v_i^T theta_hat over 2000 experiments (blue) and the prediction sigma^2 / sigma_i^2 from the corollary (orange, dashed), on a logarithmic scale."
#| fig-alt: "A chart against the index i from 1 to 128 on a logarithmic vertical axis running from about one ten-thousandth to a few trillion. An orange dashed curve rises steadily, slowly at first and then ever faster, and blue dots sit on it all the way along."
#| code-fold: true
#| code-summary: "Plot the variances"
fig = go.Figure()
fig.add_trace(go.Scatter(
x=index, y=var_theory, mode="lines", name="predicted σ²/σ<sub>i</sub>²",
line=dict(color=ORANGE, width=2.5, dash="dash"),
hovertemplate="i = %{x}<br>predicted = %{y:.2e}<extra></extra>",
))
fig.add_trace(go.Scatter(
x=index, y=coords.var(axis=0), mode="markers", name=f"sample variance, {M} experiments",
marker=dict(color=BLUE, size=7, line=dict(color="white", width=1)),
hovertemplate="i = %{x}<br>sample = %{y:.2e}<extra></extra>",
))
style(fig, height=440)
fig.update_xaxes(title="i")
fig.update_yaxes(title="variance of <b>v</b><sub>i</sub><sup>T</sup>θ̂", type="log",
exponentformat="power")
fig.show()
```
The prediction holds across more than sixteen orders of magnitude.
Summing over $i$ gives the expected squared error, and nearly all of it comes from the last few directions.
```{python}
sq_errors = np.sum(errors**2, axis=1)
print(f"average of ||theta_hat - theta||^2 = {sq_errors.mean():.2e}")
print(f"predicted = {mse_theory:.2e}")
print(f"share of the prediction from v_128 = {var_theory[-1] / mse_theory:.2f}")
print(f"share from v_119, ..., v_128 = {var_theory[-10:].sum() / mse_theory:.4f}")
```
### Unbiased, but not usefully
Least squares is still unbiased, so the average of the estimates should approach $\bftheta$.
It does, but slowly.
By [averaging away independent noise](../../lectures/lecture-02.qmd#prp-average-variance), the average of $M$ independent estimates has expected squared error $\Ebb\|\widehat{\bftheta}-\bftheta\|_2^2/M$.
```{python}
mean_theta = thetas.mean(axis=0)
needed = mse_theory / (d * 0.01**2) # M with E||average - theta||^2 = d (0.01)^2
print(f"||average of the {M} estimates - theta|| = {np.linalg.norm(mean_theta - theta):.1e}")
print(f"predicted, sqrt(E||theta_hat - theta||^2 / M) = {np.sqrt(mse_theory / M):.1e}")
print(f"experiments needed for an error of 0.01 per pixel: {needed:.0e}")
print(f"at a million experiments per second, that takes {needed / 1e6 / 3.15e7:.0f} years")
```
Even the average of $2000$ experiments is off by tens of thousands, against a $\bftheta$ of size about $5$.
Unbiasedness is a statement about averaging infinitely many experiments; in practice we get one photograph.
## Giving up a little bias
Least squares fails because it insists on being unbiased, and so it divides by every $\sigma_i$, however small.
[Ridge regression](../../lectures/lecture-10.qmd#def-ridge) accepts some bias in exchange for much less variance.
By [ridge regression filters the singular values](../../lectures/lecture-10.qmd#thm-ridge-svd),
$$
\widehat{\bftheta}_\lambda=\sum_{i=1}^d\frac{\sigma_i^2}{\sigma_i^2+\lambda}\cdot\frac{\bfu_i^\top\bfy}{\sigma_i}\bfv_i.
$$
Directions with $\sigma_i^2\gg\lambda$, where the data still see the signal, are kept as they are.
Directions with $\sigma_i^2\ll\lambda$, where the data are pure noise, are suppressed instead of amplified.
The same works for the whole photograph.
Flattening $\bfTheta$ into a vector in $\Rbb^{d^2}$ turns $\bfTheta\mapsto\bfX\bfTheta\bfX^\top$ into a $d^2\times d^2$ matrix.
Its singular vectors are the flattened matrices $\bfv_i\bfv_j^\top$, with singular values $\sigma_i\sigma_j$, since $\bfX(\bfv_i\bfv_j^\top)\bfX^\top=\sigma_i\sigma_j\,\bfu_i\bfu_j^\top$ and these matrices are orthonormal in the Frobenius inner product.
So ridge regression for the photograph is
$$
\widehat{\bfTheta}_\lambda=\sum_{i,j=1}^d\frac{(\sigma_i\sigma_j)^2}{(\sigma_i\sigma_j)^2+\lambda}\cdot\frac{\bfu_i^\top\bfY\bfu_j}{\sigma_i\sigma_j}\,\bfv_i\bfv_j^\top.
$$
```{python}
S2 = np.outer(s, s) # entry (i, j) is the singular value sigma_i sigma_j
coeffs = U.T @ blurred @ U # entry (i, j) is u_i^T Y u_j
def ridge(lam):
return V @ (coeffs * S2 / (S2**2 + lam)) @ V.T
lams = np.concatenate([[0.0], 10.0 ** np.arange(-9, -0.5, 0.5)]) # lambda = 0 is least squares
ridges = [ridge(lam) for lam in lams]
rms = np.array([np.sqrt(np.mean((R - photo)**2)) for R in ridges])
best = int(np.argmin(rms))
print(f"best lambda = {lams[best]:.0e}, with RMS error {rms[best]:.3f} per pixel")
print(f"RMS error of the blurred photo itself: {np.sqrt(np.mean((blurred - photo)**2)):.3f}")
```
```{python}
#| label: fig-ridge
#| fig-cap: "The sharp photograph, the blurred photograph, and the ridge regression estimate. Drag the slider to change lambda; lambda = 0 is least squares."
#| fig-alt: "An interactive figure with three greyscale images side by side and a slider for lambda. The first two images are the sharp and the blurred photograph. The third is the ridge estimate: static for lambda near zero, a sharper version of the photograph for lambda near one hundredth, and increasingly blurry for larger lambda."
#| code-fold: true
#| code-summary: "Build the ridge slider"
# The sharp and blurred photographs are drawn once; only the third image changes with lambda.
fig = go.Figure(go.Image(source=png_uri(side_by_side([photo, blurred])), hoverinfo="skip"))
for k, R in enumerate(ridges):
fig.add_trace(go.Image(source=png_uri(R), x0=2 * (d + GAP), visible=k == best,
hoverinfo="skip"))
fig.update_layout(sliders=[dict(
active=best, currentvalue=dict(prefix="λ = ", font=dict(color=INK)),
pad=dict(t=30), len=0.9, x=0.05,
steps=[
dict(method="update", label="0" if lam == 0 else f"{lam:.0e}",
args=[{"visible": [True] + [j == k for j in range(len(lams))]},
{"title.text": f"RMS error of the estimate: {rms[k]:.3f} per pixel"}])
for k, lam in enumerate(lams)
],
)])
style(fig, height=430, showlegend=False,
title=f"RMS error of the estimate: {rms[best]:.3f} per pixel")
label_images(fig, ["sharp", "blurred, with noise", "ridge estimate"])
fig.update_xaxes(range=[-0.5, 3 * d + 2 * GAP - 0.5])
fig.update_layout(margin=dict(l=20, r=20, t=50, b=100))
fig.show()
```
For tiny $\lambda$, ridge regression behaves like least squares and amplifies the noise; that is the variance.
For large $\lambda$, it suppresses directions the data could have recovered, and the estimate stays blurry; that is the bias.
In between, the estimate is sharper than the blurred photograph: the legs of the tripod and the edges of the coat come back into focus.
We chose $\lambda$ here by comparing with the sharp photograph, which in practice we would not have, so this is the best case.
The supplementary notes on [choosing hyperparameters](../../notes/choosing-hyperparameters.qmd) pick $\lambda$ from the blurred photograph alone, by cross-validation, for this same problem.
## A teaser: total variation
Even with the best $\lambda$, ridge regression cannot recover the finest detail: along directions where the blur pushed the signal below the noise, that detail is gone for good, and ridge regression simply declines to make it up.
The penalty $\lambda\|\bfTheta\|_F^2$ knows nothing about photographs, except that their brightnesses should not be large.
We know more.
This photograph is mostly flat patches (the coat, the sky, the grass) separated by sharp edges.
The *total variation* of $\bfTheta$,
$$
\operatorname{TV}(\bfTheta)=\sum_{p,q=1}^d\sqrt{(\Theta_{p+1,q}-\Theta_{p,q})^2+(\Theta_{p,q+1}-\Theta_{p,q})^2},
$$
adds up the size of the jump from each pixel to its neighbours below and to the right, where differences past the last row or column count as zero.
It is small for a picture made of a few flat patches, however sharp the edges between them, and large for a noisy one.
So we estimate
$$
\widehat{\bfTheta}_{\mathrm{TV}}=\argmin_{\bfTheta}\ \tfrac12\|\bfX\bfTheta\bfX^\top-\bfY\|_F^2+\lambda\operatorname{TV}(\bfTheta).
$$
Like the lasso at the end of [Lecture 10](../../lectures/lecture-10.qmd), this penalty is a sum of sizes rather than of squares, and it prefers many of the differences to be exactly zero.
Unlike ridge regression, it has no formula in terms of the SVD.
We compute it with an iterative algorithm, which is beyond this course; the code is folded below.
```{python}
#| code-fold: true
#| code-summary: "A total variation solver"
def gradient(T):
"""The differences from each pixel to the next one down and the next one right."""
down, right = np.zeros_like(T), np.zeros_like(T)
down[:-1] = T[1:] - T[:-1]
right[:, :-1] = T[:, 1:] - T[:, :-1]
return down, right
def divergence(down, right):
"""Minus the adjoint (transpose) of gradient."""
out = np.zeros_like(down)
out[:-1] += down[:-1]
out[1:] -= down[:-1]
out[:, :-1] += right[:, :-1]
out[:, 1:] -= right[:, :-1]
return out
def total_variation(Xr, Yr, lam, iters=300, tau=10.0):
"""Minimize ||Xr T X^T - Yr||_F^2 / 2 + lam TV(T) by the primal-dual method of Chambolle and Pock.
The data term only enters through the step "solve (I + tau A^T A) T = rhs" for the map
A: T -> Xr T X^T, which the SVDs of Xr and X make diagonal, as for ridge regression.
"""
_, s1, V1t = np.linalg.svd(Xr) # V1 is d x d, so pad s1 with zeros if Xr has fewer rows
V1 = V1t.T
S = np.outer(np.concatenate([s1, np.zeros(d - len(s1))]), s)
data_term = Xr.T @ Yr @ X
step = 1 / (8 * tau) # the dual step: tau * step * ||gradient||^2 <= 1
T = np.zeros((d, d))
T_bar = T.copy()
p_down, p_right = np.zeros((d, d)), np.zeros((d, d))
for _ in range(iters):
down, right = gradient(T_bar)
p_down, p_right = p_down + step * down, p_right + step * right
shrink = np.maximum(1, np.sqrt(p_down**2 + p_right**2) / lam)
p_down, p_right = p_down / shrink, p_right / shrink
rhs = T + tau * divergence(p_down, p_right) + tau * data_term
T_new = V1 @ ((V1.T @ rhs @ V) / (1 + tau * S**2)) @ V.T
T_bar = 2 * T_new - T
T = T_new
return T
```
As for ridge regression, we choose $\lambda$ here by comparing with the sharp photograph.
```{python}
tv_lams = 10.0 ** np.arange(-5, -1.5, 0.5)
tvs = [total_variation(X, blurred, lam) for lam in tv_lams]
tv_rms = np.array([np.sqrt(np.mean((T - photo)**2)) for T in tvs])
tv_best = int(np.argmin(tv_rms))
tv_estimate = tvs[tv_best]
ridge_estimate = ridges[best]
print(f"best lambda for total variation = {tv_lams[tv_best]:.0e}")
print(f"RMS error per pixel, ridge regression: {rms[best]:.3f}")
print(f"RMS error per pixel, total variation: {tv_rms[tv_best]:.3f}")
```
```{python}
#| label: fig-tv
#| fig-cap: "The sharp photograph, the blurred photograph, and the two estimates, each with the lambda that brings it closest to the sharp photograph: ridge regression and total variation."
#| fig-alt: "Four greyscale images side by side. The first is the sharp photograph of the cameraman and the second is the blurred one. The third, the ridge estimate, is sharper than the blurred photograph but covered in a mottled ripple pattern. The fourth, the total variation estimate, has clean flat regions and crisp edges on the coat, the camera, and the tripod, but the grass and sky are flattened into blotchy patches, and the picture looks cartoon-like."
#| code-fold: true
#| code-summary: "Plot the estimates"
fig = go.Figure(go.Image(source=png_uri(side_by_side([photo, blurred, ridge_estimate, tv_estimate])),
hoverinfo="skip"))
style(fig, height=280, showlegend=False)
label_images(fig, ["sharp", "blurred, with noise", "ridge regression", "total variation"])
fig.update_layout(margin=dict(l=20, r=20, t=20, b=40))
fig.show()
```
```{python}
#| label: fig-tv-row
#| fig-cap: "Row 50 of the sharp photograph (orange, dashed), of the ridge estimate (grey), and of the total variation estimate (blue)."
#| fig-alt: "A line chart of brightness against pixel position from 0 to 127. The orange dashed true row has flat stretches and sharp jumps. The grey ridge estimate follows it but wiggles throughout and overshoots at each jump. The blue total variation estimate is a staircase of flat steps that jumps where the true row does, but it cuts off the top of the narrow peak near pixel 45."
#| code-fold: true
#| code-summary: "Plot row 50 of the estimates"
fig = go.Figure()
fig.add_trace(go.Scatter(
x=pixels, y=ridge_estimate[row], mode="lines", name="ridge regression",
line=dict(color=INK, width=1.5),
hovertemplate="pixel %{x}<br>ridge = %{y:.3f}<extra></extra>",
))
fig.add_trace(go.Scatter(
x=pixels, y=tv_estimate[row], mode="lines", name="total variation",
line=dict(color=BLUE, width=2),
hovertemplate="pixel %{x}<br>TV = %{y:.3f}<extra></extra>",
))
fig.add_trace(go.Scatter(
x=pixels, y=theta, mode="lines", name="true row θ",
line=dict(color=ORANGE, width=2, dash="dash"),
hovertemplate="pixel %{x}<br>θ = %{y:.3f}<extra></extra>",
))
style(fig)
fig.update_xaxes(title="pixel")
fig.update_yaxes(title="brightness")
fig.show()
```
Total variation is not unbiased either, and its bias is of a different kind.
Where the photograph really is flat patches and sharp edges, such as the coat and the tripod, it recovers edges that ridge regression leaves soft, because it is willing to put back detail the data could not see as long as that detail is a clean jump.
Where the photograph is textured, such as the grass, it flattens the texture into patches, and it shaves the top off narrow features, such as the bright peak near pixel $45$ in @fig-tv-row.
The choice of penalty is a choice of what we believe the answer looks like, and a good belief can buy back much of what the blur destroyed.