Lab 8 - Unsupervised learning.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Lab 8: Unsupervised Learning from Unlabelled Data\n",
    "\n",
    "## A: [*Principal component analysis* ](#partA); B: [*Clustering*](#partB); C: [*Autoencoder*](#partC)\n",
    "\n",
    "[**Haiping Lu**](http://staffwww.dcs.shef.ac.uk/people/H.Lu/) -  [COM4509/6509 MLAI2020](https://github.com/maalvarezl/MLAI) @ The University of Sheffield\n",
    "\n",
    "**Sources**: Part A is based on the [Dimensionality Reduction: Latent Variable Modelling notebook](https://nbviewer.jupyter.org/github/lawrennd/mlai2015/blob/master/week8.ipynb) by Neil Lawrence [MLAI2015]. Part B covers clustering briefly and uses Scikit-learn examples. Part C is based on the [Autoencoder notebook by Lisa Zhang from APS360, University of Toronto](https://www.cs.toronto.edu/~lczhang/aps360_20191/lec/w05/autoencoder.html) and Lab 2 of my [SimplyDeep](https://github.com/haipinglu/SimplyDeep/)  notebooks.\n",
    "\n",
    "There are *five* questions in this notebook."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Why unsupervised learning?\n",
    "So far we have focused mainly on [supervised learning](https://en.wikipedia.org/wiki/Supervised_learning) problems including regression and classification, where training data samples are all labelled. Now we are going to turn to a different form of learning, commonly known as [unsupervised learning](https://en.wikipedia.org/wiki/Unsupervised_learning). In unsupervised learning, our data are not labelled, but we want models that give us a better understanding of the data, which is called [exploratory data analysis](https://en.wikipedia.org/wiki/Exploratory_data_analysis) in some contexts. Unsupervised learning is a fundamental problem in machine learning and the focus is to learn from unlabelled data. Unsupervised feature learning or representation learning can find wide usage in various applications for extracting useful information from often abundant unlabelled data."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# <a id='partA'></a>A: Principal Component Analysis (PCA)\n",
    "\n",
    "#### PCA example in Scikit-learn\n",
    "Run the Scikit-learn example [Face recognition example using eigenfaces and SVMs](https://scikit-learn.org/stable/auto_examples/applications/plot_face_recognition.html). More examples on PCA are at (the bottom of) the [PCA documentation](https://scikit-learn.org/stable/modules/generated/sklearn.decomposition.PCA.html).\n",
    "\n",
    "The following discusses the background and history of PCA that is complementary to the lecture materials. The latent factor perspective will be useful to understanding generative models in the next session. If you are interetsted in the applications (how to use) PCA only, you may *safely skip* the remaining part of this section (A1 and A2), including the two questions.\n",
    "## A1. Latent variable and latent factor analysis\n",
    "\n",
    "This section discusses the context where PCA was first proposed.\n",
    "### Latent variables\n",
    "\n",
    "Latent means hidden. [Latent variables](https://en.wikipedia.org/wiki/Latent_variable) are hidden, unobservable variables. Recall that [Artificial neural networks](https://en.wikipedia.org/wiki/Artificial_neural_network) have three types of nodes (a.k.a. units) called artificial neurons: input, **hidden**, and output neurons.\n",
    "\n",
    "Latent variables are also called **latent factors** (see below) when the focus is more on **analysing and interpreting** the data. In the following, let us look at the [factor analysis](https://en.wikipedia.org/wiki/Factor_analysis) problem first."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Factor analysis model as a multi-output regression problem\n",
    "\n",
    "Note that in the following, $\\mathbf{y}$ is observed and known while $\\mathbf{x}$ is hidden, unknown, and to be estimated. This can be confusing because in the lecture, the input is $\\mathbf{x}$ and the output (projections) is $\\mathbf{y}$ so please mind the difference.\n",
    "\n",
    "If we are given a high dimensional vector of features (perhaps questionaire answers) associated with an individual, $\\mathbf{y}$, we assume that these factors are actually **generated** from a low dimensional vector latent traits, or latent variables $\\mathbf{x}$, which determine the personality.\n",
    "$$\n",
    "\\mathbf{y} = \\mathbf{f}(\\mathbf{x}) + \\boldsymbol{\\epsilon}\n",
    "$$\n",
    "where $\\mathbf{f}(\\mathbf{x})$ is a *vector valued* function that is dependent on the latent traits and $\\boldsymbol{\\epsilon}$ is some corrupting noise. For simplicity, we assume that the function is given by a *linear* relationship,\n",
    "$$\n",
    "\\mathbf{f}(\\mathbf{x}) = \\mathbf{W}\\mathbf{x}\n",
    "$$\n",
    "where we have introduced a matrix $\\mathbf{W}$ that is sometimes referred to as the *factor loadings* but we also immediately see is related to our *multivariate linear regression* models. That is because our vector valued function is of the form\n",
    "$$\n",
    "\\mathbf{f}(\\mathbf{x}) = \\begin{bmatrix} f_1(\\mathbf{x}) \\\\ f_2(\\mathbf{x}) \\\\ \\vdots \\\\ f_p(\\mathbf{x})\\end{bmatrix}\n",
    "$$\n",
    "where there are $p$ **features** associated with the individual. If we consider any of these functions individually we have a prediction function that looks like a regression model, \n",
    "$$\n",
    "f_j(\\mathbf{x}) = \\mathbf{w}_{i, :} \\mathbf{x}, \n",
    "$$\n",
    "for each element of the vector valued function, where $\\mathbf{w}_{i, :}$ is the $i$th row of the matrix $\\mathbf{W}$. In the context of regression, each column of $\\mathbf{W}$ is a vector of *regression weights*. This is a **multiple input and multiple output regression**. Our inputs (or covariates) have dimensionality greater than 1 and our outputs (or response variables) also have dimensionality greater than one, different from predicting a single output value in our previous sessions. Just as in a standard regression, we are assuming that we don't observe the function directly (note that this *also* makes the function a *type* of latent variable), but we observe some corrupted variant of the function, where the corruption is given by $\\boldsymbol{\\epsilon}$. Just as in linear regression we can assume that this corruption is given by Gaussian noise, where the noise for the $j$th element of $\\mathbf{y}$ is by,\n",
    "$$\n",
    "\\epsilon_j \\sim \\mathcal{N}(0, \\sigma^2_j).\n",
    "$$\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Of course, just as in a regression problem we also need to make an assumption across the individual data points to form our full likelihood. Our data set now consists of many observations of $\\mathbf{y}$ for different individuals. We store these observations in a *design matrix*, $\\mathbf{Y}$, where each *row* of $\\mathbf{Y}$ contains the observation for one individual. To emphasize that $\\mathbf{y}$ is a vector derived from a row of $\\mathbf{Y}$, we represent the observation of the features associated with the $i$th individual by $\\mathbf{y}_{i, :}$, and place each individual in our data matrix,\n",
    "$$\n",
    "\\mathbf{Y} = \\begin{bmatrix} \\mathbf{y}_{1, :}^\\top \\\\ \\mathbf{y}_{2, :}^\\top \\\\ \\vdots \\\\ \\mathbf{y}_{n, :}^\\top\\end{bmatrix},\n",
    "$$\n",
    "where we have $n$ data points. Our data matrix therefore has $n$ rows and $p$ columns. The point to notice here is that each data obsesrvation appears as a row vector in the design matrix (thus the transpose operation inside the brackets). Our prediction functions are now actually a *matrix value* function, \n",
    "$$\n",
    "\\mathbf{F} = \\mathbf{X}\\mathbf{W}^\\top,\n",
    "$$\n",
    "where for each matrix the data points are in the rows and the data features are in the columns. This implies that if we have $q$ inputs to the function we have $\\mathbf{F}\\in \\Re^{n\\times p}$, $\\mathbf{W} \\in \\Re^{p \\times q}$ and $\\mathbf{X} \\in \\Re^{n\\times q}$. Note that $\\mathbf{Y}$ equals to $\\mathbf{F}$ + noise. \n",
    "\n",
    "$n$: number of data point\n",
    "\n",
    "$p$: input feature dimension\n",
    "\n",
    "$q$: latent feature dimension\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "### <a id='Q1'></a>Question 1\n",
    "\n",
    "Show that, given all the definitions above, if,\n",
    "$$\n",
    "\\mathbf{F} = \\mathbf{X}\\mathbf{W}^\\top\n",
    "$$\n",
    "and the elements of the vector valued function $\\mathbf{F}$ are given by \n",
    "$$\n",
    "f_{i, j} = f_j(\\mathbf{x}_{i, :}),\n",
    "$$\n",
    "where $\\mathbf{x}_{i, :}$ is the $i$th row of the latent variables, $\\mathbf{X}$, then show that\n",
    "$$f_j(\\mathbf{x}_{i, :}) = \\mathbf{x}_{i, :} \\mathbf{w}_{j, :} ^\\top,$$\n",
    "where $\\mathbf{w}_{j, :}$ is the $j$th row of $\\mathbf{W}$.\n",
    "\n",
    "**Answer**: \n",
    "\n",
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Gaussian modelling of the factor analysis model\n",
    "\n",
    "The difference between this factor analysis model and a multiple output regression is that in the regression case we are provided with the covariates $\\mathbf{X}$, here they are *latent variables*. These variables are unknown.\n",
    "\n",
    "Just as we have done in the past for unknowns, we now treat them with a probability distribution. In *factor analysis* we assume that the latent variables have a Gaussian density which is independent across both across the latent variables associated with the different data points, and across those associated with different data features, so we have,\n",
    "$$\n",
    "x_{i,j} \\sim \\mathcal{N}(0, 1),\n",
    "$$\n",
    "and we can write the density governing the latent variable associated with a single point as,\n",
    "$$\n",
    "\\mathbf{x}_{i, :} \\sim \\mathcal{N}(\\mathbf{0}, \\mathbf{I}).\n",
    "$$\n",
    "If we consider the values of the function for the $i$th data point as\n",
    "$$\n",
    "\\mathbf{f}_{i, :} = \\mathbf{f}(\\mathbf{x}_{i, :}) = \\mathbf{W}\\mathbf{x}_{i, :} \n",
    "$$\n",
    "then we can use the rules for multivariate Gaussian relationships to write that\n",
    "$$\n",
    "\\mathbf{f}_{i, :} \\sim \\mathcal{N}(\\mathbf{0}, \\mathbf{W}\\mathbf{W}^\\top)\n",
    "$$\n",
    "which implies that the distribution for $\\mathbf{y}_{i, :}$ is given by\n",
    "$$\n",
    "\\mathbf{y}_{i, :} = \\sim \\mathcal{N}(\\mathbf{0}, \\mathbf{W}\\mathbf{W}^\\top + \\boldsymbol{\\Sigma})\n",
    "$$\n",
    "where $\\boldsymbol{\\Sigma}$ the covariance of the noise variable, $\\epsilon_{i, :}$ which for factor analysis is a diagonal matrix (because we have assumed that the noise was **independent** across the features),\n",
    "$$\n",
    "\\boldsymbol{\\Sigma} = \\begin{bmatrix}\\sigma^2_{1} & 0 & 0 & 0\\\\\n",
    "                                     0 & \\sigma^2_{2} & 0 & 0\\\\\n",
    "                                     0 & 0 & \\ddots & 0\\\\\n",
    "                                     0 & 0 & 0 & \\sigma^2_p\\end{bmatrix}.\n",
    "$$\n",
    "For completeness, we could also add in a *mean* for the data vector $\\boldsymbol{\\mu}$, \n",
    "$$\n",
    "\\mathbf{y}_{i, :} = \\mathbf{W} \\mathbf{x}_{i, :} + \\boldsymbol{\\mu} + \\boldsymbol{\\epsilon}_{i, :}\n",
    "$$\n",
    "which would give our marginal distribution for $\\mathbf{y}_{i, :}$ a mean $\\boldsymbol{\\mu}$. However, the maximum likelihood solution for $\\boldsymbol{\\mu}$ turns out to equal the empirical mean of the data,\n",
    "$$\n",
    "\\hat{\\boldsymbol{\\mu}} = \\frac{1}{n} \\sum_{i=1}^n \\mathbf{y}_{i, :},\n",
    "$$\n",
    "*regardless* of the form of the covariance, $\\mathbf{C} = \\mathbf{W}\\mathbf{W}^\\top + \\boldsymbol{\\Sigma}$. As a result it is very common to simply preprocess the data and ensure it is zero mean. \n",
    "\n",
    "The prior density over latent variables is independent, and the likelihood is independent, that means that the marginal likelihood here is also independent over the data points. \n",
    "                                     \n",
    "Factor analysis was developed mainly in psychology and the social sciences for understanding personality and intelligence. [Charles Spearman](http://en.wikipedia.org/wiki/Charles_Spearman) was concerned with the measurements of \"the abilities of man\" and is credited with the earliest version of factor analysis.                                    \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## A2. PCA background\n",
    "This section discusses how PCA was first proposed and the mathematical background for PCA.\n",
    "### Background on PCA and factor analysis (history)\n",
    "\n",
    "In 1933 [Harold Hotelling](http://en.wikipedia.org/wiki/Harold_Hotelling) published on *principal component analysis* the first mention of this approach. Hotelling's inspiration was to provide mathematical foundation for factor analysis methods that were by then widely used within psychology and the social sciences. His model was a factor analysis model, but he considered the noiseless 'limit' of the model. In other words he took $\\sigma^2_i \\rightarrow 0$ so that he had\n",
    "$$\n",
    "\\mathbf{y}_{i, :} \\sim \\lim_{\\sigma^2 \\rightarrow 0} \\mathcal{N}(\\mathbf{0}, \\mathbf{W}\\mathbf{W}^\\top + \\sigma^2 \\mathbf{I}).\n",
    "$$\n",
    "The paper had two unfortunate effects. \n",
    "\n",
    "  * Firstly, the resulting model is no longer valid probablistically, because the covariance of this Gaussian is 'degenerate'. Because $\\mathbf{W}\\mathbf{W}^\\top$ has rank of at most $q$ where $q<p$ (due to the dimensionality reduction) the determinant of the covariance is zero, meaning the inverse doesn't exist so the density,\n",
    "$$\n",
    "p(\\mathbf{y}_{i, :}|\\mathbf{W}) = \\lim_{\\sigma^2 \\rightarrow 0} \\frac{1}{(2\\pi)^\\frac{p}{2} |\\mathbf{W}\\mathbf{W}^\\top + \\sigma^2 \\mathbf{I}|} \\exp\\left(-\\frac{1}{2}\\mathbf{y}_{i, :}\\left[\\mathbf{W}\\mathbf{W}^\\top+ \\sigma^2 \\mathbf{I}\\right]^{-1}\\mathbf{y}_{i, :}\\right),\n",
    "$$\n",
    "is *not* valid for $q<p$ (where $\\mathbf{W}\\in \\Re^{p\\times q}$). This mathematical consequence is a probability density which has no 'support' in large regions of the space for $\\mathbf{y}_{i, :}$. There are regions for which the probability of $\\mathbf{y}_{i, :}$ is zero. These are any regions that lie off the hyperplane defined by mapping from $\\mathbf{x}$ to $\\mathbf{y}$ with the matrix $\\mathbf{W}$. In factor analysis, the noise corruption, $\\boldsymbol{\\epsilon}$, allows for points to be found away from the hyperplane. In Hotelling's PCA, the noise variance is zero, so there is only support for points that fall precisely on the hyperplane. \n",
    "  * Secondly, Hotelling explicity chose to rename factor analysis as principal component analysis, arguing that the factors social scientist sought were different in nature to the concept of a mathematical factor. This was unfortunate because the factor loadings, $\\mathbf{W}$ can also be seen as factors in the mathematical sense because the model Hotelling defined is a Gaussian model with covariance given by $\\mathbf{C} = \\mathbf{W}\\mathbf{W}^\\top$ so $\\mathbf{W}$ is a *factor* of the covariance in the mathematical sense, as well as a factor loading. \n",
    "\n",
    "However, the paper had one great advantage over standard approaches to factor analysis. Despite the fact that the model was a special case that is subsumed by the more general approach of factor analysis, it is this special case that leads to a particular algorithm, namely that the factor loadings (or principal components as Hotelling referred to them) are given by an [eigendecomposition](https://en.wikipedia.org/wiki/Eigendecomposition_of_a_matrix) of the empirical covariance matrix. \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Eigendecomposition\n",
    "\n",
    "Eigenvalue problems are widespreads in physics and mathematics, they are often written as a matrix/vector equation but we prefer to write them as a full matrix equation. In an eigenvalue problem, we are looking to find a matrix of eigenvectors, $\\mathbf{U}$ and a *diagonal* matrix of eigenvalues, $\\boldsymbol{\\Lambda}$ that satisfy the *matrix* equation\n",
    "$$\n",
    "\\mathbf{A}\\mathbf{U} = \\mathbf{U}\\boldsymbol{\\Lambda},\n",
    "$$\n",
    "where $\\mathbf{A}$ is our matrix of interest. This equation is not trivially solvable through matrix inverse because matrix multiplication is not [commutative](http://en.wikipedia.org/wiki/Commutative_property), so premultiplying by $\\mathbf{U}^{-1}$ gives\n",
    "$$\n",
    "\\mathbf{U}^{-1}\\mathbf{A}\\mathbf{U} = \\boldsymbol{\\Lambda}, \n",
    "$$\n",
    "where $\\boldsymbol{\\Lambda}$ is a *diagonal* matrix, so the eigenvectors can be used to *diagonalise* the matrix. When performing the eigendecomposition on a Gaussian covariances, diagonalisation is very important because it returns the covariance to a form where there is no correlation between points. \n",
    "\n",
    "#### Positive definite\n",
    "\n",
    "We are interested in the case where $\\mathbf{A}$ is a covariance matrix, which implies it is *positive definite*. A positive definite matrix is one for which the inner product,\n",
    "$$\n",
    "\\mathbf{w}^\\top \\mathbf{C}\\mathbf{w}\n",
    "$$\n",
    "is positive for *all* values of the vector $\\mathbf{w}$ other than the zero vector. One way of creating a positive definite matrix is to assume that the symmetric and positive definite matrix $\\mathbf{C}\\in \\Re^{p\\times p}$ is factorised into, $\\mathbf{A}\\in \\Re^{p\\times p}$, a *full rank* matrix, so that\n",
    "$$\n",
    "\\mathbf{C} = \\mathbf{A}^\\top \\mathbf{A}.\n",
    "$$\n",
    "This ensures that $\\mathbf{C}$ must be positive definite because\n",
    "$$\n",
    "\\mathbf{w}^\\top \\mathbf{C}\\mathbf{w}=\\mathbf{w}^\\top \\mathbf{A}^\\top\\mathbf{A}\\mathbf{w} \n",
    "$$\n",
    "and if we now define a new *vector* $\\mathbf{b}$ as\n",
    "$$\n",
    "\\mathbf{b} = \\mathbf{A}\\mathbf{w}\n",
    "$$\n",
    "we can now rewrite as\n",
    "$$\n",
    "\\mathbf{w}^\\top \\mathbf{C}\\mathbf{w} = \\mathbf{b}^\\top\\mathbf{b} = \\sum_{i} b_i^2\n",
    "$$\n",
    "which, since it is a sum of squares, is positive or zero. The constraint that $\\mathbf{A}$ must be *full rank* ensures that there is no vector $\\mathbf{w}$, other than the zero vector, which causes the vector $\\mathbf{b}$ to be all zeros. \n",
    "\n",
    "You may want to review what is the rank of a matrix and what full rank means [here](https://en.wikipedia.org/wiki/Rank_(linear_algebra))."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "### <a id='Q2'></a>Question 2\n",
    "\n",
    "If $\\mathbf{C}=\\mathbf{A}^\\top \\mathbf{A}$ then express $c_{i,j}$, the value of the element at the $i$th row and the $j$th column of $\\mathbf{C}$, in terms of the columns of $\\mathbf{A}$. Use this to show that (i) the matrix is symmetric and (ii) the matrix has positive elements along its diagonal, assuming $\\mathbf{A}$ is of full rank.\n",
    "\n",
    "**Answer**: \n",
    "\n",
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Eigenvectors of a symmetric matrix\n",
    "\n",
    "[Symmetric matrices](https://en.wikipedia.org/wiki/Symmetric_matrix) have *orthogonal* eigenvectors, which are often normalised to unit norm to obtain *orthonormal* eigenvectors. This means that $\\mathbf{U}$ is an [orthogonal matrix](http://en.wikipedia.org/wiki/Orthogonal_matrix), $\\mathbf{U}^\\top\\mathbf{U} = \\mathbf{I}$. This implies that $\\mathbf{u}_{:, i} ^\\top \\mathbf{u}_{:, j}$ is equal to 0 if $i\\neq j$ and 1 if $i=j$. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### PCA for high-dimensional data  \n",
    "For $p$>>$n$ (generally the case in practice), you need to consider how to do the larger eigenvalue probleme efficiently without large demands on computer memory.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "When the data is quite high dimensional, solving the eigenvalue problem in the high dimensional space can take some time. At this point we turn to a neat trick, you don't have to solve the full eigenvalue problem in the $p\\times p$ covariance, you can choose instead to solve the related eigenvalue problem in the $n \\times n$ space, and in this case $n=200$ which is much smaller than $p$.\n",
    "\n",
    "The original eigenvalue problem has the form\n",
    "$$\n",
    "\\mathbf{Y}^\\top\\mathbf{Y} \\mathbf{u} = \\lambda\\mathbf{u}\n",
    "$$\n",
    "Let us look at a related eigenvalue problem.\n",
    "$$\n",
    "\\mathbf{Y}\\mathbf{Y}^\\top \\mathbf{v} = \\sigma\\mathbf{v}\n",
    "$$\n",
    "Let us multiply the left hand side by $\\mathbf{Y}^\\top$, then we have\n",
    "$$\n",
    "\\mathbf{Y}^\\top\\mathbf{Y}(\\mathbf{Y}^\\top \\mathbf{v}) = \\sigma(\\mathbf{Y}^\\top\\mathbf{v}).\n",
    "$$\n",
    "Note that $(\\mathbf{Y}^\\top\\mathbf{v})$ is an eigenvector of $\\mathbf{Y}^\\top\\mathbf{Y}$. Thus, we can obtain the eigenvector of $\\mathbf{Y}^\\top\\mathbf{Y}$ ($p\\times p$) via obtaining the eigenvector of $\\mathbf{Y}\\mathbf{Y}^\\top$ ($n\\times n$) and then multiply it by $\\mathbf{Y}^\\top$.\n",
    "\n",
    "<!--But if we premultiply by $\\mathbf{Y}$ then we can solve,\n",
    "$$\n",
    "\\mathbf{Y}\\mathbf{Y}^\\top\\mathbf{Y} \\mathbf{U} = \\mathbf{Y}\\mathbf{U}\\boldsymbol{\\Lambda}.\n",
    "$$\n",
    "Now let us multiply the right hand side by \\Lambda^{-\\frac{1}{2}}, we have\n",
    "$$\n",
    "\\mathbf{Y}\\mathbf{Y}^\\top(\\mathbf{Y} \\mathbf{U}\\Lambda^{-\\frac{1}{2}}) = (\\mathbf{Y}\\mathbf{U}\\Lambda^{-\\frac{1}{2}})\\boldsymbol{\\Lambda}.\n",
    "$$\n",
    "Here, note here we used the property that if $\\mathbf{A}$ and $\\mathbf{B}$ are diagonal, $\\mathbf{A}\\mathbf{B}=\\mathbf{B}\\mathbf{A}$.\n",
    "Let us define a new matrix\n",
    "$$\n",
    "\\mathbf{U}^\\prime = \\mathbf{Y} \\mathbf{U} \\Lambda^{-\\frac{1}{2}}\n",
    "$$\n",
    "We can verify that $\\mathbf{U}^\\prime$ is an orthorormal matrix because\n",
    "$$\n",
    "\\left.\\mathbf{U}^\\prime\\right.^\\top\\mathbf{U}^\\prime = \\Lambda^{-\\frac{1}{2}}\\mathbf{U}^\\top\\mathbf{Y}^\\top\\mathbf{Y} \\mathbf{U} \\Lambda^{-\\frac{1}{2}}\n",
    "$$\n",
    "and since $\\mathbf{U}$ diagonalises $\\mathbf{Y}^\\top\\mathbf{Y}$, \n",
    "$$\n",
    "\\mathbf{U}^\\top\\mathbf{Y}^\\top\\mathbf{Y} \\mathbf{U} = \\Lambda\n",
    "$$\n",
    "then \n",
    "$$\n",
    "\\left.\\mathbf{U}^\\prime\\right.^\\top\\mathbf{U}^\\prime = \\mathbf{I}\n",
    "$$\n",
    "Using this new variable, we have \n",
    "$$\n",
    "\\mathbf{Y}\\mathbf{Y}^\\top\\mathbf{U}^\\prime = \\mathbf{U}^\\prime\\boldsymbol{\\Lambda}.\n",
    "$$\n",
    "Thus, $\\mathbf{U}^\\prime$ is the (normalised) eigenvector of $\\mathbf{Y}\\mathbf{Y}^\\top$ ($n\\times n$). We can obtain $\\mathbf{U}^\\prime$ \n",
    "end -->"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# <a id='partB'></a>B: Clustering\n",
    "\n",
    "Clustering associates each data point, $\\mathbf{y}_{i, :}$ with one of $k$ different discrete groups. For example, clustering animals into discrete groups or clustering diseases into different sub-types. This is a fundamental unsupervised learning task and also essential in [data mining](https://en.wikipedia.org/wiki/Data_mining).\n",
    "\n",
    "### [$k$-means clustering](https://en.wikipedia.org/wiki/K-means_clustering)\n",
    "\n",
    "#### The $k$-means algorithm\n",
    "The simplest and most popular clustering algorithm. \n",
    "* *Require*: Set $k$ and a stopping criterion\n",
    "    1. Initialize cluster centres as randomly selected data points.\n",
    "    2. Assign each data point to *nearest* cluster centre (centroid).\n",
    "    3. Update each cluster centre by setting it to the mean of assigned data points.\n",
    "    4. Repeat 2 and 3 until the stopping criterion reached (e.g., cluster allocations do not change).\n",
    "    \n",
    "There is only one hyperparameter ($k$) to set.\n",
    "\n",
    "#### The objective function of $k$-means\n",
    "\n",
    "The standard $k$-means algorithm minimizes the objective (compactness)\n",
    "  $$\n",
    "    E=\\sum_{j=1}^k \\sum_{i\\ \\text{allocated to}\\ j}  \\left(\\mathbf{y}_{i, :} - \\boldsymbol{\\mu}_{j, :}\\right)^\\top\\left(\\mathbf{y}_{i, :} - \\boldsymbol{\\mu}_{j, :}\\right)\n",
    "  $$\n",
    "  *i.e.* it minimizes the sum of **Euclidean squared distances** betwen points and their associated centres. This is equivalent to minimizing the trace of the pooled within-cluster covariance/scatter matrix. The minimum is *not* guaranteed to be *global* or *unique*. Typically, multiple versions are run to selecte the best with a predefined metric. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### $k$-means example in Scikit-learn\n",
    "Run the Scikit-learn example [Demonstration of $k$-means assumptions](https://scikit-learn.org/stable/auto_examples/cluster/plot_kmeans_assumptions.html). More examples on $k$-means are at (the bottom of) the [$k$-means documentation](http://scikit-learn.org/stable/modules/generated/sklearn.cluster.KMeans.html)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### More clustering\n",
    "\n",
    "We only touch clustering lightly in this session. To learning more about clustering, refer to the Wikipedia [Cluster Analysis](https://en.wikipedia.org/wiki/Cluster_analysis) page and the excellent Scikit-learn [Clustering](https://scikit-learn.org/stable/modules/clustering.html) module documentation and examples, e.g. for [Spectral clustering](https://scikit-learn.org/stable/modules/clustering.html#spectral-clustering).\n",
    "\n",
    "Run the Scikit-learn example [Comparing different clustering algorithms on toy datasets](https://scikit-learn.org/stable/auto_examples/cluster/plot_cluster_comparison.html) to generate the following\n",
    "![compare_clustering](https://scikit-learn.org/stable/_images/sphx_glr_plot_cluster_comparison_001.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# <a id='partC'></a>C: Autoencoder with PyTorch"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, we perform unsupervised feature learning via [autoencoder](https://en.wikipedia.org/wiki/Autoencoder) in PyTorch. We are going to build an autoencoder to learn a low-dimensional representation of some specific images, in a particular dataset.  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Libaries\n",
    "\n",
    "Get ready by importing the standard APIs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import torch\n",
    "import torch.nn as nn\n",
    "import torch.nn.functional as F\n",
    "import torch.optim as optim\n",
    "import matplotlib.pyplot as plt\n",
    "from torchvision import datasets, transforms\n",
    "\n",
    "#Set the random seed for reproducibility \n",
    "torch.manual_seed(2020) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Data\n",
    "Let us work with the popular [MNIST dataset](https://en.wikipedia.org/wiki/MNIST_database) with handwritten digits. We will work on a small subset here for efficiency."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "mnist_data = datasets.MNIST('data', train=True, download=True, transform=transforms.ToTensor())\n",
    "print(len(mnist_data))\n",
    "mnist_data = list(mnist_data)[:4096]\n",
    "print(len(mnist_data))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Define the NN architecture\n",
    "Define an autoencoder class consisting of an **encoder** followed by a **decoder** below.\n",
    "<img src=\"https://miro.medium.com/max/3148/1*44eDEuZBEsmG_TCAKRI3Kw@2x.png\" style=\"height:5cm;\"/>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "class Autoencoder(nn.Module):\n",
    "    def __init__(self):\n",
    "        super(Autoencoder, self).__init__()\n",
    "        self.encoder = nn.Sequential(\n",
    "            # 1 input image channel, 16 output channel, 3x3 square convolution\n",
    "            nn.Conv2d(1, 16, 3, stride=2, padding=1),\n",
    "            nn.ReLU(),\n",
    "            nn.Conv2d(16, 32, 3, stride=2, padding=1),\n",
    "            nn.ReLU(),\n",
    "            nn.Conv2d(32, 64, 7)\n",
    "        )\n",
    "        self.decoder = nn.Sequential(\n",
    "            nn.ConvTranspose2d(64, 32, 7),\n",
    "            nn.ReLU(),\n",
    "            nn.ConvTranspose2d(32, 16, 3, stride=2, padding=1, output_padding=1),\n",
    "            nn.ReLU(),\n",
    "            nn.ConvTranspose2d(16, 1, 3, stride=2, padding=1, output_padding=1),\n",
    "            nn.Sigmoid()  #to range [0, 1]\n",
    "        )\n",
    "\n",
    "    def forward(self, x):\n",
    "        x = self.encoder(x)\n",
    "        x = self.decoder(x)\n",
    "        return x"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Again, `__init__()` defines the layers, using `nn.Sequential` this time.  `forward()` defines the *forward pass* that transform the input to the output. `backward()` is automatically defined using `autograd`. `ReLu()` and `Sigmoid()` are [rectified linear unit (ReLU)](https://en.wikipedia.org/wiki/Rectifier_(neural_networks)) and [Sigmoid function](https://en.wikipedia.org/wiki/Sigmoid_function), two popular **activation function** that performs a *nonlinear* transformation/mapping of an input variable (element-wise operation).\n",
    "\n",
    "Here, we have both convolution layers `Conv2d()` and transpose convolution layers `ConvTranspose2d()`, with nice illustrations at [Convolution arithmetic](https://github.com/vdumoulin/conv_arithmetic). The basic ones are reproduced below where blue maps indicate inputs, and cyan maps indicate outputs.\n",
    "\n",
    "<table>\n",
    "    <tr>\n",
    "    <td  style=\"text-align: left\"> Convolution with no padding, no strides.      <img src=\"https://raw.githubusercontent.com/vdumoulin/conv_arithmetic/master/gif/no_padding_no_strides.gif\" alt=\"Drawing\" style=\"width: 250px;\"/> </td>\n",
    "    <td  style=\"text-align: left\"> Transpose convolution with No padding, no strides.<img src=\"https://github.com/vdumoulin/conv_arithmetic/raw/master/gif/no_padding_no_strides_transposed.gif\" alt=\"Drawing\" style=\"width: 250px;\"/> </td>\n",
    "</tr>\n",
    "</table>\n",
    "\n",
    "#### Inspect the NN architecture\n",
    "\n",
    "Now let's take a look at the autoencoder built."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "myAE=Autoencoder()\n",
    "print(myAE)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let us check the (randomly initialised) parameters of this NN. Below, we check the first 2D convolution and the ReLu activiation function. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "params = list(myAE.parameters())\n",
    "print(len(params))\n",
    "print(params[0].size())  # First Conv2d's .weight\n",
    "print(params[1].size())  # First Conv2d's .bias\n",
    "print(params[1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Train the NN\n",
    "Next, we will feed data in this autoencoder to train it, i.e., learn its parameters so that the reconstruction error (the `loss`) is minimised, using the mean square error (MSE) and `Adam` optimiser. The dataset is loaded in batches to train the model. One `epoch` means one cycle through the full training dataset. The `outputs` at the end of each epoch save the orignal image and the reconstructed (decoded) image pairs for later inspection. Again, the steps are \n",
    "* Define the optimisation criteria and optimisation method.\n",
    "* Iterate through the whole dataset in batches, for a number of `epochs` till a maximum specified or a convergence criteria (e.g., successive change of loss < 0.000001)\n",
    "* In each batch processing, we \n",
    "    * do a forward pass\n",
    "    * compute the loss\n",
    "    * backpropagate the loss via `autograd`\n",
    "    * update the parameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "#Training (optimisation) parameters\n",
    "batch_size=64\n",
    "learning_rate=1e-3\n",
    "max_epochs = 20\n",
    "\n",
    "#Choose mean square error loss\n",
    "criterion = nn.MSELoss() \n",
    "#Choose the Adam optimiser\n",
    "optimizer = torch.optim.Adam(myAE.parameters(), lr=learning_rate, weight_decay=1e-5)\n",
    "#Specify how the data will be loaded in batches (with random shuffling)\n",
    "train_loader = torch.utils.data.DataLoader(mnist_data, batch_size=batch_size, shuffle=True)\n",
    "#Storage\n",
    "outputs = []\n",
    "\n",
    "#Start training\n",
    "for epoch in range(max_epochs):\n",
    "    for data in train_loader:\n",
    "        img, label = data\n",
    "        optimizer.zero_grad()\n",
    "        recon = myAE(img)\n",
    "        loss = criterion(recon, img)\n",
    "        loss.backward()\n",
    "        optimizer.step()            \n",
    "    if (epoch % 2) == 0:\n",
    "        print('Epoch:{}, Loss:{:.4f}'.format(epoch+1, float(loss)))\n",
    "    outputs.append((epoch, img, recon),)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can see the loss is being well minimised. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "### <a id='Q2'></a>Question 3\n",
    "\n",
    "When we choose to work on a subset above `mnist_data = list(mnist_data)[:4096]`, why do we choose the number 4096, e.g. rather than 4000? What's special about 4096?\n",
    "\n",
    "**Answer**: \n",
    "\n",
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can inspect the reconstructed images at different epochs below. Note here we use `.detach()` because the gradients are not needed for inspection purpose here. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "numImgs=12;\n",
    "for k in range(0, max_epochs, 9):\n",
    "    plt.figure(figsize=(numImgs, 2))\n",
    "    imgs = outputs[k][1].detach().numpy()    \n",
    "    recon = outputs[k][2].detach().numpy()\n",
    "    print('Epoch:', k+1)\n",
    "    for i, item in enumerate(imgs):\n",
    "        if i >= numImgs: break\n",
    "        plt.subplot(2, numImgs, i+1)\n",
    "        plt.imshow(item[0])\n",
    "        \n",
    "    for i, item in enumerate(recon):\n",
    "        if i >= numImgs: break\n",
    "        plt.subplot(2, numImgs, numImgs+i+1)\n",
    "        plt.imshow(item[0])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Generate synthesised images\n",
    "We can **interpolate** between two images via the learned embeddding. Let us pick the first `5` (on my machine) and third image `6`  (on my machine) from the first epoch and obtain the low-dimensional embedding using the learned encoder. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "epochIndex=0;\n",
    "img1Index=0;\n",
    "img2Index=2;\n",
    "\n",
    "imgs = outputs[epochIndex][1].detach().numpy()\n",
    "x1 = outputs[epochIndex][1][img1Index,:,:,:];# first image\n",
    "x2 = outputs[epochIndex][1][img2Index,:,:,:] # second image\n",
    "x = torch.stack([x1,x2])     # stack them together so we only call `encoder` once\n",
    "embedding = myAE.encoder(x)\n",
    "e1 = embedding[0] # embedding of first image\n",
    "e2 = embedding[1] # embedding of second image\n",
    "print(e1.size())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the embedding space, we do a linear interpolation between the two embeddings and then decode these interpolated embeddings into images."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "embedding_values = []\n",
    "for i in range(0, 10):\n",
    "    e = e1 * (i/10) + e2 * (10-i)/10\n",
    "    embedding_values.append(e)\n",
    "embedding_values = torch.stack(embedding_values)\n",
    "\n",
    "recons = myAE.decoder(embedding_values)\n",
    "plt.figure(figsize=(10, 2))\n",
    "for i, recon in enumerate(recons.detach().numpy()):\n",
    "    plt.subplot(2,10,i+1)\n",
    "    plt.imshow(recon[0])\n",
    "plt.subplot(2,10,11)\n",
    "plt.imshow(imgs[img2Index][0])\n",
    "plt.subplot(2,10,20)\n",
    "plt.imshow(imgs[img1Index][0])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "### <a id='Q3'></a>Question 4\n",
    "\n",
    "Apply PCA from `sklearn` to the training data from the MNIST training dataset above. Print the top 30 eigenvalues."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Answer: "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "### <a id='Q3'></a>Question 5\n",
    "\n",
    "Plot the **cumulative** variances captured by the top 30 PCs (plot 30 values in total, e.g., the cumulative variance for the top 5 PCs is the summation of variance captured by the top 5 PCs)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Answer: "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Additional ideas to explore\n",
    "\n",
    "* Code PCA using ``torch.nn`` and compare it with the close-form solution via eigendecomposition (e.g. that in `sklearn`).\n",
    "* Try out different optimisers or different loss function (the L1loss, MAE) and compare the results.\n",
    "* Change the architecture of autoencoder (e.g., depth, other layers such as max pooling, different activation functions) to compare the results.\n",
    "* Repeat the above on a subset from the CIFAR10 dataset. For example, you can interpolate a cat and a dog."
   ]
  }
 ],
 "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.8.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}