{
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "# LU and upper echelon form"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 5,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "import numpy as np\n",
        "import scipy.linalg as la"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Part I: The Problem"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "---------------\n",
        "Here's a matrix:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 6,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "A = np.array([[0.,  0., -1., -1.,  0.,  0.,  0., -1.,  0.],\n",
        "         [0.,  0.,  0.,  0.,  0., -1.,  0.,  1.,  0.],\n",
        "         [0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.],\n",
        "         [-1.,  0.,  1.,  0.,  0.,  0., -1.,  0.,  0.],\n",
        "         [1., -1.,  0.,  1., -1.,  1.,  0.,  0., -1.],\n",
        "         [0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  1.],\n",
        "         [0.,  0.,  0.,  0.,  0.,  0.,  1.,  0.,  0.],\n",
        "         [0.,  1.,  0.,  0.,  0.,  0.,  0.,  0.,  0.],\n",
        "         [0.,  0.,  0.,  0.,  1.,  0.,  0.,  0.,  0.]])\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 12,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "P, L, U = la.lu(A)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Check that it is actually a factorization:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 13,
      "metadata": {
        "collapsed": false
      },
      "outputs": [
        {
          "data": {
            "text/plain": [
              "0.0"
            ]
          },
          "execution_count": 13,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "la.norm(A-P.dot(L).dot(U))"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Now look at `U`:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 10,
      "metadata": {
        "collapsed": false
      },
      "outputs": [
        {
          "data": {
            "text/plain": [
              "array([[-1.,  0.,  1.,  0.,  0.,  0., -1.,  0.,  0.],\n",
              "       [ 0., -1.,  1.,  1., -1.,  1., -1.,  0., -1.],\n",
              "       [ 0.,  0., -1., -1.,  0.,  0.,  0., -1.,  0.],\n",
              "       [ 0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.],\n",
              "       [ 0.,  0.,  0.,  0., -1.,  1., -1., -1., -1.],\n",
              "       [ 0.,  0.,  0.,  0.,  0., -1.,  0.,  1.,  0.],\n",
              "       [ 0.,  0.,  0.,  0.,  0.,  0.,  1.,  0.,  0.],\n",
              "       [ 0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  1.],\n",
              "       [ 0.,  0.,  0.,  0.,  0.,  0.,  0.,  0., -1.]])"
            ]
          },
          "execution_count": 10,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "U"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "* What do you notice about the last two rows?\n",
        "* Would this be allowed if `U` were upper echelon?\n",
        "----"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "# Part II: Getting to Echelon Form"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 19,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "def swap_rows(mat, i, j):\n",
        "    temp = mat[i].copy()\n",
        "    mat[i] = mat[j]\n",
        "    mat[j] = temp"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 36,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "def m_echelon(A):\n",
        "    m, n = A.shape\n",
        "    M = np.eye(m)\n",
        "    U = A.copy()\n",
        "\n",
        "    row = 0 \n",
        "    for col in range(min(m, n)):\n",
        "        piv_row = row + np.argmax(np.abs(U[row:, col]))\n",
        "                 \n",
        "        if abs(U[piv_row, col]) == 0:\n",
        "            # column is exactly zero\n",
        "            continue\n",
        "                 \n",
        "        swap_rows(U, row, piv_row)\n",
        "        swap_rows(M, row, piv_row)\n",
        "                 \n",
        "        for el_row in range(row + 1, m):\n",
        "            fac = -U[el_row, col]/U[row, col]\n",
        "            U[el_row] += fac*U[row]\n",
        "            M[el_row] += fac*M[row]\n",
        "\n",
        "        row += 1 \n",
        "\n",
        "    return M, U"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Compute M and U, and check that $MA=U$:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 29,
      "metadata": {
        "collapsed": false
      },
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "0.0\n"
          ]
        }
      ],
      "source": [
        "M, U = m_echelon(A)\n",
        "\n",
        "diff = M.dot(A)-U\n",
        "\n",
        "print(la.norm(diff))"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Let's see if $U$ is actually in echelon form:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 25,
      "metadata": {
        "collapsed": false
      },
      "outputs": [
        {
          "data": {
            "text/plain": [
              "array([[-1.,  0.,  1.,  0.,  0.,  0., -1.,  0.,  0.],\n",
              "       [ 0., -1.,  1.,  1., -1.,  1., -1.,  0., -1.],\n",
              "       [ 0.,  0., -1., -1.,  0.,  0.,  0., -1.,  0.],\n",
              "       [ 0.,  0.,  0.,  0., -1.,  1., -1., -1., -1.],\n",
              "       [ 0.,  0.,  0.,  0.,  0., -1.,  0.,  1.,  0.],\n",
              "       [ 0.,  0.,  0.,  0.,  0.,  0.,  1.,  0.,  0.],\n",
              "       [ 0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  1.],\n",
              "       [ 0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.],\n",
              "       [ 0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.]])"
            ]
          },
          "execution_count": 25,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "U"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "And what does $M$ look like?"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 27,
      "metadata": {
        "collapsed": false
      },
      "outputs": [
        {
          "data": {
            "text/plain": [
              "array([[ 0.,  0.,  0.,  1.,  0.,  0.,  0.,  0.,  0.],\n",
              "       [ 0.,  0.,  0.,  1.,  1.,  0.,  0.,  0.,  0.],\n",
              "       [ 1.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.],\n",
              "       [ 1.,  0.,  0.,  1.,  1.,  0.,  0.,  1.,  0.],\n",
              "       [ 0.,  1.,  0.,  0.,  0.,  0.,  0.,  0.,  0.],\n",
              "       [ 0.,  0.,  0.,  0.,  0.,  0.,  1.,  0.,  0.],\n",
              "       [ 0.,  0.,  0.,  0.,  0.,  1.,  0.,  0.,  0.],\n",
              "       [ 0.,  0.,  1.,  0.,  0.,  0.,  0.,  0.,  0.],\n",
              "       [ 1.,  1.,  0.,  1.,  1.,  1.,  1.,  1.,  1.]])"
            ]
          },
          "execution_count": 27,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "M"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "... not much structure here."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "--------------\n",
        "But we can still have *something* a little like LU:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 34,
      "metadata": {
        "collapsed": false
      },
      "outputs": [
        {
          "data": {
            "text/plain": [
              "array([[ 0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.],\n",
              "       [ 0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.],\n",
              "       [ 0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.],\n",
              "       [ 0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.],\n",
              "       [ 0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.],\n",
              "       [ 0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.],\n",
              "       [ 0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.],\n",
              "       [ 0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.],\n",
              "       [ 0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.]])"
            ]
          },
          "execution_count": 34,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "A - la.inv(M).dot(U)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": []
    }
  ],
  "metadata": {},
  "nbformat": 4,
  "nbformat_minor": 0
}