{
 "nbformat": 4,
 "nbformat_minor": 5,
 "metadata": {
  "kernelspec": {
   "name": "python3",
   "display_name": "Python 3",
   "language": "python"
  },
  "language_info": {
   "name": "python"
  },
  "colab": {
   "name": "eigenvalue.ipynb"
  }
 },
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "# eigenvalue \u2014 Python demo\n\nNumerical companion to the entry [eigenvalue](https://dictionaryofml.org/terms/eigenvalue.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/eigenvalue.py`](https://dictionaryofml.org/terms/eigenvalue.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(), \"eigenvalue.py\")\nos.makedirs(\"pythondemos\", exist_ok=True)"
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "\"\"\"\neigenvalue.py \u2014 numerical companion to the glossary entry 'eigenvalue'.\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]   lambda is an eigenvalue of a square matrix A iff A x = lambda x\n          for some nonzero vector x: every (lambda, x) pair returned by\n          np.linalg.eig satisfies the defining equation, the eigenvectors\n          are nonzero, and applying A to an eigenvector only rescales it\n          (direction preserved \u2014 the content of the entry's figure).\n[P-conv]  eigenvalues decide convergence of iterative methods that\n          repeatedly apply an affine update w -> M w + b, here GD for\n          linear regression with M = I - (2 eta/m) X^T X: the error\n          norm is bounded by rho^t times the initial error, with rho\n          the largest eigenvalue magnitude of M; a stepsize below\n          m/lambda_max(X^T X) gives rho < 1 and convergence, a larger\n          one gives rho > 1 and a growing error; each eigenvalue of M\n          is the image of an eigenvalue of X^T X under\n          lambda -> 1 - (2 eta/m) lambda.\n[P-graph] the second-smallest eigenvalue lambda_2 of a graph Laplacian\n          measures connectivity, growing as edges are added to the same\n          six nodes (the entry's figure): lambda_1 = 0 always;\n          lambda_2 = 0 for two separate clusters, lambda_2 ~ 0.44 once\n          a connecting edge is added, and lambda_2 = 6 for the complete\n          graph; the signs of the entries of the eigenvector for\n          lambda_2 (the Fiedler vector) recover the two clusters \u2014\n          spectral clustering.\n\nOutputs\n-------\neigenvalue.png : preview figure (checking only).\n\nData generated by pythondemos/eigenvalue.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}\")"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "**[P-def]** lambda is an eigenvalue of a square matrix A iff A x = lambda x for some nonzero vector x: every (lambda, x) pair returned by np.linalg.eig satisfies the defining equation, the eigenvectors are nonzero, and applying A to an eigenvector only rescales it (direction preserved \u2014 the content of the entry's figure)."
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "print(\"[P-def] A x = lambda x for every eigenpair\")\nA = np.array([[2.0, 1.0], [1.0, 3.0]])         # symmetric -> real eigenvalues\nlam, V = np.linalg.eig(A)\nfor i in range(2):\n    x = V[:, i]\n    check(f\"pair {i}: ||A x - lambda x|| < 1e-12 \"\n          f\"(lambda = {lam[i]:.4f})\",\n          np.linalg.norm(A @ x - lam[i] * x) < 1e-12)\n    check(f\"pair {i}: eigenvector is nonzero\", np.linalg.norm(x) > 0)\n    # direction preserved: A x is collinear with x\n    cos = abs(x @ (A @ x)) / (np.linalg.norm(x) * np.linalg.norm(A @ x))\n    check(f\"pair {i}: A x collinear with x (|cos| = 1)\",\n          abs(cos - 1) < 1e-12)\n# a non-eigenvector is NOT mapped to a multiple of itself\nu = np.array([1.0, 0.0])\ncos_u = abs(u @ (A @ u)) / (np.linalg.norm(u) * np.linalg.norm(A @ u))\ncheck(\"generic vector changes direction under A (|cos| < 1)\",\n      cos_u < 1 - 1e-6)"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "**[P-conv]** eigenvalues decide convergence of iterative methods that repeatedly apply an affine update w -> M w + b, here GD for linear regression with M = I - (2 eta/m) X^T X: the error norm is bounded by rho^t times the initial error, with rho the largest eigenvalue magnitude of M; a stepsize below m/lambda_max(X^T X) gives rho < 1 and convergence, a larger one gives rho > 1 and a growing error; each eigenvalue of M is the image of an eigenvalue of X^T X under lambda -> 1 - (2 eta/m) lambda."
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "print(\"[P-conv] largest eigenvalue magnitude of the update matrix \"\n      \"decides convergence\")\nm_tr, d = 40, 2\nX = rng.standard_normal((m_tr, d)) * np.array([1.0, 3.0])\ny_tr = X @ np.array([1.0, -2.0]) + 0.1 * rng.standard_normal(m_tr)\nQ = X.T @ X\nlam_Q = np.linalg.eigvalsh(Q)\nw_hat = np.linalg.solve(Q, X.T @ y_tr)\n\n\ndef run_gd(eta, n_iter):\n    M = np.eye(d) - (2 * eta / m_tr) * Q       # affine update w -> M w + b\n    b = (2 * eta / m_tr) * (X.T @ y_tr)\n    rho = np.max(np.abs(np.linalg.eigvalsh(M)))\n    w = np.zeros(d)\n    errs = np.empty(n_iter)\n    for t in range(n_iter):\n        errs[t] = np.linalg.norm(w - w_hat)\n        w = M @ w + b\n    return rho, errs\n\n\neta_good = m_tr / (lam_Q[-1] + lam_Q[0])       # below m/lambda_max\nrho_g, err_g = run_gd(eta_good, 300)\ncheck(\"stepsize below m/lambda_max: rho < 1\", rho_g < 1)\ncheck(\"error norm bounded by rho^t times the initial error\",\n      np.all(err_g[:20] <= rho_g ** np.arange(20) * err_g[0] * (1 + 1e-9)))\ncheck(\"error converges to 0\", err_g[-1] < 1e-8 * err_g[0])\neta_bad = 1.05 * m_tr / lam_Q[-1]              # above m/lambda_max\nrho_b, err_b = run_gd(eta_bad, 300)\ncheck(\"stepsize above m/lambda_max: rho > 1 and the error grows\",\n      rho_b > 1 and err_b[-1] > err_b[0])\nlam_M = np.linalg.eigvalsh(np.eye(d) - (2 * eta_good / m_tr) * Q)\ncheck(\"eigenvalues of M are 1 - (2 eta/m) lambda for eigenvalues \"\n      \"lambda of X^T X\",\n      np.allclose(np.sort(lam_M),\n                  np.sort(1 - (2 * eta_good / m_tr) * lam_Q)))"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "**[P-graph]** the second-smallest eigenvalue lambda_2 of a graph Laplacian measures connectivity, growing as edges are added to the same six nodes (the entry's figure): lambda_1 = 0 always; lambda_2 = 0 for two separate clusters, lambda_2 ~ 0.44 once a connecting edge is added, and lambda_2 = 6 for the complete graph; the signs of the entries of the eigenvector for lambda_2 (the Fiedler vector) recover the two clusters \u2014 spectral clustering."
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "print(\"[P-graph] second-smallest Laplacian eigenvalue measures connectivity\")\n\n\ndef laplacian(edge_list, n):\n    L = np.zeros((n, n))\n    for i, j in edge_list:\n        L[i, j] -= 1.0\n        L[j, i] -= 1.0\n        L[i, i] += 1.0\n        L[j, j] += 1.0\n    return L\n\n\nn_nodes = 6\n# two clusters of three nodes each, every pair inside a cluster linked\ntwo_clusters = [(0, 1), (1, 2), (0, 2), (3, 4), (4, 5), (3, 5)]\nwith_bridge = two_clusters + [(2, 3)]          # one edge between the clusters\ncomplete = [(i, j) for i in range(n_nodes)     # every pair of nodes linked\n            for j in range(i + 1, n_nodes)]\nlam_disc = np.linalg.eigvalsh(laplacian(two_clusters, n_nodes))\nlam_conn, V_conn = np.linalg.eigh(laplacian(with_bridge, n_nodes))\nlam_comp = np.linalg.eigvalsh(laplacian(complete, n_nodes))\ncheck(\"smallest eigenvalue lambda_1 = 0 for all three graphs\",\n      abs(lam_disc[0]) < 1e-12 and abs(lam_conn[0]) < 1e-12\n      and abs(lam_comp[0]) < 1e-12)\ncheck(\"disconnected graph (no edge between clusters): lambda_2 = 0\",\n      abs(lam_disc[1]) < 1e-12)\ncheck(\"connected graph (one edge between clusters): lambda_2 ~ 0.44\",\n      lam_conn[1] > 1e-8 and abs(lam_conn[1] - 0.4384) < 1e-3)\ncheck(\"complete graph: lambda_2 = 6, the largest possible value\",\n      np.isclose(lam_comp[1], n_nodes))\nfiedler = V_conn[:, 1]\nside = fiedler > 0\ncheck(\"signs of the Fiedler vector recover the two clusters\",\n      len(set(side[:3])) == 1 and len(set(side[3:])) == 1\n      and side[0] != side[3])\n\n# ------------------------------------------------------------ preview\nfig, ax = plt.subplots(1, 3, figsize=(12.6, 4.0))\nfor i, c in zip(range(2), (\"C0\", \"C1\")):\n    x = V[:, i]\n    ax[0].arrow(0, 0, *x, head_width=0.06, color=c,\n                length_includes_head=True)\n    ax[0].arrow(0, 0, *(A @ x), head_width=0.06, color=c, alpha=0.4,\n                length_includes_head=True)\n    ax[0].annotate(f\"$\\\\lambda_{i+1}={lam[i]:.2f}$\", xy=A @ x)\nax[0].arrow(0, 0, *u, head_width=0.06, color=\"k\", length_includes_head=True)\nax[0].arrow(0, 0, *(A @ u), head_width=0.06, color=\"k\", alpha=0.35,\n            length_includes_head=True)\nax[0].set_aspect(\"equal\")\nax[0].set_xlabel(\"$x_1$\")\nax[0].set_ylabel(\"$x_2$\")\nax[0].set_title(\"[P-def] eigenvectors keep direction\")\nt_show = np.arange(40)\nax[1].semilogy(t_show, err_g[:40], \"o-\", ms=3,\n               label=f\"$\\\\rho = {rho_g:.2f} < 1$\")\nax[1].semilogy(t_show, rho_g ** t_show * err_g[0], \"k--\",\n               label=\"$\\\\rho^t \\\\cdot$ initial error\")\nax[1].semilogy(t_show, err_b[:40], \"s:\", ms=3,\n               label=f\"$\\\\rho = {rho_b:.2f} > 1$\")\nax[1].set_xlabel(\"iteration $t$\")\nax[1].set_ylabel(\"error norm\")\nax[1].set_title(\"[P-conv] update-matrix eigenvalues\")\nax[1].legend(frameon=False)\nidx = np.arange(1, n_nodes + 1)\nax[2].bar(idx - 0.26, lam_disc, width=0.24, color=\"C0\", hatch=\"//\",\n          label=\"two clusters\")\nax[2].bar(idx, lam_conn, width=0.24, color=\"C1\",\n          label=\"one connecting edge\")\nax[2].bar(idx + 0.26, lam_comp, width=0.24, color=\"C2\", hatch=\"xx\",\n          label=\"complete graph\")\nax[2].set_xlabel(\"eigenvalue index $i$\")\nax[2].set_ylabel(\"Laplacian eigenvalue $\\\\lambda_i$\")\nax[2].set_title(\"[P-graph] $\\\\lambda_2$ grows with connectivity\")\nax[2].legend(frameon=False)\nfig.tight_layout()\nfig.savefig(OUT_DIR / \"eigenvalue.png\", dpi=110)\nprint(f\"\\n{sum(ok for _, ok in report)}/{len(report)} checks passed\")\nassert all(ok for _, ok in report)"
  }
 ]
}