Theme

Blog · Classical ML ·

Support vectors: why only a handful of your data points matter

The support vector machine sells itself on accuracy. Its real trick is sparsity: after training you can throw almost all your data away, and the Karush–Kuhn–Tucker conditions say exactly which points you have to keep.

  • Interactive
  • svm
  • kernel-methods
  • classification
  • libsvm
  • iris

The report this comes out of opens Part A by contrasting two kinds of model. A parametric one fits a fixed set of weights and then the training data can be discarded; a kernel method keeps the training data around, because prediction is written in terms of it. That sounds like a straight loss: you traded a few dozen numbers for the whole dataset.

The support vector machine is the answer to that complaint, and it is a much better answer than “accuracy”. After training an SVM you do throw almost all your data away. Not as a heuristic, not as a compression step you bolt on afterwards. It falls out of the optimality conditions. Most of your points end up with a coefficient of exactly zero, and a point with a coefficient of zero contributes nothing to any future prediction. On the easiest of the three Iris problems below, with a linear kernel and C=64C = 64, 4 of 100 points survive. The other 96 could be deleted and the classifier would be bit-for-bit identical.

This post is the derivation of why that happens, the solver I have since written to check it, and a playground where you can drag a point around and watch it stop mattering.

The margin

Write the decision boundary as

y(x)=wTϕ(x)+by(\mathbf{x}) = \mathbf{w}^{\mathsf T}\boldsymbol\phi(\mathbf{x}) + b

with ϕ\boldsymbol\phi some feature map and targets tn{1,+1}t_n \in \lbrace -1, +1 \rbrace. If every point is classified correctly then tny(xn)>0t_n\,y(\mathbf{x}_n) > 0 for all nn, and the perpendicular distance from a point to the boundary is

tn(wTϕ(xn)+b)w\frac{t_n \left( \mathbf{w}^{\mathsf T}\boldsymbol\phi(\mathbf{x}_n) + b \right)}{\lVert \mathbf{w} \rVert}

The margin is the smallest of those distances. Maximising it is the whole idea: of all the boundaries that separate the classes, take the one that sits in the middle of the widest empty corridor, because it is the one least likely to be on the wrong side of a point you have not seen yet.

marginy = −1y = 0y = 1(a) the widest corridormargin(b) also correct, and worse
Both boundaries get every point right. The one on the right is more likely to be wrong about the next point, because there is less room either side of it. This is my redraw of Fig. 16 on p. 17 of the report, which was a hand sketch.

Here is the trick that turns “maximise the smallest distance” into something a solver can eat. Scaling w\mathbf{w} and bb by the same constant does not move the boundary and does not change any perpendicular distance: the scale cancels in the ratio above. So you are free to choose the scale, and the convenient choice is to make y(x)=1y(\mathbf{x}) = 1 at the closest correctly classified point. Every point then satisfies

tn(wTϕ(xn)+b)1t_n \left( \mathbf{w}^{\mathsf T}\boldsymbol\phi(\mathbf{x}_n) + b \right) \ge 1

and the margin is exactly 1/w1/\lVert\mathbf{w}\rVert. Maximising it is minimising 12w2\tfrac{1}{2}\lVert\mathbf{w}\rVert^2, which is a quadratic objective under linear constraints: a convex problem with one optimum and no local minima to get stuck in.

The dual, where the data disappears into a kernel

Attach a Lagrange multiplier an0a_n \ge 0 to each constraint and eliminate w\mathbf{w} and bb, and the problem turns inside out. Instead of minimising over the weights you maximise over the multipliers:

n=1Nan12n=1Nm=1Nanamtntmk(xn,xm)\sum_{n=1}^{N} a_n - \frac{1}{2}\sum_{n=1}^{N}\sum_{m=1}^{N} a_n a_m t_n t_m\, k(\mathbf{x}_n, \mathbf{x}_m)

subject to an0a_n \ge 0 and nantn=0\sum_n a_n t_n = 0. Two things changed, and both matter.

The first is that ϕ\boldsymbol\phi is gone. The features only ever appear as inner products ϕ(xn)Tϕ(xm)\boldsymbol\phi(\mathbf{x}_n)^{\mathsf T}\boldsymbol\phi(\mathbf{x}_m), which you can name k(xn,xm)k(\mathbf{x}_n, \mathbf{x}_m) and compute directly. You never have to build ϕ(x)\boldsymbol\phi(\mathbf{x}), which is what lets a radial basis kernel work in an infinite-dimensional feature space on a laptop from 2018.

The second is that predictions are now written in terms of the data:

y(x)=nantnk(x,xn)+by(\mathbf{x}) = \sum_{n} a_n t_n\, k(\mathbf{x}, \mathbf{x}_n) + b

which looks like exactly the disaster the introduction warned about. Every prediction is a sum over the whole training set. Except that it is not, and the reason is the next section.

Karush–Kuhn–Tucker, and where the sparsity comes from

A constrained optimum has to satisfy complementary slackness: for each constraint, either the multiplier is zero or the constraint is tight. Here that reads

an(tny(xn)1)=0a_n \left( t_n y(\mathbf{x}_n) - 1 \right) = 0

so for every single point, one of two things is true. Either tny(xn)=1t_n y(\mathbf{x}_n) = 1 (the point sits exactly on one of the two margin lines) or an=0a_n = 0, and the point contributes nothing at all to the sum above.

That is the entire sparsity argument. It is not a bound, an approximation or an empirical observation; it is the optimality condition. Points that sit on the margin are the support vectors, and they are the only rows of your dataset the trained model refers to. Everything else can be deleted.

It is also a strange thing to be told. The natural mental model of a classifier is that it somehow averages over all the data; here, a point sitting comfortably in the middle of its own class has precisely zero influence, and could be moved anywhere on its own side of the margin without the boundary twitching. That claim is much easier to believe once you have grabbed a point and tried it.

InteractiveDrag a point and watch it stop mattering
4 of these 16 points determine the boundary.The 12 faded ones have α = 0.

With JavaScript on, this becomes a canvas you can click to add points to, drag them around, and delete from, with the decision boundary, both margins and the support vectors updating live, a choice of all four LIBSVM kernels, a slider for C, and a miniature version of the report’s own (C, γ) grid search.

Two buttons under the plot make the point directly. Wander a non-support vector picks the most comfortably classified point on the plot and walks it in a circle, refusing any step that would take it inside the margin, and reports the largest change in y(x)y(\mathbf{x}) anywhere on the canvas. It stays at 0.0000. Nudge a support vector pushes a ringed point across the margin instead, and the same number jumps immediately.

That the first number is exactly zero rather than merely small is worth a footnote. There is nothing random in the solver: every iteration takes the maximal violating pair, so the same points always produce the same sequence of updates, and a point with an=0a_n = 0 that stays outside the margin is never the violator, never enters an update, and never changes anyone else’s gradient. The same points give the same (a,b)(\mathbf{a}, b) to the last bit.

The first version of this solver did not have that property. It picked the second index at random, the way the widely copied “simplified SMO” does, and with a loose stopping rule it would sometimes settle on a different point of a degenerate dual, so the number under the button flickered in the third decimal place instead of sitting at zero. Chasing that down is what sent me back to svm.cpp.

Soft margins: C is the price of a violation

Everything above assumes the classes can be separated. Real data usually cannot, and even when it can, insisting on it is a good way to let one mislabelled point ruin the boundary. The fix is to allow violations and charge for them: slack variables ζn\zeta_n let points sit inside the margin or on the wrong side, and the dual comes out with the identical objective, subject only to

0anC,n=1Nantn=00 \le a_n \le C, \qquad \sum_{n=1}^{N} a_n t_n = 0

A single box constraint. CC is the price of one unit of violation: large CC means the solver would rather have a thin margin than let anything inside it, small CC means it will happily give ground.

The interesting failure is at the small end, and the report’s results table contains a lovely example of it. Here is what my solver does on the full four-feature versicolor vs. virginica problem, min-max normalised, linear kernel:

Csupport vectorsat the bound α = Ctraining accuracyrange of y(x) over the data
0.03100 / 10010093.0%−0.212 to +0.212
1.0060 / 1005793.0%−1.949 to +1.802
64.0017 / 1001394.0%−5.368 to +4.011

Read the first row twice. Every point is a support vector, and every one of them is pinned at the upper bound, which is the exact opposite of the property we came here for. The last column says why. The margin lines are where y(x)=±1y(\mathbf{x}) = \pm 1, and the decision function does not reach ±0.22\pm 0.22 anywhere in the training set, so every single point is inside the margin. The margin has swallowed the whole dataset. That is what “a violation is nearly free” buys you: with all 100 multipliers at CC, w=nantnxn\mathbf{w} = \sum_n a_n t_n \mathbf{x}_n has collapsed to CC times the difference of the two class means, and the data has stopped choosing the direction of the boundary at all. What is left is a nearest-class-mean classifier wearing an SVM’s clothes. On this pair that still scores 93% (the class means are far enough apart) but nothing about the answer came from the margin, and on a problem where the class means are less informative it falls apart completely, which is what the report’s numbers show.

The four kernels

The choice on offer is whatever LIBSVM offers, which the report lists and svm.cpp implements in four lines (SVM/libsvm-3.22/svm.cpp:231-247):

double kernel_linear(int i, int j) const
{ return dot(x[i],x[j]); }
double kernel_poly(int i, int j) const
{ return powi(gamma*dot(x[i],x[j])+coef0,degree); }
double kernel_rbf(int i, int j) const
{ return exp(-gamma*(x_square[i]+x_square[j]-2*dot(x[i],x[j]))); }
double kernel_sigmoid(int i, int j) const
{ return tanh(gamma*dot(x[i],x[j])+coef0); }

so xx\mathbf{x}'\cdot\mathbf{x}, (γxx+C0)d(\gamma\,\mathbf{x}'\cdot\mathbf{x} + C_0)^d, exp(γxx2)\exp(-\gamma\lVert\mathbf{x}' - \mathbf{x}\rVert^2) and tanh(γxx+C0)\tanh(\gamma\,\mathbf{x}'\cdot\mathbf{x} + C_0). The widget’s kernel selector is those four and its sliders are log2\log_2 of the same γ\gamma, C0C_0 and CC the grid search stepped through, so a setting you find by hand is a LIBSVM command line: the export button writes it out.

The sigmoid is the odd one out and it shows. For most (γ,C0)(\gamma, C_0) it is not a valid kernel at all: the Gram matrix is not positive semi-definite, the dual is not concave, and what the solver returns is whatever the quadratic programme happened to reach. It is the kernel that needs the largest CC before it does anything sensible on Iris, and the only one that is still at chance at C=0.25C = 0.25 on the hard pair.

My solver is not LIBSVM, and here is how it differs

The original experiment called LIBSVM. The widget cannot, so I wrote the solver: a cut-down Sequential Minimal Optimisation, about 130 lines, in src/widgets/support-vector-machines/svm.ts. SMO’s idea is that with the equality constraint nantn=0\sum_n a_n t_n = 0 you cannot move one multiplier alone, but you can move two and solve that two-variable subproblem in closed form. Repeat until the worst remaining KKT violation is under a tolerance, the same 10310^{-3} LIBSVM defaults to.

Where mine is a toy compared to the vendored LIBSVM 3.22:

  • Working-set selection. Each iteration I take the maximal violating pair: the ii that maximises yifi-y_i\nabla f_i over the up-set and the jj that minimises it over the low-set. LIBSVM’s Solver::select_working_set (svm.cpp:786) picks the same ii and then chooses jj to minimise the resulting decrease in the objective, second-order information rather than first, and materially fewer iterations for the same answer.
  • Shrinking. LIBSVM periodically removes variables it believes are stuck at a bound from the active set (do_shrinking, svm.cpp:905). Mine scans all NN on every iteration.
  • Caching. LIBSVM keeps an LRU cache of kernel columns (class Cache, svm.cpp:67) because on a real dataset the Gram matrix does not fit in memory. I just compute the whole N×NN \times N matrix, which is fine for the couple of hundred points a browser widget ever holds and hopeless at any real scale.
  • What is LIBSVM’s. The two-variable subproblem, its box clipping, and the threshold are lifted straight from svm.cpp. I wrote Platt’s b1/b2b_1/b_2 rule for the threshold first, and it disagrees badly whenever no support vector is free of its bound, which is precisely what a very small CC produces, and precisely the regime this post is about. calculate_rho (svm.cpp:966) takes the mean gradient over the free variables, or the midpoint of the interval the bounded ones leave when there are none.

The optimum is the same convex problem, so a converged cut-down SMO and a converged LIBSVM agree. Two checks that it does.

A case you can do on paper. Two points, x1=(0,0)\mathbf{x}_1 = (0,0) with t=1t = -1 and x2=(1,0)\mathbf{x}_2 = (1,0) with t=+1t = +1, linear kernel, CC large. The dual reduces to 2aa2/22a - a^2/2, maximised at a=2a = 2, giving w=(2,0)\mathbf{w} = (2, 0), b=1b = -1 and a margin of 1/w=0.51/\lVert\mathbf{w}\rVert = 0.5 either side. My solver returns a=(2,2)a = (2, 2), b=1b = -1 and y(0.5,0)=0y(0.5, 0) = 0, exactly. Set C=0.5C = 0.5 and both multipliers clamp to 0.50.5, as the box constraint requires.

The report’s own experiment. Re-running the protocol below (50 random 60/40 splits, min-max normalisation fitted on the training half, all four features) with the converged parameters from the report’s Table 6:

ProblemKernelCminereport (Table 5)
versicolor vs virginicaLinear1.0094.35%94.50%
versicolor vs virginicaRadial basis1.0094.10%93.40%
versicolor vs virginicaLinear64.0093.05%93.65%
versicolor vs setosaLinear0.0386.75%85.35%
versicolor vs setosaRadial basis0.25100.00%100.00%
virginica vs setosaLinear1.00100.00%100.00%

Tuning it in 2018: coarse to fine, in powers of two

Following Hsu et al., the search is a grid in log2\log_2 of every parameter, repeatedly refined. Twelve rounds; each round evaluates the grid, keeps the best point, and rebuilds a grid around it at half the step (SVM/PbQ1v2.py:59-73):

for stepPow in range(12):
    scores,values = zip(*_pool.map(getscorepartial,product(*params)))
    bvi=np.argmax(scores)
    besterrparams=values[bvi]
    step1 = nextstep1
    nextstep1 = step1/2
    step2 = nextstep2
    nextstep2 = step2/2
    params[1]=[besterrparams[1]]
    params[2]=np.arange(besterrparams[2]-step1,besterrparams[2]+step1+nextstep1,nextstep1)
    params[3]=np.arange(besterrparams[3]-step2,besterrparams[3]+step2+nextstep2,nextstep2)

dd is pinned after the first round because it has to stay an integer; γ\gamma and C0C_0 keep halving. CC is not refined at all: it is swept over thirteen fixed powers of two so its effect can be read separately, which is what makes Table 5 and the results widget possible. A multiprocessing.Pool(40) evaluates the grid, one process per grid cell.

The genuinely interesting design decision is the scoring function, and the report is honest about why it is not the textbook one. vv-fold cross-validation splits the training data into vv parts and holds each out in turn. With 60 training points, five folds leaves validation sets of a dozen, on a problem where most sensible parameter settings get all of them right, so the search is handed a pile of ties at zero error and no way to choose. The replacement is a batch of random resplits, which produces a finer-grained score (the report says 40 of them at 70/30, and the code says otherwise on both counts, as the list below records) plus an early-abandon rule so hopeless settings do not eat the budget (SVM/irusUtil.py:70-82):

def getVfoldCrossValidation(param, XTrain, yTrain, n, splitRatio=0.7, prewnPercentage = 85):
    errorAccumulator = 0
    ri=np.arange(len(yTrain))

    for i in range(1,n+1):
        trainx, trainy, testx, testy = IrusDataReader.ShuffelAndSplit(XTrain,yTrain,ri,splitRatio)
        prob = svmutil.svm_problem(trainy,trainx.tolist())
        m = svmutil.svm_train(prob, param, '-q')
        p_label, p_acc, p_val = svmutil.svm_predict(testy,testx.tolist(), m, '-q')
        errorAccumulator +=p_acc[0]
        if i > 3 and errorAccumulator < prewnPercentage*i:
            return n*errorAccumulator/i
    return errorAccumulator

After three splits, if the running mean accuracy is under 85%, stop and extrapolate. The second motivation given for it is a good one that I had forgotten: badly chosen parameters are also the slowest to fit, because the solver grinds against a problem it cannot satisfy and often hits the iteration cap. Abandoning them early buys more than the fits it skips.

The grid-search panel in the widget above is that procedure in miniature: a 13 × 9 grid of (C,γ)(C, \gamma) for the radial basis kernel over whatever points are on the plot, each cell scored by five random 70/30 resplits, with the 85% rule as a checkbox so you can watch how many fits it saves. It is a coarse grid run once rather than twelve refinements, and five resplits rather than forty, because it has to finish between two animation frames. The one optimisation that makes it possible is worth stating: the kernel matrix does not depend on CC, so the 585 fits share 45 Gram matrices.

The results

The three problems are the three pairs of Iris classes. Setosa is linearly separable from both of the others (you can see the gap in any projection) and versicolor and virginica are not:

Scatter plot of petal length against petal width for Iris versicolor in blue and Iris virginica in red. The two clouds meet and interleave along a diagonal band around petal length 5 centimetres, with roughly half a dozen points of each class inside the other's territory.

Petal length (x3x_3) against petal width (x4x_4) for the hard pair, one of the 18 scatter plots SVM/PlotClassesv2.py wrote, report Fig. 21, p. 22. The overlap in the middle is the whole reason this pair tops out around 94.5% and the other two reach 100%. Load this pair into the playground above (it is a preset, with all six feature pairs selectable) and no value of C gets you a clean split.

The widget below is the report’s own numbers: mean test accuracy over 50 random 60/40 splits, for every one of the thirteen CC values, read out of SVM/OutputStats.npz. Table 5 in the report prints four of those thirteen columns; this is all of them, with the interquartile band across the 50 trials.

InteractiveThe report's Iris results, all thirteen values of C
Mean test accuracy, versicolor vs virginica, from SVM/OutputStats.npz.
CLinearPolynomialRadial basisSigmoid
0.01647.4093.4547.4546.85
0.0347.4093.4047.4046.75
0.0667.8593.3074.1547.00
0.1289.3593.1593.1547.05
0.2593.1093.5093.6047.40
1.0094.5093.1093.4091.55
4.0094.7093.3093.9593.55
64.0093.6592.6592.9593.85

With JavaScript on, this is a chart of all thirteen C values for any of the three class pairs, with quartiles and a readout of the exact numbers.

Table 5 as the report prints it (p. 25), transcribed from SVM/Output/MeanTabel.tex:

ProblemCLinearPolynomialRadial basisSigmoid
Versicolor vs Virginica0.0347.4093.4047.4046.75
0.2593.1093.5093.6047.40
1.0094.5093.1093.4091.55
64.0093.6592.6592.9593.85
Versicolor vs Setosa0.0385.3599.8077.7553.35
0.25100.0099.45100.00100.00
1.00100.0099.50100.00100.00
64.00100.00100.00100.00100.00
Virginica vs Setosa0.0392.9599.3078.4053.05
0.25100.0098.65100.0099.65
1.00100.0099.05100.00100.00
64.00100.00100.00100.00100.00

And Table 6 (p. 26), the parameters that got the lowest validation error, from SVM/Output/ParamsTabel.tex:

ProblemKerneldγC₀C
Versicolor vs VirginicaLinear8.00
Polynomial2.001.1840.910.06
Radial basis0.624.00
Sigmoid0.240.1364.00
Versicolor vs SetosaLinear0.12
Polynomial2.000.250.024.00
Radial basis0.250.25
Sigmoid0.090.091.00
Virginica vs SetosaLinear0.12
Polynomial2.000.280.051.00
Radial basis0.250.25
Sigmoid0.090.091.00

Three things come out of those two tables.

The separable pairs are boring, and that is the finding. Once C0.25C \ge 0.25, every kernel scores 100% on both setosa problems. The gap between setosa and everything else is wide enough that no point ever needs to fall inside the margin, so the extra capacity of a polynomial or radial basis kernel buys nothing, and the converged parameters in Table 6 agree, settling on small γ\gamma and small C0C_0, which is to say on the simplest function in the family.

The hard pair rewards a middle CC. Versicolor vs virginica peaks at 94.50% with a linear kernel at C=1C = 1 in Table 5, and the full sweep (which the report does not print but the .npz does) peaks a little higher still, at 94.70% at C=4C = 4, before drifting back down to 93.65% at C=64C = 64. Too small and the boundary dissolves; too large and it contorts itself around individual overlapping points and generalises worse. The report says this in one sentence in §4.4 and it is the correct sentence.

Table 6’s CC values are not to be trusted, and the report says so. The selection rule takes the first CC reaching the highest validation score, and with only five resplits of a 60-point training set on those two problems there are ties everywhere. That is how the linear kernel ends up recorded at C=0.12C = 0.12 on both setosa pairs when the full sweep is already at 100% mean test accuracy by C=0.06C = 0.06, and on the versicolor–setosa pair C=0.12C = 0.12 is, on my own solver, the last value at which not one support vector is free of its bound. The report flags this itself in §4.4: “these C values are relatively small and may be sub-optimal”. More data, or more resplits, would break the ties.

What I would change

The sparsity argument holds up, and I still think it is the most interesting thing about the model, more interesting than the accuracies, which on Iris are a foregone conclusion. What I would change is around it.

The grid search is the weak part. Twelve rounds of coarse-to-fine, on three problems, four kernels and thirteen values of CC, across 50 trials, with a 40-process pool, is a lot of compute spent to distinguish 93.4% from 93.6% on 100 flowers, and a good chunk of it, as noted above, was spent refining parameters the linear kernel ignores. The honest version of this experiment is smaller: one coarse grid, and the compute saved spent on more trials so the error bars mean something.

I would also not evaluate a sigmoid kernel without saying what it is doing. It is in the table because LIBSVM offers it, not because there was a reason to think tanh\tanh of an inner product was a sensible similarity for flower measurements, and its results are the least interpretable in the set.

The thread this sits on is worth naming, though. The previous post built a Gaussian process, which is the other kernel method in the same assignment: it also writes predictions as a weighted sum of kernel evaluations against the training data, and it also never builds ϕ\boldsymbol\phi explicitly. The difference is exactly the one this post is about. The GP keeps every training point in its predictive mean and pays O(N3)O(N^3) for the privilege; the SVM’s optimality conditions hand back a subset and let you throw the rest away. Kernel methods trade dimensions for data points, and the SVM is the one that then gives most of the data points back.