{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "2b940af6",
   "metadata": {},
   "source": [
    "# Lesson 30: Projection — The Central Idea\n",
    "\n",
    "This begins Part 3 of the course: deep learning for computer vision. Before writing a single neural network, it's worth pausing to reflect upon an idea that rarely gets much attention despite its fundamental importance. **Projection** is a simple idea: map points from one space into another, usually a simpler one, by a linear combination of their coordinates. Projection quietly ran through much of Parts 1 and 2 — many of the techniques so far have secretly been projections. Projection also underlies everything from here on — neural networks turn out to be nothing more than *learned*, *composed* projections."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "825f6d6c",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ad1bf246",
   "metadata": {},
   "source": [
    "## Two projections you've already built\n",
    "\n",
    "**PCA (Lesson 6).** Given a cloud of points, the eigenvectors of their covariance matrix gave directions to project onto. Projecting onto the top eigenvector, $y = v^\\top x$, collapses each 2D point to a single number — the coordinate along the direction of greatest spread."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5f48a370",
   "metadata": {},
   "outputs": [],
   "source": [
    "rng = np.random.default_rng(1)\n",
    "cloud = rng.multivariate_normal([0, 0], [[3, 1.5], [1.5, 1]], 100)\n",
    "\n",
    "cov = np.cov(cloud.T)\n",
    "eigvals, eigvecs = np.linalg.eigh(cov)\n",
    "principal_direction = eigvecs[:, -1]  # eigenvector of the largest eigenvalue\n",
    "\n",
    "projected = cloud @ principal_direction\n",
    "\n",
    "fig, axes = plt.subplots(1, 2, figsize=(8, 3.5))\n",
    "axes[0].scatter(cloud[:, 0], cloud[:, 1], s=15, alpha=0.6)\n",
    "axes[0].plot([0, 3 * principal_direction[0]], [0, 3 * principal_direction[1]], color='red', linewidth=2)\n",
    "axes[0].set_aspect('equal')\n",
    "axes[0].set_title('2D cloud + principal direction')\n",
    "axes[1].scatter(projected, np.zeros_like(projected), s=15, alpha=0.6)\n",
    "axes[1].set_yticks([])\n",
    "axes[1].set_title('Projected onto that direction (1D)')\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "86b6f6cb",
   "metadata": {},
   "source": [
    "**Camera projection (Lessons 24, 25).** $P = K[R|t]$ maps a 3D world point to a 2D pixel. Every row of that matrix, before the final homogeneous divide, is itself a linear projection: a dot product between the point's coordinates and a fixed direction, plus an offset."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b2e3a6ab",
   "metadata": {},
   "source": [
    "Both examples share the same mechanics: pick a direction $w$ (and maybe an offset $b$), then compute $y = w^\\top x + b$. What differs is *why* the direction was chosen: PCA picks $w$ to maximize the spread of the projected data (an unsupervised, purely geometric objective). A camera's rows are fixed by its physical geometry. Neither one is chosen to solve a classification or recognition problem — but nothing stops us from choosing $w$ for exactly that purpose."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e28c8e9e",
   "metadata": {},
   "source": [
    "## Matrix multiplication: many projections at once\n",
    "\n",
    "Everything above computed one number $y = w^\\top x$ from one direction $w$. Stack $m$ different directions as the rows of a matrix $W$ (shape $m \\times n$), and $y = Wx$ computes all $m$ projections at once — one dot product per row. This is exactly what matrix multiplication *is*: a bank of simultaneous projections. It also means that several seemingly different linear operations from earlier in this course are secretly the same operation, just with a different, structured choice of $W$:\n",
    "\n",
    "- **Convolution** (Lesson 10): each output value is a dot product between a small neighborhood and a kernel — a projection. Sliding the kernel across the whole signal is a matrix multiply by a sparse matrix whose rows are shifted copies of the kernel (a Toeplitz matrix).\n",
    "- **The Fourier transform** (Lesson 15): each output value is a dot product between the *whole* signal and one sinusoid at a fixed frequency — a projection onto that frequency. The full transform is a matrix multiply by a fixed matrix of sines and cosines.\n",
    "- **A homography or camera projection matrix** (Lessons 24, 25): each output coordinate is a dot product between a homogeneous point and one row of $H$ or $P$.\n",
    "\n",
    "We verify the first two directly: build the matrix by hand, and check it reproduces `np.convolve`/`np.fft.fft` exactly."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3132ab86",
   "metadata": {},
   "outputs": [],
   "source": [
    "signal = rng.normal(size=10)\n",
    "kernel = np.array([0.25, 0.5, 0.25])  # a small smoothing kernel, as in Lesson 10\n",
    "\n",
    "# a matrix whose rows are shifted copies of the kernel -- a Toeplitz matrix\n",
    "n, k = len(signal), len(kernel)\n",
    "W_conv = np.zeros((n - k + 1, n))\n",
    "for i in range(n - k + 1):\n",
    "    W_conv[i, i:i + k] = kernel\n",
    "\n",
    "via_matmul = W_conv @ signal\n",
    "via_convolve = np.convolve(signal, kernel, mode='valid')\n",
    "print(f'max diff, convolution as a matrix multiply: {np.abs(via_matmul - via_convolve).max():.2e}')\n",
    "\n",
    "N = 8\n",
    "idx = np.arange(N)\n",
    "W_dft = np.exp(-2j * np.pi * np.outer(idx, idx) / N)  # row k: the frequency-k sinusoid, sampled at every point\n",
    "\n",
    "x = rng.normal(size=N)\n",
    "via_matmul = W_dft @ x\n",
    "via_fft = np.fft.fft(x)\n",
    "print(f'max diff, DFT as a matrix multiply:          {np.abs(via_matmul - via_fft).max():.2e}')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8eea35f6",
   "metadata": {},
   "source": [
    "A fully-connected neural network layer is nothing more than this same $y = Wx + b$, with $W$'s rows *learned* to solve a task instead of fixed by hand-derived structure (a kernel shape, a bank of sinusoids, a camera's geometry)."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4f0b8ec4",
   "metadata": {},
   "source": [
    "## A projection for classification\n",
    "\n",
    "Suppose instead of \"maximize spread,\" the goal is \"separate two classes.\" The same formula $y = w^\\top x + b$ still applies — project each point onto a direction $w$, and classify by the *sign* of the resulting score. This is, in its entirety, a single artificial neuron with no activation function: the linear core that every neural network layer is built from."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9e568a7c",
   "metadata": {},
   "outputs": [],
   "source": [
    "class_a = rng.normal(loc=[-2, -1], scale=0.8, size=(60, 2))\n",
    "class_b = rng.normal(loc=[2, 1.5], scale=0.8, size=(60, 2))\n",
    "\n",
    "# a principled hand-picked direction: point from one class's mean toward the other's\n",
    "w = class_b.mean(axis=0) - class_a.mean(axis=0)\n",
    "w = w / np.linalg.norm(w)\n",
    "midpoint = (class_a.mean(axis=0) + class_b.mean(axis=0)) / 2\n",
    "b = -w @ midpoint  # threshold: the boundary passes through the midpoint between the classes\n",
    "\n",
    "score_a = class_a @ w + b\n",
    "score_b = class_b @ w + b\n",
    "accuracy = (np.sum(score_a < 0) + np.sum(score_b > 0)) / (len(score_a) + len(score_b))\n",
    "print(f'w = {np.round(w, 3)}, b = {b:.3f}')\n",
    "print(f'classification accuracy: {accuracy:.1%}')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "10fa3af0",
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_projection_classifier(ax, class_a, class_b, w, b, title):\n",
    "    ax.scatter(*class_a.T, s=15, label='class A')\n",
    "    ax.scatter(*class_b.T, s=15, label='class B')\n",
    "\n",
    "    # decision boundary: the line {x : w.x + b = 0}, drawn through its closest point to the origin\n",
    "    foot = -b * w\n",
    "    perp = np.array([-w[1], w[0]])\n",
    "    p1, p2 = foot + 5 * perp, foot - 5 * perp\n",
    "    ax.plot([p1[0], p2[0]], [p1[1], p2[1]], color='black', linewidth=1.5, label='decision boundary')\n",
    "    ax.arrow(*foot, *w, head_width=0.15, color='red', length_includes_head=True, label='w')\n",
    "\n",
    "    ax.set_aspect('equal')\n",
    "    ax.set_title(title, fontsize=9)\n",
    "\n",
    "fig, axes = plt.subplots(1, 2, figsize=(9, 4))\n",
    "plot_projection_classifier(axes[0], class_a, class_b, w, b, f'Decision boundary\\naccuracy={accuracy:.0%}')\n",
    "axes[0].legend(fontsize=7)\n",
    "\n",
    "axes[1].scatter(score_a, np.zeros_like(score_a), s=15, label='class A')\n",
    "axes[1].scatter(score_b, np.zeros_like(score_b), s=15, label='class B')\n",
    "axes[1].axvline(0, color='black', linewidth=1.5, label='threshold')\n",
    "axes[1].set_yticks([])\n",
    "axes[1].set_title('The same points, projected to 1D')\n",
    "axes[1].legend(fontsize=7)\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7a029d47",
   "metadata": {},
   "source": [
    "Two views of the exact same operation: on the left, $w$ is a direction in the original 2D space and the decision boundary is the line perpendicular to it; on the right, every point has been projected down onto that direction, and classification is just thresholding a single number at zero. The 2D picture is more intuitive, but the 1D picture is what actually generalizes &mdash; in 100 dimensions there's no picture to draw of the \"boundary,\" but \"project to a single number, then threshold\" still works exactly the same way."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3ec5a67c",
   "metadata": {},
   "source": [
    "## What a single projection cannot do\n",
    "\n",
    "A single linear projection can only produce a *straight-line* decision boundary (a hyperplane, in higher dimensions). Some datasets have no straight-line separator at all, no matter how $w$ and $b$ are chosen &mdash; the classic example is one class surrounding the other."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1937d148",
   "metadata": {},
   "outputs": [],
   "source": [
    "theta_in = rng.uniform(0, 2 * np.pi, 60)\n",
    "inner = np.stack([0.5 * np.cos(theta_in), 0.5 * np.sin(theta_in)], axis=1) + rng.normal(0, 0.1, (60, 2))\n",
    "theta_out = rng.uniform(0, 2 * np.pi, 60)\n",
    "outer = np.stack([2.0 * np.cos(theta_out), 2.0 * np.sin(theta_out)], axis=1) + rng.normal(0, 0.15, (60, 2))\n",
    "\n",
    "# search many directions and thresholds for the best possible LINEAR separator\n",
    "best_acc, best_w, best_b = 0, None, None\n",
    "for _ in range(2000):\n",
    "    w_try = rng.normal(size=2)\n",
    "    w_try /= np.linalg.norm(w_try)\n",
    "    for b_try in np.linspace(-3, 3, 61):\n",
    "        s_in, s_out = inner @ w_try + b_try, outer @ w_try + b_try\n",
    "        acc = max((s_in < 0).sum() + (s_out > 0).sum(), (s_in > 0).sum() + (s_out < 0).sum()) / 120\n",
    "        if acc > best_acc:\n",
    "            best_acc, best_w, best_b = acc, w_try, b_try\n",
    "\n",
    "print(f'best achievable accuracy with ANY single linear projection: {best_acc:.1%}')\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(4.5, 4.5))\n",
    "plot_projection_classifier(ax, inner, outer, best_w, best_b, f'Best possible linear boundary\\naccuracy={best_acc:.0%}')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "11deba7d",
   "metadata": {},
   "source": "No rotation or shift of a single straight line can separate a ring from the disk it surrounds &mdash; the best any linear projection can manage is a mediocre compromise. This is the wall every purely linear method runs into, and it's exactly what motivates the next two lessons: first, *learning* $w$ and $b$ automatically instead of hand-picking them (Lesson 31), and then, more importantly, *composing several projections with nonlinearities in between* (Lesson 32), which, given enough capacity, can separate regions of essentially any shape."
  },
  {
   "cell_type": "markdown",
   "id": "8aee3469",
   "metadata": {},
   "source": [
    "### Exercises\n",
    "\n",
    "1. In the PCA recap, project `cloud` onto the *second* (smaller) eigenvector instead of the principal one. How does the spread of the resulting 1D projection compare to projecting onto the principal direction, and why does that match what the eigenvalue itself tells you (Lesson 6)?\n",
    "2. Modify `class_a`/`class_b`'s means and spread so the classes overlap more. At what point does the mean-difference direction `w` stop achieving high accuracy, and can you find a *better* `w` than the mean-difference one by hand for that harder case?\n",
    "3. For the ring-and-disk dataset, instead of a linear projection, try classifying by a simple nonlinear function of the points: the distance from the origin, $r = \\sqrt{x_1^2+x_2^2}$, thresholded at some value. What accuracy does this one nonlinear feature achieve, and what does that suggest about *what kind* of projection would actually solve this problem?"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "name": "python",
   "version": "3.x"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}