{
 "nbformat": 4,
 "nbformat_minor": 5,
 "metadata": {
  "kernelspec": {
   "name": "python3",
   "display_name": "Python 3",
   "language": "python"
  },
  "language_info": {
   "name": "python"
  },
  "colab": {
   "name": "pca.ipynb"
  }
 },
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "# principal component analysis \u2014 Python demo\n\nNumerical companion to the entry [principal component analysis](https://dictionaryofml.org/terms/pca.html) of the [Dictionary of Applied Machine Learning](https://dictionaryofml.org/): it recomputes what the entry states and prints one line per check.\n\nNumerical companion to the glossary entry 'pca'. The GeoSphere Austria weather station Krems (station id 3805) records eight measurements per day: minimum, maximum and mean air temperature, precipitation, sunshine duration, relative humidity, air pressure and wind speed. This script downloads the records for 2024 from the GeoSphere data hub (dataset klima-v2-1d) and writes them to pca_weather.csv. A value of -1 for precipitation marks a trace of rain too small to record and is set to 0.\n\nRequires NumPy and Matplotlib only, and uses fixed seeds, so the printed numbers reproduce exactly. Generated from [`pythondemos/pca.py`](https://dictionaryofml.org/terms/pca.py); CC BY 4.0."
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "# Notebook shim: the script resolves output paths relative to __file__,\n# which a notebook kernel does not define; everything lands in the\n# working directory instead.\nimport os\n__file__ = os.path.join(os.getcwd(), \"pca.py\")\nos.makedirs(\"pythondemos\", exist_ok=True)"
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "\"\"\"Principal component analysis on a year of days at Krems: the principal\ndirections are the eigenvectors of the sample covariance matrix, they\nmaximize the variance of the projection, and the reconstruction error\nthey leave is the sum of the eigenvalues they drop.\n\nPurpose\n-------\nNumerical companion to the glossary entry 'pca'.  The GeoSphere Austria\nweather station Krems (station id 3805) records eight measurements per\nday: minimum, maximum and mean air temperature, precipitation, sunshine\nduration, relative humidity, air pressure and wind speed.  This script\ndownloads the records for 2024 from the GeoSphere data hub (dataset\nklima-v2-1d) and writes them to pca_weather.csv.  A value of -1 for\nprecipitation marks a trace of rain too small to record and is set to 0.\n\nTwo views of the same data.  The picture uses the two temperatures in\ndegrees Celsius, centered but not rescaled, so the principal directions\ncan be drawn on the scatterplot in the units of the measurement.  The\nidentities use all eight measurements, each scaled to zero sample mean\nand unit sample variance because temperatures, precipitation and\npressure carry different units.\n\nThe demo checks the entry's claims: (1) the first principal direction\nmaximizes the variance of the projection, beating two thousand random\nunit vectors, and the directions are orthonormal; (2) the encoder W\nbuilt from the top eigenvectors minimizes the reconstruction error, and\nthat minimum equals the number of data points times the sum of the\ndropped eigenvalues; (3) the decoder that minimizes the error for such\nan encoder is its transpose, so the decoder is determined once the\nencoder is fixed; (4) PCA is ERM with the squared error loss, its\nobjective value matching the average reconstruction error; and (5) the\nsmallest eigenvalue, four orders of magnitude below the largest, names\na redundant measurement: the mean temperature of a day is the midpoint\nof its minimum and maximum.\n\nDeterministic: the data are a fixed archive year, the decomposition is\nan eigenvalue decomposition, and the random comparison vectors are\ndrawn with a fixed seed.  Self-contained: numpy + matplotlib only\n(stdlib urllib for the download).\n\nBlocks\n------\n[B-data]      Download the 366 days with eight measurements each; center\n              the two temperatures for the picture and scale all eight\n              for the identities.\n[B-directions] Eigenvalue decomposition of the sample covariance matrix\n              of the two temperatures: check the eigenvectors are\n              orthonormal, that the first maximizes the variance of the\n              projection against two thousand random unit vectors, and\n              that the eigenvalues are those variances.\n[B-reconstruct] All eight measurements, code size two: check the\n              reconstruction error of the PCA encoder equals the number\n              of data points times the sum of the dropped eigenvalues,\n              and that no random encoder reconstructs better.\n[B-decoder]   Which of the two maps may be fixed first: for a decoder\n              with orthonormal columns the best encoder is its transpose,\n              while a generic encoder with orthonormal rows does NOT have\n              its transpose as the best decoder; and restricting the\n              decoder to orthonormal columns costs nothing.\n[B-trace]     The reconstruction is orthogonal to the error it leaves, so\n              the reconstruction error is m times a gap between two traces;\n              the second trace is maximized by the top eigenvectors.\n[B-erm]       PCA as ERM: the average squared error loss of the learned\n              pair equals the reconstruction error divided by the number\n              of data points, and it falls as the code size grows.\n[B-redundant] The smallest eigenvalue is four orders of magnitude below\n              the largest, and its direction loads on the three\n              temperatures alone: the mean temperature of a day is the\n              midpoint of its minimum and maximum, to within the 0.05\n              degree the archive records. PCA finds the redundancy\n              without being told to look for one.\n[B-plot]      Preview: the days in the two temperatures with the two\n              principal directions, and the eigenvalue spectrum of the\n              eight measurements.\n\nOutputs\n-------\npca_weather.csv    : date and the eight measurements, 366 days of 2024\npca_points.csv     : x1, x2 -- the centered temperatures of each day\npca_axis1.csv,\npca_axis2.csv      : the two principal directions as segments, scaled by\n                     the square root of their eigenvalue\npca_spectrum.csv   : index, eigenvalue, cumulative share of the variance\npca.png            : preview (checking only)\n\"\"\"\n\nimport json\nimport urllib.request\nfrom pathlib import Path\n\nimport numpy as np\nimport matplotlib\n\nmatplotlib.use(\"Agg\")\nimport matplotlib.pyplot as plt\n\nOUT_DIR = Path(__file__).parent\n\nreport = []\n\n\ndef check(name, ok):\n    report.append((name, bool(ok)))\n    print(f\"  [{'ok' if ok else 'FAIL'}] {name}\")"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "**[B-data]** Download the 366 days with eight measurements each; center the two temperatures for the picture and scale all eight for the identities."
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "PARAMS = [\"tlmin\", \"tlmax\", \"tl_mittel\", \"rr\", \"so_h\", \"rf_mittel\",\n          \"p_mittel\", \"vv_mittel\"]\nURL = (\"https://dataset.api.hub.geosphere.at/v1/station/historical/\"\n       f\"klima-v2-1d?parameters={','.join(PARAMS)}&station_ids=3805\"\n       \"&start=2024-01-01&end=2024-12-31\")\nwith urllib.request.urlopen(URL, timeout=120) as resp:\n    payload = json.load(resp)\nparams = payload[\"features\"][0][\"properties\"][\"parameters\"]\nstamps = [t[:10] for t in payload[\"timestamps\"]]\nM = np.stack([np.array(params[p][\"data\"], dtype=float) for p in PARAMS], 1)\nM[:, 3] = np.maximum(M[:, 3], 0.0)             # -1 marks a trace of rain\nwith open(OUT_DIR / \"pca_weather.csv\", \"w\") as f:\n    f.write(\"date,\" + \",\".join(PARAMS) + \"\\n\")\n    for day, row in zip(stamps, M):\n        f.write(day + \",\" + \",\".join(f\"{v:g}\" for v in row) + \"\\n\")\nT = M[:, :2] - M[:, :2].mean(axis=0)           # centered, in degrees\nX = (M - M.mean(axis=0)) / M.std(axis=0)       # scaled, all eight\nm, d = X.shape\ncheck(\"[B-data] 366 days with eight measurements each\", (m, d) == (366, 8))\ncheck(\"[B-data] the record matches the archive (Feb 1: -3.8 to 10.4 deg)\",\n      stamps[31] == \"2024-02-01\" and np.allclose(M[31, :2], [-3.8, 10.4]))\ncheck(\"[B-data] the two temperatures are centered\", np.allclose(T.mean(axis=0), 0))\n\n\ndef principal(A):\n    \"\"\"Eigenvalues and eigenvectors of the sample covariance matrix of A,\n    in decreasing order of eigenvalue.\"\"\"\n    Q = A.T @ A / len(A)\n    vals, vecs = np.linalg.eigh(Q)\n    order = np.argsort(vals)[::-1]\n    return vals[order], vecs[:, order]"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "**[B-directions]** Eigenvalue decomposition of the sample covariance matrix of the two temperatures: check the eigenvectors are orthonormal, that the first maximizes the variance of the projection against two thousand random unit vectors, and that the eigenvalues are those variances."
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "lam2, U2 = principal(T)\nrng = np.random.default_rng(0)\nangles = rng.uniform(0, 2 * np.pi, 2000)\ndirs = np.stack([np.cos(angles), np.sin(angles)], 1)\nvar_random = ((T @ dirs.T) ** 2).mean(axis=0)\nvar_first = float(((T @ U2[:, 0]) ** 2).mean())\nprint(f\"  two temperatures: eigenvalues {lam2[0]:.2f} and {lam2[1]:.2f} \"\n      f\"(squared degrees); best random direction {var_random.max():.2f}\")\ncheck(\"[B-directions] the principal directions are orthonormal\",\n      np.allclose(U2.T @ U2, np.eye(2)))\ncheck(\"[B-directions] the first maximizes the variance of the projection \"\n      \"(2000 random unit vectors)\", var_first >= var_random.max() - 1e-9)\ncheck(\"[B-directions] the eigenvalues are the variances of the projections\",\n      np.allclose([((T @ U2[:, j]) ** 2).mean() for j in (0, 1)], lam2))"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "**[B-reconstruct]** All eight measurements, code size two: check the reconstruction error of the PCA encoder equals the number of data points times the sum of the dropped eigenvalues, and that no random encoder reconstructs better."
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "lam, U = principal(X)\nCODE = 2\nW = U[:, :CODE].T                              # encoder, CODE x d\n\n\ndef min_reconstruction_error(W, A):\n    \"\"\"min over decoders R of sum_r || a^(r) - R W a^(r) ||^2.\"\"\"\n    Z = A @ W.T\n    R = np.linalg.lstsq(Z, A, rcond=None)[0].T\n    return float(((A - Z @ R.T) ** 2).sum())\n\n\nerr_pca = min_reconstruction_error(W, X)\ndropped = m * lam[CODE:].sum()\nerr_random = [min_reconstruction_error(rng.standard_normal((CODE, d)), X)\n              for _ in range(50)]\nprint(f\"  eight measurements: reconstruction error {err_pca:.1f}, sum of the \"\n      f\"dropped eigenvalues times m {dropped:.1f}; best of 50 random \"\n      f\"encoders {min(err_random):.1f}\")\ncheck(\"[B-reconstruct] the error equals m times the sum of the dropped \"\n      \"eigenvalues\", np.isclose(err_pca, dropped))\ncheck(\"[B-reconstruct] no random encoder reconstructs better\",\n      err_pca < min(err_random))"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "**[B-decoder]** Which of the two maps may be fixed first: for a decoder with orthonormal columns the best encoder is its transpose, while a generic encoder with orthonormal rows does NOT have its transpose as the best decoder; and restricting the decoder to orthonormal columns costs nothing."
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "def best_code(R, A):\n    \"\"\"argmin_z || a - R z ||^2 for every row a of A, solved independently.\"\"\"\n    return np.linalg.lstsq(R, A.T, rcond=None)[0].T\n\n\northo = [np.linalg.qr(rng.standard_normal((d, CODE)))[0] for _ in range(200)]\ngap_dec = max(abs(float(((X - (X @ R) @ R.T) ** 2).sum())\n                  - float(((X - best_code(R, X) @ R.T) ** 2).sum()))\n              for R in ortho)\ncheck(\"[B-decoder] for a decoder with orthonormal columns the best encoder \"\n      \"is its transpose\", gap_dec < 1e-8)\n\nW_gen = np.linalg.qr(rng.standard_normal((d, CODE)))[0].T\nZ_gen = X @ W_gen.T\nR_ls = np.linalg.lstsq(Z_gen, X, rcond=None)[0].T\nerr_ls = float(((X - Z_gen @ R_ls.T) ** 2).sum())\nerr_t = float(((X - Z_gen @ W_gen) ** 2).sum())\nprint(f\"  generic orthonormal-row encoder: its best decoder leaves \"\n      f\"{err_ls:.1f}, its transpose {err_t:.1f}\")\ncheck(\"[B-decoder] a generic orthonormal-row encoder does NOT have its \"\n      \"transpose as the best decoder\", err_t > err_ls + 1.0)\nZ = X @ W.T\nR_hat = np.linalg.lstsq(Z, X, rcond=None)[0].T\ncheck(\"[B-decoder] the eigenvector encoder does, its rows spanning an \"\n      \"invariant subspace\", np.allclose(R_hat, W.T, atol=1e-8))\nfree = min(min_reconstruction_error(rng.standard_normal((CODE, d)), X)\n           for _ in range(300))\ncheck(\"[B-decoder] no unconstrained pair beats the orthonormal-decoder \"\n      \"optimum\", free > err_pca)"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "**[B-trace]** The reconstruction is orthogonal to the error it leaves, so the reconstruction error is m times a gap between two traces; the second trace is maximized by the top eigenvectors."
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "Q = X.T @ X / m\nR_top = U[:, :CODE]\npyth = float(np.abs((X ** 2).sum(1) - ((X @ R_top) ** 2).sum(1)\n                    - ((X - (X @ R_top) @ R_top.T) ** 2).sum(1)).max())\ncheck(\"[B-trace] the reconstruction is orthogonal to the error it leaves\",\n      pyth < 1e-9)\nerr_trace = m * (np.trace(Q) - np.trace(R_top.T @ Q @ R_top))\ncheck(\"[B-trace] the reconstruction error is m times the gap between the \"\n      \"trace of the sample covariance matrix and the projected trace\",\n      np.isclose(err_pca, err_trace))\nbest_tr = max(float(np.trace(R.T @ Q @ R)) for R in ortho)\nprint(f\"  trace: best of {len(ortho)} random orthonormal decoders \"\n      f\"{best_tr:.4f}, sum of the top {CODE} eigenvalues \"\n      f\"{lam[:CODE].sum():.4f}\")\ncheck(\"[B-trace] no random orthonormal decoder reaches the sum of the top \"\n      \"eigenvalues\", best_tr < lam[:CODE].sum())"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "**[B-erm]** PCA as ERM: the average squared error loss of the learned pair equals the reconstruction error divided by the number of data points, and it falls as the code size grows."
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "avg_loss = float(((X - Z @ W) ** 2).sum(axis=1).mean())\ncheck(\"[B-erm] the average squared error loss equals the reconstruction \"\n      \"error per data point\", np.isclose(avg_loss, err_pca / m))\ncurve = [(k, m * lam[k:].sum() / m) for k in range(1, d + 1)]\ncheck(\"[B-erm] the average loss falls as the code size grows\",\n      all(curve[i][1] > curve[i + 1][1] for i in range(len(curve) - 1)))\nshare = np.cumsum(lam) / lam.sum()\nprint(\"  average squared error loss by code size: \"\n      + \", \".join(f\"{k}: {v:.2f}\" for k, v in curve[:4])\n      + f\"; two components carry {100 * share[1]:.0f}% of the variance\")"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "**[B-redundant]** The smallest eigenvalue is four orders of magnitude below the largest, and its direction loads on the three temperatures alone: the mean temperature of a day is the midpoint of its minimum and maximum, to within the 0.05 degree the archive records. PCA finds the redundancy without being told to look for one."
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "small = lam[-1]\nload = U[:, -1] / np.abs(U[:, -1]).max()\nnamed = {PARAMS[j]: float(load[j]) for j in range(d) if abs(load[j]) > 0.1}\nmid = (M[:, 0] + M[:, 1]) / 2.0\ngap = float(np.abs(M[:, 2] - mid).max())\nprint(f\"  smallest eigenvalue {small:.2e} against the largest {lam[0]:.2f}; \"\n      f\"its direction loads on {', '.join(named)}; the mean temperature \"\n      f\"differs from the midpoint of minimum and maximum by at most \"\n      f\"{gap:.2f} deg\")\ncheck(\"[B-redundant] the smallest eigenvalue is four orders of magnitude \"\n      \"below the largest\", small < 1e-4 * lam[0])\ncheck(\"[B-redundant] its direction loads only on the three temperatures\",\n      set(named) == {\"tlmin\", \"tlmax\", \"tl_mittel\"})\ncheck(\"[B-redundant] the mean temperature is the midpoint of the minimum and \"\n      \"the maximum, to within the recording step\",\n      gap <= 0.05 + 1e-9)\n\n# ---- CSVs for the entry's figure\nnp.savetxt(OUT_DIR / \"pca_points.csv\", T, delimiter=\",\", header=\"x1,x2\",\n           comments=\"\", fmt=\"%.2f\")\nfor j in (0, 1):\n    tip = np.sqrt(lam2[j]) * U2[:, j] * (1.0 if U2[0, j] >= 0 else -1.0)\n    with open(OUT_DIR / f\"pca_axis{j + 1}.csv\", \"w\") as f:\n        f.write(\"x1,x2\\n0.000,0.000\\n\" + f\"{tip[0]:.3f},{tip[1]:.3f}\\n\")\nwith open(OUT_DIR / \"pca_spectrum.csv\", \"w\") as f:\n    f.write(\"index,eigenvalue,share\\n\")\n    for k in range(d):\n        f.write(f\"{k + 1},{lam[k]:.4f},{share[k]:.4f}\\n\")"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "**[B-plot]** Preview: the days in the two temperatures with the two principal directions, and the eigenvalue spectrum of the eight measurements."
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(10, 4.2))\nax1.scatter(T[:, 0], T[:, 1], s=10, color=\"0.55\", label=\"day of 2024\")\nfor j, style in ((0, \"-\"), (1, \"--\")):\n    tip = np.sqrt(lam2[j]) * U2[:, j] * (1.0 if U2[0, j] >= 0 else -1.0)\n    ax1.annotate(\"\", xy=tip, xytext=(0, 0),\n                 arrowprops=dict(arrowstyle=\"->\", lw=2, ls=style, color=\"black\"))\n    ax1.annotate(f\"$u^{{({j + 1})}}$\", xy=tip * 1.12, fontsize=11)\nax1.set_aspect(\"equal\")\nax1.set_xlabel(\"centered minimum temperature in deg C\")\nax1.set_ylabel(\"centered maximum temperature in deg C\")\nax1.set_title(\"Principal directions of the 366 days, scaled by sqrt(eigenvalue)\")\nax1.legend(frameon=False, fontsize=8)\nax2.bar(np.arange(1, d + 1), lam, color=\"0.6\", edgecolor=\"black\")\nax2.plot(np.arange(1, d + 1), share * lam.max(), \"ko--\", label=\"cumulative share\")\nax2.set_xlabel(\"index of the eigenvalue\")\nax2.set_ylabel(\"eigenvalue of the sample covariance matrix\")\nax2.set_title(\"Spectrum of the eight scaled measurements\")\nax2.legend(frameon=False, fontsize=8)\nfig.tight_layout()\nfig.savefig(OUT_DIR / \"pca.png\", dpi=110)\n\nn_ok = sum(ok for _, ok in report)\nprint(f\"\\n{n_ok}/{len(report)} checks pass\")\nif n_ok != len(report):\n    raise SystemExit(1)"
  }
 ]
}