week_1_LU_pivoting.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# I. $LU$ factorization of a square matrix"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Consider a simple naive implementation of the LU decomposition. \n",
    "\n",
    "Note that we're using the `numpy` arrays to represent matrices [do **not** use `np.matrix`]."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "def diy_lu(a):\n",
    "    \"\"\"Construct the LU decomposition of the input matrix.\n",
    "    \n",
    "    Naive LU decomposition: work column by column, accumulate elementary triangular matrices.\n",
    "    No pivoting.\n",
    "    \"\"\"\n",
    "    N = a.shape[0]\n",
    "    \n",
    "    u = a.copy()\n",
    "    L = np.eye(N)\n",
    "    for j in range(N-1):\n",
    "        lam = np.eye(N)\n",
    "        gamma = u[j+1:, j] / u[j, j]\n",
    "        lam[j+1:, j] = -gamma\n",
    "        u = lam @ u\n",
    "\n",
    "        lam[j+1:, j] = gamma\n",
    "        L = L @ lam\n",
    "    return L, u"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "6"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Now, generate a full rank matrix and test the naive implementation\n",
    "\n",
    "import numpy as np\n",
    "\n",
    "N = 6\n",
    "a = np.zeros((N, N), dtype=float)\n",
    "for i in range(N):\n",
    "    for j in range(N):\n",
    "        a[i, j] = 3. / (0.6*i*j + 1)\n",
    "\n",
    "np.linalg.matrix_rank(a)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Tweak the printing of floating-point numbers, for clarity\n",
    "np.set_printoptions(precision=3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 1.     0.     0.     0.     0.     0.   ]\n",
      " [ 1.     1.     0.     0.     0.     0.   ]\n",
      " [ 1.     1.455  1.     0.     0.     0.   ]\n",
      " [ 1.     1.714  1.742  1.     0.     0.   ]\n",
      " [ 1.     1.882  2.276  2.039  1.     0.   ]\n",
      " [ 1.     2.     2.671  2.944  2.354  1.   ]] \n",
      "\n",
      "[[  3.000e+00   3.000e+00   3.000e+00   3.000e+00   3.000e+00   3.000e+00]\n",
      " [  0.000e+00  -1.125e+00  -1.636e+00  -1.929e+00  -2.118e+00  -2.250e+00]\n",
      " [  0.000e+00   0.000e+00   2.625e-01   4.574e-01   5.975e-01   7.013e-01]\n",
      " [  0.000e+00   2.220e-16   0.000e+00  -2.197e-02  -4.480e-02  -6.469e-02]\n",
      " [  0.000e+00  -4.528e-16   0.000e+00   6.939e-18   8.080e-04   1.902e-03]\n",
      " [  0.000e+00   4.123e-16   0.000e+00  -1.634e-17   0.000e+00  -1.585e-05]] \n",
      "\n",
      "[[  0.000e+00   0.000e+00   0.000e+00   0.000e+00   0.000e+00   0.000e+00]\n",
      " [  0.000e+00   0.000e+00   0.000e+00   0.000e+00   0.000e+00   0.000e+00]\n",
      " [  0.000e+00   0.000e+00   0.000e+00   2.220e-16  -1.110e-16  -1.665e-16]\n",
      " [  0.000e+00   0.000e+00   2.220e-16  -5.551e-17  -1.665e-16  -1.665e-16]\n",
      " [  0.000e+00   0.000e+00  -1.110e-16  -1.665e-16   1.665e-16   5.551e-17]\n",
      " [  0.000e+00   0.000e+00  -1.665e-16  -1.665e-16   5.551e-17   0.000e+00]]\n"
     ]
    }
   ],
   "source": [
    "L, u = diy_lu(a)\n",
    "\n",
    "print(L, \"\\n\")\n",
    "print(u, \"\\n\")\n",
    "\n",
    "# Quick sanity check: L times U must equal the original matrix, up to floating-point errors.\n",
    "print(L@u - a)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# II. The need for pivoting"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's tweak the matrix a little bit, we only change a single element:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "a1 = a.copy()\n",
    "a1[1, 1] = 3"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Resulting matix still has full rank, but the naive LU routine breaks down."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "6"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.linalg.matrix_rank(a1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[nan nan nan nan nan nan]\n",
      " [nan nan nan nan nan nan]\n",
      " [nan nan nan nan nan nan]\n",
      " [nan nan nan nan nan nan]\n",
      " [nan nan nan nan nan nan]\n",
      " [nan nan nan nan nan nan]] [[nan nan nan nan nan nan]\n",
      " [nan nan nan nan nan nan]\n",
      " [nan nan nan nan nan nan]\n",
      " [nan nan nan nan nan nan]\n",
      " [nan nan nan nan nan nan]\n",
      " [nan nan nan nan nan nan]]\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/maxkway/anaconda3/lib/python3.7/site-packages/ipykernel_launcher.py:15: RuntimeWarning: divide by zero encountered in true_divide\n",
      "  from ipykernel import kernelapp as app\n",
      "/home/maxkway/anaconda3/lib/python3.7/site-packages/ipykernel_launcher.py:15: RuntimeWarning: invalid value encountered in true_divide\n",
      "  from ipykernel import kernelapp as app\n"
     ]
    }
   ],
   "source": [
    "l, u = diy_lu(a1)\n",
    "\n",
    "print(l, u)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Test II.1\n",
    "\n",
    "For a naive LU decomposition to work, all leading minors of a matrix should be non-zero. Check if this requirement is satisfied for the two matrices `a` and `a1`.\n",
    "\n",
    "(20% of the grade)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "def check_all_minors(a, eps=1e-30):\n",
    "    N = a.shape[0]\n",
    "    minors = np.zeros((N, 1), dtype=float)\n",
    "    for j in range(N):\n",
    "        minors[j] = np.linalg.det(a[:(j+1),:(j+1)])\n",
    "    return np.abs(minors) > eps"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ True]\n",
      " [ True]\n",
      " [ True]\n",
      " [ True]\n",
      " [ True]\n",
      " [ True]]\n",
      "[[ True]\n",
      " [False]\n",
      " [ True]\n",
      " [ True]\n",
      " [ True]\n",
      " [ True]]\n"
     ]
    }
   ],
   "source": [
    "print(check_all_minors(a))\n",
    "print(check_all_minors(a1))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Test II.2\n",
    "\n",
    "Modify the `diy_lu` routine to implement column pivoting. Keep track of pivots, you can either construct a permutation matrix, or a swap array (your choice).\n",
    "\n",
    "(40% of the grade)\n",
    "\n",
    "Implement a function to reconstruct the original matrix from a decompositon. Test your routines on the matrices `a` and `a1`.\n",
    "\n",
    "(40% of the grade)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 169,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "def lu(a):\n",
    "    N = a.shape[0]\n",
    "    u = a.copy()\n",
    "    \n",
    "    L = np.eye(N)\n",
    "    l = np.eye(N)\n",
    "    P = np.eye(N)\n",
    "    p_revert = np.eye(N)\n",
    "    for j in range(N-1):\n",
    "        \n",
    "        i = np.argmax(np.abs(u[j:,j])) #finding pivot\n",
    "        i = i + j\n",
    "        \n",
    "        p = np.eye(N) # creating per-matrix\n",
    "        p[j,j] = 0\n",
    "        p[i,i] = 0 \n",
    "        p[i,j] = 1\n",
    "        p[j,i] = 1\n",
    "        \n",
    "        u = p @ u\n",
    "        \n",
    "        lam = np.eye(N)\n",
    "        gamma = u[j+1:, j] / u[j, j]\n",
    "        lam[j+1:, j] = -gamma\n",
    "        u = lam @ u\n",
    "\n",
    "        lam[j+1:, j] = gamma\n",
    "        l = l @ p @ lam\n",
    "        P = p @ P\n",
    "        p_revert = p_revert @ p\n",
    "    L = P @ l\n",
    "    return L, u, p_revert"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 154,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "def reconstruct(L, u, p_revert):\n",
    "    return P_revert @ L @ u # or did I misubderstend what is needed"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 167,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 0.000e+00  0.000e+00  0.000e+00  0.000e+00  0.000e+00  0.000e+00]\n",
      " [ 0.000e+00  0.000e+00  0.000e+00  0.000e+00  0.000e+00  0.000e+00]\n",
      " [ 0.000e+00  0.000e+00  0.000e+00  2.220e-16 -1.110e-16 -1.665e-16]\n",
      " [ 0.000e+00  0.000e+00  2.220e-16 -5.551e-17 -1.665e-16 -1.665e-16]\n",
      " [ 0.000e+00  0.000e+00 -1.110e-16 -1.665e-16  1.665e-16  5.551e-17]\n",
      " [ 0.000e+00  0.000e+00 -1.665e-16 -1.665e-16  5.551e-17  0.000e+00]]\n"
     ]
    }
   ],
   "source": [
    "L, u, p_revert = lu(a)\n",
    "t = reconstruct(L, u, p_revert)\n",
    "print(t - a)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 163,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 0.000e+00  0.000e+00  0.000e+00  0.000e+00  0.000e+00  0.000e+00]\n",
      " [ 0.000e+00  0.000e+00  0.000e+00  0.000e+00  0.000e+00  0.000e+00]\n",
      " [ 0.000e+00  0.000e+00  0.000e+00  2.220e-16 -1.110e-16 -1.665e-16]\n",
      " [ 0.000e+00  0.000e+00  2.220e-16 -5.551e-17 -1.110e-16 -1.665e-16]\n",
      " [ 0.000e+00  0.000e+00 -1.110e-16 -1.665e-16  1.665e-16  5.551e-17]\n",
      " [ 0.000e+00  0.000e+00 -1.665e-16 -1.665e-16  5.551e-17  0.000e+00]]\n"
     ]
    }
   ],
   "source": [
    "L, u, p_revert = lu(a1)\n",
    "t = reconstruct(L, u, p_revert)\n",
    "print(t - a1)"
   ]
  }
 ],
 "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.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}