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, 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 later rows, but without introducing fractions") print("3. Once a pivot has cleared out later rows, scale it to be 1") print("4. Once the matrix is in echelon form, 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,3) n = random.randint(m+1,m+2) 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("===================")