Lab 3 - Linear Regression.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Linear Algebra and Linear Regression\n",
    "\n",
    "### Modified by Mauricio Álvarez, 13th October 2019\n",
    "\n",
    "### 13th October 2015 Neil Lawrence\n",
    "\n",
    "\n",
    "## Sum of Squares Error\n",
    "\n",
    "Last week we considered a cost function for minimization of the error. We considered items (films) and users and assumed that each movie rating, $y_{i,j}$ could be summarised by an inner product between a vector associated with the item, $\\mathbf{v}_j$ and one associated with the user $\n",
    "\\mathbf{u}_i$. We justified the inner product as a measure of similarity in the space of 'movie subjects', where both the users and the items lived, giving the analogy of a library.\n",
    "\n",
    "To make predictions we encouraged the similarity to be high if the movie rating was high using the quadratic error function,\n",
    "$$\n",
    "E_{i,j}(\\mathbf{u}_i, \\mathbf{v}_j) = \\left(\\mathbf{u}_i^\\top \\mathbf{v}_j - y_{i,j}\\right)^2,\n",
    "$$\n",
    "which we then summed across all the observations to form the total error\n",
    "$$\n",
    "E(\\mathbf{U}, \\mathbf{V}) = \\sum_{i,j}s_{i,j}\\left(\\mathbf{u}_i^\\top \\mathbf{v}_j - y_{i,j}\\right)^2,\n",
    "$$\n",
    "where $s_{i,j}$ is an indicator variable which is set to 1 if the rating of movie $j$ by user $i$ is provided in our data set. \n",
    "\n",
    "This is known as a sum of squares error. Minimizing it was first proposed by [Legendre](http://en.wikipedia.org/wiki/Adrien-Marie_Legendre) in 1805. His book, which was on the orbit of comets, is available on google books, we can take a look at the relevant page by calling the code below. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pods\n",
    "pods.notebook.display_google_book(id='spcAAAAAMAAJ', page=72) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Of course, the main text is in French, but the key part we are interested in can be roughly translated as\n",
    "\n",
    "\"In most matters where we take measures data through observation, the most accurate results they can offer, it is almost always leads to a system of equations of the form\n",
    "$$E = a + bx + cy + fz + etc .$$\n",
    "where a, b, c, f etc are the known coefficients and  x , y, z etc are unknown and must be determined by the condition that the value of E is reduced, for each equation, to an amount or zero or very small.\"\n",
    "\n",
    "He continues\n",
    "\n",
    "\"Of all the principles that we can offer for this item, I think it is not broader, more accurate, nor easier than the one we have used in previous research application, and that is to make the minimum sum of the squares of the errors. By this means, it is between the errors a kind of balance that prevents extreme to prevail, is very specific to make known the state of the closest to the truth system. The sum of the squares of the errors $E^2 + \\left.E^\\prime\\right.^2 + \\left.E^{\\prime\\prime}\\right.^2 + etc$ being\n",
    "\\begin{align*}   &(a + bx + cy + fz + etc)^2 \\\\\n",
    "+ &(a^\\prime + b^\\prime x + c^\\prime y + f^\\prime z + etc ) ^2\\\\\n",
    "+ &(a^{\\prime\\prime} + b^{\\prime\\prime}x  + c^{\\prime\\prime}y +  f^{\\prime\\prime}z + etc )^2 \\\\\n",
    "+ & etc\n",
    "\\end{align*}\n",
    "if we wanted a minimum, by varying x alone, we will have the equation ...\"\n",
    "\n",
    "This is the earliest know printed version of the problem of least squares. The notation, however, is a little awkward for mordern eyes. In particular Legendre doesn't make use of the sum sign,\n",
    "$$\n",
    "\\sum_{i=1}^3 z_i = z_1 + z_2 + z_3\n",
    "$$\n",
    "nor does he make use of the inner product. \n",
    "\n",
    "In our notation, if we were to do linear regression, we would need to subsititue:\n",
    "\\begin{align*}\n",
    "a &\\leftarrow y_1-c, \\\\ a^\\prime &\\leftarrow y_2-c,\\\\ a^{\\prime\\prime} &\\leftarrow y_3 -c,\\\\ \n",
    "\\text{etc.} \n",
    "\\end{align*}\n",
    "to introduce the data observations $\\{y_i\\}_{i=1}^{n}$ alongside $c$, the offset. We would then introduce the input locations\n",
    "\\begin{align*}\n",
    "b & \\leftarrow x_1,\\\\\n",
    "b^\\prime & \\leftarrow x_2,\\\\\n",
    "b^{\\prime\\prime} & \\leftarrow x_3,\\\\\n",
    "\\text{etc.}\n",
    "\\end{align*}\n",
    "and finally the gradient of the function\n",
    "$$x \\leftarrow -m.$$\n",
    "The remaining coefficients ($c$ and $f$) would then be zero. That would give us \n",
    "\\begin{align*}   &(y_1 - (mx_1+c))^2 \\\\\n",
    "+ &(y_2 -(mx_2 + c))^2\\\\\n",
    "+ &(y_3 -(mx_3 + c))^2 \\\\\n",
    "+ & \\text{etc.}\n",
    "\\end{align*}\n",
    "which we would write in the modern notation for sums as\n",
    "$$\n",
    "\\sum_{i=1}^n (y_i-(mx_i + c))^2\n",
    "$$\n",
    "which is recognised as the sum of squares error for a linear regression.\n",
    "\n",
    "This shows the advantage of modern [summation operator](http://en.wikipedia.org/wiki/Summation), $\\sum$,  in keeping our mathematical notation compact. Whilst it may look more complicated the first time you see it, understanding the mathematical rules that go around it, allows us to go much further with the notation.\n",
    "\n",
    "Inner products (or [dot products](http://en.wikipedia.org/wiki/Dot_product)) are similar. They allow us to write\n",
    "$$\n",
    "\\sum_{i=1}^q u_i v_i\n",
    "$$\n",
    "in a more compact notation,\n",
    "$\n",
    "\\mathbf{u}\\cdot\\mathbf{v}.\n",
    "$\n",
    "\n",
    "Here we are using bold face to represent vectors, and we assume that the individual elements of a vector $\\mathbf{z}$ are given as a series of scalars\n",
    "$$\n",
    "\\mathbf{z} = \\begin{bmatrix} z_1\\\\ z_2\\\\ \\vdots\\\\ z_n \\end{bmatrix}\n",
    "$$\n",
    "which are each indexed by their position in the vector.\n",
    "\n",
    "## Linear Algebra\n",
    "\n",
    "Linear algebra provides a very similar role, when we introduce [linear algebra](http://en.wikipedia.org/wiki/Linear_algebra), it is because we are faced with a large number of addition and multiplication operations. These operations need to be done together and would be very tedious to write down as a group. So the first reason we reach for linear algebra is for a more compact representation of our mathematical formulae. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Running Example: Olympic Marathons\n",
    "\n",
    "Now we will load in the Olympic marathon data. This is data of the olympic marath times for the men's marathon from the first olympics in 1896 up until the London 2012 olympics."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data = pods.datasets.olympic_marathon_men()\n",
    "x = data['X']\n",
    "y = data['Y']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can see what these values are by typing:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(x)\n",
    "print(y)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that they are not `pandas` data frames for this example, they are just arrays of dimensionality $n\\times 1$, where $n$ is the number of data."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The aim of this lab is to have you coding linear regression in python. We will do it in two ways, once using iterative updates (coordinate ascent) and then using linear algebra. The linear algebra approach will not only work much better, it is easy to extend to multiple input linear regression and *non-linear* regression using basis functions."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Plotting the Data\n",
    "\n",
    "You can make a plot of $y$ vs $x$ with the following command:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline \n",
    "import pylab as plt\n",
    "\n",
    "plt.plot(x, y, 'rx')\n",
    "plt.xlabel('year')\n",
    "plt.ylabel('pace in min/km')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Maximum Likelihood: Iterative Solution"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we will take the maximum likelihood approach we derived in the lecture to fit a line, $y_i=mx_i + c$, to the data you've plotted. We are trying to minimize the error function:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$E(m, c) =  \\sum_{i=1}^n(y_i-mx_i-c)^2$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "with respect to $m$ and $c$. We can start with an initial guess for $m$, "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "m = -0.4\n",
    "c = 80 "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Then we use the maximum likelihood update to find an estimate for the offset, $c$.\n",
    "\n",
    "### Coordinate Descent\n",
    "\n",
    "In the movie recommender system example, we minimised the objective function by steepest descent based gradient methods. Our updates required us to compute the gradient at the position we were located, then to update the gradient according to the direction of steepest descent. This time, we will take another approach. It is known as *coordinate descent*. In coordinate descent, we choose to move one parameter at a time. Ideally, we design an algorithm that at each step moves the parameter to its minimum value. At each step we choose to move the individual parameter to its minimum.\n",
    "\n",
    "To find the minimum, we look for the point in the curve where the gradient is zero. This can be found by taking the gradient of $E(m,c)$ with respect to the parameter. \n",
    "\n",
    "#### Update for Offset\n",
    "\n",
    "Let's consider the parameter $c$ first. The gradient goes nicely through the summation operator, and we obtain\n",
    "$$\n",
    "\\frac{\\text{d}E(m,c)}{\\text{d}c} = -\\sum_{i=1}^n 2(y_i-mx_i-c).\n",
    "$$\n",
    "Now we want the point that is a minimum. A minimum is an example of a [*stationary point*](http://en.wikipedia.org/wiki/Stationary_point), the stationary points are those points of the function where the gradient is zero. They are found by solving the equation for $\\frac{\\text{d}E(m,c)}{\\text{d}c} = 0$. Substituting in to our gradient, we can obtain the following equation, \n",
    "$$\n",
    "0 = -\\sum_{i=1}^n 2(y_i-mx_i-c)\n",
    "$$\n",
    "which can be reorganised as follows,\n",
    "$$\n",
    "c^* = \\frac{\\sum_{i=1}^n(y_i-m^*x_i)}{n}.\n",
    "$$\n",
    "The fact that the stationary point is easily extracted in this manner implies that the solution is *unique*. There is only one stationary point for this system. Traditionally when trying to determine the type of stationary point we have encountered we now compute the *second derivative*,\n",
    "$$\n",
    "\\frac{\\text{d}^2E(m,c)}{\\text{d}c^2} = 2n.\n",
    "$$\n",
    "The second derivative is positive, which in turn implies that we have found a minimum of the function. This means that setting $c$ in this way will take us to the lowest point along that axes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# set c to the minimum\n",
    "c = (y - m*x).mean()\n",
    "print(c)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Update for Slope\n",
    "\n",
    "Now we have the offset set to the minimum value, in coordinate descent, the next step is to optimise another parameter. Only one further parameter remains. That is the slope of the system. \n",
    "\n",
    "Now we can turn our attention to the slope. We once again peform the same set of computations to find the minima. We end up with an update equation of the following form.\n",
    "\n",
    "$$m^* = \\frac{\\sum_{i=1}^n (y_i - c)x_i}{\\sum_{i=1}^n x_i^2}$$\n",
    "\n",
    "Communication of mathematics in data science is an essential skill, in a moment, you will be asked to rederive the equation above. Before we do that, however, we will briefly review how to write mathematics in the notebook.\n",
    "\n",
    "### $\\LaTeX$ for Maths\n",
    "\n",
    "These cells use [Markdown format](http://en.wikipedia.org/wiki/Markdown). You can include maths in your markdown using [$\\LaTeX$ syntax](http://en.wikipedia.org/wiki/LaTeX), all you have to do is write your answer inside dollar signs, as follows:\n",
    "\n",
    "To write a fraction, we write `$\\frac{a}{b}$`, and it will display like this $\\frac{a}{b}$. To write a subscript we write `$a_b$` which will appear as $a_b$. To write a superscript (for example in a polynomial) we write `$a^b$` which will appear as $a^b$. There are lots of other macros as well, for example we can do greek letters such as `$\\alpha, \\beta, \\gamma$` rendering as $\\alpha, \\beta, \\gamma$. And we can do sum and intergral signs as `$\\sum \\int \\int$`.\n",
    "\n",
    "You can combine many of these operations together for composing expressions."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Question 1 \n",
    "\n",
    "Convert the following python code expressions into $\\LaTeX$j, writing your answers below. In each case write your answer as a single equality (i.e. your maths should only contain one expression, not several lines of expressions). For the purposes of your $\\LaTeX$ please assume that `x` and `w` are $n$ dimensional vectors. \n",
    "\n",
    "(a) \n",
    "``` python\n",
    "f = x.sum()\n",
    "```\n",
    "\n",
    "(b) \n",
    "``` python \n",
    "m = x.mean()\n",
    "```\n",
    "\n",
    "(c) \n",
    "``` python\n",
    "g = (x*w).sum()\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Question 1 Answer\n",
    "\n",
    "Write your answer to the question in this box."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Gradient With Respect to the Slope\n",
    "Now that you've had a little training in writing maths with $\\LaTeX$, we will be able to use it to answer questions. The next thing we are going to do is a little differentiation practice. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Question 2\n",
    "\n",
    "Derive the the gradient of the objective function with respect to the slope, $m$. Rearrange it to show that the update equation written above does find the stationary points of the objective function. By computing its derivative show that it's a minimum."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Question 2 Answer\n",
    "\n",
    "Write your answer to the question in this box."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "m = ((y - c)*x).sum()/(x**2).sum()\n",
    "print(m)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can have a look at how good our fit is by computing the prediction across the input space. First create a vector of 'test points',"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "x_test = np.linspace(1890, 2020, 130)[:, None]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now use this vector to compute some test predictions,"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "f_test = m*x_test + c"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now plot those test predictions with a blue line on the same plot as the data,"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(x_test, f_test, 'b-')\n",
    "plt.plot(x, y, 'rx')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The fit isn't very good, we need to iterate between these parameter updates in a loop to improve the fit, we have to do this several times,"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for i in np.arange(10):\n",
    "    m = ((y - c)*x).sum()/(x*x).sum()\n",
    "    c = (y-m*x).sum()/y.shape[0]\n",
    "print(m)\n",
    "print(c)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And let's try plotting the result again"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "f_test = m*x_test + c\n",
    "plt.plot(x_test, f_test, 'b-')\n",
    "plt.plot(x, y, 'rx')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Clearly we need more iterations than 10! In the next question you will add more iterations and report on the error as optimisation proceeds. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Question 3\n",
    "\n",
    "There is a problem here, we seem to need many interations to get to a good solution. Let's explore what's going on. Write code which alternates between updates of `c` and `m`. Include the following features in your code.\n",
    "\n",
    "(a) Initialise with `m=-0.4` and `c=80`.\n",
    "(b) Every 10 iterations compute the value of the objective function for the training data and print it to the screen (you'll find hints on this in the lab from last week.\n",
    "(c) Cause the code to stop running when the error change over less than 10 iterations is smaller than $1\\times10^{-4}$. This is known as a stopping criterion.\n",
    "\n",
    "Why do we need so many iterations to get to the solution?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Question 3 Answer Code\n",
    "# Write code for you answer to this question in this box"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Question 3 Answer Here\n",
    "\n",
    "Write your answer to the question in this box."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Multiple Input Solution with Linear Algebra\n",
    "\n",
    "You've now seen how slow it can be to perform a coordinate ascent on a system. Another approach to solving the system (which is not always possible, particularly in *non-linear* systems) is to go direct to the minimum. To do this we need to introduce *linear algebra*. We will represent all our errors and functions in the form of linear algebra. \n",
    "\n",
    "As we mentioned above, linear algebra is just a shorthand for performing lots of multiplications and additions simultaneously. What does it have to do with our system then? Well the first thing to note is that the linear function we were trying to fit has the following form:\n",
    "$$\n",
    "f(x) = mx + c\n",
    "$$\n",
    "the classical form for a straight line. From a linear algebraic perspective we are looking for multiplications and additions. We are also looking to separate our parameters from our data. The data is the *givens* remember, in French the word is données literally translated means *givens* that's great, because we don't need to change the data, what we need to change are the parameters (or variables) of the model. In this function the data comes in through $x$, and the parameters are $m$ and $c$. \n",
    "\n",
    "What we'd like to create is a vector of parameters and a vector of data. Then we could represent the system with vectors that represent the data, and vectors that represent the parameters. \n",
    "\n",
    "We look to turn the multiplications and additions into a linear algebraic form, we have one multiplication ($m\\times c$) and one addition ($mx + c$). But we can turn this into a inner product by writing it in the following way,\n",
    "$$\n",
    "f(x) = m \\times x + c \\times 1,\n",
    "$$\n",
    "in other words we've extracted the unit value, from the offset, $c$. We can think of this unit value like an extra item of data, because it is always given to us, and it is always set to 1 (unlike regular data, which is likely to vary!). We can therefore write each input data location, $\\mathbf{x}$, as a vector\n",
    "$$\n",
    "\\mathbf{x} = \\begin{bmatrix} 1\\\\ x\\end{bmatrix}.\n",
    "$$\n",
    "\n",
    "Now we choose to also turn our parameters into a vector. The parameter vector will be defined to contain \n",
    "$$\n",
    "\\mathbf{w} = \\begin{bmatrix} c \\\\ m\\end{bmatrix}\n",
    "$$\n",
    "because if we now take the inner product between these to vectors we recover\n",
    "$$\n",
    "\\mathbf{x}\\cdot\\mathbf{w} = 1 \\times c + x \\times m = mx + c\n",
    "$$\n",
    "In `numpy` we can define this vector as follows"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# define the vector w\n",
    "w = np.zeros(shape=(2, 1))\n",
    "w[0] = m\n",
    "w[1] = c"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This gives us the equivalence between original operation and an operation in vector space. Whilst the notation here isn't a lot shorter, the beauty is that we will be able to add as many features as we like and still keep the same representation. In general, we are now moving to a system where each of our predictions is given by an inner product. When we want to represent a linear product in linear algebra, we tend to do it with the transpose operation, so since we have $\\mathbf{a}\\cdot\\mathbf{b} = \\mathbf{a}^\\top\\mathbf{b}$ we can write\n",
    "$$\n",
    "f(\\mathbf{x}_i) = \\mathbf{x}_i^\\top\\mathbf{w}.\n",
    "$$\n",
    "Where we've assumed that each data point, $\\mathbf{x}_i$, is now written by appending a 1 onto the original vector\n",
    "$$\n",
    "\\mathbf{x}_i = \n",
    "\\begin{bmatrix} \n",
    "1 \\\\\n",
    "x_i\n",
    "\\end{bmatrix}\n",
    "$$\n",
    "\n",
    "## Design Matrix\n",
    "\n",
    "We can do this for the entire data set to form a [*design matrix*](http://en.wikipedia.org/wiki/Design_matrix) $\\mathbf{X}$,\n",
    "\n",
    "$$\\mathbf{X} = \\begin{bmatrix} \n",
    "\\mathbf{x}_1^\\top \\\\\\ \n",
    "\\mathbf{x}_2^\\top \\\\\\ \n",
    "\\vdots \\\\\\\n",
    "\\mathbf{x}_n^\\top\n",
    "\\end{bmatrix} = \\begin{bmatrix}\n",
    "1 & x_1 \\\\\\\n",
    "1 & x_2 \\\\\\\n",
    "\\vdots & \\vdots \\\\\\\n",
    "1 & x_n \n",
    "\\end{bmatrix},$$\n",
    "\n",
    "which in `numpy` can be done with the following commands:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "X = np.hstack((np.ones_like(x), x))\n",
    "print(X)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Writing the Objective with Linear Algebra\n",
    "\n",
    "When we think of the objective function, we can think of it as the errors where the error is defined in a similar way to what it was in Legendre's day $y_i - f(\\mathbf{x}_i)$, in statistics these errors are also sometimes called [*residuals*](http://en.wikipedia.org/wiki/Errors_and_residuals_in_statistics). So we can think as the objective and the prediction function as two separate parts, first we have,\n",
    "$$\n",
    "E(\\mathbf{w}) = \\sum_{i=1}^n (y_i - f(\\mathbf{x}_i; \\mathbf{w}))^2,\n",
    "$$\n",
    "where we've made the function $f(\\cdot)$'s dependence on the parameters $\\mathbf{w}$ explicit in this equation. Then we have the definition of the function itself,\n",
    "$$\n",
    "f(\\mathbf{x}_i; \\mathbf{w}) = \\mathbf{x}_i^\\top \\mathbf{w}.\n",
    "$$\n",
    "Let's look again at these two equations and see if we can identify any inner products. The first equation is a sum of squares, which is promising. Any sum of squares can be represented by an inner product,\n",
    "$$\n",
    "a = \\sum_{i=1}^{k} b^2_i = \\mathbf{b}^\\top\\mathbf{b},\n",
    "$$\n",
    "so if we wish to represent $E(\\mathbf{w})$ in this way, all we need to do is convert the sum operator to an inner product. We can get a vector from that sum operator by placing both $y_i$ and $f(\\mathbf{x}_i; \\mathbf{w})$ into vectors, which we do by defining \n",
    "$$\n",
    "\\mathbf{y} = \\begin{bmatrix}y_1\\\\y_2\\\\ \\vdots \\\\ y_n\\end{bmatrix}\n",
    "$$\n",
    "and defining\n",
    "$$\n",
    "\\mathbf{f}(\\mathbf{X}; \\mathbf{w}) = \\begin{bmatrix}f(\\mathbf{x}_1; \\mathbf{w})\\\\f(\\mathbf{x}_2; \\mathbf{w})\\\\ \\vdots \\\\ f(\\mathbf{x}_n; \\mathbf{w})\\end{bmatrix}.\n",
    "$$\n",
    "The second of these is actually a vector-valued function. This term may appear intimidating, but the idea is straightforward. A vector valued function is simply a vector whose elements are themselves defined as *functions*, i.e. it is a vector of functions, rather than a vector of scalars. The idea is so straightforward, that we are going to ignore it for the moment, and barely use it in the derivation. But it will reappear later when we introduce *basis functions*. So we will, for the moment, ignore the dependence of $\\mathbf{f}$ on $\\mathbf{w}$ and $\\mathbf{X}$ and simply summarise it by a vector of numbers\n",
    "$$\n",
    "\\mathbf{f} = \\begin{bmatrix}f_1\\\\f_2\\\\ \\vdots \\\\ f_n\\end{bmatrix}.\n",
    "$$\n",
    "This allows us to write our objective in the folowing, linear algebraic form,\n",
    "$$\n",
    "E(\\mathbf{w}) = (\\mathbf{y} - \\mathbf{f})^\\top(\\mathbf{y} - \\mathbf{f})\n",
    "$$\n",
    "from the rules of inner products.\n",
    "\n",
    "But what of our matrix $\\mathbf{X}$ of input data? At this point, we need to dust off [*matrix-vector multiplication*](http://en.wikipedia.org/wiki/Matrix_multiplication). Matrix multiplication is simply a convenient way of performing many inner products together, and it's exactly what we need to summarise the operation\n",
    "$$\n",
    "f_i = \\mathbf{x}_i^\\top\\mathbf{w}.\n",
    "$$\n",
    "This operation tells us that each element of the vector $\\mathbf{f}$ (our vector valued function) is given by an inner product between $\\mathbf{x}_i$ and $\\mathbf{w}$. In other words it is a series of inner products. Let's look at the definition of matrix multiplication, it takes the form\n",
    "$$\n",
    "\\mathbf{c} = \\mathbf{B}\\mathbf{a}\n",
    "$$\n",
    "where $\\mathbf{c}$ might be a $k$ dimensional vector (which we can intepret as a $k\\times 1$ dimensional matrix), and $\\mathbf{B}$ is a $k\\times k$ dimensional matrix and $\\mathbf{a}$ is a $k$ dimensional vector ($k\\times 1$ dimensional matrix). \n",
    "\n",
    "The result of this multiplication is of the form\n",
    "$$\n",
    "\\begin{bmatrix}c_1\\\\c_2 \\\\ \\vdots \\\\ a_k\\end{bmatrix} = \n",
    "\\begin{bmatrix} b_{1,1} & b_{1, 2} & \\dots & b_{1, k} \\\\\n",
    "b_{2, 1} & b_{2, 2} & \\dots & b_{2, k} \\\\\n",
    "\\vdots & \\vdots & \\ddots & \\vdots \\\\\n",
    "b_{k, 1} & b_{k, 2} & \\dots & b_{k, k} \\end{bmatrix} \\begin{bmatrix}a_1\\\\a_2 \\\\ \\vdots\\\\ c_k\\end{bmatrix} = \\begin{bmatrix} b_{1, 1}a_1 + b_{1, 2}a_2 + \\dots + b_{1, k}a_k\\\\\n",
    "b_{2, 1}a_1 + b_{2, 2}a_2 + \\dots + b_{2, k}a_k \\\\ \n",
    "\\vdots\\\\ \n",
    "b_{k, 1}a_1 + b_{k, 2}a_2 + \\dots + b_{k, k}a_k\\end{bmatrix}\n",
    "$$\n",
    "so we see that each element of the result, $\\mathbf{a}$ is simply the inner product between each *row* of $\\mathbf{B}$ and the vector $\\mathbf{c}$. Because we have defined each element of $\\mathbf{f}$ to be given by the inner product between each *row* of the design matrix and the vector $\\mathbf{w}$ we now can write the full operation in one matrix multiplication,\n",
    "$$\n",
    "\\mathbf{f} = \\mathbf{X}\\mathbf{w}.\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "f = np.dot(X, w) # np.dot does matrix multiplication in python"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Combining this result with our objective function,\n",
    "$$\n",
    "E(\\mathbf{w}) = (\\mathbf{y} - \\mathbf{f})^\\top(\\mathbf{y} - \\mathbf{f})\n",
    "$$\n",
    "we find we have defined the *model* with two equations. One equation tells us the form of our predictive function and how it depends on its parameters, the other tells us the form of our objective function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "resid = (y-f)\n",
    "E = np.dot(resid.T, resid) # matrix multiplication on a single vector is equivalent to a dot product.\n",
    "print(\"Error function is:\", E)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Question 4\n",
    "\n",
    "The prediction for our movie recommender system had the form\n",
    "$$\n",
    "f_{i,j} = \\mathbf{u}_i^\\top \\mathbf{v}_j\n",
    "$$\n",
    "and the objective function was then\n",
    "$$\n",
    "E = \\sum_{i,j} s_{i,j}(y_{i,j} - f_{i, j})^2\n",
    "$$\n",
    "Try writing this down in matrix and vector form. How many of the terms can you do? For each variable and parameter carefully think about whether it should be represented as a matrix or vector. Do as many of the terms as you can. Use $\\LaTeX$ to give your answers and give the *dimensions* of any matrices you create."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Question 4 Answer\n",
    "\n",
    "Write your answer to the question in this box."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Objective Optimisation\n",
    "\n",
    "Our *model* has now been defined with two equations, the prediction function and the objective function. Next we will use multivariate calculus to define an *algorithm* to fit the model. The separation between model and algorithm is important and is often overlooked. Our model contains a function that shows how it will be used for prediction, and a function that describes the objective function we need to optimise to obtain a good set of parameters. \n",
    "\n",
    "The linear regression model we have described is still the same as the one we fitted above with a coordinate ascent algorithm. We have only played with the notation to obtain the same model in a matrix and vector notation. However, we will now fit this model with a different algorithm, one that is much faster. It is such a widely used algorithm that from the end user's perspective it doesn't even look like an algorithm, it just appears to be a single operation (or function). However, underneath the computer calls an algorithm to find the solution. Further, the algorithm we obtain is very widely used, and because of this it turns out to be highly optimised.\n",
    "\n",
    "Once again we are going to try and find the stationary points of our objective by finding the *stationary points*. However, the stationary points of a multivariate function, are a little bit more complext to find. Once again we need to find the point at which the derivative is zero, but now we need to use  *multivariate calculus* to find it. This involves learning a few additional rules of differentiation (that allow you to do the derivatives of a function with respect to  vector), but in the end it makes things quite a bit easier. We define vectorial derivatives as follows,\n",
    "$$\n",
    "\\frac{\\text{d}E(\\mathbf{w})}{\\text{d}\\mathbf{w}} = \\begin{bmatrix}\\frac{\\partial E(\\mathbf{w})}{\\partial w_1}\\\\\\frac{\\partial E(\\mathbf{w})}{\\partial w_2}\\end{bmatrix}.\n",
    "$$\n",
    "where $\\frac{\\partial E(\\mathbf{w})}{\\partial w_1}$ is the [partial derivative](http://en.wikipedia.org/wiki/Partial_derivative) of the error function with respect to $w_1$.\n",
    "\n",
    "Differentiation through multiplications and additions is relatively straightforward, and since linear algebra is just multiplication and addition, then its rules of diffentiation are quite straightforward too, but slightly more complex than regular derivatives. \n",
    "\n",
    "### Matrix Differentiation\n",
    "\n",
    "We will need two rules of differentiation. The first is diffentiation of an inner product. By remebering that the inner product is made up of multiplication and addition, we can hope that its derivative is quite straightforward, and so it proves to be. We can start by thinking about the definition of the inner product,\n",
    "$$\n",
    "\\mathbf{a}^\\top\\mathbf{z} = \\sum_{i} a_i z_i,\n",
    "$$\n",
    "which if we were to take the derivative with respect to $z_k$ would simply return the gradient of the one term in the sum for which the derivative was non zero, that of $a_k$, so we know that \n",
    "$$\n",
    "\\frac{\\text{d}}{\\text{d}z_k} \\mathbf{a}^\\top \\mathbf{z} = a_k\n",
    "$$\n",
    "and by our definition of multivariate derivatives we can simply stack all the partial derivatives of this form in a vector to obtain the result that\n",
    "$$\n",
    "\\frac{\\text{d}}{\\text{d}\\mathbf{z}} \\mathbf{a}^\\top \\mathbf{z} = \\mathbf{a}.\n",
    "$$\n",
    "The second rule that's required is differentiation of a 'matrix quadratic'. A scalar quadratic in $z$ with coefficient $c$ has the form $cz^2$. If $\\mathbf{z}$ is a $k\\times 1$ vector and $\\mathbf{C}$ is a $k \\times k$ *matrix* of coefficients then the matrix quadratic form is written as $\\mathbf{z}^\\top \\mathbf{C}\\mathbf{z}$, which is itself a *scalar* quantity, but it is a function of a *vector*. \n",
    "\n",
    "#### Matching Dimensions in Matrix Multiplications\n",
    "\n",
    "There's a trick for telling that it's a scalar result. When you are doing maths with matrices, it's always worth pausing to perform a quick sanity check on the dimensions. Matrix multplication only works when the dimensions match. To be precise, the 'inner' dimension of the matrix must match. What is the inner dimension. If we multiply two matrices $\\mathbf{A}$ and $\\mathbf{B}$, the first of which has $k$ rows and $\\ell$ columns and the second of which has $p$ rows and $q$ columns, then we can check whether the multiplication works by writing the dimensionalities next to each other,\n",
    "$$\n",
    "\\mathbf{A} \\mathbf{B} \\rightarrow (k \\times \\underbrace{\\ell)(p}_\\text{inner dimensions} \\times q) \\rightarrow (k\\times q).\n",
    "$$\n",
    "The inner dimensions are the two inside dimensions, $\\ell$ and $p$. The multiplication will only work if $\\ell=p$. The result of the multiplication will then be a $k\\times q$ matrix: this dimensionality comes from the 'outer dimensions'. Note that matrix multiplication is not [*commutative*](http://en.wikipedia.org/wiki/Commutative_property). And if you change the order of the multiplication, \n",
    "$$\n",
    "\\mathbf{B} \\mathbf{A} \\rightarrow (\\ell \\times \\underbrace{k)(q}_\\text{inner dimensions} \\times p) \\rightarrow (\\ell \\times p).\n",
    "$$\n",
    "firstly it may no longer even work, because now the condition is that $k=q$, and secondly the result could be of a different dimensionality. An exception is if the matrices are square matrices (e.g. same number of rows as columns) and they are both *symmetric*. A symmetric matrix is one for which $\\mathbf{A}=\\mathbf{A}^\\top$, or equivalently, $a_{i,j} = a_{j,i}$ for all $i$ and $j$.  \n",
    "\n",
    "You will need to get used to working with matrices and vectors applying and developing new machine learning techniques. You should have come across them before, but you may not have used them as extensively as we will now do in this course. You should get used to using this trick to check your work and ensure you know what the dimension of an output matrix should be. For our matrix quadratic form, it turns out that we can see it as a special type of inner product.\n",
    "$$\n",
    "\\mathbf{z}^\\top\\mathbf{C}\\mathbf{z} \\rightarrow (1\\times \\underbrace{k) (k}_\\text{inner dimensions}\\times k) (k\\times 1) \\rightarrow \\mathbf{b}^\\top\\mathbf{z}\n",
    "$$\n",
    "where $\\mathbf{b} = \\mathbf{C}\\mathbf{z}$ so therefore the result is a scalar,\n",
    "$$\n",
    "\\mathbf{b}^\\top\\mathbf{z} \\rightarrow (1\\times \\underbrace{k) (k}_\\text{inner dimensions}\\times 1) \\rightarrow (1\\times 1)\n",
    "$$\n",
    "where a $(1\\times 1)$ matrix is recognised as a scalar.\n",
    "\n",
    "This implies that we should be able to differentiate this form, and indeed the rule for its differentiation is slightly more complex than the inner product, but still quite simple,\n",
    "$$\n",
    "\\frac{\\text{d}}{\\text{d}\\mathbf{z}} \\mathbf{z}^\\top\\mathbf{C}\\mathbf{z}= \\mathbf{C}\\mathbf{z} + \\mathbf{C}^\\top \\mathbf{z}.\n",
    "$$\n",
    "Note that in the special case where $\\mathbf{C}$ is symmetric then we have $\\mathbf{C} = \\mathbf{C}^\\top$ and the derivative simplifies to \n",
    "$$\n",
    "\\frac{\\text{d}}{\\text{d}\\mathbf{z}} \\mathbf{z}^\\top\\mathbf{C}\\mathbf{z}= 2\\mathbf{C}\\mathbf{z}.\n",
    "$$\n",
    "### Differentiating the Objective\n",
    "\n",
    "First, we need to compute the full objective by substituting our prediction function into the objective function to obtain the objective in terms of $\\mathbf{w}$. Doing this we obtain\n",
    "$$\n",
    "E(\\mathbf{w})= (\\mathbf{y} - \\mathbf{X}\\mathbf{w})^\\top (\\mathbf{y} - \\mathbf{X}\\mathbf{w}).\n",
    "$$\n",
    "We now need to differentiate this *quadratic form* to find the minimum. We differentiate with respect to the *vector* $\\mathbf{w}$. But before we do that, we'll expand the brackets in the quadratic form to obtain a series of scalar terms. The rules for bracket expansion across the vectors are similar to those for the scalar system giving,\n",
    "$$\n",
    "(\\mathbf{a} - \\mathbf{b})^\\top (\\mathbf{c} - \\mathbf{d}) = \\mathbf{a}^\\top \\mathbf{c} - \\mathbf{a}^\\top \\mathbf{d} - \\mathbf{b}^\\top \\mathbf{c} + \\mathbf{b}^\\top \\mathbf{d}\n",
    "$$\n",
    "which substituting for $\\mathbf{a} = \\mathbf{c} = \\mathbf{y}$ and $\\mathbf{b}=\\mathbf{d} = \\mathbf{X}\\mathbf{w}$ gives\n",
    "$$\n",
    "E(\\mathbf{w})= \\mathbf{y}^\\top\\mathbf{y} - 2\\mathbf{y}^\\top\\mathbf{X}\\mathbf{w} + \\mathbf{w}^\\top\\mathbf{X}^\\top\\mathbf{X}\\mathbf{w}\n",
    "$$\n",
    "where we used the fact that $\\mathbf{y}^\\top\\mathbf{X}\\mathbf{w}= \\mathbf{w}^\\top\\mathbf{X}^\\top\\mathbf{y}$. Now we can use our rules of differentiation to compute the derivative of this form, which is,\n",
    "$$\n",
    "\\frac{\\text{d}}{\\text{d}\\mathbf{w}}E(\\mathbf{w})=- 2\\mathbf{X}^\\top \\mathbf{y} + 2\\mathbf{X}^\\top\\mathbf{X}\\mathbf{w},\n",
    "$$\n",
    "where we have exploited the fact that $\\mathbf{X}^\\top\\mathbf{X}$ is symmetric to obtain this result."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Question 5\n",
    "\n",
    "Use the equivalence between our vector and our matrix formulations of linear regression, alongside our definition of vector derivates, to match the gradients we've computed directly for $\\frac{\\text{d}E(c, m)}{\\text{d}c}$ and $\\frac{\\text{d}E(c, m)}{\\text{d}m}$ to those for $\\frac{\\text{d}E(\\mathbf{w})}{\\text{d}\\mathbf{w}}$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Question 5 Answer\n",
    "\n",
    "Write your answer to the question in this box."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Update Equation for Global Optimum\n",
    "\n",
    "Once again, we need to find the minimum of our objective function. Using our likelihood for multiple input regression we can now minimize for our parameter vector $\\mathbf{w}$. Firstly, just as in the single input case, we seek stationary points by finding parameter vectors that solve for when the gradients are zero,\n",
    "$$\n",
    "\\mathbf{0}=- 2\\mathbf{X}^\\top \\mathbf{y} + 2\\mathbf{X}^\\top\\mathbf{X}\\mathbf{w},\n",
    "$$\n",
    "where $\\mathbf{0}$ is a *vector* of zeros. Rearranging this equation we find the solution to be\n",
    "$$\n",
    "\\mathbf{w} = \\left[\\mathbf{X}^\\top \\mathbf{X}\\right]^{-1} \\mathbf{X}^\\top \\mathbf{y}\n",
    "$$ \n",
    "where $\\mathbf{A}^{-1}$ denotes [*matrix inverse*](http://en.wikipedia.org/wiki/Invertible_matrix).\n",
    "\n",
    "### Solving the Multivariate System\n",
    "\n",
    "The solution for $\\mathbf{w}$ is given in terms of a matrix inverse, but computation of a matrix inverse requires, in itself, an algorithm to resolve it. You'll know this if you had to invert, by hand, a $3\\times 3$ matrix in high school. From a numerical stability perspective, it is also best not to compute the matrix inverse directly, but rather to ask the computer to *solve* the  system of linear equations given by\n",
    "$$\\mathbf{X}^\\top\\mathbf{X} \\mathbf{w} = \\mathbf{X}^\\top\\mathbf{y}$$\n",
    "for $\\mathbf{w}$. This can be done in `numpy` using the command"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.linalg.solve?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "so we can obtain the solution using"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "w = np.linalg.solve(np.dot(X.T, X), np.dot(X.T, y))\n",
    "print(w)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can map it back to the liner regression and plot the fit as follows"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "m = w[1]; c=w[0]\n",
    "f_test = m*x_test + c\n",
    "print(m)\n",
    "print(c)\n",
    "plt.plot(x_test, f_test, 'b-')\n",
    "plt.plot(x, y, 'rx')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Multivariate Linear Regression\n",
    "\n",
    "A major advantage of the new system is that we can build a linear regression on a multivariate system. The matrix calculus didn't specify what the length of the vector $\\mathbf{x}$ should be, or equivalently the size of the design matrix. \n",
    "\n",
    "### Movie Body Count Data\n",
    "\n",
    "Let's load back in the movie body count data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data = pods.datasets.movie_body_count()\n",
    "movies = data['Y']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's remind ourselves of the features we've been provided with."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(', '.join(movies.columns))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we will build a design matrix based on the numeric features: year, Body_Count, Length_Minutes in an effort to predict the rating. We build the design matrix as follows:\n",
    "\n",
    "## Relation to Single Input System\n",
    "\n",
    "Bias as an additional feature."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "select_features = ['Year', 'Body_Count', 'Length_Minutes']\n",
    "X = movies.loc[:, select_features]\n",
    "X['Eins'] = 1 # add a column for the offset\n",
    "y = movies[['IMDB_Rating']]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now let's perform a linear regression. But this time, we will create a pandas data frame for the result so we can store it in a form that we can visualise easily."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "w = pd.DataFrame(data=np.linalg.solve(np.dot(X.T, X), np.dot(X.T, y)),  # solve linear regression here\n",
    "                 index = X.columns,  # columns of X become rows of w\n",
    "                 columns=['regression_coefficient']) # the column of X is the value of regression coefficient\n",
    "w"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can check the residuals to see how good our estimates are"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "(y - np.dot(X, w)).hist()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Which shows our model *hasn't* yet done a great job of representation, because the spread of values is large. We can check what the rating is dominated by in terms of regression coefficients."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "w"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Although we have to be a little careful about interpretation because our input values live on different scales, however it looks like we are dominated by the bias, with a small negative effect for later films (but bear in mind the years are large, so this effect is probably larger than it looks) and a positive effect for length. So it looks like long earlier films generally do better, but the residuals are so high that we probably haven't modelled the system very well."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from IPython.display import YouTubeVideo\n",
    "YouTubeVideo('ui-uNlFHoms')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from IPython.display import YouTubeVideo\n",
    "YouTubeVideo('78YNphT90-k')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Solution with QR Decomposition\n",
    "\n",
    "Performing a solve instead of a matrix inverse is the more numerically stable approach, but we can do even better. A [QR-decomposition](http://en.wikipedia.org/wiki/QR_decomposition) of a matrix factorises it into a matrix which is an orthogonal matrix $\\mathbf{Q}$, so that $\\mathbf{Q}^\\top \\mathbf{Q} = \\mathbf{I}$. And a matrix which is upper triangular, $\\mathbf{R}$. \n",
    "$$\n",
    "\\mathbf{X}^\\top \\mathbf{X} \\boldsymbol{\\beta} = \\mathbf{X}^\\top \\mathbf{y}\n",
    "$$\n",
    "$$\n",
    "(\\mathbf{Q}\\mathbf{R})^\\top (\\mathbf{Q}\\mathbf{R})\\boldsymbol{\\beta} = (\\mathbf{Q}\\mathbf{R})^\\top \\mathbf{y}\n",
    "$$\n",
    "$$\n",
    "\\mathbf{R}^\\top (\\mathbf{Q}^\\top \\mathbf{Q}) \\mathbf{R} \\boldsymbol{\\beta} = \\mathbf{R}^\\top \\mathbf{Q}^\\top \\mathbf{y}\n",
    "$$\n",
    "$$\n",
    "\\mathbf{R}^\\top \\mathbf{R} \\boldsymbol{\\beta} = \\mathbf{R}^\\top \\mathbf{Q}^\\top \\mathbf{y}\n",
    "$$\n",
    "$$\n",
    "\\mathbf{R} \\boldsymbol{\\beta} = \\mathbf{Q}^\\top \\mathbf{y}\n",
    "$$\n",
    "This is a more numerically stable solution because it removes the need to compute $\\mathbf{X}^\\top\\mathbf{X}$ as an intermediate. Computing $\\mathbf{X}^\\top\\mathbf{X}$ is a bad idea because it involves squaring all the elements of $\\mathbf{X}$ and thereby potentially reducing the numerical precision with which we can represent the solution. Operating on $\\mathbf{X}$ directly preserves the numerical precision of the model.\n",
    "\n",
    "This can be more particularly seen when we begin to work with *basis functions* in the next week. Some systems that can be resolved with the QR decomposition can not be resolved by using solve directly."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import scipy as sp\n",
    "Q, R = np.linalg.qr(X)\n",
    "w = sp.linalg.solve_triangular(R, np.dot(Q.T, y)) \n",
    "w = pd.DataFrame(w, index=X.columns)\n",
    "w"
   ]
  }
 ],
 "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.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}