{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Dimensionality Reduction: Latent Variable Modelling\n",
    "\n",
    "### 17th November 2015 Neil D. Lawrence\n",
    "\n",
    "So far in our classes we have focussed mainly on regression problems, which are examples of supervised learning. We have considered the relationship between the likelihood and the objective function and we have shown how we can find paramters by maximizing the likelihood (equivalent to minimizing the objective function) and in the last session we saw how we can *marginalize* the parameters in a process known as Bayesian inference.\n",
    "\n",
    "Now we are going to turn to a different form of learning, commonly known as *unsupervised* learning. In unsupervised learning our data isn't necessarily labelled in any form, but we want models that give us a better understanding of the data. We've actually seen an example of this already with [*matrix factorization* for collaborative filtering](./week2.ipynb),  which we introduces in the context of *objective functions*. Now we will introduce a more probabilistic approach to such models, specifically we are interested in *latent variable* modelling.\n",
    "\n",
    "## Latent Variables\n",
    "\n",
    "Latent means hidden, and hidden variables are simply *unobservable* variables. The idea of a latent variable is crucial to the concept of artificial intelligence, machine learning and experimental design. A latent variable could take many forms. We might observe a man walking along a road with a large bag of clothes and we might *infer* that the man is walking to the laundrette. Our observations are a highly complex data space, the response in our eyes is processed through our visual cortex, the combination of the indidivuals limb movememnts and the direction they are walking in all conflate in our heads to cause us to infer that (perhaps) the individual is going to the laundrette. We don't *know* that the man is walking to the laundrette, but we have a model of the world that suggests that it's a likely outcome for the very complex data. In some ways the latent variable can be seen as a *compression* of this very complex scene. If I were writing a book, I might write that \"A man tripped over whilst walking to the laundrette\". In the reader's mind an image of a man, perhaps laden with dirty clothes, may occur. All these ideas come from our expectations of the world around us. We can make further inference about the man, some of it perhaps plausible others less so. The man may be going to the laundrette because his washing machine is broken, or because he doesn't have a large enough flat to have a washing machine, or because he's carrying a duvet, or because he doesn't like ironing. All of these may *increase* in probability given our observation, but they are still *latent* variables. Unless we follow the man back to his appartment, or start making other enquirires about the man, we don't know the true answer. \n",
    "\n",
    "It's clear that to do inference about any complex system we *must* include latent variables. Latent variables are extremely powerful. In robotics, they are used to represent the *state* of the robot. The state of the robot may include its position (in x, y coordinates) its speed, its direction of facing. How are *these* variables unknown to the robot? Well the robot only posesses *sensors*, it can make observations of the nearest object in a certain direction, and it may have a map of its environment. If we represent the state of the robot as its position on a map, it may be uncertain of that position. If you go walking or running in the hills around Sheffield, you can take a very high quality ordnance survey map with you. However, unless you are a really excellent orienteer, when you are far from any given landmark, you will probably be *uncertain* about your true position on the map. These states are also latent variables. \n",
    "\n",
    "In statistical analysis of experiments you try to control for each aspect of the experiment, in particular by *randomization*. So if I'm interested in the ability of a particular fertilizer to improve the yield of a particular plant I may design an experiment where I apply the fertilizer to some plants (the treatment group) and withold the fertilizer from others (the control group). I then test to see whether the yield from the treatment group is better (or worse) than the control group. I may find that I have an excellent yield for the treatment group. However, what if I'd (unknowlingly) planted all my treatment plants in a sunny part of the field, and all the control plants in a shady part of the field. That would also be a latent variable, in this case known as a *confounder*. In statistical experimental design *randomization* is used to attempt to eliminate the correlated effects of these confounders: you aim to ensure that if these confounders *do* exist their effects are not correlated with treatment and contorl. This is known as a [randomized control trial](http://en.wikipedia.org/wiki/Randomized_controlled_trial). \n",
    "\n",
    "Greek philosophers worried a great deal about what was knowable and what was unknowable. Adherents of [philosophical Skeptisism](http://en.wikipedia.org/wiki/Skepticism) were inspired by the idea that since your senses sometimes give you contradictory information, they cannot be trusted, and in extreme cases they chose to *ignore* their senses. This is an acknowledgement that very often the true state of the world cannot be known with precision. Unfortunately, these philosophers didn't have a good understanding of probability, so they were unable to encapsulate their ideas through a *degree* of belief.\n",
    "\n",
    "We often use language to express the compression of a complex behavior or patterns in a simpler way, for example we talk about motives as a useful distallation for a perhaps very complex patter of behavior. In physics we use principles of causation and simple laws to describe the world around us. Such motives or underlying principles are difficult to observe directly, our conclusions about them emerge over a period of time by observing indirect consequences of the latent variables.\n",
    "\n",
    "Epistemic uncertainty allows us to deal with these worries by associating our degree of belief about the state of the world with a probaiblity distribution. This core idea underpins state space modelling, probabilistic graphical models and the wider field of latent variable modelling. In this session we are going to explore the idea in a simple linear system and see how it relates to *factor analysis* and *principal component analysis*.\n",
    "\n",
    "## Your Personality\n",
    "\n",
    "At the beginning of the 20th century there was a great deal of interest amoungst psychologists in formalizing patterns of thought. The approach they used became known as factor analysis. The principle is that we observe a potentially high dimensional vector of characteristics about an individual. To formalize this, social scientists designed questionaires. We can envisage many questions that we may ask, but the assumption is that underlying these questions there are only a few traits that dictate the behavior. These models are known as latent trait models and the analysis is sometimes known as factor analysis. The idea is that there are a few characteristic traits that we are looking to discern. These traits or factors can be extracted by assimilating the high dimensional characteristics of the individual into a few latent factors. \n",
    "\n",
    "### Factor Analysis Model\n",
    "\n",
    "This causes us to consider a model as follows, 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, 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 from the [previous session on linear regression](./week3.ipynb). 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}_{j, :}^\\top \\mathbf{x},\n",
    "$$\n",
    "for each element of the vector valued function, where $\\mathbf{w}_{:, j}$ is the $j$th column of the matrix $\\mathbf{W}$. In that context 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. 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",
    "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 diffetent 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}$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Assignment 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",
    "$$\n",
    "f_j(\\mathbf{x}_{i, :}) = \\mathbf{w}_{j, :}^\\top \\mathbf{x}_{i, :}\n",
    "$$\n",
    "\n",
    "*10 marks*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Question 1 Answer\n",
    "\n",
    "Write your answer to the question in this box."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Latent Variables\n",
    "\n",
    "The difference between this 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",
    "\\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. We will follow that convention for this session.\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",
    "                                   \n",
    "## Principal Component Analysis\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. 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}|^{-1}} \\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. 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 *eigenvalue decomposition* of the empirical covariance matrix.\n",
    "\n",
    "## Eigenvalue Decomposition\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 you 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 your 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 we remember that $\\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"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Assignment 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.\n",
    "\n",
    "*15 marks*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Question 2 Answer\n",
    "\n",
    "Write your answer to the question in this box."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Eigenvectors of a Symmetric Matric\n",
    "\n",
    "Symmetric matrices have *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": [
    "## Probabilistic PCA\n",
    "\n",
    "In 1997 [Tipping and Bishop](http://research.microsoft.com/pubs/67218/bishop-ppca-jrss.pdf)  and [Roweis](https://www.cs.nyu.edu/~roweis/papers/empca.pdf) independently revisited Hotelling's model and considered the case where the noise variance was finite, but *shared* across all output dimensons. Their model can be thought of as a factor analysis where\n",
    "$$\n",
    "\\boldsymbol{\\Sigma} = \\sigma^2 \\mathbf{I}.\n",
    "$$\n",
    "This leads to a marginal likelihood of the form\n",
    "$$\n",
    "p(\\mathbf{Y}|\\mathbf{W}, \\sigma^2) = \\prod_{i=1}^n\\mathcal{N}(\\mathbf{y}_{i, :} | \\mathbf{0}, \\mathbf{W}\\mathbf{W}^\\top + \\sigma^2 \\mathbf{I})\n",
    "$$\n",
    "where the limit of $\\sigma^2\\rightarrow 0$ is *not* taken. This defines a proper probabilistic model. Tippping and Bishop then went on to prove that the *maximum likelihood* solution of this model with respect to $\\mathbf{W}$ is given by an eigenvalue problem. In the probabilistic PCA case the eigenvalues and eigenvectors are given as follows. \n",
    "$$\n",
    "\\mathbf{W} = \\mathbf{U}\\mathbf{L} \\mathbf{R}^\\top\n",
    "$$\n",
    "where $\\mathbf{U}$ is the eigenvectors of the empirical covariance matrix \n",
    "$$\n",
    "\\mathbf{S} = \\sum_{i=1}^n (\\mathbf{y}_{i, :} - \\boldsymbol{\\mu})(\\mathbf{y}_{i, :} - \\boldsymbol{\\mu})^\\top,\n",
    "$$ \n",
    "which can be written $\\mathbf{S} = \\frac{1}{n} \\mathbf{Y}^\\top\\mathbf{Y}$ if the data is zero mean. The matrix $\\mathbf{L}$ is diagonal and is dependent on the *eigenvalues* of $\\mathbf{S}$, $\\boldsymbol{\\Lambda}$. If the $i$th diagonal element of this matrix is given by $\\lambda_i$ then the corresponding element of $\\mathbf{L}$ is \n",
    "$$\n",
    "\\ell_i = \\sqrt{\\lambda_i - \\sigma^2}\n",
    "$$\n",
    "where $\\sigma^2$ is the noise variance. Note that if $\\sigma^2$ is larger than any particular eigenvalue, then that eigenvalue (along with its corresponding eigenvector) is *discarded* from the solution. \n",
    "\n",
    "### Python Implementation of Probabilistic PCA\n",
    "\n",
    "We will now implement this algorithm in python. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# probabilistic PCA algorithm\n",
    "def ppca(Y, q):\n",
    "    # remove mean\n",
    "    Y_cent = Y - Y.mean(0)\n",
    "    import numpy as np\n",
    "\n",
    "    # Comute covariance\n",
    "    S = np.dot(Y_cent.T, Y_cent)/Y.shape[0]\n",
    "    lambd, U = np.linalg.eig(S)\n",
    "\n",
    "    # Choose number of eigenvectors\n",
    "    sigma2 = np.sum(lambd[q:])/(Y.shape[1]-q)\n",
    "    l = np.sqrt(lambd[:q]-sigma2)\n",
    "    W = U[:, :q]*l[None, :]\n",
    "    return W, sigma2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In practice we may not wish to compute the eigenvectors of the covariance matrix directly. This is because it requires us to estimate the covariance, which involves a sum of squares term, before estimating the eigenvectors. We can estimate the eigenvectors directly either through [QR decomposition](http://en.wikipedia.org/wiki/QR_decomposition) or [singular value decomposition](http://en.wikipedia.org/wiki/Singular_value_decomposition). We saw a similar issue arise when [computing the weights in a regression problem](./week3.ipynb), where we also wished to avoid computation of $\\mathbf{X}^\\top\\mathbf{X}$ (or in the case of [nonlinear regression with basis functions](./week4.ipynb) $\\boldsymbol{\\Phi}^\\top\\boldsymbol{\\Phi}$). "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Posterior for Principal Component Analysis\n",
    "\n",
    "Under the latent variable model justification for principal component analysis, we are normally interested in inferring something about the latent variables given the data. This is the distribution,\n",
    "$$\n",
    "p(\\mathbf{x}_{i, :} | \\mathbf{y}_{i, :})\n",
    "$$\n",
    "for any given data point. Determining this density turns out to be very similar to the approach for determining the Bayesian posterior of $\\mathbf{w}$ in Bayesian linear regression, only this time we place the prior density over $\\mathbf{x}_{i, :}$ instead of $\\mathbf{w}$. The posterior is proportional to the joint density as follows,\n",
    "$$\n",
    "p(\\mathbf{x}_{i, :} | \\mathbf{y}_{i, :}) \\propto p(\\mathbf{y}_{i, :}|\\mathbf{W}, \\mathbf{x}_{i, :}, \\sigma^2) p(\\mathbf{x}_{i, :})\n",
    "$$\n",
    "And as in the Bayesian linear regression case we first consider the log posterior,\n",
    "$$\n",
    "\\log p(\\mathbf{x}_{i, :} | \\mathbf{y}_{i, :}) = \\log p(\\mathbf{y}_{i, :}|\\mathbf{W}, \\mathbf{x}_{i, :}, \\sigma^2) + \\log p(\\mathbf{x}_{i, :}) + \\text{const}\n",
    "$$\n",
    "where the constant is not dependent on $\\mathbf{x}$. As before we collect the quadratic terms in $\\mathbf{x}_{i, :}$ and we assemble them into a Gaussian density over $\\mathbf{x}$.\n",
    "$$\n",
    "\\log p(\\mathbf{x}_{i, :} | \\mathbf{y}_{i, :}) = -\\frac{1}{2\\sigma^2} (\\mathbf{y}_{i, :} - \\mathbf{W}\\mathbf{x}_{i, :})^\\top(\\mathbf{y}_{i, :} - \\mathbf{W}\\mathbf{x}_{i, :}) - \\frac{1}{2} \\mathbf{x}_{i, :}^\\top \\mathbf{x}_{i, :} + \\text{const}\n",
    "$$\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Assignment Question 3\n",
    "\n",
    "Multiply out the terms in the brackets. Then collect the quadratic term and the linear terms together. Show that the posterior has the form\n",
    "$$\n",
    "\\mathbf{x}_{i, :} | \\mathbf{W} \\sim \\mathcal{N}( \\boldsymbol{\\mu}_x, \\mathbf{C}_x)\n",
    "$$\n",
    "where \n",
    "$$\n",
    "\\mathbf{C}_x = \\left(\\sigma^{-2} \\mathbf{W}^\\top\\mathbf{W} + \\mathbf{I}\\right)^{-1}\n",
    "$$\n",
    "and \n",
    "$$\n",
    "\\boldsymbol{\\mu}_x = \\mathbf{C}_x \\sigma^{-2}\\mathbf{W}^\\top \\mathbf{y}_{i, :} \n",
    "$$\n",
    "Compare this to the posterior for the Bayesian linear regression from last week, do they have similar forms? What matches and what differs?\n",
    "\n",
    "*30 marks*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Question 3 Answer\n",
    "\n",
    "Write your answer to the question in this box. Use latex to write your derivations, show each step of your work.\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Python Implementation of the Posterior\n",
    "\n",
    "Now let's implement the system in code."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Assignment Question 4\n",
    "\n",
    "Use the values for $\\mathbf{W}$ and $\\sigma^2$ you have computed, along with the data set $\\mathbf{Y}$ to compute the posterior density over $\\mathbf{X}$. Write a function of the form\n",
    "```python\n",
    "mu_x, C_x = posterior(Y, W, sigma2)\n",
    "```\n",
    "where `mu_x` and `C_x` are the posterior mean and posterior covariance for the given $\\mathbf{Y}$. \n",
    "\n",
    "Don't forget to subtract the mean of the data `Y` inside your function before computing the posterior: remember we assumed at the beginning of our analysis that the data had been centred (i.e. the mean was removed).\n",
    "\n",
    "\n",
    "*20 marks*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Question 4 Answer Code\n",
    "# Write code for you answer to this question in this box\n",
    "# Do not delete these comments, otherwise you will get zero for this answer.\n",
    "# Make sure your code has run and the answer is correct *before* submitting your notebook for marking.\n",
    "import numpy as np\n",
    "import scipy as sp\n",
    "def posterior(Y, W, sigma2):\n",
    "    Y_cent = Y - Y.mean(0)\n",
    "    # Compute posterior over X\n",
    "    C_x = \n",
    "    mu_x = \n",
    "    return mu_x, C_x\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Numerically Stable and Efficient Version\n",
    "\n",
    "Just as we saw for [linear regression](./week3.ipynb) and [regression with basis functions](./week4) computation of a matrix such as $\\mathbf{Y}^\\top\\mathbf{Y}$ (or its centred version) can be a bad idea in terms of loss of numerical accuracy. Fortunately, we can find the eigenvalues and eigenvectors of the matrix $\\mathbf{Y}^\\top\\mathbf{Y}$ without direct computation of the matrix. This can be done with the [*singular value decomposition*](http://en.wikipedia.org/wiki/Singular_value_decomposition). The singular value decompsition takes a matrix, $\\mathbf{Z}$ and represents it in the form,\n",
    "$$\n",
    "\\mathbf{Z} = \\mathbf{U}\\boldsymbol{\\Lambda}\\mathbf{V}^\\top\n",
    "$$\n",
    "where $\\mathbf{U}$ is a matrix of orthogonal vectors in the columns, meaning $\\mathbf{U}^\\top\\mathbf{U} = \\mathbf{I}$. It has the same number of rows and columns as $\\mathbf{Z}$. The matrices $\\mathbf{\\Lambda}$ and $\\mathbf{V}$ are both square with dimensionality given by the number of columns of $\\mathbf{Z}$. The matrix $\\mathbf{\\Lambda}$ is *diagonal* and $\\mathbf{V}$ is an orthogonal matrix so $\\mathbf{V}^\\top\\mathbf{V} = \\mathbf{V}\\mathbf{V}^\\top = \\mathbf{I}$. \n",
    "\n",
    "The eigenvalues of the matrix $\\mathbf{Y}^\\top\\mathbf{Y}$ are then given by the singular values of the matrix $\\mathbf{Y}^\\top$ squared and the eigenvectors are given by $\\mathbf{U}$.\n",
    "\n",
    "### Solution for $\\mathbf{W}$\n",
    "\n",
    "Given the singular value decomposition of $\\mathbf{Y}$ then we have\n",
    "$$\n",
    "\\mathbf{W} = \\mathbf{U}\\mathbf{L}\\mathbf{R}^\\top\n",
    "$$\n",
    "where $\\mathbf{R}$ is an arbitrary rotation matrix. This implies that the posterior is given by\n",
    "$$\n",
    "\\mathbf{C}_x = \\left[\\sigma^{-2}\\mathbf{R}\\mathbf{L}^2\\mathbf{R}^\\top + \\mathbf{I}\\right]^{-1}\n",
    "$$\n",
    "because $\\mathbf{U}^\\top \\mathbf{U} = \\mathbf{I}$. Since, by convention, we normally take $\\mathbf{R} = \\mathbf{I}$ to ensure that the principal components are orthonormal we can write\n",
    "$$\n",
    "\\mathbf{C}_x = \\left[\\sigma^{-2}\\mathbf{L}^2 + \\mathbf{I}\\right]^{-1}\n",
    "$$\n",
    "which implies that $\\mathbf{C}_x$ is actually diagonal with elements given by\n",
    "$$\n",
    "c_i = \\frac{\\sigma^2}{\\sigma^2 + \\ell^2_i}\n",
    "$$\n",
    "and allows us to write\n",
    "$$\n",
    "\\boldsymbol{\\mu}_x = [\\mathbf{L}^2 + \\sigma^2 \\mathbf{I}]^{-1} \\mathbf{L} \\mathbf{U}^\\top \\mathbf{y}_{i, :}\n",
    "$$\n",
    "$$\n",
    "\\boldsymbol{\\mu}_x = \\mathbf{D}\\mathbf{U}^\\top \\mathbf{y}_{i, :}\n",
    "$$\n",
    "where $\\mathbf{D}$ is a diagonal matrix with diagonal elements given by $d_{i} = \\frac{\\ell_i}{\\sigma^2 + \\ell_i^2}$.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# probabilistic PCA algorithm using SVD\n",
    "import scipy as sp\n",
    "import numpy as np\n",
    "def ppca(Y, q, center=True):\n",
    "    \"\"\"Probabilistic PCA through singular value decomposition\"\"\"\n",
    "    # remove mean\n",
    "    if center:\n",
    "        Y_cent = Y - Y.mean(0)\n",
    "    else:\n",
    "        Y_cent = Y\n",
    "        \n",
    "    # Comute singluar values, discard 'R' as we will assume orthogonal\n",
    "    U, sqlambd, _ = sp.linalg.svd(Y_cent.T,full_matrices=False)\n",
    "    lambd = (sqlambd**2)/Y.shape[0]\n",
    "    # Compute residual and extract eigenvectors\n",
    "    sigma2 = np.sum(lambd[q:])/(Y.shape[1]-q)\n",
    "    ell = np.sqrt(lambd[:q]-sigma2)\n",
    "    return U[:, :q], ell, sigma2\n",
    "\n",
    "def posterior(Y, U, ell, sigma2, center=True):\n",
    "    \"\"\"Posterior computation for the latent variables given the eigendecomposition.\"\"\"\n",
    "    if center:\n",
    "        Y_cent = Y - Y.mean(0)\n",
    "    else:\n",
    "        Y_cent = Y\n",
    "    C_x = np.diag(sigma2/(sigma2+ell**2))\n",
    "    d = ell/(sigma2+ell**2)\n",
    "    mu_x = np.dot(Y_cent, U)*d[None, :]\n",
    "    return mu_x, C_x"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Examples\n",
    "\n",
    "For our first example we'll consider some motion capture data of a man breaking into a run. [Motion capture data](http://en.wikipedia.org/wiki/Motion_capture) involves capturing a 3-d point cloud to represent a character, often by an underlying skeleton. For this data set, from Ohio State University, we have 54 frame of motion capture, each frame containing 102 values, which are the 3-d locations of 34 different points from the subjects skeleton."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pods\n",
    "data = pods.datasets.osu_run1()\n",
    "Y = data['Y']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Once the data is loaded in we can examine the first two principal components as follows,"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "q = 2\n",
    "U, ell, sigma2 = ppca(Y, q)\n",
    "mu_x, C_x = posterior(Y, U, ell, sigma2)\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "plt.plot(mu_x[:, 0], mu_x[:, 1], 'rx-')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here because the data is a time course, we have connected points that are neighbouring in time. This highlights the form of the run, which involves 3 paces. This projects in our low dimensional space to 3 loops. We can examin how much residual variance there is in the system by looking at `sigma2`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(sigma2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Robot Navigation Example\n",
    "\n",
    "In the next example we will load in data from a robot navigation problem. The data consists of wireless access point strengths as recorded by a robot performing a loop around the University of Washington's Computer Science department in Seattle. The robot records all the wireless access points it can cache and stores their signal strength."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data = pods.datasets.robot_wireless()\n",
    "Y = data['Y']\n",
    "Y.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There are 215 observations of 30 different access points. In this case the model is suggesting that the access point signal strength should be linearly dependent on the location in the map. In other words we are expecting the access point strength for the $j$th access point at robot position $x_{i, :}$ to be represented by $y_{i, j} = \\mathbf{w}_{j, :}^\\top \\mathbf{x}_{i, :} +  \\epsilon_{i,j}$.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "q = 2\n",
    "U, ell, sigma2 = ppca(Y, q)\n",
    "mu_x, C_x = posterior(Y, U, ell, sigma2)\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "fig, ax = plt.subplots(figsize=(8,8))\n",
    "ax.plot(mu_x[:, 0], mu_x[:, 1], 'rx-')\n",
    "ax.set_title('Latent Variable: Robot Inferred Locations')\n",
    "fig, ax = plt.subplots(figsize=(8,8))\n",
    "W = U*ell[None, :]\n",
    "ax.plot(W[:, 0], W[:, 1], 'bo')\n",
    "ax.set_title('Access Point Inferred Locations')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "U, ell, sigma2 = ppca(Y.T, q)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Relationship to Matrix Factorization\n",
    "\n",
    "We can use the robot naviation example to realise that PCA (and factor analysis) are very reminiscient of the [*matrix factorization* example](./week2.ipynb) that we used for introducing objective functions. In that system we used slightly different notation, $\\mathbf{u}_{i, :}$ for *user* location in our metaphorical library and $\\mathbf{v}_{j, :}$ for *item* location in our metaphorical library. To see how these systems are somewhat analagous, now let us think about the user as the robot and the items as the wifi access points. We can plot the relative location of both. This process is known as \"SLAM\": simultaneous *localisation* and *mapping*. A latent variable model of the type we have developed is one way of performing SLAM. We have an estimate of the *landmarks* in the system (in this case WIFI access points) and we have an estimate of the robot position. These are analagous to the estimate of the user's position and the estimate of the items positions in the library. In the matrix factorisation example users are informing us what items they are 'close' to by expressing their preferences, in the robot localization example the robot is informing us what access point it is close to by measuring signal strength.\n",
    "\n",
    "From a personal perspective, I find this analogy quite comforting. I think it is very arguable that one of the mechanisms through which we (as humans) may have developed higher reasoning is through the need to navigate around our environment, identifying landmarks and associating them with our search for food. If such a system were to exist, the idea that it could be readily adapted to other domains such as categorising the nature of the different foodstuffs we were able to forage is intriguing. \n",
    "\n",
    "From an algorithmic perspective, we also can now realise that matrix factorization and latent variable modelling are effectively the same thing. The only difference is the objective function and our probabilistic (or lack of probabilistic) treatment of the variables. But the prediction functoin for both systems,\n",
    "$$\n",
    "f_{i, j} = \\mathbf{u}_{i, :}^\\top \\mathbf{v}_{j, :} \n",
    "$$\n",
    "for matrix factorization or\n",
    "$$\n",
    "f_{i, j} = \\mathbf{x}_{i, :}^\\top \\mathbf{w}_{j, :} \n",
    "$$\n",
    "for probabilistic PCA and factor analysis are the same.\n",
    "\n",
    "## Other Interpretations of PCA: Separating Model and Algorithm\n",
    "\n",
    "Since Hotelling first introduced his perspective on factor analysis as PCA there has been somewhat of a conflation of the idea of the model underlying PCA (for which it was very clear that Hotelling was inspired by Factor Analysis) and the algorithm that is used to fit that model: the eigenvalues and eigenvectors of the covariance matrix. The eigenvectors of an ellipsoid have been known since the middle of the 19th century as the principal axes of the elipsoid, and they arise through the following additional ideas: seeking the orthogonal directions of *maximum variance* in a dataset. Pearson in 1901 arrived at the same algorithm driven by a desire to seek a *symmetric regression* between two covariate/response variables $x$ and $y$. He is, therefore, often credited with the invention of principal component analysis, but to me this seems disengenous. His aim was very different from Hotellings, it was just happened that the optimal solution for his model was coincident with that of Hotelling. The approach is also known as the [Karhunen Loeve Transform](http://en.wikipedia.org/wiki/Karhunen%E2%80%93Lo%C3%A8ve_theorem)  in stochastic process theory and in classical multidimensional scaling the same operation can be shown to be minimising a particular objective function based on interpoint distances in the data and the latent space (see the section on Classical Multidimensional Scaling in [Mardia, Kent and Bibby](http://store.elsevier.com/Multivariate-Analysis/Kanti-Mardia/isbn-9780124712522/)). One of my own contributions to machine learning was deriving yet another model whose linear variant was solved by finding the principal subspace of the covariance matrix (an approach I termed dual probabilistic PCA or probabilistic principal coordinate analysis). Finally, the approach is sometimes referred to simply as singular value decomposition (SVD). The singular value decomposition of a data set has the following form,\n",
    "$$\n",
    "\\mathbf{Y} = \\mathbf{V} \\boldsymbol{\\Lambda} \\mathbf{U}^\\top\n",
    "$$\n",
    "where $\\mathbf{V}\\in\\Re^{n\\times n}$ and $\\mathbf{U}^\\in \\Re^{p\\times p}$ are square orthogonal matrices and $\\mathbf{\\Lambda}^{n \\times p}$ is zero apart from its first $p$ diagonal entries. Singularvalue decomposition gives a diagonalisation of the covariance matrix, because under the SVD we have\n",
    "$$\n",
    "\\mathbf{Y}^\\top\\mathbf{Y} = \\mathbf{U}\\boldsymbol{\\Lambda}\\mathbf{V}^\\top\\mathbf{V} \\boldsymbol{\\Lambda} \\mathbf{U}^\\top = \\mathbf{U}\\boldsymbol{\\Lambda}^2 \\mathbf{U}^\\top\n",
    "$$\n",
    "where $\\boldsymbol{\\Lambda}^2$ is now the eigenvalues of the covariane matrix and $\\mathbf{U}$ are the eigenvectors. So performing the SVD can simply be seen as another approach to determining the principal components.\n",
    "\n",
    "### Separating Model and Algorithm\n",
    "\n",
    "I've given a fair amount of personal thought to this situation and my own opinion that this confusion about method arises because of a conflation of model and algorithm. The model of Hotelling, that which he termed principal component analysis, was really a variant of factor analysis, and it was unfortunate that he chose to rename it. However, the algorithm he derived was a very convenient way of optimising a (simplified) factor analysis, and it's therefore become very popular. The algorithm is also the optimal solution for many other models of the data, even some which might seem initally to be unrelated (e.g. seeking directions of maximum variance). It is only through the mathematics of this linear system (which also contains some intersting symmetries) that all these ides become related. However, as soon as we choose to non-linearise the system (e.g. through basis functions) we find that each of the non-linear intepretations we can derive for the different models each leads to a very different algorithm (if such an algorithm is possible). For example [principal curves](http://web.stanford.edu/~hastie/Papers/Principal_Curves.pdf) of Hastie and Stuezle attempt to non-linearise the maximum variance interpretation, [kernel PCA](http://en.wikipedia.org/wiki/Kernel_principal_component_analysis) of Schoelkopf, Smola and Mueller uses basis functions to form the eigenvalue problem in a nonlinear space, and my own work in this area [non-linearises the dual probabilistic PCA](http://jmlr.org/papers/volume6/lawrence05a/lawrence05a.pdf). \n",
    "\n",
    "My conclusion is that when you are doing machine learning you should always have it clear in your mind what your *model* is and what your *algorithm* is. You can recognise your model because it normally contains a prediction function and an objective function. The algorithm on the other hand is the sequence of steps you implement on the computer to solve for the parameters of this model. For efficient implementation, we often modify our model to allow for faster algorithms, and this is a perfectly valid pragmatist's approach, so conflation of model and algorithm is not always a bad thing. But for clarity of thinking and understanding it is necessary to maintain the separation and to maintain a handle on when and why we perform the conflation.\n",
    "\n",
    "## PCA in Practice\n",
    "\n",
    "Principal component analysis is so effective in practice that there has almost developed a mini-industry in renaming the method itself (which is ironic, given its origin). In particular [Latent Semantic Indexing](http://en.wikipedia.org/wiki/Latent_semantic_indexing) in text processing is simply PCA on a particular representation of the term frequencies of the document. There is a particular fad to rename the eigenvectors after the nature of the data you are examining, perhaps initially triggered by [Turk and Pentland's](http://www.face-rec.org/algorithms/PCA/jcn.pdf) paper on eigenfaces, but also with [eigenvoices](https://wiki.inf.ed.ac.uk/twiki/pub/CSTR/ListenSemester1_2007_8/kuhn-junqua-eigenvoice-icslp1998.pdf) and [eigengenes](http://www.biomedcentral.com/1752-0509/1/54). This seems to be an instantiation of a wider, and hopefully subconcious, tendency in academia to attempt to differentiate one idea from the same idea in related fields in order to emphasise the novelty. The unfortunate result is somewhat of a confusing literature for relatively simple model. My recommendations would be as follows. If you have multivariate data, applying some form of principal component would seem to be a very good idea as a first step. Even if you intend to later perform classification or regression on your data, it can give you understanding of the structure of the underlying data and help you to develop your intuitions about the nature of your data. Intelligent plotting and interaction with your data is always a good think, and for high dimensional data that means that you need some way of visualisation, PCA is typically a good starting point.\n",
    "\n",
    "## Marginal Likelihood\n",
    "\n",
    "We have developed the posterior density over the latent variables given the data and the parameters, and due to symmetries in the underlying prediction function, it has a very similar form to its sister density, the posterior of the weights given the data from Bayesian regression. Two key differences are as follows. If we were to do a Bayesian multiple output regression we would find that the marginal likelihood of the data is independent across the features and correlated across the data,\n",
    "$$\n",
    "p(\\mathbf{Y}|\\mathbf{X}) = \\prod_{j=1}^p \\mathcal{N}(\\mathbf{y}_{:, j}|\\mathbf{0}, \\alpha\\mathbf{X}\\mathbf{X}^\\top + \\sigma^2 \\mathbf{I})\n",
    "$$\n",
    "where $\\mathbf{y}_{:, j}$ is a column of the data matrix and the independence is across the *features*, in probabilistic PCA the marginal likelihood has the form,\n",
    "$$\n",
    "p(\\mathbf{Y}|\\mathbf{W}) = \\prod_{i=1}^n \\mathcal{N}(\\mathbf{y}_{i, :}|\\mathbf{0}, \\mathbf{W}\\mathbf{W}^\\top + \\sigma^2 \\mathbf{I})\n",
    "$$\n",
    "where $\\mathbf{y}_{i, :}$ is a row of the data matrix $\\mathbf{Y}$ and the independence is across the data points.\n",
    "\n",
    "## Computation of the Log Likelihood\n",
    "\n",
    "The quality of the model can be assessed using the log likelihood of this Gaussian form.\n",
    "$$\n",
    "\\log p(\\mathbf{Y}|\\mathbf{W}) = -\\frac{n}{2} \\log \\left| \\mathbf{W}\\mathbf{W}^\\top + \\sigma^2 \\mathbf{I} \\right| -\\frac{1}{2} \\sum_{i=1}^n \\mathbf{y}_{i, :}^\\top \\left(\\mathbf{W}\\mathbf{W}^\\top + \\sigma^2 \\mathbf{I}\\right)^{-1} \\mathbf{y}_{i, :} +\\text{const}\n",
    "$$\n",
    "but this can be computed more rapidly by exploiting the low rank form of the covariance covariance, $\\mathbf{C} = \\mathbf{W}\\mathbf{W}^\\top + \\sigma^2 \\mathbf{I}$ and the fact that $\\mathbf{W} = \\mathbf{U}\\mathbf{L}\\mathbf{R}^\\top$. \n",
    "\n",
    "Specifically, we first use the decomposition of $\\mathbf{W}$ to write:\n",
    "$$\n",
    "-\\frac{n}{2} \\log \\left| \\mathbf{W}\\mathbf{W}^\\top + \\sigma^2 \\mathbf{I} \\right| = -\\frac{n}{2} \\sum_{i=1}^q \\log (\\ell_i^2 + \\sigma^2) - \\frac{n(p-q)}{2}\\log \\sigma^2 ,\n",
    "$$\n",
    "where $\\ell_i$ is the $i$th diagonal element of $\\mathbf{L}$. Next, we use the [Woodbury matrix identity](http://en.wikipedia.org/wiki/Woodbury_matrix_identity) which allows us to write the inverse as a quantity which contains another inverse in a smaller matrix:\n",
    "$$\n",
    "(\\sigma^2 \\mathbf{I} + \\mathbf{W}\\mathbf{W}^\\top)^{-1} = \\sigma^{-2}\\mathbf{I}-\\sigma^{-4}\\mathbf{W}{\\underbrace{(\\mathbf{I}+\\sigma^{-2}\\mathbf{W}^\\top\\mathbf{W})}_{\\mathbf{C}_x}}^{-1}\\mathbf{W}^\\top\n",
    "$$\n",
    "So, it turns out that the original inversion of the $p \\times p$ matrix can be done by forming a quantity which contains the inversion of a $q \\times q$ matrix which, moreover, turns out to be the quantity $\\mathbf{C}_x$ of the posterior.\n",
    "\n",
    "Now, we put everything together to obtain:\n",
    "$$\n",
    "\\log p(\\mathbf{Y}|\\mathbf{W}) = -\\frac{n}{2} \\sum_{i=1}^q \\log (\\ell_i^2 + \\sigma^2) - \\frac{n(p-q)}{2}\\log \\sigma^2 - \\frac{1}{2} \\text{tr} \\left( \\mathbf{Y}^\\top \\left( \\sigma^{-2}\\mathbf{I}-\\sigma^{-4}\\mathbf{W} \\mathbf{C}_x \\mathbf{W}^\\top \\right) \\mathbf{Y}  \\right) + \\text{const}\n",
    "$$\n",
    "where we used the fact that a scalar sum can be written as $\\sum_{i=1}^n \\mathbf{y}_{i,:}^\\top \\mathbf{K} \\mathbf{y}_{i,:} = \\text{tr}\\left( \\mathbf{Y}^\\top \\mathbf{K} \\mathbf{Y} \\right)$, for any matrix $\\mathbf{K}$ of appropriate dimensions. We now use the properties of the trace $\\text{tr}(\\mathbf{A}+\\mathbf{B})=\\text{tr}(\\mathbf{A})+\\text{tr}(\\mathbf{B})$ and $\\text{tr}(c \\mathbf{A}) = c \\text{tr}(\\mathbf{A})$, where $c$ is a scalar and $\\mathbf{A},\\mathbf{B}$ matrices of compatible sizes. Therefore, the final log likelihood takes the form:\n",
    "$$\n",
    "\\log p(\\mathbf{Y}|\\mathbf{W}) = -\\frac{n}{2} \\sum_{i=1}^q \\log (\\ell_i^2 + \\sigma^2) - \\frac{n(p-q)}{2}\\log \\sigma^2 - \\frac{\\sigma^{-2}}{2} \\text{tr}(\\mathbf{Y}^\\top \\mathbf{Y}) +\\frac{\\sigma^{-4}}{2} \\text{tr}(\\mathbf{B}\\mathbf{C}_x\\mathbf{B}^\\top) + \\text{const}\n",
    "$$\n",
    "where we also defined $\\mathbf{B}=\\mathbf{Y}^\\top\\mathbf{W}$.\n",
    "\n",
    "Finally, notice that $\\text{tr}(\\mathbf{Y}\\mathbf{Y}^\\top)=\\text{tr}(\\mathbf{Y}^\\top\\mathbf{Y})$ can be computed faster as the sum of all the elements of $\\mathbf{Y}\\circ\\mathbf{Y}$, where $\\circ$ denotes the element-wise (or [Hadamard](http://en.wikipedia.org/wiki/Hadamard_product_(matrices)) product.\n",
    "\n",
    "\n",
    "\n",
    "## Reconstruction of the Data\n",
    "\n",
    "Given any posterior projection of a data point, we can replot the original data as a function of the input space. \n",
    "\n",
    "We will now try to reconstruct the motion capture figure form some different places in the latent plot."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Assignment Question 5\n",
    "\n",
    "Project the motion capture data onto its principal components, and then use the *mean posterior estimate* to reconstruct the data from the latent variables at the data points. Use two latent dimensions. What is the sum of squares error for the reconstruction? \n",
    "\n",
    "*25 marks*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# Question 5 Answer Code\n",
    "# Write code for you answer to this question in this box\n",
    "# Do not delete these comments, otherwise you will get zero for this answer.\n",
    "# Make sure your code has run and the answer is correct *before* submitting your notebook for marking.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Other Data Sets to Explore\n",
    "\n",
    "Below there are a few other data sets from `pods` you might want to explore with PCA. Both of them have $p$>$n$ so you need to consider how to do the larger eigenvalue probleme efficiently without large demands on computer memory.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The data is actually quite high dimensional, and 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} = \\mathbf{U}\\boldsymbol{\\Lambda}\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",
    "but it turns out that we can write\n",
    "$$\n",
    "\\mathbf{U}^\\prime = \\mathbf{Y} \\mathbf{U} \\Lambda^{\\frac{1}{2}}\n",
    "$$\n",
    "where $\\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}\\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}\\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"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Olivetti Faces\n",
    "\n",
    "You too can create your own eigenfaces. In this example we load in the 'Olivetti Face' data set, a small data set of 200 faces from the [Olivetti Research Laboratory](http://en.wikipedia.org/wiki/Olivetti_Research_Laboratory). Below we load in the data and display an image of the second face in the data set (i.e., indexed by 1). "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data = pods.datasets.olivetti_glasses()\n",
    "Y = data['X'] # this data set is set up for classification, but we will model the outputs.\n",
    "lbls = data['Y']\n",
    "import matplotlib.cm as cm # import color map\n",
    "display_index = 0\n",
    "plt.imshow(np.reshape(Y[display_index, :].flatten(), (64, 64)).T, cmap = cm.Greys_r)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that to display the face we had to reshape the appropriate row of the data matrix. This is because the images are turned into vectors by stacking columns of the image on top of each other to form a vector. The operation\n",
    "```python \n",
    "im = np.reshape(Y[1, :].flatten(), (64, 64)).T\n",
    "```\n",
    "recovers the original image into a matrix `im` by using the `np.reshape` function to return the vector to a matrix.\n",
    "\n",
    "### Visualizing the Eigenvectors\n",
    "\n",
    "Each retained eigenvector is stored in the $j$th column of $\\mathbf{U}$. Each of these eigenvectors is associated with particular directions of variation in the original data. Principal component analysis implies that we can reconstruct any face by using a weighted sum of these eigenvectors where the weights for each face are given by the relevant vector of the latent variables, $\\mathbf{x}_{i, :}$ and the diagonal elements of the matrix $\\mathbf{L}$. We can visualize the eigenvectors $\\mathbf{U}$ as images by performing the same reshape operation we used to recover the image associated with a data point above. Below we do this for the first nine eigenvectors of the Olivetti Faces data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "width=3\n",
    "height=3\n",
    "q = width*height\n",
    "fig, ax = plt.subplots(width,height,figsize=(12,12))\n",
    "\n",
    "U, ell, sigma2 = ppca(Y, q)\n",
    "lat = 0\n",
    "for i in range(width):\n",
    "    for j in range(height):\n",
    "        ax[i, j].imshow(np.reshape(U[:, lat].flatten(), (64, 64)).T, cmap = cm.Greys_r)\n",
    "        ax[i, j].set_title('Principle Component ' + str(lat+1))\n",
    "        lat += 1\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Reconstruction\n",
    "\n",
    "We can now attempt to reconstruct a given training point from these eigenvectors. As we mentioned above, the reconstruction is dependent on the value of the latent variable and the weights from the matrix $\\mathbf{L}$. First let's compute the value of the latent variables for the point we want to construct. Then we'll use them to compute the weightings of each of the eigenvectors."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mu_x, C_x = posterior(Y, U, ell, sigma2)\n",
    "reconstruction_weights = mu_x[display_index, :]*ell\n",
    "print(reconstruction_weights)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This vector of reconstruction weights is applied to the 'template images' given by the eigenvectors to give us our reconstruction. Below we weight these templates and combine to form the reconstructed image, and show the comparison to the original image."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(1, 2,figsize=(12,6))\n",
    "ax[0].imshow(np.reshape(Y[display_index, :].flatten(), (64, 64)).T, cmap = cm.Greys_r)\n",
    "ax[0].set_title('Original Example no ' + str(display_index))\n",
    "ax[1].imshow(np.reshape(np.dot(U,reconstruction_weights) + Y.mean(axis=0)[None, :], (64, 64)).T, cmap = cm.Greys_r)\n",
    "ax[1].set_title('Reconstruction of Example from ' + str(len(reconstruction_weights)) + ' Latent Variables')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The quality of the reconstruction is a bit blurry, it can be improved by increasing the number of template images used (i.e. increasing the *latent dimensionality*). "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Gene Expression\n",
    "\n",
    "Each of the cells in your body stores your entire genetic code in your DNA, but at any one moment it is only 'expressing' a small portion of that code. Your cells are mainly constructed of protein, and these proteins are manufactured by first transcribing the DNA to RNA and then translating the RNA to protein. If the DNA is the cells hard drive, then one role of the RNA is to act like a cache of data that is being read from the hard drive at any one time. Gene expression arrays allow us to measure the quantity of different types of RNA in the cell, effectively analyzing what's in the cache (although we have to destroy the cell or the tissue to access it). A gene expression experiment often consists of a time course or a series of experiments that characterise the gene expression of cells at any given time.\n",
    "\n",
    "We will now load in one of the earliest gene expression data sets from a [1998 paper by Spellman et al.](http://www.ncbi.nlm.nih.gov/pubmed/9843569), it consists of gene expression measurements of over six thousand genes in a range of conditions for brewer's yeast. The experiment was designed for understanding the cell cycle of the genes. The expectation is that there should be oscillating signals inside the cell.\n",
    "\n",
    "First we extract the principale components of the gene expression."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# load in data and replace missing values with zero\n",
    "data=pods.datasets.spellman_yeast_cdc15()\n",
    "Y = data['Y'].fillna(0)\n",
    "q = 5\n",
    "U, ell, sigma2 = ppca(Y, q)\n",
    "mu_x, C_x = posterior(Y, U, ell, sigma2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, looking through, we find that there is indeed a latent variable that appears to oscilate at approximately the right frequency. The 4th latent dimension (`index=3`) can be plotted across the time course as follows."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(mu_x[:, 3])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To reveal an oscillating shape. We can see which genes correspond to this shape by examining the associated column of $\\mathbf{U}$. Let's augment our data matrix with this score."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gene_list = Y.T\n",
    "gene_list['oscilation'] = np.sqrt(U[:, 3]**2)\n",
    "gene_list.sort(columns='oscilation', ascending=False).index[:4]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can look up the first three genes in this list which now ranks the genes according to how strongly they respond to the fourth latent dimension. [The NCBI gene database](http://www.ncbi.nlm.nih.gov/gene/) allows us to search for the function of these genes. Looking at the function of the four genes that respond most strongly to the third latent variable they are all genes that encode [histone](http://en.wikipedia.org/wiki/Histone) proteins. The histone is the support scaffold for the DNA that ensures it folds correctly within the cell creating the nucleosomes. It seems to make sense that production of histone proteins should be strongly correlated with the cell cycle, as when the cell divides it needs to create a large number of histone proteins for creating the duplicated nucleosomes. The relevant links giving the descriptions of each gene given here: [YDR224C](http://www.ncbi.nlm.nih.gov/gene/851810), [YDR225W](http://www.ncbi.nlm.nih.gov/gene/851811), [YBL003C](http://www.ncbi.nlm.nih.gov/gene/852283) and [YNL030W](http://www.ncbi.nlm.nih.gov/gene/855701)."
   ]
  }
 ],
 "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.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}