mathematics/code/pca_solution.ipynb

{
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "gaylQbZDqdYp"
      },
      "source": [
        "## Learning objectives\n",
        "1. Write code that implements PCA.\n",
        "2. Write code that implements PCA for high-dimensional datasets"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "vtOueCzkqdYq"
      },
      "source": [
        "Let's first import the packages we need for this week."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "3t7Hz7ZFqdYq"
      },
      "outputs": [],
      "source": [
        "# PACKAGE: DO NOT EDIT THIS CELL\n",
        "import numpy as np\n",
        "import timeit\n",
        "import matplotlib as mpl\n",
        "mpl.use('Agg')\n",
        "import matplotlib.pyplot as plt\n",
        "plt.style.use('fivethirtyeight')\n",
        "from ipywidgets import interact\n",
        "from sklearn.datasets import fetch_mldata"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "DdiCa79CqdYr"
      },
      "outputs": [],
      "source": [
        "# PACKAGE: DO NOT EDIT THIS CELL\n",
        "from sklearn.datasets import fetch_openml\n",
        "images, labels = fetch_openml('mnist_784', version=1, return_X_y=True)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "rI-tBSg_qdYs"
      },
      "outputs": [],
      "source": [
        "%matplotlib inline"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "QgIdwNcfqdYs"
      },
      "source": [
        "Now, let's plot a digit from the dataset:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "hfArXLuCqdYs",
        "outputId": "05f44485-83a1-4bb2-e992-ed8f7c48bdb3"
      },
      "outputs": [
        {
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAPgAAAD1CAYAAAB9TzjVAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAP+klEQVR4nO3df4wUdZrH8c9ziH/IrgPcushxKuumg4dEOV3wTo1ikJUlEkVRl0R0spz6h6xszpBz0WTRDYSc6CZEvRB/IBLl3LiiaMyNZEQnGxMcbFH5Ea91Y1hwIirM8EvjAc/9MTW93cP0t3t6arqbL+9X0pmqerqqH4r5TFVXdVeZuwtAnP6u3g0AGDwEHIgYAQciRsCBiBFwIGKnDNaCu7q6ODwP1FBTU5P1njagLbiZTTezT8zsUzO7byDLApC+qgNuZkMkPS7pF5LGS5pjZuPTagxACty9qoekf5XUUjD+W0m/7Rnv7Oz0nockl+Tt7e354UZ70Bu9nYh9Feasr5wOZBd9jKS/FozvSqYBaBBW7UdVzewmSde4+78l43MlTXb3X0vFB9lyuVwKrQLoLZPJ5If7OsjGLjq90dsJ3Ndg7qK3S8qY2U/M7FRJv5S0fgDLA5Cyqs+Du/sRM5svqUXSEEnPuPu21DoDMGAD+qCLu78h6Y2UegGQMj6qCkSMgAMRI+BAxAg4EDECDkSMgAMRI+BAxAg4EDECDkSMgAMRI+BAxAg4EDECDkSMgAMRI+BAxAg4EDECDkSMgAMRI+BAxAg4EDECDkSMgAMRI+BAxAg4EDECDkSMgAMRI+BAxAg4EDECDkRsQHcXReMZMmRIsN7U1JT6a44cOTI/PH/+/JLPO+2004LLGTduXLB+9913B+vLly8/btoLL7wgSZozZ05w3u+++y5YX7ZsWbD+4IMPBuv1MqCAm9nnkg5IOirpiLv/LI2mAKQjjS34Ve7+dQrLAZAy3oMDERtowF3Sm2b2vpndmUZDANJj7l79zGb/4O5fmNmPJW2Q9Gt3b5Okrq6u/IJzudyAGwVwvEwmkx9uamqy3vUBBbxoQWaLJR109+VSccCHDx8uSWpvb9ekSZNSeb20xdJbrY+it7S06JprrsmPN9JR9Ewmk9+4NNJR9DR/1zo7O/PDfQW86l10MxtmZj/sGZb0c0lbq10egPQN5Cj6KEnrzKxnOS+4+/+k0tUJ7uyzzw7WTz311GD90ksvPW7abbfdlh++/PLLS87bs7dUyo033his91c2m9VXX32VyrJ27doVrK9YsSJYnzVrVtF4NpvVLbfcIkk6cOBAcN4PP/wwWH/nnXeC9UZVdcDd/S+SLkyxFwAp4zQZEDECDkSMgAMRI+BAxAg4EDG+LlqFiRMnButvvfVWsN7fD5tks1mtWrWqX/M0omPHjgXrDzzwQLB+8ODBYP35558vGl+0aJFmz54tSero6AjOu2/fvmD9k08+CdYbFVtwIGIEHIgYAQciRsCBiBFwIGIEHIgYAQcixnnwKuzcuTNY/+abb4L1wbh0cVo2bdoUrBdeYECSzjjjDLW0tOTHr7rqqpLzfv/998Flr1mzpoIOK7do0SKtW7cu1WWeaNiCAxEj4EDECDgQMQIORIyAAxEj4EDECDgQMc6DV2Hv3r3B+sKFC4P1a6+9Nlj/4IMPisabm5t1zz335MfLXT44ZMuWLcH6tGnTgvVDhw4Vjbe3t2vGjBn58fPPP7/kvAsWLKigQ6SJLTgQMQIORIyAAxEj4EDECDgQMQIORIyAAxHjPPggeOWVV4L1ctdN732r2+bmZj3xxBP58QsvLH1T13nz5gWXvXz58mC993nu/tq2bVvJ2p133jmgZaP/ym7BzewZM9tjZlsLpo00sw1mlkt+jhjcNgFUo5Jd9GclTe817T5Jre6ekdSajANoMGUD7u5tknp/NvM6SauT4dWSrk+5LwApMHcv/ySzsZJed/cJyXinuw8vqO9z96Ld9K6urvyCc7lcWv0CKJDJZPLDTU1N1rtek4NskyZNktT9xYSe4UZTy95OP/30YL33Qbb33ntPkydPzo+vXLmy5LzlDrLdeuutwfratWuD9d74P+2/NPvqfRHM3qo9TfalmY2WpOTnniqXA2AQVRvw9ZJuT4Zvl/RqOu0ASFPZXXQzWytpiqQfmdkuSb+TtEzSH81snqSdkm4azCZjs3///n7PU3ispKurq+rXvuOOO4L1F198MVgvd49vNJayAXf3OSVKU1PuBUDK+KgqEDECDkSMgAMRI+BAxAg4EDG+LnoCWrx4ccnaxRdfHJz3yiuvDNavvvrqYP3NN98M1tFY2IIDESPgQMQIOBAxAg5EjIADESPgQMQIOBAxzoOfgEKXNi73ddBsNhusP/nkk8H6xo0bj5v27LPP5oc3b95cct7HH388uOxKLh+G/mELDkSMgAMRI+BAxAg4EDECDkSMgAMRI+BAxDgPHpnPPvssWG9ubg7WV61aFazPnTu3aDybzRZN610vNGzYsOCyn3vuuWC9o6MjWMfx2IIDESPgQMQIOBAxAg5EjIADESPgQMQIOBAxzoOfZNatWxes53K5YP3RRx8tGh8xYoRaW1vz41Onlr7p7NKlS4PLPuecc4L1JUuWBOu7d+8O1k9GZbfgZvaMme0xs60F0xab2W4z25I8ZgxumwCqUcku+rOSpvcx/Q/uPjF5vJFuWwDSUDbg7t4maW8NegGQMqvkOlhmNlbS6+4+IRlfLKlZ0n5JmyXd6+77Cufp6urKL7jc+zoA1clkMvnhpqYm612vNuCjJH0tySX9XtJod/9V4TyFAR8+fLgkqb29XZMmTer/v6IG6K3bhAkTgvW+DrLt2/e3v+2hg2zlrFy5Mljv70G2Rv0/TbOvzs7O/HBfAa/qNJm7f+nuR939mKQnJU2uukMAg6aqgJvZ6ILRWZK2lnougPopex7czNZKmiLpR2a2S9LvJE0xs4nq3kX/XNJdg9gjamjr1vDf6ptvvrlofMOGDUXTZs6cWXLect81v+uu8K9R4fvNvkybNi1YPxmVDbi7z+lj8tOD0AuAlPFRVSBiBByIGAEHIkbAgYgRcCBifF0U/VL4yam+pq1Zs6bkvE899VRw2aecEv51vOKKK4L1KVOmlJz29ttvB+eNFVtwIGIEHIgYAQciRsCBiBFwIGIEHIgYAQcixnlwFLnggguC9dmzZx837aGHHsoPh65UUu48dznbt28P1tva2iqadjJhCw5EjIADESPgQMQIOBAxAg5EjIADESPgQMQ4Dx6ZcePGBevz588P1m+44YZg/cwzzywaz2azuv/++ytrroyjR48G6x0dHcH6sWPHKpp2MmELDkSMgAMRI+BAxAg4EDECDkSMgAMRI+BAxDgP3oB6n2vuPW3OnL5u+Nqt3HnusWPHVt3XQG3evDlYX7JkSbC+fv36NNs5KZTdgpvZWWa20cx2mNk2M1uQTB9pZhvMLJf8HDH47QLoj0p20Y9Iutfd/0nSv0i628zGS7pPUqu7ZyS1JuMAGkjZgLt7h7tnk+EDknZIGiPpOkmrk6etlnT9YDUJoDrm7pU/2WyspDZJEyTtdPfhBbV97p7fTe/q6sovOJfLpdErgF4ymUx+uKmpyXrXKz7IZmY/kPQnSb9x9/1mxy2rpJ4L8bW3twcvyldPjdRb74Nsr732mmbOnJkfb6SDbNlsVhdddFFFz631QbZG+j8tlGZffd0MslBFp8nMbKi6w/28u7+cTP7SzEYn9dGS9gygTwCDoOwW3Lo31U9L2uHujxaU1ku6XdKy5Oerg9LhCWjUqFHB+vjx44P1xx57rGj88OHDam1tzY+fd9551Tc3QJs2bSoaHzp0aNG0hx9+uOS8r74a/hU52b/aORgq2UW/TNJcSR+b2ZZk2iJ1B/uPZjZP0k5JNw1OiwCqVTbg7v5nSaXecE9Ntx0AaeKjqkDECDgQMQIORIyAAxEj4EDE+LpoCSNHjixZW7lyZXDeiRMnBuvnnntuv3rJZrOpnft+9913g/VHHnkkWG9paSkab2tr09SpfzuZ8u2331bfHFLHFhyIGAEHIkbAgYgRcCBiBByIGAEHIkbAgYhFex78kksuCdYXLlx43LSXXnopPzx58uSS844ZM6b6xlJw+PDhkrUVK1YE5126dGmwfujQoX73w7nvxsUWHIgYAQciRsCBiBFwIGIEHIgYAQciRsCBiEV7HnzWrFn9qmez2bLzVGr79u3B+uuvvx6sHzlypGh8+vTpReevQ9/ZLnenC5xc2IIDESPgQMQIOBAxAg5EjIADESPgQMQIOBAzdw8+JJ0laaOkHZK2SVqQTF8sabekLcljRuF8nZ2d3vOQ5JK8vb09P9xoD3qjtxOxr8Kc9ZXfSj7ockTSve6eNbMfSnrfzDYktT+4+/IKlgGgDiq5P3iHpI5k+ICZ7ZBU30uaAKhIv96Dm9lYSf8saVMyab6ZfWRmz5jZiJR7AzBAlryfLv9Esx9IekfSEnd/2cxGSfpa3e8Ffi9ptLv/quf5XV1d+QXncrlUmwbQLZPJ5IebmprsuCeUO8iW/AEYKqlF0r+XqI+VtJWDbPTWSI9G7a2WB9nK7qKbmUl6WtIOd3+0YProgqfNkrS13LIA1FYlR9EvkzRX0sdmtiWZtkjSHDObqO6/JJ9LumtQOgRQtUqOov9Z0vH79tIb6bcDIE18kg2IGAEHIkbAgYgRcCBiBByIGAEHIkbAgYgRcCBiBByIGAEHIkbAgYgRcCBiBByIWMVXdOmvwiu6ABh8fV3RhS04EDECDkRs0HbRAdQfW3AgYjULuJlNN7NPzOxTM7uvVq9bCTP73Mw+NrMtZra5zr08Y2Z7zGxrwbSRZrbBzHLJz7pcg75Eb4vNbHey7raY2Yw69HWWmW00sx1mts3MFiTT677eAr3VZL3VZBfdzIZI+l9J0yTtktQuaY67bx/0F6+AmX0u6Wfu/nUD9HKFpIOSnnP3Ccm0/5S0192XJX8cR7j7fzRIb4slHaznLaySK/yOLry9lqTrJTWrzust0NvNqsF6q9UWfLKkT939L+7+vaT/lnRdjV77hOLubZL29pp8naTVyfBqdf+C1FyJ3urO3TvcPZsMH1D3jTLHqAHWW6C3mqhVwMdI+mvB+C411v3NXNKbZva+md1Z72b6MCq5R1zPveJ+XOd+emuYW1j1ur1WQ623etz6q1YB7+uyy410+P4yd79I0i8k3Z3siqIy/yXpp5ImqvsmlY/Uq5Hk9lp/kvQbd99frz760kdvNVlvtQr4LnXfZ7zHP0r6okavXZa7f5H83CNpnbrfUjSSL3vuJJP83FPnfvLc/Ut3P+ruxyQ9qTqtOzMbqu4APe/uLyeTG2K99dVbrdZbrQLeLiljZj8xs1Ml/VLS+hq9dpCZDUsOfsjMhkn6uRrvNkzrJd2eDN8u6dU69lKkEW5hVer2WmqA9Vb3W39VcvPBNB6SZqj7SPpnku6v1etW0Ne5kj5MHtvq3ZuktereZfs/de/5zJP095JaJeWSnyMbqLc1kj6W9JG6AzW6Dn1dru63fB9J2pI8ZjTCegv0VpP1xifZgIjxSTYgYgQciBgBByJGwIGIEXAgYgQciBgBByJGwIGI/T/Mu8Gck0mjpgAAAABJRU5ErkJggg==\n",
            "text/plain": [
              "<Figure size 288x288 with 1 Axes>"
            ]
          },
          "metadata": {
            "needs_background": "light"
          },
          "output_type": "display_data"
        }
      ],
      "source": [
        "plt.figure(figsize=(4,4))\n",
        "plt.imshow(images[0].reshape(28,28), cmap='gray');"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "9dRRsmjhqdYt"
      },
      "source": [
        "Before we implement PCA, we will need to do some data preprocessing. In this assessment, some of them\n",
        "will be implemented by you, others we will take care of. However, when you are working on real world problems, you will need to do all these steps by yourself.\n",
        "\n",
        "The preprocessing steps we will do are\n",
        "1. Convert unsigned interger 8 (uint8) encoding of pixels to a floating point number between 0 and 1.\n",
        "2. Subtract from each image the mean $\\boldsymbol \\mu$.\n",
        "3. Scale each dimension of each image by $\\frac{1}{\\sigma}$ where $\\sigma$ is the stardard deviation.\n",
        "\n",
        "The steps above ensure that our images will have zero mean and one variance. These preprocessing\n",
        "steps are also known as [Data Normalization or Feature Scaling](https://en.wikipedia.org/wiki/Feature_scaling)."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "AwBKtnb7qdYu"
      },
      "source": [
        "## 1. PCA\n",
        "\n",
        "Now we will implement PCA. Before we do that, let's pause for a moment and\n",
        "think about the steps for performing PCA. Assume that we are performing PCA on\n",
        "some dataset $\\boldsymbol X$ for $M$ principal components.\n",
        "We then need to perform the following steps, which we break into parts:\n",
        "\n",
        "1. Data normalization (`normalize`).\n",
        "2. Find eigenvalues and corresponding eigenvectors for the covariance matrix $S$.\n",
        "   Sort by the largest eigenvalues and the corresponding eigenvectors (`eig`).\n",
        "\n",
        "After these steps, we can then compute the projection and reconstruction of the data onto the spaced spanned by the top $n$ eigenvectors."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "scrolled": true,
        "id": "X9yyVjF2qdYu"
      },
      "outputs": [],
      "source": [
        "def normalize(X):\n",
        "    \"\"\"Normalize the given dataset X\n",
        "    Args:\n",
        "        X: ndarray, dataset\n",
        "\n",
        "    Returns:\n",
        "        (Xbar, mean, std): tuple of ndarray, Xbar is the normalized dataset\n",
        "        with mean 0 and standard deviation 1; mean and std are the\n",
        "        mean and standard deviation respectively.\n",
        "\n",
        "    Note:\n",
        "        You will encounter dimensions where the standard deviation is\n",
        "        zero, for those when you do normalization the normalized data\n",
        "        will be NaN. Handle this by setting using `std = 1` for those\n",
        "        dimensions when doing normalization.\n",
        "    \"\"\"\n",
        "    mu = np.mean(X, axis=0)\n",
        "    std = np.std(X, axis=0)\n",
        "    std_filled = std.copy()\n",
        "    std_filled[std==0] = 1.\n",
        "    Xbar = ((X-mu)/std_filled)\n",
        "    return Xbar, mu, std\n",
        "\n",
        "def eig(S):\n",
        "    \"\"\"Compute the eigenvalues and corresponding eigenvectors\n",
        "        for the covariance matrix S.\n",
        "    Args:\n",
        "        S: ndarray, covariance matrix\n",
        "\n",
        "    Returns:\n",
        "        (eigvals, eigvecs): ndarray, the eigenvalues and eigenvectors\n",
        "\n",
        "    Note:\n",
        "        the eigenvals and eigenvecs should be sorted in descending\n",
        "        order of the eigen values\n",
        "    \"\"\"\n",
        "    eigvals, eigvecs = np.linalg.eig(S)\n",
        "    k = np.argsort(eigvals)[::-1]\n",
        "    return eigvals[k], eigvecs[:,k]\n",
        "\n",
        "def projection_matrix(B):\n",
        "    \"\"\"Compute the projection matrix onto the space spanned by `B`\n",
        "    Args:\n",
        "        B: ndarray of dimension (D, M), the basis for the subspace\n",
        "\n",
        "    Returns:\n",
        "        P: the projection matrix\n",
        "    \"\"\"\n",
        "    return (B @ np.linalg.inv(B.T @ B) @ B.T)\n",
        "\n",
        "def PCA(X, num_components):\n",
        "    \"\"\"\n",
        "    Args:\n",
        "        X: ndarray of size (N, D), where D is the dimension of the data,\n",
        "           and N is the number of datapoints\n",
        "        num_components: the number of principal components to use.\n",
        "    Returns:\n",
        "        X_reconstruct: ndarray of the reconstruction\n",
        "        of X from the first `num_components` principal components.\n",
        "    \"\"\"\n",
        "    # first perform normalization on the digits so that they have zero mean and unit variance\n",
        "    # Then compute the data covariance matrix S\n",
        "    S = 1.0/len(X) * np.dot(X.T, X)\n",
        "\n",
        "    # Next find eigenvalues and corresponding eigenvectors for S\n",
        "    eig_vals, eig_vecs = eig(S)\n",
        "\n",
        "    # find indices for the largest eigenvalues, use them to sort the eigenvalues and\n",
        "    # corresponding eigenvectors. Take a look at the documenation fo `np.argsort`\n",
        "    # (https://docs.scipy.org/doc/numpy/reference/generated/numpy.argsort.html),\n",
        "    # which might be useful\n",
        "    eig_vals, eig_vecs = eig_vals[:num_components], eig_vecs[:, :num_components]\n",
        "\n",
        "    # dimensionality reduction of the original data\n",
        "    B = np.real(eig_vecs)\n",
        "    # Z = X.T.dot(W)\n",
        "    # reconstruct the images from the lower dimensional representation\n",
        "    reconst = (projection_matrix(B) @ X.T)\n",
        "    return reconst.T\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "lFvHDFK_qdYv"
      },
      "outputs": [],
      "source": [
        "## Some preprocessing of the data\n",
        "NUM_DATAPOINTS = 1000\n",
        "X = (images.reshape(-1, 28 * 28)[:NUM_DATAPOINTS]) / 255.\n",
        "Xbar, mu, std = normalize(X)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "_DcYB7TbqdYv",
        "outputId": "c34fb5c3-babc-4dd7-e42b-9b7f199ffce4"
      },
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "4.8273595646655134e-24\n",
            "5.586566850065507e-25\n",
            "5.530125631353103e-25\n",
            "3.316249326213407e-23\n",
            "1.2015015841330959e-24\n",
            "3.350219923629842e-24\n",
            "1.9496488523062087e-24\n",
            "1.5512187910326764e-24\n",
            "2.0501953750890233e-23\n",
            "1.4003549686952187e-23\n",
            "1.3753195113297822e-22\n",
            "1.3808977146200646e-23\n",
            "1.1964993443957294e-21\n",
            "9.801058172699404e-23\n",
            "8.209428275187677e-23\n",
            "2.3283754737388765e-23\n",
            "1.9457132032354168e-23\n",
            "7.155036400524889e-24\n",
            "4.6872886881379156e-24\n"
          ]
        }
      ],
      "source": [
        "for num_component in range(1, 20):\n",
        "    from sklearn.decomposition import PCA as SKPCA\n",
        "    # We can compute a standard solution given by scikit-learn's implementation of PCA\n",
        "    pca = SKPCA(n_components=num_component, svd_solver='full')\n",
        "    sklearn_reconst = pca.inverse_transform(pca.fit_transform(Xbar))\n",
        "    reconst = PCA(Xbar, num_component)\n",
        "    np.testing.assert_almost_equal(reconst, sklearn_reconst)\n",
        "    print(np.square(reconst - sklearn_reconst).sum())"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "8H-6S_KEqdYw"
      },
      "source": [
        "The greater number of of principal components we use, the smaller will our reconstruction\n",
        "error be. Now, let's answer the following question:\n",
        "\n",
        "\n",
        "> How many principal components do we need\n",
        "> in order to reach a Mean Squared Error (MSE) of less than $100$ for our dataset?\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "kywvHGQAqdYw"
      },
      "source": [
        "We have provided a function in the next cell that computes the mean squared error (MSE), which will be useful for answering the question above."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "TyDfTEO9qdYw"
      },
      "outputs": [],
      "source": [
        "def mse(predict, actual):\n",
        "    \"\"\"Helper function for computing the mean squared error (MSE)\"\"\"\n",
        "    return np.square(predict - actual).sum(axis=1).mean()"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "pLag5ZdgqdYw"
      },
      "outputs": [],
      "source": [
        "loss = []\n",
        "reconstructions = []\n",
        "\n",
        "# iterate over different numbers of principal components, and compute the MSE\n",
        "for num_component in range(1, 100):\n",
        "    reconst = PCA(Xbar, num_component)\n",
        "    error = mse(reconst, Xbar)\n",
        "    reconstructions.append(reconst)\n",
        "    # print('n = {:d}, reconstruction_error = {:f}'.format(num_component, error))\n",
        "    loss.append((num_component, error))\n",
        "\n",
        "reconstructions = np.asarray(reconstructions)\n",
        "reconstructions = reconstructions * std + mu # \"unnormalize\" the reconstructed image\n",
        "loss = np.asarray(loss)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "MPDPmIzQqdYw",
        "outputId": "872b706d-cb89-4c45-9a6c-0df64d001ab4"
      },
      "outputs": [
        {
          "data": {
            "text/html": [
              "<div>\n",
              "<style scoped>\n",
              "    .dataframe tbody tr th:only-of-type {\n",
              "        vertical-align: middle;\n",
              "    }\n",
              "\n",
              "    .dataframe tbody tr th {\n",
              "        vertical-align: top;\n",
              "    }\n",
              "\n",
              "    .dataframe thead th {\n",
              "        text-align: right;\n",
              "    }\n",
              "</style>\n",
              "<table border=\"1\" class=\"dataframe\">\n",
              "  <thead>\n",
              "    <tr style=\"text-align: right;\">\n",
              "      <th></th>\n",
              "      <th>0</th>\n",
              "      <th>1</th>\n",
              "    </tr>\n",
              "  </thead>\n",
              "  <tbody>\n",
              "    <tr>\n",
              "      <th>0</th>\n",
              "      <td>1.0</td>\n",
              "      <td>569.447737</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>1</th>\n",
              "      <td>2.0</td>\n",
              "      <td>536.059608</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>2</th>\n",
              "      <td>3.0</td>\n",
              "      <td>508.250286</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>3</th>\n",
              "      <td>4.0</td>\n",
              "      <td>487.018907</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>4</th>\n",
              "      <td>5.0</td>\n",
              "      <td>467.571610</td>\n",
              "    </tr>\n",
              "  </tbody>\n",
              "</table>\n",
              "</div>"
            ],
            "text/plain": [
              "     0           1\n",
              "0  1.0  569.447737\n",
              "1  2.0  536.059608\n",
              "2  3.0  508.250286\n",
              "3  4.0  487.018907\n",
              "4  5.0  467.571610"
            ]
          },
          "execution_count": 10,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "import pandas as pd\n",
        "# create a table showing the number of principal components and MSE\n",
        "pd.DataFrame(loss).head()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "zA46y9lZqdYw"
      },
      "source": [
        "We can also put these numbers into perspective by plotting them."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "7kJJ1nYnqdYw",
        "outputId": "7e31609a-3fab-4104-b236-88f15b3f322b"
      },
      "outputs": [
        {
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAa0AAAEXCAYAAAAeBxskAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3deXhU5f3//+c7yWQhOyTsICBBQRFUUFzqvqC1Ra229tO69Ke1/dYu1tpqW9taW612Ubtv2lZt646CSt1QVFxBBUQRwr6TQFayL/fvj3MCM5NJMhOSyfZ6XNdcmXPPmfu8z5mTec855z73bc45RERE+oKEng5AREQkWkpaIiLSZyhpiYhIn6GkJSIifYaSloiI9BlKWiIi0mcoaUm/YGaLzOyeno4jnJl9w8y2mlmzmd3cDfVvNLOburjObtmW3RGrDDxKWnFiZv8yM2dmj0d47Xz/tcaw8k+b2WIzKzGzKjNba2b/MbMs//Vx/vsiPa6P17pJZGY2Ergb+AUwCvh1NyxmJnBXF9d5IXBdF9cpXcz/Pri5p+OIt6SeDmCA2Qx8ysyGOed2BZVfDWwCRrcUmNlpwFzgZ/7rdcBE4HwgJazeOcA7YWUVXRv6wGRmyc65+k6+fQLeD8P5zrkdXRjWvricc8VdWS+Ac66kq+sU6TLOOT3i8AD+BbwIvArcEFQ+FmgAfgI0BpXfDSztoM5xgANOjCGOM4EmYExY+eeAWiDHn/4BsB4vWRYDzwFp7dS7EbgF+C1QAuzCO7JIDJpnEXBP2PtuAjZG2E7fALYCe4F7gADwVbzkXgr8DUgOq/sfwO3AbrykfU94zH69H/vrWgj8EEgKW4+fA38C9gBL2lnnc4F3/W1U5L8n3X/tZv+zCX6Ma2fb3erHW+HHfweQ0FFcfvlNsXwO/nzXAB8Fxf5YW59TNNsWb79a5C+zHHgFOCbCet7U1vb05zkYeNSvpxpYAZwXzTbv7v3Hr+N2YBtQ72+//wuL3wFfAx4AKoEtwPfC5kny948NePvhh8BXYqnHj7fV/uXHeKe/7nXADuCheH/fdeejxwMYKI+gf6Yv4n1Zml9+C/AscAWhSesGoCz8Hz+sznHEnrQS/B36+2HlTwMP+88v9P9pP4WXVKcD19Jx0ioFbgQK8JJgI/CloHkWEV3SKgfuAyYDn/b/sRcA9wNTgPOAGuD/hdVdAfzdf9+n8L7Ufhc0z814X1oXAOPxvgA3Az8LW48Kf95JwJQ21vcIf/3u8pd3jl/XA/7rGf52dMCRwHDCEkeEZd4CHAJcClQB13UUF5GTVkefw0/xvsy/7td1VFgdIZ9TlNv2AuBiv77D8L7wS4AhYbG1mbT8bbQL7//kRLwENgc4N5ptHof951d4Pxha1vMHQDNwetA8zl+HL/vxf9MvOzUsxhXAWXj74efw/tevjLYeYDBe0vu1v92GA4l4p3W3Aqfg/e/OBK7t6e+/rnz0eAAD5cH+pJXq7/in+jvZVrwvtysITVqDgPn+jroDeBL4VtiXwDj/9Wq8L6Hgx3HtxHI78FHQ9FC8o71P+tPfBtYAgRjWbyPeabDgsmeBB4OmFxFd0ioi9FfwM3i/flOCyubR+uhgI6FHdi2nVdP97VkNzA5b/mVAWdh6LIxifR8A3gkrm4P3JXaQP32K//mMjmLbvRZWdhuwtaO4iJy02vwc/G1RA1zfTjwhn1NH27aNOhLwkucX2oo1wnt+Buxsp85otnl37j91wNfClv8E8FLQtCMo0fllHwO/8J+P9+M9NGyeHwPLoq3Hn14L3Bw2z2+Bl/B/FPfHhxpixJlzrhbvn+/LwCfxThU8FWG+aufcp/F28u8D2/2/q81sctjsX8I7Ggp+vN9OGPcBk81spj/9ebxE+pw//QjeaYZNfgOSS80sM4rVWxY2vQ0YFsX7wq1yodeRdgKrnXN1YWVDw973jnOuKWj6dSAZ75fqYUAa8LiZ7W15AH8Fss0sP7ieKGI8DO9Ub7BXAMP7NR+rN8OmXwdGtTS6iSEuaP9zOAzvh9PzMcbX3rbFzMab2QN+44AKvKOWbOCgGJZxNPCGc66qjdej3ebdsf9M9J9HWv5hYWXtbf8ZfrxLw/bDH+AdGUdbT1v+CUwF1prZX8zsM2aW3MF7+hQ1xOgZf8VLKmOBfzrnGsws4ozOuY14vx7/ZWY/xDsC+h5eomqxzTm3NtqFO+dWmdlSvKOMJf7f/zrnGv3Xt5nZoXhHg6cBPwLuMLNjnXNb2qk6vMGCI7SFajPeP2ywQIR6GiLUE6msox9dwctqmfdivG0YLrjxQVtfmuFcjOWxiLRDRBtXR59DS9mBCI/vabyjmWvwrr/UA4vxvuhj0VFc0Wzz7th/2lq+RShrb/u3/D0e78i/vbqj+RxDZ3BumZmNx7vGeCrekdfPzGyWc65fNM7SkVYPcM6twksWx+Od+4/2faVE/oXYGfcDl5jZNLxrGveFLavOOfesc+57eL/cBuG1XDwQRcDIsLKjDrDOYDPNLDFo+ji8f/x1eBe7a4EJzrm1ER5NkSpsx4fAyWFlJ+N9sXzUidhnhU0fB2zvhi+aj/C2w9kxvq/NbWtmQ/COdG53zj3nnGtZRqz76bvACWaW3sbrXb3Nw7W3/6zFOz0YvvyT/Lii9a7/d2yEfXBdjPHW411iCOGc2+uce8I59028I7vJEeLus3Sk1XPOBlJdG82L/fsvMvDOx2/0n18OHI736ynYYDMbHlZW5ZyrbGf5DwK/wb8o7JxbHrTsK/F+0LyDd4H4dCCTA/9ieBH4s5l9FngPuAj4hL+MrjAE+KOZ/RavufnPgL+3nG4ys9uA2/yj2hfw9v+pwJHOuRtiXNavgPfM7E68lmjjgN8D/3HObe5E7NP9z/y/eF8038JrdNGlnHN7zew3wM1mVoO3HdLwGjv8op23trlt/XqKgS+b2Tp/3l/iXTuLxZ+ArwDzzOwneKfEDwOanHP/o+u3ebiO9p/f4R21FOOdursY75ramdEuwDm31sz+AfzdzL6Hd1o4He/UaL5z7o4Y4t2Al+TH4h21lQDfwdtuy/yyz+O1Fo50dqFPUtLqIc65alqfHgj2Cl4T3X/itQzaC6wGvuic+0/YvPMivP+PeK3D2lr+bjN7Bu/oKfxG5FK/7Jd494StB652zi1sJ95o3IeXdP+Ad9roP8Dv8E5PdoXH8JoHt5yWehTvVCoAzrmfmdl2vObQv8b7Ul2Dl7hj4pxbYWafxvtiuwbvGs5jtN6W0fo93vWfpXgt5P5M19803OJHeEnmm/4ySml9rSZcm9vWOddsZhfjfZYr8Fpo/gCv2X7UnHM7zOxE/30L8E4dF+Jdy+2ObR6u3f0H7/aIZrzbUfLxjr6+2In/i6vxkssP8ZJjBd7R2h9irOcneJcaVuNdpxzv13Ud3vWxBGAV8Bnn3OoY6+61Wppdi0gPMbONeK31ft7TsURiZouAtc65q3o6lu4yENaxv9A1LRER6TOUtEREpM/Q6UEREekzdKQlIiJ9Rp9tPVheXq5DRBGRfiw7O7vVDd460hIRkT5DSUtERPqMAZ+0CgsLe0UdXVVPb6mjq+pRLL27jq6qR7H07jq6sp4DNeCTloiI9B1KWiIi0mcoaYmISJ+hpCUiIn2GkpaIiPQZSloiItJnDNik1dDsWFJUz7+2JHHXivbGShQRkd6iz3bjdCDWlDVw6lPFVDU6IJmhRXu5dmoG/oi2IiLSSw3II63xWUkEp6eimmbWlDf2WDwiIhKdAZm0AgnGrGHJIWWLd9b1UDQiIhKtAZm0AD4xIiVk+rUd9T0UiYiIRGvAJq0Th4cmrcU769CAmCIivduATVrThgTIDOy/srW7tpnVuq4lItKrDdiklZRgHBd2Xeu1HbquJSLSmw3YpAWRTxGKiEjvNaCTVnhjjMU76mnWdS0RkV5rQCetqYMDpCfuT1J76pr5uEzXtUREeqsBnbSSEowjs5pDyhbrupaISK81oJMWwNHZTSHTr+m6lohIr6WklROatF7fqetaIiK91YBPWpPSHVnJ++/XKqlrZsWehh6MSERE2jLgk1aiwUlhTd//t6W2h6IREZH2DPikBTB7bGrI9LNKWiIivZKSFnD26NSQoUqW72lgW1VTm/OLiEjPUNIC8tMSmZkf2qXTczraEhHpdZS0fK1PEdb0UCQiItIWJS3f7DGhSeuVHXVUNTS3MbeIiPQEJS3f5JwkxmYk7puua4JF23WjsYhIb6Kk5TMzzgk72lLTdxGR3kVJK8g5Yde1nt9aq94xRER6ESWtIMcPSyEraDTjoppm3tut3jFERHoLJa0gyYnGaaNCj7ae2aRWhCIivUXckpaZbTSzD8xsmZkt9csGm9kLZlbo/831y83Mfmdma81shZkdFa84w08RPrSumsZmnSIUEekN4n2kdapzbrpzboY/fSOw0DlXACz0pwHOAQr8x9XAn+MV4LljU8lI2n+KcEd1M89vVYMMEZHeoKdPD84B7vOf3wecH1R+v/O8BeSY2Yh4BJQZSOAzE9JCyu5bUx2PRYuISAfimbQc8LyZvWtmV/tlw5xzOwD8v0P98lHAlqD3bvXL4uLySekh0y9srVVfhCIivYC5ODXpNrORzrntZjYUeAH4BjDfOZcTNE+pcy7XzJ4BfuGcW+yXLwS+55x7t2Xe8vLyfYEXFhZ2aazOwReWpVJYtT+nf2VsPVeNbezS5YiISKiCgoJ9z7Ozsy389aR4BeKc2+7/LTKzJ4BjgF1mNsI5t8M//Vfkz74VGBP09tHA9rbqDl7JWBUWFkZ8/9VNe/nuW+X7pp/Zk8Ztpw4jMaHVNmyzjq6KpS/WoVh6fyz9bX0US/fV0ZX1HKi4nB40s3Qzy2x5DpwFrATmA5f7s10OzPOfzwcu81sRzgLKW04jxsvFEwaRlrg/QW2tauJldeskItKj4nVNaxiw2MyWA+8AzzjnngVuB840s0LgTH8aYAGwHlgL/B34Wpzi3CcnJYE540Kbv9+3pireYYiISJC4nB50zq0HpkUo3wOcHqHcAdfEIbR2XX5IOg+t239z8f8217KjuokRgxLbeZeIiHSXnm7y3qvNGprMIdn783qjg3s/1tGWiEhPUdJqh5nx5cmhzd//tbqK2kb1kCEi0hOUtDpwycRBZCXvb5Cxu7aZxzfoZmMRkZ6gpNWBjEAClxaEHm399aMq4nV/m4iI7KekFYUvT04n+O6sFSUNvLmrvsfiEREZqJS0ojAuM4nZYaMa/3XV3h6KRkRk4FLSitJXp2SETD+9qZYte9Wtk4hIPClpRemkEclMztnf/L3JwZ8+1NGWiEg8KWlFycz4StjR1j9XV7GjWr2/i4jEi5JWDC45eBAjB+3fZLVNcOfyyh6MSERkYFHSikFqkvGdaZkhZfetqWJnbeue30VEpOspacXo0oJ0xmTs73uwvhnu3RLowYhERAYOJa0YJSca3w072nqqKJENFWpJKCLS3ZS0OuHzEwcxPnP/0VaTM+5YVtGDEYmIDAxKWp0QSDBumJ4VUvbwuhpW7FEvGSIi3UlJq5MunpAWMmyJA37wTrn6JBQR6UZKWp2UmGD8dGbo0dbinfU8tam2hyISEen/lLQOwNmjUzltZEpI2Y+WlFPXpKMtEZHuoKR1AMyMW4/JJpH9SWrT3ib+rO6dRES6hZLWAZqcG+DCEaHN3X+9vJJd6t5JRKTLKWl1gavHNpATNLrx3kbHj5eW92BEIiL9k5JWF8gJwI1Htm4C/8r2uh6KSESkf1LS6iJXHprOYblJIWXfebNMjTJERLqQklYXCSQYdx2fQ3DXuWsrGrn7A/UCLyLSVZS0utAxQ1O44pBBIWV3rqhkXbn6JRQR6QpKWl3sJ0dnk5+6f7PWNcG33yyjWT1liIgcMCWtLpaTksCtx2SHlL26o46/rarqoYhERPoPJa1ucPGENE4J6ynjJ0vLWVXa0EMRiYj0D0pa3cDM+P0JOWQF3btV1wRXv1pKvVoTioh0mpJWNxmTkcRvZuWElH1Q0sBt72vcLRGRzlLS6kYXHzyIz4xPCyn77Qd7eW2HbjoWEekMJa1u9pvjchg1aP8oxw64+tUS9tSqb0IRkVgpaXWznJQE/vSJ3JCbjndUN/O1xWUaMFJEJEZxTVpmlmhm75vZ0/70eDN728wKzexhM0v2y1P86bX+6+PiGWdXO3lkCtdOzQgpe25LLX/5SM3gRURiEe8jrW8Bq4Km7wDucs4VAKXAlX75lUCpc24icJc/X5/2g6OymJkfCCn78dJylu2u76GIRET6nrglLTMbDXwSuMefNuA04DF/lvuA8/3nc/xp/NdP9+fvswIJxj0nDw5pBt/QDJe9XMJuXd8SEYmKxeu6ipk9BvwCyASuB64A3vKPpjCzMcD/nHOHm9lKYLZzbqv/2jrgWOfc7pb6ysvL9wVeWFgYl3XoCgt3J3Ljx6E3Hk/PauKPh9eRrCuMIjLAFRQU7HuenZ3d6mAlKbygO5jZeUCRc+5dMzulpTjCrC6K11oJXslYFRYWHtD7Y62joAAKKePej/dfz1pWkcgfivL49rDdTJoUv1i6sw7F0vtj6W/ro1i6r46urOdAxeu3/QnAp81sI/AQ3mnBu4EcM2tJnKOB7f7zrcAYAP/1bKAkTrF2u9uPzeYTw5NDyh5cW8392+LyG0JEpM+KS9Jyzn3fOTfaOTcOuAR4yTn3BeBl4CJ/tsuBef7z+f40/usvuX7UPjyQYNx/2hAOzkoMKf/jxgDzN9b0UFQiIr1fT19FuQG4zszWAkOAe/3ye4Ehfvl1wI09FF+3yU1J4OEzhpAd1DDDYVz9aglLi9WiUEQkkrgnLefcIufcef7z9c65Y5xzE51zFzvn6vzyWn96ov/6+njHGQ8TswPcf+pgkoKu4NU2wedf3MOmSg0cKSISrqePtAa8k0emcufxoR3rFtc287kX91BW19xDUYmI9E5KWr3AZZPS+XZYjxkflzXyhZf2UN2oxCUi0kJJq5f40dFZnJEXekrw9Z31XPZSCXUag0tEBFDS6jUSzPhJQT3H5Ic2hX9xWx1XvVJCY7MSl4iIklYvkpoID585hCm5ofdrPbWplmsWl9Lcf1r9i4h0ipJWL5ObksCTZ+cxMSs0cT28roZvv6HhTERkYFPS6oWGpiUyb3YeYzNCbz6+b001N7xdrsQlIgOWklYvNSo9kXln5zFiUOhH9LdVVfx4aYUSl4gMSEpavdj4rCTmnZ1Hfmrox/T7lXu59b1KJS4RGXCUtHq5STkB5s3OY3BK6Ef16xWV3Pa+EpeIDCxKWn3AlNwAc88K7acQ4FfLK/nFssoeikpEJP6UtPqI6XnJzD0rj6xAaOL65bJKfvF+RQ9FJSISX0pafcjR+cnMPbt14rpjWSW3vKtWhSLS/ylp9TEz8pN5/Kw8MsMS150r9nLj2+Wo4wwR6c86TFpmdnHY9CFh09d2dVDSvplDI58q/OuqKm5bm0yTMpeI9FPRHGndGzb9Ztj0LV0Ui8Rg5tBk5s3OIzclNHHN25XEl18tpV6d7IpIPxRN0rIYpyVOjsxL5unZ+QxNC/0Y526o4fML91DVoGFNRKR/iSZphf9k72ha4uiwwQEWnJPHqEGhXT4t3FbHBc/toVQDSYpIPxJVQwzzJJhZYqRp6VkTswP875N5HJwV+nG8U1zPuQuK2VbV1EORiYh0rWiSVgbQCDQA9UBO0HQDkN5t0UnUxmYk8ey5+RySHnpktaqskbOeLuaj0oYeikxEpOtEk7TGAxOCHuMjPJdeID8tkb9MreX4YaEDSW6rbmL2gmIW76zrochERLpGh0nLObcp0gOoCHouvURGEjx+Vh6fHJsaUl5R77jwud08tr66hyITETlw0dyndZmZnR00PcPMtgC7zWx1+H1b0vPSkoz7Tx3MVYeGnrmtb4arXinl5+9VaBRkEemTojk9+B1gZ9D034AXgSP8v7/qhrjkACUmGL+alc2Pj85q9dqvl1dyxcslahIvIn1ONElrLPABgJmNAaYC33HOfQjcCBzbfeHJgTAzrjsik798IpfksE96/qZaZi/YzZa9jT0TnIhIJ0STtBqBliv7xwMfO+dK/OlqIK07ApOuc8nEQcyfnUde2GCSH5Q0cOpTxbyuBhoi0kdEk7ReAW41syOAbwBPBb12KKGnDqWXmjUshZc+lc9huUkh5btrm5nz7G7uWbVXvcSLSK8XTdL6FnAk8DrekdUdQa9dCjzbDXFJNxibkcRzn8xv1bKw0cH1b5VzzeIyqht1nUtEeq+kjmchEbgCr49BB2SbWbb/2p+6KS7pJhmBBB44bTC/XFbJ7WGjHv93bTXL99RzywSjoIfiExFpTzRJayOh/QuGd5Dr8BKb9BEJZtx4ZBaHDw7w1VdL2du4/+P9sLSRy5al8uesGj49TpcrRaR3ieb04AqgELgJGAcEwh7Jbb5TerXzDkrjhfPyKcgO/e1S1WRc9nIJN75dpiFORKRXiaZHjOnARcBgYDGwALgESHbONTnn1BtrHzY5N8BLn8rngghHVX/5qIrZC4rZWKlm8SLSO0TVy7tzbqVz7rt4fQ3eCZwH7DCzo6J5v5mlmtk7ZrbczD40s5/65ePN7G0zKzSzh80s2S9P8afX+q+P68zKSXQyAwn845Rc7jg2m0DYHvHe7gZOml/EvI01PROciEiQqJJWkALgZOA44H2gNMr31QGnOeemAdOB2WY2C68l4l3OuQK/riv9+a8ESp1zE4G7CG2xKN3AzPjKlAz+d24+I1JCWxBW1Dsuf7mEb75eql40RKRHRdP34GAzu8bM3gGeBPYCJznnTnXObYhmIc6z159suRbmgNOAx/zy+4Dz/edz/Gn81083M42QHAcz8pP595G1nBvWLB7g/jXVnDy/mGW763sgMhGR6I60tgNfx0tY1wBvARPN7LSWRzQLMrNEM1sGFAEvAOuAMudcywWTrcAo//koYAuA/3o5MCS6VZIDlZUE/zltMLcd0/p04dqKRs58ppg7V1TS1KxGGiISX9ZRLwhmtpHQJu/hnHMu6jG1zCwHeAL4MfBP/xRgS7+GC5xzU83sQ+Bs59xW/7V1wDHOuT0t9ZSXl++LqbCwMNrFS4xW7zV+uDqFTTWtf99Mz2ri5kn1jEpV8hKRrlFQsP8u0ezs7FZn2Dq8T8s5N64rA3LOlZnZImAWkGNmSf7R1Gi8ozrwjrrGAFvNLAnIBkoi1QehKxmrwsLCA3p/V9XRm2IJrqMAOPXwZn7wTjn3rQkdi2tZRSKXLh/Ebcdk88WCQYSfwe0t66NYuq8OxdL7Y+lN69MVYm2I0Slmlu8fYWFmacAZwCrgZbzm9ACXA/P85/P9afzXX3LqGK/HpAcS+O0Judx/6mAGp4TuMpUNjm+8XsYlL+5hZ7XufhCR7hWXpAWMAF42sxXAEuAF59zTwA3AdWa2Fu+a1b3+/PcCQ/zy6/CGQJEe9ulxabxx/lDOHJXS6rXnttYx64ldPLquWh3viki3iaYbpwPmnFuB1+luePl64JgI5bXAxXEITWI0fFAij5w5hHs/ruLHSyuoDuoCqqze8eVXS3l0fTW/OS6nB6MUkf4qXkda0o+YGVdNzmDxnKHMGtq6F6/nt9Yx64kiHtqepBaGItKllLSk0yZkJfHMOXn8bEYWKWFdJlc1On6zPpnTn9Z9XSLSdZS05IAkJhjfmJrJa58eynHDWh91LdvTwGlPF/Pdt8oor1dvGiJyYJS0pEtMygnwzDl53HVcDlmB0KbvzQ7+vqqKmXN38dBaNdQQkc5T0pIuk2DGlw5N5+0Lh0XsNb6oppmvvlbKOQt2s2KPThmKSOyUtKTLjRiUyD9PHcxvD6tlXGbr8UHfKqrnlKeKufb1UoprdG+XiERPSUu6zfG5zbx5/jC+Oy2zVUONZgf/WlPN0Y/v4vcfVFKnwSZFJApKWtKt0pKMHx6VxVvnD2P2mNY9x1c0OH60tIJZT+ziqU01ut4lIu1S0pK4GJ+VxENnDOGRM4YwMav1Pe0bKpu49KUSznt2N8t1vUtE2qCkJXF11phU3jh/KLcek01Wcush0l7fWc8p84u5+tUSNu9tjFCDiAxkSloSd8mJxjWHZfD+Z4Zx5aHpJITlLgc8sq6GGY/v4u4NAUrrdH+XiHiUtKTHDElN5DfH5fD6nKGcHqET3vpm+M+2ANMe28lvlldS1aDkJTLQKWlJj5ucG+Dxs/J49MwhTMlpfb2rot7xs/cqOPLxXfzto73Uq6WhyIClpCW9xpmjU3ltzlD+cGIOIwe13jWLapr53tvlzJi7i/8WVqkzXpEBSElLepXEBOOLBem8+5nh3Hx0FpmJrRPT5r1NfG1xGcc9WcSTG2poVjN5kQFDSUt6pbQk49ojMpk3s4bvHJHBoKTWLQ3XlDdyxaISPjGviKd1j5fIgKCkJb1aZhL86Ohs3v/MMK6enE4gwh77YWkjX3yphFOfKmbBZiUvkf5MSUv6hGGDEvnlrByWXjiMLxQMatVMHrxhUP5vYQmfmF/MvI06bSjSHylpSZ9yUGYSfzwxl7cvGMpFE9KIkLtYWdLA5S+XcMKTRTy6rppGNdgQ6TeUtKRPKsgOcM/Jg3n9/KF86qDWfRoCrCpr5MuvlnLM3F3M25mopvIi/YCSlvRpU3IDPHDaEBbPGcoF4yIfea2vbOLna1OY/thO/rCykkrdpCzSZylpSb9w+OAA/zx1MG9eMJTPTkiLeM1re3UzNy2pYOojO/n5exUUaSwvkT5HSUv6lUNzAvzt5MEsvXAYlxYMIkJLecrqHb9eXsnUR3dy7eulrC1viH+gItIpSlrSL03ISuL3J+by/kXD+MrkdFISWl/PqmvyBqKcObeI/1u4hzd31am5vEgvp6Ql/dqYjCTumJXD0zNr+N70THIiDIfigAWbazlnwW7OfKaYueuraVCLQ5FeSUlLBoScAPzgyCxWfnY4tx+bzZiMxIjzLS1u4P97pZQjHt3Jr5dXsrtW171EehMlLRlQMgIJfHWKN5bXPSfnMm1IIOJ8O6qb+fl7FUx5eCdfebWEJUX16MyhSM9T0pIBKSnBuGjCIBZ9Kp/5s/M4a3Tr8aUjArcAABreSURBVLzAG9Pr4XU1nPlMMZcuS+W+1VXsVZN5kR6jpCUDmplx0ogUHjkzj7cvGMqVh6ZH7JwXYHVVAt96o4zJD+/k+jfLWFmiVoci8aakJeI7JCfAb47L4aPPDudnM7MYlxn5uldlg+Oej6s4cV4Rpz9VxH2rq3TDskicKGmJhMlJSeAbh2fy3meG8diZQzh7TGrEnjYA3t3dwLfeKOPQh3by9cWlvFOkZvMi3an12OYiAkCCGWeMTuWM0alsqmzk7re28syeFIpqWh9VVTU6/l1Yzb8LqzkkO4kvThrEJQcPIj8t8tGaiHSOjrREonBQZhJfG9fAyouH889TcjllZOSGGwCryxv50ZIKJj+8k0tf2sPzW2rV07xIF4lL0jKzMWb2spmtMrMPzexbfvlgM3vBzAr9v7l+uZnZ78xsrZmtMLOj4hGnSEeSE40Lxg/iybPzWHbRMK6flsmIQZH/jRodPLWpls++uIepj+7klnfL1WWUyAGK15FWI/Ad59xkYBZwjZlNAW4EFjrnCoCF/jTAOUCB/7ga+HOc4hSJ2rjMJG46KosPLh7Ow2cM4byxqRH7OgTvvq87V+xlxtwirlqRwn2rqyivV+MNkVjFJWk553Y4597zn1cCq4BRwBzgPn+2+4Dz/edzgPud5y0gx8xGxCNWkVglJRhnj0nl36cP4aPPDeeWGVkUZLd9uXh5RSLfeqOMQx7awZWLSnhhq04fikTL4t3SyczGAa8ChwObnXM5Qa+VOudyzexp4Hbn3GK/fCFwg3Nuacu85eXl+wIvLCyMU/Qi0XEOVlQmMH9XEi/uTqS6qa32h57BAcfs/EbOHdrIpHSHtT+7SL9VUFCw73l2dnar/4S4th40swzgceBa51yFtf2fGemFNrNr8ErGqrCw8IDe31V19KZY+tv69FQsk4CLgL0NzczbWMO/C6t5c1d9xHlLGoz/bg/w3+0BDs1J4nMHD+KiCWmMyYj8LzrQt61iiW8dXVnPgYpb60EzC+AlrP845+b6xbtaTvv5f4v88q3AmKC3jwa2xytWka6UEUjgCwXp/O/cfN7/zDCuGtPA2DY67AX4uKyRn75bwdRHd3HOgmLu/XivOu4V8cWr9aAB9wKrnHN3Br00H7jcf345MC+o/DK/FeEsoNw5tyMesYp0p/FZSXzloAaWXTSMp8/J49KCQWQF2j4X+Oauer7zZjmHPLSTzzy/m38XVlFWpwYcMnDF6/TgCcClwAdmtswv+wFwO/CImV0JbAYu9l9bAJwLrAWqgS/FKU6RuEgw48ThKZw4PIVfzsrhf5treHhdNS9uq6MpwonwJgcLt9WxcFsd304o49jsFC5NqOacMalkJet2Sxk44pK0/AYVbf2cPD3C/A64pluDEukl0pKMCycM4sIJg9hd28Tc9V4Ce3d35Hu6GpphcWkii18tJSURzhiVygXj0zh7TCqZASUw6d/UjZNIL5KXmsjVUzK4ekoGGyoaeXxDDY+tr+bjssaI89c1wTOba3lmcy0piXD6qFTmjEtj9phUsnUEJv2QkpZILzU+K4nrp2Vy/bRMPiptYO6GGp7YUM26isiNMuqaYMHmWhZsriU5AU4dlcqcg1I5d2waOSlKYNI/KGmJ9AFTcgNMyQ3wwyMzWVHSwD/e286islQ27Y2cwOqb4bkttTy3pZYkK+MTI1L41EFpnDs2leGD1Imv9F1KWiJ9iJkxbUgyXx/XwN0TD2LZngae3FDDkxtr2kxgjQ5e3l7Hy9vruO5NmJkf4NyxaRzujJ6/60YkNkpaIn2UmXFkXjJH5iVz84wslu9pYN7GGuZtrGF9Zdv3dS0pbmBJcQOQxsFrd3L2mFRmj0njuGHJBBLUFYf0bkpaIv2AmTE9L5npecn8+OgsVpY2Mm9DDU9tqmF1eeRGHADrKpr404dV/OnDKrKSjTNGpXL2mFTOGJXCkFSdRpTeR0lLpJ8xM6YODjB1cICbjs5iTVkDT2+u5alNNbzfRjN6gIp6x9wNNczdUEOCwYy8ZM4cncKZo1M5YkiABHWIKL2AkpZIPzcpJ8B1OQGuOyKTbVVN/G9zDc9sruW1HbU0usiJqNnBO8X1vFNcz63vVzIsLYHTRqVy5qgUThuVqtaI0mOUtEQGkFHpiVw1OYOrJmfw/qpCNqeN5rkttTy/tZbdtW13D7WrppkH11bz4NpqEgxm5idzxqgUzhidyiCNqiJxpKQlMkBlJMGccWnMGZdGU7Pjvd0NPLe1lme31LKypO3TiM0O3i6q5+0i7ygsJymNM3aUcOrIFE4dmcrIdF0Lk+6jpCUiJCYYM4cmM3NoMjcdlcX2qiZe3FbL81tqWbS9jr2NbR9OlTUaj62v4bH1NQAcmpPEKSNTOGVkCicMT1HXUtKllLREpJWR6YlcNimdyyalU9/keKuonhe31vLi1lo+aqNLqRYflzXycVkjf/moiiSDGfnJnOwnsRn5alYvB0ZJS0TalZxonDQihZNGpHDLzGy2VTWxcFstC7fV8vL2Oirq2z4Ka3TwVlE9bxXVc8eySjKSjOOHJ3PyyFROHpHClNwktUqUmChpiUhMRgUdhTU2O558fz2FCfm8tK2Wd3c30NxOw4y9jY7nt9bx/NY6AAanJHDi8GQOSUziwrwGDs1Jop0RzUWUtESk85ISjGlZzVxUkMX3j8yirK6ZxTvreMXvNmptRfunEkvqmpm/qRZI5lfri8hLTfDHGUvmxBEpHJKtJCahlLREpMvkpCRw3kFpnHdQGgBb9zbyyg4vib2yo45dNe2Pury7tpknN3p9KQLkpSZw/LBkThjuNerQ6URR0hKRbjM6I4kvFCTxhYJ0nHN8XOYlsUXb63hjV/vXw8BLYvM31fpHY5CbYhw3LIXjhyVz/LAUjhgSIEkNOwYUJS0RiQszY3JugMm5Ab46JYOmZseKkgZe21HHs+tKWVGZ1G7TeoDSOrdvzDCA9CSvqf5xw5I5blgKuW33Eyz9hJKWiPSIxIT9vdSfk7qT8QePYfkeL4m9vrOON3fVd5jEqhodi7Z7R25QSaKlceTaIo4blsKxQ5OZNSyZPHX8268oaYlIr5CUYBydn8zR+clce0Qmjc2O5XsaeGNnHYt31fNmFKcTm5yxtLiBpcUN/N4vOzgrkWOHpjBrWDIz85M5JEfXxfoyJS0R6ZWCk9g3pkJTs2NlaQNv+gnsjZ31FLfTX2KLdRVNrKuo5r9rqwHISjZm5nsJ7JihyRyVl6wOgPsQJS0R6RMSE7xRm6cNSearUzJwzrGuopE3dtXvS2Qb2xn8skVFvWPhtjoWbvPuFTNgUnYSM4YmMyMvmaE1xvhmpwYevZSSloj0SWbGxOwAE7MDXDYpHYA3PyykKH00b+ys4+2iej4oaaCpg17oHbC6vJHV5Y38p7AaSGPQyh1MGxJgRn4yR+clc1R+gDHpibpnrBdQ0hKRfiMvGY7ze64H2NvQzLvFDbxdVMeSIm98sPIOrosBVDc6/+itfl9ZfmoCR+Unc1RewEtkeQEGq5FH3ClpiUi/lRFI4OSRKZw8MgWAZudYU97IkqJ6lhTXs7SonlVljUQzJFhxbTPPbanluS21+8oOykjkKD+BHZmfzBGDA2Ql6/pYd1LSEpEBI8GMQ3MCHJoT4FL/lGJ5fTPLdteztLiBJcX1vLOzhpKG6E4DbtrbxKa9NTzh9+BhQEF2EtPzAoxqTuKsrDqmDg6QoeFZuoySlogMaNnJCV6v8yNTAVizppDUkeN5t9hLZO/trmf5ngaqO7hnDLzrY2vKG1lT3ggkc9eG3fsaekzLCzB9SDLThwSYOiSgccY6SUlLRCSIGYzNSGJsRhIXjPfKGpu9Lqje213P+7vrebe4gY9KG4gij4U09Hhk3f4jsglZiUwdnMwRQwIcMTjA1MEBhg3SNbKOKGmJiHQgKcE4fHCAwwfvb6lY0+hYWeIdib23u55luxtYUx7d9TFHy/1j+zsHBhiWluAlsCFeEsuoNiY0OxLV/H4fJS0RkU5I8/s9nDk0eV9ZZUMzK/Z4iWzxxhLW1aV2ODxLsF01zbywrY4X/HvIII1BK3YwJTeJw3MD+xLnYYMH7ulFJS0RkS6SGUjYN4zK7JSdFBQcREV9MytKGljmXxtbvqeBwiiPyMBrft/SNVWwcZmJHJbrJbDDcwMclhtgXGZivz8qU9ISEelGWcktA1um7Curamjmo9JGVpTUs2JPAytKvGtkdTH0Ur+xsomNlU08s3l/E/y0ROOQnCQm5waY4v9NrTMmOtdvboyOS9Iys38A5wFFzrnD/bLBwMPAOGAj8FnnXKl5W/a3wLlANXCFc+69eMQpIhIP6YGEVqcWG5u9e8g+KGngAz+RLS+upbwx+mRT0+RYtqeBZXuCj8rSyFq+g8k5AQ7NSeLQnABTcr2/Q9MS+lwyi9eR1r+APwD3B5XdCCx0zt1uZjf60zcA5wAF/uNY4M/+XxGRfispwZiSG2BKboDPHeyVrVlTSMaoCawsaWBlaYP3t6SBtRWNNEd7fhGvv8W3i+p5u6g+pDw3peW+taR9fw/JCTCsFyezuCQt59yrZjYurHgOcIr//D5gEV7SmgPc75xzwFtmlmNmI5xzO+IRq4hIb2EGI9MTGZmeyFljUveVVzc2s7qskZUlDXxY2sCHJQ18VNrInrqOe70PVlrXursqgOxkL5lNyk7yHjkBkmt7SUtG51xcHninAVcGTZeFvV7q/30aODGofCEwI7y+srIy1/Jw0OpRdffd+16vuvvuiPO0PILrapw2rc356i6/fN98FYsWtVtnxaJF++atu/zyNudrnDZt33xLlixpt86+tk5Llixp8/Ppi+vUsj79ZZ32HnpoyPJ7cp0+vP9+rVM7+16067Tzl3e6Jz/a6X7x+lb3j6/f2u683PXhvsfS0VPanO+vsy7aN9+s6x/p9nUKKXetc0lvbIgRKY27WCsp2rWL3YWFAOTt2sW4duYt9OcDmFxXR3ob85WXl7PJn3fQ5s1MaafOLZs3U53u1XRQeTn5bcxXW1cXsvz29MV1KiwsZEY7dfaldQqOoT+sU/jye3Sd2qkPBvY6hX9HtLdONaXFjKrawqhkyMvc2+7yUxIcdc2xHTXVd3Ag97WXdlJQtZ5zhzZ12TqFM/9optv5pwefDmqIsRo4xTm3w8xGAIucc4eY2V/95w+GzxdcX3l5eZcEXlhYSEFBQY/X0Zti6W/ro1i6rw7F0vtjaauOZufYvLeJj8saWF3WyKrSBlaXN7KmrJGqaLr6aMOPj87iuiMyDyTkfbKzs1tl1Z480poPXA7c7v+dF1T+dTN7CK8BRrmuZ4mIdK0EM8ZlJjEuM4nZY/aXNzvH1qomVpd5fSiuKfN6+vi4pI7SKDoSnpjVvWklXk3eH8RrdJFnZluBn+Alq0fM7EpgM3CxP/sCvObua/GavH8pHjGKiIiXzFr6Xjxz9P7ywsJChoyZQKHfj+Jav2PgwvIGNlY27Rtsc1JOP0hazrnPt/HS6RHmdcA13RuRiIjEanBqIsemJnLssJSQ8vomx8ZKL4lNyOwHSUtERPqv5ERjUk6ASTmBbl/WwOxxUURE+iQlLRER6TOUtEREpM9Q0hIRkT5DSUtERPqMuPWI0dW6qkcMERHpnSL1iKEjLRER6TOUtEREpM/os6cHRURk4NGRloiI9BkDNmmZ2T/MrMjMVh5gPTlm9piZfWxmq8zsuM4s28wuNrMPzazZzKIaXqatdTCzb5jZar++X3ZQxxgze9mP/UMz+1as8bRVRydiSTWzd8xsuT//T/3yr5vZWjNzZpbXyTrMzG41szV+nN9srx7/PYlm9r6ZPR1rHB3UE1MsZrbRzD4ws2VmttQv68z+0qoevzyWz6jV/t7JWCL+30Qbi5kd4q9Hy6PCzK6NNZa26unEdvm2P99KM3vQ3w9j2l/aqKMz++23/Do+DFqXWLdLqzpi3SbdJl4jF/e2B3AScBRBoyl3sp77gKv858lATmeWDUwGDgEWEWGk5hjqORV4EUjxp4d2UMcI4Cj/eSawBpgSSzzt1BFrLAZk+M8DwNvALOBIvJGvNwJ5nazjS8D9QEI0sfjzXAf8F28cOGKJo4N6Yool0vI6ub9EqifWz6jV/t7JWCLVE1MsQXUlAjuBgzoTSxv1RB0LMArYAKT5048AV8S437ZVR6z7yuHASrwxKJP8dSiIZbu0U0enPp+ufgzYDnOdc6+aNzBlp5lZFl7iuMKvsx6o78yynXOr/DqjXn4b6/D/gNudc3X+PEUd1LED2OE/rzSzVcAo59wL0cbTVh3Al2OMxQEtw60G/Idzzr0fQywR68DbLv/nnGuOJhYzGw18ErgVL+kQSxzt1RNrLJF0Zn9pQ9T7Szv7e1kssbRVj5nFtO8GOR1Y55zbFLSMKN8auR4z+1WMsSQBaWbWgPdlv70T+0urOoCfE9u+Mhl4yzlX7S/7FeAC59wvY4glYh14Awx35vPpUgP29GAXmQAUA//0T//cY2ZtjcQdL5OAT5jZ22b2ipnNjPaNfgI8Eu/opFPC6og5FvNOpS0DioAXnHMxx9JGHQcDnzOzpWb2PzPraDjYu4HvAR0MMN6hSPXEGosDnjezd83s6gOIJVI9sXxGXbW/t1VPZ/fdS4AHOxFHe/VEHYtzbhvwa7xxAXfgDVz7fCwLbqeOWPeVlcBJZjbEzAbhjU04poP3RFtHp79bupKS1oFJwjs992fn3JFAFXBjz4ZEEpCLd0rsu3gDbXb488rMMoDHgWudcxWdWXCEOmKOxTnX5JybDowGjjGzw2ONo406UoBa59wM4O/AP9pZj/OAIufcu7EuO8p6oo7Fd4Jz7ijgHOAaMzupkyFFqieWz6ir9ve26ol5fzGzZODTwKOdiKO9eqKOxcxygTnAeGAkkG5mX4xx+W3VEdO+4h+B3wG8ADwLLAcaY4mlnTo69d3S1ZS0DsxWYGvQ0cBjeP+MPWkrMNd53sH7hd9R44UAXrL5j3NubmcW2kYdMcfSwjlXhnf+fXZn4olQx1Y/PoAngCPaeesJwKfNbCPwEHCamf27EyG0VU8sseCc2+7/LfLnP6YTsbRVTyyfUVft723V05n95RzgPefcrk7E0V49scRyBrDBOVfsnGsA5gLHx7j8tuqIaV8BcM7d65w7yjl3ElACFMYYS1t1dPr/uSspaR0A59xOYIuZHeIXnQ581IMhATwJnAZgZpPwLnLvbmtm/5fSvcAq59ydnVlgO3XEGku+meX4z9Pw/pE/jjGWturYFwtwMl5jkYicc993zo12zo3DO2X0knMupl/OHdQTdSxmlm5mmS3PgbPwTt/EpJ16ov6Mump/b6eemPYX3+fpmlOD4fXEEstmYJaZDfL/F04HVsW4/LbqiHpfaWFmQ/2/Y4EL6cT2aaOOznw+Xc/1QOuP3vDwP4QdQAPeL4grO1nPdGApsML/UHM7s2y8C51bgTpgF/BcJ+tJBv6N94X0HnBaB3WciHetYwWwzH+cG0s87dQRayxHAO/79awEfuyXf9OPpRHv4vQ9nagjB3gG+AB4E5gW5ed7Cvtb/UUdRwf1RB0L3vWf5f7jQ+CHfnlM+0s79cT6GbXa3zu570aqJ9ZYBgF7gOygss7EEqmeWGP5Kd6Po5XAA3in9WLaX9qoI+b9FngN70fAcuD0Tu4vkeqIaZt010M9YoiISJ+h04MiItJnKGmJiEifoaQlIiJ9hpKWiIj0GUpaIiLSZyhpiYhIn6GkJSJdyrwhUM7o6Tikf1LSEhGRPkNJS/ol/9f+9Wa2wszKzexh8wbVu8LMFofN68xsov/8X2b2J79H7b1m9rqZDTezu82s1LxBC4+MYvljzGyumRWb2R4z+4NfnmBmN5nZJvMG8LzfzLL918b5sXzJzLb4y/uqmc3016OspR5//iv8+H7vr+PHZnZ60OsjzWy+mZWYNxjhl4Neu9nMHvGXX2neoH4zwt77uB//BgsafLC995rZA8BY4Cl/+33P3+7/9rdDmZktMbNhsX+qIkpa0r99Fq+z3PF43TtdEcP7bsLrDLQOr/uc9/zpx4B2+2g0s0TgaWAT3iCAo/A6zMWP4Qq8AfUmABnAH8KqOBZv0L3P4Q1t8kO8PhQPAz5rZieHzbvej+0nwFwzG+y/9iBe1z0jgYuA24KTGl6v5g/hdRU0vyUOM0sAnsLrwmcUXj9415rZ2R291zl3KV4/ep9yzmU4bxyny4FsvOEthgBfBWra2YQibVLSkv7sd8657c65Erwv4elRvu8J59y7zrlavJ61a51z9zvnmoCH8cYLa88xeIniu865KudcrXOu5ejuC8Cdzrn1zrm9wPeBS8wseEDWn/nveR5v2I4HnXNFzhtz6bWw5RcBdzvnGpxzDwOrgU+a2Ri8PiFv8OtaBtwDXBr03sXOuQX+ej0ATPPLZwL5zrlbnHP1zrn1eMNiXBLFeyNpwEtWE503bMy7rpPD34gM2JGLZUDYGfS8Gi+RRCN4mIuaCNMZHbx/DLDJORdpHKOReEdgLTbh/R8Gny6LZfnbXGgHopv8ZYwESpxzlWGvzQiaDt8+qX7yPAgYaWZlQa8n4iXMdt/bxjo/gLdNHvJ74P83Xoe9DRHmFWmXjrRkoKnC69UbADMb3g3L2AKMDTt6arEdLym0GIvXC3hnx4Ma5Q9lEVzfdv8xuGU4kqDXtkVR5xa8sZ1ygh6Zzrlzo4wppBdu/yjwp865KXhjRJ0HXBZlXSIhlLRkoFkOHGZm080sFbi5G5bxDt6QMbf741ilmtkJ/msPAt82s/HmjfR8G/BwG0co0RgKfNPMAmZ2MTAZWOCc2wK8AfzCX/4ReEPX/CfK+CvM7AYzSzOzRDM73KIfXn0X3vU6AMzsVDOb6l/rq8A7XdgU/SqK7KekJQOKc24NcAvwIt5orIvbf0enltEEfAqYiNcoYSteowrwhkt/AHgV2ADUAt84gMW9jddoYzdwK3CRc26P/9rn8RqCbMe7NvcT59wLMcQ/3Y9xN971sOwoY/oFcJPfUvB6YDheA5YKvIENX8E7RSgSM42nJdJHmdkVwFXOuRN7OhaReNGRloiI9BlqPSjSCWY2Fm848kimOOc2xzMekYFCpwdFRKTP0OlBERHpM5S0RESkz1DSEhGRPkNJS0RE+gwlLRER6TOUtEREpM/4/wGrm7tHwVcg0gAAAABJRU5ErkJggg==\n",
            "text/plain": [
              "<Figure size 432x288 with 1 Axes>"
            ]
          },
          "metadata": {
            "needs_background": "light"
          },
          "output_type": "display_data"
        }
      ],
      "source": [
        "fig, ax = plt.subplots()\n",
        "ax.plot(loss[:,0], loss[:,1]);\n",
        "ax.axhline(100, linestyle='--', color='r', linewidth=2)\n",
        "ax.xaxis.set_ticks(np.arange(1, 100, 5));\n",
        "ax.set(xlabel='num_components', ylabel='MSE', title='MSE vs number of principal components');"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "0dMUwF0PqdYx"
      },
      "source": [
        "But _numbers dont't tell us everything_! Just what does it mean _qualitatively_ for the loss to decrease from around\n",
        "$450.0$ to less than $100.0$?\n",
        "\n",
        "Let's find out! In the next cell, we draw the the leftmost image is the original dight. Then we show the reconstruction of the image on the right, in descending number of principal components used."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "scrolled": true,
        "id": "KibETu40qdYx",
        "outputId": "f00129db-efdc-43ce-86db-fdf6f16cb8df",
        "colab": {
          "referenced_widgets": [
            "6ba659149ab54409b2b5c07c4137d932"
          ]
        }
      },
      "outputs": [
        {
          "data": {
            "application/vnd.jupyter.widget-view+json": {
              "model_id": "6ba659149ab54409b2b5c07c4137d932",
              "version_major": 2,
              "version_minor": 0
            },
            "text/plain": [
              "interactive(children=(IntSlider(value=500, description='image_idx', max=1000), Output()), _dom_classes=('widge…"
            ]
          },
          "metadata": {},
          "output_type": "display_data"
        }
      ],
      "source": [
        "@interact(image_idx=(0, 1000))\n",
        "def show_num_components_reconst(image_idx):\n",
        "    fig, ax = plt.subplots(figsize=(20., 20.))\n",
        "    actual = X[image_idx]\n",
        "    # concatenate the actual and reconstructed images as large image before plotting it\n",
        "    x = np.concatenate([actual[np.newaxis, :], reconstructions[:, image_idx]])\n",
        "    ax.imshow(np.hstack(x.reshape(-1, 28, 28)[np.arange(10)]),\n",
        "              cmap='gray');\n",
        "    ax.axvline(28, color='orange', linewidth=2)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "eFnr3VN-qdYx"
      },
      "source": [
        "We can also browse throught the reconstructions for other digits. Once again, `interact` becomes handy for visualing the reconstruction."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "scrolled": true,
        "id": "K6gWnyjLqdYx",
        "outputId": "2a0cfe5a-231f-4584-c18c-26cb44366d14",
        "colab": {
          "referenced_widgets": [
            "4996f7863094434cb948efd878147c8a"
          ]
        }
      },
      "outputs": [
        {
          "data": {
            "application/vnd.jupyter.widget-view+json": {
              "model_id": "4996f7863094434cb948efd878147c8a",
              "version_major": 2,
              "version_minor": 0
            },
            "text/plain": [
              "interactive(children=(IntSlider(value=1, description='i', max=10), Output()), _dom_classes=('widget-interact',…"
            ]
          },
          "metadata": {},
          "output_type": "display_data"
        }
      ],
      "source": [
        "@interact(i=(0, 10))\n",
        "def show_pca_digits(i=1):\n",
        "    \"\"\"Show the i th digit and its reconstruction\"\"\"\n",
        "    plt.figure(figsize=(4,4))\n",
        "    actual_sample = X[i].reshape(28,28)\n",
        "    reconst_sample = (reconst[i, :] * std + mu).reshape(28, 28)\n",
        "    plt.imshow(np.hstack([actual_sample, reconst_sample]), cmap='gray')\n",
        "    plt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "7Z-QtjrwqdYx"
      },
      "source": [
        "## 2. PCA for high-dimensional datasets\n",
        "\n",
        "Sometimes, the dimensionality of our dataset may be larger than the number of samples we\n",
        "have. Then it might be inefficient to perform PCA with our implementation above. Instead,\n",
        "we can implement PCA in a more efficient manner, which we call \"PCA for high dimensional data\" (PCA_high_dim).\n",
        "\n",
        "Below are the steps for performing PCA for high dimensional dataset\n",
        "1. Compute the matrix $\\boldsymbol X\\boldsymbol X^T$ (a $N$ by $N$ matrix with $N \\ll D$)\n",
        "2. Compute eigenvalues $\\lambda$s and eigenvectors $V$ for $\\boldsymbol X\\boldsymbol X^T$\n",
        "3. Compute the eigenvectors for the original covariance matrix as $\\boldsymbol X^T\\boldsymbol V$. Choose the eigenvectors associated with the M largest eigenvalues to be the basis of the principal subspace $U$.\n",
        "4. Compute the orthogonal projection of the data onto the subspace spanned by columns of $\\boldsymbol U$."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "DYUUrYkbqdYx"
      },
      "outputs": [],
      "source": [
        "### PCA for high-dimensional datasets\n",
        "\n",
        "def PCA_high_dim(X, n_components):\n",
        "    \"\"\"Compute PCA for small sample size but high-dimensional features.\n",
        "    Args:\n",
        "        X: ndarray of size (N, D), where D is the dimension of the sample,\n",
        "           and N is the number of samples\n",
        "        num_components: the number of principal components to use.\n",
        "    Returns:\n",
        "        X_reconstruct: (N, D) ndarray. the reconstruction\n",
        "        of X from the first `num_components` pricipal components.\n",
        "    \"\"\"\n",
        "    N, D = X.shape\n",
        "    M = np.dot(X, X.T) / N\n",
        "    eig_vals, eig_vecs = eig(M)\n",
        "    eig_vals, eig_vecs = eig_vals[:n_components], eig_vecs[:, :n_components]\n",
        "    U = (X.T @ (eig_vecs))\n",
        "    answer = np.zeros((N, D))\n",
        "    answer = ((U @ np.linalg.inv(U.T @ U) @ U.T) @ X.T).T\n",
        "    return answer"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "NHVpJi76qdYx"
      },
      "source": [
        "Given the same dataset, `PCA_high_dim` and `PCA` should give the same output.\n",
        "Assuming we have implemented `PCA`, correctly, we can then use `PCA` to test the correctness\n",
        "of `PCA_high_dim`. Given the same dataset, `PCA` and `PCA_high_dim` should give identical results.\n",
        "\n",
        "We can use this __invariant__\n",
        "to test our implementation of PCA_high_dim, assuming that we have correctly implemented `PCA`."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "frgHUXHpqdYx"
      },
      "outputs": [],
      "source": [
        "np.testing.assert_almost_equal(PCA(Xbar, 2), PCA_high_dim(Xbar, 2))"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "UVog5hiaqdYx"
      },
      "source": [
        "Now let's compare the running time between `PCA` and `PCA_high_dim`.\n",
        "\n",
        "__Tips__ for running benchmarks or computationally expensive code:\n",
        "\n",
        "When you have some computation that takes up a non-negligible amount of time. Try separating\n",
        "the code that produces output from the code that analyzes the result (e.g. plot the results, compute statistics of the results). In this way, you don't have to recompute when you want to produce more analysis."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "0I1EN5khqdYx"
      },
      "source": [
        "The next cell includes a function that records the time taken for executing a function `f` by repeating it for `repeat` number of times. You do not need to modify the function but you can use it to compare the running time for functions which you are interested in knowing the running time."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "nxOjHt1_qdYx"
      },
      "outputs": [],
      "source": [
        "def time(f, repeat=10):\n",
        "    times = []\n",
        "    for _ in range(repeat):\n",
        "        start = timeit.default_timer()\n",
        "        f()\n",
        "        stop = timeit.default_timer()\n",
        "        times.append(stop-start)\n",
        "    return np.mean(times), np.std(times)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "ZiRcQmZ_qdYy"
      },
      "source": [
        "We first benchmark the time taken to compute $\\boldsymbol X^T\\boldsymbol X$ and $\\boldsymbol X\\boldsymbol X^T$. Jupyter's magic command `%time` is quite handy."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "f3fCr7D1qdYy"
      },
      "source": [
        "The next cell finds the running time for computing $\\boldsymbol X^T\\boldsymbol X$ and $\\boldsymbol X\\boldsymbol X^T$ for different dimensions of $\\boldsymbol X$."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "IRTIWRAmqdYy"
      },
      "outputs": [],
      "source": [
        "times_mm0 = []\n",
        "times_mm1 = []\n",
        "\n",
        "# iterate over datasets of different size\n",
        "for datasetsize in np.arange(4, 784, step=20):\n",
        "    XX = Xbar[:datasetsize] # select the first `datasetsize` samples in the dataset\n",
        "    # record the running time for computing X.T @ X\n",
        "    mu, sigma = time(lambda : XX.T @ XX)\n",
        "    times_mm0.append((datasetsize, mu, sigma))\n",
        "\n",
        "    # record the running time for computing X @ X.T\n",
        "    mu, sigma = time(lambda : XX @ XX.T)\n",
        "    times_mm1.append((datasetsize, mu, sigma))\n",
        "\n",
        "times_mm0 = np.asarray(times_mm0)\n",
        "times_mm1 = np.asarray(times_mm1)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "uNcr3vFiqdYy"
      },
      "source": [
        "Having recorded the running time for computing `X @ X.T` and `X @ X.T`, we can plot them."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "_ETiC2BkqdYy"
      },
      "outputs": [],
      "source": [
        "fig, ax = plt.subplots()\n",
        "ax.set(xlabel='size of dataset', ylabel='running time')\n",
        "bar = ax.errorbar(times_mm0[:, 0], times_mm0[:, 1], times_mm0[:, 2], label=\"$X^T X$ (PCA)\", linewidth=2)\n",
        "ax.errorbar(times_mm1[:, 0], times_mm1[:, 1], times_mm1[:, 2], label=\"$X X^T$ (PCA_high_dim)\", linewidth=2)\n",
        "ax.legend();"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "ExfhkGCEqdYy"
      },
      "source": [
        "Alternatively, use the `time` magic command for benchmarking functions."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "IfisU19mqdYy"
      },
      "outputs": [],
      "source": [
        "%time Xbar.T @ Xbar\n",
        "%time Xbar @ Xbar.T\n",
        "pass # Put this here so that our output does not show result of computing `Xbar @ Xbar.T`"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "422ln_qbqdYy"
      },
      "source": [
        "Next we benchmark PCA, PCA_high_dim."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "0fqIASTtqdY6"
      },
      "outputs": [],
      "source": [
        "times0 = []\n",
        "times1 = []\n",
        "\n",
        "# iterate over datasets of different size\n",
        "for datasetsize in np.arange(4, 784, step=100):\n",
        "    XX = Xbar[:datasetsize]\n",
        "    npc = 2\n",
        "    mu, sigma = time(lambda : PCA(XX, npc), repeat=10)\n",
        "    times0.append((datasetsize, mu, sigma))\n",
        "\n",
        "    mu, sigma = time(lambda : PCA_high_dim(XX, npc), repeat=10)\n",
        "    times1.append((datasetsize, mu, sigma))\n",
        "\n",
        "times0 = np.asarray(times0)\n",
        "times1 = np.asarray(times1)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "0iqm75HzqdY6"
      },
      "source": [
        "Let's plot the running time. Spend some time and think about what this plot means. We mentioned in lectures that PCA_high_dim are advantageous when\n",
        "we have dataset size $N$ < data dimension $M$. Although our plot does not for the two running time does not intersect exactly at $N = M$, it does show the trend."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "scrolled": true,
        "id": "reqG7wgbqdY6"
      },
      "outputs": [],
      "source": [
        "fig, ax = plt.subplots()\n",
        "ax.set(xlabel='number of datapoints', ylabel='run time')\n",
        "ax.errorbar(times0[:, 0], times0[:, 1], times0[:, 2], label=\"PCA\", linewidth=2)\n",
        "ax.errorbar(times1[:, 0], times1[:, 1], times1[:, 2], label=\"PCA_high_dim\", linewidth=2)\n",
        "ax.legend();"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "KyrDVruxqdY6"
      },
      "source": [
        "Again, with the magic command `time`."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "qya81hl5qdY6"
      },
      "outputs": [],
      "source": [
        "%time PCA(Xbar, 2)\n",
        "%time PCA_high_dim(Xbar, 2)\n",
        "pass"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "poCm_3GdqdY7"
      },
      "outputs": [],
      "source": []
    }
  ],
  "metadata": {
    "coursera": {
      "course_slug": "mathematics-machine-learning-pca",
      "graded_item_id": "CXC11",
      "launcher_item_id": "ub5A7"
    },
    "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.7"
    },
    "colab": {
      "provenance": []
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}