{ "metadata": { "name": "" }, "nbformat": 3, "nbformat_minor": 0, "worksheets": [ { "cells": [ { "cell_type": "heading", "level": 1, "metadata": {}, "source": [ "Gaussian Elimination - I" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We need the following functions which we wrote in the last class." ] }, { "cell_type": "code", "collapsed": false, "input": [ "def read_matrix(filename) :\n", " fl = open(filename, 'r')\n", " matrix = []\n", " for line in fl :\n", " row = []\n", " words = line.split()\n", " for word in words :\n", " row.append(float(word))\n", " matrix.append(row)\n", " return matrix" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 1 }, { "cell_type": "code", "collapsed": false, "input": [ "A = read_matrix(\"files/A.txt\")\n", "print A" ], "language": "python", "metadata": {}, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ "[[1.0, 2.0, 3.0], [2.0, 1.0, 2.0], [3.0, 2.0, 1.0]]\n" ] } ], "prompt_number": 2 }, { "cell_type": "code", "collapsed": false, "input": [ "def print_matrix(matrix, myformat) :\n", " for row in matrix :\n", " for no in row :\n", " print myformat % no,\n", " print" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 3 }, { "cell_type": "code", "collapsed": false, "input": [ "print_matrix(A, \"%10.0f\")" ], "language": "python", "metadata": {}, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ " 1 2 3\n", " 2 1 2\n", " 3 2 1\n" ] } ], "prompt_number": 4 }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Row operation 1 : Interchanging two rows\n", "As one can guess, this will be a function which takes a matrix as an input, and row numbers of two rows; like `elem1(A, 0, 2)`. However, before we do that, it might be useful to define a function called `copy_row(source)`." ] }, { "cell_type": "code", "collapsed": false, "input": [ "def copy_row(src) :\n", " dest = []\n", " for entry in src :\n", " dest.append(entry)\n", " return dest\n", "\n", "def elem1(matrix, row1, row2) :\n", " if row1 < 0 or row2 < 0 or row1 >= len(matrix) or row2 >= len(matrix) :\n", " print \"Row out of range.\"\n", " else :\n", " temp_row = copy_row(matrix[row1])\n", " matrix[row1] = copy_row(matrix[row2])\n", " matrix[row2] = copy_row(temp_row)\n", " return matrix" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 5 }, { "cell_type": "markdown", "metadata": {}, "source": [ "###Row operation 2 : multiplying a row by a non-zero scalar\n", "The functionw will take a matrix, and row number and a scalar as input. We'll also check if the scalar is non-zero." ] }, { "cell_type": "code", "collapsed": false, "input": [ "def elem2(matrix, row, scalar) :\n", " if row < 0 or row >= len(matrix) :\n", " print \"row %d out of range.\" % row\n", " elif scalar == 0 :\n", " print \"scalar has to be non-zero.\"\n", " else :\n", " for i in range(len(matrix[row])) :\n", " matrix[row][i] *= scalar\n", " return matrix" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 6 }, { "cell_type": "code", "collapsed": false, "input": [ "print_matrix(A, \"%2.0f\")\n", "print \"--------------------------\"\n", "print_matrix(elem2(A, 0, 4), \"%2.0f\")\n", "print \"-------------------------\"\n", "print_matrix(A, \"%2.0f\")" ], "language": "python", "metadata": {}, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ " 1 2 3\n", " 2 1 2\n", " 3 2 1\n", "--------------------------\n", " 4 8 12\n", " 2 1 2\n", " 3 2 1\n", "-------------------------\n", " 4 8 12\n", " 2 1 2\n", " 3 2 1\n" ] } ], "prompt_number": 7 }, { "cell_type": "markdown", "metadata": {}, "source": [ "####row operation 3 : replacing a row by the sum of that row and a multiple of another row.\n", "Here we need 4 inputs to the function:\n", "+ The matrix\n", "+ The row to be changed\n", "+ The row whose multiple will be added to the first row,\n", "+ a scalar : the multiplication factor.\n", "\n", "Again for this we use some helper functions which are readily available for vectors." ] }, { "cell_type": "code", "collapsed": false, "input": [ "def add_vects(lst1, lst2) :\n", " return [a + b for (a, b) in zip(lst1, lst2)]\n", "\n", "def scalar_mult(a, lst) :\n", " return [a * no for no in lst]\n", "\n", "\n", "def elem3(matrix, row_2b_changed, row_used_4_change, scalar) :\n", " temp_row = copy_row(scalar_mult(scalar, matrix[row_used_4_change]))\n", " matrix[row_2b_changed] = add_vects(matrix[row_2b_changed], temp_row)\n", " return matrix" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 8 }, { "cell_type": "code", "collapsed": false, "input": [ "print_matrix(A, \"%2.0f\")\n", "print \"-\"*20\n", "print_matrix(elem3(A, 0, 1, -2), \"%2.0f\")\n", "print \"-\"*21\n", "print_matrix(A, \"%2.0f\")" ], "language": "python", "metadata": {}, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ " 4 8 12\n", " 2 1 2\n", " 3 2 1\n", "--------------------\n", " 0 6 8\n", " 2 1 2\n", " 3 2 1\n", "---------------------\n", " 0 6 8\n", " 2 1 2\n", " 3 2 1\n" ] } ], "prompt_number": 9 }, { "cell_type": "heading", "level": 2, "metadata": {}, "source": [ "Sweeping a column" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "If we recall, sweeping a column requires a matrix and an entry (position in the matrix) which we call pivot. We shall define the pivot as the 2-tuple (row_no, col_no)." ] }, { "cell_type": "code", "collapsed": false, "input": [ "def sweep(matrix, pivot) :\n", " (r, c) = pivot\n", " if r < 0 or r >= len(matrix) :\n", " print \"row out of range\"\n", " elif c < 0 or c >= len(matrix[r]) :\n", " print \"column out of range\"\n", " elif matrix[r][c] == 0:\n", " # pivot cannot be zero.\n", " print(\"pivot cannot be zero.\")\n", " else :\n", " # print_matrix(matrix, \"%5.2f\")\n", " pivot = matrix[r][c]\n", " \n", " # Step 1 : multiply the row containing the pivot by 1/pivot to make the pivot 1.\n", " matrix = elem2(matrix, r, 1.0/pivot)\n", " # print_matrix(matrix, \"%5.2f\")\n", " \n", " # Step 2 : for row != that of pivot, subtract matrix[row][c] times the row containing pivot\n", " for i in range(len(matrix)) :\n", " if i != r :\n", " matrix = elem3(matrix, i, r, -matrix[i][c])\n", " # print_matrix(matrix, \"%5.2f\")\n", " return matrix" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 10 }, { "cell_type": "heading", "level": 2, "metadata": {}, "source": [ "Gaussian elimination" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Gaussian elemination is an algorithm to reduce a matrix to its reduced echelon form. Here *reduction* means performing row operations till the final matrix satisfies the definition of a reduced echelon form. For your convenience let us recall the definition.\n", "\n", "An $m \\times n$ matrix is said to be in echelon form, if it has $r$, $0 \\leq r \\leq m$ non-zero rows and\n", "1. All the non-zero rows are on the top.\n", "1. For $1 \\leq i \\leq r$, $p_i$ denotes the column containing the *first non-zero entry* of the $i$-th row, then\n", " $p_1 < p_2 < \\dotsb < p_r$.\n", "1. $a_{i p_i} = 1$.\n", "The matrix is said to be in reduced echelon form if in addition to being in the echelon form, the $p_i$'th columns have all but one zeroes. The non-zero entry is $a_{i p_i}$." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "####Example :\n", "The matrix\n", "\n", "0 1 4 5\n", "\n", "1 4 3 3\n", "\n", "0 0 0 5\n", "\n", "is not even in echelon form as $p_1 = 2$ and $p_2 = 1$ and hence $p_1 < p_2$ is not satisfied. Note, for the sake of giving an example, here $p_3 = 4$.\n", "\n", "2 1 4 5\n", "\n", "1 4 3 3\n", "\n", "0 0 0 5\n", "\n", "is also not in echelon form because of the same reason.\n", "\n", "1 4 5 3\n", "\n", "0 3 8 2\n", "\n", "0 0 5 6\n", "\n", "does satisy the conditions on $p_i$s but $a_{2 p_2} \\neq 1$ and hence is not in echelon form.\n", "\n", "1 4 0 3\n", "\n", "0 1 0 2\n", "\n", "0 0 1 6\n", "\n", "is in echelon form but not in reduced echelon form as $p_2$-th column has extra non-zero entries (namely 4).\n", "\n", "1 0 0 3\n", "\n", "0 1 0 2\n", "\n", "0 0 1 6\n", "\n", "is in reduced echelon form.\n", "\n", "\n", "#### The algorithm to reduce to echelon form.\n", "1. First find the *first* non-zero column. Suppose the the column number is $p_1$ and the first non-zero entry is in the $i$-th row.\n", "1. Switch the 1st row with the $i$-th row.\n", "1. Sweep the $p_1$-th column\n", "1. Repeat this for the smaller matrix spanning row 2 - last row and column $p_1 + 1$ to the last column." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### First let us write a function to find the first column with a non-zero entry.\n", "It will return a tuple (row, column) for the non-zero entry. If the matrix is zero. It will return (-1, -1)." ] }, { "cell_type": "code", "collapsed": false, "input": [ "def first_non_zero(matrix) :\n", " found_entry = False\n", " column_scanning = 0\n", " if matrix == [] or matrix == [[]] :\n", " print \"The matrix is empty.\"\n", " row = -1\n", " col = -1\n", " else :\n", " while found_entry == False and column_scanning < len(matrix[0]) :\n", " # len(matrix[0]) = number of columns of the matrix.\n", " row_scanning = 0\n", " while found_entry == False and row_scanning < len(matrix) :\n", " #print \"r, c : \", row_scanning, column_scanning\n", " if matrix[row_scanning][column_scanning] != 0 :\n", " found_entry = True\n", " row = row_scanning\n", " col = column_scanning\n", " else :\n", " row_scanning += 1\n", " column_scanning += 1\n", " if found_entry == False :\n", " row = -1\n", " col = -1\n", " print \"The matrix is the zero matrix.\"\n", " return (row, col)" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 11 }, { "cell_type": "code", "collapsed": false, "input": [ "# To test it out :\n", "print first_non_zero([[0, 0, 0, 0], [0, 2, 1, 0], [0, 0, 5, 1], [0, 5, 5, 6]])\n", "print first_non_zero([[0, 0, 0, 0], [0, 0, 0, 0]])\n", "print first_non_zero([[], []])" ], "language": "python", "metadata": {}, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ "(1, 1)\n", "The matrix is the zero matrix.\n", "(-1, -1)\n", "The matrix is the zero matrix.\n", "(-1, -1)\n" ] } ], "prompt_number": 12 }, { "cell_type": "heading", "level": 2, "metadata": {}, "source": [ "Digression : recursion" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In python a function can call itself. It has some uses.\n", "\n", "In math some things are defined recursively.\n", "\n", "#### Example \n", "*Fibonacci sequences* is defined to be sequence $a_1, a_2, \\dotsc$ where\n", "\n", "+ $a_1 = 1$,\n", "+ $a_2 = 1$,\n", "+ for an other $n \\geq 3$, $a_n = a_{n-1} + a_{n-2}$.\n", "\n", "As demonstration let us define a function `fib(n)` which *return*s the $n$-th entry of the Fibonacci sequence, first without recursion, second with recursion." ] }, { "cell_type": "code", "collapsed": false, "input": [ "def fib_nonrec(n) :\n", " if n != int(n) :\n", " print \"Please input an positive integer.\"\n", " retval = None\n", " elif n <= 0 :\n", " print \"Please input a natural number.\"\n", " retval = None\n", " elif n == 1 or n == 2 :\n", " retval = 1\n", " else :\n", " i = 2\n", " fi = 1\n", " fi_1 = 1\n", " \n", " while i < n :\n", " i += 1\n", " fi_2 = fi_1\n", " fi_1 = fi\n", " fi = fi_1 + fi_2\n", " retval = fi\n", " return retval" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 13 }, { "cell_type": "code", "collapsed": false, "input": [ "for i in range(11): \n", " print fib_nonrec(i)\n", "print fib_nonrec(4.5)" ], "language": "python", "metadata": {}, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ "Please input a natural number.\n", "None\n", "1\n", "1\n", "2\n", "3\n", "5\n", "8\n", "13\n", "21\n", "34\n", "55\n", "Please input an positive integer.\n", "None\n" ] } ], "prompt_number": 14 }, { "cell_type": "code", "collapsed": false, "input": [ "def fib_rec(n) :\n", " if n != int(n) or n <= 0 :\n", " print \"Please input a natural number.\"\n", " retval = None\n", " elif n == 1 or n == 2 :\n", " retval = 1\n", " else :\n", " retval = fib_rec(n-1) + fib_rec(n-2)\n", " return retval" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 15 }, { "cell_type": "code", "collapsed": false, "input": [ "for i in range(11) :\n", " print fib_rec(i)\n", "print fib_rec(4.5)" ], "language": "python", "metadata": {}, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ "Please input a natural number.\n", "None\n", "1\n", "1\n", "2\n", "3\n", "5\n", "8\n", "13\n", "21\n", "34\n", "55\n", "Please input a natural number.\n", "None\n" ] } ], "prompt_number": 16 }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Another example\n", "One can define the factorial of a number to be fact(0) = 1, fact(n) = n * fact(n-1)" ] }, { "cell_type": "code", "collapsed": false, "input": [ "def my_fact(n) :\n", " if n != int(n) or n < 0 :\n", " print \"Expected a non-negative integer.\"\n", " retval = None\n", " elif n == 0 :\n", " retval = 1\n", " else :\n", " retval = n * my_fact(n-1)\n", " return retval" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 17 }, { "cell_type": "code", "collapsed": false, "input": [ "print my_fact(100.0), my_fact(-100), my_fact(1.0/100), my_fact(100)" ], "language": "python", "metadata": {}, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ "9.33262154439e+157 Expected a non-negative integer.\n", "None Expected a non-negative integer.\n", "None 93326215443944152681699238856266700490715968264381621468592963895217599993229915608941463976156518286253697920827223758251185210916864000000000000000000000000\n" ] } ], "prompt_number": 18 }, { "cell_type": "heading", "level": 3, "metadata": {}, "source": [ "Now we are ready to do Gaussian elimination" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We shall first do it using recursion. For that we need a function to find the first non-zero column, `c`, and the first row which has a non-zero element in column `c`." ] }, { "cell_type": "code", "collapsed": false, "input": [ "def left_most_non_zero(M) :\n", " \"\"\"Finds the left most non-zero entry in M. Returns (-1, -1) for zero and empty matrices.\"\"\"\n", " \n", " if M == [[]] or M == []:\n", " print \"The matrix is empty.\"\n", " r = -1\n", " c = -1\n", " else :\n", " cell_found = False\n", " present_col = 0\n", " no_of_cols = len(M[0])\n", " no_of_rows = len(M)\n", " \n", " while present_col < no_of_cols and cell_found == False :\n", " present_row = 0\n", " while present_row < no_of_rows and cell_found == False :\n", " if M[present_row][present_col] != 0 :\n", " cell_found = True\n", " r = present_row\n", " c = present_col\n", " present_row += 1\n", " present_col += 1\n", " \n", " if cell_found == False :\n", " r = -1\n", " c = -1\n", " return (r, c)\n", " \n", " " ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 19 }, { "cell_type": "code", "collapsed": false, "input": [ "my_matrix = [[0, 0, 1], [0, 0, 2], [0, 1, 3]]\n", "print left_most_non_zero(my_matrix)" ], "language": "python", "metadata": {}, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ "(2, 1)\n" ] } ], "prompt_number": 20 }, { "cell_type": "markdown", "metadata": {}, "source": [ "Here is the algorithm.\n", "\n", "1. First find the first non-zero column, say $c$, and find a row which has a non-zero entry on that column. Let that row by $r$.\n", "1. Interchange rows $0$ and $r$. (elementary operation 1).\n", "1. Sweep respect to pivot $(0, c)$.\n", "1. Repeat the procedure on the submatrix whose first row is the elements of the 1st row from $c+1$ till the end till the last row from $c+1$ to the last column.\n", "\n", "Thus *to use recursion* we shall write a function which given a matrix $M$ and a coordinate $(r, c)$ *return*s a submatrix as above. We also need a function which will help us to copy the submatrix at the correct position." ] }, { "cell_type": "code", "collapsed": false, "input": [ "def extract_sub_matrix(M, t) :\n", " \"\"\"M : matrix; t = (r, c) is the tuple such that every entry of the submatrix has\n", " coordinates (i, j) with i > r and j > c.\"\"\"\n", " \n", " r = t[0]\n", " c = t[1]\n", " \n", " submat = []\n", " no_rows_M = len(M)\n", " if no_rows_M == 0 :\n", " print \"The matrix is empty.\"\n", " else :\n", " no_cols_M = len(M[0])\n", " for i in range(r+1, no_rows_M) :\n", " current_row = []\n", " for j in range(c+1, no_cols_M) :\n", " current_row.append(M[i][j])\n", " submat.append(current_row)\n", " return submat" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 21 }, { "cell_type": "code", "collapsed": false, "input": [ "new_matrix=[[1,2,3,4],[5,6,7,8],[9,10,11,12]]\n", "print_matrix(new_matrix, \"%3.0f\")\n", "print\n", "print_matrix(extract_sub_matrix(new_matrix, (-1,-1)), \"%3d\")\n" ], "language": "python", "metadata": {}, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ " 1 2 3 4\n", " 5 6 7 8\n", " 9 10 11 12\n", "\n", " 1 2 3 4\n", " 5 6 7 8\n", " 9 10 11 12\n" ] } ], "prompt_number": 22 }, { "cell_type": "code", "collapsed": false, "input": [ "def copy_back_submatrix(M, S) :\n", " \"\"\"Copies submatrix S into M to the lower right of t.\"\"\"\n", " no_rows_M = len(M)\n", " no_rows_S = len(S)\n", "\n", " if no_rows_M == 0 :\n", " print \"M is anyway empty.\"\n", " retM = []\n", " elif no_rows_S == 0 :\n", " print \"Nothing to copy: S is empty\"\n", " retM = M\n", " else :\n", " no_cols_M = len(M[0])\n", " no_cols_S = len(S[0])\n", " \n", " r = no_rows_M - no_rows_S\n", " c = no_cols_M - no_cols_S\n", " if r < 0 or c < 0 :\n", " print \"Matrix S is too big.\"\n", " retM = S\n", " else :\n", " retM = []\n", " for i in range(no_rows_M) :\n", " present_row = []\n", " for j in range(no_cols_M) :\n", " if i < r or j < c :\n", " present_row.append(M[i][j])\n", " else :\n", " present_row.append(S[i-r][j-c])\n", " retM.append(present_row)\n", " return retM\n", " " ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 23 }, { "cell_type": "code", "collapsed": false, "input": [ "print_matrix(new_matrix, \"%3d\")\n", "print \"-\"*15\n", "print_matrix(copy_back_submatrix(new_matrix, [[-2], [-3]]), \"%3d\")" ], "language": "python", "metadata": {}, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ " 1 2 3 4\n", " 5 6 7 8\n", " 9 10 11 12\n", "---------------\n", " 1 2 3 4\n", " 5 6 7 -2\n", " 9 10 11 -3\n" ] } ], "prompt_number": 24 }, { "cell_type": "code", "collapsed": false, "input": [ "def red_to_ech_form(M) :\n", " \"\"\"Reduce matrix M to echelon form.\"\"\"\n", "\n", " # First find the non-zero column.\n", " t = left_most_non_zero(M)\n", " r = t[0]\n", " c = t[1]\n", " \n", " if r == -1 or c == -1 :\n", " retM = M\n", " else :\n", " # Swap row r and row 0 and sweep\n", " retM = M\n", " retM = elem1(retM, r, 0)\n", " retM = sweep(retM, (0, c))\n", " \n", " Mprime = extract_sub_matrix(retM, (0, c))\n", " Mprimeech = red_to_ech_form(Mprime)\n", " retM = copy_back_submatrix(retM, Mprimeech)\n", " \n", " return retM" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 25 }, { "cell_type": "code", "collapsed": false, "input": [ "another_matrix = [[0,.5,1,0],[2,11,0,0],[1,0,0,9]]\n", "print_matrix(another_matrix, \"%10.5f\")\n", "print \"-\"*50\n", "print_matrix(red_to_ech_form(another_matrix), \"%10.5f\")" ], "language": "python", "metadata": {}, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ " 0.00000 0.50000 1.00000 0.00000\n", " 2.00000 11.00000 0.00000 0.00000\n", " 1.00000 0.00000 0.00000 9.00000\n", "--------------------------------------------------\n", "The matrix is empty.\n", "Nothing to copy: S is empty\n", " 1.00000 5.50000 0.00000 0.00000\n", " 0.00000 1.00000 2.00000 0.00000\n", " 0.00000 0.00000 1.00000 0.81818\n" ] } ], "prompt_number": 26 } ], "metadata": {} } ] }