Create a Program to Create PLU Decomposition of Matrices in Python Assignment Solution

Instructions

Objective

Write a python assignment program to create PLU decomposition of matrices.

Requirements and Specifications

Source Code

## CS 132 Homework 4 A ### Due Wednesday August 4th at Midnight (1 minute after 11:59pm) in Gradescope (with grace period of 6 hours) ### Homeworks may be submitted up to 24 hours late with a 10% penalty (same grace period) Please read through the entire notebook, reading the expository material and then do the problems, following the instuctions to try various features; among the exposition of various features of Python there are three problems which will be graded. All problems are worth 10 points. You will need to complete this notebook and then convert to a PDF file in order to submit to Gradescope. Instructions are given here: https://www.cs.bu.edu/fac/snyder/cs132/HWSubmissionInstructions.html # Here are some imports which will be used in the code in the rest of the lab # Imports used for the code in CS 237 import numpy as np # arrays and functions which operate on arrays import matplotlib.pyplot as plt # normal plotting # NOTE: You may not use any other libraries than those listed here without explicit permission. # Gaussian Elimination # number of digits of precision to print out prec = 4 ######################################################################################################## # Returns the Row Echelon (upper triangular) form def forwardElimination(A,traceLevel=0): A = (np.copy(A)).astype(float) if (traceLevel > 0): print("Running Forward Elimination on:\n") print(np.round(A, decimals=4),"\n") print() (numRows,numCols) = A.shape # r = row we are currently working on pivot value is A[r][c] r = 0 for c in range(numCols): # solve for variable in column c # find row in column c with first non-zero element, and exchange with row r r1 = r while(r1 < numRows): if (not np.isclose(A[r1][c],0.0)): # A[r1][c] is first non-zero element at or below r in column c break r1 += 1 if(r1 == numRows): # all zeros below r in this column continue # go on to the next column, but still working on same row if(r1 != r): # exchange rows r1 and r A[[r1,r],:] = A[[r,r1],:] if (traceLevel == 2): print("Exchange R" + str(r1+1) + " and R" + str(r+1) + "\n") print(np.round(A, decimals=4)) print() # now use pivot A[r][c] to eliminate all vars in this column below row r for r2 in range(r+1,numRows): rowReduce(A,r,c,r2,traceLevel) r += 1 if (r >= numRows): break return A # for pivot A[r][c], eliminate variable in location A[r2][c] in row r2 using row operations def rowReduce(A,r,c,r2,traceLevel=0): if(not np.isclose(A[r2][c],0.0)): factor = -A[r2][c]/A[r][c] A[r2] += factor * A[r] if(traceLevel == 2): print("R" + str(r2+1) + " += " + str(np.around(factor,prec)) + "*R" + str(r+1) + "\n") print(np.round(A, decimals=4)) print() # Take a Row Echelon Form and return a Reduced Row Echelon Form def backwardSubstitute(A,augmented=True,traceLevel=0): numRows,numCols = A.shape if (A.dtype != 'float64'): A = A.astype(float) # now back-substitute the variables from bottom row to top if (traceLevel > 0): print("Creating Reduced Row Echelon Form...\n") for r in range(numRows): # find variable in this row for c in range(numCols): if(not np.isclose(A[r][c],0.0)): break if ((augmented and c >= numCols-1) or (c >= numCols)): # inconsistent or redundant row continue # A[r][c] is variable to eliminate factor = A[r][c] if (np.isclose(factor,0.0)): # already eliminated continue if(not np.isclose(factor,1.0)): A[r] *= 1/factor if (traceLevel == 2): print("R" + str(r+1) + " = R" + str(r+1) + "/" + str(np.around(factor,prec)) + "\n") print(np.round(A, decimals=4)) print() for r2 in range(r): rowReduce(A,r,c,r2,traceLevel) return A # try to find row of all zeros except for last column, in augmented matrix. def noSolution(A): numRows,numCols = A.shape for r in range(numRows-1,-1,-1): # start from bottom, since inconsistent rows end up there for c in range(numCols): if(not np.isclose(A[r][c],0.0)): # found first non-0 in this row if(c == numCols-1): return True else: break return False ######################################################################################################## # Runs GE and returns a reduced row echelon form # If b == None assumes that A is NOT an augmented matrix, and runs GE and returns Reduced Row Echelon Form # If b is a column matrix then adjoins it to A and runs GE; # Always returns the Reduced Row Echelon Form # If inconsistent then also prints out "Inconsistent!" # If b is a length n array instead of a 1 x n array (column vector) # b will be converted to a column vector, for convenience. # traceLevel 0 will not print anything out during the run # traceLevel 1 will print out various stages of the process and the intermediate matrices # traceLevel 2 will also print out the row operations used at each step. # If you want to produce an Echelon Form (NOT reduced), then use forwardElimination instead. # See examples for more explanation of how to use this def GaussianElimination(A,b=None, traceLevel = 0): if( type(b) != type(None)): if( (A.shape)[0] == 1 ): b = np.array([b]).T Ab = (np.hstack((A.copy(),b))).astype(float) else: Ab = A.copy().astype(float) if( traceLevel > 0 and type(b) == type(None)): print("Creating Reduced Row Echelon Form:\n") print(np.round(Ab, decimals=4)) print() elif( traceLevel > 0 ): print("Running Gaussian Elimination on augmented matrix:\n") print(np.round(Ab, decimals=4)) print() B = forwardElimination(Ab,traceLevel) if( traceLevel > 0 ): print("Echelon Form:\n") print( np.round(B, decimals=4) + np.zeros(B.shape),"\n") # adding -0.0 + 0 gives 0.0 print() if ( type(b) != type(None) and noSolution(B) ): print("No solution!") C = backwardSubstitute(B,type(b)!=type(None),traceLevel) if( traceLevel > 0 ): print("Reduced Row Echelon Form:\n") print( np.round(C, decimals=4) + np.zeros(C.shape),"\n") # adding -0.0 + 0 gives 0.0 print() return C ######################################################################################################## def GE(A,b=None, traceLevel = 0): return GaussianElimination(A,b,traceLevel) ## Problem 1 In this problem you will write code to calculate the PLU decomposition of a matrix, AND the determinant in the case that the matrix is square, as discussed in lecture. Please review lectures 8B and 9 for a complete discussion of what you need to do; a brief summary is as follows. Gaussian Elimination can be implemented completely by matrix multiplication using elementary matrices to perform row operations; since we are only concerned with creating the echelon form (not the reduced form, with pivots = 1.0 and back substitution), we only need to perform forward elimination, using only row exchanges and row reductions (adding a multiple of one row to another). Therefore if we denote the echelon form (an upper-triangular matrix) by $U$, a row exchange by $P$, a row reduction by $R$, and either one of these by $E$, Gaussian Elimination with $i$ row exchanges and $j$ row reductions can be expressed as $$U\ =\ E_{i+j} E_{i+j-1} \ldots E_{2} E_{1} A$$ As explained in lecture, the row exchanges commute with the row reductions, so this can be written as $$U\ =\ R_{j}\ldots R_{1} P_{i} \ldots P_{1} A$$ or, since elementary matrices are invertible, as $$A\ =\ P_{1}\ldots P_{i} R^{-1}_{1} \ldots R^{-1}_{j} U.$$ Note that each $P_i$ is its own inverse, and to invert a row reduction, instead of adding a multiple of one row to another, simply *subtract*. The PLU decomposition is therefore $$P=P_{1}\ldots P_{i},$$ $$L=R^{-1}_{1} \ldots R^{-1}_{j},$$ and $U$ is the echelon form produced by forward elimination. In order to calculate $P$ and $L$, initialize each of these to $I_m$, where $m$ is the number of rows in $A$ ($A$ need not be square, but $P$ and $L$ will be square). Then, for each row operation, update $P$ and $L$ by *right* multiplying by the inverse of the appropriate elementary matrix: $$ P = P @ P_k\quad\quad\quad \text{ or }\quad\quad\quad L = L @ R^{-1}_k.$$ Again, review the lectures to understand this process so you can add code to the following, which is simply the code for `forwardElimination` from the previous code cell, with comment hints about what to do. You should not have to change anything except where "your code here" is indicated. Tests are given below. # Calculate the PLU decomposition of a matrix by keeping track of # the row ops and accumulating the appropriate matrices. # While we are at it, if the matrix is square, we calculate the determinant by # taking the product of the diagonal of U, and determine # the sign by counting the number of row exchanges # NOTE: i and j are matrix rows, R1 is row 0 def PLU_Decomposition(A,traceRowOps=False): A = (np.copy(A)).astype(float) if (traceRowOps): print("\nRunning Forward Elimination on:\n") print(np.round(A, decimals=4),"\n") print() if (traceRowOps): print("Creating Echelon Form...\n") (numRows,numCols) = A.shape # Now initialize the permutation matrix P, the lower triangular L, and also a counter for the number of # row exchanges. P and L will be square matrices with the same number of rows as A # YOUR CODE HERE P = np.eye(numRows) # just to get it to compile, you'll have to change these two L = np.eye(numRows) numExchs = 0 # number of row exchanges # r = row we are currently working on pivot value is A[r][c] r = 0 for c in range(numCols): # solve for variable in column c # find row in column c with first non-zero element, and exchange with row r r1 = r while (r1 < numRows): if (not np.isclose(A[r1][c],0.0)): # A[r1][c] is first non-zero element at or below r in column c break r1 += 1 if (r1 == numRows): # all zeros below r in this column continue # go on to the next column, but still working on same row if (r1 != r): # Exchange rows r1 and r tmp = A[r1].copy() A[r1] = A[r] A[r] = tmp # YOUR CODE HERE to count the number of row exchanges AND update P numExchs += 1 tmp = P[r1].copy() P[r1] = P[r] P[r] = tmp if (traceRowOps): print("Exchange R" + str(r1+1) + " and R" + str(r+1) + "\n") print(np.round(A, decimals=4)) print() # now use pivot A[r][c] to eliminate all vars in this column below row r for r2 in range(r+1,numRows): if(not np.isclose(A[r2][c],0.0)): factor = -A[r2][c]/A[r][c] A[r2] += factor * A[r] # YOUR CODE HERE to update L with inverse of the row reduction just done L[r2,c] = -factor if(traceRowOps): print("R" + str(r2+1) + " += " + str(np.around(factor,prec)) + "*R" + str(r+1) + "\n") print(np.round(A, decimals=4)) print() r += 1 if (r >= numRows): break # Note: A has actually been turned into U at this point U = A # Now calculate the determinant, if A is square # Determinant of U is PRODUCT of the diagonal elements # but the sign changes each time there is a row exchange. # If matrix is not square set d to None d = None # YOUR CODE HERE if numRows == numCols: # A is square d = U.diagonal().prod() * (-1)**numExchs return (P,L,U,d) A.size # (A) def test(A): (P,L,U,d) = PLU_Decomposition(A) print("A:\n",np.around(A,4),'\n') print("P:\n",np.around(P,4),'\n') print("L:\n",np.around(L,4),'\n') print("U:\n",np.around(U,4),'\n') if(type(d) != type(None)): print("d:\n",np.around(d,4),'\n') if(P.size > 0): print("P@L@U:\n",np.around(P@L@U,4)) # Test it! np.isclose(A,P@L@U).all() A = np.array([[1,2], [3,4]]) test(A) # (B) B = np.array([[1,-1,0,0], [-1,2,-1,0],[0,-1,2,-1],[0,0,-1,2]]) test(B) # (C) C = np.array([[1,2,3], [4,5,6],[7,8,9]]) test(C) #(D) D = np.array([[1,1,1], [1,2,2],[1,2,3]]) test(D) # (E) E = np.random.random((51,72)) test(E) ## Problem 2 You must use your code for `analyzeMatrix` from HW 02 for this problem. In the next cell, paste in your solution and then complete the following code template to calculate the eigenvalues of a 2x2 matrix. Further hints are written in the template. # Paste your code for analyzeMatrix here def analyzeMatrix(A): B = GaussianElimination(A) (numRows, numCols) = B.shape r = 0 pCols = [] basis = [] nbasis = [] rank = 0 for r in range(numRows): for c in range(numCols): if not np.isclose(B[r][c], 0.0): pCols.append(c) break for x in pCols: basis.append(np.vstack(A[:,x].copy().astype(float))) for y in range(numCols): if y not in pCols: nbasis.append(np.vstack(B[:,y]).copy().astype(float)) count = len(basis) for z in range(len(nbasis)): nbasis[z] = nbasis[z]*-1 if count < numCols: nbasis[z][count] = 1 count += 1 #print("rank is", len(basis)) #print("nullity is", len(nbasis)) #print("basis is", basis) #print("nbasis is", nbasis) return (len(basis), basis, len(nbasis), nbasis) # Solution: Calculate the eigenvalues of a 2 x 2 matrix # Use the following two function to determine the characteristic polynomial # compute the trace of a 2x2 matrix def trace(A): return A[0,0] + A[1,1] # compute the determinant of a 2x2 matrix def det(A): return A[0,0]*A[1,1] - A[1,0]*A[0,1] def eigen(A): # use the quadratic formula to calculate the roots of the characteristic polynomial a = 1 # just to get it to compile -- your code here b = -(A[0,0]+A[1,1]) # just to get it to compile -- your code here c = A[0,0]*A[1,1] - A[0,1]*A[1,0] # just to get it to compile -- your code here # d = discriminant (expression inside the square root) d = b**2 - 4*a*c # just to get it to compile -- your code here # determine whether there are only complex eigenvalues, 1 (duplicated) real eigenvalue, or 2 distinct real eigenvalues # then use analyzeMatrix to calculate the eigenbasis for each real eigenvalue. # Print out which case (based on discriminant), and any real eigenvalues with associated eigenspace # as shown below in Part A # Remember to use np.isclose(d,0.0) instead of d == 0.0. # Your code here iscomplex = False if b**2 - 4*a*c < 0: iscomplex = True if iscomplex: l1 = np.complex(-b/(2*a), -np.sqrt(np.abs(b**2-4*a*c)/(2*a))) l2 = np.complex(-b/(2*a), np.sqrt(np.abs(b**2-4*a*c)/(2*a))) else: l1 = (-b - np.sqrt(b**2 - 4*a*c))/(2*a) l2 = (-b + np.sqrt(b**2 - 4*a*c))/(2*a) (n_basis, basis, n_nbasis, nbasis) = analyzeMatrix(A) n = 1 if l1 != l2: n = 2 # check if real or complex eig_type = "real" if iscomplex: eig_type = "complex" if n == 1: print("one {} eigenvalue".format(eig_type)) else: print("two {} eigenvalues".format(eig_type)) print() print("lambda_1 = {:.1f}".format(l1)) print(basis) print() if n > 1: print("lambda_2 = {:.1f}".format(l2)) print(nbasis) ### (A) Test your code with the following matrix; you should design your code so that it prints out something like this: [[1 0] [0 0]] two distinct real eigenvalues lambda_1 = 1.0 [[1.] [0.]] lambda_2 = 0.0 [[0.] [1.]] Out[3]: (array([1., 0.]), array([[1., 0.], [0., 1.]])) A = np.array( [[1,0],[0,0]] ) # Has two distinct real eigenvalues print(A,'\n') eigen(A) # test it; remember that eigenvectors may be scaled differently and # this function returns eigenvectors as columns in a matrix np.linalg.eig(A) # (B) Test on the following matrix in the same way B = np.array( [[1,1],[1,1]] ) # Has two distinct real eigenvalues print(B,'\n') eigen(B) # test it; remember that eigenvectors may be scaled differently and # this function returns eigenvectors as columns in a matrix np.linalg.eig(B) # (C) Test on the following matrix in the same way C = np.array( [[-1,1],[0,-1]] ) # Should have only 1 eigenvalue, with 1D eigenspace print(C,'\n') eigen(C) # test it; remember that eigenvectors may be scaled differently and # this function returns eigenvectors as columns in a matrix np.linalg.eig(C) # (D) Test on the following matrix in the same way D = np.array( [[1,0],[0,1]] ) # Should have only 1 eigenvalue, with 2D eigenspace print(D,'\n') eigen(D) # test it; remember that eigenvectors may be scaled differently and # this function returns eigenvectors as columns in a matrix np.linalg.eig(D) # (E) Test on the following matrix in the same way E = np.array( [[1,3],[-1,-2]] ) # Has complex eigenvalues print(E,'\n') eigen(E) # test it; remember that eigenvectors may be scaled differently and # this function returns eigenvectors as columns in a matrix np.linalg.eig(E)