MathWiki

3. Gaussian elimination algorithm

Gaussian elimination algorithm allows us to solve systems of linear equations.

Let us have a system of linear equations:

$$ \begin{pmatrix} a_{11}&a_{12}&\dots &a_{1n}\\ a_{21}&a_{22}&\dots &a_{2n}\\ \vdots & \vdots &\ddots&\vdots\\ a_{m1}&a_{m2}&\dots &a_{mn} \end{pmatrix} \begin{pmatrix} x_1\\ x_2\\ \vdots\\ x_m \end{pmatrix}= \begin{pmatrix} b_1\\ b_2\\ \vdots\\ b_m \end{pmatrix} $$

and its augmented matrix:

$$ \left(\begin{array}{rrrr|r} a_{11}&a_{12}&\dots &a_{1n}&b_1\\ a_{21}&a_{22}&\dots &a_{2n}&b_2\\ \vdots & \vdots &\ddots&\vdots&\vdots\\ a_{m1}&a_{m2}&\dots &a_{mn}&b_m \end{array}\right). $$

The following operations with a matrix are called elementary row operations:
  1. Switching two rows.
  2. Multiplying a row by a nonzero constant.
  3. Adding to some row a nonzero constant times another row.

The matrix $\begin{pmatrix}-2g & -2h & -2i\\a+d & b+e & c+f\\a & b & c\end{pmatrix}$ has been obtained from the matrix $\begin{pmatrix}a & b & c\\d & e & f\\g & h & i\end{pmatrix}$ by elementary row operations.
Name the elementary row operations used.

Answer

Write an IPython function for (a) adding to a constant $l$ times row $i$ a constant $k$ times row $j$ in a given matrix, return the new matrix; (b) switching rows $i$ and $j$ in a given matrix, return the new matrix.

Solution

Solution 

from numpy import array
 
def add_row_multiples(M, i, j, l, k):      # M is matrix, i is row multiplied by l and add k times row j  
    if (i > len(M) or j > len(M)):         # check do the rows i and j exist in the matrix M
        return 'The row(s) do not exist'
    else:
        R = array(M)                       # the result matrix obtains the value of the initial matrix
        R[i,:] = l*M[i,:] + k*M[j,:]       # row i in R obtains the value l times row i added k times row j
    return R

from numpy import array
 
def swap_rows(M, i, j):                     # M is matrix, rows i and j are swapped
    if (i > len(M) or j > len(M)):          # check whether the rows i and j exist in the matrix M
        return 'The rows do not exist'
    else:
        R = array(M)                        # the result matrix obtains the value of the initial matrix 
        R[[i,j],:]=R[[j,i],:]               # matrix rows i and j are swapped
    return R


In terms of the augmented matrix, the first elementary row operation means that we can reorder equations within the system.
The second elementary row operation means that we can multiply an equation with some constant, i.e.

$ \begin{align*} a_{i1}x_1+\cdots+a_{in}x_n=\beta\qquad\Longleftrightarrow\qquad \gamma a_{i1}x_1+\cdots+\gamma a_{in}x_n=\gamma\beta_i\enspace. \end{align*} $

The third elementary row operation means adding two equations:

$ \begin{align*} \begin{aligned} a_{k1}x_1+\cdots+a_{kn}x_n=\beta_k\\ a_{l1}x_1+\cdots+a_{ln}x_n=\beta_l\\ \end{aligned} \qquad\Longleftrightarrow\qquad (a_{k1}+\gamma a_{l1})x_1+\cdots+(a_{kn}+\gamma a_{ln})x_n=\beta_k+\beta_l.\\ \end{align*} $

Theorem. The elementary row operations deliver an augmented matrix for a system of equations which has the same solution set as the original system.

For example, let

$ \begin{align*} \left(\begin{array}{rrr|r} 1 & 3 & 6 & 25\\ 2 & 7 & 14 & 58\\ 0 & 2 & 5 & 19\\ \end{array}\right) \quad\overset{R_2-2\cdot R_1}\longrightarrow\quad \left(\begin{array}{rrr|r} 1 & 3 & 6 & 25\\ 0 & 1 & 2 & 8\\ 0 & 2 & 5 & 19\\ \end{array}\right)\enspace, \end{align*} $

then

$ \begin{align*} \left\{ \begin{aligned} 1x + 3y + 6z &= 25\\ 2x + 7y + 14z &= 58\\ 2y + 5z &= 19 \end{aligned} \right. \end{align*}\quad $ and $ \quad\begin{align*} \left\{ \begin{aligned} 1x + 3y + 6z &= 25\\ 1y + 2z &= 8\\ 2y + 5z &= 19 \end{aligned} \right. \end{align*} $

have one and the same solution ($x=1,\ y=6,\ z=3$).

3.1 Simple Gaussian elimination algorithm

Let $m=n$ (i.e. the number of equations equals to the number of unknown) and let the system of linear equations have precisely one solution, that is, there is only one set of values $x_1=c_1$, $x_2=c_2$, …, $x_m=c_m$ for which all the equation in the system are true.

Then, the idea of the Gaussian elimination algorithm is

$$ \left(\begin{array}{rrrr|r} a_{11}&a_{12}&\dots &a_{1m}&b_1\\ a_{21}&a_{22}&\dots &a_{2m}&b_2\\ \vdots & \vdots &\ddots&\vdots&\vdots\\ a_{m1}&a_{m2}&\dots &a_{mm}&b_m \end{array}\right)\qquad\overset{\text{elementary row operations}}{\longrightarrow \cdots\longrightarrow}\qquad \left(\begin{array}{rrrr|r} 1&0&\dots &0&c_1\\ 0&1&\dots &0&c_2\\ \vdots & \vdots &\ddots&\vdots&\vdots\\ 0&0&\dots &1&c_m \end{array}\right), $$

i.e. the augmented matrix of a system of linear equations is converted into such diagonal form by using the elementary row operations.

Note that the diagonal form corresponds to the following system of linear equations

$ \begin{align*} \left\{ \begin{aligned} x_1&=c_1\\ x_2&=c_1\\ \cdots&\\ x_n&=c_m \end{aligned}\right. \end{align*} $

and thus, the last column in the diagonal form obtained from the initial augmented matrix by elementary row operations, i.e. $\begin{pmatrix}c_1&\cdots &c_n\end{pmatrix}^T$, is the solution of the original system (see Theorem).

Let us now use Gaussian elimination algorithm for solving the system

$ \begin{align*} \left\{ \begin{aligned} x_2-10x_3&=0\\ 3x_1-x_2-5x_3&=9\\ -2x_1+x_2+3x_3&=-6 \end{aligned}\right. \end{align*}. $

The augmented matrix of this system is

$$ \left(\begin{array}{rrr|r} 0&1&-10&0\\ 3&-1&-5&9\\ -2&1&3&-6 \end{array}\right). $$

STEP 1. (a) Switch the first row of the augmented matrix with rows below till the first element in the first column is nonzero:

$$ \left(\begin{array}{rrr|r} 0&1&-10&0\\ 3&-1&-5&9\\ -2&1&3&-6 \end{array}\right)\quad \overset{\text{switch}\ R_1\ \text{and}\ R_2}{\longrightarrow}\quad \left(\begin{array}{rrr|r} 3&-1&-5&9\\ 0&1&-10&0\\ -2&1&3&-6 \end{array}\right). $$

(b) By using the first element of the first column we zero out the rest of the entries in the first column:

$$ \left(\begin{array}{rrr|r} \color{blue}3&-1&-5&9\\ \color{red}0&1&-10&0\\ \color{red}{-2}&1&3&-6 \end{array}\right) \quad\overset{3R_3+2R_1}{\longrightarrow}\quad \left(\begin{array}{rrr|r} \color{blue}3&-1&-5&9\\ \color{red}0&1&-10&0\\ \color{red}0&1&-1&0 \end{array}\right) $$

Verify the step in IPython by using the functions you created for elementary row operations (e.g., swap_rows and add_row_multiples).

Solution

Solution 

from numpy import array
 
def swap_rows(M, i, j):                                # function for swapping rows i and j in matrix M
    if (i > len(M) or j > len(M)):               
        return 'The rows do not exist'
    else:
        R = array(M)                                    
        R[[i,j],:] = R[[j,i],:]                  
    return R                                               
 
def add_row_multiples(M, i, j, l, k):                  # function for adding l times row i and k times row j
    if (i > len(M) or j > len(M)):               
        return 'The rows do not exist'
    else:
        R = array(M)                                      
        R[i,:] = l*M[i,:] + k*M[j,:]              
    return R
 
 
M = array([[0,1,-10,0],[3,-1,-5,9],[-2,1,3,-6]])       # initial matrix
i = 1
R = array(M)                                           # the result matrix
while (R[0,0] == 0 and i < len(R)):                    # swap rows till the element a_11 is nonzero
    R = swap_rows(R,0,i)
    i = i + 1
print R                                                # matrix with swapped rows
 
for j in range(1, len(R)):                             # loop for making elements a_i1 zero from the row two
    if R[j,0] <> 0:
        l = R[0,0]
        k = R[j,0]
        if cmp(l,0) <> cmp(k,0):                       # choosing constants for row multiplication
            R = add_row_multiples(R,j,0,abs(l),abs(k))
        else: R = add_row_multiples(R, j, 0, l, -k)
print R                                                # matrix where 1st column is zeroed out below 1st row


Solve the task in IPython for $\mathbb{Z}_{11}$.

STEP 2. (a) Now switch the second row with rows below till the second element in the second column is nonzero, i.e. apply STEP 1 (a) for the submatrix obtained by ignoring the first row and the first column:

$$ \left(\begin{array}{rrr|r} \color{lightgray}3&\color{lightgray}{-1}&\color{lightgray}{-5}&\color{lightgray}9\\ \color{lightgray}0&1&-10&0\\ \color{lightgray}0&1&-1&0 \end{array}\right). $$

In our case the second element of the second row is already nonzero, thus no switching is needed.

(b) By using the second element of the second column we zero out the rest of the entries in the second column:

$$ \left(\begin{array}{rrr|r} 3&\color{red}{-1}&-5&9\\ 0&\color{blue}1&-10&0\\ 0&\color{red}1&-1&0 \end{array}\right) \quad\overset{R_1+R_2}{\longrightarrow}\quad \left(\begin{array}{rrr|r} 3&\color{red}{0}&-15&9\\ 0&\color{blue}1&-10&0\\ 0&\color{red}1&-1&0 \end{array}\right) \quad\overset{R_3-R_2}{\longrightarrow}\quad \left(\begin{array}{rrr|r} 3&\color{red}0&-15&9\\ 0&\color{blue}1&-10&0\\ 0&\color{red}0&9&0 \end{array}\right). $$

Verify the step in IPython (like in step 1).

Solution

Solution 

from numpy import array
 
def swap_rows(M, i, j):                              # function for swapping rows i and j in matrix M
    if (i > len(M) or j > len(M)):               
        return 'The rows do not exist'
    else:
        R = array(M)                                    
        R[[i,j],:] = R[[j,i],:]                  
    return R                                               
 
def add_row_multiples(M, i, j, l, k):                # function for adding l times row i and k times row j
    if (i > len(M) or j > len(M)):               
        return 'The rows do not exist'
    else:
        R = array(M)                                      
        R[i,:] = l*M[i,:] + k*M[j,:]              
    return R
 
 
M = array([[3,-1,-5,9],[0,1,-10,0],[0,1,-1,0]])           # initial matrix
i = 2
R = array(M)                                              # the result matrix
while (R[1,1] == 0 and i < len(R)):                       # swap rows till the element a_22 is nonzero
    R = swap_rows(R,1,i)
    i = i + 1
print R                                                   # matrix with swapped rows
 
for j in range(len(R)):                                   # loop for making elements a_i2, except i=2, zero
    if (R[j,1] <> 0 and j <> 1):
        l = R[1,1]
        k = R[j,1]
        if cmp(l,0) <> cmp(k,0):                          # choosing constants for row multiplication
            R = add_row_multiples(R, j, 1, abs(l), abs(k))
        else: R = add_row_multiples(R, j, 1, l, -k)
        print R                                           # matrix printed if row i is changed


Solve the task in IPython for $\mathbb{Z}_{11}$.

STEP 3. The third element in the third row is nonzero (if it were zero, the system would have more then one solution, see General Gaussian elimination algorithm) and can be used to zero out the rest of the elements in the third column:

$$ \left(\begin{array}{rrr|r} 3&0&\color{red}{-15}&9\\ 0&1&\color{red}{-10}&0\\ 0&0&\color{blue}{9}&0 \end{array}\right) \quad\overset{9R_1+15R_3}{\longrightarrow}\quad \left(\begin{array}{rrr|r} 27&0&\color{red}{0}&81\\ 0&1&\color{red}{-10}&0\\ 0&0&\color{blue}{9}&0 \end{array}\right) \quad\overset{9R_2+10R_3}{\longrightarrow}\quad \left(\begin{array}{rrr|r} 27&0&\color{red}0&81\\ 0&9&\color{red}0&0\\ 0&0&\color{blue}{9}&0 \end{array}\right). $$

Verify the step in IPython (like in step 2).

Solution

Solution 

from numpy import array
 
def swap_rows(M, i, j):                              # function for swapping rows i and j in matrix M
    if (i > len(M) or j > len(M)):               
        return 'The rows do not exist'
    else:
        R = array(M)                                    
        R[[i,j],:] = R[[j,i],:]                  
    return R                                               
 
def add_row_multiples(M, i, j, l, k):                # function for adding l times row i and k times row j
    if (i > len(M) or j > len(M)):               
        return 'The rows do not exist'
    else:
        R = array(M)                                      
        R[i,:] = l*M[i,:] + k*M[j,:]              
    return R
 
 
M = array([[3,0,-15,9],[0,1,-10,0],[0,0,9,0]])            # initial matrix
i = 3
R = array(M)                                              # the result matrix
while (R[2,2] == 0 and i < len(R)):                       # swap rows till the element a_33 is nonzero
    R = swap_rows(R,1,i)
    i = i + 1
print R                                                   # matrix with swapped rows
 
for j in range(len(R)):                                   # loop for making elements a_i3, except i=3, zero
    if (R[j,2] <> 0 and j <> 2):
        l = R[2,2]
        k = R[j,2]
        if cmp(l,0) <> cmp(k,0):                          # choosing constants for row multiplication
            R = add_row_multiples(R, j, 2, abs(l), abs(k))
        else: R = add_row_multiples(R, j, 2, l, -k)
        print R                                           # matrix printed for every i if a_i3 is changed


Solve the task in IPython for $\mathbb{Z}_{11}$.

THE RESULT. It remains to multiply each row by the inverse element of the value of the leading entry:

$$ \left(\begin{array}{rrr|r} \color{red}{27}&0&0&81\\ 0&\color{red}{9}&0&0\\ 0&0&\color{red}{9}&0 \end{array}\right) \quad\overset{\frac{1}{27}R_1,\ \frac{1}{9}R_2,\ \frac{1}{9}R_3}{\longrightarrow}\quad \left(\begin{array}{rrr|r} \color{blue}{1}&0&0&3\\ 0&\color{blue}{1}&0&0\\ 0&0&\color{blue}{1}&0 \end{array}\right), $$

and thus, $x_1=3,\ x_2=0,\ x_3=0$.

In general, the $i^{th}$ step of the simple Gaussian elimination algorithm can be described as follows:

(a) if $a_{ii}=0$ and $i>m$ switch the $i^{th}$ row with a row $j$, where $j>i$ till $a_{ii}$ is nonzero;

(b) use $a_{ii}$ to zero out all of the other elements in the $i^{th}$ column.

Write a function in IPython for performing the $i^{th}$ step of the algorithm.

Solution

Solution 

from numpy import array
 
def swap_rows(M, i, j):                              # function for swapping rows i and j in matrix M
    if (i > len(M) or j > len(M)):               
        return 'The rows do not exist'
    else:
        R = array(M)                                    
        R[[i,j],:] = R[[j,i],:]                  
    return R                                               
 
def add_row_multiples(M, i, j, l, k):                # function for adding l times row i and k times row j
    if (i > len(M) or j > len(M)):               
        return 'The rows do not exist'
    else:
        R = array(M)                                      
        R[i,:] = l*M[i,:] + k*M[j,:]              
    return R
 
#############################################################################################################
 
def gaussian_step(M, i):
    if (i < 1 or i >len(M)):
        return 'No such step in this matrix'
 
    R = array(M)                                          # the result matrix
    s = i
    while (R[i-1,i-1] == 0 and s < len(R)):               # swap rows till the element a_ii is nonzero
        R = swap_rows(R, i-1, s)
        s = s + 1
        print R                                           # print matrix with swapped rows
 
    if R[i-1,i-1] == 0:
        print 'Error: the system does not have one solution'
        return numpy.zeros(shape=(len(R),len(R.T)))       # if we continue, we have row multiplication with zero
 
    for j in range(len(R)):                               # loop for making elements in row i, except a_ii, zero
        if (R[j,i-1] <> 0 and j <> (i - 1)):
            l = R[i-1,i-1]
            k = R[j,i-1]
            if cmp(l,0) <> cmp(k,0):                      # choosing constants for row multiplication
                R = add_row_multiples(R, j, i-1, abs(l), abs(k))
            else : R = add_row_multiples(R, j, i-1, l, -k)
            print R                                       # matrix printed every time it changes
    return R


Write a function in IPython for solving a system of linear equation with the simple Gaussian elimination algorithm (mind the assumptions for the algorithm; use the function created in the previous exercise (e.g., gaussian_step)).

Solution

Solution 

from numpy import array
 
def swap_rows(M, i, j):                              # function for swapping rows i and j in matrix M
    if (i > len(M) or j > len(M)):               
        return 'The rows do not exist'
    else:
        R = array(M)                                    
        R[[i,j],:] = R[[j,i],:]                  
    return R                                               
 
def add_row_multiples(M, i, j, l, k):                # function for adding l times row i and k times row j
    if (i > len(M) or j > len(M)):               
        return 'The rows do not exist'
    else:
        R = array(M)                                      
        R[i,:] = l*M[i,:] + k*M[j,:]              
    return R
 
def gaussian_step(M, i):                             # function for step i in gaussian elimination algorithm
    if (i < 1 or i >len(M)):
        return 'No such step in this matrix'
    R = array(M)
    s = i
    while (R[i-1,i-1] == 0 and s < len(R)):
        R = swap_rows(R, i-1, s)
        s = s + 1
        print R
    if R[i-1,i-1] == 0:
        print 'Error: the system does not have one solution'
        return numpy.zeros(shape=(len(R),len(R.T)))
    for j in range(len(R)):
        if (R[j,i-1] <> 0 and j <> (i - 1)):
            l = R[i-1,i-1]
            k = R[j,i-1]
            if cmp(l,0) <> cmp(k,0):
                R = add_row_multiples(R, j, i-1, abs(l), abs(k))
            else : R = add_row_multiples(R, j, i-1, l, -k)
            print R
    return R
 
############################################################################################################
 
def simple_gaussian(M): 
    R = array(M)                                  # the result matrix
    for i in range(len(R)):
        R = gaussian_step(R,i+1)
        if (R == numpy.zeros(shape=(len(R),len(R.T)))).all() :
            print 'Error: not solvable with this method'
            return
 
    T = R.astype(numpy.float)                    # needed when initial matrix is int, but solution is float
    for j in range(len(R)):
        T[j,:] = (float(1)/R[j,j])*R[j,:]        # multiplying matrix rows so that on diagonal we get 1
        print T
    return T[:,len(R)]                           # returns last column of the matrix, i.e. the solution


Find a solution to the equation $A\mathbf{x}=\mathbf{b}$ when

$$\begin{align*}A=\begin{pmatrix}1 & 1 & 1\\1 & 2 & 4\\1 & 3 & 9\end{pmatrix}\end{align*}$$ and $\mathbf{b}=(1\ 1\ 1)^T$ over $\mathbb{Z}_2$, $\mathbb{Z}_3$ and $\mathbb{Z}_5$.

3.2 General Gaussian elimination algorithm

Let us now use the simple Gaussian elimination algorithm for the system

$ \begin{align*} \left\{ \begin{aligned} 2x_1-3x_2+5x_3&=1\\ 4x_1-6x_2+2x_3&=2\\ 2x_1-3x_2-11x_3&=1\\ \end{aligned}\right.\enspace. \end{align*} $

The first step is

$$ \left(\begin{array}{rrr|r} 2&-3&5&1\\ 4&-6&2&2\\ 2&-3&-11&1 \end{array}\right) \quad\overset{R_2-2R_1}{\longrightarrow}\quad \left(\begin{array}{rrr|r} 2&-3&5&1\\ 0&0&-8&0\\ 2&-3&-11&1 \end{array}\right) \quad\overset{R_3-R_2}{\longrightarrow}\quad \left(\begin{array}{rrr|r} 2&-3&5&1\\ 0&0&-8&0\\ 0&0&-16&0 \end{array}\right). $$

The second step can not be performed, since $a_{22}$ is zero and it can not be changed by switching second row with any row below:

$$ \left(\begin{array}{rrr|r} 2&-3&5&1\\ 0&\color{red}0&-8&0\\ 0&\color{red}0&-16&0 \end{array}\right). $$

So lets skip the second step.

The third step is again standard:

$$ \left(\begin{array}{rrr|r} 2&-3&5&1\\ 0&0&-8&0\\ 0&0&-16&0 \end{array}\right) \quad\overset{16R_1+5R_3}{\longrightarrow}\quad \left(\begin{array}{rrr|r} 32&-48&0&16\\ 0&0&-8&0\\ 0&0&-16&0 \end{array}\right) \quad\overset{2R_2-R_3}{\longrightarrow}\quad \left(\begin{array}{rrr|r} 32&-48&0&16\\ 0&0&0&0\\ 0&0&-16&0 \end{array}\right). $$

Note, that the equation $0x_1+0x_2+0x_3=0$ holds always and hence, does not change the solution of the system. Thus, we may disregard the row and thus, have obtain the system

$ \begin{align*} \left\{ \begin{aligned} 32x_1-48x_2&=16\\ -16x_3&=0\\ \end{aligned}\right.\enspace, \end{align*} $

which is equivalent to

$ \begin{align*} \left\{ \begin{aligned} x_1&=\frac{16}{32}+\frac{48}{32}x_2\\ x_3&=0\\ \end{aligned}\right.\enspace. \end{align*} $

Note, that $x_3$ has a concrete value, $x_1$ has a criteria for evaluating it, and $x_2$ has no boundaries for valuations ($x_2$ is a free variable in this system). Clearly, there isn't just one set of values $x_1=c_1,\ x_2=c_2,\ x_3=c_3$ which satisfy the system of equations, that is, the system has no unique solution. In such a case, the solution is described as a set

$$\left\{\left(\frac{16}{32}+\frac{48}{32}c_2,\ c_2,\ 0\right): c_2\in\mathbb{R}\right\}.$$

By giving $c_2$ a concrete value, for example, let us take $c_2=0$, we obtain one solution of the system: $(\frac{16}{32},0,0)$, but it is not the only solution. Actually, the system has infinitely many solutions, since we can choose infinitely many different values for $c_2$.

Use the IPython function you created for simple Gaussian elimination algoritm (e.g., simple_gaussian) for solving the system

$ \begin{align*} \left\{ \begin{aligned} 2x_1-3x_2+5x_3&=1\\ 4x_1-6x_2+2x_3&=2\\ 2x_1-3x_2-11x_3&=1\\ \end{aligned}\right.\enspace. \end{align*} $

Analyse the result: what should be changed in order to obtain the final augmented matrix from above.

Solution

Solution 

from numpy import array
 
def swap_rows(M, i, j):                         # function for swapping rows i and j in matrix M
    if (i > len(M) or j > len(M)):               
        return 'The rows do not exist'
    else:
        R = array(M)                                    
        R[[i,j],:] = R[[j,i],:]                  
    return R                                               
 
def add_row_multiples(M, i, j, l, k):           # function for adding l times row i and k times row j
    if (i > len(M) or j > len(M)):               
        return 'The rows do not exist'
    else:
        R = array(M)                                      
        R[i,:] = l*M[i,:] + k*M[j,:]              
    return R
 
def changed_gaussian_step(M, i):                # function for step i in simple gaussian elimination
    if (i < 1 or i >len(M)):
        return 'No such step in this matrix'
    R = array(M)
    if any(R[range(i-1,len(R)),i-1] <> 0):      # here is a change 
        s = i
        while (R[i-1,i-1] == 0 and s < len(R)):               
            R = swap_rows(R, i-1, s)
            s = s + 1
            print R
        for j in range(len(R)):
            if (R[j,i-1] <> 0 and j <> (i - 1)):
                l = R[i-1,i-1]
                k = R[j,i-1]
                if cmp(l,0) <> cmp(k,0):
                    R = add_row_multiples(R, j, i-1, abs(l), abs(k))
                else : R = add_row_multiples(R, j, i-1, l, -k)
                print R
    else : print R
    return R
 
###########################################################################################################
 
# it is possible to get the final augmented matrix with the changes 
 
def changed_simple_gaussian(M): 
    R = array(M)
    for i in range(len(R)):
        R = changed_gaussian_step(R,i+1)                        # here is a change
        if (R == numpy.zeros(shape=(len(R),len(R.T)))).all() :
            print 'Error: not solvable with this method'
            return
 
    T = R.astype(numpy.float)
    for j in range(len(T)):
        if R[j,j] <> 0: T[j,:] = (float(1)/R[j,j])*R[j,:]       # here is a change
        print T
    return T
 
M = array([[2,-3,5,1],[4,-6,2,2],[2,-3,-11,1]])
changed_simple_gaussian(M)

Now apply the changed simple Gaussian elimination to the system

$ \begin{align*} \left\{ \begin{aligned} 2x_1-3x_2-x_3+5x_4&=1\\ 4x_1-6x_2-2x_3+2x_4&=2\\ 2x_1-3x_2-x_3-11x_4&=1\\ \end{aligned}\right.\enspace. \end{align*} $

Do you get the correct answer? Answer


As you can see, the simple Gaussian elimination algorithm does not work in such a case. Next we will generalize the simple Gaussian elimination algorithm in order to be able to solve any (solvable) system.

So, let us have a system of linear equation without any restrictions to $m$ or $n$ or number of solutions

$$ \begin{pmatrix} a_{11}&a_{12}&\dots &a_{1n}\\ a_{21}&a_{22}&\dots &a_{2n}\\ \vdots & \vdots &\ddots&\vdots\\ a_{m1}&a_{m2}&\dots &a_{mn} \end{pmatrix} \begin{pmatrix} x_1\\ x_2\\ \vdots\\ x_m \end{pmatrix}= \begin{pmatrix} b_1\\ b_2\\ \vdots\\ b_m \end{pmatrix}. $$

The goal of the general Gaussian elimination algorithm is to give to the augmented matrix of the linear system the following form


by using elementary row operations and leaving out rows consisting only of zeros ($r\leq m$).

In the obtained matrix, if $r=n$ and $a_{ii}=1$ for every $i=1,\dots,n$, then the system has a unique solution. The solution is the last column of the augmented matrix, see simple Gaussian elimination algorithm.

If $r<n$, then the system has more than one solution. In the case of solving such a system over $\mathbb{R}$ (i.e. the solution of the system is searched among real numbers), there are usually infinitely many solutions (see the example) which can not all be listed. There exists theory for giving the general formula for the solution, however, we will not go into it and just give one particular solution of the system.
In the case of solving such a system over a finite field, e.g., $\mathbb{Z}_n$, where $n$ is a prime number, all the solution can be listed ($\mathbb{Z}_n$ has finite number of elements).

If for any row $i$ the only nonzero element is in the last column of the augmented matrix, i.e. in row $i$ only $c_i\neq 0$, then the system has no solution at all. Because then the $i^{th}$ row re-presence the equation

$0x_{1}+0x_{2}+\dots+0x_{n}=c_{i}\quad\Longleftrightarrow\quad 0=c_{i}$

which is never true ($c_{i}$ is nonzero).

Which of the following augmented matrices are the results of general Gaussian elimination algorithm. If some matrix is the result of general Gaussian elimination algorithm, then decide weather it has unique solution or not.

$\left(\begin{array}{rrrrr|r}1&0&0&5&8&0\\0&0&1&2&7&0\\0&0&0&0&0&1\end{array}\right)$ $\left(\begin{array}{rrrrr|r}1&0&6&5&8&2\\0&0&2&2&7&3\end{array}\right)$
Answer Answer
$\left(\begin{array}{rrr|r}1&0&0&3\\0&1&0&-2\\0&0&1&10\\0&0&0&0\end{array}\right)$ $\left(\begin{array}{rrrrr|r}1&0&0&0&0&3\\0&0&9&1&0&-3\\0&1&-2&0&0&2\\0&0&0&0&1&5\end{array}\right)$
Answer Answer 
Note that this time we do not move within the augmented matrix along the diagonal elements a_ii.
However, the simple Gaussian elimination is a special case of the general Gaussian elimination:
just assume, that all the blocks of arbitrary numbers are with 0 width.
The general Gaussian elimination step:

row $i$ and column $j$

if there exists row $k$, $k\leq j$ such that $a_{kj}\neq 0$ then switch rows till $a_{ij}$ is nonzero and by using the new $a_{ij}$ zero out all the rest of the elements in column $j$; next step row $i+1$ and column $j+1$

else next step row $i$ and column $j+1$.

Write a function in IPython for the general Gaussian elimination step.

Solution

Solution 

from numpy import array
 
def swap_rows(M, i, j):                          # function for swapping rows i and j
    if (i > len(M) or j > len(M)):               
        return 'The rows do not exist'
    else:
        R = array(M)                                    
        R[[i,j],:] = R[[j,i],:]                  
    return R                                               
 
def add_row_multiples(M, i, j, l, k):            # function for adding l times row i and k times row j
    if (i > len(M) or j > len(M)):               
        return 'The rows do not exist'
    else:
        R = array(M)                                      
        R[i,:] = l*M[i,:] + k*M[j,:]              
    return R
 
###########################################################################################################
 
def general_gaussian_step(M, i, j):                   # i is row index and j is column index: i,j=0,1,...
    print M
 
    if (i < 0 or i > len(M)-1 or j < 0 or j > len(M.T)-2):
        return 'No such step in this matrix'
 
    R = array(M)                                      # the result matrix
    if any(R[range(i,len(R)),j] <> 0):                # if nonzero elements in column j starting from row i
        s = i + 1
        while (R[i,j] == 0 and s < len(R)):           # switched rows till the element a_ij is nonzero
            R = swap_rows(R, i, s)
            s = s + 1
            print R
        for t in range(len(R)):                       # zeros the elements below and above
            if (R[t,j] <> 0 and t <> i):
                l = R[i,j]
                k = R[t,j]
                if cmp(l,0) <> cmp(k,0):
                    R = add_row_multiples(R, t, i, abs(l), abs(k))
                else : R = add_row_multiples(R, t, i, l, -k)
                print R
    return R


Write a function in IPython for solving a system of linear equations (no restriction to the system here) with the Gaussian elimination algorithm.

Solution

Solution 

from numpy import array
 
def swap_rows(M, i, j):                              # function for swapping rows i and j
    if (i > len(M) or j > len(M)):               
        return 'The rows do not exist'
    else:
        R = M                                    
        R[[i,j],:] = R[[j,i],:]                  
    return R                                               
 
def add_row_multiples(M, i, j, l, k):                # function for adding l times row i and k times row j
    if (i > len(M) or j > len(M)):               
        return 'The rows do not exist'
    else:
        R = M                                      
        R[i,:] = l*M[i,:] + k*M[j,:]              
    return R
 
def general_gaussian_step(M, i, j):                  # function for general gaussian step, row i and col j
    if (i < 0 or i > len(M)-1 or j < 0 or j > len(M.T)-2):
        return 'No such step in this matrix'
 
    R = array(M) 
    if any(R[range(i,len(R)),j] <> 0):
        s = i + 1
        while (R[i,j] == 0 and s < len(R)): 
            R = swap_rows(R, i, s)
            s = s + 1
            print R
        for t in range(len(R)): 
            if (R[t,j] <> 0 and t <> i):
                l = R[i,j]
                k = R[t,j]
                if cmp(l,0) <> cmp(k,0):
                    R = add_row_multiples(R, t, i, abs(l), abs(k))
                else : R = add_row_multiples(R, t, i, l, -k)
                print R
    return R
 
#############################################################################################################
 
def general_gaussian(M): 
    R = array(M)                                      # the result matrix
    i = 0                                             # row index
    j = 0                                             # column index
    while (j < len(M.T)-1 and i < len(R)):
        R = general_gaussian_step(R,i,j)              # general gaussian step for row i and column j
        if R[i,j] <> 0: i = i+1                       # criteria for changing row index
        j = j+1                                       # in every step changes column index
 
    rowind = array(range(len(R)))                     # vector of row indexes in R
    R = R[(rowind < i) | (R[:,len(R.T)-1] <> 0),:]    # leaving out rows consisting only of zeros
    print R
 
    if (all(R[len(R)-1,range(len(R.T)-1)] == 0)):     # enough to check the last row for contradiction
        print 'The system has no solution'
        return
 
    R = R.astype(numpy.float)                              
    t = 0                                             # row index
    s = 0                                             # column index
    while (s < len(R.T) and t < len(R)):
        if R[t,s] <> 0:                               # if the leading element is not zero
            R[t,:] = (float(1)/R[t,s])*R[t,:]         # the leading el. one by mult. the row with the inverse 
            print R
            t = t+1
        s = s+1
 
    if len(R) == len(R.T)-1:                          # precisely one solution, return the last column
        print R[:,len(R.T)-1]
        return
    else : 
        print 'The system has more then one solution'
        S = [0]*(len(R.T)-1)                          # we take the free variables equal to 0
        for m in range(len(R)):
            for n in range(len(R.T)):
                if R[m,n] == 1:                       #if leading el., the variable equals to constant term 
                    S[n] = R[m,len(R.T)-1]
                    break
        print 'One of the solutions is', S
        return


Find all solutions of the system $A\mathbf{x}=\mathbf{b}$ where $$A=\begin{pmatrix}1&1&0\\0&0&1\\0&0&2\end{pmatrix}$$ and $\mathbf{b}=\begin{pmatrix}1&2&4\end{pmatrix}^T$ over $\mathbb{Z}_3$.

4. Solving systems of linear equations via invertible matrices

linear_algebra/03_gaussian_elimination.txt · Last modified: 2014/12/14 00:12 by jaan