Sunday, November 04, 2012

Gaussian Elimination

Gaussian Elimination is an elementary transformation that converts a matrix into a triangle, or row-reduced echelon form (RREF). It forms the basis of a number of operations in linear algebra to solve systems of equations, invert matrices, and minimize systems of equations among other things (I'll cover these in later posts). The Gaussian Elimination algorithm itself is straight-forward (you probably learnt it in high school). Given a system of equations, e.g.
  a + 2b +  c + 4d = 13
 2a +      4c + 3d = 28
 4a + 2b + 2c +  d = 20
-3a +  b + 3c + 2d = 6
We can form an augmented matrix to represent it, and use Gaussian elimination to solve it. The goal is to produce a triangle-matrix representation, so that we can solve the equations by back-substitution. In other words, we want to (eventually) have one row represent each variable, and for all other rows, that variable should be zero. (i.e. solved). Gaussian elimination takes us part of the way there by giving us a set of equations with a starting point which we can then later solve.

Representing the above equations as a matrix, we have:
    1     2     1     4    13 
    2     0     4     3    28 
    4     2     2     1    20 
   -3     1     3     2     6 
The first step is to select a pivot row, which we can use to eliminate/reduce the other rows. When we eliminate the other rows, we want the that variables value to be 0. In this example, we pick the first row, and then subtract that twice from the row below, to ensure that the row below will have zero a's.
    1     2     1     4    13 <-- pivot
    0    -4     2    -5     2 <-- subtract pivot row * 2
    4     2     2     1    20 
   -3     1     3     2     6 
Likewise, four times the third tow, and negative three times the final row.
    1     2     1     4    13 <-- pivot
    0    -4     2    -5     2 
    0    -6    -2   -15   -32 <-- subtract pivot row * 4
    0     7     6    14    45 <-- subtract pivot row * -3
Great. Our first variable (a) has been eliminated. We now repeat this step, starting from the second row, with the variable 'b'. We don't want to use the first row, as we want to preserve that row's representation of the 'a' variable.
    1     2     1     4    13 
    0    -4     2    -5     2 <-- pivot
    0     0    -5  -7.5   -35 <-- subtract pivot row * 1.5
    0     0   9.5   5.2    48 <-- subtract pivot row * -1.75
Now, we repeat the process again, starting from the third row.
    1     2     1     4    13 
    0    -4     2    -5     2 
    0     0    -5  -7.5   -35 <-- pivot
    0     0     0    -9   -18 <-- subtract pivot row * -1.9
Done. In pseudo-code/words, the algorithm is:
  • For each row (except the last), select a pivot. (In my example, I just take the first available row each time)
  • For each row that is below the pivot, calculate the number of times we need to subtract the row (i.e. divide)
  • For each element in this row, subtract the corresponding element in the pivot row, multiplied by the value we calculated above.
The code to achieve this is:
//input a m (col) by n (row) matrix ('mat')
    //p is the pivot - which row we will use to eliminate
    for (int p=0;p<n-1;p++) { //pivot through all the rows
        for (int r=p+1; r < n; r++) { //for each row that isn't the pivot
            float multiple = mat[r][p] / mat[p][p]; //how many multiples of the pivot row do we need (to eliminate this row)?
            for (int c = 0; c<m; c++) { //for each element in this row
                mat[r][c] = mat[r][c] - mat[p][c]*multiple; //subtract the pivot row element (multiple times)
            }
             
        }
    }
(full code here) Next time, we continue on to solve the equations - available here! (2,4,-1,3)

No comments: