Support vector machines

What it is and when you reach for it

A support vector machine is a binary classifier that picks, among all hyperplanes separating two classes, the one with the maximum margin: the largest distance to the nearest training point on either side. It sits between logistic regression and neural networks in the family tree. Like logistic regression it fits a linear decision function w·x + b, and with the kernel trick it becomes nonlinear without ever learning features the way a network does; the features are fixed by the choice of kernel and only their weights are learned. You reach for an SVM today when the dataset is small to medium, the feature space is already good (classic text features, molecular descriptors, hand-crafted signals), and you want a strong decision boundary with almost no architecture decisions. You reach for its ideas far more often than for the model itself: hinge loss, margins, and support vectors reappear throughout modern metric learning and contrastive training, which is why the derivation is still worth owning even in a deep learning world.

The math

The max-margin objective

Take labels y_i ∈ {−1, +1} and a decision function f(x) = w·x + b, predicting the sign of f. The distance from a point x_i to the hyperplane f(x) = 0 is |w·x_i + b| / ‖w‖, and the point is correctly classified when y_i(w·x_i + b) > 0. The scale of (w, b) is arbitrary: doubling both doubles every score and moves no boundary. The classical move is to spend that free scale by requiring the closest points to satisfy y_i(w·x_i + b) = 1, at which point the geometric margin equals 1/‖w‖ exactly. Maximizing the margin then becomes minimizing ‖w‖, and for differentiability we minimize ½‖w‖² subject to y_i(w·x_i + b) ≥ 1 for every i. That is the hard-margin SVM, and it is infeasible the moment the classes overlap, so the soft-margin version introduces slack variables ξ_i ≥ 0 measuring how far each point falls short of its required margin, and pays for them linearly:

minimize    ½‖w‖² + C Σᵢ ξᵢ
subject to  yᵢ(w·xᵢ + b) ≥ 1 − ξᵢ,   ξᵢ ≥ 0

C large  → violations expensive → narrow margin, fits harder
C small  → violations cheap     → wide margin, regularizes harder

Hinge loss and the unconstrained form

The constrained problem hides a plain unconstrained one. At the optimum each slack takes the smallest value the constraints allow, which is ξ_i = max(0, 1 − y_i(w·x_i + b)). Substituting that in eliminates the constraints entirely and leaves

L(w, b) = (1/n) Σᵢ max(0, 1 − yᵢ(w·xᵢ + b)) + (λ/2)‖w‖²

hinge(m) = max(0, 1 − m)   where m = yᵢ f(xᵢ) is the margin

   loss
    │\
    │ \            zero loss once the point is
    │  \           beyond the margin, not merely
    │   \          on the right side
    └────\──────────  m
    0    1

with λ = 1/(nC) absorbing the constant. This is exactly hinge loss plus L2 regularization, the same shape as every regularized empirical risk you have ever minimized, which means ordinary gradient methods apply. The hinge is the whole personality of the SVM: it charges nothing for points classified with margin at least 1, so those points exert zero force on the solution, and the boundary is determined entirely by the points that touch or violate the margin, the support vectors. Compare the logistic loss log(1 + e^−m), which never reaches zero and therefore lets every point, however confidently classified, keep tugging on the weights forever.

Subgradient descent on the primal, with one worked step

The hinge has a corner at m = 1, so the objective is convex but not differentiable. Convexity saves us: at every point there is a subgradient, any vector that supports the function from below, and subgradient descent with a decaying step converges on convex problems. For a single example the subgradient with respect to w is −y_ix_i when the margin is below 1 and 0 when it is above; at exactly 1 any value in between is valid, and choosing 0 there is the convention autograd frameworks use for max. Adding the regularizer, one step on a batch B is

g_w = λw − (1/|B|) Σ_{i ∈ B, yᵢf(xᵢ) < 1} yᵢ xᵢ
g_b =     − (1/|B|) Σ_{i ∈ B, yᵢf(xᵢ) < 1} yᵢ
w ← w − η g_w,   b ← b − η g_b

A concrete step, small enough to check by hand. Take w = (1, 0), b = 0, λ = 0.1, learning rate η = 0.1, and the single positive example x = (0.5, 1), y = +1. The score is w·x + b = 0.5, so the margin is 0.5 < 1 and the hinge loss is 0.5. The subgradient is g_w = λw − yx = (0.1, 0) − (0.5, 1) = (−0.4, −1) and g_b = −y = −1. The update gives w = (1.04, 0.1) and b = 0.1, after which the score on the same point is 1.04·0.5 + 0.1·1 + 0.1 = 0.72: the violating point pulled the hyperplane toward classifying it with more room, exactly as it should. This mean-over-violators update with λ-shrinkage is the core of Pegasos, the classic primal SVM solver, and it is precisely what autograd computes for us below.

The dual and the kernel trick, at concept level

Lagrangian duality turns the constrained problem into one over a multiplier α_i ≥ 0 per training point: maximize Σ_i α_i − ½ Σ_iΣ_j α_iα_j y_iy_j (x_i·x_j) subject to 0 ≤ α_i ≤ C and Σ α_iy_i = 0. Two structural facts make this worth knowing even if you never solve the dual yourself. First, the optimal weights are a weighted sum of training points, w = Σ α_iy_ix_i, and complementary slackness forces α_i = 0 for every point strictly outside the margin, so the sum runs only over support vectors and the model is sparse in the data. Second, and more consequential, the training points enter the dual only through pairwise dot products x_i·x_j, so replacing every dot product with a kernel function K(x_i, x_j) = φ(x_i)·φ(x_j) trains a linear SVM in the feature space of φ without ever computing φ, which may be enormous or infinite-dimensional. The decision function follows the same pattern: f(x) = Σ_i α_iy_iK(x_i, x) + b, a similarity-weighted vote of the support vectors.

The RBF kernel on actual numbers

The workhorse kernel is the radial basis function, K(x, z) = exp(−γ‖x − z‖²), whose implicit feature space is infinite-dimensional. It reads as a similarity that decays with distance at a rate set by γ. Take γ = 0.5 and three points in the plane: a = (0, 0), b = (1, 1), c = (3, 4). Then ‖a − b‖² = 2, so K(a, b) = e⁻¹ ≈ 0.368; ‖a − c‖² = 25, so K(a, c) = e^−12.5 ≈ 3.7 × 10⁻⁶; and K(a, a) = e⁰ = 1 always. Point b is a moderately similar neighbor of a, point c is effectively invisible to it. Now the decision function makes physical sense: a query point is scored by summing α_iy_i over support vectors, each discounted by this similarity, so only nearby support vectors vote and the boundary can bend around any cluster shape. The same numbers show the failure mode: raise γ to 5 and K(a, b) drops to e⁻¹⁰ ≈ 4.5 × 10⁻⁵, meaning even close neighbors stop voting for each other, every training point becomes its own island, and the model memorizes. γ is a bandwidth, and tuning it is not optional.

Implementation, twice

The primal linear SVM is a five-line loss plus a gradient loop, and both frameworks make the subgradient handling invisible: the derivative they assign to max(0, ·) at the corner is a valid subgradient, so plain autograd descent is correct. The PyTorch version leans on the optimizer; the JAX version jits a pure update step. Both operate on the whole batch as one matrix expression, no per-sample loops.

import torch

def train_linear_svm(X, y, lam=1e-2, lr=0.1, epochs=500):
    """Primal linear SVM: hinge loss + L2 by (sub)gradient descent.

    X: (n, d) float tensor.  y: (n,) tensor with values in {-1, +1}.
    Returns (w, b, losses).  Only w is regularized, never b:
    shrinking the bias would pull the hyperplane toward the origin
    for no statistical reason.
    """
    n, d = X.shape
    w = torch.zeros(d, requires_grad=True)
    b = torch.zeros((), requires_grad=True)
    opt = torch.optim.SGD([w, b], lr=lr)
    losses = []
    for _ in range(epochs):
        opt.zero_grad()
        margins = y * (X @ w + b)                  # (n,)
        # clamp's gradient at the corner (margin exactly 1) is 0,
        # which is a valid subgradient of the hinge.
        hinge = torch.clamp(1.0 - margins, min=0.0).mean()
        loss = hinge + 0.5 * lam * w.dot(w)
        loss.backward()
        opt.step()
        losses.append(loss.item())
    return w.detach(), b.detach(), losses

import jax
import jax.numpy as jnp

def svm_loss(params, X, y, lam):
    """Hinge + L2.  y in {-1, +1}.  Only w is regularized, not b."""
    w, b = params
    margins = y * (X @ w + b)                      # (n,)
    hinge = jnp.maximum(0.0, 1.0 - margins).mean()
    return hinge + 0.5 * lam * jnp.dot(w, w)

@jax.jit
def step(params, X, y, lam, lr):
    # value_and_grad gives the loss and a valid subgradient of the
    # hinge (0 at the corner) in one traced pass.
    loss, grads = jax.value_and_grad(svm_loss)(params, X, y, lam)
    params = jax.tree.map(lambda p, g: p - lr * g, params, grads)
    return params, loss

def train_linear_svm(X, y, lam=1e-2, lr=0.1, epochs=500):
    params = (jnp.zeros(X.shape[1]), jnp.array(0.0))
    losses = []
    for _ in range(epochs):
        params, loss = step(params, X, y, lam, lr)
        losses.append(float(loss))
    w, b = params
    return w, b, losses

The kernel machine needs no new training code here because prediction is where the structure lives: given multipliers α_i from any dual solver (scikit-learn will hand you its own, see the verification recipe below), the decision function is one Gram matrix times one vector. Both versions compute the squared distances with the expansion ‖a − b‖² = ‖a‖² + ‖b‖² − 2a·b so the whole Gram matrix is a single matmul plus broadcasting, and both clamp tiny negative distances that the cancellation can produce in float32.

import torch

def rbf_gram(A, B, gamma):
    """K[i, j] = exp(-gamma * ||A_i - B_j||^2), one matmul, no loops."""
    sq = (A * A).sum(1)[:, None] + (B * B).sum(1)[None, :] - 2.0 * A @ B.T
    return torch.exp(-gamma * sq.clamp(min=0.0))   # clamp: float32 cancellation

def kernel_decision(X_test, X_sv, alpha_y, b, gamma):
    """f(x) = sum_i alpha_i y_i K(x_i, x) + b over support vectors only.

    alpha_y: (n_sv,) the products alpha_i * y_i (sklearn's dual_coef_).
    """
    K = rbf_gram(X_test, X_sv, gamma)              # (m, n_sv)
    return K @ alpha_y + b

def kernel_predict(X_test, X_sv, alpha_y, b, gamma):
    return torch.sign(kernel_decision(X_test, X_sv, alpha_y, b, gamma))

import jax
import jax.numpy as jnp

def rbf_gram(A, B, gamma):
    """K[i, j] = exp(-gamma * ||A_i - B_j||^2), one matmul, no loops."""
    sq = (A * A).sum(1)[:, None] + (B * B).sum(1)[None, :] - 2.0 * A @ B.T
    return jnp.exp(-gamma * jnp.maximum(sq, 0.0))  # float32 cancellation guard

@jax.jit
def kernel_decision(X_test, X_sv, alpha_y, b, gamma):
    """f(x) = sum_i alpha_i y_i K(x_i, x) + b over support vectors only.

    alpha_y: (n_sv,) the products alpha_i * y_i (sklearn's dual_coef_).
    """
    K = rbf_gram(X_test, X_sv, gamma)              # (m, n_sv)
    return K @ alpha_y + b

def kernel_predict(X_test, X_sv, alpha_y, b, gamma):
    return jnp.sign(kernel_decision(X_test, X_sv, alpha_y, b, gamma))

Using it on a real shape of problem

A realistic linear-SVM shape is a few thousand points in a few dozen dimensions with partial overlap. The snippet below builds two Gaussian blobs whose means are close enough that a few percent of points land on the wrong side, then trains the primal solver from the previous section. Expect the loss to drop steeply for the first tens of epochs and then flatten to a positive floor: hinge loss does not go to zero on overlapping data, because the overlapping points buy permanent slack. Accuracy should land in the mid-to-high 90s on this configuration, and the fraction of points with margin below 1, the active set, should shrink to roughly the overlap fraction. Exact numbers vary with the seed and machine, but the shapes of those curves should not.

import torch

torch.manual_seed(0)
n, d = 2000, 20
X_pos = torch.randn(n // 2, d) + 0.45     # overlapping blobs
X_neg = torch.randn(n // 2, d) - 0.45
X = torch.cat([X_pos, X_neg])
y = torch.cat([torch.ones(n // 2), -torch.ones(n // 2)])

w, b, losses = train_linear_svm(X, y, lam=1e-2, lr=0.1, epochs=300)

acc = ((X @ w + b).sign() == y).float().mean()
active = ((y * (X @ w + b)) < 1.0).float().mean()
print(f"loss {losses[0]:.3f} -> {losses[-1]:.3f}, "
      f"acc {acc:.3f}, margin violations {active:.3f}")
# loss falls fast then plateaus above zero; acc ~0.95+ (seed-dependent)

import jax
import jax.numpy as jnp

key = jax.random.PRNGKey(0)
k1, k2 = jax.random.split(key)
n, d = 2000, 20
X_pos = jax.random.normal(k1, (n // 2, d)) + 0.45  # overlapping blobs
X_neg = jax.random.normal(k2, (n // 2, d)) - 0.45
X = jnp.concatenate([X_pos, X_neg])
y = jnp.concatenate([jnp.ones(n // 2), -jnp.ones(n // 2)])

w, b, losses = train_linear_svm(X, y, lam=1e-2, lr=0.1, epochs=300)

acc = jnp.mean(jnp.sign(X @ w + b) == y)
active = jnp.mean(y * (X @ w + b) < 1.0)
print(f"loss {losses[0]:.3f} -> {losses[-1]:.3f}, "
      f"acc {acc:.3f}, margin violations {active:.3f}")
# loss falls fast then plateaus above zero; acc ~0.95+ (seed-dependent)

Applications

The SVM's reign was text. Through the 2000s, linear SVMs on sparse bag-of-words or tf-idf features were the standard for spam filtering, topic classification, and sentiment analysis, following Joachims' demonstration that high-dimensional sparse text is close to linearly separable and margins handle it beautifully; LIBLINEAR made training on millions of documents a matter of seconds. Computer vision leaned on SVMs almost as heavily before deep learning: the Dalal and Triggs pedestrian detector was HOG features under a linear SVM, deformable part models scored parts with SVMs, and even at the dawn of the deep era the original R-CNN classified its CNN features with per-class SVMs rather than a softmax head. Bioinformatics used kernel SVMs on string and graph kernels for protein and gene classification, and RankSVM turned the margin objective into one of the first practical learning-to-rank methods.

The model faded; the margin did not. Metric learning is the clearest heir: the triplet loss that trained FaceNet is max(0, d(a, p) − d(a, n) + m), a hinge demanding that the anchor-positive distance beat the anchor-negative distance by a margin m, and it inherits the SVM's key behavior that satisfied triplets contribute zero gradient, so training is driven by hard examples, the support vectors of embedding space. Face recognition heads like ArcFace build an explicit angular margin into softmax for the same reason the SVM wants a cushion: separation with slack generalizes better than bare separation. Contrastive representation learning keeps the same shape in its margin-based losses, and even where InfoNCE replaces the hinge with a softmax, the practice of hard-negative mining is support-vector thinking under a new name. The SVM survives less as a model than as a design principle: charge nothing for easy examples, spend the entire gradient budget on the hard ones at the boundary.

Against the real libraries

The lineage of production SVMs runs through two libraries from Chih-Jen Lin's group at National Taiwan University. LIBSVM solves the kernel dual with sequential minimal optimization, an algorithm that repeatedly picks the pair of multipliers most violating the optimality conditions and solves for that pair analytically, with shrinking heuristics that freeze bounded multipliers and a cache for kernel rows. It is the engine behind scikit-learn's SVC. LIBLINEAR is the linear-only sibling: it drops the kernel machinery and solves the linear dual by coordinate descent (or the primal by trust-region Newton), one multiplier at a time with a closed-form update, which scales linearly in the data and is why LinearSVC trains on millions of sparse documents while SVC, whose SMO is between quadratic and cubic in n, becomes impractical somewhere around 10⁴ to 10⁵ samples. The third option in scikit-learn, SGDClassifier with hinge loss, is the same primal subgradient method implemented on this page, plus learning-rate schedules and averaging.

What the libraries add over the reference implementation is mostly exactness and scale, not different math: SMO and coordinate descent converge to the true optimum with certificates from the duality gap, handle the box constraints natively, and exploit sparsity in ways a dense matmul cannot. The from-scratch primal solver is genuinely enough when the problem is linear, dense, and moderate-sized, or when you want the SVM loss inside a larger differentiable system where a dual solver cannot go. To verify it, train sklearn's LinearSVC(loss="hinge", C=1/(lam*n)) on the same data; the objectives differ only by the constant factor λ, so the minimizers coincide, and the learned w and b should agree to two or three decimals once both are converged (LIBLINEAR also regularizes the intercept, a small systematic difference worth knowing about). For the kernel path, fit SVC(kernel="rbf"), pull out clf.support_vectors_, clf.dual_coef_ (which is exactly the α_iy_i vector), and clf.intercept_, feed them to the kernel_decision function above, and check it matches clf.decision_function to about 10⁻⁵: that single test confirms the Gram computation, the sign conventions, and the bias handling all at once.

And when should you not use an SVM at all? Logistic regression beats it in practice more often than the textbooks suggest. If you need probabilities, logistic regression produces calibrated ones natively, while SVC's probability=True bolts on Platt scaling via an internal five-fold cross-validation that is expensive and can disagree with the decision function near the boundary. On large sparse problems the accuracy of a well-regularized logistic model and a linear SVM is usually indistinguishable, and the logistic loss is smooth, so quasi-Newton solvers converge fast and reliably. The honest summary: kernels on small weird data and margin intuition are the SVM's remaining edges; for a plain linear classifier at scale, logistic regression with the machinery on its own page is the default.

Traps and misconceptions

SVM scores are not probabilities. The decision function is a signed distance in feature space, and its magnitude is set by the margin convention, not by likelihood. Thresholding it works; feeding it to anything expecting a probability does not, and Platt scaling is a patch, not a property of the model.

Skipping feature scaling, especially with RBF. The RBF kernel is a function of Euclidean distance, so a feature measured in thousands drowns one measured in units, and a single γ cannot be right for both. Standardize features first; with kernels this is not a nicety but a correctness issue, and the worked numbers above show how fast e^−γd² collapses when distances inflate.

"Max margin means it cannot overfit." The margin argument holds for a fixed feature space, but C and γ move the effective capacity enormously: large C with large γ will carve an island around every training point and generalize terribly. Grid-search both on a log scale; the good region is usually a diagonal band, not a point.

Regularizing the bias. The λ term belongs on w only. Penalizing b drags the hyperplane toward the origin, which encodes nothing about the data; it merely biases the intercept. The implementations above exclude b deliberately, and this is also the largest source of small discrepancies against LIBLINEAR, which by default folds the intercept into the regularized weights.

Thinking the kernel trick computes features. Nothing is ever lifted into the feature space; the trick is that the optimization and the prediction only ever needed inner products, and the kernel supplies those directly. This is also why kernel SVMs must keep their support vectors around forever: the model is those points, and prediction cost grows with their number, unlike the linear case where everything collapses into one weight vector.