Gaussian elimination (incomplete)

gauss.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Gauss–Jordan elimination in Python\n",
    "\n",
    "Note that we use sympy for its `Matrix` object and its primitive operations (`row_insert`, `col_delete`, `row(i)`, etc.). We also import some extra classes for type annotations. But we are strictly rolling our own row operations!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from sympy import Matrix, Eq, Number, init_printing\n",
    "init_printing()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We define a couple helper functions to make it easier to work with sympy matrices, which have strange mutation semantics (the `_del` family mutates while the `_insert` family does not).\n",
    "\n",
    "All our functions are strictly functional from the perspective of the callee; that is, none of the functions mutate their arguments.\n",
    "\n",
    "Where appropriate, we use PEP 3107 function annotations as documentation to indicate recommended argument type. Of course, Python does not perform type checking on function calls."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "p = lambda *x: print(*x, sep=\"\\n -> \", end=\"\\n\\n\")\n",
    "\n",
    "def rup(m: Matrix, i: int, new: Matrix) -> Matrix:\n",
    "    '''Update row with different values.'''\n",
    "    m = Matrix(m)\n",
    "    m.row_del(i)\n",
    "    return m.row_insert(i, new)\n",
    "\n",
    "def iter_rows(m: Matrix) -> (list):\n",
    "    rows, _ = m.shape\n",
    "    for i in range(rows):\n",
    "        yield list(m.row(i))\n",
    "\n",
    "def iter_cols(m: Matrix) -> (list):\n",
    "    _, cols = m.shape\n",
    "    for i in range(cols):\n",
    "        yield list(m.col(i))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We define test matrices `consistent` and `inconsistent` so we can play with the elementary row functions as we tweak the algorithm."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "consistent = Matrix(\n",
    "   [[  0,  1,  5, -2 ],\n",
    "    [  1,  0, -3,  8 ],\n",
    "    [  2,  2,  9,  7 ]]\n",
    ")\n",
    "\n",
    "inconsistent = Matrix(\n",
    "   [[  1, -5,  4, -3 ],\n",
    "    [  2, -7,  3, -2 ],\n",
    "    [ -2,  1,  7, -1 ]]\n",
    ")\n",
    "\n",
    "ref = Matrix(\n",
    "   [[  5,  8,  3, -1 ],\n",
    "    [  0,  2,  0,  8 ],\n",
    "    [  0,  0,  1,  1 ]]\n",
    ")\n",
    "\n",
    "rref = Matrix(\n",
    "   [[  1,  0,  0,  2 ],\n",
    "    [  0,  1,  0,  2 ],\n",
    "    [  0,  0,  1,  2 ]]\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Elementary Row Operations\n",
    "\n",
    "In Gauss–Jordan elimination, we have three _elementary row operations_ at our disposal.\n",
    "\n",
    "1. `rswap` swaps two rows `i` and `j`\n",
    "2. `rci` multiplies a row `i` by a constant `c`\n",
    "3. `rcij` adds the product of a row `i` and a constant `c` to a different row `j`\n",
    "\n",
    "The awful names for these functions come from their equivalents in the HP-48G RPL language and represent the order of the function arguments."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def rswap(m: Matrix, i: int, j: int) -> Matrix:\n",
    "    '''Swap two rows.'''\n",
    "\n",
    "    # get frozen versions of current i and j rows\n",
    "    row_i = m.row(i)\n",
    "    row_j = m.row(j)\n",
    "\n",
    "    m = rup(m, i, row_j)\n",
    "    m = rup(m, j, row_i)\n",
    "\n",
    "    return m\n",
    "\n",
    "def rci(m: Matrix, c: Number, i: int) -> Matrix:\n",
    "    '''Multiply a row by a constant.'''\n",
    "    row_i = m.row(i)\n",
    "    return rup(m, i, c * row_i)\n",
    "\n",
    "def rcij(m: Matrix, c: Number, i: int, j: int) -> Matrix:\n",
    "    '''Adds the product of a row i and a constant c to a different row j.'''\n",
    "    row_i = m.row(i)\n",
    "    row_j = m.row(j)\n",
    "    return rup(m, j, c * row_i + row_j)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "m\n",
      " -> Matrix([[0, 1, 5, -2], [1, 0, -3, 8], [2, 2, 9, 7]])\n",
      "\n",
      "rswap(m, 0, 1)\n",
      " -> Matrix([[1, 0, -3, 8], [0, 1, 5, -2], [2, 2, 9, 7]])\n",
      "\n",
      "rci(m, 100, 0)\n",
      " -> Matrix([[0, 100, 500, -200], [1, 0, -3, 8], [2, 2, 9, 7]])\n",
      "\n",
      "rcij(m, 100, 0, 1)\n",
      " -> Matrix([[0, 1, 5, -2], [1, 100, 497, -192], [2, 2, 9, 7]])\n",
      "\n"
     ]
    }
   ],
   "source": [
    "p(\"m\", consistent)\n",
    "p(\"rswap(m, 0, 1)\", rswap(consistent, 0, 1))\n",
    "p(\"rci(m, 100, 0)\", rci(consistent, 100, 0))\n",
    "p(\"rcij(m, 100, 0, 1)\", rcij(consistent, 100, 0, 1))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "Now that elementary row operations are defined, we need to be able to programmatically determine whether a matrix is in row-echelon or reduced-row-echelon form."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def is_zeros(row: list) -> bool:\n",
    "    return all([ Eq(x, 0) for x in row ])\n",
    "\n",
    "def zeros_at_bottom(m: Matrix) -> bool:\n",
    "    '''Once we find a row with all zeros, we expect\n",
    "    subsequent rows to be all zeros iff all zero-rows\n",
    "    are beneath all non-zero rows.'''\n",
    "\n",
    "    found_zero_row = False\n",
    "    for row in iter_rows(m):\n",
    "        if is_zeros(row):\n",
    "            found_zero_row = True\n",
    "        elif found_zero_row:\n",
    "            return False\n",
    "    return True\n",
    "\n",
    "def first_nonzero_element(row: list) -> (int, Number) or (None, None):\n",
    "    for e, x in enumerate(row):\n",
    "        if not Eq(x, 0):\n",
    "            return (e, x)\n",
    "    return (None, None)\n",
    "\n",
    "def is_ref(m: Matrix, reduced=False) -> bool:\n",
    "    '''Is the first nonzero element of each nonzero\n",
    "    row in a column greater than the column of\n",
    "    the first nonzero element of the previous row?\n",
    "    Reduced row-echelon form additionally requires\n",
    "    that every leading coefficient is 1 and the only\n",
    "    nonzero entry in the column.'''\n",
    "\n",
    "    if not zeros_at_bottom(m):\n",
    "        return False\n",
    "\n",
    "    # index of last leading entry\n",
    "    last_e = -1\n",
    "\n",
    "    for i, row in enumerate(iter_rows(m)):\n",
    "\n",
    "        if is_zeros(row):\n",
    "            continue\n",
    "        e, x = first_nonzero_element(row)\n",
    "        if e <= last_e:\n",
    "            return False\n",
    "\n",
    "        # reduced row-echelon form additionally requires\n",
    "        if reduced:\n",
    "\n",
    "            # that the first nonzero element is 1\n",
    "            if not Eq(x, 1):\n",
    "                return False\n",
    "\n",
    "            # that all other elements in the first nonzero\n",
    "            # element's column are zero\n",
    "            col = list(m.col(e))\n",
    "            col.pop(i)\n",
    "            if not is_zeros(col):\n",
    "                return False\n",
    "\n",
    "        last_e = e\n",
    "    return True\n",
    "\n",
    "def is_rref(m: Matrix) -> bool:\n",
    "    return is_ref(m, reduced=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "generic matrix is ref?\n",
      " -> False\n",
      "\n",
      "generic matrix is rref?\n",
      " -> False\n",
      "\n",
      "ref matrix is ref?\n",
      " -> True\n",
      "\n",
      "ref matrix is rref?\n",
      " -> False\n",
      "\n",
      "rref matrix is ref?\n",
      " -> True\n",
      "\n",
      "rref matrix is rref?\n",
      " -> True\n",
      "\n"
     ]
    }
   ],
   "source": [
    "p(\"generic matrix is ref?\", is_ref(consistent))\n",
    "p(\"generic matrix is rref?\", is_rref(consistent))\n",
    "p(\"ref matrix is ref?\", is_ref(ref))\n",
    "p(\"ref matrix is rref?\", is_rref(ref))\n",
    "p(\"rref matrix is ref?\", is_ref(rref))\n",
    "p(\"rref matrix is rref?\", is_rref(rref))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Elimination algorithm helpers"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def to_bottom(m: Matrix, i: int) -> Matrix:\n",
    "    m = Matrix(m)\n",
    "    \n",
    "    to_insert = m.row(i)\n",
    "    m.row_del(i)\n",
    "    \n",
    "    rows, _ = m.shape\n",
    "    return m.row_insert(rows, to_insert)\n",
    "\n",
    "def which_zero_rows(m: Matrix) -> [int]:\n",
    "    return [ i for i, row in iter_rows(m) if is_zeros(row) ]\n",
    "\n",
    "def zeros_to_bottom(m: Matrix) -> (int, Matrix):\n",
    "    for row in which_zero_rows(m):\n",
    "        m = to_bottom(m, row)\n",
    "    return m\n",
    "\n",
    "def list_swap(l: list, i: int, j: int) -> list:\n",
    "    l = list(l)\n",
    "\n",
    "    # get frozen versions of current i and j rows\n",
    "    new_i = l[j]\n",
    "    new_j = l[i]\n",
    "    \n",
    "    l[i].pop()\n",
    "    l.insert(i, new_i)\n",
    "    l[j].pop()\n",
    "    l.insert(j, new_j)\n",
    "\n",
    "    return l\n",
    "\n",
    "def rsort(m: Matrix) -> Matrix:\n",
    "    '''Each row is sorted by the column position of its first nonzero element.'''\n",
    "    rows = [ (first_nonzero_element(row)[0], row) for row in iter_rows(m) ]\n",
    "    rows.sort(key=lambda x: x[0])\n",
    "    rows = [ row[1] for row in rows ]\n",
    "    return Matrix(rows)\n",
    "\n",
    "def rnormalize(m: Matrix) -> Matrix:\n",
    "    '''Makes all first nonzero elements equal one.'''\n",
    "    for i, row in enumerate(iter_rows(m)):\n",
    "        _, x = first_nonzero_element(row)\n",
    "        m = rci(m, 1 / x, i)\n",
    "    return m"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def gauss(m: Matrix) -> Matrix:\n",
    "    if not zeros_at_bottom(m):\n",
    "        m = zeros_to_bottom(m)\n",
    "\n",
    "    # for each first nonzero element in a row, make it the only nonzero element\n",
    "    # in the column.\n",
    "    rows, cols = m.shape\n",
    "    for i in range(rows):\n",
    "        row = m.row(i)\n",
    "\n",
    "        # get first nonzero element in each row\n",
    "        e, pivot = first_nonzero_element(m.row(i))\n",
    "\n",
    "        # we're done if all we have left is zero rows\n",
    "        if e is None:\n",
    "            break\n",
    "\n",
    "        # if the first nonzero element is the last element, there's no solution\n",
    "        if e >= cols - 1:\n",
    "            return None\n",
    "\n",
    "        # zero out all other numbers in the first nonzero element's column\n",
    "        for j in range(rows):\n",
    "\n",
    "            # don't change the first nonzero element's row, also don't change if\n",
    "            # it's already zero\n",
    "            if j == i or m.row(j)[e] == 0:\n",
    "                continue\n",
    "            \n",
    "            # find out what constant c to apply to row i such that\n",
    "            # m_j,e + m_i,e * c = 0\n",
    "            # m_j,e = -m_i,e * c\n",
    "            # m_j,e / -m_i,e = c\n",
    "            c = -1 * m.row(j)[e] / m.row(i)[e]\n",
    "            m = rcij(m, c, i, j)\n",
    "    m = rsort(m)\n",
    "    m = rnormalize(m)\n",
    "    return m"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$$\\left[\\begin{matrix}0 & 1 & 5 & -2\\\\1 & 0 & -3 & 8\\\\2 & 2 & 9 & 7\\end{matrix}\\right]$$"
      ],
      "text/plain": [
       "⎡0  1  5   -2⎤\n",
       "⎢            ⎥\n",
       "⎢1  0  -3  8 ⎥\n",
       "⎢            ⎥\n",
       "⎣2  2  9   7 ⎦"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "consistent"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$$\\left[\\begin{matrix}1 & 0 & 0 & 5\\\\0 & 1 & 0 & 3\\\\0 & 0 & 1 & -1\\end{matrix}\\right]$$"
      ],
      "text/plain": [
       "⎡1  0  0  5 ⎤\n",
       "⎢           ⎥\n",
       "⎢0  1  0  3 ⎥\n",
       "⎢           ⎥\n",
       "⎣0  0  1  -1⎦"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "s = gauss(consistent)\n",
    "s"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$$\\left[\\begin{matrix}1 & -5 & 4 & -3\\\\2 & -7 & 3 & -2\\\\-2 & 1 & 7 & -1\\end{matrix}\\right]$$"
      ],
      "text/plain": [
       "⎡1   -5  4  -3⎤\n",
       "⎢             ⎥\n",
       "⎢2   -7  3  -2⎥\n",
       "⎢             ⎥\n",
       "⎣-2  1   7  -1⎦"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "inconsistent"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "is_rref(gauss(consistent))\n",
      " -> True\n",
      "\n",
      "gauss(inconsistent) is None\n",
      " -> True\n",
      "\n"
     ]
    }
   ],
   "source": [
    "p(\"is_rref(gauss(consistent))\", is_rref(s))\n",
    "p(\"gauss(inconsistent) is None\", gauss(inconsistent) is None)"
   ]
  }
 ],
 "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.5.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}