B-plot writes when the script runsDictionary of Applied Machine Learning · robust principal component analysis (robust PCA)
Numerical companion to the entry robust principal component analysis (robust PCA): it recomputes what the entry states and prints one line per check
Numerical 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.
Run it without installing anything:uv run https://dictionaryofml.org/terms/rpca.py
uv downloads this script and the pinned NumPy and Matplotlib it needs, then runs it; the script fetches any input file it uses. To keep the output files, download rpca.py into a folder and run uv run rpca.py there. With NumPy and Matplotlib already installed, python3 rpca.py, from any directory — it writes its output files into the current directory. Fixed seeds, so the printed numbers reproduce exactly. Download rpca.py · Notebook · Open in Colab
One cell per block of the script: the code, and what that code printed when it last ran here
"""Robust PCA on a year of hourly temperatures arranged as a matrix of
days by hours: the daily cycle makes the matrix nearly low rank, injected
sensor faults are sparse, and principal component pursuit separates the
two where PCA cannot.
Purpose
-------
Numerical 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.
Principal component pursuit recovers both parts by minimizing the
nuclear norm of the low-rank part plus a weighted l1 norm of the sparse
part under the constraint that the two sum to the observed matrix
(Candes, Li, Ma and Wright). The demo solves it by the alternating
direction method of multipliers with singular value thresholding for
the low-rank part and soft thresholding for the sparse part.
The demo checks the entry's claims: the clean matrix is nearly low
rank (three singular values carry at least 95% of its energy); the
recovered low-rank part is within a few percent of the clean matrix,
whereas the best rank-three approximation of the corrupted matrix by
PCA is much further from it; the entries of the recovered sparse part
at the scale of a fault are the injected faults, and removing them
leaves a matrix within a few percent of the clean one; and the recovered
low-rank part has rank well below the number of hours. The clean
matrix is only approximately low rank, so the pursuit also assigns the
small remainder of the daily curves to the sparse part, which is why
the support is read at the fault scale.
Deterministic: fixed seed for the faults; the readings are fetched from
the public archive. Self-contained: numpy + matplotlib only.
Blocks
------
[B-matrix] The 366 x 24 matrix of hourly readings; check that three
singular values carry at least 95% of its energy.
[B-faults] Inject faults into 5% of the entries; the corrupted matrix.
[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.
[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.
[B-plot] Write one day's curves and the singular values for the
entry's figure, plus the preview.
Outputs
-------
rpca_temperature.csv : hour, temperature -- the readings
rpca_day.csv : hour, clean, corrupted, recovered, pca -- one day with an
injected fault: the clean reading, the corrupted entry,
the row of the recovered low-rank part, and the row of
the rank-three PCA approximation of the corrupted matrix
rpca_singular.csv : index, clean, corrupted, recovered -- the singular
values of the three matrices
rpca.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
ARCHIVE = "https://dataset.api.hub.geosphere.at/v1/station/historical/"
STATION = 3805
report = [] # collects (check name, pass/fail) pairs
def check(name, ok): # records and prints one verification
report.append((name, bool(ok)))
print(f" [{'ok' if ok else 'FAIL'}] {name}")
matrix of 366 days x 24 hours; three singular values carry 99.5% of the energy, the largest three are [1398.0, 158.0, 91.0]
The 366 x 24 matrix of hourly readings; check that three singular values carry at least 95% of its energy.
def fetch(resource, parameter, start, end):
"""One parameter of one station from the GeoSphere Austria archive."""
url = (f"{ARCHIVE}{resource}?parameters={parameter}"
f"&station_ids={STATION}&start={start}&end={end}")
with urllib.request.urlopen(url, timeout=180) as resp:
payload = json.load(resp)
values = payload["features"][0]["properties"]["parameters"]
key = list(values.keys())[0]
return (np.array(values[key]["data"], dtype=float),
[t[:13] for t in payload["timestamps"]])
temp, stamps = fetch("klima-v2-1h", "tl", "2024-01-01T00:00", "2024-12-31T23:00")
temp = np.where(np.isnan(temp), np.nanmean(temp), temp) # a few gaps
with open(OUT_DIR / "rpca_temperature.csv", "w") as f:
f.write("hour,temperature\n")
for h, v in zip(stamps, temp):
f.write(f"{h},{v:.1f}\n")
L0 = temp.reshape(-1, 24) # one row per day, one column per hour
n_days, n_hours = L0.shape
sv_clean = np.linalg.svd(L0, compute_uv=False)
energy3 = float((sv_clean[:3] ** 2).sum() / (sv_clean ** 2).sum())
print(f" matrix of {n_days} days x {n_hours} hours; three singular values carry "
f"{100 * energy3:.1f}% of the energy, the largest three are {np.round(sv_clean[:3], 0).tolist()}")
check("[B-matrix] three singular values carry at least 95% of the energy of the clean matrix",
energy3 >= 0.95)
[ok] [B-matrix] three singular values carry at least 95% of the energy of the clean matrix 441 of 8784 entries corrupted by 20-40 degrees
Inject faults into 5% of the entries; the corrupted matrix.
rng = np.random.default_rng(0)
FAULT_SHARE = 0.05
mask = rng.uniform(size=L0.shape) < FAULT_SHARE
S0 = np.zeros_like(L0)
S0[mask] = rng.choice([-1.0, 1.0], size=int(mask.sum())) * rng.uniform(20.0, 40.0, size=int(mask.sum()))
Mobs = L0 + S0 # the corrupted matrix
print(f" {int(mask.sum())} of {L0.size} entries corrupted by 20-40 degrees")
check("[B-faults] the faults touch about 5% of the entries",
0.04 <= mask.mean() <= 0.06)
[ok] [B-faults] the faults touch about 5% of the entries pursuit with lambda = 1/sqrt(366) = 0.052: 3000 iterations; relative error of the low-rank part 5.89%; sparse part: 100% of the flagged entries are faults, 100% of the faults are flagged; rank of the low-rank part 9 of 24; removing the found faults leaves a matrix within 1.66% of the clean one
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.
def soft(A, tau):
return np.sign(A) * np.maximum(np.abs(A) - tau, 0.0)
def svt(A, tau):
"""Singular value thresholding: shrink the singular values by tau."""
U, s, Vt = np.linalg.svd(A, full_matrices=False)
return (U * np.maximum(s - tau, 0.0)) @ Vt
def pursuit(M, lam, iters=3000, tol=1e-7):
"""min ||L||_* + lam ||S||_1 s.t. L + S = M, by ADMM."""
mu = 0.25 * M.size / np.abs(M).sum() # weight of the augmented Lagrangian, as in Candes et al.
S = np.zeros_like(M)
Y = np.zeros_like(M) # scaled dual variable
for it in range(iters):
L = svt(M - S + Y / mu, 1.0 / mu)
S = soft(M - L + Y / mu, lam / mu)
residual = M - L - S
Y = Y + mu * residual
if np.linalg.norm(residual) <= tol * np.linalg.norm(M):
break
return L, S, it + 1
LAM = 1.0 / np.sqrt(max(n_days, n_hours))
L_hat, S_hat, n_iter = pursuit(Mobs, LAM)
rel_L = np.linalg.norm(L_hat - L0) / np.linalg.norm(L0)
support_hat = np.abs(S_hat) > 10.0 # entries at the scale of a fault
tp = int((support_hat & mask).sum())
flagged_are_faults = tp / max(int(support_hat.sum()), 1) # share of flagged entries that are faults
faults_flagged = tp / int(mask.sum()) # share of faults that are flagged
sv_hat = np.linalg.svd(L_hat, compute_uv=False)
rank_hat = int((sv_hat > 1e-3 * sv_hat[0]).sum())
cleaned = Mobs - np.where(support_hat, S_hat, 0.0) # the faults the pursuit found, removed
rel_cleaned = np.linalg.norm(cleaned - L0) / np.linalg.norm(L0)
print(f" pursuit with lambda = 1/sqrt({max(n_days, n_hours)}) = {LAM:.3f}: {n_iter} iterations; "
f"relative error of the low-rank part {100 * rel_L:.2f}%; sparse part: {100 * flagged_are_faults:.0f}% "
f"of the flagged entries are faults, {100 * faults_flagged:.0f}% of the faults are flagged; rank of the low-rank part {rank_hat} "
f"of {n_hours}; removing the found faults leaves a matrix within {100 * rel_cleaned:.2f}% of the clean one")
check("[B-pursuit] the two parts sum to the observed matrix",
np.linalg.norm(L_hat + S_hat - Mobs) <= 1e-5 * np.linalg.norm(Mobs))
check("[B-pursuit] the recovered low-rank part is within 8% of the clean matrix",
rel_L < 0.08)
check("[B-pursuit] the entries of the sparse part at the scale of a fault are the faults "
"(at least 90% of them are faults, and at least 90% of the faults are among them)",
flagged_are_faults >= 0.9 and faults_flagged >= 0.9)
check("[B-pursuit] removing the faults the pursuit found leaves a matrix within 2.5% of the "
"clean one", rel_cleaned < 0.025)
check("[B-pursuit] the recovered low-rank part has rank at most half the number of hours",
rank_hat <= n_hours // 2)
[ok] [B-pursuit] the two parts sum to the observed matrix [ok] [B-pursuit] the recovered low-rank part is within 8% of the clean matrix [ok] [B-pursuit] the entries of the sparse part at the scale of a fault are the faults (at least 90% of them are faults, and at least 90% of the faults are among them) [ok] [B-pursuit] removing the faults the pursuit found leaves a matrix within 2.5% of the clean one [ok] [B-pursuit] the recovered low-rank part has rank at most half the number of hours rank-3 PCA approximation of the corrupted matrix: relative error 19.21% against 5.89% for the pursuit
PCA of the corrupted matrix: check that its rank-three approximation is much further from the clean matrix than the recovered low-rank part.
U, s, Vt = np.linalg.svd(Mobs, full_matrices=False)
L_pca = (U[:, :3] * s[:3]) @ Vt[:3] # best rank-3 approximation of Mobs
rel_pca = np.linalg.norm(L_pca - L0) / np.linalg.norm(L0)
sv_obs = np.linalg.svd(Mobs, compute_uv=False)
print(f" rank-3 PCA approximation of the corrupted matrix: relative error {100 * rel_pca:.2f}% "
f"against {100 * rel_L:.2f}% for the pursuit")
check("[B-pca] the rank-three PCA approximation of the corrupted matrix is at least "
"three times further from the clean matrix than the recovered low-rank part",
rel_pca >= 3.0 * rel_L)
[ok] [B-pca] the rank-three PCA approximation of the corrupted matrix is at least three times further from the clean matrix than the recovered low-rank part
Write one day's curves and the singular values for the entry's figure, plus the preview.
faulty_days = np.flatnonzero(mask.sum(axis=1) >= 2)
day = int(faulty_days[len(faulty_days) // 2]) # a day in the middle of the year with two faults
with open(OUT_DIR / "rpca_day.csv", "w") as f:
f.write("hour,clean,corrupted,recovered,pca\n")
for h in range(n_hours):
f.write(f"{h},{L0[day, h]:.2f},{Mobs[day, h]:.2f},{L_hat[day, h]:.2f},{L_pca[day, h]:.2f}\n")
with open(OUT_DIR / "rpca_singular.csv", "w") as f:
f.write("index,clean,corrupted,recovered\n")
for i in range(n_hours):
f.write(f"{i + 1},{sv_clean[i]:.3f},{sv_obs[i]:.3f},{sv_hat[i]:.3f}\n")
fig, axes = plt.subplots(1, 3, figsize=(15, 3.6))
ax = axes[0]
hours = np.arange(n_hours)
ax.plot(hours, Mobs[day], "x", color="0.5", markersize=6, label="corrupted entries")
ax.plot(hours, L0[day], "-", color="black", label="clean readings")
ax.plot(hours, L_hat[day], "--", color="black", label="recovered low-rank part")
ax.plot(hours, L_pca[day], ":", color="0.5", label="rank-3 PCA of corrupted")
ax.set_xlabel("hour of the day")
ax.set_ylabel("temperature (°C)")
ax.set_title(f"Day {day + 1} of 2024 with two injected faults")
ax.legend(frameon=False, fontsize=8)
ax = axes[1]
ax.semilogy(np.arange(1, n_hours + 1), sv_clean, "o-", color="black", markersize=3, label="clean matrix")
ax.semilogy(np.arange(1, n_hours + 1), sv_obs, "s--", color="0.5", markersize=3, label="corrupted matrix")
ax.semilogy(np.arange(1, n_hours + 1), np.maximum(sv_hat, 1e-3), "^:", color="black", markersize=3,
label="recovered low-rank part")
ax.set_xlabel("index")
ax.set_ylabel("singular value")
ax.set_title("Singular values")
ax.legend(frameon=False, fontsize=8)
ax = axes[2]
ax.imshow(np.abs(S_hat.T), aspect="auto", cmap="gray_r", vmin=0, vmax=40)
ax.set_xlabel("day of the year")
ax.set_ylabel("hour of the day")
ax.set_title("Recovered sparse part (absolute value)")
fig.tight_layout()
fig.savefig(OUT_DIR / "rpca.png", dpi=150)
plt.close(fig)
check("[B-plot] the CSV files were written",
(OUT_DIR / "rpca_day.csv").exists() and (OUT_DIR / "rpca_singular.csv").exists())
passed = sum(1 for _, ok in report if ok)
print(f"\n{passed}/{len(report)} checks pass")
[ok] [B-plot] the CSV files were written 9/9 checks pass
B-plot writes when the script runs