Back substitution

Picking up from the last post on elimination, we represented a system of equations in matrix form as $Ax = b$:

$$ \begin{equation} \begin{bmatrix} 1 & 2 & 3 \\ 2 & 2 & 6 \\ 4 & 5 & 6 \end{bmatrix} % \begin{bmatrix} x_1 \\ x_2 \\ x_3 \end{bmatrix} % =% \begin{bmatrix} 0 \\ -2 \\ 5 \end{bmatrix} \end{equation} $$

We converted this to an equivalent system of equations $Ux = \hat b$ (Note that the same elimination steps applied to $A$ are also applied to the rows of $b$):

$$ \begin{equation} \begin{bmatrix} 1 & 2 & 3 \\ 0 & -2 & 0 \\ 0 & 0 & -6 \end{bmatrix} % \begin{bmatrix} x_1 \\ x_2 \\ x_3 \end{bmatrix} % = \begin{bmatrix} 0 \\ -2 \\ 8 \end{bmatrix} \end{equation} $$

Now we can systematically solve for the components of $x$ with ease. Clearly, $x_3$ must be $\frac{8}{-6} = -1.33$.

Solving for $x_2$ is easy in this example, since there is no $x_3$ coefficient, but if there was we could multiply $x_3$ by the coefficient to arrive at the equation:

$$ \begin{equation} 0 - 2x_2 + (z \times -1.33) = -2 \end{equation} $$

Where $z$ stands for some non-zero coefficient.

Rearranging:

$$ \begin{equation} 0 - 2x_2 = -2 - (z \times -1.33) \end{equation} $$

We then have $-2x_2 = something$ and can solve for $x_2$ easily. In our actual case, we have $-2x_2 = -2$ and so clearly $x_2 = 1$.

The last row is more interesting:

$$ \begin{equation} x_1 + 2x_2 + 3x_3 = 0 \end{equation} $$

Substituting what we know already:

$$ \begin{equation} x_1 + 2 - 4 = 0 \end{equation} $$

And so:

$$ \begin{equation} x_1 = 2 \end{equation} $$

We see that our solutions work for the original system of equations:

$$ \begin{equation} \begin{bmatrix} 1 & 2 & 3 \\ 2 & 2 & 6 \\ 4 & 5 & 6 \end{bmatrix} % \begin{bmatrix} 2 \\ 1 \\ -1.33 \end{bmatrix} % =% \begin{bmatrix} 0 \\ -2 \\ 5 \end{bmatrix} \end{equation} $$

Code implementation

First we run elimination on the augmented matrix $[A | b]$:

def backsub(self, b):
    augmented = cat(self, b, axis=1)
    _, _, _, U, _, _ = augmented.elimination()

Next we initiate a list where we will store the coefficients as they are solved. We loop over the row indices in reverse order. For each step (apart from the first), we take the previous coefficient and multiply the corresponding column of U, then subtract the result from the current $b$. To do this, we use an elimination matrix on the right side of $U$ (to manipulate the columns rather than the rows):

coeff = []
for idx in range(-1, -(U.size(0)+1), -1):
    if idx < -1:
        E = eye([U.size(0)+1, U.size(1)])
        E.data[idx][U.size(1)-1] = -1*(coeff[-1])
        U = U.multiply(E)

With the rearranged row we can attempt to solve for the unknown. There are a number of possibilities. If we have a singular matrix, elimination will have produced a row of zeros and in this case we can't solve for a non-zero in $b$. If $b$ happens to also contain a zero at this position, then any coefficient will multiply to produce a zero and we have an infinite number of solutions. We default to picking a coefficient of $1$ in this case. Otherwise, we solve for the coefficient by dividing the $b$ term with term in the row:

row = U.data[idx]
# check solution possibilities
if row[idx-1] == 0 and row[-1] != 0:
   print('No solution!')
   return None
elif row[idx-1] == 0 and row[-1] == 0:
   print('Infinite solutions!')
   coeff.append(1)
else:
    coeff.append(row[-1]/row[idx-1])

Since we solved for the coefficients in reverse order, we reverse our coefficient list, and return them as a column vector:

coeffs = list(reversed(coeff))
return Mat([coeffs]).tr()

Demo

We create a matrix, call the backsub method and print the result. We also print the result of the multiplication $Ax$ to confirm that we get $b$:

import linalg as la

A = la.Mat([[1, 2, 3],
            [2, 2, 6],
            [4, 5, 6]])

b = la.Mat([[0],
            [-2],
            [5]])

x = A.backsub(b)

la.print_mat(x, 2)
la.print_mat(A.multiply(x))

Outputs:

>>> la.print_mat(x, 2)
[2.0]
[1.0]
[-1.33]

>>> la.print_mat(A.multiply(x))
[0.0]
[-2.0]
[5.0]

< $EA = U$

Pivots, rank, singularity, determinant >

back to project main page
back to home

---