33 - Support Vector Machine: the maximum margin and the kernel trick

Two neighbourhoods are separated by an empty plot of land, and you have to draw a road that divides them. You can run it very close to the houses of the first neighbourhood or to those of the second, or down the middle, as far as possible from both. If someone later builds a house slightly out of place, the central road is the one that still separates the two neighbourhoods well. A Support Vector Machine (SVM) does exactly this: among all the boundaries that separate the classes, it picks the one with the widest margin.

Then it adds two powerful ideas: accepting a few errors when the classes are not separable, and the kernel trick, which allows curved boundaries without ever computing coordinates in gigantic spaces. In this lesson we derive the mathematics, implement a linear SVM by hand, and see where the library behaves treacherously.

Three mathematical tools, explained before we use them

• Dot product w · x = w₁x₁ + w₂x₂ + ... + w_nx_n. It is a number that says "how much x points in the direction of w". If w is a unit vector, w · x is the length of the shadow of x along w. With fixed weights w, it is also a score: the weighted sum of the features.
• Norm ‖w‖ = √(w₁² + ... + w_n²): the length of the vector. Dividing a vector by its norm makes it a unit vector.
• Convex problem: an optimisation problem whose function to minimise is bowl-shaped, with no secondary dips. It has a single bottom and any descent finds it.

The three tools are needed at three different moments. The dot product tells you which side of the boundary a point is on, the norm turns that score into a real distance, and convexity guarantees that training has a single right answer. If you have to remember one, remember the norm: the whole story of the margin is a division by ‖w‖. A kernel, which appears in the title, is instead a function that measures the "similarity" between two points and lets the SVM draw curved boundaries: we define it further down, once the margin is clear.

The hyperplane and the distance from the boundary

In two dimensions a linear boundary is a line; in n dimensions it is a hyperplane, the set of points x with

w · x + b = 0

where w is the weight vector (perpendicular to the hyperplane) and b the intercept. A point lies on one side or the other depending on the sign of w · x + b. The distance of a point x from the hyperplane is

d(x) = |w · x + b| / ‖w‖

The factor ‖w‖ appears because multiplying w and b by a constant does not change the hyperplane, but it would change the value of w · x + b: dividing by ‖w‖ makes the distance independent of this arbitrary choice.

In the drawing the black line is the boundary; the two dashed lines are the edges of the margin; the circled points are the support vectors, the ones that "touch" the edges. If you remove a far-away point, the boundary does not move: it depends only on the points closest to it.

The margin and the optimisation problem

Let us call the two classes y = +1 and y = −1. A correct classifier satisfies y_i (w · x_i + b) > 0 for every point. Since w and b can be scaled freely, we fix the scale so that the points closest to the boundary have exactly y_i (w · x_i + b) = 1. These points are called support vectors. Their distance from the hyperplane is 1/‖w‖ on each side, so the total margin is

margin = 2 / ‖w‖

Maximising the margin is equivalent to minimising ‖w‖, or (for convenience in the derivation) ½‖w‖². The hard margin problem is:

minimise    ½ ‖w‖²
subject to  y_i (w · x_i + b) ≥ 1    for every i

It is a convex problem: it has a single minimum and no misleading local minima. But it works only if the classes are perfectly separable by a hyperplane, and a single misplaced point (noise, a labelling error) is enough to make it impossible.

Soft margin: the parameter C

Some points are allowed to violate the margin, at the price of a penalty. For each point we introduce a slack variable ξ_i ≥ 0, the measure of the violation:

minimise    ½ ‖w‖² + C · Σ_i ξ_i
subject to  y_i (w · x_i + b) ≥ 1 − ξ_i ,   ξ_i ≥ 0

For a given (w, b), the smallest admissible value of ξ_i is max(0, 1 − y_i (w · x_i + b)). Substituting it, the problem becomes unconstrained and reads as "regularisation + error":

minimise    ½ ‖w‖² + C · Σ_i max(0, 1 − y_i (w · x_i + b))
                ↑ wide margin        ↑ hinge loss: errors and violations

The function max(0, 1 − m) with m = y(w·x+b) is called the hinge loss: it costs zero if the point is correctly classified and outside the margin (m ≥ 1), and it grows linearly as the point enters the margin or lands on the wrong side. The parameter C sets the trade-off:

• Large C: every violation costs a lot, the model tries to classify even the noisy points correctly, narrow margin, risk of overfitting.
• Small C: violations cost little, wide margin, the model tolerates errors, risk of underfitting.

The two items are not two models but the two ends of the same knob, and neither value is right in itself: it depends on how much noise there is in the labels. C rewards accuracy on the training set, while 1/C rewards the width of the margin. Further down, a measured table shows what happens to the support vectors and to the margin when C goes from 0.001 to 1000, on data that really do overlap.

A linear SVM by hand

The unconstrained form can be minimised with gradient descent (you will find it in full in Gradient descent; here the idea is enough: you take a step in the direction that reduces the error). The hinge loss is not differentiable at the junction point, but it has a simple subgradient: where the point violates the margin (y(w·x+b) < 1) the contribution to the gradient is −y·x, elsewhere it is zero. Dividing the objective by C·n we get (λ/2)‖w‖² + the mean of the hinge, with λ = 1/(C·n):

The term λ·w comes from the regularisation, and at every step it acts as a small force pulling w back towards zero. The second term appears only for the "violating" points, that is those inside the margin or on the wrong side, and pushes w towards their label. Written out in full:

∂/∂w = λ·w − (1/n) · Σ_{violating i} y_i x_i
∂/∂b =     − (1/n) · Σ_{violating i} y_i

The listing translates these two lines into NumPy: viola is the mask of the points that violate the margin, gw and gb are the two derivatives written above, and the step lr decreases as 1/√t. At the end it compares the result with scikit-learn's SVC(kernel="linear"), which solves the same problem by another route.

import numpy as np
from sklearn.datasets import make_blobs
from sklearn.svm import SVC

X, y01 = make_blobs(n_samples=200, centers=[(-1.5, -1), (1.5, 1)],
                    cluster_std=1.1, random_state=1)
y = 2 * y01 - 1                       # labels in {-1, +1}

def svm_sgd(X, y, C=1.0, epochs=2000, lr0=0.5):
    n, d = X.shape
    lam = 1.0 / (C * n)
    w, b = np.zeros(d), 0.0
    for t in range(1, epochs + 1):
        viola = y * (X @ w + b) < 1                    # points inside the margin or misclassified
        gw = lam * w - (y[viola, None] * X[viola]).sum(axis=0) / n
        gb = -y[viola].sum() / n
        lr = lr0 / np.sqrt(t)                           # decreasing step
        w, b = w - lr * gw, b - lr * gb
    return w, b

w, b = svm_sgd(X, y)
sk = SVC(kernel="linear", C=1.0).fit(X, y)
w_sk = sk.coef_[0]

cos = w @ w_sk / (np.linalg.norm(w) * np.linalg.norm(w_sk))
print("hand-made / sklearn accuracy: %.3f / %.3f" % ((np.sign(X @ w + b) == y).mean(), sk.score(X, y)))
print("hand-made w:", np.round(w, 3), " sklearn w:", np.round(w_sk, 3))
print("cosine between the two directions: %.4f" % cos)
print("hand-made margin %.3f   sklearn %.3f" % (2 / np.linalg.norm(w), 2 / np.linalg.norm(w_sk)))
print("support vectors per class (sklearn):", sk.n_support_, "out of", len(X), "points")

hand-made / sklearn accuracy: 0.950 / 0.950
hand-made w: [1.31  0.889]  sklearn w: [1.415 0.914]
cosine between the two directions: 0.9997
hand-made margin 1.263   sklearn 1.187
support vectors per class (sklearn): [14 15] out of 200 points

What the code does. make_blobs generates two Gaussian clouds of 100 points each around the centres (−1.5, −1) and (1.5, 1), with standard deviation 1.1: they overlap a little, so the problem is not perfectly separable and the soft margin is needed. y = 2*y01 - 1 converts the labels from {0,1} to {−1,+1}, as the hinge loss requires. In svm_sgd the loop repeats 2000 times: it computes the functional margin of every point y * (X @ w + b), finds the points that violate it (value below 1) in the boolean vector viola, computes the subgradient only on those and takes a step of size lr, decreasing as 1/√t (a shrinking step makes the search settle instead of oscillating). In the comparison, SVC(kernel="linear") exposes coef_ (the vector w) and intercept_ (b).

The direction of w coincides almost perfectly (cosine 0.9997) and the accuracy is identical; the norm is slightly different because gradient descent with a decreasing step is an approximation, whereas the library solves the problem in dual form with much greater precision. The margin is 1.26 against 1.19. Notice also that the support vectors are 29 out of 200 points: only they determine the boundary. You can remove any other point and the boundary does not move. This is the property that gives the algorithm its name.

The kernel trick

If the data are not separable by a hyperplane in the original space (two interleaved moons, concentric circles), the idea is to transform the features with a map φ into a larger space where they become separable, and to use a hyperplane there. A concrete example: in 2D you cannot separate with a line the points inside a circle from those outside; but if you add the coordinate x₁² + x₂² the circle becomes a threshold on a single coordinate.

The cost of this idea looks prohibitive: φ can send the data into thousands (or infinitely many) dimensions. The trick lies in the dual form of the problem, in which the data appear only through dot products between pairs of points:

f(x) = Σ_i α_i y_i ⟨φ(x_i), φ(x)⟩ + b      (α_i > 0 only for the support vectors)

If there is a function K(x, z) = ⟨φ(x), φ(z)⟩ that can be computed directly from the original points, φ never has to be built. An example you can check by hand: in two dimensions, the polynomial kernel K(x, z) = (x · z)² corresponds to φ(x) = (x₁², √2·x₁x₂, x₂²). Indeed:

⟨φ(x), φ(z)⟩ = x₁²z₁² + 2·x₁x₂z₁z₂ + x₂²z₂²
             = (x₁z₁ + x₂z₂)²  =  (x · z)²

Let us check the identity numerically on two hand-picked points, x = (1, 2) and z = (3, −1): on one side the dot product of the images phi(x) and phi(z) in three dimensions, on the other the kernel (x·z)² computed directly in the two original dimensions.

phi = lambda v: np.array([v[0] ** 2, np.sqrt(2) * v[0] * v[1], v[1] ** 2])
x, z = np.array([1.0, 2.0]), np.array([3.0, -1.0])
print("kernel (x.z)^2 computed in the original space:", (x @ z) ** 2)
print("dot product after the phi map, in 3D:         %.6f" % (phi(x) @ phi(z)))

kernel (x.z)^2 computed in the original space: 1.0
dot product after the phi map, in 3D:         1.000000

The code checks the identity numerically: phi turns a 2D point into three numbers (x₁², √2·x₁x₂, x₂²). The dot product of the two transformed points is 1.000000, and the kernel (x·z)² computed without transforming anything is 1.0 (here x·z = 3 − 2 = 1). This is the heart of the trick: the same result, with a computation in 2 dimensions instead of 3. For kernels such as the RBF the map φ has infinitely many dimensions and direct computation is the only practicable route.

The most widely used kernel is the RBF (Gaussian): K(x, z) = exp(−γ‖x − z‖²). It corresponds to a φ with infinitely many dimensions, and is interpreted as a "similarity": it equals 1 for coincident points and tends to 0 as the distance grows. The parameter γ sets the radius of influence of each point: large γ means local influence (a very jagged boundary), small γ wide influence (a smooth, almost linear boundary).

Comparing kernels on non-linear data

from sklearn.datasets import make_moons
from sklearn.model_selection import train_test_split
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import StandardScaler

Xm, ym = make_moons(n_samples=600, noise=0.25, random_state=0)
Xt, Xe, yt, ye = train_test_split(Xm, ym, test_size=0.3, stratify=ym, random_state=0)

for kernel in ("linear", "poly", "rbf"):
    m = make_pipeline(StandardScaler(), SVC(kernel=kernel, C=1.0)).fit(Xt, yt)
    print("%-6s train %.3f  test %.3f  support vectors: %d" % (
        kernel, m.score(Xt, yt), m.score(Xe, ye), m[-1].n_support_.sum()))

linear train 0.886  test 0.856  support vectors: 126
poly   train 0.876  test 0.839  support vectors: 135
rbf    train 0.936  test 0.917  support vectors: 103

What the code does. make_moons produces two interleaved half-circles with noise 0.25 (600 points, 420 of them for training). The Pipeline chains StandardScaler (which brings every feature to mean 0 and standard deviation 1, using only the training set) and SVC; the attribute m[-1] takes the last step of the pipeline in order to read n_support_, the number of support vectors per class.

The RBF kernel beats the linear one on the test set (0.917 against 0.856): the curved boundary follows the shape of the two moons. The polynomial kernel (degree 3 by default) does not help here: a more expressive kernel is not automatically better. Note also that the number of support vectors is high (103 out of 420 training points with the RBF, about one in four): with overlapping classes many points end up inside the margin.

Overfitting with C and gamma: the training accuracy that lies

With large C and large γ the SVM can "encircle" every single training point. Look at what happens when we fix C = 1000 and vary γ:

for gamma in (0.1, 1, 10, 100, 1000):
    m = make_pipeline(StandardScaler(), SVC(kernel="rbf", C=1000, gamma=gamma)).fit(Xt, yt)
    print("gamma %6.1f  train %.3f  test %.3f  support vectors: %d" % (
        gamma, m.score(Xt, yt), m.score(Xe, ye), m[-1].n_support_.sum()))

gamma    0.1  train 0.957  test 0.917  support vectors: 66
gamma    1.0  train 0.964  test 0.911  support vectors: 50
gamma   10.0  train 1.000  test 0.867  support vectors: 109
gamma  100.0  train 1.000  test 0.856  support vectors: 378
gamma 1000.0  train 1.000  test 0.667  support vectors: 419

What the code does: it keeps C fixed at 1000 and varies γ over five orders of magnitude, printing training accuracy, test accuracy and the number of support vectors. The test set here serves only to show the phenomenon: in a real project you would choose γ with cross-validation on the training set.

Training accuracy reaches 1.000 already at γ = 10, while test accuracy drops from 0.917 (γ = 0.1) to 0.667 (γ = 1000): the model has memorised the noise. The number of support vectors is a useful alarm bell: with γ = 1000 there are 419 out of 420, that is almost every point is a support vector and the model is memorising. The two hyperparameters must be chosen with cross-validation (GridSearchCV), searching on a logarithmic scale (0.001, 0.01, 0.1, 1, 10...), never with the test set.

The feature-scale trap

The SVM with an RBF kernel is based on Euclidean distances: a feature with values in the thousands dominates all the others. To show it, we multiply by 1000 one of the informative features of a synthetic dataset:

from sklearn.datasets import make_classification

Xc, yc = make_classification(n_samples=600, n_features=5, n_informative=3,
                             n_redundant=0, random_state=3)
Xc[:, 0] *= 1000                                   # one feature with a huge scale
Xa, Xb, ya, yb = train_test_split(Xc, yc, test_size=0.3, stratify=yc, random_state=0)

print("without scaling: test %.3f" % SVC().fit(Xa, ya).score(Xb, yb))
print("with scaling:    test %.3f" % make_pipeline(StandardScaler(), SVC()).fit(Xa, ya).score(Xb, yb))

without scaling: test 0.817
with scaling:    test 0.950

What the code does: make_classification with 3 informative features out of 5; then Xc[:, 0] *= 1000 multiplies the first column by a thousand, simulating for example an amount in euros next to percentages. The two evaluation lines train the same SVC with and without the StandardScaler.

Same model, same data: standardisation alone moves the accuracy from 0.817 to 0.950. The scaler goes inside the Pipeline (otherwise you would compute the mean and standard deviation on the whole dataset and cause data leakage, as the process of an ML project explains).

Trade-offs

• Pros: effective in high-dimensional spaces, excellent with little data, a convex problem (no local minima), kernels allow non-linear boundaries, small memory footprint (only the support vectors are kept).
• Cons: training scales badly with the number of samples (the scikit-learn documentation gives a scaling between n_features·n_samples² and n_features·n_samples³ for the libsvm-based implementation): with hundreds of thousands of rows it becomes slow; it requires scaling; delicate hyperparameters; it does not natively produce probabilities (probability=True adds a calibration with five-fold cross-validation, which the documentation describes as expensive).
• Practical choice: small or medium datasets with many features → RBF SVM is an excellent candidate; huge datasets → LinearSVC or SGDClassifier, or tree-based models.

Pros and cons are two faces of the same structure: the solution depends only on the support vectors, so it is compact and robust, but to find them the solver has to compare pairs of samples, and the cost grows quickly. The practical rule is therefore one of scale. Up to a few thousand or tens of thousands of rows the RBF is a good candidate; beyond that, the choice is linear or tree-based. If you need probabilities, factor in the cost of calibration.

Real traps

• Forgetting the StandardScaler (see above).
• Choosing C and γ by looking at the test set: use GridSearchCV on the training set.
• Labels in {0,1} in the hand-made implementation: the hinge loss requires {−1,+1}. Measured by replacing y with y01 in the call to svm_sgd: no error, w = (0.131, 0.009), and the accuracy against the true labels drops to 0.500, that is chance.
• Imbalanced classes: use class_weight="balanced" and evaluate with precision/recall.
• Expecting calibrated probabilities from decision_function: it is a signed distance, not a probability.

The five traps do not weigh the same. The first two silently distort the measurement: the model looks good because it was evaluated badly. The third produces an absurd model without any error. The last two concern how to read the result, not how to obtain it. If you have to check one thing before handing over a model, check that the scaler is inside the pipeline.

The parameter C up close: margin, support vectors and accuracy

In the derivation you saw that C sets the trade-off between a wide margin and few errors. We measure it on a problem with overlapping classes (two close Gaussian clouds, hence not separable): for each C we train a linear SVM and read the number of support vectors, the width of the margin 2/‖w‖ computed from coef_, and the accuracy on training and test sets (200 points each, so a single point is worth 0.005). At the end we print the three loss functions, 0-1 error, hinge and logistic, computed for seven values of the functional margin m.

from sklearn.datasets import make_blobs, load_digits, make_classification
from sklearn.model_selection import train_test_split

X_ov, y_ov = make_blobs(n_samples=400, centers=[(-1, -0.7), (1, 0.7)], cluster_std=1.3, random_state=2)
X_ov_tr, X_ov_te, y_ov_tr, y_ov_te = train_test_split(X_ov, y_ov, test_size=0.5, stratify=y_ov, random_state=0)

print("C         support vectors   margin width   train acc.   test acc.")
for C in (0.001, 0.01, 0.1, 1, 10, 1000):
    m = SVC(kernel="linear", C=C).fit(X_ov_tr, y_ov_tr)
    margin = 2 / np.linalg.norm(m.coef_[0])
    print("%-8g  %12d      %10.3f      %9.3f    %9.3f" % (C, m.n_support_.sum(), margin, m.score(X_ov_tr, y_ov_tr), m.score(X_ov_te, y_ov_te)))

# the three loss functions, as a function of the functional margin m = y (w.x + b)
margins = np.array([-2, -1, -0.5, 0, 0.5, 1, 2])
print("\nm              ", margins)
print("0-1 error      ", (margins <= 0).astype(int))
print("hinge (SVM)    ", np.round(np.maximum(0, 1 - margins), 3))
print("logistic      ", np.round(np.log(1 + np.exp(-margins)), 3))

C         support vectors   margin width   train acc.   test acc.
0.001              200           7.963          0.805        0.780
0.01               123           3.573          0.805        0.795
0.1                 90           2.315          0.805        0.805
1                   82           2.040          0.825        0.800
10                  82           2.041          0.825        0.800
1000                82           2.040          0.825        0.800

m               [-2.  -1.  -0.5  0.   0.5  1.   2. ]
0-1 error       [1 1 1 1 0 0 0]
hinge (SVM)     [3.  2.  1.5 1.  0.5 0.  0. ]
logistic       [2.127 1.313 0.974 0.693 0.474 0.313 0.127]

The table tells of three phases. With C = 0.001 the penalty for violations is so low that all 200 training points become support vectors and the margin widens to 7.96. With C = 0.1 the vectors drop to 90, the margin to 2.32 and the test accuracy reaches the maximum of the table, 0.805. From C = 1 upwards the solution no longer changes (82 vectors, margin 2.04): in this experiment beyond that threshold the term C·Σξ dominates and the violations are already at the lowest possible level. Mind the scale of the differences, though: between 0.780 and 0.805 on the test set there are five points out of 200, so I would not read "C = 0.1 is the best" as a discovery, but as an indication that very small values underfit.

The rows below show why the SVM is called that. The hinge equals zero for m ≥ 1 (the point is outside the margin and produces no gradient), whereas the logistic loss, which is the one of logistic regression, never reaches zero: even a distant, well-classified point keeps contributing a little. It is the same property of the support vectors seen from the loss side.

More than two classes: one-vs-one, one-vs-rest

The SVM was born for two classes. For K classes several binary classifiers are combined. In the one-vs-rest strategy (OvR) K of them are trained, each "class k against the rest", and the class with the highest score is chosen. In the one-vs-one strategy (OvO) K(K − 1)/2 are trained, one for every pair of classes, and the decision is made by majority vote. With the 10 digits of Digits: 10 classifiers against 45. We measure accuracy and training time, adding as a reference a linear SVM (LinearSVC), which uses no kernel.

The listing does two things: it compares the three multiclass strategies on the digits and then measures how much the fit time grows when the samples double. The time is read with time.time() around the fit alone; what matters is the ratio with the previous row (t_rbf / previous), not the absolute value.

import time
from sklearn.svm import LinearSVC
from sklearn.multiclass import OneVsRestClassifier

X_d, y_d = load_digits(return_X_y=True)
X_d_tr, X_d_te, y_d_tr, y_d_te = train_test_split(X_d, y_d, test_size=0.3, stratify=y_d, random_state=0)

models = {"SVC RBF, one-vs-one (default)": SVC(kernel="rbf"),
          "SVC RBF, one-vs-rest": OneVsRestClassifier(SVC(kernel="rbf")),
          "LinearSVC (linear, one-vs-rest)": LinearSVC(max_iter=5000)}
for name, clf in models.items():
    pipe = make_pipeline(StandardScaler(), clf)
    t0 = time.time(); pipe.fit(X_d_tr, y_d_tr); elapsed = time.time() - t0
    print("%-52s accuracy %.3f   training time %.2f s" % (name, pipe.score(X_d_te, y_d_te), elapsed))

print("\ntraining time as the number of samples grows (20 features):")
print("n          SVC (RBF)     LinearSVC")
previous = None
for n in (1000, 2000, 4000, 8000):
    X_n, y_n = make_classification(n_samples=n, n_features=20, n_informative=8, random_state=0)
    X_n = StandardScaler().fit_transform(X_n)
    t0 = time.time(); SVC().fit(X_n, y_n); t_rbf = time.time() - t0
    t0 = time.time(); LinearSVC(max_iter=5000).fit(X_n, y_n); t_lin = time.time() - t0
    print("%-9d  %8.3f s    %8.3f s    %s" % (n, t_rbf, t_lin, "" if previous is None else "(x%.1f versus n/2)" % (t_rbf / previous)))
    previous = t_rbf

# Output (one run; the times change at every run)
# SVC RBF, one-vs-one (default)                        accuracy 0.981   training time 0.08 s
# SVC RBF, one-vs-rest                                 accuracy 0.980   training time 0.22 s
# LinearSVC (linear, one-vs-rest)                      accuracy 0.948   training time 0.45 s
#
# training time as the number of samples grows (20 features):
# n          SVC (RBF)     LinearSVC
# 1000          0.022 s       0.003 s    
# 2000          0.079 s       0.004 s    (x3.6 versus n/2)
# 4000          0.246 s       0.010 s    (x3.1 versus n/2)
# 8000          0.501 s       0.016 s    (x2.0 versus n/2)

The first comparison trains three variants on the digits (10 classes, 1,797 images of 8×8), always with standardisation inside the Pipeline. The accuracy is 0.98 for the SVM with the RBF kernel, both with the one-vs-one strategy (the default choice: the SVC documentation indicates it as the multiclass scheme) and with OneVsRestClassifier; LinearSVC reaches 0.95. The kernel is worth it: there is a three-point gap with respect to the linear model. On timing, in two consecutive runs one-vs-one took 0.07-0.09 s against 0.16-0.21 s for one-vs-rest (the run shown above gave 0.08 s and 0.22 s): the order was the same both times, but these are tenths of a second on a single machine and must not drive the choice. The plausible explanation is that the 45 one-vs-one classifiers each see only two classes, hence little data, and the cost of a kernel SVM grows more than linearly.

The second experiment is the one that matters. Going from 1,000 to 8,000 samples (eight times the data) the kernel SVM took, over six separate runs, from 17 to 30 times as long, that is an exponent between 1.4 and 1.6; LinearSVC from 5 to 9 times. Author's measurements, 2026-11-12, 12 logical threads, Python 3.14.4, scikit-learn 1.8.0: they vary from one run to the next, and at these values of n they do not yet show the asymptotic law. The SVC documentation says that the training time grows at least quadratically with the samples and that it can be impracticable beyond a few tens of thousands, and recommends LinearSVC or SGDClassifier, possibly after a kernel approximation. Do not extrapolate from 8,000 to 800,000 with my exponent: trust the documented limit.

Inside the dual solution: the decision as a sum over the support vectors

In the derivation we stated that the decision function is written f(x) = Σ_i α_i y_i K(x_i, x) + b, with α_i different from zero only for the support vectors. Let us check it by rebuilding by hand the predictions of an SVC from its internal attributes: support_vectors_ (the points x_i), dual_coef_ (the products α_i y_i) and intercept_ (b). For the RBF kernel we compute K(x_i, x) = exp(−γ‖x_i − x‖²), with γ = 0.5 fixed in the listing (gamma=0.5, in place of the default 'scale') so that the hand formula uses the same number as the library.

X_m, y_m = make_blobs(n_samples=300, centers=[(-1, -0.7), (1, 0.7)], cluster_std=1.1, random_state=5)
rbf = SVC(kernel="rbf", C=1.0, gamma=0.5).fit(X_m, y_m)

SV = rbf.support_vectors_                                         # (n_sv, 2): the support vectors x_i
alpha_y = rbf.dual_coef_[0]                                       # (n_sv,): the products alpha_i * y_i

def decision_by_hand(x_new):
    d2 = ((SV[None, :, :] - x_new[:, None, :]) ** 2).sum(axis=2)  # (n_new, n_sv): squared distances
    return (np.exp(-0.5 * d2) @ alpha_y) + rbf.intercept_[0]      # sum_i alpha_i y_i K(x_i, x) + b

x_test = np.array([[0.0, 0.0], [2.0, 1.5], [-2.0, -1.0], [0.5, -1.0]])
print("support vectors: %d out of %d points (%.0f%%)" % (len(SV), len(X_m), 100 * len(SV) / len(X_m)))
print("decision by hand:      ", np.round(decision_by_hand(x_test), 4))
print("decision_function:     ", np.round(rbf.decision_function(x_test), 4))
print("max |difference| over the whole training set: %.2e" % np.abs(decision_by_hand(X_m) - rbf.decision_function(X_m)).max())

# in the linear case the sum can be compressed into a single vector w = sum_i alpha_i y_i x_i
lin = SVC(kernel="linear", C=1.0).fit(X_m, y_m)
w_from_dual = lin.dual_coef_[0] @ lin.support_vectors_
print("w rebuilt from the dual form: %s   scikit-learn coef_: %s" % (np.round(w_from_dual, 4), np.round(lin.coef_[0], 4)))

support vectors: 108 out of 300 points (36%)
decision by hand:       [-0.1481  1.0233 -1.0086 -0.049 ]
decision_function:      [-0.1481  1.0233 -1.0086 -0.049 ]
max |difference| over the whole training set: 0.00e+00
w rebuilt from the dual form: [1.405  0.7192]   scikit-learn coef_: [1.405  0.7192]

The listing has three parts. It trains an SVM with an RBF kernel and extracts the three ingredients of the dual formula. decision_by_hand computes, for each new point, the squared distances from all the support vectors using broadcasting (shape (n_new, n_sv)), turns them into RBF similarities with exp(−γ d²), weights them by α_i y_i and sums, adding the bias. The values coincide with scikit-learn's decision_function over the whole training set, with a maximum difference of zero in this run (at most of the order of 10⁻¹⁵): the dual formula is not an abstraction, it is exactly what the library computes. This has a concrete consequence: the cost of prediction is proportional to the number of support vectors, because every new point has to be compared with each of them; with a linear kernel, instead, the sum compresses into a single vector w = Σ α_i y_i x_i (the last line shows it: it coincides with coef_) and prediction costs a single multiplication.

It follows that a model with too many support vectors is slow and bulky in production, and that a small C, which produces many of them, is not "free" in terms of inference cost.

How to choose the kernel and the hyperparameters: a procedure

1. Standardise (inside the pipeline). Without it, it is the most frequent cause of poor results.
2. Start from the linear one (LinearSVC or SVC(kernel="linear")): fast, and it tells you at once whether the data are nearly separable with a line.
3. Move to the RBF if the linear one is insufficient. Search C and γ on a logarithmic grid (for example C in 10⁻²...10³ and γ in 10⁻⁴...10¹) with cross-validation on the training set only.
4. Check the number of support vectors: if it is almost the whole training set, the model is memorising (γ too large) or the problem is very noisy.
5. Compare with simpler models and with several seeds: a logistic regression or a forest often does as well with less tuning.

The order of the five steps is not random. Standardisation comes first because it changes the distances the RBF works on; the linear model comes before the RBF because it is the term of comparison, and if it is enough you stop. The grid comes only afterwards, because searching C and γ on a model that does not need them wastes time. Two additions: look at the plateau around the maximum of the grid, because an isolated peak is a fragile choice; and remember that prediction costs as much as the number of support vectors. If you need probabilities, the documentation points to probability=True (Platt scaling, with internal cross-validation) or CalibratedClassifierCV, to be validated separately.

Three experiments to try

Exercise 1. Repeat the experiment on C with separable data (cluster_std=0.6): does the margin keep varying for large C? Exercise 2. Compute by hand the number of binary classifiers of OvO for 100 classes: how much does it grow compared with OvR? Exercise 3. Use SVC(kernel="rbf", probability=True) and measure how much slower training becomes: it is the cost of the internal calibration of the probabilities.

If you want to practise on the same topic with a different domain, telling extra virgin and lampante olive oil apart from two chemical indices, there is an exercise skeleton in Python with three parts to complete (the hyperplane from the support vectors, the kernel, the scale). It uses only NumPy and, as long as the points to complete are empty, it prints which exercises are missing.

What you can do now

You can derive the margin 2/‖w‖ and explain what changes with C, you can check by hand that the kernel is a dot product in a space you do not build, you can recognise from the number of support vectors a model that is memorising, and you can choose between linear SVM, RBF and alternatives according to the number of rows. In the next lesson, KNN, the boundary disappears altogether: you classify by looking at the nearest neighbours.

Sources

Previous lesson in the course: Naive Bayes: Bayes' theorem with a naive assumption. Other lessons in the series that link to this one: Constrained Optimisation: Lagrange Multipliers, Duality and the KKT Conditions.