Support Vector Machines and Kernel Methods

A linear support vector classifier balances a wide margin around its decision boundary against penalties for points on the wrong side of that margin. On separable data, the hard-margin version seeks the widest separating strip. With overlap or noise, a soft margin permits violations. The fitted decision rule can be expressed using weighted training points called support vectors. Kernels then let the same calculation describe nonlinear boundaries.

The margin

Write the score as \(s(x)=w^\top x+b\): \(w\) controls the boundary’s orientation and \(b\) its offset. Use labels \(y_i\in\{-1,+1\}\) in the mathematics; the code passes 0 and 1, which the classifier handles internally. Starting from a separating boundary on separable data, multiply both \(w\) and \(b\) by the same positive constant to make the smallest signed score \(y_i s(x_i)\) equal to 1. The boundary \(s(x)=0\) is unchanged. The planes \(s(x)=+1\) and \(s(x)=-1\) are each \(1/\|w\|\) from it, so their separation is \(2/\|w\|\). For example, \(w=(2,0)\), \(b=0\) gives planes at \(x_1=\pm0.5\), a strip of width 1. Maximizing that width under \(y_i s(x_i)\geq1\) is equivalent to minimizing \(\tfrac12\|w\|^2\).

The soft-margin constraints are \(y_i s(x_i)\geq1-\xi_i\), with \(\xi_i\geq0\). Minimize \(\tfrac12\|w\|^2+C\sum_i\xi_i\), where \(C>0\) controls the penalty for violations. At the optimum, \(\xi_i=\max(0,1-y_i s(x_i))\): this is hinge loss. A signed score of 1.4 has zero loss, 0.4 has loss 0.6 despite a correct class prediction, and −0.2 has loss 1.2 and the wrong sign. Thus margin violations and classification errors are different quantities. The code’s margin is the width between score planes, even when points lie within or beyond them; it is not an empty gap between the classes.

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

X, y = make_classification(n_samples=600, n_features=4, n_informative=4,
                           n_redundant=0, class_sep=2.0, flip_y=0.01,
                           random_state=0)

for C in (0.001, 0.01, 0.1, 1.0, 10.0, 1000.0):
    m = SVC(kernel="linear", C=C).fit(X, y)
    w = np.linalg.norm(m.coef_[0])
    print(f"C={C:8.3f}  |w| {w:7.4f}  margin {2 / w:7.4f}"
          f"  support vectors {len(m.support_):4d}  train acc {m.score(X, y):.4f}")
# C=   0.001  |w|  0.3315  margin  6.0335  support vectors  466  train acc 0.7650
# C=   0.010  |w|  0.6531  margin  3.0624  support vectors  342  train acc 0.7983
# C=   0.100  |w|  0.8401  margin  2.3808  support vectors  301  train acc 0.7950
# C=   1.000  |w|  0.8670  margin  2.3069  support vectors  294  train acc 0.7917
# C=  10.000  |w|  0.8683  margin  2.3033  support vectors  294  train acc 0.7933
# C=1000.000  |w|  0.8698  margin  2.2994  support vectors  294  train acc 0.7933

At \(C=0.001\), the strip is 6.0335 wide and the model uses 466 support vectors. Support vectors can lie on a margin plane, inside the strip, or on the wrong side of the decision boundary; their count is not the number strictly inside the strip. Increasing \(C\) to 1 reduces the strip width to 2.3069 and the support count to 294. The count remains 294 in the next two rows, but the coefficients and accuracy still change. These rows show a relatively stable fit at large \(C\), not that the remaining errors cannot be corrected at any price. Training accuracy itself need not improve monotonically, because the optimized loss is hinge loss rather than the error count.

For a fixed dataset and kernel, larger \(C\) weakens the norm penalty relative to hinge loss. Dividing the objective by \(Cn\) gives average hinge loss plus \(\|w\|^2/(2Cn)\). In the convention “average loss plus \(\lambda\|w\|^2/2\),” this means \(\lambda=1/(Cn)\). The sample count and the loss normalization matter when comparing libraries or training sizes. These examples use SVC(kernel="linear"); LinearSVC has different defaults, including squared hinge loss.

Support vectors and refitting

The dual optimization assigns a multiplier \(0\leq\alpha_i\leq C\) to each observation. At an exact optimum, \(w=\sum_i\alpha_i y_i x_i\) and \(\sum_i\alpha_i y_i=0\). Points with \(\alpha_i>0\) are support vectors. The optimality conditions imply that a point with signed score strictly greater than 1 has \(\alpha_i=0\), while \(0<\alpha_i<C\) places a point on a margin plane. Points with \(\alpha_i=C\) may be on or beyond that margin. A point on the plane need not have a positive multiplier in a degenerate solution.

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

X, y = make_classification(n_samples=600, n_features=4, n_informative=4,
                           n_redundant=0, class_sep=2.0, flip_y=0.01,
                           random_state=0)
m = SVC(kernel="linear", C=1.0).fit(X, y)
sv = m.support_
m2 = SVC(kernel="linear", C=1.0).fit(X[sv], y[sv])       # refit on support vectors only

print(f"support vectors {len(sv)} of {len(X)} ({len(sv) / len(X):.1%})")
print(f"all data  coef {np.round(m.coef_[0], 5)}  b {m.intercept_[0]:.5f}")
print(f"SVs only  coef {np.round(m2.coef_[0], 5)}  b {m2.intercept_[0]:.5f}")
print(f"max coefficient difference {np.abs(m.coef_[0] - m2.coef_[0]).max():.2e}")
print(f"predictions differing: {(m.predict(X) != m2.predict(X)).sum()} of {len(X)}")

rng = np.random.default_rng(0)
non_sv = np.setdiff1d(np.arange(len(X)), sv)
keep = np.setdiff1d(np.arange(len(X)),
                    rng.choice(non_sv, size=len(non_sv) // 2, replace=False))
m3 = SVC(kernel="linear", C=1.0).fit(X[keep], y[keep])
print(f"drop half the non-SVs: max coef diff {np.abs(m.coef_[0] - m3.coef_[0]).max():.2e}"
      f"   predictions differing {(m.predict(X) != m3.predict(X)).sum()}")
mt = SVC(kernel="linear", C=1.0, tol=1e-8).fit(X, y)
mt2 = SVC(kernel="linear", C=1.0, tol=1e-8).fit(X[mt.support_], y[mt.support_])
print(f"tight tolerance: max coef diff {np.max(np.abs(mt.coef_ - mt2.coef_)):.2e}"
      f"   predictions differing {(mt.predict(X) != mt2.predict(X)).sum()}")
# support vectors 294 of 600 (49.0%)
# all data  coef [0.28612 0.39012 0.59976 0.39732]  b 1.06210
# SVs only  coef [0.28715 0.39047 0.60088 0.39834]  b 1.06498
# max coefficient difference 1.12e-03
# predictions differing: 1 of 600
# drop half the non-SVs: max coef diff 9.13e-06   predictions differing 0
# tight tolerance: max coef diff 1.10e-08   predictions differing 0

With the default stopping tolerance, removing half the non-support vectors changes coefficients by at most \(9.13\times10^{-6}\) and changes none of the 600 training predictions. Removing all of them gives a coefficient difference of \(1.12\times10^{-3}\) and one changed prediction. Repeating both fits at tol=1e-8 reduces the coefficient difference to about \(1.10\times10^{-8}\) and gives no changed predictions here. These are numerical comparisons, not evidence that deleting an exact zero-multiplier point must move the optimum.

For the exact optimization problem, retaining all positive-multiplier points and keeping the same \(C\), features, and kernel preserves an optimal solution: the remaining dual multipliers still certify it. The norm penalty fixes a unique \(w\), but the intercept can be non-unique in some problems. Numerical solvers can also return slightly different fits. Recomputing a scaler, data-dependent kernel parameters, or class weights after deletion would change the problem. In a kernel model, storing \(s\) dense support vectors takes roughly \(O(sd)\) space, and prediction sums \(s\) kernel contributions. A linear model can instead combine them into \(w\) and predict without storing those rows. The support fraction depends on \(C\), the kernel, and data geometry; it is not a test for separability.

The kernel trick

The decision score can be written as \(s(x)=\sum_{i:\alpha_i>0}\alpha_i y_i K(x_i,x)+b\). For a linear model, \(K(x_i,x)=x_i^\top x\). A valid inner-product kernel satisfies \(K(x,z)=\phi(x)^\top\phi(z)\) for a feature map \(\phi\). Its Gram matrix, whose entry \(K_{ij}\) compares training rows \(i\) and \(j\), must be symmetric and positive semidefinite for every finite input set: \(c^\top Kc\geq0\) for every vector \(c\). An arbitrary similarity function need not have this property. Substituting a valid kernel in the dual fits a linear separator in feature space while evaluating similarities from the original inputs.

import numpy as np
from sklearn.datasets import make_classification

X, _ = make_classification(n_samples=200, n_features=4, n_informative=4,
                           n_redundant=0, random_state=0)

K = (1 + X @ X.T) ** 2                          # the degree-2 polynomial kernel

def phi(x):                                     # the feature map it corresponds to
    d = len(x)
    return np.array([1.0] + [np.sqrt(2) * v for v in x] + [v * v for v in x]
                    + [np.sqrt(2) * x[i] * x[j]
                       for i in range(d) for j in range(i + 1, d)])

Phi = np.array([phi(row) for row in X])
print(f"max |K - Phi Phi^T| {np.abs(K - Phi @ Phi.T).max():.2e}")
print(f"feature-map dimension {Phi.shape[1]}   kernel matrix {K.shape}")
# max |K - Phi Phi^T| 4.55e-13
# feature-map dimension 15   kernel matrix (200, 200)

For two inputs, \((1+x^\top z)^2=1+2\sum_jx_jz_j+\sum_jx_j^2z_j^2+2\sum_{i<j}x_ix_jz_iz_j\). The factors \(\sqrt2\) in phi produce the coefficients 2 when mapped coordinates are multiplied. With four original features, there are 1 constant, 4 linear, 4 square, and 6 cross terms: 15 coordinates. The maximum discrepancy in the printed check is \(4.55\times10^{-13}\), consistent with floating-point arithmetic. The expansion supplies the identity; checking 200 inputs illustrates it.

An inhomogeneous degree-10 polynomial kernel on 1,000 input features has \(\binom{1010}{10}\approx2.91\times10^{23}\) monomial coordinates. Evaluating \((1+x^\top z)^{10}\) needs a length-1,000 dot product and a power, avoiding that explicit expansion. The RBF kernel \(K(x,z)=\exp(-\gamma\|x-z\|^2)\), for \(\gamma>0\), has an infinite-dimensional feature representation on Euclidean space. Finite approximations such as random Fourier features are also possible. Gamma controls the distance scale, while \(C\) controls the norm-versus-loss tradeoff in that kernel space; both influence the fitted boundary.

A fully materialized Gram matrix contains \(n^2\) entries: at 100,000 rows, float64 entries alone occupy about 80 GB. SVC can compute entries as needed and retain a bounded kernel cache, so it does not always allocate the full matrix. Kernel training can nevertheless become expensive as sample count grows, and runtime depends on the data, kernel, tolerance, and cache. There is no universal 100,000-row cutoff. For larger problems, compare linear solvers or finite kernel approximations. The SVM user guide describes the solver and caching considerations.

Exercises

1. Gamma and the RBF distance scale. At fixed \(C=1\), sweep six gamma values from 0.01 to 1000, spanning five orders of magnitude, on a circular boundary. Report training accuracy, cross-validated accuracy, and the support count.

Compare both ends of the sweep. Can a high support fraction appear with weak training performance as well as perfect training accuracy?

import numpy as np
from sklearn.svm import SVC
from sklearn.model_selection import cross_val_score

rng = np.random.default_rng(0)
n = 400
angle = rng.uniform(0, 2 * np.pi, n)
radius = rng.uniform(0, 2, n)
X = np.column_stack([radius * np.cos(angle), radius * np.sin(angle)])
y = (radius > 1).astype(int)

for g in (0.01, 0.1, 1.0, 10.0, 100.0, 1000.0):
    m = SVC(kernel="rbf", gamma=g, C=1.0)
    train = m.fit(X, y).score(X, y)
    print(f"gamma {g:8.2f}  train {train:.4f}"
          f"  cv {cross_val_score(m, X, y, cv=5).mean():.4f}"
          f"  support vectors {len(m.support_):4d}")
# gamma     0.01  train 0.6675  cv 0.6550  support vectors  400
# gamma     0.10  train 0.9700  cv 0.9650  support vectors  180
# gamma     1.00  train 0.9900  cv 0.9825  support vectors   77
# gamma    10.00  train 0.9900  cv 0.9700  support vectors  145
# gamma   100.00  train 1.0000  cv 0.9450  support vectors  346
# gamma  1000.00  train 1.0000  cv 0.7450  support vectors  393

2. Changing feature units. Generate two independent features that enter the label rule symmetrically, rescale only the second, and compare an RBF SVM with and without standardization.

Check whether a unit change alters raw-model accuracy and whether fitting a scaler inside each training fold removes that change.

import numpy as np
from sklearn.svm import SVC
from sklearn.pipeline import Pipeline
from sklearn.preprocessing import StandardScaler
from sklearn.model_selection import cross_val_score

rng = np.random.default_rng(0)
n = 1000
a = rng.normal(size=n)
b = rng.normal(size=n)
y = (a + b > 0).astype(int)

for label, X in (("both in units",  np.column_stack([a, b])),
                 ("b x 1000",       np.column_stack([a, b * 1000])),
                 ("b x 0.001",      np.column_stack([a, b * 0.001]))):
    raw = cross_val_score(SVC(kernel="rbf"), X, y, cv=5).mean()
    scaled = cross_val_score(Pipeline([("s", StandardScaler()),
                                       ("m", SVC(kernel="rbf"))]), X, y, cv=5).mean()
    print(f"{label:16s} raw {raw:.4f}   scaled {scaled:.4f}")
# both in units    raw 0.9890   scaled 0.9890
# b x 1000         raw 0.7780   scaled 0.9890
# b x 0.001        raw 0.7450   scaled 0.9890

3. What SVR’s epsilon does. Sweep epsilon from 0 to 4 on the noisy sine example below. Report the support fraction and cross-validated mean absolute error (MAE) for each value.

Find where fewer support vectors come with a material increase in held-out error.

import numpy as np
from sklearn.svm import SVR
from sklearn.model_selection import cross_val_score

rng = np.random.default_rng(1)
n = 500
X = rng.uniform(-3, 3, (n, 1))
y = np.sin(X.ravel()) * 3 + rng.normal(size=n)          # noise sd = 1

for eps in (0.0, 0.1, 0.5, 1.0, 2.0, 4.0):
    m = SVR(kernel="rbf", C=10.0, epsilon=eps)
    mae = -cross_val_score(m, X, y, cv=5,
                           scoring="neg_mean_absolute_error").mean()
    m.fit(X, y)
    print(f"epsilon {eps:4.1f}  support vectors {len(m.support_):4d}"
          f" ({len(m.support_) / n:5.1%})   cv MAE {mae:.4f}")
# epsilon  0.0  support vectors  500 (100.0%)   cv MAE 0.8676
# epsilon  0.1  support vectors  463 (92.6%)   cv MAE 0.8656
# epsilon  0.5  support vectors  336 (67.2%)   cv MAE 0.8620
# epsilon  1.0  support vectors  175 (35.0%)   cv MAE 0.8541
# epsilon  2.0  support vectors   29 ( 5.8%)   cv MAE 0.8665
# epsilon  4.0  support vectors    3 ( 0.6%)   cv MAE 1.3474