A matrix factorization writes \(A\) as a product of simpler matrices (triangular, orthogonal, diagonal, or all-nonnegative). Each standard factorization is also the solution to an optimization problem: minimize an error subject to a structural constraint.
The sections below state each problem, then apply it to synthetic matrices and three real datasets.
import io
import time
import matplotlib.pyplot as plt
import matplotlib.ticker as mticker
import numpy as np
import pandas as pd
import pyreadr
import requests
import seaborn as sns
from numpy.linalg import cond, eigh, svd
from scipy.linalg import cho_factor, cho_solve, cholesky, qr, solve, solve_triangular
from sklearn.decomposition import NMF
sns.set_theme(style="whitegrid", rc={"axes.edgecolor": "0.85"})
ACCENT = "#4A3AA7"
GREY = "#999999"
rng = np.random.default_rng(7)Ordinary least squares minimizes
\[ \min_\beta \|y-X\beta\|_2^2. \]
Calculus gives the closed form
\[ \hat\beta=(X^TX)^{-1}X^Ty. \]
The condition number \(\kappa(A)\) measures input–output error amplification. For \(X^TX\),
\[ \kappa(X^TX)=\kappa(X)^2 \]
because eigenvalues of \(X^TX\) are squared singular values of \(X\). A solve via the normal equations amplifies floating-point error roughly as \(\kappa(X)^2\varepsilon\) instead of \(\kappa(X)\varepsilon\), where \(\varepsilon\) is machine epsilon.
Any full-column-rank \(X\) admits \(X=QR\) with orthonormal \(Q\) (\(\kappa(Q)=1\)) and upper-triangular \(R\). Substituting into the normal equations:
\[ R^TQ^TQR\beta=R^TQ^Ty \;\Rightarrow\; R\beta=Q^Ty. \]
Back-substitution on \(R\beta=Q^Ty\) costs \(O(n^2)\) and is conditioned on \(\kappa(X)\), not \(\kappa(X)^2\), because \(X^TX\) is never formed.
The demo builds \(X\) with known \(\kappa(X)=10^6\): random orthonormal \(U_0,V_0\), singular values decaying geometrically from \(1\) to \(10^{-6}\), then \(X=(U_0 \operatorname{diag}(\sigma))V_0^T\). Real design matrices become ill-conditioned from nearly dependent columns (redundant sensors, correlated markers). The response \(y=X\beta_{\text{true}}\) is noiseless so recovery error is pure floating-point conditioning.
# Build X with a controlled condition number instead of a raw high-degree
# Vandermonde basis, whose columns are so nearly dependent that beta itself
# becomes ill-posed, not just numerically delicate -- singular values decaying
# geometrically isolate the conditioning effect the derivation is about.
n, p = 200, 12
U0, _ = np.linalg.qr(rng.normal(size=(n, p)))
V0, _ = np.linalg.qr(rng.normal(size=(p, p)))
singvals = np.logspace(0, -6, p) # cond(X) ~ 1e6, cond(X.T @ X) ~ 1e12
X = (U0 * singvals) @ V0.T
beta_true = rng.normal(size=p)
y = X @ beta_true # noiseless: recovery error below is pure floating-point conditioning
print(f"cond(X) = {cond(X):.3e}")
print(f"cond(X.T@X) = {cond(X.T @ X):.3e}")
beta_normal = solve(X.T @ X, X.T @ y)
Q, R = qr(X, mode="economic")
beta_qr = solve_triangular(R, Q.T @ y)
print(f"||beta_normal - beta_true|| = {np.linalg.norm(beta_normal - beta_true):.2e}")
print(f"||beta_qr - beta_true|| = {np.linalg.norm(beta_qr - beta_true):.2e}")cond(X) = 1.000e+06
cond(X.T@X) = 1.000e+12
||beta_normal - beta_true|| = 2.70e-05
||beta_qr - beta_true|| = 5.63e-12On identical data, the normal-equations solve recovers \(\beta\) to about five decimal places; QR to about eleven.
Given \(A\), a rank-\(k\) summary minimizes Frobenius error among matrices of rank at most \(k\). The Frobenius norm is the Euclidean norm on vectorized entries.
The SVD writes
\[ A=U\Sigma V^T \]
with orthonormal \(U,V\) and diagonal \(\Sigma\) of singular values \(\sigma_1\ge\sigma_2\ge\dots\). Truncation at rank \(k\) gives
\[ A_k=\sum_{i\le k}\sigma_i u_i v_i^T. \]
For any \(B\) with \(\operatorname{rank}(B)\le k\),
\[ \|A-B\|_F \ge \sqrt{\sum_{i>k}\sigma_i^2}, \]
with equality at \(A_k\). The proof uses the Courant–Fischer min-max characterization of singular values.
A = rng.normal(size=(60, 40))
U, s, Vt = svd(A, full_matrices=False)
k = 10
A_k = (U[:, :k] * s[:k]) @ Vt[:k]
residual_frobenius = np.linalg.norm(A - A_k, "fro")
tail_bound = np.sqrt(np.sum(s[k:] ** 2))
print(f"||A - A_k||_F = {residual_frobenius:.6f}")
print(f"sqrt(sum tail^2) = {tail_bound:.6f}")||A - A_k||_F = 32.942888
sqrt(sum tail^2) = 32.942888The residual equals the tail bound to printed precision.
For symmetric \(A\), maximize projected variance subject to unit length:
\[ \max_{\|v\|=1} v^TAv. \]
The Lagrangian \(\mathcal L(v,\lambda)=v^TAv-\lambda(v^Tv-1)\) yields stationarity
\[ Av=\lambda v. \]
Eigenvectors are stationary points of the Rayleigh quotient \(v^TAv/v^Tv\); eigenvalues are the values attained there. Repeating on orthogonal complements gives \(A=Q\Lambda Q^T\).
A covariance matrix is \(\Sigma=X_c^TX_c\) for centered \(X_c\). Eigendecomposing \(\Sigma\) squares conditioning: eigenvalues of \(\Sigma\) equal squared singular values of \(X_c\). SVD of \(X_c\) returns the same spectrum at \(\kappa(X_c)\).
Xc = X - X.mean(axis=0)
Sigma = Xc.T @ Xc
eigvals, eigvecs = eigh(Sigma)
_, s_centered, _ = svd(Xc, full_matrices=False)
print(f"top eigenvalues: {np.sort(eigvals)[::-1][:4].round(3)}")
print(f"top squared singular vals: {(s_centered[:4] ** 2).round(3)}")
print(f"cond(Sigma) = {cond(Sigma):.3e} cond(Xc) = {cond(Xc):.3e}")top eigenvalues: [0.989 0.08 0.007 0.001]
top squared singular vals: [0.989 0.08 0.007 0.001]
cond(Sigma) = 9.911e+11 cond(Xc) = 9.955e+05Top eigenvalues match squared singular values; \(\kappa(\Sigma)\approx\kappa(X_c)^2\).
SVD allows cancellation via negative entries. NMF constrains both factors to be nonnegative:
\[ \min_{W,H\ge0} \|V-WH\|_F^2. \]
Lee and Seung’s multiplicative updates
\[ H \leftarrow H \odot \frac{W^TV}{W^TWH}, \qquad W \leftarrow W \odot \frac{VH^T}{WHH^T} \]
decrease the objective while preserving nonnegativity. The problem is biconvex, not jointly convex; updates converge to a stationary point depending on initialization. nndsvda initializes from the SVD for determinism, not global optimality.
V_demo = rng.gamma(shape=2.0, scale=1.0, size=(50, 30))
demo_model = NMF(n_components=5, init="nndsvda", max_iter=500, random_state=7)
W_demo = demo_model.fit_transform(V_demo)
H_demo = demo_model.components_
print(f"reconstruction error = {demo_model.reconstruction_err_:.4f}")
print(f"W, H nonnegative: {(W_demo >= 0).all()}, {(H_demo >= 0).all()}")reconstruction error = 41.9530
W, H nonnegative: True, TrueGamma-distributed random data stands in for count/intensity matrices (word counts, probe intensities).
For symmetric positive definite (SPD) \(A\) (equivalently \(x^TAx>0\) for all nonzero \(x\)),
\[ A=LL^T \]
with lower-triangular \(L\) and positive diagonal. Factorization costs \(O(n^3/3)\); each subsequent solve is two triangular substitutions at \(O(n^2)\). For \(m\) right-hand sides,
\[ O(n^3+mn^2). \]
Explicit inversion has similar asymptotics and worse conditioning.
Sampling: if \(z\sim\mathcal N(0,I)\) then \(x=Lz\sim\mathcal N(0,\Sigma)\) when \(A=\Sigma\).
sizes = [50, 100, 200, 400, 800]
factor_times, inv_times = [], []
for size in sizes:
Rm = rng.normal(size=(size, size))
Am = Rm @ Rm.T + size * np.eye(size)
Bm = rng.normal(size=(size, 5))
t0 = time.perf_counter()
c, low = cho_factor(Am)
cho_solve((c, low), Bm)
factor_times.append(time.perf_counter() - t0)
t0 = time.perf_counter()
np.linalg.inv(Am) @ Bm
inv_times.append(time.perf_counter() - t0)
fig, ax = plt.subplots(figsize=(6, 4))
ax.loglog(sizes, factor_times, "o-", color=ACCENT, label="cho_factor + cho_solve")
ax.loglog(sizes, inv_times, "o-", color=GREY, label="explicit inverse")
ax.xaxis.set_major_locator(mticker.FixedLocator(sizes))
ax.xaxis.set_major_formatter(mticker.FixedFormatter([str(s) for s in sizes]))
ax.xaxis.set_minor_locator(mticker.NullLocator())
ax.set_xlabel("matrix size n")
ax.set_ylabel("wall-clock seconds")
ax.set_title("Factor once, solve many: Cholesky vs. explicit inverse")
ax.legend()
plt.show()
realAWSCloudwatch category.(This section fetches data at render time; project renders respect freeze: auto.)
base = "https://raw.githubusercontent.com/numenta/NAB/master"
ec2_ids = ["53ea38", "24ae8d", "5f5533", "77c1ca", "825cc2", "ac20cd", "c6585a", "fe7f93"]
labels = requests.get(f"{base}/labels/combined_labels.json", timeout=30).json()
columns, anomaly_idx = {}, {}
for id_ in ec2_ids:
key = f"realAWSCloudwatch/ec2_cpu_utilization_{id_}.csv"
df = pd.read_csv(io.StringIO(requests.get(f"{base}/data/{key}", timeout=30).text))
df["timestamp"] = pd.to_datetime(df["timestamp"])
columns[id_] = df["value"].to_numpy()
label_times = pd.to_datetime(labels.get(key, []))
anomaly_idx[id_] = df.index[df["timestamp"].isin(label_times)].to_numpy()
# stacked by position, not wall-clock: fleet-wide regime structure, not calendar sync
telemetry = pd.DataFrame(columns)
print(f"real telemetry matrix: {telemetry.shape[0]} timesteps x {telemetry.shape[1]} real EC2 instances")real telemetry matrix: 4032 timesteps x 8 real EC2 instancesColumns are stacked by timestep index (not wall-clock alignment).
X_tel = telemetry.to_numpy()
Xc_tel = X_tel - X_tel.mean(axis=0)
U_tel, s_tel, Vt_tel = svd(Xc_tel, full_matrices=False)
fig, ax = plt.subplots(figsize=(6, 4))
ax.semilogy(np.arange(1, len(s_tel) + 1), s_tel, "o-", color=ACCENT)
ax.set_xlabel("component")
ax.set_ylabel("singular value (log scale)")
ax.set_title("Singular values of 8 real EC2 CPU-utilization traces")
plt.show()
fig, ax = plt.subplots(figsize=(8, 2.8))
sns.heatmap(U_tel[:, :5].T, cmap="RdBu_r", center=0, cbar_kws={"label": "loading"}, ax=ax)
ax.set_xlabel("timestep")
ax.set_ylabel("left singular vector")
ax.set_title("Top 5 temporal regime vectors across the real fleet")
plt.show()
k_tel = 4
X_k_tel = (U_tel[:, :k_tel] * s_tel[:k_tel]) @ Vt_tel[:k_tel]
residual = np.abs(Xc_tel - X_k_tel)
target = "fe7f93"
j = ec2_ids.index(target)
fig, ax = plt.subplots(figsize=(8, 4))
ax.plot(residual[:, j], color=ACCENT, lw=1)
for idx in anomaly_idx[target]:
ax.axvline(idx, color="crimson", ls="--", alpha=0.7)
ax.set_title(f"Rank-{k_tel} residual, instance {target} (dashed = real labeled anomaly)")
ax.set_xlabel("timestep")
ax.set_ylabel("|residual|")
plt.show()
Top left singular vectors encode shared fleet regimes; rank-\(k\) residual energy peaks near human-labelled anomalies.
wheat dataset in BGLR-R; CIMMYT Global Wheat Program.Linkage disequilibrium — correlated inheritance blocks — is collinearity under a genetic name.
resp = requests.get("https://raw.githubusercontent.com/gdlc/BGLR-R/master/data/wheat.RData", timeout=60)
with open("/tmp/wheat.RData", "wb") as f:
f.write(resp.content)
wheat = pyreadr.read_r("/tmp/wheat.RData")
X_markers = wheat["wheat.X"].to_numpy()
grain_yield = wheat["wheat.Y"].iloc[:, 0].to_numpy() # real yield, environment 1
print(f"real markers: {X_markers.shape[0]} wheat lines x {X_markers.shape[1]} real DArT markers")real markers: 599 wheat lines x 1279 real DArT markersmarker_var = X_markers.var(axis=0)
ranked_markers = np.argsort(marker_var)[::-1]
n_lines = X_markers.shape[0]
perm = rng.permutation(n_lines)
train_idx, test_idx = perm[:450], perm[450:]
y_train, y_test = grain_yield[train_idx], grain_yield[test_idx]
marker_counts = [50, 150, 250, 350, 400, 430, 449]
corr_qr, cond_x, cond_xtx = [], [], []
for p_count in marker_counts:
Xp = X_markers[:, ranked_markers[:p_count]].astype(float)
Xp -= Xp.mean(axis=0)
X_train, X_test = Xp[train_idx], Xp[test_idx]
Q_tr, R_tr = qr(X_train, mode="economic")
beta_qr = solve_triangular(R_tr, Q_tr.T @ y_train)
corr_qr.append(np.corrcoef(X_test @ beta_qr, y_test)[0, 1])
cond_x.append(cond(X_train))
cond_xtx.append(cond(X_train.T @ X_train))
print("markers: ", marker_counts)
print("held-out corr:", [f"{c:.3f}" for c in corr_qr])
print("cond(X): ", [f"{c:.2e}" for c in cond_x])
print("cond(X.T@X): ", [f"{c:.2e}" for c in cond_xtx])
X_sub = X_markers[:, ranked_markers[:350]].astype(float)
X_sub -= X_sub.mean(axis=0)markers: [50, 150, 250, 350, 400, 430, 449]
held-out corr: ['0.200', '0.267', '0.181', '0.090', '0.106', '0.114', '0.050']
cond(X): ['2.28e+01', '8.56e+01', '2.37e+02', '5.25e+02', '1.13e+03', '2.21e+03', '1.89e+04']
cond(X.T@X): ['5.19e+02', '7.33e+03', '5.62e+04', '2.75e+05', '1.27e+06', '4.88e+06', '3.57e+08']fig, axes = plt.subplots(1, 3, figsize=(14, 4))
Q20, R20 = qr(X_sub[:, :20], mode="economic")
sns.heatmap(np.triu(R20), cmap="RdBu_r", center=0, ax=axes[0], cbar_kws={"label": "R value"})
axes[0].set_title("Upper-triangular R\n(20 real markers)")
axes[1].plot(marker_counts, corr_qr, "o-", color=ACCENT)
axes[1].axhline(0, color=GREY, lw=1)
axes[1].set_xlabel("markers used (of 450 training lines)")
axes[1].set_ylabel("held-out yield correlation")
axes[1].set_title("Plain OLS overfits as p -> n")
axes[2].semilogy(marker_counts, cond_x, "o-", color=ACCENT, label="cond(X)")
axes[2].semilogy(marker_counts, cond_xtx, "o-", color=GREY, label="cond(X'X)")
axes[2].axhline(1 / np.finfo(float).eps, color="crimson", ls="--", lw=1, label="float64 ceiling")
axes[2].set_xlabel("markers used")
axes[2].set_ylabel("condition number")
axes[2].set_title("Normal eqs burns precision twice as fast")
axes[2].legend(fontsize=8)
plt.tight_layout()
plt.show()
Held-out correlation peaks near 150 markers then falls as \(p\to n\). \(\kappa(X^TX)\approx\kappa(X)^2\) at every marker count. Normal equations lose precision near \(\kappa(X)\sim10^8\); QR remains usable to \(\kappa(X)\sim10^{16}\).
Genomic relationship matrix \(G\propto XX^T\) is SPD; Cholesky samples correlated breeding values:
G = (X_sub @ X_sub.T) / X_sub.shape[1]
G += 1e-6 * np.eye(G.shape[0]) # ridge for numerical PD
L = cholesky(G, lower=True)
breeding_values = L @ rng.normal(size=G.shape[0])
print(f"sampled breeding values: mean={breeding_values.mean():.3f}, var={breeding_values.var():.3f}")sampled breeding values: mean=0.000, var=0.137Leading eigenvectors of \(G\) encode population structure (EIGENSTRAT-style covariates in GWAS):
eigvals_G, eigvecs_G = eigh(G)
order = np.argsort(eigvals_G)[::-1]
pc1, pc2 = eigvecs_G[:, order[0]], eigvecs_G[:, order[1]]
fig, ax = plt.subplots(figsize=(5, 5))
ax.scatter(pc1, pc2, s=16, color=ACCENT, alpha=0.7)
ax.set_xlabel("PC1 (top eigenvector of G)")
ax.set_ylabel("PC2")
ax.set_title("Population structure, 599 real wheat lines")
plt.show()
wheat.sets is a cross-validation fold, not a population label; points are uncoloured.
gbase = "https://raw.githubusercontent.com/dharsandip/Classification_of_Cancer_by_Gene_Expression_Dataset/master"
train_raw = pd.read_csv(f"{gbase}/data_set_ALL_AML_train.csv")
indep_raw = pd.read_csv(f"{gbase}/data_set_ALL_AML_independent.csv")
actual = pd.read_csv(f"{gbase}/actual.csv")
def expression_matrix(raw):
id_cols = [c for c in raw.columns if not c.startswith("call") and c not in ("Gene Description", "Gene Accession Number")]
return raw[id_cols].to_numpy(dtype=float).T, [int(c) for c in id_cols]
X_train_expr, train_patients = expression_matrix(train_raw)
X_test_expr, test_patients = expression_matrix(indep_raw)
expr = np.vstack([X_train_expr, X_test_expr]) # 72 real patients x 7129 real genes
patient_order = train_patients + test_patients
diagnosis = actual.set_index("patient").loc[patient_order, "cancer"].to_numpy()
print(f"real expression matrix: {expr.shape[0]} patients x {expr.shape[1]} genes "
f"({(diagnosis == 'ALL').sum()} ALL / {(diagnosis == 'AML').sum()} AML, real diagnoses)")real expression matrix: 72 patients x 7129 genes (47 ALL / 25 AML, real diagnoses)expr_log = np.log2(np.clip(expr, 1.0, None))
gene_var = expr_log.var(axis=0)
top_genes = np.argsort(gene_var)[::-1][:1500]
V_leuk = expr_log[:, top_genes]
leuk_model = NMF(n_components=2, init="nndsvda", max_iter=1000, random_state=7)
W_leuk = leuk_model.fit_transform(V_leuk)
H_leuk = leuk_model.components_
dominant = np.argmax(W_leuk, axis=1)
print(pd.crosstab(diagnosis, dominant, rownames=["real diagnosis"], colnames=["dominant metagene"]))dominant metagene 0 1
real diagnosis
ALL 26 21
AML 25 0fig, axes = plt.subplots(1, 2, figsize=(11, 4.5))
colors = np.where(diagnosis == "ALL", ACCENT, "#E07B39")
axes[0].scatter(W_leuk[:, 0], W_leuk[:, 1], c=colors, s=30, alpha=0.85)
axes[0].set_xlabel("metagene 1 weight")
axes[0].set_ylabel("metagene 2 weight")
axes[0].set_title("NMF weights, colored by real diagnosis")
sns.heatmap(H_leuk[:, :150], cmap="viridis", cbar_kws={"label": "gene weight"}, ax=axes[1])
axes[1].set_xlabel("gene (top-variance subset, first 150 shown)")
axes[1].set_ylabel("metagene")
axes[1].set_title("Learned H: two real metagenes")
plt.tight_layout()
plt.show()
One metagene dominates 21 patients, all ALL; the other splits the remainder without clean ALL/AML separation. NMF isolates a pure-ALL subgroup, not a full binary partition.
| Need | Factorization | Constraint traded |
|---|---|---|
| Stable least squares | QR | Avoid \(X^TX\) |
| Best rank-\(k\) compression | Truncated SVD | Exact fit |
| Symmetric spectrum / PCA | Eigendecomposition or SVD of \(X_c\) | Second orthogonal basis |
| Interpretable parts | NMF | Global optimum |
| Repeated SPD solves / sampling | Cholesky | Explicit inverse |
wheat dataset — CIMMYT genomic selection panel.