{
 "nbformat": 4,
 "nbformat_minor": 5,
 "metadata": {
  "kernelspec": {
   "name": "python3",
   "display_name": "Python 3",
   "language": "python"
  },
  "language_info": {
   "name": "python"
  },
  "colab": {
   "name": "covmtx.ipynb"
  }
 },
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "# covariance matrix \u2014 Python demo\n\nNumerical companion to the entry [covariance matrix](https://dictionaryofml.org/terms/covmtx.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/covmtx.py`](https://dictionaryofml.org/terms/covmtx.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(), \"covmtx.py\")\nos.makedirs(\"pythondemos\", exist_ok=True)"
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "\"\"\"\ncovmtx.py \u2014 numerical companion to the glossary entry 'covariance matrix'.\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 covariance matrix C = E{(x - E x)(x - E x)^T}: the empirical\n         outer-product average of iid draws from a random vector with\n         analytic covariance A A^T recovers that matrix; its entry (j, j')\n         is the covariance of entries j and j' (checked against np.cov),\n         and its diagonal holds the per-entry variances.\n[P-mvn]  The basic ML picture under a Gaussian model of (feature, label):\n         the principal axes of the constant-density ellipse are the\n         eigenvectors of the 2 x 2 covariance matrix \u2014 the draws'\n         coordinates along them are uncorrelated with variances equal to\n         the eigenvalues \u2014 and the smallest-risk hypothesis map is\n         linear with slope cov(x, y)/var(x).\n[P-prec] The diagonal of the inverse covariance matrix (the precision\n         matrix): 1/(C^{-1})_{jj} equals the conditional variance of\n         entry j given the remaining entries (the variance of the error\n         of the best linear prediction of entry j from the rest), and on\n         iid Gaussian draws the Bayes estimator E{x_j | rest} attains\n         this value as its squared-error risk \u2014 the baseline that no\n         prediction can beat; a constant prediction does worse.\n\nOutputs\n-------\ncovmtx.png : preview figure (checking only).\n\nData generated by pythondemos/covmtx.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]** The covariance matrix C = E{(x - E x)(x - E x)^T}: the empirical outer-product average of iid draws from a random vector with analytic covariance A A^T recovers that matrix; its entry (j, j') is the covariance of entries j and j' (checked against np.cov), and its diagonal holds the per-entry variances."
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "print(\"[P-def] C = E{(x - Ex)(x - Ex)^T} \u2014 empirical vs analytic\")\nA = np.array([[1.0, 0.0, 0.0], [0.5, 0.8, 0.0], [-0.2, 0.3, 1.1]])\nmu = np.array([2.0, -1.0, 0.5])\nC = A @ A.T                                    # analytic covariance\nm = 10**6\nx = rng.standard_normal((m, 3)) @ A.T + mu\nxc = x - x.mean(axis=0)                        # centered\nC_emp = (xc[:, :, None] * xc[:, None, :]).mean(axis=0)\ncheck(\"empirical outer-product average recovers C (|err| < 5e-3)\",\n      np.max(np.abs(C_emp - C)) < 5e-3)\ncheck(\"matches np.cov (ddof=0)\",\n      np.allclose(C_emp, np.cov(x.T, ddof=0), atol=1e-9))\ncovs = np.array([[np.mean(xc[:, j] * xc[:, k]) for k in range(3)]\n                 for j in range(3)])\ncheck(\"entry (j, j') is the covariance of entries j and j'\",\n      np.allclose(covs, C_emp, atol=1e-12))\ncheck(\"diagonal entries are the per-entry variances\",\n      np.allclose(np.diag(C_emp),\n                  [np.mean(xc[:, j] ** 2) for j in range(3)], atol=1e-12))"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "**[P-mvn]** The basic ML picture under a Gaussian model of (feature, label): the principal axes of the constant-density ellipse are the eigenvectors of the 2 x 2 covariance matrix \u2014 the draws' coordinates along them are uncorrelated with variances equal to the eigenvalues \u2014 and the smallest-risk hypothesis map is linear with slope cov(x, y)/var(x)."
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "print(\"[P-mvn] EVD of the covariance matrix gives the principal axes \"\n      \"of the Gaussian density\")\nC2 = np.array([[1.0, 0.8], [0.8, 1.0]])        # (feature, label) covariance\nmu2 = np.array([2.5, 3.0])\nlam2, V2 = np.linalg.eigh(C2)\ncheck(\"eigenvalues of the 2 x 2 covariance matrix are 0.2 and 1.8\",\n      np.allclose(lam2, [0.2, 1.8]))\nz = rng.standard_normal((10**6, 2)) @ np.linalg.cholesky(C2).T + mu2\nproj = (z - z.mean(axis=0)) @ V2               # coordinates on the axes\ncheck(\"coordinates along the eigenvectors are uncorrelated\",\n      abs(np.mean(proj[:, 0] * proj[:, 1])) < 5e-3)\ncheck(\"variance along each principal axis equals the eigenvalue\",\n      np.allclose(proj.var(axis=0), lam2, atol=5e-3))\nslope_hat = np.polyfit(z[:, 0], z[:, 1], 1)[0]\ncheck(\"smallest-risk linear slope matches cov(x, y)/var(x) = 0.8\",\n      abs(slope_hat - C2[0, 1] / C2[0, 0]) < 5e-3)"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "**[P-prec]** The diagonal of the inverse covariance matrix (the precision matrix): 1/(C^{-1})_{jj} equals the conditional variance of entry j given the remaining entries (the variance of the error of the best linear prediction of entry j from the rest), and on iid Gaussian draws the Bayes estimator E{x_j | rest} attains this value as its squared-error risk \u2014 the baseline that no prediction can beat; a constant prediction does worse."
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "print(\"[P-prec] diagonal of C^{-1}: conditional variance and \"\n      \"squared-error baseline\")\nP = np.linalg.inv(C)\nmse_bayes = np.empty(3)\nfor j in range(3):\n    rest = [k for k in range(3) if k != j]\n    # conditional variance of entry j given the rest = variance of the\n    # error of the best linear prediction of entry j from the rest\n    w = np.linalg.solve(C[np.ix_(rest, rest)], C[rest, j])\n    cond_var = C[j, j] - C[rest, j] @ w\n    check(f\"entry {j}: 1/(C^-1)_jj equals the conditional variance\",\n          np.isclose(1 / P[j, j], cond_var))\n    # Bayes estimator E{x_j | rest} is linear in the Gaussian case; its\n    # squared-error risk on iid Gaussian draws attains the baseline\n    pred = mu[j] + (x[:, rest] - mu[rest]) @ w\n    mse_bayes[j] = np.mean((x[:, j] - pred) ** 2)\n    check(f\"entry {j}: risk of the Bayes estimator matches 1/(C^-1)_jj\",\n          abs(mse_bayes[j] - 1 / P[j, j]) < 1e-2)\nmse_const = np.mean((x[:, 0] - mu[0]) ** 2)    # constant prediction\ncheck(\"constant prediction of entry 0 has larger risk than the baseline\",\n      mse_const > 1 / P[0, 0])\n\n# ------------------------------------------------------------ preview\nfig, ax = plt.subplots(1, 4, figsize=(13.6, 3.0))\nfor a, M, t in ((ax[0], C, \"analytic C = A A^T\"),\n                (ax[1], C_emp, \"empirical (m = 1e6)\")):\n    im = a.imshow(M, cmap=\"gray\")\n    a.set_title(t)\n    fig.colorbar(im, ax=a, shrink=0.8)\nzs = z[:400]\nax[2].plot(zs[:, 0], zs[:, 1], \".\", color=\"0.6\", ms=2)\nang = np.linspace(0, 2 * np.pi, 100)\nell = mu2[:, None] + V2 @ (np.sqrt(lam2)[:, None]\n                           * np.vstack([np.cos(ang), np.sin(ang)]))\nax[2].plot(ell[0], ell[1], \"k-\", lw=1.5)\nfor lam_i, v_i in zip(lam2, V2.T):\n    ax[2].annotate(\"\", xy=mu2 + np.sqrt(lam_i) * v_i, xytext=mu2,\n                   arrowprops=dict(arrowstyle=\"->\", lw=1.5))\nxg = np.linspace(z[:, 0].min(), z[:, 0].max(), 2)\nax[2].plot(xg, C2[0, 1] / C2[0, 0] * (xg - mu2[0]) + mu2[1], \"--\",\n           color=\"C1\", label=\"hypothesis map\")\nax[2].set_xlabel(\"feature x\")\nax[2].set_ylabel(\"label y\")\nax[2].set_title(\"[P-mvn] principal axes from the EVD\")\nax[2].legend(frameon=False)\njj = np.arange(3)\nax[3].bar(jj - 0.18, np.diag(C), width=0.36, color=\"C0\", hatch=\"//\",\n          label=\"variance $C_{j,j}$\")\nax[3].bar(jj + 0.18, 1 / np.diag(P), width=0.36, color=\"C1\",\n          label=\"baseline $1/(C^{-1})_{j,j}$\")\nax[3].set_xlabel(\"entry j\")\nax[3].set_ylabel(\"squared-error risk\")\nax[3].set_title(\"[P-prec] baseline from $C^{-1}$\")\nax[3].legend(frameon=False)\nfig.suptitle(\"covariance matrix, Gaussian principal axes, and the \"\n             \"$C^{-1}$ baseline\")\nfig.tight_layout()\nfig.savefig(OUT_DIR / \"covmtx.png\", dpi=110)\nprint(f\"\\n{sum(ok for _, ok in report)}/{len(report)} checks passed\")\nassert all(ok for _, ok in report)"
  }
 ]
}