{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Before you turn this problem in, make sure everything runs as expected. First, **restart the kernel** (in the menubar, select Kernel$\\rightarrow$Restart) and then **run all cells** (in the menubar, select Cell$\\rightarrow$Run All).\n",
    "\n",
    "Make sure you fill in any place that says `YOUR CODE HERE` or \"YOUR ANSWER HERE\", as well as your name below.\n",
    "\n",
    "Rename this problem sheet as follows:\n",
    "\n",
    "    ps{number of lab}_{your user name}_problem{number of problem sheet in this lab}\n",
    "    \n",
    "for example\n",
    "    \n",
    "    ps2_blja_problem1\n",
    "\n",
    "Submit your homework within one week until next Monday, 9 a.m."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "NAME = \"\"\n",
    "EMAIL = \"\"\n",
    "USERNAME = \"\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Part A - PCA for dimension reduction on the MNIST data set\n",
    "\n",
    "As we have already mentioned, principal component analysis can be used to decrease the computational cost of different learning procedures by reducing the number of variables in our model.\n",
    "Therefore, we now use a slightly larger data set.\n",
    "The **MNIST** data set is widely used for testing.\n",
    "It contains 70,000 grey-valued images of size 28x28, each assigned with a digit from 0 to 9.\n",
    "\n",
    "**Task (1 point)**: Download the `mnist_784.csv` from the class web page.\n",
    "Adapt the following code and read the `csv` file into a `pandas.DataFrame` named `df`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "deletable": false,
    "nbgrader": {
     "cell_type": "code",
     "checksum": "e7e8ee8ed72c8f4277c7aab06b70706d",
     "grade": false,
     "grade_id": "cell-3556285c140a4301",
     "locked": false,
     "schema_version": 3,
     "solution": true,
     "task": false
    }
   },
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import numpy as np\n",
    "# YOUR CODE HERE\n",
    "raise NotImplementedError()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "deletable": false,
    "editable": false,
    "nbgrader": {
     "cell_type": "code",
     "checksum": "fa787d85c172382e5ea6b42ddbc116d7",
     "grade": true,
     "grade_id": "cell-3e6a7fdfdaca7707",
     "locked": true,
     "points": 1,
     "schema_version": 3,
     "solution": false,
     "task": false
    }
   },
   "outputs": [],
   "source": [
    "assert type(df) == pd.core.frame.DataFrame\n",
    "assert abs(df['pixel300'].mean() - 94.15655714285714) < 1e-8\n",
    "assert df.shape == (70000, 785)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following cell splits the data into $60000$ training and $10000$ test data.\n",
    "It should be executable once you read the `csv` file correctly."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ntrain = 60000\n",
    "Xtrain = df.iloc[:ntrain,:-1].values.astype(float)\n",
    "ytrain = df.iloc[:ntrain,-1].values\n",
    "Xtest = df.iloc[ntrain:,:-1].values.astype(float)\n",
    "ytest = df.iloc[ntrain:,-1].values\n",
    "\n",
    "# Number of pixels along on axis\n",
    "npxl = np.sqrt(Xtrain.shape[1]).astype(int)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following cell plots the first numbers in the data set.\n",
    "Images can be plottet using the function `plt.imshow()`.\n",
    "\n",
    "**Task (no points)**: Execute the code from below to plot the first 16 samples"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "plt.rcParams['figure.figsize'] = (15,10)\n",
    "npics = 4\n",
    "fig, ax = plt.subplots(npics,npics)\n",
    "ax=ax.flatten()\n",
    "\n",
    "for i in range(npics**2):\n",
    "    x = Xtrain[i,:]\n",
    "    ax[i].imshow(x.reshape((npxl,npxl)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Task (1 points)**: Train and fit a `StandardScaler` using your **training data**.\n",
    "Store the `StandardScaler` object in the variable `stdScaler`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "deletable": false,
    "nbgrader": {
     "cell_type": "code",
     "checksum": "b7e5e50441172a9e4f22c599d7cb7054",
     "grade": false,
     "grade_id": "cell-e12c4546ea0d2858",
     "locked": false,
     "schema_version": 3,
     "solution": true,
     "task": false
    }
   },
   "outputs": [],
   "source": [
    "# YOUR CODE HERE\n",
    "raise NotImplementedError()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "deletable": false,
    "editable": false,
    "nbgrader": {
     "cell_type": "code",
     "checksum": "76b6f34834ad4f0c299062a96a6b2632",
     "grade": true,
     "grade_id": "cell-d19a9adbc87a776f",
     "locked": true,
     "points": 1,
     "schema_version": 3,
     "solution": false,
     "task": false
    }
   },
   "outputs": [],
   "source": [
    "assert abs(stdScaler.mean_.mean() - 33.318421449829934) < 1e-8\n",
    "assert abs(stdScaler.var_.mean() - 4373.017134018651) < 1e-8"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Task (1 point)**: *Transform* both, your training and test set using `stdScaler`.\n",
    "Store the scaled training set as `X`, and the scaled test set as `Xtest_scaled`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "deletable": false,
    "nbgrader": {
     "cell_type": "code",
     "checksum": "22538487dae540f792806fd1a023e380",
     "grade": false,
     "grade_id": "cell-624f1ad1dfdf888c",
     "locked": false,
     "schema_version": 3,
     "solution": true,
     "task": false
    }
   },
   "outputs": [],
   "source": [
    "# YOUR CODE HERE\n",
    "raise NotImplementedError()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "deletable": false,
    "editable": false,
    "nbgrader": {
     "cell_type": "code",
     "checksum": "655b720fd032fe835beccfc8609aa3c5",
     "grade": true,
     "grade_id": "cell-8f0c40aaf5430dda",
     "locked": true,
     "points": 1,
     "schema_version": 3,
     "solution": false,
     "task": false
    }
   },
   "outputs": [],
   "source": [
    "assert abs(X.mean()) < 1e-10\n",
    "assert abs(Xtest_scaled.mean() - 0.00249569730891525) < 1e-10\n",
    "assert abs(Xtest_scaled.var() - 0.916126707981634) < 1e-8"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Task (no points)**: Now plot the scaled numbers again. Try to explain, why some numbers appear lighter and some darker."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(npics,npics)\n",
    "ax=ax.flatten()\n",
    "\n",
    "for i in range(npics**2):\n",
    "    x = X[i,:]\n",
    "    ax[i].imshow(x.reshape((npxl,npxl)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Task (1 point)**: Perform a linear discriminant analysis on your training data set `Xtrain`.\n",
    "Measure the time your computer needs to perform this task.\n",
    "You can do this easily by using the *magic command* `% time` in front of your `*.fit(Xtrain, ytrain)`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "deletable": false,
    "nbgrader": {
     "cell_type": "code",
     "checksum": "96d75341e92cc32fd38a0f24a778f14c",
     "grade": false,
     "grade_id": "cell-6eda76ed48f21ac4",
     "locked": false,
     "schema_version": 3,
     "solution": true,
     "task": false
    }
   },
   "outputs": [],
   "source": [
    "# YOUR CODE HERE\n",
    "raise NotImplementedError()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "deletable": false,
    "editable": false,
    "nbgrader": {
     "cell_type": "code",
     "checksum": "c9af1304c4708b187ef7b61f2a3fd458",
     "grade": true,
     "grade_id": "cell-3e8eed48e99da50c",
     "locked": true,
     "points": 1,
     "schema_version": 3,
     "solution": false,
     "task": false
    }
   },
   "outputs": [],
   "source": [
    "assert abs(lda.coef_.mean() - 9.726994039367636e-05) < 1e-8\n",
    "assert abs(lda.explained_variance_ratio_.mean() - 0.1111111111111111) < 1e-8"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Task (1 point)**: What is the proportion of correct classications on your (scaled) test set? Store your answer in `full_score`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "deletable": false,
    "nbgrader": {
     "cell_type": "code",
     "checksum": "cd2c0ea7ef7b096b29489dab486c1c95",
     "grade": false,
     "grade_id": "cell-4d102227c3620f13",
     "locked": false,
     "schema_version": 3,
     "solution": true,
     "task": false
    }
   },
   "outputs": [],
   "source": [
    "# YOUR CODE HERE\n",
    "raise NotImplementedError()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "deletable": false,
    "editable": false,
    "nbgrader": {
     "cell_type": "code",
     "checksum": "4af8747517cb019fbe46e3c91c4c3c28",
     "grade": true,
     "grade_id": "cell-06ac7425bffe0c4f",
     "locked": true,
     "points": 1,
     "schema_version": 3,
     "solution": false,
     "task": false
    }
   },
   "outputs": [],
   "source": [
    "assert abs(full_score - 0.873) < 1e-8"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we want to perform a truncated principal component analysis of our scaled data `X`.\n",
    "Depending wheather the optional parameter `n_components` is an integer larger or equal to 1, or a float betwean 0 and 1, the behaviour is different:\n",
    "- setting the option `n_components = 0.9` lets the algorithm choose the number of components, such that these principal components declare 90\\% of the variability in the data\n",
    "- setting the option `n_components` to an integer larger or equal to 1 specifies the number of principal components directly\n",
    "\n",
    "**Task (1 point)**: Perform a principal component analysis on the training data `X` and store the learned model in the variable `pca`.\n",
    "\n",
    "**Caution**: You can find out the number of principal components in the model using the attribute `*.n_components_`.\n",
    "There is also the attribute `*.n_components` which is the value of the option that you specified."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "deletable": false,
    "nbgrader": {
     "cell_type": "code",
     "checksum": "b43a04ea386e3e79473af6ec76a92015",
     "grade": false,
     "grade_id": "cell-3db0124e3d51243c",
     "locked": false,
     "schema_version": 3,
     "solution": true,
     "task": false
    }
   },
   "outputs": [],
   "source": [
    "# YOUR CODE HERE\n",
    "raise NotImplementedError()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "deletable": false,
    "editable": false,
    "nbgrader": {
     "cell_type": "code",
     "checksum": "19370fc7a4d7bb5ead2d566b99d41257",
     "grade": true,
     "grade_id": "cell-d3ad652d69ad10a9",
     "locked": true,
     "points": 1,
     "schema_version": 3,
     "solution": false,
     "task": false
    }
   },
   "outputs": [],
   "source": [
    "assert abs(pca.components_.mean() - 0.0006007088858304203) < 1e-8"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Task (1 point)**:\n",
    "Store the principal components as a `numpy.ndarray` named `pc`, and the number of principal components as `n_pc`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "deletable": false,
    "nbgrader": {
     "cell_type": "code",
     "checksum": "cbeaa7060700759784484a4c2284d8ed",
     "grade": false,
     "grade_id": "cell-0ecc5d7d8de9c518",
     "locked": false,
     "schema_version": 3,
     "solution": true,
     "task": false
    }
   },
   "outputs": [],
   "source": [
    "# YOUR CODE HERE\n",
    "raise NotImplementedError()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "deletable": false,
    "editable": false,
    "nbgrader": {
     "cell_type": "code",
     "checksum": "cf6598a2f7235b0c8ce5c485de493fbc",
     "grade": true,
     "grade_id": "cell-6f3cee455fe6169b",
     "locked": true,
     "points": 1,
     "schema_version": 3,
     "solution": false,
     "task": false
    }
   },
   "outputs": [],
   "source": [
    "assert n_pc == 236\n",
    "assert abs(pc[:,100].max() - 63.234284245188896) < 1e-8"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Task (1 point)**: Perform now an LDA using these principal components. Measure the time that is necessary for this operation, and compare the score on your test set with the score obtained by using the full model.\n",
    "Store the proportion of correct classifications on the (scaled) test set in the variable `red_score`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "deletable": false,
    "nbgrader": {
     "cell_type": "code",
     "checksum": "0b01f41a1f1e99dc1b4dd041a3db8f98",
     "grade": false,
     "grade_id": "cell-94cdc08542d78aae",
     "locked": false,
     "schema_version": 3,
     "solution": true,
     "task": false
    }
   },
   "outputs": [],
   "source": [
    "# YOUR CODE HERE\n",
    "raise NotImplementedError()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "deletable": false,
    "editable": false,
    "nbgrader": {
     "cell_type": "code",
     "checksum": "52fb145b982548f91ee9140b9df2cb44",
     "grade": true,
     "grade_id": "cell-db2543fa51a2d935",
     "locked": true,
     "points": 1,
     "schema_version": 3,
     "solution": false,
     "task": false
    }
   },
   "outputs": [],
   "source": [
    "assert abs(red_score - 0.8717) < 1e-8"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You should observe that the computing time has been quartered approximately.\n",
    "The score of 87.2\\% is comparable to the previous reached 87.3\\% in the full model."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### A closer look at PCA, and its connection to SVD\n",
    "As you all know from the lecture, principal component anaylsis (PCA) is strongly connected to a (truncated) singular value decomposition (SVD).\n",
    "\n",
    "Assume, that the matrix $X \\in \\mathbb{R}^{n \\times p}$ is scaled, i.e. has column-mean zero and variance one.\n",
    "Then the covariance matrix $C$ is given by $C = X^T X$.\n",
    "This is a symmetric matrix, and thus can be diagonalized into\n",
    "\n",
    "$$ C = V \\Lambda V^T $$\n",
    "\n",
    "with $V \\in \\mathbb{R}^{p \\times p}$ the matrix of eigenvectors of $C$ and $\\Lambda = \\text{diag}(\\lambda_1, \\ldots, \\lambda_p)$ the matrix of eigenvalues on the diagonal.\n",
    "The columns of the matrix $V$ are called principal directions of the data, and projections of the data on the principal directions are called *principal components*.\n",
    "Thus, the $j$-th column of the matrix $XV$ is called the $j$-th principal component.\n",
    "\n",
    "#### Full SVD\n",
    "\n",
    "If we perform a full SVD of the matrix $X \\in \\mathbb{R}^{n \\times p}$, we obtain the decomposition\n",
    "\n",
    "$$ X = U \\Sigma V^T $$\n",
    "\n",
    "with a unitary matrix $U \\in \\mathbb{R}^{n \\times n}$ (left-singular vectors),\n",
    "a \"diagonal\" matrix $\\Sigma \\in \\mathbb{R}^{n \\times p}$ (singular values $\\Sigma_{i,i}$, rest zero),\n",
    "and a unitary matrix $V \\in \\mathbb{R}^{p \\times p}$ (right-singular vectors).\n",
    "\n",
    "Thus, we see that\n",
    "\n",
    "$$ C = X^T X = (U \\Sigma V^T)^T (U \\Sigma V^T) = V \\Sigma^2 V^T. $$\n",
    "\n",
    "In other words, this says that the right-singular vectors are principal directions, and the principal components are given by\n",
    "\n",
    "$$ XV = (U \\Sigma V^T) V = U \\Sigma. $$\n",
    "\n",
    "#### Truncated SVD\n",
    "\n",
    "In comparison to a full SVD, a truncated SVD **approximates** the matrix by restricting only on the $k$ largest singular values, i.e.,\n",
    "\n",
    "$$ X \\approx U_k \\Sigma_k V_k^T $$\n",
    "\n",
    "with matrices \n",
    "$U_k \\in \\mathbb{R}^{n \\times k}$,\n",
    "$\\Sigma_k \\in \\mathbb{R}^{k \\times k}$,\n",
    "$V_k \\in \\mathbb{R}^{p \\times k}$.\n",
    "\n",
    "Thus, we get the first $k$ principal components by\n",
    "\n",
    "$$ X V_k = U_k \\Sigma_k \\in \\mathbb{R}^{n \\times k}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Task (no points)**: Please read the documentation of the sklearn function `PCA` [here](https://scikit-learn.org/stable/modules/generated/sklearn.decomposition.PCA.html)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You should learn, that `pca.components_` is equivalent to the matrix $V_k^T$.\n",
    "We assign the matrix `V` by `pca.components_.T` and check the size of the matrix.\n",
    "\n",
    "The rest of this notebook is prepared for you.\n",
    "Please make sure, that you understand each step.\n",
    "Ask questions or have a look into the documentation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "V = pca.components_.T\n",
    "print(V.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Can we resemble the principal components? Yes, we can! Remember, it's the matrix product of $X$ and $V_k$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.abs(X.dot(V)-pc).max()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, we can assign the matrix with the singular values.\n",
    "The method `pca.singular_values_` returns a vector, but we can easily store it as a diagonal matrix $\\Sigma$.\n",
    "This is also true for $\\Sigma^{-1}$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Sigma = np.diag(pca.singular_values_)\n",
    "SigmaInv = np.diag(1./pca.singular_values_)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Is the size correct?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Sigma.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And therefore, we can also compute the matrix $U_k$ by simple matrix multiplication.\n",
    "This can, as you know, done either by the `numpy.ndarray` method `.dot()`, or by the operator `@`.\n",
    "\n",
    "**Remember**: The operator `*` performs an element-wise multiplication."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "U = X @ V @ SigmaInv"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Check the size."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "U.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It is an interesting exercise to plot the first principal directions, i.e., the first columns of the matrix $V_k$.\n",
    "What do you observe?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "\n",
    "npics = 4\n",
    "fig, ax = plt.subplots(npics,npics)\n",
    "ax = ax.flatten()\n",
    "for i in range(npics**2):\n",
    "    x = V[:,i]\n",
    "    ax[i].imshow(x.reshape((npxl,npxl)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, we want to compare some numbers and their corresponding approximations using PCA.\n",
    "\n",
    "We can compute the approximations simply by setting \n",
    "\n",
    "$$ X_{\\text{approx}} = U_k \\Sigma_k V_k^T $$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Xapprox = U @ Sigma @ V.T"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, plotting is easy. You can play with $m$ to display different samples in the data set."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "m = 0\n",
    "fig, ax = plt.subplots(1,2)\n",
    "ax[0].imshow(X[m,:].reshape((npxl,npxl)))\n",
    "ax[0].set_title('Original (%d pixel)' %  (npxl**2))\n",
    "ax[1].imshow(Xapprox[m,:].reshape((npxl,npxl)))\n",
    "ax[1].set_title('Approximation')\n",
    "ax[1].set_xlabel('%d principal components')\n"
   ]
  }
 ],
 "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.6.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}