Dictionary of Applied Machine Learning · interpretable machine learning

interpretable machine learning — Python demo

Numerical companion to the entry interpretable machine learning: it recomputes what the entry states and prints one line per check

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.

Run it with python3 pythondemos/interpretableml.py, from the repository root — it writes its data files under pythondemos/. Requires NumPy only, and uses fixed seeds, so the printed numbers reproduce exactly. Download interpretableml.py

The script, block by block

One cell per block of the script: the code, and what that code printed when it last ran here

setup

"""
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

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.

check("[B-reg]   dense fits training data better but generalizes worse",
      train_dense < train_sparse and test_sparse < test_dense)
  [ok] [B-reg]   dense fits training data better but generalizes worse

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.

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)
  [ok] [B-acc]   sparse hypothesis uses few features at comparable test error (3 of 20)

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.

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)
  [ok] [B-tree]  tree prediction equals the two-rule lookup on a 100 x 100 grid; 3 cells
  feature 1 moved by 0.7: prediction moves -0.5348, weight times step -0.5348

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.

# 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)
  [ok] [B-effect] a feature change moves the prediction by weight times step
  label 4 perturbed by 0.5: prediction shifts -0.0544, normal equations predict -0.0544
  [ok] [B-effect] a perturbed training label shifts the prediction by the amount the normal equations predict

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.

# 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")
  [ok] [B-lime]  weighted linear fit matches h near x0, not globally
  [ok] [B-lime]  the 5-feature local surrogate depends on 2 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.

# 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)
  [ok] [B-eerm]  the penalty shrinks the departure from the user predictions and the test error

8/8 checks pass (train MSE dense 0.136 < sparse 0.262; test MSE sparse 1.781 < dense 5.558; EERM: disc 0.807 -> 0.008, test 0.749 -> 0.058)
wrote pythondemos/interpretableml_lime.csv, pythondemos/interpretableml_eerm.csv, pythondemos/interpretableml_eermdata.csv