notebooks/L28demo_pca2.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## PCA (Part 2)\n",
    "\n",
    "CPSC 340: Machine Learning and Data Mining\n",
    "\n",
    "The University of British Columbia\n",
    "\n",
    "2017 Winter Term 2\n",
    "\n",
    "Mike Gelbart"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pickle\n",
    "import os.path\n",
    "import random\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "from sklearn.decomposition import PCA"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "We have 2414 faces, each 32x32 pixels.\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWQAAAC6CAYAAACZWDfLAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAHsRJREFUeJztnXuwXXV1x7/r3gTyuHmSGAghiSACIfLQUjpWphaoHQpB7Vilgqid6bRlqG3HtrZ2UJrW1s5YqzMdW6aWME7RVqAVEYW08ihQok4mIBAJRkLI+/1OCLn37P6xd+Ds9fuee9c5996c39XvZyaTu9f57b1/e+/fXmef9d1r/awoCgghhOg+Pd3ugBBCiBI5ZCGEyAQ5ZCGEyAQ5ZCGEyAQ5ZCGEyAQ5ZCGEyAQ55GFiZh8xs8dP0L7ea2YbzOygmV18IvYpRh8zu9XM/q3b/cgZM/s9M9tWjf1Tut2f0WLMOmQze8TM9pjZyd3uywnkcwBuLoqiryiKVd3ujIhROZHj/xpmdqRp+fpu9y93zGw8gM8DeFc19nd1u0+jxZh0yGa2EMBlAAoA17axnpnZmDzmigUAnut2J0R7VE6kryiKPgAvA1jSZLuz2/3rBDPrPYG7mwNgAn4Gxv5YdU43AlgB4A4AHx6sYfUk/RkzewLAYQBnmtlHzexHZnbAzF40s99pav9OM9toZh83s+1mtsXMPtr0+Slm9k0z229m3wdwltvf283sB2a2r/r/7a4vf21m/1c9Hd1Xbe/Oans/qL5s/DGcbGYHAfQCeNrMflLZ/8zMflIdx2oze69b77ebjnO1mb21ss81s3vMbIeZrTOzj4XOuhhNTjKzr1TX6jkz+7njH7RzvczsajNbVY2nDWZ2a9NnD5jZza7902b269Xf55rZf5vZbjNbY2bvb2p3h5n9k5l928wOAfjlwfZVrXOjma03s11mdouZvWRmV1af9TSN311m9nUzm0mO580A1lSLe83socr+xWqf+81spZld1rROr5l9suneWGlmZwx1jFlQFMWY+wdgLYCbALwNwDEAcwZp+wjKp5LzAYwDMB7A1SgdqQH4JZSO+q1V+3cC6AewtGr7a9XnM6rP/x3A1wFMBrAYwCYAj1efzQSwB8CHqn39ZrV8SlNf1lb7ngZgNYAXAFxZtf8KgGWDHEsB4E1Ny78BYC7KL9YPADgE4LSmzzYBuKQ6zjehfMLuAbASwKcAnATgTAAvAvjVbl/Xn4V/AF4CcKWz3QrglWqs9QL4WwArqs/aul7V+H1Ltd4FALYBeE/12Y0AnmhquwjAXgAnV+N5A4CPVmPxrQB2Aji/ansHgH0AfrHa9oQh9rUIwEEA76j6/TmU9+qV1ed/iPKhal61/9sAfK3FMS2sxv64JtsNAE6p+vpxAFsBTKg++xMAzwA4pxr7F1ZtBz3GHP51vQMdDOh3VBd2VrX8PIA/GqT9IwCWDrHNbwD4g6YBfcRd/O0AfqG6WY4BOLfps7/B6w75QwC+77b9JICPNPXlL5o++3sA32laXgLgqUH6WXPI5POnALy7+vvB48fk2lwK4GVn+3MM8kWgfyM6fl8Cd8j/07S8CMCRkbheAL4A4B+qv6eg/NJeUC1/BsDt1d8fAPCYW/c2AJ+u/r4DwFfa2Nenmh0sgEkAXsXrDvlHAK5o+vy06t4aR7a7EM4hkzZ7AFxY/b3m+H3g2gx6jDn8G4shiw8DWF4Uxc5q+asYImyB8lvxNczsKjNbUf1s2YvyyWRWU5NdRVH0Ny0fBtAHYDbKb9bm7a1v+nuuWz7++elNy9ua/j5ClvuGOJbm47jRzJ4ys73VcSxuOo4zAPyErLYAwNzj61TrfRJlnE50j61Nfx8GMMHMxqHN62Vml5rZw1V4Yx+A30U1JoqiOADgfgDXVc2vA3A8hr0AwKVuP9cDOLVp8/4+arkvlPfCa+2LojgMoFmMWwDgv5r29SMAA62Oixznx6tw3L5q/WmIjf2hjrGrjOt2B9rBzCYCeD+AXjM7PoBPBjDdzC4siuLpFqu+VtLOyrcy7kH58+3eoiiOmdk3UP60GYodKMMZZ6B8MgeA+U2fb0Z50ZuZD+CBwLbbwswWAPgXAFcAeLIoigEzewqvH8cGuPh2k31dURRnj3SfxKjQ7vX6KoB/BHBVURSvmNkXUH/Y+BqAT5vZ/wKYCODhpv08WhTFrwyybV8acrB9bUEZMgDw2r3b/LraBgC/VRTFE8Hjeo0qXvwJlGP/uaIoGma2B+nYf9atGjnGrjLWnpDfg/JbdBGAi6p/5wF4DKWDjXASSie+A0C/mV0F4F2RFYuiGADwnwBuNbNJZrYI9afzbwN4s5l90MzGmdkHqr5+K9i3dpiM8gbZAQBWCo+Lmz7/MoA/NrO3WcmbKif+fQD7zewTZjaxEkAWm9klo9BHMXzavV5TAOyuHOTPA/ig+/zbKB8algL4j6IoGpX9WyjH7ofMbHz17xIzO2+Qvg22r7sBLLFS5D4JwF+i/tDzzwA+U41JmNlsM3v3kGfj9f32oxz748zsUwCmNn3+ZQB/ZWZnV2P/AivfXe7kGE8oY80hfxhl7Ozloii2Hv+H8lv6+uon3qBUP9s+hlKY24NyEH2zjT7cjDKssBVlXG1Z07Z3AbgGpciwC8CfArimKbwyYhRFsRplDPpJlGGPtwB4ounzu1DGCL8K4ADKOPnM6ktlCcovs3UoRY0vo/zJJzKjg+t1E4ClZnYAZRz36257R1E+VFyJcmwctx9A+WByHcpfelsB/B3Kh5dWtNxXURTPAfh9lCL4FpRjcDuAo1WTL6K875ZX669AGS+P8CCA76AUxNejFESbwymfr/qyHMB+AP8KYGKHx3hCsSqwLYQQo4aZ9aF8o+PsoijWdbs/uTLWnpCFEGMEM1tShfYmo3zt7RmUb5mIFsghCyFGi3ejDA1sBnA2gOsK/SQfFIUshBAiE/SELIQQmSCHLIQQmdBWYsiECROKvr56ItnJJ6dvjPT21gtB9fSkft+3aWVrNBq1ZRZiefXVV4dcr1W7SMimv79/yDZA2v+BgYGkDTsX0bCR7wc7RmYzS3Ne2D69jfU/sl6rfQ5Fo9FAo9Fof8VhMmXKlGLWrFk127hx6a3hj4kdY9TmxwEbF5H12mkXWa/Tazmc0GdkXTau2XqRdsMZ194WuQd37NiBAwcOhMZ1Ww65r68PS5Ysqdne+MY3Ju2mT5+erOeZOnVqYmPtjh49WltmznHDhg2J7dChQ4ntpZdeSmzHjh1LbJ49e/YkNnYh/DGxPkyYMCGxRR3+rl31MrCHDx9O2hw8eDCxnXTSSaF9+i+sffv2JW389QBiXzyRm/rAgQNDthkNZs2ahaVLl9ZsfgwDwMSJE2vLzGmzc83a+XEwadKkpA172GHjh7Xz/WAPO9GHA98u6ggZbKz4dVkbdp+yByx2z/kxy8ZZ9KHuyJEjg24bSO/BW265JWnTCoUshBAiE+SQhRAiE9oKWfT29mLatHrGJgs9+J9748ePT9pEQwr+JwHb1sqVKxMb+1nCfqK9/PLLtWUWBoiENQD+s9DDfuKwY2J99T+roj+P2U/MV155JbH5n4psW9Hwit/nSMWZTxRsfPqfouwasZ/b7Dh9mIGFIlgYw99/ADBzZlLXPWnHtsVg1zwajvCwa87Gjz+vPiwA8PuZXaP9+/cnNn/fsGvE9smumz8/7Hj8tto5f3pCFkKITJBDFkKITJBDFkKITJBDFkKITGhb1JsyZUrNxoLtq1atqi2zQLt/pxYoX6D2bN++vbbMBDAWWGeCAhPnvDAWfVmeEXlpnG2fvQMZEYdmz56dtPHnC+Dvd7Nj8iIJEynZeqyvnogA1q26Kv39/cl4jAhSTBhl15Lht8/GxeTJk0M2n9QCAPPnz68tL1y4MGnDBEL2TrPvW+RdYoDfb+y+92OWjeEtW7YkNvaePBPnIqIaG+tM4PQ2Jpb689XOuNYTshBCZIIcshBCZIIcshBCZEJbMeSiKJJ4zP3335+087FgFv9hRGLBwyluw+KYPg7IYkksbs3iSxGisbZIIaRnnnkmtM+tW7cO3ShINJnDn0d2jVgyRDdoNBqJzsEShLxeEk3yYdc8Ej9nNl9PA+AawebNm2vLLN7t48wAj0f7sc7GABvD7Bzu3JlOL/n888/Xlpm/YGOYbT/iQ5gfYLHzSIEpdr19nD+aWAboCVkIIbJBDlkIITJBDlkIITJBDlkIITKhLWXqwIEDePTRR2u2NWvWJO2Y2OFhgkUkiM7otBoVEJuFI5p44mECYXRGiQjRAuOdzijR6fUAUtGKnUMvdnQrMWRgYCAR9Vh/fRIRE9PYOWOijrcx0Y1VM4tWMfSTKqxbty5pE60W6JMfWBsmQrO+MsHO21iCGEsCYeeMiWxePGbJNUxgZtvy4l9kogclhgghxBhEDlkIITJBDlkIITJBDlkIITKhLVFvYGAgCa5HBDwmBLH1olOcs355WCA9IrJ1Ws0MAK644oraMjuexx57LLStyFQ6TFCInouIuBg9biYO+WNnokmk4tmJoKenJxGumPATyc5kghSbtTySGRipvAbwce3XjVZcZBUEfXZgdOoqVnmNnQsvXrKZ05loyEQ9VinO999XrASAuXPnDrkekAqcrI2/T9n90Qo9IQshRCbIIQshRCbIIQshRCbIIQshRCa0Jeo1Gg0qPgxFVMCLiEhMoGLbimax+QB8VFhatGhRYrvvvvtqy3fffXfShpU8vOOOO4bsFxCbHqjT6aai67Lzw8RFD8t6YmJUNxg/fjxOPfXUmo0JUhs3bqwtMwGPTYu0du3axPbDH/6wtszO67nnnpvYzj///MS2d+/exOaFq2gmHbtXvSjFxEaWjci2xc6rb8fGE5sqid0jDz/88JDbZ0LcpZdeGtqnv74LFixI2px22mmJLYqekIUQIhPkkIUQIhPkkIUQIhPaiiGbWUdxv2i8mLWLxHhZTCu6/UhiCzvm22+/PbH5GCyr1jVz5szEdsoppyQ2VhUrEt9myQudJotEpkEH+Av6kSmculXdzdPb24sZM2bUbKy6mJ9GyE+TBPBryZIOVqxYkfSB9cszb968xMbitz4OG5m+DOD3g99WRIsB4olFPgmHJXew5AoWT3/wwQcTmz9OFkM+77zzEtu2bdsSm9cRmO9hyTVR9IQshBCZIIcshBCZIIcshBCZIIcshBCZ0JaoN23aNFx99dU1GxM2Hnnkkdoyq/DEqmkxEcC/0B5NTIhWbYuIlDfddFNiY8fk+8YEGJ8QAPAX9CP9GsnKbq1ske2zpA8vIuUi4DHGjRuXiHqbNm1K2s2ZM6e2zKqGsaSAqVOnJrZLLrmktswSJljSARNV2fn391c0oYqNKS8aMhGRJZ4woY/13/d1+vTpoW2xfixevDix+XPL7ktmYwKtFxzZtfXJI+28CKEnZCGEyAQ5ZCGEyAQ5ZCGEyAQ5ZCGEyIS2RL3p06fj2muvrdlYts/1119fW77nnnuSNnfddVdiiwgPUdGKEckWY1k8F110UWJjgoKHiS1svU4Fl06FuVZEMvVYX1k7v61o9b1u0Nvbmwg4F198cdLu7LPPri0zIY5NP8REqve97321ZSaKRbMiFy5cmNi8UMbOP+tXZOowJvyxccey62bNmjXkukzAY+eajZ9rrrkmsfljZy8UsCp9LAvvggsuqC339fUlbfx0ZWw7rdATshBCZIIcshBCZIIcshBCZIIcshBCZEJbol5vb28S/GbZLL7NOeeck7Q566yzEtuyZcsSm8+MiZZxjJbf9Hz2s59NbCzThokfkf35jDCAC0grV64ccvtM/IhMpwRwESaSURRdz/dtOGLsaMPGNcvCi4hbrHwla+e3xUSr/fv3JzYmCrPz6PfJhCUmZDEheqhtt+oD2ycrP+uFdJbtyLJZ2fhnQqU/TiY2RkvN+vPD7nF/P0R8xWvbC7cUQggxqsghCyFEJsghCyFEJsghCyFEJrQl6hVFkQT0mQiwa9eu2jITmljZussuuyyx+VKeftutiAbSfZYTEwqYaHXGGWcMue3vfe97iY0JNb6sI8DPK+tbZL2RnM+OiTcRwS66Xjfo6elJhKWIwDmcueW8jWXqMQEsOk+dv+dYhh8TpCKZeuwY2bnwGWut1vWZc0wMZH1lfiUiLrJ7hNnaEeOa8eci8jLBa2072qMQQogRRw5ZCCEyQQ5ZCCEyoa0gydGjR7F27dqajb2w7WMvu3fvptvysNjRVVddVVueO3du0mbVqlWJ7YUXXkhsl19+eWLzL41v3LgxacNieRs2bEhsy5cvry0/++yzSRsWQ/vxj3+c2EYyvjrSVeEi24rE0bytW4kiPT09oTij7180hsySOXw7Fq+cNGkS77CDxVK9VsHut0gSCJAeZ3TqsIjmAaT3Pau4yJI5ojqLH2fR2DlrFxmjvg+KIQshxBhEDlkIITJBDlkIITJBDlkIITKhLVHvyJEjePrpp2s29sL54cOHa8ssSL9gwYLEduGFFyY2X02OTZmyePHixPbkk08mtn379iW2devW1ZZvuOGGpM29996b2Pbs2ZPYXnzxxdry1q1bkzZMEGHnMCKSMLEg+tL+SApouSR4DIeI8BMRZ9j1ZUkmkQSr6PZZUokXElkfoskQft3ouItOARYROFlf2TEx8c+3Y2O/0+vNjns495aekIUQIhPkkIUQIhPkkIUQIhPkkIUQIhPaEvWOHTuWTKkUmdKHiVabNm1KbLNmzUpsPgtp7969SRsmFLAMQpbRd91119WWv/SlLyVtHn300cTWaTA/GvCPCA9RMS2aqee3FxUDoyLJUPvrFkVRJOOMHZMXllgbJkhF7pFophsT8JjNi3/R7DQmwPuqbUw4i/QB4NmHPquQ9cFXhAPiQp8//9EMvEiVRHaNohmKDD0hCyFEJsghCyFEJsghCyFEJsghCyFEJrQl6vX39ydTKDFhwAfzWXB/+vTpiY0JcaeffvqQ29qyZUtiW716dWKbMWNGYnvooYdqy48//njS5uDBg4mNiQe+jCATBZigcOjQoVC7KVOm1JZZWcfI1EOt6DTDiO0zF8EuQlEUtDylJ5KxFhGVWLvINE8AF4zY+ImIlExsZ/v04hm759lYZNuKiMJMuGfTQTFxMSIwD6ccrT//bJz7c9HOvaAnZCGEyAQ5ZCGEyAQ5ZCGEyIS2Y8jbtm2r2dgL4T7Owl709lMntcLHrBnr169PbBdddFFimzdvXmLbvHlzbZklgUQrcflzwdZj8SR2flj8zcc5o1PDdFrtLVoBi00h5BlObHu0KYoiifuxOGmk2l6nSTit+uVhcV9m88ezc+fOpE102iJvY1OtseqHLPmLxbv9tGwsnu8rSAI8ySSSmDOc+8YTqb6nGLIQQoxB5JCFECIT5JCFECIT5JCFECIT2hL1Go1GErD2lZqANNjOXhpnYiATP7yQwoL7ixYtCtnYFFE+WWTZsmVJGwYTqTqt8tSpOBQVi6LtvCDC2jABLzKlT6SaWbeSSVi1NzY+fSLCcBJiIpX1ojaGF/HWrl2btGHjlQlxZ555Zm156tSpSRsm6vnp3gA+frx/mD17dtImkoABcMExUiUxavMiHjseb5OoJ4QQYxA5ZCGEyAQ5ZCGEyAQ5ZCGEyIS2RD0gDVgzkc1nqLEgPcvG6evrG3L/LEB+1llnJTY2HRTL+vMV5m6++eakDZvWiQkKXjxgWTwsk6hTmCgTzUKKZM51OjVTtJ3vQ6fV5oZLURShLEt/zdk57HR6oKiAGp3qyd9fbLwymJjpBcIjR44kbVi2IJt2iWXqeducOXOSNtGpktg9F6n2xmDb8teE9UHV3oQQ4qcAOWQhhMgEOWQhhMgEOWQhhMiEthSmgYGBJADPSkf6TB4mKDChjwmEPiDO1mNB+meffTaxsSmiPKeeempiu/POOxPbbbfdlti++93v1paj2T9MQIpM3xPJbATiWYURIS6Slcf6Npyss9Gm0WgkIhgb1/442TWKlhn1xx4R5lq1Y+PHZ8syoYwJzGxqNS/OsfuZTbHEbOwe99mybGomdl4johsQK7/J7qXIPcj259eTqCeEEGMQOWQhhMgEOWQhhMgEOWQhhMiEtkS9np6epLwdy4hjgXvPjh07EhvLEvIiGxMwVqxYkdgOHjyY2JjI4LP8HnjggaTNrbfemtgioli0PCM7bia4+O1FS22yvkYyBpn4MZzynkP1q5sinxdnmKDmxSbWJoo/djYG2PZZOzbOpkyZUluOzDUH8Ow6v30murH7krXz/QLS+5IdDztutk+WMejPNesXg903fowy4U+ZekII8VOAHLIQQmSCHLIQQmRCWzFkVhWLVVDzVdv8i98An/Jl/vz5ic1PEbV+/fqkDZsqhsWcFi5cmNiWL19eW2bVqBiRuFA0IYPF99jL975dNF4cbeePKVoVK1qBLFeKokjON3vh348ppkmwcRepCMauN6uqFo0r+/g/m9ooen19QlW0CmA0iciPH3buo/H6SLJINKGH9dVvi103TeEkhBA/BcghCyFEJsghCyFEJsghCyFEJgx7PiEWgPeBbhZEZ1XbWIB806ZNteXVq1eH+sBsa9asSWw+cB8Vo5iw4W3suNl6rK8Mf37Yi/EMJq4wEaPTKZQ6XS+Xam9M1GMikj/f7LgjUz8BqWDHhDkm6rFtsXW9LSqKRZJ8olXootN9RaaXYscYve995b6o2NhptTffV4l6QggxBpFDFkKITJBDFkKITJBDFkKITGhL1Gs0GonQwConzZw5s7bMAt9M/Ni5c2di85l5UVEjihcxokLZSApSkcprQHoeI9PVAFyciExjw44xWgEu0qZTMXCkYeOajQN/naLV/Nh18lXJWHXCjRs3Jjafucq2xfbJBMLotfT3ePQeYbB9+vPKts8yDdm22BRs/gUCdt2iIrf3P5FzL1FPCCHGIHLIQgiRCXLIQgiRCXLIQgiRCW2X3/QCUWQKk23btiVtmBi4ffv2xOZFjGj5ShZInzp1amLzAgITV6KZdBGRKlq6kOGPkx1j1BYtQehh5z+S3djpNE8ngoGBgaTEJCut6YUlJuhEM/W8EM3EOnbfbNmyJbExwc6fW3Zto2KyHyvsGCNlNVvhxwbzDUzUmzZtWmJjfYtMEcWuN8Nfc5a16Kewa6cUrZ6QhRAiE+SQhRAiE+SQhRAiE+SQhRAiE9rO1POiFwus+3n2WFCbzV3HbJGAOBPFmGjFxA8v5kQy2AAu9EVEjOjcYwwvRkQz3SKlQoHY3GbRvvrzH81q6wb9/f3YvXt3zTZp0qSkHRObPOyYWCYpG58eJjS94Q1v6Gj7bFvRse6PKVpqs1PRkI07tk92Pdg+I+c6OheiF/HYepEXEVqhJ2QhhMgEOWQhhMgEOWQhhMiEYSeG+JeggTTOEo0ht9pnM8OJm0ZiotFqadG4rGckp04aTgw2sk+2fXYtI8kukeSabsWUmTayd+/epF0khhyNu/vzyLY9b968xMbuNxbH9Ptk22fXjV1fv312jCx2y/bJEjd8X6PjjsXFI2MxMj0XEJt6iyWS+bETTSwD9IQshBDZIIcshBCZIIcshBCZIIcshBCZ0JaoB6TBdVbxKsJITukTrTQVrQrXKZEKW8MRJXOl00pfkRf2TwSs2lvk2jHRhyWURJJiWDUzJlqxioXsPPprwq4R63+n1Q7ZeGXHxM7PxIkTa8vR5KmoAO9hxx2ZZgtIX0bYt29f0sZX5GP7a4WekIUQIhPkkIUQIhPkkIUQIhPkkIUQIhPaFvW8GBEJoo+kmBYVA6MBfy+4RIUm1o9INbbh9N/botM1RYlc22i1t0iFMJ8d1U5G00jSaDRoBpzHX18m1syYMSOxRaYkigpZ0apqkWvJxgoT4vx1YWOTCZBerGvVzo+DaPXGaNW5SKYh2xYTOH1VQL/caltR9IQshBCZIIcshBCZIIcshBCZIIcshBCZ0Lao5wP6kQA8C6J3Oq0QI5qVN3ny5MTmM5+2b9+etOk0o2w4WXmjXYqy0+mmokKiFzbY9YgKhKONmSXXgGVp+Ww+BhN0+vr6EpsfUxHhD+BZf8zmzzcTIKMlM/19w+5nJtax/keOKZr5ydpFMu7Y/cyuLSvB6sVfVsrTn+t27uU87gghhBByyEIIkQtyyEIIkQnDnsIpkugQTYboNG7KYmHsBX2Gn7KbxdAiSQOMaAw2OpVOJCYdjct2moQRjeVFpsaKtDkRFEURmlLJJxiwSl/svPppf4A0rjxt2rSkTbTCWWSssBgvg8W7/fajCUNsWqQIUZ2C9YPdq/4eZ3FfNqWcX4/B+hCZwq4VekIWQohMkEMWQohMkEMWQohMkEMWQohMaEvUM7MRE2KYCMAEEd8u+tI4E1xYu4i4FZmCB4j1lREV8Hy74QiETIzw60b7FbFFqu+NdjJMK8wsEYZZAoMXa9g5jCZN+Hbs2FkyB5vCiVVVY0K3J5p4Erku0enK2DnzYmn0vDKxlIlzPumDrccSgdj590JlVGyMoidkIYTIBDlkIYTIBDlkIYTIBDlkIYTIBGtHSDGzHQDWj153xM84C4qimH2id6pxLUaZ8LhuyyELIYQYPRSyEEKITJBDFkKITJBDFkKITJBDFkKITJBDFkKITJBDFkKITJBDFkKITJBDFkKITJBDFkKITPh/kVpQ3skJzh8AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x114edadd8>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "with open(os.path.join(\"img\",\"faces.pkl\"), \"rb\") as f:\n",
    "    faces = pickle.load(f)\n",
    "\n",
    "# helper functions\n",
    "def display_face(f, cmap='Greys_r', vmin=None, vmax=None):\n",
    "    plt.imshow(np.reshape(f, (32,32)).T, cmap=cmap, vmin=vmin, vmax=vmax)\n",
    "    plt.tick_params(axis='both', which='both', bottom='off', left='off', labelbottom='off', labelleft='off')\n",
    "    plt.grid(False)    \n",
    "def reconstruct(model, X):\n",
    "    Z = model.transform(X)\n",
    "    try:\n",
    "        return model.inverse_transform(Z)\n",
    "    except AttributeError: # no inverse_transform function\n",
    "        W = model.components_\n",
    "        return Z@W\n",
    "    # Note on the above: you can't just do Z@W for PCA because of the centring of the data, \n",
    "    # which needs to be undone for the reconstruction\n",
    "\n",
    "\n",
    "plt.figure()\n",
    "plt.subplot(1,2,1)\n",
    "display_face(random.choice(faces))\n",
    "plt.title(\"A random face\")\n",
    "plt.subplot(1,2,2)\n",
    "display_face(np.mean(faces,axis=0))\n",
    "plt.title(\"The average face\")\n",
    "\n",
    "N,D = faces.shape\n",
    "print(\"We have %d faces, each 32x32 pixels.\" % N)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Testing it out with scikit-learn's PCA"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "PCA(copy=True, iterated_power='auto', n_components=None, random_state=None,\n",
       "  svd_solver='auto', tol=0.0, whiten=False)"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pca = PCA()\n",
    "pca.fit(faces)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "With  1 components, variance explained = 36%\n",
      "With  2 components, variance explained = 71%\n",
      "With  3 components, variance explained = 76%\n",
      "With  4 components, variance explained = 79%\n",
      "With  5 components, variance explained = 80%\n",
      "With  6 components, variance explained = 82%\n",
      "With  7 components, variance explained = 83%\n",
      "With  8 components, variance explained = 84%\n",
      "With  9 components, variance explained = 85%\n",
      "With 10 components, variance explained = 85%\n",
      "With 11 components, variance explained = 86%\n",
      "With 12 components, variance explained = 86%\n",
      "With 13 components, variance explained = 87%\n",
      "With 14 components, variance explained = 87%\n",
      "With 15 components, variance explained = 88%\n",
      "With 16 components, variance explained = 88%\n",
      "With 17 components, variance explained = 88%\n",
      "With 18 components, variance explained = 89%\n",
      "With 19 components, variance explained = 89%\n",
      "With 20 components, variance explained = 89%\n"
     ]
    }
   ],
   "source": [
    "cumulative_var_expl = np.cumsum(pca.explained_variance_ratio_)\n",
    "for i in range(20):\n",
    "    print(\"With %2d components, variance explained = %.0f%%\" % (i+1, cumulative_var_expl[i]*100))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Our own implementation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {},
   "outputs": [],
   "source": [
    "pixel_means = np.mean(faces,axis=0)\n",
    "faces_centered = faces - pixel_means\n",
    "u,s,W = np.linalg.svd(faces_centered, full_matrices=False)\n",
    "Z = u@np.diag(s)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(1024,)"
      ]
     },
     "execution_count": 45,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pixel_means.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(2414, 1024)"
      ]
     },
     "execution_count": 46,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Z.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(1024, 1024)"
      ]
     },
     "execution_count": 47,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "W.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "With  1 components, variance explained = 36%\n",
      "With  2 components, variance explained = 71%\n",
      "With  3 components, variance explained = 76%\n",
      "With  4 components, variance explained = 79%\n",
      "With  5 components, variance explained = 80%\n",
      "With  6 components, variance explained = 82%\n",
      "With  7 components, variance explained = 83%\n",
      "With  8 components, variance explained = 84%\n",
      "With  9 components, variance explained = 85%\n",
      "With 10 components, variance explained = 85%\n",
      "With 11 components, variance explained = 86%\n",
      "With 12 components, variance explained = 86%\n",
      "With 13 components, variance explained = 87%\n",
      "With 14 components, variance explained = 87%\n",
      "With 15 components, variance explained = 88%\n",
      "With 16 components, variance explained = 88%\n",
      "With 17 components, variance explained = 88%\n",
      "With 18 components, variance explained = 89%\n",
      "With 19 components, variance explained = 89%\n"
     ]
    }
   ],
   "source": [
    "for i in range(1,20):\n",
    "    reconstr = Z[:,:i]@W[:i]\n",
    "    var_expl = 1-np.sum((reconstr-faces_centered)**2)/np.sum(faces_centered**2)\n",
    "    print(\"With %2d components, variance explained = %.0f%%\" % (i, var_expl*100))"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.8.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}