{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# CSC338. Homework 4\n",
"\n",
"Due Date: Wednesday Feburary 5, 9pm\n",
"\n",
"Please see the guidelines at https://www.cs.toronto.edu/~lczhang/338/homework.html\n",
"\n",
"### What to Hand In\n",
"\n",
"Please hand in 2 files:\n",
"\n",
"- Python File containing all your code, named `hw4.py`.\n",
"- PDF file named `hw4_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.\n",
"\n",
"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.\n",
"\n",
"**Make sure to remove or comment out all matplotlib or other expensive code before submitting your homework!**\n",
"\n",
"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."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"import math\n",
"import numpy as np"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Question 1\n",
"\n",
"For this question, we will again start from code from tutorial 3."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Code from tutorial 3\n",
"\n",
"def backward_substitution(A, b):\n",
" \"\"\"Return a vector x with np.matmul(A, x) == b, where \n",
" * A is an nxn numpy matrix that is upper-triangular and non-singular\n",
" * b is an nx1 numpy vector\n",
" \"\"\"\n",
" n = A.shape[0]\n",
" x = np.zeros_like(b, dtype=np.float)\n",
" for i in range(n-1, -1, -1):\n",
" s = 0\n",
" for j in range(n-1, i, -1):\n",
" s += A[i,j] * x[j]\n",
" x[i] = (b[i] - s) / A[i,i]\n",
" return x\n",
"\n",
"def eliminate(A, b, k):\n",
" \"\"\"Eliminate the k-th row of A, in the system np.matmul(A, x) == b,\n",
" so that A[i, k] = 0 for i < k. The elimination is done in place.\"\"\"\n",
" n = A.shape[0]\n",
" for i in range(k + 1, n):\n",
" m = A[i, k] / A[k, k]\n",
" for j in range(k, n):\n",
" A[i, j] = A[i, j] - m * A[k, j]\n",
" b[i] = b[i] - m * b[k]\n",
"\n",
"def gauss_elimination(A, b):\n",
" \"\"\"Return a vector x with np.matmul(A, x) == b using\n",
" the Gauss Elimination algorithm, without partial pivoting.\"\"\"\n",
" for k in range(A.shape[0] - 1):\n",
" eliminate(A, b, k)\n",
" x = backward_substitution(A, b)\n",
" return x"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Part (a) -- 1 pt\n",
"\n",
"Solve the system $Ax = b$ for the values of $A$ and $b$ below using\n",
"the function `gauss_elimination` from last time.\n",
"Save the solution you obtain in the variable `soln_nopivot`."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"e = pow(2, -100)\n",
"A = np.array([[e, 1],\n",
" [1, 1]])\n",
"b = np.array([1 + e, 2])\n",
"\n",
"soln_nopivot = None"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Part (b) -- 2 pt\n",
"\n",
"Write a helper function `partial_pivot` that performs partial pivoting\n",
"on `A` at column `k`, so that the function \n",
"`gauss_elimination_partial_pivot` performs Gauss Elimination with Partial Pivoting."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"def partial_pivot(A, b, k):\n",
" \"\"\"Perform partial pivoting for column k. That is, swap row k\n",
" with row j > k so that the new element at A[k,k] is the largest\n",
" amongst all other values in column k below the diagonal.\n",
" \n",
" This function should modify A and b in place.\n",
" \"\"\"\n",
" # TODO\n",
"\n",
"def gauss_elimination_partial_pivot(A, b):\n",
" \"\"\"Return a vector x with np.matmul(A, x) == b using\n",
" the Gauss Elimination algorithm, with partial pivoting.\"\"\"\n",
" for k in range(A.shape[0] - 1):\n",
" partial_pivot(A, b, k)\n",
" eliminate(A, b, k)\n",
" x = backward_substitution(A, b)\n",
" return x"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Part (c) -- 1 pt\n",
"\n",
"Solve the system $Ax = b$ for the values of $A$ and $b$ below using `gauss_elimination_partial_pivot`. Save the solution you obtain in the variable `soln_pivot`."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"e = pow(2, -100)\n",
"A = np.array([[e, 1],\n",
" [1, 1]])\n",
"b = np.array([1 + e, 2])\n",
"\n",
"soln_pivot = None"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Part (d) -- 2 pt\n",
"\n",
"Do your answers in parts (a) and (d) match? If not, which is the correct\n",
"answer? Include your explanation in the PDF write-up.\n",
"\n",
"## Question 2\n",
"\n",
"### Part (a) -- 3 pt\n",
"\n",
"Consider the following matrices `M1`, `M2`, and `M3`. Compute each of their $L_1$, $L_2$, and $L_\\infty$ norms. Save the results in the variables below.\n",
"\n",
"For the $L_2$ norm you may find the function `np.linalg.norm` helpful. You can compute the $L_1$ and $L_\\infty$ norms either by hand or write a function."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"M1 = np.array([[3., 0.],\n",
" [-4., 2.]])\n",
"\n",
"M2 = np.array([[2., -2., 0.3],\n",
" [0.5, 1., 0.9],\n",
" [-4., -2., 5]])\n",
"\n",
"M3 = np.array([[0.2, -0.2],\n",
" [1.0, 0.2]])\n",
"\n",
"# fill in these answers\n",
"M1_l_1 = 0\n",
"M1_l_2 = 0\n",
"M1_l_infty = 0\n",
"M2_l_1 = 0\n",
"M2_l_2 = 0\n",
"M2_l_infty = 0\n",
"M3_l_1 = 0\n",
"M3_l_2 = 0\n",
"M3_l_infty = 0"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Part (b) -- 4 pt\n",
"\n",
"Show that an induced matrix norm $|| \\cdot ||$ has the property\n",
"$$||A + B|| \\leq ||A|| + ||B||$$\n",
"\n",
"### Part (c) -- 4 pt\n",
"\n",
"Is it true that for a vector, $||v||_\\infty \\le ||v||_1$? What about for a matrix: is it true that for a matrix, $||M||_\\infty \\le ||M||_1$? Include your solution and justification in your pdf writeup.\n",
"\n",
"## Question 3\n",
"\n",
"### Part (a) [3 pt]\n",
"\n",
"Write a function `matrix_condition_number` that computes the condition\n",
"number of a $2 \\times 2$ matrix. Use the $$L_1$$ matrix norm."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"def matrix_condition_number(M):\n",
" \"\"\"\n",
" Returns the condition number of the 2x2 matrix M.\n",
" Use the $L_1$ matrix norm.\n",
"\n",
" Precondition: M.shape == [2, 2] \n",
" M is non-singular\n",
"\n",
" >>> matrix_condition_number(np.array([[1., 0.], [0., 1.]]))\n",
" 1\n",
" \"\"\""
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Part (b) [2 pt]\n",
"\n",
"Classify each of the following matrices `A1`, `A2`, `A3` and `A4`,\n",
"as well-conditioned or ill-conditioned.\n",
"\n",
"You may do this question either by hand, or by using the function above.\n",
"\n",
"Save the classifications in a Python array called `conditioning`.\n",
"Each item of the array\n",
"should be either the string \"well\" or the string \"ill\"."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"A1 = np.array([[3, 0], \n",
" [0, pow(3, -30)]])\n",
"A2 = np.array([[pow(3, 20), 0],\n",
" [0, pow(5, 50)]])\n",
"A3 = np.array([[pow(4, -40), 0],\n",
" [0, pow(2, -80)]])\n",
"A4 = np.array([[pow(5, -17), pow(5, -16)],\n",
" [pow(5, -18), pow(5, -17)]])\n",
"\n",
"# whether A1, A2, A3, A4 is \"well\" or \"ill\" conditioned\n",
"conditioning = [None, None, None, None]"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Part (c) [2 pt]\n",
"\n",
"It should be immediate obvious that the matrix\n",
"\n",
"$$\\begin{bmatrix}100 & 2 \\\\ 201 & 4\\end{bmatrix}$$\n",
"\n",
"is ill-conditioned. Explain why. Include your answer in your pdf write-up.\n",
"\n",
"\n",
"### Part (d) [2 pt]\n",
"\n",
"Suppose that $A$ and $B$ are two $n \\times n$ matrices, and both are well-conditioned. Is $A(B^{-1})$ also well-conditioned? Why or why not? Include your answer and justifications in your pdf write-up. Be specific.\n",
"\n",
"### Part (e) [4 pt]\n",
"\n",
"Describe an efficient algorithm to compute $d^T B^T A^{-1} B d$\n",
"\n",
"Where:\n",
"\n",
"* $A$ is an invertible $n \\times n$ matrix, \n",
"* $B$ is an $n \\times n$ matrix, and \n",
"* $d$ is an $n \\times 1$ vectors\n",
"\n",
"Be clear and specific. Include your strategy in your pdf write-up."
]
}
],
"metadata": {},
"nbformat": 4,
"nbformat_minor": 2
}