CSC338. Homework 3¶
Due Date: Wednesday January 29, 9pm
Please see the guidelines at https://www.cs.toronto.edu/~lczhang/338/homework.html
What to Hand In¶
Please hand in 2 files:
• Python File containing all your code, named hw3.py.
• PDF file named hw3_written.pdf containing your solutions to the written parts of the assignment. Your solution can be hand-written, but must be legible. Graders may deduct marks for illegible or poorly presented solutions.
If you are using Jupyter Notebook to complete the work, your notebook can be exported as a .py file (File -> Download As -> Python). Your code will be auto-graded using Python 3.6, so please make sure that your code runs. There will be a 20% penalty if you need a remark due to small issues that renders your code untestable.
Make sure to remove or comment out all matplotlib or other expensive code before submitting your homework!
Submit the assignment on MarkUs by 9pm on the due date. See the syllabus for the course policy regarding late assignments. All assignments must be done individually.
In [ ]:
import math
import numpy as np
Question 1¶
Part (a) — 4 pt¶
Solve the following system $A{\bf x} = {\bf b}$ using Gauss Elimination by hand. Show all your steps. Show all your steps in your pdf writeup.
\begin{align*} A = \begin{bmatrix} 1 & 0 & -2 \\ -1 & 1 & 1 \\ 0 & 2 & -1 \end{bmatrix}, b = \begin{bmatrix} -1 \\ 2 \\ 3 \\ \end{bmatrix}, \end{align*}
Part (b) — 2 pt¶
Using your work from part (a), find the LU factorization of $A$. Show all your steps in your pdf writeup.
Question 2¶
For this question, you can refer to and use any of the code that we wrote together in tutorial 3:
In [ ]:
# Code from tutorial 3
def backward_substitution(A, b):
“””Return a vector x with np.matmul(A, x) == b, where
* A is an nxn numpy matrix that is upper-triangular and non-singular
* b is an nx1 numpy vector
“””
n = A.shape[0]
x = np.zeros_like(b, dtype=np.float)
for i in range(n-1, -1, -1):
s = 0
for j in range(n-1, i, -1):
s += A[i,j] * x[j]
x[i] = (b[i] – s) / A[i,i]
return x
def eliminate(A, b, k):
“””Eliminate the k-th row of A, in the system np.matmul(A, x) == b,
so that A[i, k] = 0 for i < k. The elimination is done in place."""
n = A.shape[0]
for i in range(k + 1, n):
m = A[i, k] / A[k, k]
for j in range(k, n):
A[i, j] = A[i, j] - m * A[k, j]
b[i] = b[i] - m * b[k]
def gauss_elimination(A, b):
"""Return a vector x with np.matmul(A, x) == b using
the Gauss Elimination algorithm, without partial pivoting."""
for k in range(A.shape[0] - 1):
eliminate(A, b, k)
x = backward_substitution(A, b)
return x
Part (a) -- 3 pt¶
Write a function forward_substitution that solves the lower-triangular system $Ax = b$.
Hint: This function should be very similar to the backward_substitution function from the tutorial.
In [ ]:
def forward_substitution(A, b):
"""Return a vector x with np.matmul(A, x) == b, where
* A is an nxn numpy matrix that is lower-triangular and non-singular
* b is a nx1 numpy vector
>>> A = np.array([[2., 0.],
[1., -2.]])
>>> b = np.array([1., 2.])
>>> forward_substitution(A, b)
array([ 0.5 , -0.75])
“””
Part (b) — 4 pt¶
Write a function elementary_elimination_matrix that returns the $k$-th elementary elimination matrix ($M_k$ in your notes).
You may assume that A[i,j] = 0 for i > j, j < k - 1. (The subtraction in $k-1$ is because in Python, indices begin at 0.)
In [ ]:
def elementary_elimination_matrix(A, k):
"""Return the elements below the k-th diagonal of the
k-th elementary elimination matrix.
(Do not use partial pivoting, since we haven't
introduced the idea yet.)
Precondition: A is an nxn numpy matrix, non-singular
A[i,j] = 0 for i > j, j < k - 1
As always, these examples are for your understanding only.
The actual Python output might differ slightly.
>>> A = np.array([[2., 0., 1.],
[1., 1., 0.],
[2., 1., 2.]])
>>> elementary_elimination_matrix(A, 1)
np.array([[1., 0., 0.],
[-0.5, 1., 0.],
[-1., 0., 1.]])
>>> A = np.array([[2., 0., 1.],
[0., 1., 0.],
[0., 1., 2.]])
>>> elementary_elimination_matrix(A, 2)
np.array([[1., 0., 0.],
[0., 1., 0.],
[0., -1., 1.]])
“””
Part (c) — 4pt¶
Write a function lu_factorize that factors a matrix A into its upper and lower triangular components. Use the function elementary_elimination_matrix as a helper.
In [ ]:
def lu_factorize(A):
“””Return two matrices L and U, where
* L is lower triangular
* U is upper triangular
* and np.matmul(L, U) == A
>>> A = np.array([[2., 0., 1.],
[1., 1., 0.],
[2., 1., 2.]])
>>> L, U = lu_factorize(A)
>>> L
array([[1. , 0. , 0. ],
[0.5, 1. , 0. ],
[1. , 1. , 1. ]])
>>> U
array([[ 2. , 0. , 1. ],
[ 0. , 1. , -0.5],
[ 0. , 0. , 1.5]])
“””
Part (d) — 2pt¶
Write a function solve_lu that solves a linear system Ax=b by
* factoring $A = LU$ (using the `lu_factorize` function)
* solving $Ly = b$ using forward substitution (using the `forward_substitution` function)
* solving $Ux = y$ using backward substitution (using the `backward_substitution` function)
In [ ]:
def solve_lu(A, b):
“””Return a vector x with np.matmul(A, x) == b using
LU factorization. (Do not use partial pivoting, since we haven’t
introduced the idea yet.)
“””
Part (e) — 5 pt¶
Write a function invert_matrix that takes an $n \times n$ matrix A, and computes its inverse by solving $n$ systems of linear equations of the form $Ax = b$.
You can use the functions that you wrote earlier, but note that you will be graded on efficiency. In particular, avoid writing code that repeats a computation unnecessarily.
In [ ]:
def invert_matrix(A):
“””Return the inverse of the nxn matrix A by
solving n systems of linear equations of the form
Ax = b.
>>> A = np.array([[0.5, 0.],
[-1., 2.]])
>>> invert_matrix(A)
array([[ 2. , -0. ],
[ 1. , 0.5]])
“””
Question 3¶
Part (a) — 4 pt¶
Prove that the matrix
\begin{align*} A = \begin{bmatrix} 0 & 1 \\ 1 & 0 \\ \end{bmatrix} \end{align*}
has no $LU$ facotirzation. That is, there are no lower triangular matrix $L$ and upper triangular matrix $U$ such that $A = LU$.
Include your proof in your pdf writeup.
Part (b) — 2 pt¶
Show that for an elementary matrix $M_k = I – {\bf m}{\bf e_k^T}$, we have $M_k^{-1} = I + {\bf m}{\bf e_k^T}$. (Property 4 of the elementary elimination matrix from your lecture notes)
Include your solution in your pdf writeup.