python built-in function to do matrix reduction

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

6 Answers

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]) 

see

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.

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.

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:

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 xrange(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: #pivot index pid = argmax( abs(A[i:,pcol]) ) #Row exchange A[i,:],A[pid+i,:] = A[pid+i,:].copy(),A[i,:].copy() pcol += 1 if pcol >= n : break pivot = float(A[i,pcol]) for j in xrange(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 some simple tests

print "/*--------------------------------------/" print "/ Simple TEST /" print "/--------------------------------------*/" A = array([[1,2,3],[4,5,6],[-7,8,9]]) print "A:\n",R R = rref(A,precision=6) 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) /*--------------------------------------/ / Simple TEST / /--------------------------------------*/ A: [[ 1. 0. 1.] [-0. 1. 1.] [ 0. 0. 0.] [ 0. 0. 0.]] 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.]] print "/*--------------------------------------/" print "/ Not Invertable TEST /" print "/--------------------------------------*/" A = array([ [2,2,4, 4], [3,1,6, 2], [5,3,10,6]]) print "A:\n",A R = rref(A,precision=2) print "R:\n",R print print "A^{T}:\n",A.T R = rref(A.T,precision=10) print "R:\n",R /*--------------------------------------/ / 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.]] 

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 

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)