Partial pivoting code in python

Question

import numpy as np

def eliminateCol(A, k=0):
    # guarantee that A is float type
    if A.dtype.kind != 'f' and A.dtype.kind != 'c':  return None
    # load dimensions
    dim0,dim1 = np.shape(A)
    # initialize outputs
    swap = np.identity(dim0,dtype=float)
    oper = np.identity(dim0,dtype=float)
    # k < dim1
    k = k%dim1
    # if the column k is zero column, return
    if A[:,k].any() == False:   return swap,oper

    i0 = 0
    
    if k != 0:
        # find a row i0 whose leftmost k-1 entries are all zero
        while i0 < dim0:
            if not A[i0, 0:k].any():
                break
            i0 += 1
    
    # find a nonzero element A[i, k] to use as a pivot
    i = i0
    while i < dim0:
        if A[i, k] != 0:
            break
        i += 1
    
    if i == dim0:
        return swap, oper
    
    # Partial Pivoting: ensures that the largest element of the column is placed on the diagonal
    max_index = np.argmax(np.abs(A[i0:, k])) + i0
    A[[i0, max_index], ] = A[[max_index, i0], ]; swap[[i0, max_index], ] = swap[[max_index, i0], ]
    
    # row i0 is the pivot row and A[j, k] is zero for i0 < j <= i
    for i in range(i0 + 1, dim0):
        if A[i, k] != 0:
            quot = A[i, k] / A[i0, k]
            A[i, :] -= quot * A[i0, :]
            oper[i, i0] = -quot
    
    return swap, oper

def eliminateRow(A, k=-1):
    # guarantee that A is float type
    if A.dtype.kind != 'f' and A.dtype.kind != 'c':  return None
    # load dimensions
    dim0,dim1 = np.shape(A)
    # initialize outputs
    oper = np.identity(dim0,dtype=float)
    # k < dim0
    k = k%dim0
    # if the row k is zero column, return
    if k == 0 or A[k,:].any() == False:    return oper

    # find a nonzero element A[k,j] to use as a pivot
    for j in range(dim1):
        if A[k, j] != 0:    break
        
    for i in range(k - 1, -1, -1):
        if A[i, j] != 0:
            quot = A[i, j] / A[k, j]
            A[i, :] -= quot * A[k, :]
            oper[i, k] = -quot
    return oper  

def eliminateScale(A):
    # guarantee that A is float type
    if A.dtype.kind != 'f' and A.dtype.kind != 'c':  return None
    # load dimensions
    dim0,dim1 = np.shape(A)
    # initialize outputs
    oper = np.identity(dim0,dtype=float)
    # find pivots
    i = 0
    while i < dim0:
        j = 0
        while j < dim1:
            if A[i,j] == 0:     j+= 1; continue
            else:
                oper[i,i] = 1/A[i,j]
                A[i,j:]   = A[i,j:]/A[i,j]
                break
            j+= 1
        i+= 1
    return  oper

def get_reduced_echelon_form (A):
    
    for k in range(A.shape[1]):
        eliminateCol(A, k)
        eliminateRow(A, k)

    # Finally, scale the matrix
    eliminateScale(A)
    return A

def get_echelon_form(A):
    dim0, dim1 = np.shape(A)

    for k in range(min(dim0, dim1)):
        # Eliminate below diagonal in column k
        swap, oper = eliminateCol(A, k)

    return A

# Example Usage:
B = np.array([[0.02, 0.01, 0.0, 0.0, 0.02],
              [1.0, 2.0, 1.0, 0.0, 1.0],
              [0.0, 1.0, 2.0, 1.0, 4.0],
              [0.0, 0.0, 100.0, 200.0, 800.0]], dtype=float)
A = np.array([[1, 2, 3],
              [4, 5, 6],
              [7, 8, 9],
              [10, 11, 12]], dtype=float)

echelon = get_echelon_form(A.copy())
reduced = get_reduced_echelon_form(A.copy())

print(B)
print(echelon)
print(reduced)

This is my code to find echelon and reduced echelon form of a matrix using partial pivoting. Please help me to any bugs or shortcomings.(please ignore using for loop instead of while loop)

At the result of execution, correct reduced row echelon form is found for zero column matrix, zero row matrix, and a example usage matrix for partial pivoting. I am looking for a bug or shortcoming

Answer 1

lint

def eliminateCol(A, k=0):

You are sharing code, seeking to collaborate with the broader python community. One of the community guidelines, pep8, asks that you spell it eliminate_col. (I won't pursue A, as I understand that math conventions can also be helpful.)

It would have been a kindness to the reader to annotate A as being of type np.ndarray.

Please run your source code through black before asking others to read it.

exceptions

This is just Wrong:

    if A.dtype.kind != 'f' and A.dtype.kind != 'c':  return None

I read the signature. You made a promise. By the time eliminate_col() returns, it shall have eliminated a column. But in this case it fails to make good on that promise. That is incorrect behavior.

Also, please don't condense two lines down to one. And especially avoid abominations like this: if A[i,j] == 0: j+= 1; continue

This code should check the mandatory pre-condition, and throw when caller messed up:

    if A.dtype.kind not in ('f', 'c'):
        raise ValueError(f"unsupported dtype: {A.dtype.kind}")

Skipping down to the return statement, it's unclear how the (swap, oper) tuple satisfies your original promise to eliminate a column. You elected not to offer a """docstring""" for this function which explains what it does. In the review context you elected not to Cite Your Reference(s), so it's unclear what math theorems you are relying on here. All of which makes it much more difficult for the reader to examine this and confidently proclaim, "oh, yes, this definitely does the right thing with all input matrices."

Absent a citation, I arbitrarily choose to impose this one for purposes of reviewing the code.

All right, back to checking pre-conditions. I really do appreciate that the .shape tuple unpack strictly enforces that the ndarray passed in shall be a 2-D array, or we die with a helpful diagnostic.

extract helper

    # find a nonzero element A[i, k] to use as a pivot

It seems pretty clear here that the code is asking you to extract a helper function, so you can write

    i = find_nonzero_pivot(A, k)

div by zero

        if A[i, k] != 0:
            quot = A[i, k] / A[i0, k]

I don't understand that logic. It would be a bad thing for the quotient to be zero?!?

It seems like testing the denominator for zeroness is the more pressing concern.

DRY

I don't understand eliminateRow. It is related to eliminateCol, but different, in ways the comments fail to motivate and which the missing docstring fails to describe. You have not made it easy to trace from elements of the cited reference to elements of the implementation.

I'm willing to say that this function will return oper. More than that? It's up to you to offer a more eloquent argument for a stronger post-condition.

validation

You did not describe why sympy's rref() is unsuitable for your use case. More worryingly, you did not offer a test suite that {compares, contrasts} how this implementation and others available on pypi will handle matrices of interest. One way to instill confidence in the correctness of an implementation is to show how it behaves similarly to pre-existing implementations that folks have placed their trust in.

---

I am happy for you that the example reduced echelon form meets your needs. I would not be willing to call this code in a production setting where maintenance engineers might be called on to diagnose or repair this code.

Answer 2

There are a few points to consider and potential improvements:

Return of Elimination Functions: Your eliminateCol and eliminateRow functions return transformation matrices, but these are not used in your main functions get_echelon_form and get_reduced_echelon_form. If you're only interested in the transformed matrix and not the individual transformation matrices, this is fine. However, if you want to keep track of the transformations (which could be useful for finding inverses or for other linear algebra applications), you should modify your main functions to apply these transformation matrices.

In-Place Modifications: You're modifying the matrix A in place in your functions. This is efficient in terms of memory usage, but it means the original matrix is lost. If you need to keep the original matrix for any reason, you should work on a copy of the matrix instead.

Floating-Point Operations: When dealing with floating-point numbers, direct comparisons (like if A[i, k] != 0:) can sometimes be problematic due to precision issues. It's often safer to check if the number is "close enough" to zero, using a small tolerance value.

Efficiency: Using nested loops can be inefficient for large matrices. While optimization might not be your primary concern, be aware that this implementation might not scale well for very large matrices.