Support Vector Machines and the Kernel Trick, in depth: maximum margin, the dual, a hand-solved XOR example, a working SMO solver and practical tuning

Maximum margin from first principles

Take labelled points x_i with labels y_i in {-1, +1}. A linear classifier predicts the sign of f(x) = w·x + b. If the data are separable, many hyperplanes separate them. Scale w and b so the closest points satisfy y_i f(x_i) = 1. The distance from the hyperplane to those points is then 1/||w||, so the margin, the empty band between the classes, has width 2/||w||. Maximising it is the same as minimising ||w||²/2 subject to y_i(w·x_i + b) ≥ 1 for every i. That is the hard-margin primal: a quadratic objective with linear constraints, so it is convex and has one global optimum.

Real data overlap, so the soft-margin SVM adds a slack ξ_i ≥ 0 to each constraint, y_i(w·x_i + b) ≥ 1 - ξ_i, and pays C ∑ ξ_i in the objective. At the optimum each slack equals max(0, 1 - y_i f(x_i)), the hinge loss, so the soft-margin SVM is regularised hinge-loss minimisation. C sets the trade: large C punishes every violation and fits the training data tightly; small C accepts violations in exchange for a wider, simpler margin. Points outside the margin contribute nothing, which is why only a few matter.

The dual and the support vectors

Introduce a multiplier α_i ≥ 0 for each margin constraint and eliminate w and b. Setting derivatives of the Lagrangian to zero gives w = ∑ α_i y_i x_i and ∑ α_i y_i = 0. Substituting back gives the dual problem: maximise ∑ α_i - (1/2) ∑_i ∑_j α_i α_j y_i y_j (x_i·x_j) subject to 0 ≤ α_i ≤ C and ∑ α_i y_i = 0. The derivation is a standard application of Lagrangian duality, and because the primal is convex with linear constraints, strong duality holds and the two problems have the same optimum.

The KKT conditions explain the name. A point with α_i = 0 is outside the margin and has no influence. A point with α_i strictly between 0 and C lies exactly on the margin, y_i f(x_i) = 1, and can be used to compute b. A point with α_i = C is inside the margin or misclassified. Only points with α_i > 0 are support vectors, and the prediction is f(x) = ∑ α_i y_i (x_i·x) + b over those points alone. Delete every other training point, retrain, and you get the same model.

The kernel trick

Look at where the data appear in the dual and in the prediction: only inside dot products. Suppose you map each point through a feature map φ into a bigger space and run the linear SVM there. You would need only φ(x)·φ(z). A kernel is a function K(x, z) that equals that dot product for some φ, computed without building φ. Substitute K for every dot product and the algorithm is unchanged, but the boundary is linear in feature space and curved in input space.

For two-dimensional x, the polynomial kernel K(x, z) = (x·z + 1)² expands to exactly φ(x)·φ(z) with φ(x) = (1, √2 x₁, √2 x₂, x₁², x₂², √2 x₁ x₂), a six-dimensional space that includes the product x₁x₂. Computing K costs one dot product in two dimensions. For degree d in n dimensions the explicit feature space grows combinatorially while K still costs O(n). The RBF kernel exp(-γ ||x - z||²) corresponds to an infinite-dimensional φ, so it can only be used through the kernel.

A valid kernel (Mercer's condition) is symmetric and yields positive semi-definite Gram matrices, which keeps the dual concave.

Worked example: XOR by hand

XOR is the classic dataset no linear classifier can separate: (1, 1) and (-1, -1) are class +1, (1, -1) and (-1, 1) are class -1. Use K(x, z) = (x·z + 1)² and a very large C so the margin is hard.

Build the Gram matrix. Each point with itself: x·x = 2, so K = 9. Opposite corners: x·z = -2, so K = 1. Adjacent corners: x·z = 0, so K = 1. By symmetry all four multipliers are equal, say a, which also satisfies ∑ α_i y_i = 0. For any point, ∑_j y_i y_j K_ij = 9 + 1 - 1 - 1 = 8, so the dual objective is 4a - (1/2)(4)(8)a² = 4a - 16a², maximised at a = 1/8. The bias is 0 by symmetry.

Now the decision function. With s = x₁ + x₂ and d = x₁ - x₂, the two positive terms give (1/8)(2s² + 2) and the two negative terms give (1/8)(2d² + 2), so f(x) = (s² - d²)/4 = x₁ x₂. The SVM has discovered that the product feature, one of the six coordinates of φ, separates XOR. The boundary is the two axes, the margins are the hyperbolas x₁x₂ = ±1, and all four points sit on them. ||w||² = 32a² = 1/2, so the margin width in feature space is 2/||w|| = 2√2.

import numpy as np
from sklearn.svm import SVC

X = np.array([[1, 1], [-1, -1], [1, -1], [-1, 1]], float)
y = np.array([1, 1, -1, -1])
m = SVC(kernel="poly", degree=2, gamma=1, coef0=1, C=1e6).fit(X, y)
print(m.dual_coef_)          # y_i * alpha_i: [[-0.125 -0.125  0.125  0.125]]
print(m.intercept_)          # [0.]
print(m.decision_function([[2, 3], [0.5, -0.5]]))   # [ 6.   -0.25] = x1 * x2

Note that scikit-learn's dual_coef_ stores y_iα_i, not α_i, and support_ lists indices in its own order.

Solving the dual with SMO

General-purpose quadratic programming solvers handle the dual but scale badly. Sequential Minimal Optimisation, introduced by John Platt in 1998 and the basis of LIBSVM, the library behind scikit-learn's SVC, uses a simple observation: the equality constraint ∑ α_i y_i = 0 means you cannot change one multiplier alone, but you can change two. For a pair, the problem is one-dimensional along a line, has a closed-form solution, and is clipped to the box [0, C]. SMO repeatedly picks a pair that violates the KKT conditions, solves it exactly, and updates the bias, until no violations larger than a tolerance remain.

import numpy as np

def rbf(A, B, gamma):
    d = (A**2).sum(1)[:, None] + (B**2).sum(1)[None, :] - 2 * A @ B.T
    return np.exp(-gamma * np.maximum(d, 0))

def smo(K, y, C, tol=1e-3, max_iter=100_000, seed=0):
    """Simplified SMO for the dual SVM. K: n x n Gram matrix, y in {-1, +1}."""
    n = len(y); a = np.zeros(n); b = 0.0
    f = np.zeros(n)                       # f[i] = sum_j a_j y_j K[j, i]
    rng = np.random.default_rng(seed)

    def step(i, j):
        nonlocal b
        if i == j: return False
        Ei, Ej = f[i] + b - y[i], f[j] + b - y[j]
        ai, aj = a[i], a[j]
        if y[i] != y[j]: L, H = max(0, aj - ai), min(C, C + aj - ai)
        else:            L, H = max(0, ai + aj - C), min(C, ai + aj)
        eta = K[i, i] + K[j, j] - 2 * K[i, j]          # curvature along the pair
        if L >= H or eta <= 1e-12: return False
        a[j] = np.clip(aj + y[j] * (Ei - Ej) / eta, L, H)
        if abs(a[j] - aj) < 1e-10: return False
        a[i] = ai + y[i] * y[j] * (aj - a[j])          # keeps sum(a * y) == 0
        di, dj = (a[i] - ai) * y[i], (a[j] - aj) * y[j]
        b1 = b - Ei - di * K[i, i] - dj * K[i, j]
        b2 = b - Ej - di * K[i, j] - dj * K[j, j]
        b = b1 if 0 < a[i] < C else b2 if 0 < a[j] < C else (b1 + b2) / 2
        f[:] += di * K[i] + dj * K[j]
        return True

    for _ in range(max_iter):
        E = f + b - y
        viol = ((y * E < -tol) & (a < C)) | ((y * E > tol) & (a > 0))
        if not viol.any(): return a, b                 # KKT holds within tol
        moved = False
        for i in rng.permutation(np.flatnonzero(viol)):
            j = int(np.argmax(np.abs(E - E[i])))       # largest step heuristic
            if step(i, j) or any(step(i, k) for k in rng.permutation(n)[:20]):
                moved = True; break
        if not moved: break
    return a, b

Predict with sign(rbf(X_new, X, gamma) @ (a * y) + b). On the XOR Gram matrix this solver returns all four multipliers at 0.125 and b = 0. On 200 random two-dimensional points labelled by whether they fall outside a circle, with RBF γ = 1 and C = 1, it found the same 60 support vectors as scikit-learn, a bias of 1.0742 against 1.0749, and identical predictions on 2,000 test points. The fallback over random partners matters: a version without it stopped early with 15 support vectors and agreed on only about 91 percent of points. Stopping rules, not update formulas, are where SVM solvers usually go wrong.

Using SVMs in practice

In practice you use LIBSVM through scikit-learn and spend your effort on preprocessing and two hyperparameters. The RBF kernel compares Euclidean distances, so a feature measured in thousands drowns one measured in fractions; scale features, inside the pipeline so cross-validation does not leak test statistics. Then search C and γ on a log grid. The default gamma='scale' uses 1 / (n_features × X.var()), a reasonable centre for the grid.

from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import StandardScaler
from sklearn.svm import SVC
from sklearn.model_selection import GridSearchCV, StratifiedKFold

pipe = make_pipeline(StandardScaler(), SVC(kernel="rbf", cache_size=1000))
grid = {"svc__C": [0.1, 1, 10, 100], "svc__gamma": ["scale", 0.001, 0.01, 0.1, 1]}
search = GridSearchCV(pipe, grid, cv=StratifiedKFold(5, shuffle=True, random_state=0),
                      scoring="roc_auc", n_jobs=-1).fit(X_train, y_train)
print(search.best_params_, search.best_score_)
print("support vectors:", search.best_estimator_[-1].n_support_)

Read the support vector count: if nearly every point is one, prediction will be slow, since it costs one kernel evaluation per support vector. Small γ gives smooth boundaries that approach a linear model; large γ makes each point an island and memorises the training set.

Cost is the main limit. Training needs kernel values for many pairs, so memory and time grow roughly between n² and n³; LIBSVM caches kernel rows within cache_size megabytes to soften it. Beyond a few tens of thousands of samples, switch approach: LinearSVC for linear problems, which uses squared hinge loss by default, or an approximate feature map such as Nystroem or RBFSampler followed by a linear model such as SGDClassifier. These trade some accuracy for training that scales with n.

Failure modes

• Unscaled features. The most common cause of a poor SVM. Scale inside the pipeline.
• Probabilities from decision values. SVC scores are distances, not probabilities. probability=True fits Platt scaling with internal five-fold cross-validation, which is slow, and its predict_proba can disagree with predict. Calibrate separately if you need probabilities; see model calibration.
• Class imbalance. The hinge loss treats each error equally, so the minority class gets ignored. Use class_weight='balanced' or tuned weights, and score with a metric that sees the minority class.
• Huge C with noisy labels. Mislabelled points become support vectors at the bound and bend the boundary around themselves.
• Invalid kernels. A custom similarity that is not positive semi-definite breaks the convexity the solver relies on; check the eigenvalues of a sample Gram matrix.
• Multiclass cost. SVC trains one-versus-one classifiers internally, so cost grows with the number of class pairs.

Trade-offs

SVMs suit small to medium dense datasets, and sparse text with a linear kernel, and their convex training is reproducible. Against gradient boosted trees, they lose on tabular data with mixed feature types and missing values, and they need scaling that trees do not. Against regularised linear models, an RBF SVM captures interactions without feature engineering but costs far more to train and predict. Pick an SVM when n is moderate, features are numeric and scaled, and you want a strong, deterministic baseline.

What to do next

1. Reproduce the XOR example by hand, then confirm it with the scikit-learn snippet.
2. Run the SMO code against SVC on a synthetic dataset and compare support vectors and predictions.
3. Put StandardScaler and SVC in one pipeline and grid-search C and gamma on a log scale.
4. Check the support vector count and the training time before choosing the model.
5. For more than a few tens of thousands of rows, benchmark LinearSVC or Nystroem plus a linear model.
6. Calibrate scores separately if a downstream system needs probabilities.