What it actually does
Of all the hyperplanes that separate two classes, take the one that sits as far as possible from the nearest points on either side. That clearance is the margin, and maximizing it is a bet worth stating plainly: a boundary with room around it is more likely to survive new data than one that squeaks between training points. Two things follow. Only the points on or inside the margin — the support vectors — affect the answer, so most of the training set could be deleted without changing the model. And because the solution depends on the data only through inner products, replacing that inner product with a kernel produces a nonlinear boundary without ever computing the coordinates it lives in.
- The classes are separable, or nearly so, in the feature space the kernel implies — the kernel choice IS the modelling assumption
- Features are scaled comparably, since both the margin and the RBF kernel are defined through distances
- A geometric margin is a meaningful notion for the problem; where classes overlap heavily, the boundary is arbitrary and a probabilistic model is the better fit
- The two error types are equally costly unless class weights say otherwise — the objective has no notion of a base rate
- The sample is small enough that an n-by-n kernel matrix is affordable, which is a constraint on the model, not on the hardware
Margin objective
Make the margin as wide as possible, and charge C for every unit by which a point intrudes into it. Large C makes violations expensive, so the boundary contorts to avoid them; small C buys a wider, smoother margin that tolerates mistakes. Written as a loss the constraint becomes a hinge — max(0, 1 - y·f(x)) — and the defining property of that hinge is that a point classified correctly and comfortably contributes exactly zero. That is why the solution is sparse in the training data, and why this is a margin objective rather than a loss on probabilities: it is indifferent to how right it already is.
- inverse margin width — minimizing it is literally maximizing the gap
- slack: how far example i intrudes past its margin, zero for a comfortably correct point
- the price of a unit of intrusion, and the entire bias-variance dial
- the feature map, which is never evaluated explicitly — only its inner product is
- label in {-1, +1}; the margin is defined on signs, so 0/1 labels break the objective silently
How it is minimized
Sequential minimal optimization (SMO) on the dual QP; stochastic subgradient descent on the primal for the linear case
- dual multiplier for example i; non-zero exactly for the support vectors
- the kernel — an inner product in feature space, computed without visiting that space
- box constraint; a multiplier pinned at C marks a margin violator
- the equality constraint that forces SMO to move multipliers in pairs
The primal carries a d-dimensional weight vector and n slacks; the dual carries n multipliers and touches the data only through K(x_i, x_j). That substitution is the kernel trick, and it is the whole reason the dual is what gets solved: in the primal you cannot represent an infinite-dimensional feature map, and in the dual you never need to. SMO then exploits the equality constraint — you cannot move one multiplier without violating it, so the smallest legal working set is two, and a two-variable box-constrained QP has an analytic solution. That is why an SVM solver is a loop over pairs rather than a call to a general-purpose QP library.
C
Cost of a margin violation. Small C underfits toward a wide smooth margin, large C toward memorization (log grid from 1e-3 to 1e3, cross-validated jointly with gamma)
kernel
The model assumption. Linear for high-dimensional sparse data, RBF as the nonlinear default, polynomial rarely, custom where the domain supplies a similarity
gamma (RBF)
Inverse width of the kernel. Too large and every point becomes its own support vector (scale heuristic 1/(d · var(X)) as the centre of a log grid)
epsilon (SVR only)
Width of the insensitive tube — errors smaller than epsilon cost nothing, which is what makes SVR sparse
nu (one-class)
Upper bound on the outlier fraction and lower bound on the support-vector fraction — an interpretable dial, unlike C
class_weight
Rescales C per class; the standard handling for imbalance, since the margin itself has no notion of a base rate
The dual is a convex QP with a unique optimum, so any correct solver reaches the same answer and there are no seeds or local minima to worry about. The failure modes are about cost and calibration rather than correctness. Cost: training scales roughly between O(n^2) and O(n^3) and the kernel matrix is O(n^2) memory, so an SVM stops being trainable somewhere around 10^5 samples — the single most important practical fact about the model, and the reason it lost the tabular default to gradient boosting. Calibration: the decision function is a signed distance, not a probability, and converting it requires Platt scaling fitted on held-out data. And with an RBF kernel a large gamma memorizes: every point becomes its own support vector, training accuracy hits one, and the model is pure variance.
Training: O(n^2 d) to O(n^3 d) depending on how much of the kernel matrix is cached, with O(n^2) memory for the full Gram. Inference: O(n_sv · d) per query for a kernel SVM, so latency grows with the number of support vectors — and on noisy data that count grows roughly linearly with n, meaning more training data makes the model slower to serve as well as slower to fit. A linear SVM collapses to one dot product, O(d), which is why the linear and kernel variants are operationally different models rather than settings of one.
How, where, and why it is used
Time-Series Forecasting
AdaptedSupport vector regression on a lagged design matrix: the same supervised table any regressor would use, fitted with an epsilon-insensitive tube so residuals smaller than epsilon cost nothing. That tube is the reason to reach for it — it produces a forecaster that ignores small noise entirely and is therefore less prone to chasing it than a squared-error model.
- Short-horizon nonlinear forecasting on small samples, where a flexible model would overfit and an RBF kernel with a wide margin does not
- Financial and energy series in the pre-boosting era, where SVR was the standard nonlinear benchmark and remains a reference point in the literature
- Problems with a few thousand observations and a nonlinear driver relationship, which is the regime where the O(n^2) cost is not yet a problem
The honest framing is historical and still useful: SVR is a genuinely strong small-sample nonlinear forecaster, and it has been largely displaced by gradient boosting on tabular lag features, which trains faster, handles missing values, and does not carry a quadratic memory cost. Reach for it when n is small and the relationship is smooth and nonlinear; do not reach for it on a million rows, because it cannot be fitted there at all.
- Scale features and the target — an RBF kernel is a function of squared distance, so unscaled inputs make gamma meaningless
- Lags, rolling statistics and Fourier seasonal terms, exactly as for any supervised forecaster
- Difference or detrend first: SVR cannot extrapolate a trend, since its prediction is a bounded combination of kernel evaluations
- Choose epsilon relative to the noise scale of the target, not as a fixed default
Rolling-origin backtesting with C, gamma and epsilon tuned inside each fold — tuning them once on the full series is a leak, and with three interacting hyperparameters it is a large one. Score MASE against seasonal-naive.
- Scaling the target on the full series before splitting, which leaks the future through the scaler
- Trend: the model cannot predict outside the range of the training targets, so a trending series produces forecasts that flatten
- Tuning C and gamma independently — they interact strongly, and a 1-D sweep on each finds a worse optimum than a coarse 2-D grid
Anomaly Detection
Primary fitOne-class SVM turns the margin objective inside out: with no labels, find the smallest region containing most of the data — a hyperplane separating the data from the origin in feature space, or equivalently (SVDD) the smallest enclosing sphere. The parameter nu directly bounds the fraction of training points allowed outside, so the contamination rate is something you set rather than something you discover.
- Novelty detection where a clean reference period exists and anomalies are unlabelled and rare
- Industrial process and equipment monitoring with a well-characterized normal operating envelope
- Network intrusion detection against a profile of normal traffic
- Screening inputs to a deployed model for being outside its training distribution
It gives a closed, nonlinear boundary around normal behaviour with a directly interpretable contamination parameter, and the decision is a geometric one rather than a density estimate — which is what makes it work in moderate dimensions where density estimation has already failed. The reasons against it are concrete: it is O(n^2), it is sensitive to gamma in a way that is hard to tune without labels, and on large or heavily contaminated data an isolation forest is both faster and more robust.
- Scale on a confirmed-clean window and persist those statistics — the boundary is defined in scaled space
- Reduce dimension before fitting; the enclosing boundary loosens quickly as dimension grows
- Set nu from the believed contamination rate rather than tuning it for a target alert volume, which conflates the model with the threshold
Precision@k and PR-AUC against confirmed incidents, holding out entire anomaly episodes rather than individual points. Validate gamma by checking the support-vector fraction: a fit where most points are support vectors has wrapped itself around the training data and will flag everything new.
- Contaminated training data pulls the boundary outward until it encloses the anomalies
- Gamma too large produces a boundary that hugs every training point — perfect on train, useless on anything unseen
- Quadratic training cost makes periodic refits on a growing reference window quietly unaffordable
Optimization
ViableThe SVM is a constrained convex program studied as one: a quadratic objective under box constraints and one equality constraint. Working through it is where the applied vocabulary of optimization is actually acquired — Lagrangian duality, complementary slackness, KKT conditions as a convergence test rather than a formality, and active-set reasoning, since the support vectors are precisely the active constraints.
- The canonical worked example of duality in applied machine learning — the primal-to-dual transformation is what makes the kernel possible
- Decomposition methods: SMO is an active-set method with a working set of two, chosen because the equality constraint forbids one
- KKT violation as a principled stopping rule, in contrast to "the objective stopped moving"
It is the smallest model in this section where the constraint structure, rather than the objective, determines the algorithm. You cannot pick a working set of one, because the equality constraint couples the multipliers; you can solve two analytically; therefore the solver is a loop over pairs. That chain of reasoning — from constraint geometry to algorithm — is the transferable content, and it is why the classical-ml verdict for this domain names the SVM alongside Lasso.
- Scale features so the box constraints and the kernel width are comparable across dimensions
- Use class weights to rescale C per class rather than resampling, which changes the problem rather than the objective
Check the KKT conditions at the returned solution: every multiplier at zero should correspond to a point outside the margin, every multiplier at C to a violator, and every interior multiplier to a point exactly on it. That is an exact correctness test, unlike an iteration budget.
- Treating a stalled objective as convergence when KKT violations remain on points the working-set heuristic has not revisited
- Ignoring the O(n^2) kernel matrix in the problem formulation, so the model is chosen before anyone checks it can be fitted
- Reading the dual multipliers as importances — they are constraint activities, not effect sizes
Computer Vision
AdaptedHistorically: hand-engineered descriptors — SIFT, HOG, bag-of-visual-words — fed to an SVM, which was the state of the art in object recognition until 2012. Currently: a linear SVM as a classification head on frozen features from a pretrained network, which is the standard way to get a strong classifier from a few hundred labelled images.
- Small-data image classification on top of a frozen pretrained backbone, where fine-tuning would overfit
- Pedestrian and object detection with HOG features, still deployed in embedded and low-power settings
- Linear probing — the standard diagnostic for how linearly separable a learned representation is
On raw pixels it is genuinely superseded, and the honest statement is that learned features beat designed ones decisively once data and compute allowed. What survives is the regime where the SVM was always strongest: few labelled examples, high-dimensional inputs, and a fixed representation. A linear SVM on frozen embeddings trains in seconds on a few hundred images and reliably beats a fine-tuned network at that sample size, which is a narrow but real and frequently encountered case.
- L2-normalize descriptors and embeddings before fitting; margin geometry is meaningless across inconsistent vector magnitudes
- Prefer a linear kernel on high-dimensional embeddings — RBF adds cost and rarely accuracy once the features are already learned
- For classical descriptors, pool to a fixed-length representation before the SVM, which has no notion of spatial structure
Per-class accuracy and confusion structure rather than a single top-1 number, since one-vs-rest decomposition hides which pairs are actually confused. For linear probing, report the frozen-feature baseline alongside the fine-tuned result or the number means nothing.
- Multiclass by one-vs-rest with uncalibrated decision values, so scores from different binary problems are compared on incomparable scales
- Quadratic training cost, which makes even a modest image dataset infeasible without subsampling
- Assuming an SVM adds representational power to a backbone; it adds a decision rule, and the features set the ceiling
Natural Language
ViableA linear SVM on tf-idf features is the classical strong baseline for text classification: documents are very high-dimensional and very sparse, which is exactly the regime where a max-margin linear boundary generalizes well and a kernel adds nothing. Training on sparse features is cheap enough that the usual scaling objection does not apply.
- Topic and intent classification with a few thousand labelled documents
- Spam, abuse and content-policy classifiers where the deliverable must be inspectable
- The baseline any transformer-based classifier is expected to beat before its cost is justified
With d far larger than n, a max-margin boundary is exactly the right inductive bias: there are many separating hyperplanes and the widest-margin one generalizes best. It is also close to free to train and trivial to serve. What it cannot do is anything requiring word order, negation, or semantics beyond the vocabulary — for which a pretrained language model is the answer, and the honest question is whether the accuracy gain justifies the operational cost on the task at hand.
- tf-idf with sublinear term frequency and L2-normalized rows
- Character n-grams for noisy, multilingual or adversarial text where tokenization is unreliable
- Keep the sparse representation end to end — densifying a tf-idf matrix is what makes people conclude, wrongly, that SVMs cannot handle text
Macro-F1 rather than accuracy when classes are imbalanced, on a temporal split where language drifts. Inspect the top-weighted features per class — with a linear kernel they are directly readable, and this is where label noise and leakage usually surface.
- Fitting the vectorizer on the full corpus before splitting, which leaks vocabulary and document frequencies
- Reading decision-function values as confidences without Platt scaling
- Losing sparsity through a dense transformation, which turns a cheap fit into an infeasible one
What it costs to run
Superlinear in samples and quadratic in memory: the kernel matrix alone is n^2 doubles, which is roughly 80 GB at n = 100,000. This decides whether an SVM is a candidate at all, and it decides it before any accuracy comparison is run. A linear SVM on sparse features escapes this entirely and trains in seconds.
A linear SVM is one dot product. A kernel SVM evaluates the kernel against every support vector, so latency is proportional to the support-vector count — which is itself proportional to how much the classes overlap. A model that retained 40% of its training set as support vectors is a slow model, and the training log says so before production does.
Infrequent by necessity rather than preference: a full refit is expensive and there is no incremental update that preserves the exact solution. Where data arrives continuously, the primal form fitted by stochastic subgradient descent is the variant that can keep up, at the cost of the kernel.
- Track the support-vector fraction across refits — a rising fraction means the classes are overlapping more, and it predicts an inference-latency regression before it predicts an accuracy one
- Monitor the distribution of decision values near zero; a thickening band at the boundary is drift arriving before any accuracy metric moves
- If Platt scaling is deployed, recalibrate it on its own schedule — the sigmoid drifts well before the boundary does
- Persist and monitor the feature scaler as a first-class artefact, since the RBF kernel width is defined in scaled units
- Labels must be -1 and +1; passing 0/1 does not error, it silently optimizes a different objective
- Feature scaling is mandatory, and gamma is expressed in squared scaled units — changing the scaler invalidates every tuned gamma
- decision_function returns a signed distance, not a probability; anything downstream that thresholds it as a confidence is wrong unless Platt scaling was fitted on held-out data
- Multiclass is one-vs-rest or one-vs-one wrapped around binary SVMs, so multiclass "probabilities" are two transformations away from the objective that was actually optimized
- The support vectors are literal training rows shipped inside the model, which makes deletion requests and data licensing a model problem, exactly as with kNN
Pros and cons, in context
Margin maximization generalizes well from small samples
Decisive when n is small and d is large — text, genomics, embeddings — the regime where a flexible model overfits and this one does not. The advantage fades as n grows, because at that point the data itself controls variance and cheaper models catch up.
The kernel trick buys nonlinearity without an explicit feature map
You get to work in an effectively infinite-dimensional space at the cost of an n-by-n matrix. An excellent trade while n is small and a disqualifying one when it is not — the same condition, stated twice, governs when this model is appropriate.
The solution is sparse in the training points
Only support vectors survive into the model, which makes it compact when classes are well separated. On noisy, overlapping data the support-vector count balloons and the sparsity advantage inverts into an inference-cost problem.
Convex objective with a unique optimum
No seeds, no local minima, reproducible fits across runs and machines — which matters more in audited pipelines than the last point of accuracy, and is the same argument that keeps linear models deployed.
Does not scale: quadratic memory and superlinear time in the sample count
A hard ceiling near 10^5 samples that no amount of tuning moves. This is the single fact that decides most SVM-versus-boosting arguments before accuracy is discussed, and it should be checked first.
No native probabilities
The decision function is a signed distance. Getting calibrated probabilities requires Platt scaling on held-out data — an extra fitting step that is routinely skipped, after which people read distances as confidences and threshold them.
C and gamma interact strongly and must be tuned jointly
The search is genuinely two-dimensional, and a default RBF fit at C = 1 is usually a poor showing rather than a fair baseline. Comparisons against SVMs are frequently unfair for exactly this reason.
Superseded by gradient boosting on tabular data
Boosting generally wins on accuracy, training time, missing-value handling, and interpretability of feature importance. The SVM keeps the high-dimensional small-sample regime, and honestly not much outside it.
Intuitive, idiomatic, optimized
The same algorithm in three languages, each progressing through Kent Beck's three stages. Every optimization names the lever it pulls and the trade-off it accepts — cost profiles are illustrative, not measured benchmarks.
The math, transcribed literally. No dependencies.
"""Soft-margin SVM by subgradient descent - the primal, transcribed.
Minimize (1/2)||w||^2 + C * sum_i max(0, 1 - y_i (w . x_i + b)).
The hinge is not differentiable at the margin, so this takes a subgradient: a
point inside the margin contributes -C * y_i * x_i, and a point safely outside
contributes nothing at all. That single if-statement is the whole difference
between this and a least-squares fit.
"""
def fit(X, y, C=1.0, lr=0.001, epochs=1_000):
"""y must be in {-1, +1}: the margin is defined on signs, not on 0/1."""
n = len(X)
d = len(X[0])
w = [0.0] * d
b = 0.0
for _ in range(epochs):
grad_w = [0.0] * d
grad_b = 0.0
for i in range(n):
# margin = y_i (w . x_i + b)
score = b
for j in range(d):
score += w[j] * X[i][j]
margin = y[i] * score
if margin < 1.0: # hinge is active
for j in range(d):
grad_w[j] -= C * y[i] * X[i][j]
grad_b -= C * y[i]
# The L2 term contributes w itself; the hinge term was accumulated above.
for j in range(d):
w[j] -= lr * (w[j] + grad_w[j])
b -= lr * grad_b
return w, b
def decision_function(w, b, x):
score = b
for j in range(len(w)):
score += w[j] * x[j]
return scoreO(n*d) per epoch in interpreter loops. Linear only: the primal carries an explicit w, so there is nowhere for a kernel to go.
How this language actually wants it written.
"""Kernel SVM by sequential minimal optimization - typed, validated, dual."""
from collections.abc import Callable
from dataclasses import dataclass
import numpy as np
from numpy.typing import NDArray
Vector = NDArray[np.float64]
Matrix = NDArray[np.float64]
Kernel = Callable[[Matrix, Matrix], Matrix]
def linear_kernel() -> Kernel:
return lambda a, b: a @ b.T
def rbf_kernel(gamma: float) -> Kernel:
"""K(a, b) = exp(-gamma ||a - b||^2), via the squared-norm expansion."""
def kernel(a: Matrix, b: Matrix) -> Matrix:
squared = (
np.einsum("ij,ij->i", a, a)[:, None]
- 2.0 * (a @ b.T)
+ np.einsum("ij,ij->i", b, b)[None, :]
)
# Cancellation in the expansion can produce small negatives; clip before
# exponentiating rather than letting them become spurious similarities.
return np.exp(-gamma * np.maximum(squared, 0.0))
return kernel
@dataclass(frozen=True)
class SvmModel:
"""Only the support vectors survive - every other alpha is exactly zero."""
support_vectors: Matrix
dual_coefficients: Vector # alpha_i * y_i, one per support vector
intercept: float
kernel: Kernel
@property
def n_support(self) -> int:
return int(self.support_vectors.shape[0])
def decision_function(self, X: Matrix) -> Vector:
"""A signed distance to the boundary. NOT a probability."""
if X.ndim != 2 or X.shape[1] != self.support_vectors.shape[1]:
raise ValueError(
f"expected (n, {self.support_vectors.shape[1]}) input, got {X.shape}"
)
return self.kernel(X, self.support_vectors) @ self.dual_coefficients + self.intercept
def predict(self, X: Matrix) -> Vector:
return np.sign(self.decision_function(X))
def _bounds(y_i: float, y_j: float, a_i: float, a_j: float, C: float) -> tuple[float, float]:
"""The segment alpha_j may move along without breaking sum(alpha y) = 0."""
if y_i != y_j:
return max(0.0, a_j - a_i), min(C, C + a_j - a_i)
return max(0.0, a_i + a_j - C), min(C, a_i + a_j)
def fit(
X: Matrix,
y: Vector,
C: float = 1.0,
kernel: Kernel | None = None,
max_passes: int = 10,
tol: float = 1e-3,
) -> SvmModel:
"""Simplified SMO on the dual. Raises ValueError on malformed input."""
if X.ndim != 2:
raise ValueError(f"X must be 2-D, got shape {X.shape}")
if X.shape[0] != y.shape[0]:
raise ValueError(f"X has {X.shape[0]} rows but y has {y.shape[0]}")
if not np.isin(y, (-1.0, 1.0)).all():
raise ValueError("labels must be -1 or +1; 0/1 optimizes a different objective")
if C <= 0.0:
raise ValueError(f"C must be positive, got {C}")
# None rather than a default kernel object: a mutable default would be
# shared across every call.
active_kernel = linear_kernel() if kernel is None else kernel
n = X.shape[0]
alpha = np.zeros(n, dtype=np.float64)
intercept = 0.0
rng = np.random.default_rng(0)
passes = 0
while passes < max_passes:
changed = 0
for i in range(n):
# Kernel rows are computed on demand: the full Gram is O(n^2) and
# is the thing that puts a ceiling on this model.
row_i = active_kernel(X[i : i + 1], X)[0]
error_i = float(row_i @ (alpha * y)) + intercept - y[i]
violates_kkt = (y[i] * error_i < -tol and alpha[i] < C) or (
y[i] * error_i > tol and alpha[i] > 0.0
)
if not violates_kkt:
continue
j = int(rng.integers(n - 1))
j = j + 1 if j >= i else j
row_j = active_kernel(X[j : j + 1], X)[0]
error_j = float(row_j @ (alpha * y)) + intercept - y[j]
old_i, old_j = float(alpha[i]), float(alpha[j])
low, high = _bounds(float(y[i]), float(y[j]), old_i, old_j, C)
if high - low < 1e-12:
continue
# Second derivative of the objective along the alpha_j direction.
eta = 2.0 * row_i[j] - row_i[i] - row_j[j]
if eta >= 0.0:
continue
alpha[j] = float(np.clip(old_j - y[j] * (error_i - error_j) / eta, low, high))
if abs(alpha[j] - old_j) < 1e-5:
continue
alpha[i] = old_i + y[i] * y[j] * (old_j - alpha[j])
# Restore the KKT condition on b using whichever multiplier is
# unbound, since an unbound point sits exactly on the margin.
b_i = intercept - error_i - y[i] * (alpha[i] - old_i) * row_i[i] - y[j] * (
alpha[j] - old_j
) * row_i[j]
b_j = intercept - error_j - y[i] * (alpha[i] - old_i) * row_i[j] - y[j] * (
alpha[j] - old_j
) * row_j[j]
if 0.0 < alpha[i] < C:
intercept = float(b_i)
elif 0.0 < alpha[j] < C:
intercept = float(b_j)
else:
intercept = float((b_i + b_j) / 2.0)
changed += 1
passes = passes + 1 if changed == 0 else 0
support = alpha > 1e-8
return SvmModel(
support_vectors=X[support],
dual_coefficients=alpha[support] * y[support],
intercept=intercept,
kernel=active_kernel,
)What changed: This is a change of formulation, not of style, and it is the point of the entry: subgradient descent on the primal carries an explicit weight vector, so it can never use a kernel. The dual touches the data only through inner products, which is what makes a nonlinear boundary expressible at all — so the algorithm becomes SMO over pairs of multipliers, paired because the equality constraint forbids moving one alone. Around that, the idiomatic changes: the kernel is an injected callable rather than hard-coded arithmetic, the fitted model is a frozen dataclass retaining only support vectors, labels are validated as ±1 (the failure that otherwise optimizes a different objective in silence), and the KKT test is a guard clause rather than nested conditionals.
O(n*d) per kernel row, two rows per candidate pair. Memory stays O(n) because the Gram is never materialized.
Faster, and the justification for each change.
"""Kernel SVM - precomputed Gram, maintained error cache, in-place updates."""
from collections.abc import Callable
import numpy as np
from numpy.typing import NDArray
Vector = NDArray[np.float64]
Matrix = NDArray[np.float64]
Kernel = Callable[[Matrix, Matrix], Matrix]
def fit_cached(
X: Matrix,
y: Vector,
kernel: Kernel,
C: float = 1.0,
max_passes: int = 10,
tol: float = 1e-3,
) -> tuple[Vector, float]:
"""SMO with the two structural accelerations every real solver uses.
First, the Gram matrix is computed once as a single kernel call over the
whole dataset - one GEMM inside the kernel rather than 2n row evaluations
per pass. The cost is O(n^2) memory, and that is not an implementation
detail: it is the reason SVMs stop being trainable somewhere near 100k
samples, and it should be checked before the model is chosen.
Second, an error cache. E_i depends on every alpha, so recomputing it from
a Gram row is O(n) per candidate; but a pair update changes exactly two
alphas, so every cached error can be corrected with one vector expression.
"""
design = np.ascontiguousarray(X, dtype=np.float64)
labels = np.ascontiguousarray(y, dtype=np.float64)
n = design.shape[0]
gram = np.ascontiguousarray(kernel(design, design), dtype=np.float64)
alpha = np.zeros(n, dtype=np.float64) # allocated once
errors = -labels.copy() # E_i at alpha = 0, b = 0
scratch = np.empty(n, dtype=np.float64) # reused per update
intercept = 0.0
rng = np.random.default_rng(0)
passes = 0
while passes < max_passes:
changed = 0
for i in range(n):
if not (
(labels[i] * errors[i] < -tol and alpha[i] < C)
or (labels[i] * errors[i] > tol and alpha[i] > 0.0)
):
continue
# Second-choice heuristic: the pair with the largest error gap makes
# the largest step, so the working set is chosen rather than sampled.
j = int(np.argmax(np.abs(errors - errors[i])))
if j == i:
j = int(rng.integers(n - 1))
j = j + 1 if j >= i else j
old_i, old_j = float(alpha[i]), float(alpha[j])
if labels[i] != labels[j]:
low, high = max(0.0, old_j - old_i), min(C, C + old_j - old_i)
else:
low, high = max(0.0, old_i + old_j - C), min(C, old_i + old_j)
if high - low < 1e-12:
continue
eta = 2.0 * gram[i, j] - gram[i, i] - gram[j, j]
if eta >= 0.0:
continue
new_j = min(max(old_j - labels[j] * (errors[i] - errors[j]) / eta, low), high)
if abs(new_j - old_j) < 1e-5:
continue
new_i = old_i + labels[i] * labels[j] * (old_j - new_j)
delta_i = labels[i] * (new_i - old_i)
delta_j = labels[j] * (new_j - old_j)
b_i = intercept - errors[i] - delta_i * gram[i, i] - delta_j * gram[i, j]
b_j = intercept - errors[j] - delta_i * gram[i, j] - delta_j * gram[j, j]
if 0.0 < new_i < C:
new_intercept = float(b_i)
elif 0.0 < new_j < C:
new_intercept = float(b_j)
else:
new_intercept = float((b_i + b_j) / 2.0)
# errors += d_i * K[i] + d_j * K[j] + (b_new - b_old), fused into
# the reusable scratch buffer so no temporary row is materialized.
np.multiply(gram[i], delta_i, out=scratch)
scratch += delta_j * gram[j]
scratch += new_intercept - intercept
errors += scratch
alpha[i], alpha[j] = new_i, new_j
intercept = new_intercept
changed += 1
passes = passes + 1 if changed == 0 else 0
return alpha, interceptWhat changed: Two accelerations, both structural. The kernel matrix is evaluated once for the whole dataset instead of two rows per candidate pair, which moves the work into one BLAS-backed call — and buys it with O(n^2) memory, the constraint that defines where this model can be used at all. Then the error cache: E_i depends on every multiplier, so recomputing it per candidate was the dominant cost, but a pair update perturbs exactly two multipliers, so every cached error can be corrected by one fused vector expression. Working-set selection also stops being random and picks the maximum-error-gap partner, which is what turns simplified SMO into something that converges in a sensible number of passes.
Vectorize to NumPy/BLAS (@, np.dot, einsum)
The Gram matrix is one kernel call over the full dataset, which for linear and RBF kernels bottoms out in a single GEMM rather than n separate row evaluations per pass.
Trade-off: O(n^2) memory — roughly 80 GB at n = 100,000 — which is precisely the ceiling that rules SVMs out on large datasets. This optimization buys speed by spending the resource the model is already short of.
Eliminate Python-level loops over samples
Every cached error is corrected by one vector expression after a pair update, replacing an O(n) per-candidate recomputation from a kernel row.
Trade-off: The cache is incrementally maintained, so floating-point error accumulates over many updates and a long run needs a periodic full refresh — and a stale cache produces a wrong answer rather than a slow one.
Pre-allocate output arrays and use in-place operations
The multipliers, error vector and update scratch are allocated once; the error correction writes through the scratch buffer rather than materializing a temporary row per update.
Trade-off: The in-place chain is order-dependent and the scratch buffer is shared state, so reordering the update lines corrupts the errors silently instead of raising.
O(n^2 d) once for the Gram, then O(n) per multiplier-pair update. Illustrative, not a measured benchmark.