Kernelizing Your SVM for Non-linear Data

The Challenge: Why We Need a New Approach

Our previous SVM implementation solved the primal problem, where we directly optimized for the weight vector $\mathbf{w}$ and bias $b$. The decision function was simply ${\mathbf{w}}^T\mathbf{x}+b$. This approach is intuitive, but it's fundamentally tied to the original feature space of $\mathbf{x}$.

Kernels, on the other hand, operate by implicitly mapping our data $\mathbf{x}$ into a higher-dimensional feature space $\Phi (\mathbf{x})$. To work in this new space, we would need to find a weight vector ${\mathbf{w}}_{\varphi }$ in that same high-dimensional space. The "kernel trick" is that we can avoid this explicit, costly transformation by using a different mathematical formulation of the SVM problem: the dual problem.

The dual formulation rephrases the optimization problem so that it only depends on the dot products of data points, which we can then replace with our kernel function.

To get a quick and clear visual intuition for this process, let's watch a short video.

The SVM Dual Formulation

Instead of optimizing for $\mathbf{w}$, the dual problem involves optimizing a set of variables called Lagrange multipliers, denoted by ${\alpha }_i$, one for each training sample.

The dual objective function we want to maximize is:
$L_D(\alpha )={\sum }_{i=1}^n{\alpha }_i-\frac{1}{2}{\sum }_{i=1}^n{\sum }_{j=1}^n{\alpha }_i{\alpha }_j{y}_i{y}_{j}K({\mathbf{x}}_i,{\mathbf{x}}_j)$
This is subject to the constraints $0\le {\alpha }_i\le C$.

The key things to notice are:

1. The weight vector $\mathbf{w}$ is gone!
2. The data points ${\mathbf{x}}_i$ and ${\mathbf{x}}_j$ only appear inside the kernel function $K({\mathbf{x}}_i,{\mathbf{x}}_j)$. This is where we plug in our polynomial or RBF kernel.

Once we find the optimal $\alpha$ values, the prediction for a new data point $\mathbf{z}$ is given by:
$\text{prediction}=\text{sign}\left({\sum }_{i=1}^n{\alpha }_i{y}_{i}K({\mathbf{x}}_i,\mathbf{z})+b\right)$

Data points for which ${\alpha }_i>0$ are the support vectors. The decision boundary is determined entirely by them.

Implementation: Building a Kernelized SVM from Scratch

Now, let's get to the fun part: implementation. We will build a new SVM class that solves the dual problem using gradient ascent. This will be a major upgrade to your previous linear SVM.

We will follow the excellent guide from the article "Support Vector Machines for Beginners – Training Algorithms," which details this exact process.

Inspired by that guide, let's structure our new KernelSVM class.

Step 1: __init__ and Kernel Functions

First, we'll set up the class to accept a kernel as a parameter. We'll also bring in the kernel functions you wrote in the last lesson. Since we're building a new class, let's put them inside for self-containment.

import numpy as np

class KernelSVM:
    def __init__(self, kernel='rbf', C=1.0, gamma=0.1, degree=3, coef=1.0):
        self.C = C # Regularization parameter
        self.kernel_name = kernel
        self.gamma = gamma   # For RBF kernel
        self.degree = degree # For polynomial kernel
        self.coef = coef     # For polynomial kernel
        
        if self.kernel_name == 'rbf':
            self.kernel = self._rbf_kernel
        elif self.kernel_name == 'poly':
            self.kernel = self._polynomial_kernel
        else: # Default to linear kernel
            self.kernel = self._linear_kernel

        self.alphas = None
        self.b = 0
        self.X_train = None
        self.y_train = None

    def _linear_kernel(self, X1, X2):
        return X1 @ X2.T

    def _polynomial_kernel(self, X1, X2):
        return (X1 @ X2.T + self.coef) ** self.degree

    def _rbf_kernel(self, X1, X2):
        # A pure NumPy implementation without SciPy
        # This is for self-containment, though cdist is often faster
        sq_dists = np.sum(X1**2, axis=1).reshape(-1, 1) + np.sum(X2**2, axis=1) - 2 * X1 @ X2.T
        return np.exp(-self.gamma * sq_dists)

Step 2: The fit() Method

This is where the core logic of the dual formulation resides. We'll optimize the alphas using gradient ascent.

1. Initialization: We initialize an alpha vector with one entry per training sample.
2. Kernel Matrix: We pre-compute a matrix y_iy_jk_ij = np.outer(y, y) * self.kernel(X, X). This contains the term $y_i{y}_{j}K({\mathbf{x}}_i,{\mathbf{x}}_j)$ for all pairs and is constant, so we only compute it once.
3. Training Loop:
   • Gradient: We compute the gradient for all alphas in a vectorized way: gradient = 1 - y_iy_jk_ij.dot(self.alphas).
   • Update: We update alphas using gradient ascent: self.alphas += lr * gradient.
   • Clipping: We enforce the constraint $0\le {\alpha }_i\le C$ by clipping the values.
4. Bias Calculation: After the loop, we find the support vectors (where 0 < alpha < C) and use them to compute a stable value for the bias b.

# Add this method to your KernelSVM class

def fit(self, X, y, lr=1e-3, epochs=100):
    n_samples, n_features = X.shape
    self.X_train = X
    self.y_train = y
    
    # 1. Initialize alphas
    self.alphas = np.zeros(n_samples)
    self.b = 0
    
    # 2. Pre-compute the kernel-related matrix
    kernel_matrix = self.kernel(X, X)
    y_iy_jk_ij = np.outer(y, y) * kernel_matrix

    # 3. Training loop
    for epoch in range(epochs):
        # Vectorized gradient calculation
        gradient = np.ones(n_samples) - y_iy_jk_ij.dot(self.alphas)
        
        # Update alphas
        self.alphas += lr * gradient
        
        # Clipping to enforce 0 <= alpha <= C
        self.alphas[self.alphas < 0] = 0
        self.alphas[self.alphas > self.C] = self.C

    # 4. Calculate bias b
    # Find support vectors (more robust to use a small tolerance)
    sv_indices = np.where((self.alphas > 1e-5) & (self.alphas < self.C - 1e-5))[0]
    
    if len(sv_indices) > 0:
        # Calculate b for each support vector and average
        # b = y_k - sum_i(alpha_i * y_i * K(x_i, x_k))
        b_k = y[sv_indices] - (self.alphas * y).dot(kernel_matrix[:, sv_indices])
        self.b = np.mean(b_k)
    else:
        # Fallback if no SVs are found in the margin (can happen with some datasets/params)
        print("Warning: Could not find support vectors in the margin to compute bias. Bias remains 0.")
        self.b = 0

Step 3: The predict() Method

Prediction now relies on the learned alphas and the kernel function, not a weight vector w. For a new point $\mathbf{z}$, we compute the kernel between $\mathbf{z}$ and all the training points ${\mathbf{x}}_i$, then weigh the results by ${\alpha }_i{y}_i$.

# Add this method to your KernelSVM class

def predict(self, X):
    # K(all_training_data, new_data)
    kernel_matrix = self.kernel(self.X_train, X)
    
    # decision_function = sum(alpha_i * y_i * K(x_i, z)) + b
    decision_function = (self.alphas * self.y_train).dot(kernel_matrix) + self.b
    
    return np.sign(decision_function)

Seeing it in Action

Now you have a complete KernelSVM! Let's test it on a classic non-linear dataset, like the "moons" dataset from scikit-learn.

from sklearn.datasets import make_moons
from sklearn.model_selection import train_test_split
import matplotlib.pyplot as plt

# Generate data
X, y = make_moons(n_samples=200, noise=0.1, random_state=42)
y[y == 0] = -1 # SVMs typically use labels {-1, 1}

X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=42)

# --- Train with RBF Kernel ---
svm_rbf = KernelSVM(kernel='rbf', C=1.0, gamma=10)
svm_rbf.fit(X_train, y_train, lr=1e-2, epochs=200)

# --- Plot decision boundary (helper function) ---
def plot_decision_boundary(svm, X, y):
    x_min, x_max = X[:, 0].min() - 0.1, X[:, 0].max() + 0.1
    y_min, y_max = X[:, 1].min() - 0.1, X[:, 1].max() + 0.1
    xx, yy = np.meshgrid(np.arange(x_min, x_max, 0.01),
                       np.arange(y_min, y_max, 0.01))
    
    grid_points = np.c_[xx.ravel(), yy.ravel()]
    Z = svm.predict(grid_points)
    Z = Z.reshape(xx.shape)
    
    plt.contourf(xx, yy, Z, alpha=0.4)
    plt.scatter(X[:, 0], X[:, 1], c=y, s=20, edgecolors='k')
    plt.title(f"SVM with {svm.kernel_name} kernel")
    plt.show()

plot_decision_boundary(svm_rbf, X_train, y_train)

If you run this code, you should see a beautiful, non-linear boundary that successfully separates the two moon shapes! Try experimenting with a polynomial kernel as well.

Conclusion

Congratulations! You have successfully upgraded your SVM to handle non-linear data by implementing the dual formulation with kernels. This is a significant milestone and represents one of the most elegant ideas in classical machine learning.

Key Takeaways:

• Dual Formulation: To use kernels, we switched from the primal problem (optimizing $\mathbf{w}$) to the dual problem (optimizing $\alpha$).
• Optimization: We used gradient ascent to find the optimal alpha values, which represent the importance of each data point.
• Sparsity: Most alphas become zero, meaning the model relies only on a subset of data (the support vectors) for prediction, making it efficient.
• Integration: You successfully modified your fit and predict methods to use a generic kernel function, making your SVM versatile.

Preview of the Next Lesson:
We are now moving to a completely different family of algorithms. So far, our models (Linear/Logistic Regression, SVMs) have been geometric in nature, finding lines or planes to separate data. In the next module, we will begin our study of Decision Trees. These models make predictions by learning a hierarchy of simple if/else rules, creating a tree-like structure that is highly interpretable. Our first step will be to implement the logic for finding the best rule, or "split," in a dataset.