6.7 Support vector machines

Separate two classes with the widest possible street. Only the cases on the edge of the street, the support vectors, decide where it goes. A kernel bends the street when a straight one will not do.

import matplotlib.pyplot as plt
import numpy as np
from matplotlib.axes import Axes
from sklearn import datasets
from sklearn.linear_model import LogisticRegression
from sklearn.model_selection import StratifiedKFold, cross_val_score
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import StandardScaler
from sklearn.svm import SVC

from statly import style

style.use()

iris = datasets.load_iris(as_frame=True).frame
PETALS = ["petal length (cm)", "petal width (cm)"]
cancer = datasets.load_breast_cancer()  # 569 breast tumour samples, 30 measurements
CLASS_COLORS = np.array([style.BLUE, style.CORAL])


def plot_svm(
    ax: Axes,
    model: SVC,
    points: np.ndarray,
    labels: np.ndarray,
    limits: tuple[tuple[float, float], tuple[float, float]],
    title: str,
) -> None:
    """Classes, the decision boundary (solid), the street's edges (dashed) and support vectors."""
    (x0, x1), (y0, y1) = limits
    gx, gy = np.meshgrid(np.linspace(x0, x1, 300), np.linspace(y0, y1, 300))
    score = model.decision_function(np.c_[gx.ravel(), gy.ravel()]).reshape(gx.shape)
    ax.contourf(gx, gy, score, levels=[-1e9, 0, 1e9], colors=list(CLASS_COLORS), alpha=0.08)
    ax.contour(
        gx,
        gy,
        score,
        levels=[-1, 0, 1],
        colors=style.INK,
        linestyles=["--", "-", "--"],
        linewidths=[1, 1.6, 1],
    )
    ax.scatter(points[:, 0], points[:, 1], c=CLASS_COLORS[labels], s=14, alpha=0.8)
    support = model.support_vectors_
    ax.scatter(
        support[:, 0],
        support[:, 1],
        s=70,
        facecolors="none",
        edgecolors=style.INK,
        linewidths=1,
    )
    ax.set(title=title, xlim=(x0, x1), ylim=(y0, y1))

6.7.1 The widest street

Setosa and versicolor irises are easy to tell apart by their petals: a straight line separates them with room to spare. But which line? Infinitely many lines separate them. The maximal margin classifier picks the one with the widest empty street around it.

pair = iris[iris.target < 2]
points = pair[PETALS].to_numpy()
labels = pair.target.to_numpy()
hard = SVC(kernel="linear", C=1e6).fit(
    points, labels
)  # a huge C: no flower may enter the street

fig, ax = plt.subplots(figsize=(6, 3.6))
plot_svm(ax, hard, points, labels, ((0.8, 5.4), (0, 1.9)), "The widest street")
ax.set(xlabel=PETALS[0], ylabel=PETALS[1])
plt.show()

Figure 6.7.1: Setosa (blue) and versicolor (red) by petal size. The solid line is the middle of the widest street that separates them; the dashed lines are its edges. The circled flowers on the edges are the support vectors: move any other flower (without crossing the street) and nothing changes.

The street is set by the few circled flowers on its edges. That is the main idea of the method, and its name: the cases on the edge support the street, and the others do not matter. A classifier built from a few borderline cases is robust to what happens far from the boundary.

In symbols: code the two classes as $y_i = +1$ and $-1$. A straight boundary is the set of points where $f(x) = \beta_0 + \beta_1 x_1 + \beta_2 x_2 = 0$; the class is the sign of $f(x)$. Scale the coefficients so that the edges of the street are $f(x) = +1$ and $f(x) = -1$. The street’s width is then $2 / \|\beta\|$, where $\|\beta\| = \sqrt{\beta_1^2 + \beta_2^2}$, and the widest street solves

$$\text{minimise } \|\beta\|^2 \quad \text{subject to } y_i f(x_i) \ge 1 \text{ for every case}$$

every case on the right side, at least on the edge of the street.

6.7.2 When the classes overlap

Versicolor and virginica overlap: no street is empty. The support vector classifier allows some flowers inside the street, or even on the wrong side, at a price. A knob, C, sets the price. A high C tolerates few violations and gives a narrow street held by a few flowers; a low C tolerates many and gives a wide street held by many.

pair = iris[iris.target > 0]
points = pair[PETALS].to_numpy()
labels = (pair.target == 2).astype(int).to_numpy()
LIMITS = ((2.8, 7.2), (0.8, 2.7))

fig, axes = style.subplots(2, height=3.4, width=9.5, sharey=True)
for ax, (c, title) in zip(axes, ((0.05, "Low price: C = 0.05"), (100, "High price: C = 100"))):
    model = SVC(kernel="linear", C=c).fit(points, labels)
    plot_svm(
        ax, model, points, labels, LIMITS, f"{title}, {len(model.support_)} support vectors"
    )
    ax.set_xlabel(PETALS[0])
axes[0].set_ylabel(PETALS[1])
fig.tight_layout()
plt.show()

Figure 6.7.2: Versicolor (blue) and virginica (red). Left: a low price for violations gives a wide street with many support vectors (circled). Right: a high price gives a narrow street with few. The boundary moves less than the street.

C is the flexibility knob of this method. A high C bends to every training flower near the boundary: low bias, high variance. A low C averages over many: the reverse. As always, choose it by cross-validation.

The same thing as a loss

The support vector classifier can be written like every other model in this part: as a loss on the training cases plus a penalty. Its loss is the hinge loss, $\max(0,\, 1 - y_i f(x_i))$: zero for a case safely on the right side of the street, and growing linearly the further a case strays into the street or across it. The penalty is $\lambda \|\beta\|^2$, the ridge penalty of Chapter 6.4, with λ playing the role of 1/C. Logistic regression minimises almost the same thing with a smooth loss instead:

margin = np.linspace(-3, 3, 400)

fig, ax = plt.subplots(figsize=(8, 3.2))
ax.plot(margin, np.maximum(0, 1 - margin), color=style.BLUE, lw=2, label="hinge loss (SVM)")
ax.plot(
    margin,
    np.log(1 + np.exp(-margin)) / np.log(2),
    color=style.CORAL,
    lw=2,
    label="logistic loss (scaled)",
)
ax.plot(
    margin,
    (margin < 0).astype(float),
    color=style.GREY,
    ls="--",
    lw=1.5,
    label="0-1 loss (errors)",
)
ax.axvline(1, color=style.INK, ls=":", lw=1)
ax.text(1.05, 3.3, "edge of\nthe street", fontsize=8.5, color=style.INK)
ax.legend(fontsize=8, loc="upper center")
ax.set(xlabel="y · f(x): how far on the right side", ylabel="loss for the case", ylim=(-0.1, 4))
plt.show()

Figure 6.7.3: Loss for one case as a function of y·f(x), how far it lies on the right side of the boundary (positive) or the wrong side (negative). The hinge loss of the support vector machine (blue) is exactly zero once a case is beyond the edge of the street, which is why only the cases near or inside the street matter. The logistic loss (red) never quite reaches zero, so every case has a little say. The 0-1 loss (grey) counts errors but gives no direction to improve.

The shapes are so close that the two methods usually agree, which the comparison in Section 6.7.4 confirms.

6.7.3 Kernels: bending the street

Some classes cannot be separated by any straight line, however wide the street. One fix is to add new columns, such as squares and products of the measurements, and find a straight street in that bigger space; seen back in the original space, it is curved. A kernel does this without ever making the columns: it only needs a way to measure how similar two cases are.

The solution of the street problem turns out to depend on the cases only through their inner products, $\langle x_i, x_{i'} \rangle = \sum_j x_{ij} x_{i'j}$, a measure of similarity. Replace the inner product by another similarity, a kernel $K(x_i, x_{i'})$, and the street bends:

$$f(x) = \beta_0 + \sum_{i \in S} \alpha_i\, K(x, x_i)$$

a weighted sum over the support vectors S only. Two common kernels:

• the polynomial kernel of degree d, $K = (1 + \langle x_i, x_{i'} \rangle)^d$, equivalent to adding all products of the measurements up to degree d;
• the radial kernel, $K = \exp(-\gamma \sum_j (x_{ij} - x_{i'j})^2)$, which calls two cases similar when they are close. Its knob γ sets how close “close” is.

ring_rng = np.random.default_rng(3)
angle = ring_rng.uniform(0, 2 * np.pi, 300)
radius = np.where(
    np.arange(300) < 150, ring_rng.uniform(0, 1.0, 300), ring_rng.uniform(1.3, 2.2, 300)
)
ring = np.c_[radius * np.cos(angle), radius * np.sin(angle)] + ring_rng.normal(
    0, 0.18, (300, 2)
)
ring_labels = (np.arange(300) >= 150).astype(int)
RING_LIMITS = ((-2.8, 2.8), (-2.8, 2.8))
KERNELS = [
    ("Straight", SVC(kernel="linear", C=1)),
    ("Radial, γ = 1", SVC(kernel="rbf", gamma=1, C=1)),
    ("Radial, γ = 20", SVC(kernel="rbf", gamma=20, C=1)),
]

fig, axes = style.subplots(3, height=3.3, width=10.5, sharey=True)
for ax, (title, model) in zip(axes, KERNELS):
    model.fit(ring, ring_labels)
    errors = (model.predict(ring) != ring_labels).mean()
    plot_svm(ax, model, ring, ring_labels, RING_LIMITS, f"{title}: {errors:.0%} wrong")
    ax.set(aspect="equal")
fig.tight_layout()
plt.show()

Figure 6.7.4: Simulated data: one class inside a ring of the other. Left: a straight street fails. Middle: a radial kernel with a moderate γ draws a round boundary. Right: a very large γ draws a jagged boundary that hugs single cases, and almost every case becomes a support vector: overfitting.

The γ = 20 machine makes almost no errors on the training cases and would be poor on new ones; the near-perfect training score is the warning sign, not the goal (Chapter 6.3). C and γ are tuned together by cross-validation.

6.7.4 Support vector machines in practice

FOLDS = StratifiedKFold(10, shuffle=True, random_state=1)
MODELS = [
    ("logistic regression", LogisticRegression(max_iter=2000)),
    ("linear SVM", SVC(kernel="linear", C=0.1)),
    ("radial SVM", SVC(kernel="rbf", C=3, gamma="scale")),
]
accuracy = [
    cross_val_score(
        make_pipeline(StandardScaler(), model), cancer.data, cancer.target, cv=FOLDS
    )
    for _, model in MODELS
]


def plot_accuracy(ax: Axes) -> None:
    """Accuracy in each fold, with the mean, per model."""
    for i, scores in enumerate(accuracy):
        ax.plot(scores, np.full(len(scores), i), "o", ms=5, color=style.BLUE, alpha=0.6)
        ax.plot([scores.mean()] * 2, [i - 0.3, i + 0.3], color=style.CORAL, lw=2.5)
        ax.text(1.008, i, f"{scores.mean():.1%}", va="center", fontsize=9)
    ax.set_yticks(range(len(MODELS)), [name for name, _ in MODELS])
    ax.xaxis.set_major_formatter(style.PERCENT)
    ax.invert_yaxis()
    ax.set_xticks([0.92, 0.94, 0.96, 0.98, 1.0])
    ax.set(xlabel="accuracy in each fold", xlim=(0.92, 1.03))


fig, ax = plt.subplots(figsize=(8, 2.5))
plot_accuracy(ax)
plt.show()

Figure 6.7.5: Cross-validated accuracy on the breast tumour data (all 30 measurements, standardised), in ten folds. A linear support vector classifier, a radial one and logistic regression perform almost identically.

That is typical. With a straight boundary, support vector machines and logistic regression are close relatives and usually agree. The radial kernel pays off when the boundary is truly curved and there are enough cases to find it. Some practical points:

• Standardise the measurements first: the street’s width is measured in their units.
• More than two classes: one-versus-one fits a machine for every pair of classes and lets them vote; one-versus-rest fits one machine per class against all the others and picks the class with the largest score.
• The same ideas give support vector regression, which ignores residuals smaller than a set width (a street around the line) and penalises larger ones linearly.
• A support vector machine gives a score, not a probability. If you need probabilities, logistic regression gives them directly.

6.7.5 In one picture

• A support vector machine finds the widest street between two classes: width $2/\|\beta\|$.
• Equivalently, it minimises the hinge loss plus a ridge penalty; logistic regression uses a smooth cousin of the same loss.
• Only the support vectors, the cases on or inside the street, decide where it goes.
• C prices violations: high C, narrow street, more variance.
• Kernels replace the inner product by another similarity and bend the street; the radial kernel’s γ sets how local it is.
• With straight boundaries, SVMs and logistic regression usually agree.

Next: how to explain what a fitted model does.