{
 "nbformat": 4,
 "nbformat_minor": 5,
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "name": "python",
   "version": "3.11"
  },
  "colab": {
   "provenance": []
  }
 },
 "cells": [
  {
   "cell_type": "markdown",
   "id": "cell-0000",
   "metadata": {},
   "source": [
    "# PCA: Dimensionality Reduction\n",
    "\n",
    "Squeeze 30 features into 2 with principal component analysis — choose the number of components with explained variance, and use PCA as a preprocessing step.\n",
    "\n",
    "*Part of the free [Machine Learning](https://ramadnsyh.dev/courses/machine-learning) course by [Muhammad Ramadiansyah](https://ramadnsyh.dev). This notebook is generated from the interactive lesson — [read it online](https://ramadnsyh.dev/courses/machine-learning/pca).*"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cell-0001",
   "metadata": {},
   "source": [
    "Clustering grouped rows; this lesson compresses **columns**. Real datasets\n",
    "routinely have dozens or hundreds of features, and that width causes real\n",
    "problems: you can't plot it, distances get less meaningful, and models\n",
    "overfit. Principal component analysis (PCA) is the classic fix — it finds a\n",
    "small set of new axes that keep as much of the data's variation as possible.\n",
    "You'll see how it works geometrically, how to choose the number of\n",
    "components, and how to drop it into a scikit-learn pipeline."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cell-0002",
   "metadata": {},
   "source": [
    "## The curse of dimensionality\n",
    "\n",
    "Every extra feature adds a dimension to the space your data lives in, and\n",
    "high-dimensional space behaves badly. The volume grows exponentially, so a\n",
    "fixed number of rows spreads thinner and thinner — with 30 features, 569\n",
    "patients (the breast-cancer dataset we'll use) barely sketch the space.\n",
    "Distances lose contrast: in very high dimensions, the nearest and farthest\n",
    "neighbors of a point end up almost equally far away, which quietly degrades\n",
    "anything built on distance — KNN, K-Means, SVMs with RBF kernels. And more\n",
    "features means more parameters, which means more ways to memorize noise.\n",
    "\n",
    "Dimensionality reduction attacks this from three angles at once:\n",
    "\n",
    "- **Visualization** — compress to 2 or 3 dimensions so you can actually look\n",
    "  at the data\n",
    "- **Compression** — store or transmit the same information in fewer numbers\n",
    "- **Feature extraction** — build a few dense, informative features out of\n",
    "  many redundant ones, and feed those to a model\n",
    "\n",
    "The good news: real features are rarely independent. A tumor's radius,\n",
    "perimeter, and area are three columns telling one story. When features are\n",
    "correlated, the data doesn't fill the 30-dimensional space — it hugs a much\n",
    "lower-dimensional sheet inside it, and PCA finds that sheet."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cell-0003",
   "metadata": {},
   "source": [
    "## Directions of maximal variance\n",
    "\n",
    "PCA's rule is simple: find the direction along which the data varies the\n",
    "most. That's **principal component 1 (PC1)**. Then find the direction at a\n",
    "right angle to it with the most remaining variance — **PC2** — and so on.\n",
    "Projecting the data onto the first few components keeps the spread (the\n",
    "information) and throws away the flat, noisy directions.\n",
    "\n",
    "Try it below. Rotate the candidate axis through the cloud and watch the\n",
    "projected variance change, then project onto PC1 and see how much of the\n",
    "2-D structure survives in 1-D:"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cell-0004",
   "metadata": {},
   "source": [
    "> 🎛️ **Interactive demo** — this section has a hands-on visualization in the web version of this lesson: [open it here](https://ramadnsyh.dev/courses/machine-learning/pca)."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cell-0005",
   "metadata": {},
   "source": [
    "Notice that the best axis isn't either original feature — it's a *diagonal\n",
    "combination* of them. That's the general pattern: each principal component\n",
    "is a weighted mix of all the original features, and the components are\n",
    "computed in one shot via singular value decomposition (SVD), the same matrix\n",
    "factorization trick that powers half of classical ML. You'll also notice\n",
    "what gets lost: the small wiggles perpendicular to PC1. PCA bets that\n",
    "low-variance directions are noise. Usually a good bet — not always."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cell-0006",
   "metadata": {},
   "source": [
    "## Scale first, always\n",
    "\n",
    "PCA chases variance, and variance has units. In the breast-cancer data,\n",
    "`mean area` lives in the hundreds while `mean smoothness` lives around 0.1 —\n",
    "unscaled, the \"biggest\" direction is just whichever feature has the biggest\n",
    "numbers. Watch how badly this skews things:"
   ]
  },
  {
   "cell_type": "code",
   "id": "cell-0007",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": [
    "from sklearn.datasets import load_breast_cancer\n",
    "from sklearn.preprocessing import StandardScaler\n",
    "from sklearn.decomposition import PCA\n",
    "\n",
    "X, y = load_breast_cancer(return_X_y=True)\n",
    "print(f\"Data: {X.shape[0]} patients x {X.shape[1]} features\")\n",
    "\n",
    "pca_raw = PCA().fit(X)\n",
    "pca_scaled = PCA().fit(StandardScaler().fit_transform(X))\n",
    "\n",
    "print(f\"\\\\nPC1 share of variance, unscaled: {pca_raw.explained_variance_ratio_[0]:.1%}\")\n",
    "print(f\"PC1 share of variance, scaled  : {pca_scaled.explained_variance_ratio_[0]:.1%}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cell-0008",
   "metadata": {},
   "source": [
    "Unscaled, PC1 \"explains\" over 98% of the variance — but it's just pointing\n",
    "at the large-unit features, not at structure. After standardizing, PC1\n",
    "carries a believable 44%, and the rest is spread across many components."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cell-0009",
   "metadata": {},
   "source": [
    "> **StandardScaler before PCA**\n",
    "> \n",
    "> Put a `StandardScaler` in front of `PCA` every time your features have\n",
    "> different units — which is essentially every real dataset. In a `Pipeline`,\n",
    "> that's one extra line."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cell-0010",
   "metadata": {},
   "source": [
    "## Choosing the number of components\n",
    "\n",
    "Each fitted PCA exposes `explained_variance_ratio_` — the fraction of total\n",
    "variance each component captures. Its **cumulative sum** is the standard\n",
    "tool for choosing `n_components`: plot it and read off how many components\n",
    "you need to keep, say, 95% of the variance."
   ]
  },
  {
   "cell_type": "code",
   "id": "cell-0011",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from sklearn.datasets import load_breast_cancer\n",
    "from sklearn.preprocessing import StandardScaler\n",
    "from sklearn.decomposition import PCA\n",
    "\n",
    "X, y = load_breast_cancer(return_X_y=True)\n",
    "X_scaled = StandardScaler().fit_transform(X)\n",
    "\n",
    "pca = PCA().fit(X_scaled)\n",
    "cumvar = np.cumsum(pca.explained_variance_ratio_)\n",
    "\n",
    "plt.figure(figsize=(7, 4))\n",
    "plt.plot(range(1, 31), cumvar, \"o-\")\n",
    "plt.axhline(0.95, color=\"gray\", ls=\"--\")\n",
    "plt.xlabel(\"Number of components\")\n",
    "plt.ylabel(\"Cumulative explained variance\")\n",
    "plt.title(\"30 features, but how much real information?\")\n",
    "plt.show()\n",
    "\n",
    "n95 = int(np.argmax(cumvar >= 0.95)) + 1\n",
    "print(f\"Components for 95% of the variance: {n95}\")\n",
    "print(f\"First 2 components already hold {cumvar[1]:.1%}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cell-0012",
   "metadata": {},
   "source": [
    "Ten components carry 95% of the variance of thirty features — a third of\n",
    "the width for almost all of the information. That redundancy is exactly the\n",
    "radius/perimeter/area correlation showing up in the math. A convenient\n",
    "shortcut: `PCA(n_components=0.95)` picks the count for you automatically."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cell-0013",
   "metadata": {},
   "source": [
    "## Thirty dimensions on one screen\n",
    "\n",
    "The most immediate payoff is visualization. Project the standardized data\n",
    "onto PC1 and PC2 and color each patient by diagnosis — a plot that would be\n",
    "impossible with the raw 30 columns:"
   ]
  },
  {
   "cell_type": "code",
   "id": "cell-0014",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "from sklearn.datasets import load_breast_cancer\n",
    "from sklearn.preprocessing import StandardScaler\n",
    "from sklearn.decomposition import PCA\n",
    "\n",
    "data = load_breast_cancer()\n",
    "X_scaled = StandardScaler().fit_transform(data.data)\n",
    "\n",
    "pca = PCA(n_components=2)\n",
    "X2 = pca.fit_transform(X_scaled)\n",
    "\n",
    "plt.figure(figsize=(7, 5))\n",
    "for label, color in [(0, \"tab:red\"), (1, \"tab:blue\")]:\n",
    "    mask = data.target == label\n",
    "    plt.scatter(X2[mask, 0], X2[mask, 1], s=14, c=color,\n",
    "                alpha=0.6, label=data.target_names[label])\n",
    "plt.xlabel(f\"PC1 ({pca.explained_variance_ratio_[0]:.0%} of variance)\")\n",
    "plt.ylabel(f\"PC2 ({pca.explained_variance_ratio_[1]:.0%} of variance)\")\n",
    "plt.legend()\n",
    "plt.title(\"Breast-cancer data in PCA space\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cell-0015",
   "metadata": {},
   "source": [
    "Remember: PCA never saw the diagnosis labels. It compressed 30 measurements\n",
    "into 2 numbers per patient using variance alone — and the malignant and\n",
    "benign groups separate almost cleanly anyway. That's strong evidence the\n",
    "labels are learnable, spotted before training a single classifier."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cell-0016",
   "metadata": {},
   "source": [
    "## PCA as a preprocessing step\n",
    "\n",
    "Because `PCA` is a transformer, it slots straight into a `Pipeline` between\n",
    "the scaler and the model. Fewer, decorrelated inputs can mean faster\n",
    "training and sometimes less overfitting:"
   ]
  },
  {
   "cell_type": "code",
   "id": "cell-0017",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": [
    "from sklearn.datasets import load_breast_cancer\n",
    "from sklearn.pipeline import make_pipeline\n",
    "from sklearn.preprocessing import StandardScaler\n",
    "from sklearn.decomposition import PCA\n",
    "from sklearn.svm import SVC\n",
    "from sklearn.model_selection import cross_val_score\n",
    "\n",
    "X, y = load_breast_cancer(return_X_y=True)\n",
    "\n",
    "models = [\n",
    "    (\"SVC on all 30 features\", make_pipeline(StandardScaler(), SVC())),\n",
    "    (\"SVC on 10 PCA components\", make_pipeline(StandardScaler(), PCA(n_components=10), SVC())),\n",
    "    (\"SVC on 2 PCA components\", make_pipeline(StandardScaler(), PCA(n_components=2), SVC())),\n",
    "]\n",
    "\n",
    "for name, model in models:\n",
    "    scores = cross_val_score(model, X, y, cv=5)\n",
    "    print(f\"{name:26s} accuracy = {scores.mean():.3f}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cell-0018",
   "metadata": {},
   "source": [
    "Ten components match the full model's accuracy with a third of the inputs —\n",
    "and even two components stay remarkably close, which the scatter plot above\n",
    "already predicted. In a real project you'd treat `n_components` as a\n",
    "hyperparameter and let `GridSearchCV` tune it along with the classifier's\n",
    "settings, since the whole pipeline cross-validates as one unit."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cell-0019",
   "metadata": {},
   "source": [
    "## What PCA can't do\n",
    "\n",
    "PCA has two honest limitations. First, it's **linear**: components are\n",
    "straight axes, so if your data curls along a spiral or an S-shaped surface,\n",
    "no rotation captures it — that's what nonlinear methods like t-SNE (next\n",
    "lesson) and kernel PCA are for. Second, components trade interpretability\n",
    "for compactness: \"0.22 times mean radius minus 0.10 times smoothness plus\n",
    "28 more terms\" is much harder to explain to a stakeholder than any original\n",
    "column. And one subtle trap: PCA is unsupervised, so it keeps high-variance\n",
    "directions whether or not they help your prediction task. Almost always they\n",
    "do — but if the signal lives in a low-variance direction, PCA will happily\n",
    "throw it away."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cell-0020",
   "metadata": {},
   "source": [
    "### 🏋️ Exercise — Compress the digits\n",
    "\n",
    "= 0.90) + 1.\",\n",
    "    \"load_digits has 64 features (8x8 pixel images) — expect far fewer components to reach 90%.\",\n",
    "    \"Wrap StandardScaler, PCA, and LogisticRegression(max_iter=2000) in make_pipeline for the comparison.\",\n",
    "  ]}\n",
    "  solution={`\n",
    "import numpy as np\n",
    "from sklearn.datasets import load_digits\n",
    "from sklearn.preprocessing import StandardScaler\n",
    "from sklearn.decomposition import PCA\n",
    "from sklearn.pipeline import make_pipeline\n",
    "from sklearn.linear_model import LogisticRegression\n",
    "from sklearn.model_selection import cross_val_score\n",
    "\n",
    "X, y = load_digits(return_X_y=True)\n",
    "X_scaled = StandardScaler().fit_transform(X)\n",
    "\n",
    "cumvar = PCA().fit(X_scaled).explained_variance_ratio_.cumsum()\n",
    "n90 = int(np.argmax(cumvar >= 0.90)) + 1\n",
    "print(f\"Components for 90% variance: {n90} of {X.shape[1]}\")\n",
    "\n",
    "full = make_pipeline(StandardScaler(),\n",
    "                     LogisticRegression(max_iter=2000))\n",
    "reduced = make_pipeline(StandardScaler(), PCA(n_components=n90),\n",
    "                        LogisticRegression(max_iter=2000))\n",
    "\n",
    "print(f\"all 64 features : {cross_val_score(full, X, y, cv=3).mean():.3f}\")\n",
    "print(f\"{n90} components   : {cross_val_score(reduced, X, y, cv=3).mean():.3f}\")\n",
    "\n",
    "# Roughly 30 components carry 90% of the variance of 64 pixels,\n",
    "# and the classifier loses almost nothing.\n",
    "`}\n",
    ">\n",
    "Repeat the workflow on `load_digits` (64 pixel features). Scale the data,\n",
    "find how many components you need for **90%** of the variance, then compare\n",
    "a logistic-regression pipeline on all 64 features against one that reduces\n",
    "to that number of components first (3-fold cross-validation). How much\n",
    "accuracy does the compression cost?"
   ]
  },
  {
   "cell_type": "code",
   "id": "cell-0021",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": [
    "# Your solution here"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cell-0022",
   "metadata": {},
   "source": [
    "<details><summary>✅ Show solution</summary>\n",
    "\n",
    "```python\n",
    "import numpy as np\n",
    "from sklearn.datasets import load_digits\n",
    "from sklearn.preprocessing import StandardScaler\n",
    "from sklearn.decomposition import PCA\n",
    "from sklearn.pipeline import make_pipeline\n",
    "from sklearn.linear_model import LogisticRegression\n",
    "from sklearn.model_selection import cross_val_score\n",
    "\n",
    "X, y = load_digits(return_X_y=True)\n",
    "X_scaled = StandardScaler().fit_transform(X)\n",
    "\n",
    "cumvar = PCA().fit(X_scaled).explained_variance_ratio_.cumsum()\n",
    "n90 = int(np.argmax(cumvar >= 0.90)) + 1\n",
    "print(f\"Components for 90% variance: {n90} of {X.shape[1]}\")\n",
    "\n",
    "full = make_pipeline(StandardScaler(),\n",
    "                     LogisticRegression(max_iter=2000))\n",
    "reduced = make_pipeline(StandardScaler(), PCA(n_components=n90),\n",
    "                        LogisticRegression(max_iter=2000))\n",
    "\n",
    "print(f\"all 64 features : {cross_val_score(full, X, y, cv=3).mean():.3f}\")\n",
    "print(f\"{n90} components   : {cross_val_score(reduced, X, y, cv=3).mean():.3f}\")\n",
    "\n",
    "# Roughly 30 components carry 90% of the variance of 64 pixels,\n",
    "# and the classifier loses almost nothing.\n",
    "```\n",
    "\n",
    "</details>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cell-0023",
   "metadata": {},
   "source": [
    "Next up: when straight axes aren't enough — t-SNE for visualizing nonlinear\n",
    "structure, and topic modeling for finding themes in text."
   ]
  }
 ]
}