{
 "nbformat": 4,
 "nbformat_minor": 5,
 "metadata": {
  "kernelspec": {
   "name": "python3",
   "display_name": "Python 3",
   "language": "python"
  },
  "language_info": {
   "name": "python"
  },
  "colab": {
   "name": "dimred.ipynb"
  }
 },
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "# dimensionality reduction \u2014 Python demo\n\nNumerical companion to the entry [dimensionality reduction](https://dictionaryofml.org/terms/dimred.html) of the [Dictionary of Applied Machine Learning](https://dictionaryofml.org/): it recomputes what the entry states and prints one line per check.\n\nChecks the three benefits the entry attributes to using fewer features, on synthetic image-like data: the statistical benefit (less overfitting of linear regression on the learned features than on the raw ones), the computational benefit (the matrix that linear regression inverts shrinks with the number of features), and visualization (two learned features place the data points of two digits in separate regions of a scatterplot). It also checks that a random projection approximately preserves distances. Self-contained (numpy and matplotlib only), fixed seed.\n\nRequires NumPy and Matplotlib only, and uses fixed seeds, so the printed numbers reproduce exactly. Generated from [`pythondemos/dimred.py`](https://dictionaryofml.org/terms/dimred.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(), \"dimred.py\")\nos.makedirs(\"pythondemos\", exist_ok=True)"
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "\"\"\"\ndimred.py \u2014 numerical companion to the glossary entry 'dimensionality\nreduction'.\n\nPurpose\n-------\nChecks the three benefits the entry attributes to using fewer features,\non synthetic image-like data: the statistical benefit (less overfitting\nof linear regression on the learned features than on the raw ones), the\ncomputational benefit (the matrix that linear regression inverts shrinks\nwith the number of features), and visualization (two learned features\nplace the data points of two digits in separate regions of a\nscatterplot).  It also checks that a random projection approximately\npreserves distances.  Self-contained (numpy and matplotlib only), fixed\nseed.\n\nSetup\n-----\nEach data point is a vector of d = 50 grayscale values: one of two digit\ntemplates, which differ along one direction, plus a random variation\nalong two further directions plus pixel noise.  The label is a linear\nfunction of the two variation coordinates plus noise.  Training set\nm = 60 data points, and 400 further data points to measure the error on\nnew data points.  The\nlearned transformation is PCA with d' = 2 features; the random projection\nuses d' = 20 features and a matrix of independent Gaussian entries.\n\nBlocks\n------\n[B-pca]   PCA with d' = 2 reconstructs the raw features with an average\n          squared reconstruction error below 20 percent of their total\n          variance, since the data vary mainly along the digit direction\n          and the larger of the two variation directions.\n[B-stat]  Linear regression on the d = 50 raw features reaches an error\n          near zero on the training set but an error more than ten times\n          larger on new data points; on the two PCA features the error on\n          new data points is less than twice the error on the training\n          set and smaller than with the raw features.\n[B-comp]  The matrix inverted by linear regression has 50 x 50 entries\n          with the raw features and 2 x 2 with the PCA features.\n[B-vis]   In the scatterplot of the two PCA features, the data points of\n          the two digits occupy regions whose centers are far apart\n          relative to the spread within each digit.\n[B-rp]    A random projection to d' = 20 features changes every pairwise\n          distance between the 60 training data points by less than a\n          factor of two, and the average distance by less than 15 percent.\n\nOutputs\n-------\ndimred_scatter.csv : z1, z2, digit for the 60 training data points.\ndimred_errors.csv  : features, err_train, err_new for the raw and the\n                     PCA features.\ndimred.png         : matplotlib preview (checking only).\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)\n\nd, d_pca, d_rp = 50, 2, 20\nm_tr, m_va = 60, 400\ntemplates = rng.uniform(0.0, 1.0, (2, d))          # two digit templates\nQ = np.linalg.qr(rng.standard_normal((d, 3)))[0]     # three directions at right angles\ndirections = Q[:, :2]                                # two variation directions\ntemplates[1] = templates[0] + 3.0 * Q[:, 2]          # digits differ along the third\n\n\ndef make(m):\n    digit = rng.integers(0, 2, m)\n    coords = rng.standard_normal((m, 2)) * np.array([1.5, 0.4])\n    X = templates[digit] + coords @ directions.T + 0.1 * rng.standard_normal((m, d))\n    y = 1.0 * coords[:, 0] - 0.5 * coords[:, 1] + 0.1 * rng.standard_normal(m)\n    return X, y, digit\n\n\nX_tr, y_tr, dig_tr = make(m_tr)\nX_va, y_va, dig_va = make(m_va)"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "**[B-pca]** PCA with d' = 2 reconstructs the raw features with an average squared reconstruction error below 20 percent of their total variance, since the data vary mainly along the digit direction and the larger of the two variation directions."
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "mu = X_tr.mean(axis=0)\nC = (X_tr - mu).T @ (X_tr - mu) / m_tr\nevals, evecs = np.linalg.eigh(C)\nW = evecs[:, ::-1][:, :d_pca]                        # d x d' transformation\nZ_tr = (X_tr - mu) @ W; Z_va = (X_va - mu) @ W\nrecon = mu + Z_tr @ W.T\nrec_err = float(np.mean(np.sum((X_tr - recon) ** 2, axis=1)))\ntotal_var = float(np.trace(C))\ncheck(f\"[B-pca]   average squared reconstruction error {rec_err:.3f} = \"\n      f\"{100 * rec_err / total_var:.1f} percent of the total variance \"\n      f\"{total_var:.3f}\", rec_err / total_var < 0.2)"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "**[B-stat]** Linear regression on the d = 50 raw features reaches an error near zero on the training set but an error more than ten times larger on new data points; on the two PCA features the error on new data points is less than twice the error on the training set and smaller than with the raw features."
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "def linreg(Xa, ya, Xb, yb):\n    A = np.c_[Xa, np.ones(len(Xa))]; B = np.c_[Xb, np.ones(len(Xb))]\n    w = np.linalg.lstsq(A, ya, rcond=None)[0]\n    return float(np.mean((ya - A @ w) ** 2)), float(np.mean((yb - B @ w) ** 2))\n\n\ntr_raw, va_raw = linreg(X_tr, y_tr, X_va, y_va)\ntr_pca, va_pca = linreg(Z_tr, y_tr, Z_va, y_va)\ncheck(f\"[B-stat]  raw features: error {tr_raw:.4f} on the training set, \"\n      f\"{va_raw:.3f} on new data points; PCA features: {tr_pca:.4f} and {va_pca:.4f}\",\n      va_raw > 10 * tr_raw and va_pca < 2 * tr_pca and va_pca < va_raw)"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "**[B-comp]** The matrix inverted by linear regression has 50 x 50 entries with the raw features and 2 x 2 with the PCA features."
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "size_raw = (X_tr.T @ X_tr).shape; size_pca = (Z_tr.T @ Z_tr).shape\ncheck(f\"[B-comp]  matrix to invert: {size_raw[0]} x {size_raw[1]} entries \"\n      f\"with raw features, {size_pca[0]} x {size_pca[1]} with PCA features\",\n      size_raw == (d, d) and size_pca == (d_pca, d_pca))"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "**[B-vis]** In the scatterplot of the two PCA features, the data points of the two digits occupy regions whose centers are far apart relative to the spread within each digit."
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "c0, c1 = Z_tr[dig_tr == 0].mean(axis=0), Z_tr[dig_tr == 1].mean(axis=0)\nu = (c1 - c0) / np.linalg.norm(c1 - c0)             # direction between the centers\nspread = max(((Z_tr[dig_tr == 0] - c0) @ u).std(), ((Z_tr[dig_tr == 1] - c1) @ u).std())\nsep = float(np.linalg.norm(c0 - c1))\ncheck(f\"[B-vis]   centers of the two digits {sep:.2f} apart in the \"\n      f\"scatterplot, spread of a digit along that direction at most {spread:.2f}\",\n      sep > 3 * spread)"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "**[B-rp]** A random projection to d' = 20 features changes every pairwise distance between the 60 training data points by less than a factor of two, and the average distance by less than 15 percent."
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "R = rng.standard_normal((d, d_rp)) / np.sqrt(d_rp)\nP_tr = X_tr @ R\n\n\ndef pairwise(A):\n    diff = A[:, None, :] - A[None, :, :]\n    return np.sqrt((diff ** 2).sum(axis=2))[np.triu_indices(len(A), 1)]\n\n\nratio = pairwise(P_tr) / pairwise(X_tr)\ncheck(f\"[B-rp]    random projection to {d_rp} features: distance ratios \"\n      f\"in [{ratio.min():.2f}, {ratio.max():.2f}], mean {ratio.mean():.3f}\",\n      ratio.min() > 0.5 and ratio.max() < 2.0 and abs(ratio.mean() - 1) < 0.15)\n\n# ---------------------------------------------------------------- CSV\nwith open(OUT_DIR / \"dimred_scatter.csv\", \"w\") as fh:\n    fh.write(\"z1,z2,digit\\n\")\n    for z, g in zip(Z_tr, dig_tr):\n        fh.write(f\"{z[0]:.4f},{z[1]:.4f},{g}\\n\")\nwith open(OUT_DIR / \"dimred_errors.csv\", \"w\") as fh:\n    fh.write(\"features,nfeatures,err_train,err_new\\n\")\n    fh.write(f\"raw,{d},{tr_raw:.4f},{va_raw:.4f}\\n\")\n    fh.write(f\"pca,{d_pca},{tr_pca:.4f},{va_pca:.4f}\\n\")\n\n# -------------------------------------------------------------- preview\nfig, (ax1, ax2) = plt.subplots(1, 2, figsize=(9.2, 3.9))\nfor g, mk, lab in ((0, \"o\", \"digit A\"), (1, \"^\", \"digit B\")):\n    sel = dig_tr == g\n    ax1.plot(Z_tr[sel, 0], Z_tr[sel, 1], mk, color=\"k\", mfc=\"k\" if g == 0 else \"none\",\n             ms=5, label=lab)\nax1.set_xlabel(\"learned feature $z_1$\"); ax1.set_ylabel(\"learned feature $z_2$\")\nax1.set_title(\"two learned features separate the digits\")\nax1.legend(frameon=False, fontsize=8)\nxs = np.arange(2)\nax2.plot(xs, [tr_raw, tr_pca], \"ko-\", ms=7, lw=1.0, label=\"error on the training set\")\nax2.plot(xs, [va_raw, va_pca], \"k^--\", mfc=\"none\", ms=8, lw=1.0,\n         label=\"error on new data points\")\nax2.set_xticks(xs); ax2.set_xticklabels([f\"{d} raw features\", f\"{d_pca} PCA features\"])\nax2.set_xlim(-0.5, 1.5)\nax2.set_ylabel(\"average squared error\"); ax2.set_xlabel(\"features used by linear regression\")\nax2.set_title(\"fewer features, less overfitting\")\nax2.legend(frameon=False, fontsize=8)\nfig.tight_layout()\nfig.savefig(OUT_DIR / \"dimred.png\", dpi=110)\n\nn_ok = sum(ok for _, ok in report)\nprint(f\"\\n{n_ok}/{len(report)} checks pass\")\nprint(f\"wrote {OUT_DIR / 'dimred_scatter.csv'}, {OUT_DIR / 'dimred_errors.csv'}, \"\n      f\"{OUT_DIR / 'dimred.png'}\")\nif n_ok != len(report):\n    raise SystemExit(1)"
  }
 ]
}