## Sunday, November 04, 2012

### 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.

Peterdell said...

Great Explanation !!

Greetings !

Peterdell said...

Great Explanation !!

Greetings !