{ "cells": [ { "cell_type": "markdown", "id": "e8c4a792", "metadata": {}, "source": [ "# CS6140 Machine Learning — Fall 2026\n", "# HW2 Starter — Gradient Descent, L1/Lasso, Perceptron, Boosting for Classification, Active Learning\n", "\n", "This notebook implements `HW2_26F.html`. Where a well-known library offers the same functionality\n", "you'll run it first as a sanity-check baseline, then implement the same thing from scratch and\n", "evaluate it exactly the same way, so your two numbers should land right next to each other once\n", "everything is filled in.\n", "\n", "**Problems 1-5 (required):**\n", "1. Gradient descent: L2-regularized linear regression + logistic regression (Housing, Spambase)\n", "2. L1 / Lasso regularization and feature selection (Spambase, Housing, 20 Newsgroups)\n", "3. Perceptron from scratch (linearly-separable toy data)\n", "4. Gradient boosting for classification, reusing the HW1 decision tree (Spambase; multiclass\n", " Newsgroups stretch)\n", "5. Active learning by uncertainty sampling, reusing the classifiers from Problems 1 and 4\n", "\n", "**Problem 6 (optional, no credit):**\n", "6. Lasso from scratch by coordinate descent, verified against sklearn's `Lasso`\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` or produce obviously wrong output until you fill it in — that's expected, not\n", "a bug. Work top to bottom: Problem 4's boosting reuses the `DecisionTree` class you fill in earlier\n", "in that same problem; Problem 5 reuses the logistic regression from Problem 1 and the\n", "gradient-boosted classifier from Problem 4, so finish those `TODO`s first. Problem 2 is entirely\n", "library code (no from-scratch component) and Problem 6 is a standalone optional problem that only\n", "needs Problem 1's Spambase data and the `Lasso` import from Problem 2. Once a \"library baseline\"\n", "cell's TODOs are filled in and run, the printed number is a target: your matching \"from scratch\"\n", "implementation just below it should land on (or very near) the same value.\n", "\n", "Datasets are read from the shared course `data/` folder (paths are relative to\n", "`2_GD_REG_pton_NN/code_hw/`, 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": "048071d0", "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, confusion_matrix, ConfusionMatrixDisplay, mean_squared_error as mse_loss\n", "\n", "DATA = \"../../data\"\n", "np.random.seed(42)\n" ] }, { "cell_type": "markdown", "id": "4e153ccf", "metadata": {}, "source": [ "## Data loading\\n\\nHousing comes pre-split (regression). Spambase is one file, so we split it ourselves -- same convention as HW1. Both are standardized using *training* statistics only (per the assignment's normalization note), applied to train and test alike." ] }, { "cell_type": "code", "execution_count": null, "id": "3bcfb989", "metadata": {}, "outputs": [], "source": [ "# ---- Housing (regression): already split into train/test by the course ----\n", "train_housing = np.loadtxt(f\"{DATA}/housing_train.txt\")\n", "test_housing = np.loadtxt(f\"{DATA}/housing_test.txt\")\n", "\n", "X_train_housing, y_train_housing = train_housing[:, :-1], train_housing[:, -1]\n", "X_test_housing, y_test_housing = test_housing[:, :-1], test_housing[:, -1]\n", "\n", "housing_scaler = StandardScaler().fit(X_train_housing)\n", "X_train_housing = housing_scaler.transform(X_train_housing)\n", "X_test_housing = housing_scaler.transform(X_test_housing)\n", "\n", "print(\"Housing train:\", X_train_housing.shape, \" test:\", X_test_housing.shape)\n" ] }, { "cell_type": "code", "execution_count": null, "id": "7240b9e7", "metadata": {}, "outputs": [], "source": [ "# ---- Spambase (binary classification): one file, we split it ----\n", "spam_raw = np.loadtxt(f\"{DATA}/spambase/spambase.data\", delimiter=\",\")\n", "X_spam, y_spam = spam_raw[:, :-1], spam_raw[:, -1]\n", "\n", "X_train_spam, X_test_spam, y_train_spam, y_test_spam = train_test_split(\n", " X_spam, y_spam, test_size=0.2, random_state=42)\n", "\n", "spam_scaler = StandardScaler().fit(X_train_spam)\n", "X_train_spam = spam_scaler.transform(X_train_spam)\n", "X_test_spam = spam_scaler.transform(X_test_spam)\n", "\n", "print(\"Spambase train:\", X_train_spam.shape, \" test:\", X_test_spam.shape)\n" ] }, { "cell_type": "markdown", "id": "e83de52d", "metadata": {}, "source": [ "---\n", "## Problem 1 — Gradient Descent: Linear (L2-regularized) + Logistic Regression\n", "\n", "Both models share the same two helpers: fold a bias column into `X`, and run batch gradient\n", "descent from a fixed random initialization. `MyOLSGD` also supports an L2 penalty (a genuine\n", "`+ 2*lambda*theta` term in the gradient, bias excluded) so it matches the assignment's\n", "\"L2-regularized linear regression by GD\" -- not just plain OLS.\n" ] }, { "cell_type": "code", "execution_count": null, "id": "49548fe8", "metadata": {}, "outputs": [], "source": [ "class MyOLSGD:\n", " \"\"\"Linear regression fit by batch gradient descent, with an optional L2 (ridge) penalty.\"\"\"\n", "\n", " def __init__(self, lr=1e-3, epochs=2000, lmbda=0.0, bias=True):\n", " self.lr = lr\n", " self.epochs = epochs\n", " self.lmbda = lmbda\n", " self.bias = bias\n", " self.theta = None\n", "\n", " def _add_bias(self, X):\n", " if X.ndim == 1:\n", " return np.append(X, 1.0) if self.bias else X\n", " n = X.shape[0]\n", " return np.concatenate([X, np.ones((n, 1))], axis=1) if self.bias else X\n", "\n", " def fit(self, X, y):\n", " Xb = self._add_bias(X)\n", " n, m = Xb.shape\n", " rng = np.random.RandomState(0)\n", " self.theta = rng.normal(size=m)\n", "\n", " for _ in range(self.epochs):\n", " reg = 2 * self.... ### TODO_STUDENT ### implement the L2 penalty gradient term, 2*lambda*theta\n", " if self.bias:\n", " reg[-1] = ... ### TODO_STUDENT ### don't regularize the bias/intercept term\n", " grad = ... ### TODO_STUDENT ### implement the (regularized) mean-squared-error gradient w.r.t. theta\n", " self.theta = self.theta - self.... ### TODO_STUDENT ### apply the gradient descent update step\n", " return self\n", "\n", " def decision_function(self, X):\n", " \"\"\"Raw (unthresholded) linear score -- the regression output itself.\"\"\"\n", " return self._add_bias(X) @ self.... ### TODO_STUDENT ### compute the linear score as X @ theta\n", "\n", " def predict(self, X):\n", " return self.decision_function(X)\n" ] }, { "cell_type": "code", "execution_count": null, "id": "703b7f45", "metadata": {}, "outputs": [], "source": [ "# Housing: L2-regularized linear regression by GD\n", "housing_gd = MyOLSGD(lr=0.05, epochs=2000, lmbda=1.0).fit(X_train_housing, y_train_housing)\n", "mse_train = mse_loss(y_train_housing, housing_gd.predict(X_train_housing))\n", "mse_test = mse_loss(y_test_housing, housing_gd.predict(X_test_housing))\n", "print(f\"[scratch] Housing linear GD (L2, lambda=1): train MSE = {mse_train:.3f}, test MSE = {mse_test:.3f}\")\n", "\n", "# Spambase: same model, real-valued output thresholded at 0.5 for classification\n", "spam_gd = MyOLSGD(lr=0.05, epochs=2000, lmbda=1.0).fit(X_train_spam, y_train_spam)\n", "acc_train = accuracy_score(y_train_spam, spam_gd.predict(X_train_spam) >= 0.5)\n", "acc_test = accuracy_score(y_test_spam, spam_gd.predict(X_test_spam) >= 0.5)\n", "print(f\"[scratch] Spambase linear GD (L2, threshold=0.5): train acc = {acc_train:.3f}, test acc = {acc_test:.3f}\")\n" ] }, { "cell_type": "markdown", "id": "b3efbb4b", "metadata": {}, "source": [ "### Logistic regression by GD" ] }, { "cell_type": "code", "execution_count": null, "id": "7d6b291f", "metadata": {}, "outputs": [], "source": [ "class MyLogReg:\n", " \"\"\"Logistic regression fit by batch gradient descent on the (unregularized) log-loss.\"\"\"\n", "\n", " def __init__(self, lr=1e-2, epochs=5000, bias=True):\n", " self.lr = lr\n", " self.epochs = epochs\n", " self.bias = bias\n", " self.theta = None\n", "\n", " def _add_bias(self, X):\n", " if X.ndim == 1:\n", " return np.append(X, 1.0) if self.bias else X\n", " n = X.shape[0]\n", " return np.concatenate([X, np.ones((n, 1))], axis=1) if self.bias else X\n", "\n", " def sigmoid(self, z):\n", " return ... ### TODO_STUDENT ### implement the sigmoid function\n", "\n", " def fit(self, X, y):\n", " Xb = self._add_bias(X)\n", " n, m = Xb.shape\n", " rng = np.random.RandomState(0)\n", " self.theta = rng.normal(size=m)\n", "\n", " for _ in range(self.epochs):\n", " p = self.... ### TODO_STUDENT ### compute the predicted probability p = sigmoid(X @ theta)\n", " grad = (1 / n) * Xb.... ### TODO_STUDENT ### implement the log-loss gradient w.r.t. theta\n", " self.theta = self.theta - self.... ### TODO_STUDENT ### apply the gradient descent update step\n", " return self\n", "\n", " def predict_proba(self, X):\n", " \"\"\"p(y=1|x) -- the probability of the positive class.\"\"\"\n", " return ... ### TODO_STUDENT ### compute p(y=1|x) = sigmoid(X @ theta)\n", "\n", " def predict(self, X, threshold=0.5):\n", " return ... ### TODO_STUDENT ### threshold the predicted probability into a 0/1 class label\n" ] }, { "cell_type": "code", "execution_count": null, "id": "8941f9e2", "metadata": {}, "outputs": [], "source": [ "spam_lr = MyLogReg(lr=0.1, epochs=5000).fit(X_train_spam, y_train_spam)\n", "acc_train = accuracy_score(y_train_spam, spam_lr.predict(X_train_spam))\n", "acc_test = accuracy_score(y_test_spam, spam_lr.predict(X_test_spam))\n", "print(f\"[scratch] Spambase logistic GD: train acc = {acc_train:.3f}, test acc = {acc_test:.3f}\")\n" ] }, { "cell_type": "markdown", "id": "50fda079", "metadata": {}, "source": [ "### Library baseline: sklearn `Ridge` / `LogisticRegression`\n", "\n", "Sanity check, not an exact-match requirement -- GD should reach a *comparable* answer to the\n", "closed-form/library solver, not necessarily an identical one.\n" ] }, { "cell_type": "code", "execution_count": null, "id": "97a8f684", "metadata": {}, "outputs": [], "source": [ "from sklearn.linear_model import Ridge, LogisticRegression as SkLogisticRegression\n", "\n", "# sklearn's Ridge minimizes the *summed* squared error + alpha*||w||^2, while our GD gradient comes\n", "# from the *mean* squared error + lmbda*||w||^2 -- same lambda=1 means a very different regularization\n", "# strength unless we correct for that factor of n. alpha = lmbda * n_train makes the two penalties\n", "# equivalent (sum = n * mean), so this is the actually-matching library baseline.\n", "sk_ridge_house = Ridge... ### TODO_STUDENT ### fit sklearn's Ridge (alpha = lambda * n_train) as the library baseline on Housing\n", "mse_train_lib = mse_loss(y_train_housing, sk_ridge_house.predict(X_train_housing))\n", "mse_test_lib = mse_loss(y_test_housing, sk_ridge_house.predict(X_test_housing))\n", "print(f\"[library] Housing Ridge(alpha=lambda*n_train): train MSE = {mse_train_lib:.3f}, test MSE = {mse_test_lib:.3f}\")\n", "\n", "sk_ridge_spam = Ridge... ### TODO_STUDENT ### fit sklearn's Ridge (alpha = lambda * n_train) as the library baseline on Spambase\n", "acc_train_lib = accuracy_score(y_train_spam, sk_ridge_spam.predict(X_train_spam) >= 0.5)\n", "acc_test_lib = accuracy_score(y_test_spam, sk_ridge_spam.predict(X_test_spam) >= 0.5)\n", "print(f\"[library] Spambase Ridge(alpha=lambda*n_train, threshold=0.5): train acc = {acc_train_lib:.3f}, test acc = {acc_test_lib:.3f}\")\n", "\n", "sk_logreg_spam = SkLogisticRegression... ### TODO_STUDENT ### fit sklearn's LogisticRegression(penalty=None) as the library baseline on Spambase\n", "acc_train_lib = accuracy_score(y_train_spam, sk_logreg_spam.predict(X_train_spam))\n", "acc_test_lib = accuracy_score(y_test_spam, sk_logreg_spam.predict(X_test_spam))\n", "print(f\"[library] Spambase LogisticRegression (unregularized): train acc = {acc_train_lib:.3f}, test acc = {acc_test_lib:.3f}\")\n" ] }, { "cell_type": "markdown", "id": "406b61ce", "metadata": {}, "source": [ "### Confusion matrices and ROC/AUC, linear vs. logistic (Spambase)" ] }, { "cell_type": "code", "execution_count": null, "id": "90fb3914", "metadata": {}, "outputs": [], "source": [ "fig, axes = plt.subplots(1, 2, figsize=(10, 4))\n", "\n", "gd_cm = confusion_matrix(y_test_spam, spam_gd.predict(X_test_spam) >= 0.5)\n", "ConfusionMatrixDisplay(gd_cm, display_labels=[0, 1]).plot(ax=axes[0], colorbar=False)\n", "axes[0].set_title(\"Linear regression (thresholded), test\")\n", "\n", "lr_cm = confusion_matrix(y_test_spam, spam_lr.predict(X_test_spam))\n", "ConfusionMatrixDisplay(lr_cm, display_labels=[0, 1]).plot(ax=axes[1], colorbar=False)\n", "axes[1].set_title(\"Logistic regression, test\")\n", "\n", "plt.tight_layout()\n", "plt.show()\n" ] }, { "cell_type": "code", "execution_count": null, "id": "1135f407", "metadata": {}, "outputs": [], "source": [ "def roc_curve_manual(y_true, scores, n_thresholds=200):\n", " \"\"\"Sweep thresholds across the scores' own range (works for both bounded probabilities and\n", " unbounded linear-regression outputs), then sort by FPR so trapezoid integration is valid.\"\"\"\n", " thresholds = np.linspace(scores.min(), scores.max(), n_thresholds)\n", " tprs, fprs = [], []\n", " for t in thresholds:\n", " y_pred = (scores >= t).astype(int)\n", " tprs.append(np.sum((y_true == 1) & (y_pred == 1)) / np.sum(y_true == 1))\n", " fprs.append(np.sum((y_true == 0) & (y_pred == 1)) / np.sum(y_true == 0))\n", " order = np.argsort(fprs)\n", " fprs, tprs = np.array(fprs)[order], np.array(tprs)[order]\n", " auc = np.trapezoid(tprs, fprs)\n", " return fprs, tprs, auc\n", "\n", "fpr_gd, tpr_gd, auc_gd = roc_curve_manual(y_test_spam, spam_gd.decision_function(X_test_spam))\n", "fpr_lr, tpr_lr, auc_lr = roc_curve_manual(y_test_spam, spam_lr.predict_proba(X_test_spam))\n", "\n", "plt.figure(figsize=(5, 5))\n", "plt.plot(fpr_gd, tpr_gd, color=\"#2a78d6\", label=f\"linear regression (AUC={auc_gd:.3f})\")\n", "plt.plot(fpr_lr, tpr_lr, color=\"#e34948\", label=f\"logistic regression (AUC={auc_lr:.3f})\")\n", "plt.plot([0, 1], [0, 1], \"k--\", linewidth=1, label=\"chance\")\n", "plt.xlabel(\"False Positive Rate\")\n", "plt.ylabel(\"True Positive Rate\")\n", "plt.title(\"Spambase: ROC, linear vs. logistic regression\")\n", "plt.legend()\n", "plt.show()\n" ] }, { "cell_type": "markdown", "id": "56fa39a0", "metadata": {}, "source": [ "---\n", "## Problem 2 — L1 / Lasso Regularization and Feature Selection\n", "\n", "### (a) Spambase — what the assignment actually asks for\n", "\n", "Sweep the L1 penalty, threshold the (real-valued) Lasso output at 0.5 for classification, and\n", "look at two things together: does test accuracy hold up as the penalty grows, and how many\n", "features does Lasso actually keep (nonzero coefficients) at each penalty level?\n" ] }, { "cell_type": "code", "execution_count": null, "id": "7ced5ca8", "metadata": {}, "outputs": [], "source": [ "from sklearn.linear_model import Lasso\n", "\n", "alphas = np.logspace(-4, 0, 30)\n", "train_accs, test_accs, n_selected = [], [], []\n", "for a in alphas:\n", " lasso_spam = Lasso... ### TODO_STUDENT ### fit sklearn's Lasso(alpha=a) on Spambase at this penalty value\n", " train_accs.append(...) ### TODO_STUDENT ### threshold the Lasso prediction at 0.5 to get a class label, then compute train accuracy\n", " test_accs.append(...) ### TODO_STUDENT ### threshold the Lasso prediction at 0.5 to get a class label, then compute test accuracy\n", " n_selected.append(...) ### TODO_STUDENT ### count the number of nonzero (selected) coefficients at this penalty value\n", "\n", "fig, axes = plt.subplots(1, 2, figsize=(11, 4))\n", "axes[0].plot(alphas, train_accs, color=\"#2a78d6\", marker=\"o\", markersize=4, label=\"train\")\n", "axes[0].plot(alphas, test_accs, color=\"#e34948\", marker=\"o\", markersize=4, label=\"test\")\n", "axes[0].set_xscale(\"log\")\n", "axes[0].set_xlabel(\"L1 penalty (alpha)\")\n", "axes[0].set_ylabel(\"Accuracy (threshold=0.5)\")\n", "axes[0].set_title(\"Spambase: Lasso penalty vs. accuracy\")\n", "axes[0].legend()\n", "\n", "axes[1].plot(alphas, n_selected, color=\"#2a78d6\", marker=\"o\", markersize=4)\n", "axes[1].set_xscale(\"log\")\n", "axes[1].set_xlabel(\"L1 penalty (alpha)\")\n", "axes[1].set_ylabel(\"Nonzero coefficients (of 57)\")\n", "axes[1].set_title(\"Spambase: Lasso feature selection vs. penalty\")\n", "\n", "plt.tight_layout()\n", "plt.show()\n" ] }, { "cell_type": "markdown", "id": "f93b875d", "metadata": {}, "source": [ "### (b) Housing — penalty vs. train/test MSE" ] }, { "cell_type": "code", "execution_count": null, "id": "0f79b7d3", "metadata": {}, "outputs": [], "source": [ "lasso_lambdas = np.arange(0.01, 2, 0.01)\n", "train_mses, test_mses = [], []\n", "for l in lasso_lambdas:\n", " lasso_house = Lasso(alpha=l, max_iter=10000).fit(X_train_housing, y_train_housing)\n", " train_mses.append(mse_loss(y_train_housing, lasso_house.predict(X_train_housing)))\n", " test_mses.append(mse_loss(y_test_housing, lasso_house.predict(X_test_housing)))\n", "\n", "plt.figure(figsize=(6, 4))\n", "plt.plot(lasso_lambdas, train_mses, color=\"#2a78d6\", label=\"train MSE\")\n", "plt.plot(lasso_lambdas, test_mses, color=\"#e34948\", label=\"test MSE\")\n", "plt.xlabel(\"Lambda (alpha)\")\n", "plt.ylabel(\"MSE\")\n", "plt.title(\"Housing: Lasso penalty vs. train/test MSE\")\n", "plt.legend()\n", "plt.show()\n" ] }, { "cell_type": "markdown", "id": "809dd500", "metadata": {}, "source": [ "### (c) 20 Newsgroups (8-class subset) — L1 feature selection on text\n", "\n", "L1 logistic regression (one-vs-rest) ranks features by average |coefficient| across the 8\n", "one-vs-rest classifiers; we keep the top 200 and refit an L2 logistic model on just those, then\n", "report per-class and overall accuracy. (The course also has a local, pre-extracted\n", "`data/20newsgroup.zip` / `data/8newsgroup.zip` for this -- either source works fine here, it isn't\n", "an important detail for this HW, so we use the standard sklearn fetch for convenience.)\n" ] }, { "cell_type": "code", "execution_count": null, "id": "9fcba886", "metadata": {}, "outputs": [], "source": [ "from sklearn.datasets import fetch_20newsgroups\n", "from sklearn.feature_extraction.text import TfidfVectorizer\n", "from sklearn.multiclass import OneVsRestClassifier\n", "from sklearn.linear_model import LogisticRegression\n", "\n", "categories = [\n", " \"rec.sport.hockey\", \"sci.space\", \"comp.graphics\", \"talk.politics.mideast\",\n", " \"alt.atheism\", \"soc.religion.christian\", \"misc.forsale\", \"rec.autos\",\n", "]\n", "newsgroups = fetch_20newsgroups(subset=\"test\", categories=categories, remove=(\"headers\", \"footers\", \"quotes\"))\n", "\n", "vectorizer = TfidfVectorizer(stop_words=\"english\", max_features=3000)\n", "X_ng = vectorizer.fit_transform(newsgroups.data).toarray()\n", "y_ng = newsgroups.target\n", "\n", "X_train_ng, X_test_ng, y_train_ng, y_test_ng = train_test_split(X_ng, y_ng, test_size=0.2, random_state=42)\n", "print(\"Newsgroups (8-class):\", X_train_ng.shape, X_test_ng.shape)\n" ] }, { "cell_type": "code", "execution_count": null, "id": "25744c40", "metadata": {}, "outputs": [], "source": [ "l1_log_reg = OneVsRestClassifier(LogisticRegression(penalty=\"l1\", solver=\"liblinear\", max_iter=2000))\n", "l1_log_reg.fit(X_train_ng, y_train_ng)\n", "\n", "n_classes, n_features = len(l1_log_reg.estimators_), X_train_ng.shape[1]\n", "subclassifier_coefs = np.zeros((n_classes, n_features))\n", "for i, subclassifier in enumerate(l1_log_reg.estimators_):\n", " subclassifier_coefs[i, :] = np.abs(subclassifier.coef_)\n", "\n", "average_coefs = np.mean(subclassifier_coefs, axis=0)\n", "top_200 = np.argsort(-average_coefs)[:200]\n", "\n", "X_train_ng200 = X_train_ng[:, top_200]\n", "X_test_ng200 = X_test_ng[:, top_200]\n", "print(\"Selected feature matrix:\", X_train_ng200.shape)\n" ] }, { "cell_type": "code", "execution_count": null, "id": "297243c4", "metadata": {}, "outputs": [], "source": [ "l2_log_reg = OneVsRestClassifier(LogisticRegression(penalty=\"l2\", max_iter=2000))\n", "l2_log_reg.fit(X_train_ng200, y_train_ng)\n", "\n", "y_train_pred = l2_log_reg.predict(X_train_ng200)\n", "y_test_pred = l2_log_reg.predict(X_test_ng200)\n", "\n", "for k, name in enumerate(newsgroups.target_names):\n", " mask_train, mask_test = y_train_ng == k, y_test_ng == k\n", " print(f\" {name:24s} train acc = {np.mean(y_train_pred[mask_train] == k):.3f}, \"\n", " f\"test acc = {np.mean(y_test_pred[mask_test] == k):.3f}\")\n", "\n", "print(f\"Overall: train acc = {accuracy_score(y_train_ng, y_train_pred):.3f}, \"\n", " f\"test acc = {accuracy_score(y_test_ng, y_test_pred):.3f}\")\n" ] }, { "cell_type": "markdown", "id": "2933cc94", "metadata": {}, "source": [ "---\n", "## Problem 3 — Perceptron from Scratch\n", "\n", "Mistake-driven updates on a linearly-separable toy dataset (`perceptronData.txt`, 4 features + a\n", "+-1 label). We track the number of mistakes made per full pass over the data; it should hit 0 once\n", "the weight vector separates the classes.\n" ] }, { "cell_type": "code", "execution_count": null, "id": "2cb1c059", "metadata": {}, "outputs": [], "source": [ "perceptron_data = np.loadtxt(f\"{DATA}/perceptronData.txt\")\n", "X_perceptron, y_perceptron = perceptron_data[:, :-1], perceptron_data[:, -1].astype(int)\n", "print(\"Perceptron data:\", X_perceptron.shape, y_perceptron.shape)\n" ] }, { "cell_type": "code", "execution_count": null, "id": "b66c8b8a", "metadata": {}, "outputs": [], "source": [ "class MyPerceptron:\n", " \"\"\"Classic mistake-driven perceptron: on a misclassified point, nudge theta toward its label.\"\"\"\n", "\n", " def __init__(self, lr=1.0, epochs=100, bias=True):\n", " self.lr = lr\n", " self.epochs = epochs\n", " self.bias = bias\n", " self.theta = None\n", "\n", " def _add_bias(self, X):\n", " if X.ndim == 1:\n", " return np.insert(X, 0, 1.0) if self.bias else X\n", " n = X.shape[0]\n", " return np.concatenate([np.ones((n, 1)), X], axis=1) if self.bias else X\n", "\n", " def predict(self, X):\n", " return np.... ### TODO_STUDENT ### classify by the sign of theta . x\n", "\n", " def fit(self, X, y):\n", " Xb = self._add_bias(X)\n", " n, m = Xb.shape\n", " rng = np.random.RandomState(0)\n", " self.theta = rng.normal(size=m)\n", "\n", " self.mistakes_per_epoch = []\n", " for epoch in range(self.epochs):\n", " mistakes = 0\n", " for xi, yi in zip(Xb, y):\n", " if ...: ### TODO_STUDENT ### detect a mistake: the point is misclassified if y*(theta.x) <= 0\n", " mistakes += 1\n", " self.theta = self.theta + self.... ### TODO_STUDENT ### apply the perceptron mistake-driven update rule\n", " self.mistakes_per_epoch.append(mistakes)\n", " if mistakes == 0:\n", " break\n", " return self\n" ] }, { "cell_type": "code", "execution_count": null, "id": "80f2d151", "metadata": {}, "outputs": [], "source": [ "perceptron = MyPerceptron(lr=1.0, epochs=100).fit(X_perceptron, y_perceptron)\n", "\n", "for epoch, mistakes in enumerate(perceptron.mistakes_per_epoch):\n", " print(f\" epoch {epoch + 1}: mistakes = {mistakes}\")\n", "print(\"Converged after\", len(perceptron.mistakes_per_epoch), \"epochs\")\n", "print(\"Final weight vector (bias first):\", perceptron.theta)\n" ] }, { "cell_type": "markdown", "id": "5a936c3a", "metadata": {}, "source": [ "### Library baseline: sklearn `Perceptron`\n", "\n", "The data is linearly separable, so both perceptrons should reach 0 training errors; we compare how\n", "many epochs each needed, and whether the two separating hyperplanes point the same way.\n" ] }, { "cell_type": "code", "execution_count": null, "id": "5d07424c", "metadata": {}, "outputs": [], "source": [ "from sklearn.linear_model import Perceptron as SkPerceptron\n", "\n", "sk_perceptron = SkPerceptron... ### TODO_STUDENT ### fit sklearn's Perceptron as the library baseline\n", "sk_train_acc = accuracy_score(y_perceptron, sk_perceptron.predict(X_perceptron))\n", "scratch_train_acc = accuracy_score(y_perceptron, perceptron.predict(X_perceptron))\n", "\n", "print(f\"[library] sklearn Perceptron: train acc = {sk_train_acc:.3f}, converged in {sk_perceptron.n_iter_} epochs\")\n", "print(f\"[scratch] our perceptron: train acc = {scratch_train_acc:.3f}, \"\n", " f\"converged in {len(perceptron.mistakes_per_epoch)} epochs\")\n", "\n", "# Both encode the same separating hyperplane w.x + b = 0, just with the bias in a different slot\n", "# (sklearn: [w1..w4, b]; ours: [b, w1..w4]) -- compare orientation via cosine similarity, not raw values.\n", "sk_normal = np.append(sk_perceptron.coef_[0], sk_perceptron.intercept_[0])\n", "scratch_normal = perceptron.theta[[1, 2, 3, 4, 0]]\n", "cos_sim = (sk_normal @ scratch_normal) / (np.linalg.norm(sk_normal) * np.linalg.norm(scratch_normal))\n", "print(f\"Cosine similarity between the two separating hyperplanes: {cos_sim:.3f}\")\n" ] }, { "cell_type": "markdown", "id": "5dc4a286", "metadata": {}, "source": [ "---\n", "## Problem 4 — Gradient Boosting for Classification\n", "\n", "### From-scratch weak learner: the HW1 `DecisionTree`\n", "\n", "Reused verbatim from HW1 Problem 2 (`mode=\"regression\"`, so it fits real-valued residuals rather\n", "than class labels).\n" ] }, { "cell_type": "code", "execution_count": null, "id": "2b015ccf", "metadata": {}, "outputs": [], "source": [ "class TreeNode:\n", " def __init__(self, feature=None, threshold=None, left=None, right=None, value=None):\n", " self.feature = feature\n", " self.threshold = threshold\n", " self.left = left\n", " self.right = right\n", " self.value = value # only set on leaves\n", "\n", " def is_leaf(self):\n", " return self.value is not None\n", "\n", "\n", "def variance(y):\n", " if len(y) == 0:\n", " return 0.0\n", " return np.... ### TODO_STUDENT ### implement the variance (mean squared deviation from the mean) formula\n", "\n", "\n", "def candidate_thresholds(column):\n", " values = np.unique(column)\n", " return ... ### TODO_STUDENT ### compute the midpoints between consecutive sorted unique values as candidate split thresholds\n", "\n", "\n", "class DecisionTree:\n", " \"\"\"Binary regression tree with threshold splits, grown by variance reduction. (Reused from\n", " HW1 Problem 2 -- only the regression mode is needed here.)\"\"\"\n", "\n", " def __init__(self, mode=\"regression\", max_depth=5, min_samples_split=2):\n", " self.mode = mode\n", " self.max_depth = max_depth\n", " self.min_samples_split = min_samples_split\n", "\n", " def fit(self, X, y):\n", " self.root = self._grow(X, y, depth=0)\n", " return self\n", "\n", " def _leaf_value(self, y):\n", " return np.... ### TODO_STUDENT ### return the mean target value among this leaf's points\n", "\n", " def _grow(self, X, y, depth):\n", " n, m = X.shape\n", " if ...: ### TODO_STUDENT ### decide when to stop splitting and turn this node into a leaf\n", " return TreeNode(value=self._leaf_value(y))\n", "\n", " parent_impurity = variance(y)\n", " best_gain, best_feature, best_threshold = 0.0, None, None\n", "\n", " for feature in range(m):\n", " for threshold in candidate_thresholds(X[:, feature]):\n", " left_mask = ... ### TODO_STUDENT ### define which points fall left vs. right of this candidate threshold\n", " right_mask = ~left_mask\n", " if ...: ### TODO_STUDENT ### skip candidate thresholds that would leave one side empty\n", " continue\n", " y_left, y_right = y[left_mask], y[right_mask]\n", " weighted_child_impurity = ... ### TODO_STUDENT ### weight each child's variance by the fraction of the parent's points it received\n", " gain = ... ### TODO_STUDENT ### compute the variance-reduction gain of this candidate split\n", " if gain > best_gain:\n", " best_gain, best_feature, best_threshold = ... ### TODO_STUDENT ### keep track of the best (feature, threshold) split seen so far\n", "\n", " if best_feature is None:\n", " return TreeNode(value=self._leaf_value(y))\n", "\n", " left_mask = ... ### TODO_STUDENT ### recompute the partition using the winning (best_feature, best_threshold) before recursing\n", " left = self._grow... ### TODO_STUDENT ### recursively grow the left subtree on the left split\n", " right = self._grow... ### TODO_STUDENT ### recursively grow the right subtree on the right split\n", " return TreeNode(feature=best_feature, threshold=best_threshold, left=left, right=right)\n", "\n", " def _predict_one(self, x, node):\n", " if node.is_leaf():\n", " return node.value\n", " branch = ... ### TODO_STUDENT ### walk left or right at this node based on the split threshold\n", " return self._predict_one(x, branch)\n", "\n", " def predict(self, X):\n", " return np.array([self._predict_one(x, self.root) for x in X])\n" ] }, { "cell_type": "markdown", "id": "5d026d81", "metadata": {}, "source": [ "### Core: binary gradient boosting (Spambase)\n", "\n", "Maintain an accumulated score F(x); the probability is p(x) = sigmoid(F(x)). Each round, fit a\n", "shallow regression tree to the pseudo-residuals r = y - p(x), and add it to F.\n" ] }, { "cell_type": "code", "execution_count": null, "id": "aa17dbd1", "metadata": {}, "outputs": [], "source": [ "class GradientBoostBinaryClassifier:\n", " def __init__(self, n_rounds=60, tree_depth=2):\n", " self.n_rounds = n_rounds\n", " self.tree_depth = tree_depth\n", "\n", " def sigmoid(self, z):\n", " return ... ### TODO_STUDENT ### implement the sigmoid function\n", "\n", " def fit(self, X, y):\n", " self.trees = []\n", " F = np.zeros(X.shape[0])\n", " for _ in range(self.n_rounds):\n", " residual = ... ### TODO_STUDENT ### compute the pseudo-residual r = y - sigmoid(F(x))\n", " tree = DecisionTree... ### TODO_STUDENT ### fit a shallow regression tree to the current pseudo-residuals\n", " self.trees.append(tree)\n", " F = ... ### TODO_STUDENT ### accumulate this round's tree prediction into the running score F\n", " return self\n", "\n", " def decision_function(self, X):\n", " F = np.zeros(X.shape[0])\n", " for tree in self.trees:\n", " F = ... ### TODO_STUDENT ### accumulate every tree's prediction into the total score F\n", " return F\n", "\n", " def predict_proba(self, X):\n", " \"\"\"p(y=1|x) = sigmoid(F(x)).\"\"\"\n", " return ... ### TODO_STUDENT ### convert the accumulated score F(x) into a probability via sigmoid\n", "\n", " def predict(self, X, threshold=0.5):\n", " return ... ### TODO_STUDENT ### threshold the predicted probability into a 0/1 class label\n", "\n", " def staged_decision_function(self, X):\n", " \"\"\"Yield F(x) after each additional tree (for the T-sweep plot).\"\"\"\n", " F = np.zeros(X.shape[0])\n", " for tree in self.trees:\n", " F = ... ### TODO_STUDENT ### accumulate this round's tree prediction into the running score F\n", " yield F\n" ] }, { "cell_type": "code", "execution_count": null, "id": "5261dda5", "metadata": {}, "outputs": [], "source": [ "n_rounds = 60\n", "gb_spam = GradientBoostBinaryClassifier(n_rounds=n_rounds, tree_depth=2).fit(X_train_spam, y_train_spam)\n", "\n", "train_err, test_err = [], []\n", "for F_train, F_test in zip(gb_spam.staged_decision_function(X_train_spam),\n", " gb_spam.staged_decision_function(X_test_spam)):\n", " train_err.append(1 - accuracy_score(y_train_spam, gb_spam.sigmoid(F_train) >= 0.5))\n", " test_err.append(1 - accuracy_score(y_test_spam, gb_spam.sigmoid(F_test) >= 0.5))\n", "\n", "print(f\"[scratch] Spambase boosting, final (T={n_rounds}): train err = {train_err[-1]:.3f}, test err = {test_err[-1]:.3f}\")\n" ] }, { "cell_type": "markdown", "id": "1b3919e5", "metadata": {}, "source": [ "### Library baseline: sklearn `GradientBoostingClassifier`\n", "\n", "Same tree depth (2) and an unshrunk step (`learning_rate=1.0`, matching our from-scratch version's\n", "full step each round) so the two curves are directly comparable round-for-round.\n" ] }, { "cell_type": "code", "execution_count": null, "id": "0e2fb0d7", "metadata": {}, "outputs": [], "source": [ "from sklearn.ensemble import GradientBoostingClassifier\n", "\n", "sk_gbc_spam = GradientBoostingClassifier(\n", " loss=\"log_loss\", n_estimators=n_rounds, max_depth=2, learning_rate=1.0, random_state=0\n", ")... ### TODO_STUDENT ### fit sklearn's GradientBoostingClassifier(loss='log_loss') as the library baseline\n", "\n", "lib_train_err = [1 - accuracy_score(y_train_spam, pred) for pred in sk_gbc_spam.staged_predict(X_train_spam)]\n", "lib_test_err = [1 - accuracy_score(y_test_spam, pred) for pred in sk_gbc_spam.staged_predict(X_test_spam)]\n", "print(f\"[library] Spambase boosting, final (T={n_rounds}): train err = {lib_train_err[-1]:.3f}, test err = {lib_test_err[-1]:.3f}\")\n", "\n", "plt.figure(figsize=(6, 4))\n", "rounds = range(1, n_rounds + 1)\n", "plt.plot(rounds, lib_train_err, \"b--\", label=\"library train error\")\n", "plt.plot(rounds, lib_test_err, \"r--\", label=\"library test error\")\n", "plt.plot(rounds, train_err, \"b-\", label=\"scratch train error\")\n", "plt.plot(rounds, test_err, \"r-\", label=\"scratch test error\")\n", "plt.xlabel(\"Number of boosting rounds T\")\n", "plt.ylabel(\"Error rate\")\n", "plt.title(\"Gradient boosting for classification: Spambase, library vs. scratch\")\n", "plt.legend()\n", "plt.show()\n" ] }, { "cell_type": "markdown", "id": "4efa5aa5", "metadata": {}, "source": [ "### Stretch: multiclass boosting on 20 Newsgroups (8-class)\n", "\n", "Same pseudo-residual idea generalized to K classes via softmax over K accumulated scores, one\n", "weak learner per class per round. Reuses the 8-class Newsgroups data already loaded in Problem 2c\n", "-- no new fetch, no new training data. The weak learner here is sklearn's\n", "`DecisionTreeRegressor(max_depth=1)` rather than our from-scratch `DecisionTree`: with ~3000 TF-IDF\n", "features, the from-scratch tree's plain-Python double loop over every feature and every candidate\n", "threshold would be far too slow to run `n_estimators` times per class.\n" ] }, { "cell_type": "code", "execution_count": null, "id": "7cce4b7b", "metadata": {}, "outputs": [], "source": [ "from sklearn.tree import DecisionTreeRegressor\n", "\n", "class MyGradientBoostClassifier:\n", " def __init__(self, n_estimators=100, lr=0.5):\n", " self.n_estimators = n_estimators\n", " self.lr = lr\n", "\n", " def fit(self, X, y):\n", " n = X.shape[0]\n", " self.n_classes = len(np.unique(y))\n", " F = np.zeros((n, self.n_classes))\n", " self.weak_learners = [[DecisionTreeRegressor(max_depth=1) for _ in range(self.n_classes)]\n", " for _ in range(self.n_estimators)]\n", "\n", " for t in range(self.n_estimators):\n", " softmax = ... ### TODO_STUDENT ### convert the K accumulated scores into class probabilities via softmax\n", " for k in range(self.n_classes):\n", " residual_k = ... ### TODO_STUDENT ### compute class k's pseudo-residual: (y == k) - p_k\n", " self.weak_learners[t][k]... ### TODO_STUDENT ### fit class k's weak learner to its pseudo-residual\n", " F[:, k] += ... ### TODO_STUDENT ### accumulate this round's prediction into class k's running score\n", " return self\n", "\n", " def decision_function(self, X):\n", " n = X.shape[0]\n", " F = np.zeros((n, self.n_classes))\n", " for t in range(self.n_estimators):\n", " for k in range(self.n_classes):\n", " F[:, k] += ... ### TODO_STUDENT ### accumulate every round's prediction into class k's total score\n", " return F\n", "\n", " def predict(self, X):\n", " F = self.decision_function(X)\n", " return ... ### TODO_STUDENT ### predict the class with the highest accumulated score\n" ] }, { "cell_type": "code", "execution_count": null, "id": "cb60329e", "metadata": {}, "outputs": [], "source": [ "gb_clf_ng = MyGradientBoostClassifier(n_estimators=150, lr=0.6).fit(X_train_ng, y_train_ng)\n", "\n", "acc_train = accuracy_score(y_train_ng, gb_clf_ng.predict(X_train_ng))\n", "acc_test = accuracy_score(y_test_ng, gb_clf_ng.predict(X_test_ng))\n", "print(f\"[stretch] Newsgroups (8-class) boosting: train acc = {acc_train:.3f}, test acc = {acc_test:.3f}\")\n" ] }, { "cell_type": "markdown", "id": "99eb8b21", "metadata": {}, "source": [ "### Library baseline: sklearn `GradientBoostingClassifier` (multiclass)\n", "\n", "`GradientBoostingClassifier` supports multiclass natively (it also fits one weak learner per class\n", "per round internally), so it's a direct library counterpart for this stretch too -- same weak-learner\n", "depth (1) and learning rate (0.6) as our from-scratch version, for a like-for-like comparison.\n" ] }, { "cell_type": "code", "execution_count": null, "id": "49266ec7", "metadata": {}, "outputs": [], "source": [ "sk_gbc_ng = GradientBoostingClassifier(\n", " loss=\"log_loss\", n_estimators=150, max_depth=1, learning_rate=0.6, random_state=0\n", ")... ### TODO_STUDENT ### fit sklearn's GradientBoostingClassifier as the multiclass library baseline\n", "\n", "acc_train_lib = accuracy_score(y_train_ng, sk_gbc_ng.predict(X_train_ng))\n", "acc_test_lib = accuracy_score(y_test_ng, sk_gbc_ng.predict(X_test_ng))\n", "print(f\"[library] Newsgroups (8-class) boosting: train acc = {acc_train_lib:.3f}, test acc = {acc_test_lib:.3f}\")\n" ] }, { "cell_type": "markdown", "id": "a5ca9c86", "metadata": {}, "source": [ "---\n", "## Problem 5 — Active Learning (Uncertainty Sampling)\n", "\n", "One pool-based loop, reused for two base classifiers built above. Both classifiers expose a\n", "`predict_proba` that returns p(y=1|x), so the same \"distance from 0.5\" uncertainty measure works\n", "for both:\n", "\n", "- **(a) gradient-boosted trees** (Problem 4 core) -- uncertainty = |sigmoid(F(x)) - 0.5|\n", "- **(b) logistic regression** (Problem 1) -- uncertainty = |p(y|x) - 0.5|\n", "\n", "Each round we retrain the classifier from scratch on the currently-labeled set, so the boosting\n", "model here uses fewer rounds (15, vs. 60 above) purely to keep the repeated retraining fast --\n", "the active-learning *loop* is the point being tested, not a fully-tuned boosted model.\n" ] }, { "cell_type": "code", "execution_count": null, "id": "7c938e5b", "metadata": {}, "outputs": [], "source": [ "def margin_uncertainty(model, X):\n", " return np.... ### TODO_STUDENT ### define uncertainty as distance from the decision threshold 0.5\n", "\n", "\n", "def active_learning_curve(model_factory, X_pool, y_pool, X_test, y_test,\n", " n_seed=60, n_query=20, n_rounds=15, strategy=\"uncertainty\", seed=0):\n", " rng = np.random.RandomState(seed)\n", " order = rng.permutation(len(X_pool))\n", " labeled = list(order[:n_seed])\n", " unlabeled = list(order[n_seed:])\n", "\n", " n_labeled_list, accs = [], []\n", " for _ in range(n_rounds):\n", " model = model_factory... ### TODO_STUDENT ### retrain the classifier from scratch on the currently-labeled set\n", " accs.append(accuracy_score(y_test, model.predict(X_test)))\n", " n_labeled_list.append(len(labeled))\n", "\n", " if not unlabeled:\n", " break\n", " if strategy == \"uncertainty\":\n", " scores = margin_uncertainty(model, X_pool[unlabeled])\n", " chosen = list... ### TODO_STUDENT ### select the n_query least-certain (smallest-margin) unlabeled points\n", " else:\n", " chosen = list(rng.choice(len(unlabeled), size=min(n_query, len(unlabeled)), replace=False))\n", "\n", " chosen_global = [unlabeled[i] for i in chosen]\n", " labeled.extend(chosen_global)\n", " chosen_set = set(chosen_global)\n", " unlabeled = [i for i in unlabeled if i not in chosen_set]\n", "\n", " return n_labeled_list, accs\n" ] }, { "cell_type": "code", "execution_count": null, "id": "2914f55d", "metadata": {}, "outputs": [], "source": [ "gb_factory = lambda: GradientBoostBinaryClassifier(n_rounds=15, tree_depth=2)\n", "logreg_factory = lambda: MyLogReg(lr=0.1, epochs=3000)\n", "\n", "n_gb_active, acc_gb_active = active_learning_curve(gb_factory, X_train_spam, y_train_spam, X_test_spam, y_test_spam, strategy=\"uncertainty\")\n", "n_gb_random, acc_gb_random = active_learning_curve(gb_factory, X_train_spam, y_train_spam, X_test_spam, y_test_spam, strategy=\"random\")\n", "\n", "n_lr_active, acc_lr_active = active_learning_curve(logreg_factory, X_train_spam, y_train_spam, X_test_spam, y_test_spam, strategy=\"uncertainty\")\n", "n_lr_random, acc_lr_random = active_learning_curve(logreg_factory, X_train_spam, y_train_spam, X_test_spam, y_test_spam, strategy=\"random\")\n", "\n", "fig, axes = plt.subplots(1, 2, figsize=(11, 4))\n", "axes[0].plot(n_gb_active, acc_gb_active, color=\"#2a78d6\", marker=\"o\", markersize=4, label=\"uncertainty sampling\")\n", "axes[0].plot(n_gb_random, acc_gb_random, color=\"#e34948\", marker=\"o\", markersize=4, label=\"random baseline\")\n", "axes[0].set_xlabel(\"Number of labeled points\")\n", "axes[0].set_ylabel(\"Test accuracy\")\n", "axes[0].set_title(\"Active learning: gradient-boosted trees\")\n", "axes[0].legend()\n", "\n", "axes[1].plot(n_lr_active, acc_lr_active, color=\"#2a78d6\", marker=\"o\", markersize=4, label=\"uncertainty sampling\")\n", "axes[1].plot(n_lr_random, acc_lr_random, color=\"#e34948\", marker=\"o\", markersize=4, label=\"random baseline\")\n", "axes[1].set_xlabel(\"Number of labeled points\")\n", "axes[1].set_ylabel(\"Test accuracy\")\n", "axes[1].set_title(\"Active learning: logistic regression\")\n", "axes[1].legend()\n", "\n", "plt.tight_layout()\n", "plt.show()\n" ] }, { "cell_type": "markdown", "id": "414d0eaf", "metadata": {}, "source": [ "### Which classifier benefits more from active selection?\n", "\n", "*(Fill in after running: compare how far above the random baseline each uncertainty-sampling curve\n", "sits, and at what labeling budget the gap is largest.)*\n" ] }, { "cell_type": "markdown", "id": "d04b0928", "metadata": {}, "source": [ "---\n", "## Problem 6 (optional) — Lasso from Scratch (Coordinate Descent)\n", "\n", "L1 isn't differentiable at 0, so Lasso isn't fit by plain gradient descent -- the standard\n", "from-scratch approach is **coordinate descent**: cycle through features j = 1..m, and at each one\n", "compute the partial residual r_j = y - Xtheta + theta_j*x_j (i.e. the residual if feature j's\n", "contribution were removed), its correlation with feature j (rho_j = x_j . r_j), and update\n", "\n", "theta_j <- soft_threshold(rho_j / ||x_j||^2, alpha_eff / ||x_j||^2), soft_threshold(z,t) = sign(z)*max(|z|-t, 0)\n", "\n", "Sweep until the largest coordinate change in a sweep drops below a tolerance. One normalization\n", "detail: sklearn's `Lasso(alpha=...)` minimizes `(1/(2n))*sum-squared-error + alpha*||w||_1`, while\n", "the soft-threshold update above comes from the *un-normalized* `(1/2)*sum-squared-error +\n", "alpha_eff*||w||_1` -- so to match a given sklearn `alpha`, we solve with `alpha_eff = alpha * n`\n", "internally (same normalization fix as Problem 1's Ridge-vs-GD comparison).\n" ] }, { "cell_type": "code", "execution_count": null, "id": "ef0b18e6", "metadata": {}, "outputs": [], "source": [ "class LassoCoordinateDescent:\n", " \"\"\"Lasso fit by cyclic coordinate descent. `alpha` uses sklearn's convention; internally we\n", " solve the equivalent un-normalized problem the soft-threshold update assumes.\"\"\"\n", "\n", " def __init__(self, alpha=1.0, max_sweeps=2000, tol=1e-6):\n", " self.alpha = alpha\n", " self.max_sweeps = max_sweeps\n", " self.tol = tol\n", "\n", " def soft_threshold(self, z, t):\n", " return np.... ### TODO_STUDENT ### implement the soft-thresholding operator, sign(z)*max(|z|-t, 0)\n", "\n", " def fit(self, X, y):\n", " n, m = X.shape\n", " self.y_mean_ = y.mean()\n", " y_centered = y - self.y_mean_ # X is already mean-0 (standardized), so centering y alone\n", " effective_alpha = self.alpha * n # gives an intercept-free fit equivalent to sklearn's\n", " col_norm_sq = np.sum(X ** 2, axis=0)\n", "\n", " theta = np.zeros(m)\n", " pred = X @ theta # kept in sync with theta as we go, so each update is O(n) not O(nm)\n", "\n", " for sweep in range(self.max_sweeps):\n", " max_change = 0.0\n", " for j in range(m):\n", " if col_norm_sq[j] == 0:\n", " continue\n", " r_j = ... ### TODO_STUDENT ### compute the partial residual with feature j's own contribution removed\n", " rho_j = ... ### TODO_STUDENT ### compute feature j's correlation with the partial residual\n", " theta_j_new = self.soft_threshold... ### TODO_STUDENT ### apply the coordinate-descent soft-threshold update for theta_j\n", " delta = theta_j_new - theta[j]\n", " if delta != 0.0:\n", " pred = ... ### TODO_STUDENT ### keep the running prediction in sync with the updated theta_j\n", " max_change = max(max_change, abs(delta))\n", " theta[j] = theta_j_new\n", " if ...: ### TODO_STUDENT ### check for convergence: stop once no coordinate changes by more than the tolerance\n", " break\n", "\n", " self.n_sweeps_ = sweep + 1\n", " self.theta = theta\n", " return self\n", "\n", " def predict(self, X):\n", " return self.y_mean_ + X @ self.theta\n" ] }, { "cell_type": "markdown", "id": "d8c42c34", "metadata": {}, "source": [ "### Verification against sklearn's `Lasso`, at a handful of alphas spanning Problem 2a's sweep" ] }, { "cell_type": "code", "execution_count": null, "id": "7f5cd347", "metadata": {}, "outputs": [], "source": [ "test_alphas = [0.0001, 0.001, 0.01, 0.1, 1.0]\n", "\n", "for a in test_alphas:\n", " sk_lasso = Lasso(alpha=a, max_iter=10000).fit(X_train_spam, y_train_spam)\n", " my_lasso = LassoCoordinateDescent(alpha=a).fit(X_train_spam, y_train_spam)\n", "\n", " coef_diff = np.linalg.norm(sk_lasso.coef_ - my_lasso.theta)\n", " acc_sk = accuracy_score(y_test_spam, sk_lasso.predict(X_test_spam) >= 0.5)\n", " acc_mine = accuracy_score(y_test_spam, my_lasso.predict(X_test_spam) >= 0.5)\n", " nnz_sk = int(np.sum(np.abs(sk_lasso.coef_) > 1e-8))\n", " nnz_mine = int(np.sum(np.abs(my_lasso.theta) > 1e-8))\n", "\n", " print(f\"alpha={a:<8} coef diff (L2 norm) = {coef_diff:.6f}, sweeps = {my_lasso.n_sweeps_:>4} | \"\n", " f\"test acc sklearn={acc_sk:.3f} scratch={acc_mine:.3f} | \"\n", " f\"nonzero coefs sklearn={nnz_sk:>2} scratch={nnz_mine:>2}\")\n" ] } ], "metadata": {}, "nbformat": 4, "nbformat_minor": 5 }