"""Autoencoders on a year of days at Krems and on a curved toy dataset: a linear autoencoder finds the subspace PCA finds, a smaller code costs reconstruction error, and a nonlinear autoencoder follows a curve that no linear one can. Purpose ------- Numerical companion to the glossary entry 'autoencoder'. An autoencoder learns an encoder and a decoder together, judged by how well the decoder rebuilds a data point from the code the encoder produced. Nothing in that criterion needs a label: what the decoder has to reproduce is the data point itself. Two datasets. The GeoSphere Austria weather station Krems (station id 3805) records eight measurements per day; this script downloads the records for 2024 from the GeoSphere data hub (dataset klima-v2-1d), writes them to autoencoder_weather.csv, and scales each measurement to zero sample mean and unit sample variance because temperatures, precipitation and pressure carry different units. The curved dataset is synthetic on purpose: points along a parabola in the plane, where the structure a nonlinear autoencoder can follow and a linear one cannot is known in advance. The demo checks the entry's claims: (1) a linear autoencoder trained by gradient descent on the weather data reaches the reconstruction error PCA attains in closed form, and reconstructs into the same subspace, the largest principal angle between the decoder's range and the principal subspace staying below a hundredth of a degree; (2) the code size controls what can be rebuilt, the error falling as the code grows and matching the sum of the eigenvalues the code drops; (3) on the curved dataset a nonlinear autoencoder with a code of one number beats the best linear one by more than a factor of ten, because a line cannot follow a parabola; (4) the hand-written gradient used for the training agrees with a finite-difference gradient, so the fits rest on a gradient that was checked rather than assumed. Deterministic: the weather data are a fixed archive year, PCA is an eigenvalue decomposition, and every initialization is drawn from a fixed seed. Self-contained: numpy + matplotlib only (stdlib urllib for the download). Blocks ------ [B-data] Download the 366 days with eight measurements each and scale them; build the curved dataset. [B-linear] Train a linear autoencoder (encoder and decoder both matrices) by gradient descent with code size two: check its reconstruction error reaches the PCA minimum, that the decoder reconstructs into the principal subspace, and that its reconstruction map agrees with the PCA projector on the data. The comparison is made on the decoder because the encoder's row space is not determined: the data carry a direction of almost no variance. [B-code] The code size is what forces the choice: the reconstruction error of the best linear autoencoder equals the sum of the dropped eigenvalues, and falls as the code grows. [B-gradcheck] The hand-written gradient of the nonlinear autoencoder agrees with a finite-difference gradient to eight digits. [B-nonlinear] On the curved dataset, a nonlinear autoencoder with one code number beats the best linear one by more than a factor of ten. [B-plot] Preview: the curved dataset with both reconstructions, and the reconstruction error against the code size. Outputs ------- autoencoder_weather.csv : date and the eight measurements autoencoder_curve.csv : x1, x2 -- the curved dataset autoencoder_linear.csv : x1, x2 -- its linear (PCA) reconstruction autoencoder_nonlinear.csv : x1, x2 -- its nonlinear reconstruction autoencoder_codesize.csv : code, error -- error against code size autoencoder.png : preview (checking only) """ import json import urllib.request from pathlib import Path import numpy as np import matplotlib matplotlib.use("Agg") import matplotlib.pyplot as plt OUT_DIR = Path(__file__).parent report = [] def check(name, ok): report.append((name, bool(ok))) print(f" [{'ok' if ok else 'FAIL'}] {name}") # ---- [B-data] the weather measurements and the curved dataset PARAMS = ["tlmin", "tlmax", "tl_mittel", "rr", "so_h", "rf_mittel", "p_mittel", "vv_mittel"] URL = ("https://dataset.api.hub.geosphere.at/v1/station/historical/" f"klima-v2-1d?parameters={','.join(PARAMS)}&station_ids=3805" "&start=2024-01-01&end=2024-12-31") with urllib.request.urlopen(URL, timeout=120) as resp: payload = json.load(resp) params = payload["features"][0]["properties"]["parameters"] stamps = [t[:10] for t in payload["timestamps"]] M = np.stack([np.array(params[p]["data"], dtype=float) for p in PARAMS], 1) M[:, 3] = np.maximum(M[:, 3], 0.0) # -1 marks a trace of rain with open(OUT_DIR / "autoencoder_weather.csv", "w") as f: f.write("date," + ",".join(PARAMS) + "\n") for day, row in zip(stamps, M): f.write(day + "," + ",".join(f"{v:g}" for v in row) + "\n") X = (M - M.mean(axis=0)) / M.std(axis=0) m, d = X.shape check("[B-data] 366 days with eight measurements each", (m, d) == (366, 8)) rng = np.random.default_rng(0) t = np.linspace(-1.5, 1.5, 300) C = np.stack([t, t ** 2], 1) + 0.03 * rng.standard_normal((300, 2)) C = C - C.mean(axis=0) # centered, for a fair line check("[B-data] the curved dataset lies along a parabola", len(C) == 300) def eig_spectrum(A): vals, vecs = np.linalg.eigh(A.T @ A / len(A)) order = np.argsort(vals)[::-1] return vals[order], vecs[:, order] lam, U = eig_spectrum(X) def recon_error(A, W, R): """sum_r || a^(r) - R W a^(r) ||^2 for an encoder W and a decoder R.""" return float(((A - A @ W.T @ R.T) ** 2).sum()) # ---- [B-linear] a linear autoencoder reaches the PCA optimum def train_linear(A, code, steps=4000, lr=0.02, seed=1): g = np.random.default_rng(seed) n, p = A.shape W = 0.1 * g.standard_normal((code, p)) R = 0.1 * g.standard_normal((p, code)) for _ in range(steps): Z = A @ W.T E = A - Z @ R.T # residual, n x p gR = -2.0 * E.T @ Z / n gW = -2.0 * (E @ R).T @ A / n R -= lr * gR W -= lr * gW return W, R CODE = 2 W_ae, R_ae = train_linear(X, CODE) err_ae = recon_error(X, W_ae, R_ae) err_pca = m * lam[CODE:].sum() print(f" linear autoencoder: reconstruction error {err_ae:.1f} against the " f"PCA minimum {err_pca:.1f}") check("[B-linear] the trained linear autoencoder reaches the PCA minimum " "(within one percent)", err_ae < 1.01 * err_pca) def principal_angles(B1, B2): """Angles in degrees between the subspaces spanned by the columns.""" Q1 = np.linalg.qr(B1)[0] Q2 = np.linalg.qr(B2)[0] s = np.clip(np.linalg.svd(Q1.T @ Q2, compute_uv=False), -1.0, 1.0) return np.degrees(np.arccos(s)) # The subspace the autoencoder RECONSTRUCTS into is the decoder's column # space, and that is what PCA's principal subspace is compared with. The # encoder's row space is not determined here: the measurements carry a # direction of almost no variance (see below), and an encoder may read # along it without changing any reconstruction. ang = principal_angles(R_ae, U[:, :CODE]) ang_enc = principal_angles(W_ae.T, U[:, :CODE]) print(f" largest principal angle to the PCA subspace: decoder " f"{ang.max():.4f} degrees, encoder {ang_enc.max():.2f} degrees") check("[B-linear] the decoder spans the subspace PCA spans (largest " "principal angle below a hundredth of a degree)", ang.max() < 0.01) check("[B-linear] the reconstruction map agrees with the PCA projector on " "the data", np.allclose(X @ (R_ae @ W_ae).T, X @ (U[:, :CODE] @ U[:, :CODE].T), atol=1e-3)) # ---- [B-code] the code size decides what can be rebuilt curve = [] for k in range(1, 6): Wk, Rk = train_linear(X, k, seed=2 + k) curve.append((k, recon_error(X, Wk, Rk) / m, m * lam[k:].sum() / m)) print(" error per data point by code size: " + ", ".join(f"{k}: {e:.2f} (PCA {p:.2f})" for k, e, p in curve)) check("[B-code] the error falls as the code grows", all(curve[i][1] > curve[i + 1][1] for i in range(len(curve) - 1))) check("[B-code] each equals the sum of the dropped eigenvalues (within two " "percent)", all(e < 1.02 * p + 1e-9 for _, e, p in curve)) with open(OUT_DIR / "autoencoder_codesize.csv", "w") as f: f.write("code,error,pca\n") for k, e, p in curve: f.write(f"{k},{e:.4f},{p:.4f}\n") # ---- the nonlinear autoencoder: one hidden layer on each side def init_nonlinear(p, hidden, seed=3): g = np.random.default_rng(seed) s = 0.8 return {"A1": s * g.standard_normal((p, hidden)), "b1": np.zeros(hidden), "a2": s * g.standard_normal((hidden, 1)), "c2": np.zeros(1), "a3": s * g.standard_normal((1, hidden)), "b3": np.zeros(hidden), "A4": s * g.standard_normal((hidden, p)), "b4": np.zeros(p)} def forward(P, A): H1 = np.tanh(A @ P["A1"] + P["b1"]) # encoder hidden Z = H1 @ P["a2"] + P["c2"] # the code, one number H2 = np.tanh(Z @ P["a3"] + P["b3"]) # decoder hidden Xh = H2 @ P["A4"] + P["b4"] # reconstruction return H1, Z, H2, Xh def loss_and_grad(P, A): n = len(A) H1, Z, H2, Xh = forward(P, A) E = A - Xh loss = float((E ** 2).sum() / n) dXh = -2.0 * E / n g = {"A4": H2.T @ dXh, "b4": dXh.sum(0)} dH2 = dXh @ P["A4"].T dpre2 = dH2 * (1.0 - H2 ** 2) g["a3"] = Z.T @ dpre2 g["b3"] = dpre2.sum(0) dZ = dpre2 @ P["a3"].T g["a2"] = H1.T @ dZ g["c2"] = dZ.sum(0) dH1 = dZ @ P["a2"].T dpre1 = dH1 * (1.0 - H1 ** 2) g["A1"] = A.T @ dpre1 g["b1"] = dpre1.sum(0) return loss, g # ---- [B-gradcheck] the hand-written gradient against finite differences P0 = init_nonlinear(2, 12) _, g0 = loss_and_grad(P0, C) worst = 0.0 probe = np.random.default_rng(4) for name in P0: flat = P0[name].ravel() for _ in range(3): i = int(probe.integers(flat.size)) eps, keep = 1e-6, flat[i] flat[i] = keep + eps lp = loss_and_grad(P0, C)[0] flat[i] = keep - eps lm = loss_and_grad(P0, C)[0] flat[i] = keep num = (lp - lm) / (2 * eps) worst = max(worst, abs(num - g0[name].ravel()[i])) print(f" gradient check: largest difference to finite differences {worst:.2e}") check("[B-gradcheck] the hand-written gradient agrees with finite differences", worst < 1e-8) def train_nonlinear(A, steps=20000, lr=0.05, seed=3, hidden=12): P = init_nonlinear(A.shape[1], hidden, seed) vel = {k: np.zeros_like(v) for k, v in P.items()} for _ in range(steps): _, g = loss_and_grad(P, A) for k in P: vel[k] = 0.9 * vel[k] - lr * g[k] P[k] = P[k] + vel[k] return P # ---- [B-nonlinear] a curve no line can follow lamC, UC = eig_spectrum(C) w_lin = UC[:, :1].T lin_recon = C @ w_lin.T @ w_lin err_lin = float(((C - lin_recon) ** 2).sum() / len(C)) P = train_nonlinear(C) nonlin_recon = forward(P, C)[3] err_non = float(((C - nonlin_recon) ** 2).sum() / len(C)) print(f" curved dataset, code of one number: linear error {err_lin:.4f}, " f"nonlinear error {err_non:.4f} (factor {err_lin / err_non:.1f})") check("[B-nonlinear] the nonlinear autoencoder beats the best linear one by " "more than a factor of ten", err_lin > 10.0 * err_non) np.savetxt(OUT_DIR / "autoencoder_curve.csv", C, delimiter=",", header="x1,x2", comments="", fmt="%.4f") np.savetxt(OUT_DIR / "autoencoder_linear.csv", lin_recon, delimiter=",", header="x1,x2", comments="", fmt="%.4f") order = np.argsort(nonlin_recon[:, 0]) np.savetxt(OUT_DIR / "autoencoder_nonlinear.csv", nonlin_recon[order], delimiter=",", header="x1,x2", comments="", fmt="%.4f") # ---- [B-plot] preview fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(10, 4.2)) ax1.scatter(C[:, 0], C[:, 1], s=9, color="0.6", label="data point") ax1.plot(lin_recon[np.argsort(lin_recon[:, 0]), 0], lin_recon[np.argsort(lin_recon[:, 0]), 1], "k--", lw=2, label="linear autoencoder (the PCA line)") ax1.plot(nonlin_recon[order, 0], nonlin_recon[order, 1], "k-", lw=2, label="nonlinear autoencoder") ax1.set_aspect("equal") ax1.set_xlabel("first feature $x_1$") ax1.set_ylabel("second feature $x_2$") ax1.set_title("One code number: a line cannot follow a curve") ax1.legend(frameon=False, fontsize=8) ks = [k for k, _, _ in curve] ax2.plot(ks, [e for _, e, _ in curve], "ko-", label="trained linear autoencoder") ax2.plot(ks, [p for _, _, p in curve], "s--", color="0.5", label="PCA (sum of the dropped eigenvalues)") ax2.set_xlabel("code size") ax2.set_ylabel("reconstruction error per data point") ax2.set_title("The code size decides what can be rebuilt") ax2.legend(frameon=False, fontsize=8) fig.tight_layout() fig.savefig(OUT_DIR / "autoencoder.png", dpi=110) n_ok = sum(ok for _, ok in report) print(f"\n{n_ok}/{len(report)} checks pass") if n_ok != len(report): raise SystemExit(1)