Support Vector Machines & Kernel Methods

Max-Margin and the Hinge Loss

The previous example changed the penalty; this one changes the loss. A support vector machine replaces squared error with the hinge loss and seeks the separating hyperplane with the widest margin. The hinge is exactly zero once a point sits on the correct side of the margin and grows linearly inside it, so most of the data contributes nothing at all to the fit: on the two-blob toy only 18 of 240 points satisfy $y_i f(x_i) < 1$. The other 222 could be moved anywhere on their own side of the margin without shifting the boundary by a hair — which is what support vector means.

PegasosSVM optimises that objective by stochastic sub-gradient descent with a decaying step $\eta_t = 1/(\lambda t)$. It reaches AUC 0.704 on credit default against scikit-learn's exact LinearSVC at 0.711 — the same hinge objective, SGD against an exact solver.

Does max-margin actually help?

That comparison is worth stating carefully, because the SVM does not win here. The exact linear SVM (0.711) ties an ordinary logistic regression (0.715); the 0.011 the from-scratch version trails is the optimiser, not the loss. R tells the same story more sharply — e1071 scores 0.699 against a logistic fit to the very same 3,000 rows at 0.713. Subsampling is not the excuse: the logistic barely moves between 3,000 rows and 21,000 (0.713 to 0.712). On heavily overlapping, near-linearly-separable-in-no-sense data, max-margin buys nothing over maximum likelihood.

The kernel trick

The kernel trick makes the linear machine nonlinear for free. Replace every inner product $x\cdot x'$ with $k(x,x')=\langle\phi(x),\phi(x')\rangle$ and the method operates in a rich implicit feature space $\phi$ that is never constructed. On the two-moons data the linear SVM is helpless and the RBF kernel separates cleanly. On the real credit data the kernel lifts the SVM only back to the logistic baseline (0.715) and stops — there is little smooth nonlinear structure in a noisy, 22%-imbalanced problem, and the curved boundary costs $O(n^2)$ to find out.

For regression the kernel attaches to ridge instead, and the payoff is real: kernel ridge and SVR cut California test RMSE from the linear model's 0.737 to 0.647 and 0.644. The from-scratch closed form matches scikit-learn to 3×10⁻¹⁴, since it solves the identical linear system.

Kernel Ridge Is a Gaussian Process

Which sets up the result the whole subsection is built around. A Gaussian process with kernel $k$ and noise variance $\sigma^2$ has posterior mean $\bar f(x_*) = k(x_*,X)[K+\sigma^2 I]^{-1}y$ — exactly kernel ridge regression with $\lambda=\sigma^2$. Verified against scikit-learn's GP to 3.6×10⁻¹⁴ and against kernlab::gausspr in R to 5.5×10⁻¹⁴. The frequentist kernel machine and the Bayesian nonparametric one are not cousins; they are the same estimator.

Then why does the GP score better?

That identity leaves something to explain, and chasing it down is the most useful thing in this example. Kernel ridge scores 0.647 above, while the Gaussian-process example reports 0.611 — how can the same estimator give two answers? It cannot. Choosing the same two numbers by 5-fold cross-validation on the training subsample alone, with the test set never touched, gives $\gamma=0.125,\ \alpha=0.1$ and test RMSE 0.612. The gap was never the estimator; it was that the GP selected its length-scale and noise by maximising the marginal likelihood while kernel ridge had them handed to it. Left at $\gamma=1/p$ and $\alpha=1$, the same machine gives up 0.036 RMSE. What the Bayesian formulation buys here is not a better predictor but hyperparameters read off the data rather than guessed — and a posterior variance on top.

A scaling default worth knowing

A second discrepancy has a plainer cause and is worth knowing about. R's SVR scores 0.601 where Python's scores 0.644 — outside anything resampling explains. The reason is a default: e1071::svm scales its inputs unless told otherwise, and with scale=FALSE the identical call returns 0.656. This is not cosmetic. An isotropic RBF kernel measures distance in whatever units it is handed, so rescaling the features is a change of kernel. AveOccup is what makes it bite: its full-sample SD is set by a handful of extreme blocks, leaving the standardised subsample nearly constant on that axis (SD 0.11) until e1071 restandardises it back to 1.

Where this sits

The threads meet here. Penalties are priors, so ridge is a Gaussian prior and the Variable Selection arc is the fully Bayesian version of the lasso; kernels are covariance functions, so kernel ridge is a Gaussian process and the smoothing penalty in penalised splines is another face of the same object. And the standing practical caveat: kernel machines share the GP's $O(n^2\!-\!n^3)$ cost, which forces the subsampling used throughout and is the recurring reason gradient boosting, not the kernel machine, is the default for large tabular data.

Notebooks

Downloads

Kernel-Machine Module — Source Code

"""From-scratch support-vector and kernel machines.

Two load-bearing algorithms, built by hand and validated against scikit-learn:

* PegasosSVM -- the soft-margin linear support-vector classifier trained by
  stochastic sub-gradient descent on the hinge-loss objective
      (lambda/2)||w||^2 + (1/n) sum_i max(0, 1 - y_i (w.x_i)),
  Shalev-Shwartz et al., "Pegasos" (ICML 2007).  The max-margin idea and the
  hinge loss are the whole story of the SVM; this is the optimiser.

* KernelRidge -- kernel ridge regression in closed form,
      alpha = (K + n*lambda*I)^{-1} y,   f(x) = sum_i alpha_i k(x, x_i),
  which solves the SAME linear system as scikit-learn's KernelRidge (so the
  fits match to numerical tolerance) and, with an RBF kernel and noise
  variance = n*lambda, equals the POSTERIOR MEAN of a Gaussian process --
  the exact bridge to the Bayesian-nonparametric arc demonstrated in the
  notebook.

The kernel trick: replace every inner product x.x' by k(x,x') = <phi(x),phi(x')>
and a linear method becomes nonlinear without ever forming phi.  rbf_kernel and
poly_kernel below are the two used in the notebook.
"""
import numpy as np


def rbf_kernel(A, B, gamma):
    """Gaussian / squared-exponential kernel exp(-gamma ||a-b||^2)."""
    a2 = np.sum(A**2, axis=1)[:, None]
    b2 = np.sum(B**2, axis=1)[None, :]
    sq = np.maximum(a2 + b2 - 2 * A @ B.T, 0.0)
    return np.exp(-gamma * sq)


def poly_kernel(A, B, degree=3, coef0=1.0, gamma=1.0):
    """Polynomial kernel (gamma <a,b> + coef0)^degree."""
    return (gamma * (A @ B.T) + coef0) ** degree


class PegasosSVM:
    """Linear soft-margin SVM via the Pegasos stochastic sub-gradient method.

    Labels are taken in {0,1} and mapped internally to {-1,+1}.  `lam` is the
    L2 regularisation strength (larger = wider margin, more regularised).
    decision_function returns the signed distance w.x + b; predict thresholds
    it at 0.
    """

    def __init__(self, lam=1e-4, n_epochs=20, seed=0):
        self.lam = lam
        self.n_epochs = n_epochs
        self.seed = seed

    def fit(self, X, y):
        X = np.asarray(X, float); n, d = X.shape
        yy = np.where(np.asarray(y) > 0, 1.0, -1.0)
        Xb = np.hstack([X, np.ones((n, 1))])           # absorb bias as an extra feature
        w = np.zeros(d + 1)
        rng = np.random.default_rng(self.seed)
        t = 0
        for _ in range(self.n_epochs):
            for i in rng.permutation(n):
                t += 1
                eta = 1.0 / (self.lam * t)             # decaying step size
                if yy[i] * (Xb[i] @ w) < 1:            # inside the margin -> hinge active
                    w = (1 - eta * self.lam) * w + eta * yy[i] * Xb[i]
                else:
                    w = (1 - eta * self.lam) * w
        self.coef_ = w[:d]; self.intercept_ = w[d]
        return self

    def decision_function(self, X):
        return np.asarray(X, float) @ self.coef_ + self.intercept_

    def predict(self, X):
        return (self.decision_function(X) > 0).astype(int)


class KernelRidge:
    """Kernel ridge regression in closed form (= GP posterior mean).

    kernel: 'rbf' or 'poly'.  alpha is the ridge penalty (scikit-learn's alpha).
    Fits dual_ = (K + alpha*I)^{-1} y and predicts k(x, X) @ dual_ -- no centering
    of y, matching scikit-learn's KernelRidge, so standardize the target outside
    the class if you want a data-mean prior rather than a zero-mean one.
    """

    def __init__(self, kernel="rbf", gamma=1.0, alpha=1.0, degree=3, coef0=1.0):
        self.kernel = kernel; self.gamma = gamma; self.alpha = alpha
        self.degree = degree; self.coef0 = coef0

    def _K(self, A, B):
        if self.kernel == "rbf":
            return rbf_kernel(A, B, self.gamma)
        return poly_kernel(A, B, self.degree, self.coef0, self.gamma)

    def fit(self, X, y):
        # matches scikit-learn KernelRidge exactly (no centering): dual = (K+alpha*I)^{-1} y.
        # With an RBF kernel and alpha = noise variance this is the posterior mean of a
        # zero-mean Gaussian process -- standardize y outside for accuracy / a data-mean prior.
        self.X_ = np.asarray(X, float); y = np.asarray(y, float)
        K = self._K(self.X_, self.X_)
        n = K.shape[0]
        self.dual_ = np.linalg.solve(K + self.alpha * np.eye(n), y)
        return self

    def predict(self, X):
        return self._K(np.asarray(X, float), self.X_) @ self.dual_

References

• Cortes, C. & Vapnik, V. (1995). Support-vector networks. Machine Learning 20(3), 273–297. — the soft-margin SVM
• Shalev-Shwartz, S., Singer, Y., Srebro, N. & Cotter, A. (2011). Pegasos: primal estimated sub-gradient solver for SVM. Mathematical Programming 127(1), 3–30. — the optimiser implemented here
• Boser, B., Guyon, I. & Vapnik, V. (1992). A training algorithm for optimal margin classifiers. COLT, 144–152. — the kernel trick
• Schölkopf, B. & Smola, A. (2002). Learning with Kernels. MIT Press. — kernels, the representer theorem, and kernel ridge
• Rasmussen, C. E. & Williams, C. K. I. (2006). Gaussian Processes for Machine Learning. MIT Press. — §6.2 on the kernel-ridge/GP equivalence used above
• Chang, C.-C. & Lin, C.-J. (2011). LIBSVM: a library for support vector machines. ACM TIST 2(3). — the solver behind both scikit-learn's SVC and R's e1071