PY3-solving-linear-systems.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## <center> Solving Linear Systems Using `Numpy` </center>##"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Linear Systems ###\n",
    "\n",
    "An $m\\times n$ [linear system of equations](https://en.wikipedia.org/wiki/System_of_linear_equations) is a collection of linear equations:\n",
    "\n",
    "$$\n",
    "       \\begin{eqnarray*}\n",
    "        a_{11}x_1 + a_{12}x_2 + \\cdots + a_{1n}x_n &=& b_1 \\\\\n",
    "        a_{21}x_1 + a_{22}x_2 + \\cdots + a_{2n}x_n &=& b_2 \\\\\n",
    "        \\vdots \\hspace{0.5in} &  & \\\\\n",
    "        a_{m1}x_1 + a_{m2}x_2 + \\cdots + a_{mn}x_n &=& b_m \\\\\n",
    "       \\end{eqnarray*}\n",
    "$$\n",
    "\n",
    "In matrix notation, a linear system  is $A\\vec{x} = \\vec{b}$, where\n",
    "\n",
    "$$\n",
    "  A = \\begin{bmatrix}\n",
    "                a_{11} & a_{12} & \\cdots & a_{1n}  \\\\\n",
    "                a_{21} & a_{22} & \\cdots & a_{2n} \\\\\n",
    "                \\vdots &        &       &         \\\\\n",
    "                a_{m1} & a_{m2} & \\cdots & a_{mn} \n",
    "      \\end{bmatrix},\\ \\ \\ \\ \n",
    "\\vec{x} = \\left[ \\begin{array}{c}\n",
    "                    x_1 \\\\ x_2 \\\\ \\vdots \\\\ x_n\n",
    "                    \\end{array}\n",
    "             \\right], \\ \\ \\ \\ \n",
    "\\vec{b} = \\begin{bmatrix} b_1 \\\\ b_2 \\\\ b_3 \\end{bmatrix}\n",
    "$$\n",
    "\n",
    "The corresponding  $m\\times (n+1)$ augmented matrix is\n",
    "$$ [\\mathrm{A}|\\vec{b}]\\ = \\left[ \\begin{array}{cccc|c}\n",
    "                            a_{11} & a_{12} & \\cdots & a_{1n} & b_1 \\\\\n",
    "                            a_{21} & a_{22} & \\cdots & a_{2n} & b_2\\\\\n",
    "                            \\vdots &        &                 &    \\\\\n",
    "                            a_{m1} & a_{m2} & \\cdots & a_{mn} & b_m\\\\\n",
    "                           \\end{array}\n",
    "                    \\right]\n",
    "$$   \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Solving Linear Systems with `numpy.linalg.solve` or `numpy.linalg.inv` ###\n",
    "\n",
    "\n",
    "If a coefficient matrix $A$ is square and nonsingular (i.e., $\\mathrm{det}\\,A\\neq 0$), we can use the function [numpy.linalg.solve](https://docs.scipy.org/doc/numpy/reference/generated/numpy.linalg.solve.html#numpy.linalg.solve) or [numpy.linalg.inv](https://docs.scipy.org/doc/numpy/reference/generated/numpy.linalg.inv.html#numpy.linalg.inv). (For more  `numpy.linalg` functions [see](https://docs.scipy.org/doc/numpy/reference/routines.linalg.html#module-numpy.linalg) this.)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The functions return a solution of the system $A\\vec{x}=\\vec{b}$. For example,"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "A = np.array([[1,1],[1,-1]])\n",
    "print(A)\n",
    "print(\"\")\n",
    "print(\"det(A) = \", np.linalg.det(A))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Since $\\mathrm{det}\\,A\\neq 0$, the linear system has a unique solution and we may use `numpy.linalg.solve` to find the solution. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "b1 = np.array([2,0])\n",
    "print(b1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x1 = np.linalg.solve(A,b1)\n",
    "print(x1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that the output  is returned as a 1D NumPy array when the vector  (the right hand side) is entered as a 1D NumPy array. If we input  as a 2D NumPy array, then the output is a 2D NumPy array. For example:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "A = np.array([[1,1],[1,-1]])\n",
    "b2 = np.array([2,0]).reshape(2,1)\n",
    "x2 = np.linalg.solve(A,b2)\n",
    "print(x2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If we want to use `numpy.linalg.inv`,"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "A_inv = np.linalg.inv(A)\n",
    "print(A_inv)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And multiply $A^{-1}\\vec{b}$ to find $\\vec{x}$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = A_inv @ b1\n",
    "print(x)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We get the same result."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Inverse or Solve ?** It's not a good idea to use the inverse $A^{-1}$ to solve $A\\vec{x} = \\vec{b}$ if $A$ is large. It's too computationally expensive. Let's create a large random matrix $A$ and a vector $\\vec{b}$, and compute the solution $\\vec{x}$ in two ways: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "N = 1000\n",
    "B = np.random.rand(N,N)\n",
    "b3 = np.random.rand(N,1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Check the first entries of $B$ and $\\vec{b}_3$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(B[:3,:3])\n",
    "print(\" \")\n",
    "print(b2[:3,:])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, we compare the speed of `numpy.linalg.solve` with `numpy.linalg.inv`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%timeit\n",
    "x = np.linalg.solve(B,b3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%timeit\n",
    "x = np.linalg.inv(B) @ b3"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Solving with `scipy.linalg.solve` is about three times faster. (The exact speed may be different depending on your machine.)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "** Remark:** If $A$ is not square or nonsingular, then there are other `numpy.linalg` functions to be employed. For example, `np.linalg.lstsq`."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Gaussian Elimination###\n",
    "\n",
    "The general procedure to solve a linear system of equation is called [Gaussian elimination](https://en.wikipedia.org/wiki/Gaussian_elimination). This method is for any dimensions of matrices, including non-square matrices. The idea is to perform [elementary row operations](https://en.wikipedia.org/wiki/Elementary_matrix#Elementary_row_operations) to reduce the system to its row echelon form and then solve. Elementary row operations are\n",
    "  1. Interchange rows $i$ and $j$ ($R_i\\leftrightarrow R_j$).\n",
    "  2. Multiply row $i$ by a *nonzero* scalar $k$ ($kR_i\\rightarrow R_i$).\n",
    "  3. Add a *nonzero* scalar $k$ times row $i$ to row $j$ ($kR_i+R_j\\rightarrow R_j$), where $i\\neq j$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Elementary Matrices###  \n",
    "\n",
    "Each of the elementary row operations is the result of matrix multiplication by an elementary matrix (on the left).\n",
    "\n",
    "  [1] For $R_i\\leftrightarrow R_j$ in an $m\\times n$ matrix $A$, we multiply $A$ by an $m\\times m$ matrix $E$ where $E$ is equal to the identity matrix $I_m$ except $E_{ii}=E_{jj}=0,$ and $E_{ij}=E_{ji}=1$. (Equivalently, we can interchange $i$-th row and $j$-th row of the identity matrix $I$ to get $E$.) For example, if $A$ is $3\\times 4$ and we want to switch $R_1$ and $R_2$, then\n",
    "  \n",
    "  $$\n",
    "   \\begin{bmatrix} 1 & 0 & 0 \\\\ 0 & 0 & 1\\\\ 0 & 1 & 0 \\end{bmatrix}\n",
    "  $$\n",
    "\n",
    "Let's verify it:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "A1 = np.array([[1,1,2, -1],[0, -1,3,1],[1,0,5,2]])\n",
    "print(\"A is \")\n",
    "print(A)\n",
    "print(\" \")\n",
    "\n",
    "E1 = np.array([[1,0,0],[0,0,1],[0,1,0]])\n",
    "print(\"E1 is\")\n",
    "print(E1)\n",
    "print(\" \")\n",
    "\n",
    "A1_new = E1 @ A1\n",
    "print(\"A1_new is \")\n",
    "print(A1_new)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "  [2] For $kR_i\\rightarrow R_i$ in an $m\\times n$ matrix $A$, we multiply $A$ by the matrix $E$ where $E$  is equal to the identity matrix $I_m$ except $E_{ii}=k$. For example, if $A$ is 4 by 2 and we want to multiply row 2 by -2 then\n",
    "  \n",
    "   $$\n",
    "   \\begin{bmatrix} 1 & 0 & 0 & 0\\\\ 0 & -2 & 0 &0\\\\ 0 & 0 & 1 &0 \\\\\n",
    "       0&  0 &0 &1 \\end{bmatrix}\n",
    "  $$\n",
    "\n",
    "Let's verify it:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "A2 = np.array([[2, -1],[0, -1],[3,1],[1,0]])\n",
    "print(\"A2 is \")\n",
    "print(A2)\n",
    "print(\" \")\n",
    "\n",
    "E2 = np.array([[1,0,0,0],[0,-2,0,0],[0,0,1,0],[0,0,0,1]])\n",
    "print(\"E2 is\")\n",
    "print(E2)\n",
    "print(\" \")\n",
    "\n",
    "A2_new = E2 @ A2\n",
    "print(\"A_new2 is \")\n",
    "print(A2_new)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "  [3]  For $kR_i+R_j\\rightarrow R_j$ in an $m\\times n$ matrix $A$, we multiply $A$ by the matrix $E$ where $E$  is equal to the identity matrix $I_m$ except $E_{ji}=k$. For example, if $A$ is 3 by 3 and we want to add 3 times $R_3$ to row $R_1$ (i.e., $3R_3+R_1\\rightarrow R_1$), then\n",
    "  \n",
    "   $$\n",
    "   \\begin{bmatrix} 1 & 0 & 3 \\\\ 0 & 1 & 0\\\\ 0 & 0 & 1 \\end{bmatrix}\n",
    "  $$\n",
    "\n",
    "Let's verify it:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "A3 = np.array([[1,1,2],[-1,3,1],[0,5,2]])\n",
    "print(\"A3 is \")\n",
    "print(A3)\n",
    "print(\" \")\n",
    "\n",
    "E3 = np.array([[1,0,3],[0,1,0],[0,0,1]])\n",
    "print(\"E3 is\")\n",
    "print(E3)\n",
    "print(\" \")\n",
    "\n",
    "A3_new = E3 @ A3\n",
    "print(\"A3_new is \")\n",
    "print(A3_new)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Implementation ###\n",
    "\n",
    "Let's write function to implement the elementary row operations. \n",
    "\n",
    "  [1] First of all, let's write a function called `switch_rows` which takes 3 input parameters $A, i$ and $j$ and returns the matrix that results from interchanging rows $i$ and $j$ in the matrix $A$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def switch_rows(A,i,j):\n",
    "    \"Switch rows 1 and j in matrix A.\"\n",
    "    m = A.shape[0]\n",
    "    E = np.eye(m)\n",
    "    E[i,i] = 0\n",
    "    E[j,j] = 0\n",
    "    E[i,j] = 1\n",
    "    E[j,i] = 1\n",
    "    return E @ A"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "A4 = np.array([[1,1,1],[1,-1,0]])\n",
    "print(A4)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "switch_rows(A4,0,1)  # Note that the index in python starts from 0 but not 1."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "  [2] Second, let's write a function called `scale_row` which takes 3 input parameters $A, k,$ and $i$ and returns the matrix that results from  $k$ times row $i$ in the matrix $A$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def scale_row(A,k,i):\n",
    "    \"Multiply row i by k in matrix A\"\n",
    "    m = A.shape[0]\n",
    "    E = np.eye(m)\n",
    "    E[i,i] = k\n",
    "    return E @ A"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "A5 = np.array([[3,1],[-2,7]])\n",
    "print(A5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "scale_row(A5,3,1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "  [3] Finally, let's write a function called `add_rows` which takes input parameters $A, k, i$ and $j$ and returns the NumPy array resulting from adding $k$ times row $i$  to row $j$  in the matrix $A$. Note that we expect $i\\neq j$, so have to put a warning message."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def add_row(A,k,i,j):\n",
    "    \"Add k times row i to row j in matrix A.\"\n",
    "    m = A.shape[0]\n",
    "    E = np.eye(m)\n",
    "    if i==j:\n",
    "        print(\"Warning: i=j. i and j must be different.\")\n",
    "        return None\n",
    "    else:\n",
    "        E[j,i] = k\n",
    "        return E @ A"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "A6 = np.array([[1,1],[3,2]])\n",
    "print(A6)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "add_row(A6,2,0,0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercises ##\n",
    "\n",
    "  1. You want to make a certain kind of tropical punch, using bananas, oranges, and papayas. Suppose you don't know how many of each to put in the punch, but you know that there are seven pieces of fruit in the mix, and there are twice as many oranges as bananas. You also know that the seven pieces of fruit cost 15.25 Dhs, where bananas cost 1.5 Dhs each, oranges cost 2.25 Dhs, and papayas cost 3.75 Dhs each. \n",
    "  \n",
    "    (a) Construct a system of equation describing the problem.    \n",
    "    (b) Solve the system using `numpy.linalg.solve`.   \n",
    "    (c) How many bananas, oranges and papayas would you need to make this punch? \n",
    "    \n",
    "  2. Consider the *tridiagonal* linear system\n",
    " \n",
    " $$ \\begin{eqnarray*}\n",
    "       3x_1 - x_2 &=&  2\\\\\n",
    "       -x_1+3x_2-x_3 & = & 1\\\\\n",
    "           &\\vdots & \\\\\n",
    "       -x_{n-2}+3x_{n-1}-x_{n} &=& 1\\\\\n",
    "       -x_{n-1}+3x_n &=& 2\n",
    "      \\end{eqnarray*}\n",
    "  $$\n",
    "  \n",
    "    (a) Compare the speed of `numpy.linalg.solve` with `numpy.linalg.inv` when solving the linear system for $n=10$.    \n",
    "    (b) Compare the speed of `numpy.linalg.solve` with `numpy.linalg.inv` when solving the linear system for $n=100$. \n",
    "    \n",
    "  3. Use the functions defined in this notebook (`switch_rows`, `scale_row`, and `add_row` to find the general solution of the system\n",
    "   $$ \\begin{eqnarray*}\n",
    "       3x_2 - 6x_3 + 6x_4 +4x_5 &=& -5\\\\\n",
    "       3x_1 - 7x_2 + 8x_3 - 5x_4 + 8x_5 &=& 9\\\\\n",
    "       3x_1 - 9x_2 + 12x_3 - 9x_4 +6x_5 &=& 15\n",
    "      \\end{eqnarray*}\n",
    "  $$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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
}