# coding: utf-8
# # 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.