week_3/3.2_dimensionality_reduction/dimensionality_reduction.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Lecture 3.2: Dimensionality Reduction\n",
    "\n",
    "This lecture, we are going to reduce the dimensionality of a dataset using PCA.\n",
    "\n",
    "**Learning goals:**\n",
    "\n",
    "- Use dimensionality reduction to compress images\n",
    "- Implement PCA on an emoji dataset\n",
    "- Analyse the cumulative explained variance of a PCA model\n",
    "- Use dimensionality reduction for data visualization\n",
    "\n",
    "## 1. Introduction\n",
    "\n",
    "In the lecture slides, we have seen a toy example of PCA which projected a 2D dataset down to 1D. But we have also learned that dimensionality reduction shines when applied to datasets with many dimensions. We saw how pictures were good examples, since they have as many dimensions as they have pixels.\n",
    "\n",
    "Let's try use PCA to reduce the dimensionality of an emoji picture dataset.\n",
    "\n",
    "## 2. Data Munging\n",
    "\n",
    "First let's load the images using [pillow](https://pillow.readthedocs.io/en/stable/), like in lecture 2.5:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-02-08T22:13:34.850375Z",
     "start_time": "2022-02-08T22:13:33.938414Z"
    }
   },
   "outputs": [],
   "source": [
    "from PIL import Image\n",
    "import glob\n",
    "\n",
    "def open_image(path):\n",
    "    with open(path, 'rb') as f:\n",
    "        img = Image.open(f)\n",
    "        img.load()\n",
    "    return img\n",
    "\n",
    "paths = glob.glob('emojis/*.png')\n",
    "\n",
    "images = [open_image(p) for p in paths]\n",
    "images[458]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The images are in grayscale to avoid having even more dimensions with three RGB channels. Speaking of which, how many pixels are there per image?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-02-08T22:13:34.854750Z",
     "start_time": "2022-02-08T22:13:34.852175Z"
    }
   },
   "outputs": [],
   "source": [
    "print(len(images))\n",
    "print(images[0].size)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We have 1067 emojis 🎊. $64\\times64$ pixels means $4096$ dimensional vectors. That's a lot! Our goal is to reduce those vectors' dimensions using PCA.\n",
    "\n",
    "The [sklearn PCA class](https://scikit-learn.org/stable/modules/generated/sklearn.decomposition.PCA.html) expects data in the form of NumPy arrays, where each row is one __example__, and each column is one _feature_ (like our tabular data in lectures 1.2 & 1.3). \n",
    "\n",
    "Let's convert these pillow `Image`s to NumPy `ndarrays`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-02-08T22:13:37.520703Z",
     "start_time": "2022-02-08T22:13:37.502685Z"
    }
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "arrays = [np.asarray(im) for im in images]\n",
    "arrays[0].shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Remember how images are just arrays of numbers? We have loaded our grayscale emojis into $64 \\times 64$ _matrices_. But we want $4096 \\times 1$ _vectors_."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For this, we can use the NumPy [`.reshape()`](https://docs.scipy.org/doc/numpy/reference/generated/numpy.reshape.html) function:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-02-08T22:13:39.524656Z",
     "start_time": "2022-02-08T22:13:39.518641Z"
    }
   },
   "outputs": [],
   "source": [
    "vectors = [arr.reshape((4096,)) for arr in arrays]\n",
    "vectors[0].shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We have our rows, but they are \"stuck\" in a list. We want them to be in a large matrix! We can stack this list of vectors into one big `ndarray` with [`np.stack`](https://docs.scipy.org/doc/numpy/reference/generated/numpy.stack.html), leaving us with a $1067 \\times 4096$ matrix where each row is a vector corresponding to one emoji:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-02-08T22:13:40.130137Z",
     "start_time": "2022-02-08T22:13:40.123592Z"
    }
   },
   "outputs": [],
   "source": [
    "data = np.stack(vectors)\n",
    "data.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-02-08T22:13:40.555395Z",
     "start_time": "2022-02-08T22:13:40.551949Z"
    }
   },
   "outputs": [],
   "source": [
    "data[0]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3. Training\n",
    "\n",
    "We are now ready to _train_ our PCA model with our `data`. sklearn makes this super easy:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-02-08T22:13:43.483173Z",
     "start_time": "2022-02-08T22:13:43.049100Z"
    }
   },
   "outputs": [],
   "source": [
    "from sklearn.decomposition import PCA\n",
    "pca = PCA()\n",
    "pca.fit(data)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Just like k-Means (lecture 3.1), PCA is a _learning_ algorithm. Therefore, we must _train_ it, then _use_ it. But before we use our model to reduce dimensions, we can investigate the _fitting_ process by having a look at the [eigenvectors](https://web.mit.edu/be.400/www/SVD/Singular_Value_Decomposition.htm). Since we are dealing with image data, these can be visualised as follows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-02-08T22:13:45.297882Z",
     "start_time": "2022-02-08T22:13:44.682193Z"
    }
   },
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "\n",
    "fig, axes = plt.subplots(3, 8, figsize=(9, 4),\n",
    "                         subplot_kw={'xticks':[], 'yticks':[]},\n",
    "                         gridspec_kw=dict(hspace=0.1, wspace=0.1))\n",
    "for i, ax in enumerate(axes.flat):\n",
    "    ax.imshow(pca.components_[i].reshape(64, 64), cmap='bone')\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Spooky images! 👻 These are representations of the optimal vector basis that our PCA model has inferred. The most \"important\" eigenvectors, i.e those that explain most of the data's variance, are found in the top left. Notice how the most \"important\" eigenvectors first focus on general lightness/darkness, then refine common orientations and patterns found in emojis.\n",
    "\n",
    "But enough analysis, let's reduce some dimensions! And how many dimensions do we want to reduce to exactly? There is a common method to help us decide, which plots the cumulative explained variance of the eigenvectors vs the number of principal components used. This gives us a rough measure of how much information is lost in the compression. More details [here](https://stats.stackexchange.com/questions/22569/pca-and-proportion-of-variance-explained). \n",
    "\n",
    "Let's construct a `PCA` model with no arguments, which means that it will conserve _all_ principal components. This allows us to plot the explained proportion of variance graph:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-02-08T20:57:04.009354Z",
     "start_time": "2022-02-08T20:57:03.829712Z"
    }
   },
   "outputs": [],
   "source": [
    "plt.plot(np.cumsum(pca.explained_variance_ratio_))\n",
    "plt.xlabel('number of components')\n",
    "plt.ylabel('cumulative explained variance');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Notice how ~ 200 components still \"explains\" more than 90% of the variance of the data, despite being 20 times less dimensions than the original images! Let's project our emojis down to 200 dimensions:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 4. Dimensionality Reduction"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-02-08T20:57:16.988103Z",
     "start_time": "2022-02-08T20:57:16.589399Z"
    }
   },
   "outputs": [],
   "source": [
    "from sklearn.decomposition import PCA\n",
    "pca = PCA(n_components=200)\n",
    "pca.fit(data)\n",
    "components = pca.transform(data)\n",
    "len(components[0])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As promised, our data is reduced to 200 dimensions! However this exercise is futile if the compression destroyed too information... Luckily we are dealing with image data, so we can visually investigate:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-02-08T21:08:44.838656Z",
     "start_time": "2022-02-08T21:08:44.343217Z"
    }
   },
   "outputs": [],
   "source": [
    "# inverse transform to project back to 4096 dimensions\n",
    "projected = pca.inverse_transform(components)\n",
    "\n",
    "# plot the images\n",
    "fig, ax = plt.subplots(2, 10, figsize=(10, 2.5),\n",
    "                       subplot_kw={'xticks':[], 'yticks':[]},\n",
    "                       gridspec_kw=dict(hspace=0.1, wspace=0.1))\n",
    "for i in range(10):\n",
    "    ax[0, i].imshow(data[i].reshape(64, 64), cmap='binary_r')\n",
    "    ax[1, i].imshow(projected[i].reshape(64, 64), cmap='binary_r')\n",
    "    \n",
    "ax[0, 0].set_ylabel('full-dim\\ninput')\n",
    "ax[1, 0].set_ylabel('200-dim\\nreconstruction');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Those images got quite blurry! Notice how we can still identify the emojis despite having removed ~ 3800 dimensions. This is like deleting 95% of the columns of a large dataframe, and still having access to most of the data! 🤔 Think of how difficult it would have been to write a _rule-based algorithm_ to achieve the same task.\n",
    "\n",
    "ℹ️ Also, note that this compression was not content-aware, and had nothing to do with repetitions of pixels. This was pure mathematics using the power of singular value decomposition and eigenvector bases. Impressive! 🧙‍♂️\n",
    "\n",
    "💪 Play around with different values of `n_components` to get a feel of how much information is lost. What's the minimum value which still allows you to recognise the Vulcan salutation?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 5. Visualization"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Earlier, we learned that dimensionality reduction can be used in data visualization to help graph high dimensional data. Our emoji dataset is _very_ high dimensional, and as we've seen in the cumulative explained variance graph, 2 dimensions won't have much explanatory power at all. Nevertheless, we can hope to see trends in a 2D graphs that might correspond to those first light vs darkness eigenvectors.\n",
    "\n",
    "For this, we create a `.plot_components()` function to visualise emoji thumbnails on the 2D projection:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-02-07T21:22:46.062593Z",
     "start_time": "2022-02-07T21:22:46.055675Z"
    }
   },
   "outputs": [],
   "source": [
    "from matplotlib import offsetbox\n",
    "\n",
    "def plot_components(data, model, images=None, ax=None,\n",
    "                    thumb_frac=0.05, cmap='gray'):\n",
    "    ax = ax or plt.gca()\n",
    "    \n",
    "    proj = model.fit_transform(data)\n",
    "    ax.plot(proj[:, 0], proj[:, 1], '.k')\n",
    "    \n",
    "    if images is not None:\n",
    "        min_dist_2 = (thumb_frac * max(proj.max(0) - proj.min(0))) ** 2\n",
    "        shown_images = np.array([2 * proj.max(0)])\n",
    "        for i in range(data.shape[0]):\n",
    "            dist = np.sum((proj[i] - shown_images) ** 2, 1)\n",
    "            if np.min(dist) < min_dist_2:\n",
    "                # don't show points that are too close\n",
    "                continue\n",
    "            shown_images = np.vstack([shown_images, proj[i]])\n",
    "            imagebox = offsetbox.AnnotationBbox(\n",
    "                offsetbox.OffsetImage(images[i], cmap=cmap),\n",
    "                                      proj[i])\n",
    "            ax.add_artist(imagebox)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-02-07T21:23:15.966016Z",
     "start_time": "2022-02-07T21:23:14.820259Z"
    }
   },
   "outputs": [],
   "source": [
    "small_images = [im[::4, ::4] for im in arrays]\n",
    "fig, ax = plt.subplots(figsize=(10, 10))\n",
    "plot_components(data,\n",
    "                model=PCA(n_components=2),\n",
    "                images=small_images, thumb_frac=0.02)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The result is pretty good considering we just projected 4096 dimensions down to 2! Notice how the only \"variations\" that were conserved by PCA were light/dark, thin/large. This must be because these features explain most of the variation in the data. \n",
    "\n",
    "🧠🧠 Can you think of ways this could be combined with clustering (lecture 3.1) to visualize high dimensional unlabeled datasets?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Summary\n",
    "\n",
    "Today was our introduction to **dimensionality reduction**. We learned how methods like **PCA** can project a dataset down to **fewer dimensions**. We explained how PCA eliminates the **principal components** which least explain the data. Finally, we implemented PCA to reduce the dimensions of an emoji image dataset, which we used both for **compression**, and **data visualization**.\n",
    "\n",
    "\n",
    "# Resources\n",
    "\n",
    "\n",
    "### Core Resources\n",
    "\n",
    "- [Python data science handbook](https://jakevdp.github.io/PythonDataScienceHandbook/05.09-principal-component-analysis.html)  \n",
    "Great blog post on the mechanisms and uses of PCA\n",
    "- [Coursera - PCA Algorithm](https://www.coursera.org/lecture/machine-learning/principal-component-analysis-algorithm-ZYIPa)  \n",
    "Excellent video from the legendary Andrew Ng, going more in depth into the mathematics of the PCA algorithm\n",
    "- [Tutorial on PCA](https://arxiv.org/pdf/1404.1100.pdf)  \n",
    "Comprehensive tutorial from Jonathon Shlens from Google Research, intuitive mathematical explanation of PCA\n",
    "\n",
    "### Additional Resources\n",
    "\n",
    "- [PCA 4 dummies](https://georgemdallas.wordpress.com/2013/10/30/principal-component-analysis-4-dummies-eigenvectors-eigenvalues-and-dimension-reduction/)  \n",
    "Short but effective visualisations to explain PCA"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.4"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}