class Matrix:
#### Linear Algebra Methods ####
def scale(self,i,r,**arguments):
if "verbose" in arguments and arguments["verbose"]:
if abs(r)>1:
print("Scaling row " + str(i+1) + " by " +str (r) )
elif abs(r)<1:
print("Scaling row " + str(i+1) + " by 1/"+ str(1/r))
elif r == -1:
print("Negating row " + str(i+1) )
elif r == 1:
print("Leaving row " + str(i+1) + " unchanged")
for j in range(self.n):
self.Array[i][j] = r * self.Array[i][j]
def interchange(self,i1,i2,**arguments):
if "verbose" in arguments and arguments["verbose"]:
print("Interchanging rows " + str( i1 +1) + " and " + str( i2+1))
for j in range(self.n):
temp = self.Array[i2][j]
self.Array[i2][j] = self.Array[i1][j]
self.Array[i1][j] = temp
def replacement(self,i1, i2, r,**arguments):
if "verbose" in arguments and arguments["verbose"]:
print("Replacing row " + str(i1+1) + " with " +str ( r) + " times row " + str( i2 +1) )
for j in range(self.n):
self.Array[i1][j] = self.Array[i1][j] + r*self.Array[i2][j]
#### Support Methods ####
def __init__(self,Array):
self.Array = Array
self.m = len(Array)
if self.m == 0:
print("ERROR: Empty array")
self.n = len(Array[0])
for i in range(len(Array)):
if len(Array[i]) != self.n:
print("ERROR: inconsistent array dimension")
def elt(self,i,j):
if (i<0) or (i>self.m) or (j<0) or (j>self.n):
print("ERROR: invalid dimension", i, " vs ,", j, "in")
self.print()
return self.Array[i][j]
def row_lnze_index(self,i):
''' Returns index of row LNZE if it exists.
Otherwise, returns -1'''
for j in range(self.n):
if self.elt(i,j) != 0: return j
return -1
def print(self):
# print("[",end="")
# for i in range(self.m):
# print(self.Array[i],end="")
# if i != self.m-1: print(",")
# print("]",end="")
for i in range(self.m):
print(self.Array[i])
print(".")
# def print(self):
# print("[",end="")
# for i in range(self.m):
# print(self.Array[i],end="")
# if i != self.m-1: print(",")
# print("]",end="")
def print_latex(self):
print(" \\begin{bmatrix}")
for i in range(self.m):
delim = " & "
print( delim.join(map(str,self.Array[i])) + " \\\\" )
# for j in range(self.n):
# print(self.Array[i][j],end=" ")
# if j != self.n-1:
# print("& ",end="")
# elif i!= self.m-1:
# print( "\\\\ \n",end="" )
print("\n\\end{bmatrix} \n")
import random
def random_matrix_of_rank(m,n,rank):
if (1>m) or (1>n):
print("Error: invalid dimensions", m, n)
return Matrix([[]])
if (rank>n) or (rank>m):
print("Error: rank cannot exceed m and n", m, n)
return Matrix([[]])
elif 1>rank:
return Matrix([ [0 for i in range(n) ] for j in range(m) ] )
## Pick the columns to contain pivots
## Always put a pivot in the first column
if rank==1:
pivot_cols=[]
else:
pivot_cols = random.sample( [i for i in range(1,n)] ,rank-1 )
pivot_cols.append(0)
pivot_cols.sort()
## Intialize an empty array
values = [0, 0, 0, 0, 0, 0, 1, 1, 1, -1,-1, 2,-2,3,-3]
Array = [ [ 0 for i in range(n) ] for j in range(m) ]
for row,col in enumerate(pivot_cols):
if (row>m) or (col>n):
print("Error")
return Matrix([[]])
Array[row][col] = 1
for row , lnze in enumerate(pivot_cols):
# lnze = pivot_cols[row]
for col in range(lnze+1,n):
if not(col in pivot_cols):
Array[row][col] = random.choice(values)
M = Matrix(Array)
num_ops = 2*(n+m+2)
row_operations = random.choices(["replacement","interchange","scaling"],k=num_ops)
for operation in row_operations:
if operation == "replacement":
r = random.choice([-1,1])
i1,i2 = random.sample([i for i in range(m)],k=2)
M.replacement(i1,i2,r)
elif operation == "interchange":
i1,i2 = random.sample([i for i in range(m)],k=2)
M.interchange(i1,i2)
elif operation == "scaling":
i= random.choice([i for i in range(m)])
r = random.choice([-2,-1,1,2])
M.scale(i,r)
return M
def RREF_classic( A: Matrix , **arguments):
''' Row reduces the matrix using classical Gaussian elimination'''
if "verbose" in arguments and arguments["verbose"]:
print("Reducing matrix using classical Gaussian elimination")
A.print()
for i in range(A.m):
# Flip rows as needed to ensure we have the leftmost LNZE
for i2 in range(i+1,A.m):
if ( -1 < A.row_lnze_index(i2) < A.row_lnze_index(i) ) or ( A.row_lnze_index(i) == -1 ) :
A.interchange(i,i2,**arguments)
A.print()
# Store column of current LNZE
j = A.row_lnze_index(i)
if j>-1: # If there is a nonzero pivot
# Convert the LNZE to 1
A.scale(i,1/A.elt(i,j),**arguments)
A.print()
# Use LNZE to clear out all other rows
for i2 in range(A.m):
if i2 != i and A.elt(i2,j)!=0 :
r = - A.elt(i2,j)/A.elt(i,j)
A.replacement(i2,i,r,**arguments)
A.print()
def gcd(a, b):
'''Calculate the Greatest Common Divisor of a and b.
FROM https://stackoverflow.com/questions/11175131/code-for-greatest-common-divisor-in-python#:~:text=The%20greatest%20common%20divisor%20(GCD)%20of%20a%20and%20b%20is
Unless b==0, the result will have the same sign as b (so that when
b is divided by it, the result comes out positive).
'''
while b:
a, b = b, a%b
return a
def RREF_lcm( A: Matrix , **arguments):
''' Uses Gaussian elimination modified to avoid division
by scaling the row to be replaced so it is a multiple
of the row being used.
'''
if "verbose" in arguments and arguments["verbose"]:
print("Reducing matrix using GCM's")
A.print()
for i in range(A.m):
# Flip rows as needed to ensure we have the leftmost LNZE
for i2 in range(i+1,A.m):
if ( -1 < A.row_lnze_index(i2) < A.row_lnze_index(i) ) or ( A.row_lnze_index(i) == -1 ) :
A.interchange(i,i2,**arguments)
A.print()
# Store column of current LNZE
j = A.row_lnze_index(i)
if j >-1 :# If there is a nonzero pivot
# Scale row to be replaced to be a multiple of the row used.
for i2 in range(A.m):
if i2 != i and A.elt(i2,j)!=0 :
target_lnze = A.elt(i2,j)
A.scale(i2, A.elt(i,j)/gcd(A.elt(i,j),target_lnze),**arguments)
A.print()
A.replacement(i2, i, -target_lnze/gcd(A.elt(i,j),target_lnze), **arguments)
A.print()
for i in range(A.m):
j = A.row_lnze_index(i)
if j>-1 :
if A.elt(i,j)!=0 :
A.scale(i,1/(A.elt(i,j)),**arguments)
A.print()
def RREF_euclidean( A: Matrix , **arguments):
''' Row reduces the matrix using the Euclidean algorithm to get Echelon Form.'''
if "verbose" in arguments and arguments["verbose"]:
print("Reducing matrix using Euclidean algorithm:")
print("At each stage")
print("1. Sort columns by pivot index and size")
print("2. Use the smallest magnitude pivot to reduce its position in all other rows, but without introducing fractions")
print("3. Repeat until in Echelon form.")
print("4. Scale all rows to get 1, then clear up")
print(".")
A.print()
current_row = 0 # Start with row zero
done_with_row = False
while current_row < A.m:
# Rearrange rows by LNZE
for i in range(A.m):
# Flip rows as needed to ensure we have the leftmost LNZE
for i2 in range(i+1,A.m):
if ( -1 < A.row_lnze_index(i2) < A.row_lnze_index(i) ) :
A.interchange(i,i2,**arguments)
A.print()
elif ( (A.row_lnze_index(i) == -1) and (A.row_lnze_index(i2)>-1) ):
A.interchange(i,i2,**arguments)
A.print()
# Sort rows by ABS VALUE of the LNZE itself
for i in range(A.m):
for i2 in range(i+1,A.m):
if (A.row_lnze_index(i) != -1 ) and (A.row_lnze_index(i) == A.row_lnze_index(i2)):
if ( abs(A.elt(i,A.row_lnze_index(i))) > abs(A.elt(i2,A.row_lnze_index(i)))):
A.interchange(i,i2,**arguments)
A.print()
# Store column of current LNZE
j = A.row_lnze_index(current_row)
if j == -1:
## Since we have sorted the rows, if the current row has no LNZE,
## we are at the bottom of the matrix
break
elif j > -1:
# Use the current nonzero pivot to reduce (and maybe zero out)
# its value in later rows
done_with_row = True
for i2 in range(current_row+1,A.m):
if i2 != current_row and A.elt(i2,j)!=0 :
remainder = A.elt(i2,j) % A.elt(current_row,j)
if( remainder ) == 0:
# print("divisible")
r = -A.elt(i2,j)/A.elt(current_row,j)
else:
done_with_row = False
r = -(A.elt(i2,j)-remainder)/A.elt(current_row,j)
A.replacement(i2,current_row,r,**arguments)
A.print()
if done_with_row :
# print("Rescaling row")
A.scale( current_row, 1/A.elt(current_row,j) , **arguments)
A.print()
current_row = current_row +1
print("Matrix is in Echelon form")
print("Now finding Reduced Echelon form\n")
for current_row in range(A.m):
j = A.row_lnze_index(current_row)
if j>-1:
for i2 in range(current_row):
if A.elt(i2,j) != 0:
A.replacement(i2,current_row, -A.elt(i2,j),**arguments)
A.print()
print("Matrix is now in Reduced Echelon Form")
m = random.randint(2,4)
n = random.randint(m+1,8)
rank = min( random.randint(int(n/2),n), m, n)
A = random_matrix_of_rank(m,n,rank)
B = Matrix([[A.Array[i][j] for j in range(n)] for i in range(m)])
C = Matrix([[A.Array[i][j] for j in range(n)] for i in range(m)])
A.print_latex()
print("===================")
RREF_classic(A, verbose=True)
print("===================")
RREF_lcm(B, verbose=True)
print("===================")
RREF_euclidean(C, verbose=True)
print("===================")
A.print_latex()
B.print_latex()
C.print_latex()