labs/linear_equations/notebook.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "# Linear Equations\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Outline \n",
    "\n",
    "1. Setup\n",
    "2. Special cases\n",
    "3. L-U Factorization\n",
    "4. Pivoting\n",
    "5. Benchmarking\n",
    "6. Ill conditioning\n",
    "7. Iterative methods\n",
    "8. Resources\n",
    "\n",
    "\n",
    "## Setup"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The linear equation is the most elementary problem that arises in computational economic analysis. In a linear equation, an $n \\times n$ matrix $A$ and an n-vector $b$ are given, and one must compute the $n$-vector $x$ that satisfies\n",
    "$Ax = b$.\n",
    "\n",
    "Linear equations arise naturally in many  economic applications such as\n",
    "\n",
    "- Linear multicommodity market equilibrium models\n",
    "- Finite-state financial market models \n",
    "- Markov chain models \n",
    "- Ordinary least squares\n",
    "\n",
    "They more commonly arise indirectly from numerical solution to nonlinear and functional equations:\n",
    "\n",
    "- Nonlinear multicommodity market models\n",
    "- Multiperson static game models\n",
    "- Dynamic optimization models\n",
    "- Rational expectations models \n",
    "\n",
    "Applications often require the repeated solution of very large linear equation systems. In these situations, issues regarding speed, storage requirements, and preciseness of the solution of such equations can arise."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import numpy as np\n",
    "\n",
    "from linear_algorithms import backward_substitution\n",
    "from linear_algorithms import forward_substitution\n",
    "from linear_algorithms import gauss_seidel\n",
    "from linear_algorithms import solve\n",
    "\n",
    "from linear_plots import plot_operation_count\n",
    "\n",
    "from temfpy.linear_equations import get_ill_cond_lin_eq\n",
    "from linear_problems import get_inverse_demand_problem"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "## Special cases\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can start with some special cases to develop a basic understanding for the core building blocks for more complicated settings. Let's start with the case of a lower triangular matrix $A$, where we can solve the linear equation by a simple backward or forward substitution. Let's consider the following setup.\n",
    "\n",
    "$$\n",
    "A  =  \\begin{bmatrix}\n",
    "   a_{11} & 0      & 0 \\\\\n",
    "   a_{21} & a_{22} & 0 \\\\\n",
    "   a_{31}   & a_{32} & a_{33} \\\\\n",
    "\\end{bmatrix}\n",
    "$$\n",
    "\n",
    "Consider an algorithmic implementation of forward-substitution as an example.\n",
    "\n",
    "$$\n",
    "x_i = \\left ( b_i - \\sum_{j=1}^{i-1} a_{ij}x_j \\right )/a_{ii}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "??forward_substitution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def test_problem():\n",
    "    A = np.tril(np.random.normal(size=(3, 3)))\n",
    "    x = np.random.normal(size=3)\n",
    "    b = np.matmul(A, x)\n",
    "    return A, b, x\n",
    "\n",
    "\n",
    "for _ in range(10):\n",
    "    A, b, x_true = test_problem()\n",
    "    x_solve = forward_substitution(A, b)\n",
    "    np.testing.assert_almost_equal(x_solve, x_true)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### _Questions_ \n",
    "     \n",
    "* How can we make the test code more generic and sample test problems of different dimensions?\n",
    "* Is there a way to control the randomness in the test function?\n",
    "* Is there software out there that allows to automate parts of the testing?\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If we have an upper triangular matrix, we can use backward substitution to solve the linear system"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "??backward_substitution"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### _Exercise_ \n",
    "\n",
    "\n",
    "* Implement the same testing setup as above the backward-substitution function."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can now build on these two functions to tackle more complex tasks. This is a good example on how to develop scientific software step-by-step ensuring that each component is well tested before integrating into more involved settings."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "## L-U Factorization\n",
    "        \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Most linear equations encountered in practice, however, do not have a triangular $A$ matrix. Doolittle and Crout have shown that any matrix $A$ can be decomposed into the product of a (row-permuted) lower and upper triangular matrix $L$ and $U$, respectively $A=L \\times U$ using **Gaussian elimination**. We will not look into the Gaussian elimination algorithm, but there is an example application in our textbook where you can follow along step by step. The L-U algorithm is designed to decompose the $A$ matrix into the product of lower and upper triangular matrices, allowing the linear equation to be solved using a combination of backward and forward substitution. \n",
    "\n",
    "Here are the two core steps:\n",
    "\n",
    "* Factorization phase\n",
    "\n",
    "\\begin{align*}\n",
    "A = LU\n",
    "\\end{align*}\n",
    "\n",
    "* Solution phase:\n",
    "\n",
    "\\begin{align*}\n",
    "Ax = (LU)x=L(Ux) = b, \n",
    "\\end{align*}\n",
    "\n",
    "where we solve $Ly = b$ using forward-substitution and $Ux=y$ using backward-substitution.\n",
    "\n",
    "Adding to this the two building blocks we developed earlier `forward_substitution` and `backward_substitution`, we can now write a quite generic function to solve systems of linear equations. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "??solve"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's see if this is actually working."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "A = np.array([[3, 1], [1, 2]])\n",
    "x_true = np.array([9, 8])\n",
    "b = A @ x_true\n",
    "x_solve = solve(A, b)\n",
    "np.testing.assert_almost_equal(x_true, x_solve)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Pivoting\n",
    "\n",
    "\n",
    "Rounding error can cause serious error when solving linear equations. Let's consider the following example, where $\\epsilon$ is a tiny number.\n",
    "\n",
    "$$\n",
    "\\begin{bmatrix} \\epsilon & 1\\\\ 1 & 1 \\end{bmatrix} \\times \\left[ \\begin{array}{c} x_1 \\\\ x_2 \\end{array} \\right] = \\left[\\begin{array}{c} 1 \\\\ 2 \\end{array} \\right] \n",
    "$$\n",
    "\n",
    "\n",
    "It is easy to verify that the right solution is\n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    "x_1 & = \\frac{1}{1 - \\epsilon} \\\\\n",
    "x_2 & = \\frac{1 - 2 \\epsilon}{1 - \\epsilon}\n",
    "\\end{aligned}\n",
    "$$\n",
    "\n",
    "and thus $x_1$ is slightly more than one and $x_2$ is slightly less than one. To solve the system using Gaussian elimination we need to add $-1/\\epsilon$ times the first row to the second row. We end up with\n",
    "\n",
    "\n",
    "$$\n",
    "\\begin{bmatrix} \\epsilon & 1 \\\\ 0 & 1 - \\frac{1}{\\epsilon} \\end{bmatrix} \\times \\left[ \\begin{array}{c} x_1 \\\\ x_2 \\end{array} \\right] = \\left[ \\begin{array}{c} 1 \\\\ 2 - \\frac{1}{\\epsilon} \\end{array} \\right],\n",
    "$$\n",
    "\n",
    "which we can then solve recursively.\n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    "x_2 & = \\frac{2 - 1/\\epsilon}{1 - 1/\\epsilon} \\\\\n",
    "x_1 & = \\frac{1 - x_2}{\\epsilon}\n",
    "\\end{aligned}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's translate this into code."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "eps = 1e-17\n",
    "A = np.array([[eps, 1], [1, 1]])\n",
    "b = np.array([1, 2])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can now use our solution algorithm."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0., 1.])"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "solve(A, b)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "But we have to realize that the results are grossly inaccurate.  What happened? Is there any hope to apply `numpy`'s routine?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([1., 1.])"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.linalg.solve(A, b)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This algorithm  does automatically check whether such rounding errors can be avoided by simply changing the order of rows. This is called **pivoting** and changes the recursive solution to\n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    "x_2 & = \\frac{1 - 2\\epsilon}{1 - \\epsilon} \\\\\n",
    "x_1 & = 2 - x_2\n",
    "\\end{aligned}\n",
    "$$\n",
    "\n",
    "which can be solved more accurately. Our implementation also solves the modified problem well. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([1., 1.])"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A = np.array([[1, 1], [eps, 1]])\n",
    "b = np.array([2, 1])\n",
    "solve(A, b)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Building your own numerical routines is usually the only way to really understand the algorithms and learn about all the potential pitfalls. However, the default should be to rely on battle-tested production code. For linear algebra there are numerous well established libraries available. Building your own numerical routines is usually the only way to really understand the algorithms and learn about all the potential pitfalls. However, the default should be to rely on battle-tested production code. For linear algebra there are numerous well established libraries available."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "## Benchmarking\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "How does solving a system of linear equations by an $L-U$ decomposition compare to other alternatives of solving the system of linear equations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX0AAAEJCAYAAAB4yveGAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAAqL0lEQVR4nO3dd5xU5d338c9vZvsuuJQFqYKKBRVBV4JKFMWCiqKJCkYNSXhi7vvBxFtNbIkllieaGBI1GkNsGCuW3GCNii12UWygRlSElV6WbWybuZ4/zllmFpa2Mztndub7fr3mNee6zplzfjvKd89ep5lzDhERyQ6hoAsQEZHUUeiLiGQRhb6ISBZR6IuIZBGFvohIFskJuoBtGTdunHv22WeDLkNEpLOxtjrTfk9/9erVQZcgIpIx0j70RUQkeRT6IiJZRKEvIpJFFPoiIllEoS8ikkUU+iIiWUShLyKSbj56BN69E6LRpK867S/OEhHJKtXL4ekLoX49fPIYfP9O6NonaavXnn47lZSUtGovWrSIfffdt1XfVVddxY033pjKskSks3v6V17gA1R9CwU7JXX1Cn0RkXSxYBZ8OjvWPvFmyCtK6iYU+iIi6WDDOnjql7H2AT+EXQ9P+mY6/Zj+oEue6rB1L7r+hA5bt4hIK//6DdSu9KZLdoajr+mQzWhPP0nM2ryh3Rb7RUQ2+vJF+OC+WHv8NCgs7ZBNKfSTpEePHqxbt65V39q1a+nZs2dAFYlIp9BQA0+cF2vvcwrs1XGjDJ1+eCddhmBKSkro06cPc+bMYezYsaxdu5Znn32W8847b9sfFpHs9eK1ULnYmy7sBsf9vkM31+lDPyh1dXX0799/Y/uCCy7g3nvvZerUqVx44YUAXHnlley2225BlSgi6W7JO/D27bH2uOuhpFeHblKh307RLVwp99JLL6W4EhHplJobYNa5gPPaux8FwyZ2+GY1pi8iEoRXb4TVn3vTeSUw/k+QghM/FPoiIqm2/BN4bVqsfdRVUDowJZtW6IuIpFKkGWafC9Fmrz1gFJRPSdnmFfoiIqn01m2wdJ43Hc6Dk26BUOqiWKEvIpIqa76El66LtQ+/GMr2SGkJCn0RkVRwzrsIq7nea/feDw5N/XU82wx9M7vLzFaa2Sdxfd3N7Hkz+8J/7xY371IzW2hmn5vZsXH9B5rZx/68m62T359g01sri4hs1fszYNG/vWkLw4S/QDg35WVsz57+PcC4TfouAeY454YAc/w2ZjYUmATs43/mNjML+5/5K3AOMMR/bbpOiROJRIIuQUSSpWopPHd5rH3Iz6Hv8EBK2WboO+deBdZu0j0BmOFPzwBOjut/yDnX4Jz7GlgIjDSzPkBX59ybzjkH3Bv3mU7t5ZdfZsyYMZx66qnstddenHnmmTjneOaZZzj99NNbLXfiiScC8Nxzz3HwwQdzwAEHcNppp1FTUwPAoEGDuPrqqxk9ejSPPPIIN998M0OHDmXYsGFMmjQJgNraWn7yk59w0EEHMWLECGbNmpX6H1pEtp9z8OQF0FDltbvvBmMuCayc9l6R29s5twzAObfMzFquG+4HvBW3XIXf1+RPb9rfJjM7B++vAgYO3Ma5q1cl96kyrde9frsWmzdvHvPnz6dv374ceuihvP766xx99NH87Gc/o7a2luLiYh5++GEmTpzI6tWrufbaa3nhhRcoLi7mhhtuYNq0aVxxxRUAFBQU8NprrwHQt29fvv76a/Lz86msrATguuuu48gjj+Suu+6isrKSkSNHctRRR1FcXNwhX4GIJGj+4/CfZ2Ltk26B3MLAykn2gdy2xundVvrb5Jyb7pwrd86Vl5WVJa24jjJy5Ej69+9PKBRi+PDhLFq0iJycHMaNG8cTTzxBc3MzTz31FBMmTOCtt95iwYIFHHrooQwfPpwZM2bwzTffbFzXxImxy7CHDRvGmWeeyX333UdOjvf7+bnnnuP6669n+PDhjBkzhvr6ehYvXpzyn1lEtkPtGnj6oli7fAoMOjS4emj/nv4KM+vj7+X3Afw7/1MBDIhbrj+w1O/v30Z/RsjPz984HQ6HaW72LrqYOHEit956K927d+eggw6iS5cuOOc4+uijefDBB9tcV/we+1NPPcWrr77K7Nmzueaaa5g/fz7OOR577DH23HPPjv2hRCRx/7oU6lZ70137eVfeBqy9oT8bmAxc77/Piut/wMymAX3xDti+45yLmFm1mY0C3gZ+CNySUOUttnMIJghjxoxhypQp/P3vf9+4Bz9q1CimTp3KwoUL2X333amrq6OiooI99mh9rm40GmXJkiUcccQRjB49mgceeICamhqOPfZYbrnlFm655RbMjHnz5jFixIggfjwR2Zr/PAcfPRxrj/8TFHQNrh7f9pyy+SDwJrCnmVWY2RS8sD/azL4AjvbbOOfmAzOBBcCzwFTnXMtpKP8N3IF3cPdL4BkyXDgcZvz48TzzzDOMHz8egLKyMu655x7OOOMMhg0bxqhRo/jss882+2wkEuGss85iv/32Y8SIEZx//vmUlpZy+eWX09TUxLBhw9h33325/PLLN/usiASsvgqePD/W3u902OPYLS+fQuadTJO+ysvL3dy5c4MuQ0Rk+z11Ibx7hzdd1AOmvgvFPVJdRZvXQumKXBGRZPrmjVjgg/ckrNQH/hYp9EVEkqVpg/9gFN8e42Df7wdXTxsU+iIiyfLKDbD2S286rwucMC0lD0bZEQp9EZFkWPoBvH5zrH3M1bDTFq9BDYxCX0QkUZEm78EoLScr7jIaDvhRoCVtiUJfRCRRb9wMyz/2pnMK4KSbU/pglB2RnlWJiHQWq7+Al2+ItY+4DHrsFlw926DQFxFpr2gUZv8cIg1eu89wGDU10JK2RaEvItJec++ExW9606Ec/8Eo7b27TWoo9EVE2qNyCbxwVaw9+nzYeb/AytleCn0RkR3lnHdvnUbvAUj03AMO+1WwNW0nhb6IyI76aCYsfN5vGJz0F8jJ3+pH0oVCX0RkR9SsgmcvjrVHngMDvxNcPTtIoS8isiOeuQg2rPOmdxoAY68Itp4dpNAXEdlenz3lPfO2xYl/hvySwMppD4W+iMj22FDp3Se/xf4/gN2PCqyc9lLoi4hsj+evgOpl3nRxGRx7XbD1tJNCX0RkW756Bd6fEWsffyMUdQ+ungQo9EVEtqaxDp74Ray913gYOiG4ehKk0BcR2ZqXroN1i7zp/J3ghD+m3YNRdoRCX0RkS759D966LdY+9jrosnNw9SSBQl9EpC3NjTDr5+CiXnvw4TDirGBrSgKFvohIW17/M6yc703nFsGJN3XqYZ0WCn0RkU2t/Axe+X2sfeRvoPvg4OpJIoW+iEi8aMR73m20yWv3OxC+81/B1pRECn0RkXjvTIeKd73pUK53B81QONiakkihLyLSYt0imHN1rH3YL6H30MDK6QgJhb6ZnW9m883sEzN70MwKzKy7mT1vZl/4793ilr/UzBaa2edmdmzi5YuIJEk04j3vtqnOa/caCqMvCLamDtDu0DezfsAvgHLn3L5AGJgEXALMcc4NAeb4bcxsqD9/H2AccJuZZc7fTCLSub14DXz9qjdtIf/BKHnB1tQBEh3eyQEKzSwHKAKWAhOAlptUzABO9qcnAA855xqcc18DC4GRCW5fRCRxnzwOr/0p1h59AfQ/MLh6OlC7Q9859y1wI7AYWAasd849B/R2zi3zl1kG9PI/0g9YEreKCr9vM2Z2jpnNNbO5q1atam+JIiLbtvxjmDU11h5yDBxxWXD1dLBEhne64e29Dwb6AsVmtrXL1dq6qsG1taBzbrpzrtw5V15WVtbeEkVEtq5uLTx0Zmwcv/tu8L2/Z9TZOptKZHjnKOBr59wq51wT8DhwCLDCzPoA+O8r/eUrgAFxn++PNxwkIpJ6kWZ49MdQ+Y3XziuBSQ9AYWmgZXW0REJ/MTDKzIrMzICxwKfAbGCyv8xkYJY/PRuYZGb5ZjYYGAK8k8D2RUTab85V8NXLsfYpf4NeewVVTcrktPeDzrm3zexR4H2gGZgHTAdKgJlmNgXvF8Np/vLzzWwmsMBffqpzLpJg/SIiO+6jR+CNW2Ltwy+GvccHV08KmXNtDqunjfLycjd37tygyxCRTLHsQ7jzGGiu99p7HOcN64Qy7lrVNu8Ol3E/pYjIFtWu9g7ctgR+zz3ge9MzMfC3KHt+UhHJbpEmeORHsN4/czy/q7eHX9A10LJSTaEvItnhucth0b/9hnmnZvYcEmhJQVDoi0jm++BBePuvsfYRl8Ge44KrJ0AKfRHJbN++D0+cF2vvNR6++8vg6gmYQl9EMlfNSnj4LIg0eO2yveCU27PqwO2msvcnF5HMFmmCmZOh6luvnb+Td+A2v0uwdQVMoS8imenZS2HxG37D4NQ7ocdugZaUDhT6IpJ53v8HvPv3WHvsFTDk6ODqSSMKfRHJLBVz4am4J14NPRlGnx9YOelGoS8imaN6hX/gttFr99oHJtwK1uYdCbKSQl9EMkNzI8w8G6qXee2CUph0P+SXBFpWulHoi0hmeOYiWPK2N20hOPUu6D442JrSkEJfRDq/uXfDe3fH2kf9FnYfG1w9aUyhLyKd2+K34elfxdr7ngqH/Dy4etKcQl9EOq+qZd44frTJa++8H5x0iw7cboVCX0Q6p+YG70ydmhVeu7A7TLwf8oqCrSvNKfRFpPNxzjsX/1v/qXoWhtPugW67BFpWZ6DQF5HO5907YN59sfYx18KuhwdXTyei0BeRzuWbN+DZS2LtYRNh1H8HV08no9AXkc5jfQXM/CFEm712n/3hxJt04HYHKPRFpHNo2uAduK1d5bWLenoHbnMLg62rk1Hoi0j6cw6evACWzvPaoRw4fQaUDgi2rk5IoS8i6e/tv8GHD8Tax/4OBo0Orp5OTKEvIunt61fhX5fF2sPPgpE/Da6eTk6hLyLpq3IxPPIjcBGv3e9AOOGPOnCbAIW+iKSnxjp46EyoW+O1i3vB6f+A3IJg6+rkEgp9Mys1s0fN7DMz+9TMDjaz7mb2vJl94b93i1v+UjNbaGafm9mxiZcvIhnJOXjiPFj+kdcO5cDp98JO/YKtKwMkuqd/E/Csc24vYH/gU+ASYI5zbggwx29jZkOBScA+wDjgNjMLJ7h9EclEb94KH8+MtY/7PexycHD1ZJB2h76ZdQUOA+4EcM41OucqgQnADH+xGcDJ/vQE4CHnXINz7mtgITCyvdsXkQz15Uvw/OWx9gGTofwnwdWTYRLZ098VWAXcbWbzzOwOMysGejvnlgH477385fsBS+I+X+H3bcbMzjGzuWY2d9WqVQmUKCKdyrpF8OiPwUW9dv+RcPwfdOA2iRIJ/RzgAOCvzrkRQC3+UM4WtPVfzbW1oHNuunOu3DlXXlZWlkCJItJpNNZ6B243rPPaJTt74/g5+cHWlWESCf0KoMI55z+UkkfxfgmsMLM+AP77yrjl4y+f6w8sTWD7IpIpnINZU2HFJ147lAsT/wFd+wRbVwZqd+g755YDS8xsT79rLLAAmA1M9vsmA7P86dnAJDPLN7PBwBDgnfZuX0QyyKt/gPn/jLVP+CMM0CG/jpCT4Od/DtxvZnnAV8CP8X6RzDSzKcBi4DQA59x8M5uJ94uhGZjqXMsVFyKStV75A7x0XaxdPgUOnLzl5SUh5lybw+ppo7y83M2dOzfoMkQk2ZyDl38Hr9wQ6xt8GJz5GOTkBVdX5mjz6Heie/oiIjvOOZhzNbw2Lda36xiY9KACv4Mp9EUktZyD534Db/4l1rf7UTDxPt0bPwUU+iKSOs55jzp8+/ZY3x7jdGpmCin0RSQ1olF4+pcw985Y317j4dS7NaSTQgp9Eel40Sg8eR68f2+sb+jJ8P07IJwbWFnZSKEvIh0rGoFZ57Z+8tV+p8HJt0NYEZRq+sZFpONEmuF//ws+fiTWt/8ZMOFWCOkmu0FQ6ItIx4g0weM/bX2l7Yiz4cSbFPgBUuiLSPI1N8JjP4FPn4j1lf8Ejv8jhPTAviAp9EUkuZobYOZk+M8zsb6RP4PjbtAtktOAQl9EkqepHh4+CxY+H+s7+Fw45loFfppQ6ItIcjTWwUM/gK9eivWNPh/GXqnATyMKfRFJXGMtPDARFv071nfYRXDEZQr8NKPQF5HENFTD/afD4jdifUf8Gg6/KLiaZIsU+iLSfvXr4b5ToSLueUhjr4TvXhBcTbJVCn0RaZ8NlXDf9+Db92J9x1wHh5wbWEmybQp9EdlxdWvhHyfDsg9jfcf9Hr7zs8BKku2j0BeRHVO7Bu6dACs+jvWdMA0OmhJcTbLdFPoisv1qVsG9J8HKBX6HwUk3wwE/DLQs2X4KfRHZPtXLYcZJsPpzv8Pg5Ntg+A8CLUt2jEJfRLatainMOBHWLPTaFoJTpsOw04KtS3aYQl9Etq5yiRf467722haGU++EfU4Jti5pF4W+iGzZum9gxnioXOy1Q7lw2t2w94nB1iXtptAXkbat/QruORGqKrx2OM97gPmexwVblyREoS8im1u90NvDr17mtcP5MOl+GHJ0sHVJwhT6ItLaqs+9MfyaFV47pwDOeBB2OzLYuiQpFPoiErNigRf4dau9dm4R/OBhGHxYsHVJ0iT83DIzC5vZPDN70m93N7PnzewL/71b3LKXmtlCM/vczI5NdNsikkTLPoJ7TogFfl4JnPWYAj/DJONhlecBn8a1LwHmOOeGAHP8NmY2FJgE7AOMA24zMz0dWSQdLJ3n7eFvWOu187vC2f+EXQ4Jti5JuoRC38z6AycAd8R1TwBm+NMzgJPj+h9yzjU4574GFgIjE9m+iCRBxVyYMQHqK712wU5w9v/CAP3zzESJ7un/GbgIiMb19XbOLQPw33v5/f2AJXHLVfh9mzGzc8xsrpnNXbVqVYIlisgWffki3HsyNKz32oXd4Iezof+BgZYlHafdoW9m44GVzrn3trmw/5E2+lxbCzrnpjvnyp1z5WVlZe0tUUS2JNIMc66Bf3wPGqu9vqIeMPkJ6Ds80NKkYyVy9s6hwElmdjxQAHQ1s/uAFWbWxzm3zMz6ACv95SuAAXGf7w8sTWD7ItIeVUvhsf8D37we6yvp7Q3p9B4aWFmSGu3e03fOXeqc6++cG4R3gPZF59xZwGxgsr/YZGCWPz0bmGRm+WY2GBgCvIOIpM4Xz8Pto1sH/uDD4Wf/VuBniY44T/96YKaZTQEWA6cBOOfmm9lMYAHQDEx1zkU6YPsisqlIE7x4Dbx+U6zPQjDmMu95tiGdSJctzLk2h9XTRnl5uZs7d27QZYh0XpVL4LEpsOTtWF+XPvD9O2DQ6ODqko7W1nFUXZErktE+exr+979jp2MC7H4UnPI3KO4ZWFkSHIW+SCZqboQXroK3bo31WRjGXg6HnAehZFyXKZ2RQl8k06xbBI/8GJa+H+vr2g9OvQsGjgqsLEkPCn2RTLJgNsw6N3axFcAe4+Dkv0JR9+DqkrSh0BfJBE318Pzl8M70WF8oB476LRw8FazNY3qShRT6Ip3dmi/hkR/B8o9ifaUD4dR7dDsF2YxCX6Qz++QxmH1e7FYKAHuNhwm3QmFpYGVJ+lLoi3RGTRvg2UvgvXtifeE8OOY6GPlTDefIFin0RTqb1V94wzkrPon1dRsMp90NfUcEVpZ0Dgp9kc7kw4fhyfOhqTbWt88pcOLNUNA1uLqk01Doi3QGjbXw9EXwwX2xvnA+HHc9HPhjDefIdlPoi6S7lZ96wzmrPov19dgdTrsHdt4vqKqkk1Loi6Qr5+CD++GpX0Lzhlj/sIlwwjTILwmuNum0FPoi6aihBp66AD56ONaXUwjH/wFGnKXhHGk3hb5Iuln+CTwyGdYsjPWV7eUN5/TaO7CyJDMo9EXShXPw3t3wzCUQaYj1Dz8Ljv895BUHV5tkDIW+SDqor4InzoP5j8f6coth/DTYf1JwdUnGUeiLBG3pB97ZOeu+jvX12scbzinbI6CiJFMp9EWC4hy883d47tcQaYz1H/hjGPc7yC0MrjbJWAp9kSCs/xaeuQg+ezLWl9cFTvwz7HdqYGVJ5lPoi6TShkp47U/w9u3QXB/r33mYN5zTY7egKpMsodAXSYXmBm8o5983woZ1reeNPAeOvgZyC4KpTbKKQl+kI0Wj8PFMePE6WL+49by+I+Doq2HwYcHUJllJoS/SURbOgeevhBUft+7vNgjGXgFDT4FQKJDSJHsp9EWSbekH8MKV8NXLrfuLesLhF8OBP4KcvAAKE1HoiyTPukUw5xr45NHW/blFcPC5cMjPdc97CZxCXyRRtWvg1T/Au3dAtCnWb2E44Icw5hLosnNw9YnEaXfom9kA4F5gZyAKTHfO3WRm3YGHgUHAIuB059w6/zOXAlOACPAL59y/EqpeJEiNdfDWbfD6TdBQ1XreXuNh7JW6olbSTiJ7+s3Ahc65982sC/CemT0P/AiY45y73swuAS4BLjazocAkYB+gL/CCme3hnIsk9iOIpFik2bvP/cu/g+plrecNGOWdkTPwO8HUJrIN7Q5959wyYJk/XW1mnwL9gAnAGH+xGcDLwMV+/0POuQbgazNbCIwE3mxvDSIp5Rx8/gy8cBWs/rz1vJ57wFFXwZ7H6173ktaSMqZvZoOAEcDbQG//FwLOuWVm1stfrB/wVtzHKvy+ttZ3DnAOwMCBA5NRokhilrwDz18BizfZRynZGY641Lv9cViHyCT9Jfx/qZmVAI8B/+Ocq7It7+W0NcO1taBzbjowHaC8vLzNZURSYvUXMOe38OkTrfvzusDo82DU/9V97qVTSSj0zSwXL/Dvd8613Ah8hZn18ffy+wAr/f4KYEDcx/sDSxPZvkiHqV4OL18P798L8YedQrlw0BQ47FdQ3DO4+kTaKZGzdwy4E/jUOTctbtZsYDJwvf8+K67/ATObhncgdwjwTnu3L9IhGqrh9Zvhzb9AU13reft+H468HLoPDqY2kSRIZE//UOBs4GMz+8Dvuwwv7Gea2RRgMXAagHNuvpnNBBbgnfkzVWfuSNpoboT37oFXboC61a3nDT4Mjvot9DsgkNJEksmcS+8h8/Lycjd37tygy5BM5RzM/yfMubr1k6sAeu/rhf3uY3VGjnRGbf5Pq9MNJHt9/W/vjJyl77fu32kAHPkb2O903RBNMo5CX7LPsg/hxWvhi+da9xeUwmG/hIN+qnvbS8ZS6Et2qF0NHz8CHzwAyz9qPS+cD6P+C0afD4XdgqlPJEUU+pK5mhvhi3/BBw9679HmTRYwGP4DOOIy2Kl/ICWKpJpCXzKLc7DsAy/oP34ENqzdfJlwPuw9Hr57IfTeJ+UligRJoS+ZoXo5fDQTPnwQVi5oe5kB34H9z4B9ToHC0pSWJ5IuFPrSeTXVw+dPe+P0X84BF918ma79Yf9JXtj33D31NYqkGYW+dC7OQcVc79bG8x+H+vWbL5NbBHuf5I3XD/quTrsUiaPQl85hfQV8+JA3fLNmYdvL7DLaC/qhJ0F+l9TWJ9JJKPQlfTXWeXe3/PAB+OoV2rwpa7dBsP8PYP+J3rRIhnDO0dAcpSA3nNT1KvQlvTgH37zhBf38WdBYvfkyeV1gn5O9vfqBB+sWCdLp1DU2s6KqgeXr61lZXc+KqnqWr29gRXU9K9bXe+9VDRyxZxl/O7s8qdtW6Et6WLcoNnyzblEbCxjsOsYL+r3GQ15RausT2Q6NzVFW1TSwosoP76p6VlQ3tAryFevrqW7Y9JqRtq2oakh6jQp9CU5DNSyY5Z1T/81rbS/TYwgMPwOGTdQFVBKYaNSxpraRFVXenvny9Q1x016Yr6yuZ3VNY1K3W13flNT1gUJfUi0ahUWvekH/6ezN71kPULCTd+/6/X8A/cs1fCMdpjkSZU1tI6uqG1hV08Dq6gZWVjewssoL8uVV9aysqmdldQPN0eTdkTgvHKJX13x6dy1g564F9Oqaz85dC+gdN92rawEl+cmPaIW+dLz133rPll3yNnz2NFRVbL6MhWD3o7zz6fc8Xjc8k3aLRh1r67wgX13T0Ordm47NW1vXSDLvLh8y6Fnihbn3yt/k3Xt1K8plK4+W7VAKfUmuaBRWfeaF/OK3vNf6xVtevtdQL+iHnQ5ddk5dndKpOOdYv6Fp4x55WwHeMm9tbSORJO6VtygtyqV3l9Z75ZuGec+SPHLC6X1diEJfEtNUD0vnxUJ+yVttXzAVr7A77Head1C2z/4avslSzZEo6+qaWFvbyJpaL6xXbxxmafTe4/bUmyLJD3Iz6F6UR1mXfMq65NOzJJ+eJXmtgrxl+CXZp04GRaEvO6ZurTdMs/hNWPy29wCSyDYOXuUWeWPzAw+GgaO8i6hy8lJTr6RMQ3PEC/CaRtbWNvph3shaP9DX1DSyrq6lr5H1G5qSOrQSr7Qol54l+ZSV5NOzi/fuhXos4MtK8ulenP575smm0Jctcw4qv/GHafw9+VWfbftzxWVeuLeE/M7DIJzb8fVK0jjnqG2MsLYmthce/1qzaajXNFLb2LGPvO6Sn9NmgPf02y3TPUryyM/JjL3yjqDQl5hIM6yc3zrkq5dt+3M9do8L+YOh+64askkj0aijur6Zyg2NVNY1Ubmhico6b097XW0TlRsaWRcX5C1h3tjcxg3sksgMSgtz6V6cR4/ifLoV524W4GVxIZ8pwytBU+hns8Za7+ZlLSFf8S401mz9M6Ec6DPcD/lRMGAUlJSlpNxsF4k6quubqKxrYl1dI5Ubmlhf5wW4F+Stp9fHhXsHHNfcTE7I6FacR4/iPLr7L286n+4leXQv8vtKvPfSwtysG1pJBwr9bFKzMnZGzeI3vWfFum38SZ7XBQaMjA3V9DtQV8MmoOV+KlX1TVTXN7PeD+51dbG98PVbCPGq+o4bA29LQW6IHsX5rQK8W6swbwlwb5muBTmBnYYo20+hn2k2rIPKxbDuG288vmV6zRew9qttf75Ln9gwzcBR3pOlQvqzukVzJEpNQ/PGwK6ub94Y4FV+u7q+KdZX39IXm98Y6dhhk7Z0yc9hp6JcuhXlUVqUy06FuZQW5VJa6LW7FeVttjdelKd4yET6r9rZNNb6gb7YC/WN4f4NrFsMDds4XXJTvYZ6T5RqCfnSgRk7Ht8ciVLbGKG2odkP7iaqNjRvFtBVGzYNbG+56vqmDj9YuS1dC3IobRXceXQryqW0MJedirwhk9KilpfX7lqYS66GUcSn0E83zQ1QuQQqF7W9x163uv3rDud5wzMtB137HwRF3ZNWerI1RaIbA7quMUJNQzO1G18Rahv9eQ1x8xr9ef7nahtj8xs6+MDk9soLh+hSkEPXwtxWIR4f3N2Kvb3wnfz+0iJv+ERj4JIohX6qRZqh6tu4vfTFraerl9HmfeO3V06ht7deOhC67QKlu8Smy/ZO+u0NolFvjHpDU8R7NUaoj5ve0OS3/ek6f0/bC+UIdX5w124W7JFAhkG2xQxK8nPoWpDrBXdBLl0Lc+hS4AV4l7h2y/yWgG9p6ywUCZJCPxGRZu9+7w0trxrvPb6vvsq710zLHvv6b7d98HRrQrne3SZbBfqg2HRJLxzQ0BylMRKlsTnqTTdHqV/dyIamDdT7Abx5SEc3C+k2263mpV8wb4kZlOTlUJQfpji/JbBbwrjtgO6ySV9JXg6hUGYOf0l2SHnom9k44CYgDNzhnLs+pQVEI95piS0B3Sqk2+prO8xdQw3WvCH55RGiOq8Xa/P6sDanN6ty+rAi3Jvl1ptloV6scKXUR4zGyiiNq71gb2iK0hippLF5LY1+2GeCkEFxfg4l+TkUt7zywhv7ivLCW53X+rNhCnPDOrtEsl5KQ9/MwsCtwNFABfCumc12zi1I5nY++ueN5Cx6hbxILfnNteRF6yiI1JIfraXA1SdlG4lEx0pXyhJXxhJXRoUrY4nrtbG93PWgqX5r/1l28EBtChTkhijM9UK1IC+8cbowL0xB7ubtkvwwRXmxQC7K98Pb72tp5+eEFNIiSZbqPf2RwELn3FcAZvYQMAFIaujXL/mAketfTeYq2xRxRi2FVFNIjSuklgJqXCE1fruGQmooYKXr5oe7F/INdPx9Z/LCIfJzQuTFvQpz40LYD+eC3DCFedsR2q2Wjy2TnxPScIdIJ5Lq0O8HLIlrVwDf2XQhMzsHOAdg4MCBO7yR5nDxVufXuIK4YC7wA9sL6eqN4V1ENYXUtiwbH+R+X1OogJxwiJxQiJwcIydk5IRChENGbtj89xC5fgD3zgkxIKcljMPkhb0wzs+JC+hwiPzckD8vvFlwb1w2HN7Y3nR+Xlh7yCLStlSHfltJtNmpKs656cB0gPLy8h0+laXkO2fx0pKDiOZ1IZpbjMvrgssvgbyuuLxicnPCGwM5HDIKw0ZJKES/kJET9oI7x59uWSY3FCIcbgl2L9AVrCLS2aQ69CuAAXHt/sDSZG9kWPl3ofy7yV6tiEinl+orPd4FhpjZYDPLAyYBs1Ncg4hI1krpnr5zrtnMzgX+hXfK5l3OufmprEFEJJul/Dx959zTwNOp3q6IiKR+eEdERAKk0BcRySIKfRGRLKLQFxHJIuZS+fy1djCzVcA3QdeRoJ5AAjfCzyj6LlrT99Gavo+YRL+L1c65cZt2pn3oZwIzm+ucKw+6jnSg76I1fR+t6fuI6ajvQsM7IiJZRKEvIpJFFPqpMT3oAtKIvovW9H20pu8jpkO+C43pi4hkEe3pi4hkEYW+iEgWUeh3EDMbYGYvmdmnZjbfzM4LuqZ0YGZhM5tnZk8GXUvQzKzUzB41s8/8/08ODrqmoJjZ+f6/k0/M7EEzKwi6plQys7vMbKWZfRLX193MnjezL/z3bsnYlkK/4zQDFzrn9gZGAVPNbGjANaWD84BPgy4iTdwEPOuc2wvYnyz9XsysH/ALoNw5ty/ebdcnBVtVyt0DbHoh1SXAHOfcEGCO306YQr+DOOeWOefe96er8f5B9wu2qmCZWX/gBOCOoGsJmpl1BQ4D7gRwzjU65yoDLSpYOUChmeUARXTAE/XSmXPuVWDtJt0TgBn+9Azg5GRsS6GfAmY2CBgBvB1wKUH7M3AREA24jnSwK7AKuNsf7rrDzIqDLioIzrlvgRuBxcAyYL1z7rlgq0oLvZ1zy8DbiQR6JWOlCv0OZmYlwGPA/zjnqoKuJyhmNh5Y6Zx7L+ha0kQOcADwV+fcCKCWJP353tn4Y9UTgMFAX6DYzM4KtqrMpdDvQGaWixf49zvnHg+6noAdCpxkZouAh4Ajzey+YEsKVAVQ4Zxr+evvUbxfAtnoKOBr59wq51wT8DhwSMA1pYMVZtYHwH9fmYyVKvQ7iJkZ3njtp865aUHXEzTn3KXOuf7OuUF4B+ledM5l7d6cc245sMTM9vS7xgILAiwpSIuBUWZW5P+7GUuWHtTexGxgsj89GZiVjJWm/Bm5WeRQ4GzgYzP7wO+7zH9GsAjAz4H7zSwP+Ar4ccD1BMI597aZPQq8j3fW2zyy7HYMZvYgMAboaWYVwJXA9cBMM5uC94vxtKRsS7dhEBHJHhreERHJIgp9EZEsotAXEckiCn0RkSyi0BcRySIKfcloZna7mV0edB0tzGygmdWYWTjoWiQ76ZRN6dT8K3x7453fHcG7wOleYLpzTvf4EdmE9vQlE5zonOsC7IJ3QcvF+HevFJHWFPqSMZxz651zs4GJwGQz29fM7jGzawHMbIyZVZjZRf4DK5aZ2clmdryZ/cfM1prZZS3rM7OQmV1iZl+a2Rozm2lm3f15g8zMmdlkM1tsZqvN7Ndxnx1pZnPNrMrMVpjZtE0+l+O3+5rZbH/bC83sp3HruMrf5r1mVu0/ZKQ8Nd+mZCqFvmQc59w7eDc0+24bs3cGCvCebXAF8HfgLOBAf/krzGxXf9lf4N3D/HC8uz+uA27dZH2jgT3x7hdzhZnt7fffBNzknOsK7AbM3EK5D/q19gVOBf6fmY2Nm38S3g3qSvHuxfKXrf7wItug0JdMtRTo3kZ/E3CdfzfHh4CeeOFc7ZybD8wHhvnL/gz4tXOuwjnXAFwFnNqyl+77rXNug3PuQ+BDvCdgtWxndzPr6Zyrcc69tWkhZjYA75fGxc65eufcB3gPmDk7brHXnHNPO+ciwD/i1i/SLgp9yVT92PxJRABr/AAF2OC/r4ibvwEo8ad3Af5pZpVmVol358cI3oHjFsvjpuviPjsF2AP4zMze9Z8nsKm+wFr/yWotvqH1E9Y2XX/BJr90RHaIQl8yjpkdhBecryW4qiXAcc650rhXgf+kp61yzn3hnDsD72lHNwCPtvFkrKVAdzPrEtc3ENjm+kXaS6EvGcPMuvp71A8B9znnPk5wlbcD15nZLv76y8xswnbWcpaZlfmnjVb63ZH4ZZxzS4A3gN+ZWYGZDcP7C+H+BOsW2SL9mSiZ4Akza8Z79u4CYBpeYCfqJsCA58ysL96Tix5m+x5mMQ6YZmZFeEM2k5xz9d4zQlo5w691Kd6B4iudc88noXaRNuniLBGRLKLhHRGRLKLQFxHJIgp9EZEsotAXEckiCn0RkSyi0BcRySIKfRGRLKLQFxHJIv8fu+MWctvV0jsAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plot_operation_count()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The right setup for your numerical needs depends on your particular problem. For example, this trade-off looks very different if you have to solve numerous linear equations that only differ in $b$ but not $A$. In this case you only need to compute the inverse once."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### _Exercise_ \n",
    "   \n",
    "\n",
    "* Set up a benchmarking exercise that compares the time to solution for the two approaches for $m=\\{1, 100\\}$ and $n = \\{50, 100\\}$, where $n$ denotes the number of linear equations and $m$ the number of repeated solutions."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "## Ill Conditioning\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Some linear equations are inherently difficult to solve accurately on a computer. This difficulty occurs when the A matrix is structured in such a way that a small perturbation $\\delta b$ in the data vector $b$ induces a large change $\\delta x$ in the solution vector $x$. In such cases the linear equation or, more generally, the $A$ matrix is said to be **ill conditioned**.\n",
    "\n",
    "One measure of ill conditioning in a linear equation Ax = b is the “elasticity” of the solution vector $x$ with respect to the data vector $b$ \n",
    "\n",
    "$$\n",
    "\\epsilon = \\sup_{||\\delta  b|| > 0} \\frac{||\\delta x|| / ||x||}{||\\delta b|| / ||b||}\n",
    "$$\n",
    "\n",
    "The elasticity gives the maximum percentage change in the size of the solution vector $x$ induced by a $1$ percent change in the size of the data vector $b$. If the elasticity is large, then small errors in the computer representation of the data vector $b$ can produce large errors in the computed solution vector x. Equivalently, the computed solution $x$ will have far fewer significant digits than the data vector $b$.\n",
    "\n",
    "In practice, the elasticity is estimated using the condition number of the matrix $A$, which for invertible $A$ is defined by $\\kappa \\equiv ||A|| \\cdot ||A^{-1} ||$. The condition number is always greater than or equal to one. Numerical analysts often use the rough rule of thumb that for each power of $10$ in the condition number, one significant digit is lost in the computed solution vector $x$. Thus, if $A$ has a condition number of $1,000$, the computed solution vector $x$ will have about three fewer significant digits than the data vector $b$.\n",
    "\n",
    "Let's look at an example, where the solution vector is all ones but but the linear equation is notoriously ill-conditioned."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "??get_ill_cond_lin_eq"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "How does the solution error depend on the condition number in this setting."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "rslt = dict(((\"Condition\", []), (\"Error\", []), (\"Dimension\", [])))\n",
    "grid = [5, 10, 15, 25]\n",
    "\n",
    "for n in grid:\n",
    "    A, b, x_true = get_ill_cond_lin_eq(n)\n",
    "    x_solve = np.linalg.solve(A, b)\n",
    "    rslt[\"Condition\"].append(np.linalg.cond(A))\n",
    "    rslt[\"Error\"].append(np.linalg.norm(x_solve - x_true, 1))\n",
    "    rslt[\"Dimension\"].append(n)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>Condition</th>\n",
       "      <th>Error</th>\n",
       "      <th>Dimension</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>2.616969e+04</td>\n",
       "      <td>1.278000e+03</td>\n",
       "      <td>5</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>2.106258e+12</td>\n",
       "      <td>1.762345e+09</td>\n",
       "      <td>10</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>2.582411e+21</td>\n",
       "      <td>4.947153e+16</td>\n",
       "      <td>15</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>2.035776e+22</td>\n",
       "      <td>2.383979e+20</td>\n",
       "      <td>25</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "      Condition         Error  Dimension\n",
       "0  2.616969e+04  1.278000e+03          5\n",
       "1  2.106258e+12  1.762345e+09         10\n",
       "2  2.582411e+21  4.947153e+16         15\n",
       "3  2.035776e+22  2.383979e+20         25"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pd.DataFrame.from_dict(rslt)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There is a more general lesson here. Always be skeptical about the quality of your numerical results. See these two papers for an exploratory analysis of econometric software packages and nonliner optimization. Yes, they are rather old and much progress has been made, but the general points remain valid.\n",
    "\n",
    "* McCullough, B. D., and Vinod, H. D. (2003). [Verifying the Solution from a Nonlinear Solver: A Case Study](https://www.aeaweb.org/articles?id=10.1257/000282803322157133). *American Economic Review*, 93 (3): 873-892.\n",
    "\n",
    "* McCullough, B., D., and Vinod H. D. (1999). [The Numerical Reliability of Econometric Software](https://www.aeaweb.org/articles?id=10.1257/jel.37.2.633). *Journal of Economic Literature*, 37 (2): 633-665. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### _Exercise_\n",
    "\n",
    "Let's consider the following example as well, which is taken from Johansson (2015, p.131).\n",
    "\n",
    "$$\n",
    "\\left[\\begin{array}{c} 1 \\\\ 2 \\end{array}\\right] = \\begin{bmatrix} 1 & \\sqrt{p}\\\\ 1 & \\frac{1}{\\sqrt{p}} \\end{bmatrix} \\times \\left[ \\begin{array}{c} x_1 \\\\ x_2 \\end{array} \\right]\n",
    "$$\n",
    "\n",
    "This system is singular for $p=1$ and for $p$ in the vicinity of one is ill-conditioned.\n",
    "\n",
    "* Create two plots that show the condition number and the error between the analytic and numerical solution."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "## Iterative methods\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Algorithms based on Gaussian elimination are called exact or, more properly, direct methods because they would generate exact solutions for the linear equation $Ax = b$ after a finite number of operations, if not for rounding error. Such methods are ideal for moderately sized linear equations but may be impractical for large ones. Other methods, called iterative methods, can often be used to solve large linear equations more efficiently if the $A$ matrix is sparse, that is, if $A$ is composed mostly of zero entries. Iterative methods are designed to generate a sequence of increasingly accurate approximations to the solution of a linear equation, but they generally do not yield an exact solution after a prescribed number of steps, even in theory.\n",
    "\n",
    "\n",
    "The most widely used iterative methods for solving a linear equation $Ax = b$ are developed by choosing an easily invertible matrix $Q$ and writing the linear equation in the equivalent form\n",
    "\n",
    "$$\n",
    "Qx = b + (Q - A)x\n",
    "$$\n",
    "\n",
    "or \n",
    "\n",
    "$$\n",
    "x = Q^{-1} b + (I - Q^{-1} A)x\n",
    "$$\n",
    "\n",
    "This form of the linear equation suggests the iteration rule \n",
    "\n",
    "$$\n",
    "x^{k+1}\\leftarrow Q^{-1} b + (I - Q^{-1} A)x^{k}\n",
    "$$\n",
    "\n",
    "which, if convergent, must converge to a solution of the linear equation. Ideally, the so-called splitting matrix $Q$ will satisfy two criteria. First, $Q^{-1}b$ and $Q^{-1} A$ should be relatively easy to compute. This criterion is met if $Q$ is either diagonal or triangular. There are two popular approaches:\n",
    "\n",
    "* The **Gauss-Seidel** method sets $Q$ equal to the upper triangular matrix formed from the upper triangular elements of $A$.\n",
    "\n",
    "* The **Gauss-Jacobi** method sets $Q$ equal to the diagonal matrix formed from the diagonal entries of $A$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "??gauss_seidel"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### _Exercise_ \n",
    "\n",
    "* Implement the Gauss-Jacobi method."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "from linear_solutions_tests import gauss_jacobi  # noqa: E402"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's conclude with an economic application as outlined in Judd (1998). Suppose we  have the following inverse demand function $p = 10 - q$ and the following supply curve $q = p / 2 +1$. Equilibrium is where supply equals demand and thus we need to solve the following linear system.\n",
    "\n",
    "$$\n",
    "\\left[\\begin{array}{c} 10 \\\\ -2 \\end{array}\\right] = \\begin{bmatrix} 1 & 1\\\\ 1 & -2\\end{bmatrix} \\times \\left[ \\begin{array}{c} p \\\\ q \\end{array} \\right]\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "A, b, x_true = get_inverse_demand_problem()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can compre the two solution approaches and make sure that they in fact give the same result."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "x_seidel = gauss_seidel(A, b)\n",
    "x_jacobi = gauss_jacobi(A, b)\n",
    "np.testing.assert_almost_equal(x_seidel, x_jacobi)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "## Resources\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Software\n",
    "\n",
    "* **The PARDISO Solver Project**: https://www.pardiso-project.org\n",
    "\n",
    "* **LAPACK — Linear Algebra PACKage**: http://www.netlib.org/lapack\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### References\n",
    "\n",
    "- Robert Johansson. *Numerical Python: scientific computing and data science applications with NumPy, SciPy and Matplotlib*. Apress, 2018.\n",
    "\n",
    "- William H Press, Brian P Flannery, Saul A Teukolsky, and William T Vetterling. *Numerical recipes: The art of scientific computing*. Cambridge University Press, 1986.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}