linear equations (+ HW).ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Solving linear equations"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# I. $LU$ factorization of a square matrix"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let us consider a decomposition of a square $n \\times n$ matrix A as follows:\n",
    "$$A = L \\cdot U, \\; \\mbox{where} \\; A = \\begin{pmatrix} \n",
    "                                            a_{11} & a_{12} & a_{13} & \\ldots & a_{1n} \\\\\n",
    "                                            a_{21} & a_{22} & a_{23} & \\ldots & a_{2n} \\\\\n",
    "                                            a_{31} & a_{32} & a_{33} & \\ldots & a_{3n} \\\\\n",
    "                                            \\vdots & \\vdots & \\vdots & \\ddots & \\vdots \\\\\n",
    "                                            a_{n1} & a_{n2} & a_{n3} & \\ldots & a_{nn} \\\\\n",
    "                                        \\end{pmatrix}\n",
    "                               , \\; L = \\begin{pmatrix} \n",
    "                                            1 & 0 & 0 & \\ldots & 0 \\\\\n",
    "                                            * & 1 & 0 & \\ldots & 0 \\\\\n",
    "                                            * & * & 1 & \\ldots & 0 \\\\\n",
    "                                            \\vdots & \\vdots & \\vdots & \\ddots & \\vdots \\\\\n",
    "                                            * & * & * & \\ldots & 1 \\\\\n",
    "                                        \\end{pmatrix}\n",
    "                               , \\; U = \\begin{pmatrix} \n",
    "                                            a_{11} & * & * & \\ldots & * \\\\\n",
    "                                            0 & a_{22} & * & \\ldots & * \\\\\n",
    "                                            0 & 0 & a_{33} & \\ldots & * \\\\\n",
    "                                            \\vdots & \\vdots & \\vdots & \\ddots & \\vdots \\\\\n",
    "                                            0 & 0 & 0 & \\ldots & a_{nn} \\\\\n",
    "                                        \\end{pmatrix} .$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's start with Gaussian elimination. When we are working with the first column, we combine the first row and the second row multiplied by coefficient $$\\gamma_{21} = \\cfrac{a_{21}}{a{11}};$$ then the first row and the third row multiplied by coefficient $$\\gamma_{31} = \\cfrac{a_{31}}{a{11}},$$ and so on."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Hereby, to eliminate all elements below $a_{11}$ we need to multiply matrix A by matrix $$\\Lambda_1 = \\begin{pmatrix} \n",
    "                                                        1 & 0 & 0 & \\ldots & 0 \\\\\n",
    "                                                        -\\gamma_{21} & 1 & 0 & \\ldots & 0 \\\\\n",
    "                                                        -\\gamma_{31} & 0 & 1 & \\ldots & 0 \\\\\n",
    "                                                        \\vdots & \\vdots & \\vdots & \\ddots & \\vdots \\\\\n",
    "                                                        -\\gamma_{n1} & 0 & 0 & \\ldots & 1 \\\\\n",
    "                                                    \\end{pmatrix}$$\n",
    "(prove it by trying to act with $\\Lambda_1$ on the first column of matrix A)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Likewise, we can construct matrix $\\Lambda_2$ as $$\\Lambda_2 = \\begin{pmatrix} \n",
    "                                                                1 & 0 & 0 & \\ldots & 0 \\\\\n",
    "                                                                0 & 1 & 0 & \\ldots & 0 \\\\\n",
    "                                                                0 & -\\gamma_{32} & 1 & \\ldots & 0 \\\\\n",
    "                                                                \\vdots & \\vdots & \\vdots & \\ddots & \\vdots \\\\\n",
    "                                                                0 & -\\gamma_{n2} & 0 & \\ldots & 1 \\\\\n",
    "                                                            \\end{pmatrix}.$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, we will get the upper triangular matrix $$U = \\Lambda_n \\cdot \\Lambda_{n-1} \\cdot \\ldots \\cdot \\Lambda_2 \\cdot \\Lambda_1 A. $$\n",
    "Hence the lower triangular matrix $L = \\Lambda_1^{-1} \\cdot \\Lambda_2^{-1} \\cdot \\ldots \\cdot \\Lambda_{n-1}^{-1} \\cdot \\Lambda_n^{-1}.$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "One can show that, for example, $$\\Lambda_1^{-1} = \\begin{pmatrix} \n",
    "                                                        1 & 0 & 0 & \\ldots & 0 \\\\\n",
    "                                                        \\gamma_{21} & 1 & 0 & \\ldots & 0 \\\\\n",
    "                                                        \\gamma_{31} & 0 & 1 & \\ldots & 0 \\\\\n",
    "                                                        \\vdots & \\vdots & \\vdots & \\ddots & \\vdots \\\\\n",
    "                                                        \\gamma_{n1} & 0 & 0 & \\ldots & 1 \\\\\n",
    "                                                    \\end{pmatrix}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that we're using the `numpy` arrays to represent matrices [do **not** use `np.matrix`]."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "def diy_lu_ext(a):\n",
    "    \"\"\"\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",
    "    #Initializing the factors\n",
    "    u = a.copy()\n",
    "    L = np.eye(N)\n",
    "    \n",
    "    for j in range(N-1):\n",
    "        lam = np.eye(N)\n",
    "        \n",
    "        #Creating the vector of gammas\n",
    "        gamma = np.zeros(N-j-1)\n",
    "        for i in range(N-j-1):\n",
    "            gamma[i] = u[j+1+i, j]/u[j, j]\n",
    "        \n",
    "        #Creating matrix \\Lambda_i\n",
    "        for i in range(N-j-1):\n",
    "            lam[j+1+i, j] = -gamma[i]\n",
    "        \n",
    "        #Acting with \\Lambda_i on A to get U\n",
    "        u_new = np.zeros((N, N))\n",
    "        for ind_i in range(N):\n",
    "            for ind_j in range(N):\n",
    "                for ind_k in range(N):\n",
    "                    u_new[ind_i, ind_j] += lam[ind_i, ind_k] * u[ind_k, ind_j]\n",
    "        u = u_new.copy()\n",
    "        \n",
    "        #Creating matrix \\Lambda_i^{-1}\n",
    "        for i in range(N-j-1):\n",
    "            lam[j+1+i, j] = gamma[i]\n",
    "            \n",
    "        #Multiplying L and \\Lambda_i^{-1} o get new L\n",
    "        L_new = np.zeros((N, N))\n",
    "        for ind_i in range(N):\n",
    "            for ind_j in range(N):\n",
    "                for ind_k in range(N):\n",
    "                    L_new[ind_i, ind_j] += L[ind_i, ind_k] * lam[ind_k, ind_j]\n",
    "        L = L_new.copy()\n",
    "        \n",
    "    return L, u"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "6"
      ]
     },
     "execution_count": 43,
     "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": 44,
   "metadata": {
    "scrolled": false
   },
   "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 -1.110e-16  1.110e-16  1.110e-16 -5.551e-17]\n",
      " [ 0.000e+00  0.000e+00  3.331e-16 -2.220e-16 -5.551e-17  0.000e+00]\n",
      " [ 0.000e+00  0.000e+00  0.000e+00 -1.110e-16 -1.665e-16  0.000e+00]\n",
      " [ 0.000e+00  0.000e+00 -1.110e-16 -2.776e-16  1.110e-16  1.110e-16]]\n"
     ]
    }
   ],
   "source": [
    "L, u = diy_lu_ext(a)\n",
    "print(L@u - a)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "LU can be programmed in a more simple way by using the perks of `numpy`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "def diy_lu(a):\n",
    "    \"\"\"\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",
    "        \n",
    "        #Creating the vector of gammas\n",
    "        gamma = u[j+1:, j] / u[j, j]\n",
    "        \n",
    "        #Creating matrix \\Lambda_i\n",
    "        lam[j+1:, j] = -gamma\n",
    "        \n",
    "        #Acting with \\Lambda_i on A to get U\n",
    "        u = lam @ u\n",
    "        \n",
    "        #Creating matrix \\Lambda_i^{-1}\n",
    "        lam[j+1:, j] = gamma\n",
    "            \n",
    "        #Multiplying L and \\Lambda_i^{-1} o get new L\n",
    "        L = L @ lam\n",
    "    return L, u"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Tweak the printing of floating-point numbers, for clarity\n",
    "np.set_printoptions(precision=3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "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 -1.110e-16  1.110e-16  1.110e-16 -5.551e-17]\n",
      " [ 0.000e+00  0.000e+00  3.331e-16 -2.220e-16 -5.551e-17  0.000e+00]\n",
      " [ 0.000e+00  0.000e+00  0.000e+00 -1.110e-16 -1.665e-16  0.000e+00]\n",
      " [ 0.000e+00  0.000e+00 -1.110e-16 -2.776e-16  1.110e-16  1.110e-16]]\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": 48,
   "metadata": {},
   "outputs": [],
   "source": [
    "a1 = a.copy()\n",
    "a1[1, 1] = 3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Resulting matrix still has full rank, but the naive LU routine breaks down."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "6"
      ]
     },
     "execution_count": 49,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.linalg.matrix_rank(a1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "metadata": {
    "scrolled": false
   },
   "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": [
      "/usr/lib/python3/dist-packages/ipykernel_launcher.py:18: RuntimeWarning: divide by zero encountered in true_divide\n",
      "/usr/lib/python3/dist-packages/ipykernel_launcher.py:24: RuntimeWarning: invalid value encountered in matmul\n",
      "/usr/lib/python3/dist-packages/ipykernel_launcher.py:30: RuntimeWarning: invalid value encountered in matmul\n",
      "/usr/lib/python3/dist-packages/ipykernel_launcher.py:18: RuntimeWarning: invalid value encountered in true_divide\n"
     ]
    }
   ],
   "source": [
    "l, u = diy_lu(a1)\n",
    "\n",
    "print(l, u)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": []
  },
  {
   "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`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Non-zero of A-matrix principle minors: True\n",
      "Non-zero of A1-matrix principle minors: False\n"
     ]
    }
   ],
   "source": [
    "def check_principle_minors(matrix):\n",
    "    return all([np.linalg.det(matrix[:i, :i]) for i in range(1, matrix.shape[0])])\n",
    "\n",
    "\n",
    "print(f\"Non-zero of A-matrix principle minors: {check_principle_minors(a)}\")\n",
    "print(f\"Non-zero of A1-matrix principle minors: {check_principle_minors(a1)}\")\n"
   ]
  },
  {
   "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",
    "Implement a function to reconstruct the original matrix from a decompositon. Test your routines on the matrices `a` and `a1`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 52,
   "metadata": {},
   "outputs": [],
   "source": [
    "def diy_lu_pivot(a):\n",
    "    \n",
    "    N = a.shape[0]\n",
    "    \n",
    "    P = np.arange(N)\n",
    "    \n",
    "    u = a.copy()\n",
    "    L = np.eye(N)\n",
    "    for j in range(N - 1):\n",
    "        lam = np.eye(N)\n",
    "        max_row_index = np.argmax(abs(u[j:, j])) + j\n",
    "        P[[j, max_row_index]] = P[[max_row_index, j]]\n",
    "        u[[j, max_row_index]] = u[[max_row_index, j]]\n",
    "        gamma = u[j+1:, j] / u[j, j]\n",
    "        lam[j+1:, j] = -gamma\n",
    "        u = lam @ u\n",
    "        lam[j+1:, j] = gamma\n",
    "        L = L @ lam\n",
    "    \n",
    "    return L, u, P\n",
    "\n",
    "L1, u1, P1 = diy_lu_pivot(a1)\n",
    "L, u, P = diy_lu_pivot(a)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 3.     3.     3.     3.     3.     3.   ]\n",
      " [ 3.     0.75   0.429  0.3    0.231  0.188]\n",
      " [ 3.     1.364 -0.506 -0.892 -1.132 -1.295]\n",
      " [ 3.     1.071  0.652  0.424  0.292  0.206]\n",
      " [ 3.     0.882  0.517  0.366  0.284  0.232]\n",
      " [ 3.     3.     2.752  2.661  2.605  2.567]] \n",
      "\n",
      "[[3.    3.    3.    3.    3.    3.   ]\n",
      " [3.    0.75  0.429 0.3   0.231 0.188]\n",
      " [3.    1.364 0.779 0.458 0.253 0.111]\n",
      " [3.    1.071 0.652 0.473 0.375 0.313]\n",
      " [3.    0.882 0.517 0.366 0.283 0.23 ]\n",
      " [3.    1.875 1.467 1.262 1.138 1.055]]\n"
     ]
    }
   ],
   "source": [
    "def original_matrix_reconstruction(l, u):\n",
    "    return l @ u\n",
    "\n",
    "print(original_matrix_reconstruction(L1, u1), \"\\n\")\n",
    "print(original_matrix_reconstruction(L, u))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Sum all elements in matrix `L1` and `u1` separately (for Google Form)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "15.166249240140564\n",
      "-3.456659008328081\n"
     ]
    }
   ],
   "source": [
    "print(np.sum(L1))\n",
    "print(np.sum(u1))"
   ]
  }
 ],
 "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.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}