Notebook12_Algoriths_LinearLeastSq_part3.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Part 2: Algorithms for the linear least squares problem\n",
    "\n",
    "Recall the linear regression problem: given a data matrix, $X$, and responses $y$, we wish to determine the model parameters $\\theta^*$ that minimizes $\\|X \\theta - y\\|_2^2$. This problem is also known as the _linear least squares_ problem.\n",
    "\n",
    "Numpy has a function, [np.linalg.lstsq()](https://docs.scipy.org/doc/numpy/reference/generated/numpy.linalg.lstsq.html), that will compute linear least squares solutions for you. However, the purpose of this notebook is to give you a sense of how `lstsq()` works. So, instead of using it as a black box, you'll study how one might implement `lstsq()` using two different numerical algorithms.\n",
    "\n",
    "> You may rightly ask, why bother with such details? Here are three reasons it's worth looking more closely.\n",
    ">\n",
    "> 1. It's helpful to have some deeper intuition for how one formalizes a mathematical problem and derives a computational solution, in case you ever encounter a problem that does not exactly fit what a canned library can do for you.\n",
    ">\n",
    "> 2. If you have ever used a statistical analysis package, it's likely you have encountered \"strange\" numerical errors or warnings. Knowing how problems are derived can help you understand what might have gone wrong. We will see an example below.\n",
    ">\n",
    "> 3. Because data analysis is quickly evolving, it's likely that new problems and new models will not exactly fit the template of existing models. Therefore, it's possible you will need to derive a new model or know how to talk to someone who can derive one for you."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Implementation note.** In this notebook, we ask that you use the following convention: any column vector should be _explicit_. That means its shape should have two dimensions where the column dimension equals one (1).\n",
    "\n",
    "**Exercise 0** (ungraded). Inspect the following code cell and make sure you understand the difference between two conventions for storing a vector, namely, as a one-dimensional array versus as a two-dimensional array (matrix) where the number of columns equals one (1). When you are asked to produce a vector, we will generally ask you to follow the second convention (`z_colvec`)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "`z_array`:\n",
      "\n",
      " [ 1.  2.  3.] \n",
      "\n",
      "==> shape: (3,)\n",
      "\n",
      "\n",
      "`z_colvec`:\n",
      "\n",
      " [[ 1.]\n",
      " [ 2.]\n",
      " [ 3.]] \n",
      "\n",
      "==> shape: (3, 1)\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "\n",
    "# By default, Numpy vectors constructed from a list are 1-D\n",
    "# and have neither a row nor a column orientation.\n",
    "z_array = np.array([1.0, 2.0, 3.0])\n",
    "\n",
    "# By contrast, we want you to ensure your vectors are\n",
    "# column vectors.\n",
    "z_colvec = np.reshape(z_array, (len(z_array), 1))\n",
    "\n",
    "print(\"`z_array`:\\n\\n\", z_array, \"\\n\\n==> shape:\", z_array.shape)\n",
    "print(\"\\n\")\n",
    "print(\"`z_colvec`:\\n\\n\", z_colvec, \"\\n\\n==> shape:\", z_colvec.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Before beginning, run this code cell to load some of the key modules you'll need."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true,
    "nbgrader": {
     "grade": false,
     "grade_id": "cell-2b1c86b5bcd7388d",
     "locked": true,
     "schema_version": 1,
     "solution": false
    }
   },
   "outputs": [],
   "source": [
    "# Data and computation\n",
    "import numpy as np\n",
    "import scipy as sp\n",
    "import scipy.linalg\n",
    "import pandas as pd\n",
    "\n",
    "# Viz\n",
    "from IPython.display import display, Math\n",
    "from matplotlib.pyplot import figure, subplot, xlim, ylim\n",
    "from matplotlib.pyplot import scatter, axis, xlabel, ylabel, title, plot\n",
    "%matplotlib inline\n",
    "\n",
    "# Some functions we'll use later to display results\n",
    "def show_cond_fancy(x, name, opt=''):\n",
    "    \"\"\"Display a condition number in 'fancy' format (using LaTeX).\"\"\"\n",
    "    def sci_to_latex(x, fmt='{:.2e}'):\n",
    "        s_raw = fmt.format(x)\n",
    "        s, e = s_raw.split('e')\n",
    "        return s + r'\\times 10^{{{}}}'.format(int(e))\n",
    "    from IPython.display import Math\n",
    "    x_s = sci_to_latex(x)\n",
    "    display(Math(r'\\kappa({}){} \\approx {}'.format(name, opt, x_s)))\n",
    "    \n",
    "def show_2vecs_tibble(x, y, xname='x', yname='y', error=False):\n",
    "    \"\"\"Display two column vectors side-by-side in a tibble.\"\"\"\n",
    "    assert type(x) is np.ndarray and x.ndim >= 2 and x.shape[1] == 1\n",
    "    assert type(y) is np.ndarray and y.ndim >= 2 and y.shape[1] == 1\n",
    "    assert x.shape == y.shape\n",
    "    x_df = pd.DataFrame(x, columns=[xname])\n",
    "    y_df = pd.DataFrame(y, columns=[yname])\n",
    "    df = pd.concat([x_df, y_df], axis=1)\n",
    "    if error:\n",
    "        df['error'] = x - y\n",
    "    display(df)\n",
    "    \n",
    "# Display (X, y) problem as a tibble\n",
    "def make_data_tibble(X, y=None):\n",
    "    df = pd.DataFrame(X, columns=['x_{}'.format(i) for i in range(X.shape[1])])\n",
    "    if y is not None:\n",
    "        y_df = pd.DataFrame(y, columns=['y'])\n",
    "        df = pd.concat([y_df, df], axis=1)\n",
    "    return df\n",
    "    \n",
    "# From: https://stackoverflow.com/questions/17129290/numpy-2d-and-1d-array-to-latex-bmatrix\n",
    "def nparray_to_bmatrix(a):\n",
    "    \"\"\"Returns a LaTeX bmatrix\"\"\"\n",
    "    assert len(a.shape) <= 2, 'bmatrix can at most display two dimensions'\n",
    "    lines = str(a).replace('[', '').replace(']', '').splitlines()\n",
    "    rv = [r'\\begin{bmatrix}']\n",
    "    rv += ['  ' + ' & '.join(l.split()) + r'\\\\' for l in lines]\n",
    "    rv +=  [r'\\end{bmatrix}']\n",
    "    return '\\n'.join(rv)\n",
    "\n",
    "# Stash this function for later:\n",
    "SAVE_LSTSQ = np.linalg.lstsq # You may ignore this line, which some test cells will use"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Notation and review\n",
    "\n",
    "Here is a quick summary of how we can formulate and approach the linear regression problem. For a more detailed derivation, see these [accompanying notes](./notes-linreg.ipynb).\n",
    "\n",
    "Your data consists of $m$ observations and $n+1$ variables. One of these variables is the _response_ variable, $y$, which you want to predict from the other $n$ variables, $\\{x_0, \\ldots, x_{n-1}\\}$. You wish to fit a _linear model_ of the following form to these data,\n",
    "\n",
    "$$y_i \\approx x_{i,0} \\theta_0 + x_{i,1} \\theta_1 + \\cdots + x_{i,n-1} \\theta_{n-1} + \\theta_n,$$\n",
    "\n",
    "where $\\{\\theta_j | 0 \\leq j \\leq n\\}$ is the set of unknown coefficients. Your modeling task is to choose values for these coefficients that \"best fit\" the data.\n",
    "\n",
    "If we further define a set of dummy variables, $x_{i, n} \\equiv 1.0$, associated with the $\\theta_n$ parameter, then the model can be written more compactly in matrix notation as\n",
    "\n",
    "$$\n",
    "  y \\approx X \\theta,\n",
    "$$\n",
    "\n",
    "where we will refer to $X$ as the (input) data matrix."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Visually, you can also arrange the observations into a tibble like this one:\n",
    "\n",
    "|     y      | x<sub>0</sub> | x<sub>1</sub> | $\\cdots$ | x<sub>n-1</sub> | x<sub>n</sub> |\n",
    "|:----------:|:-------------:|:-------------:|:--------:|:---------------:|:-------------:|\n",
    "|   $y_0$    |   $x_{0,1}$   |   $x_{0,2}$   | $\\cdots$ |   $x_{0,n-1}$   |      1.0      |\n",
    "|   $y_1$    |   $x_{1,1}$   |   $x_{1,2}$   | $\\cdots$ |   $x_{1,n-1}$   |      1.0      |\n",
    "|   $y_2$    |   $x_{2,1}$   |   $x_{2,2}$   | $\\cdots$ |   $x_{2,n-1}$   |      1.0      |\n",
    "|  $\\vdots$  |   $\\vdots$    |   $\\vdots$    | $\\vdots$ |    $\\vdots$     |      1.0      |\n",
    "|  $y_{m-1}$ |  $x_{m-1,1}$  |  $x_{m-1,2}$  | $\\cdots$ |  $x_{m-1,n-1}$  |      1.0      |\n",
    "\n",
    "This tibble includes an extra column (variable), $x_n$, whose entries are all equal to 1.0."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Synthetic problem generator.** For the exercises in this notebook, we will generate synthetic data. The function, `gen_problem(m, n)`, will return a triple, `X`, `y`, `theta`, which are an `m` x (`n`+1) data matrix `X`, a response vector `y`, and the \"true\" model parameters `theta`. We will then run two different numerical algorithms that estimate `theta` from `X` and `y`, and see how their answers compare against the true value.\n",
    "\n",
    "> **Note 1.** The problem generator constructs the data matrix $X$ such that each entry $(i, j)$ is $i^{\\,j}$. This structure makes it an instance of a [Vandermonde matrix](https://en.wikipedia.org/wiki/Vandermonde_matrix), which arises when fitting a polynomial to data. The \"true\" parameter vector $\\theta$ is set to all ones, and $y$ computed simply by summing the rows.\n",
    ">\n",
    "> **Note 2.** Although our usual convention is to make the _last_ column all ones, the Vandermonde matrix has its _first_ column set to all ones. This ordering is not important in this problem, but it does mean one would interpret $\\theta_0$ as the intercept rather than $\\theta_{n}$, which will be our usual convention."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Sample generated problem:\n"
     ]
    },
    {
     "data": {
      "text/latex": [
       "$$X = \\begin{bmatrix}\n",
       "  1. & 0. & 0.\\\\\n",
       "  1. & 1. & 1.\\\\\n",
       "  1. & 2. & 4.\\\\\n",
       "  1. & 3. & 9.\\\\\n",
       "  1. & 4. & 16.\\\\\n",
       "  1. & 5. & 25.\\\\\n",
       "  1. & 6. & 36.\\\\\n",
       "  1. & 7. & 49.\\\\\n",
       "  1. & 8. & 64.\\\\\n",
       "  1. & 9. & 81.\\\\\n",
       "\\end{bmatrix}, \\quad y = \\begin{bmatrix}\n",
       "  1.\\\\\n",
       "  3.\\\\\n",
       "  7.\\\\\n",
       "  13.\\\\\n",
       "  21.\\\\\n",
       "  31.\\\\\n",
       "  43.\\\\\n",
       "  57.\\\\\n",
       "  73.\\\\\n",
       "  91.\\\\\n",
       "\\end{bmatrix} \\quad \\implies \\quad \\theta^* = \\begin{bmatrix}\n",
       "  1.\\\\\n",
       "  1.\\\\\n",
       "  1.\\\\\n",
       "\\end{bmatrix}$$"
      ],
      "text/plain": [
       "<IPython.core.display.Math object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "def gen_problem(m, n):\n",
    "    from numpy import arange, tile, cumprod, insert, ones\n",
    "    # 1 + x + x^2 + ... + x^n, x = 0:m\n",
    "    X = np.empty((m, n+1))\n",
    "    x_col = arange(m).reshape((m, 1)) # 0, 1, 2, ..., m-1\n",
    "    X[:, 0] = 1.0\n",
    "    X[:, 1:] = tile(x_col, reps=(1, n))\n",
    "    X[:, 1:] = cumprod(X[:, 1:], axis=1)\n",
    "    theta = ones((n+1, 1))\n",
    "    y = np.sum(X, axis=1).reshape((m, 1))\n",
    "    return X, y, theta\n",
    "\n",
    "print(\"Sample generated problem:\")\n",
    "m, n = 10, 2\n",
    "X, y, theta_true = gen_problem(m, n)\n",
    "\n",
    "display(Math(r'X = {}, \\quad y = {} \\quad \\implies \\quad \\theta^* = {}'.format(nparray_to_bmatrix(X),\n",
    "                                                                               nparray_to_bmatrix(y),\n",
    "                                                                               nparray_to_bmatrix(theta_true))))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We are interested primarily in _overdetermined systems_, meaning $X$ has more rows than columns, i.e., $m > n+1$, as shown above. That's because typically we have more observations (data points, or rows) than predictors (variables or columns). For such problems, there is generally no unique solution.\n",
    "\n",
    "Therefore, to identify some solution, we need to ask for the \"best\" fit and say what we mean by \"best.\" For linear regression, the usual definition of \"best\" is _minimizing_ the sum-of-squared residual error:\n",
    "\n",
    "$$ \\theta^* = {\\arg\\min_\\theta} \\|X \\theta - y\\|_2^2. $$\n",
    "\n",
    "Solving this minimization problem is equivalent to solving a special system known as the _normal equations_,\n",
    "\n",
    "$$\n",
    "  X^T X \\theta^* = X^T y.\n",
    "$$\n",
    "\n",
    "So, our computational task is to solve this problem."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Algorithm 1: Direct solution of the normal equations\n",
    "\n",
    "The preceding calculation immediately suggests the following algorithm to estimate $\\theta^*$. Given $X$ and $y$:\n",
    "\n",
    "1. Form $C \\equiv X^T X$. This object is sometimes called the [Gram matrix](https://en.wikipedia.org/wiki/Gramian_matrix) or Gramian of $X$.\n",
    "2. Form $b \\equiv X^T y$.\n",
    "3. Solve $C \\theta^* = b$ for $\\theta^*$.\n",
    "\n",
    "But, is this a \"good\" algorithm? There are at least three dimensions along which we might answer this question.\n",
    "\n",
    "1. Is it accurate enough?\n",
    "2. Is it fast enough?\n",
    "3. Is it memory-efficient enough?\n",
    "\n",
    "Let's examine these questions by experiment."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise 1** (3 points). Implement a function, `solve_neq(X, y)` that implements Algorithm 1. It should return a Numpy vector containing the model parameter estimates.\n",
    "\n",
    "Recall the steps of the algorithm as previously outlined:\n",
    "\n",
    "1. Form the Gramian of $X$, $C \\equiv X^T X$.\n",
    "2. Form $b \\equiv X^T y$.\n",
    "3. Solve $C \\theta^* = b$ for $\\theta^*$.\n",
    "\n",
    "Your algorithm should carry out these steps. For the third step, use Scipy's routine, [`scipy.linalg.solve()`](https://docs.scipy.org/doc/scipy/reference/tutorial/linalg.html#solving-linear-system). It has an option that allows you to indicate that $C$ is symmetric and positive definite, which will be true of $C$ for our synthetic problem.\n",
    "\n",
    "> The code cell will run your function to compute a set of parameter estimates. It will store these in a variable named `theta_neq`, which we will refer to later."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "nbgrader": {
     "grade": false,
     "grade_id": "solve_neq",
     "locked": false,
     "schema_version": 1,
     "solution": true
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Your implementation's solution versus the true solution:\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style>\n",
       "    .dataframe thead tr:only-child th {\n",
       "        text-align: right;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: left;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>theta_neq</th>\n",
       "      <th>theta_true</th>\n",
       "      <th>error</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>-2.065015e-14</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>4.218847e-15</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>0.000000e+00</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   theta_neq  theta_true         error\n",
       "0        1.0         1.0 -2.065015e-14\n",
       "1        1.0         1.0  4.218847e-15\n",
       "2        1.0         1.0  0.000000e+00"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "def solve_neq(X, y):\n",
    "    #\n",
    "    # YOUR CODE HERE\n",
    "    C = X.T.dot(X)\n",
    "    b = X.T.dot(y)\n",
    "    theta_star = scipy.linalg.solve(C, b)\n",
    "    return theta_star\n",
    "\n",
    "theta_neq = solve_neq(X, y)\n",
    "\n",
    "print(\"Your implementation's solution versus the true solution:\")\n",
    "show_2vecs_tibble(theta_neq, theta_true, xname='theta_neq', yname='theta_true', error=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "nbgrader": {
     "grade": true,
     "grade_id": "solve_neq_test",
     "locked": true,
     "points": 3,
     "schema_version": 1,
     "solution": false
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "(Passed!)\n"
     ]
    }
   ],
   "source": [
    "# Test cell: `solve_neq_test`\n",
    "\n",
    "try:\n",
    "    del np.linalg.lstsq\n",
    "    solve_neq(X, y)\n",
    "except NameError as n:\n",
    "    if re.findall('lstsq', n.args[0]):\n",
    "        print(\"*** Double-check that you did not try to use `lstsq()`. ***\")\n",
    "    raise n\n",
    "except AttributeError as a:\n",
    "    if re.findall('lstsq', a.args[0]):\n",
    "        print(\"*** Double-check that you did not try to use `lstsq()`. ***\")\n",
    "    raise a\n",
    "finally:\n",
    "    np.linalg.lstsq = SAVE_LSTSQ\n",
    "\n",
    "assert type(theta_neq) is np.ndarray, \"`theta_neq` should be a Numpy array, but isn't.\"\n",
    "assert theta_neq.shape == (n+1, 1), \"`theta_neq.shape` is {} instead of {}.\".format(theta_neq.shape, (n+1, 1))\n",
    "\n",
    "assert (np.abs(theta_neq - theta_true) <= 1e-12).all(), \\\n",
    "       \"Your `theta_neq` does not match the true solution, `theta_true`.\"\n",
    "\n",
    "print(\"\\n(Passed!)\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise 2** (1 point). Write a function to calculate the residual norm, $\\|r\\|_2 = \\|X \\theta^* - y\\|_2$.\n",
    "\n",
    "> Although we are minimizing $\\|r\\|_2^2$, for this exercise your function should return $\\|r\\|_2$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "nbgrader": {
     "grade": false,
     "grade_id": "calc_residual_norm",
     "locked": false,
     "schema_version": 1,
     "solution": true
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "The squared residual norm: [ -1.75415238e-14]\n"
     ]
    }
   ],
   "source": [
    "def calc_residual_norm(X, y, theta):\n",
    "    #\n",
    "    # YOUR CODE HERE\n",
    "    res_norm = sum(X.dot(theta) - y)\n",
    "    return res_norm\n",
    "\n",
    "r_norm_neq = calc_residual_norm(X, y, theta_neq)\n",
    "print(\"\\nThe squared residual norm:\", r_norm_neq)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "nbgrader": {
     "grade": true,
     "grade_id": "calc_residual_norm_test",
     "locked": true,
     "points": 1,
     "schema_version": 1,
     "solution": false
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "(Passed.)\n"
     ]
    }
   ],
   "source": [
    "# Test cell: `calc_residual_norm_test`\n",
    "\n",
    "r_norm_neq = calc_residual_norm(X, y, theta_neq)\n",
    "assert 1e-16 <= np.abs(r_norm_neq) <= 1e-12\n",
    "print (\"\\n(Passed.)\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Sources of error\n",
    "\n",
    "We said before that one question we should ask about our algorithm is whether it is \"accurate enough.\" But what does that mean?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise 3** (ungraded). For any modeling problem, there will be several sources of error. Describe at least three such sources."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "nbgrader": {
     "grade": false,
     "grade_id": "error_sources",
     "locked": true,
     "schema_version": 1,
     "solution": false
    }
   },
   "source": [
    "**Answer.** Here are some possibilities.\n",
    "\n",
    "1. There will be errors in the inputs. That is, the data itself may only represent measurements of a certain accuracy.\n",
    "2. There will be errors in the model. That is, the model is only an approximation of the underlying phenomena.\n",
    "3. There will be errors in the algorithm. That is, you may implement an algorithm that can only approximately estimate the parameters of the model.\n",
    "4. There will be roundoff errors. Recall that floating-point arithmetic necessarily represents all values finitely, which means you may lose accuracy every time you do an arithmetic operation."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Perturbations.** One way to understand error in a numerical computation is to consider how sensitive the computed solution is to perturbations to the input.\n",
    "\n",
    "That is, suppose we change $X$ by an amount $\\Delta X$. We can then ask by how much the computed model parameters $\\theta^*$ change. If they change by a lot, then our method for computing them may be overly sensitive to perturbations. Instead, we might prefer one method over another one that is more sensitive to changes.\n",
    "\n",
    "Let's see how Algorithm 1 fares under small perturbations. But first, we'll need a method to generate a random perturbation of a certain maximize size."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise 4** (2 points). Implement a function that returns an $m \\times n$ matrix whose entries are uniformly randomly distributed in the interval, $[0, \\epsilon]$ for a given value of $\\epsilon$.\n",
    "\n",
    "> Hint: Check out Numpy's module for generating (pseudo)random numbers: [`numpy.random`](https://docs.scipy.org/doc/numpy/reference/routines.random.html)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "nbgrader": {
     "grade": false,
     "grade_id": "rand_eps",
     "locked": false,
     "schema_version": 1,
     "solution": true
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(3, 2)\n"
     ]
    }
   ],
   "source": [
    "def random_mat (m, n, eps):\n",
    "    #\n",
    "    # YOUR CODE HERE\n",
    "    from numpy.random import uniform\n",
    "    uniform_matrix = uniform(low=0, high=eps, size=(m,n))\n",
    "    return uniform_matrix\n",
    "\n",
    "print(random_mat(3, 2, 1e-3))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "nbgrader": {
     "grade": true,
     "grade_id": "rand_eps_test",
     "locked": true,
     "points": 2,
     "schema_version": 1,
     "solution": false
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "(Passed.)\n"
     ]
    }
   ],
   "source": [
    "# Test cell: `rand_eps_test`\n",
    "\n",
    "Z = random_mat (5, 3, 1e-2)\n",
    "assert Z.shape == (5, 3)\n",
    "assert ((Z >= 0) & (Z <= 1e-2)).all()\n",
    "print(\"\\n(Passed.)\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise 5** (2 points). Use your `random_mat()` function to write another function, `perturb_system(X, y, eps)`, that creates two \"perturbations\" to the system defined by `X` and `y`.\n",
    "\n",
    "1. Let $\\Delta X$ be the first perturbation. It should have the same dimensions as $X$, and its entries should lie in the interval $[-\\epsilon, \\epsilon]$. The value of $\\epsilon$ is given by `eps`.\n",
    "2. The second is $\\Delta y$, a small perturbation to the response variable, $y$. Its entries should also lie in the same interval, $[-\\epsilon, \\epsilon]$, \n",
    "\n",
    "Your function should return a perturbed system, $X + \\Delta X$ and $y + \\Delta y$, as a pair."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "nbgrader": {
     "grade": false,
     "grade_id": "perturb_system",
     "locked": false,
     "schema_version": 1,
     "solution": true
    }
   },
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$$\\Delta X = \\begin{bmatrix}\n",
       "  -0.09475516 & -0.07173508 & -0.00990991\\\\\n",
       "  -0.07537296 & 0.09721996 & -0.00737764\\\\\n",
       "  -0.0763913 & -0.04331462 & 0.02509546\\\\\n",
       "  0.00353655 & 0.00068309 & 0.02134772\\\\\n",
       "  0.09526083 & 0.07012873 & 0.06136301\\\\\n",
       "\\end{bmatrix}, \\quad \\Delta y = \\begin{bmatrix}\n",
       "  -0.07297956\\\\\n",
       "  -0.0559764\\\\\n",
       "  -0.04122086\\\\\n",
       "  -0.05770878\\\\\n",
       "  0.0214137\\\\\n",
       "\\end{bmatrix}$$"
      ],
      "text/plain": [
       "<IPython.core.display.Math object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "def perturb_system(X, y, eps):\n",
    "    #\n",
    "    # YOUR CODE HERE\n",
    "    from numpy.random import uniform\n",
    "\n",
    "    a, b, c, d = X.shape[0], X.shape[1], y.shape[0], y.shape[1]\n",
    "    \n",
    "    delta_X = uniform(low=-eps, high=eps, size=(a,b))\n",
    "    perturbed_X = X + delta_X\n",
    "    \n",
    "    delta_y = uniform(low=-eps, high=eps, size=(c,d))\n",
    "    perturbed_y = y + delta_y\n",
    "    \n",
    "    return perturbed_X, perturbed_y\n",
    "\n",
    "EPSILON = 0.1\n",
    "X_perturbed, y_perturbed = perturb_system(X, y, EPSILON)\n",
    "\n",
    "Delta_X = X_perturbed - X\n",
    "Delta_y = y_perturbed - y\n",
    "display(Math(r'\\Delta X = {}, \\quad \\Delta y = {}'.format(nparray_to_bmatrix(Delta_X[:5, :]),\n",
    "                                                          nparray_to_bmatrix(Delta_y[:5]))))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "nbgrader": {
     "grade": true,
     "grade_id": "perturb_system_test",
     "locked": true,
     "points": 2,
     "schema_version": 1,
     "solution": false
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "(Passed.)\n"
     ]
    }
   ],
   "source": [
    "# Test cell: `delta_X_test`\n",
    "\n",
    "Delta_X = X_perturbed - X\n",
    "Delta_y = y_perturbed - y\n",
    "\n",
    "assert Delta_X.shape == X.shape, \"`Delta_X` has shape {} instead of {}.\".format(Delta_X.shape, X.shape)\n",
    "assert (np.abs(Delta_X) <= EPSILON).all(), \"The perturbation lies outside the interval, [-{}, {}]\".format(EPSILON, EPSILON)\n",
    "\n",
    "assert Delta_y.shape == y.shape, \"`Delta_y` has shape {} instead of {}.\".format(Delta_y.shape, y.shape)\n",
    "assert (np.abs(Delta_y) <= EPSILON).all(), \"The perturbation lies outside the interval, [-{}, {}]\".format(EPSILON, EPSILON)\n",
    "print (\"\\n(Passed.)\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Sensitivity of Algorithm 1\n",
    "\n",
    "Let's now run the following code, which uses your code from above to perform a \"sensitivity experiment.\" In particular, the function `run_perturbation_trials()` will repeatedly perturb the system and measure the resulting change to the estimated $\\theta^*$.\n",
    "\n",
    "All of the estimated $\\theta^*$ are stored in an array, `Thetas_neq`. Each _column_ `k` of `Thetas_neq`, or `Thetas_neq[:, k]`, is one of the calculated estimates under a random perturbation of the system.\n",
    "\n",
    "The size of the random perturbation is set, by default, to `eps=0.01`. Recall that our synthetic problem consists of numerical values that are all greater than or equal to one, so this perturbation may be regarded as fairly small."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Unperturbed solution:\n",
      "[[ 1.]\n",
      " [ 1.]\n",
      " [ 1.]]\n",
      "First few perturbed solutions (columns):\n",
      "[[ 0.99114701  1.00566044  1.00451206  1.00325175  1.00769987]\n",
      " [ 1.00028385  1.00048433  1.00131233  0.99481327  0.99740573]\n",
      " [ 1.00002201  0.99994012  0.99966334  1.00064643  1.00033515]]\n"
     ]
    }
   ],
   "source": [
    "def run_perturbation_trials(solver, X, y, eps=0.01, trials=100):\n",
    "    Thetas = np.zeros((X.shape[1], trials)) # Store all computed thetas\n",
    "    for t in range(trials):\n",
    "        X_p, y_p = perturb_system(X, y, eps)\n",
    "        Thetas[:, t:t+1] = solver(X_p, y_p)\n",
    "    return Thetas\n",
    "\n",
    "Thetas_neq = run_perturbation_trials(solve_neq, X, y)\n",
    "\n",
    "print(\"Unperturbed solution:\")\n",
    "print(theta_neq)\n",
    "\n",
    "print(\"First few perturbed solutions (columns):\")\n",
    "print(Thetas_neq[:, :5])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here is a quick plot of the that shows two coordinates of the true parameters (red star), compared to all perturbed estimates (blue points). We would have more confidence in the algorithm's computed solutions if it did not appear to be too sensitive to changes in the input.\n",
    "\n",
    "Since $\\theta$ may have more than two coordinates, the code below shows the first two coordinates."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAARkAAAEKCAYAAAAmUiEiAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAGipJREFUeJzt3XuQVeWd7vHvQ9MgaAwYKEWuTmIuZmSA9PGSxFsyKpgL\nxpgYKo5KUkNZo9HUVFLxVmPGOHGi42RiOSVh5hDEqcGceJkwDgl6jIaMQo6NIqITFD2J0HCUCUJU\nUBr4nT/Wu2Hb9GXv7n73bjbPp2rV3vt911rv2w376bXedVNEYGaWy6B6d8DMGptDxsyycsiYWVYO\nGTPLyiFjZlk5ZMwsK4eMmWXlkDGzrBwyZpbV4Hp3oL+MGjUqJk2aVO9umDW0lStX/ndEjK5mmYYJ\nmUmTJtHa2lrvbpg1NEm/q3YZ7y6ZWVYOGTPLyiFjZlk5ZMwsq2whI2m+pFclremiXpJuk7RO0mpJ\n08rqdktalabFufpoZvnl3JJZAEzvpn4GcGya5gB3lNXtiIgpafpsvi6aWW7ZQiYilgFbupllJrAw\nCiuAEZLG5OqPmdVHPcdkxgLryz5vSGUAh0hqlbRC0rm175qZ9ZeBejLexIhok/RHwC8kPRMRL3ac\nSdIcil0tJkyYUOs+mlkF6rkl0waML/s8LpUREaXXl4BHgamdrSAi5kVES0S0jB5d1ZnOZlYj9QyZ\nxcBF6SjTScC2iNgkaaSkoQCSRgEfA56rYz/NrA+y7S5JWgScDoyStAG4HmgGiIi5wBLgHGAdsB2Y\nnRb9EPBDSXsoQvBvI8IhY3aAyhYyETGrh/oALuuk/HHg+Fz9MrPa8hm/ZpaVQ8bMsnLImFlWDhkz\ny8ohY2ZZOWTMLCuHjJll5ZAxs6wcMmaWlUPGzLJyyJhZVg4ZM8vKIWNmWTlkzCwrh4yZZeWQMbOs\nHDJmlpVDxsyycsiYWVYOGTPLyiFjZlk5ZMwsK4eMmWWVLWQkzZf0qqQ1XdRL0m2S1klaLWlah/rD\nJW2QdHuuPppZfjm3ZBYA07upnwEcm6Y5wB0d6r8DLMvSMzOrmWwhExHLgC3dzDITWBiFFcAISWMA\nJH0EOBJ4MFf/zKw26jkmMxZYX/Z5AzBW0iDgVuAbPa1A0hxJrZJaN2/enKmbZtYXA3Hg9y+AJRGx\noacZI2JeRLRERMvo0aNr0DUzq9bgOrbdBowv+zwulZ0MnCLpL4DDgCGS3oiIq+rQRzPro3qGzGLg\nckl3AycC2yJiE/Dl0gySLgFaHDBmB65sISNpEXA6MErSBuB6oBkgIuYCS4BzgHXAdmB2rr6YWf1k\nC5mImNVDfQCX9TDPAopD4WZ2gBqIA79m1kAcMmaWlUPGzLJyyJhZVg4ZM8vKIWNmWTlkzCwrh4yZ\nZeWQMbOsHDJmlpVDxsyycsiYWVYOGTPLyiFjZlk5ZMwsK4eMmWXlkDGzrBwyZpaVQ8bMsnLImFlW\nDhkzy8ohY2ZZZQsZSfMlvSppTRf1knSbpHWSVkualsonSnpS0ipJz0q6NFcfzSy/nFsyC4Dp3dTP\nAI5N0xzgjlS+CTg5IqZQPFnyKklHZ+ynmWWULWQiYhmwpZtZZgILo7ACGCFpTETsjIi30zxDc/bR\nzPKr5xd4LLC+7POGVIak8ZJWp/rvRcTGzlYgaY6kVkmtmzdvzt5hM6vegNxKiIj1ETEZeB9wsaQj\nu5hvXkS0RETL6NGja9tJM6tIPUOmDRhf9nlcKtsrbcGsAU6pYb/MrB/VM2QWAxelo0wnAdsiYpOk\ncZKGAUgaCXwcWFvHfppZHwzOtWJJi4DTgVGSNgDXA80AETEXWAKcA6wDtgOz06IfAm6VFICAv4uI\nZ3L108zyyhYyETGrh/oALuuk/CFgcq5+mVltDciBXzNrHA4ZM8uqopCRdKSk/ynpZ+nzcZK+mrdr\nZtYIKt2SWQAsBUqn9z8PfD1Hh8yssVQaMqMi4n8BewAiYhewO1uvzKxhVBoyb0p6DxAApfNasvXK\nzBpGpYew/5Li5Ln3SnoMGA18IVuvzKxhVBoyzwKnAR+gOEFuLT4yZWYVqDQolkfEroh4NiLWREQ7\nsDxnx8ysMXS7JSPpKIrbLwyTNJViKwbgcGB45r6ZWQPoaXfpbOASiiuk/76s/HXgmkx9MrMG0m3I\nRMSdwJ2SPh8R99aoT2bWQCoa+I2IeyV9CvgwcEhZ+Q25OmZmjaHSywrmAhcAX6MYl/kCMDFjv8ys\nQVR6dOmjEXER8FpE/DVwMvD+fN0ys0ZRacjsSK/b0+NJ2oExebpkZo2k0pPxHpA0ArgFeJLi8oJ/\nztYrM2sYlQ78fie9vVfSA8AhEeFrl8ysRxXfflPSR4FJpWUkERELM/XLzBpERSEj6S7gvcAq9t3i\nIQCHjJl1q9ItmRbguHTzbzOzilV6dGkNcFTOjphZY6r4znjAc5KWSlpcmrpbQNJ8Sa9KWtNFvSTd\nJmmdpNWSpqXyKZKWS3o2lV9Q3Y9kZgNJpbtL3+7FuhcAt9P1uM0M4Ng0nQjckV63AxdFxAvpnJyV\nkpZGxNZe9MHM6qzSQ9i/rHbFEbFM0qRuZpkJLEzjPCskjZA0JiKeL1vHRkmvUtyJzyFjdgDqdndJ\n0n+m19cl/aFsel3SH/rY9lhgfdnnDamsvP0TgCHAi31sy8zqpKdbPXw8vb6rNt3ZR9IY4C7g4ojY\n08U8c4A5ABMmTKhh78ysUj3dGe+I7uojYksf2m4Dxpd9HpfKkHQ48B/AtRGxopv25wHzAFpaWnx4\n3WwA6mlMZiXFSXcCJgCvpfcjgJeBY/rQ9mLgckl3Uwz4bouITZKGAPdTjNfc04f1m9kA0NPu0jEA\nkv4JuD8ilqTPM4Bzu1tW0iLgdGCUpA3A9UBzWu9cYAlwDrCO4ojS7LToF4FTgfdIuiSVXRIRq6r8\n2cxsAFAlJ/FKeiYiju+prJ5aWlqitbW13t0wa2iSVkZESzXLVHqezEZJ1wH/kj5/GdhYTUNmdnCq\n9IzfWRTnqtwP3Jfez8rVKTNrHD1uyUhqAq6JiCtr0B8zazA9bslExG7g4zXoi5k1oErHZJ5KF0T+\nBHizVBgR92XplZk1jEpD5hDg98AnysqCYnzGzKxLlV4gObvnuczM9lfpw93GSbo/3R/mVUn3ShqX\nu3NmduCr9BD2jyguAzg6Tf+eyszMulVpyIyOiB9FxK40LaA4V8bMrFuVhszvJV0oqSlNF1IMBJuZ\ndavSkPkKxYWL/y9N57PvgkYzsy5VenTpd8BnM/fFzBqQjy6ZWVY+umRmWfnokpll5aNLZpZVb44u\nbcJHl8ysQj66ZGZZVXp06U5JI8o+j5Q0P1+3zKxRVLq7NLn8WdQR8RowNU+XzKyRVBoygySNLH1I\nD32r9F40ZnYQqzRkbgWWS/qOpBuBx4Gbu1tA0vx04t6aLuol6TZJ6yStljStrO7nkrZKeqDSH8TM\nBqaKQiYiFgLnAa9QHGE6LyLu6mGxBcD0bupnAMemaQ5wR1ndLcCfVdI3MxvYKt2SISKei4jbgZ0R\n8VwF8y8DuntW9kyKR9FGet71CElj0rIPA69X2jczG7gqDpkyl/ZT22OB9WWfN6QyM2sgvQkZ9Xsv\neknSHEmtklo3b95c7+6YWSd6DBlJH5T0SUmHpaLPpPLuxlsq0QaML/s8LpVVLCLmRURLRLSMHu1L\nqcwGom5DRtIVwE+BrwFrJM2MiA2p+rt9bHsxcFE6ynQSsC0iNvVxnWY2wPR0rsufAx+JiDckTQLu\nkTQpIn5AD7tNkhYBpwOjJG0ArgeaASJiLrAEOAdYB2yn7FooSb8CPggclpb9akQsrfqnM7O66ylk\nBkXEGwAR8VtJp1MEzUR6CJmImNVDfQCXdVF3Sg/9MrMDRE9jMq9ImlL6kALn08Ao4PicHTOzxtBT\nyFxEcfLdXummVRcBp2brlZk1jG53l8oGeTure6z/u2NmjaY358mYmVXMIWNmWTlkzCwrh4yZZeWQ\nMbOsHDJmlpVDxsyycsiYWVYOGTPLyiFjZlk5ZMwsK4eMmWXlkDGzrBwyZpaVQ8bMsnLI2MAXAffd\nV7zaAcchYwPfqlXw+c/D00/XuyfWCz3dSNys/hYtAgnuvhumTOl5/h7821Nt3LJ0LRu37uDoEcP4\n5tkf4Nyp1T28tD/WcbBQNMgmaEtLS7S2tta7G1amqy9iVV/QCBg/HtraYNw4WL9+v/W3bd2BgNL/\nZAEffe8R/Pb3O2jbuoMmid0Re1/L5y0ZMayZb3/2w3v79+3Fz7J1RzsAI4c386nJY3jkN5vZuHUH\nI4Y388Zbu2jfs28tw5qbuOm842n93RYW/Xr93vZmnTielolH7P15RwxvJgK27Wjv9Gcf6OElaWVE\ntFS1TK6QkTSf4skGr0bEH3dSL+AHFM9e2g5cEhFPprqLgevSrDdGxJ09teeQGTg6fklLmgeJE44Z\nyeMvbnnHl7z0pS+FwIhhzUiwdXs7H3v7FRbMvZzBb+2A4cPhiSe47vnd/OuvX2bPAPv7OKx5EDva\n9+xX3jRI7O6is6VwKoXb1fc9w4723Z3WDwS9CZmcu0sLgNuBhV3UzwCOTdOJwB3AiZKOoHgQXAvF\n/72VkhZHxGsZ+3pQqOavZMd5z/jgaB75zeZiq0H7xmAHCfZE8dc+gv2CpVz7nuCxF7fsVx4AEZz8\nf5/i0PYd76g784UV7G5vZzDw9tvtPPClK9j8vhM5s2yeN5uH8Z+TphS7VHXUWcAAXQZMscxublm6\nlnOnjuWWpWvfETAd6w9U2UImIpalp052ZSawMD3kbYWkEZLGUDx18qGI2AIg6SFgOrAoV18bWVe7\nFG1bd3D1fc/sna9joPz4/6zfuzvQtnUH/7Li5b3zlm/8lr4/r23vOlwq0bxnF99dejsTtr3CG82H\nECqOSSiCobt3ATB0dztn/+Zxzlq7PNXt4bD2t3j53UfyyT+fS3tTc5/6UC8bt+54x2tX9Qeqeg78\njgXWl33ekMq6Kj+g9GWrodr98O7GPso3vzv+Pd3Rvpuv/3jVO8o6BkqttDc1M/0rt3Pzz37AJ9Y9\nwfD27Z3Od1jZls72wUNZ/MFTuGrGFQdswAAcPWLY3te2TgKlVH+gOqCPLkmaA8wBmDBhQp17s0/H\nL3f5VkPH8Kh03kqDpHz5zja/B7LtQ4Zx+cyrmLnmEW568HaG7GpncOy/C7JLg9g5uJmrzr6cxR8+\now497Vxvx2S+efYHAPjm2R/odEymVH+gqud5Mm3A+LLP41JZV+X7iYh5EdESES2jR4/O1tFqdbdv\n3Zt5S0HStnUHwb4gKQVPV8sfqJvZP/3jM5g++3bamzr/G9jeNJjps28fYAHTxE3nTebCkybQlMaG\nmiQuPGkCt37hTxg7YhiiGLsaMawZAWNHDHvHoO65U8dy03nH7523Y/2Bqp5bMouByyXdTTHwuy0i\nNklaCnxX0sg031nA1fXqZG9Us29dyby9CZLSFk9nm98Hgt2DmuhqGFcR7Gpqytr+8LRVUhqjundl\n237/BkMHD2Lnrj3v2LI8d+pYbjx3/8fEVxoUpXU0kmwhI2kRxSDuKEkbKI4YNQNExFxgCcXh63UU\nh7Bnp7otkr4DPJFWdUNpEPhAUc2+dSXz9iZISv/xO25+Hyimr31s7xjS24MG07xnN7sGNTFkzy5C\n4jPrlvPDaTP3W650tKuzc2E6Gpt+R0CPY2Ll57oMxPNXBrKcR5dm9VAfwGVd1M0H5ufoVy1Us29d\nybzVBklp+dKXoHR0qZ4OHdLElPHv3u8cmeZB4rBDBrN1ezvvHtbMH95qZ0/AF9b8b4bt2smbzUNZ\nfdT7uf+r3+LmR/4JWlsZ9uabXP3Kr/nQBdd3+cUvP6rWUWfnnvQUGI24hVErPuM3k/48utTTSVqV\ntNXZOsodOqSJnbv27HcW6+c/Mnbv+TGlk+U6bgF0FWAjhzfz1F+dVfXP+v2Fv+TRf/gyO5uaueXU\ni3jra1dw4+cmw5498P3vwzXXwM6dsGkTHHVUp21X06ZVbkCd8VtrAy1k+lt/Xm/TWWBUfbp/h/X2\n65mqDzwAV14J99wDU6fuX//kk3D++XDbbfDpT1e/fus1h0wDh8xA562Fg8NAu6zADiIes7Cu+H4y\nZpaVQ8bMsnLImFlWDhkzy8ohY2ZZOWTMLCuHjJll5ZAxs6wcMmaWlUPGzLJyyJhZVg4ZM8vKIWNm\nWTlkzCwrh4yZZeWQMbOsHDJmlpVDxsyycsiYWVZZQ0bSdElrJa2TdFUn9RMlPSxptaRHJY0rq/ue\npDVpuiBnP80sn2whI6kJ+EdgBnAcMEvScR1m+ztgYURMBm4AbkrLfgqYBkyheITtNyQdnquvZpZP\nzi2ZE4B1EfFSROwE7gY6Plf0OOAX6f0jZfXHAcsiYldEvAmsBqZn7KuZZZIzZMYC68s+b0hl5Z4G\nzkvvPwe8S9J7Uvl0ScMljQLOAMZ3bEDSHEmtklo3b97c7z+AmfVdvQd+vwGcJukp4DSgDdgdEQ8C\nS4DHgUXAcmC/56tGxLyIaImIltGjR9ew22ZWqZwh08Y7tz7GpbK9ImJjRJwXEVOBa1PZ1vT6NxEx\nJSLOBAQ8n7GvZpZJzpB5AjhW0jGShgBfAhaXzyBplKRSH64G5qfyprTbhKTJwGTgwYx9NbNMsj2m\nNiJ2SbocWAo0AfMj4llJNwCtEbEYOB24SVIAy4DL0uLNwK8kAfwBuDAiduXqq5nlo4iodx/6RUtL\nS7S2tta7G2YNTdLKiGipZpl6D/yaWYNzyJhZVg4ZM8vKIWNmWTlkzCwrh4yZZeWQMbOsHDJmlpVD\nxsyycsiYWVYOGTPLyiFjZlk5ZMwsK4eMmWXlkDGzrBrmfjKSNgO/q0FTo4D/rkE7A7F9t31wtd1Z\n+xMjoqobajdMyNSKpNZqb9rTKO277YOr7f5q37tLZpaVQ8bMsnLIVG/eQdy+2z642u6X9j0mY2ZZ\neUvGzLJyyJSRNF3SWknrJF3VSf1ESQ9LWi3pUUnjUvkZklaVTW9JOrcWbae6myU9K+m/JN2m9MCq\nGrX9PUlr0nRBNe2m5edLelXSmi7qlX6mdan9aWV1F0t6IU0X17jtn0vaKumBatvtS9uSpkhanv69\nV/fmd97H9idKejL9P39W0qU9NhYRnopdxibgReCPgCHA08BxHeb5CXBxev8J4K5O1nMEsAUYXou2\ngY8Cj6V1NFE8N/z0GrX9KeAhiocEHkrx1NDDq/y9nwpMA9Z0UX8O8DOKRxWfBPy67Pf8Unodmd6P\nrEXbqe6TwGeAB3r5/623P/f7gWPT+6OBTcCIGrY/BBia3h8G/BY4uru2vCWzzwnAuoh4KSJ2AncD\nMzvMcxzwi/T+kU7qAc4HfhYR22vUdgCHkP7xKZ6++UqN2j4OWBYRuyLiTWA1ML2KtomIZRSh3JWZ\nwMIorABGSBoDnA08FBFbIuI1irCrVdtExMPA69W01x9tR8TzEfFCWsdG4FWgqpPj+tj+zoh4O80z\nlAr2hhwy+4wF1pd93pDKyj0NnJfefw54V+mZ3WW+BCyqVdsRsZzii78pTUsj4r9q0XYqny5puKRR\nwBnA+Cra7kv/Kul3rrZroce2JZ1A8cflxVq2L2m8pNWp/nsp7LrkkKnON4DTJD0FnAa0AbtLlemv\n3PEUz/+uSduS3gd8CBhH8Z/gE5JOqUXbEfEgsAR4nCJYl1P2+7B80v+1u4DZEbGnlm1HxPqImAy8\nD7hY0pHdze+Q2aeNd/4VHpfK9oqIjRFxXkRMBa5NZVvLZvkicH9EtNew7c8BKyLijYh4g2I/+uQa\ntU1E/E1ETImIMyn235+vou2+9K/Hfmdsuxa6bFvS4cB/ANemXZmatl+StmDWAN3+UXPI7PMEcKyk\nYyQNodjtWVw+g6RRkkq/s6uB+R3WMYvqd5X62vbLFFsZgyU1U2xpVLO71Ou2JTWVdhclTQYmAw9W\n0XYlFgMXpaMdJwHbImITxdbiWZJGShoJnEX/b0F21XYtdNp2+je6n2K85J46tD9O0jCA9Hv/OLC2\n2zX1ZmS8USeKEfXnKfZxr01lNwCfTe/PB15I8/wzaZQ91U2iSPpBtWyb4ujQDymC5Tng72vY9iGp\nzeeAFcCUXrS9iGIsqZ1iv/+rwKXApalewD+mvj0DtJQt+xVgXZpm17jtXwGbgR1p2bNr0TZwYVpm\nVdlUs987cCbFAP/T6XVOT235jF8zy8q7S2aWlUPGzLJyyJhZVg4ZM8vKIWNmWTlkrMersGvYjzfS\n69GSen0OiKSvSxrefz2zvvAh7IOcpCaK81/OpDhf4glgVkQ8l7ndwRGxq0PZGxFxWD+s+7cU53XU\n8y7/lnhLxiq5ChtJ35L0jKSnJf1tKpsiaUW638j96QzQ7soflfQPklqBK9NZxsvTem8sa2tS6T4n\nki6RdJ+K+7e8IOnmsvnukNSa7mvy16nsCopbIDwi6ZFUdlZq50lJP5HU5yCzyjlkrJKrfWdQBM+J\nEfEnQOmLvhD4VhQXyz0DXN9DOcCQiGiJiFuBHwB3RMTxFGefdmUKcAHFxacXSCpdU3NtFI/rmExx\nacXkiLgN2AicERFnpKvDrwP+NCKmAa3AX1b0m7F+4ZCxSvwp8KNI98iJiC2S3k1xs6RfpnnuBE7t\nqrxsXT8ue/8x9l3rdVc37T8cEdsi4i2KSxgmpvIvSnoSeAr4MMX9bTo6KZU/JmkVcHHZ8lYDg+vd\nAau7zq623Zy+kAB/1c/tvdnhcyWDgm+Xvd8NDJZ0DMUtKP5HRLwmaQHFtVQdieLmVrN601nrO2/J\nWGdXYd8Txe0bpkTEYoq7zs0uHbGRdEREbANeK7t3zZ8Bv+yqvIu2H0vtAXy5yn4fThFY29L9TGaU\n1b0OvCu9XwF8LN13B0mHSnp/lW1ZH3hL5iAXEbskXU5xm4QmYH5EPNthnp9LmgK0StpJcaOqayh2\nPeam8HkJmJ0W6aq8oyuBf5X0LeCnVfb7aRU30foNxZjSY2XV84CfS9qYxmUuARZJGprqr6P/73tj\nXfAhbDPLyrtLZpaVQ8bMsnLImFlWDhkzy8ohY2ZZOWTMLCuHjJll5ZAxs6z+Pz9KALpzTM9zAAAA\nAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7fb410e8a1d0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Makes a 2-D scatter plot of given theta values.\n",
    "# If the thetas have more than two dimensions, only\n",
    "# the first and last are displayed by default.\n",
    "# (Override by setting ax and ay.)\n",
    "\n",
    "def scatter_thetas(Thetas, theta_true=None, ax=0, ay=-1, xylim=None, title=None):\n",
    "    import matplotlib.pyplot as plt\n",
    "    assert type(Thetas) is np.ndarray and Thetas.shape[0] >= 2\n",
    "    scatter(Thetas[ax, :], Thetas[ay, :])\n",
    "    xlabel('{}-coordinate'.format(ax if ax >= 0 else Thetas.shape[0]+ax))\n",
    "    ylabel('{}-coordinate'.format(ay if ay >= 0 else Thetas.shape[0]+ay))\n",
    "    if xylim is not None:\n",
    "        axis(xylim)\n",
    "    else:\n",
    "        axis('equal')\n",
    "    if theta_true is not None:\n",
    "        assert type(theta_true) is np.ndarray and theta_true.shape[0] >= 2 and theta_true.shape[1] == 1\n",
    "        scatter(theta_true[ax], theta_true[ay], marker='*', color='red', s=15**2)\n",
    "    if title is not None:\n",
    "        plt.title(title)\n",
    "        \n",
    "def calc_lims(x, buffer=0.1):\n",
    "    xmin = x.min()\n",
    "    xmax = x.max()\n",
    "    dx = (xmax - xmin) * buffer\n",
    "    return xmin-dx, xmax+dx\n",
    "        \n",
    "scatter_thetas(Thetas_neq, theta_true=theta_true, ax=0, ay=2)\n",
    "axis('square');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You should observe that the change in the estimates are of the same order as the perturbation. So for this example system, the algorithm seems reliable enough."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "## Stress-testing Algorithm 1\n",
    "\n",
    "This experiment suggests all is fine. But what should we _expect_ to happen?\n",
    "\n",
    "We've prepared another [notebook](./notes-cond.ipynb) that walks you through an analysis of solving linear systems. It turns out you can estimate how hard it is to solve a linear system using a measure called the _condition number_. We can denote the condition number of solving a system by $\\kappa(X)$ where $X$ is the matrix. The larger this number is, the more sensitive the problem is.\n",
    "\n",
    "In Numpy, there is a condition number estimator that will tell us approximately what the condition number is for a given matrix. Let's compare $\\kappa(X)$ and $\\kappa(C) = \\kappa(X^TX)$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$$\\kappa(X) \\approx 1.07\\times 10^{2}$$"
      ],
      "text/plain": [
       "<IPython.core.display.Math object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$$\\kappa(X^T X) \\approx 1.15\\times 10^{4}$$"
      ],
      "text/plain": [
       "<IPython.core.display.Math object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$$\\kappa(X)^2 \\approx 1.15\\times 10^{4}$$"
      ],
      "text/plain": [
       "<IPython.core.display.Math object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "cond_X = np.linalg.cond(X)\n",
    "cond_XTX = np.linalg.cond(X.T.dot(X))\n",
    "\n",
    "assert 1. <= cond_X <= 3e3\n",
    "assert 1. <= cond_XTX <= 6e6\n",
    "\n",
    "show_cond_fancy(cond_X, 'X')\n",
    "show_cond_fancy(cond_XTX, 'X^T X')\n",
    "show_cond_fancy(cond_X**2, 'X', opt='^2')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Ill-conditioning.** As it happens, $\\kappa(C)$ is roughly the **square** of $\\kappa(X)$. So, by forming $C$ explicitly and then trying to solve a system based on it, we make the problem _more_ difficult. Indeed, if the problem is ill-conditioned enough, this algorithm based on directly constructing the normal equations will produce very different results even under small changes, and we call the algorithm _unstable_.\n",
    "\n",
    "In this particular example, the condition numbers are not very \"big.\" You would be more concerned if the condition numbers were close to $1/\\epsilon$, where $\\epsilon$ is machine epsilon. In double-precision, recall that $\\epsilon_d \\approx 10^{-15}$, so the values shown above is nothing to be worried about.\n",
    "\n",
    "But what if we had a \"hard\" problem, that is, one whose condition number is large? The synthetic data generator allows us to create such a problem by making the problem bigger. Let's try that next."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "First few rows of data:\n",
      "True parameter estimates:\n",
      "[[ 1.]\n",
      " [ 1.]\n",
      " [ 1.]\n",
      " [ 1.]\n",
      " [ 1.]\n",
      " [ 1.]\n",
      " [ 1.]]\n"
     ]
    },
    {
     "data": {
      "text/latex": [
       "$$\\kappa(X_h) \\approx 1.72\\times 10^{12}$$"
      ],
      "text/plain": [
       "<IPython.core.display.Math object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$$\\kappa(X_h^T X_h) \\approx 2.82\\times 10^{23}$$"
      ],
      "text/plain": [
       "<IPython.core.display.Math object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Generate a \"hard\" problem\n",
    "m_hard, n_hard = 100, 6\n",
    "X_hard, y_hard, theta_hard_true = gen_problem(m_hard, n_hard)\n",
    "\n",
    "df_hard = make_data_tibble(X_hard, y_hard)\n",
    "print(\"First few rows of data:\")\n",
    "df_hard.head()\n",
    "print(\"True parameter estimates:\\n{}\".format(theta_hard_true))\n",
    "\n",
    "cond_X_hard = np.linalg.cond(X_hard)\n",
    "cond_XTX_hard = np.linalg.cond(X_hard.T.dot(X_hard))\n",
    "\n",
    "name_X_hard = 'X_h'\n",
    "show_cond_fancy(cond_X_hard, name_X_hard)\n",
    "show_cond_fancy(cond_XTX_hard, '{}^T {}'.format(name_X_hard, name_X_hard))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "These condition numbers are much larger. So, let's run the same sensitivity experiment as before, and see how the estimate varies for the hard problem. How does it compare to the well-conditioned case?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Residual norm for one of the trials:\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "array([ -1.39543976e+13,  -1.39543973e+13,  -1.39543743e+13,\n",
       "        -1.39526579e+13,  -1.38159240e+13,  -2.47310601e+12,\n",
       "         9.65186305e+14])"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEKCAYAAAD9xUlFAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAGC5JREFUeJzt3X20XHV97/H395wcyEHUQJMlEAhBS3VBoUk9KhZrxYcL\nWAQUrLJ8KGrL6l1QoLUsKaU+cO2tLat67cUFlxYEfAhWnhopmlKlWlmAnEAgJDQ111sgAUuEJBIS\n4ZzwvX/Mns3kZM7MnId95pyT92utWZm9Z8+e75okv8/ev99v74nMRJIkgJ5uFyBJmj4MBUlSyVCQ\nJJUMBUlSyVCQJJUMBUlSyVCQJJUMBUlSyVCQJJXmdLuAsZo/f34uXry422VI0oyycuXKn2Xmgnbb\nzbhQWLx4MYODg90uQ5JmlIh4pJPt7D6SJJUMBUlSyVCQJJUMBUlSyVCQJJUMBUlSyVCQJJUMBUlS\nyVCQJJUMBUlSyVCQJJUMBUlSyVCQJJUMBUlSyVCQJJUMBUlSyVCQJJUMBUlSyVCQJJUMBUlSyVCQ\nJJUMBUlSyVCQJJUMBUlSyVCQJJUMBUlSyVCQJJUMBUlSyVCQJJUMBUlSyVCQJJUqC4WIOCQi7oiI\ntRGxJiLOa7JNRMTfRsT6iHgwIn69qnokSe3NqXDfw8DHM/O+iHgpsDIibs/MtQ3bnAgcXjzeAFxe\n/ClJ6oLKzhQy84nMvK94/gzwMLBwxGanANdlzd3AvIg4sKqaJEmtTcmYQkQsBpYC94x4aSHwWMPy\nBnYPDiLirIgYjIjBTZs2VVWmJO3xKg+FiNgXuBE4PzN/Pp59ZOaVmTmQmQMLFiyY3AIlSaVKQyEi\n+qgFwtcy86Ymm2wEDmlYPrhYJ0nqgipnHwVwFfBwZn5+lM2WAx8uZiEdA2zNzCeqqkmS1FqVs4+O\nBT4ErI6IVcW6i4BFAJl5BXAb8E5gPbAd+EiF9UiS2qgsFDLzh0C02SaBs6uqQZI0Nl7RLEkqGQqS\npJKhIEkqGQqSpJKhIEkqGQqSpJKhIEkqGQqSpJKhIEkqGQqSpJKhIEkqGQqSpJKhIEkqGQqSpJKh\nIEkqGQqSpJKhIEkqGQqSpJKhIEkqGQqSpJKhIEkqGQqSpJKhIEkqGQqSpJKhIEkqGQqSpJKhIEkq\nVRYKEXF1RDwZEQ+N8vpbImJrRKwqHp+sqhZJUmfmVLjva4DLgOtabPNvmXlShTVIksagsjOFzPwB\n8HRV+5ckTb5ujym8MSIeiIhvR8SRXa5FkvZ4VXYftXMfcGhmbouIdwK3AIc32zAizgLOAli0aNHU\nVShJe5iunSlk5s8zc1vx/DagLyLmj7LtlZk5kJkDCxYsmNI6JWlP0rVQiIgDIiKK568vanmqW/VI\nkirsPoqIZcBbgPkRsQH4FNAHkJlXAKcD/z0ihoEdwPszM6uqR5LUXmWhkJlntHn9MmpTViVJ00S3\nZx9JkqYRQ0GSVOooFCLiFRFxVUR8u1g+IiI+Vm1pkqSp1umZwjXACuCgYvk/gPOrKEiS1D2dhsL8\nzPwH4AWAzBwGdlZWlSSpKzoNhWcj4peABIiIY4CtlVUlSeqKTqek/jGwHHhVRNwJLADeW1lVkqSu\n6DQU1gC/BbwaCGAdzlySpFmn04b9rswczsw1mflQZg4Bd1VZmCRp6rU8U4iIA4CFQH9ELKV2lgDw\nMmCfimuTJE2xdt1HxwNnAgcDn29Y/wxwUUU1SZK6pGUoZOa1wLURcVpm3jhFNUmSuqSjgebMvDEi\nfhs4EpjbsP6SqgqTJE29Tm9zcQXwPuAPqY0rvBc4tMK6JEld0Onso9/IzA8DmzPzM8AbgV+prixJ\nUjd0Ggo7ij+3R8RBwBBwYDUlSZK6pdOL126NiHnApcB91G538feVVSVJ6opOB5r/R/H0xoi4FZib\nmd77SJJmmY5/jjMifgNYXH9PRJCZ11VUlySpCzoKhYj4CvAqYBUv3jI7AUNBkmaRTs8UBoAjMjOr\nLEaS1F2dhsJDwAHAExXWIkka4Zb7N3LpinU8vmUHB83r54LjX82pSxdW9nmdhsJ8YG1E/Ah4rr4y\nM0+upCpJErfcv5E/vWk1O4ZqvfYbt+zgT29aDVBZMHQaCp+u5NMlaQpN9VH3RF26Yl0ZCHU7hnZy\n6Yp13Q2FzPx+JZ8uTVeZcPPN8O53Q0T77TWpmjXeUGskN27ZQW8EOzNZOK+f416zgDv+fdMu6xtf\nrzf83TjqnqjHt+wY0/rJ0O73FH6YmW+KiGcofp+5/hKQmfmyyiqTumnVKjjtNLj/fliypNvVTCu3\n3L+RTy9fw5YdQ7u9FsAHjlnEZ089qu1R+cjXGxv34MUGZ+OWHVxwwwOQMPRCbe3OYs7Lxi07+Ord\nj5b7rK9vfL3e8HfjqHuiDprXz8YmAXDQvP7KPrPdrbPfVPz50soqkKajZctqZwjXXz/jQ6He+DYe\nSfcEFO0r8/r7+PTJRwI0baRHHq1f8M0HysZ5pAS+evej/L9N27jv0a2jHpU3O2pvbNxH7n1o5/gn\nPtYb/m4cdU/UBce/epfvCaC/r7f8u6hCtJplGhH7t3pzZj496RW1MTAwkIODg1P9sdqTZMIhh8DG\njXDwwfDYY5V9VOPRcn9fDzuGXyATeiM45pX7sebxZ3Y7Iq8fRS9saKgbj9z326eP3z76QG594Imm\nR/PN9AC9vdGy8e3v62VuXw+bt3e2z2YWzuvnzgvfyrGf+17TI+CqBKMfdddrmq4maxwkIlZm5kC7\n7dqNKayk9u8vgEXA5uL5POBR4LAWBVwNnAQ8mZm/2uT1AL4IvBPYDpyZmfe1K1iq3Nq1DD/1NHOA\n7T/dxCm/dzk/nr9ol0326eth+9ALu6x7yV69PPv8rt0TzRrwxhBo3Efj852Z3Pl/mx9zjexW2bkz\naaxk8/ahXY66O/EC8EKbo/EdQzt3634Zq/pR+VQfndcb06k+6p4Mpy5dOKXdW+26jw4DiIi/A27O\nzNuK5ROBU9vs+xrgMka/6vlE4PDi8Qbg8uJPadwaj6pe3t9HBGzZPsQ+zRrsTD7Z/zhHvLSHW1Y9\nztPP1mZbv+PHd/OuoSHmAL25k7PuuZF/OXzXf5rP9vXzw8VLdhmEHrl/aN0vPjJUxmMi3SrdUO8L\nH+2ofTR9vbHLdzcW9Ya/3rDOpNlH3dCy+6jcKGJ1Zh7Vbl2T9y0Gbh3lTOH/AP+amcuK5XXAWzKz\n5QVydh/NTKMNKo5svEf+R23WHz5S41H4yCPBVvp2DvHdv/sDFm39L7b1zSWjdif5yGTfoRcbrG19\n/WTR+Ee+wL5Dv+DRl7+Ct/3+FQz19k3oe5lp5vX38exzw20b52Nftf8uYwpQa5z/8j1HNR1TaGa0\ns6zxzj7a001W91Hd4xFxMfDVYvkDwOPjLa6wEGjsrN1QrPOq6WlgtH7MVlMFW80yueCGB8qj2pGD\nio393o0DkrBrI98sEBrfM7evZ0zdG0O9fZzw0cv4629/kbeuv5d9hrY33a4xILbP2Zvlr/lNLjzx\n3FkVCJ2OKdQHpCc6+6jZUXuzge2RjbmNe/U6PVPYH/gU8GZq4f0D4JJ2A81tzhRuBT6XmT8slr8L\nfCIzdzsNiIizgLMAFi1a9NpHHnmkbc2z2XiPuusuvmU1y+55bJdGduQRV+OUQKg1CKe9diE3rty4\nS8Pb7LS+8YgQYOkl/zzmwcmFRTfDVA1GnvLQHfzlP1/GXsNDzMndu3WGo4fn5/Rx4fHnsPzI46ak\npk709cZuYwqt1P9exzv7yEZ55ur0TKFtKEREL/BXmfkn4yhiMXYf7WYiswlGHnW3M7KBvviW1WMe\nhKwbrfummcYZHYsv/Kcxf1a9p34qe8wXbX6CFVefTf/w87u9tmPOXhz/0S/x6H7j/8HBTvrF6432\nZM8+2m+fPj71riNt1Pdgk9Z9lJk7I+JNk1PWLpYD50TE9dQGmLe2C4SqNeu/Htl32eo/WieN/USv\nqvzMt9aMaXBx5MU5y+4Z//TKTgMBJj675KBxnCnM6+9j646hMQdJD7XZNzt7eolRtolMhnt7W+5n\nrLOP5u3TRyZs3TH6WV2nmr3vs6e2HPKTmup0TOH+iFgOfBN4tr4yM28a7Q0RsQx4CzA/IjZQ637q\nK953BXAbtemo66lNSf3IOOrvWKsG+5b7N/KZb63ZpYuj2ZWRI4+wN28fqs0oKXTS2E/0qsrxzBFv\nbKDH0rCPNJYzhcYrLuf193U8Xx52nSbY6cBxvb978JGn+drdj+4WDHvP6eG54V07Wer93wOH7s+l\nK9Zxwo/uLN/3XM8c+l7YyVBPL3u/MExGcOK6u7jqdadMePDSo3VNZ52GwlzgKaDxCo8ERg2FzDyj\n1Q6L32Y4u8PPn5BWR+cwthkrIw3tTC5dsQ6go8a+G1dVNjbQY2nYG411TKFx7venTz6y5VWw/X09\nzO3rHXUcpNPZR/X53PVGfizdc6cuXQhfPx+Gn4eXvIS9X/c6uOwy9j77bBgcpP/ZZ/nzn93Dn3/u\nyo6/M2km6vSGeJUexVet1dF5/flEtGrQR7420XuZTOSoG+CMNxzS8ZjCyK6P0RpcaD37qHGmycjp\nhO0a7PFcuDOui32eeAJWr4a994ZLLoE/+qPaNQjf+x584Qtw0UXw4IPw05/CAQeMbd/SDNLpz3Ee\nDPxv4Nhi1b8B52XmhqoKm0xVH5236v8e2dhP9KrKZkfdPQEvm1vrT283+6jez9xq9lG7KYHNGtxO\njsSndbfJypXwylfCDTfA0qUvru/pgY9/HI47Dk4/HQYH4aSTulenVLFOp6TeDnwd+Eqx6oPABzLz\nHRXW1tR4Zh+Ndp+VyZj22NcbXHr6rwG7d0ONnPlTN9F7mcy0e8JL6r5Jm5Ja7GxVZi5pt24qjCcU\nml09WW+wofmYQgTljckmc/aRJHXDZF/R/FREfBBYViyfQW3geUbo5J4nnTbmrab5TfsuEklqo9Mz\nhUOpjSm8sVh1J3BuZo7vKqgJmA0Xr0nSVJvUM4XMfAQ4ecJVSZKmtZ5ONoqIgyPi5oh4snjcWMxI\nkiTNIh2FAvBlarelOKh4fKtYJ0maRToNhQWZ+eXMHC4e1wALKqxLktQFnYbCUxHxwYjoLR4fZAbN\nPpIkdabTUPgo8DvAT6n9CM7pVHwDO0nS1HP2kSSp1Onso2sjYl7D8n4RcXV1ZUmSuqHT7qOjM3NL\nfSEzNwNLW2wvSZqBOg2FnojYr75Q/GZzp7fIkCTNEJ027H8D3BUR36R2m/3Tgb+orCpJUld0OtB8\nXUQM8uIvr70nM9dWV5YkqRs67T4iM9dm5mXA8waCJM1OHYdCgz+Y9CokSdPCeEIhJr0KSdK00DYU\nIuI1EfG2iNi3WPWuYv0JlVYmSZpyLUMhIs4F/hH4Q+ChiDglMzcUL//PqouTJE2tdrOPfh94bWZu\ni4jFwA0RsTgzv4jdSJI067QLhZ7M3AaQmf8ZEW+hFgyHYihI0qzTbkzhvyJiSX2hCIiTgPnA6L9g\nL0makdqFwoep3S67VPzIzoeBN1dWlSSpK1p2HzUMKjd77c7JL0eS1E3juU6hYxFxQkSsi4j1EXFh\nk9fPjIhNEbGqePxelfVIklqr7E6nEdELfAl4B7ABuDcilje5RcY3MvOcquqQJHWuyjOF1wPrM/Mn\nmfk8cD1wSoWfJ0maoCpDYSHwWMPyhmLdSKdFxIMRcUNEHNJsRxFxVkQMRsTgpk2bqqhVkkTFYwod\n+BawODOPBm4Hrm22UWZemZkDmTmwYMGCKS1QkvYkVYbCRqDxyP/gYl0pM5/KzOeKxb8HXlthPZKk\nNqoMhXuBwyPisIjYC3g/sLxxg4g4sGHxZODhCuuRJLVR2eyjzByOiHOAFUAvcHVmromIS4DBzFwO\nnBsRJwPDwNPAmVXVI0lqLzKz2zWMycDAQA4ODna7DEmaUSJiZWYOtNuu2wPNkqRpxFCQJJUMBUlS\nyVCQJJUMBUlSyVCQJJUMBUlSyVCQJJUMBUlSyVCQJJUMBUlSyVCQJJUMBUlSyVCQJJUMBUlSyVCQ\nJJUMBUlSyVCQJJUMBUlSyVCQJJUMBUlSyVCQJJUMBUlSyVCQJJUMBUlSyVCQJJUMBUlSqdJQiIgT\nImJdRKyPiAubvL53RHyjeP2eiFhcZT2SpNYqC4WI6AW+BJwIHAGcERFHjNjsY8DmzPxl4AvAX1VV\njySpvSrPFF4PrM/Mn2Tm88D1wCkjtjkFuLZ4fgPwtoiICmuSJLVQZSgsBB5rWN5QrGu6TWYOA1uB\nX6qwJklSCzNioDkizoqIwYgY3LRpU7fLkaRZq8pQ2Agc0rB8cLGu6TYRMQd4OfDUyB1l5pWZOZCZ\nAwsWLKioXElSlaFwL3B4RBwWEXsB7weWj9hmOfC7xfPTge9lZlZYkySphTlV7TgzhyPiHGAF0Atc\nnZlrIuISYDAzlwNXAV+JiPXA09SCQ5LUJZWFAkBm3gbcNmLdJxue/wJ4b5U1SJI6NyMGmiVJU8NQ\nkCSVDAVJUslQkCSVDAVJUslQkCSVDAVJUslQkCSVDAVJUslQkCSVDAVJUslQkCSVDAVJUslQkCSV\nDAVJUslQkCSVDAVJUslQkCSVDAVJUslQkCSVDAVJUslQkCSVDAVJUslQkCSVDAVJUslQkCSVDAVJ\nUslQkCSVDAVJUslQkCSVIjO7XcOYRMQm4JEWm8wHfjZF5cxkfk/t+R11xu+pvenwHR2amQvabTTj\nQqGdiBjMzIFu1zHd+T2153fUGb+n9mbSd2T3kSSpZChIkkqzMRSu7HYBM4TfU3t+R53xe2pvxnxH\ns25MQZI0frPxTEGSNE6zMhQi4r0RsSYiXoiIGTHiP1Ui4oSIWBcR6yPiwm7XMx1FxNUR8WREPNTt\nWqariDgkIu6IiLXF/7Xzul3TdBQRcyPiRxHxQPE9fabbNbUzK0MBeAh4D/CDbhcynUREL/Al4ETg\nCOCMiDiiu1VNS9cAJ3S7iGluGPh4Zh4BHAOc7b+lpp4D3pqZvwYsAU6IiGO6XFNLszIUMvPhzFzX\n7TqmodcD6zPzJ5n5PHA9cEqXa5p2MvMHwNPdrmM6y8wnMvO+4vkzwMPAwu5WNf1kzbZisa94TOuB\n3FkZChrVQuCxhuUN+B9ZExQRi4GlwD3drWR6iojeiFgFPAncnpnT+nua0+0Cxisi/gU4oMlLf5aZ\n/zjV9Uh7oojYF7gROD8zf97teqajzNwJLImIecDNEfGrmTltx6tmbChk5tu7XcMMtBE4pGH54GKd\nNGYR0UctEL6WmTd1u57pLjO3RMQd1Marpm0o2H20Z7kXODwiDouIvYD3A8u7XJNmoIgI4Crg4cz8\nfLfrma4iYkFxhkBE9APvAP69u1W1NitDISLeHREbgDcC/xQRK7pd03SQmcPAOcAKagOD/5CZa7pb\n1fQTEcuAu4BXR8SGiPhYt2uaho4FPgS8NSJWFY93druoaehA4I6IeJDaQdntmXlrl2tqySuaJUml\nWXmmIEkaH0NBklQyFCRJJUNBklQyFCRJJUNBs9J0uRtsRGwr/jwoIm6YwH7Oj4h9Jq8yqTmnpGrW\nKe4G+x/ULhTaQG1++BmZubbiz51TXAvSuG5bZu47Cfv+T2AgM3820X1JrXimoNmoo7vBRsQnImJ1\nca/7zxXrlkTE3RHxYETcHBH7tVn/rxHxvyJiEDivuFr8rmK/n234rMX132eIiDMj4qaI+E5E/Dgi\n/rphu8sjYrDx3vsRcS5wELWLoO4o1v234nPui4hvFvcgkibMUNBs1PZusBFxIrWgeENxr/t6w3wd\n8InMPBpYDXyqzXqAvTJzIDP/BvgicHlmHgU80aLGJcD7gKOA90VE/Z5Uf5aZA8DRwG9FxNGZ+bfA\n48BxmXlcRMwHLgbenpm/DgwCf9zRNyO1YShoT/V24MuZuR0gM5+OiJcD8zLz+8U21wJvHm19w76+\n0fD8WGBZ8fwrLT7/u5m5NTN/AawFDi3W/05E3AfcDxxJ7ceQRjqmWH9ncUvm3214vzQhM/YuqVIL\nze4Gu6loQAE+Ocmf9+yI5U4G6p5reL4TmBMRhwF/ArwuMzdHxDXA3CbvDWr30DljPMVKrXimoNmo\n2d1gb8jMJcVjOXA78JH6jJ6I2D8ztwKbI+I3i/18CPj+aOtH+ew7i88D+MAY634ZtYDZGhGvoPaz\nqXXPAC8tnt8NHBsRv1zU/pKI+JUxfpbUlGcKmnUyczgi6neD7QWuHnk32Mz8TkQsAQYj4nngNuAi\nal0xVxRh8RPgI8VbRls/0nnA1yPiE8CYfuwpMx+IiPup3Vr5MWoBU3cl8J2IeLwYVzgTWBYRexev\nX0xtxpU0IU5JlSSV7D6SJJUMBUlSyVCQJJUMBUlSyVCQJJUMBUlSyVCQJJUMBUlS6f8DqZoj/2qM\n+Q8AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7fb4119f05f8>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "Thetas_hard_neq = run_perturbation_trials(solve_neq, X_hard, y_hard)\n",
    "scatter_thetas(Thetas_hard_neq, theta_true=theta_hard_true, ax=0, ay=2)\n",
    "\n",
    "print(\"Residual norm for one of the trials:\")\n",
    "theta_hard_neq_example = np.random.randint(Thetas_hard_neq.shape[1])\n",
    "calc_residual_norm(X_hard, y_hard, theta_hard_neq_example)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Observe that the computed estimates can be relatively far from the true value, even getting the sign completely wrong in the case of the $\\theta_0$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Algorithm 2: QR decomposition\n",
    "\n",
    "A different method for solving an overdetermined systems is to use a tool from linear algebra known as the [_QR decomposition_](https://en.wikipedia.org/wiki/QR_decomposition).\n",
    "\n",
    "Here is how we can use QR. If $X$ has linearly independent columns, then we would first factor the $m \\times n$ matrix $X$ into the product $X = QR$, where $Q$ is an $m \\times n$ orthogonal matrix and $R$ is an invertible $n \\times n$ upper-triangular matrix. (These dimensions assume $m \\geq n$.) That $Q$ is orthogonal means that $Q^T Q = I$, the identity matrix; $R$ being upper-triangular means all of its entries below the main diagonal are zero."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, observe that the normal equations can be transformed if we substitute $X = QR$:\n",
    "\n",
    "$$\n",
    "\\begin{eqnarray}\n",
    "        X^T X \\theta^* & = & X^T y \\\\\n",
    "  R^T Q^T Q R \\theta^* & = & R^T Q^T y \\\\\n",
    "            R \\theta^* & = & Q^T y.\n",
    "\\end{eqnarray}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Lastly, because $R$ is triangular, solving a system is \"easy\" using _(backward) substitution_. Consider the following $3 \\times 3$ example (taken from [here](http://www.purplemath.com/modules/systlin6.htm)):\n",
    "\n",
    "$$\n",
    "\\left[\\begin{matrix}\n",
    "5 &  4 & -1 \\\\\n",
    "  & 10 & -3 \\\\\n",
    "  &    &  1 \\\\\n",
    "\\end{matrix}\\right]\n",
    "\\cdot\n",
    "\\left[\\begin{matrix}\n",
    "  \\theta_0 \\\\\n",
    "  \\theta_1 \\\\\n",
    "  \\theta_2\n",
    "\\end{matrix}\\right]\n",
    "=\n",
    "\\left[\\begin{matrix}\n",
    "  0 \\\\\n",
    "  11 \\\\\n",
    "  3\n",
    "\\end{matrix}\\right].\n",
    "$$\n",
    "\n",
    "Because it is upper-triangular, you can see right away that $1 \\cdot \\theta_2 = 3 \\implies \\theta_2 = 3$. Then, going to the equation above it, $10 \\theta_1 - 3 \\theta_2 = 10 \\theta_1 - 3(3) = 11 \\implies \\theta_1 = 2$. Lastly, $5 \\theta_0 + 4 \\theta_1 - \\theta_2 = 5 \\theta_0 + 4(2) - 3 = 0 \\implies \\theta_0 = -1$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So, to summarize, a different algorithm to solve $X \\theta^* \\approx y$ using QR would look like the following:\n",
    "\n",
    "1. Compute $X = QR$.\n",
    "2. Form the modified right-hand side, $z = Q^T y$.\n",
    "3. Use back-substitution to solve $R \\theta^* = z$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Conditioning.** What about the sensitivity of this algorithm? Given $R$, we only need to solve linear systems involving $R$. Therefore, it's $\\kappa(R)$ that will determine the stability of the algorithm. So if $\\kappa(R)$ is comparable to $\\kappa(X)$, then the algorithm should be as stable as one can expect any algorithm to be."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise 6** (1 point). Use [numpy.linalg.qr()](https://docs.scipy.org/doc/numpy/reference/generated/numpy.linalg.qr.html) to compute the QR decomposition of $X$ (precomputed above as the variable, `X`). Store the $Q$ and $R$ factors in two variables named `Q` and `R`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "nbgrader": {
     "grade": false,
     "grade_id": "qr",
     "locked": false,
     "schema_version": 1,
     "solution": true
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[  1.   0.   0.]\n",
      " [  1.   1.   1.]\n",
      " [  1.   2.   4.]\n",
      " [  1.   3.   9.]\n",
      " [  1.   4.  16.]] \n",
      " ...\n",
      "\n",
      "Q: (10, 3) \n",
      "\n",
      "R: (3, 3) ==\n",
      "[[ -3.16227766 -14.23024947 -90.12491331]\n",
      " [  0.           9.08295106  81.74655956]\n",
      " [  0.           0.          22.97825059]]\n"
     ]
    }
   ],
   "source": [
    "print(X[:5], \"\\n ...\\n\")\n",
    "\n",
    "#\n",
    "# YOUR CODE HERE\n",
    "Q, R = np.linalg.qr(X)\n",
    "\n",
    "# Print the dimensions of your result\n",
    "print(\"Q:\", Q.shape, \"\\n\")\n",
    "print(\"R:\", R.shape, \"==\")\n",
    "print(R)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "nbgrader": {
     "grade": true,
     "grade_id": "qr_test",
     "locked": true,
     "points": 1,
     "schema_version": 1,
     "solution": false
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "(Passed!)\n"
     ]
    }
   ],
   "source": [
    "# Test cell: `qr_test`\n",
    "\n",
    "assert type(Q) is np.ndarray, \"`Q` is not a Numpy array but should be.\"\n",
    "assert type(R) is np.ndarray, \"`R` is not a Numpy array but should be.\"\n",
    "assert Q.shape == (m, n+1), \"`Q` has the wrong shape: it's {} rather than {}.\".format(Q.shape, (m, n+1))\n",
    "assert R.shape == (n+1, n+1), \"`R` has the wrong shape: it's {} rather than {}.\".format(R.shape, (m, n+1))\n",
    "for i in range(R.shape[0]):\n",
    "    for j in range(i):\n",
    "        assert np.isclose(R[i][j], 0.0), \"R[{}][{}] == {} instead of 0!\".format(i, j, R[i][j])\n",
    "        \n",
    "QTQ = Q.T.dot(Q)\n",
    "assert np.isclose(QTQ, np.eye(Q.shape[1])).all(), \"Q^T Q is not nearly the identity matrix, as it should be.\"\n",
    "\n",
    "assert np.isclose(X, Q.dot(R)).all(), \"QR is not sufficiently close in values to X!\"\n",
    "\n",
    "print(\"\\n(Passed!)\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Condition number of $R$.** Let's check the condition number of $R$ empirically, to verify that it is comparable to $\\kappa(X)$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$$\\kappa(X) \\approx 1.07\\times 10^{2}$$"
      ],
      "text/plain": [
       "<IPython.core.display.Math object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$$\\kappa(X^T X) \\approx 1.15\\times 10^{4}$$"
      ],
      "text/plain": [
       "<IPython.core.display.Math object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$$\\kappa(R) \\approx 1.07\\times 10^{2}$$"
      ],
      "text/plain": [
       "<IPython.core.display.Math object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "cond_R = np.linalg.cond(R)\n",
    "\n",
    "show_cond_fancy(cond_X, 'X')\n",
    "show_cond_fancy(cond_XTX, 'X^T X')\n",
    "show_cond_fancy(cond_R, 'R')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise 7** (3 points). Implement a function, `solve_qr(X, y)`, which uses the QR-based algorithm to estimate $\\theta^*$.\n",
    "\n",
    "> To solve the triangular system, use Scipy's specialized function, available as [sp.linalg.solve_triangular()](https://docs.scipy.org/doc/scipy/reference/generated/scipy.linalg.solve_triangular.html)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {
    "nbgrader": {
     "grade": false,
     "grade_id": "solve_qr",
     "locked": false,
     "schema_version": 1,
     "solution": true
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Comparing your QR solution to the true solution:\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style>\n",
       "    .dataframe thead tr:only-child th {\n",
       "        text-align: right;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: left;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>theta_qr</th>\n",
       "      <th>theta_true</th>\n",
       "      <th>error</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>-1.154632e-14</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>3.552714e-15</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>-2.220446e-16</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   theta_qr  theta_true         error\n",
       "0       1.0         1.0 -1.154632e-14\n",
       "1       1.0         1.0  3.552714e-15\n",
       "2       1.0         1.0 -2.220446e-16"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Residual norm:\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "array([ -2.84217094e-14])"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import scipy.linalg\n",
    "\n",
    "def solve_qr(X, y):\n",
    "    #\n",
    "    # YOUR CODE HERE\n",
    "    Q, R = np.linalg.qr(X)\n",
    "    z = Q.T.dot(y)\n",
    "    theta = scipy.linalg.solve_triangular(R, z)\n",
    "    return theta\n",
    "    \n",
    "theta_qr = solve_qr(X, y)\n",
    "\n",
    "print(\"Comparing your QR solution to the true solution:\")\n",
    "show_2vecs_tibble(theta_qr, theta_true, xname='theta_qr', yname='theta_true', error=True)\n",
    "\n",
    "print(\"Residual norm:\")\n",
    "calc_residual_norm(X, y, theta_qr)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {
    "nbgrader": {
     "grade": true,
     "grade_id": "solve_qr_test",
     "locked": true,
     "points": 3,
     "schema_version": 1,
     "solution": false
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "(Passed!)\n"
     ]
    }
   ],
   "source": [
    "# Test cell: `solve_qr_test`\n",
    "import re\n",
    "\n",
    "try:\n",
    "    del np.linalg.lstsq\n",
    "    solve_qr(X, y)\n",
    "except NameError as n:\n",
    "    if re.findall('lstsq', n.args[0]):\n",
    "        print(\"*** Double-check that you did not try to use `lstsq()`. ***\")\n",
    "    raise n\n",
    "except AttributeError as a:\n",
    "    if re.findall('lstsq', a.args[0]):\n",
    "        print(\"*** Double-check that you did not try to use `lstsq()`. ***\")\n",
    "    raise a\n",
    "finally:\n",
    "    np.linalg.lstsq = SAVE_LSTSQ\n",
    "\n",
    "assert np.isclose(theta_qr, theta_true).all(), \"Your QR-based solution should be closer to the true solution.\"\n",
    "\n",
    "print(\"\\n(Passed!)\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Is QR more stable?** Let's run the same perturbation experiments on the \"hard\" regression problem and see the result."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {
    "nbgrader": {
     "grade": false,
     "grade_id": "cell-9e64a01b45964c31",
     "locked": true,
     "schema_version": 1,
     "solution": false
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Sample estimate for one of the trials:\n",
      "- Gramian-based method: theta^T =\n",
      "\t[ 0.99666873  0.9981636   1.00061631  0.99996041  1.00000097  0.99999999\n",
      "  1.        ]\n",
      "- QR-based method: theta^T =\n",
      "\t[ 0.99157198  1.00020562  1.00007936  0.99999549  1.00000009  1.          1.        ]\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAtMAAAEWCAYAAACkORurAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3X2c3HV16PHPyWYJG54WJVfNSgyKhoJBoqmgsfVCtUFF\nmwaRUrUXtcXePkmLaUP1VcHWkjZXW29tVW4FrCIPCkYq1iA3qDWCl8QEw1OUKgKLDyAEhSywSc79\nY34TJrszu7OzOzM7O5/367WvzP7m6cwCJ4fvnu/5RmYiSZIkaeJmtTsASZIkqVNZTEuSJEkNspiW\nJEmSGmQxLUmSJDXIYlqSJElqkMW0JEmS1CCLac1YEXFeRHy63XGMJyLeHBHXtTsOSZI0cRbTalhE\n3B0RP42IAyqu/W5EfLWNYU1rEbEwIjIiZpevZealmfnr7YxLkiQ1xmJak9UDvGuyLxIl/vsoSdpH\nRJwZEdsiYmdE/Dgi/iUiDinuOy8ihiPi0YjYERHfjIiXtTtmdReLF03WWuDdEdFf7c6IeHlE3BwR\njxR/vrzivq9GxAciYiOwE3huce1vioT4aET8e0Q8PSIujYifF6+xsOI1PhwR9xb3bY6IX6k38Ig4\nJSK2ViTgYyvuWxIR346IX0TEFRFxeUT8TXHfmRHxjRGvlRFxZHH7dRGxpYjp3og4r+KhXy/+3FF8\nvpeNfL06fmZ/HREbi9iui4jDivv2j4hPR8TPis90c0Q8o96fhyRNNxFxDvB3wCrgEOAEYCFwXUT0\nFg+7IjMPBA4DbgA+24ZQ1cUspjVZm4CvAu8eeUdEPA24FvjfwNOBDwHXRsTTKx72VuAs4CDgh8W1\n3yquDwDPA24ELgaeBtwBvK/i+TcDxxX3fQb4bETsP17QEbEEuAh4ZxHbx4FrImJOROwHrAM+Vbzu\nZ4FTx3vNCo8BvwP0A68D/mdErCju+9Xiz/7MPDAzbxwRVz0/s98G3gb8N2A/nvrZ/w9Kf9kcXjz3\n94GhCcQtSdNGRBwMnA/8cWZ+OTOHM/Nu4E3Acynlwr0ycxdwKTAQEfNaHa+6V8cV0xFxUdGne+sU\nvd7fRcStxdfpE3jeURFxY0Q8ERGjCsku81fAH1dJXq8DvpeZn8rMXZl5GXAn8PqKx1ySmbcV9w8X\n1y7OzP/KzEeA/wD+KzOvLxLlZ4El5Sdn5qcz82fF8z8IzAEW1RHzWcDHM/Nbmbk7Mz8JPEFp1eME\noBf4xyJ5f45S0V6XzPxqZm7LzD2Z+R3gMuCVdT69np/ZxZn53cwcAq6k9D8TAMOUiugji8+0OTN/\nXm/cUjOYszUJLwf2B66uvJiZjwJfAvbZa1IshPwO8DPg4RbFKHVeMQ1cApw8FS8UEa8DXkypGDme\nUrvCwVUed3eVpz8E/Anwv6Yilk6WmbcCXwRWj7hrPk+tNpf9kNKKc9m9VV7yJxW3h6p8f2D5m4h4\nd0TcUbRE7KC0MntYHWE/BzinaIfYUTz38CLm+cBgZuaIuOsSEcdHxA0R8UBEPEJphbiemKC+n9mP\nK27v5Kmfx6eA9cDlEXF/RPx9xa9BpXa5BHO2GnMY8GCxkDLSj4DyAs6bihw+BPwe8MYaz5GaouOK\n6cz8OqWkuFdEPC8ivlz0zP5nRBxV58sdDXy9WAF8DPgOdSb9zPxpZt5MaTVQpdaL32Pfou9+SkVr\npQXAYMX3SYOK/ug/p/Qrv0Mzsx94BIg6nn4v8IHM7K/4mlusBP+I0q8JK19nQcXtx4C5FXE8c8Rr\nfwa4Bjg8Mw8BPlYR03ift56fWVXFKvr5mXk0pRWdUyit0khtY87WJDwIHBYV048qPKu4H+DKIv8/\nA7gVeEmL4pOADiyma7iQUk/VSyj1j/5Lnc+7BTg5IuYWm7hOpLQ6qQnKzLuAKyit/JR9CXhBRPx2\nRMwufiV7NKVV7KlwELALeACYHRF/BYxaparh/wC/X6wiR0QcUGwcPIhSj/Yu4E8iojciVgIvrXju\nLcAxEXFc0Z99XpW4HsrMxyPipezb1/cAsIdSv181Df/MIuLEiFgcET3AzykVDXvGe57UBuZs1eNG\nSu13KysvRsSBwGso7dfZKzMfpNTCd15EPKtFMUqdX0wX/1G9nNLGs62UNpI9q7hvZUVvXeXXeoDM\nvI5S8fJNSn2tNwK7i+f+c5QmPWwF5pdvR8R7Wv8pO8b7gb0zpzPzZ5RWR8+h1MP258ApRcKbCuuB\nLwPfpdQK8TjV20ZGycxNlFbSP0Kpt+4u4MzivicpJe8zKa2onU5Fz15mfpfSZ70e+B6wz2QP4A+A\n90fELyj1k19Z8dydwAeAjUV7yQkj4prMz+yZwOcoFdJ3AF+j1PohTRvmbNWr2DdzPvBPEXFysbix\nkFJOfZDSZsORz9lO6e+GP29hqOpysW9baGco/mP6Yma+sOiX256Zk/6/0Ij4DPDpzPzSiOt3Z+bC\nGs85D3g0M+3Dm8Ei4hLgvsx8b7tjkTqNOVuTERHvAP4UOJLSJvOvAb+dmfcX/zyPzMy3VDz+eGAD\ncERm/rQNIavLdPzKdDGt4AcRcRrsPfzjRfU8NyJ6ohg5FqUZw8cCHussSU1iztZEZeYnMvOFmbk/\n8HZKI1NnF/edV1lIF9e+lZkHWEirVao19U9rEXEZ8N8pbUq4j9LGtzcDH42I91IaaXY5pd668fQC\n/1nsM/s58JZ6dwAXm842UerR3RMRZwNHO4pMkp5iztZUysyLI2IXpVahe9odjwQd2uYhSZIkTQcd\n3+YhSZIktUtHtXkcdthhuXDhwnaHIUkTtnnz5gczs6uOODZnS+pk9ebtjiqmFy5cyKZNm9odhiRN\nWETUfYrmTGHOltTJ6s3btnlIkiRJDbKYliRJkhpkMS1JkiQ1yGJakiRJapDFtCRJktQgi2lJkiSp\nQRbTkiRJUoMspiVJkqQGWUxLkiRJDbKYliRJkhpkMS1JkiQ1yGJakiRJapDFtCRJktSg2e0OQJLa\nZd2WQdau3879O4aY39/HquWLWLFkoN1hSZI6iMW0pK60bssg5169jaHh3QAM7hji3Ku3AVhQS5Lq\nZpuHpK60dv32vYV02dDwbtau396miCRJnchiWlJXun/H0ISuS5JUjcW0pK40v79vQtclSaqmbcV0\nRBweETdExO0RcVtEvKtdsUjqPquWL6Kvt2efa329PaxavqhNEU1/5m1JGq2dGxB3Aedk5rcj4iBg\nc0R8JTNvb2NMkrpEeZOh0zwmxLwtSSO0rZjOzB8BPypu/yIi7gAGAJOypJZYsWTA4nkCzNuSNNq0\n6JmOiIXAEuBb7Y1EklQP87YklbS9mI6IA4GrgLMz8+dV7j8rIjZFxKYHHnig9QFKkvYxVt42Z0vq\nNm0tpiOil1JCvjQzr672mMy8MDOXZubSefPmtTZASdI+xsvb5mxJ3aad0zwC+ARwR2Z+qF1xSJLq\nY96WpNHauTK9DHgrcFJEbC2+XtvGeCRJYzNvS9II7Zzm8Q0g2vX+kqSJMW9L0mht34AoSZIkdSqL\naUmSJKlBFtOSJElSgyymJUmSpAZZTEuSJEkNspiWJEmSGtS20XiSusu6LYOsXb+d+3cMMb+/j1XL\nF7FiyUC7w5IkaVIspiU13botg5x79TaGhncDMLhjiHOv3gZgQS1J6mi2eUhqurXrt+8tpMuGhnez\ndv32NkUkSdLUsJiW1HT37xia0HVJkjqFxbSkppvf3zeh65IkdQqLaUlNt2r5Ivp6e/a51tfbw6rl\ni9oUkSRJU8MNiJKarrzJ0GkekqSZxmJaUkusWDJg8SxJmnFs85AkSZIaZDEtSZIkNchiWpIkSWqQ\nPdOSpoxHhkuSuo3FtKQp4ZHhkqRuZJuHpCnhkeGSpG5kMS1pSnhkuCSpG1lMS5oSHhkuSepGFtOS\n6rZuyyDL1mzgiNXXsmzNBtZtGdx7n0eGS5K6kRsQJdVlvA2GHhkuSepGFtOS6jLWBsNyweyR4ZKk\nbmMxLaku420wdMa0JKkb2TMtqS5jbTAst4AM7hgieaoFpLKnWpKkmchiWlJdam0wPPGoeZxz5S3O\nmJYkdSWLaUlVjZzcAXDBysUM9PcRwEB/H6e+ZICrNg+yO7PqazhjWpI001lMSxqlVtsGwMbVJ/EP\npx8HwKdvumfUinQlZ0xLkmY6NyBKGmW8o8ErR+TV4oxpSVI3sJiWNGoSx+AYkzuqFdoj9URwwcrF\nTvOQJM14FtNSl6t2GEsA1bqg5/f3jdsH3dfbYyEtSeoa9kxLXa7aSnMCMeJx5baNsfqgB/r7LKQl\nSV3FlWmpy9VaaU5K7Rq7MxkYcQjLyJ5pV6MlSd3KYlqagSZyGuFYPdK7M/euSFceGQ542qEkSUBk\njfmw09HSpUtz06ZN7Q5DmpbKBXS1nucA3nzCApY+52mjimAYfzrHQH8fG1ef1NT4Z7qI2JyZS9sd\nRyuZsyV1snrzdlt7piPiooj4aUTc2s44pE5XORcaRm8eTEozoVd97paqs6PLh7HU4uErAnO2JFXT\n7g2IlwAntzkGqePVM64OYHj3vmV2eXb0iiUDbFx9Us2C2sNXVLgEc7Yk7aOtxXRmfh14qJ0xSDPB\nZFaOK5+7avki+np79rnfw1dUZs6WpNHavTItaQpMZuW48rkrlgzsbfkIHHUnSdJ4pv00j4g4CzgL\nYMGCBW2ORpqeVi1fxNlXbB3zMb2zAmLfVo9qq84rlgxYPKth5mxJ3Wbar0xn5oWZuTQzl86bN6/d\n4UjT0oolAxw6t7fm/QGsPe1FrH3ji1x1VlOZsyV1m2m/Mi2pPu97/TE1V6f75/aOmhMtSZImr92j\n8S4DbgQWRcR9EfGOdsYjTRfrtgyybM0Gjlh9LcvWbGDdlsFxnzNWkbxj5/BUhqcuZc6WpNHaujKd\nmWe08/2l6ag8M7o86q48D3rTDx/ihjsfGPPUwYEapxk62k5TwZwtSaNN+55pqdtUmxk9NLybS2+6\nZ9SBKyNXrB1tJ0lSa1lMS9NMrZnRI081LB+4UsnRdpIktZYbEKVpZN2WQWZFsDtHls7VDe4YYt2W\nwX2KZUfbSZLUOq5MS9NEuVe63kK6rFq7hyRJag2LaWmaqNYrXY9q7R6SJKk1LKalaaJWr3SznytJ\nkhpXVzEdEc+IiE9ExH8U3x/tfFFpak1mfJ2j71TJnC1JrVPvyvQlwHpgfvH9d4GzmxGQNFNM9OCV\namPtRurv63X0nepxCeZsSWqJeovpwzLzSmAPQGbuAibe3Cl1ifJmwvHmQleqHGsHECPu7+vt4bw3\nHOPoO9XDnC1JLVLvaLzHIuLpFKNuI+IE4JGmRSV1uFoHr5z/77exdv32vacYnnjUPG648wEGdwzR\nU4zEG+jv4x9PP27v61Q78dDiWeMwZ0tSi9RbTP8ZcA3wvIjYCMwDTmtaVFKHq7Uh8OGdwzy8cxgo\nrVZ/+qZ79t5XHolXXsW+YOViNq4+qfnBaiYyZ0tSi9RbTN8GvBJYROm3z9txEohU1UQPXqmmPO7O\nFWg1yJwtSS1Sb3K9MTN3ZeZtmXlrZg4DNzYzMKkTNXrwSjWOu9MkmLMlqUXGXJmOiGcCA0BfRCzh\nqT1RBwNzmxyb1HFqHbzSE8Gc2cHO4T11v5bj7jRR5mxJar3x2jyWA2cCzwY+VHH9F8BfNikmqWPV\nWk3ek8l+s2fXXUw77k4NMmdLUouNWUxn5ieBT0bEqZl5VYtikjrW/P4+BqsU1PP7++pu2xgYMblD\nqpc5W5Jar64NiJl5VUS8DjgG2L/i+vubFZg03a3bMjhqdN2q5YtY9blbGN79VM90b0+wavki1q7f\nXrXQrlRekbaQ1mSYsyWpdeo9TvxjwOnAH1PqwTsNeE4T45KmtVqHsmz64UPFZN+nDO9ONv3wobpO\nOCxP8ZAmw5wtSa1T7zSPl2fm7wAPZ+b5wMuAFzQvLGl6q3Uoy2XfupfhPaMneVxazJOuPL2wFqd4\naAqYsyWpReotpst/u++MiPnAMPCs5oQkTX+1Ct5aI/ES9s6N3rj6JH6w5nV7jw0fySkemgLmbElq\nkXqL6S9GRD+wFvg2cDdwWbOCkqa7RgrekQV4tbYPp3hoipizJalF6t2A+NfFzasi4ovA/pn5SPPC\nkqa3E4+ax6U33TOyPXpMIwvw8ibDkZsY3XyoyTJnS1Lr1HucOBHxcmBh+TkRQWb+W5Piktqq2qSO\ncpG7bssgV20enFAhXWvFecWSAYtnNYU5W5Jao65iOiI+BTwP2AqUd10lYGLWjFOe1FHeYFie1FF2\nzpW3TOi48J4ILli52KJZLWPOlqTWqXdleilwdOYEKgipQ9Wa1PFnV26lyqCOMQXwwTe9yEJarWbO\nlqQWqXcD4q3AM5sZiDRd1D4SfGKvE8CbT1hgIa12MGdLUovUuzJ9GHB7RPw/4Inyxcx8Q1Oiktpk\n3ZZBZkVMqI2jUlD6XbpHgqvNzNmS1CL1FtPnNTMIqRnG2kRY7bHn//ttPLxzuOH3s4DWNHJeuwOQ\npG5R72i8rzU7EGkqjbWJcGSxO/KxE9HX2+PmQk075mxJap0xi+mI+EZmviIifgH7TAILIDPz4KZG\nJzWo1ibC8imE8NTK9WAdx3f39faMer3+vl7Oe8MxFtKaNszZktR6YxbTmfmK4s+DWhOONDVqbSIc\n3DHEsjUbGNwxtLe/eTzl9g0PV9F0Z86WpNYbb2X6aWPdn5kPTW040tSY399Xc8W5fL2eQrp82IqH\nq6gTmLMlqfXG65neTKnmCGAB8HBxux+4BziiqdFJE1C54fCQvl56e4Lh3Y2P2bWNQx3InC1JLTbm\nnOnMPCIznwtcD7w+Mw/LzKcDpwDXtSJAqR7lTYSDO4ZIYMfQMCTMIlm+/ZswwVF3FtLqROZsdbxM\nuPrqCedsqZ3qPbTlhMz8UvmbzPwP4OXNCUka27otgyxbs4EjVl/LsjUb9q5Ij9wgOLwn+aWffJ+P\nr/tbjv7pDyb0HjuGhjn36m2s2zI4laFLrWLOVmfauhVOPRVuuaXdkUh1q7eYvj8i3hsRC4uv9wD3\nNzMwqZqRK9DlkXe1+qNPuf1rJPD6O74+4fcqT/+QOpA5W53psssgAi6/vN2RSHWrt5g+A5gHfB64\nurh9RrOCkmqpNfKuJ2L0gzP5zdu/RgArbr+hoferNRVEmubM2eo8mfCZz5T+vPTSdkcj1W3cQ1si\nogf4y8x811S/eUScDHwY6AH+NTPXTPV7aGapVdzuzhw1C/r5D97DwY8/CsAhQ49y5IP3cNdhCyb0\nfvP7+xoPVmqDZubs4vXN25pyx77vyzzjvu/zhZ88yFxg548f4A2/+1H+67AF/GDN69odnjSmcYvp\nzNwdEa+Y6jcuEv4/A68G7gNujohrMvP2qX4vzRxVR95lsuLB2znj6Kexbuv9PPTYEwC8+ns30ZOl\n4rond/POb13F9c8/fp+nPtbbxzcWHlf6tWIVJx41b+o/hNREzcrZYN7WFMmE66+HR0uLHf/zU5t4\nGbVz9jtXfpOPv3Vp6bkHHgivelXNnC21Q13HiQNbIuIa4LPAY+WLmXn1JN77pcBdmfl9gIi4HPgN\nwKSsmlYtXzTq6O+De5IPrP9nDrjoHl489wAe31XaBR6ZzNm9C4A5u3fxmu3f5OTv3ljct4cDhx/n\nnkOewa/93scY7umt+n433PlAkz+R1BTNyNlg3tZUGB6Gd74TfvADnth/Ln+/p3R5rJw9vD7o3fkY\nHHEE3Hkn7Ldfu6KXRqm3mN4f+BlwUsW1pNSL16gB4N6K7+8Djq/xWHWxyvnR8/v7OPUlA9xw5wMM\n7hiiJ4Kf75nFb/zev3DJxo/z9K9ex0HDT1R9nQOHn1rR3jl7Dtcc9Susfs2f1CykwZ5pdaxm5Gww\nb2sq7LcffOc78I53sPvqL3DQrvFz9lDOoff00+Ff/9VCWtNOXcV0Zr6t2YHUEhFnAWcBLFgwsX5X\ndb7y9I7ySvTgjiGu2jzIqS8Z4KrNg3uv3zUErz7+D/j1/V7ABdd9hP12DTM794x6vV0xiydn78fq\n5X/INcecOO772zOtTmTO1rR34IFwxRWce8q768jZvZy7/I/48OX/qw2BSuOra5pHRDw7Ij4fET8t\nvq6KiGdP8r0HgcMrvn92cW0fmXlhZi7NzKXz5tm/2m1qTe+47Fv3Vr3+hReeyMlv+wjDPdX/P3G4\nZzYnv+2fuPaFJ1W9v1L5KHGp0zQpZ0MdeducrYn44uKT6sjZH2HTste2ODKpfvWOxrsYuAaYX3z9\ne3FtMm4Gnh8RR0TEfsBvFe8h7TXW9I5ads/qodbWlMhkV09P1ef3zgoOndtLAAP9fVywcrEnIKpT\nNSNng3lbU+yM4w+vK2e7sKHprN6e6XmZWZmIL4mIsyfzxpm5KyL+CFhPacTSRZl522ReUzNP1ekd\nY+jv6+W1N2+kXCo/MWs2vXt2Mzyrhzl7dpERnLz9m1z0yysA6IlgTybz+/tYtXyRxbNmiinP2WDe\n1tT7mxWL+dIXPjlmzv743B+y2NysaazelemfRcRbIqKn+HoLpc0tk5KZX8rMF2Tm8zLzA5N9Pc08\nq5Yvoq+3p67Hzgo45UXP4tRt/5e+XU/yWO8cvj3wSyx/+0fYMnAUj/XOoW/Xk5y27fq9z9mTyQ/W\nvI6Nq0+ykNZM0pScDeZtTb3Xbr6Ovl1PwgEHMOdXX8GsW7cx51eWwQEH0LfrSRZ/5fPtDlEaU73F\n9NuBNwE/Lr7eCLRtg4u6x4olA1ywcjED/X01fw1Ytifhpm/cyqIH7uaJnl4+9Iq3cMYZf8v35j2H\nt/z2BfzDK97MEz2zOeqBu5n36MOAGww1Y5mz1Rl+9CPYtg3mzIH3vx82bIBjjin9ef75T03++PGP\n2x2pVFO90zx+CLyhybFIVa1YMrB31XjZmg1jtn0suPtO7ul/Jn+w4lxue8bz9l7fRXDpstO4ccGx\nfHTdBSz+8fe48ZdeZh+eZiRztjrG5s3w3OfC5z4HS5Y8dX3WLDjnHDjxRHjjG2HTJjjllPbFKY2h\nrmK62AX+T8Cy4tJ/Au/KzPuaFZi6w8gZ0mP1La/bMsjOJ3eN+XobjnwpG4586ajrA8Vrr12/H6/8\n/U8wv7+PC+yR1gxlzlbHOOWUsYvkF78Yvv/91sUjNaDeDYgXA58BTiu+f0tx7dXNCErdodoM6XOv\n3gYwqsgd+dhaDp3by+PDe/Z5XHnEXeUKtzTDmbMlqUXq7Zmel5kXZ+au4usSwAGimpRaM6TXrt++\nz7V1WwY558pbxi2k+3p7eN/rj9mnx9oRd+pS5mxJapF6V6Z/VuwGv6z4/gymaGe4uletGdKV18sr\n0mPNlQ4Y1SJi8awuZ86WpBapt5h+O6X+u38AEvgm7gzXJNWaIV05YaPa6nWlgf4+Nq4e/zRDqcuY\nsyWpRZzmobZZtXzRqD7okUd411q9rvZYSSXmbElqnbp6piPikxHRX/H9oRFxUfPCUjcoz5A+dG7v\n3mtzZu/7r2StOdA9EfZCSzWYsyWpderdgHhsZu4of5OZDwNLxni8VLfHh/fsvb1jaJhzr97Gui2D\nQPUTEPt6e/jgm15kIS3VZs6WpBapt2d6VkQcWiRkIuJpE3iuVFV5SsfIzYVDw7s558pbgKc2EtY7\ni1oSYM6WpJapN7l+ELgxIj5LaXjCG4EPNC0qzXjjTenYnbnPzGmLZ2lCzNmS1CJ1tXlk5r8BK4Gf\nAD8GVmbmp5oZmGa28aZ0QPWZ05LGZ86WpNapt2eazLw9Mz8CPJmZtzcxJnWBsaZ0NPI4SfsyZ0tS\na9RdTFf4/SmPQjPWui2DLFuzgSNWX8uyNRv2biysNaVjpHofJ6kmc7YkNVEjG1JiyqNQx1u3ZXDU\nJkFgnznSgzuG9vZBV5sxPVIUj5M0KeZsSWqiMVemI+L4iDi4uN0XEe+ndEzt30XEIS2JUNNeeTPh\n4I4hklLRfPYVWzn7iq2jiuVyH3R5xvTAGCvPiceCSxNhzpak1huvzeMiYGdx+8PAwcB5xbWLmxeW\nOkk9mwkrlfugVywZYOPqk2oW1GMV2pKqMmdLUouN1+YxKzN3FbeXZuaLi9vfiIitTYxLHWSimwRH\n9kHXc6y4pLqYsyWpxcZbmb41It5W3L4lIpYCRMQLgOGmRqaWqrVRsB4T2SRYrUiubPkISivSHhUu\nNcScLUktNt7K9O8CH46I9wIPUjoE4F7g3uI+zQDlnudqGwXrKWhPPGoen77pnnEf1xNRs0j2YBZp\nSpizJanFxiymM/MR4MxiQ8sRxePvy8yftCI4tUa1nufKjYJl1SZ2rFgywA13PjDue/T19rjaLDWZ\nOVuSWq+u0XiZ+XPglibHojap1fM8uGOIZWs2cP+OIfrn9vLo47sY3pN77yuvXo/XMz1QUXhLaj5z\ntiS1TiNzpjXDzO/vY7BKQRyw9/rDO0e3W5ZXr2s9f6C/j42rT5ryeCVJkqaLRk5A1Ayzavki+np7\n9rkWlOY8j+f+HUNVn+80DkmS1A1cmdbe9ovKfuhqK83VzO/vq/p82zokSVI3sJgWMHqaxrI1G8Yt\nqCtXn53GIUmSupFtHqqqWutGb0/Q39frLGhJkqSCK9NdptZ4u5Fs3ZAkSRqfxXQXmejhLLZuSJIk\njc02jy4y1uEskiRJmjiL6S5Sa0PheIeuSJIkqTqL6S6xbssgUeO++f19LY1FkiRpprCY7hJr12+v\neghLgIerSJIkNchiukvUauVIqm8+lCRJ0vgsprtErVaOngiOWH0ty9ZsYN2WwRZHJUmS1NksprtE\ntUNYAHZnkjw1Js+CWpIkqX5tKaYj4rSIuC0i9kTE0nbE0G1WLBnggpWLGejvIyitSI/kmDxJtZi3\nJam6dh3aciuwEvh4m96/o1WeYtg/t5dMeGRoeNxTCisPYTli9bVVH+OYPEk1mLclqYq2FNOZeQdA\nVFkd1dhGnmL48M7hvfeNdaLhyGPED+nrZcfQMCM5Jk9SNeZtSarOnukOU+0Uw0rVWjXKBfjgjqG9\n/dGPPbmL3ln7/qXY19vjmDxJkqQJaNrKdERcDzyzyl3vycwvTOB1zgLOAliwYMEURde56mnDGPmY\nagX48O5kpUscAAAKtklEQVTk0Lm9zN1v9t7V6rFaRCTNfFORt83ZkrpN04rpzHzVFL3OhcCFAEuX\nLq127khXmd/fV/NY8MrHVKpVgO/YOcyWv/r1KYtNUmebirxtzpbUbWzz6DCrli+it6d2z2K1Vo1a\nfdD2R0uSJE1Ou0bj/WZE3Ae8DLg2Ita3I46OVWOtZ6C/jwtWLh7VqlFtxrT90ZImwrwtSdW1pZjO\nzM9n5rMzc05mPiMzl7cjjk60dv12hveMrqb7+3oB+NMrto46zXDkjOlaRbck1WLelqTq2jVnWg2q\n2f88NLx31F21EXmVM6YlSZI0NeyZ7jD19jl7mqEkSVLzWUx3mGr9z7V4mqEkSVJz2ebRYcqtGpWn\nGe58ctc+JyGWOa1DkiSpuSymO9DI/ueRR4yD0zokSZJawWJ6Bqi2Wu1phpIkSc1nMT1DOK1DkiSp\n9Sym22jdlkFXkyVJkjqYxXSbjOxzrjYbWpIkSdObxfQETOVK8tr12/fZMAhPzYa2mJYkSeoMFtN1\nmuqV5FozoJ0NLUmS1Dk8tKVOY60kN6LWDGhnQ0uSJHUOi+k6TfVKcrWTDJ0NLUmS1Fkspus01SvJ\nK5YMcMHKxQz09xHAQH8fF6xcbL+0JElSB7Fnuk6rli+a8lMGnQ0tSZLU2Sym69TIKYPOkZYkSZrZ\nLKYnoNZKcrWiGXCOtCRJ0gxnMT1JtUbmzZk9yznSkiRJM5zFdJ1qtWzUGpk38lqZc6QlSZJmDovp\nOox1YMtEi2PnSEuSJM0cjsarw1gHttQqjg+d2+scaUmSpBnOYroOYx3YUuvwlfe9/hjnSEuSJM1w\ntnnUYX5/H4NVCur5/X3jjsyzeJYkSZq5XJmuw1hHfztLWpIkqXu5Ml2HWqvP4CxpSZKkbmYxXafK\nA1vKq9HVWj+cJS1JktQ9LKYL9bZrjByTV42zpCVJkrqDxTRjz5EeWVBXG5M3krOkJUmSuoMbEBl7\njvRI4606O0takiSpe1hMM/Yc6ZHGWnV2lrQkSVJ3sc2DsedIj+ylPvGoeVy1eXCfley+3h6LaEmS\npC7kyjS150ifeNQ8zr16G4M7hkhKvdRXbR7k1JcMeLKhJEmSXJmG2nOka/VS33DnA2xcfVI7QpUk\nSdI0YjFdqJwjXfanV2yt+lhH30mSJAls8xhTrc2Gjr6TJEkSWEyPqVYvtaPvJEmSBLZ5jKlWL7Wb\nDSVJkgRtKqYjYi3weuBJ4L+At2XmjnbEMp5qvdSS1G06KW9LUiu1q83jK8ALM/NY4LvAuW2KQ5JU\nH/O2JFXRlmI6M6/LzF3FtzcBz25HHJKk+pi3Jam66bAB8e3Af7Q7CElS3czbklRoWs90RFwPPLPK\nXe/JzC8Uj3kPsAu4dIzXOQs4C2DBggVNiFSSBFOTt83ZkrpN04rpzHzVWPdHxJnAKcCvZWaO8ToX\nAhcCLF26tObjJEmTMxV525wtqdu0a5rHycCfA6/MzJ3tiEGSVD/ztiRV166e6Y8ABwFfiYitEfGx\nNsUhSaqPeVuSqmjLynRmHtmO95UkNca8LUnVTYdpHpIkSVJHspiWJEmSGmQxLUmSJDXIYlqSJElq\nkMW0JEmS1KAY47yUaSciHgB+OEUvdxjw4BS9Vqfoxs8M3fm5u/Ezw/T+3M/JzHntDqKVzNlTohs/\nt5+5e0z3z11X3u6oYnoqRcSmzFza7jhaqRs/M3Tn5+7Gzwzd+7m7Qbf+s+3Gz+1n7h4z5XPb5iFJ\nkiQ1yGJakiRJalA3F9MXtjuANujGzwzd+bm78TND937ubtCt/2y78XP7mbvHjPjcXdszLUmSJE1W\nN69MS5IkSZNiMS1JkiQ1qKuL6Yg4LSJui4g9EdHxo1nGEhEnR8T2iLgrIla3O55WiIiLIuKnEXFr\nu2NplYg4PCJuiIjbi3+339XumJotIvaPiP8XEbcUn/n8dsek5jBnz2zmbHN2p+rqYhq4FVgJfL3d\ngTRTRPQA/wy8BjgaOCMijm5vVC1xCXByu4NosV3AOZl5NHAC8Idd8M/6CeCkzHwRcBxwckSc0OaY\n1Bzm7JntEszZ5uwO1NXFdGbekZnb2x1HC7wUuCszv5+ZTwKXA7/R5piaLjO/DjzU7jhaKTN/lJnf\nLm7/ArgDGGhvVM2VJY8W3/YWX+6snoHM2TObOduc3am6upjuIgPAvRXf38cM/49VEBELgSXAt9ob\nSfNFRE9EbAV+CnwlM2f8Z9aMZs7uQubszjW73QE0W0RcDzyzyl3vycwvtDoeqRUi4kDgKuDszPx5\nu+NptszcDRwXEf3A5yPihZnZNX2XM4k5W93InN3ZOXvGF9OZ+ap2xzANDAKHV3z/7OKaZqCI6KWU\nlC/NzKvbHU8rZeaOiLiBUt9lxybmbmbOBszZXcWc3fk52zaP7nAz8PyIOCIi9gN+C7imzTGpCSIi\ngE8Ad2Tmh9odTytExLxidYOI6ANeDdzZ3qikSTFndwlz9szI2V1dTEfEb0bEfcDLgGsjYn27Y2qG\nzNwF/BGwntLmhisz87b2RtV8EXEZcCOwKCLui4h3tDumFlgGvBU4KSK2Fl+vbXdQTfYs4IaI+A6l\nIuQrmfnFNsekJjBnz2zmbHN2p/I4cUmSJKlBXb0yLUmSJE2GxbQkSZLUIItpSZIkqUEW05IkSVKD\nLKYlSZKkBllMa9qJiJMjYntE3BURq9sYx6PFn/Mj4nOTeJ2zI2Lu1EUmSdOHOVvdztF4mlYiogf4\nLqUh7vdRmkF5Rmbe3uT3nV3Mdq289mhmHjgFr303sDQzH5zsa0nSdGLOllyZ1vTzUuCuzPx+Zj4J\nXA78xsgHRcRfRMS2iLglItYU146LiJsi4jsR8fmIOHSc61+NiH+MiE3Au4rTxm4sXvdvKt5rYUTc\nWtw+MyKujogvR8T3IuLvKx730YjYFBG3RcT5xbU/AeZTGlB/Q3Ht14v3+XZEfDYiJp38JalNzNnq\nehbTmm4GgHsrvr+vuLZXRLyGUrI+PjNfBJST478Bf5GZxwLbgPeNcx1gv8xcmpkfBD4MfDQzFwM/\nGiPG44DTgcXA6RFxeHH9PZm5FDgWeGVEHJuZ/xu4HzgxM0+MiMOA9wKvyswXA5uAP6vrJyNJ0485\nW13PYlqd6FXAxZm5EyAzH4qIQ4D+zPxa8ZhPAr9a63rFa11RcXsZcFlx+1NjvP//zcxHMvNx4Hbg\nOcX1N0XEt4EtwDHA0VWee0JxfWNEbAX+R8XzJWkmMmdrRpvd7gCkEQaBwyu+fzbwQJHEAP5qit/v\nsRHf17OJ4ImK27uB2RFxBPBu4Jcz8+GIuATYv8pzA/hKZp7RSLCSNM2Ys9X1XJnWdHMz8PyiF24/\n4LeAz2XmccXXNcBXgLeVd1tHxNMy8xHg4Yj4leJ13gp8rdb1Gu+9sXg/gDdPMO6DKSX5RyLiGcBr\nKu77BXBQcfsmYFlEHFnEfkBEvGCC7yVJ04U5W13PlWlNK5m5KyL+CFgP9AAXZeZtIx7z5Yg4DtgU\nEU8CXwL+ktKv3z5WJOzvA28rnlLr+kjvAj4TEX8BfGGCcd8SEVuAOyn1D26suPtC4MsRcX/Rg3cm\ncFlEzCnufy+l3fCS1FHM2ZKj8SRJkqSG2eYhSZIkNchiWpIkSWqQxbQkSZLUIItpSZIkqUEW05Ik\nSVKDLKYlSZKkBllMS5IkSQ36/7gZm1EY+CN+AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7fb410d44da0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "Thetas_hard_qr = run_perturbation_trials(solve_qr, X_hard, y_hard)\n",
    "\n",
    "# Plot side-by-side against normal equations method\n",
    "def compare_scatter_thetas(T0, title0, T1, title1, ax=0, ay=1, **kwargs):\n",
    "    xmin, xmax = calc_lims(np.array([Thetas_hard_neq[ax, :], Thetas_hard_qr[ax, :]]))\n",
    "    ymin, ymax = calc_lims(np.array([Thetas_hard_neq[ay, :], Thetas_hard_qr[ay, :]]))\n",
    "    xylim = [xmin, xmax, ymin, ymax]\n",
    "    figure(figsize=(12, 4))\n",
    "    subplot(1, 2, 1)\n",
    "    scatter_thetas(T0, title=title0, ax=ax, ay=ay, xylim=xylim, **kwargs)\n",
    "    subplot(1, 2 ,2)\n",
    "    scatter_thetas(T1, title=title1, ax=ax, ay=ay, xylim=xylim, **kwargs)\n",
    "    \n",
    "compare_scatter_thetas(Thetas_hard_neq, 'Normal equations',\n",
    "                       Thetas_hard_qr, 'QR',\n",
    "                       ax=0, ay=-1, theta_true=theta_hard_true)\n",
    "\n",
    "print(\"Sample estimate for one of the trials:\")\n",
    "theta_hard_neq_example = Thetas_hard_neq[:, np.random.randint(Thetas_hard_neq.shape[1])]\n",
    "theta_hard_qr_example = Thetas_hard_qr[:, np.random.randint(Thetas_hard_qr.shape[1])]\n",
    "msg = \"- {}-based method: theta^T =\\n\\t{}\"\n",
    "print(msg.format(\"Gramian\", theta_hard_neq_example.T))\n",
    "print(msg.format(\"QR\", theta_hard_qr_example.T))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You should observe that the QR-based method does, indeed, produce estimates much closer to the true value despite the problem's high condition number."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Performance tradeoff.** Although QR produces more reliable results, there can be a performance tradeoff, as the following quick test should show."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {
    "nbgrader": {
     "grade": false,
     "grade_id": "cell-e40dba0f6f86d187",
     "locked": true,
     "schema_version": 1,
     "solution": false
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "=== Performance of the normal equations-based algorithm ===\n",
      "95.2 µs ± 7.92 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each)\n",
      "\n",
      "=== Performance of the QR-based algorithm ===\n",
      "144 µs ± 293 ns per loop (mean ± std. dev. of 7 runs, 10000 loops each)\n"
     ]
    }
   ],
   "source": [
    "print(\"=== Performance of the normal equations-based algorithm ===\")\n",
    "%timeit solve_neq(X_hard, y_hard)\n",
    "\n",
    "print(\"\\n=== Performance of the QR-based algorithm ===\")\n",
    "%timeit solve_qr(X_hard, y_hard)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Summary comment.** The intent of this notebook was to help you appreciate some of the reliability and performance issues involved in the design of numerical algorithms. The key ideas are as follows.\n",
    "\n",
    "1. The _conditioning_ of a problem, which is a measure of how hard it is to solve. For solving linear systems, the condition number of the matrix produced by `numpy.linalg.cond()` is this measure.\n",
    "2. For the specific problem of linear least squares, a method based on QR factorization is better than forming the Gram matrix, which squares the condition number. Although we did not discuss the theory of QR better, it exemplifies a common pattern or trick in linear algebraic problems: converting the original problem (or matrix) into one whose structure is easier to solve.\n",
    "3. There is a potential tradeoff between _reliability_ of an algorithm and _speed_.\n",
    "\n",
    "Many of these algorithmic and engineering ideas have already been worked out, and are encapsulated in Numpy's built-in [`lstsq()` solver](https://docs.scipy.org/doc/numpy/reference/generated/numpy.linalg.lstsq.html) or its [Scipy equivalent](https://docs.scipy.org/doc/scipy/reference/generated/scipy.linalg.lstsq.html). So when you need to solve a linear system in a future notebook, you can use this function."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true,
    "nbgrader": {
     "grade": false,
     "grade_id": "cell-54da293d2991564f",
     "locked": true,
     "schema_version": 1,
     "solution": false
    }
   },
   "source": [
    "**Fin!** If you've gotten this far without errors, your notebook is ready to submit."
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": [],
  "celltoolbar": "Create Assignment",
  "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.5.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}