Support Vector Machine

1import numpy as np
2
3class LinearSVM:
4    """Hard/soft-margin linear SVM using simplified SMO (Sequential Minimal Optimization)."""
5
6    def __init__(self, C=1.0, tol=1e-3, max_passes=100):
7        self.C = C
8        self.tol = tol
9        self.max_passes = max_passes
10        self.alphas = None
11        self.b = None
12        self.X = None
13        self.y = None
14
15    def _kernel(self, x1, x2):
16        return x1 @ x2  # Linear kernel: K(x,z) = x·z
17
18    def fit(self, X, y):
19        """Simplified SMO for linear kernel SVM."""
20        n, d = X.shape
21        self.X = X.copy()
22        self.y = y.astype(float)  # y ∈ {-1, +1}
23
24        alphas = np.zeros(n)
25        b = 0.0
26        passes = 0
27
28        # Precompute kernel matrix K[i,j] = X[i] · X[j]
29        K = X @ X.T  # (n, n)
30
31        while passes < self.max_passes:
32            num_changed = 0
33
34            for i in range(n):
35                # Compute decision function for sample i
36                f_i = float(np.sum(alphas * self.y * K[:, i])) + b
37                E_i = f_i - self.y[i]
38
39                # Check KKT violation
40                if (self.y[i] * E_i < -self.tol and alphas[i] < self.C) or                    (self.y[i] * E_i > self.tol and alphas[i] > 0):
41
42                    # Select j ≠ i randomly (simplified heuristic)
43                    j = np.random.choice([k for k in range(n) if k != i])
44                    f_j = float(np.sum(alphas * self.y * K[:, j])) + b
45                    E_j = f_j - self.y[j]
46
47                    alpha_i_old, alpha_j_old = alphas[i], alphas[j]
48
49                    # Compute bounds L and H for alpha_j
50                    if self.y[i] == self.y[j]:
51                        L = max(0, alphas[i] + alphas[j] - self.C)
52                        H = min(self.C, alphas[i] + alphas[j])
53                    else:
54                        L = max(0, alphas[j] - alphas[i])
55                        H = min(self.C, self.C + alphas[j] - alphas[i])
56
57                    if L >= H:
58                        continue
59
60                    # Compute eta = 2·K(i,j) - K(i,i) - K(j,j)
61                    eta = 2 * K[i, j] - K[i, i] - K[j, j]
62                    if eta >= 0:
63                        continue
64
65                    # Update alpha_j
66                    alphas[j] -= self.y[j] * (E_i - E_j) / eta
67                    alphas[j] = np.clip(alphas[j], L, H)
68
69                    if abs(alphas[j] - alpha_j_old) < 1e-5:
70                        continue
71
72                    # Update alpha_i (maintains Σ αᵢyᵢ = 0)
73                    alphas[i] += self.y[i] * self.y[j] * (alpha_j_old - alphas[j])
74
75                    # Update bias b
76                    b1 = b - E_i - self.y[i]*(alphas[i]-alpha_i_old)*K[i,i]                                    - self.y[j]*(alphas[j]-alpha_j_old)*K[i,j]
77                    b2 = b - E_j - self.y[i]*(alphas[i]-alpha_i_old)*K[i,j]                                    - self.y[j]*(alphas[j]-alpha_j_old)*K[j,j]
78
79                    if 0 < alphas[i] < self.C:
80                        b = b1
81                    elif 0 < alphas[j] < self.C:
82                        b = b2
83                    else:
84                        b = (b1 + b2) / 2.0
85
86                    num_changed += 1
87
88            passes = passes + 1 if num_changed == 0 else 0
89
90        self.alphas = alphas
91        self.b = b
92
93        # Recover w for linear kernel: w = Σ αᵢ yᵢ xᵢ
94        self.w = (alphas * self.y) @ X
95        sv_mask = alphas > 1e-5
96        print(f"Support vectors: {sv_mask.sum()} / {n} ({100*sv_mask.mean():.1f}%)")
97        return self
98
99    def decision_function(self, X):
100        return X @ self.w + self.b
101
102    def predict(self, X):
103        return np.sign(self.decision_function(X)).astype(int)
104
105    def score(self, X, y):
106        return np.mean(self.predict(X) == y)
107
108
109# ── Demo ──────────────────────────────────────────────────────────────────────
110from sklearn.datasets import make_classification
111from sklearn.model_selection import train_test_split
112from sklearn.preprocessing import StandardScaler
113
114X, y = make_classification(n_samples=300, n_features=2, n_redundant=0,
115                            random_state=42, n_clusters_per_class=1)
116y = np.where(y == 0, -1, 1)  # Convert to {-1, +1}
117
118X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=42)
119
120scaler = StandardScaler()
121X_train_s = scaler.fit_transform(X_train)
122X_test_s  = scaler.transform(X_test)
123
124svm = LinearSVM(C=1.0, max_passes=50)
125svm.fit(X_train_s, y_train)
126print(f"Train accuracy: {svm.score(X_train_s, y_train):.4f}")
127print(f"Test accuracy:  {svm.score(X_test_s, y_test):.4f}")
128print(f"Weight vector:  {svm.w.round(4)}")
129print(f"Bias:           {svm.b:.4f}")