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