"""Bishop (2006), Chapter 7: Sparse Kernel Machines. Reproduces the numerical figures of Chapter 7 of C. M. Bishop, Pattern Recognition and Machine Learning (Springer, 2006), pp. 325-358, written from scratch with numpy and matplotlib (no optimisation library is used at all). Every algorithm is a small function meant to be read next to the lecture notes (ch07-sparse-kernel-machines.html): * Section 7.1, support vector machines. smo_solve is one sequential-minimal-optimisation solver (Platt 1999, with the second-order working-set rule of Fan, Chen and Lin 2005) for the dual quadratic programme min 1/2 a'Qa + p'a s.t. y'a = 0, 0 <= a <= C. It is used four ways: - svc_fit: hard-margin (C = infinity, dual (7.32)-(7.34)) and soft-margin (C finite, dual (7.43)-(7.45)) classification, decision function (7.13) and (7.37), b from (7.37) averaged over the unbounded support vectors; - nu_svc_fit: the nu-SVM (7.49)-(7.52), obtained as the C-SVM whose support-vector fraction is nu (Scholkopf et al. 2000: the two formulations have the same solutions); - svr_fit: epsilon-insensitive regression, dual (7.61)-(7.63) with coefficients (a_n - a^_n) as in (7.64), written as one 2N-variable problem with y = (+1,...,+1,-1,...,-1); - nu_svr_fit: the nu-SVR (7.68)-(7.72), again as the epsilon-SVR whose support-vector fraction is nu. The primal and dual objectives, the margin 1/||w|| and the KKT residuals are returned so that checks_ch07.py can verify them. * Section 7.2, relevance vector machines. rvm_regression is the evidence-maximising RVM: the posterior (7.81)-(7.83), the re-estimation equations (7.87)-(7.89) (MacKay updates) or the EM updates, pruning of basis functions whose alpha_i diverges; log_evidence_direct (7.85)-(7.86) and log_evidence_woodbury (the same number from the M x M matrices, (3.86)) are two routes to the marginal likelihood. rvm_fast_regression is the sequential sparse Bayesian learning algorithm of Section 7.2.2 (sparsity s_i and quality q_i, (7.97)-(7.101), Faul and Tipping 2002), whose log marginal likelihood never decreases. rvm_classification is Section 7.2.3: IRLS for the mode (7.110)-(7.112), the Laplace covariance (7.113), the evidence (7.114) and the update (7.116). Figures: 7.1, 7.2, 7.4, 7.5, 7.6, 7.8, 7.9, 7.10, 7.11, 7.12 (7.3 and 7.7 are pure sketches of slack variables and are not reproduced). Run `python3 ch07_sparse_kernel_machines.py --list`, `python3 ch07_sparse_kernel_machines.py` (all figures, a few seconds) or `--fig 7.2 7.12`. Files go to ../img/ch07/. What differs from the book, and the settings chosen where the book gives none: * The book's synthetic data (Appendix A) are not distributed with it. Three stand-ins are used, all generated with seeded numpy generators: - Figures 7.1 and 7.2: a small separable two-class set (Fig. 7.2: a blue Gaussian blob inside a red ring, 14 + 18 points, rejection-sampled so that the classes are separable; Gaussian kernel exp(-gamma ||x - x'||^2), gamma = 0.8, hard margin as the figure shows no slack). - Figures 7.4 and 7.12: the Ripley-type data of Appendix A (each class an equal mixture of two Gaussians, centres (-0.7, 0.3), (0.3, 0.3) and (-0.3, 0.7), (0.4, 0.7), standard deviation 0.25, which overlap heavily), 100 points per class, shifted and scaled to zero mean and unit standard deviation per axis (the book's plot has axes about -2.3 to 2.6). Fig. 7.4: nu-SVM with nu = 0.5 and Gaussian kernel gamma = 0.45 (the book's gamma; its nu is not stated). Fig. 7.12: RVM with the Gaussian kernel exp(-gamma||x-x'||^2), gamma = 0.45 as well. - Figures 7.8 and 7.9: ten points t = sin(2 pi x) + N(0, 0.2^2) at x = 0, 1/9, ..., 1 (seed 3), Gaussian kernel exp(-gamma (x-x')^2), gamma = 15 for the nu-SVR (nu = 0.7, C = 5) and the RVM (the book states only 'the same Gaussian kernel functions'). * The nu-SVM is solved as a C-SVM with C found by bisection so that the support-vector fraction equals nu; at the solution the margin-error fraction is at most nu, as (7.50)-(7.52) require. * Figure 7.10 is a schematic in the book; here it is computed with a chosen t and phi. * Figure 7.11 plots (7.97) for (q^2, s) = (4, 1) and (1, 2). """ import numpy as np import matplotlib.pyplot as plt import common # --------------------------------------------------------------------------- kernels class Kernel: """k(x, x') for the two kernels used here: 'linear' x'x' and 'rbf' exp(-gamma ||x - x'||^2) (6.23).""" def __init__(self, kind="rbf", gamma=1.0): self.kind, self.gamma = kind, gamma def __call__(self, A, B): A, B = np.atleast_2d(A), np.atleast_2d(B) if self.kind == "linear": return A @ B.T d2 = (A ** 2).sum(1)[:, None] + (B ** 2).sum(1)[None, :] - 2 * A @ B.T return np.exp(-self.gamma * np.maximum(d2, 0.0)) # --------------------------------------------------------------------------- SMO def smo_solve(Q, p, y, C, tol=1e-10, max_iter=200000, a0=None): """Solve min_a 1/2 a'Qa + p'a s.t. y'a = 0, 0 <= a_i <= C (y_i = +-1) by SMO. Each iteration picks the maximal violating pair i in I_up (largest -y_i G_i) and the partner j in I_low that gives the largest decrease of the objective (second-order rule), solves the two-variable sub-problem exactly and clips it to the box. G = Qa + p is updated in O(N). Stops when the KKT gap max_{I_up}(-yG) - min_{I_low}(-yG) falls below tol. Returns a, G, iterations. C may be np.inf (hard margin). """ n = len(p) a = np.zeros(n) if a0 is None else np.array(a0, float) G = Q @ a + p TAU = 1e-12 for it in range(max_iter): g = -y * G up = ((y > 0) & (a < C)) | ((y < 0) & (a > 0)) low = ((y > 0) & (a > 0)) | ((y < 0) & (a < C)) gu = np.where(up, g, -np.inf) i = int(np.argmax(gu)) Gmax = gu[i] gl = np.where(low, g, np.inf) if Gmax - gl.min() < tol: return a, G, it bt = Gmax - g cand = low & (g < Gmax) aij = Q[i, i] + np.diag(Q) - 2 * y[i] * y * Q[i] aij = np.where(aij <= 0, TAU, aij) obj = np.where(cand, -(bt ** 2) / aij, np.inf) j = int(np.argmin(obj)) if y[i] != y[j]: quad = Q[i, i] + Q[j, j] + 2 * Q[i, j] quad = quad if quad > 0 else TAU delta = (-G[i] - G[j]) / quad diff = a[i] - a[j] ai, aj = a[i] + delta, a[j] + delta if diff > 0: if aj < 0: aj, ai = 0.0, diff else: if ai < 0: ai, aj = 0.0, -diff if diff > 0: if ai > C: ai, aj = C, C - diff else: if aj > C: aj, ai = C, C + diff else: quad = Q[i, i] + Q[j, j] - 2 * Q[i, j] quad = quad if quad > 0 else TAU delta = (G[i] - G[j]) / quad s = a[i] + a[j] ai, aj = a[i] - delta, a[j] + delta if s > C: if ai > C: ai, aj = C, s - C else: if aj < 0: aj, ai = 0.0, s if s > C: if aj > C: aj, ai = C, s - C else: if ai < 0: ai, aj = 0.0, s dai, daj = ai - a[i], aj - a[j] a[i], a[j] = ai, aj G += Q[:, i] * dai + Q[:, j] * daj return a, G, max_iter def _bias(a, G, y, C, eps=1e-9): """b from the unbounded multipliers, b = -y_n G_n (7.37 / 7.64); midpoint of the bounds if none.""" free = (a > eps) & (a < C - eps) if free.any(): return float(np.mean(-y[free] * G[free])) ub, lb = np.inf, -np.inf yG = y * G for k in range(len(a)): if a[k] >= C - eps: if y[k] < 0: ub = min(ub, yG[k]) else: lb = max(lb, yG[k]) else: if y[k] > 0: ub = min(ub, yG[k]) else: lb = max(lb, yG[k]) return float(-(ub + lb) / 2) # --------------------------------------------------------------------------- SVM classification class SVC: """A fitted support vector classifier: f(x) = sum_n coef_n k(x, x_n) + b, coef_n = a_n t_n (7.13), (7.37).""" def __init__(self, X, t, kernel, C, a, G, iters): self.X, self.t, self.kernel, self.C, self.a, self.iters = X, t, kernel, C, a, iters self.coef = a * t self.b = _bias(a, G, t.astype(float), C) K = kernel(X, X) self.K = K self.wnorm2 = float(self.coef @ K @ self.coef) # ||w||^2 = sum_n sum_m a_n a_m t_n t_m k(x_n, x_m) self.f_train = K @ self.coef + self.b self.xi = np.maximum(0.0, 1 - t * self.f_train) # slack variables (7.38) self.dual = float(a.sum() - 0.5 * self.wnorm2) # (7.32) / (7.43) self.primal = 0.5 * self.wnorm2 + (C * self.xi.sum() if np.isfinite(C) else 0.0) # (7.40) self.sv = a > 1e-8 * (1.0 if not np.isfinite(C) else min(1.0, C)) self.margin = 1.0 / np.sqrt(self.wnorm2) def decision_function(self, Z): return self.kernel(np.atleast_2d(Z), self.X) @ self.coef + self.b def predict(self, Z): return np.where(self.decision_function(Z) >= 0, 1, -1) def svc_fit(X, t, kernel=None, C=np.inf, tol=1e-10): """Maximum-margin classifier, targets t in {-1,+1}. C = inf is the hard margin of Section 7.1 (dual (7.32)-(7.34)); finite C is the soft margin of Section 7.1.1 (dual (7.43)-(7.45)).""" X = np.asarray(X, float) t = np.asarray(t, float) kernel = kernel or Kernel("linear") K = kernel(X, X) Q = (t[:, None] * t[None, :]) * K a, G, it = smo_solve(Q, -np.ones(len(t)), t, C, tol=tol) return SVC(X, t, kernel, C, a, G, it) def nu_svc_fit(X, t, nu, kernel=None, tol=1e-7): """nu-SVM (Section 7.1.1, (7.49)-(7.52)): the C-SVM whose fraction of support vectors is nu. nu is an upper bound on the fraction of margin errors and a lower bound on that of support vectors. C is found by bisection on ln C (the support-vector fraction decreases with C).""" lo, hi = -6.0, 5.0 # for nu below the smallest attainable fraction (about 0.5 on heavily overlapping data) the search ends at C = e^5 best = None for _ in range(40): mid = 0.5 * (lo + hi) m = svc_fit(X, t, kernel, np.exp(mid), tol=tol) frac = m.sv.mean() best = m if frac > nu: lo = mid else: hi = mid if abs(frac - nu) < 0.5 / len(t): break return best # --------------------------------------------------------------------------- SVM regression class SVR: """f(x) = sum_n (a_n - a^_n) k(x, x_n) + b (7.64), with the epsilon-insensitive tube (7.51-7.54).""" def __init__(self, X, t, kernel, C, eps, a, G, iters): n = len(t) self.X, self.t, self.kernel, self.C, self.eps, self.a, self.iters = X, t, kernel, C, eps, a, iters self.a_up, self.a_lo = a[:n], a[n:] self.coef = self.a_up - self.a_lo yy = np.r_[np.ones(n), -np.ones(n)] self.b = _bias(a, G, yy, C) K = kernel(X, X) self.f_train = K @ self.coef + self.b r = t - self.f_train self.xi = np.maximum(0.0, r - eps) # xi_n (7.53): t_n above the tube, t_n > y_n + eps; its multiplier is a_n self.xi_hat = np.maximum(0.0, -r - eps) # xi^_n (7.54): t_n below the tube, t_n < y_n - eps; its multiplier is a^_n self.wnorm2 = float(self.coef @ K @ self.coef) self.primal = self.C * float((self.xi + self.xi_hat).sum()) + 0.5 * self.wnorm2 # (7.55) # dual (7.61): -1/2 sum sum (a-a^)(a-a^) k - eps sum (a + a^) + sum (a - a^) t self.dual = float(-0.5 * self.wnorm2 - eps * (self.a_up + self.a_lo).sum() + (self.coef * t).sum()) self.sv = (self.a_up > 1e-8) | (self.a_lo > 1e-8) def predict(self, Z): return self.kernel(np.atleast_2d(Z), self.X) @ self.coef + self.b def svr_fit(X, t, kernel, C, eps, tol=1e-10): """epsilon-SVR: dual (7.61)-(7.63) as the 2N-variable problem min 1/2 z'Qz + p'z, y'z = 0, 0 <= z <= C with z = (a, a^), y = (+1, -1), Q = (y y') (x) K, p = (eps - t, eps + t).""" X = np.asarray(X, float) t = np.asarray(t, float) n = len(t) K = kernel(X, X) y = np.r_[np.ones(n), -np.ones(n)] Q = (y[:, None] * y[None, :]) * np.block([[K, K], [K, K]]) p = np.r_[eps - t, eps + t] a, G, it = smo_solve(Q, p, y, C, tol=tol) return SVR(X, t, kernel, C, eps, a, G, it) def nu_svr_fit(X, t, kernel, C, nu, tol=1e-9): """nu-SVR (7.68)-(7.72): the epsilon-SVR whose support-vector fraction is nu (at most nu N points outside the tube, at least nu N support vectors). epsilon is found by bisection.""" lo, hi = 0.0, float(np.abs(t - t.mean()).max()) + 1e-9 best = None for _ in range(60): mid = 0.5 * (lo + hi) m = svr_fit(X, t, kernel, C, mid, tol=tol) frac = m.sv.mean() best = m if frac > nu: lo = mid else: hi = mid if abs(frac - nu) < 0.5 / len(t) and abs(hi - lo) < 1e-3: break return best # --------------------------------------------------------------------------- RVM regression def design_matrix(X, kernel, centres=None, bias=True): """Phi_ni = k(x_n, x_i) (7.78), with a constant first column for the bias w_0.""" Phi = kernel(X, X if centres is None else centres) return np.hstack([np.ones((len(Phi), 1)), Phi]) if bias else Phi def rvm_posterior(Phi, t, alpha, beta): """Posterior mean m = beta Sigma Phi't (7.82) and covariance Sigma = (A + beta Phi'Phi)^-1 (7.83).""" S = np.linalg.inv(np.diag(alpha) + beta * Phi.T @ Phi) return beta * S @ Phi.T @ t, S def log_evidence_direct(Phi, t, alpha, beta): """ln p(t | alpha, beta) = -1/2 [N ln 2 pi + ln|C| + t'C^-1 t], C = beta^-1 I + Phi A^-1 Phi' (7.85)-(7.86).""" N = len(t) C = np.eye(N) / beta + (Phi / alpha) @ Phi.T sign, logdet = np.linalg.slogdet(C) return -0.5 * (N * np.log(2 * np.pi) + logdet + t @ np.linalg.solve(C, t)) def log_evidence_woodbury(Phi, t, alpha, beta): """The same number from the M x M posterior: ln p = N/2 ln beta - N/2 ln 2pi - E(m) - 1/2 ln|A^-1... (3.86) with a diagonal A: -E(m) + 1/2 sum ln alpha_i + N/2 ln beta - 1/2 ln|Sigma^-1| - N/2 ln 2 pi, E(m) = beta/2 ||t - Phi m||^2 + 1/2 m'Am.""" N = len(t) m, S = rvm_posterior(Phi, t, alpha, beta) E = 0.5 * beta * np.sum((t - Phi @ m) ** 2) + 0.5 * np.sum(alpha * m ** 2) return (-E + 0.5 * np.log(alpha).sum() + 0.5 * N * np.log(beta) - 0.5 * np.linalg.slogdet(np.linalg.inv(S))[1] - 0.5 * N * np.log(2 * np.pi)) def rvm_regression(Phi, t, n_iter=2000, alpha0=1.0, beta0=None, update="mackay", prune=1e9, tol=1e-9): """RVM regression by evidence maximisation (Section 7.2.1). update = 'mackay': alpha_i = gamma_i / m_i^2 (7.87), beta^-1 = ||t - Phi m||^2 / (N - sum gamma_i) (7.88), gamma_i = 1 - alpha_i Sigma_ii (7.89). update = 'em': alpha_i = 1/(m_i^2 + Sigma_ii), and beta^-1 = (||t - Phi m||^2 + beta^-1 sum gamma_i) / N, the EM updates of Section 9.3.4, which cannot decrease the evidence. Basis functions with alpha_i > prune are removed. Returns a dict with the active index set, m, Sigma, alpha, beta and the trace of ln p(t|alpha, beta) after every iteration. """ N, M = Phi.shape act = np.arange(M) alpha = np.full(M, float(alpha0)) beta = 1.0 / (0.1 * np.var(t)) if beta0 is None else beta0 trace = [] for it in range(n_iter): P = Phi[:, act] m, S = rvm_posterior(P, t, alpha[act], beta) trace.append(log_evidence_direct(P, t, alpha[act], beta)) gam = 1 - alpha[act] * np.diag(S) resid = np.sum((t - P @ m) ** 2) if update == "mackay": a_new = gam / np.maximum(m ** 2, 1e-300) b_inv = resid / max(N - gam.sum(), 1e-12) else: a_new = 1.0 / (m ** 2 + np.diag(S)) b_inv = (resid + gam.sum() / beta) / N delta = np.max(np.abs(np.log(a_new[a_new < prune]) - np.log(alpha[act][a_new < prune]))) if (a_new < prune).any() else 0 alpha_act = a_new keep = alpha_act < prune alpha[act] = alpha_act beta = 1.0 / b_inv act = act[keep] if it > 2 and delta < tol: break P = Phi[:, act] m, S = rvm_posterior(P, t, alpha[act], beta) trace.append(log_evidence_direct(P, t, alpha[act], beta)) return dict(active=act, m=m, Sigma=S, alpha=alpha[act], beta=beta, trace=np.array(trace), iters=it + 1) def rvm_predict(Phi_star, fit): """Predictive mean m'phi(x) (7.90) and variance beta^-1 + phi' Sigma phi (7.91).""" P = Phi_star[:, fit["active"]] return P @ fit["m"], 1.0 / fit["beta"] + np.einsum("ij,jk,ik->i", P, fit["Sigma"], P) def rvm_fast_regression(Phi, t, beta, n_iter=500, tol=1e-10): """Sequential sparse Bayesian learning (Section 7.2.2) with beta fixed. Start with the basis vector of largest |phi't|/|phi| and alpha from (7.101); then cycle through all candidates: add (alpha_i = s_i^2/(q_i^2 - s_i), (7.101)), re-estimate or delete (alpha_i = inf) according to whether q_i^2 > s_i, using Q_i, S_i from (7.102)-(7.107). Each action maximises the marginal likelihood over one alpha_i exactly, so ln p(t|alpha) never decreases. Returns the active set, alpha, trace of ln p after each step. """ N, M = Phi.shape alpha = np.full(M, np.inf) i0 = int(np.argmax((Phi.T @ t) ** 2 / np.sum(Phi ** 2, 0))) nrm = np.sum(Phi[:, i0] ** 2) alpha[i0] = nrm / ((Phi[:, i0] @ t) ** 2 / nrm - 1 / beta + 1e-300) if (Phi[:, i0] @ t) ** 2 / nrm > 1 / beta else 1.0 trace = [] for sweep in range(n_iter): changed = False for i in range(M): act = np.isfinite(alpha) P = Phi[:, act] Cinv_P = None C = np.eye(N) / beta + (P / alpha[act]) @ P.T Cinv = np.linalg.inv(C) phi = Phi[:, i] Si = phi @ Cinv @ phi Qi = phi @ Cinv @ t if np.isfinite(alpha[i]): s = alpha[i] * Si / (alpha[i] - Si) q = alpha[i] * Qi / (alpha[i] - Si) else: s, q = Si, Qi if q ** 2 > s: new = s ** 2 / (q ** 2 - s) # (7.101) else: new = np.inf if new != alpha[i] and not (np.isinf(new) and np.isinf(alpha[i])): if np.isfinite(new) and np.isfinite(alpha[i]) and abs(np.log(new / alpha[i])) < tol: continue if np.isinf(new) and act.sum() == 1: continue alpha[i] = new changed = True act2 = np.isfinite(alpha) trace.append(log_evidence_direct(Phi[:, act2], t, alpha[act2], beta)) if not changed: break act = np.where(np.isfinite(alpha))[0] m, S = rvm_posterior(Phi[:, act], t, alpha[act], beta) return dict(active=act, m=m, Sigma=S, alpha=alpha[act], beta=beta, trace=np.array(trace)) # --------------------------------------------------------------------------- RVM classification def _bern_logpost(w, Phi, t, alpha): a = Phi @ w # ln p(t|w) = sum t a - ln(1 + e^a) (stable) return float(np.sum(t * a - np.logaddexp(0, a)) - 0.5 * np.sum(alpha * w ** 2)) # (7.109) up to a constant def logpost_grad(w, Phi, t, alpha): """grad ln p(w|t, alpha) = Phi'(t - y) - A w (7.110).""" y = common.sigmoid(Phi @ w) return Phi.T @ (t - y) - alpha * w def logpost_hess(w, Phi, t, alpha): """Hessian = -(Phi' B Phi + A), B = diag(y_n (1 - y_n)) (7.111).""" y = common.sigmoid(Phi @ w) return -(Phi.T @ (Phi * (y * (1 - y))[:, None]) + np.diag(alpha)) def laplace_mode(Phi, t, alpha, w0=None, n_iter=100, tol=1e-12): """IRLS (Newton-Raphson with step halving) for the mode w* of (7.109); returns w*, Sigma (7.113).""" w = np.zeros(Phi.shape[1]) if w0 is None else w0.copy() f = _bern_logpost(w, Phi, t, alpha) for _ in range(n_iter): g = logpost_grad(w, Phi, t, alpha) H = logpost_hess(w, Phi, t, alpha) step = np.linalg.solve(H, g) s = 1.0 while s > 1e-10: wn = w - s * step fn = _bern_logpost(wn, Phi, t, alpha) if fn >= f - 1e-14: break s *= 0.5 w, fold, f = wn, f, fn if np.max(np.abs(g)) < tol or abs(f - fold) < 1e-15: break H = logpost_hess(w, Phi, t, alpha) return w, np.linalg.inv(-H) def laplace_log_evidence(Phi, t, alpha, w, Sigma): """ln p(t|alpha) = ln p(t|w*) + ln p(w*|alpha) + M/2 ln 2 pi + 1/2 ln|Sigma| (7.114).""" M = len(w) a = Phi @ w lik = np.sum(t * a - np.logaddexp(0, a)) lprior = 0.5 * np.log(alpha).sum() - 0.5 * M * np.log(2 * np.pi) - 0.5 * np.sum(alpha * w ** 2) return float(lik + lprior + 0.5 * M * np.log(2 * np.pi) + 0.5 * np.linalg.slogdet(Sigma)[1]) def rvm_classification(Phi, t, n_iter=500, alpha0=1.0, prune=1e9, tol=1e-6): """RVM classification (Section 7.2.3), targets t in {0, 1}: IRLS for w* given alpha, then alpha_i = gamma_i / w*_i^2 (7.116), gamma_i = 1 - alpha_i Sigma_ii, pruning alpha_i > prune; repeat.""" N, M = Phi.shape act = np.arange(M) alpha = np.full(M, float(alpha0)) w = np.zeros(M) trace = [] for it in range(n_iter): P = Phi[:, act] w_act, S = laplace_mode(P, t, alpha[act], w0=w[act]) w[:] = 0 w[act] = w_act trace.append(laplace_log_evidence(P, t, alpha[act], w_act, S)) gam = 1 - alpha[act] * np.diag(S) a_new = gam / np.maximum(w_act ** 2, 1e-300) delta = np.max(np.abs(np.log(a_new[a_new < prune]) - np.log(alpha[act][a_new < prune]))) if (a_new < prune).any() else 0 alpha[act] = a_new act = act[a_new < prune] if it > 3 and delta < tol: break P = Phi[:, act] w_act, S = laplace_mode(P, t, alpha[act], w0=w[act]) trace.append(laplace_log_evidence(P, t, alpha[act], w_act, S)) return dict(active=act, w=w_act, Sigma=S, alpha=alpha[act], trace=np.array(trace), iters=it + 1) def rvm_class_prob(Phi_star, fit, probit_approx=True): """p(C1|x) = sigma(kappa(sigma_a^2) mu_a) with the probit approximation (4.153) to the predictive (4.151).""" P = Phi_star[:, fit["active"]] mu = P @ fit["w"] if not probit_approx: return common.sigmoid(mu) var = np.einsum("ij,jk,ik->i", P, fit["Sigma"], P) return common.sigmoid(mu / np.sqrt(1 + np.pi * var / 8)) # --------------------------------------------------------------------------- data def ripley_like(n_per_class=100, seed=11): """Two heavily overlapping classes, each a 50:50 mixture of two spherical Gaussians (sd 0.25), standardised. Stand-in for the Ripley synthetic data of Appendix A. Returns X (2n x 2), t in {-1, +1} (blue = -1, red = +1).""" rng = np.random.default_rng(seed) c0 = np.array([[-0.7, 0.3], [0.3, 0.3]]) c1 = np.array([[-0.3, 0.7], [0.4, 0.7]]) X, t = [], [] for lab, cs in ((-1, c0), (1, c1)): k = rng.integers(0, 2, n_per_class) X.append(cs[k] + 0.25 * rng.normal(size=(n_per_class, 2))) t += [lab] * n_per_class X = np.vstack(X) X = X[:, ::-1] # axes swapped so that the blue class lies on the left, as in the book's plots X = (X - X.mean(0)) / X.std(0) return X, np.array(t) def sinusoid_data(n=10, noise=0.2, seed=3): """t = sin(2 pi x) + N(0, noise^2) at x = 0, 1/(n-1), ..., 1 (the synthetic regression set of Appendix A).""" rng = np.random.default_rng(seed) x = np.linspace(0, 1, n) return x[:, None], np.sin(2 * np.pi * x) + noise * rng.normal(size=n) def blob_ring_data(seed=5, n_blue=14, n_red=18): """A blue Gaussian blob surrounded by a red ring, rejection-sampled to be separable by a 0.55 gap (Fig. 7.2 stand-in).""" rng = np.random.default_rng(seed) blue = rng.normal(0, 0.45, size=(n_blue, 2)) red = [] while len(red) < n_red: r = rng.uniform(1.0, 2.0) th = rng.uniform(0, 2 * np.pi) p = np.array([r * np.cos(th) * 1.25, r * np.sin(th)]) if np.min(np.linalg.norm(blue - p, axis=1)) > 0.8 and all(np.linalg.norm(p - q) > 0.25 for q in red): red.append(p) X = np.vstack([blue, np.array(red)]) return X, np.r_[-np.ones(n_blue), np.ones(n_red)] # --------------------------------------------------------------------------- plotting helpers BLUE_X, RED_X = "#0000ff", "#ff0000" CIRC = "#00e000" def _scatter_classes(ax, X, t, ms=28): ax.scatter(*X[t > 0].T, marker="x", c=RED_X, s=ms, linewidths=1.6, zorder=3) ax.scatter(*X[t < 0].T, marker="x", c=BLUE_X, s=ms, linewidths=1.6, zorder=3) def _circle_sv(ax, X, mask, s=130): ax.scatter(*X[mask].T, s=s, facecolors="none", edgecolors=CIRC, linewidths=2.0, zorder=4) def _grid(X, pad=0.5, n=160): x0, x1 = X[:, 0].min() - pad, X[:, 0].max() + pad y0, y1 = X[:, 1].min() - pad, X[:, 1].max() + pad gx, gy = np.meshgrid(np.linspace(x0, x1, n), np.linspace(y0, y1, n)) return gx, gy, np.c_[gx.ravel(), gy.ravel()] # --------------------------------------------------------------------------- figures def fig7_1(out_dir=None): """Figure 7.1: the margin (left: an arbitrary separating line) and the maximum-margin boundary with its support vectors (right).""" rng = np.random.default_rng(10) Xb = rng.normal([-1.0, -0.8], 0.55, size=(12, 2)) Xr = rng.normal([1.0, 0.9], 0.55, size=(12, 2)) X = np.vstack([Xb, Xr]) t = np.r_[-np.ones(12), np.ones(12)] m = svc_fit(X, t) w = m.coef @ X fig, axes = plt.subplots(1, 2, figsize=(7.4, 3.5), sharex=True, sharey=True) xs = np.array([-3.3, 3.3]) # (a) an arbitrary separating line through the middle gap, tilted, with its own margin w_a = np.array([0.55, 1.0]) b_a = -(w_a @ (X[t < 0].mean(0) + X[t > 0].mean(0)) / 2) + 0.15 for ax in axes: _scatter_classes(ax, X, t) ax.set_xlim(-3, 3) ax.set_ylim(-2.7, 2.9) ax.set_xticks([]) ax.set_yticks([]) f = (X @ w_a + b_a) / np.linalg.norm(w_a) dist = np.abs(f) ax = axes[0] ax.plot(xs, -(w_a[0] * xs + b_a) / w_a[1], "k-", lw=1.4) k = np.argmin(dist) d = dist[k] n_a = w_a / np.linalg.norm(w_a) for sgn in (+1, -1): ax.plot(xs, -(w_a[0] * xs + b_a - sgn * d * np.linalg.norm(w_a)) / w_a[1], "k--", lw=0.9) ax.set_title("(a) a separating line and its margin %.2f" % d, fontsize=9.5) ax = axes[1] ax.plot(xs, -(w[0] * xs + m.b) / w[1], "k-", lw=1.4) for sgn in (+1, -1): ax.plot(xs, -(w[0] * xs + m.b - sgn) / w[1], "k--", lw=0.9) _circle_sv(ax, X, m.sv) ax.set_title("(b) maximum margin %.2f, %d support vectors" % (m.margin, m.sv.sum()), fontsize=9.5) fig.tight_layout() common.savefig(fig, 7, "fig7-1", out_dir) def fig7_2(out_dir=None): """Figure 7.2: hard-margin SVM with a Gaussian kernel on synthetic two-class data: contours of y(x), boundary, margins, support vectors.""" X, t = blob_ring_data() m = svc_fit(X, t, Kernel("rbf", 0.8)) gx, gy, G = _grid(X, 0.7) Z = m.decision_function(G).reshape(gx.shape) fig, ax = plt.subplots(figsize=(4.9, 3.7)) levels = [-3, -2, -1, 0, 1, 2, 3] ax.contour(gx, gy, Z, levels=[l for l in levels if l not in (-1, 0, 1)], colors="k", linewidths=0.6, linestyles="solid") ax.contour(gx, gy, Z, levels=[-1, 1], colors="k", linewidths=1.5, linestyles="solid") ax.contour(gx, gy, Z, levels=[0], colors="k", linewidths=3.0) _scatter_classes(ax, X, t) _circle_sv(ax, X, m.sv) ax.set_xticks([]) ax.set_yticks([]) fig.tight_layout() print(" Fig 7.2: N = %d, support vectors = %d, hard margin, gamma = 0.8, max alpha = %.3g" % (len(t), m.sv.sum(), m.a.max())) common.savefig(fig, 7, "fig7-2", out_dir) def fig7_4(out_dir=None): """Figure 7.4: nu-SVM (nu = 0.5, Gaussian kernel gamma = 0.45) on overlapping classes; support vectors circled.""" X, t = ripley_like() m = nu_svc_fit(X, t, 0.5, Kernel("rbf", 0.45)) gx, gy, G = _grid(X, 0.4) Z = m.decision_function(G).reshape(gx.shape) fig, ax = plt.subplots(figsize=(5.2, 3.8)) ax.contour(gx, gy, Z, levels=[0], colors="k", linewidths=1.5) _scatter_classes(ax, X, t, 22) _circle_sv(ax, X, m.sv, 90) ax.set_xlim(gx.min(), gx.max()) ax.set_ylim(gy.min(), gy.max()) fig.tight_layout() me = (m.xi > 1e-8).mean() print(" Fig 7.4: C = %.3g, SV fraction %.3f (nu = 0.5), margin-error fraction %.3f, training error %.3f" % (m.C, m.sv.mean(), me, np.mean(np.where(m.f_train >= 0, 1, -1) != t))) common.savefig(fig, 7, "fig7-4", out_dir) def fig7_5(out_dir=None): """Figure 7.5: the hinge error (7.45-style, blue), logistic error ln(1 + e^-z)/ln 2 (red), misclassification (black) and squared error (green) against z = t y.""" z = np.linspace(-2, 2, 400) fig, ax = plt.subplots(figsize=(4.6, 3.4)) ax.plot(z, np.maximum(0, 1 - z), color="#0000ff", lw=2, label="hinge $E_{SV}$") ax.plot(z, np.logaddexp(0, -z) / np.log(2), color="#ff0000", lw=2, label="logistic (rescaled)") ax.plot(z[z < 0], np.ones((z < 0).sum()), "k-", lw=1.2, label="misclassification") ax.plot(z[z >= 0], np.zeros((z >= 0).sum()), "k-", lw=1.2) ax.plot(z, (z - 1) ** 2, color="#00c000", lw=2, label="squared") ax.set_ylim(-0.1, 4.0) ax.set_xlim(-2, 2) ax.set_xlabel("$z = t\\,y(x)$") ax.set_ylabel("$E(z)$") ax.axhline(0, color="k", lw=0.6) ax.axvline(0, color="k", lw=0.6) ax.legend(frameon=False, fontsize=8) fig.tight_layout() common.savefig(fig, 7, "fig7-5", out_dir) def fig7_6(out_dir=None): """Figure 7.6: the epsilon-insensitive error (7.51, red) and the quadratic error (green).""" z = np.linspace(-1.5, 1.5, 400) eps = 0.5 fig, ax = plt.subplots(figsize=(4.6, 3.0)) ax.plot(z, np.maximum(0, np.abs(z) - eps), color="#ff0000", lw=2, label="$\\epsilon$-insensitive, $\\epsilon=0.5$") ax.plot(z, 0.5 * z ** 2, color="#00c000", lw=2, label="quadratic") ax.set_xticks([-eps, 0, eps]) ax.set_xticklabels(["$-\\epsilon$", "0", "$\\epsilon$"]) ax.set_ylim(-0.05, 1.0) ax.set_xlabel("$z = y(x) - t$") ax.set_ylabel("$E(z)$") ax.legend(frameon=False, fontsize=8, loc="upper center") fig.tight_layout() common.savefig(fig, 7, "fig7-6", out_dir) def _sin_axes(ax): ax.set_xlim(-0.1, 1.1) ax.set_ylim(-1.6, 1.6) ax.set_xticks([0, 1]) ax.set_yticks([-1, 0, 1]) ax.set_xlabel("$x$") ax.set_ylabel("$t$", rotation=0) def fig7_8(out_dir=None): """Figure 7.8: nu-SVR (nu = 0.7, C = 5, gamma = 15) on the sinusoidal data: curve, epsilon-tube, support vectors circled.""" X, t = sinusoid_data() m = nu_svr_fit(X, t, Kernel("rbf", 15.0), 5.0, 0.7) xs = np.linspace(0, 1, 300)[:, None] f = m.predict(xs) fig, ax = plt.subplots(figsize=(4.8, 3.5)) ax.fill_between(xs[:, 0], f - m.eps, f + m.eps, color="#ffd6d6", lw=0) ax.plot(xs[:, 0], f, color="#ff0000", lw=2) ax.scatter(X[:, 0], t, s=45, facecolors="none", edgecolors=CIRC, linewidths=2, zorder=3) ax.scatter(X[m.sv, 0], t[m.sv], s=190, facecolors="none", edgecolors="#0000ff", linewidths=2, zorder=3) _sin_axes(ax) fig.tight_layout() print(" Fig 7.8: epsilon = %.3f, support vectors = %d of %d (book: 7), points outside the tube = %d" % (m.eps, m.sv.sum(), len(t), int(((m.xi + m.xi_hat) > 1e-8).sum()))) common.savefig(fig, 7, "fig7-8", out_dir) def fig7_9(out_dir=None): """Figure 7.9: RVM regression on the same data and kernel: predictive mean, one-sd band, relevance vectors circled.""" X, t = sinusoid_data() k = Kernel("rbf", 15.0) Phi = design_matrix(X, k) fit = rvm_regression(Phi, t) xs = np.linspace(0, 1, 300)[:, None] mu, var = rvm_predict(design_matrix(xs, k, centres=X), fit) sd = np.sqrt(var) rv = fit["active"][fit["active"] > 0] - 1 # data points that are relevance vectors (column 0 is the bias) fig, ax = plt.subplots(figsize=(4.8, 3.5)) ax.fill_between(xs[:, 0], mu - sd, mu + sd, color="#ffd6d6", lw=0) ax.plot(xs[:, 0], mu, color="#ff0000", lw=2) ax.scatter(X[:, 0], t, s=45, facecolors="none", edgecolors=CIRC, linewidths=2, zorder=3) ax.scatter(X[rv, 0], t[rv], s=190, facecolors="none", edgecolors="#0000ff", linewidths=2, zorder=3) _sin_axes(ax) fig.tight_layout() print(" Fig 7.9: relevance vectors = %d (book: 3; the bias column is %s), 1/beta^0.5 = %.3f, %d iterations" % (len(rv), "kept" if 0 in fit["active"] else "pruned", fit["beta"] ** -0.5, fit["iters"])) common.savefig(fig, 7, "fig7-9", out_dir) def fig7_10(out_dir=None): """Figure 7.10: the mechanism of sparsity, N = 2, one badly aligned basis vector phi: unit-Mahalanobis ellipses of C = beta^-1 I (alpha = inf, left) and C = beta^-1 I + alpha^-1 phi phi' with the same |C| (right), with ln N(t|0, C).""" t = np.array([1.35, 1.05]) phi = np.array([1.0, -0.55]) sig2 = t @ t / 2 # most probable beta^-1 when alpha = inf: |t|^2 / N # right panel: finite alpha chosen so that |C| is the same (beta^-1 changes), alpha^-1 |phi|^2 = 1.2 ainv = 1.2 / (phi @ phi) # |C| = s (s + ainv |phi|^2) = sig2^2 -> s c = ainv * (phi @ phi) s = (-c + np.sqrt(c * c + 4 * sig2 ** 2)) / 2 Cl = sig2 * np.eye(2) Cr = s * np.eye(2) + ainv * np.outer(phi, phi) fig, axes = plt.subplots(1, 2, figsize=(7.2, 3.5)) th = np.linspace(0, 2 * np.pi, 300) circ = np.vstack([np.cos(th), np.sin(th)]) for ax, C, noise, lab in ((axes[0], Cl, sig2, "(a) $\\alpha=\\infty$"), (axes[1], Cr, s, "(b) finite $\\alpha$, same $|\\mathbf{C}|$")): L = np.linalg.cholesky(C) e = L @ circ ax.plot(e[0], e[1], color="#ff0000", lw=2) ax.plot(np.sqrt(noise) * circ[0], np.sqrt(noise) * circ[1], color="#00c000", lw=1.6, ls="--") ax.plot(*t, "kx", ms=8, mew=2) ax.text(t[0] + 0.07, t[1] + 0.05, "$\\mathbf{t}$") ld = -0.5 * (2 * np.log(2 * np.pi) + np.linalg.slogdet(C)[1] + t @ np.linalg.solve(C, t)) ax.set_title("%s: $\\ln p(\\mathbf{t})=%.3f$" % (lab, ld), fontsize=9.5) ax.axhline(0, color="k", lw=0.6) ax.axvline(0, color="k", lw=0.6) ax.set_xlim(-2.4, 2.4) ax.set_ylim(-2.4, 2.4) ax.set_aspect("equal") ax.set_xlabel("$t_1$") ax.set_ylabel("$t_2$", rotation=0) fig.tight_layout() print(" Fig 7.10: |C| equal: %.5f vs %.5f; ln p(t) %.3f (alpha = inf) > %.3f (finite alpha)" % (np.linalg.det(Cl), np.linalg.det(Cr), -0.5 * (2 * np.log(2 * np.pi) + np.linalg.slogdet(Cl)[1] + t @ np.linalg.solve(Cl, t)), -0.5 * (2 * np.log(2 * np.pi) + np.linalg.slogdet(Cr)[1] + t @ np.linalg.solve(Cr, t)))) common.savefig(fig, 7, "fig7-10", out_dir) def lam(alpha, q2, s): """lambda(alpha_i) = 1/2 [ln alpha_i - ln(alpha_i + s_i) + q_i^2 / (alpha_i + s_i)] (7.97).""" return 0.5 * (np.log(alpha) - np.log(alpha + s) + q2 / (alpha + s)) def fig7_11(out_dir=None): """Figure 7.11: lambda(alpha_i) against ln alpha_i for q^2 = 4, s = 1 (finite maximum at s^2/(q^2-s), left) and q^2 = 1, s = 2 (maximum at infinity, right).""" la = np.linspace(-5, 5, 400) fig, axes = plt.subplots(1, 2, figsize=(7.2, 3.1), sharey=True) for ax, (q2, s), lab in zip(axes, ((4.0, 1.0), (1.0, 2.0)), "ab"): ax.plot(la, lam(np.exp(la), q2, s), color="#ff0000", lw=2) if q2 > s: a_star = s ** 2 / (q2 - s) ax.plot(np.log(a_star), lam(a_star, q2, s), "ko", ms=5) ax.set_title("(%s) $q^2=%g > s=%g$: $\\alpha^\\star=s^2/(q^2-s)=%.3f$" % (lab, q2, s, a_star), fontsize=9) else: ax.set_title("(%s) $q^2=%g < s=%g$: maximum at $\\alpha=\\infty$" % (lab, q2, s), fontsize=9) ax.set_xlabel("$\\ln\\alpha_i$") ax.set_ylim(-4, 2) ax.set_xlim(-5, 5) axes[0].set_ylabel("$\\lambda(\\alpha_i)$") fig.tight_layout() common.savefig(fig, 7, "fig7-11", out_dir) def fig7_12(out_dir=None): """Figure 7.12: RVM classification on the same data as Fig. 7.4: relevance vectors and decision boundary (left), posterior probability (right).""" X, tpm = ripley_like() t = (tpm > 0).astype(float) k = Kernel("rbf", 0.45) Phi = design_matrix(X, k) fit = rvm_classification(Phi, t) gx, gy, G = _grid(X, 0.4) P = common.sigmoid # noqa pr = rvm_class_prob(design_matrix(G, k, centres=X), fit).reshape(gx.shape) rv = fit["active"][fit["active"] > 0] - 1 fig, axes = plt.subplots(1, 2, figsize=(9.4, 3.6)) ax = axes[0] ax.contour(gx, gy, pr, levels=[0.5], colors="k", linewidths=1.5) _scatter_classes(ax, X, tpm, 22) _circle_sv(ax, X[rv], np.ones(len(rv), bool), 90) ax = axes[1] ax.imshow(pr, extent=(gx.min(), gx.max(), gy.min(), gy.max()), origin="lower", cmap="bwr", vmin=0, vmax=1, alpha=0.6, aspect="auto", rasterized=True) ax.contour(gx, gy, pr, levels=[0.5], colors="k", linewidths=1.5) _scatter_classes(ax, X, tpm, 18) for ax in axes: ax.set_xlim(gx.min(), gx.max()) ax.set_ylim(gy.min(), gy.max()) fig.tight_layout() err = np.mean((rvm_class_prob(Phi, fit, False) > 0.5) != (t > 0.5)) print(" Fig 7.12: relevance vectors = %d of %d (nu-SVM of Fig 7.4 used %d), training error %.3f, %d outer iterations" % (len(rv), len(t), nu_svc_fit(X, tpm, 0.5, k).sv.sum(), err, fit["iters"])) common.savefig(fig, 7, "fig7-12", out_dir) FIGURES = {"7.1": fig7_1, "7.2": fig7_2, "7.4": fig7_4, "7.5": fig7_5, "7.6": fig7_6, "7.8": fig7_8, "7.9": fig7_9, "7.10": fig7_10, "7.11": fig7_11, "7.12": fig7_12} if __name__ == "__main__": common.run(FIGURES, 7)