{
 "nbformat": 4,
 "nbformat_minor": 5,
 "metadata": {
  "kernelspec": {
   "name": "python3",
   "display_name": "Python 3",
   "language": "python"
  },
  "language_info": {
   "name": "python"
  },
  "colab": {
   "name": "svm.ipynb"
  }
 },
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "# support vector machine (SVM) \u2014 Python demo\n\nNumerical companion to the entry [support vector machine (SVM)](https://dictionaryofml.org/terms/svm.html) of the [Dictionary of Applied Machine Learning](https://dictionaryofml.org/): it recomputes what the entry states and prints one line per check.\n\nA small, hand-crafted non-separable toy trainset in R^2 whose soft-margin SVM solution leaves exactly one data point misclassified *beyond* the opposite margin, i.e. with functional margin y (w^T x + b) < -1 (a bounded support vector, hinge loss > 2). The same dataset drives the entry's non-separable figure (Fig. 2). Self-contained (numpy/matplotlib only), fixed seed; the SVM is solved by a pure-numpy simplified SMO so no external solver is needed.\n\nRequires NumPy and Matplotlib only, and uses fixed seeds, so the printed numbers reproduce exactly. Generated from [`pythondemos/svm.py`](https://dictionaryofml.org/terms/svm.py); CC BY 4.0."
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "# Notebook shim: the script resolves output paths relative to __file__,\n# which a notebook kernel does not define; everything lands in the\n# working directory instead.\nimport os\n__file__ = os.path.join(os.getcwd(), \"svm.py\")\nos.makedirs(\"pythondemos\", exist_ok=True)"
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "\"\"\"\nsvm.py \u2014 numerical companion to the glossary entry\n'support vector machine (SVM)'.\n\nPurpose\n-------\nA small, hand-crafted non-separable toy trainset in R^2 whose soft-margin\nSVM solution leaves exactly one data point misclassified *beyond* the\nopposite margin, i.e. with functional margin y (w^T x + b) < -1 (a bounded\nsupport vector, hinge loss > 2). The same dataset drives the entry's\nnon-separable figure (Fig. 2). Self-contained (numpy/matplotlib only),\nfixed seed; the SVM is solved by a pure-numpy simplified SMO so no external\nsolver is needed.\n\nBlocks\n------\n[P-nonsep] Two symmetric sets of points (class +1 right, class -1 left) plus one\n           outlier labelled +1 planted deep in the -1 region. The soft-margin\n           SVM keeps the vertical boundary between the two classes (correcting the\n           outlier would misclassify the whole -1 class) and yields\n           w_hat = (0.5, 0), b_hat = 0: boundary x1 = 0, margins x1 = +-2,\n           margin 1/||w_hat|| = 2.\n[P-sv]     The four points at x1 = +-2 are margin support vectors\n           (0 < alpha < C, functional margin exactly 1); the outlier is a\n           bounded support vector (alpha = C); the two points at x1 = +-3 are\n           not support vectors (alpha = 0) and can be removed without changing\n           w_hat.\n[P-alpha]   Overly large regularization: ||w|| <= rho/(2 alpha),\n           constant majority predictions, minority class all\n           misclassified support vectors.\n[P-viol]   The outlier x = (-4, 0), y = +1 has y (w^T x + b) = -2 < -1, hinge\n           loss xi = 3.\n\nOutputs\n-------\nsvm_points.csv : the 7 data points with columns x1,x2,label,cls,fmargin\n                 (cls in {pos,neg,svpos,svneg,viol}) read by the entry's\n                 pgfplots figure.\nsvm.png        : matplotlib preview of that figure (checking only).\n\"\"\"\n\nimport numpy as np\nimport matplotlib\n\nmatplotlib.use(\"Agg\")\nimport matplotlib.pyplot as plt\n\nfrom pathlib import Path\n\nOUT_DIR = Path(__file__).parent\n\nreport = []\n\n\ndef check(name, ok):\n    report.append((name, bool(ok)))\n    print(f\"  [{'ok' if ok else 'FAIL'}] {name}\")\n\n\n# --------------------------------------------------------------------------\n# Pure-numpy soft-margin SVM (simplified SMO, Platt 1998 / CS229 notes).\n# Solves the C-SVM dual   max sum_i a_i - 1/2 sum_ij a_i a_j y_i y_j <x_i,x_j>\n# s.t. 0 <= a_i <= C, sum_i a_i y_i = 0.  Output: w = sum_i a_i y_i x_i, b, a.\n# --------------------------------------------------------------------------\ndef smo(X, y, C, tol=1e-6, max_passes=500, seed=0):\n    rng = np.random.default_rng(seed)\n    m = X.shape[0]\n    K = X @ X.T\n    a = np.zeros(m)\n    b = 0.0\n    passes = 0\n\n    def score(i):\n        return (a * y) @ K[i] + b\n\n    while passes < max_passes:\n        changed = 0\n        for i in range(m):\n            Ei = score(i) - y[i]\n            if (y[i] * Ei < -tol and a[i] < C) or (y[i] * Ei > tol and a[i] > 0):\n                j = int(rng.integers(m - 1))\n                j = j + (j >= i)                       # uniform j != i\n                Ej = score(j) - y[j]\n                ai, aj = a[i], a[j]\n                if y[i] != y[j]:\n                    L, H = max(0.0, aj - ai), min(C, C + aj - ai)\n                else:\n                    L, H = max(0.0, ai + aj - C), min(C, ai + aj)\n                if L == H:\n                    continue\n                eta = 2 * K[i, j] - K[i, i] - K[j, j]\n                if eta >= 0:\n                    continue\n                a[j] = np.clip(aj - y[j] * (Ei - Ej) / eta, L, H)\n                if abs(a[j] - aj) < 1e-12:\n                    continue\n                a[i] = ai + y[i] * y[j] * (aj - a[j])\n                b1 = b - Ei - y[i] * (a[i] - ai) * K[i, i] - y[j] * (a[j] - aj) * K[i, j]\n                b2 = b - Ej - y[i] * (a[i] - ai) * K[i, j] - y[j] * (a[j] - aj) * K[j, j]\n                b = b1 if 0 < a[i] < C else b2 if 0 < a[j] < C else 0.5 * (b1 + b2)\n                changed += 1\n        passes = passes + 1 if changed == 0 else 0\n    w = (a * y) @ X\n    return w, b, a\n\n\n# --------------------------------------------------------------------------\n# Toy dataset in R^2\n# --------------------------------------------------------------------------\nX = np.array([\n    [ 2.0,  1.0], [ 2.0, -1.0], [ 3.0, 0.0],   # class +1 points (right)\n    [-2.0,  1.0], [-2.0, -1.0], [-3.0, 0.0],   # class -1 points (left)\n    [-4.0,  0.0],                              # outlier labelled +1, deep in -1 region\n])\ny = np.array([1.0, 1.0, 1.0, -1.0, -1.0, -1.0, 1.0])\nOUT = 6                                         # index of the outlier\nC = 1.0                                         # C = 1/(2 m alpha); m=7 -> alpha = 1/14\n\nw, b, alpha = smo(X, y, C)\nfmargin = y * (X @ w + b)                        # functional margins y_i (w^T x_i + b)\nxi = np.maximum(0.0, 1.0 - fmargin)              # hinge losses\n\nprint(f\"w_hat = {np.round(w,4).tolist()}, b_hat = {round(b,4)}, \"\n      f\"margin 1/||w|| = {1/np.linalg.norm(w):.3f}\")\nprint(f\"functional margins: {np.round(fmargin,3).tolist()}\")\nprint(f\"alpha (0<alpha<C: margin SV; alpha=C: bounded SV): {np.round(alpha,3).tolist()}\")"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "**[P-nonsep]** Two symmetric sets of points (class +1 right, class -1 left) plus one outlier labelled +1 planted deep in the -1 region. The soft-margin SVM keeps the vertical boundary between the two classes (correcting the outlier would misclassify the whole -1 class) and yields w_hat = (0.5, 0), b_hat = 0: boundary x1 = 0, margins x1 = +-2, margin 1/||w_hat|| = 2."
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "check(\"w_hat = (0.5, 0), b_hat = 0 (vertical boundary x1 = 0)\",\n      np.allclose(w, [0.5, 0.0], atol=1e-3) and abs(b) < 1e-3)\ncheck(\"margin 1/||w_hat|| = 2\", abs(1 / np.linalg.norm(w) - 2.0) < 1e-3)"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "**[P-sv]** The four points at x1 = +-2 are margin support vectors (0 < alpha < C, functional margin exactly 1); the outlier is a bounded support vector (alpha = C); the two points at x1 = +-3 are not support vectors (alpha = 0) and can be removed without changing w_hat."
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "margin_sv = np.array([0, 1, 3, 4])\nnon_sv = np.array([2, 5])\ncheck(\"the four points at x1 = +-2 are margin support vectors (0 < alpha < C)\",\n      np.all((alpha[margin_sv] > 1e-6) & (alpha[margin_sv] < C - 1e-6)))\ncheck(\"the two points at x1 = +-3 are not support vectors (alpha = 0)\",\n      np.all(alpha[non_sv] < 1e-6))\ncheck(\"the outlier is a bounded support vector (alpha = C)\", abs(alpha[OUT] - C) < 1e-6)\n# removing a non-support vector leaves w_hat unchanged\nkeep = [i for i in range(len(y)) if i != non_sv[0]]\nw2, b2, _ = smo(X[keep], y[keep], C)\ncheck(\"removing a non-support vector leaves w_hat, b_hat unchanged\",\n      np.allclose(w2, w, atol=1e-3) and abs(b2 - b) < 1e-3)"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "**[P-viol]** The outlier x = (-4, 0), y = +1 has y (w^T x + b) = -2 < -1, hinge loss xi = 3."
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "check(\"outlier x=(-4,0), y=+1 has y(w^T x + b) < -1 (bounded, hinge loss > 2)\",\n      fmargin[OUT] < -1 and xi[OUT] > 2)\nprint(f\"outlier: y(w^T x + b) = {fmargin[OUT]:+.3f}, hinge loss xi = {xi[OUT]:.3f}\")\n\n# --------------------------------------------------------------------------\n# CSV for the entry's pgfplots figure\n# --------------------------------------------------------------------------\ncls = np.array([\"pos\", \"pos\", \"pos\", \"neg\", \"neg\", \"neg\", \"viol\"], dtype=object)\ncls[margin_sv[:2]] = \"svpos\"        # +1 margin SVs  (2, +-1)\ncls[margin_sv[2:]] = \"svneg\"        # -1 margin SVs  (-2, +-1)\nwith open(OUT_DIR / \"svm_points.csv\", \"w\") as fh:\n    fh.write(\"x1,x2,label,cls,fmargin\\n\")\n    for (x1, x2), yi, ci, fm in zip(X, y, cls, fmargin):\n        fh.write(f\"{x1:g},{x2:g},{int(yi):+d},{ci},{fm:g}\\n\")\nprint(\"wrote svm_points.csv\")\n\n# --------------------------------------------------------------------------\n# matplotlib preview (checking only)\n# --------------------------------------------------------------------------\nfig, ax = plt.subplots(figsize=(6.4, 3.2))\nfor xv, lab in [(0, r\"$0$\"), (-2, r\"$-1$\"), (2, r\"$+1$\")]:\n    ax.axvline(xv, ls=\"-\" if xv == 0 else \"--\", color=\"0.4\", lw=1.2)\n    ax.text(xv, 1.55, lab, ha=\"center\", fontsize=9)\nax.text(0, 1.9, r\"$\\hat w^\\top x + \\hat b$\", ha=\"center\", fontsize=9)\npos = y > 0\nax.scatter(X[pos & (cls != \"viol\") & (alpha < C), 0],\n           X[pos & (cls != \"viol\") & (alpha < C), 1], marker=\"o\", c=\"tab:blue\", s=40, label=\"+1\")\nax.scatter(X[~pos, 0], X[~pos, 1], marker=\"s\", facecolors=\"none\", edgecolors=\"k\", s=40, label=\"-1\")\nsv = (alpha > 1e-6) & (alpha < C - 1e-6)\nax.scatter(X[sv, 0], X[sv, 1], marker=\"o\", facecolors=\"none\", edgecolors=\"r\", s=180, lw=2)\nax.scatter(X[OUT, 0], X[OUT, 1], marker=\"o\", c=\"tab:blue\", edgecolors=\"r\", s=120, lw=2)\n# xi of the outlier: gap to its own (+1) margin line at x1 = +2\nax.annotate(\"\", xy=(2, 0), xytext=(-4, 0),\n            arrowprops=dict(arrowstyle=\"<->\", color=\"r\", lw=1))\nax.text(-1, 0.15, r\"$\\xi = 3$\", color=\"r\", ha=\"center\", fontsize=9)\nax.text(-4, -0.35, \"misclassified\\noutlier\", color=\"r\", ha=\"center\", fontsize=8)\nax.text(3, 0.3, \"not a\\nsupport vector\", ha=\"center\", fontsize=8)\nax.set_xlim(-4.8, 3.8)\nax.set_ylim(-1.6, 2.1)\nax.set_xlabel(r\"$x_1$\")\nax.set_ylabel(r\"$x_2$\")\nax.set_title(\"Data generated by pythondemos/svm.py\")\nfig.tight_layout()\n# --------------------------------------------------------------------------"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "**[P-alpha]** Overly large regularization: ||w|| <= rho/(2 alpha), constant majority predictions, minority class all misclassified support vectors."
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "# underfit: ||w_hat|| <= rho/(2 alpha) (rho = largest feature norm), the\n# predictions collapse to the constant majority class sign(b_hat), and\n# every minority-class data point becomes a misclassified support vector\n# (sparsity of the expansion is lost). Note: with the unpenalized offset\n# b, NOT every data point becomes a support vector (majority points can\n# sit at margin slightly above one) \u2014 only the homogeneous b = 0 SVM has\n# an all-support-vector threshold.\n# --------------------------------------------------------------------------\nrho = np.max(np.linalg.norm(X, axis=1))\nnorms_alpha = []\nfor alpha_reg in (1.0, 10.0, 50.0):\n    C_a = 1.0 / (2 * len(y) * alpha_reg)\n    w_a, b_a, _ = smo(X, y, C_a, max_passes=2000)\n    norms_alpha.append(np.linalg.norm(w_a))\n    check(f\"[P-alpha] ||w_hat|| <= rho/(2 alpha) at alpha = {alpha_reg:g}\",\n          np.linalg.norm(w_a) <= rho / (2 * alpha_reg) + 1e-9)\ncheck(\"[P-alpha] ||w_hat|| shrinks as alpha grows\",\n      norms_alpha[0] > norms_alpha[1] > norms_alpha[2])\nmarg_a = y * (X @ w_a + b_a)                      # alpha = 50 solution\nminority = y == -1                                 # 3 of 7 labels\ncheck(\"[P-alpha] predictions collapse to the constant majority class\",\n      np.all(np.sign(X @ w_a + b_a) == 1.0))\ncheck(\"[P-alpha] every minority-class point is a misclassified \"\n      \"support vector\", np.all(marg_a[minority] < 0))\ncheck(\"[P-alpha] sparsity lost: more support vectors than at the \"\n      \"default C\", (marg_a <= 1 + 1e-6).sum() > 5)\n\n# --------------------------------------------------------------------------"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "**[P-rerm]** the primal RERM view: (1/m) sum_i hinge_i + lam ||w||^2 with"
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "# lam = 1/(2 C m) is equivalent to the dual that SMO solves \u2014 subgradient\n# descent on this nonsmooth convex objective (Pegasos-style step 1/(lam t))\n# reaches the same solution.\n# --------------------------------------------------------------------------\nm_svm = len(y)\nlam = 1.0 / (2.0 * C * m_svm)\ndef primal(wv, bv):\n    return (np.mean(np.maximum(0.0, 1.0 - y * (X @ wv + bv)))\n            + lam * wv @ wv)\nwp = np.zeros(2); bp = 0.0\nfor t in range(1, 20001):\n    act = y * (X @ wp + bp) < 1.0                    # margin violators\n    gw = 2.0 * lam * wp - (y[act, None] * X[act]).sum(0) / m_svm\n    gb = -y[act].sum() / m_svm\n    step = 1.0 / (lam * t)\n    wp -= step * gw; bp -= step * gb\ncheck(\"[P-rerm] subgradient descent on the primal reaches the dual \"\n      \"(SMO) objective value\",\n      abs(primal(wp, bp) - primal(w, b)) < 5e-3)\ncheck(\"[P-rerm] primal and dual parameter vectors agree\",\n      np.linalg.norm(wp - w) < 0.05 and abs(bp - b) < 0.1)\n\n# --------------------------------------------------------------------------"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "**[P-kernel]** kernel extension via new features: the polynomial features"
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "# phi(x) = (x1^2, sqrt(2) x1 x2, x2^2) linearize a circular pattern\n# that no linear classifier on the raw features separates, and its Gram\n# matrix equals the polynomial kernel (x^T x')^2 \u2014 the kernel trick.\n# --------------------------------------------------------------------------\nrng_k = np.random.default_rng(7)\nr_in = rng_k.uniform(0.0, 0.8, 20); a_in = rng_k.uniform(0, 2*np.pi, 20)\nr_out = rng_k.uniform(1.4, 2.0, 20); a_out = rng_k.uniform(0, 2*np.pi, 20)\nXc = np.vstack([np.c_[r_in*np.cos(a_in), r_in*np.sin(a_in)],\n                np.c_[r_out*np.cos(a_out), r_out*np.sin(a_out)]])\nyc = np.concatenate([np.ones(20), -np.ones(20)])\nphi = lambda Z: np.c_[Z[:, 0]**2, np.sqrt(2)*Z[:, 0]*Z[:, 1], Z[:, 1]**2]\nw_lin, b_lin, _ = smo(Xc, yc, 10.0)\nacc_lin = np.mean(np.sign(Xc @ w_lin + b_lin) == yc)\nw_phi, b_phi, _ = smo(phi(Xc), yc, 10.0)\nacc_phi = np.mean(np.sign(phi(Xc) @ w_phi + b_phi) == yc)\nprint(f\"[P-kernel] accuracy raw features {acc_lin:.2f} vs polynomial \"\n      f\"features {acc_phi:.2f}\")\ncheck(\"[P-kernel] no linear classifier separates the circles \"\n      \"(raw accuracy well below 1)\", acc_lin < 0.8)\ncheck(\"[P-kernel] the same SVM on mapped features separates them\",\n      acc_phi == 1.0)\ncheck(\"[P-kernel] Gram matrix of phi equals the polynomial kernel \"\n      \"(x^T x')^2 (kernel trick)\",\n      np.allclose(phi(Xc) @ phi(Xc).T, (Xc @ Xc.T) ** 2, atol=1e-8))\n\nfig.savefig(OUT_DIR / \"svm.png\", dpi=110)\nprint(\"wrote svm.png\")\n\nn_ok = sum(ok for _, ok in report)\nprint(f\"\\n{n_ok}/{len(report)} checks passed\")\nassert n_ok == len(report), \"some checks FAILED\""
  }
 ]
}