{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Name: **Your name here** \n",
"UID: **Your student ID num here**"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Homework 2: More Linear Algebra "
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"# Setup the environment - do not modify this cell, but run it before anything else\n",
"import numpy as np\n",
"from scipy.signal import convolve2d\n",
"from scipy.linalg import hilbert\n",
"from numpy.random import randn, normal\n",
"from numpy.linalg import norm, inv\n",
"from numpy.fft import fft2, ifft2\n",
"import urllib\n",
"import matplotlib.pyplot as plt\n",
"np.random.seed(0)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Problem 1 - condition number\n",
"Run the following code. "
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Recovery error (clean) = 1.63e-06\n",
"measurement error = 1.79e-07\n",
"Recovery error (noisy) = 178.930\n"
]
}
],
"source": [
"# Do not modify this block!\n",
"# Create a linear system Ax=b\n",
"n=8\n",
"A = hilbert(n) # construct an n by n matrix\n",
"x = randn(n,1) # construct n by 1 signal \n",
"b = A@x # Note that '@' is matrix multiplication in python, while '*' denotes entry-wise multiplication\n",
"\n",
"# Solve the system to recover x\n",
"x_recovered = inv(A)@b \n",
"print('Recovery error (clean) = %0.3g'%norm(x_recovered-x))\n",
"\n",
"# Add some noise\n",
"b_noise = b+randn(n,1)*0.0000001\n",
"print('measurement error = %0.3g'%norm(b_noise-b))\n",
"\n",
"# Solve the noisy system\n",
"x_noise = inv(A)@b_noise\n",
"print('Recovery error (noisy) = %0.3f'%norm(x_noise-x))\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Why did this happen? Write a few sentences.\n",
"*Put your work here*"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"In general, suppose that I want to solve $Ax=b$ to find $x$, but I have noisy measurements $\\hat b.$ To do this, I compute $\\hat x = A^{-1} \\hat b.$ \n",
"**Prove the following**\n",
"$$\\frac{\\|x-\\hat x\\|}{\\|x\\|} \\le \\kappa \\frac{\\|b-\\hat b\\|}{\\|b\\|}, $$\n",
"where $\\kappa$ is the condition number of $A.$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Proof\n",
"*Put your work here*"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Problem 2 - Adjoints\n",
"Suppose you have two functions, `A` an `At` that each implement linear operators. Write a *randomized* method for checking whether `At` is the adjoint of `A`. Your test should directly verify the definition of the adjoint\n",
" $$\\langle A(x), y \\rangle = \\langle x, At(y)\\rangle$$\n",
"where $\\langle \\cdot,\\cdot \\rangle$ denotes the Hermitian inner product. When $x$ and $h$ happen to be column vectors, this condition becomes\n",
" $$A(x)^H y = x^H At(y)$$\n",
" where $x^H$ is the Hermitian transpose of $x.$\n",
" The arguments of the checker method are the functions `A` and `At`, and a tuple containing the dimensions of the argument to `A`. The method returns `True` if the methods are adjoints of one another, and `False` otherwise.\n",
" \n",
" Your method must work for inputs $x$ of any dimension and shape.\n",
" Note: If you choose to use the numpy transpose operator in your solution, remember that this built-in operator is not the Hermitian transpose - It's just a regular transpose without taking the conjugate."
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [],
"source": [
"def check_adjoint(A,At,dims):\n",
" # start with this line - create a random input for A()\n",
" x = normal(size=dims)+1j*normal(size=dims)\n",
" \n",
" # FILL IN METHOD BODY HERE\n",
" \n",
" # End with something like this (change these lines if you want)...\n",
" if relative_error < 1e-10:\n",
" print('Adjoint Test Passed, rel_diff = %s'%rel_error)\n",
" return True\n",
" else:\n",
" print('Adjoint Test Failed, rel_diff = %s'%rel_error)\n",
" return False"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"After filling in the body of your method, run this unit test to make sure it's ok."
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [
{
"ename": "NameError",
"evalue": "name 'relative_error' is not defined",
"output_type": "error",
"traceback": [
"\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
"\u001b[0;31mNameError\u001b[0m Traceback (most recent call last)",
"\u001b[0;32m\u001b[0m in \u001b[0;36m\u001b[0;34m\u001b[0m\n\u001b[1;32m 5\u001b[0m \u001b[0mA\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;32mlambda\u001b[0m \u001b[0mx\u001b[0m\u001b[0;34m:\u001b[0m \u001b[0mfft2\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mx\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 6\u001b[0m \u001b[0mAt\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;32mlambda\u001b[0m \u001b[0mx\u001b[0m\u001b[0;34m:\u001b[0m \u001b[0mifft2\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mx\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 7\u001b[0;31m \u001b[0mresult\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mcheck_adjoint\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mA\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mAt\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;36m100\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;36m200\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 8\u001b[0m \u001b[0;32massert\u001b[0m \u001b[0;34m(\u001b[0m\u001b[0;32mnot\u001b[0m \u001b[0mresult\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m\"Adjoint test should have failed, but succeeded! Double check your solution!\"\u001b[0m \u001b[0;31m# Throw an exception if the result is unexpected\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 9\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
"\u001b[0;32m\u001b[0m in \u001b[0;36mcheck_adjoint\u001b[0;34m(A, At, dims)\u001b[0m\n\u001b[1;32m 6\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 7\u001b[0m \u001b[0;31m# End with something like this (change these lines if you want)...\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 8\u001b[0;31m \u001b[0;32mif\u001b[0m \u001b[0mrelative_error\u001b[0m \u001b[0;34m<\u001b[0m \u001b[0;36m1e-10\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 9\u001b[0m \u001b[0mprint\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m'Adjoint Test Passed, rel_diff = %s'\u001b[0m\u001b[0;34m%\u001b[0m\u001b[0mrel_error\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 10\u001b[0m \u001b[0;32mreturn\u001b[0m \u001b[0;32mTrue\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
"\u001b[0;31mNameError\u001b[0m: name 'relative_error' is not defined"
]
}
],
"source": [
"# Adjoint unit test - DO NOT MODIFY THIS BLOCK\n",
"# This method will throw a nasty exception if your code doesn't perform as expected\n",
"\n",
"# Test #1: This test should fail because the standard DFT is not self adjoint\n",
"A = lambda x: fft2(x)\n",
"At = lambda x: ifft2(x)\n",
"result = check_adjoint(A,At,(100,200))\n",
"assert (not result), \"Adjoint test should have failed, but succeeded! Double check your solution!\" # Throw an exception if the result is unexpected\n",
"\n",
"# Test #2: This test should pass though, because F^H = conj(F). \n",
"A = lambda x: fft2(x)\n",
"At = lambda x: np.conj(fft2(np.conj(x)))\n",
"result = check_adjoint(A,At,(100,200))\n",
"assert result, \"Adjoint test should have succeeded, but failed! Double check your solution!\" # Throw an exception if the result is unexpected\n",
"\n",
"print(\"Tests PASSED! You're on your way to understanding linear operators!\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Problem 3 - Convolutions\n",
"When computing total variation, you need to produce the image gradient, which containts the horizontal and vertical differences between adjacent pixels in a 2d array.\n",
"Choose the kernels below so that the methods `gradh` and `gradv` produce differences (discrete gradients) in the horizontal and vertical directions. Remember that convolve2d assumes the middle element of the kernel (array index 1 for a kernel of length 3) is the center of the kernel. This differs from standard convolutions, in which array index 0 is the center of the kernel.\n",
"\n",
"Output i,j of the horizontal gradients should contain `x[i,j+1]-x[i,j]`, while the veritical graident should contain `x[i+1,j]-x[i,j]`."
]
},
{
"cell_type": "code",
"execution_count": 71,
"metadata": {},
"outputs": [],
"source": [
"# Put your kernels here! \n",
"kernel_h = # Complete this line of code by defining a 3-element, 2d array\n",
"kernel_v = # Complete this line of code by defining a 3-element, 2d array"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"** Now, run the cell below. It will create the gradient operators using your kernels, and unit test them. Do NOT modify any of the code in the cell below.**"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Do not modify ANYTHING in this cell. \n",
"def gradh(x):\n",
" \"\"\"Discrete gradient/difference in horizontal direction\"\"\"\n",
" return convolve2d(x,kernel_h, mode='same', boundary='wrap')\n",
"def gradv(x):\n",
" \"\"\"Discrete gradient/difference in vertical direction\"\"\"\n",
" return convolve2d(x,kernel_v, mode='same', boundary='wrap')\n",
"def grad2d(x):\n",
" \"\"\"The full gradient operator: compute both x and y differences and return them all. The x and y \n",
" differences are stacked so that rval[0] is a 2D array of x differences, and rval[1] is the y differences.\"\"\"\n",
" return np.stack([gradh(x),gradv(x)])\n",
"\n",
"# Perform unit tests - this will throw exceptions if your method is screwed up!\n",
"x = randn(10,20)\n",
"ghx = gradh(x)\n",
"assert ghx[0,0] == x[0,1] - x[0,0], 'Failed test 1'\n",
"assert ghx[0,-1] == x[0,0] - x[0,-1], 'Failed test 2'\n",
"assert ghx[1,1] == x[1,2] - x[1,1], 'Failed test 3'\n",
"gvx = gradv(x)\n",
"assert gvx[0,0] == x[1,0] - x[0,0], 'Failed test 4'\n",
"assert gvx[-1,0] == x[0,0] - x[-1,0], 'Failed test 5'\n",
"assert gvx[1,1] == x[2,1] - x[1,1], 'Failed test 6'\n",
"print('TESTS PASSED! YOU ROCK!')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"**Now, implement the adjoint/transpose of these operators!** No looping allowed! Your implementation of `gradht` and `gradvt` must call `convolve2d` exactly once. Ideally, you'll only write 1 line of code per line."
]
},
{
"cell_type": "code",
"execution_count": 73,
"metadata": {},
"outputs": [],
"source": [
"# Fill in the implementations of these functions\n",
"def gradht(x):\n",
" \"\"\"Adjoint of gradh\"\"\"\n",
" # Your work here!\n",
" return # Put a return value here\n",
"def gradvt(x):\n",
" \"\"\"Adjoint of gradv\"\"\"\n",
" # Your work here!\n",
" return # Put a return value here\n",
"def divergence2d(x):\n",
" \"The methods is the adjoint of grad2d.\"\n",
" # Your work here!\n",
" return # Put a return value here"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"**After writing the adjoint routines, run the unit test below!**"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Do not modify ANYTHING in this block!\n",
"is_pass = check_adjoint(gradh, gradht,(10,20))\n",
"assert is_pass, 'Your gradht method is not the adjoint of gradh.'\n",
"\n",
"is_pass = check_adjoint(gradv, gradvt,(10,20))\n",
"assert is_pass, 'Your gradvt method is not the adjoint of gradv.'\n",
"\n",
"is_pass = check_adjoint(grad2d, divergence2d,(10,20))\n",
"assert is_pass, 'Your divergence2d method is not the adjoint of grad2d.'\n",
"\n",
"print(\"Unit tests PASSED! You're getting really good at this!\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Problem 4 - FFT \n",
"Now, re-implement these methods using the FFT!\n",
"Your code must call `np.fft.fft2()` and `np.fft.ifft2(x)`, and you cannot call convolve2d. No loops allowed! Remember, when you convolve things using the FFT, you're relying on the convolution theorem. This theorem assumes the center of the kernel is at index 0. Also, we're using the 2D DFT here rather than the 1D DFT from your last homework. All the basic ideas still apply. "
]
},
{
"cell_type": "code",
"execution_count": 75,
"metadata": {},
"outputs": [],
"source": [
"def gradh_fft(x):\n",
" \"\"\"Discrete gradient/difference in horizontal direction\"\"\"\n",
" # Fill in method body\n",
" return # Fill in return value\n",
"def gradv_fft(x):\n",
" \"\"\"Discrete gradient/difference in vertical direction\"\"\"\n",
" # Fill in method body\n",
" return # Fill in return value"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"**Now, run the unit tests below!**"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Do not modify ANYTHING in this cell\n",
"# create random array\n",
"x = randn(100,200)\n",
"\n",
"# verify that gradh_fft = gradh\n",
"h_error = norm(gradh_fft(x)-gradh(x))\n",
"print('Horizontal error = ', h_error)\n",
"assert h_error<1e-10, 'Horizontal FFT gradient is incorrect!'\n",
"\n",
"# verify that gradv_fft = gradv\n",
"v_error = norm(gradv_fft(x)-gradv(x))\n",
"print('Vertical error = ',v_error)\n",
"assert v_error<1e-10, 'Vertical FFT gradient is incorrect!'\n",
"\n",
"print(\"Tests PASSED! Wow - you're a linear algebra GENIUS!\")\n",
"\n",
"f = urllib.request.urlopen(\"https://www.cs.umd.edu/~tomg/img/important_memes/good_job_cat.png\")\n",
"a = plt.imread(f)\n",
"fig = plt.imshow(a)\n",
"fig.axes.get_xaxis().set_visible(False)\n",
"fig.axes.get_yaxis().set_visible(False)\n",
"plt.show()"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.7.6"
}
},
"nbformat": 4,
"nbformat_minor": 2
}