Exercise 2 Solution.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercise 2 - Implementing the LU Decomposition"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Your exercise today is straight forward. Write a function `my_lu(A)`, which takes a matrix $A$ and returns the factors $L$ and $U$ of the LU decomposition without pivoting.\n",
    "\n",
    "* Try to vectorize the code, that means the inner-most operation of subtracting a multiple of one row from another row should be performed by operating on the whole slice of the corresponding matrix. Overall, your code should only require a double-for loop.\n",
    "* Test your function by computing $\\|A-LU\\|$. If this value is in the order of machine precision for various examples your algorithm should be correct.\n",
    "* Compare the speed of your LU decomposition with the speed of the `scipy.linalg.lu` function.\n",
    "* You may want to use the following skeleton to write your LU function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "def lu(A):\n",
    "    \n",
    "    import numpy as np\n",
    "    \n",
    "    # Return an error if matrix is not square\n",
    "    if not A.shape[0]==A.shape[1]:\n",
    "        raise ValueError(\"Input matrix must be square\")\n",
    "        \n",
    "    n = A.shape[0] # The number of rows/columns of A\n",
    "    \n",
    "    L = np.zeros((n,n),dtype='float64') # Reserve space for L\n",
    "    U = np.zeros((n,n),dtype='float64') # Reserve space for U\n",
    "    U[:] = A # Copy A into U as we do not want to modify A\n",
    "    np.fill_diagonal(L,1) # fill the diagonal of L with 1\n",
    "    \n",
    "    for i in range(n-1):\n",
    "        # The outer iteration \n",
    "        for j in range(i+1,n):\n",
    "            L[j,i] = U[j,i]/U[i,i]\n",
    "            U[j,i:] = U[j,i:]-L[j,i]*U[i,i:]\n",
    "            U[j,i] = 0 # Explicitly set to zero, not really necessary later on\n",
    "    return (L,U)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Implementation test\n",
      "The residual is 1.9058066906814195e-12\n",
      "===========\n",
      "53.8 ms ± 3.97 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)\n",
      "578 µs ± 65.8 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)\n",
      "Time for my LU decomposition is 0.05074618400012696 s\n",
      "Time for Scipy LU decomposition is 0.0004466583260000334 s\n",
      "Speed ratio: 113.61298121222791\n"
     ]
    }
   ],
   "source": [
    "# Test implementation\n",
    "import numpy as np\n",
    "from scipy.linalg import lu as scipy_lu\n",
    "\n",
    "n = 200\n",
    "A = np.random.rand(n,n)\n",
    "L,U = lu(A)\n",
    "residual = np.linalg.norm(A-np.dot(L,U))\n",
    "print(\"Implementation test\")\n",
    "print(\"The residual is {0}\".format(residual))\n",
    "print(\"===========\")\n",
    "\n",
    "# Now do some timing comparisons\n",
    "\n",
    "runtime_my_lu = %timeit -o lu(A)\n",
    "runtime_scipy_lu = %timeit -o scipy_lu(A)\n",
    "\n",
    "\n",
    "print(\"Time for my LU decomposition is {0} s\".format(runtime_my_lu.best))\n",
    "\n",
    "print(\"Time for Scipy LU decomposition is {0} s\".format(runtime_scipy_lu.best))\n",
    "\n",
    "print(\"Speed ratio: {0}\".format(runtime_my_lu.best / runtime_scipy_lu.best))\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "On my computer for larger dimensions Scipy is almost 100 times faster than our own implementation of the LU decomposition. We will later learn more details of why this is the case."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Note**: The above way of measuring times is complicated but very accurate. A simpler but slightly less accurate way of measuring execution times is as follows:\n",
    "\n",
    "    import time\n",
    "    start_time = time.time()\n",
    "    # your code\n",
    "    elapsed_time = time.time() - start_time"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python [conda env:dev]",
   "language": "python",
   "name": "conda-env-dev-py"
  },
  "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.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}