Support Vector Machines & Kernels

Suppose your data is linearly separable: there exists a hyperplane that puts all positive examples on one side and all negative examples on the other. There are infinitely many such hyperplanes. Which one should you choose?

Imagine you are drawing a road between two countries. You want the road to be as wide as possible so that small measurement errors, new arrivals near the border, and future data points are least likely to end up on the wrong side. This is the maximum-margin principle.

Formal setup: each training example is a pair (xᵢ, yᵢ) where xᵢ ∈ ℝᵈ and yᵢ ∈ {-1, +1}. A hyperplane is defined by a weight vector w ∈ ℝᵈ and bias b ∈ ℝ as {x : w·x + b = 0}. The functional margin of example i is yᵢ(w·xᵢ + b). For a correctly classified point, this is positive.

The geometric margin of a point is its Euclidean distance to the hyperplane: yᵢ(w·xᵢ + b) / ‖w‖. The margin of the dataset is the minimum geometric margin over all examples.

The SVM seeks the hyperplane that maximizes this minimum distance. The examples that sit exactly on the margin boundary — the ones that would destroy the margin if removed — are called support vectors. Everything else is irrelevant: if you deleted non-support-vector training examples, you would get exactly the same hyperplane. This sparsity is a key feature, not a coincidence.

import numpy as np

def geometric_margin(w, b, X, y):
    """Compute per-point geometric margins."""
    return y * (X @ w + b) / np.linalg.norm(w)

# For a 2D toy example:
w = np.array([1.0, 1.0])
b = -1.0
X = np.array([[0, 2], [2, 0], [-1, -1], [-2, -2]])
y = np.array([1, 1, -1, -1])

margins = geometric_margin(w, b, X, y)
print(margins)          # positives > 0, negatives > 0 if correct
print(margins.min())    # this is the dataset margin

To maximize the geometric margin, we can fix ‖w‖ = 1 and maximize minᵢ yᵢ(w·xᵢ + b). Equivalently — and more conveniently — we can require that the support vectors satisfy yᵢ(w·xᵢ + b) = 1 (we can always rescale w to achieve this) and then notice that the margin equals 2/‖w‖. Maximizing the margin is therefore equivalent to minimizing ‖w‖, or equivalently ½‖w‖².

This gives the hard-margin SVM primal:

minimize    ½ ‖w‖²
subject to  yᵢ(w·xᵢ + b) ≥ 1   for all i

This is a convex quadratic program with linear constraints. It has a unique global optimum. The objective is strictly convex (the Hessian is the identity), so there are no local minima to worry about.

The geometric picture: the two margin hyperplanes are w·x + b = +1 and w·x + b = -1. The decision boundary is w·x + b = 0, equidistant between them. The support vectors are the points that satisfy yᵢ(w·xᵢ + b) = 1 with equality — they lie on the margin lines.

Note the elegant separation: minimizing ½‖w‖² controls the complexity of the classifier (sharper w → narrower margin → more capacity), while the constraints enforce correctness. This is regularization and fit in a single clean formulation.

from scipy.optimize import minimize
import numpy as np

def svm_primal(X, y):
    """
    Solve hard-margin SVM via scipy (for illustration only;
    real code uses a QP solver like CVXOPT or sklearn).
    """
    n, d = X.shape

    def objective(params):
        w = params[:d]
        return 0.5 * np.dot(w, w)

    def constraint_i(params, i):
        w, b = params[:d], params[d]
        return y[i] * (X[i] @ w + b) - 1

    params0 = np.zeros(d + 1)
    constraints = [
        {'type': 'ineq', 'fun': constraint_i, 'args': (i,)}
        for i in range(n)
    ]
    result = minimize(objective, params0,
                      constraints=constraints, method='SLSQP')
    w_opt = result.x[:d]
    b_opt = result.x[d]
    return w_opt, b_opt

The primal is a perfectly good QP, but its solution reveals something stunning when we look at it through the lens of Lagrange duality — and this is where the kernel trick becomes possible.

Attach a non-negative Lagrange multiplier αᵢ ≥ 0 to each constraint. The Lagrangian is:

L(w, b, α) = ½‖w‖² − Σᵢ αᵢ [yᵢ(w·xᵢ + b) − 1]

Setting ∂L/∂w = 0 gives w = Σᵢ αᵢ yᵢ xᵢ. The optimal w is a linear combination of the training examples, weighted by α. Setting ∂L/∂b = 0 gives Σᵢ αᵢ yᵢ = 0.

Substituting back, the dual problem is:

maximize    Σᵢ αᵢ − ½ Σᵢ Σⱼ αᵢ αⱼ yᵢ yⱼ (xᵢ · xⱼ)
subject to  αᵢ ≥ 0 for all i
            Σᵢ αᵢ yᵢ = 0

The objective depends only on dot products xᵢ·xⱼ. This is the crucial observation.

KKT complementary slackness tells us αᵢ [yᵢ(w·xᵢ + b) − 1] = 0 for all i. This means: either αᵢ = 0 (the point is not a support vector, it contributes nothing) or yᵢ(w·xᵢ + b) = 1 (it lies exactly on the margin and IS a support vector). Most αᵢ will be exactly zero — hence the sparsity.

At test time, the prediction for a new point x is sign(w·x + b) = sign(Σᵢ αᵢ yᵢ (xᵢ·x) + b). Again, only dot products appear.

# Dual SVM with CVXOPT (the standard QP library for SVMs)
import numpy as np
import cvxopt

def svm_dual(X, y):
    n = X.shape[0]
    # Gram matrix: Q[i,j] = y_i * y_j * (x_i . x_j)
    Q = np.outer(y, y) * (X @ X.T)

    # cvxopt minimizes 0.5 x'Px + q'x
    # We negate to turn maximization into minimization
    P = cvxopt.matrix(Q, tc='d')
    q = cvxopt.matrix(-np.ones(n), tc='d')
    G = cvxopt.matrix(-np.eye(n), tc='d')   # -alpha <= 0
    h = cvxopt.matrix(np.zeros(n), tc='d')
    A = cvxopt.matrix(y.reshape(1, -1), tc='d')  # sum alpha_i y_i = 0
    b_eq = cvxopt.matrix(np.zeros(1), tc='d')

    sol = cvxopt.solvers.qp(P, q, G, h, A, b_eq)
    alphas = np.array(sol['x']).flatten()

    # Support vectors: alpha > threshold
    sv = alphas > 1e-5
    w = (alphas[sv] * y[sv] @ X[sv])  # w = sum alpha_i y_i x_i
    b = np.mean(y[sv] - X[sv] @ w)
    return w, b, alphas

Here is the insight that transforms SVMs from a clever linear classifier into a universal learning machine.

The dual objective and prediction rule depend on the data only through dot products xᵢ·xⱼ. What if instead of working with raw xᵢ, we first applied a feature map φ: ℝᵈ → ℝᴴ (where H can be enormous or even infinite) and replaced every dot product xᵢ·xⱼ with φ(xᵢ)·φ(xⱼ)?

If we can compute φ(xᵢ)·φ(xⱼ) directly — without explicitly computing φ — we get all the benefit of the high-dimensional feature space at the cost of the original computation. A kernel function K(x, z) = φ(x)·φ(z) is exactly this shortcut.

Mercer's theorem tells us when a symmetric function K is a valid kernel: if and only if the matrix [K(xᵢ, xⱼ)] is positive semi-definite for any finite set of points. Positive semi-definiteness guarantees the dual QP remains convex.

Common kernels:

• Linear: K(x,z) = x·z (no feature map; same as primal)
• Polynomial: K(x,z) = (x·z + c)ᵈ (maps to all monomials up to degree d)
• RBF / Gaussian: K(x,z) = exp(−‖x−z‖²/2σ²) (infinite-dimensional feature space!)
• Sigmoid: K(x,z) = tanh(κ x·z + θ) (not always PSD — use with care)

The RBF kernel is the most popular. Its feature space corresponds to expanding the Gaussian in a Taylor series: you get an infinite sum of polynomial feature spaces, all weighted by e^{−‖x‖²/2σ²}. Yet evaluating K takes a single call to exp.

import numpy as np

def rbf_kernel(X1, X2, sigma=1.0):
    """
    Gram matrix for RBF kernel.
    K[i,j] = exp(-||x_i - x_j||^2 / (2*sigma^2))
    """
    # Efficient via broadcast: ||a-b||^2 = ||a||^2 + ||b||^2 - 2 a.b
    sq_dists = (
        np.sum(X1**2, axis=1, keepdims=True)
        + np.sum(X2**2, axis=1)
        - 2 * X1 @ X2.T
    )
    return np.exp(-sq_dists / (2 * sigma**2))

def poly_kernel(X1, X2, degree=3, c=1.0):
    return (X1 @ X2.T + c) ** degree

# Kernel SVM prediction (after solving dual with kernel Gram matrix):
def predict_kernel(X_train, y_train, alphas, b, X_test, kernel_fn):
    K = kernel_fn(X_test, X_train)   # shape: (n_test, n_train)
    scores = K @ (alphas * y_train) + b
    return np.sign(scores)

The Radial Basis Function (Gaussian) kernel deserves special attention because it is the default choice for nonlinear SVMs and its behavior is subtle.

What it measures: K(x,z) = exp(−‖x−z‖²/2σ²) lies in [0,1]. K = 1 iff x = z; K → 0 as the points grow far apart. It is a measure of similarity: nearby points have high kernel value, distant points have kernel value near zero.

The bandwidth σ: This is the most important hyperparameter. Small σ means the Gaussian drops off quickly — only very nearby training points influence a test prediction. This creates a decision boundary that wiggles tightly around each training point (high variance, risk of overfitting). Large σ means the Gaussian is nearly flat — every training point influences every prediction. This approximates a linear classifier (high bias).

Equivalently, γ = 1/2σ² is often used as the kernel hyperparameter. Large γ = small σ = narrow Gaussian = high complexity.

Infinite feature space: The Taylor expansion of the RBF kernel is: K(x,z) = e^{−‖x‖²/2σ²} · e^{−‖z‖²/2σ²} · e^{x·z/σ²} = Σₖ (1/k!) (x·z/σ²)ᵏ · (Gaussian normalization terms)

This is a weighted sum of polynomial kernels of all degrees. The SVM in this feature space is simultaneously fitting polynomials of every degree, with higher degrees exponentially down-weighted by the 1/k! factor.

Choosing σ: Use cross-validation. A practical grid: σ ∈ {0.01, 0.1, 1, 10, 100} (or log-uniform) combined with C ∈ {0.01, 0.1, 1, 10, 100}. sklearn's GridSearchCV or RandomizedSearchCV handle this automatically.

from sklearn.svm import SVC
from sklearn.model_selection import GridSearchCV
from sklearn.preprocessing import StandardScaler
from sklearn.pipeline import Pipeline

def train_rbf_svm(X_train, y_train):
    # IMPORTANT: always scale features before SVM
    pipe = Pipeline([
        ('scaler', StandardScaler()),
        ('svm', SVC(kernel='rbf'))
    ])

    param_grid = {
        'svm__C':     [0.01, 0.1, 1, 10, 100],
        'svm__gamma': [0.001, 0.01, 0.1, 1, 'scale']
    }

    gs = GridSearchCV(pipe, param_grid, cv=5,
                      scoring='accuracy', n_jobs=-1)
    gs.fit(X_train, y_train)
    print(f"Best params: {gs.best_params_}")
    print(f"Best CV accuracy: {gs.best_score_:.3f}")
    return gs.best_estimator_

The hard-margin SVM requires the data to be linearly separable (in the original or kernel feature space). Real data is messy: classes overlap, there are outliers, and even with a kernel, perfect separation may require wildly overfitted boundaries.

The soft-margin SVM introduces slack variables ξᵢ ≥ 0, one per training example. A slack ξᵢ = 0 means the point is on the correct side of its margin (or beyond). ξᵢ ∈ (0, 1] means the point is inside the margin but still correctly classified. ξᵢ > 1 means the point is misclassified.

The modified primal is:

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

The hyperparameter C > 0 controls the bias-variance tradeoff:

• Large C: the penalty for misclassification is high. The optimizer fits the training data tightly, possibly accepting a narrow margin. Risk of overfitting.
• Small C: misclassifications are cheap. The optimizer prefers a wide margin, accepting more training errors. Risk of underfitting.

The dual of the soft-margin SVM is almost identical to the hard-margin dual, but with a box constraint: 0 ≤ αᵢ ≤ C (instead of αᵢ ≥ 0). The kernel trick applies identically.

The soft-margin SVM is equivalent to minimizing the hinge loss Σᵢ max(0, 1 − yᵢ f(xᵢ)) plus an L2 regularizer ½‖w‖². This framing connects SVMs directly to the broader family of regularized empirical risk minimizers.

import numpy as np

# Hinge loss framing (equivalent to soft-margin SVM)
def hinge_loss(w, b, X, y):
    margins = y * (X @ w + b)
    return np.maximum(0, 1 - margins).mean()

# The full soft-margin objective:
def soft_svm_objective(w, b, X, y, C):
    reg = 0.5 * np.dot(w, w)
    loss = hinge_loss(w, b, X, y)
    return reg + C * loss * len(y)  # note: sklearn uses sum not mean

# Practical usage:
from sklearn.svm import SVC, LinearSVC

# Kernel SVM (solved via dual QP)
clf_kernel = SVC(C=1.0, kernel='rbf', gamma='scale')

# Linear SVM (primal, via liblinear — much faster for text)
clf_linear = LinearSVC(C=1.0, max_iter=5000)

# C=1 is the default; always cross-validate

After solving the dual, three types of training examples emerge, classified by their Lagrange multiplier αᵢ:

1. αᵢ = 0: Non-support vectors. These points are correctly classified with margin > 1 (slack ξᵢ = 0). They contribute nothing to the decision boundary. Remove them and retrain — you get the same classifier.

2. 0 < αᵢ < C: These are the margin support vectors. They sit exactly on the margin boundary (yᵢ f(xᵢ) = 1, ξᵢ = 0). They define the hyperplane completely.

3. αᵢ = C: These are the bounded support vectors. They are either inside the margin or misclassified (ξᵢ > 0). They are also support vectors but pinned at the box constraint boundary.

This sparsity is a crucial property. In practice, the number of support vectors grows with dataset noise and the C parameter, but for clean separable data with a good kernel, a small fraction of training points become support vectors.

Practical implications:

• Prediction cost is O(n_sv · d) per test point (or O(n_sv) kernel evaluations). For large n_sv this can be slow at inference time.
• Examining the support vectors gives you insight: which examples are the hardest to classify? Which sit on the decision boundary?
• b (the bias) is recovered from any margin support vector: b = yᵢ − w·xᵢ for any i with 0 < αᵢ < C. In practice, average over all such support vectors for numerical stability.

from sklearn.svm import SVC
import numpy as np

clf = SVC(C=1.0, kernel='rbf', gamma=0.5)
clf.fit(X_train, y_train)

# Support vector indices into training set
print("Number of SVs:", clf.n_support_)      # [n_sv_class0, n_sv_class1]
print("SV indices:", clf.support_)            # indices in X_train
print("SV vectors:", clf.support_vectors_.shape)

# Dual coefficients: alpha_i * y_i for each support vector
# clf.dual_coef_ has shape (n_classes-1, n_sv)
print("Dual coefs:", clf.dual_coef_.shape)

# Decision function f(x) = sum_i (alpha_i * y_i) * K(sv_i, x) + b
scores = clf.decision_function(X_test)  # raw margins, not probabilities
preds  = np.sign(scores)                # +1 or -1

The kernel trick is not just a computational shortcut — it has a deep mathematical home: Reproducing Kernel Hilbert Spaces (RKHS).

Given a positive-definite kernel K, there exists a unique Hilbert space ℋ_K of functions f: X → ℝ such that:

1. For every x, the function K(·, x) is in ℋ_K.
2. The reproducing property holds: ⟨f, K(·, x)⟩_ℋ = f(x) for all f ∈ ℋ_K, x ∈ X.

The reproducing property is remarkable: evaluating a function f at a point x is the same as taking an inner product with K(·, x) in ℋ_K. The kernel function is the feature map.

The Representer Theorem (Kimeldorf & Wahba, 1971) is the theorem that makes kernels practical: any minimizer of a regularized empirical risk of the form min_{f ∈ ℋ_K} Σᵢ L(yᵢ, f(xᵢ)) + λ ‖f‖²_ℋ

has a finite-dimensional representation f*(x) = Σᵢ αᵢ K(xᵢ, x). Even though ℋ_K may be infinite-dimensional, the optimal solution lives in the span of the n training kernels. This is why the SVM dual has n variables (one αᵢ per training point) rather than infinitely many.

The SVM is a special case: the loss is hinge loss, the regularizer is ½‖f‖²_ℋ, and the representer theorem guarantees the dual form we derived. The RKHS perspective unifies SVMs, ridge regression (kernel ridge), Gaussian processes, and spline smoothing.

# Demonstrating the representer theorem:
# The learned function f*(x) = sum_i alpha_i * K(x_i, x)
# Kernel ridge regression (closed form; same RKHS framework as SVM)

import numpy as np

def kernel_ridge_regression(X_train, y_train, X_test,
                             kernel_fn, lam=1.0):
    """
    Solve (K + lambda*I) alpha = y, then predict K_test @ alpha.
    This is the 'kernelized' version of ridge regression.
    The representer theorem guarantees f*(x) = sum alpha_i K(x_i, x).
    """
    K = kernel_fn(X_train, X_train)     # n x n Gram matrix
    n = K.shape[0]
    alpha = np.linalg.solve(K + lam * np.eye(n), y_train)

    K_test = kernel_fn(X_test, X_train) # n_test x n
    return K_test @ alpha               # predictions

The vanilla SVM is a binary classifier. Real problems often have K > 2 classes. Two strategies dominate:

One-vs-Rest (OvR): Train K binary SVMs, where SVM k treats class k as positive and all other classes as negative. At test time, run all K classifiers and pick the class with the highest decision function score. This is sklearn's default for multiclass SVMs.

One-vs-One (OvO): Train K(K−1)/2 binary SVMs, one for each pair of classes. At test time, run all classifiers and use majority voting. Each SVM is trained on a smaller dataset (only two classes), which is faster per SVM but you have quadratically many of them. sklearn's SVC uses OvO by default.

Multiclass SVM (Crammer & Singer): Extend the margin directly to K classes with a single joint optimization. More principled but harder to implement and rarely needed in practice.

Which to use: OvR is faster at training and inference for large K. OvO can be better when classes are well-separated. For K ≤ 10 the difference is small; cross-validate.

In either strategy, the kernel trick applies identically to each binary subproblem.

from sklearn.svm import SVC
from sklearn.multiclass import OneVsRestClassifier, OneVsOneClassifier

# sklearn SVC defaults: OvO with kernel SVM
clf_ovo = SVC(kernel='rbf', C=1.0, decision_function_shape='ovo')

# OvR explicitly:
clf_ovr = OneVsRestClassifier(SVC(kernel='rbf', C=1.0))

# Linear SVC (uses OvR, much faster for text/high-d data):
from sklearn.svm import LinearSVC
clf_linear = LinearSVC(C=1.0, multi_class='ovr')  # or 'crammer_singer'

SVMs are not always the right tool. Here is a principled guide:

Use SVMs when:

• Dataset is small to medium (n < 50k rows). For large n, QP solve is O(n²)–O(n³) and memory becomes an issue (the n×n Gram matrix).
• Features are continuous and roughly standardized (after scaling).
• You need strong theoretical guarantees (unique global optimum, margin-based generalization bounds).
• The decision boundary is nonlinear but you have domain knowledge about a good kernel.
• Text classification with tf-idf features (linear SVM is often the best classifier here — bag-of-words spaces are nearly linearly separable).
• Number of features is large relative to n (e.g., genomics, text): linear SVMs work well in this regime.

Prefer neural networks when:

• You have raw perceptual data (images, audio, raw text) where feature engineering is hard.
• n > 100k and you have GPU resources.
• Transfer learning applies (use a pretrained backbone).

Prefer gradient boosting (XGBoost/LightGBM) when:

• Data is tabular with mixed feature types.
• n is in the tens of thousands range.
• You need fast training and prediction.

Always:

1. Scale your features (StandardScaler or MinMaxScaler) before SVM — kernel values depend on distances.
2. Cross-validate both C and γ (for RBF). This is non-negotiable.
3. Check if the linear kernel already gives good results before jumping to RBF.

# A practical SVM workflow
from sklearn.pipeline import Pipeline
from sklearn.preprocessing import StandardScaler
from sklearn.svm import SVC
from sklearn.model_selection import RandomizedSearchCV
from scipy.stats import loguniform
import numpy as np

pipe = Pipeline([
    ('scale', StandardScaler()),
    ('clf',   SVC(kernel='rbf', probability=False))
])

# Log-uniform search is much better than a grid for C and gamma
param_dist = {
    'clf__C':     loguniform(1e-2, 1e3),
    'clf__gamma': loguniform(1e-4, 1e1)
}

search = RandomizedSearchCV(
    pipe, param_dist, n_iter=50, cv=5,
    scoring='f1_macro', n_jobs=-1, random_state=42
)
search.fit(X_train, y_train)
print(search.best_params_, search.best_score_)