Reduced Row Echelon form script doesn't work in specific cases

Question

I'm trying to write a script that puts a matrix in reduced row echelon form, however my script (and other python scripts that I have found on the internet) seem to fail in specific cases.

For example, my script produces correct results on the following matrices:

In: [[-1  1]
     [-1  2]
     [-3  2]]

Out: [[1. 0.]
      [0. 1.]
      [0. 0.]]



In:  [[-1  1  0  1]
      [-2 -3 -1 -2]
      [-3 -1 -2 -1]]

Out: [[ 1.         0.          0.         -0.33333333]
     [ 0.          1.          0.          0.66666667]
     [-0.         -0.          1.          0.66666667]]

But it fails on matrices such as this one:

In: [[ 0 -3  1 -1]
     [-2  1  0  3]]



 Out: [[ 1.  -0.5 -0.  -1.5]
       [ 0.   0.   0.   0. ]]

Where the correct answer should be:

[[ 1. 0. -0.16666, -1.33333
 [ 0. 1. -0.33333, 0.33333 ]]

The code for the rref is the following:

import numpy as np

def rref(M):
        lead = 0
        rowCount = len(M)
        columnCount = len(M[0])
        for r in range(rowCount):
                if lead >= columnCount:
                        return
                i = r
                while (M[i][lead] == 0):
                        i += 1
                        if (i == rowCount):
                                i = r
                                lead += 1
                                if (columnCount == lead):
                                        return
                M[r],M[i] = M[i],M[r]

                if (M[r][lead] != 0):
                        M[r] = M[r]/M[r][lead]
                for i in range (rowCount):
                        if (i != r):
                                M[i] = M[i] - M[i][lead]*M[r]
                lead += 1

a = np.array([[0,-3,1,-1],[-2,1,0,3]])

#prints original matrix
print(a)

#ARRAY MUST BE FLOATS FOR RREF FUNCTION

a = a.astype(np.float)
rref(a)

#prints rref matrix
print(a)

This code is an improvement on the also non-functional code found here: https://rosettacode.org/wiki/Reduced_row_echelon_form#Python

I am aware that sympy has a built-in rref function and I know how to use it, but I need to implement my own rref code as I will later need to keep track of the row swaps/operations performed on the original matrix.

Answer 1

Your rows didn't swap correctly.

M[r],M[i] = M[i],M[r]

results in

[-2.  1.  0.  3.]
[-2.  1.  0.  3.]

The proper way to swap is

M[[r,i]] = M[[i,r]]