{ "metadata": { "name": "" }, "nbformat": 3, "nbformat_minor": 0, "worksheets": [ { "cells": [ { "cell_type": "heading", "level": 1, "metadata": {}, "source": [ "Matrices, Inverses, Equations, Solution" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Today I'll be using lists as I am still not sure if numpy and scipy are installed in the lab." ] }, { "cell_type": "heading", "level": 2, "metadata": {}, "source": [ "Matrices, reading from a file." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "It will be painful to input a 16x16 matrix whenever we need to. So we shall load matrices from files. We shall use the split function (*method*) we learnt before. We shall write a function which takes the name of file as input and returns a list of lists of numbers." ] }, { "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": "markdown", "metadata": {}, "source": [ "The output isn't really very pretty. Let us make it pretty by" ] }, { "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": "heading", "level": 2, "metadata": {}, "source": [ "Coding elementary row operations" ] }, { "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": [ "Note lists are also passed by reference. So modifying matrix actually changes the actual matrix. You can see it below." ] }, { "cell_type": "code", "collapsed": false, "input": [ "print_matrix(elem1(A, 1, 2), \"%2.0f\")\n", "print \"-------------------------\"\n", "print_matrix(A, \"%2.0f\")" ], "language": "python", "metadata": {}, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ " 1 2 3\n", " 3 2 1\n", " 2 1 2\n", "-------------------------\n", " 1 2 3\n", " 3 2 1\n", " 2 1 2\n" ] } ], "prompt_number": 6 }, { "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": 7 }, { "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", " 3 2 1\n", " 2 1 2\n", "--------------------------\n", " 4 8 12\n", " 3 2 1\n", " 2 1 2\n", "-------------------------\n", " 4 8 12\n", " 3 2 1\n", " 2 1 2\n" ] } ], "prompt_number": 8 }, { "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": 9 }, { "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", " 3 2 1\n", " 2 1 2\n", "--------------------\n", "-2 4 10\n", " 3 2 1\n", " 2 1 2\n", "---------------------\n", "-2 4 10\n", " 3 2 1\n", " 2 1 2\n" ] } ], "prompt_number": 10 }, { "cell_type": "heading", "level": 2, "metadata": {}, "source": [ "Sweeping a column" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let us read another matrix $B$. This is the same one as in the class last week." ] }, { "cell_type": "code", "collapsed": false, "input": [ "B = read_matrix(\"files/B.txt\")\n", "print_matrix(B, \"%2.0f\")" ], "language": "python", "metadata": {}, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ " 0 2 3 2\n", " 1 4 5 5\n", " 3 0 2 6\n" ] } ], "prompt_number": 11 }, { "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": 12 }, { "cell_type": "code", "collapsed": false, "input": [ "print_matrix(sweep(B, (1, 2)), \"%6.2f\")" ], "language": "python", "metadata": {}, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ " -0.60 -0.40 0.00 -1.00\n", " 0.20 0.80 1.00 1.00\n", " 2.60 -1.60 0.00 4.00\n" ] } ], "prompt_number": 13 }, { "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": 14 }, { "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": 15 }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now we are ready to do Gaussian elimination" ] }, { "cell_type": "code", "collapsed": false, "input": [ "def gauss_elim_prelim(matrix) :\n", " row = 0\n", " row_offset = 0\n", " col_offset = 0\n", " \n", " (r, c) = first_non_zero(matrix)\n", " if r == -1 :\n", " print \"Nothing to do: The matrix is zero.\"\n", " else :\n", " matrix = elem1(matrix, 0, r)\n", " matrix = sweep(matrix, (0, c))\n", " print_matrix(matrix, \"%6.2f\")" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 16 }, { "cell_type": "code", "collapsed": false, "input": [ "gauss_elim_prelim([[0,0,2,2],\\\n", " [1,4,4,4],\\\n", " [0,1,3,3]])" ], "language": "python", "metadata": {}, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ " 1.00 4.00 4.00 4.00\n", " 0.00 0.00 2.00 2.00\n", " 0.00 1.00 3.00 3.00\n" ] } ], "prompt_number": 17 }, { "cell_type": "code", "collapsed": false, "input": [ "def gauss_elim(matrix) :\n", " row = 0\n", " row_offset = 0\n", " col_offset = 0\n", " \n", " while row < len(matrix) :\n", " submatrix = []\n", " for i in range(len(matrix) - row_offset) :\n", " row_submatrix = []\n", " for j in range(len(matrix[0]) - col_offset) :\n", " row_submatrix.append(matrix[i + row_offset][j + col_offset])\n", " submatrix.append(row_submatrix)\n", " #print_matrix(submatrix, \"%6.2f\")\n", " print row, row_offset, col_offset, submatrix, matrix\n", " print \"-\"*20\n", " \n", " (r, c) = first_non_zero(submatrix)\n", " if r == -1 :\n", " print \"Nothing to do: The (sub)matrix is zero.\"\n", " else :\n", " r += row_offset\n", " c += col_offset\n", " matrix = elem1(matrix, row_offset, r)\n", " matrix = sweep(matrix, (row_offset, c))\n", " row_offset += 1\n", " col_offset = c+1\n", " row += 1\n", " return matrix" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 18 }, { "cell_type": "code", "collapsed": false, "input": [ "read_matrix(\"files/B.txt\")\n", "print_matrix(gauss_elim(B), \"%6.2f\")\n", "print_matrix(gauss_elim([[0,0,0,0], [0,1,1,0], [0,0,1,0]]), \"%6.2f\")" ], "language": "python", "metadata": {}, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ "0 0 0 [[-0.6000000000000001, -0.40000000000000036, 0.0, -1.0], [0.2, 0.8, 1.0, 1.0], [2.6, -1.6, 0.0, 4.0]] [[-0.6000000000000001, -0.40000000000000036, 0.0, -1.0], [0.2, 0.8, 1.0, 1.0], [2.6, -1.6, 0.0, 4.0]]\n", "--------------------\n", "1 1 1 [[0.6666666666666666, 1.0, 0.6666666666666667], [-3.333333333333335, 0.0, -0.33333333333333304]] [[1.0, 0.6666666666666672, -0.0, 1.6666666666666665], [0.0, 0.6666666666666666, 1.0, 0.6666666666666667], [0.0, -3.333333333333335, 0.0, -0.33333333333333304]]\n", "--------------------\n", "2 2 2 [[5.000000000000002, 3.0000000000000018]] [[1.0, 0.0, -1.0000000000000009, 0.9999999999999993], [0.0, 1.0, 1.5, 1.0], [0.0, 0.0, 5.000000000000002, 3.0000000000000018]]\n", "--------------------\n", " 1.00 0.00 0.00 1.60\n", " 0.00 1.00 0.00 0.10\n", " 0.00 0.00 1.00 0.60\n", "0 0 0 [[0, 0, 0, 0], [0, 1, 1, 0], [0, 0, 1, 0]] [[0, 0, 0, 0], [0, 1, 1, 0], [0, 0, 1, 0]]\n", "--------------------\n", "1 1 2 [[0.0, 0.0], [1.0, 0.0]] [[0.0, 1.0, 1.0, 0.0], [0.0, 0.0, 0.0, 0.0], [0.0, 0.0, 1.0, 0.0]]\n", "--------------------\n", "2 2 3 [[0.0]] [[0.0, 1.0, 0.0, 0.0], [0.0, 0.0, 1.0, 0.0], [0.0, 0.0, 0.0, 0.0]]\n", "--------------------\n", "The matrix is the zero matrix.\n", "Nothing to do: The (sub)matrix is zero.\n", " 0.00 1.00 0.00 0.00\n", " 0.00 0.00 1.00 0.00\n", " 0.00 0.00 0.00 0.00\n" ] } ], "prompt_number": 19 } ], "metadata": {} } ] }