{
"cells": [
{
"cell_type": "markdown",
"id": "004ac186",
"metadata": {},
"source": [
"# CS6140 Machine Learning — Fall 2026\n",
"# HW3 Starter — Kernels, SVM, Kernel Ridge Regression, PCA, Kernel PCA\n",
"\n",
"This notebook implements `HW3_26F.html`. Where a library offers the same functionality you'll run\n",
"it first as a sanity-check baseline, then implement the same thing from scratch and evaluate it\n",
"exactly the same way, so your two numbers should land right next to each other once everything is\n",
"filled in.\n",
"\n",
"**Problems 1-5 (required):**\n",
"1. SVM via a library, several kernels: Spambase (binary), Digits/MNIST (multiclass)\n",
"2. SVM from scratch: **Simplified SMO** (CS229 handout, `materials/smo.pdf`), validated on a\n",
" linear kernel then run on RBF (the required non-linear case) — TwoSpirals\n",
"3. Kernel ridge regression from scratch (linear / RBF / polynomial kernels): WAVE, CRESCENT\n",
"4. PCA from scratch, feature-corruption narrative: polluted Spambase + logistic regression\n",
"5. Kernel PCA from scratch + linear classifier: TwoSpirals, ThreeCircles\n",
"\n",
"**How to work through this notebook.** Everything runs except the lines marked\n",
"`### TODO_STUDENT ###` — those have their code replaced with a `...` placeholder and a one-line\n",
"description of what to implement there. Read the surrounding (working) code and the markdown above\n",
"each section for context, then replace each `...` with real code. A cell with a `...` in it will\n",
"raise a `SyntaxError` until you fill it in — that's expected, not a bug. A few TODOs are\n",
"mirror-image pairs (Problem 2's L/H bounds and its b1/b2 threshold formulas): one side of each\n",
"pair is given to you as a fully worked example, and the comment on the other side tells you it\n",
"mirrors the given one — use it as your template. Work top to bottom: Problem 2's\n",
"`linear_gram`/`rbf_gram`/`simplified_smo`/`predict_smo` are reused in both Step 1 (linear) and\n",
"Step 2 (RBF); Problem 4's `fit_eval_logreg` helper is reused by Problem 5. Every hyperparameter or\n",
"dimension choice (Problem 1's kernel grid, Problem 3's kernel/regularization, Problem 5's `T`) is\n",
"selected on a **validation** split and confirmed once on the held-out **test** split — never pick\n",
"a configuration by looking at test performance. Once a \"library baseline\" cell's TODOs are filled\n",
"in and run, the printed number is a target: your matching \"from scratch\" implementation just below\n",
"it should land on (or very near) the same value.\n",
"\n",
"Datasets are read from the shared course `data/` folder (paths relative to\n",
"`6_SVM_kernels/code_HW6/`, i.e. two levels up).\n",
"\n",
"**Requirements.**\n",
"1. Fill in every `TODO_STUDENT` blank according to the algorithm steps covered in lecture and the\n",
" linked materials, and make sure the notebook runs top to bottom without errors.\n",
"2. Understand the *whole* notebook — including the provided/given code, not just the blanks you\n",
" filled in — well enough to explain any part of it during office hours.\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "fed55f52",
"metadata": {},
"outputs": [],
"source": [
"import numpy as np\n",
"import matplotlib.pyplot as plt\n",
"from sklearn.preprocessing import StandardScaler\n",
"from sklearn.model_selection import train_test_split\n",
"from sklearn.metrics import accuracy_score, mean_squared_error\n",
"from sklearn.svm import SVC\n",
"from sklearn.multiclass import OneVsRestClassifier\n",
"from sklearn.decomposition import PCA, KernelPCA\n",
"from sklearn.kernel_ridge import KernelRidge\n",
"from sklearn.linear_model import LogisticRegression\n",
"\n",
"DATA = \"../../data\"\n",
"np.random.seed(42)\n"
]
},
{
"cell_type": "markdown",
"id": "739541dc",
"metadata": {},
"source": [
"## Data loading"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "84bc0838",
"metadata": {},
"outputs": [],
"source": [
"# ---- Spambase (binary classification): one file, split into train/val/test ----\n",
"# Validation is for model/hyperparameter selection (Problem 1's kernel grid); test is reserved\n",
"# for a single final reported number, never used to pick a configuration.\n",
"spam_raw = np.loadtxt(f\"{DATA}/spambase/spambase.data\", delimiter=\",\")\n",
"X_spam, y_spam = spam_raw[:, :-1], spam_raw[:, -1] # y in {0,1}\n",
"\n",
"X_train_spam, X_temp_spam, y_train_spam, y_temp_spam = train_test_split(\n",
" X_spam, y_spam, test_size=0.4, random_state=42, stratify=y_spam)\n",
"X_val_spam, X_test_spam, y_val_spam, y_test_spam = train_test_split(\n",
" X_temp_spam, y_temp_spam, test_size=0.5, random_state=42, stratify=y_temp_spam)\n",
"\n",
"spam_scaler = StandardScaler().fit(X_train_spam)\n",
"X_train_spam = spam_scaler.transform(X_train_spam)\n",
"X_val_spam = spam_scaler.transform(X_val_spam)\n",
"X_test_spam = spam_scaler.transform(X_test_spam)\n",
"\n",
"print(\"Spambase train:\", X_train_spam.shape, \" val:\", X_val_spam.shape, \" test:\", X_test_spam.shape)\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "010b0140",
"metadata": {},
"outputs": [],
"source": [
"# ---- MNIST/Digits (Haar features): pre-extracted 200-dim features, subsampled like HW1 ----\n",
"X_train_digits_full = np.loadtxt(f\"{DATA}/mnist_haar_bingyu/training_image.txt\", delimiter=\",\")\n",
"y_train_digits_full = np.loadtxt(f\"{DATA}/mnist_haar_bingyu/training_label.txt\", delimiter=\",\")\n",
"X_test_digits_full = np.loadtxt(f\"{DATA}/mnist_haar_bingyu/testing_image.txt\", delimiter=\",\")\n",
"y_test_digits_full = np.loadtxt(f\"{DATA}/mnist_haar_bingyu/testing_label.txt\", delimiter=\",\")\n",
"\n",
"n_train_digits, n_test_digits = 3000, 1000\n",
"rng = np.random.RandomState(0)\n",
"train_idx = rng.choice(len(X_train_digits_full), size=n_train_digits, replace=False)\n",
"test_idx = rng.choice(len(X_test_digits_full), size=n_test_digits, replace=False)\n",
"\n",
"X_train_digits = X_train_digits_full[train_idx]\n",
"y_train_digits = y_train_digits_full[train_idx]\n",
"X_test_digits = X_test_digits_full[test_idx]\n",
"y_test_digits = y_test_digits_full[test_idx]\n",
"\n",
"digits_scaler = StandardScaler().fit(X_train_digits)\n",
"X_train_digits = digits_scaler.transform(X_train_digits)\n",
"X_test_digits = digits_scaler.transform(X_test_digits)\n",
"\n",
"print(\"Digits (Problem 1b) train:\", X_train_digits.shape, \" test:\", X_test_digits.shape)\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "31170ffb",
"metadata": {},
"outputs": [],
"source": [
"# ---- WAVE and CRESCENT: 1-D input / 1-D output regression toy sets, one file each ----\n",
"# Same train/val/test convention as Spambase: val picks Problem 3's kernel/hyperparameters,\n",
"# test is reported once for the chosen configuration.\n",
"def load_1d_regression(name, val_size=0.2, test_size=0.2, seed=42):\n",
" x = np.loadtxt(f\"{DATA}/{name}X.txt\").reshape(-1, 1)\n",
" y = np.loadtxt(f\"{DATA}/{name}Y.txt\")\n",
" x_train, x_temp, y_train, y_temp = train_test_split(\n",
" x, y, test_size=val_size + test_size, random_state=seed)\n",
" x_val, x_test, y_val, y_test = train_test_split(\n",
" x_temp, y_temp, test_size=test_size / (val_size + test_size), random_state=seed)\n",
" return x_train, x_val, x_test, y_train, y_val, y_test\n",
"\n",
"X_train_wave, X_val_wave, X_test_wave, y_train_wave, y_val_wave, y_test_wave = load_1d_regression(\"wave\")\n",
"X_train_crescent, X_val_crescent, X_test_crescent, y_train_crescent, y_val_crescent, y_test_crescent = load_1d_regression(\"crescent\")\n",
"\n",
"print(\"WAVE train:\", X_train_wave.shape, \" val:\", X_val_wave.shape, \" test:\", X_test_wave.shape)\n",
"print(\"CRESCENT train:\", X_train_crescent.shape, \" val:\", X_val_crescent.shape, \" test:\", X_test_crescent.shape)\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "a332ff2b",
"metadata": {},
"outputs": [],
"source": [
"# ---- Polluted Spambase (original features + junk/duplicated ones), already train/test split ----\n",
"X_train_pol = np.loadtxt(f\"{DATA}/spam_polluted/train_feature.txt\")\n",
"y_train_pol = np.loadtxt(f\"{DATA}/spam_polluted/train_label.txt\")\n",
"X_test_pol = np.loadtxt(f\"{DATA}/spam_polluted/test_feature.txt\")\n",
"y_test_pol = np.loadtxt(f\"{DATA}/spam_polluted/test_label.txt\")\n",
"\n",
"pol_scaler = StandardScaler().fit(X_train_pol)\n",
"X_train_pol_s = pol_scaler.transform(X_train_pol)\n",
"X_test_pol_s = pol_scaler.transform(X_test_pol)\n",
"\n",
"print(\"Polluted Spambase train:\", X_train_pol.shape, \" test:\", X_test_pol.shape)\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "8554cb99",
"metadata": {},
"outputs": [],
"source": [
"# ---- TwoSpirals (binary, label in {-1,+1}) / ThreeCircles: 2-D points ----\n",
"# ThreeCircles is NOT binary despite the shared file format -- it's a balanced 3-class dataset,\n",
"# label in {-1, 0, +1} (333/333/334 points), so its chance level is ~33%, not ~50%.\n",
"# Same train/val/test convention as the other datasets: val is for Problem 5's T selection.\n",
"two_spirals = np.loadtxt(f\"{DATA}/TwoSpirals/twoSpirals.txt\")\n",
"three_circles = np.loadtxt(f\"{DATA}/TwoSpirals/threecircles.txt\", delimiter=\",\")\n",
"\n",
"X_train_spirals, X_temp_spirals, y_train_spirals, y_temp_spirals = train_test_split(\n",
" two_spirals[:, :2], two_spirals[:, 2], test_size=0.4, random_state=42, stratify=two_spirals[:, 2])\n",
"X_val_spirals, X_test_spirals, y_val_spirals, y_test_spirals = train_test_split(\n",
" X_temp_spirals, y_temp_spirals, test_size=0.5, random_state=42, stratify=y_temp_spirals)\n",
"\n",
"X_train_circles, X_temp_circles, y_train_circles, y_temp_circles = train_test_split(\n",
" three_circles[:, :2], three_circles[:, 2], test_size=0.4, random_state=42, stratify=three_circles[:, 2])\n",
"X_val_circles, X_test_circles, y_val_circles, y_test_circles = train_test_split(\n",
" X_temp_circles, y_temp_circles, test_size=0.5, random_state=42, stratify=y_temp_circles)\n",
"\n",
"print(\"TwoSpirals (2 classes) train:\", X_train_spirals.shape, \" val:\", X_val_spirals.shape, \" test:\", X_test_spirals.shape)\n",
"print(\"ThreeCircles (3 classes) train:\", X_train_circles.shape, \" val:\", X_val_circles.shape, \" test:\", X_test_circles.shape)\n"
]
},
{
"cell_type": "markdown",
"id": "06464364",
"metadata": {},
"source": [
"---\n",
"## Problem 1 — SVM (library), several kernels\n",
"\n",
"### Part A — Spambase (binary)\n",
"Try linear, polynomial, and RBF kernels; a small grid for each kernel's own hyperparameter.\n",
"Configurations are compared on the **validation** set; the test set is reported only once, for\n",
"the configurations finally chosen, so it never influences which hyperparameters get picked.\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "69b75e1c",
"metadata": {},
"outputs": [],
"source": [
"svm_grid_results = []\n",
"for kernel, param_grid in [\n",
" (\"linear\", [{\"C\": c} for c in [0.1, 1.0, 10.0]]),\n",
" (\"poly\", [{\"C\": 1.0, \"degree\": d, \"coef0\": c0} for d in [2, 3] for c0 in [0.0, 1.0]]),\n",
" (\"rbf\", [{\"C\": 1.0, \"gamma\": g} for g in [0.001, 0.01, 0.05, 0.1]]),\n",
"]:\n",
" for params in param_grid:\n",
" clf = SVC... ### TODO_STUDENT ### fit an SVC with this kernel and hyperparameters as the library SVM baseline\n",
" acc = clf.score(X_val_spam, y_val_spam)\n",
" svm_grid_results.append((kernel, params, acc, clf.support_.shape[0]))\n",
" print(f\"kernel={kernel:<7} params={params} val acc={acc:.4f} #SV={clf.support_.shape[0]}\")\n",
"\n",
"best_kernel, best_params, best_val_acc, _ = max(svm_grid_results, key=lambda r: r[2])\n",
"print(f\"\\nBest on validation: kernel={best_kernel}, params={best_params}, val acc={best_val_acc:.4f}\")\n"
]
},
{
"cell_type": "markdown",
"id": "def5e9ce",
"metadata": {},
"source": [
"**Best linear and RBF configurations on Spambase** (selected on validation, confirmed once on test). Problem 2's from-scratch SMO runs on a different dataset, TwoSpirals — see that problem for why."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "83cc9ecb",
"metadata": {},
"outputs": [],
"source": [
"best_linear = max((r for r in svm_grid_results if r[0] == \"linear\"), key=lambda r: r[2])\n",
"_, best_linear_params, best_linear_val_acc, best_linear_nsv = best_linear\n",
"BEST_LINEAR_C = best_linear_params[\"C\"]\n",
"\n",
"linear_baseline = SVC... ### TODO_STUDENT ### fit a linear-kernel SVC baseline at the best C found above\n",
"print(f\"Library linear-SVM baseline (selected on validation): C={BEST_LINEAR_C}\")\n",
"print(f\" train acc={linear_baseline.score(X_train_spam, y_train_spam):.4f} \"\n",
" f\"val acc={linear_baseline.score(X_val_spam, y_val_spam):.4f} \"\n",
" f\"test acc={linear_baseline.score(X_test_spam, y_test_spam):.4f} #SV={linear_baseline.support_.shape[0]}\")\n",
"\n",
"best_rbf = max((r for r in svm_grid_results if r[0] == \"rbf\"), key=lambda r: r[2])\n",
"_, best_rbf_params, best_rbf_val_acc, best_rbf_nsv = best_rbf\n",
"BEST_GAMMA = best_rbf_params[\"gamma\"]\n",
"BEST_C = best_rbf_params[\"C\"]\n",
"\n",
"rbf_baseline = SVC... ### TODO_STUDENT ### fit an RBF-kernel SVC baseline at the best (C, gamma) found above\n",
"print(f\"\\nLibrary RBF-SVM baseline (selected on validation): C={BEST_C}, gamma={BEST_GAMMA}\")\n",
"print(f\" train acc={rbf_baseline.score(X_train_spam, y_train_spam):.4f} \"\n",
" f\"val acc={rbf_baseline.score(X_val_spam, y_val_spam):.4f} \"\n",
" f\"test acc={rbf_baseline.score(X_test_spam, y_test_spam):.4f} #SV={rbf_baseline.support_.shape[0]}\")\n"
]
},
{
"cell_type": "markdown",
"id": "28bf3074",
"metadata": {},
"source": [
"### Part B — Digits/MNIST (multiclass, required)\\n\\n`sklearn.svm.SVC` handles multiclass natively (one-vs-one internally, 45 binary classifiers for 10 classes). We also build an explicit one-vs-rest wrapper (10 classifiers) with `OneVsRestClassifier` and compare the two schemes."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "83ce9f55",
"metadata": {},
"outputs": [],
"source": [
"svc_multiclass = SVC(kernel=\"rbf\", C=1.0, gamma=0.01, decision_function_shape=\"ovo\")\n",
"svc_multiclass... ### TODO_STUDENT ### fit sklearn's native multiclass SVC (one-vs-one internally) on Digits\n",
"acc_ovo = svc_multiclass.score(X_test_digits, y_test_digits)\n",
"print(f\"Native multiclass SVC (one-vs-one, 45 classifiers): test acc={acc_ovo:.4f}\")\n",
"\n",
"ovr = OneVsRestClassifier(SVC(kernel=\"rbf\", C=1.0, gamma=0.01))\n",
"ovr... ### TODO_STUDENT ### fit an explicit one-vs-rest multiclass wrapper on Digits\n",
"acc_ovr = ovr.score(X_test_digits, y_test_digits)\n",
"print(f\"Explicit one-vs-rest wrapper (10 classifiers): test acc={acc_ovr:.4f}\")\n"
]
},
{
"cell_type": "markdown",
"id": "423627e0",
"metadata": {},
"source": [
"---\n",
"## Problem 2 — SVM from scratch: Simplified SMO\n",
"\n",
"Implements the **CS229 Simplified SMO algorithm** (`materials/smo.pdf`, Section 3): iterate over\n",
"all αi, pick a random αj whenever αi violates the\n",
"KKT conditions, jointly optimize the pair, stop after `max_passes` sweeps with no change (capped by\n",
"`max_sweeps`, see below). This is deliberately *not* full Platt SMO (no error cache, no\n",
"entire-set/non-bound alternation, no second-choice heuristic) — that fuller version is one of the\n",
"optional extensions below.\n",
"\n",
"The solver itself never looks at raw features, only a precomputed kernel (Gram) matrix `K` — so\n",
"the exact same function works for any kernel; only the Gram matrix passed in changes.\n",
"\n",
"**Dataset choice: TwoSpirals, not Spambase.** On Spambase a linear kernel already gets 0.929 test\n",
"accuracy vs. RBF's 0.926 (Problem 1) — kernels aren't buying anything there, so it's a poor\n",
"demonstration of *why* the kernel trick matters. **TwoSpirals** (reused from Problem 5) makes the\n",
"point directly: it is not linearly separable by construction, so a linear-kernel SVM is stuck near\n",
"chance while an RBF-kernel SVM — the exact same solver, just a different Gram matrix — solves it\n",
"almost perfectly.\n",
"\n",
"**Step 1 — linear kernel, a correctness check (not expected to classify well here).** With a linear\n",
"kernel the dual solution gives back an explicit primal weight vector,\n",
"`w = ∑αiyixi`, which can be compared directly against\n",
"`SVC(kernel='linear').coef_` — a stronger correctness check than matching accuracy alone, though\n",
"still only a *capped-sweep* approximation on this non-separable data: expect the **direction** to\n",
"match well (cosine similarity ≈ 0.99) while the **magnitude** (relative L2 norm) does not,\n",
"since `max_sweeps` stops the run before full convergence. Matching accuracy and direction here still\n",
"confirms the optimizer is doing the right thing; the poor accuracy itself is the kernel's fault, not\n",
"the solver's.\n",
"\n",
"**Step 2 — RBF kernel, the required non-linear case** — identical solver code, RBF Gram matrix\n",
"substituted in, benchmarked against a library RBF-SVM on the same data.\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "7a9698b3",
"metadata": {},
"outputs": [],
"source": [
"def linear_gram(A, B):\n",
" return A @ B.T # given -- the linear kernel is just the dot product, nothing to derive\n",
"\n",
"\n",
"def rbf_gram(A, B, gamma):\n",
" sq_a = np.sum(A**2, axis=1)[:, None]\n",
" sq_b = np.sum(B**2, axis=1)[None, :]\n",
" sq_dist = sq_a + sq_b - 2 * A @ B.T\n",
" np.maximum(sq_dist, 0, out=sq_dist)\n",
" return ... ### TODO_STUDENT ### implement the RBF/Gaussian kernel formula given the squared distances\n",
"\n",
"\n",
"def simplified_smo(K, y, C=1.0, tol=1e-3, max_passes=5, max_sweeps=300, verbose=False):\n",
" \"\"\"CS229 'Simplified SMO' (materials/smo.pdf, eqs 10-19), given a precomputed kernel\n",
" matrix K. Kernel-agnostic: swap in a linear or RBF Gram matrix, nothing else changes.\n",
"\n",
" The handout's stopping rule (`max_passes` consecutive sweeps with no change) assumes\n",
" near-separable data and isn't guaranteed to trigger on noisy, non-separable real data --\n",
" on this dataset the linear kernel in particular keeps making small updates near the margin\n",
" for a very long time (still adjusting alphas at sweep 300 on TwoSpirals). `max_sweeps` is a\n",
" hard cap so runtime stays bounded regardless. Empirically the RBF kernel's *accuracy* already\n",
" matches the library's exactly well before the cap; the linear kernel's classification accuracy\n",
" and weight-vector *direction* (cosine similarity) also match closely, but its weight vector has\n",
" not fully converged in magnitude by sweep 300 -- see Step 1's printed relative difference.\n",
" \"\"\"\n",
" m = K.shape[0]\n",
" alpha = np.zeros(m)\n",
" b = 0.0\n",
" passes = 0\n",
"\n",
" def f(i):\n",
" return ... ### TODO_STUDENT ### implement the SVM decision function f(x_i) = sum_k alpha_k y_k K(x_k, x_i) + b\n",
"\n",
" sweep = 0\n",
" while passes < max_passes and sweep < max_sweeps:\n",
" sweep += 1\n",
" num_changed = 0\n",
" for i in range(m):\n",
" Ei = f(i) - y[i] # given, worked example -- E_j below (same formula, index j) is the one to fill in\n",
" if ...: ### TODO_STUDENT ### check whether alpha_i violates the KKT conditions (smo.pdf eqs 6-8)\n",
" j = i\n",
" while j == i:\n",
" j = np.random.randint(0, m) # given -- mechanical index bookkeeping, not part of the KKT/SMO math\n",
" Ej = ... ### TODO_STUDENT ### compute the error E_j = f(x_j) - y_j (mirrors E_i above)\n",
" alpha_i_old, alpha_j_old = alpha[i], alpha[j]\n",
"\n",
" if y[i] != y[j]:\n",
" # given, worked example -- the y_i == y_j case below is the mirror-image derivation to fill in\n",
" L = max(0.0, alpha[j] - alpha[i])\n",
" H = min(C, C + alpha[j] - alpha[i])\n",
" else:\n",
" L = ... ### TODO_STUDENT ### compute the lower bound L when y_i == y_j (smo.pdf eq 11) -- mirrors the y_i != y_j case above\n",
" H = ... ### TODO_STUDENT ### compute the upper bound H when y_i == y_j (smo.pdf eq 11) -- mirrors the y_i != y_j case above\n",
" if L == H:\n",
" continue\n",
"\n",
" eta = ... ### TODO_STUDENT ### compute eta, the curvature of the objective along the constraint line (smo.pdf eq 14)\n",
" if eta >= 0:\n",
" continue\n",
"\n",
" alpha[j] = ... ### TODO_STUDENT ### compute the unclipped update for alpha_j (smo.pdf eq 12)\n",
" alpha[j] = min(H, max(L, alpha[j])) # given -- mechanical clamp into [L, H] once L, H are known (smo.pdf eq 15)\n",
" if abs(alpha[j] - alpha_j_old) < 1e-5:\n",
" alpha[j] = alpha_j_old # given -- revert: no meaningful update, so undo it rather than\n",
" continue # leave alpha_j drifted with no compensating alpha_i change\n",
"\n",
" alpha[i] = ... ### TODO_STUDENT ### update alpha_i to offset the change in alpha_j (smo.pdf eq 16)\n",
"\n",
" # b1 is given, worked example -- b2 (the mirror-image derivation, eq 18) is the one to fill in\n",
" b1 = (b - Ei - y[i] * (alpha[i] - alpha_i_old) * K[i, i]\n",
" - y[j] * (alpha[j] - alpha_j_old) * K[i, j])\n",
" b2 = ... ### TODO_STUDENT ### compute the b threshold implied by alpha_j (smo.pdf eq 18) -- mirrors b1 above\n",
" if 0 < alpha[i] < C:\n",
" b = b1\n",
" elif 0 < alpha[j] < C:\n",
" b = b2\n",
" else:\n",
" b = ... ### TODO_STUDENT ### pick b1, b2, or their average per the three KKT cases (smo.pdf eq 19)\n",
" num_changed += 1\n",
" if verbose and sweep % 20 == 0:\n",
" print(f\" sweep {sweep}: alphas changed this sweep = {num_changed}\")\n",
" passes = passes + 1 if num_changed == 0 else 0\n",
"\n",
" sv = alpha > 1e-6\n",
" return alpha, b, sv\n",
"\n",
"\n",
"def predict_smo(X_train, y_train, alpha, b, gram_fn, X_query):\n",
" sv = alpha > 1e-6\n",
" K = gram_fn(X_train[sv], X_query)\n",
" scores = (alpha[sv] * y_train[sv]) @ K + b # given -- same decision-function formula as f(i) above, restricted to support vectors\n",
" return np.where(scores >= 0, 1.0, -1.0) # given -- threshold into +-1 (>=0, not np.sign, so an exact-zero score is +1 per convention, not the undefined 0)\n"
]
},
{
"cell_type": "markdown",
"id": "ecc06dde",
"metadata": {},
"source": [
"### Step 1 — linear kernel: validate against `SVC(kernel='linear')` by comparing weight vectors"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "fac36f1f",
"metadata": {},
"outputs": [],
"source": [
"import time\n",
"\n",
"SMO_C = 1.0\n",
"linear_baseline_spirals = SVC... ### TODO_STUDENT ### fit a linear-kernel SVC baseline on TwoSpirals\n",
"\n",
"K_lin_train = linear_gram(X_train_spirals, X_train_spirals)\n",
"t0 = time.time()\n",
"alpha_lin, b_lin, sv_lin = simplified_smo(\n",
" K_lin_train, y_train_spirals, C=SMO_C, tol=1e-3, max_passes=5, verbose=True)\n",
"t_lin = time.time() - t0\n",
"\n",
"w_scratch = ... ### TODO_STUDENT ### recover the primal weight vector w = sum(alpha_i y_i x_i) from the dual solution (only valid for a linear kernel)\n",
"w_lib = linear_baseline_spirals.coef_.ravel()\n",
"\n",
"w_rel_diff = np.linalg.norm(w_scratch - w_lib) / np.linalg.norm(w_lib)\n",
"cos_sim = np.dot(w_scratch, w_lib) / (np.linalg.norm(w_scratch) * np.linalg.norm(w_lib))\n",
"\n",
"pred_train_lin = predict_smo(X_train_spirals, y_train_spirals, alpha_lin, b_lin, linear_gram, X_train_spirals)\n",
"pred_test_lin = predict_smo(X_train_spirals, y_train_spirals, alpha_lin, b_lin, linear_gram, X_test_spirals)\n",
"acc_train_lin = np.mean(pred_train_lin == y_train_spirals)\n",
"acc_test_lin = np.mean(pred_test_lin == y_test_spirals)\n",
"\n",
"print(f\"\\n[Simplified SMO, linear kernel, TwoSpirals] C={SMO_C} #SV={sv_lin.sum()} train time={t_lin:.1f}s\")\n",
"print(f\" ||w_scratch - w_lib|| / ||w_lib|| = {w_rel_diff:.4f} cosine(w_scratch, w_lib) = {cos_sim:.4f}\")\n",
"print(f\" b_scratch={b_lin:.4f} b_lib={linear_baseline_spirals.intercept_[0]:.4f}\")\n",
"print(f\" train acc={acc_train_lin:.4f} test acc={acc_test_lin:.4f} (near chance, as expected -- \"\n",
" f\"the kernel, not the optimizer, is the bottleneck here)\")\n",
"print(f\"\\n[Library SVC linear baseline, TwoSpirals] C={SMO_C} \"\n",
" f\"train acc={linear_baseline_spirals.score(X_train_spirals, y_train_spirals):.4f} \"\n",
" f\"test acc={linear_baseline_spirals.score(X_test_spirals, y_test_spirals):.4f}\")\n"
]
},
{
"cell_type": "markdown",
"id": "4fccc67f",
"metadata": {},
"source": [
"### Step 2 — RBF kernel: same solver, RBF Gram matrix, the required non-linear case"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "861be21e",
"metadata": {},
"outputs": [],
"source": [
"SMO_GAMMA = 1.0\n",
"rbf_baseline_spirals = SVC... ### TODO_STUDENT ### fit an RBF-kernel SVC baseline on TwoSpirals\n",
"\n",
"K_rbf_train = rbf_gram(X_train_spirals, X_train_spirals, SMO_GAMMA)\n",
"t0 = time.time()\n",
"alpha_rbf, b_rbf, sv_rbf = simplified_smo(\n",
" K_rbf_train, y_train_spirals, C=SMO_C, tol=1e-3, max_passes=5, verbose=True)\n",
"t_rbf = time.time() - t0\n",
"\n",
"rbf_gram_fixed = lambda A, B: rbf_gram(A, B, SMO_GAMMA)\n",
"pred_train_rbf = predict_smo(X_train_spirals, y_train_spirals, alpha_rbf, b_rbf, rbf_gram_fixed, X_train_spirals)\n",
"pred_test_rbf = predict_smo(X_train_spirals, y_train_spirals, alpha_rbf, b_rbf, rbf_gram_fixed, X_test_spirals)\n",
"acc_train_rbf = np.mean(pred_train_rbf == y_train_spirals)\n",
"acc_test_rbf = np.mean(pred_test_rbf == y_test_spirals)\n",
"\n",
"print(f\"\\n[Simplified SMO, RBF kernel, TwoSpirals] C={SMO_C}, gamma={SMO_GAMMA} #SV={sv_rbf.sum()} \"\n",
" f\"train time={t_rbf:.1f}s\")\n",
"print(f\" train acc={acc_train_rbf:.4f} test acc={acc_test_rbf:.4f}\")\n",
"\n",
"print(f\"\\n[Library SVC RBF baseline, TwoSpirals] C={SMO_C}, gamma={SMO_GAMMA} \"\n",
" f\"train acc={rbf_baseline_spirals.score(X_train_spirals, y_train_spirals):.4f} \"\n",
" f\"test acc={rbf_baseline_spirals.score(X_test_spirals, y_test_spirals):.4f}\")\n"
]
},
{
"cell_type": "markdown",
"id": "865ba354",
"metadata": {},
"source": [
"**The point:** identical solver, identical data, one Gram matrix swapped in for another — linear stays near chance, RBF solves it almost perfectly."
]
},
{
"cell_type": "markdown",
"id": "2d663bb1",
"metadata": {},
"source": [
"**Optional extensions (not implemented here):** running this same solver on Spambase or\n",
"ThreeCircles; SMO with the *full* Platt heuristics (error cache, entire-set/non-bound alternation,\n",
"second-choice pair selection — see `6_SVM_kernels/code_HW6/bingyu/SVM/SVM.py` and `SMO.java` for a\n",
"reference implementation of that fuller version); running this SMO on the Digits dataset via a\n",
"multiclass wrapper.\n"
]
},
{
"cell_type": "markdown",
"id": "9de7a45c",
"metadata": {},
"source": [
"---\n",
"## Problem 3 — Kernel Ridge Regression from scratch\n",
"\n",
"Dual solution: α = (K + λI)-1y, prediction f(x) = ∑αiK(xi,x).\n",
"The kernel evaluation may use a library convention; the wrapper (`fit`/`predict`, solving the\n",
"linear system) is from scratch. Hyperparameters are selected on the **validation** set; test MSE\n",
"is reported once, for the chosen configuration, and separately checked against sklearn's\n",
"`KernelRidge` at the same (kernel, hyperparameters) as a library baseline.\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "7a6a560f",
"metadata": {},
"outputs": [],
"source": [
"class KernelRidgeRegression:\n",
" def __init__(self, kernel=\"rbf\", alpha=1.0, gamma=1.0, degree=3, coef0=1.0):\n",
" self.kernel = kernel\n",
" self.alpha = alpha\n",
" self.gamma = gamma\n",
" self.degree = degree\n",
" self.coef0 = coef0\n",
"\n",
" def _kernel_matrix(self, A, B):\n",
" if self.kernel == \"linear\":\n",
" return A @ B.T # given -- the linear kernel is just the dot product, nothing to derive\n",
" elif self.kernel == \"rbf\":\n",
" sq = np.sum(A**2, 1)[:, None] + np.sum(B**2, 1)[None, :] - 2 * A @ B.T\n",
" return ... ### TODO_STUDENT ### implement the RBF kernel given the squared distances\n",
" elif self.kernel == \"poly\":\n",
" return ... ### TODO_STUDENT ### implement the polynomial kernel (gamma* + coef0)^degree\n",
" raise ValueError(self.kernel)\n",
"\n",
" def fit(self, X, y):\n",
" self.X_train_ = X\n",
" K = self._kernel_matrix(X, X)\n",
" n = K.shape[0]\n",
" self.dual_coef_ = np.linalg.... ### TODO_STUDENT ### solve the kernel ridge dual system (K + lambda*I) @ dual_coef = y\n",
" return self\n",
"\n",
" def predict(self, X):\n",
" K = self._kernel_matrix(X, self.X_train_)\n",
" return ... ### TODO_STUDENT ### predict as a weighted sum of kernel evaluations, K @ dual_coef\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "a466fb60",
"metadata": {},
"outputs": [],
"source": [
"def tune_and_fit(X_train, y_train, X_val, y_val, kernel):\n",
" grids = {\n",
" \"linear\": [{\"alpha\": a} for a in [0.001, 0.01, 0.1, 1.0]],\n",
" \"rbf\": [{\"alpha\": a, \"gamma\": g}\n",
" for a in [0.001, 0.01, 0.1] for g in [0.5, 1.0, 2.0, 5.0, 10.0]],\n",
" \"poly\": [{\"alpha\": a, \"degree\": d}\n",
" for a in [0.001, 0.01, 0.1] for d in [2, 3, 5]],\n",
" }\n",
" best = None\n",
" for params in grids[kernel]:\n",
" model = KernelRidgeRegression(kernel=kernel, gamma=params.get(\"gamma\", 1.0),\n",
" alpha=params[\"alpha\"], degree=params.get(\"degree\", 3)).fit(X_train, y_train)\n",
" mse_val = mean_squared_error(y_val, model.predict(X_val))\n",
" if best is None or mse_val < best[1]:\n",
" best = (params, mse_val, model)\n",
" return best # (params, val_mse, fitted_model)\n",
"\n",
"\n",
"kr_data = {\n",
" \"WAVE\": (X_train_wave, X_val_wave, X_test_wave, y_train_wave, y_val_wave, y_test_wave),\n",
" \"CRESCENT\": (X_train_crescent, X_val_crescent, X_test_crescent, y_train_crescent, y_val_crescent, y_test_crescent),\n",
"}\n",
"kr_results = {}\n",
"for name, (Xtr, Xv, Xte, ytr, yv, yte) in kr_data.items():\n",
" print(f\"=== {name} ===\")\n",
" kr_results[name] = {}\n",
" for kernel in [\"linear\", \"rbf\", \"poly\"]:\n",
" params, mse_val, model = tune_and_fit(Xtr, ytr, Xv, yv, kernel)\n",
" mse_train = mean_squared_error(ytr, model.predict(Xtr))\n",
" mse_test = mean_squared_error(yte, model.predict(Xte))\n",
" kr_results[name][kernel] = (params, mse_train, mse_val, mse_test, model)\n",
" print(f\" {kernel:<7} params={params} (selected on val) \"\n",
" f\"train MSE={mse_train:.4f} val MSE={mse_val:.4f} test MSE={mse_test:.4f}\")\n"
]
},
{
"cell_type": "markdown",
"id": "96bada65",
"metadata": {},
"source": [
"**Library baseline: `sklearn.kernel_ridge.KernelRidge`**, at each dataset's validation-selected (kernel, hyperparameters), to confirm the from-scratch dual solve."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "ce9f92a5",
"metadata": {},
"outputs": [],
"source": [
"print(\"sklearn KernelRidge baseline, at each dataset's validation-selected (kernel, hyperparameters):\")\n",
"for name, (Xtr, Xv, Xte, ytr, yv, yte) in kr_data.items():\n",
" best_kernel = min(kr_results[name], key=lambda k: kr_results[name][k][2])\n",
" params, _, _, mse_test_scratch, model = kr_results[name][best_kernel]\n",
" sk_kwargs = {\"alpha\": params[\"alpha\"]}\n",
" if \"gamma\" in params:\n",
" sk_kwargs[\"gamma\"] = params[\"gamma\"]\n",
" if \"degree\" in params:\n",
" sk_kwargs[\"degree\"] = params[\"degree\"]\n",
" sk_model = KernelRidge... ### TODO_STUDENT ### fit sklearn's KernelRidge as a library baseline at the same (kernel, hyperparameters)\n",
" mse_test_sklearn = mean_squared_error(yte, sk_model.predict(Xte))\n",
" max_abs_diff = np.max(np.abs(sk_model.predict(Xte) - model.predict(Xte)))\n",
" print(f\" {name}: kernel={best_kernel} scratch test MSE={mse_test_scratch:.4f} \"\n",
" f\"sklearn test MSE={mse_test_sklearn:.4f} max|scratch-sklearn| on test predictions={max_abs_diff:.2e}\")\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "875eee9a",
"metadata": {},
"outputs": [],
"source": [
"fig, axes = plt.subplots(1, 2, figsize=(11, 4.5))\n",
"for ax, name in zip(axes, [\"WAVE\", \"CRESCENT\"]):\n",
" Xtr, Xv, Xte, ytr, yv, yte = kr_data[name]\n",
" order = np.argsort(Xte[:, 0])\n",
" ax.scatter(Xtr[:, 0], ytr, s=10, alpha=0.3, label=\"train\")\n",
" ax.scatter(Xte[:, 0], yte, s=14, alpha=0.6, label=\"test (true)\")\n",
" best_kernel = min(kr_results[name], key=lambda k: kr_results[name][k][2])\n",
" model = kr_results[name][best_kernel][4]\n",
" ax.plot(Xte[order, 0], model.predict(Xte)[order], \"r-\", lw=2,\n",
" label=f\"predicted ({best_kernel})\")\n",
" ax.set_title(f\"{name}: best kernel (on val) = {best_kernel}\")\n",
" ax.legend()\n",
"plt.tight_layout()\n",
"plt.show()\n"
]
},
{
"cell_type": "markdown",
"id": "ac62b208",
"metadata": {},
"source": [
"**Which kernel does better, and why.** WAVE oscillates several times over `x in [-1, 1]` — a\n",
"genuinely non-linear, locally-varying function, so the linear kernel underfits badly while RBF\n",
"(local, bandwidth-tunable) and polynomial (if given a high-enough degree) track it much more\n",
"closely. CRESCENT is a single smooth hump — close to a low-degree polynomial already — so a\n",
"low-degree polynomial or a wide-bandwidth RBF both fit it well, and the gap to the linear kernel\n",
"(which still can't capture the curvature) is smaller than for WAVE but still visible.\n"
]
},
{
"cell_type": "markdown",
"id": "86cfb917",
"metadata": {},
"source": [
"---\n",
"## Problem 4 — PCA from scratch (feature corruption & rescue)\n",
"\n",
"Downstream classifier: `sklearn.linear_model.LogisticRegression` (per the assignment note, reusing\n",
"HW2's classifier choice, not from-scratch here — the from-scratch piece in this problem is PCA).\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "dbc98cdc",
"metadata": {},
"outputs": [],
"source": [
"def fit_eval_logreg(X_train, y_train, X_test, y_test, max_iter=2000):\n",
" clf = LogisticRegression... ### TODO_STUDENT ### fit sklearn's LogisticRegression as the downstream classifier (library baseline, reused across (A)/(B)/(C) and Problem 5)\n",
" return accuracy_score(y_train, clf.predict(X_train)), accuracy_score(y_test, clf.predict(X_test))\n",
"\n",
"\n",
"# (0) control: same classifier on the ORIGINAL, unpolluted 57-feature Spambase (loaded in setup)\n",
"# -- this is the \"before\" number the (A) vs (B)/(C) comparison below is a drop/rescue relative to.\n",
"acc_train_clean, acc_test_clean = fit_eval_logreg(X_train_spam, y_train_spam, X_test_spam, y_test_spam)\n",
"print(f\"(0) Unpolluted Spambase, 57 dims (control): \"\n",
" f\"train acc={acc_train_clean:.4f} test acc={acc_test_clean:.4f}\")\n",
"\n",
"# (A) raw polluted features\n",
"acc_train_raw, acc_test_raw = fit_eval_logreg(X_train_pol_s, y_train_pol, X_test_pol_s, y_test_pol)\n",
"print(f\"(A) Raw polluted features ({X_train_pol_s.shape[1]} dims): \"\n",
" f\"train acc={acc_train_raw:.4f} test acc={acc_test_raw:.4f}\")\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "2724deb6",
"metadata": {},
"outputs": [],
"source": [
"# (B) library PCA -> ~100 features (fit on train only, applied to test)\n",
"# svd_solver=\"full\" (not the \"auto\" default, which picks randomized/approximate SVD for a\n",
"# matrix this shape): keeps this an apples-to-apples comparison with the from-scratch SVD\n",
"# in (C), which is also exact.\n",
"N_COMPONENTS = 100\n",
"lib_pca = PCA... ### TODO_STUDENT ### fit sklearn's PCA as the library baseline, reducing to N_COMPONENTS dimensions\n",
"X_train_pca_lib = lib_pca.transform(X_train_pol_s)\n",
"X_test_pca_lib = lib_pca.transform(X_test_pol_s)\n",
"\n",
"acc_train_lib, acc_test_lib = fit_eval_logreg(X_train_pca_lib, y_train_pol, X_test_pca_lib, y_test_pol)\n",
"print(f\"(B) sklearn PCA -> {N_COMPONENTS} dims: \"\n",
" f\"train acc={acc_train_lib:.4f} test acc={acc_test_lib:.4f}\")\n",
"print(f\" explained variance ratio (cumulative): {lib_pca.explained_variance_ratio_.sum():.4f}\")\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "70edee47",
"metadata": {},
"outputs": [],
"source": [
"# (C) own PCA\n",
"class MyPCA:\n",
" def __init__(self, n_components):\n",
" self.n_components = n_components\n",
"\n",
" def fit(self, X):\n",
" self.mean_ = X.mean(axis=0)\n",
" Xc = X - self.mean_\n",
" # SVD of the centered data is numerically nicer than eigendecomposing the covariance\n",
" # matrix directly, and gives the same principal components.\n",
" U, S, Vt = np.linalg.... ### TODO_STUDENT ### compute the SVD of the centered data to get the principal directions\n",
" self.components_ = ... ### TODO_STUDENT ### keep only the top n_components right-singular vectors (the principal axes)\n",
" self.explained_variance_ = (S[: self.n_components] ** 2) / (X.shape[0] - 1)\n",
" return self\n",
"\n",
" def transform(self, X):\n",
" return ... ### TODO_STUDENT ### project centered data onto the stored principal components\n",
"\n",
"\n",
"my_pca = MyPCA(n_components=N_COMPONENTS).fit(X_train_pol_s)\n",
"X_train_pca_mine = my_pca.transform(X_train_pol_s)\n",
"X_test_pca_mine = my_pca.transform(X_test_pol_s)\n",
"\n",
"acc_train_mine, acc_test_mine = fit_eval_logreg(X_train_pca_mine, y_train_pol, X_test_pca_mine, y_test_pol)\n",
"print(f\"(C) own PCA -> {N_COMPONENTS} dims: \"\n",
" f\"train acc={acc_train_mine:.4f} test acc={acc_test_mine:.4f}\")\n",
"\n",
"# sanity check: own PCA vs library PCA should span the same subspace (components can differ by sign)\n",
"per_component_corr = [abs(np.corrcoef(X_train_pca_lib[:, k], X_train_pca_mine[:, k])[0, 1])\n",
" for k in range(N_COMPONENTS)]\n",
"print(f\" |correlation| between library- and own-PCA components: \"\n",
" f\"min={min(per_component_corr):.4f} mean={np.mean(per_component_corr):.4f}\")\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "71a3351f",
"metadata": {},
"outputs": [],
"source": [
"print(f\"{'method':<28}{'train acc':>12}{'test acc':>12}\")\n",
"print(f\"{'(0) unpolluted (control)':<28}{acc_train_clean:>12.4f}{acc_test_clean:>12.4f}\")\n",
"print(f\"{'(A) raw polluted':<28}{acc_train_raw:>12.4f}{acc_test_raw:>12.4f}\")\n",
"print(f\"{'(B) library PCA':<28}{acc_train_lib:>12.4f}{acc_test_lib:>12.4f}\")\n",
"print(f\"{'(C) own PCA':<28}{acc_train_mine:>12.4f}{acc_test_mine:>12.4f}\")\n"
]
},
{
"cell_type": "markdown",
"id": "26bdc26d",
"metadata": {},
"source": [
"**Explanation.** The polluted features add many junk/duplicated dimensions on top of the\n",
"original 57 Spambase features; a linear classifier trained directly on all of them overfits the\n",
"noise, dropping test accuracy relative to the (0) unpolluted control above. PCA finds the ~100\n",
"directions of highest variance in the *training* data (fit once, then applied to test with the\n",
"stored mean/components — no leakage), which are dominated by the signal from the original\n",
"features, so the downstream classifier recovers most of the accuracy lost to the junk dimensions\n",
"— back close to, and sometimes above, the (0) control. The own-PCA implementation matches the\n",
"library one almost exactly (components agree up to sign, which a linear classifier is invariant\n",
"to).\n"
]
},
{
"cell_type": "markdown",
"id": "e39bb106",
"metadata": {},
"source": [
"---\n",
"## Problem 5 — Kernel PCA + linear classifier\n",
"\n",
"**Note:** TwoSpirals is binary (`{-1,+1}`, chance ≈ 50%). ThreeCircles, despite the shared\n",
"2-D-points file format, is a balanced **3-class** dataset (`{-1,0,+1}`, chance ≈ 33%) — see\n",
"the Data loading note above. Keep that in mind when reading the raw-classifier accuracies below.\n",
"\n",
"### (A) A linear classifier fails on the raw 2-D data\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "4705232b",
"metadata": {},
"outputs": [],
"source": [
"for name, (Xtr, Xte, ytr, yte) in {\n",
" \"TwoSpirals\": (X_train_spirals, X_test_spirals, y_train_spirals, y_test_spirals),\n",
" \"ThreeCircles\": (X_train_circles, X_test_circles, y_train_circles, y_test_circles),\n",
"}.items():\n",
" acc_train, acc_test = fit_eval_logreg(Xtr, ytr, Xte, yte)\n",
" n_classes = len(np.unique(ytr))\n",
" print(f\"{name} ({n_classes} classes, chance {'~50%' if n_classes == 2 else '~33%'}): \"\n",
" f\"raw-2D logistic regression train acc={acc_train:.4f} test acc={acc_test:.4f}\")\n"
]
},
{
"cell_type": "markdown",
"id": "eb450488",
"metadata": {},
"source": [
"### (B) Kernel PCA from scratch (Gaussian kernel)\n",
"\n",
"Build the Gram matrix, center it in feature space, eigendecompose, and project new points using\n",
"the standard out-of-sample kernel-PCA formula.\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "593f0195",
"metadata": {},
"outputs": [],
"source": [
"class MyKernelPCA:\n",
" def __init__(self, n_components, gamma=1.0):\n",
" self.n_components = n_components\n",
" self.gamma = gamma\n",
"\n",
" @staticmethod\n",
" def _rbf(A, B, gamma):\n",
" sq = np.sum(A**2, 1)[:, None] + np.sum(B**2, 1)[None, :] - 2 * A @ B.T\n",
" return ... ### TODO_STUDENT ### implement the RBF kernel given the squared distances\n",
"\n",
" def fit(self, X):\n",
" self.X_train_ = X\n",
" n = X.shape[0]\n",
" K = self._rbf(X, X, self.gamma) # given -- just calling the kernel already implemented above on the training data\n",
" one_n = np.full((n, n), 1.0 / n)\n",
" K_centered = ... ### TODO_STUDENT ### center the kernel matrix in feature space (the kernel-PCA centering formula)\n",
"\n",
" eigvals, eigvecs = np.linalg.... ### TODO_STUDENT ### eigendecompose the centered kernel matrix\n",
" order = np.argsort(eigvals)[::-1][: self.n_components] # given -- sort descending, keep the top n_components indices\n",
" eigvals, eigvecs = eigvals[order], eigvecs[:, order]\n",
" eigvals = np.maximum(eigvals, 1e-12)\n",
"\n",
" self.eigvals_ = eigvals\n",
" self.alphas_ = ... ### TODO_STUDENT ### normalize each eigenvector by sqrt of its eigenvalue, so transform() is a single matmul\n",
" self.K_train_ = K\n",
" return self\n",
"\n",
" def transform(self, X):\n",
" n_train = self.X_train_.shape[0]\n",
" K_test = self._rbf(X, self.X_train_, self.gamma) # [n_query, n_train]\n",
" row_mean_train = self.K_train_.mean(axis=0, keepdims=True) # [1, n_train]\n",
" total_mean_train = self.K_train_.mean()\n",
" row_mean_test = K_test.mean(axis=1, keepdims=True) # [n_query, 1]\n",
" K_test_centered = ... ### TODO_STUDENT ### center the test-time kernel values using the stored training kernel statistics\n",
" return ... ### TODO_STUDENT ### project the centered test kernel values onto the stored eigenvectors\n"
]
},
{
"cell_type": "markdown",
"id": "2a4c8c2d",
"metadata": {},
"source": [
"`T` is chosen on the **validation** set (smallest `T` reaching a \"good performance\" bar of 90%\n",
"validation accuracy, falling back to the largest `T` tried if none clear it), then confirmed once\n",
"on the held-out test set — same protocol as Problems 1 and 3.\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "a5235ce1",
"metadata": {},
"outputs": [],
"source": [
"GAMMA_KPCA = 1.0\n",
"T_VALUES = [2, 5, 10, 20, 50]\n",
"GOOD_ACC_THRESHOLD = 0.9\n",
"\n",
"kpca_data = {\n",
" \"TwoSpirals\": (X_train_spirals, X_val_spirals, X_test_spirals, y_train_spirals, y_val_spirals, y_test_spirals),\n",
" \"ThreeCircles\": (X_train_circles, X_val_circles, X_test_circles, y_train_circles, y_val_circles, y_test_circles),\n",
"}\n",
"kpca_results = {}\n",
"chosen_T = {}\n",
"for name, (Xtr, Xv, Xte, ytr, yv, yte) in kpca_data.items():\n",
" print(f\"=== {name} ===\")\n",
" kpca_results[name] = []\n",
" for T in T_VALUES:\n",
" kpca = MyKernelPCA(n_components=T, gamma=GAMMA_KPCA).fit(Xtr)\n",
" Xtr_kpca, Xv_kpca = kpca.transform(Xtr), kpca.transform(Xv)\n",
" acc_train, acc_val = fit_eval_logreg(Xtr_kpca, ytr, Xv_kpca, yv)\n",
" kpca_results[name].append((T, acc_train, acc_val))\n",
" print(f\" T={T:<3} train acc={acc_train:.4f} val acc={acc_val:.4f}\")\n",
" good_Ts = [T for T, _, acc_val in kpca_results[name] if acc_val >= GOOD_ACC_THRESHOLD]\n",
" chosen_T[name] = good_Ts[0] if good_Ts else T_VALUES[-1]\n",
" print(f\" -> smallest T with val acc >= {GOOD_ACC_THRESHOLD:.0%}: T={chosen_T[name]}\")\n",
"\n",
"print(\"\\nFinal held-out test confirmation, at each dataset's chosen T:\")\n",
"test_acc_at_chosen_T = {}\n",
"for name, (Xtr, Xv, Xte, ytr, yv, yte) in kpca_data.items():\n",
" T = chosen_T[name]\n",
" kpca = MyKernelPCA(n_components=T, gamma=GAMMA_KPCA).fit(Xtr)\n",
" _, acc_test = fit_eval_logreg(kpca.transform(Xtr), ytr, kpca.transform(Xte), yte)\n",
" test_acc_at_chosen_T[name] = acc_test\n",
" print(f\" {name}: T={T} test acc={acc_test:.4f}\")\n"
]
},
{
"cell_type": "markdown",
"id": "fabf2093",
"metadata": {},
"source": [
"**Library baseline: `sklearn.decomposition.KernelPCA`**, at each dataset's chosen `T`."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "d15b5ce3",
"metadata": {},
"outputs": [],
"source": [
"print(\"sklearn KernelPCA baseline, at each dataset's chosen T:\")\n",
"for name, (Xtr, Xv, Xte, ytr, yv, yte) in kpca_data.items():\n",
" T = chosen_T[name]\n",
" sk_kpca = KernelPCA... ### TODO_STUDENT ### fit sklearn's KernelPCA as a library baseline at the chosen T\n",
" Xtr_sk, Xte_sk = sk_kpca.transform(Xtr), sk_kpca.transform(Xte)\n",
" _, acc_test_sk = fit_eval_logreg(Xtr_sk, ytr, Xte_sk, yte)\n",
" print(f\" {name}: T={T} scratch test acc={test_acc_at_chosen_T[name]:.4f} \"\n",
" f\"sklearn KernelPCA -> logistic test acc={acc_test_sk:.4f}\")\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "70145823",
"metadata": {},
"outputs": [],
"source": [
"fig, axes = plt.subplots(1, 2, figsize=(11, 4))\n",
"for ax, name in zip(axes, [\"TwoSpirals\", \"ThreeCircles\"]):\n",
" Ts = [r[0] for r in kpca_results[name]]\n",
" val_accs = [r[2] for r in kpca_results[name]]\n",
" ax.plot(Ts, val_accs, \"o-\")\n",
" ax.axhline(GOOD_ACC_THRESHOLD, color=\"gray\", ls=\"--\", lw=1, label=\"good-performance bar\")\n",
" ax.axvline(chosen_T[name], color=\"r\", ls=\":\", lw=1, label=f\"chosen T={chosen_T[name]}\")\n",
" ax.set_xlabel(\"T (kernel-PCA dimensions)\")\n",
" ax.set_ylabel(\"validation accuracy\")\n",
" ax.set_title(name)\n",
" ax.set_ylim(0.2, 1.02)\n",
" ax.legend()\n",
"plt.tight_layout()\n",
"plt.show()\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "d13bece8",
"metadata": {},
"outputs": [],
"source": [
"print(f'How large must T be? (>= {GOOD_ACC_THRESHOLD:.0%} validation accuracy counts as \"good\")')\n",
"for name in [\"TwoSpirals\", \"ThreeCircles\"]:\n",
" print(f\" {name}: smallest T reaching the bar = {chosen_T[name]}\")\n",
"print()\n",
"print(\"ThreeCircles' three classes separate almost immediately once the RBF kernel unfolds the \"\n",
" \"circle structure; TwoSpirals' tighter interleaved geometry needs far more of the centered \"\n",
" \"kernel matrix's top eigenvectors before a linear classifier can separate it well -- these \"\n",
" \"two datasets are not equally easy, despite both being 2-D and both being 'unfolded' by the \"\n",
" \"same kernel PCA machinery.\")\n"
]
}
],
"metadata": {},
"nbformat": 4,
"nbformat_minor": 5
}