""" interpretableml.py — numerical companion to the glossary entry 'interpretable machine learning (interpretable ML)'. Purpose ------- Backs the entry's two empirical claims on synthetic tabular data: restricting the hypothesis space to interpretable (sparse) hypotheses acts as regularization, and the resulting interpretable hypothesis achieves accuracy comparable to an opaque (dense) one. Also verifies that the entry's triage decision tree is an explicit, traceable function. Self-contained (numpy only), fixed seeds. Blocks ------ [B-reg] Regularization by pruning the model: on m = 30 training points with d = 20 features (3 informative), the dense least-squares hypothesis attains lower TRAINING error than the sparse one (it also fits the 17 noise features), but its TEST error is higher — restricting the hypothesis space to sparse, interpretable hypotheses reduces overfitting. [B-acc] Accuracy parity: the sparse hypothesis uses at most 6 of the 20 features (a human can read its few terms) and its test error is at most that of the dense hypothesis. [B-tree] The entry's triage tree is its own explanation: the tree prediction (temperature > 39 -> urgent; else heart rate > 120 -> urgent; else routine) coincides with the explicit two-rule lookup on a grid of 10000 (temperature, heart rate) pairs, and it partitions the feature space into exactly three axis-aligned cells. [B-effect] The two effects an interpretable hypothesis makes traceable, both exact for a linear model fitted by the normal equations: moving feature j by Delta moves the prediction by w_j * Delta, and perturbing training label r by delta shifts a prediction by delta * x^T (X^T X)^{-1} x^(r) -- computable without refitting, which refitting then confirms. [B-lime] LIME idea: a proximity-weighted linear fit approximates a nonlinear learned hypothesis near the data point x0 = 0.8 (max error < 0.3 within |x - x0| < 0.3) but not globally (error more than 4x larger far from x0). In d = 5 features, the local linear surrogate of a hypothesis that depends only on the first two features has near-zero coefficients for the other three: the surrogate depends on few features. [B-eerm] EERM idea (Zhang et al., 2024): the penalty charges the departure of the learned hypothesis from the predictions one user supplies. Those predictions are not elicited one data point at a time here; they are the predictions of a simpler proxy model the user considers interpretable, a least-squares line. Adding lambda * mean((h(x) - userpred(x))^2) to least-squares training of a degree-10 polynomial pulls the learned hypothesis toward that line; the EERM fit has (i) smaller departure from the user predictions and (ii) smaller test error than the unregularized fit, at slightly larger training error. The sparse hypothesis is the relaxed Lasso under an orthonormal design: the Lasso there has the closed form of soft-thresholding the least-squares coefficients (Hastie et al., 2009, Sect. 3.4.2), it is used to select the features, and least squares is refit on the selected support; no iterative solver is needed. Outputs ------- interpretableml_lime.csv : x, h (learned hypothesis), g (local linear surrogate; nan outside the window |x - x0| <= 0.6) for the LIME figure. interpretableml_eerm.csv : x, hopaque (unregularized polynomial), hsurr (the user predictions, from the proxy line), heerm (EERM-regularized polynomial) on a grid, for the EERM figure. interpretableml_eermdata.csv : x, y training points of the EERM figure. (The decision-tree figure is schematic TikZ and needs no data.) """ import numpy as np # the only dependency 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}") rng = np.random.default_rng(4) # synthetic tabular regression: d = 20 features, only 3 informative d, m_train = 20, 30 w_true = np.zeros(d) w_true[:3] = (3.0, -2.0, 1.5) # the interpretable ground truth sigma = 0.5 # label noise level # orthonormal design: X^T X = I, so least squares and Lasso are closed form X, _ = np.linalg.qr(rng.standard_normal((m_train, d))) y = X @ w_true + sigma * rng.standard_normal(m_train) w_dense = X.T @ y # dense least-squares hypothesis lam = 1.0 # Lasso penalty weight # Lasso under orthonormal design = soft-thresholding of w_dense; the # Lasso here SELECTS the features, and the sparse hypothesis refits # least squares on the selected support (relaxed Lasso) — under the # orthonormal design that refit just keeps the selected entries of # w_dense support = np.abs(w_dense) > lam w_sparse = np.where(support, w_dense, 0.0) # training errors: dense fits the noise features, so it fits better train_dense = np.mean((y - X @ w_dense) ** 2) train_sparse = np.mean((y - X @ w_sparse) ** 2) # expected test error for x ~ N(0, I): ||w - w_true||^2 + sigma^2 test_dense = np.sum((w_dense - w_true) ** 2) + sigma ** 2 test_sparse = np.sum((w_sparse - w_true) ** 2) + sigma ** 2 # ----------------------------------------------------------------- [B-reg] check("[B-reg] dense fits training data better but generalizes worse", train_dense < train_sparse and test_sparse < test_dense) # ----------------------------------------------------------------- [B-acc] n_used = int(np.count_nonzero(w_sparse)) check("[B-acc] sparse hypothesis uses few features at comparable " f"test error ({n_used} of {d})", n_used <= 6 and test_sparse <= test_dense) # ---------------------------------------------------------------- [B-tree] def tree(temp, hr): # the entry's triage decision tree if temp > 39.0: return "urgent" if hr > 120.0: return "urgent" return "routine" def two_rules(temp, hr): # the explicit lookup a human reads off return "urgent" if (temp > 39.0 or hr > 120.0) else "routine" temps = np.linspace(35.0, 42.0, 100) hrs = np.linspace(50.0, 180.0, 100) agree = all(tree(t, h) == two_rules(t, h) for t in temps for h in hrs) # the tree partitions the plane into 3 axis-aligned cells: # {temp > 39}, {temp <= 39, hr > 120}, {temp <= 39, hr <= 120} cells = {(t > 39.0, t <= 39.0 and h > 120.0) for t in temps for h in hrs} check("[B-tree] tree prediction equals the two-rule lookup on a " "100 x 100 grid; 3 cells", agree and len(cells) == 3) # -------------------------------------------------------------- [B-effect] # The two effects the entry claims a user of a linear model can work out: # that of a change to one feature, and that of a change to one training # label. # Both are exact here, not approximate, which is what "traceable" means. rng_s = np.random.default_rng(11) m_s, d_s = 12, 3 Xs = rng_s.standard_normal((m_s, d_s)) ys = Xs @ np.array([1.5, -0.8, 0.3]) + 0.2 * rng_s.standard_normal(m_s) G = Xs.T @ Xs # the normal equations w_hat = np.linalg.solve(G, Xs.T @ ys) # (i) feature change: moving feature j by Delta moves the prediction by # exactly w_j * Delta, whatever the other features hold x0 = rng_s.standard_normal(d_s) j, delta_x = 1, 0.7 x1 = x0.copy() x1[j] += delta_x moved = float(x1 @ w_hat - x0 @ w_hat) print(f" feature {j} moved by {delta_x}: prediction moves {moved:+.4f}, " f"weight times step {w_hat[j] * delta_x:+.4f}") check("[B-effect] a feature change moves the prediction by weight times step", abs(moved - w_hat[j] * delta_x) < 1e-12) # (ii) label change: the normal equations make the map from training set to # hypothesis LINEAR in the labels, so perturbing label r by delta shifts a # prediction by exactly delta * x^T (X^T X)^{-1} x^(r). The user can read # that off without refitting; refitting here only confirms it. r, delta_y = 4, 0.5 predicted = float(delta_y * x0 @ np.linalg.solve(G, Xs[r])) ys_pert = ys.copy() ys_pert[r] += delta_y w_pert = np.linalg.solve(G, Xs.T @ ys_pert) actual = float(x0 @ w_pert - x0 @ w_hat) print(f" label {r} perturbed by {delta_y}: prediction shifts {actual:+.4f}, " f"normal equations predict {predicted:+.4f}") check("[B-effect] a perturbed training label shifts the prediction by the " "amount the normal equations predict", abs(actual - predicted) < 1e-12) # ---------------------------------------------------------------- [B-lime] # LIME idea: a proximity-weighted linear fit approximates the learned # hypothesis near the data point x0, using few features. rng_l = np.random.default_rng(5) def h_hat(x): # stand-in for an opaque learned hypothesis return np.sin(3.0 * x) + 0.5 * x ** 2 x0 = 0.8 # the data point to be explained xs = np.linspace(-2.0, 2.0, 201) wts = np.exp(-((xs - x0) ** 2) / (2 * 0.25 ** 2)) # proximity weights at x0 A = np.c_[np.ones_like(xs), xs] # linear design (1, x) sw = np.sqrt(wts) coef, *_ = np.linalg.lstsq(A * sw[:, None], h_hat(xs) * sw, rcond=None) g = A @ coef # local linear surrogate near = np.abs(xs - x0) < 0.3 # around x0 ... far = np.abs(xs + 1.5) < 0.3 # ... vs far away err_near = np.max(np.abs(h_hat(xs) - g)[near]) err_far = np.max(np.abs(h_hat(xs) - g)[far]) check("[B-lime] weighted linear fit matches h near x0, not globally", err_near < 0.3 and err_far > 4.0 * err_near) # few features: in d = 5, a hypothesis that depends on features 0 and 1 # only yields a local surrogate with near-zero remaining coefficients Z = rng_l.standard_normal((400, 5)) * 0.3 # perturbations around a point z0 = np.array([0.8, -0.2, 0.5, 0.1, -0.4]) Zp = z0 + Z hz = np.sin(3.0 * Zp[:, 0]) + 0.5 * Zp[:, 1] ** 2 # depends on features 0, 1 wz = np.exp(-np.sum(Z ** 2, axis=1) / (2 * 0.3 ** 2)) Az = np.c_[np.ones(len(Zp)), Zp] swz = np.sqrt(wz) cz, *_ = np.linalg.lstsq(Az * swz[:, None], hz * swz, rcond=None) used = np.abs(cz[1:]) > 0.1 # features the surrogate uses check("[B-lime] the 5-feature local surrogate depends on 2 features", used[0] and used[1] and not used[2:].any()) # write the LIME figure data (surrogate only inside the window) gwin = np.where(np.abs(xs - x0) <= 0.6, g, np.nan) with open("pythondemos/interpretableml_lime.csv", "w") as fh: fh.write("x,h,g\n") for xi, hi, gi in zip(xs, h_hat(xs), gwin): fh.write(f"{xi:.4f},{hi:.4f},{gi:.4f}\n") # ---------------------------------------------------------------- [B-eerm] # EERM (Zhang et al., 2024): regularize the training of a high-capacity # model by penalizing the departure from the predictions one user supplies. # The user is characterized by those predictions; here they come from a # proxy model the user considers interpretable, rather than being elicited # one data point at a time. rng_e2 = np.random.default_rng(9) m2 = 14 xt = np.sort(rng_e2.uniform(-1.0, 1.0, m2)) # training inputs f_true = lambda x: 1.2 * x + 0.4 * np.sin(2.0 * np.pi * x) # noqa: E731 yt = f_true(xt) + 0.4 * rng_e2.standard_normal(m2) deg = 10 # high-capacity polynomial model Phi = np.vander(xt, deg + 1, increasing=True) # train design xg = np.linspace(-1.0, 1.0, 201) # unlabeled grid for the penalty Psi = np.vander(xg, deg + 1, increasing=True) # grid design # the user predictions, taken from a proxy model the user can follow: # the least-squares line through the data cl, *_ = np.linalg.lstsq(np.c_[np.ones(m2), xt], yt, rcond=None) surr = cl[0] + cl[1] * xg # user predictions on the grid w_plain, *_ = np.linalg.lstsq(Phi, yt, rcond=None) # unregularized training lam2 = 1.0 # EERM penalty weight # EERM objective (1/m)||y - Phi w||^2 + lam (1/N)||Psi w - surr||^2 as one # stacked least-squares problem Astack = np.vstack([Phi / np.sqrt(m2), np.sqrt(lam2 / len(xg)) * Psi]) bstack = np.concatenate([yt / np.sqrt(m2), np.sqrt(lam2 / len(xg)) * surr]) w_eerm, *_ = np.linalg.lstsq(Astack, bstack, rcond=None) h_plain, h_eerm = Psi @ w_plain, Psi @ w_eerm disc_plain = np.mean((h_plain - surr) ** 2) # departure from the user disc_eerm = np.mean((h_eerm - surr) ** 2) test_plain = np.mean((h_plain - f_true(xg)) ** 2) # error vs the true mean test_eerm = np.mean((h_eerm - f_true(xg)) ** 2) tr_plain = np.mean((yt - Phi @ w_plain) ** 2) tr_eerm = np.mean((yt - Phi @ w_eerm) ** 2) check("[B-eerm] the penalty shrinks the departure from the user " "predictions and the test error", disc_eerm < disc_plain and test_eerm < test_plain and tr_plain <= tr_eerm) with open("pythondemos/interpretableml_eerm.csv", "w") as fh: fh.write("x,hopaque,hsurr,heerm\n") for r in zip(xg, h_plain, surr, h_eerm): fh.write(",".join(f"{v:.4f}" for v in r) + "\n") with open("pythondemos/interpretableml_eermdata.csv", "w") as fh: fh.write("x,y\n") for xi, yi in zip(xt, yt): fh.write(f"{xi:.4f},{yi:.4f}\n") n_ok = sum(ok for _, ok in report) print(f"\n{n_ok}/{len(report)} checks pass " f"(train MSE dense {train_dense:.3f} < sparse {train_sparse:.3f}; " f"test MSE sparse {test_sparse:.3f} < dense {test_dense:.3f}; " f"EERM: disc {disc_plain:.3f} -> {disc_eerm:.3f}, " f"test {test_plain:.3f} -> {test_eerm:.3f})") print("wrote pythondemos/interpretableml_lime.csv, " "pythondemos/interpretableml_eerm.csv, " "pythondemos/interpretableml_eermdata.csv") if n_ok != len(report): raise SystemExit(1)