17 Autoencoder — Exercises¶
Test your understanding of autoencoders: reconstruction gradients, the PCA connection, the VAE KL term, and the reparameterization trick.
Prerequisites. Read theory.md and work through first_principles.ipynb before attempting these.
import random
import matplotlib.pyplot as plt
import numpy as np
%load_ext autoreload
%autoreload 2
SEED = 42
random.seed(SEED)
rng = np.random.default_rng(SEED)
plt.rcParams["figure.dpi"] = 90
Exercise 1 — Hand Derivation: Gradients Through a 1-Hidden-Unit Linear Autoencoder¶
Consider the smallest possible autoencoder: input $x \in \mathbb{R}^2$, a single latent unit ($k = 1$), and no biases or activations:
$$z = w_e^\top x, \qquad \hat{x} = z\, w_d, \qquad L = \lVert x - \hat{x} \rVert_2^2$$
with encoder weights $w_e \in \mathbb{R}^2$ and decoder weights $w_d \in \mathbb{R}^2$.
Task: Using the chain rule, derive closed-form expressions for $\nabla_{w_d} L$ and $\nabla_{w_e} L$ in terms of the residual $r = x - \hat{x}$. Then evaluate both gradients by hand at
$$x = \begin{bmatrix} 1 \\ 2 \end{bmatrix}, \qquad w_e = \begin{bmatrix} 0.5 \\ -0.5 \end{bmatrix}, \qquad w_d = \begin{bmatrix} 1 \\ 0.5 \end{bmatrix}$$
Questions:
- Compute $z$, $\hat{x}$, $r$, and $L$.
- Derive and evaluate $\nabla_{w_d} L$.
- Derive and evaluate $\nabla_{w_e} L$. Which factor plays the role of the "backpropagated error" arriving at the latent unit?
Expected results:
| Quantity | Value |
|---|---|
| $z$ | $-0.5$ |
| $\hat{x}$ | $(-0.5,\ -0.25)$ |
| $r$ | $(1.5,\ 2.25)$ |
| $L$ | $7.3125$ |
| $\nabla_{w_d} L$ | $(1.5,\ 2.25)$ |
| $\nabla_{w_e} L$ | $(-5.25,\ -10.5)$ |
Solution 1¶
Forward pass. $z = 0.5 \cdot 1 + (-0.5) \cdot 2 = -0.5$, so $\hat{x} = -0.5 \cdot (1,\ 0.5) = (-0.5,\ -0.25)$, the residual is $r = x - \hat{x} = (1.5,\ 2.25)$, and $L = 1.5^2 + 2.25^2 = 7.3125$.
Decoder gradient. Only $\hat{x} = z\, w_d$ depends on $w_d$ ($z$ is a function of $w_e$ alone). By the chain rule,
$$\nabla_{w_d} L = \frac{\partial L}{\partial \hat{x}} \frac{\partial \hat{x}}{\partial w_d} = (-2r) \cdot z = -2 z\, r$$
Evaluating: $-2 \cdot (-0.5) \cdot (1.5,\ 2.25) = (1.5,\ 2.25)$.
Encoder gradient. $w_e$ affects $L$ only through the scalar $z$. First backpropagate to $z$ (sum over output components, chain rule):
$$\frac{\partial L}{\partial z} = \sum_{j=1}^{2} \frac{\partial L}{\partial \hat{x}_j} \frac{\partial \hat{x}_j}{\partial z} = -2 r^\top w_d$$
This scalar is the backpropagated error at the latent unit. Then, since $\partial z / \partial w_e = x$,
$$\nabla_{w_e} L = \frac{\partial L}{\partial z} \cdot x = -2 (w_d^\top r)\, x$$
Evaluating: $w_d^\top r = 1 \cdot 1.5 + 0.5 \cdot 2.25 = 2.625$, so $\nabla_{w_e} L = -2 \cdot 2.625 \cdot (1,\ 2) = (-5.25,\ -10.5)$.
Result: $\nabla_{w_d} L = -2 z\, r$ and $\nabla_{w_e} L = -2 (w_d^\top r)\, x$.
# Verify the hand derivation numerically with central finite differences
x = np.array([1.0, 2.0])
w_e = np.array([0.5, -0.5])
w_d = np.array([1.0, 0.5])
def recon_loss(w_e, w_d, x):
z = w_e @ x
x_hat = z * w_d
return np.sum((x - x_hat) ** 2)
# Analytic gradients from the solution
z = w_e @ x
r = x - z * w_d
grad_wd = -2.0 * z * r
grad_we = -2.0 * (w_d @ r) * x
print(f"z = {z}, x_hat = {z * w_d}, r = {r}, L = {recon_loss(w_e, w_d, x)}")
print(f"analytic grad_wd = {grad_wd}")
print(f"analytic grad_we = {grad_we}")
# Central finite differences
h = 1e-6
num_wd = np.zeros(2)
num_we = np.zeros(2)
for j in range(2):
e = np.zeros(2)
e[j] = h
num_wd[j] = (recon_loss(w_e, w_d + e, x) - recon_loss(w_e, w_d - e, x)) / (2 * h)
num_we[j] = (recon_loss(w_e + e, w_d, x) - recon_loss(w_e - e, w_d, x)) / (2 * h)
assert np.allclose(grad_wd, num_wd, atol=1e-6), f"decoder gradient mismatch: {grad_wd} vs {num_wd}"
assert np.allclose(grad_we, num_we, atol=1e-6), f"encoder gradient mismatch: {grad_we} vs {num_we}"
# Match the expected hand-computed values
assert np.isclose(recon_loss(w_e, w_d, x), 7.3125, atol=1e-12)
assert np.allclose(grad_wd, [1.5, 2.25], atol=1e-12)
assert np.allclose(grad_we, [-5.25, -10.5], atol=1e-12)
print("\nAnalytic gradients match finite differences and hand values. All checks passed.")
z = -0.5, x_hat = [-0.5 -0.25], r = [1.5 2.25], L = 7.3125 analytic grad_wd = [1.5 2.25] analytic grad_we = [ -5.25 -10.5 ] Analytic gradients match finite differences and hand values. All checks passed.
Exercise 2 — Coding: Linear Autoencoder Recovers the Top Principal Component¶
Theory §3.2 (Baldi & Hornik, 1989) says a linear autoencoder with a $k$-dimensional bottleneck spans the same subspace as the top-$k$ principal components. Verify the $k = 1$ case by training one and comparing against PCA computed via SVD.
Task: Build a linear autoencoder $4 \to 1 \to 4$ from two Dense layers
of the repo library (ml_first_principles.nn_core) and train it with
full-batch gradient descent on the centered dataset generated below
(a noisy 1D line embedded in $\mathbb{R}^4$).
Requirements:
- Train with the MSE gradient $\dfrac{\partial L}{\partial \hat{X}} = \dfrac{2(\hat{X} - X)}{nd}$,
backpropagating decoder $\to$ encoder (
Dense.backwardapplies the update). - After training, compare the normalized decoder column with the top right
singular vector $v_1$ of the centered data: $\vert \cos \angle(w_d, v_1) \vert$
must equal $1$ within
atol=1e-4. - The autoencoder reconstruction MSE must match the rank-1 PCA reconstruction
MSE within
atol=1e-4.
Hint: lr=0.2 and 6000 steps converge comfortably. The sign of $w_d$ is not
identifiable (flipping the signs of $w_e, w_d$ together leaves $\hat{X}$
unchanged), hence the absolute value in the cosine check.
from ml_first_principles.nn_core import Dense
# Dataset: noisy 1D line embedded in R^4, mean-centered
n, d = 60, 4
direction = np.array([2.0, -1.0, 0.5, 1.5])
direction /= np.linalg.norm(direction)
t_latent = rng.normal(size=(n, 1)) * 2.0
X_data = t_latent * direction + 0.05 * rng.normal(size=(n, d))
X_c = X_data - X_data.mean(axis=0)
# PCA reference via SVD (stable, per repo standards)
U, S, Vt = np.linalg.svd(X_c, full_matrices=False)
v1 = Vt[0]
X_hat_pca = (X_c @ v1[:, None]) @ v1[None, :]
mse_pca = float(np.mean((X_c - X_hat_pca) ** 2))
print(f"PCA rank-1 reconstruction MSE: {mse_pca:.6f}")
# TODO: build encoder = Dense(d, 1, ...), decoder = Dense(1, d, ...)
# TODO: full-batch GD loop:
# Z = encoder.forward(X_c); X_hat = decoder.forward(Z)
# grad = 2 * (X_hat - X_c) / (n * d)
# encoder.backward(decoder.backward(grad, lr), lr)
# TODO: check cosine alignment with v1 and MSE against mse_pca
PCA rank-1 reconstruction MSE: 0.001536
Solution 2¶
encoder = Dense(d, 1, random_state=SEED)
decoder = Dense(1, d, random_state=SEED + 1)
lr = 0.2
n_steps = 6000
for step in range(n_steps):
Z = encoder.forward(X_c)
X_hat = decoder.forward(Z)
grad = 2.0 * (X_hat - X_c) / (n * d)
encoder.backward(decoder.backward(grad, lr), lr)
X_hat = decoder.forward(encoder.forward(X_c))
mse_ae = float(np.mean((X_c - X_hat) ** 2))
# The decoder column spans the learned 1D subspace
w_dir = decoder.weights.ravel()
w_dir = w_dir / np.linalg.norm(w_dir)
cos_align = abs(float(w_dir @ v1))
print(f"Linear AE reconstruction MSE: {mse_ae:.6f} (PCA: {mse_pca:.6f})")
print(f"cos angle(decoder direction, top PC) = {cos_align:.8f}")
# Deterministic checks
assert np.allclose(cos_align, 1.0, atol=1e-4), f"subspace mismatch: cos = {cos_align}"
assert np.allclose(mse_ae, mse_pca, atol=1e-4), f"MSE gap too large: {mse_ae - mse_pca:.2e}"
print("Linear AE recovers the top principal component. All checks passed.")
Linear AE reconstruction MSE: 0.001541 (PCA: 0.001536) cos angle(decoder direction, top PC) = 1.00000000 Linear AE recovers the top principal component. All checks passed.
Exercise 3 — Hand Calculation: Closed-Form VAE KL Divergence¶
Theory §6.4 gives, for encoder posterior $q = \mathcal{N}(\mu, \operatorname{diag}(\sigma^2))$ and prior $p = \mathcal{N}(0, I)$ in $k$ dimensions:
$$D_{\text{KL}}(q \,\Vert\, p) = \frac{1}{2}\sum_{j=1}^{k}\left(\mu_j^2 + \sigma_j^2 - \log\sigma_j^2 - 1\right)$$
Task: For a $k = 2$ latent code with
$$\mu = (0.5,\ -1.0), \qquad \log\sigma^2 = \bigl(0,\ \log 0.25\bigr)$$
compute $D_{\text{KL}}(q \,\Vert\, p)$ by hand, term by term.
Questions:
- Which dimension contributes more KL, and why?
- What values of $(\mu_j, \sigma_j^2)$ make dimension $j$'s contribution exactly zero?
- What does the answer to question 2 imply about the failure mode where the KL term drops to zero for every dimension (posterior collapse, theory §7.3)?
Expected result: $D_{\text{KL}} = \frac{1}{2}(0.25 + 1.636294) \approx 0.943147$.
Solution 3¶
Per-dimension terms $\tfrac{1}{2}(\mu_j^2 + \sigma_j^2 - \log\sigma_j^2 - 1)$:
- Dimension 1: $\mu_1^2 = 0.25$, $\sigma_1^2 = 1$, $\log\sigma_1^2 = 0$, so the term is $\tfrac{1}{2}(0.25 + 1 - 0 - 1) = 0.125$.
- Dimension 2: $\mu_2^2 = 1$, $\sigma_2^2 = 0.25$, $\log\sigma_2^2 = \log 0.25 \approx -1.386294$, so the term is $\tfrac{1}{2}(1 + 0.25 + 1.386294 - 1) = \tfrac{1}{2}(1.636294) \approx 0.818147$.
$$D_{\text{KL}} = 0.125 + 0.818147 \approx 0.943147$$
Answers.
- Dimension 2 dominates: its mean is farther from 0 and its variance $0.25$ is far from 1 — the function $\sigma^2 - \log\sigma^2 - 1$ is zero at $\sigma^2 = 1$ and grows in both directions.
- The contribution is zero iff $\mu_j = 0$ and $\sigma_j^2 = 1$, i.e. the posterior marginal equals the prior for that dimension.
- If the KL is zero for every dimension, $q(z \mid x) = p(z)$ for all $x$: the code $z$ is statistically independent of the input and carries no information — exactly the posterior-collapse failure of theory §7.3.
Result: $D_{\text{KL}}(q \,\Vert\, p) \approx 0.943147$.
# Verify: closed form vs Monte Carlo estimate of E_q[log q(z) - log p(z)]
mu = np.array([0.5, -1.0])
log_var = np.array([0.0, np.log(0.25)])
sigma2 = np.exp(log_var)
kl_closed = 0.5 * np.sum(mu**2 + sigma2 - log_var - 1)
print(f"Closed-form KL: {kl_closed:.6f}")
assert np.isclose(kl_closed, 0.943147, atol=1e-5), "hand calculation mismatch"
# Monte Carlo check with a local seeded generator
rng_mc = np.random.default_rng(SEED)
n_mc = 200_000
z_samples = mu + np.sqrt(sigma2) * rng_mc.normal(size=(n_mc, 2))
log_q = -0.5 * np.sum((z_samples - mu) ** 2 / sigma2 + np.log(2 * np.pi * sigma2), axis=1)
log_p = -0.5 * np.sum(z_samples**2 + np.log(2 * np.pi), axis=1)
kl_mc = float(np.mean(log_q - log_p))
print(f"Monte Carlo KL ({n_mc} samples): {kl_mc:.6f}")
assert np.allclose(kl_mc, kl_closed, atol=0.01), f"MC estimate off: {kl_mc} vs {kl_closed}"
print("Closed form matches the Monte Carlo estimate. All checks passed.")
Closed-form KL: 0.943147 Monte Carlo KL (200000 samples): 0.939282 Closed form matches the Monte Carlo estimate. All checks passed.
Exercise 4 — Coding: Reparameterization Trick¶
The VAE needs gradients of $\mathbb{E}_{z \sim q_\theta(z \mid x)}[\cdot]$ with respect to $\theta$, but sampling is not differentiable. Theory §6.5 rewrites the sample as a deterministic function of the parameters plus external noise:
$$z = \mu + \sigma \odot \varepsilon, \qquad \varepsilon \sim \mathcal{N}(0, I), \qquad \sigma = \exp\bigl(\tfrac{1}{2}\log\sigma^2\bigr)$$
Task: Implement reparameterize(mu, log_var, rng) taking a local
np.random.default_rng generator (repo standard: never mutate NumPy global
state).
Requirements:
- With $\mu = (1, -2)$ and $\sigma = (2, 0.3)$, the empirical mean and
standard deviation over 50,000 draws must match $\mu$ and $\sigma$ within
atol=0.03. - As $\log\sigma^2 \to -\infty$ the sample must collapse to $\mu$ exactly
(check with
atol=1e-10) — the trick degrades gracefully to a deterministic autoencoder. - In one sentence: why can gradients flow through $\mu$ and $\sigma$ here but not through a direct call to a Gaussian sampler parameterized by them?
# Solution 4
def reparameterize(mu, log_var, rng):
"""Sample z = mu + sigma * eps with eps ~ N(0, I) from a local generator."""
std = np.exp(0.5 * log_var)
eps = rng.normal(size=np.shape(mu))
return mu + std * eps
mu_r = np.array([1.0, -2.0])
log_var_r = np.array([np.log(4.0), np.log(0.09)]) # sigma = (2, 0.3)
# Requirement 1: sample statistics match (mu, sigma)
rng_rep = np.random.default_rng(SEED)
draws = np.stack([reparameterize(mu_r, log_var_r, rng_rep) for _ in range(50_000)])
print(f"empirical mean: {draws.mean(axis=0)} (target {mu_r})")
print(f"empirical std: {draws.std(axis=0)} (target {np.exp(0.5 * log_var_r)})")
assert np.allclose(draws.mean(axis=0), mu_r, atol=0.03), "mean mismatch"
assert np.allclose(draws.std(axis=0), np.exp(0.5 * log_var_r), atol=0.03), "std mismatch"
# Requirement 2: sigma -> 0 collapses to a deterministic code
z_det = reparameterize(mu_r, np.full(2, -60.0), np.random.default_rng(0))
assert np.allclose(z_det, mu_r, atol=1e-10), "should collapse to mu when sigma -> 0"
print("sigma -> 0 gives z == mu (deterministic autoencoder limit)")
print("All reparameterization checks passed.")
# Requirement 3 (answer): z is now a *deterministic, differentiable* function of
# (mu, sigma) with all randomness in the parameter-free noise eps, so
# dz/dmu = 1 and dz/dsigma = eps exist; sampling directly from N(mu, sigma^2)
# has no such differentiable path from parameters to sample.
empirical mean: [ 0.98777821 -2.00070649] (target [ 1. -2.]) empirical std: [1.99751551 0.30258633] (target [2. 0.3]) sigma -> 0 gives z == mu (deterministic autoencoder limit) All reparameterization checks passed.
Exercise 5 — Failure Analysis: Undercomplete vs Overcomplete (Identity Collapse)¶
Theory §2.3 and §7.1 warn that an overcomplete autoencoder ($k \geq d$) with no regularization can learn a useless identity map.
Questions:
- For a linear autoencoder with $k = d$ and no constraints, exhibit weights that achieve exactly zero reconstruction loss on any dataset. Why does zero loss not imply useful features here?
- An undercomplete autoencoder ($k < d$) cannot pull this trick. What does its reconstruction error on arbitrary (off-manifold) inputs tell you that the overcomplete one's cannot?
- Name two mechanisms from theory §4–§6 that make an overcomplete code useful again, and state in one sentence what each one penalizes or corrupts to rule out the identity solution.
Your answers:
1: ...
2: ...
3: ...
Solution 5¶
- Take $W_e = I_d$, $W_d = I_d$ (biases zero). Then $\hat{x} = W_d W_e x = x$ for every $x \in \mathbb{R}^d$, so the loss is exactly zero on any dataset — including pure noise never seen in training. Zero loss is achieved without the code depending on the data distribution at all: the latent representation is just a copy, it compresses nothing, and it ranks an off-manifold garbage input as perfectly "reconstructable". The training objective is satisfied while the actual goal (learning the data manifold) is not.
- An undercomplete map $\hat{x} = W_d W_e x$ has rank at most $k < d$, so it must discard $d - k$ directions. To minimize loss it is forced to keep the directions where the training data actually varies (Exercise 2: the principal subspace). Consequently its reconstruction error is small only near the data manifold and large off it — which is why reconstruction error from an undercomplete (or otherwise regularized) autoencoder works as an anomaly score, while the identity-collapsed model scores everything as normal.
- Any two of:
- Denoising (§4): corrupts the input ($\tilde{x} = x + \varepsilon$) but scores against the clean target, so copying the input reproduces the noise and is no longer optimal — the map must point back toward the manifold.
- Sparsity (§5): an L1 or KL penalty on activations makes "all units active copying all coordinates" expensive, forcing few active units per input.
- KL to the prior (§6): pulls $q(z \mid x)$ toward $\mathcal{N}(0, I)$, charging a price (in nats) for information carried by the code, so only structure that pays for itself in reconstruction survives.
# Numerical demo: identity collapse scores garbage inputs as perfect
d5 = 3
W_e_id = np.eye(d5) # overcomplete "solution": k = d, pure copy
W_d_id = np.eye(d5)
# Arbitrary off-manifold inputs (pure noise, local generator)
rng_noise = np.random.default_rng(SEED)
X_garbage = rng_noise.normal(size=(10, d5))
X_rec_over = (X_garbage @ W_e_id) @ W_d_id
mse_over = float(np.mean((X_garbage - X_rec_over) ** 2))
assert np.allclose(X_rec_over, X_garbage, atol=1e-12), "identity map should be exact"
# Undercomplete k = 1: rank-1 map cannot reproduce arbitrary 3D inputs
v = np.array([1.0, 0.0, 0.0]) # any unit vector; rank-1 projector v v^T
X_rec_under = (X_garbage @ v[:, None]) @ v[None, :]
mse_under = float(np.mean((X_garbage - X_rec_under) ** 2))
print(f"Overcomplete identity map, MSE on pure noise: {mse_over:.2e} (perfect -> useless as anomaly score)")
print(f"Undercomplete rank-1 map, MSE on pure noise: {mse_under:.4f} (large -> flags off-manifold inputs)")
assert mse_under > 0.1, "rank-deficient map must fail on arbitrary inputs"
print("Identity collapse demonstrated: zero loss with zero learning.")
Overcomplete identity map, MSE on pure noise: 0.00e+00 (perfect -> useless as anomaly score) Undercomplete rank-1 map, MSE on pure noise: 0.4496 (large -> flags off-manifold inputs) Identity collapse demonstrated: zero loss with zero learning.