07 - Support Vector Machines: Kernel Trick and Non-Linear SVM

You have two crescent-shaped clouds of points, one nested inside the curve of the other, and a linear classifier separates them with a straight line that, on the crescents of the example further down, gets 16 points out of 100 wrong (test accuracy 0.840). A Support Vector Machine (SVM) solves the problem in two moves: among the infinitely many lines that separate the classes, it picks the one farthest from the points of both; then, with the kernel trick, it computes that line in a space where the crescents become separable, without ever building that space.

The price is two numbers to choose, C and gamma, and a training time that grows at least quadratically with the number of samples. This article is the seventh in the ML Fundamentals series; you need logistic regression and cross-validation, and after trees and forests it is another non-linear model you meet in the series.

Every formula you read below was checked with a number: the margin of five points computed by hand, the equality between the primal and dual objectives, the polynomial kernel compared with the explicit features, the effect of C and gamma, the measured training times. Where the experiment contradicted a sentence that circulates in textbooks or in earlier versions of this text, the experiment wins and I say so.

The margin, computed by hand on five points

A hyperplane in two dimensions is a line, described by a vector $w$ (the direction perpendicular to the line) and a number $b$: the points with $w \cdot x + b = 0$ lie on the line, those with a positive sign on one side, those with a negative sign on the other. With the labels written as $y_i \in \{-1, +1\}$, a point is classified correctly when $y_i (w \cdot x_i + b) > 0$.

The margin is the width of the empty strip around the line: if you require $y_i (w \cdot x_i + b) \geq 1$ for every point, the strip has width $2 / \|w\|$, so maximizing the margin is equivalent to minimizing $\tfrac{1}{2}\|w\|^2$. The points that touch the edge of the strip are the support vectors.

Take two negative points, $(0,0)$ and $(0,1)$, and three positive ones, $(2,0)$, $(2,1)$ and $(5,3)$. By eye the best line is $x_1 = 1$, that is $w = (1, 0)$ and $b = -1$: the point $(2,0)$ gives $1 \cdot 2 - 1 = 1$, the point $(0,0)$ gives $-(0 - 1) = 1$, and the margin is $2 / 1 = 2$.

The listing asks libsvm (the library underneath SVC) for the same answer. I set C=1e6: a huge C makes the error on the margin so expensive that the model behaves like a hard-margin SVM, which suits separable data like this.

# Five points in the plane: the margin is computed by hand, then compared with libsvm
import numpy as np
from sklearn.svm import SVC

X = np.array([[0, 0], [0, 1], [2, 0], [2, 1], [5, 3]], dtype=float)
y = np.array([-1, -1, 1, 1, 1])

clf = SVC(kernel="linear", C=1e6).fit(X, y)   # huge C: almost hard margin
w = clf.coef_[0]
b = clf.intercept_[0]
print("w =", np.round(w, 6), " b =", round(b, 6))
print("margin 2/||w|| =", round(2 / np.linalg.norm(w), 6))
print("support vector indices:", clf.support_.tolist())
print("alpha * y of the support vectors:", np.round(clf.dual_coef_[0], 6).tolist())

# w rebuilt from the dual: sum of alpha_i y_i x_i
w_dual = clf.dual_coef_[0] @ clf.support_vectors_
print("w from the dual =", np.round(w_dual, 6))
# dual constraint: sum of alpha_i y_i = 0
print("sum of alpha_i y_i =", round(float(clf.dual_coef_[0].sum()), 9))
# y_i f(x_i) for every point: 1 on the support vectors, larger elsewhere
print("y f(x) =", np.round(y * clf.decision_function(X), 4).tolist())

# Output:
# w = [1. 0.]  b = -1.0
# margin 2/||w|| = 2.0
# support vector indices: [1, 3]
# alpha * y of the support vectors: [-0.5, 0.5]
# w from the dual = [1. 0.]
# sum of alpha_i y_i = 0.0
# y f(x) = [1.0, 1.0, 1.0, 1.0, 4.0]

The model finds exactly $w = (1, 0)$, $b = -1$ and margin 2, like the hand computation. The y f(x) line shows that four points out of five are worth 1 (they sit on the edge), while $(5,3)$ is worth 4 and is far away. There is a surprise, though: only two support vectors are returned, indices 1 and 3, that is $(0,1)$ and $(2,1)$.

The other two points on the edge have coefficient zero. This happens because the dual solution is not unique when several points sit on the margin: two coefficients $\alpha = 0.5$ are enough to produce $w = 0.5\,(2,1) - 0.5\,(0,1) = (1, 0)$. So in scikit-learn «support vector» means «point with a non-zero coefficient», which is a subset of the points that touch the margin.

Soft margin, C and the dual problem

With real data the classes overlap and no line separates everything. The soft margin lets every point violate the margin by an amount $\xi_i \geq 0$ (the slack) and makes it pay for that: you minimize $\tfrac{1}{2}\|w\|^2 + C \sum_i \xi_i$.

The optimal slack of a point is its hinge loss, $\max(0,\, 1 - y_i (w \cdot x_i + b))$: zero if the point is outside the margin, growing linearly if it violates it. C is the price of the violation: a large C punishes the error and narrows the margin, a small C accepts more violations in exchange for a wide margin.

The same problem can be written in the Lagrange multipliers $\alpha_i$ (the dual): you maximize $\sum_i \alpha_i - \tfrac{1}{2}\sum_{i,j} \alpha_i \alpha_j y_i y_j\, x_i \cdot x_j$ with $0 \leq \alpha_i \leq C$ and $\sum_i \alpha_i y_i = 0$. The data appear only inside the dot products $x_i \cdot x_j$, and that is the entry point of the kernel.

Because the problem is convex, at the optimum the two objectives have the same value. The listing checks it on two noisy moons, computing the primal objective from the weights and the slacks and the dual one from the coefficients returned by dual_coef_, which holds the products $\alpha_i y_i$.

# Strong duality: on the soft margin the primal and dual objectives coincide at the optimum
import numpy as np
from sklearn.datasets import make_moons
from sklearn.svm import SVC
from sklearn.metrics.pairwise import linear_kernel

X, y01 = make_moons(n_samples=200, noise=0.35, random_state=0)
y = 2 * y01 - 1          # labels in {-1, +1}
C = 1.0
clf = SVC(kernel="linear", C=C, tol=1e-8).fit(X, y)

w, b = clf.coef_[0], clf.intercept_[0]
xi = np.maximum(0.0, 1.0 - y * (X @ w + b))               # slack = hinge loss
primal = 0.5 * w @ w + C * xi.sum()

a = clf.dual_coef_[0]                                     # alpha_i * y_i
Ksv = linear_kernel(clf.support_vectors_)
dual = np.abs(a).sum() - 0.5 * a @ Ksv @ a

print(f"support vectors: {len(clf.support_)} out of {len(X)}")
print(f"alpha = C (violate the margin): {int((np.abs(a) > C - 1e-6).sum())}")
print(f"0 < alpha < C (on the margin):  {int((np.abs(a) < C - 1e-6).sum())}")
print(f"primal objective = {primal:.6f}")
print(f"dual objective   = {dual:.6f}")
print(f"difference       = {abs(primal - dual):.2e}")

# Output:
# support vectors: 86 out of 200
# alpha = C (violate the margin): 83
# 0 < alpha < C (on the margin):  3
# primal objective = 82.499680
# dual objective   = 82.499680
# difference       = 1.16e-07

The two values (82.499680) differ by 1e-07, that is by the solver tolerance: scikit-learn's dual and the primal you wrote by hand describe the same optimum. Counting the coefficients shows the three regimes of the dual:

• $\alpha_i = 0$: the point is outside the margin and does not enter the model (114 points out of 200);
• $0 < \alpha_i < C$: the point sits exactly on the edge of the margin (3 points);
• $\alpha_i = C$: the point violates the margin or is misclassified (83 points).

With high noise almost all the support vectors are of the third type: the solution is determined by the points inside the margin or on the wrong side, not by those on the edge. For prediction the 86 support vectors are enough, $f(x) = \sum_i \alpha_i y_i\, x_i \cdot x + b$; the other 114 points can be thrown away and the model does not change.

The kernel trick, checked with numbers

The crescents cannot be separated by a line, but they can if you add to the coordinates $(x_1, x_2)$ their squares and their product. The kernel is a function $K(x, z)$ that directly returns the dot product between the «enriched» versions $\varphi(x) \cdot \varphi(z)$ without computing $\varphi$.

Since the dual uses the data only through $x_i \cdot x_j$, it is enough to replace it with $K(x_i, x_j)$ and the SVM works in the enriched space. The smallest case you can check: in two dimensions $(x \cdot z)^2 = \varphi(x) \cdot \varphi(z)$ with $\varphi(x) = (x_1^2,\ \sqrt{2}\,x_1 x_2,\ x_2^2)$.

The values in parentheses are the SVC defaults in the scikit-learn 1.9 documentation. The default kernel is RBF; gamma='scale' uses $1 / (n_{feature} \cdot \mathrm{Var}(X))$, where $\mathrm{Var}(X)$ is the variance of all the values of the matrix taken together.

The listing below checks four things: the identity of the polynomial kernel with the explicit features, the hand-written RBF formula, the positive semidefiniteness of the RBF Gram matrix (the condition for a kernel to correspond to a real dot product, by Mercer's theorem) and a counterexample for the sigmoid.

# The kernel trick checked with numbers
import numpy as np
from sklearn.datasets import make_moons
from sklearn.metrics.pairwise import rbf_kernel, polynomial_kernel

rng = np.random.default_rng(0)
X = rng.normal(size=(6, 2))

# 1) Homogeneous degree-2 polynomial kernel in 2 dimensions:
#    (x . z)^2 = phi(x) . phi(z)  with  phi(x) = (x1^2, sqrt(2) x1 x2, x2^2)
def phi(x):
    return np.array([x[0] ** 2, np.sqrt(2) * x[0] * x[1], x[1] ** 2])

K_kernel = (X @ X.T) ** 2
F = np.array([phi(x) for x in X])
K_explicit = F @ F.T
print("poly degree 2: max |K - phi.phi| =", f"{np.abs(K_kernel - K_explicit).max():.2e}")
K_sk = polynomial_kernel(X, degree=2, gamma=1.0, coef0=0.0)
print("matches sklearn polynomial_kernel:", np.allclose(K_kernel, K_sk))

# 2) RBF: hand-written formula against sklearn
gamma = 0.7
D2 = ((X[:, None, :] - X[None, :, :]) ** 2).sum(-1)
K_rbf = np.exp(-gamma * D2)
print("hand-written RBF = sklearn rbf_kernel:", np.allclose(K_rbf, rbf_kernel(X, gamma=gamma)))

# 3) Mercer in practice: the RBF Gram matrix on 300 points is positive semidefinite
P, _ = make_moons(n_samples=300, noise=0.3, random_state=0)
G = rbf_kernel(P, gamma=1.0)
lam = np.linalg.eigvalsh(G)
print(f"smallest eigenvalue of the RBF Gram matrix: {lam.min():.2e}  (>= 0 up to numerical error)")

# 4) The sigmoid kernel does not guarantee semidefiniteness: a numerical counterexample
S = np.tanh(2.0 * (X @ X.T) - 1.0)
print(f"smallest eigenvalue of the sigmoid Gram matrix (gamma=2, coef0=-1): {np.linalg.eigvalsh(S).min():.3f}")

# Output:
# poly degree 2: max |K - phi.phi| = 8.88e-16
# matches sklearn polynomial_kernel: True
# hand-written RBF = sklearn rbf_kernel: True
# smallest eigenvalue of the RBF Gram matrix: -1.17e-14  (>= 0 up to numerical error)
# smallest eigenvalue of the sigmoid Gram matrix (gamma=2, coef0=-1): -3.434

The polynomial kernel and the dot product of the explicit features coincide to within 1e-15, that is to rounding error. For the RBF you cannot write $\varphi$ (the space has infinite dimension), but the Gram matrix on 300 points has a smallest eigenvalue of the order of -1e-14: no negative eigenvalue beyond numerical noise, as the theory guarantees for every $\gamma > 0$.

The sigmoid, on the other hand, with these parameters and these six points, produces an eigenvalue of -3.43: its matrix is not positive semidefinite and does not correspond to any dot product. It is not an implementation defect, it is the reason why in practice you start from linear and RBF.

What C and gamma do to the boundary

With the RBF kernel gamma sets the scale of the similarity: $K(x, z)$ is close to 1 for nearby points and close to 0 for distant ones, and a large gamma narrows the «nearby». A small gamma gives smooth, almost linear boundaries; a large gamma gives a boundary that wraps around the individual points.

To see it I trained the SVM on the same set of 200 training points (400 samples, noise 0.30, half and half) with several pairs of values and measured the accuracy on the 200 points held out. The listing also prints the number of support vectors, which is the quickest signal of how the model is doing.

# Effect of C and gamma on the boundary, on two noisy moons (synthetic data, fixed seed)
import numpy as np
from sklearn.datasets import make_moons
from sklearn.model_selection import train_test_split
from sklearn.svm import SVC

X, y = make_moons(n_samples=400, noise=0.30, random_state=1)
Xtr, Xte, ytr, yte = train_test_split(X, y, test_size=0.5, random_state=1, stratify=y)

print("kernel   C       gamma   acc_train acc_test  n_SV  (out of", len(Xtr), "train)")
rows = [("linear", 1, None), ("rbf", 1, 0.01), ("rbf", 1, 1.0), ("rbf", 1, 100.0),
        ("rbf", 0.01, 1.0), ("rbf", 1000, 1.0), ("rbf", 1000, 100.0)]
for kern, C, g in rows:
    kw = {} if g is None else {"gamma": g}
    m = SVC(kernel=kern, C=C, **kw).fit(Xtr, ytr)
    gs = "-" if g is None else g
    print(f"{kern:7s} {C:7g} {str(gs):7s} {m.score(Xtr, ytr):9.3f} {m.score(Xte, yte):8.3f} {m.n_support_.sum():5d}")

# Output:
# kernel   C       gamma   acc_train acc_test  n_SV  (out of 200 train)
# linear        1 -           0.860    0.840    72
# rbf           1 0.01        0.830    0.825   144
# rbf           1 1.0         0.905    0.895    68
# rbf           1 100.0       0.950    0.830   185
# rbf        0.01 1.0         0.845    0.830   200
# rbf        1000 1.0         0.915    0.900    50
# rbf        1000 100.0       1.000    0.805   170

Read it row by row, remembering that with 200 test points one point is worth half a percentage point and differences of 1-2 points are noise. With gamma=0.01 the boundary is too rigid (0.830 on training, 0.825 on test, 144 support vectors): it is underfitting, at the level of the line (0.840, difference within the noise). With gamma=1 and C=1 you get 0.905 and 0.895 with 68 support vectors, a good balance.

With gamma=100 the training score rises to 0.950 but the test score falls to 0.830: overfitting. The last row takes it to the extreme, C=1000 and gamma=100: perfect training, 1.000, and test at 0.805, the worst value in the table. The number of support vectors is high in both wrong directions: 185 with a too high gamma, 200 out of 200 with C=0.01 (no point manages to stay outside the very wide margin).

That is why a high fraction of support vectors is a symptom, not a diagnosis: it tells you that something is wrong, not whether it is underfitting or overfitting.

Scaling is not optional

The Euclidean distance sums the squares of the differences in every feature, so the feature measured with the largest numbers rules. The breast cancer dataset included in scikit-learn (load_breast_cancer, 569 samples and 30 features) is a textbook case: the variance of the 30 features goes from 7e-06 to 3e+05, more than ten orders of magnitude.

The comparison below uses the same RBF SVM with the same 5 folds, with and without StandardScaler. The scaler sits inside the Pipeline, and that is not a detail: if you applied it to the whole dataset before the cross-validation, the mean and standard deviation of the test folds would enter the computation: a small leak of information from test to training.

# Without scaling, the Euclidean distance is decided by the feature with the largest unit
import numpy as np
from sklearn.datasets import load_breast_cancer
from sklearn.model_selection import StratifiedKFold, cross_val_score
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import StandardScaler
from sklearn.svm import SVC

X, y = load_breast_cancer(return_X_y=True)
cv = StratifiedKFold(n_splits=5, shuffle=True, random_state=0)

print("variance per feature: min = %.2e, max = %.2e" % (X.var(axis=0).min(), X.var(axis=0).max()))
without = cross_val_score(SVC(kernel="rbf"), X, y, cv=cv)
with_sc = cross_val_score(make_pipeline(StandardScaler(), SVC(kernel="rbf")), X, y, cv=cv)
print(f"RBF without scaling:  accuracy {without.mean():.4f} (std {without.std():.4f})")
print(f"RBF with StandardScaler in the pipeline: accuracy {with_sc.mean():.4f} (std {with_sc.std():.4f})")

# gamma='scale' = 1 / (n_features * X.var()): recomputed by hand
m = SVC(kernel="rbf").fit(X, y)
print("gamma used by sklearn:", m._gamma, "  formula:", 1.0 / (X.shape[1] * X.var()))

# Output:
# variance per feature: min = 6.99e-06, max = 3.24e+05
# RBF without scaling:  accuracy 0.9210 (std 0.0308)
# RBF with StandardScaler in the pipeline: accuracy 0.9771 (std 0.0070)
# gamma used by sklearn: 6.395533747973492e-07   formula: 6.395533747973492e-07

Without scaling the accuracy is 0.921 with a standard deviation of 0.031 across folds; with scaling 0.977 with a standard deviation of 0.007. The gain of more than five points does not come from a better model, but from the fact that the model without scaling looks almost only at the few features with huge values.

The last line confirms the documentation's formula: gamma='scale' is $1 / (30 \cdot \mathrm{Var}(X))$, about 6.4e-07 on the raw data. Note also that _gamma is a private attribute, used here only for the check. The practical guide by the libsvm authors recommends scaling every attribute to [0, 1] or [-1, +1], with the same coefficients on training and test.

Choosing C and gamma with a grid

No value of C and gamma is right in absolute terms: they depend on the scale of the data and on the noise. The libsvm authors recommend a grid search with cross-validation over values that grow exponentially, because the effect of these parameters is multiplicative.

The listing splits the dataset into 398 training rows and 171 test rows, searches the grid with 5 stratified folds on the training set only, and then measures once on the test set. That test set has never been seen by the search: it is the only honest number, whereas the cross-validation score is optimistic because it is the maximum over sixteen combinations.

# Choosing C and gamma with a logarithmic grid; scaling stays INSIDE the validation
import numpy as np
from sklearn.datasets import load_breast_cancer
from sklearn.model_selection import GridSearchCV, StratifiedKFold, train_test_split
from sklearn.pipeline import Pipeline
from sklearn.preprocessing import StandardScaler
from sklearn.svm import SVC

X, y = load_breast_cancer(return_X_y=True)
Xtr, Xte, ytr, yte = train_test_split(X, y, test_size=0.3, random_state=0, stratify=y)

pipe = Pipeline([("scale", StandardScaler()), ("svc", SVC(kernel="rbf"))])
grid = {"svc__C": [0.1, 1, 10, 100], "svc__gamma": [0.0003, 0.003, 0.03, 0.3]}
cv = StratifiedKFold(n_splits=5, shuffle=True, random_state=0)
gs = GridSearchCV(pipe, grid, cv=cv, scoring="accuracy").fit(Xtr, ytr)

print("best parameters:", gs.best_params_)
print(f"CV accuracy (398 training rows): {gs.best_score_:.4f}")
print(f"accuracy on the held-out test set (171 rows): {gs.score(Xte, yte):.4f}")
svc = gs.best_estimator_.named_steps["svc"]
print(f"support vectors: {svc.n_support_.sum()} out of {len(Xtr)} = {svc.n_support_.sum() / len(Xtr):.1%}")

# CV table: mean for each pair (rows = C, columns = gamma)
means = gs.cv_results_["mean_test_score"].reshape(4, 4)
print("\nC \\ gamma   " + "  ".join(f"{g:>7g}" for g in grid["svc__gamma"]))
for C, row in zip(grid["svc__C"], means):
    print(f"{C:>9g}   " + "  ".join(f"{v:7.3f}" for v in row))

# Output:
# best parameters: {'svc__C': 100, 'svc__gamma': 0.003}
# CV accuracy (398 training rows): 0.9874
# accuracy on the held-out test set (171 rows): 0.9591
# support vectors: 35 out of 398 = 8.8%
#
# C \ gamma    0.0003    0.003     0.03      0.3
#       0.1     0.628    0.897    0.950    0.628
#         1     0.900    0.965    0.985    0.937
#        10     0.965    0.980    0.980    0.937
#       100     0.985    0.987    0.977    0.937

The search picks C=100 and gamma=0.003 with 0.9874 in cross-validation and 0.9591 on the test set: the difference of almost three points depends partly on having chosen the maximum among many trials and partly on having only 171 test rows (one error more or less moves the value by 0.6 points). The table shows a diagonal ridge: high values of C go with low gamma and vice versa (0.985 with C=1 and gamma=0.03, 0.985 with C=100 and gamma=0.0003).

The two knobs partly compensate for the same effect, so several different pairs give an almost equivalent model. The chosen model uses 35 support vectors out of 398, 8.8%. At the edges of the table the values fall to 0.628, which is the frequency of the most numerous class in the training set: a model that always predicts the same class would do the same.

More than two classes

The SVM was born for two classes. With $K$ classes SVC trains $K(K-1)/2$ classifiers, one for every pair, and the prediction is a majority vote (the one-vs-one strategy). The scikit-learn documentation says that the solver always uses this strategy and that decision_function_shape changes only the shape of the output of decision_function: 'ovo' returns one score per pair, 'ovr' (the default) turns it into one score per class.

The one-vs-rest strategy is something else: $K$ models, each trained on all the samples, which you get with OneVsRestClassifier. The listing uses the digits of load_digits (1,797 8x8 images, 10 classes).

# SVC is always one-vs-one; decision_function_shape changes only the SHAPE of the output
import numpy as np
from sklearn.datasets import load_digits
from sklearn.model_selection import train_test_split
from sklearn.multiclass import OneVsRestClassifier
from sklearn.svm import SVC

X, y = load_digits(return_X_y=True)
Xtr, Xte, ytr, yte = train_test_split(X / 16.0, y, test_size=0.3, random_state=0, stratify=y)

ovo = SVC(kernel="rbf", gamma=0.05, C=10, decision_function_shape="ovo").fit(Xtr, ytr)
ovr_shape = SVC(kernel="rbf", gamma=0.05, C=10, decision_function_shape="ovr").fit(Xtr, ytr)
print("classes K =", len(ovo.classes_), " pairs K(K-1)/2 =", len(ovo.classes_) * (len(ovo.classes_) - 1) // 2)
print("decision_function (ovo):", ovo.decision_function(Xte).shape)
print("decision_function (ovr):", ovr_shape.decision_function(Xte).shape)
print("same predictions:", bool((ovo.predict(Xte) == ovr_shape.predict(Xte)).all()))
print(f"SVC accuracy (one-vs-one): {ovo.score(Xte, yte):.4f}")

# Explicit one-vs-rest: K models, each trained on ALL the samples
ovr = OneVsRestClassifier(SVC(kernel="rbf", gamma=0.05, C=10)).fit(Xtr, ytr)
print("models trained by OneVsRestClassifier:", len(ovr.estimators_))
print(f"one-vs-rest accuracy: {ovr.score(Xte, yte):.4f}")

# Output:
# classes K = 10  pairs K(K-1)/2 = 45
# decision_function (ovo): (540, 45)
# decision_function (ovr): (540, 10)
# same predictions: True
# SVC accuracy (one-vs-one): 0.9926
# models trained by OneVsRestClassifier: 10
# one-vs-rest accuracy: 0.9907

With 10 classes the pairwise classifiers are 45, and decision_function has 45 columns with 'ovo' and 10 with 'ovr'. The predictions are identical in both cases: the parameter does not change the model. The accuracy is 0.9926 for one-vs-one and 0.9907 for one-vs-rest, a difference of one sample out of 540: neither strategy wins here.

There is a difference in cost, though: each pairwise classifier sees only the samples of two classes, while each of the ten one-vs-rest models sees all the samples, so with many classes the pairwise strategy builds many small models.

When the SVM does not scale

The SVC documentation says that the training time grows at least quadratically with the number of samples and that it can become impractical beyond a few tens of thousands, and recommends LinearSVC, SGDClassifier or a Nystroem transformation (which approximates the RBF kernel with a fixed number of explicit features) before a linear model.

The listing measures both routes on a synthetic problem with 20 features, with training samples from 3,000 to 48,000 and 4,000 test points. These are timings: they change from one run to the next, because they depend on the machine and on what else is running, and that is why I ran the listing five times and report ranges. The label of the last printed line, esponente stimato, means «estimated exponent»: the listing is the same one I ran, so I left it as it is.

# Training time: RBF SVC against Nystroem + LinearSVC (TIMES: they change from one run to the next)
import time
import numpy as np
from sklearn.datasets import make_classification
from sklearn.kernel_approximation import Nystroem
from sklearn.model_selection import train_test_split
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import StandardScaler
from sklearn.svm import SVC, LinearSVC

X, y = make_classification(n_samples=52000, n_features=20, n_informative=8,
                           flip_y=0.05, class_sep=1.0, random_state=0)
Xtr, Xte, ytr, yte = train_test_split(X, y, test_size=4000, random_state=0, stratify=y)
sc = StandardScaler().fit(Xtr)
Xtr, Xte = sc.transform(Xtr), sc.transform(Xte)

print(f"{'n train':>8} {'SVC t(s)':>9} {'SVC acc':>8} {'Nys+Lin t(s)':>13} {'acc':>6}")
tempi = []
for n in [3000, 6000, 12000, 24000, 48000]:
    t0 = time.perf_counter()
    m = SVC(kernel="rbf", gamma="scale").fit(Xtr[:n], ytr[:n])
    t_svc = time.perf_counter() - t0
    t0 = time.perf_counter()
    a = make_pipeline(Nystroem(gamma=1 / (20 * 1.0), n_components=300, random_state=0),
                      LinearSVC(dual="auto", max_iter=5000)).fit(Xtr[:n], ytr[:n])
    t_nys = time.perf_counter() - t0
    tempi.append(t_svc)
    print(f"{n:8d} {t_svc:9.2f} {m.score(Xte, yte):8.3f} {t_nys:13.2f} {a.score(Xte, yte):6.3f}")
esp = np.polyfit(np.log([3000, 6000, 12000, 24000, 48000]), np.log(tempi), 1)[0]
print(f"esponente stimato t ~ n^p per SVC: p = {esp:.2f}")

# Output:
#  n train  SVC t(s)  SVC acc  Nys+Lin t(s)    acc
#     3000      0.10    0.947          0.14  0.934
#     6000      0.34    0.951          0.31  0.941
#    12000      1.33    0.955          0.62  0.943
#    24000      4.52    0.959          1.29  0.943
#    48000     24.53    0.963          2.71  0.947
# esponente stimato t ~ n^p per SVC: p = 1.95

The accuracy is deterministic, the times are not. On an AMD Ryzen 5 5600H with 2 threads, in five runs (with other processes active on the machine), SVC took 0.08-0.17 s at 3,000 samples and 25-46 s at 48,000: sixteen times the samples, from about 240 to 580 times the time, with an estimated exponent between 1.95 and 2.22, that is roughly quadratic. The documentation says the time grows at least quadratically and that the complexity of the QP solver lies between quadratic and cubic; an estimate from five points has an uncertainty of that order.

Nystroem with 300 components followed by LinearSVC is slower or equal at 3,000-6,000 samples (0.10-0.34 s), then wins: 2.6-3.3 s at 48,000, from nine to sixteen times faster, but losing about 1.6 points of accuracy (0.947 against 0.963). So it is not true that the approximation «always speeds things up by 10-50 times»: below tens of thousands of samples the exact SVM is competitive, and the advantage appears where the exact SVM becomes slow. This synthetic problem is not a general benchmark, it serves to show the shape of the curve.

Probabilities from an SVM: Platt scaling

An SVM produces a signed score, the distance from the hyperplane in the kernel space, not a probability. Platt scaling trains a one-variable logistic regression on the score, using out-of-sample scores, and turns it into a probability. In SVC this was done by the probability=True parameter, which in scikit-learn 1.9 is deprecated and will be removed in 1.11: the documentation points to CalibratedClassifierCV(SVC(), ensemble=False).

The same documentation warns that probability=True slows down fit because it internally uses a 5-fold cross-validation and that predict_proba can be inconsistent with predict. The listing checks the warning, calibrates with the recommended method and compares the Brier score with an unfitted sigmoid.

# Probabilities from an SVM: Platt scaling with CalibratedClassifierCV (SVC(probability=True) is deprecated in 1.9)
import warnings
import numpy as np
import sklearn
from sklearn.calibration import CalibratedClassifierCV
from sklearn.datasets import load_breast_cancer
from sklearn.metrics import brier_score_loss
from sklearn.model_selection import train_test_split
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import StandardScaler
from sklearn.svm import SVC

print("scikit-learn", sklearn.__version__)
X, y = load_breast_cancer(return_X_y=True)
Xtr, Xte, ytr, yte = train_test_split(X, y, test_size=0.3, random_state=0, stratify=y)
base = make_pipeline(StandardScaler(), SVC(kernel="rbf", C=10, gamma=0.003))

# 1) the old way: what the installed version says
with warnings.catch_warnings(record=True) as w:
    warnings.simplefilter("always")
    make_pipeline(StandardScaler(), SVC(probability=True)).fit(Xtr, ytr)
print("warning:", str(w[0].message)[:110] if w else "none")

# 2) the way the documentation indicates: Platt scaling (sigmoid) on out-of-sample scores
cal = CalibratedClassifierCV(base, method="sigmoid", cv=5, ensemble=False).fit(Xtr, ytr)
p = cal.predict_proba(Xte)[:, 1]
print(f"Brier score on the test set (lower = better): {brier_score_loss(yte, p):.4f}")

# 3) the naive alternative: sigmoid of the raw score, not fitted to the data
score = base.fit(Xtr, ytr).decision_function(Xte)
p_naive = 1 / (1 + np.exp(-score))
print(f"Brier with 1/(1+exp(-score)) not fitted:     {brier_score_loss(yte, p_naive):.4f}")
print(f"agreement between predict and the 0.5 threshold on predict_proba: {((p > 0.5) == base.predict(Xte)).mean():.3f}")

# Output:
# scikit-learn 1.9.0
# warning: The `probability` parameter was deprecated in 1.9 and will be removed in version 1.11. Use `CalibratedClassifi
# Brier score on the test set (lower = better): 0.0315
# Brier with 1/(1+exp(-score)) not fitted:     0.0421
# agreement between predict and the 0.5 threshold on predict_proba: 0.994

The Brier score is the mean squared error between the predicted probability and the 0/1 outcome, and lower is better. The calibration fitted to the data gives 0.0315, the naive sigmoid $1/(1+e^{-s})$ applied to the raw score gives 0.0421: the shape of the transformation must be learned, not assumed.

Note also that the 0.5 threshold on the calibrated probability and predict agree in 99.4% of cases, not 100%: one sample out of 171 changes class, and that is exactly the inconsistency the documentation points out. If you only need a confidence score, and not a probability, use decision_function and skip the calibration.

Worked exercise

Compute by hand the RBF kernel matrix with $\gamma = 0.5$ on the points $a=(0,0)$, $b=(1,1)$ and $c=(3,0)$, and say which pair is the most similar according to the kernel. The procedure: compute the squared distances, multiply them by $-\gamma$, exponentiate. Then check with the listing.

# Exercise: RBF kernel matrix with gamma = 0.5 on three points
import numpy as np

P = np.array([[0.0, 0.0], [1.0, 1.0], [3.0, 0.0]])
gamma = 0.5
d2 = ((P[:, None, :] - P[None, :, :]) ** 2).sum(-1)    # squared Euclidean distances
K = np.exp(-gamma * d2)
print("distances^2:\n", d2)
print("K:\n", np.round(K, 4))

# Output:
# distances^2:
#  [[0. 2. 9.]
#  [2. 0. 5.]
#  [9. 5. 0.]]
# K:
#  [[1.     0.3679 0.0111]
#  [0.3679 1.     0.0821]
#  [0.0111 0.0821 1.    ]]

What you can do now

After this article you can do five things, each with a numerical check at hand:

• compute by hand the margin of a set of points and verify it with SVC;
• explain why the primal and dual objectives coincide and what it means for a point to be a support vector;
• check that a kernel really is a dot product;
• choose C and gamma with a logarithmic grid inside a Pipeline with the scaler;
• recognize underfitting and overfitting from the difference between training and test, and know when to give up the exact SVM.

The five abilities hold together: the margin and the dual explain what the model does, the kernel explains how it becomes non-linear, and the last two items tell you how to use it without being fooled by the numbers of a single run or by a grid that picks the maximum among too many trials.

We did not cover support vector regression, kernels on strings or graphs, and the derivation of SMO: the scikit-learn guide is the starting point, and the Machine Learning Engineer series has an article on the maximum margin and the kernel trick. In the next article we move from supervised to unsupervised learning with K-Means and hierarchical clustering; the evaluation of models, with the metrics that here we used only as accuracy, is in the article on metrics. If the features are many, PCA is a way to reduce them before a kernel model, and as for any other step the effect must be measured; taking an SVM to production, with the scaling checks that have already appeared here, is the subject of deployment.