{
 "nbformat": 4,
 "nbformat_minor": 5,
 "metadata": {
  "kernelspec": {
   "name": "python3",
   "display_name": "Python 3",
   "language": "python"
  },
  "language_info": {
   "name": "python"
  },
  "colab": {
   "name": "rpca.ipynb"
  }
 },
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "# robust principal component analysis (robust PCA) \u2014 Python demo\n\nNumerical companion to the entry [robust principal component analysis (robust PCA)](https://dictionaryofml.org/terms/rpca.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 'rpca'. The hourly air temperature at Krems an der Donau in 2024 (GeoSphere Austria archive, station 3805) fills a matrix with one row per day and one column per hour, 366 x 24. Every row is a daily cycle on top of the day's level, so a few singular values carry almost all of the energy: the matrix is nearly low rank. A sensor fault is then injected into 5% of the entries (a fixed-seed choice of entries, each shifted by 20 to 40 degrees in a random direction): the corrupted matrix is the sum of a low-rank part and a sparse part.\n\nRequires NumPy and Matplotlib only, and uses fixed seeds, so the printed numbers reproduce exactly. Generated from [`pythondemos/rpca.py`](https://dictionaryofml.org/terms/rpca.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(), \"rpca.py\")\nos.makedirs(\"pythondemos\", exist_ok=True)"
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "\"\"\"Robust PCA on a year of hourly temperatures arranged as a matrix of\ndays by hours: the daily cycle makes the matrix nearly low rank, injected\nsensor faults are sparse, and principal component pursuit separates the\ntwo where PCA cannot.\n\nPurpose\n-------\nNumerical companion to the glossary entry 'rpca'.  The hourly air\ntemperature at Krems an der Donau in 2024 (GeoSphere Austria archive,\nstation 3805) fills a matrix with one row per day and one column per\nhour, 366 x 24.  Every row is a daily cycle on top of the day's level, so\na few singular values carry almost all of the energy: the matrix is\nnearly low rank.  A sensor fault is then injected into 5% of the entries\n(a fixed-seed choice of entries, each shifted by 20 to 40 degrees in a\nrandom direction): the corrupted matrix is the sum of a low-rank part\nand a sparse part.\n\nPrincipal component pursuit recovers both parts by minimizing the\nnuclear norm of the low-rank part plus a weighted l1 norm of the sparse\npart under the constraint that the two sum to the observed matrix\n(Candes, Li, Ma and Wright).  The demo solves it by the alternating\ndirection method of multipliers with singular value thresholding for\nthe low-rank part and soft thresholding for the sparse part.\n\nThe demo checks the entry's claims: the clean matrix is nearly low\nrank (three singular values carry at least 95% of its energy); the\nrecovered low-rank part is within a few percent of the clean matrix,\nwhereas the best rank-three approximation of the corrupted matrix by\nPCA is much further from it; the entries of the recovered sparse part\nat the scale of a fault are the injected faults, and removing them\nleaves a matrix within a few percent of the clean one; and the recovered\nlow-rank part has rank well below the number of hours.  The clean\nmatrix is only approximately low rank, so the pursuit also assigns the\nsmall remainder of the daily curves to the sparse part, which is why\nthe support is read at the fault scale.\n\nDeterministic: fixed seed for the faults; the readings are fetched from\nthe public archive.  Self-contained: numpy + matplotlib only.\n\nBlocks\n------\n[B-matrix]   The 366 x 24 matrix of hourly readings; check that three\n             singular values carry at least 95% of its energy.\n[B-faults]   Inject faults into 5% of the entries; the corrupted matrix.\n[B-pursuit]  Principal component pursuit by the alternating direction\n             method of multipliers: check the constraint, the recovery\n             of the clean matrix, the support of the sparse part, and\n             the rank of the low-rank part.\n[B-pca]      PCA of the corrupted matrix: check that its rank-three\n             approximation is much further from the clean matrix than\n             the recovered low-rank part.\n[B-plot]     Write one day's curves and the singular values for the\n             entry's figure, plus the preview.\n\nOutputs\n-------\nrpca_temperature.csv : hour, temperature -- the readings\nrpca_day.csv : hour, clean, corrupted, recovered, pca -- one day with an\n               injected fault: the clean reading, the corrupted entry,\n               the row of the recovered low-rank part, and the row of\n               the rank-three PCA approximation of the corrupted matrix\nrpca_singular.csv : index, clean, corrupted, recovered -- the singular\n               values of the three matrices\nrpca.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\nARCHIVE = \"https://dataset.api.hub.geosphere.at/v1/station/historical/\"\nSTATION = 3805\n\nreport = []                         # collects (check name, pass/fail) pairs\n\n\ndef check(name, ok):                # records and prints one verification\n    report.append((name, bool(ok)))\n    print(f\"  [{'ok' if ok else 'FAIL'}] {name}\")"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "**[B-matrix]** The 366 x 24 matrix of hourly readings; check that three singular values carry at least 95% of its energy."
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "def fetch(resource, parameter, start, end):\n    \"\"\"One parameter of one station from the GeoSphere Austria archive.\"\"\"\n    url = (f\"{ARCHIVE}{resource}?parameters={parameter}\"\n           f\"&station_ids={STATION}&start={start}&end={end}\")\n    with urllib.request.urlopen(url, timeout=180) as resp:\n        payload = json.load(resp)\n    values = payload[\"features\"][0][\"properties\"][\"parameters\"]\n    key = list(values.keys())[0]\n    return (np.array(values[key][\"data\"], dtype=float),\n            [t[:13] for t in payload[\"timestamps\"]])\n\n\ntemp, stamps = fetch(\"klima-v2-1h\", \"tl\", \"2024-01-01T00:00\", \"2024-12-31T23:00\")\ntemp = np.where(np.isnan(temp), np.nanmean(temp), temp)      # a few gaps\nwith open(OUT_DIR / \"rpca_temperature.csv\", \"w\") as f:\n    f.write(\"hour,temperature\\n\")\n    for h, v in zip(stamps, temp):\n        f.write(f\"{h},{v:.1f}\\n\")\nL0 = temp.reshape(-1, 24)                       # one row per day, one column per hour\nn_days, n_hours = L0.shape\nsv_clean = np.linalg.svd(L0, compute_uv=False)\nenergy3 = float((sv_clean[:3] ** 2).sum() / (sv_clean ** 2).sum())\nprint(f\"  matrix of {n_days} days x {n_hours} hours; three singular values carry \"\n      f\"{100 * energy3:.1f}% of the energy, the largest three are {np.round(sv_clean[:3], 0).tolist()}\")\ncheck(\"[B-matrix] three singular values carry at least 95% of the energy of the clean matrix\",\n      energy3 >= 0.95)"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "**[B-faults]** Inject faults into 5% of the entries; the corrupted matrix."
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "rng = np.random.default_rng(0)\nFAULT_SHARE = 0.05\nmask = rng.uniform(size=L0.shape) < FAULT_SHARE\nS0 = np.zeros_like(L0)\nS0[mask] = rng.choice([-1.0, 1.0], size=int(mask.sum())) * rng.uniform(20.0, 40.0, size=int(mask.sum()))\nMobs = L0 + S0                                  # the corrupted matrix\nprint(f\"  {int(mask.sum())} of {L0.size} entries corrupted by 20-40 degrees\")\ncheck(\"[B-faults] the faults touch about 5% of the entries\",\n      0.04 <= mask.mean() <= 0.06)"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "**[B-pursuit]** Principal component pursuit by the alternating direction method of multipliers: check the constraint, the recovery of the clean matrix, the support of the sparse part, and the rank of the low-rank part."
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "def soft(A, tau):\n    return np.sign(A) * np.maximum(np.abs(A) - tau, 0.0)\n\n\ndef svt(A, tau):\n    \"\"\"Singular value thresholding: shrink the singular values by tau.\"\"\"\n    U, s, Vt = np.linalg.svd(A, full_matrices=False)\n    return (U * np.maximum(s - tau, 0.0)) @ Vt\n\n\ndef pursuit(M, lam, iters=3000, tol=1e-7):\n    \"\"\"min ||L||_* + lam ||S||_1  s.t.  L + S = M, by ADMM.\"\"\"\n    mu = 0.25 * M.size / np.abs(M).sum()         # weight of the augmented Lagrangian, as in Candes et al.\n    S = np.zeros_like(M)\n    Y = np.zeros_like(M)                          # scaled dual variable\n    for it in range(iters):\n        L = svt(M - S + Y / mu, 1.0 / mu)\n        S = soft(M - L + Y / mu, lam / mu)\n        residual = M - L - S\n        Y = Y + mu * residual\n        if np.linalg.norm(residual) <= tol * np.linalg.norm(M):\n            break\n    return L, S, it + 1\n\n\nLAM = 1.0 / np.sqrt(max(n_days, n_hours))\nL_hat, S_hat, n_iter = pursuit(Mobs, LAM)\nrel_L = np.linalg.norm(L_hat - L0) / np.linalg.norm(L0)\nsupport_hat = np.abs(S_hat) > 10.0             # entries at the scale of a fault\ntp = int((support_hat & mask).sum())\nflagged_are_faults = tp / max(int(support_hat.sum()), 1)   # share of flagged entries that are faults\nfaults_flagged = tp / int(mask.sum())                       # share of faults that are flagged\nsv_hat = np.linalg.svd(L_hat, compute_uv=False)\nrank_hat = int((sv_hat > 1e-3 * sv_hat[0]).sum())\ncleaned = Mobs - np.where(support_hat, S_hat, 0.0)   # the faults the pursuit found, removed\nrel_cleaned = np.linalg.norm(cleaned - L0) / np.linalg.norm(L0)\nprint(f\"  pursuit with lambda = 1/sqrt({max(n_days, n_hours)}) = {LAM:.3f}: {n_iter} iterations; \"\n      f\"relative error of the low-rank part {100 * rel_L:.2f}%; sparse part: {100 * flagged_are_faults:.0f}% \"\n      f\"of the flagged entries are faults, {100 * faults_flagged:.0f}% of the faults are flagged; rank of the low-rank part {rank_hat} \"\n      f\"of {n_hours}; removing the found faults leaves a matrix within {100 * rel_cleaned:.2f}% of the clean one\")\ncheck(\"[B-pursuit] the two parts sum to the observed matrix\",\n      np.linalg.norm(L_hat + S_hat - Mobs) <= 1e-5 * np.linalg.norm(Mobs))\ncheck(\"[B-pursuit] the recovered low-rank part is within 8% of the clean matrix\",\n      rel_L < 0.08)\ncheck(\"[B-pursuit] the entries of the sparse part at the scale of a fault are the faults \"\n      \"(at least 90% of them are faults, and at least 90% of the faults are among them)\",\n      flagged_are_faults >= 0.9 and faults_flagged >= 0.9)\ncheck(\"[B-pursuit] removing the faults the pursuit found leaves a matrix within 2.5% of the \"\n      \"clean one\", rel_cleaned < 0.025)\ncheck(\"[B-pursuit] the recovered low-rank part has rank at most half the number of hours\",\n      rank_hat <= n_hours // 2)"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "**[B-pca]** PCA of the corrupted matrix: check that its rank-three approximation is much further from the clean matrix than the recovered low-rank part."
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "U, s, Vt = np.linalg.svd(Mobs, full_matrices=False)\nL_pca = (U[:, :3] * s[:3]) @ Vt[:3]              # best rank-3 approximation of Mobs\nrel_pca = np.linalg.norm(L_pca - L0) / np.linalg.norm(L0)\nsv_obs = np.linalg.svd(Mobs, compute_uv=False)\nprint(f\"  rank-3 PCA approximation of the corrupted matrix: relative error {100 * rel_pca:.2f}% \"\n      f\"against {100 * rel_L:.2f}% for the pursuit\")\ncheck(\"[B-pca] the rank-three PCA approximation of the corrupted matrix is at least \"\n      \"three times further from the clean matrix than the recovered low-rank part\",\n      rel_pca >= 3.0 * rel_L)"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "**[B-plot]** Write one day's curves and the singular values for the entry's figure, plus the preview."
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "faulty_days = np.flatnonzero(mask.sum(axis=1) >= 2)\nday = int(faulty_days[len(faulty_days) // 2])    # a day in the middle of the year with two faults\nwith open(OUT_DIR / \"rpca_day.csv\", \"w\") as f:\n    f.write(\"hour,clean,corrupted,recovered,pca\\n\")\n    for h in range(n_hours):\n        f.write(f\"{h},{L0[day, h]:.2f},{Mobs[day, h]:.2f},{L_hat[day, h]:.2f},{L_pca[day, h]:.2f}\\n\")\nwith open(OUT_DIR / \"rpca_singular.csv\", \"w\") as f:\n    f.write(\"index,clean,corrupted,recovered\\n\")\n    for i in range(n_hours):\n        f.write(f\"{i + 1},{sv_clean[i]:.3f},{sv_obs[i]:.3f},{sv_hat[i]:.3f}\\n\")\n\nfig, axes = plt.subplots(1, 3, figsize=(15, 3.6))\nax = axes[0]\nhours = np.arange(n_hours)\nax.plot(hours, Mobs[day], \"x\", color=\"0.5\", markersize=6, label=\"corrupted entries\")\nax.plot(hours, L0[day], \"-\", color=\"black\", label=\"clean readings\")\nax.plot(hours, L_hat[day], \"--\", color=\"black\", label=\"recovered low-rank part\")\nax.plot(hours, L_pca[day], \":\", color=\"0.5\", label=\"rank-3 PCA of corrupted\")\nax.set_xlabel(\"hour of the day\")\nax.set_ylabel(\"temperature (\u00b0C)\")\nax.set_title(f\"Day {day + 1} of 2024 with two injected faults\")\nax.legend(frameon=False, fontsize=8)\nax = axes[1]\nax.semilogy(np.arange(1, n_hours + 1), sv_clean, \"o-\", color=\"black\", markersize=3, label=\"clean matrix\")\nax.semilogy(np.arange(1, n_hours + 1), sv_obs, \"s--\", color=\"0.5\", markersize=3, label=\"corrupted matrix\")\nax.semilogy(np.arange(1, n_hours + 1), np.maximum(sv_hat, 1e-3), \"^:\", color=\"black\", markersize=3,\n            label=\"recovered low-rank part\")\nax.set_xlabel(\"index\")\nax.set_ylabel(\"singular value\")\nax.set_title(\"Singular values\")\nax.legend(frameon=False, fontsize=8)\nax = axes[2]\nax.imshow(np.abs(S_hat.T), aspect=\"auto\", cmap=\"gray_r\", vmin=0, vmax=40)\nax.set_xlabel(\"day of the year\")\nax.set_ylabel(\"hour of the day\")\nax.set_title(\"Recovered sparse part (absolute value)\")\nfig.tight_layout()\nfig.savefig(OUT_DIR / \"rpca.png\", dpi=150)\nplt.close(fig)\ncheck(\"[B-plot] the CSV files were written\",\n      (OUT_DIR / \"rpca_day.csv\").exists() and (OUT_DIR / \"rpca_singular.csv\").exists())\n\npassed = sum(1 for _, ok in report if ok)\nprint(f\"\\n{passed}/{len(report)} checks pass\")"
  }
 ]
}