{
 "nbformat": 4,
 "nbformat_minor": 5,
 "metadata": {
  "kernelspec": {
   "name": "python3",
   "display_name": "Python 3",
   "language": "python"
  },
  "language_info": {
   "name": "python"
  },
  "colab": {
   "name": "projection.ipynb"
  }
 },
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "# projection \u2014 Python demo\n\nNumerical companion to the entry [projection](https://dictionaryofml.org/terms/projection.html) of the [Dictionary of Applied Machine Learning](https://dictionaryofml.org/): it recomputes what the entry states and prints one line per check.\n\nOne block per paragraph of the entry (marked [P...]): each block verifies numerically what the corresponding statement asserts. Self-contained (numpy/matplotlib only), fixed seed.\n\nRequires NumPy and Matplotlib only, and uses fixed seeds, so the printed numbers reproduce exactly. Generated from [`pythondemos/projection.py`](https://dictionaryofml.org/terms/projection.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(), \"projection.py\")\nos.makedirs(\"pythondemos\", exist_ok=True)"
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "\"\"\"\nprojection.py \u2014 numerical companion to the glossary entry 'projection'.\n\nOne block per paragraph of the entry (marked [P...]): each block verifies\nnumerically what the corresponding statement asserts. Self-contained\n(numpy/matplotlib only), fixed seed.\n\nBlocks\n------\n[P-def]    The projection onto a closed set is a closest point: the exact\n           l1-ball projection (via sorting-based soft threshold) beats a\n           dense sample of other points of the set in Euclidean\n           distance; for the convex l1-ball it is unique, and the\n           segment w - proj(w) meets the face of the ball at a right\n           angle (the marker in the entry's figure). On the two-point\n           set {(-1,0), (1,0)} the closest point is not unique for\n           (0, 0.7) \u2014 nor for any other point of the vertical axis,\n           which consists of exactly the points equidistant from the\n           two elements. For a subspace, the projection map is linear\n           (orthogonal projection matrix P = B (B^T B)^{-1} B^T with\n           P^2 = P), while the l1-ball is not a subspace ((tau, 0) is\n           in the ball, 2 (tau, 0) is not) and its projection violates\n           additivity \u2014 it is not linear (the entry's closing remark).\n[P-fund]   Idempotence and self-adjointness: every projection is\n           idempotent (also the nonlinear l1-ball projection); the\n           subspace projection matrix satisfies P^2 = P and P^T = P;\n           an oblique projection (idempotent, not self-adjoint) sends\n           a vector to a point of its range that is farther away than\n           the orthogonal projection foot.\n[P-projgd] Projected GD for Lasso: gradient steps on the training error\n           followed by l1-ball projections converge to a feasible\n           iterate whose training error is (near-)optimal among\n           feasible points, enforcing ||w||_1 <= tau in every\n           iteration. The projection step is what enforces the\n           constraint: the same gradient steps without it produce\n           iterates with ||w||_1 > tau.\n\nOutputs\n-------\nprojection.png : preview figure (checking only).\n\nData generated by pythondemos/projection.py.\n\"\"\"\n\nimport numpy as np\nimport matplotlib\n\nmatplotlib.use(\"Agg\")\nimport matplotlib.pyplot as plt\nfrom pathlib import Path\n\nOUT_DIR = Path(__file__).parent\n\nrng = np.random.default_rng(42)\nreport = []\n\n\ndef check(name, ok):\n    report.append((name, bool(ok)))\n    print(f\"  [{'ok' if ok else 'FAIL'}] {name}\")\n\n\ndef proj_l1(w, tau):\n    \"\"\"Exact Euclidean projection onto the l1-ball of radius tau.\"\"\"\n    if np.abs(w).sum() <= tau:\n        return w.copy()\n    u = np.sort(np.abs(w))[::-1]\n    css = np.cumsum(u)\n    rho = np.nonzero(u * np.arange(1, len(w) + 1) > css - tau)[0][-1]\n    theta = (css[rho] - tau) / (rho + 1.0)\n    return np.sign(w) * np.maximum(np.abs(w) - theta, 0.0)"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "**[P-def]** The projection onto a closed set is a closest point: the exact l1-ball projection (via sorting-based soft threshold) beats a dense sample of other points of the set in Euclidean distance; for the convex l1-ball it is unique, and the segment w - proj(w) meets the face of the ball at a right angle (the marker in the entry's figure). On the two-point set {(-1,0), (1,0)} the closest point is not unique for (0, 0.7) \u2014 nor for any other point of the vertical axis, which consists of exactly the points equidistant from the two elements. For a subspace, the projection map is linear (orthogonal projection matrix P = B (B^T B)^{-1} B^T with P^2 = P), while the l1-ball is not a subspace ((tau, 0) is in the ball, 2 (tau, 0) is not) and its projection violates additivity \u2014 it is not linear (the entry's closing remark)."
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "print(\"[P-def] projection = closest point; linear only on subspaces\")\ntau = 1.0\nw0 = np.array([1.6, 1.1])\np = proj_l1(w0, tau)\ncheck(\"projection is feasible (||p||_1 <= tau)\",\n      np.abs(p).sum() <= tau + 1e-12)\n# closest point: sample the l1-ball densely, none is closer\nangles = rng.uniform(0, 2 * np.pi, 20000)\nradii = rng.uniform(0, 1, 20000)\nraw = np.stack([np.cos(angles), np.sin(angles)], axis=1)\nball = radii[:, None] * raw / np.abs(raw).sum(axis=1, keepdims=True)\ndists = np.linalg.norm(ball - w0, axis=1)\ncheck(\"no sampled point of the ball is closer than the projection\",\n      np.all(dists >= np.linalg.norm(p - w0) - 1e-9))\n# right angle (figure marker): w - proj(w) is orthogonal to the face\n# x_1 + x_2 = tau of the ball that contains the projection\nface = np.array([1.0, -1.0])\ncheck(\"w - proj(w) is orthogonal to the face containing proj(w)\",\n      abs((w0 - p) @ face) < 1e-9)\n# subspace: orthogonal projection matrix, linear and idempotent\nB = rng.normal(size=(4, 2))                    # columns span a 2-d subspace\nP = B @ np.linalg.solve(B.T @ B, B.T)\nu1, u2 = rng.normal(size=4), rng.normal(size=4)\ncheck(\"subspace projection is linear: P(u + u') = P u + P u'\",\n      np.allclose(P @ (u1 + u2), P @ u1 + P @ u2))\ncheck(\"idempotent: P^2 = P\", np.allclose(P @ P, P))\ncheck(\"residual orthogonal to the subspace: B^T (u - P u) = 0\",\n      np.max(np.abs(B.T @ (u1 - P @ u1))) < 1e-10)\n# uniqueness on the convex ball: every near-minimizer is near p\nnear = ball[dists <= np.linalg.norm(p - w0) + 1e-3]\ncheck(\"convex set: every near-closest point lies near the projection\",\n      np.all(np.linalg.norm(near - p, axis=1) < 0.15))\n# existence on a nonconvex closed set (two points): minimum attained,\n# but NOT unique for the midpoint\nS = np.array([[1.0, 0.0], [-1.0, 0.0]])\nmid = np.array([0.0, 0.7])\nd_mid = np.linalg.norm(S - mid, axis=1)\ncheck(\"nonconvex closed set: a closest point exists but is not unique\",\n      np.isclose(d_mid[0], d_mid[1]))\n# ... and likewise for every point of the vertical axis, which consists\n# of exactly the points equidistant from the two elements\naxis_pts = np.stack([np.zeros(5), np.linspace(-2.0, 2.0, 5)], axis=1)\ncheck(\"every point of the vertical axis is equidistant from both elements\",\n      np.allclose(np.linalg.norm(axis_pts - S[0], axis=1),\n                  np.linalg.norm(axis_pts - S[1], axis=1)))\noff_axis = np.array([0.3, 0.7])\ncheck(\"a point off the vertical axis has a unique closest element\",\n      not np.isclose(np.linalg.norm(off_axis - S[0]),\n                     np.linalg.norm(off_axis - S[1])))\n# the l1-ball is not a subspace ((tau,0) is in the ball, 2(tau,0) is\n# not), and its projection is NOT linear\nedge = np.array([tau, 0.0])\ncheck(\"l1-ball is not a subspace: (tau,0) inside, 2(tau,0) outside\",\n      np.abs(edge).sum() <= tau and np.abs(2 * edge).sum() > tau)\na1, a2 = np.array([1.5, 0.0]), np.array([0.0, 1.5])\ncheck(\"l1-ball projection violates additivity (not a linear map)\",\n      not np.allclose(proj_l1(a1 + a2, tau),\n                      proj_l1(a1, tau) + proj_l1(a2, tau)))"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "**[P-fund]** Idempotence and self-adjointness: every projection is idempotent (also the nonlinear l1-ball projection); the subspace projection matrix satisfies P^2 = P and P^T = P; an oblique projection (idempotent, not self-adjoint) sends a vector to a point of its range that is farther away than the orthogonal projection foot."
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "print(\"[P-fund] idempotent + self-adjoint characterizes orthogonal projection\")\ncheck(\"l1-ball projection is idempotent: proj(proj(w)) = proj(w)\",\n      np.allclose(proj_l1(proj_l1(w0, tau), tau), proj_l1(w0, tau)))\ncheck(\"subspace projection matrix: P^2 = P and P^T = P\",\n      np.allclose(P @ P, P) and np.allclose(P.T, P))\n# near miss: an oblique projection is idempotent but not self-adjoint,\n# and its output is not a closest point of its range (the x-axis)\nA = np.array([[1.0, 1.0], [0.0, 0.0]])\nu = np.array([0.3, 0.8])\nfoot = np.array([u[0], 0.0])\ncheck(\"oblique projection: idempotent but not self-adjoint\",\n      np.allclose(A @ A, A) and not np.allclose(A.T, A))\ncheck(\"oblique output is farther from u than the orthogonal foot\",\n      np.linalg.norm(A @ u - u) > np.linalg.norm(foot - u) + 1e-9)"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "**[P-projgd]** Projected GD for Lasso: gradient steps on the training error followed by l1-ball projections converge to a feasible iterate whose training error is (near-)optimal among feasible points, enforcing ||w||_1 <= tau in every iteration. The projection step is what enforces the constraint: the same gradient steps without it produce iterates with ||w||_1 > tau."
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "print(\"[P-projgd] projected GD solves the Lasso constraint form\")\nm, d = 40, 2\nX = rng.normal(size=(m, d))\ny = X @ np.array([1.2, -0.3]) + 0.05 * rng.normal(size=m)\nL = 2 * np.linalg.eigvalsh(X.T @ X / m).max()\nw = np.zeros(d)\nfeasible_all = True\nfor _ in range(400):\n    w = proj_l1(w - (1 / L) * (2 / m) * X.T @ (X @ w - y), tau)\n    feasible_all &= np.abs(w).sum() <= tau + 1e-10\ntrainerr = lambda v: np.mean((y - X @ v) ** 2)\n# compare against a dense sample of the feasible set\ncand = tau * ball / np.maximum(np.abs(ball).sum(axis=1, keepdims=True), 1e-12)\ntrainerrs = np.mean((y[None, :] - cand @ X.T) ** 2, axis=1)\ncheck(\"every iterate stayed feasible\", feasible_all)\n# the cost of omitting the projection: plain GD violates the constraint\nw_plain = np.zeros(d)\nfor _ in range(400):\n    w_plain = w_plain - (1 / L) * (2 / m) * X.T @ (X @ w_plain - y)\ncheck(\"without the projection step, iterates violate ||w||_1 <= tau\",\n      np.abs(w_plain).sum() > tau + 1e-10)\ncheck(\"projected-GD training error <= best sampled feasible one + 1e-3\",\n      trainerr(w) <= trainerrs.min() + 1e-3)\n\n# ------------------------------------------------------------ preview\nfig, ax = plt.subplots(figsize=(4.2, 4.0))\nsquare = np.array([[1, 0], [0, 1], [-1, 0], [0, -1], [1, 0]]) * tau\nax.plot(square[:, 0], square[:, 1], \"k-\")\nax.plot(*w0, \"ko\"); ax.annotate(\"w\", w0)\nax.plot(*p, \"rs\"); ax.annotate(\"proj(w)\", p)\nax.plot([w0[0], p[0]], [w0[1], p[1]], \"k--\")\nax.plot(*w, \"b^\"); ax.annotate(\"projected GD\", w)\nax.set_xlabel(\"$w_1$\"); ax.set_ylabel(\"$w_2$\")\nax.set_aspect(\"equal\"); ax.set_title(\"[P-def] projection onto the l1-ball\")\nfig.tight_layout()\nfig.savefig(OUT_DIR / \"projection.png\", dpi=110)\nprint(f\"\\n{sum(ok for _, ok in report)}/{len(report)} checks passed\")\nassert all(ok for _, ok in report)"
  }
 ]
}