python built-in function to do matrix reduction

Question

Does python have a built-in function that converts a matrix into row echelon form (also known as upper triangular)?

Wat

Answer 1

If you can use sympy, Matrix.rref() can do it:

In [8]: sympy.Matrix(np.random.random((4,4))).rref()
Out[8]: 
([1, 1.42711055402454e-17, 0, -1.38777878078145e-17]
[0,                  1.0, 0,  2.22044604925031e-16]
[0, -2.3388341405089e-16, 1, -2.22044604925031e-16]
[0, 3.65674099486992e-17, 0,                   1.0],
 [0, 1, 2, 3])

Fortney

How to interpret the result? – Barneybarnhart

The 4x4 matrix is your matrix in RRE (watch out for floating point precision), and the 1x4 matrix lists the indices of your pivot vars. – Arcature

by hand, ~5-10 minutes. By sympy, 1.13 ms :-3 – Effluent

Thanks for the useful answer. It is not that straightforward to convert back to a numpy array, though. I used np.array(sympy.lambdify((),result[0])()). Maybe there is a better way? – Inane

Never mind, np.array(result[0].tolist(), dtype=float) is simpler. – Inane

Answer 2

Yes. In scipy.linalg, lu does LU decomposition which will essentially get you row-echelon form.

There are other factorizations such as qr, rq, svd, and more, if you're interested.

Documentation.

Biebel

LU decomposition is not the same as row-echelon form. See math.stackexchange.com/a/1614952 – Glossography

Yes, it is possible to compute the reduced row echelon from these functions but why make the user jump through the hoops. Most sensible applications such as scilab, Matlab, Mathematica etc, have these built-in. Instead, Python uses must go through the symbolic algebra library. – Cob

"In scipy.linalg, lu does LU decomposition which will essentially get you row-echelon form": Perhaps you can provide an example showing how to use the p, l, u matrices returned by lu to obtain the RRE form. – Goglet

Answer 3

see http://mail.scipy.org/pipermail/numpy-discussion/2008-November/038705.html

Basically: don't do it.

The rref algorithm produces too much inaccuracy when implemented on a computer. So you either want to solve the problem in another way, or use symbolics like @aix suggested.

Eustasius

Just wrap your numbers in a rational datatype such as python's fraction.Fraction, which is already overloaded for the arithmetic operations. You can even vectorize the constructor with numpy, ie VecFraction = np.vectorize(fraction.Fraction) and rational_array = VecFraction(numeric_array). There's no reason a computer could not perform exact RREF. – Pubis

The link in the answer doesn't work for me, numpy-discussion.10968.n7.nabble.com/… – Ovida

"or use symbolics like @aix suggested", where's this written? – Lett

If rref is inaccurate on a computer, then why do we have the similar algorithm LU implemented on scipy for production use? – Idolist

Answer 4

I agree with a @Mile comment to @WinstonEwert answer There's no reason a computer could not perform RREF with given precision.

The realization of RREF shouldn't be very complicated, and matlab somehow managed to have this function, so numpy should have too.

I did a very simple and straightforward realization, which is very inefficient. But for simple matrices it works pretty well:

UPDATE:

Seems, @JustMe spotted a bug in this rref realization, which appeared if an input matrix has the first column of zeros. So here an updated version.

from numpy import *

def rref(mat,precision=0,GJ=False):
    m,n = mat.shape
    p,t = precision, 1e-1**precision
    A = around(mat.astype(float).copy(),decimals=p )
    if GJ:
        A = hstack((A,identity(n)))
    pcol = -1 #pivot colum
    for i in range(m):
        pcol += 1
        if pcol >= n : break
        #pivot index
        pid = argmax( abs(A[i:,pcol]) )
        #Row exchange
        A[i,:],A[pid+i,:] = A[pid+i,:].copy(),A[i,:].copy()
        #pivot with given precision
        while pcol < n and abs(A[i,pcol]) < t:
            pcol += 1
            if pcol >= n : break
            #pivot index
            pid = argmax( abs(A[i:,pcol]) )
            #Row exchange
            A[i,:],A[pid+i,:] = A[pid+i,:].copy(),A[i,:].copy()
        if pcol >= n : break
        pivot = float(A[i,pcol])
        for j in range(m):
            if j == i: continue
            mul = float(A[j,pcol])/pivot
            A[j,:] = around(A[j,:] - A[i,:]*mul,decimals=p)
        A[i,:] /= pivot
        A[i,:] = around(A[i,:],decimals=p)
        
    if GJ:
        return A[:,:n].copy(),A[:,n:].copy()
    else:
        return A   

Here are some simple tests

print("/*--------------------------------------/")
print("/             Simple TEST               /")
print("/--------------------------------------*/")
A = array([[1,2,3],[4,5,6],[-7,8,9]])
R = rref(A,precision=6)
print("A:\n",A)
print("R:\n",R)
print()
print("With GJ ")
R,E =   rref(A,precision=6,GJ=True)
print("R:\n",R)
print("E:\n",E)
print("AdotE:\n",around( dot(A,E), decimals=0))
print()
A = array([[0, 0, 1], [0, 1, 0]]) 
R = rref(A, precision=1)
print("A:\n",A)
print("R:\n",R)
print()
A = array([[1,2,2,2],[2,4,6,8],[3,6,8,10]])
R = rref(A,precision=6)
print("A:\n",A)
print("R:\n",around(R, decimals=0))
print()

print("/*--------------------------------------/")
print( "/        Not Invertable TEST            /")
print( "/--------------------------------------*/")
A = array([
    [2,2,4, 4],
    [3,1,6, 2],
    [5,3,10,6]])
R = rref(A,precision=2)
print("A:\n",A)
print("R:\n",R)
print()
print("A^{T}:\n",A.T)
R = rref(A.T,precision=10)
print("R:\n",R)

/*--------------------------------------/
/             Simple TEST               /
/--------------------------------------*/
A:
 [[ 1  2  3]
  [ 4  5  6]
  [-7  8  9]]
R:
 [[1. 0. 0.]
  [0. 1. 0.]
  [0. 0. 1.]]

With GJ 
R:
 [[1. 0. 0.]
  [0. 1. 0.]
  [0. 0. 1.]]
E:
 [[-0.071428  0.142857 -0.071429]
  [-1.857142  0.714285  0.142857]
  [ 1.595237 -0.523809 -0.071428]]
AdotE:
 [[ 1.  0.  0.]
  [ 0.  1.  0.]
  [-0.  0.  1.]]

A:
 [[0 0 1]
  [0 1 0]]
R:
 [[0. 1. 0.]
  [0. 0. 1.]]

A:
 [[ 1  2  2  2]
  [ 2  4  6  8]
  [ 3  6  8 10]]
R:
 [[ 1.  2.  0. -2.]
  [ 0.  0.  1.  2.]
  [ 0.  0.  0.  0.]]

/*--------------------------------------/
/        Not Invertable TEST            /
/--------------------------------------*/
A:
 [[ 2  2  4  4]
  [ 3  1  6  2]
  [ 5  3 10  6]]
R:
 [[ 1.  0.  2.  0.]
  [-0.  1. -0.  2.]
  [ 0.  0.  0.  0.]]

A^{T}:
 [[ 2  3  5]
  [ 2  1  3]
  [ 4  6 10]
  [ 4  2  6]]
R:
 [[ 1.  0.  1.]
  [-0.  1.  1.]
  [ 0.  0.  0.]
  [ 0.  0.  0.]]

Lafayette

the line pcol += 1 should be the first thing we do in the while ... loop for this implementation to work correctly. – Foreshore

@joni i cannot edit the answer ("too many pending edits"). If you agree with the correction, could you please fix it for future users' benefit? – Foreshore

I'm not. @JustMe could you please explane why? – Lafayette

it appears the intention was that when you exit the loop, i, pcol should point to the pivot element - the largest element in the lower portion of column pcol, brought up to the current row i. Then we should apply row operations to the same column to clear the elements below it. But as it is, we advance pcol to the next column. If this isn't convincing, A = array([[0,0, 1], [0, 1, 0]]) print(rref(A, precision=5)) gives [[0. 0. 1.] [0. 1. 0.]] which is not in row echelon form. – Foreshore

That said we should also add a check if pcol >= n : break inside the while loop after we advance pcol. – Foreshore

Answer 5

You can define it yourself:

def rref(matrix):
    A = np.array(matrix, dtype=np.float64)

    i = 0 # row
    j = 0 # column
    while True:
        # find next nonzero column
        while all(A.T[j] == 0.0):
            j += 1
            # if reached the end, break
            if j == len(A[0]) - 1 : break
        # if a_ij == 0 find first row i_>=i with a 
        # nonzero entry in column j and swap rows i and i_
        if A[i][j] == 0:
            i_ = i
            while A[i_][j] == 0:
                i_ += 1
                # if reached the end, break
                if i_ == len(A) - 1 : break
            A[[i, i_]] = A[[i_, i]]
        # divide ith row a_ij to make it a_ij == 1
        A[i] = A[i] / A[i][j]
        # eliminate all other entries in the jth column by subtracting
        # multiples of of the ith row from the others
        for i_ in range(len(A)):
            if i_ != i:
                A[i_] = A[i_] - A[i] * A[i_][j] / A[i][j]
        # if reached the end, break
        if (i == len(A) - 1) or (j == len(A[0]) - 1): break
        # otherwise, we continue
        i += 1
        j += 1

    return A

Ernestoernestus

There are a couple of mistakes with this, but I still find the process elegant: 1. if j == len(A[0]) - 1 should break out of the external while 2. you should check whether A[i][j] is not zero (and ideally not one), before normalizing 3. you should only multiply other rows if A[i][j] is not zero – Leif

Answer 6

Here's a working version that's pretty much just a numpy version of MATLAB's rref function:

def rref(A, tol=1.0e-12):
    m, n = A.shape
    i, j = 0, 0
    jb = []

    while i < m and j < n:
        # Find value and index of largest element in the remainder of column j
        k = np.argmax(np.abs(A[i:m, j])) + i
        p = np.abs(A[k, j])
        if p <= tol:
            # The column is negligible, zero it out
            A[i:m, j] = 0.0
            j += 1
        else:
            # Remember the column index
            jb.append(j)
            if i != k:
                # Swap the i-th and k-th rows
                A[[i, k], j:n] = A[[k, i], j:n]
            # Divide the pivot row i by the pivot element A[i, j]
            A[i, j:n] = A[i, j:n] / A[i, j]
            # Subtract multiples of the pivot row from all the other rows
            for k in range(m):
                if k != i:
                    A[k, j:n] -= A[k, j] * A[i, j:n]
            i += 1
            j += 1
    # Finished
    return A, jb

Example:

A = np.array([[16.0, 2, 3, 13], [5, 11, 10, 8],
              [9, 7, 6, 12], [4, 14, 15, 1]])
Areduced, jb = rref(A)
print(f"The matrix as rank {len(jb)}")
print(Areduced)
    

Punt

I tried this version, but it has some weird behavior with floating point given your example: [[ 1. 0. 0. 1.] [ 0. 1. 0. 3.] [ 0. 0. 1. -3.] [ 0. 0. 0. 0.]] vs. int [[1 0 0 0] [0 1 0 0] [0 0 1 0] [0 0 0 1]] – Unsatisfactory

Answer 7

For any matrix (having full rank or deficiency), you can loop through the rows (or columns, this is equivalent) and see which row contributes to rank. Selecting the first n rows, you get a matrix having at most rank n.

from numpy.linalg import matrix_rank

def get_linind_rows(M):
    Mrows=[]
    Mranks=[]
    for i,Mrow in enumerate(M):
        Mrows.append(Mrow)
        Mrank = matrix_rank(Mrows) # this rank can be max i+1
        #print(i,Mrow,Mrank)
        Mranks.append(Mrank)
    Mranks_diff = np.diff(Mranks,prepend=0) # see where the rank increases by 1
    return M[Mranks_diff > 0] # select these rows as they are linear independent

Audient