{
 "nbformat": 4,
 "nbformat_minor": 5,
 "metadata": {
  "kernelspec": {
   "name": "python3",
   "display_name": "Python 3",
   "language": "python"
  },
  "language_info": {
   "name": "python"
  },
  "colab": {
   "name": "hilbertspace.ipynb"
  }
 },
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "# Hilbert space \u2014 Python demo\n\nNumerical companion to the entry [Hilbert space](https://dictionaryofml.org/terms/hilbertspace.html) of the [Dictionary of Applied Machine Learning](https://dictionaryofml.org/): it recomputes what the entry states and prints one line per check.\n\nBacks the entry's completeness discussion and its two examples: it generates the Cauchy sequence shown in the entry's figure (converging to a limit that again belongs to the space), and it verifies that the expectation E{x x'} is an inner product on the finite-variance random variables of a common probability space, once random variables that are equal with probability one are identified. Self-contained (numpy only), fixed seeds.\n\nRequires NumPy and Matplotlib only, and uses fixed seeds, so the printed numbers reproduce exactly. Generated from [`pythondemos/hilbertspace.py`](https://dictionaryofml.org/terms/hilbertspace.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(), \"hilbertspace.py\")\nos.makedirs(\"pythondemos\", exist_ok=True)"
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "\"\"\"\nhilbertspace.py \u2014 numerical companion to the glossary entry 'Hilbert space'.\n\nPurpose\n-------\nBacks the entry's completeness discussion and its two examples: it\ngenerates the Cauchy sequence shown in the entry's figure (converging\nto a limit that again belongs to the space), and it verifies that the\nexpectation E{x x'} is an inner product on the finite-variance random\nvariables of a common probability space, once random variables that\nare equal with probability one are identified.  Self-contained (numpy\nonly), fixed seeds.\n\nBlocks\n------\n[B-cauchy]  The sequence w^(t) = w + 0.82^t * 2.3 * (cos(0.65 t + 2.7),\n            sin(0.65 t + 2.7)) in R^2 is a Cauchy sequence: the\n            diameter of the tail {w^(t) : t >= N} shrinks\n            monotonically to 0.  Its limit is w = (-1.2, 1.6), an\n            element of R^2 with finite induced norm \u2014 the limit stays\n            in the space (completeness of the Euclidean space).\n[B-rvspace] On a finite probability space (7 elements, one of\n            them with probability zero), random variables are vectors in\n            R^7 and the inner product <x, x'> = E{x x'} is the\n            probability-weighted inner product of the vectors.  It is\n            symmetric, bilinear and positive semi-definite; two random\n            variables that differ only on the zero-probability element\n            satisfy E{(x - x')^2} = 0 and have zero induced distance \u2014\n            they are identified (equal with probability one), and on\n            the identified space the inner product is positive\n            definite.\n[B-proj]    Orthogonality principle in the Hilbert space of\n            finite-variance random variables: the best linear\n            approximation of y by a x (in the induced norm) has\n            coefficient a* = E{x y} / E{x^2}, and the residual\n            y - a* x is orthogonal to x, E{(y - a* x) x} = 0 \u2014\n            optimal estimation is orthogonal projection.\n\nOutputs\n-------\nhilbertspace_cauchy.csv : points w^(1), ..., w^(14) of the Cauchy\n                          sequence (columns x,y) for the entry's\n                          pgfplots figure; the limit w = (-1.2, 1.6)\n                          is drawn in TikZ directly.\nhilbertspace.png        : matplotlib preview (checking only).\n\"\"\"\n\nimport numpy as np                  # the only numerical dependency\nimport matplotlib\n\nmatplotlib.use(\"Agg\")\nimport matplotlib.pyplot as plt\n\nfrom pathlib import Path\n\nOUT_DIR = Path(__file__).parent\n\nreport = []                         # collects (check name, pass/fail) pairs\n\n\ndef check(name, ok):                # records and prints one verification\n    report.append((name, bool(ok)))\n    print(f\"  [{'ok' if ok else 'FAIL'}] {name}\")"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "**[B-cauchy]** The sequence w^(t) = w + 0.82^t * 2.3 * (cos(0.65 t + 2.7), sin(0.65 t + 2.7)) in R^2 is a Cauchy sequence: the diameter of the tail {w^(t) : t >= N} shrinks monotonically to 0. Its limit is w = (-1.2, 1.6), an element of R^2 with finite induced norm \u2014 the limit stays in the space (completeness of the Euclidean space)."
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "w_lim = np.array([-1.2, 1.6])       # the limit, an element of R^2\nt = np.arange(1, 15)                # sequence indices t = 1, ..., 14\n# a wide, slowly contracting spiral: the early elements sweep through the\n# left half of the figure before closing in on the limit\npts = w_lim + (0.82 ** t)[:, None] * 2.3 * np.c_[\n    np.cos(0.65 * t + 2.7), np.sin(0.65 * t + 2.7)]\n\n# tail diameters diam{w^(t) : t >= N} shrink monotonically to zero\ndiams = [np.max([np.linalg.norm(p - q) for p in pts[N:] for q in pts[N:]])\n         for N in range(0, 10)]\ncheck(\"[B-cauchy]  tail diameters shrink monotonically (Cauchy property)\",\n      all(d1 > d2 for d1, d2 in zip(diams, diams[1:])) and diams[-1] < 0.5)\ncheck(\"[B-cauchy]  the sequence converges to its limit in the space\",\n      np.linalg.norm(pts[-1] - w_lim) < 0.2\n      and np.isfinite(np.sqrt(w_lim @ w_lim)))\n\nwith open(OUT_DIR / \"hilbertspace_cauchy.csv\", \"w\") as fh:\n    fh.write(\"x,y\\n\")               # pgfplots table for the entry's figure\n    for p in pts:\n        fh.write(f\"{p[0]:.4f},{p[1]:.4f}\\n\")"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "**[B-rvspace]** On a finite probability space (7 elements, one of them with probability zero), random variables are vectors in R^7 and the inner product <x, x'> = E{x x'} is the probability-weighted inner product of the vectors. It is symmetric, bilinear and positive semi-definite; two random variables that differ only on the zero-probability element satisfy E{(x - x')^2} = 0 and have zero induced distance \u2014 they are identified (equal with probability one), and on the identified space the inner product is positive definite."
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "rng = np.random.default_rng(7)\nK = 7                               # 7 elements, the last of probability zero\np = np.array([0.10, 0.15, 0.20, 0.25, 0.18, 0.12, 0.00])\n\n\ndef E(x):                           # expectation of a random variable\n    return p @ x\n\n\ndef ip(x, y):                       # the entry's inner product <x, y> = E{x y}\n    return E(x * y)\n\n\nx, y, z = rng.standard_normal((3, K)) * 2.0          # three random variables\na, b = rng.standard_normal(2)                         # random coefficients\n\nok_sym = np.isclose(ip(x, y), ip(y, x))               # symmetry\nok_lin = np.isclose(ip(a * x + b * z, y),             # bilinearity\n                    a * ip(x, y) + b * ip(z, y))\nok_psd = ip(x, x) >= 0                                # positive semi-definite\n# E{x y} = probability-weighted inner product of the two vectors\nok_rep = np.isclose(ip(x, y), np.sum(p * x * y))\ncheck(\"[B-rvspace] E{x y}: symmetric, bilinear, psd, weighted inner \"\n      \"product\", ok_sym and ok_lin and ok_psd and ok_rep)\n\nxp = x.copy()\nxp[-1] += 5.0                       # differs only on the probability-0 element\nok_ident = (np.isclose(ip(x - xp, x - xp), 0.0)       # E{(x - x')^2} = 0\n            and np.isclose(np.sqrt(ip(x - xp, x - xp)), 0.0))\nok_pd = ip(x, x) > 0                # positive definite after identification\ncheck(\"[B-rvspace] x and x' equal with probability one are identified\",\n      ok_ident and ok_pd)"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "**[B-proj]** Orthogonality principle in the Hilbert space of finite-variance random variables: the best linear approximation of y by a x (in the induced norm) has coefficient a* = E{x y} / E{x^2}, and the residual y - a* x is orthogonal to x, E{(y - a* x) x} = 0 \u2014 optimal estimation is orthogonal projection."
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "a_star = ip(x, y) / ip(x, x)        # best linear approximation of y by a x\nresid = y - a_star * x              # estimation error\nok_orth = np.isclose(ip(resid, x), 0.0)               # orthogonality principle\n# a* minimizes the induced norm of the residual over a grid of coefficients\ngrid = a_star + np.linspace(-1.0, 1.0, 201)\nnorms = [ip(y - a * x, y - a * x) for a in grid]\ncheck(\"[B-proj]    residual of the best linear estimator is orthogonal \"\n      \"to x\", ok_orth and np.argmin(norms) == 100)\n\n# ------------------------------------------------------------------ preview\nfig, ax = plt.subplots(figsize=(4.2, 3.4))\nax.plot(pts[:, 0], pts[:, 1], \"k.-\", ms=4, lw=0.6)\nax.annotate(\"\", xy=w_lim, xytext=(0, 0),\n            arrowprops=dict(arrowstyle=\"-|>\", lw=2))\nax.plot(*w_lim, \"k*\", ms=9)\nax.set_aspect(\"equal\")\nax.set_xlabel(\"$w_1$\"); ax.set_ylabel(\"$w_2$\")\nax.set_title(\"Cauchy sequence and its limit in $\\\\mathbb{R}^2$\")\nfig.tight_layout()\nfig.savefig(OUT_DIR / \"hilbertspace.png\", dpi=110)"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "**[B-rkhs]**"
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "# Kernel methods work in an RKHS: the Gaussian kernel's Gram matrix on\n# any point set is symmetric psd (it defines an inner product on the\n# span of kernel sections), and for f = sum_i a_i k(x_i, .) the\n# reproducing property <f, k(x_j, .)> = f(x_j) holds \u2014 the inner product\n# computed via the Gram matrix equals the pointwise evaluation computed\n# directly from the kernel function.\nrng_r = np.random.default_rng(3)\nPk = rng_r.normal(size=(8, 2))\nkfun = lambda u, v: np.exp(-np.sum((u - v) ** 2) / 2.0)\nK = np.array([[kfun(Pk[i], Pk[j]) for j in range(8)] for i in range(8)])\ncheck(\"[B-rkhs]   Gaussian Gram matrix is symmetric psd (an inner \"\n      \"product on kernel sections)\",\n      np.allclose(K, K.T) and np.linalg.eigvalsh(K).min() > -1e-10)\na_coef = rng_r.normal(size=8)\ninner_via_gram = K @ a_coef                # <f, k(x_j,.)> = (K a)_j\nf_pointwise = np.array([sum(a_coef[i] * kfun(Pk[i], Pk[j])\n                            for i in range(8)) for j in range(8)])\ncheck(\"[B-rkhs]   reproducing property <f, k(x_j,.)> = f(x_j)\",\n      np.allclose(inner_via_gram, f_pointwise))\n\nn_ok = sum(ok for _, ok in report)\nprint(f\"\\n{n_ok}/{len(report)} checks pass\")\nprint(\"wrote hilbertspace_cauchy.csv, pythondemos/hilbertspace.png\")\nif n_ok != len(report):\n    raise SystemExit(1)"
  }
 ]
}