Part VI — Nonlinear Methods · Lecture 12

Support Vector Machines: Margins, Duality, and Kernels

Keep the linear boundary; change the space it lives in. Getting there requires the maximum-margin principle, Lagrangian duality, and the observation that the dual touches the data only through inner products — which is the whole kernel trick in one sentence.

By the end you can

Derive the distance from a point to a hyperplane, and turn 'maximise the margin' into a convex QP
Derive the soft-margin problem and show it is hinge loss plus \(L_2\) regularisation
Derive the Lagrangian dual, the representer form \(\vw=\sum_i\alpha_i y_i\vx_i\), and read the KKT conditions
Explain what a support vector is and why most \(\alpha_i\) are exactly zero
Derive the kernel trick from the dual, and construct the explicit feature map of a polynomial kernel
State Mercer's condition and check it numerically on a Gram matrix

Assumed

Phase A, Lecture 8: inner products, norms, projections, orthogonality
Phase A, Lecture 12: gradients, Lagrange multipliers, constrained optimisation
Phase B, Lecture 4: loss functions; Lecture 6: L2 regularisation

1Which separating hyperplane?

Take a linearly separable two-class dataset. There are infinitely many hyperplanes that separate it perfectly, and every one of them has zero training error. The perceptron returns whichever it stumbles into. Logistic regression, as Lecture 5 showed, does something worse: with perfect separation the likelihood has no finite maximiser and \(\|\vw\|\to\infty\).

Empirical risk minimisation has run out of opinions. We need a second principle to choose among the zero-error solutions, and the one that turns out to be both geometrically natural and theoretically defensible is: pick the hyperplane that is as far as possible from the nearest point of either class.

Intuition

A boundary squeezed right up against the training data is a boundary that will change if one point moves slightly. A boundary sitting in the widest available corridor has room to be wrong about where the points are. Margin is a proxy for stability, and stability is a proxy for generalisation.

Phase C makes this precise: generalisation bounds for linear classifiers depend on the margin, not on the dimension \(d\). That is why a maximum-margin classifier in a million-dimensional space is not absurd — and it is what lets the kernel trick work in an infinite-dimensional space without immediately overfitting.

The distance from a point to a hyperplane0 / 4

\[\mathcal{H}=\{\vx : \vw^\top\vx + b = 0\}\]

The hyperplane. Note \(\vw\) is normal to it: if \(\vx_1,\vx_2\in\mathcal H\) then \(\vw^\top(\vx_1-\vx_2)=0\), so \(\vw\) is orthogonal to every direction lying inside the plane. Phase A Lecture 8's orthogonality, doing the first piece of work.

\[\vx = \vx_\perp + r\,\frac{\vw}{\|\vw\|}, \qquad \vx_\perp\in\mathcal{H}\]

Decompose any point. Write \(\vx\) as its closest point on the plane plus a displacement along the unit normal. The signed distance is \(r\), and this decomposition is unique because \(\vw/\|\vw\|\) is orthogonal to the plane.

\[\vw^\top\vx + b = \vw^\top\vx_\perp + b + r\frac{\vw^\top\vw}{\|\vw\|} = 0 + r\|\vw\|\]

Substitute and use \(\vw^\top\vx_\perp+b=0\). The first two terms vanish because \(\vx_\perp\) is on the plane; \(\vw^\top\vw=\|\vw\|^2\), so the last term is \(r\|\vw\|\).

\[\boxed{\ r = \frac{\vw^\top\vx+b}{\|\vw\|}\ }, \qquad \text{signed distance}\]

Solve for \(r\). The sign tells you which side. Multiplying by the label, the geometric margin of a correctly classified point is \(\gamma_i = \frac{y_i(\vw^\top\vx_i+b)}{\|\vw\|} > 0\). The numerator alone, \(y_i(\vw^\top\vx_i+b)\), is called the functional margin — and it can be inflated arbitrarily by scaling \((\vw,b)\), which is exactly the loophole the next derivation has to close.

From 'maximise the margin' to a convex quadratic program0 / 5

\[\max_{\vw,b}\ \ \min_{i}\ \frac{y_i(\vw^\top\vx_i+b)}{\|\vw\|}\]

Write down what we want. Maximise the distance to the closest point. As stated this is horrible: a max of a min, non-convex, and — worse — the objective is invariant under \((\vw,b)\to(c\vw,cb)\) for any \(c>0\), so the solution is a whole ray rather than a point.

\[\text{normalise: } \ \min_i\ y_i(\vw^\top\vx_i+b) = 1\]

Use the scale invariance instead of fighting it. Since we may rescale freely, choose the scale so that the closest point has functional margin exactly 1. This is not an assumption or an approximation — it picks one representative from each ray, and every hyperplane has exactly one such representative.

\[\max_{\vw,b}\ \frac{1}{\|\vw\|} \quad\text{s.t.}\quad y_i(\vw^\top\vx_i+b)\ge 1 \ \ \forall i\]

Rewrite under that normalisation. With \(\min_i y_i(\vw^\top\vx_i+b)=1\) the geometric margin is exactly \(1/\|\vw\|\). Replacing the equality by \(\ge 1\) is harmless: at the optimum it must be tight for at least one point, since otherwise you could scale \(\vw\) down and improve the objective.

\[\boxed{\ \min_{\vw,b}\ \tfrac12\|\vw\|^2 \quad\text{s.t.}\quad y_i(\vw^\top\vx_i+b)\ge 1\ \ \forall i\ }\]

Maximising \(1/\|\vw\|\) is minimising \(\|\vw\|\), and minimising \(\|\vw\|\) is minimising \(\tfrac12\|\vw\|^2\) — both by strict monotonicity, so the argmin is unchanged (Phase A Lecture 2's argmax theorem again). The \(\tfrac12\) and the square are for differentiability, nothing more. What we now have is a convex quadratic objective under linear inequality constraints: a convex QP, with a unique global minimum and no local minima to worry about. The margin is \(2/\|\vw\|\) counting both sides of the corridor.

\[\text{total corridor width} = \frac{2}{\|\vw\|}\]

Sanity check the geometry. The two supporting hyperplanes are \(\vw^\top\vx+b=\pm1\); by the distance formula they sit at \(\pm 1/\|\vw\|\) from the boundary. Small \(\|\vw\|\) means a wide corridor, so minimising \(\|\vw\|^2\) is literally maximising the margin. Notice that this is the same \(\|\vw\|^2\) that Lecture 6 introduced as a ridge penalty from a Gaussian prior — two completely different stories arriving at the same term.

Every separating hyperplane has zero training error. Only one has the widest corridor.Interactive

class +1class −1your hyperplanemaximum-margin solution

orientation of your hyperplane (degrees)-14

offset of your hyperplane0.00

show the maximum-margin solution

Rotate and shift your hyperplane. As long as it separates the classes, the training error stays at zero — ERM has nothing more to say. The readout tracks your corridor width against the widest one available. The optimum is determined by only a handful of points (ringed): move any other point and nothing changes at all.

2Soft margin: what to do when no corridor exists

Real data is not separable, and even when it is, one mislabelled point can collapse the margin to nothing. Introduce a slack \(\xi_i\ge0\) per point measuring how far it falls short of margin 1, and pay for it.

The soft-margin problem is hinge loss plus ridge0 / 4

\[\min_{\vw,b,\vxi}\ \tfrac12\|\vw\|^2 + C\sum_{i=1}^n \xi_i \quad\text{s.t.}\quad y_i(\vw^\top\vx_i+b)\ge 1-\xi_i,\ \ \xi_i\ge0\]

The soft-margin QP. \(\xi_i=0\) means the point is at or beyond margin 1; \(0<\xi_i<1\) means inside the corridor but correctly classified; \(\xi_i>1\) means misclassified. \(C\) is the price of a unit of violation.

\[\text{for fixed }(\vw,b):\quad \xi_i^\star = \max\big(0,\ 1-y_i(\vw^\top\vx_i+b)\big)\]

Optimise out the slacks. Each \(\xi_i\) appears in exactly one constraint pair and contributes \(C\xi_i\ge0\) to the objective, so we make it as small as the constraints allow: at least \(0\), and at least \(1-y_i(\vw^\top\vx_i+b)\). The larger of the two, which is the max. No calculus needed — just read the constraints.

\[\boxed{\ \min_{\vw,b}\ \underbrace{\sum_{i=1}^{n}\max\big(0,\,1-y_i f(\vx_i)\big)}_{\text{hinge loss}} \;+\; \frac{1}{2C}\|\vw\|^2\ }\]

Substitute back and divide by \(C\). The SVM is not a special kind of object: it is empirical risk minimisation with the hinge loss and an \(L_2\) penalty, exactly the template of Lectures 4 and 6. \(C\) is an inverse regularisation strength — large \(C\) means little regularisation and a hard margin; \(\lambda = \frac{1}{2C}\) in Lecture 6's notation. Everything you learnt about choosing \(\lambda\) transfers.

\[\ell_{\text{hinge}}(m)=\max(0,1-m) \quad\text{vs}\quad \ell_{\text{log}}(m)=\log(1+e^{-m}), \qquad m=y f(\vx)\]

Compare to logistic regression. Both are convex upper bounds on the 0–1 loss \(\mathbb{1}[m<0]\), and both are decreasing in the margin. The difference is one property: hinge is exactly zero for \(m\ge1\), while logistic is positive everywhere. That single fact is the reason SVMs have sparse solutions and logistic regression does not — a point past the margin contributes no gradient, so it cannot influence \(\vw\). Everything about support vectors follows from the flat piece of the hinge.

Margin losses: what each one charges for being rightInteractive

hinge max(0, 1−m)logistic log(1+e^−m)/log 2squared (1−m)²0–1 loss

margin m = y·f(x)1.60

Only the hinge is exactly zero once m ≥ 1 — that flat piece is what makes most αᵢ vanish and the solution sparse. Logistic never quite reaches zero, so every point keeps a say. Squared loss grows again for m > 1: it penalises being TOO confidently right, which is why plain least squares is a poor classifier.

3The dual, and why it is worth the trouble

We could stop here and minimise hinge-plus-ridge by subgradient descent. Almost nobody does, because the dual of this problem has a property the primal does not: the data appears only inside inner products. That is the door the kernel trick walks through, so the dual is not an optional detour — it is the whole point.

The Lagrangian dual of the soft-margin SVM0 / 6

\[\mathcal{L}(\vw,b,\vxi,\valpha,\vm) = \tfrac12\|\vw\|^2 + C\sum_i \xi_i - \sum_i \alpha_i\big[y_i(\vw^\top\vx_i+b)-1+\xi_i\big] - \sum_i \mu_i\xi_i\]

Form the Lagrangian. One multiplier \(\alpha_i\ge0\) per margin constraint, one \(\mu_i\ge0\) per non-negativity constraint on \(\xi_i\). Sign convention: constraints are written as \(\ge 0\) and subtracted, so that \(\mathcal L \le\) the primal objective for any feasible point. This is Phase A Lecture 12's Lagrange multiplier machinery, now with inequalities instead of equalities.

\[\frac{\partial\mathcal L}{\partial\vw}=\vw-\sum_i\alpha_i y_i\vx_i = \vzero \;\Longrightarrow\; \boxed{\ \vw=\sum_i \alpha_i y_i \vx_i\ }\]

Stationarity in \(\vw\). This is the single most consequential line in the lecture. The optimal weight vector is a linear combination of the training points, with the labels as signs and the multipliers as coefficients. It is a special case of the representer theorem: any \(L_2\)-penalised linear model has its solution in the span of the data. Anything orthogonal to the data adds to \(\|\vw\|^2\) and changes no prediction, so the optimum never uses it.

\[\frac{\partial\mathcal L}{\partial b}=-\sum_i\alpha_i y_i = 0, \qquad \frac{\partial\mathcal L}{\partial\xi_i}=C-\alpha_i-\mu_i=0 \;\Longrightarrow\; 0\le\alpha_i\le C\]

The other two stationarity conditions. The first is an equality constraint the dual must respect. The second, combined with \(\mu_i\ge0\), box-constrains the multipliers: \(\alpha_i\in[0,C]\). Notice how the soft-margin parameter \(C\) enters the dual — as a ceiling on how much influence any single point may have. That is a precise statement of robustness: no outlier can pull the boundary by more than \(C\).

\[\boxed{\ \max_{\valpha}\ \sum_{i=1}^n \alpha_i - \tfrac12\sum_{i=1}^n\sum_{j=1}^n \alpha_i\alpha_j y_i y_j\, \vx_i^\top\vx_j \quad\text{s.t.}\quad 0\le\alpha_i\le C,\ \ \sum_i\alpha_i y_i=0\ }\]

Substitute all three conditions back. The \(\xi\) terms cancel identically (\(C\xi_i-\alpha_i\xi_i-\mu_i\xi_i=(C-\alpha_i-\mu_i)\xi_i=0\)); the \(b\) term cancels by \(\sum_i\alpha_iy_i=0\); and \(\tfrac12\|\vw\|^2-\vw^\top\vw = -\tfrac12\|\vw\|^2\) supplies the quadratic term. What is left is a concave QP in \(\valpha\) alone. Read the boxed expression carefully: \(\vx_i\) appears only as \(\vx_i^\top\vx_j\). Nowhere else. Hold that thought.

\[\alpha_i\big[y_i(\vw^\top\vx_i+b)-1+\xi_i\big]=0, \qquad \mu_i\xi_i=0\]

Complementary slackness (KKT). For each \(i\), either the multiplier is zero or the constraint is tight. Three cases: (i) \(\alpha_i=0\) — the point is strictly outside the margin and is irrelevant; (ii) \(0<\alpha_i<C\) — then \(\mu_i>0\) so \(\xi_i=0\), and the point sits exactly on the margin, \(y_if(\vx_i)=1\); (iii) \(\alpha_i=C\) — the point is inside the margin or misclassified. Cases (ii) and (iii) are the support vectors, and only they appear in \(\vw=\sum_i\alpha_iy_i\vx_i\).

\[b = y_k - \sum_i \alpha_i y_i \vx_i^\top\vx_k \quad\text{for any } k \text{ with } 0<\alpha_k<C\]

Recovering \(b\). Case (ii) gives \(y_k f(\vx_k)=1\), and multiplying by \(y_k\) (using \(y_k^2=1\)) gives \(f(\vx_k)=y_k\). Solve for \(b\). In practice average over all such \(k\) for numerical stability. Note \(b\) is recovered from inner products too — the kernel trick will survive this step.

One thing we are quoting, not deriving

Substituting the stationarity conditions gives the dual, and weak duality (dual optimum \(\le\) primal optimum) follows in two lines from the construction. What we are using without proof is strong duality: for this problem the two optima are equal, so solving the dual solves the primal exactly. It holds because the objective is convex, the constraints are affine, and Slater's condition is satisfied.

Marked as deferred, not skipped. The general theorem belongs in a convex-optimisation treatment and Phase C covers it. Meanwhile the code block below computes both optima on a real problem and reports the gap, which comes out at \(10^{-14}\) — verification is not proof, but it is not nothing.

4The kernel trick

Here is the payoff. Suppose we want a nonlinear boundary. The old way (Lecture 2) is to choose a feature map \(\vphi:\RR^d\to\RR^{D}\) and fit a linear model in the new coordinates. That works, but \(D\) grows explosively — all monomials of degree \(\le p\) in \(d\) variables is \(\binom{d+p}{p}\) features, which for \(d=100,\ p=4\) is over four million.

But the dual never touches \(\vphi(\vx_i)\) on its own. It touches only \(\vphi(\vx_i)^\top\vphi(\vx_j)\). If we can compute that inner product without ever forming \(\vphi\), the dimension \(D\) stops mattering — it can be four million, or infinite.

Kernel

A function \(k:\RR^d\times\RR^d\to\RR\) is a (positive semi-definite, or Mercer) kernel if for every finite set \(\vx_1,\dots,\vx_n\) the Gram matrix \(\mK_{ij}=k(\vx_i,\vx_j)\) is symmetric positive semi-definite.

Mercer's theorem: \(k\) is such a kernel if and only if there exists a feature map \(\vphi\) into some inner-product space with \(k(\vx,\vz)=\vphi(\vx)^\top\vphi(\vz)\). So "PSD Gram matrix for every finite sample" and "is secretly an inner product" are the same condition — you check the easy one and get the useful one.

The polynomial kernel, with its feature map written out0 / 5

\[k(\vx,\vz) = (1+\vx^\top\vz)^2, \qquad \vx,\vz\in\RR^2\]

A kernel we can expand by hand. Computing it costs one 2-dimensional dot product, one addition and one multiplication. Whatever feature space it corresponds to, we are paying \(O(d)\), not \(O(D)\).

\[(1+x_1z_1+x_2z_2)^2 = 1 + x_1^2z_1^2 + x_2^2z_2^2 + 2x_1z_1 + 2x_2z_2 + 2x_1z_1x_2z_2\]

Expand \((a+b+c)^2=a^2+b^2+c^2+2ab+2ac+2bc\) with \(a=1,\ b=x_1z_1,\ c=x_2z_2\). Every term factors into (something in \(\vx\)) × (the same thing in \(\vz\)) — which is exactly what an inner product looks like.

\[\vphi(\vx) = \big(1,\ \sqrt2\,x_1,\ \sqrt2\,x_2,\ x_1^2,\ x_2^2,\ \sqrt2\,x_1x_2\big)^\top \in \RR^6\]

Read the feature map off the expansion. Check: \(\vphi(\vx)^\top\vphi(\vz) = 1+2x_1z_1+2x_2z_2+x_1^2z_1^2+x_2^2z_2^2+2x_1x_2z_1z_2\), matching term for term. The \(\sqrt2\) factors split the cross-term coefficient of 2 evenly between the two vectors — that is where they come from, and the code block checks this to machine precision.

\[\text{a linear boundary in } \RR^6 \ \longleftrightarrow\ \text{a general conic in } \RR^2\]

What we bought. \(\vw^\top\vphi(\vx)+b=0\) with those six coordinates is \(w_4x_1^2+w_5x_2^2+\sqrt2 w_6 x_1x_2+\dots=0\): every ellipse, hyperbola and parabola in the plane. Circles and rings become linearly separable. We got all of it from one scalar formula, without ever building a 6-vector.

\[\text{degree } p, \text{ dimension } d: \quad D=\binom{d+p}{p} \ \text{ features, } \ O(d) \ \text{ to evaluate } k\]

And the scaling is the real prize. \(d=100,\ p=4\) gives \(D=4{,}598{,}126\) features but the kernel still costs 100 multiplications. The Gaussian RBF kernel \(k(\vx,\vz)=e^{-\gamma\|\vx-\vz\|^2}\) goes further: expanding the exponential in its power series produces monomials of every degree, so \(D=\infty\). You are fitting a linear classifier in an infinite-dimensional space, in finite time, and only the margin — not \(D\) — controls whether it generalises.

Kernelising the dual: substitute and stop0 / 3

\[\max_{\valpha}\ \sum_i \alpha_i - \tfrac12\sum_{i,j}\alpha_i\alpha_j y_iy_j\,k(\vx_i,\vx_j), \qquad 0\le\alpha_i\le C,\ \ \sum_i\alpha_iy_i=0\]

Replace \(\vx_i^\top\vx_j\) with \(k(\vx_i,\vx_j)\). That is the entire kernel trick. Not an approximation, not a heuristic: since \(k(\vx_i,\vx_j)=\vphi(\vx_i)^\top\vphi(\vx_j)\) for some \(\vphi\), this is literally the dual of a linear SVM in the feature space, written in coordinates we never compute.

\[f(\vx) = \vw^\top\vphi(\vx)+b = \sum_{i\in SV}\alpha_i y_i\, k(\vx_i,\vx) + b\]

Prediction also kernelises. Substitute \(\vw=\sum_i\alpha_iy_i\vphi(\vx_i)\) from the representer step and the \(\vphi\)'s pair up into kernels. Note that \(\vw\) itself may be infinite-dimensional and is never formed; the model is the list of support vectors and their \(\alpha_i\). This is a fundamentally different kind of object from a weight vector — it is nonparametric, and its size grows with the data.

\[\text{cost: } \ O(n^2) \ \text{kernel evaluations to build } \mK, \ \ O(n^2)\!-\!O(n^3) \ \text{to solve, } \ O(n_{SV}) \ \text{per prediction}\]

And the bill. Independent of \(D\) — but quadratic to cubic in \(n\). That is the trade we made: dimension became free, sample size became expensive. It is exactly why kernel SVMs dominated in the 1990s and 2000s, when \(n\) was in the thousands, and why they largely lost to boosted trees and neural networks once \(n\) reached millions. The mathematics did not get worse; the datasets got bigger.

A linear boundary in feature space, seen from the input spaceInteractive

class +1class −1decision boundarysupport vectors

kernel: 0 linear · 1 polynomial (deg 2) · 2 polynomial (deg 3) · 3 RBF3

log₁₀ γ (RBF width)-0.3

log₁₀ C (inverse regularisation)0.00

dataset: 0 rings · 1 two moons · 2 xor blobs0

The dashed contours are f = ±1, the margin in feature space bent back into the input space. Ringed points are the support vectors — everything else could be deleted without changing the model. Push γ up and the boundary fragments into islands around individual points, the support vector count climbs toward n, and you are memorising. Push C up with γ moderate and the boundary contorts to catch every point. Both knobs overfit, in visibly different ways.

Check yourself

You fit an RBF SVM and let \(\gamma\to\infty\) with \(C\) fixed and large. What happens?

Not quite. That is the \(\gamma\to0\) limit. As \(\gamma\to0\), \(e^{-\gamma\|\vx-\vz\|^2}\approx 1-\gamma\|\vx-\vz\|^2\) and the kernel becomes affine in the inner product, so the boundary does tend to linear. Large \(\gamma\) does the opposite.

Correct. \(k(\vx_i,\vx_j)=e^{-\gamma\|\vx_i-\vx_j\|^2}\to 0\) for \(i\ne j\) and \(=1\) for \(i=j\), so \(\mK\to\mI\). The dual decouples into \(n\) independent problems with solution \(\alpha_i=\min(C,1)\), giving \(f(\vx)=\sum_i \alpha_iy_ik(\vx_i,\vx)\) — a sum of vanishingly narrow bumps. Perfect training accuracy, zero generalisation. Pure memorisation, and the figure above shows it happening.

Not quite. \(\mK\to\mI\) is the best-conditioned matrix there is. Numerically it is trivial; statistically it is a disaster. Do not confuse the two — this is a generalisation failure, not a numerical one.

Not quite. SVMs are emphatically not scale-invariant — that is why you standardise features before fitting one. And \(\gamma\) is not a rescaling: it changes the effective dimension of the feature space and the smoothness of every function the model can represent.

Practical engineering note

Using kernel SVMs in practice.

Standardise your features. The RBF kernel uses a single \(\gamma\) for a Euclidean distance, so a feature measured in thousands drowns out a feature measured in units. This is not a nicety; an unstandardised RBF SVM is often catastrophically bad.
Tune \(C\) and \(\gamma\) jointly on a log grid — they interact strongly, so tuning one at a time finds the wrong optimum. Powers of ten first, then refine.
Know the scaling wall. The Gram matrix is \(n\times n\). At \(n=10^5\) that is 80 GB in double precision. Past a few tens of thousands of rows, use a linear SVM (LIBLINEAR, \(O(nd)\)) or an explicit random-feature approximation of the RBF kernel (Rahimi–Recht), which converts the kernel back into a finite feature map you can feed to any linear solver.
SVMs do not output probabilities. \(f(\vx)\) is a signed distance, not a log-odds. Platt scaling fits a logistic model to \(f\) on held-out data if you need calibrated probabilities — note that this is an admission that logistic regression was modelling the right thing all along.
Multiclass is bolted on. The derivation above is intrinsically binary. Libraries use one-vs-one (\(K(K-1)/2\) classifiers) or one-vs-rest. Compare with softmax regression in Lecture 5, which handles \(K\) classes from first principles.

Why SVMs faded, honestly. Not because the mathematics was wrong — it is some of the most beautiful in the subject, and the representer theorem and duality are permanently useful. They faded because \(O(n^2)\) memory is fatal at modern scale, because gradient-boosted trees beat them on tabular data with less tuning, and because neural networks learn the feature map instead of requiring you to choose it. That last point is the real one: a kernel is a feature map you commit to in advance. Phase D is about learning it instead.

5Solve the dual yourself

Dual coordinate ascent, KKT verification, and the kernel trick checked numericallyPython · NumPy

import numpy as np, time
np.set_printoptions(precision=4, suppress=True)
_t0 = time.time()

def dual_svm(K, y, C, iters=400):
    """Coordinate ascent on  max sum(a) - 1/2 a'Qa  subject to  0 <= a_i <= C,
       where Q_ij = y_i y_j K_ij.  The bias is folded into the kernel (K + 1),
       which drops the sum(a_i y_i) = 0 constraint -- the standard 'penalised
       intercept' variant.  Update:  a_i <- clip(a_i + (1 - y_i f_i)/Q_ii)."""
    n = len(y)
    a = np.zeros(n); f = np.zeros(n)          # f_i = sum_j a_j y_j K_ij
    d = np.diag(K).copy(); d[d < 1e-12] = 1e-12
    for _ in range(iters):
        for i in range(n):
            new = min(C, max(0.0, a[i] + (1 - y[i]*f[i])/d[i]))
            if abs(new - a[i]) > 1e-14:
                f += (new - a[i])*y[i]*K[:, i]
                a[i] = new
    return a

# ============ 1. a tiny separable problem: check every KKT condition =========
X = np.array([[3., 4.], [1., 4.], [2., 3.], [6., -1.], [7., -1.], [5., -3.]])
y = np.array([1., 1., 1., -1., -1., -1.])
Xa = np.hstack([X, np.ones((len(X), 1))])     # bias as a constant feature
K  = Xa @ Xa.T

print("--- 1. HARD MARGIN: KKT CONDITIONS, POINT BY POINT ---")
a  = dual_svm(K, y, C=1e6, iters=3000)
wb = (a*y) @ Xa; w, b = wb[:2], wb[2]         # w = sum_i a_i y_i x_i  (representer)
marg = y*(X@w + b)
print(f"  w = {w},  b = {b:.4f},  ||w|| = {np.linalg.norm(w):.4f}")
print(f"  geometric margin 1/||w|| = {1/np.linalg.norm(w):.4f}   corridor = {2/np.linalg.norm(w):.4f}")
print("   i   alpha_i    y_i f(x_i)   status")
for i in range(len(y)):
    st = "support vector (on the margin)" if a[i] > 1e-8 else "inactive (alpha = 0)"
    print(f"  {i}   {a[i]:8.4f}   {marg[i]:9.4f}   {st}")
print("  Every alpha > 0 sits at y_i f(x_i) = 1 exactly, and every alpha = 0 sits")
print("  strictly beyond it. That is complementary slackness, verified.")

# ============ 2. strong duality: the gap should be ~0 =======================
print()
print("--- 2. STRONG DUALITY (the one thing we quoted rather than proved) ---")
for C in (1e6, 1.0, 0.1):
    a  = dual_svm(K, y, C, iters=3000)
    wb = (a*y) @ Xa
    dual   = a.sum() - 0.5*a @ ((y[:, None]*y[None, :])*K) @ a
    primal = 0.5*wb@wb + C*np.maximum(0, 1 - y*(Xa@wb)).sum()
    print(f"  C={C:<8g} dual {dual:.10f}   primal {primal:.10f}   gap {abs(primal-dual):.2e}"
          f"   #SV {int((a > 1e-8).sum())}")
print("  Gap at machine precision for every C. Solving the dual solves the primal.")

# ============ 3. the representer theorem, made concrete =====================
print()
print("--- 3. THE REPRESENTER THEOREM: w LIVES IN THE SPAN OF THE DATA ---")
a  = dual_svm(K, y, 1e6, iters=3000); wb = (a*y) @ Xa
print(f"  w (augmented) = {wb}")
# Embed the same data in R^8 by padding with directions the data never uses.
# If the theorem is right, w must be exactly zero in every padded coordinate.
Xp = np.hstack([Xa, np.zeros((len(Xa), 5))])
ap = dual_svm(Xp@Xp.T, y, 1e6, iters=3000); wp = (ap*y) @ Xp
S  = np.linalg.svd(Xp, compute_uv=False)
print(f"  same data padded into R^8 (rank {int((S > 1e-10).sum())}):")
print(f"    w on the 3 used directions   = {wp[:3]}")
print(f"    w on the 5 unused directions = {wp[3:]}   -> max |.| = {np.abs(wp[3:]).max():.2e}")
print("  Any component orthogonal to every x_i would add to ||w||^2 while changing")
print("  no prediction (w.x_i is unaffected), so the optimum never contains one.")
print("  This is the representer theorem, and it is what makes kernelisation possible:")
print("  w is always a combination of the data, so w.phi(x) is always a sum of kernels.")

# ============ 4. the kernel trick, checked against the explicit feature map ==
def phi2(X):
    """The RR^6 feature map of  k(x,z) = (1 + x.z)^2  in two dimensions."""
    x1, x2 = X[:, 0], X[:, 1]
    return np.column_stack([np.ones(len(X)), np.sqrt(2)*x1, np.sqrt(2)*x2,
                            x1**2, x2**2, np.sqrt(2)*x1*x2])

rgk = np.random.default_rng(0); Z = rgk.normal(size=(6, 2))
print()
print("--- 4. THE KERNEL TRICK IS AN IDENTITY, NOT AN APPROXIMATION ---")
print(f"  max |(1 + x.z)^2  -  phi(x).phi(z)|  =  "
      f"{np.abs((1 + Z@Z.T)**2 - phi2(Z)@phi2(Z).T).max():.2e}")
print(f"  6 explicit features vs 1 scalar formula. In d=100 at degree 4 that is")
print(f"  {int(round(np.exp(sum(np.log(np.arange(101,105))) - sum(np.log(np.arange(1,5))))))}"
      f" features vs 100 multiplications.")

def rbf(A, B, g): return np.exp(-g*((A[:, None, :] - B[None, :, :])**2).sum(-1))
W = rgk.normal(size=(80, 2))
ev = np.linalg.eigvalsh(rbf(W, W, 0.8))
print(f"  RBF Gram matrix (80x80): smallest eigenvalue {ev.min():.3e} >= 0  -> PSD,")
print(f"  so by Mercer it IS an inner product in some space. Condition number {ev.max()/max(ev.min(),1e-300):.1e}")

# ============ 5. linear vs kernel on data no line can separate ==============
def rings(n, rg):
    t = rg.uniform(0, 2*np.pi, n)
    r = np.where(rg.random(n) < 0.5, 1.0, 2.2)
    return (np.column_stack([r*np.cos(t), r*np.sin(t)]) + rg.normal(0, 0.45, (n, 2)),
            np.where(r < 1.5, 1.0, -1.0))

Xr, yr = rings(100, np.random.default_rng(1))
Xs, ys = rings(2000, np.random.default_rng(2))
Ra, Sa = np.hstack([Xr, np.ones((len(Xr), 1))]), np.hstack([Xs, np.ones((len(Xs), 1))])
print()
print("--- 5. LINEAR vs RBF ON CONCENTRIC RINGS ---")
al = dual_svm(Ra@Ra.T, yr, 1.0, iters=400)
print(f"  linear kernel        test accuracy {float((np.sign(Sa@Ra.T@(al*yr)) == ys).mean()):.4f}"
      f"   SVs {int((al > 1e-8).sum())}/100")
print("  gamma   train acc   test acc   support vectors")
for g in (0.05, 0.2, 1.0, 5.0, 50.0, 300.0):
    Kt = rbf(Xr, Xr, g); ak = dual_svm(Kt, yr, 1.0, iters=400)
    tr = float((np.sign(Kt@(ak*yr)) == yr).mean())
    te = float((np.sign(rbf(Xs, Xr, g)@(ak*yr)) == ys).mean())
    print(f"  {g:6g}    {tr:.4f}     {te:.4f}    {int((ak > 1e-8).sum()):3d}/100")
print("  Underfits at gamma=0.05, best near 0.2, and by gamma=300 the training")
print("  accuracy is perfect while every point is a support vector: pure memorisation.")

# ============ 6. Vapnik's leave-one-out bound, tested =======================
print()
print("--- 6. LOO ERROR <= (number of support vectors)/n  -- Vapnik's bound ---")
Xv, yv = rings(60, np.random.default_rng(1))
Kv = rbf(Xv, Xv, 0.5); Cv = 1.0
av = dual_svm(Kv, yv, Cv, iters=300); nsv = int((av > 1e-8).sum())
loo = 0
for i in range(len(yv)):
    m = np.ones(len(yv), bool); m[i] = False
    ai = dual_svm(Kv[np.ix_(m, m)], yv[m], Cv, iters=300)
    if np.sign(Kv[i, m] @ (ai*yv[m])) != yv[i]: loo += 1
print(f"  support vectors      {nsv}/60 = {nsv/60:.4f}")
print(f"  leave-one-out error  {loo}/60 = {loo/60:.4f}")
print(f"  bound satisfied: {loo/60 <= nsv/60}")
print("  Why it holds: deleting a NON-support vector cannot change the solution at all")
print("  (it has alpha = 0 and does not appear in w), so that point is still classified")
print("  correctly when held out. Only support vectors can possibly become LOO errors.")
print()
print(f"total run time {time.time() - _t0:.1f} s")

6Exercises

★

Three points in \(\RR^2\): \(\vx_1=(1,1),\ y_1=+1\); \(\vx_2=(2,2),\ y_2=+1\); \(\vx_3=(0,0),\ y_3=-1\).

Find the maximum-margin hyperplane by inspection and geometry — do not use the dual.
What is \(\|\vw\|\) and the corridor width?
Which points are support vectors? What must \(\alpha_2\) be, and why?
Verify \(\vw=\sum_i\alpha_iy_i\vx_i\) and \(\sum_i\alpha_iy_i=0\) by solving for the \(\alpha_i\).

1. The closest opposite-class pair is \(\vx_1=(1,1)\) and \(\vx_3=(0,0)\), at distance \(\sqrt2\). The maximum-margin hyperplane is the perpendicular bisector of that segment: it passes through the midpoint \((0.5,0.5)\) with normal along \(\vx_1-\vx_3=(1,1)\). So the boundary is \(x_1+x_2=1\).

2. Normalise so the functional margins are \(\pm1\). Try \(\vw=(c,c),\ b=-c\): \(y_1(\vw^\top\vx_1+b)=2c-c=c\) and \(y_3(\vw^\top\vx_3+b)=-(0-c)=c\). Setting \(c=1\) gives \(\vw=(1,1),\ b=-1\), \(\|\vw\|=\sqrt2\). Corridor width \(2/\|\vw\|=\sqrt2\approx1.4142\) — which equals the distance between \(\vx_1\) and \(\vx_3\), as it must, since those two points lie on opposite margin hyperplanes and the segment joining them is perpendicular to both.

3. Check \(\vx_2\): \(y_2(\vw^\top\vx_2+b)=4-1=3>1\), strictly outside the margin. By complementary slackness \(\alpha_2=0\): it is not a support vector, and deleting it changes nothing. \(\vx_1\) and \(\vx_3\) both have margin exactly 1 and are the support vectors.

4. \(\sum_i\alpha_iy_i=0\) gives \(\alpha_1-\alpha_3=0\), so \(\alpha_1=\alpha_3=\alpha\). Then \(\vw=\alpha(1,1)-\alpha(0,0)=(\alpha,\alpha)\). Matching \(\vw=(1,1)\) gives \(\alpha=1\). So \(\valpha=(1,0,1)\). Check the dual objective: \(\sum\alpha_i-\tfrac12\|\vw\|^2 = 2-\tfrac12(2)=1\), and the primal is \(\tfrac12\|\vw\|^2=1\). Equal — strong duality, on a problem you can do by hand. ∎

★★

Prove that the RBF kernel \(k(\vx,\vz)=e^{-\gamma\|\vx-\vz\|^2}\) is a valid Mercer kernel, and show its feature space is infinite-dimensional.

Hints: expand \(\|\vx-\vz\|^2\); the class of PSD kernels is closed under products, sums and positive scalar multiples, and under \(k(\vx,\vz)\mapsto g(\vx)k(\vx,\vz)g(\vz)\) for any function \(g\).

Step 1 — split off the norms. \(\|\vx-\vz\|^2=\|\vx\|^2-2\vx^\top\vz+\|\vz\|^2\), so

\[ k(\vx,\vz)=\underbrace{e^{-\gamma\|\vx\|^2}}_{g(\vx)}\; e^{2\gamma\,\vx^\top\vz}\; \underbrace{e^{-\gamma\|\vz\|^2}}_{g(\vz)}. \]

By the last closure property it suffices to show \(e^{2\gamma\vx^\top\vz}\) is PSD.

Step 2 — the exponential of a kernel is a kernel. \(e^{2\gamma\vx^\top\vz}=\sum_{m=0}^{\infty}\frac{(2\gamma)^m}{m!}(\vx^\top\vz)^m\). The linear kernel \(\vx^\top\vz\) is PSD (its Gram matrix is \(\mX\mX^\top\), and \(\vv^\top\mX\mX^\top\vv=\|\mX^\top\vv\|^2\ge0\)). Products of PSD kernels are PSD (the Gram matrix of a product is the Schur/Hadamard product of the Gram matrices, and the Schur product theorem says that is PSD), so \((\vx^\top\vz)^m\) is PSD for every \(m\). Non-negative combinations of PSD kernels are PSD, and a pointwise limit of PSD matrices is PSD (each \(\vv^\top\mK_N\vv\ge0\) and the limit of non-negatives is non-negative). Hence the series is PSD. ∎

Step 3 — infinite dimension. The term \((\vx^\top\vz)^m\) has, as in the lecture's degree-2 expansion, a feature map consisting of all degree-\(m\) monomials. The series includes every \(m\ge0\) with a strictly positive coefficient \(\frac{(2\gamma)^m}{m!}\), so the concatenated feature map contains monomials of every degree, scaled by \(g(\vx)\). That is a countably infinite set of linearly independent features.

Concretely: for distinct points the RBF Gram matrix is strictly positive definite (full rank) for every \(n\) — you can check this numerically in the code block. A finite \(D\)-dimensional feature map would force rank \(\le D\), so no finite \(D\) can work.

★★

Take the soft-margin dual and consider the effect of \(C\).

What happens to the solution as \(C\to\infty\)? As \(C\to0\)?
Show that no single training point can shift \(\vw\) by more than \(C\|\vx_i\|\), and explain why this makes the soft-margin SVM robust to outliers in a way that least squares is not.
A point has \(\alpha_i=C\) exactly. What do you know about it? What do you not know?
You double every feature (\(\vx\to2\vx\)). How must \(C\) change for the solution to be equivalent?

1. \(C\to\infty\) removes the box constraint, recovering the hard-margin problem — which is infeasible if the data is not separable, so the solver diverges or the margin collapses to nothing chasing one bad point. \(C\to0\) drives every \(\alpha_i\to0\), hence \(\vw\to\vzero\): the regularisation term \(\frac{1}{2C}\|\vw\|^2\) dominates and the model degenerates to a constant. \(C\) is \(1/(2\lambda)\) from Lecture 6, and both limits behave exactly as \(\lambda\to0\) and \(\lambda\to\infty\) there.

2. From \(\vw=\sum_i\alpha_iy_i\vx_i\), point \(i\) contributes \(\alpha_iy_i\vx_i\), whose norm is \(\alpha_i\|\vx_i\|\le C\|\vx_i\|\) since \(\alpha_i\le C\). So a single point's total influence is capped, no matter how badly it is misclassified. Compare least squares: the contribution of point \(i\) to \(\vw\) is proportional to its residual, which is unbounded — one point at \(y=10^6\) moves the fit arbitrarily far. The cap comes from the hinge loss growing only linearly in the violation, and the box constraint \(\alpha_i\le C\) is that linearity showing up in the dual. This is the same robustness argument as Huber and absolute loss in Lecture 4.

3. \(\alpha_i=C\) means \(\mu_i=C-\alpha_i=0\), so \(\xi_i\) is unconstrained by complementary slackness and the margin constraint is tight: \(y_if(\vx_i)=1-\xi_i\) with \(\xi_i\ge0\). So you know the point is at or inside the margin. You do not know whether it is correctly classified: \(\xi_i<1\) means correct but inside the corridor, \(\xi_i>1\) means misclassified, and both give \(\alpha_i=C\). To distinguish them you must evaluate \(f(\vx_i)\).

4. With \(\vx\to2\vx\), taking \(\vw\to\vw/2,\ b\to b\) leaves every \(y_if(\vx_i)\) unchanged, so the hinge term is identical while \(\|\vw\|^2\) falls by a factor of 4. To keep the same balance in \(\sum_i\ell_i+\frac{1}{2C}\|\vw\|^2\) we need \(\frac{1}{2C}\) four times larger, i.e. \(C\to C/4\). In general scaling features by \(s\) requires \(C\to C/s^2\) — which is precisely why an unstandardised SVM is so sensitive, and why the good value of \(C\) found on one preprocessing pipeline is meaningless on another.

★★★

Kernelise ridge regression and compare with the SVM.

Starting from \(\min_{\vw}\|\vy-\mX\vw\|^2+\lambda\|\vw\|^2\), prove the identity \((\mX^\top\mX+\lambda\mI)^{-1}\mX^\top = \mX^\top(\mX\mX^\top+\lambda\mI)^{-1}\).
Use it to write the prediction using only inner products, and kernelise.
When is the right-hand form cheaper?
Kernel ridge regression uses all \(n\) training points at prediction time, while a kernel SVM uses only the support vectors. Where exactly does that difference come from?

1. Start from the identity \(\mX^\top(\mX\mX^\top+\lambda\mI) = (\mX^\top\mX+\lambda\mI)\mX^\top\) — expand both sides and each equals \(\mX^\top\mX\mX^\top+\lambda\mX^\top\). Both bracketed matrices are invertible for \(\lambda>0\) (they are positive definite: \(\vv^\top(\mX^\top\mX+\lambda\mI)\vv=\|\mX\vv\|^2+\lambda\|\vv\|^2>0\)). Left-multiply by \((\mX^\top\mX+\lambda\mI)^{-1}\) and right-multiply by \((\mX\mX^\top+\lambda\mI)^{-1}\):

\[(\mX^\top\mX+\lambda\mI)^{-1}\mX^\top = \mX^\top(\mX\mX^\top+\lambda\mI)^{-1}. \ \blacksquare\]

2. The ridge solution (Lecture 6) is \(\hat\vw=(\mX^\top\mX+\lambda\mI)^{-1}\mX^\top\vy = \mX^\top(\mX\mX^\top+\lambda\mI)^{-1}\vy = \mX^\top\valpha\) with \(\valpha=(\mK+\lambda\mI)^{-1}\vy\) and \(\mK=\mX\mX^\top\). So \(\hat\vw=\sum_i\alpha_i\vx_i\) — the representer theorem again, arrived at by pure linear algebra rather than by Lagrangian duality. Prediction:

\[\hat f(\vx)=\hat\vw^\top\vx = \sum_{i=1}^n \alpha_i\,\vx_i^\top\vx \ \longrightarrow\ \sum_{i=1}^n \alpha_i\, k(\vx_i,\vx), \qquad \valpha=(\mK+\lambda\mI)^{-1}\vy.\]

3. The left form inverts a \(d\times d\) matrix, the right an \(n\times n\) one. Use the right when \(n<d\) — and always when kernelising, since \(d\) may be infinite while \(n\) is whatever you collected. Cost \(O(n^3)\) versus \(O(d^3)\).

4. From the loss. Squared loss has a non-zero derivative at every point except where the residual is exactly zero, so \(\alpha_i=r_i/\lambda\ne0\) for essentially every \(i\) and the whole training set is retained. The hinge loss is flat for \(m\ge1\), so points beyond the margin have zero derivative, contribute nothing, and get \(\alpha_i=0\) exactly. Sparsity in kernel methods comes from a loss with a flat region — nothing else. The same reasoning explains why LASSO is sparse (a kink at zero in the penalty) and ridge is not: in both cases sparsity traces back to non-differentiability or flatness at the right place, which is a satisfying thing to notice across two very different models.

End of Phase B

That is the classical canon: the learning problem itself, linear and logistic regression, regularisation, honest evaluation, PCA and clustering, generative models, trees and ensembles, and kernel methods. You can now derive every one of them from a loss and a prior, fit them from scratch in NumPy, and say what each assumes and where each fails.

Two threads have run through all twelve lectures and are worth naming before you leave. First: almost every model is a loss plus a penalty, and the interesting content is in which loss and which penalty. Second: every capacity knob you have met — polynomial degree, \(\lambda\), \(k\), tree depth, number of boosting rounds, \(C\), \(\gamma\) — is the same bias–variance dial wearing a different hat.

Phase C asks the question all of this has been quietly assuming an answer to: why does minimising training error tell us anything about unseen data at all? That is statistical learning theory — concentration inequalities, VC dimension, Rademacher complexity, and margin bounds — and it will justify the several results this phase had to quote. Phase D then stops choosing feature maps by hand.

Previous11 · Decision Trees, Random Forests, and Gradient Boosting