$EA = U$

Solving a set of simultaneous equations (when there is an unique solution) can be done by subtracting multiples of one from another such that an unknown is 'eliminated'. By repeating this procedure, we can eliminate all but one unknown from an equation, leaving a trivially easy solution of the remaining unknown. Having found the value of that unknown, we can move to trivially solving any equation involving two unknowns provided one of them is the one we have previously solved. Continuing this procedure of 'back substitution' will systematically solve for all the unknowns.

For example, the equations: $$ \begin{equation} x_1 + 2x_2 + 3x_3 = 0 \\ 2x_1 + 2x_2 + 6x_3 = -2 \\ 4x_1 + 5x_2 + 6x_3 = 5 \end{equation} $$ can be translated, after elimination, into: $$ \begin{equation} x_1 + 2x_2 + 3x_3 = 0 \\ 0 - 2x_2 + 0 = -2 \\ 0 + 0 - 6x_3 = 8 \end{equation} $$

In linear algebra, the procedure of elimination can be carried out by multiplying a matrix $A$ of coefficients with an elimination matrix $E$. The result is an upper triangular matrix $U$.

The initial problem looks like: $$ \begin{equation} Ax = b =% \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} $$ The elimination procedure looks like: $$ \begin{equation} Ux = EAx = Eb \end{equation} $$ $$ \begin{equation} Ux = % \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} $$

In the case that we find a zero sitting in a column above an unknown that we want to eliminate, we can exchange the row with the zero with some row below that doesn't have a zero in that column.

What we would like to do is encode the results of each elimination step in the matrix $E$ and to have the result $U$ coming from the multiplication of $E$ with $A$.

Code implementation

First we create a few matrices which we'll use later on. The first is the identity matrix $I$, which will become our matrix $E$, but which could have more columns than rows (this allows us to find the inverse, described in a later post). The second is a permutation matrix $P$, which allows us to exchange the rows of a matrix through multiplication:

def elimination(self):
    # we assume the matrix is invertible
    singular = False

    # size of elimination and perumtation matrices
    mat_size = [self.size(0)]*2

    # create identity matrix which we'll turn into an E matrix
    E = eye(mat_size)

    # create a permutation matrix for row exchanges
    P = eye(mat_size)

Each iteration of elimination will produce an $E$ matrix for that step. We continually multiply those $E$'s together to keep track of the overall $E$. The input matrix is copied to $U$, and will become the upper triangular matrix. We also produce a permutation matrix P, that encodes any row exchanges needed.

For each entry on the diagonal of the matrix (the 'pivot positions'), we carry out a row exchange if we find a zero:

U = dc(self)
pivot_count = 0
row_exchange_count = 0
for idx in range(U.size(0)-1):
    for sub_row in range(idx+1, U.size(0)):
        # create elimination mat
        nextE = eye(mat_size)
        nextP = eye(mat_size)

        # handle a zero in the pivot position
        if U.data[idx][pivot_count] == 0:
            row_exchange_count += 1
            # look for a non-zero value to use as the pivot
            options = [row[pivot_count] for row in U.data[sub_row:]]
            exchange = sub_row + options.index(max(options, key=abs))

            # build and apply a purmutation matrix
            nextP.data[idx][pivot_count] = 0
            nextP.data[idx][exchange] = 1
            nextP.data[exchange][exchange] = 0
            nextP.data[exchange][pivot_count] = 1
            U = nextP.multiply(U)
            P = nextP.multiply(P)

If after a possible row exchange we have a zero in the pivot position then the matrix is singular and we flag that fact. We also undo the failed permutation and carry on with elimination in the next column.

# check if the permutation avoided a zero in the pivot position
if U.data[idx][idx] == 0:
    singular = True
    # undo the row exchanges that failed
    row_exchange_count -= 1
    U = nextP.tr().multiply(U)
    P = nextP.tr().multiply(P)
    # move on to the next column
    break

If we have avoided a zero pivot, we determine what multiple of the higher row to subtract from the lower row so as to eliminate the unknown. This multiple goes into our nextE matrix and is applied to $U$:

    # determine how much to subtract to create a zero
    ratio = U.data[sub_row][pivot_count]/U.data[idx][pivot_count]
    # create the elimination matrix for this step
    nextE.data[sub_row][idx] = -ratio
    # apply the elimination step to U
    U = nextE.multiply(U)

We update the overall $E$ by multiplying it with nextE:

    # update the overall E
    E = nextE.multiply(E)
pivot_count += 1

In the case of a 1x1 matrix, the loops in the code above won't happen, and we just take the reciprocal of the number, if it's not zero, to be the pivot. After the final pivot has been placed into $U$, we just check if it was zero, and if so indicate that with a flag on returning the results:

# If self was a 1x1 matrix, the above loops didn't happen. Take the
# reciprocal of the number:
if U.size(0) == 1 and U.size(1) == 2:
    if U.ind(0,0) != 0:
        U.data[0] = [1/U.ind(0,0), 1]
    i = -1

# check if the matrix is square
if U.size(1) == U.size(0):
    # check if the permutation avoided a zero in the pivot position
    if U.data[idx+1][idx+1] == 0:
        singular = True

return P, E, self, U, singular, row_exchange_count

Demo

We create a matrix, call the elimination method and print the result:

import linalg as la

A = la.Mat([[1, 2, 3],
            [4, 8, 6],
            [7, 8, 9]])

P, E, A, U, singular, row_exchange_count = A.elimination()

print(singular)
print()
print(row_exchange_count)
print()
la.print_mat(P)
la.print_mat(E)
la.print_mat(A)
la.print_mat(U)

Outputs:

>>> print(singular)
0
>>> print(row_exchange_count)
1
>>> la.print_mat(P)
[1, 0, 0]
[0, 0, 1]
[0, 1, 0]

>>> la.print_mat(E)
[1.0, 0.0, 0.0]
[-4.0, 1.0, 0.0]
[-7.0, 0.0, 1.0]

>>> la.print_mat(A)
[1, 2, 3]
[4, 8, 6]
[7, 8, 9]

>>> la.print_mat(U)
[1.0, 2.0, 3.0]
[0.0, -6.0, -12.0]
[0.0, 0.0, -6.0]

< Dot product, length and matrix multiplication

Back substitution >

back to project main page
back to home

---