Back-substitution and inverting matricies

Matrix triangulation and back-substitution algorithms can be used in combination with gaussian elimination to solve systems of equations or to find the inverse of a matrix. I previously covered gaussian elimination, continuing on we can now solve the systems of equations using back substitution.

The matrix we had to solve was:

    1     2     1     4    13 
    0    -4     2    -5     2 
    0     0    -5  -7.5   -35 
    0     0     0    -9   -18 

First we normalise the upper-triangle matrix, by simply dividing each row such that the leading coefficient is one:

    1     2     1     4    13 
    0     1  -0.5   1.2  -0.5 
    0     0     1   1.5     7 
    0     0     0     1     2 

(this simplifies the back-substitution, but we can skip/combine this step with the back-substitution)

For back-substitution we work our way backwards from the bottom of the matrix to the top, progressively eliminating each variable. As with gaussian elimination we select a pivot row, and subtract that from the rows above it. First, we start with the last row, and subtract 1.5 times that row from the row above.

    1     2     1     4    13 
    0     1  -0.5   1.2  -0.5 
    0     0     1     0     4 <-- subtract pivot row * 1.5
    0     0     0     1     2 <-- pivot 

Similarly, we continue on to the second row, subtracting 1.2 times, and the top row, subtracting four times.

    1     2     1     0     5 <-- subtract pivot row * 4
    0     1  -0.5     0    -3 <-- subtract pivot row * 1.2
    0     0     1     0     4 
    0     0     0     1     2 <-- pivot

Again, we repeat the process for the third column:

    1     2     0     0     1 
    0     1     0     0    -1 
    0     0     1     0     4 <-- pivot
    0     0     0     1     2 

And finally, the second column:

    1     0     0     0     3 
    0     1     0     0    -1 <-- pivot
    0     0     1     0     4 
    0     0     0     1     2 

Now we have our solution to the system of equations from our original gaussian elimination problem.
a = 3, b = -1, c = 4 and d = 2.
In words/pseudo-code, the process is:

• Pivot through all the rows, starting from the bottom to the top
• For each row above the pivot, calculate how many times we need to subtract the pivot row from this row.
• For each element in the row, subtract the corresponding element from the pivot row, multiplied by the value above.

In code:

for (int p=n-1;p>0;p--) { //pivot backwards through all the rows
        for (int r=p-1;r>=0;r--) { //for each row above the pivot
            float multiple = mat[r][p] / mat[p][p]; //how many multiples of the pivot row do we need (to subtract)?
            for (int c=p-1;c<m;c++) {
                mat[r][c] = mat[r][c] - mat[p][c]*multiple; //subtract the pivot row element (multiple times)
            }
        }
    }

(complete code here)

This process can be applied to find the inverse of a general matrix. Beginning with any matrix we want to invert, we augment it with the identity matrix. For example:

    2     4    -2     1     0     0 
    4     9    -3     0     1     0 
   -2    -3     7     0     0     1 

Now we can apply gaussian elimination to generate:

    2     4    -2     1     0     0 
    0     1     1    -2     1     0 
    0     0     4     3    -1     1 

The normalise the upper triangle to get:

    1     2    -1   0.5     0     0 
    0     1     1    -2     1     0 
    0     0     1  0.75 -0.25  0.25 

And finally, back-substitution to get our solved inverse:

    1     0     0   6.8  -2.8  0.75 
    0     1     0  -2.8   1.2 -0.25 
    0     0     1  0.75 -0.25  0.25 

In this entire discussion I have left out ill-conditioned and singular matrices, but I'll leave modifying the code for that as an exercise for the reader.