week_4_self_check_solutions.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Week 4 Self Check Questions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Question 1:\n",
    "\n",
    "Write a Python function to implement the LU Decomposition with partial pivoting. You can take the following skeleton as starting point."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Residual with Scipy LU: 3.0356551825401023e-15\n",
      "Residual with partial pivot LU: 8.041205793198706e-13\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "from scipy.linalg import solve_triangular\n",
    "from scipy.linalg import lu\n",
    "\n",
    "def partial_pivot_lu(A):\n",
    "    \"\"\"Implementation of a pivoted LU decomposition.\"\"\"\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",
    "    p = np.arange(n) # The permutation vector\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",
    "        # IMPLEMENT THE ROW EXCHANGE FOR PARTIAL PIVOTING HERE.\n",
    "        # STORE THE PERMUTATION IN THE p VECTOR.\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\n",
    "    # The vector p stores the forward permutations.\n",
    "    # The following code takes an identity matrix and swaps the rows accordingly.\n",
    "    P = np.eye(n)[p,:] \n",
    "    return (P,L,U)\n",
    "\n",
    "n = 10\n",
    "\n",
    "A = np.random.randn(n, n)\n",
    "b = np.random.randn(n)\n",
    "\n",
    "# Example solution with Scipy\n",
    "# Note that Scipy computes A = PLU and not PA=LU.\n",
    "P, L, U = lu(A)\n",
    "x_compare = solve_triangular(\n",
    "    U, solve_triangular(L, P.T @ b, lower=True))\n",
    "# Print the residual\n",
    "print(f\"Residual with Scipy LU: {np.linalg.norm(b - A @ x_compare)}\")\n",
    "\n",
    "P2, L2, U2 = partial_pivot_lu(A)\n",
    "x_pivot = solve_triangular(\n",
    "    U2, solve_triangular(L2, P2 @ b, lower=True))\n",
    "print(f\"Residual with partial pivot LU: {np.linalg.norm(b - A @ x_pivot)}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Solution:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Residual with Scipy LU: 1.1123893155135927e-15\n",
      "Residual with partial pivot LU: 2.7294119545946647e-15\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "from scipy.linalg import solve_triangular\n",
    "from scipy.linalg import lu\n",
    "\n",
    "def partial_pivot_lu(A):\n",
    "    \"\"\"Implementation of a pivoted LU decomposition.\"\"\"\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",
    "    p = np.arange(n) # The permutation vector\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",
    "        # Permute the rows to ensure the element with largest magnitude is on top\n",
    "        max_index = i+np.argmax(np.abs(U[i:,i]))\n",
    "        U[[i,max_index],:] = U[[max_index,i],:]\n",
    "        L[[i,max_index],:i] = L[[max_index,i],:i]\n",
    "        p[[i,max_index]] = p[[max_index,i]]\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\n",
    "    # The vector p stores the forward permutations.\n",
    "    # The following code takes an identity matrix and swaps the rows accordingly.\n",
    "    P = np.eye(n)[p,:] \n",
    "    return (P,L,U)\n",
    "\n",
    "n = 10\n",
    "\n",
    "A = np.random.randn(n, n)\n",
    "b = np.random.randn(n)\n",
    "\n",
    "# Example solution with Scipy\n",
    "# Note that Scipy computes A = PLU and not PA = LU\n",
    "P, L, U = lu(A)\n",
    "x_compare = solve_triangular(\n",
    "    U, solve_triangular(L, P.T @ b, lower=True))\n",
    "# Print the residual\n",
    "print(f\"Residual with Scipy LU: {np.linalg.norm(b - A @ x_compare)}\")\n",
    "\n",
    "P2, L2, U2 = partial_pivot_lu(A)\n",
    "x_pivot = solve_triangular(\n",
    "    U2, solve_triangular(L2, P2 @ b, lower=True))\n",
    "print(f\"Residual with partial pivot LU: {np.linalg.norm(b - A @ x_pivot)}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Question 1:\n",
    "\n",
    "Review again why the LU Decomposition with partial pivoting is not backward stable in general. What is the condition for the LU with partial pivoting to be backward stable?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Solution: See the corresponding Jupyter Notebook for an in-depth explanation of the backward stability of LU decomposition. The LU Decomposition with partial pivoting is backward stable if the growth factor $\\rho$ does not become large."
   ]
  }
 ],
 "metadata": {
  "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
}