09 Support Vector Machines — First Principles¶

Goal: Build a linear SVM from scratch using sub-gradient descent on the hinge loss, understand margin geometry, identify support vectors, and compare with sklearn.

Prerequisites: Linear Algebra, Calculus & Optimization, 04 Logistic Regression

Theory: theory.md — read §1–§6 before this notebook.

In [1]:

import random

import matplotlib.pyplot as plt
import numpy as np

%load_ext autoreload
%autoreload 2

SEED = 42
random.seed(SEED)
rng = np.random.default_rng(SEED)

import random import matplotlib.pyplot as plt import numpy as np %load_ext autoreload %autoreload 2 SEED = 42 random.seed(SEED) rng = np.random.default_rng(SEED)

1. Problem Setup — WHY¶

Many hyperplanes can separate two linearly separable classes. SVM picks the one with the maximum margin — the widest gap between classes. This choice is geometrically motivated: a wider margin means the classifier is more robust to perturbations.

In [2]:

# Generate linearly separable 2D data
n_per_class = 50
X_pos = rng.normal(loc=[2.0, 2.0], scale=0.6, size=(n_per_class, 2))
X_neg = rng.normal(loc=[-2.0, -2.0], scale=0.6, size=(n_per_class, 2))
X_demo = np.vstack([X_pos, X_neg])
y_demo = np.array([1] * n_per_class + [-1] * n_per_class)

# Show the data and two arbitrary separating hyperplanes
fig, axes = plt.subplots(1, 2, figsize=(14, 5))

for ax in axes:
    ax.scatter(X_pos[:, 0], X_pos[:, 1], c='tab:blue', label='y = +1', edgecolors='k', s=40)
    ax.scatter(X_neg[:, 0], X_neg[:, 1], c='tab:red', label='y = −1', edgecolors='k', s=40)
    ax.set(xlabel='$x_1$', ylabel='$x_2$')
    ax.legend()
    ax.grid(True, alpha=0.3)

# Left: a bad (narrow margin) separator
x_line = np.linspace(-4, 4, 100)
axes[0].plot(x_line, -x_line + 0.8, 'g--', linewidth=2, label='Narrow margin')
axes[0].set_title('Narrow margin hyperplane')
axes[0].legend()

# Right: a better (wider margin) separator
axes[1].plot(x_line, -x_line, 'g-', linewidth=2, label='Wide margin')
axes[1].fill_between(x_line, -x_line - 0.7, -x_line + 0.7, alpha=0.15, color='green')
axes[1].set_title('Wide margin hyperplane (SVM goal)')
axes[1].legend()

plt.tight_layout()
plt.show()

# Generate linearly separable 2D data n_per_class = 50 X_pos = rng.normal(loc=[2.0, 2.0], scale=0.6, size=(n_per_class, 2)) X_neg = rng.normal(loc=[-2.0, -2.0], scale=0.6, size=(n_per_class, 2)) X_demo = np.vstack([X_pos, X_neg]) y_demo = np.array([1] * n_per_class + [-1] * n_per_class) # Show the data and two arbitrary separating hyperplanes fig, axes = plt.subplots(1, 2, figsize=(14, 5)) for ax in axes: ax.scatter(X_pos[:, 0], X_pos[:, 1], c='tab:blue', label='y = +1', edgecolors='k', s=40) ax.scatter(X_neg[:, 0], X_neg[:, 1], c='tab:red', label='y = −1', edgecolors='k', s=40) ax.set(xlabel='$x_1$', ylabel='$x_2$') ax.legend() ax.grid(True, alpha=0.3) # Left: a bad (narrow margin) separator x_line = np.linspace(-4, 4, 100) axes[0].plot(x_line, -x_line + 0.8, 'g--', linewidth=2, label='Narrow margin') axes[0].set_title('Narrow margin hyperplane') axes[0].legend() # Right: a better (wider margin) separator axes[1].plot(x_line, -x_line, 'g-', linewidth=2, label='Wide margin') axes[1].fill_between(x_line, -x_line - 0.7, -x_line + 0.7, alpha=0.15, color='green') axes[1].set_title('Wide margin hyperplane (SVM goal)') axes[1].legend() plt.tight_layout() plt.show()

2. Mathematical Core — WHAT¶

Soft-margin SVM objective (hinge loss form)¶

$$\min_{w, b} \quad \frac{1}{2}\|w\|_2^2 + C \cdot \frac{1}{n} \sum_{i=1}^n \max(0, 1 - y_i(w^T x_i + b))$$

• First term $\frac{1}{2}\|w\|^2$: regularizer — maximizes margin ($\text{margin} = 2/\|w\|$).
• Second term: average hinge loss — penalizes points inside the margin or misclassified.
• $C$: trade-off parameter. Large $C$ → narrow margin, few violations. Small $C$ → wide margin.

Subgradient w.r.t. $w$ and $b$¶

For a mini-batch $\mathcal{B}$:

$$g_w = w - \frac{C}{|\mathcal{B}|} \sum_{i \in \mathcal{B}:\, y_i f(x_i) < 1} y_i x_i, \qquad g_b = -\frac{C}{|\mathcal{B}|} \sum_{i \in \mathcal{B}:\, y_i f(x_i) < 1} y_i.$$

3. Solution Method — HOW¶

We use sub-gradient descent (SGD variant) on the unconstrained hinge-loss objective. At each step:

1. Sample a mini-batch.
2. Compute the functional margin $m_i = y_i(w^T x_i + b)$ for each sample.
3. Identify violators: samples with $m_i < 1$.
4. Compute the subgradient and update $w$, $b$.

Margin geometry visualization¶

The margin boundaries are the hyperplanes $w^T x + b = +1$ and $w^T x + b = -1$. The decision boundary is $w^T x + b = 0$. Support vectors lie on the margin boundaries.

In [3]:

def plot_svm_margin(w, b, X, y, title='SVM Decision Boundary and Margin'):
    """Visualize the decision boundary, margin bands, and support vectors."""
    fig, ax = plt.subplots(figsize=(8, 6))
    
    # Scatter points
    mask_pos = y == 1
    mask_neg = y == -1
    ax.scatter(X[mask_pos, 0], X[mask_pos, 1], c='tab:blue', label='y = +1',
               edgecolors='k', s=50, zorder=3)
    ax.scatter(X[mask_neg, 0], X[mask_neg, 1], c='tab:red', label='y = −1',
               edgecolors='k', s=50, zorder=3)
    
    # Decision boundary and margins
    x_min, x_max = X[:, 0].min() - 1, X[:, 0].max() + 1
    x_vals = np.linspace(x_min, x_max, 200)
    
    if np.abs(w[1]) > 1e-10:
        # w[0]*x1 + w[1]*x2 + b = 0  =>  x2 = -(w[0]*x1 + b) / w[1]
        boundary = -(w[0] * x_vals + b) / w[1]
        margin_pos = -(w[0] * x_vals + b - 1) / w[1]
        margin_neg = -(w[0] * x_vals + b + 1) / w[1]
        
        y_min, y_max = X[:, 1].min() - 1, X[:, 1].max() + 1
        ax.plot(x_vals, boundary, 'k-', linewidth=2, label='Decision boundary')
        ax.plot(x_vals, margin_pos, 'k--', linewidth=1, label='Margin boundary')
        ax.plot(x_vals, margin_neg, 'k--', linewidth=1)
        ax.fill_between(x_vals, margin_neg, margin_pos, alpha=0.1, color='green')
        ax.set_ylim(y_min, y_max)
    
    # Highlight support vectors (points with margin <= 1 + epsilon)
    margins = y * (X @ w + b)
    sv_mask = margins <= 1.0 + 1e-2
    ax.scatter(X[sv_mask, 0], X[sv_mask, 1], s=150, facecolors='none',
               edgecolors='gold', linewidths=2, label='Support vectors', zorder=4)
    
    margin_width = 2.0 / np.linalg.norm(w)
    ax.set(title=f'{title}\nMargin width = {margin_width:.3f}', xlabel='$x_1$', ylabel='$x_2$')
    ax.legend(loc='best')
    ax.grid(True, alpha=0.3)
    return fig

def plot_svm_margin(w, b, X, y, title='SVM Decision Boundary and Margin'): """Visualize the decision boundary, margin bands, and support vectors.""" fig, ax = plt.subplots(figsize=(8, 6)) # Scatter points mask_pos = y == 1 mask_neg = y == -1 ax.scatter(X[mask_pos, 0], X[mask_pos, 1], c='tab:blue', label='y = +1', edgecolors='k', s=50, zorder=3) ax.scatter(X[mask_neg, 0], X[mask_neg, 1], c='tab:red', label='y = −1', edgecolors='k', s=50, zorder=3) # Decision boundary and margins x_min, x_max = X[:, 0].min() - 1, X[:, 0].max() + 1 x_vals = np.linspace(x_min, x_max, 200) if np.abs(w[1]) > 1e-10: # w[0]*x1 + w[1]*x2 + b = 0 => x2 = -(w[0]*x1 + b) / w[1] boundary = -(w[0] * x_vals + b) / w[1] margin_pos = -(w[0] * x_vals + b - 1) / w[1] margin_neg = -(w[0] * x_vals + b + 1) / w[1] y_min, y_max = X[:, 1].min() - 1, X[:, 1].max() + 1 ax.plot(x_vals, boundary, 'k-', linewidth=2, label='Decision boundary') ax.plot(x_vals, margin_pos, 'k--', linewidth=1, label='Margin boundary') ax.plot(x_vals, margin_neg, 'k--', linewidth=1) ax.fill_between(x_vals, margin_neg, margin_pos, alpha=0.1, color='green') ax.set_ylim(y_min, y_max) # Highlight support vectors (points with margin <= 1 + epsilon) margins = y * (X @ w + b) sv_mask = margins <= 1.0 + 1e-2 ax.scatter(X[sv_mask, 0], X[sv_mask, 1], s=150, facecolors='none', edgecolors='gold', linewidths=2, label='Support vectors', zorder=4) margin_width = 2.0 / np.linalg.norm(w) ax.set(title=f'{title}\nMargin width = {margin_width:.3f}', xlabel='$x_1$', ylabel='$x_2$') ax.legend(loc='best') ax.grid(True, alpha=0.3) return fig

4. Implementation — BUILD¶

From-scratch Linear SVM via sub-gradient descent¶

In [4]:

class LinearSVCScratch:
    """Linear soft-margin SVM using sub-gradient descent on hinge loss.
    
    Minimizes: (1/2)||w||^2 + C * mean(max(0, 1 - y_i * (w @ x_i + b)))
    """
    
    def __init__(self, C=1.0, lr=0.01, max_iter=1000, batch_size=32, random_state=42):
        self.C = C
        self.lr = lr
        self.max_iter = max_iter
        self.batch_size = batch_size
        self.random_state = random_state
        self.w_ = None
        self.b_ = None
        self.classes_ = None
        self.loss_history_ = []
    
    def _encode_labels(self, y):
        """Map original labels to {-1, +1}."""
        self.classes_ = np.unique(y)
        if len(self.classes_) != 2:
            raise ValueError('LinearSVCScratch supports exactly two classes')
        return np.where(y == self.classes_[0], -1.0, 1.0)
    
    def _hinge_loss(self, X, y_enc):
        """Compute the full objective: 0.5*||w||^2 + C*mean(hinge)."""
        margins = y_enc * (X @ self.w_ + self.b_)
        hinge = np.maximum(0, 1 - margins)
        return 0.5 * np.dot(self.w_, self.w_) + self.C * np.mean(hinge)
    
    def fit(self, X, y):
        """Fit the linear SVM using mini-batch sub-gradient descent."""
        X = np.asarray(X, dtype=float)
        y_enc = self._encode_labels(y)
        n, p = X.shape
        
        # Initialize weights to zero
        self.w_ = np.zeros(p)
        self.b_ = 0.0
        self.loss_history_ = []
        
        rng_fit = np.random.default_rng(self.random_state)
        batch_size = min(self.batch_size, n)
        
        for epoch in range(self.max_iter):
            indices = rng_fit.permutation(n)
            for start in range(0, n, batch_size):
                batch_idx = indices[start:start + batch_size]
                X_batch = X[batch_idx]
                y_batch = y_enc[batch_idx]
                
                # Functional margins
                margins = y_batch * (X_batch @ self.w_ + self.b_)
                
                # Sub-gradient
                violations = margins < 1.0
                
                # Gradient w.r.t. w: w - C * mean(y_i * x_i for violators)
                grad_w = self.w_.copy()
                if np.any(violations):
                    grad_w -= self.C * np.mean(
                        y_batch[violations, None] * X_batch[violations], axis=0
                    )
                
                # Gradient w.r.t. b: -C * mean(y_i for violators)
                grad_b = 0.0
                if np.any(violations):
                    grad_b = -self.C * np.mean(y_batch[violations])
                
                self.w_ -= self.lr * grad_w
                self.b_ -= self.lr * grad_b
            
            self.loss_history_.append(self._hinge_loss(X, y_enc))
        
        return self
    
    def decision_function(self, X):
        """Compute w^T x + b for each sample."""
        return np.asarray(X, dtype=float) @ self.w_ + self.b_
    
    def predict(self, X):
        """Predict original class labels."""
        raw = (self.decision_function(X) >= 0.0).astype(int)
        return self.classes_[raw]
    
    def score(self, X, y):
        """Classification accuracy."""
        return np.mean(self.predict(X) == y)
    
    def support_vectors(self, X, y):
        """Identify support vectors: points with functional margin <= 1."""
        y_enc = np.where(y == self.classes_[0], -1.0, 1.0)
        margins = y_enc * (X @ self.w_ + self.b_)
        sv_mask = margins <= 1.0 + 1e-3  # small tolerance
        return X[sv_mask], np.where(sv_mask)[0]

class LinearSVCScratch: """Linear soft-margin SVM using sub-gradient descent on hinge loss. Minimizes: (1/2)||w||^2 + C * mean(max(0, 1 - y_i * (w @ x_i + b))) """ def __init__(self, C=1.0, lr=0.01, max_iter=1000, batch_size=32, random_state=42): self.C = C self.lr = lr self.max_iter = max_iter self.batch_size = batch_size self.random_state = random_state self.w_ = None self.b_ = None self.classes_ = None self.loss_history_ = [] def _encode_labels(self, y): """Map original labels to {-1, +1}.""" self.classes_ = np.unique(y) if len(self.classes_) != 2: raise ValueError('LinearSVCScratch supports exactly two classes') return np.where(y == self.classes_[0], -1.0, 1.0) def _hinge_loss(self, X, y_enc): """Compute the full objective: 0.5*||w||^2 + C*mean(hinge).""" margins = y_enc * (X @ self.w_ + self.b_) hinge = np.maximum(0, 1 - margins) return 0.5 * np.dot(self.w_, self.w_) + self.C * np.mean(hinge) def fit(self, X, y): """Fit the linear SVM using mini-batch sub-gradient descent.""" X = np.asarray(X, dtype=float) y_enc = self._encode_labels(y) n, p = X.shape # Initialize weights to zero self.w_ = np.zeros(p) self.b_ = 0.0 self.loss_history_ = [] rng_fit = np.random.default_rng(self.random_state) batch_size = min(self.batch_size, n) for epoch in range(self.max_iter): indices = rng_fit.permutation(n) for start in range(0, n, batch_size): batch_idx = indices[start:start + batch_size] X_batch = X[batch_idx] y_batch = y_enc[batch_idx] # Functional margins margins = y_batch * (X_batch @ self.w_ + self.b_) # Sub-gradient violations = margins < 1.0 # Gradient w.r.t. w: w - C * mean(y_i * x_i for violators) grad_w = self.w_.copy() if np.any(violations): grad_w -= self.C * np.mean( y_batch[violations, None] * X_batch[violations], axis=0 ) # Gradient w.r.t. b: -C * mean(y_i for violators) grad_b = 0.0 if np.any(violations): grad_b = -self.C * np.mean(y_batch[violations]) self.w_ -= self.lr * grad_w self.b_ -= self.lr * grad_b self.loss_history_.append(self._hinge_loss(X, y_enc)) return self def decision_function(self, X): """Compute w^T x + b for each sample.""" return np.asarray(X, dtype=float) @ self.w_ + self.b_ def predict(self, X): """Predict original class labels.""" raw = (self.decision_function(X) >= 0.0).astype(int) return self.classes_[raw] def score(self, X, y): """Classification accuracy.""" return np.mean(self.predict(X) == y) def support_vectors(self, X, y): """Identify support vectors: points with functional margin <= 1.""" y_enc = np.where(y == self.classes_[0], -1.0, 1.0) margins = y_enc * (X @ self.w_ + self.b_) sv_mask = margins <= 1.0 + 1e-3 # small tolerance return X[sv_mask], np.where(sv_mask)[0]

Train on the demo data¶

In [5]:

# Train our from-scratch SVM
svm_scratch = LinearSVCScratch(C=1.0, lr=0.005, max_iter=500, batch_size=32, random_state=42)
svm_scratch.fit(X_demo, y_demo)

print(f'Weights: w = [{svm_scratch.w_[0]:.4f}, {svm_scratch.w_[1]:.4f}]')
print(f'Bias:    b = {svm_scratch.b_:.4f}')
print(f'Accuracy: {svm_scratch.score(X_demo, y_demo):.4f}')
print(f'Margin width: {2.0 / np.linalg.norm(svm_scratch.w_):.4f}')

# Support vectors
sv_points, sv_idx = svm_scratch.support_vectors(X_demo, y_demo)
print(f'Number of support vectors: {len(sv_idx)} / {len(y_demo)}')

# Train our from-scratch SVM svm_scratch = LinearSVCScratch(C=1.0, lr=0.005, max_iter=500, batch_size=32, random_state=42) svm_scratch.fit(X_demo, y_demo) print(f'Weights: w = [{svm_scratch.w_[0]:.4f}, {svm_scratch.w_[1]:.4f}]') print(f'Bias: b = {svm_scratch.b_:.4f}') print(f'Accuracy: {svm_scratch.score(X_demo, y_demo):.4f}') print(f'Margin width: {2.0 / np.linalg.norm(svm_scratch.w_):.4f}') # Support vectors sv_points, sv_idx = svm_scratch.support_vectors(X_demo, y_demo) print(f'Number of support vectors: {len(sv_idx)} / {len(y_demo)}')

Weights: w = [0.4292, 0.3750]
Bias:    b = 0.1480
Accuracy: 1.0000
Margin width: 3.5091
Number of support vectors: 3 / 100

In [6]:

# Visualize decision boundary and margin
fig = plot_svm_margin(svm_scratch.w_, svm_scratch.b_, X_demo, y_demo,
                      title='From-Scratch Linear SVM')
plt.show()

# Visualize decision boundary and margin fig = plot_svm_margin(svm_scratch.w_, svm_scratch.b_, X_demo, y_demo, title='From-Scratch Linear SVM') plt.show()

In [7]:

# Learning curve
fig, ax = plt.subplots(figsize=(8, 4))
ax.plot(svm_scratch.loss_history_, color='tab:blue')
ax.set(title='SVM Hinge Loss During Training', xlabel='Epoch', ylabel='Objective')
ax.grid(True, alpha=0.3)
plt.show()

# Learning curve fig, ax = plt.subplots(figsize=(8, 4)) ax.plot(svm_scratch.loss_history_, color='tab:blue') ax.set(title='SVM Hinge Loss During Training', xlabel='Epoch', ylabel='Objective') ax.grid(True, alpha=0.3) plt.show()

5. Library Comparison¶

5.1 Compare with src/ library LinearSVC¶

In [8]:

from ml_first_principles.svm_models import LinearSVC as LibLinearSVC

svm_lib = LibLinearSVC(C=1.0, lr=0.005, max_iter=500, batch_size=32, random_state=42)
svm_lib.fit(X_demo, y_demo)

acc_scratch = svm_scratch.score(X_demo, y_demo)
acc_lib = svm_lib.score(X_demo, y_demo)

print('=== Library LinearSVC ===')
print(f'Weights: {svm_lib.coef_}')
print(f'Bias:    {svm_lib.intercept_:.4f}')
print(f'Accuracy: {acc_lib:.4f}')

print(f'\n=== Accuracy comparison ===')
print(f'Scratch:  {acc_scratch:.4f}')
print(f'Library:  {acc_lib:.4f}')

# Both should achieve perfect accuracy on this separable data
assert acc_scratch >= 0.95, f'Scratch accuracy too low: {acc_scratch}'
assert acc_lib >= 0.95, f'Library accuracy too low: {acc_lib}'
print('\nBoth achieve >= 95% accuracy ✓')

from ml_first_principles.svm_models import LinearSVC as LibLinearSVC svm_lib = LibLinearSVC(C=1.0, lr=0.005, max_iter=500, batch_size=32, random_state=42) svm_lib.fit(X_demo, y_demo) acc_scratch = svm_scratch.score(X_demo, y_demo) acc_lib = svm_lib.score(X_demo, y_demo) print('=== Library LinearSVC ===') print(f'Weights: {svm_lib.coef_}') print(f'Bias: {svm_lib.intercept_:.4f}') print(f'Accuracy: {acc_lib:.4f}') print(f'\n=== Accuracy comparison ===') print(f'Scratch: {acc_scratch:.4f}') print(f'Library: {acc_lib:.4f}') # Both should achieve perfect accuracy on this separable data assert acc_scratch >= 0.95, f'Scratch accuracy too low: {acc_scratch}' assert acc_lib >= 0.95, f'Library accuracy too low: {acc_lib}' print('\nBoth achieve >= 95% accuracy ✓')

=== Library LinearSVC ===
Weights: [0.28196178 0.30920986]
Bias:    -0.0277
Accuracy: 1.0000

=== Accuracy comparison ===
Scratch:  1.0000
Library:  1.0000

Both achieve >= 95% accuracy ✓

5.2 Compare with sklearn¶

In [9]:

from sklearn.svm import LinearSVC as SkLinearSVC, SVC as SkSVC

# sklearn LinearSVC (hinge loss + L2 regularization)
sk_linear = SkLinearSVC(C=1.0, max_iter=5000, random_state=42)
sk_linear.fit(X_demo, y_demo)
acc_sk_linear = sk_linear.score(X_demo, y_demo)

# sklearn SVC with linear kernel (uses the dual / libsvm)
sk_svc = SkSVC(kernel='linear', C=1.0, random_state=42)
sk_svc.fit(X_demo, y_demo)
acc_sk_svc = sk_svc.score(X_demo, y_demo)

print('=== sklearn Comparison ===')
print(f'sklearn LinearSVC accuracy: {acc_sk_linear:.4f}')
print(f'sklearn SVC(linear) accuracy: {acc_sk_svc:.4f}')
print(f'Our scratch accuracy:        {acc_scratch:.4f}')

# sklearn SVC exposes support vectors
print(f'\nsklearn SVC support vectors: {sk_svc.n_support_} (per class)')
print(f'Total sklearn SVs: {len(sk_svc.support_vectors_)}')
print(f'Our scratch SVs:   {len(sv_idx)}')

from sklearn.svm import LinearSVC as SkLinearSVC, SVC as SkSVC # sklearn LinearSVC (hinge loss + L2 regularization) sk_linear = SkLinearSVC(C=1.0, max_iter=5000, random_state=42) sk_linear.fit(X_demo, y_demo) acc_sk_linear = sk_linear.score(X_demo, y_demo) # sklearn SVC with linear kernel (uses the dual / libsvm) sk_svc = SkSVC(kernel='linear', C=1.0, random_state=42) sk_svc.fit(X_demo, y_demo) acc_sk_svc = sk_svc.score(X_demo, y_demo) print('=== sklearn Comparison ===') print(f'sklearn LinearSVC accuracy: {acc_sk_linear:.4f}') print(f'sklearn SVC(linear) accuracy: {acc_sk_svc:.4f}') print(f'Our scratch accuracy: {acc_scratch:.4f}') # sklearn SVC exposes support vectors print(f'\nsklearn SVC support vectors: {sk_svc.n_support_} (per class)') print(f'Total sklearn SVs: {len(sk_svc.support_vectors_)}') print(f'Our scratch SVs: {len(sv_idx)}')

=== sklearn Comparison ===
sklearn LinearSVC accuracy: 1.0000
sklearn SVC(linear) accuracy: 1.0000
Our scratch accuracy:        1.0000

sklearn SVC support vectors: [1 1] (per class)
Total sklearn SVs: 2
Our scratch SVs:   3

In [10]:

# Side-by-side: scratch vs sklearn
fig, axes = plt.subplots(1, 2, figsize=(14, 5))

for ax, (w, b_val, title_str) in zip(axes, [
    (svm_scratch.w_, svm_scratch.b_, 'From-Scratch LinearSVC'),
    (sk_svc.coef_[0], sk_svc.intercept_[0], 'sklearn SVC(linear)'),
]):
    mask_pos = y_demo == 1
    mask_neg = y_demo == -1
    ax.scatter(X_demo[mask_pos, 0], X_demo[mask_pos, 1], c='tab:blue', edgecolors='k', s=40)
    ax.scatter(X_demo[mask_neg, 0], X_demo[mask_neg, 1], c='tab:red', edgecolors='k', s=40)
    
    x_vals = np.linspace(X_demo[:, 0].min() - 1, X_demo[:, 0].max() + 1, 200)
    if np.abs(w[1]) > 1e-10:
        boundary = -(w[0] * x_vals + b_val) / w[1]
        margin_p = -(w[0] * x_vals + b_val - 1) / w[1]
        margin_n = -(w[0] * x_vals + b_val + 1) / w[1]
        ax.plot(x_vals, boundary, 'k-', linewidth=2)
        ax.plot(x_vals, margin_p, 'k--', linewidth=1)
        ax.plot(x_vals, margin_n, 'k--', linewidth=1)
        ax.fill_between(x_vals, margin_n, margin_p, alpha=0.1, color='green')
    
    margin_w = 2.0 / np.linalg.norm(w)
    ax.set(title=f'{title_str}\nMargin = {margin_w:.3f}', xlabel='$x_1$', ylabel='$x_2$')
    ax.set_ylim(X_demo[:, 1].min() - 1, X_demo[:, 1].max() + 1)
    ax.grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

# Side-by-side: scratch vs sklearn fig, axes = plt.subplots(1, 2, figsize=(14, 5)) for ax, (w, b_val, title_str) in zip(axes, [ (svm_scratch.w_, svm_scratch.b_, 'From-Scratch LinearSVC'), (sk_svc.coef_[0], sk_svc.intercept_[0], 'sklearn SVC(linear)'), ]): mask_pos = y_demo == 1 mask_neg = y_demo == -1 ax.scatter(X_demo[mask_pos, 0], X_demo[mask_pos, 1], c='tab:blue', edgecolors='k', s=40) ax.scatter(X_demo[mask_neg, 0], X_demo[mask_neg, 1], c='tab:red', edgecolors='k', s=40) x_vals = np.linspace(X_demo[:, 0].min() - 1, X_demo[:, 0].max() + 1, 200) if np.abs(w[1]) > 1e-10: boundary = -(w[0] * x_vals + b_val) / w[1] margin_p = -(w[0] * x_vals + b_val - 1) / w[1] margin_n = -(w[0] * x_vals + b_val + 1) / w[1] ax.plot(x_vals, boundary, 'k-', linewidth=2) ax.plot(x_vals, margin_p, 'k--', linewidth=1) ax.plot(x_vals, margin_n, 'k--', linewidth=1) ax.fill_between(x_vals, margin_n, margin_p, alpha=0.1, color='green') margin_w = 2.0 / np.linalg.norm(w) ax.set(title=f'{title_str}\nMargin = {margin_w:.3f}', xlabel='$x_1$', ylabel='$x_2$') ax.set_ylim(X_demo[:, 1].min() - 1, X_demo[:, 1].max() + 1) ax.grid(True, alpha=0.3) plt.tight_layout() plt.show()

6. Experiments and Failures — VERIFY¶

6.1 Effect of C on margin width¶

In [11]:

# Generate data with some overlap for soft-margin demo
X_overlap_pos = rng.normal(loc=[1.5, 1.5], scale=1.0, size=(60, 2))
X_overlap_neg = rng.normal(loc=[-1.5, -1.5], scale=1.0, size=(60, 2))
X_overlap = np.vstack([X_overlap_pos, X_overlap_neg])
y_overlap = np.array([1] * 60 + [-1] * 60)

C_values = [0.01, 0.1, 1.0, 10.0]
fig, axes = plt.subplots(1, 4, figsize=(20, 4))

for ax, C_val in zip(axes, C_values):
    svm_c = LinearSVCScratch(C=C_val, lr=0.005, max_iter=500, batch_size=32, random_state=42)
    svm_c.fit(X_overlap, y_overlap)
    
    mask_pos = y_overlap == 1
    mask_neg = y_overlap == -1
    ax.scatter(X_overlap[mask_pos, 0], X_overlap[mask_pos, 1], c='tab:blue', s=20, edgecolors='k', linewidths=0.5)
    ax.scatter(X_overlap[mask_neg, 0], X_overlap[mask_neg, 1], c='tab:red', s=20, edgecolors='k', linewidths=0.5)
    
    x_vals = np.linspace(-5, 5, 200)
    w, b_val = svm_c.w_, svm_c.b_
    if np.abs(w[1]) > 1e-10:
        boundary = -(w[0] * x_vals + b_val) / w[1]
        margin_p = -(w[0] * x_vals + b_val - 1) / w[1]
        margin_n = -(w[0] * x_vals + b_val + 1) / w[1]
        ax.plot(x_vals, boundary, 'k-', linewidth=2)
        ax.plot(x_vals, margin_p, 'k--', linewidth=1)
        ax.plot(x_vals, margin_n, 'k--', linewidth=1)
        ax.fill_between(x_vals, margin_n, margin_p, alpha=0.1, color='green')
    
    _, sv_i = svm_c.support_vectors(X_overlap, y_overlap)
    margin_w = 2.0 / np.linalg.norm(w) if np.linalg.norm(w) > 1e-10 else float('inf')
    ax.set(title=f'C = {C_val}\nMargin = {margin_w:.2f}, SVs = {len(sv_i)}',
           xlabel='$x_1$', ylabel='$x_2$')
    ax.set_xlim(-5, 5)
    ax.set_ylim(-5, 5)
    ax.grid(True, alpha=0.3)

plt.suptitle('Effect of C on Margin Width', fontsize=14, y=1.02)
plt.tight_layout()
plt.show()

# Generate data with some overlap for soft-margin demo X_overlap_pos = rng.normal(loc=[1.5, 1.5], scale=1.0, size=(60, 2)) X_overlap_neg = rng.normal(loc=[-1.5, -1.5], scale=1.0, size=(60, 2)) X_overlap = np.vstack([X_overlap_pos, X_overlap_neg]) y_overlap = np.array([1] * 60 + [-1] * 60) C_values = [0.01, 0.1, 1.0, 10.0] fig, axes = plt.subplots(1, 4, figsize=(20, 4)) for ax, C_val in zip(axes, C_values): svm_c = LinearSVCScratch(C=C_val, lr=0.005, max_iter=500, batch_size=32, random_state=42) svm_c.fit(X_overlap, y_overlap) mask_pos = y_overlap == 1 mask_neg = y_overlap == -1 ax.scatter(X_overlap[mask_pos, 0], X_overlap[mask_pos, 1], c='tab:blue', s=20, edgecolors='k', linewidths=0.5) ax.scatter(X_overlap[mask_neg, 0], X_overlap[mask_neg, 1], c='tab:red', s=20, edgecolors='k', linewidths=0.5) x_vals = np.linspace(-5, 5, 200) w, b_val = svm_c.w_, svm_c.b_ if np.abs(w[1]) > 1e-10: boundary = -(w[0] * x_vals + b_val) / w[1] margin_p = -(w[0] * x_vals + b_val - 1) / w[1] margin_n = -(w[0] * x_vals + b_val + 1) / w[1] ax.plot(x_vals, boundary, 'k-', linewidth=2) ax.plot(x_vals, margin_p, 'k--', linewidth=1) ax.plot(x_vals, margin_n, 'k--', linewidth=1) ax.fill_between(x_vals, margin_n, margin_p, alpha=0.1, color='green') _, sv_i = svm_c.support_vectors(X_overlap, y_overlap) margin_w = 2.0 / np.linalg.norm(w) if np.linalg.norm(w) > 1e-10 else float('inf') ax.set(title=f'C = {C_val}\nMargin = {margin_w:.2f}, SVs = {len(sv_i)}', xlabel='$x_1$', ylabel='$x_2$') ax.set_xlim(-5, 5) ax.set_ylim(-5, 5) ax.grid(True, alpha=0.3) plt.suptitle('Effect of C on Margin Width', fontsize=14, y=1.02) plt.tight_layout() plt.show()

6.2 Support vector identification¶

Only points on or inside the margin boundary affect the solution. Removing non-support vectors does not change the decision boundary.

In [12]:

# Train on full data
svm_full = LinearSVCScratch(C=1.0, lr=0.005, max_iter=500, batch_size=32, random_state=42)
svm_full.fit(X_demo, y_demo)

# Identify support vectors
y_enc_demo = np.where(y_demo == svm_full.classes_[0], -1.0, 1.0)
margins_full = y_enc_demo * (X_demo @ svm_full.w_ + svm_full.b_)
sv_mask_full = margins_full <= 1.0 + 1e-3
non_sv_mask = ~sv_mask_full

print(f'Total points:         {len(y_demo)}')
print(f'Support vectors:      {np.sum(sv_mask_full)}')
print(f'Non-support vectors:  {np.sum(non_sv_mask)}')
print(f'\nSupport vector indices: {np.where(sv_mask_full)[0]}')

# Train on full data svm_full = LinearSVCScratch(C=1.0, lr=0.005, max_iter=500, batch_size=32, random_state=42) svm_full.fit(X_demo, y_demo) # Identify support vectors y_enc_demo = np.where(y_demo == svm_full.classes_[0], -1.0, 1.0) margins_full = y_enc_demo * (X_demo @ svm_full.w_ + svm_full.b_) sv_mask_full = margins_full <= 1.0 + 1e-3 non_sv_mask = ~sv_mask_full print(f'Total points: {len(y_demo)}') print(f'Support vectors: {np.sum(sv_mask_full)}') print(f'Non-support vectors: {np.sum(non_sv_mask)}') print(f'\nSupport vector indices: {np.where(sv_mask_full)[0]}')

Total points:         100
Support vectors:      3
Non-support vectors:  97

Support vector indices: [ 2 47 69]

6.3 Kernel SVM demo (using sklearn)¶

Linear SVM cannot handle non-linearly separable data. The kernel trick maps data to a higher-dimensional space where a linear separator exists.

In [13]:

# Generate non-linearly separable data: concentric circles
from sklearn.datasets import make_circles

X_circles, y_circles = make_circles(n_samples=200, noise=0.1, factor=0.4, random_state=42)
y_circles_svm = np.where(y_circles == 0, -1, 1)  # convert to {-1, +1}

fig, axes = plt.subplots(1, 3, figsize=(18, 5))

# Linear kernel — fails
sk_linear_circ = SkSVC(kernel='linear', C=1.0, random_state=42)
sk_linear_circ.fit(X_circles, y_circles_svm)

# RBF kernel — succeeds
sk_rbf = SkSVC(kernel='rbf', C=10.0, gamma=2.0, random_state=42)
sk_rbf.fit(X_circles, y_circles_svm)

# Polynomial kernel
sk_poly = SkSVC(kernel='poly', degree=3, C=10.0, random_state=42)
sk_poly.fit(X_circles, y_circles_svm)

for ax, (model, title_str) in zip(axes, [
    (sk_linear_circ, f'Linear kernel (acc={sk_linear_circ.score(X_circles, y_circles_svm):.2f})'),
    (sk_rbf, f'RBF kernel (acc={sk_rbf.score(X_circles, y_circles_svm):.2f})'),
    (sk_poly, f'Poly kernel d=3 (acc={sk_poly.score(X_circles, y_circles_svm):.2f})'),
]):
    # Create mesh for contour
    xx, yy = np.meshgrid(
        np.linspace(X_circles[:, 0].min() - 0.5, X_circles[:, 0].max() + 0.5, 200),
        np.linspace(X_circles[:, 1].min() - 0.5, X_circles[:, 1].max() + 0.5, 200),
    )
    Z = model.predict(np.column_stack([xx.ravel(), yy.ravel()])).reshape(xx.shape)
    
    ax.contourf(xx, yy, Z, alpha=0.3, cmap='coolwarm')
    mask_p = y_circles_svm == 1
    mask_n = y_circles_svm == -1
    ax.scatter(X_circles[mask_p, 0], X_circles[mask_p, 1], c='tab:blue', edgecolors='k', s=30)
    ax.scatter(X_circles[mask_n, 0], X_circles[mask_n, 1], c='tab:red', edgecolors='k', s=30)
    ax.set(title=title_str, xlabel='$x_1$', ylabel='$x_2$')
    ax.grid(True, alpha=0.3)

plt.suptitle('Kernel SVM on Non-Linear Data (Concentric Circles)', fontsize=14, y=1.02)
plt.tight_layout()
plt.show()

# Generate non-linearly separable data: concentric circles from sklearn.datasets import make_circles X_circles, y_circles = make_circles(n_samples=200, noise=0.1, factor=0.4, random_state=42) y_circles_svm = np.where(y_circles == 0, -1, 1) # convert to {-1, +1} fig, axes = plt.subplots(1, 3, figsize=(18, 5)) # Linear kernel — fails sk_linear_circ = SkSVC(kernel='linear', C=1.0, random_state=42) sk_linear_circ.fit(X_circles, y_circles_svm) # RBF kernel — succeeds sk_rbf = SkSVC(kernel='rbf', C=10.0, gamma=2.0, random_state=42) sk_rbf.fit(X_circles, y_circles_svm) # Polynomial kernel sk_poly = SkSVC(kernel='poly', degree=3, C=10.0, random_state=42) sk_poly.fit(X_circles, y_circles_svm) for ax, (model, title_str) in zip(axes, [ (sk_linear_circ, f'Linear kernel (acc={sk_linear_circ.score(X_circles, y_circles_svm):.2f})'), (sk_rbf, f'RBF kernel (acc={sk_rbf.score(X_circles, y_circles_svm):.2f})'), (sk_poly, f'Poly kernel d=3 (acc={sk_poly.score(X_circles, y_circles_svm):.2f})'), ]): # Create mesh for contour xx, yy = np.meshgrid( np.linspace(X_circles[:, 0].min() - 0.5, X_circles[:, 0].max() + 0.5, 200), np.linspace(X_circles[:, 1].min() - 0.5, X_circles[:, 1].max() + 0.5, 200), ) Z = model.predict(np.column_stack([xx.ravel(), yy.ravel()])).reshape(xx.shape) ax.contourf(xx, yy, Z, alpha=0.3, cmap='coolwarm') mask_p = y_circles_svm == 1 mask_n = y_circles_svm == -1 ax.scatter(X_circles[mask_p, 0], X_circles[mask_p, 1], c='tab:blue', edgecolors='k', s=30) ax.scatter(X_circles[mask_n, 0], X_circles[mask_n, 1], c='tab:red', edgecolors='k', s=30) ax.set(title=title_str, xlabel='$x_1$', ylabel='$x_2$') ax.grid(True, alpha=0.3) plt.suptitle('Kernel SVM on Non-Linear Data (Concentric Circles)', fontsize=14, y=1.02) plt.tight_layout() plt.show()

6.4 Failure case: unscaled features¶

In [14]:

from ml_first_principles.data_utils import standardize

# Create data with very different feature scales
X_scaled = X_demo.copy()
X_unscaled = X_demo.copy()
X_unscaled[:, 1] *= 100  # Feature 2 is 100x larger

# Train on unscaled
svm_unscaled = LinearSVCScratch(C=1.0, lr=0.0001, max_iter=500, batch_size=32, random_state=42)
svm_unscaled.fit(X_unscaled, y_demo)
acc_unscaled = svm_unscaled.score(X_unscaled, y_demo)

# Train on standardized
X_std, _, _ = standardize(X_unscaled)
svm_scaled = LinearSVCScratch(C=1.0, lr=0.005, max_iter=500, batch_size=32, random_state=42)
svm_scaled.fit(X_std, y_demo)
acc_scaled = svm_scaled.score(X_std, y_demo)

print('=== Feature Scaling Impact ===')
print(f'Unscaled features — accuracy: {acc_unscaled:.4f}')
print(f'Standardized features — accuracy: {acc_scaled:.4f}')
print(f'\nLesson: Always standardize features before SVM training.')

from ml_first_principles.data_utils import standardize # Create data with very different feature scales X_scaled = X_demo.copy() X_unscaled = X_demo.copy() X_unscaled[:, 1] *= 100 # Feature 2 is 100x larger # Train on unscaled svm_unscaled = LinearSVCScratch(C=1.0, lr=0.0001, max_iter=500, batch_size=32, random_state=42) svm_unscaled.fit(X_unscaled, y_demo) acc_unscaled = svm_unscaled.score(X_unscaled, y_demo) # Train on standardized X_std, _, _ = standardize(X_unscaled) svm_scaled = LinearSVCScratch(C=1.0, lr=0.005, max_iter=500, batch_size=32, random_state=42) svm_scaled.fit(X_std, y_demo) acc_scaled = svm_scaled.score(X_std, y_demo) print('=== Feature Scaling Impact ===') print(f'Unscaled features — accuracy: {acc_unscaled:.4f}') print(f'Standardized features — accuracy: {acc_scaled:.4f}') print(f'\nLesson: Always standardize features before SVM training.')

=== Feature Scaling Impact ===
Unscaled features — accuracy: 1.0000
Standardized features — accuracy: 1.0000

Lesson: Always standardize features before SVM training.

6.5 Failure case: wrong kernel choice¶

In [15]:

# Linear kernel on circle data (wrong choice)
sk_wrong = SkSVC(kernel='linear', C=10.0, random_state=42)
sk_wrong.fit(X_circles, y_circles_svm)
acc_wrong = sk_wrong.score(X_circles, y_circles_svm)

# RBF kernel on circle data (right choice)
sk_right = SkSVC(kernel='rbf', C=10.0, gamma=2.0, random_state=42)
sk_right.fit(X_circles, y_circles_svm)
acc_right = sk_right.score(X_circles, y_circles_svm)

print('=== Kernel Choice Impact (Concentric Circles) ===')
print(f'Linear kernel accuracy: {acc_wrong:.4f}')
print(f'RBF kernel accuracy:    {acc_right:.4f}')
print(f'\nLinear kernel fails because the data is not linearly separable.')
print(f'RBF maps to infinite-dimensional space where a linear separator exists.')

# Linear kernel on circle data (wrong choice) sk_wrong = SkSVC(kernel='linear', C=10.0, random_state=42) sk_wrong.fit(X_circles, y_circles_svm) acc_wrong = sk_wrong.score(X_circles, y_circles_svm) # RBF kernel on circle data (right choice) sk_right = SkSVC(kernel='rbf', C=10.0, gamma=2.0, random_state=42) sk_right.fit(X_circles, y_circles_svm) acc_right = sk_right.score(X_circles, y_circles_svm) print('=== Kernel Choice Impact (Concentric Circles) ===') print(f'Linear kernel accuracy: {acc_wrong:.4f}') print(f'RBF kernel accuracy: {acc_right:.4f}') print(f'\nLinear kernel fails because the data is not linearly separable.') print(f'RBF maps to infinite-dimensional space where a linear separator exists.')

=== Kernel Choice Impact (Concentric Circles) ===
Linear kernel accuracy: 0.5950
RBF kernel accuracy:    0.9950

Linear kernel fails because the data is not linearly separable.
RBF maps to infinite-dimensional space where a linear separator exists.

7. Connections¶

Concept	SVM view	Related topic
Decision boundary	Maximum-margin hyperplane $w^T x + b = 0$	Logistic Regression uses same boundary, different loss
Regularization	$\frac{1}{2}\Vert w\Vert^2$ = L2 penalty	Regularization
Loss function	Hinge $\max(0, 1-m)$ vs log-loss $\log(1+e^{-m})$	Loss Functions
Geometry	Margin = $2/\Vert w\Vert $, support vectors define boundary	Geometry of ML
Non-linearity	Kernel trick	Neural Networks (learned features)

Takeaway¶

SVM finds the maximum-margin linear separator by minimizing a combination of the L2 norm of the weight vector and the hinge loss. The solution depends only on a subset of training points — the support vectors — making SVM memory-efficient at prediction time. The kernel trick extends linear SVM to non-linear boundaries by implicitly mapping data to higher-dimensional spaces. Key practical considerations: always standardize features, tune $C$ via cross-validation, and choose the kernel based on data geometry.