Removing redundant constraints from a linear system of equations

Suppose you are solving an optimization problem with linear equality constraints of the form 

$ Ax = b $

Several off-the-shelf optimizers, including those available in CVXOPT and scipy will complain if A is not full rank.  That is, if A contains redundant constraints then optimization will likely fail.  So how can we easily deal with this issue?  The code below is something you can plug into your program for a quick fix to the issue.  Simply call A, b = reduce_row_echelon(A, b) to get an equivalent set of constraints with redundancies removed. 

# Consumes a *consistent* rectangular system of equations
# Produces a new compacted system of equations with linearly dependent rows removed
# Note: this algorithm computes the row echelon form of the augmented matrix
# This algorithm should be numerically robust
def reduce_row_echelon(B, y):
    A = np.concatenate([B, y[:,np.newaxis]], axis=1)
    m,n = A.shape
    c = 0
    for r in range(m):
        # find first non-zero column
        while c+1 < n and np.all(A[r:, c] == 0):
            c += 1
        if c+1 == n:
            return A[:r,:-1], A[:r,-1]
        # find row with largest value in that column
        i_max = r+np.argmax(np.abs(A[r:,c]))
        # swap those rows
        A[[i_max, r]] = A[[r, i_max]]
        A[r,:] /= A[r,c]
        # zero out column for every row below r
        for i in range(r+1, m):
            A[i,:] -= A[i,c]*A[r,:]
        # clamp small values to 0
        A[np.abs(A) < 1e-12] = 0
    return A[:,:-1], A[:,-1]

The above code is a quick fix, but is maybe not an ideal solution. A better solution is to carefully analyze the structure of your constraints and set it up so all constraints are linearly independent. However, that may be nontrivial in some cases, so having a quick fix like the code above might be a useful starting point. It works by transforming the linear system into reduced row echelon form, and removing zero rows. Unfortunately, there is no function available in numpy to do this to my knowledge, since the primary use case of reduced row echelon form is to solve a linear system of equations, and this approach is known to suffer from numerical stability issues. I was able to find an rref function in sympy, but that is not an ideal solution because sympy is primarily a library for symbolic computation rather than numerical computation. I wrote this code several years ago for a problem I was working on for a project, and encountered a new problem today where I needed to dig it up. Since I wasn't able to find any simple and general-purpose solution for the problem of redundant constraints in my google search, I hope this code will help others facing similar issues.