02b - SolvingSystemsOfLinearEquations.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "FI4gjuRVLvEK"
   },
   "source": [
    "Credit: some material here has been adapted from [Sam Roweis](https://www.cs.nyu.edu/home/people/in_memoriam/samroweis.html)' [Linear Algebra Review](http://www.cs.ubc.ca/~murphyk/Teaching/Papers/roweis_linAlgebra.ps).\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 0,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "3dyf5H0Gjwul"
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "EHJv4OiTLvEv"
   },
   "source": [
    "# Systems of Linear Equations\n",
    "\n",
    "A central problem in linear algebra is the solution of a system of linear equations like this:\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "3x + 4y + 2z &= 12\\\\\n",
    "x + y + z &=5\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "Typically, we express this system as a single *matrix equation* in the form: $A\\mathbf{x} = \\mathbf{b}$, where $A$ is an $m$-by-$n$ matrix, $\\mathbf{x}$ is an $n$-dimensional column vector, and $\\mathbf{b}$ is an $m$-dimensional column vector. The number of unknowns is $n$ and the number of equations or constraints is $m$. Here is another simple example:\n",
    "\n",
    "$$\n",
    "\\left[ \n",
    "\\begin{array}{rr}\n",
    "2 & -1\\\\\n",
    "1 & 1\n",
    "\\end{array}\n",
    "\\right] \n",
    "\\left[\n",
    "\\begin{array}{c}\n",
    "x_1 \\\\\n",
    "x_2\n",
    "\\end{array}\n",
    "\\right]\n",
    "=\n",
    "\\left[\n",
    "\\begin{array}{c}\n",
    "1\\\\\n",
    "5\n",
    "\\end{array}\n",
    "\\right]\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "dR-vTjLdLvEw"
   },
   "source": [
    "How do we go about \"solving\" this system of equations? \n",
    "\n",
    "* Finding $\\mathbf{b}$ given $A$ and $\\mathbf{x}$ is easy - just multiply\n",
    "* For a single $\\mathbf{x}$ and $\\mathbf{b}$ there are usually a great many matrices $A$ which satisfy the equation\n",
    "* The only interesting problem left is to find $\\mathbf{x}$ given $A$ and $\\mathbf{b}$!\n",
    "\n",
    "This kind of equation is really a problem statement. It says \"we applied the function $A$ and generated the output $\\mathbf{b}$; what was the input $\\mathbf{x}$?\"\n",
    "\n",
    "The matrix $A$ is dictated to us by our problem, and represents our model of how the system we are studying converts inputs to outputs. The vector $\\mathbf{b}$ is the output that we observe (or desire) - we know it. The vector $\\mathbf{x}$ is the set of inputs that we are trying to find.\n",
    "\n",
    "There are two ways of thinking about this system of equations. One is *rowwise* as a set of $m$ equations, or constraints that correspond geometrically to $m$ intersecting constraint surfaces:\n",
    "\n",
    "$$\n",
    "\\left[\n",
    "\\begin{array}{r}\n",
    "2x_1 - x_2 &= 1\\\\\n",
    "x_1 + x_2 &= 5\n",
    "\\end{array}\n",
    "\\right]\n",
    "$$\n",
    "\n",
    "The goal is to find the point(s), for example $(x_1, x_2)$ above, which are at the intersection of all the constraint surfaces. In the example above, these surfaces are two lines in a plane. If the lines intersect then there is a solution. If they are parallel, there is not. If they are coincident then there are infinite solutions. In higher dimensions there are more geometric solutions, but in general there can be no solutions, one solution, or infinite solutions.\n",
    "\n",
    "The other way of thinking about this system is *columnwise* in which we think of the entire system as a single vector relation:\n",
    "\n",
    "$$\n",
    "x_1\n",
    "\\left[\n",
    "\\begin{array}{r}\n",
    "2\\\\\n",
    "1\n",
    "\\end{array}\n",
    "\\right]\n",
    "+\n",
    "x_2\n",
    "\\left[\n",
    "\\begin{array}{r}\n",
    "-1\\\\\n",
    "1\n",
    "\\end{array}\n",
    "\\right]\n",
    "=\n",
    "\\left[\n",
    "\\begin{array}{r}\n",
    "1\\\\\n",
    "5\n",
    "\\end{array}\n",
    "\\right]\n",
    "$$\n",
    "\n",
    "The goal here is to discover which linear combination(s) $(x_1, x_2)$, if any, of the $n$ column vectors on the left will give the one on the right.\n",
    "\n",
    "Either way, the matrix equation $A\\mathbf{x} = \\mathbf{b}$ is an almost ubiquitous problem whose solution comes up again and again in theoretical derivations and in practical data analysis problems. One of the most important applications is the minimization of quadratic energy functions: if $A$ is symmetric positive definite then the quadratic form $\\mathbf{x}^TA\\mathbf{x} - 2 \\mathbf{x}^T\\mathbf{b} + c$ is *minimized* at the point where $A\\mathbf{x}=\\mathbf{b}$. Such quadratic forms arise often in the study of linear models with Gaussian noise since the log likelihood of data under such models is always a matrix quadratic.\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "za5ON7Hm8R7P"
   },
   "source": [
    "## Matrix inversion and Cramer's rule\n",
    "\n",
    "Let's first consider the case of a single $\\mathbf{x}$. As noted above, geometrically we can think of the rows of the system as encoding constraint surfaces in which the solution vector $\\mathbf{x}$ must lie and what we are looking for is the point(s) at which these surfaces intersect. Of course, there may be no solution satisfying the equation; then we typically need some way to pick the \"best\" approximate solution. The constraints may also intersect along an entire line or surface in which case there are an infinity of solutions; once again we would like to think about which one might be best.\n",
    "\n",
    "Let's consider finding exact solutions first. The most naive thing we could do is just find the inverse of $A$ and solve the equations as follows:\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "A^{-1}A\\mathbf{x} & = A^{-1} \\mathbf{b}\\\\\n",
    "I\\mathbf{x} &= A^{-1}\\mathbf{b}\\\\\n",
    "\\mathbf{x} &= A^{-1}\\mathbf{b}\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "which is known as *Cramer's rule*.\n",
    "\n",
    "Let's apply Cramer's rule to the system of equations above:\n",
    "\n",
    "$$\n",
    "\\left[\n",
    "\\begin{array}{r}\n",
    "2x_1 - x_2 &= 1\\\\\n",
    "x_1 + x_2 &= 5\n",
    "\\end{array}\n",
    "\\right]\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 0,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "SI09BXSY8rXE"
   },
   "outputs": [],
   "source": [
    "A = np.array([[2, -1], [1, 1]])\n",
    "print(A)\n",
    "b = np.array([1, 5])\n",
    "print(b)\n",
    "x = np.linalg.inv(A).dot(b) #Cramer's rule, slow\n",
    "print(x)\n",
    "print(A.dot(np.linalg.inv(A).dot(b))-b) #check"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "wyaYq5Wp9BPP"
   },
   "source": [
    "There are several problems with this approach. Most importantly, many interesting functions are not invertible. Beyond that, matrix inversion is a difficult and potentially numerically unstable operation. **Don't invert a matrix unless you really have to!**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "OYClkK7-6TVM"
   },
   "source": [
    "## Exact solutions: ``np.linalg.solve``\n",
    "In cases where the matrix $A$ is square and of full-rank, one can use ``numpy.linalg.solve``. \n",
    "\n",
    "Let's return to the same example and solve it more efficiently."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 0,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "14yTX1Nk7AYt"
   },
   "outputs": [],
   "source": [
    "A = np.array([[2, -1], [1, 1]])\n",
    "b = np.array([1, 5])\n",
    "x = np.linalg.solve(A, b)\n",
    "print(x)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "jldnLrrc9peA"
   },
   "source": [
    "Note, though, that this exact method fails when $A$ is not square."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 0,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "mz89_xug9mHZ"
   },
   "outputs": [],
   "source": [
    "A = np.array([[1, 2, 3],[4, 5, 6]])\n",
    "print(A)\n",
    "b = np.array([[5],[6]])\n",
    "print(b)\n",
    "x = np.linalg.solve(A, b)\n",
    "print(x)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "20l4NBzU6QK7"
   },
   "source": [
    "## Least squares:  ``np.linalg.lstsq``\n",
    "\n",
    "In fact, there is a much better way to define what we want as a solution. We will say that we want a solution $\\mathbf{x}^*$ which minimizes the error:\n",
    "\n",
    "$$E = ||A\\mathbf{x}^* - \\mathbf{b}||^2 = \\mathbf{x}^TA^TA\\mathbf{x} - 2\\mathbf{x}^TA^T\\mathbf{b} + \\mathbf{b}^T\\mathbf{b}$$\n",
    "\n",
    "This problem is known as *linear least squares*, for obvious reasons. If there is an exact solution (one giving zero error) it will certainly minimize the above cost, but if there is no solution, we can still find the best possible approximation. If we take the matrix derivative of this expression, we can find the best solution:\n",
    "\n",
    "$$\\mathbf{x}^* = (A^TA)^{-1}A^T\\mathbf{b}$$\n",
    "\n",
    "which takes advantage of the fact that even if $A$ is not invertible, $A^TA$ may be.\n",
    "\n",
    "But what if the problem is degenerate? In other words, what if there is more than one exact solution or more than one inexact solution which all achieve the same minimum error. How can this occur? \n",
    "\n",
    "Imagine an equation like this:\n",
    "\n",
    "$$\\left[\n",
    "\\begin{array}{rr}\n",
    "1 &  -1\n",
    "\\end{array}\n",
    "\\right] \\mathbf{x} = 4$$\n",
    "\n",
    "in which $A = \\left[ \\begin{array}{rr}\n",
    "1 &  -1 \n",
    "\\end{array} \\right]$. This equation constrains the difference between the two elements of $\\mathbf{x}$ to be 4 but the sum of the two elements can be as large or small as we want. \n",
    "\n",
    "We can take things one step further to get around this problem. The answer is to ask for the *minimum norm* vector $\\mathbf{x}$ that still minimizes the above error. This breaks the degeneracies in both the exact and inexact cases. In terms of the cost function, this corresponds to adding an infinitesimal penalty on $\\mathbf{x}^T\\mathbf{x}$:\n",
    "\n",
    "$$\n",
    "E = \\text{lim}_{\\epsilon \\rightarrow 0} \\left[ \\mathbf{x}^TA^TA\\mathbf{x} - 2\\mathbf{x}^TA^T\\mathbf{b} + \\mathbf{b}^T\\mathbf{b} + \\epsilon \\mathbf{x}^T\\mathbf{x} \\right]\n",
    "$$\n",
    "\n",
    "and the optimal solution becomes\n",
    "\n",
    "$$\n",
    "\\mathbf{x}^* = \\lim_{\\epsilon \\rightarrow 0}  \\left[ \\left(A^TA + \\epsilon I\\right)^{-1}A^T\\mathbf{b}\\right]\n",
    "$$\n",
    "\n",
    "Now, of course actually computing these solutions efficiently and in a numerically stable way is the topic of much study in numerical methods. However, in Python you don't have to worry about any of this, you can just type ``np.linalg.lstsq(A, b)`` and rely on the NumPy maintainers to deliver an optimized solution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 0,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "CN4HfcdJLvE0"
   },
   "outputs": [],
   "source": [
    "A = np.array([[1, 2, 3],[4, 5, 6]])\n",
    "print(A)\n",
    "b = np.array([[5],[6]])\n",
    "print(b)\n",
    "# computes the least-squares solution to a linear matrix equation\n",
    "# equation may be under-, well, or over- determined\n",
    "x, residuals, rank, s =  np.linalg.lstsq(A, b, rcond=None)\n",
    "print(x)\n",
    "print(A.dot(x))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "qsbiG_u8_URV"
   },
   "source": [
    "## Reduced-row echelon form\n",
    "\n",
    "You will notice that we didn't use Gaussian elimination to reduce a matrix to reduced-row echelon form. In fact, numerical linear algebra is very different than the linear algebra we see in textbooks. Reduced-row echelon form doesn't really work when using floating-point values. NumPy, as such, does not have an implementation of Gaussian elimination.\n",
    "\n",
    "To see a Python implementation, you can head to the SymPy library which does exact computations. It has a good [tutorial on matrix manipulation](https://docs.sympy.org/latest/tutorial/matrices.html)."
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "colab": {
   "collapsed_sections": [],
   "name": "SolvingSystemsOfLinearEquations.ipynb",
   "provenance": [
    {
     "file_id": "1iExflxBw3D88L2PLjKguNvAJnise2JgZ",
     "timestamp": 1564338887194
    }
   ]
  },
  "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.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}