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): 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)]) print("===================") RREF_classic(A, verbose=True) print("===================") RREF_lcm(B, verbose=True) print("===================") RREF_euclidean(C, verbose=True) print("===================")