How to Perform Gaussian Elimination Using Pivoting in Python

Gaussian elimination is also known as the row reduction method. It is an algorithm commonly used to solve linear problems. The algorithm involves a series of row operations on a matrix of coefficients drawn from the linear equations until the matrix is reduced to echelon form.

The following three operations are performed when transforming a matrix into an echelon form. These include:

• Multiplying one row with a scalar and adding it or subtracting it from another row.
• Swapping rows.
• Multiplying rows with a scalar.

A matrix is in row echelon form if the first element in each row, also known as the leading entry, is a zero; the leading entry in each row is one column to the right of the leading entry in the previous row. Furthermore, rows with zero elements should be below elements having non-zero rows.

Row operations are largely classified into two. Forward elimination involves reducing a matrix to echelon form to determine whether a viable solution exists and if it is finite or infinite. On the other hand, back substitution further reduces the matrix to row reduced echelon form.

Gaussian Elimination With Pivoting in Python

Pivoting is the interchange of rows and columns to get the suitable pivot element. A suitable pivot element should both be non-zero and significantly large but smaller when compared to the other row entries.

Pivoting is classified into partial pivoting and complete pivoting. Under the partial pivoting algorithm, the largest element is considered the pivot element to minimize rounding off errors.

On the other hand, complete pivoting includes the interchange of rows and columns to get the best pivot element, thus increasing accuracy.

A set of equations is considered linear if no variable has an exponent of more than one. The Gaussian elimination involves a series of steps; the first step consists of creating a coefficient matrix.

A coefficient matrix is simply a matrix of coefficients drawn from a set of linear equations. The next step involves creating an augmented matrix that is then subjected to a series of operations that reduce it into echelon form.

However, in cases where the pivot element is zero or of a very small magnitude, we have to interchange the pivot row with a lower row.

While implementing the Gaussian elimination method, we need to note the array indexes. Python iterables such as lists and arrays often start at index 0 and end at index n-1.

We can then read the contents of the augmented matrix, apply the elimination method, back substitution, and finally display the answer.

from numpy import array, zeros, fabs, linalg

a = array(
    [
        [0, 6, -1, 2, 2],
        [0, 3, 4, 1, 7],
        [5, 1, 0, 3, -1],
        [3, 1, 3, 0, 2],
        [4, 4, 1, -2, 1],
    ],
    float,
)
# the b matrix constant terms of the equations
b = array([5, 7, 2, 3, 4], float)

print("Solution by NumPy:")


print(linalg.solve(a, b))

n = len(b)
x = zeros(n, float)

# first loop specifys the fixed row
for k in range(n - 1):
    if fabs(a[k, k]) < 1.0e-12:

        for i in range(k + 1, n):
            if fabs(a[i, k]) > fabs(a[k, k]):
                a[[k, i]] = a[[i, k]]
                b[[k, i]] = b[[i, k]]
                break

    # applies the elimination below the fixed row

    for i in range(k + 1, n):
        if a[i, k] == 0:
            continue

        factor = a[k, k] / a[i, k]
        for j in range(k, n):
            a[i, j] = a[k, j] - a[i, j] * factor
            # we also calculate the b vector of each row
        b[i] = b[k] - b[i] * factor
print(a)
print(b)


x[n - 1] = b[n - 1] / a[n - 1, n - 1]
for i in range(n - 2, -1, -1):
    sum_ax = 0

    for j in range(i + 1, n):
        sum_ax += a[i, j] * x[j]

    x[i] = (b[i] - sum_ax) / a[i, i]

print("The solution of the system is:")
print(x)

Output:

Solution by NumPy:
[ 0.38947368  0.49473684 -0.10877193  0.12982456  0.83157895]
[[ 5.00000000e+00  1.00000000e+00  0.00000000e+00  3.00000000e+00  -1.00000000e+00]
 [ 0.00000000e+00  3.00000000e+00  4.00000000e+00  1.00000000e+00  7.00000000e+00]
 [ 0.00000000e+00  0.00000000e+00  4.50000000e+00  0.00000000e+00  6.00000000e+00]
 [ 0.00000000e+00  4.44089210e-16  0.00000000e+00  3.52702703e+00  2.95945946e+00]
 [ 0.00000000e+00  0.00000000e+00  0.00000000e+00  0.00000000e+00  2.11354647e+00]]
[2.         7.         4.5        2.91891892 1.75758075]
The solution of the system is:
[ 0.38947368  0.49473684 -0.10877193  0.12982456  0.83157895]

In the code above, we used the array function and the fabs function provided by the NumPy library to create a matrix and read absolute values. We have also used Linalg; a NumPy sublibrary used to perform operations such as calculating eigenvalues and vectors and determinants.

from numpy import array, zeros, fabs, linalg

The variable a represents the constant coefficients of the variables, while b represents the matrix of constant terms corresponding to the equations. Generally, the matrix and the vector are drawn from the matrix representation of four linear equations that we intend to solve.

a = array(
    [
        [0, 6, -1, 2, 2],
        [0, 3, 4, 1, 7],
        [5, 1, 0, 3, -1],
        [3, 1, 3, 0, 2],
        [4, 4, 1, -2, 1],
    ],
    float,
)
# the b matrix constant terms of the equations
b = array([5, 7, 2, 3, 4], float)

The goal of elimination is to transform the matrix to have zeros in the leading diagonal.

When performing Gaussian elimination, the first row should remain fixed. However, if the pivot element in the first row is equal to zero, we need to swap the first row with any subsequent row with no zero as the leading element.

We have used the k loop to specify the fixed row in the code above. In this case, k represents the index of the fixed row which ranges from 1 to n-1.

The subsequent if statements checks if the pivot element is zero and swap it with the appropriate row. The i used represents the index of the rows below the fixed row.

for k in range(n - 1):
    if fabs(a[k, k]) < 1.0e-12:

        for i in range(k + 1, n):
            if fabs(a[i, k]) > fabs(a[k, k]):
                a[[k, i]] = a[[i, k]]
                b[[k, i]] = b[[i, k]]
                break

In the code below, we use a loop to ensure that elimination is not applied in the first row.

We have then multiplied each row below the fixed row by a factor. A factor equals the diagonal element corresponding to the fixed row divided by the first element in the corresponding row.

Next, we subtract the two rows from the fixed row to get the element in the second column under the main diagonal. The subsequent rows can also be reduced in the same procedure.

As we have mentioned, the first row will always remain unchanged. We have also calculated the factor and the b vector for each row as shown below.

for i in range(k + 1, n):
    # first row should remain intact
    if a[i, k] == 0:
        continue
    # calculating the factor
    factor = a[k, k] / a[i, k]
    for j in range(k, n):
        a[i, j] = a[k, j] - a[i, j] * factor
        # calculating the b vector
    b[i] = b[k] - b[i] * factor
print(a)
print(b)

The term back substitution refers to the fact that the values of x can be calculated from the last equation to the first equation. In this case, we leverage the corresponding values of the b scalar vector multiplied by the coefficients of x. The result is then divided by the corresponding element in the main diagonal.

We have prioritized the calculation of values of x starting from n and the remaining values of x from n-1 to the first one. On the other hand, summation should begin from j+1 since only the values of x are known to us from the previous calculations. In addition, the value of the sum is a result of the values of x multiplied by the values of a. We can now calculate the value of x since we have the value of the sum.

# calculating the value of x beginning from n-1
x[n - 1] = b[n - 1] / a[n - 1, n - 1]
for i in range(n - 2, -1, -1):
    sum_ax = 0
    # calculating the value of the sum
    for j in range(i + 1, n):
        sum_ax += a[i, j] * x[j]
    # calculating the value of x
    x[i] = (b[i] - sum_ax) / a[i, i]

print(x)

Output:

[[ 2.00000000e+00  3.00000000e+00  4.00000000e+00  1.00000000e+00
   7.00000000e+00]
 [ 0.00000000e+00  7.00000000e+00 -1.00000000e+00  3.00000000e+00
   1.00000000e+00]
 [ 0.00000000e+00  0.00000000e+00 -1.30000000e+01  2.00000000e+00
  -2.10000000e+01]
 [ 0.00000000e+00  0.00000000e+00  0.00000000e+00  2.81250000e+00
   5.81250000e+00]
 [ 0.00000000e+00  8.88178420e-16  0.00000000e+00 -4.44089210e-16
   4.95734797e+00]]
[  7.           5.         -14.           0.625        0.15371622]
The solution of the system is:
[0.02170543 0.79224806 1.05116279 0.15813953 0.03100775]

The following NumPy function should also result in the same answer.

print("The numpy solution:")
print(linalg.solve(a, b))

Output:

Solution by Nuumpy:
[0.02170543 0.79224806 1.05116279 0.15813953 0.03100775]