{
 "nbformat": 4,
 "nbformat_minor": 5,
 "metadata": {
  "kernelspec": {
   "name": "python3",
   "display_name": "Python 3",
   "language": "python"
  },
  "language_info": {
   "name": "python"
  },
  "colab": {
   "name": "convex.ipynb"
  }
 },
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "# convex \u2014 Python demo\n\nNumerical companion to the entry [convex](https://dictionaryofml.org/terms/convex.html) of the [Dictionary of Applied Machine Learning](https://dictionaryofml.org/): it recomputes what the entry states and prints one line per check.\n\nVerifies the entry's two definitions and its ML claim on concrete objects: segment membership for a convex set, the chord inequality and epigraph convexity for the average squared error loss of linear regression, and the local-equals-global-minimum property via gradient descent from many random initializations. Self-contained (numpy only), fixed seed.\n\nRequires NumPy and Matplotlib only, and uses fixed seeds, so the printed numbers reproduce exactly. Generated from [`pythondemos/convex.py`](https://dictionaryofml.org/terms/convex.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(), \"convex.py\")\nos.makedirs(\"pythondemos\", exist_ok=True)"
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "\"\"\"\nconvex.py \u2014 numerical companion to the glossary entry 'convex'.\n\nPurpose\n-------\nVerifies the entry's two definitions and its ML claim on concrete objects:\nsegment membership for a convex set, the chord inequality and epigraph\nconvexity for the average squared error loss of linear regression, and the\nlocal-equals-global-minimum property via gradient descent from many random\ninitializations.  Self-contained (numpy only), fixed seed.\n\nBlocks\n------\n[B-set]    The ellipse set C = {w : w^T diag(1, 4) w <= 4} contains\n           beta w + (1-beta) w' for 500 random pairs w, w' in C and\n           21 values beta in [0, 1].\n[B-fn]     The linear-regression objective f(w) = (1/m)||y - X w||^2\n           (m = 10, d = 2, seed 0) satisfies the chord inequality\n           f(beta w + (1-beta) w') <= beta f(w) + (1-beta) f(w')\n           for 500 random pairs and 21 values of beta.\n[B-epi]    Convex combinations of 500 random point pairs of the epigraph\n           {(w, t) : t >= f(w)} stay in the epigraph.\n[B-global] Gradient descent on f from 20 random initializations always\n           reaches the same minimizer (pairwise distance < 1e-6) with the\n           least-squares optimum as its objective value: every local\n           minimum is global.\n[B-half]   Halfspace representation: C is contained in the intersection\n           P_k of k supporting halfspaces, and the area of P_k minus C\n           shrinks as k grows (k = 8, 32, 128): the intersection of all\n           supporting halfspaces recovers C.\n\nOutputs\n-------\nconvex.png : matplotlib preview \u2014 1-D slice of f with a chord above the\n             graph (checking only; the entry's figure is schematic TikZ).\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\nrng = np.random.default_rng(0)\nbetas = np.linspace(0.0, 1.0, 21)"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "**[B-set]** The ellipse set C = {w : w^T diag(1, 4) w <= 4} contains beta w + (1-beta) w' for 500 random pairs w, w' in C and 21 values beta in [0, 1]."
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "D = np.diag([1.0, 4.0])\n\n\ndef in_set(w):\n    return np.einsum(\"...i,ij,...j->...\", w, D, w) <= 4.0 + 1e-12\n\n\npts = rng.uniform(-2.0, 2.0, (5000, 2))\npts = pts[in_set(pts)][:1000]\npairs = pts.reshape(-1, 2, 2)[:500]\nok_set = all(in_set(b * p[0] + (1 - b) * p[1])\n             for p in pairs for b in betas)\ncheck(\"[B-set]    segments between set points stay in the set\", ok_set)"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "**[B-fn]** The linear-regression objective f(w) = (1/m)||y - X w||^2 (m = 10, d = 2, seed 0) satisfies the chord inequality f(beta w + (1-beta) w') <= beta f(w) + (1-beta) f(w') for 500 random pairs and 21 values of beta."
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "m, d = 10, 2\nX = rng.standard_normal((m, d))\ny = X @ np.array([1.0, -2.0]) + 0.1 * rng.standard_normal(m)\n\n\ndef f(w):\n    r = y - X @ w\n    return float(r @ r) / m\n\n\nW = rng.standard_normal((500, 2, 2)) * 3.0\nok_fn = all(f(b * w1 + (1 - b) * w2) <= b * f(w1) + (1 - b) * f(w2) + 1e-10\n            for w1, w2 in W for b in betas)\ncheck(\"[B-fn]     chord inequality for the linreg objective\", ok_fn)"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "**[B-epi]** Convex combinations of 500 random point pairs of the epigraph {(w, t) : t >= f(w)} stay in the epigraph."
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "ok_epi = True\nfor w1, w2 in W:\n    t1 = f(w1) + abs(rng.standard_normal())      # points above the graph\n    t2 = f(w2) + abs(rng.standard_normal())\n    for b in betas:\n        wb, tb = b * w1 + (1 - b) * w2, b * t1 + (1 - b) * t2\n        if tb < f(wb) - 1e-10:\n            ok_epi = False\ncheck(\"[B-epi]    epigraph is a convex set\", ok_epi)"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "**[B-global]** Gradient descent on f from 20 random initializations always reaches the same minimizer (pairwise distance < 1e-6) with the least-squares optimum as its objective value: every local minimum is global."
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "L = 2.0 * np.linalg.eigvalsh(X.T @ X / m).max()\neta = 1.0 / L\nminimizers = []\nfor _ in range(20):\n    w = rng.standard_normal(2) * 5.0\n    for _ in range(2000):\n        w = w - eta * (2.0 / m) * X.T @ (X @ w - y)\n    minimizers.append(w)\nminimizers = np.array(minimizers)\nspread = np.max(np.linalg.norm(minimizers - minimizers[0], axis=1))\nw_star, *_ = np.linalg.lstsq(X, y, rcond=None)\nok_glob = spread < 1e-6 and abs(f(minimizers[0]) - f(w_star)) < 1e-10\ncheck(\"[B-global] GD from 20 random inits reaches the global minimum\",\n      ok_glob)"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "**[B-half]** Halfspace representation: C is contained in the intersection P_k of k supporting halfspaces, and the area of P_k minus C shrinks as k grows (k = 8, 32, 128): the intersection of all supporting halfspaces recovers C."
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "# Supporting halfspace of C = {w : w^T D w <= 4} at boundary point p:\n# normal n = D p, halfspace {w : n^T w <= n^T p}.\ndef polygon_gap_area(k, n_mc=200000):\n    \"\"\"Monte-Carlo area of P_k \\\\ C for the intersection P_k of k\n    supporting halfspaces at equally spaced boundary points.\"\"\"\n    th = np.linspace(0.0, 2.0 * np.pi, k, endpoint=False)\n    bnd = np.c_[2.0 * np.cos(th), np.sin(th)]          # boundary of C\n    N = bnd @ D                                        # outward normals\n    c = np.einsum(\"ki,ki->k\", N, bnd)\n    S = rng_mc.uniform(-2.4, 2.4, (n_mc, 2))\n    in_pk = np.all(S @ N.T <= c + 1e-12, axis=1)\n    in_c = in_set(S)\n    box = 4.8 * 4.8\n    contain_ok = not np.any(in_c & ~in_pk)             # C subset of P_k\n    return box * np.mean(in_pk & ~in_c), contain_ok\n\n\nrng_mc = np.random.default_rng(1)\ngaps, contains = zip(*(polygon_gap_area(k) for k in (8, 32, 128)))\ncheck(\"[B-half]   C inside every halfspace intersection; gap area shrinks\",\n      all(contains) and gaps[0] > gaps[1] > gaps[2] and gaps[2] < 0.02)"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "**[B-sep]**"
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "# Linear separability of a binary trainset holds exactly when the convex\n# hulls of the two classes do not intersect. Both directions are checked:\n# hulls disjoint -> a separating (w, b) exists; hulls overlapping -> none\n# does. Hull intersection is decided by the LP feasibility of\n#   sum_i a_i x_i^+ = sum_j c_j x_j^-,  a, c >= 0,  sum a = sum c = 1,\n# solved here by projected gradient on the squared distance between the\n# two hulls (zero distance <=> the hulls meet).\ndef hull_distance(P, N, iters=5000):\n    \"\"\"Distance between conv(P) and conv(N) by Frank-Wolfe with exact line\n    search on z = p - n, which stays in the (convex) Minkowski difference:\n    each step moves z toward the vertex p_i - n_j that the current z sees as\n    steepest, so ||z|| decreases monotonically to the true distance.\"\"\"\n    z = P[0] - N[0]\n    for _ in range(iters):\n        # vertex of conv(P) - conv(N) minimizing the linear approximation\n        s_vert = P[np.argmin(P @ z)] - N[np.argmax(N @ z)]\n        diff = z - s_vert\n        denom = diff @ diff\n        if denom < 1e-18:\n            break\n        gam = min(1.0, max(0.0, (z @ diff) / denom))   # exact line search\n        z = z - gam * diff\n    return float(np.linalg.norm(z))\n\n\ndef separable(P, N, iters=20000, lr=0.1):\n    \"\"\"Perceptron-style search for (w, b) with y (w^T x + b) > 0.\"\"\"\n    X = np.vstack([P, N])\n    y = np.r_[np.ones(len(P)), -np.ones(len(N))]\n    w, b = np.zeros(X.shape[1]), 0.0\n    for _ in range(iters):\n        m = y * (X @ w + b) <= 0\n        if not m.any():\n            return True\n        w += lr * (y[m] @ X[m])\n        b += lr * y[m].sum()\n    return bool(np.all(y * (X @ w + b) > 0))\n\n\nrng_sep = np.random.default_rng(3)\nP_far = rng_sep.normal(size=(12, 2)) * 0.3 + np.array([2.5, 0.0])\nN_far = rng_sep.normal(size=(12, 2)) * 0.3 + np.array([-2.5, 0.0])\nP_mix = rng_sep.normal(size=(12, 2)) * 0.9\nN_mix = rng_sep.normal(size=(12, 2)) * 0.9\ncheck(f\"[B-sep]    disjoint hulls (distance {hull_distance(P_far, N_far):.3f}) \"\n      f\"-> the trainset is linearly separable\",\n      hull_distance(P_far, N_far) > 1e-3 and separable(P_far, N_far))\ncheck(f\"[B-sep]    overlapping hulls (distance {hull_distance(P_mix, N_mix):.3f}) \"\n      f\"-> no linear classifier separates the trainset\",\n      hull_distance(P_mix, N_mix) < 1e-3 and not separable(P_mix, N_mix))"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "**[B-expfam]**"
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "# Exponential family p(x; w) ~ h(x) exp(w^T T(x) - A(w)) on x in {0,1,2}\n# with T(x) = (x, x^2), h = 1. Two claims of the entry: the log-partition\n# function A is convex in w, and its gradient is E{T(x)}, the average of\n# the sufficient statistics under p(x; w).\nTstat = np.c_[np.arange(3), np.arange(3) ** 2].astype(float)\n\n\ndef logpart(w):\n    return float(np.log(np.exp(Tstat @ w).sum()))\n\n\nrng_ef = np.random.default_rng(5)\nviol = 0\nfor _ in range(20000):                       # chord test for convexity of A\n    u, v = rng_ef.normal(size=2) * 1.5, rng_ef.normal(size=2) * 1.5\n    t = rng_ef.uniform()\n    if logpart(t * u + (1 - t) * v) > t * logpart(u) + (1 - t) * logpart(v) + 1e-12:\n        viol += 1\ncheck(f\"[B-expfam] the log-partition function A is convex \"\n      f\"({viol} chord violations in 20000 random pairs)\", viol == 0)\n\nw_ef = np.array([0.3, -0.2])\nh_ef = 1e-6\ngradA = np.array([(logpart(w_ef + h_ef * e) - logpart(w_ef - h_ef * e)) / (2 * h_ef)\n                  for e in np.eye(2)])\npw = np.exp(Tstat @ w_ef)\npw /= pw.sum()\ncheck(f\"[B-expfam] grad A equals the expected sufficient statistics \"\n      f\"({gradA.round(4)} vs {(pw @ Tstat).round(4)})\",\n      np.allclose(gradA, pw @ Tstat, atol=1e-5))\n\n# -------------------------------------------------------------- preview\nw0, w1v = np.array([-2.0, -3.0]), np.array([3.0, 1.0])\nts = np.linspace(0.0, 1.0, 100)\nfig, ax = plt.subplots(figsize=(4.6, 3.0))\nax.plot(ts, [f(t * w1v + (1 - t) * w0) for t in ts], \"k-\",\n        label=\"$f$ along the segment\")\nax.plot([0, 1], [f(w0), f(w1v)], \"k--\", label=\"chord\")\nax.set_xlabel(r\"$\\beta$\"), ax.set_ylabel(\"objective value\")\nax.set_title(\"[B-fn] the chord lies above the function\")\nax.legend(frameon=False, fontsize=8)\nfig.tight_layout()\nfig.savefig(OUT_DIR / \"convex.png\", dpi=110)\n\nn_ok = sum(ok for _, ok in report)\nprint(f\"\\n{n_ok}/{len(report)} checks pass\")\nprint(\"wrote convex.png\")\nif n_ok != len(report):\n    raise SystemExit(1)"
  }
 ]
}