Assignment 1-done.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# TMA4215 - Assignment 1, LU-factorisation with pivoting\n",
    "\n",
    "\n",
    "**Deadline:** Wednesday September 2nd, 11:59PM.\n",
    "\n",
    "**Requirements for approval:** You need to do two things: 1) Upload a Jupyter Notebook on Blackboard that contains the complete solution to the assigment and 2) Answer the Control question form of the assignment, you should make sure that at least some of the answers are correct.\n",
    "\n",
    "**Supervision** of this assignment is digital. On Thursday August 27, 1600- the TA will be available for answering questions asked on Piazza and for one-to-one video conferencing on the platform Whereby.com.\n",
    "At any other time you may ask questions on Piazza, but please be a little patient then, it may take some time before your question is answered.\n",
    "\n",
    "\n",
    "This assignment is our way of compensating for not lecturing Gaussian elimination and LU-factorisation with pivoting. Many of you have seen this algorithm before. If your background on this method is weak, you can read up on chapters 3.2-3.5 in the text book.\n",
    "\n",
    "In this assignment you are to implement and test a Python-program for doing LU-factorisation with pivoting.\n",
    "In large simulations in real life it is often recommended to use as little memory as possible, and also avoid unnecessary operations that are time consuming. From this point of departure we would like you to make a code that takes the coefficient matrix $A$ as input and returns a representation of the matrices $P$, $L$, and $U$ which solve the equation\n",
    "\n",
    "$$\n",
    "     PA = LU\n",
    "$$\n",
    "\n",
    "where $P$ is a permutaton matrix, $L$ is a lower-triangular matrix with unit diagonal, and $U$ is an upper-triangular matrix. We *represent* the matrices in question as follows: The permutation matrix $P$ is $n\\times n$, but is represented as a vector  $\\mathtt{P}$ such that row number $k$ in $P$ is the canonical unit vector $e_{\\mathtt{P}_k}$. Let us illustrate this by an example\n",
    "\n",
    "$$\n",
    "\\mathtt{P}=\n",
    "\\left[\n",
    "\\begin{array}{r} 3 \\\\ 1 \\\\ 2 \\end{array}\n",
    "\\right]\\quad\\Rightarrow\\quad\n",
    "P=\\left[\n",
    "\\begin{array}{ccc}\n",
    "0 & 0 & 1 \\\\ 1 & 0 & 0 \\\\ 0 & 1 & 0\n",
    "\\end{array}\n",
    "\\right]\n",
    "$$\n",
    "\n",
    "We stipulate that a Python function takes a two-dimensional numpy-array $\\mathtt{A}$ as input, and returns\n",
    "an *over-written* $\\mathtt{A}$ which contains $L$ and $U$ in the following sense upon return:\n",
    "\n",
    "$$\n",
    "\\begin{array}{ll}\n",
    "\\mathtt{A}[\\mathtt{P}[i],j] = L_{ij} & \\text{for}\\ i<j \\\\\n",
    "\\mathtt{A}[\\mathtt{P}[i],j] = U_{ij} & \\text{for}\\ i\\geq j\n",
    "\\end{array}\n",
    "$$\n",
    "\n",
    "That $L$ has 1 on the diagonal is always the case, so the diagonal of $L$ needs not be stored. The remaining elements of $L$ and $U$ are zero and need not be stored either. The algorithm can be formulated as follows (compare to text book):\n",
    "\n",
    "\n",
    "- Input: $A$ of size $n\\times n$\n",
    "- Initialisation\n",
    "    * Let $P_i = i,\\ i=0,\\ldots,n-1$ be a vector (array) with $n$ components\n",
    "- for $k$ **in** range(n-1):\n",
    "    1. Find index $P_\\ell$ such that $|\\mathtt{A}_{P_\\ell,k}|=\\max_{k\\leq i \\leq n-1} |\\mathtt{A}_{P_i,k}|$, i.e. scan column $k$ from the diagonal and down for the largest element in absolute value. \n",
    "    2. Swap $P_k$ by $P_\\ell$.\n",
    "    3. Find multipliers $A_{P_i,k}\\leftarrow \\frac{A_{P_i,k}}{A_{P_k,k}},\\ i=k+1,\\ldots,n-1$.\n",
    "    4. Perform elimination, i.e. $A_{P_i,j}\\leftarrow A_{P_i,j}-A_{P_i,k}\\cdot A_{P_k,j},\\ i,j=k+1,\\ldots,n-1$\n",
    "- Output: A,P\n",
    "\n",
    "**Comment:** In practice there are of course off the shelf implementations both in Python libraries and other places for solving linear systems in an optimal way with respect to accuracy and efficiency. It also happens often that many elements of the coefficient matrix are zero and that can be exploited in various different ways. The solver you make is rather general. Irrespective of whether the code you implement is standard software, it is a useful experience to have written such a program yourself at least once, so that you understand how it works, the kind of errors that may occur etc, and you will gain understanding for instance in error analysis.\n",
    "\n",
    "\n",
    "\n",
    "**Problem 1** Write a function for LU-factorisation with row-wise pivoting as indicated above.\n",
    "A template could be\n",
    "\n",
    "\n",
    "    def mylu(A):\n",
    "   \n",
    "    \n",
    "and it should return the pivot vector (permutation vector) $\\mathtt{P}$, and and over-written  version of $A$. You can also choose to copy $A$ into some other matrix $\\mathtt{LU}$ from the beginning using e.g. \n",
    "\n",
    "    LU = A.copy()\n",
    "\n",
    "and write into and return this matrix in order to save the input matrix $A$ \n",
    "\n",
    "Use the algorithm described above. See hints about indexing and useful numpy-functions that can be used below.\n",
    "\n",
    "**Problem 2** Combine your function (mylu) with the functions for forward and backward substitution given below for computing solutions to linear systems $Ax=b$. Test it out by using $A$ and $b$ from the function getAb() below.\n",
    "\n",
    "**Control question 1:** Give the permutation vector $\\mathtt{P}$ from this numerical test (multiple choice)\n",
    "\n",
    "**Control question 2:** Give the first component of the intermediate result $c$ (where $Lc=Pb$) (multiple choice)\n",
    "\n",
    "**Control question 3:** Give the last component in the final answer $x$ (where $Ax=b$) with the given example.\n",
    "\n",
    "**Approval**\n",
    "A Jupyter notebook file with code that solves Problems 1 and 2 isto be uploaded in Blackboard and you need to answer the control questions.\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Codes for forward and bacwards substitution\n",
    "To make life easier for you, we give you two functions that can be called in order to solve Ax=b after one has LU-factorised with pivoting. It presumed that a matrix LU has been computed, and a permutation vector P."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "def forward_subs(LU,P,b):\n",
    "    ''' Forward substitution algorithm\n",
    "    Input:\n",
    "        LU contains both L and U, even if only L is needed here \n",
    "        P Integer permutation vector\n",
    "        b Vector with right hand side in the problem to be solved\n",
    "    Output:\n",
    "        c The solution to the linear lower triangular system Lc=b \n",
    "    '''\n",
    "    n, m = LU.shape\n",
    "    Pb = b[P]\n",
    "    c = np.zeros(n)\n",
    "    c[0] = Pb[0]\n",
    "    for k in range(1,n):\n",
    "        c[k] = Pb[k] - LU[P[k],0:k] @ c[0:k]\n",
    "        \n",
    "    return c\n",
    "\n",
    "def backward_subs(LU,P,c):\n",
    "    ''' Backward substitution algorithm\n",
    "    Input:\n",
    "        LU contains both L and U, even though just U is needed here\n",
    "        P Integer permutation vector\n",
    "        c Vector containing right hand side, i.e. the function solves Ux=c\n",
    "    Output:\n",
    "        x Solution to the linear upper triangular system Ux = c\n",
    "    '''\n",
    "    n,m = LU.shape\n",
    "    x = np.zeros(n)\n",
    "    x[n-1] = c[n-1]/LU[P[n-1],n-1]\n",
    "    for k in range(n-1,0,-1):\n",
    "        x[k-1] = (c[k-1]-LU[P[k-1],k:] @ x[k:])/LU[P[k-1],k-1]\n",
    "        \n",
    "    return x"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Problem 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "def myLU(A):\n",
    "    LU = A.copy()\n",
    "    n = LU.shape[0]\n",
    "    P = [i for i in range(n)]\n",
    "    for k in range(n-1):\n",
    "        piv = np.argmax(abs(LU[P[k:], k])) #1\n",
    "        val = (LU[P[k:], k])[piv] #1\n",
    "        piv += k\n",
    "        P[piv], P[k] = P[k], P[piv] #switcheroo\n",
    "        \n",
    "        mults = LU[P[k+1:],k] / LU[P[k],k]\n",
    "        LU[P[k+1:], k+1:] = LU[P[k+1:],k+1:] - np.outer(mults,LU[P[k],k+1:])\n",
    "        LU[P[k+1:], k] = mults\n",
    "    return LU, P"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Function for constructing A and b in Problem 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def getAb():\n",
    "    A=np.array([[0.3050, 0.5399, 0.9831, 0.4039, 0.1962],\n",
    "                [0.2563, -0.1986, 0.7903, 0.6807, 0.5544],\n",
    "                [0.7746, 0.6253, -0.1458, 0.1704,  0.5167],\n",
    "                [0.4406, 0.9256, 0.4361, -0.2254, 0.7784],\n",
    "                [0.4568, 0.2108, 0.6006, 0.3677, -0.8922]])\n",
    "    b=np.array([0.9876,-1.231,0.0987,-0.5544,0.7712])\n",
    "    return A,b\n",
    "    \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Problem 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Control question 1:\n",
      " [2, 3, 1, 0, 4]\n",
      "Control question 2:\n",
      " 0.0987\n",
      "Control question 3:\n",
      " -1.4130390910442732\n"
     ]
    }
   ],
   "source": [
    "A, b = getAb()\n",
    "LU, P = myLU(A)\n",
    "c = forward_subs(LU,P,b)\n",
    "x = backward_subs(LU,P,c)\n",
    "print(\"Control question 1:\\n\", P)\n",
    "print(\"Control question 2:\\n\", c[0])\n",
    "print(\"Control question 3:\\n\", x[4])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Some Python and numpy stuff\n",
    "In LU-factorisation with partial (row-wise) pivoting, one needs to make use of indirect addressin of arrays via a permutation vector $\\mathtt{P}$. Below you find some examples of how things can be done."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "A=\n",
      " [[ 1  2  3  4]\n",
      " [ 5  6  7  8]\n",
      " [ 9 10 11 12]\n",
      " [13 14 15 16]]\n",
      "\n",
      "\n",
      "Permutation vector P= [2 3 0 1]\n",
      "\n",
      "\n",
      "The rows of A sorted according to P\n",
      ", written: A[P,]\n",
      "\n",
      " [[ 9 10 11 12]\n",
      " [13 14 15 16]\n",
      " [ 1  2  3  4]\n",
      " [ 5  6  7  8]]\n",
      "\n",
      "\n",
      "Column k(= 1 ) sorted according to P, written: A[P,k]\n",
      " [10 14  2  6]\n",
      "\n",
      "\n",
      " Fetch the last part of column k in A sorted according P: A[P[k:],k]\n",
      " [14  2  6]\n",
      "\n",
      "\n",
      "***\n",
      "\n",
      "\n",
      "The numpy-function np.outer can be handy in doing item 4 above\n",
      "\n",
      "Example:\n",
      " Let x,y be\n",
      "\n",
      "x= [1 2 3] \n",
      "y= [ 1  0 -1]\n",
      "\n",
      "Then np.outer(x,y)=\n",
      " [[ 1  0 -1]\n",
      " [ 2  0 -2]\n",
      " [ 3  0 -3]]\n",
      "\n",
      "\n",
      "The numpy function argmax can be useful in seeking the pivot-element (row to be swapped)\n",
      "\n",
      "Write for instance: pivot = np.argmax(abs(A[P[k:], k])+k\n",
      " yielding pivot= 3 når k= 2\n",
      "\n",
      "\n",
      " AA=\n",
      " [[0.57885785 0.28767529 0.17731517 0.81429465 0.86021489]\n",
      " [0.7778398  0.90839834 0.10378993 0.69777731 0.1597144 ]\n",
      " [0.18696517 0.89066415 0.11685609 0.81030253 0.84750771]\n",
      " [0.4146784  0.84396208 0.02987348 0.13143525 0.06514406]\n",
      " [0.9245245  0.99149909 0.70222343 0.63021531 0.33206339]]\n"
     ]
    }
   ],
   "source": [
    "\n",
    "\n",
    "A=np.array([[1,2,3,4],[5,6,7,8],[9,10,11,12],[13,14,15,16]])\n",
    "P=np.array([2,3,0,1])\n",
    "m,n = A.shape\n",
    "print('A=\\n',A)\n",
    "print('\\n\\nPermutation vector P=',P)\n",
    "k=1\n",
    "\n",
    "print('\\n\\nThe rows of A sorted according to P\\n, written: A[P,]\\n\\n', A[P,])\n",
    "\n",
    "print('\\n\\nColumn k(=',k,') sorted according to P, written: A[P,k]\\n',A[P,k])\n",
    "print('\\n\\n Fetch the last part of column k in A sorted according P: A[P[k:],k]\\n', A[P[k:],k])\n",
    "\n",
    "\n",
    "\n",
    "print('\\n\\n***\\n\\n')\n",
    "\n",
    "print('The numpy-function np.outer can be handy in doing item 4 above\\n')\n",
    "print('Example:\\n Let x,y be\\n') \n",
    "x=np.array([1,2,3]) \n",
    "y=np.array([1,0,-1])\n",
    "print('x=',x,'\\ny=',y)\n",
    "print('\\nThen np.outer(x,y)=\\n',np.outer(x,y))\n",
    "\n",
    "k=2\n",
    "print('\\n\\nThe numpy function argmax can be useful in seeking the pivot-element (row to be swapped)\\n')\n",
    "pivot = np.argmax(abs(A[P[k:], k]))+k\n",
    "print('Write for instance: pivot = np.argmax(abs(A[P[k:], k])+k\\n yielding pivot=',pivot, 'når k=',k)\n",
    "\n",
    "AA=np.random.rand(5,5)\n",
    "print('\\n\\n AA=\\n',AA)\n",
    "\n"
   ]
  }
 ],
 "metadata": {
  "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.7.4"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": true,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}