Advanced Numerical Linear Algebra/week_7_how_to_implement_linear_regression.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## How to implement Linear Regression\n",
    "\n",
    "        We will look into how we could write our own implementation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn import datasets, linear_model, metrics\n",
    "from sklearn.model_selection import train_test_split\n",
    "from sklearn.preprocessing import PolynomialFeatures\n",
    "import math, scipy , numpy as np\n",
    "from scipy import linalg"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.set_printoptions(precision=6)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "data = datasets.load_diabetes()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "feature_names=['age', 'sex', 'bmi', 'bp', 's1', 's2', 's3', 's4', 's5', 's6']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "trn,test,y_trn,y_test = train_test_split(data.data, data.target, test_size=0.2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "((353, 10), (89, 10))"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "trn.shape, test.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "def regr_metrics(act, pred):\n",
    "    return (math.sqrt(metrics.mean_squared_error(act, pred)), \n",
    "     metrics.mean_absolute_error(act, pred))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## How did sklearn do it?\n",
    "\n",
    "        By checking the source code, you can see that in the dense case, it calls scipy.linalg.lstqr , which is calling the LAPACK method:\n",
    "\n",
    "        1. gelsd : uses SVD and a divide and conquer method\n",
    "        2. gelsy : uses QR factorization\n",
    "        3. gelss : uses SVD\n",
    "\n",
    "\n",
    "    Scipy Sparse Least Squares\n",
    "    It uses an iterative method called Golub and Kahan bidiagonalization\n",
    "    Preconditioning is another way to reduce the number of iterations. If it is possible to solve a related system M*x = b, efficiently, where M approximates A in some helpful way, LSQR may converge more rapidly on the system"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### linalg.lstqr\n",
    "        The sklearn implementation handeled adding a constant term for us"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "trn_int = np.c_[trn, np.ones(trn.shape[0])]\n",
    "test_int = np.c_[test, np.ones(test.shape[0])]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "        Since linalg.lstsq lets us specify which LAPACK routine we want to use, lets try them all and do some timing comparisons:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The slowest run took 5.14 times longer than the fastest. This could mean that an intermediate result is being cached.\n",
      "607 µs ± 510 µs per loop (mean ± std. dev. of 7 runs, 1 loop each)\n"
     ]
    }
   ],
   "source": [
    "%timeit coef, _,_,_ = linalg.lstsq(trn_int, y_trn, lapack_driver=\"gelsd\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1.24 ms ± 424 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)\n"
     ]
    }
   ],
   "source": [
    "%timeit coef, _,_,_ = linalg.lstsq(trn_int, y_trn, lapack_driver=\"gelsy\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The slowest run took 4.58 times longer than the fastest. This could mean that an intermediate result is being cached.\n",
      "587 µs ± 366 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)\n"
     ]
    }
   ],
   "source": [
    "%timeit coef, _,_,_ = linalg.lstsq(trn_int, y_trn, lapack_driver=\"gelss\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Naive Solution\n",
    "    Recall that we want to find x that minimizes:\n",
    "$$ ||Ax-b||_2 $$\n",
    "    Another way to think about this is that we are interested in where vector b is closest to the subspace spanned by A. This is the projection of b onto A. Since b-Ax must be perpendicular to the subspace spanned by A, we see that\n",
    "\n",
    "$$ A^T (b - Ax) = 0 $$\n",
    "\n",
    "    This leads us to the normal equations\n",
    "$$ x = (A^T A)^{-1} A^Tb $$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "def ls_naive(A, b):\n",
    "     return np.linalg.inv(A.T @ A) @ A.T @ b"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The slowest run took 13.44 times longer than the fastest. This could mean that an intermediate result is being cached.\n",
      "128 ms ± 106 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)\n"
     ]
    }
   ],
   "source": [
    "%timeit coeffs_naive = ls_naive(trn_int, y_trn)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(56.32960590877427, 47.82386741198034)"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "coeffs_naive = ls_naive(trn_int, y_trn)\n",
    "regr_metrics(y_test, test_int @ coeffs_naive)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Normal Equations (Cholesky)\n",
    "\n",
    "    Normal Equations :\n",
    "\n",
    "$$ A^T Ax = A^T b $$\n",
    "\n",
    "If A has full rank, the pseudo-inverse $(A^TA)^{-1} A^T$ is a square, hermitian positive definite matrix. The satndard way of solving such a system is Cholesky Factorization, which dfinds iupper-triangular R, s.t $A^T A = R^TR$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "A = trn_int"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "b = y_trn"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "AtA = A.T @ A\n",
    "Atb = A.T @ b"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "R = scipy.linalg.cholesky(AtA)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 0.912 ,  0.1804,  0.1842,  0.3173,  0.2687,  0.2217, -0.0509,\n",
       "         0.1901,  0.2657,  0.3002, -0.1708],\n",
       "       [ 0.    ,  0.8748,  0.0565,  0.1561, -0.0118,  0.1097, -0.3566,\n",
       "         0.2886,  0.0754,  0.1325, -0.1085],\n",
       "       [ 0.    ,  0.    ,  0.8607,  0.307 ,  0.169 ,  0.1812, -0.2942,\n",
       "         0.2979,  0.3477,  0.3038, -0.2325],\n",
       "       [ 0.    ,  0.    ,  0.    ,  0.7696,  0.1211,  0.0471,  0.0286,\n",
       "         0.0233,  0.1789,  0.163 ,  0.0297],\n",
       "       [ 0.    ,  0.    ,  0.    ,  0.    ,  0.8123,  0.7239,  0.1359,\n",
       "         0.3783,  0.33  ,  0.1165, -1.2538],\n",
       "       [ 0.    ,  0.    ,  0.    ,  0.    ,  0.    ,  0.3846, -0.4024,\n",
       "         0.2558, -0.3441, -0.0579, -0.7216],\n",
       "       [ 0.    ,  0.    ,  0.    ,  0.    ,  0.    ,  0.    ,  0.6497,\n",
       "        -0.4917, -0.5276, -0.1386,  0.3961],\n",
       "       [ 0.    ,  0.    ,  0.    ,  0.    ,  0.    ,  0.    ,  0.    ,\n",
       "         0.2804, -0.018 ,  0.0037, -0.3866],\n",
       "       [ 0.    ,  0.    ,  0.    ,  0.    ,  0.    ,  0.    ,  0.    ,\n",
       "         0.    ,  0.2952,  0.0783, -0.3811],\n",
       "       [ 0.    ,  0.    ,  0.    ,  0.    ,  0.    ,  0.    ,  0.    ,\n",
       "         0.    ,  0.    ,  0.7238, -0.0866],\n",
       "       [ 0.    ,  0.    ,  0.    ,  0.    ,  0.    ,  0.    ,  0.    ,\n",
       "         0.    ,  0.    ,  0.    , 18.7177]])"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.set_printoptions(suppress=True, precision=4)\n",
    "R"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "5.684557367178372e-14"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.linalg.norm(AtA - R.T @ R)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$ A^T Ax = A^T b $$\n",
    "$$ R^T Rx = A^T b $$\n",
    "$$ R^T w = A^T b $$\n",
    "$$ Rx = w $$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "w = scipy.linalg.solve_triangular(R, Atb, lower=False, trans='T')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "7.277126723879252e-12"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.linalg.norm(R.T @ w - Atb)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [],
   "source": [
    "coeffs_chol = scipy.linalg.solve_triangular(R, w, lower=False)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.3048672769382318e-13"
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.linalg.norm(R @ coeffs_chol - w)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [],
   "source": [
    "def ls_chol(A, b):\n",
    "    R = scipy.linalg.cholesky(A.T @ A)\n",
    "    w = scipy.linalg.solve_triangular(R, A.T @ b, trans='T')\n",
    "    return scipy.linalg.solve_triangular(R, w)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The slowest run took 8.18 times longer than the fastest. This could mean that an intermediate result is being cached.\n",
      "120 ms ± 65.9 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)\n"
     ]
    }
   ],
   "source": [
    "%timeit coeffs_chol = ls_chol(trn_int, y_trn)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(56.3296059087742, 47.82386741198027)"
      ]
     },
     "execution_count": 28,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "coeffs_chol = ls_chol(trn_int, y_trn)\n",
    "regr_metrics(y_test, test_int @ coeffs_chol)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### QR Factorization\n",
    "$$ Ax = b $$\n",
    "\n",
    "$$ A = QR$$\n",
    "$$ QRx = b $$\n",
    "$$ Rx = Q^T b$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [],
   "source": [
    "def ls_qr(A,b):\n",
    "    Q, R = scipy.linalg.qr(A, mode='economic')\n",
    "    return scipy.linalg.solve_triangular(R, Q.T @ b)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "409 µs ± 103 µs per loop (mean ± std. dev. of 7 runs, 1 loop each)\n"
     ]
    }
   ],
   "source": [
    "%timeit coeffs_qr = ls_qr(trn_int, y_trn)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(56.32960590877418, 47.82386741198025)"
      ]
     },
     "execution_count": 31,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "coeffs_qr = ls_qr(trn_int, y_trn)\n",
    "regr_metrics(y_test, test_int @ coeffs_qr)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### SVD    \n",
    "\n",
    "$$ Ax = b $$\n",
    "$$ A = U \\sum V$$\n",
    "$$ \\sum V x = U^T b $$\n",
    "$$ \\sum w = U^T b $$\n",
    "$$ x = V^T w $$\n",
    "\n",
    "        SVD gives a pseudo-inverse"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [],
   "source": [
    "def ls_svd(A,b):\n",
    "    m, n = A.shape\n",
    "    U, sigma, Vh = scipy.linalg.svd(A, full_matrices=False, lapack_driver='gesdd')\n",
    "    w = (U.T @ b)/ sigma\n",
    "    return Vh.T @ w"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "550 µs ± 149 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)\n"
     ]
    }
   ],
   "source": [
    "%timeit coeffs_svd = ls_svd(trn_int, y_trn)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2.44 ms ± 775 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)\n"
     ]
    }
   ],
   "source": [
    "%timeit coeffs_svd = ls_svd(trn_int, y_trn)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(56.32960590877416, 47.82386741198024)"
      ]
     },
     "execution_count": 35,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "coeffs_svd = ls_svd(trn_int, y_trn)\n",
    "regr_metrics(y_test, test_int @ coeffs_svd)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Random Sketching Technique for Least Squares Regression\n",
    "\n",
    "        1. Sample a r x n random matrix\n",
    "        2. Compute SA and Sb\n",
    "        3. Find Exact Solution x to regression SA x = Sb"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Timing Comparison"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [],
   "source": [
    "import timeit\n",
    "import pandas as pd"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {},
   "outputs": [],
   "source": [
    "def scipylstq(A, b):\n",
    "    return scipy.linalg.lstsq(A,b)[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {},
   "outputs": [],
   "source": [
    "row_names = ['Normal Eqns- Naive',\n",
    "             'Normal Eqns- Cholesky', \n",
    "             'QR Factorization', \n",
    "             'SVD', \n",
    "             'Scipy lstsq']\n",
    "\n",
    "name2func = {'Normal Eqns- Naive': 'ls_naive', \n",
    "             'Normal Eqns- Cholesky': 'ls_chol', \n",
    "             'QR Factorization': 'ls_qr',\n",
    "             'SVD': 'ls_svd',\n",
    "             'Scipy lstsq': 'scipylstq'}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {},
   "outputs": [],
   "source": [
    "m_array = np.array([100, 1000, 10000])\n",
    "n_array = np.array([20, 100, 1000])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [],
   "source": [
    "index = pd.MultiIndex.from_product([m_array, n_array], names=['# rows', '# cols'])\n",
    "pd.options.display.float_format = '{:,.6f}'.format\n",
    "df = pd.DataFrame(index=row_names, columns=index)\n",
    "df_error = pd.DataFrame(index=row_names, columns=index)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# %%prun\n",
    "for m in m_array:\n",
    "    for n in n_array:\n",
    "        if m >= n:        \n",
    "            x = np.random.uniform(-10,10,n)\n",
    "            A = np.random.uniform(-40,40,[m,n])   # removed np.asfortranarray\n",
    "            b = np.matmul(A, x) + np.random.normal(0,2,m)\n",
    "            for name in row_names:\n",
    "                fcn = name2func[name]\n",
    "                t = timeit.timeit(fcn + '(A,b)', number=5, globals=globals())\n",
    "                df.set_value(name, (m,n), t)\n",
    "                coeffs = locals()[fcn](A, b)\n",
    "                reg_met = regr_metrics(b, A @ coeffs)\n",
    "                df_error.set_value(name, (m,n), reg_met[0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_error"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "store = pd.HDFStore('least_squares_results.h5')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "store['df'] = df"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Conditioning & Stability\n",
    "\n",
    "        Condition Number is a measure of how small changes to the input cause the output to change.The relative condition number is defined by \n",
    "\n",
    "$$ \\kappa = \\text{sup} \\frac{||\\delta f ||}{|| f(x)||} / \\frac{||\\delta x||}{||x||} $$\n",
    "\n",
    "\n",
    "        Conditioning : perturbation behavior of a mathematical problem\n",
    "        Stability : perturabation  behavior of an algorithm used to solve that problem on a computer\n",
    "\n",
    "\n",
    "### Conditioning Example\n",
    "        The problem of computing eigenvalues of a non-symmetric matrix is often ill-conditioned"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "metadata": {},
   "outputs": [],
   "source": [
    "A = [[1, 1000], [0, 1]]\n",
    "B = [[1, 1000], [0.001, 1]]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "metadata": {},
   "outputs": [],
   "source": [
    "wA, vrA = scipy.linalg.eig(A)\n",
    "wB, vrB = scipy.linalg.eig(B)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(array([1.+0.j, 1.+0.j]), array([2.+0.j, 0.+0.j]))"
      ]
     },
     "execution_count": 50,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "wA, wB"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Condition Number of a Matrix\n",
    "\n",
    "The product $||A|| || A^{-1}|| $ come up so often it has its own name: the condition number of A. Note that normalyy we talk about the conditioning of problems, not matrices.\n",
    "\n",
    "        The condition number of A relates to :\n",
    "        1. computing b given A and x in Ax = b\n",
    "        2. computing x given A and b in Ax = b"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Matrix Inversion is Unstable"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.linalg import hilbert"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 52,
   "metadata": {},
   "outputs": [],
   "source": [
    "n = 14\n",
    "A = hilbert(n)\n",
    "x = np.random.uniform(-10,10,n)\n",
    "b = A @ x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "metadata": {},
   "outputs": [],
   "source": [
    "A_inv = np.linalg.inv(A)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "4.7177114010624965"
      ]
     },
     "execution_count": 54,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.linalg.norm(np.eye(n) - A @ A_inv)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 55,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "9.274324459181185e+17"
      ]
     },
     "execution_count": 55,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.linalg.cond(A)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 56,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 1.    ,  0.    ,  0.    ,  0.0003, -0.0008, -0.0033,  0.0068,\n",
       "        -0.0587, -0.0236,  0.3125, -0.9121, -0.1484,  0.0065, -0.0569],\n",
       "       [-0.    ,  1.    , -0.0001,  0.0014, -0.0126,  0.0785, -0.1674,\n",
       "         0.2113, -0.2763,  0.625 , -0.9451, -0.1367,  0.1234, -0.0101],\n",
       "       [ 0.    , -0.    ,  1.    , -0.0002,  0.002 ,  0.002 ,  0.0156,\n",
       "        -0.0156,  0.3125,  0.125 , -0.2812, -0.0312,  0.    , -0.0127],\n",
       "       [ 0.    ,  0.    ,  0.    ,  0.9999,  0.0047, -0.0407,  0.0894,\n",
       "        -0.6306,  0.9308, -1.4375,  1.2176, -1.0351,  0.4388, -0.0656],\n",
       "       [ 0.    , -0.    ,  0.    , -0.    ,  1.0011,  0.0031,  0.0148,\n",
       "        -0.1498,  0.4747,  0.    ,  0.2246, -0.0547,  0.0454, -0.0039],\n",
       "       [ 0.    , -0.    ,  0.    , -0.0002,  0.0062,  0.964 ,  0.1136,\n",
       "        -0.3853,  0.6504, -0.4688,  0.5615, -0.2886,  0.2377, -0.0267],\n",
       "       [ 0.    , -0.    ,  0.    , -0.0005,  0.007 , -0.0264,  1.0804,\n",
       "        -0.2277,  0.6103, -0.2188,  0.1478, -0.1211,  0.0712, -0.0013],\n",
       "       [-0.    , -0.    ,  0.    ,  0.0002,  0.0009, -0.0055,  0.0384,\n",
       "         0.9296,  0.2239, -0.1875,  0.3086, -0.1875,  0.1137, -0.0081],\n",
       "       [ 0.    , -0.    ,  0.    , -0.0007,  0.0071, -0.0468,  0.1594,\n",
       "        -0.6134,  2.0899, -1.375 ,  1.2291, -0.8576,  0.336 , -0.0507],\n",
       "       [-0.    , -0.    ,  0.    , -0.0001,  0.0026, -0.0024,  0.032 ,\n",
       "        -0.0885,  0.1898,  1.    , -0.1539, -0.0812,  0.0166,  0.0062],\n",
       "       [ 0.    , -0.    , -0.    , -0.0001,  0.0012, -0.0099,  0.0268,\n",
       "        -0.1045,  0.1998, -0.375 ,  1.0435, -0.1308,  0.078 , -0.0024],\n",
       "       [ 0.    , -0.    ,  0.    , -0.0004,  0.0057, -0.0331,  0.114 ,\n",
       "        -0.3636,  0.9706, -0.9375,  0.6756,  0.4871,  0.2194, -0.024 ],\n",
       "       [ 0.    , -0.    ,  0.    , -0.0003,  0.0026, -0.0082, -0.0202,\n",
       "         0.1145,  0.0151,  0.5938, -0.6964,  0.3672,  0.849 ,  0.0227],\n",
       "       [-0.    ,  0.    , -0.    ,  0.0003, -0.0026,  0.0154, -0.0878,\n",
       "         0.1345, -0.3189,  0.4688, -0.9179,  0.375 , -0.0765,  1.0153]])"
      ]
     },
     "execution_count": 56,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A @ A_inv"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 57,
   "metadata": {},
   "outputs": [],
   "source": [
    "row_names = ['Normal Eqns- Naive',\n",
    "             'QR Factorization', \n",
    "             'SVD', \n",
    "             'Scipy lstsq']\n",
    "\n",
    "name2func = {'Normal Eqns- Naive': 'ls_naive', \n",
    "             'QR Factorization': 'ls_qr',\n",
    "             'SVD': 'ls_svd',\n",
    "             'Scipy lstsq': 'scipylstq'}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 58,
   "metadata": {},
   "outputs": [],
   "source": [
    "pd.options.display.float_format = '{:,.9f}'.format\n",
    "df = pd.DataFrame(index=row_names, columns=['Time', 'Error'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "        Even if A is incredibly sparse, A inverse is generally dense. For large matrices , A inverse could be so dense as to not fit in memory\n",
    "\n",
    "\n",
    "## Runtime\n",
    "1. Matrix Inversion : $2n^3$\n",
    "2. Matrix Multiplication : $n^3$\n",
    "3. Cholesky : $\\frac{1}{3} n^3$\n",
    "4. QR, Gram Schmidt : $ 2mn^2$, $m \\ge n$\n",
    "5. QR, Householder : $2mn^2 - \\frac{2}{3} n^3$\n",
    "6. Solving a triangular system : $n^2$"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3.7.0 64-bit",
   "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.7.0"
  },
  "orig_nbformat": 4,
  "vscode": {
   "interpreter": {
    "hash": "c0f2ef85a013b3e049ebe5ad4c39b42bd34e85c4b0686b7bef6ce5ff62aa91e1"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}