Python Partial pivoting in getting the reduced row echelon form of a matrix

[psie]

TL;DR Summary

I've been trying to write a small code that takes as input an augmented matrix and outputs the reduced row echelon form of the matrix. A classic exercise probably. It's important to me that the algorithm follows Gauss elimination method, and not the Gauss-Jordan method, although these are very similar. I'm not sure the code below is true to the Gaussian elimination method, maybe it isn't? Anyway, I have a question about another thing that's part of the algorithm; partial pivoting.

I assume the reader is familiar with row echelon form of a matrix and reduced row echelon form (rref) of a matrix, the latter which is unique for each matrix. For instance, the augmented matrix $$\begin{pmatrix} 1&0&1&1\\ 1&1&-1&2\\ 2&0&1&0 \\ 0&-1&1&-3\end{pmatrix}$$has rref $$\begin{pmatrix}1&0&0&1\\ 0&1&0&5\\ 0&0&1&2\\ 0&0&0&0\end{pmatrix}.$$By augmented matrix I mean that the last column of the matrix is actually not part of the rest of the matrix (the last column is the constant vector $b$ and the rest of the matrix the so-called coefficient matrix $A$ in the system $Ax=b$), but everything gets row reduced by the same operation, so hence it's convenient to put it all into one matrix.

For the distinction between Gauss elimination and the Gauss-Jordan method, I have found this website to provide a good overview. In my book for instance, they say that Gauss elimination consists of two parts (where I myself have included the partial pivoting step in part 1):

Part $1$ (forward pass):

• In the leftmost nonzero column of an $m\times n$ augmented matrix, identify in absolute value the largest row entry and then, if necessary, swap this row with the $1$st row. This is called partial pivoting and apparently avoids floating-point precision error. The row entry that is largest in absolute value and now lies in the $1$st row is called a pivot. Make the pivot equal $1$ and then row reduce all entries below the pivot to $0$.
• Then consider row $2$ and create a $1$ in the leftmost possible column by partial pivoting, but without using previous row(s). Row reduce the entries below the pivot to $0$.
• Continue with the next row and create a $1$ in the leftmost possible column by partial pivoting, but without using previous row(s). Terminate when the column or row of the previous pivot equals $\min(m,n-1)$.

Part $2$ (backward pass):

• Starting from bottom to top, find the column of the first nonzero entry in each row (this nonzero entry will be a pivot and equal $1$). Then row reduce the entries above the pivot to $0$.
• Repeat the previous step for each row and stop when you reach the first row.

On the other hand, in Gauss-Jordan elimination one transforms the (augmented) matrix into reduced row echelon form by using the first nonzero entry in each row to make zero all other entries in its column.

That's how my book describes the difference between the two methods. It's not very clear to me what the difference is actually.

Anyway, below follows an implementation of transforming an augmented matrix into rref that I've been working on, using what I think is Gauss-Jordan's method. First, part $0$ and part $1$ of the code below is a modification of some code I found online. I then added part $2$, i.e. getting the augmented matrix from row echelon form (with normalized pivots, i.e. pivots equal to $1$) to reduced row echelon form:

[CODE lang="python" title="rref of augmented matrix"]
# reduced row echelon form (rref) with partial pivot (using Gauss-Jordan's method?)

import numpy as np
nearzero = 1.0e-10

# part 0; ask for input
M = int(input("Enter the no. of rows of augmented matrix: "))
N = int(input("Enter the no. of columns of augmented matrix: "))
A = np.zeros((M,N))
print("Enter the rows of the augmented matrix; {} numbers, separated by comma (e.g. 1,2): ".format(N))
for i in range(M):
print("Type in row {}: ".format(i+1))
l=[float(x) for x in input().split(",")]
A = np.array(l)
print("The augmented matrix is: \n", A)

# part 1; forward pass (put the matrix on row echelon form with normalized pivots)
pivotRow = 0
for column in range(N-1):
if pivotRow == min(N-1, M):
break

# partial pivot
bestRow = pivotRow
for row in range(pivotRow + 1, M):
if (abs(A[row, column]) > abs(A[bestRow, column])):
bestRow = row
if bestRow != pivotRow:
A[[pivotRow, bestRow], column:] = A[[bestRow, pivotRow], column:]

# elimination (provided non-zero pivot)
if abs(A[pivotRow, column]) > nearzero:
A[pivotRow, column:] = A[pivotRow, column:] / A[pivotRow, column]
for row in range(pivotRow + 1, M):
A[row, column:] -= A[row, column] * A[pivotRow, column:]
# avoiding round-off errors
A[row, column] = 0.0
pivotRow += 1
else:
A[pivotRow,column] = 0.0

print("The row echelon form of the augmented matrix is: \n",A)

# part 2; backward pass (put the matrix on reduced row echelon form)
for i in reversed(range(M)):
for j in range(N-1):
if abs(A[i, j]) > nearzero:
for row in reversed(range(i)):
if abs(A[row, j]) > nearzero:
A[row, j:] -= A[row, j] * A[i, j:]
break

print("The rref of the augmented matrix is: \n ", A)
[/CODE]

Here's the output of the code when I input the augmented matrix above.
[CODE title="output"]
Enter the no. of rows of augmented matrix: 4
Enter the no. of columns of augmented matrix: 4
Enter the rows of the augmented matrix; 4 numbers, separated by comma (e.g. 1,2):
Type in row 1:
1,0,1,1
Type in row 2:
1,1,-1,2
Type in row 3:
2,0,1,0
Type in row 4:
0,-1,1,-3
The augmented matrix is:
[[ 1. 0. 1. 1.]
[ 1. 1. -1. 2.]
[ 2. 0. 1. 0.]
[ 0. -1. 1. -3.]]
The row echelon form of the augmented matrix is:
[[ 1. 0. 0.5 0. ]
[ 0. 1. -1.5 2. ]
[ 0. 0. 1. 2. ]
[ 0. 0. 0. 0. ]]
The rref of the augmented matrix is:
[[ 1. 0. 0. -1.]
[ 0. 1. 0. 5.]
[ 0. 0. 1. 2.]
[ 0. 0. 0. 0.]]
[/CODE]

My question is; in the forward pass, one is swapping rows to partly avoid floating-point precision error. Shouldn't this be done in the backward pass too then? For instance, if you consider the row echelon form of the matrix in the output above, i.e. the third column, we have that $-1.5$ is greater in absolute value than $1.$. Swapping their rows would not be a good idea though. Maybe I don't fully understand what partial pivoting is all about?

 

[psie]

Maybe it isn't clear what the distinction between Gaussian elimination and Gauss-Jordan's method is in the literature. And I guess a partial pivoting in the backward pass just makes no sense (maybe there isn't even a backward pass depending on which method you choose to implement).

 

[pbuk]

I'm not sure the code below is true to the Gaussian elimination method, maybe it isn't?

No, it is Gauss-Jordan but we will come back to this later.

Anyway, I have a question about another thing that's part of the algorithm; partial pivoting.
...
My question is; in the forward pass, one is swapping rows to partly avoid floating-point precision error. Shouldn't this be done in the backward pass too then? For instance, if you consider the row echelon form of the matrix in the output above, i.e. the third column, we have that $-1.5$ is greater in absolute value than $1.$.

In what you have called the "backward pass", you are trasforming the matrix which is in row echelon form to reduced row ecehlon form. If you swap rows 2 and 3 what will be the value in column 2 of row 3? Will the matrix still be in row echelon form?

For instance, the augmented matrix $$\begin{pmatrix} 1&0&1&1\\ 1&1&-1&2\\ 2&0&1&0 \\ 0&-1&1&-3\end{pmatrix}$$has rref $$\begin{pmatrix}1&0&0&1\\ 0&1&0&5\\ 0&0&1&2\\ 0&0&0&0\end{pmatrix}.$$

That is not a good example to work with because it represents a system of 4 equations in 3 unknowns.

In my book for instance, they say that Gauss elimination consists of two parts

I think that is a confusing statement. Gaussian elimination refers to transformation of a matrix to row echelon form by elementary row operatons. To obtain the solution of a system of linear equations involves a second stage: back-substitution.

In my book for instance, they say that Gauss elimination consists of two parts (where I myself have included the partial pivoting step in part 1):

Partial pivoting is not a separate step, it is simply a strategy for deciding which row operations to perform.

This is called partial pivoting and apparently avoids floating-point precision error.

Partial pivoting does not avoid floating point precision error (which we normally call 'round-off error'), only reduce it. You say 'apparently' - it is essential that you understand how it does this and you should go back and revise this topic, either in your book, the notes I linked here before, or perhaps https://en.wikipedia.org/wiki/Round-off_error

Part $1$ (forward pass):

Part 1 is Gaussian elimination (with partial pivoting and normalization).

Part $2$ (backward pass):

• Starting from bottom to top, find the column of the first nonzero entry in each row (this nonzero entry will be a pivot and equal $1$). Then row reduce the entries above the pivot to $0$.
• Repeat the previous step for each row and stop when you reach the first row.

Part 2 is Gauss-Jordan.