Bird Tracks Through Math Land: Gaussian Elimination with Back Substitution

Last time we went over some basic matrix operations and this time around we will talk about how Gaussian Elimination works starting with the concepts and working into an implementation using Haskell. Gaussian Elimination is an algorithm that can be used to solve a system of linear equations. Let's forget about matrices for a now and just get a feel for how this works with equations. Suppose you are given the following linear equations:

2x  + y  + z = 1    (L1)
4x  + 3y + z = -1   (L2)
-2x + 2y + z = 7    (L3)

and you are asked to solve for x, y and z. Gaussian Elimination gives you a systematic way of finding the solution. The first major step is to "zero out" the lower triangle of coefficients. So we want the x coefficient from L2 to be 0 and we want the x and y coefficients from L3 to be zero.

Step 1.1: We'll zero out the x coefficient in L2 by subtracting 2*L1 giving us L2'

0x + y - z = -3     (L2')

Step 1.2: We can zero out x in L3 by adding L1 giving us L3'

0x + 3y + 2z = 8    (L3')

Step 1.3: Now we can eliminate y from our modified L3' by subtracting 3*L2' giving us L3'':

0x + 0y + 5z = 17   (L3'')

You'll notice a pattern forming by step 1.3. Because we're multiplying L3' by a factor of L2' both of which already have 0x coefficients we can be sure that any choice we make to zero out the y coefficient will still leave us with a zero x coefficient in L3''.

Lets take a look at our new equations in with zero coefficients under the diagonal (triangular form)

2x + y  + z  = 1    (L1)
0x + y  - z  = -3   (L2')
0x + 0y + 5z = 17   (L3'')

OK cool, now we can use back substitution to solve for z, y and x (in that order).

Step 2.1: Solve for z using L3''

5z = 17
z = 17/5 = 3.4

Step 2.2: Solve for y using L2' and the value of z we just calculated

y - 3.4 = -3
y = 0.4

Step 2.3: Solve for z using L1 along with our calculated z and y:

2x + 0.4 + 3.4 = 1
2x = -2.8
x = -1.4

Done!

Well, sort of... There is one special case that we need to account for. Let's repeat our example after changing the x coefficient in L2 to 6.

2x  + y  + z = 1    (M1)
6x  + 3y + z = -1   (M2)
-2x + 2y + z = 7    (M3)

We can zero out the x column as before

At this point we're supposed to zero out the y coefficient in M3' by adding a multiple of M2' but that's not possible because the y coefficient in M2' is 0. The way we fix this is by having one of the equations below M2' with a non-zero y coefficient take M2's place. In this particular case M3' is our only option.

2x + y  + z  = 1    (M1)
0x + 3y + 2z = 8    (M3')
0x + 0y - 2z = -4   (M2')

And after solving with back substitution:

x = -7/6 = -1.1666
y = 4/3 = 1.333
z = 2

In this case we're done diagonalizing the linear equations (if there were more coefficients to zero out we would just keep going as before).

> module Math.BTTML.GaussianElimination (
>   solveWithGaussAndBackSub,
>   gaussianElimination,
>   backSubstituteUpper,
>   backSubstituteLower) where
>
> import Data.List
> import Math.BTTML.BasicMatrixOps

One last thing before we move into code. All of the coefficients on the left hand side will be represented as a matrix and all of the values on the right hand side will be represented as a column vector. For example, our initial system of linear equations (L1, L2 and L3) can be represented as:

> lCoefMatrix =
>   [[2.0,  1.0, 1.0],
>    [4.0,  3.0, 1.0],
>    [-2.0, 2.0, 1.0]]
>
> lRHSVector = [1.0, -1.0, 7.0]

and our M equations can be represented as

> mCoefMatrix =
>   [[2.0,  1.0, 1.0],
>    [6.0,  3.0, 1.0],
>    [-2.0, 2.0, 1.0]]
>
> mRHSVector = [1.0, -1.0, 7.0]

Now for the code to implement the concepts we just went over.

> solveWithGaussAndBackSub coefMat rhsVec =
>   let (gaussElimCoefMat, gaussElimRHSVec) = gaussianElimination coefMat rhsVec
>   in backSubstituteUpper gaussElimCoefMat gaussElimRHSVec

"zero out" the lower triangle of the given coefficient matrix using the gaussian elimination method

> gaussianElimination :: (Fractional a) => [[a]] -> [a] -> ([[a]], [a])
> gaussianElimination [] [] = ([], [])
> gaussianElimination [] _ = error "Matrix row count is greater than right hand size!"
> gaussianElimination _ [] = error "Matrix row count is less than right hand size!"
> gaussianElimination coefMat rhsVec =
>   let
>       -- work on zeroing out the leftmost coefficients in all equations
>       -- below the top (ie below topCoefs). The ensureNonZeroFirstCoef
>       -- function helps us deal with the special case where the first
>       -- coefficient is zero
>       (topCoefs:otherCoefsMat, topRHSVal:otherRHSVals) =
>           ensureNonZeroFirstCoef coefMat rhsVec
>       zeroOutFactors = map (zeroOutFactor topCoefs) otherCoefsMat
>      
>       zeroOutCoefMatAddends = zipWith (map . (*)) zeroOutFactors (repeat topCoefs)
>       zeroedOutCoefMat = otherCoefsMat `matAdd` zeroOutCoefMatAddends
>      
>       -- now calculate the right-hand-side to go along with the
>       -- zeroed out coefficients
>       zeroOutRHSAddends = map (topRHSVal *) zeroOutFactors
>       zeroedOutRHS = eqZipWith (+) otherRHSVals zeroOutRHSAddends
>      
>       -- now that the leftmost coefs are zero we need to recurse into the
>       -- lower right corner
>       lowerRightMat = map tail zeroedOutCoefMat
>       (newLowerRightMat, newLowerRHS) = gaussianElimination lowerRightMat zeroedOutRHS
>      
>       -- OK now we just have to zip it up like we need
>       newLowerMat = eqZipWith (:) (map head zeroedOutCoefMat) newLowerRightMat
>       finalMat = topCoefs:newLowerMat
>       finalRHS = topRHSVal:newLowerRHS
>   in
>       (finalMat, finalRHS)

this function is just to make sure that the left-most coefficient is non-zero. If the 1st leftmost coefficient is zero in the 1st row then that row is swapped with the 1st row that has a non-zero starting coefficient

> ensureNonZeroFirstCoef coefMat rhsVec =
>   let firstNonZeroIndex = findIndex ((/= 0) . head) coefMat
>   in case firstNonZeroIndex of
>       Just 0  -> (coefMat, rhsVec)
>       Just i  -> (zeroSwap i coefMat, zeroSwap i rhsVec)
>       Nothing -> error "Failed to find non-zero coefficient!"
>   where zeroSwap i xs =
>           let (start, finish) = splitAt i xs
>           in head finish : start ++ tail finish

the zeroOutFactor determines the factor that should be used to zero out the left-most coefficient in otherCoefs

> zeroOutFactor :: (Fractional a) => [a] -> [a] -> a
> zeroOutFactor topCoefs otherCoefs =
>   let topLeftCoef = head topCoefs
>       otherLeftCoef = head otherCoefs
>   in negate $ otherLeftCoef / topLeftCoef

Solve a matrix that is zeroed out below the diagonal using back substitution. This function delegates all the real work to backSubstituteLower.

> backSubstituteUpper upperTriangleMat rhsVec =
>   let
>       -- reverse column and row order in the matrix, then reverse the RHS
>       revMat = reverse $ map reverse upperTriangleMat
>       revRHS = reverse rhsVec
>   in
>       -- now we can back substiture as a lower triangle matrix, reverse
>       -- the result and we're done!
>       reverse $ backSubstituteLower revMat revRHS

solve a matrix that is zeroed out above the diagonal using back substitution

> backSubstituteLower lowerTriangleMat rhsVec = go [] lowerTriangleMat rhsVec
>   where
>       go _ [] [] = []
>       go prevSolutions (matHead:matTail) (rhsHead:rhsTail) =
>           let (solved, (solveNow:solveLater)) = splitAt (length prevSolutions) matHead
>               subRHS = rhsHead - (sum $ eqZipWith (*) prevSolutions solved)
>               solution = subRHS / solveNow
>           in
>               solution:(go (prevSolutions ++ [solution]) matTail rhsTail)
>
>       -- If we haven't hit one of the previous patterns that means that the
>       -- row is different than the RHS count
>       go _ _ _ =
>           error "Matrix row count doesn't match right hand side length!"

That should do it. You can test it by invoking:

> test1 = solveWithGaussAndBackSub lCoefMatrix lRHSVector

or

> test2 = solveWithGaussAndBackSub mCoefMatrix mRHSVector