Support Vector Machines

Imagine you have two classes of data points that are linearly separable. There are infinitely many lines you could draw to separate them. Which one should you choose?

SVMs choose the line that's as far as possible from the nearest points of both classes. This distance is called the margin, and the line that maximizes it is the maximum margin hyperplane.

The points that lie exactly on the margin boundaries are called support vectors. These are the critical points that completely determine the decision boundary. You could remove all other points and get the same hyperplane. This is why the algorithm is called Support Vector Machine.

Let's formalize this. We have training data {(x₁, y₁), …, (xₙ, yₙ)} where xᵢ ∈ ℝᵈ and yᵢ ∈ {−1, +1}.

We want to find a hyperplane defined by:

For a point xᵢ, the signed distance to the hyperplane is:

For the hyperplane to correctly classify all points with margin γ:

The yᵢ ensures we're measuring distance on the correct side (positive for +1 class, negative for -1 class). By convention, we scale w and b so that the closest points (support vectors) satisfy yᵢ(w · xᵢ + b) = 1.

This sets γ = 1/||w||, so maximizing the margin becomes:

The (1/2) is for mathematical convenience (makes the derivative cleaner). This is a convex quadratic programming problem with linear inequality constraints.

Real data isn't always linearly separable. We need to allow some points to be on the wrong side of the margin or even misclassified. We introduce slack variables ξᵢ ≥ 0 for each point:

Now we want to minimize margin violations while still maximizing the margin:

The parameter C controls the tradeoff between margin size and violations:

Solving the primal problem directly is hard with inequality constraints. This is where Lagrange multipliers and duality theory come in. The Lagrangian function incorporates the constraints via multipliers αᵢ ≥ 0:

1# Lagrangian for soft-margin SVM
2L(w, b, ξ, α, β) = (1/2)||w||² + C * Σξᵢ 
3                  - Σαᵢ[yᵢ(w·xᵢ + b) - 1 + ξᵢ] 
4                  - Σβᵢξᵢ

To find the optimal point, we take derivatives with respect to the primal variables (w, b, ξ) and set them to zero. These are the KKT (Karush-Kuhn-Tucker) stationarity conditions:

Substituting these back into the Lagrangian gives us the dual problem:

This is beautiful because instead of optimizing over w (d-dimensional), we're optimizing over α (n-dimensional, but sparse). The dual problem only depends on inner products xᵢ·xⱼ—this will be crucial for the kernel trick.

Now we need to solve the dual problem. Traditional quadratic programming solvers are slow. The SMO algorithm, invented by John Platt, is much faster.

The key insight: if we fix all α values except two, the problem has an analytical closed-form solution. SMO repeatedly selects pairs ofα values and optimizes them jointly.

1def smo_algorithm(X, y, C, tol=1e-3, max_passes=5):
2    """Sequential Minimal Optimization for SVM."""
3    n_samples = len(y)
4    alpha = np.zeros(n_samples)
5    b = 0
6    E = -y.copy()  # Error cache
7    
8    passes = 0
9    while passes < max_passes:
10        num_changed = 0
11        
12        for i in range(n_samples):
13            # Check KKT conditions for alpha[i]
14            Ei = E[i]
15            
16            if (y[i] * Ei < -tol and alpha[i] < C) or \
17               (y[i] * Ei > tol and alpha[i] > 0):
18                
19                # Select j != i (maximize |Ei - Ej|)
20                j = select_second_alpha(i, Ei, E)
21                Ej = E[j]
22                
23                # Save old alphas
24                alpha_i_old = alpha[i]
25                alpha_j_old = alpha[j]
26                
27                # Compute bounds L and H
28                if y[i] != y[j]:
29                    L = max(0, alpha[j] - alpha[i])
30                    H = min(C, C + alpha[j] - alpha[i])
31                else:
32                    L = max(0, alpha[i] + alpha[j] - C)
33                    H = min(C, alpha[i] + alpha[j])
34                
35                if L == H:
36                    continue
37                
38                # Compute eta (second derivative)
39                eta = 2 * K[i,j] - K[i,i] - K[j,j]
40                if eta >= 0:
41                    continue
42                
43                # Update alpha[j]
44                alpha[j] = alpha_j_old - y[j] * (Ei - Ej) / eta
45                alpha[j] = np.clip(alpha[j], L, H)
46                
47                # Update alpha[i]
48                alpha[i] = alpha_i_old + y[i]*y[j]*(alpha_j_old - alpha[j])
49                
50                # Update bias and error cache
51                b = update_bias(b, alpha, i, j, Ei, Ej, ...)
52                E = update_errors(E, alpha, X, y, b)
53                
54                num_changed += 1
55        
56        if num_changed == 0:
57            passes += 1
58        else:
59            passes = 0
60    
61    return alpha, b

The analytical solution for updating α₂ is:

After updating α₂ (clipped to bounds [L, H]), we compute α₁:

1# Update alpha_1 from the constraint
2alpha_1_new = alpha_1_old + y1 * y2 * (alpha_2_old - alpha_2_new)
3
4# Update bias using KKT conditions
5b1 = b - E1 - y1*(alpha_1_new - alpha_1_old)*K[i1,i1] \
6         - y2*(alpha_2_new - alpha_2_old)*K[i1,i2]
7b2 = b - E2 - y1*(alpha_1_new - alpha_1_old)*K[i1,i2] \
8         - y2*(alpha_2_new - alpha_2_old)*K[i2,i2]
9
10if 0 < alpha_1_new < C:
11    b_new = b1
12elif 0 < alpha_2_new < C:
13    b_new = b2
14else:
15    b_new = (b1 + b2) / 2

Here's where things get really cool. Look at the dual formulation—it only uses inner productsxᵢ·xⱼ. What if we replace those inner products with a more general function K(xᵢ, xⱼ)?

If K is a valid kernel (satisfies Mercer's condition), it corresponds to an inner product in some (possibly infinite-dimensional) feature space:

The magic: we never need to compute φ(x) explicitly! The kernel computes the inner product in the transformed space directly.

Popular kernels include:

def linear_kernel(x1, x2):
    return np.dot(x1, x2)

def poly_kernel(x1, x2, d=3, r=1):
    return (r + np.dot(x1, x2)) ** d

def rbf_kernel(x1, x2, gamma=0.1):
    return np.exp(-gamma * np.linalg.norm(x1 - x2) ** 2)

SVMs have several key hyperparameters that can play a large role in how the model trains:

A general tuning strategy:

SVMs are not the flashiest models in 2025, but they remain devastatingly useful. They shine in specific situations:

Though we've built SVMs from scratch, here are various ways to extend and optimize: