P-kernel writes when the script runsDictionary of Applied Machine Learning · support vector machine
Numerical companion to the entry support vector machine: it recomputes what the entry states and prints one line per check
A small, hand-crafted non-separable toy trainset in R^2 whose soft-margin SVM solution leaves exactly one data point misclassified *beyond* the opposite margin, i.e. with functional margin y (w^T x + b) < -1 (a bounded support vector, hinge loss > 2). The same dataset drives the entry's non-separable figure (Fig. 2). Self-contained (numpy/matplotlib only), fixed seed; the SVM is solved by a pure-numpy simplified SMO so no external solver is needed.
Run it with python3 pythondemos/svm.py, from the repository root — it writes its data files under pythondemos/. Requires NumPy and Matplotlib only, and uses fixed seeds, so the printed numbers reproduce exactly. Download svm.py
One cell per block of the script: the code, and what that code printed when it last ran here
"""
svm.py — numerical companion to the glossary entry
'support vector machine (SVM)'.
Purpose
-------
A small, hand-crafted non-separable toy trainset in R^2 whose soft-margin
SVM solution leaves exactly one data point misclassified *beyond* the
opposite margin, i.e. with functional margin y (w^T x + b) < -1 (a bounded
support vector, hinge loss > 2). The same dataset drives the entry's
non-separable figure (Fig. 2). Self-contained (numpy/matplotlib only),
fixed seed; the SVM is solved by a pure-numpy simplified SMO so no external
solver is needed.
Blocks
------
[P-nonsep] Two symmetric sets of points (class +1 right, class -1 left) plus one
outlier labelled +1 planted deep in the -1 region. The soft-margin
SVM keeps the vertical boundary between the two classes (correcting the
outlier would misclassify the whole -1 class) and yields
w_hat = (0.5, 0), b_hat = 0: boundary x1 = 0, margins x1 = +-2,
margin 1/||w_hat|| = 2.
[P-sv] The four points at x1 = +-2 are margin support vectors
(0 < alpha < C, functional margin exactly 1); the outlier is a
bounded support vector (alpha = C); the two points at x1 = +-3 are
not support vectors (alpha = 0) and can be removed without changing
w_hat.
[P-alpha] Overly large regularization: ||w|| <= rho/(2 alpha),
constant majority predictions, minority class all
misclassified support vectors.
[P-viol] The outlier x = (-4, 0), y = +1 has y (w^T x + b) = -2 < -1, hinge
loss xi = 3.
Outputs
-------
svm_points.csv : the 7 data points with columns x1,x2,label,cls,fmargin
(cls in {pos,neg,svpos,svneg,viol}) read by the entry's
pgfplots figure.
svm.png : matplotlib preview of that figure (checking only).
"""
import numpy as np
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
report = []
def check(name, ok):
report.append((name, bool(ok)))
print(f" [{'ok' if ok else 'FAIL'}] {name}")
# --------------------------------------------------------------------------
# Pure-numpy soft-margin SVM (simplified SMO, Platt 1998 / CS229 notes).
# Solves the C-SVM dual max sum_i a_i - 1/2 sum_ij a_i a_j y_i y_j <x_i,x_j>
# s.t. 0 <= a_i <= C, sum_i a_i y_i = 0. Returns w = sum_i a_i y_i x_i, b, a.
# --------------------------------------------------------------------------
def smo(X, y, C, tol=1e-6, max_passes=500, seed=0):
rng = np.random.default_rng(seed)
m = X.shape[0]
K = X @ X.T
a = np.zeros(m)
b = 0.0
passes = 0
def score(i):
return (a * y) @ K[i] + b
while passes < max_passes:
changed = 0
for i in range(m):
Ei = score(i) - y[i]
if (y[i] * Ei < -tol and a[i] < C) or (y[i] * Ei > tol and a[i] > 0):
j = int(rng.integers(m - 1))
j = j + (j >= i) # uniform j != i
Ej = score(j) - y[j]
ai, aj = a[i], a[j]
if y[i] != y[j]:
L, H = max(0.0, aj - ai), min(C, C + aj - ai)
else:
L, H = max(0.0, ai + aj - C), min(C, ai + aj)
if L == H:
continue
eta = 2 * K[i, j] - K[i, i] - K[j, j]
if eta >= 0:
continue
a[j] = np.clip(aj - y[j] * (Ei - Ej) / eta, L, H)
if abs(a[j] - aj) < 1e-12:
continue
a[i] = ai + y[i] * y[j] * (aj - a[j])
b1 = b - Ei - y[i] * (a[i] - ai) * K[i, i] - y[j] * (a[j] - aj) * K[i, j]
b2 = b - Ej - y[i] * (a[i] - ai) * K[i, j] - y[j] * (a[j] - aj) * K[j, j]
b = b1 if 0 < a[i] < C else b2 if 0 < a[j] < C else 0.5 * (b1 + b2)
changed += 1
passes = passes + 1 if changed == 0 else 0
w = (a * y) @ X
return w, b, a
# --------------------------------------------------------------------------
# Toy dataset in R^2
# --------------------------------------------------------------------------
X = np.array([
[ 2.0, 1.0], [ 2.0, -1.0], [ 3.0, 0.0], # class +1 points (right)
[-2.0, 1.0], [-2.0, -1.0], [-3.0, 0.0], # class -1 points (left)
[-4.0, 0.0], # outlier labelled +1, deep in -1 region
])
y = np.array([1.0, 1.0, 1.0, -1.0, -1.0, -1.0, 1.0])
OUT = 6 # index of the outlier
C = 1.0 # C = 1/(2 m alpha); m=7 -> alpha = 1/14
w, b, alpha = smo(X, y, C)
fmargin = y * (X @ w + b) # functional margins y_i (w^T x_i + b)
xi = np.maximum(0.0, 1.0 - fmargin) # hinge losses
print(f"w_hat = {np.round(w,4).tolist()}, b_hat = {round(b,4)}, "
f"margin 1/||w|| = {1/np.linalg.norm(w):.3f}")
print(f"functional margins: {np.round(fmargin,3).tolist()}")
print(f"alpha (0<alpha<C: margin SV; alpha=C: bounded SV): {np.round(alpha,3).tolist()}")
w_hat = [0.5, 0.0], b_hat = 0.0, margin 1/||w|| = 2.000 functional margins: [1.0, 1.0, 1.5, 1.0, 1.0, 1.5, -2.0] alpha (0<alpha<C: margin SV; alpha=C: bounded SV): [0.493, 0.132, -0.0, 0.993, 0.632, 0.0, 1.0] [ok] w_hat = (0.5, 0), b_hat = 0 (vertical boundary x1 = 0) [ok] margin 1/||w_hat|| = 2 [ok] the four points at x1 = +-2 are margin support vectors (0 < alpha < C) [ok] the two points at x1 = +-3 are not support vectors (alpha = 0) [ok] the outlier is a bounded support vector (alpha = C) [ok] removing a non-support vector leaves w_hat, b_hat unchanged [ok] outlier x=(-4,0), y=+1 has y(w^T x + b) < -1 (bounded, hinge loss > 2) outlier: y(w^T x + b) = -2.000, hinge loss xi = 3.000 wrote pythondemos/svm_points.csv
Two symmetric sets of points (class +1 right, class -1 left) plus one outlier labelled +1 planted deep in the -1 region. The soft-margin SVM keeps the vertical boundary between the two classes (correcting the outlier would misclassify the whole -1 class) and yields w_hat = (0.5, 0), b_hat = 0: boundary x1 = 0, margins x1 = +-2, margin 1/||w_hat|| = 2.
check("w_hat = (0.5, 0), b_hat = 0 (vertical boundary x1 = 0)",
np.allclose(w, [0.5, 0.0], atol=1e-3) and abs(b) < 1e-3)
check("margin 1/||w_hat|| = 2", abs(1 / np.linalg.norm(w) - 2.0) < 1e-3)
The four points at x1 = +-2 are margin support vectors (0 < alpha < C, functional margin exactly 1); the outlier is a bounded support vector (alpha = C); the two points at x1 = +-3 are not support vectors (alpha = 0) and can be removed without changing w_hat.
margin_sv = np.array([0, 1, 3, 4])
non_sv = np.array([2, 5])
check("the four points at x1 = +-2 are margin support vectors (0 < alpha < C)",
np.all((alpha[margin_sv] > 1e-6) & (alpha[margin_sv] < C - 1e-6)))
check("the two points at x1 = +-3 are not support vectors (alpha = 0)",
np.all(alpha[non_sv] < 1e-6))
check("the outlier is a bounded support vector (alpha = C)", abs(alpha[OUT] - C) < 1e-6)
# removing a non-support vector leaves w_hat unchanged
keep = [i for i in range(len(y)) if i != non_sv[0]]
w2, b2, _ = smo(X[keep], y[keep], C)
check("removing a non-support vector leaves w_hat, b_hat unchanged",
np.allclose(w2, w, atol=1e-3) and abs(b2 - b) < 1e-3)
The outlier x = (-4, 0), y = +1 has y (w^T x + b) = -2 < -1, hinge loss xi = 3.
check("outlier x=(-4,0), y=+1 has y(w^T x + b) < -1 (bounded, hinge loss > 2)",
fmargin[OUT] < -1 and xi[OUT] > 2)
print(f"outlier: y(w^T x + b) = {fmargin[OUT]:+.3f}, hinge loss xi = {xi[OUT]:.3f}")
# --------------------------------------------------------------------------
# CSV for the entry's pgfplots figure
# --------------------------------------------------------------------------
cls = np.array(["pos", "pos", "pos", "neg", "neg", "neg", "viol"], dtype=object)
cls[margin_sv[:2]] = "svpos" # +1 margin SVs (2, +-1)
cls[margin_sv[2:]] = "svneg" # -1 margin SVs (-2, +-1)
with open("pythondemos/svm_points.csv", "w") as fh:
fh.write("x1,x2,label,cls,fmargin\n")
for (x1, x2), yi, ci, fm in zip(X, y, cls, fmargin):
fh.write(f"{x1:g},{x2:g},{int(yi):+d},{ci},{fm:g}\n")
print("wrote pythondemos/svm_points.csv")
# --------------------------------------------------------------------------
# matplotlib preview (checking only)
# --------------------------------------------------------------------------
fig, ax = plt.subplots(figsize=(6.4, 3.2))
for xv, lab in [(0, r"$0$"), (-2, r"$-1$"), (2, r"$+1$")]:
ax.axvline(xv, ls="-" if xv == 0 else "--", color="0.4", lw=1.2)
ax.text(xv, 1.55, lab, ha="center", fontsize=9)
ax.text(0, 1.9, r"$\hat w^\top x + \hat b$", ha="center", fontsize=9)
pos = y > 0
ax.scatter(X[pos & (cls != "viol") & (alpha < C), 0],
X[pos & (cls != "viol") & (alpha < C), 1], marker="o", c="tab:blue", s=40, label="+1")
ax.scatter(X[~pos, 0], X[~pos, 1], marker="s", facecolors="none", edgecolors="k", s=40, label="-1")
sv = (alpha > 1e-6) & (alpha < C - 1e-6)
ax.scatter(X[sv, 0], X[sv, 1], marker="o", facecolors="none", edgecolors="r", s=180, lw=2)
ax.scatter(X[OUT, 0], X[OUT, 1], marker="o", c="tab:blue", edgecolors="r", s=120, lw=2)
# xi of the outlier: gap to its own (+1) margin line at x1 = +2
ax.annotate("", xy=(2, 0), xytext=(-4, 0),
arrowprops=dict(arrowstyle="<->", color="r", lw=1))
ax.text(-1, 0.15, r"$\xi = 3$", color="r", ha="center", fontsize=9)
ax.text(-4, -0.35, "misclassified\noutlier", color="r", ha="center", fontsize=8)
ax.text(3, 0.3, "not a\nsupport vector", ha="center", fontsize=8)
ax.set_xlim(-4.8, 3.8)
ax.set_ylim(-1.6, 2.1)
ax.set_xlabel(r"$x_1$")
ax.set_ylabel(r"$x_2$")
ax.set_title("Data generated by pythondemos/svm.py")
fig.tight_layout()
# --------------------------------------------------------------------------
Overly large regularization: ||w|| <= rho/(2 alpha), constant majority predictions, minority class all misclassified support vectors.
# underfit: ||w_hat|| <= rho/(2 alpha) (rho = largest feature norm), the
# predictions collapse to the constant majority class sign(b_hat), and
# every minority-class data point becomes a misclassified support vector
# (sparsity of the expansion is lost). Note: with the unpenalized offset
# b, NOT every data point becomes a support vector (majority points can
# sit at margin slightly above one) — only the homogeneous b = 0 SVM has
# an all-support-vector threshold.
# --------------------------------------------------------------------------
rho = np.max(np.linalg.norm(X, axis=1))
norms_alpha = []
for alpha_reg in (1.0, 10.0, 50.0):
C_a = 1.0 / (2 * len(y) * alpha_reg)
w_a, b_a, _ = smo(X, y, C_a, max_passes=2000)
norms_alpha.append(np.linalg.norm(w_a))
check(f"[P-alpha] ||w_hat|| <= rho/(2 alpha) at alpha = {alpha_reg:g}",
np.linalg.norm(w_a) <= rho / (2 * alpha_reg) + 1e-9)
check("[P-alpha] ||w_hat|| shrinks as alpha grows",
norms_alpha[0] > norms_alpha[1] > norms_alpha[2])
marg_a = y * (X @ w_a + b_a) # alpha = 50 solution
minority = y == -1 # 3 of 7 labels
check("[P-alpha] predictions collapse to the constant majority class",
np.all(np.sign(X @ w_a + b_a) == 1.0))
check("[P-alpha] every minority-class point is a misclassified "
"support vector", np.all(marg_a[minority] < 0))
check("[P-alpha] sparsity lost: more support vectors than at the "
"default C", (marg_a <= 1 + 1e-6).sum() > 5)
# --------------------------------------------------------------------------
[ok] [P-alpha] ||w_hat|| <= rho/(2 alpha) at alpha = 1 [ok] [P-alpha] ||w_hat|| <= rho/(2 alpha) at alpha = 10 [ok] [P-alpha] ||w_hat|| <= rho/(2 alpha) at alpha = 50 [ok] [P-alpha] ||w_hat|| shrinks as alpha grows [ok] [P-alpha] predictions collapse to the constant majority class [ok] [P-alpha] every minority-class point is a misclassified support vector [ok] [P-alpha] sparsity lost: more support vectors than at the default C
the primal RERM view: (1/m) sum_i hinge_i + lam ||w||^2 with
# lam = 1/(2 C m) is equivalent to the dual that SMO solves — subgradient
# descent on this nonsmooth convex objective (Pegasos-style step 1/(lam t))
# reaches the same solution.
# --------------------------------------------------------------------------
m_svm = len(y)
lam = 1.0 / (2.0 * C * m_svm)
def primal(wv, bv):
return (np.mean(np.maximum(0.0, 1.0 - y * (X @ wv + bv)))
+ lam * wv @ wv)
wp = np.zeros(2); bp = 0.0
for t in range(1, 20001):
act = y * (X @ wp + bp) < 1.0 # margin violators
gw = 2.0 * lam * wp - (y[act, None] * X[act]).sum(0) / m_svm
gb = -y[act].sum() / m_svm
step = 1.0 / (lam * t)
wp -= step * gw; bp -= step * gb
check("[P-rerm] subgradient descent on the primal reaches the dual "
"(SMO) objective value",
abs(primal(wp, bp) - primal(w, b)) < 5e-3)
check("[P-rerm] primal and dual parameter vectors agree",
np.linalg.norm(wp - w) < 0.05 and abs(bp - b) < 0.1)
# --------------------------------------------------------------------------
[ok] [P-rerm] subgradient descent on the primal reaches the dual (SMO) objective value [ok] [P-rerm] primal and dual parameter vectors agree
kernel extension via a feature map: the polynomial feature
# map phi(x) = (x1^2, sqrt(2) x1 x2, x2^2) linearizes a circular pattern
# that no linear classifier on the raw features separates, and its Gram
# matrix equals the polynomial kernel (x^T x')^2 — the kernel trick.
# --------------------------------------------------------------------------
rng_k = np.random.default_rng(7)
r_in = rng_k.uniform(0.0, 0.8, 20); a_in = rng_k.uniform(0, 2*np.pi, 20)
r_out = rng_k.uniform(1.4, 2.0, 20); a_out = rng_k.uniform(0, 2*np.pi, 20)
Xc = np.vstack([np.c_[r_in*np.cos(a_in), r_in*np.sin(a_in)],
np.c_[r_out*np.cos(a_out), r_out*np.sin(a_out)]])
yc = np.concatenate([np.ones(20), -np.ones(20)])
phi = lambda Z: np.c_[Z[:, 0]**2, np.sqrt(2)*Z[:, 0]*Z[:, 1], Z[:, 1]**2]
w_lin, b_lin, _ = smo(Xc, yc, 10.0)
acc_lin = np.mean(np.sign(Xc @ w_lin + b_lin) == yc)
w_phi, b_phi, _ = smo(phi(Xc), yc, 10.0)
acc_phi = np.mean(np.sign(phi(Xc) @ w_phi + b_phi) == yc)
print(f"[P-kernel] accuracy raw features {acc_lin:.2f} vs feature map "
f"{acc_phi:.2f}")
check("[P-kernel] no linear classifier separates the circles "
"(raw accuracy well below 1)", acc_lin < 0.8)
check("[P-kernel] the same SVM on mapped features separates them",
acc_phi == 1.0)
check("[P-kernel] Gram matrix of phi equals the polynomial kernel "
"(x^T x')^2 (kernel trick)",
np.allclose(phi(Xc) @ phi(Xc).T, (Xc @ Xc.T) ** 2, atol=1e-8))
fig.savefig("pythondemos/svm.png", dpi=110)
print("wrote pythondemos/svm.png")
n_ok = sum(ok for _, ok in report)
print(f"\n{n_ok}/{len(report)} checks passed")
assert n_ok == len(report), "some checks FAILED"
[P-kernel] accuracy raw features 0.68 vs feature map 1.00 [ok] [P-kernel] no linear classifier separates the circles (raw accuracy well below 1) [ok] [P-kernel] the same SVM on mapped features separates them [ok] [P-kernel] Gram matrix of phi equals the polynomial kernel (x^T x')^2 (kernel trick) wrote pythondemos/svm.png 19/19 checks passed
P-kernel writes when the script runs