## Instructions

**Objective**

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