ECSE 343 Assignment 0 to 4 solutions

Description

ECSE 343 Assignment 0: Floating Point and Vectorized Python

1.2 Installing Python, Libraries and Tools
The assignments and group project are to be implemented using Python 3.7 and the latest version of the
numpy and matplotlib libraries.
You are free to install and configure your development environment, including IDE choices, as you wish. One
popular, more self-contained installation package that we recommend is the Individual Edition of the
Anaconda software installer.
After following your platform specific installation instructions, the Anaconda Navigator provides a simple
graphical interface that allows you to:
define isolated development environments with the appropriate Python version (3.7)
download and install the required libraries (numpy and matplotlib), including their dependencies,
into the environment, and
[optionally] to pick between a variety of development IDEs.
1.3 Python Language and Library Usage Rules
Python is a powerful language, with many built-in features. You should feel free to explore the base
language features and apply them as a convenience, whenever necessary. A good example is that, if you
need to sort values in a list before plotting them, you should feel free to use the built-in sort function rather
than implementing your own sorting algorithm (although, that’s perfectly fine, too!):
Every time you submit a new file, your previous submission will be overwritten. We will only grade
the final submitted file, so feel free to submit often as you progress through the assignment.
ⓘ
We will, however, draw exceptions to the use of (typically external) library routines that allow you to shortcut
through the core learning objective(s) of an assignment. For example, if we ask you to develop a linear
solver and apply it to a problem, and you instead rely on calling one of numpy’s built-in solvers, you will be
penalized. When in doubt as to whether a library (or even a built-in) routine is “safe” to use in your solution,
please contact the TA.
To help, the (purposefully minimal) base code we provide you with includes a superset of all the library
imports we could imagine you using for the assignment.
As you may have noticed, this course will rely heavily on numpy — in fact, you’ll likely learn just as much
about the power (and peculiarities) of the Python programming language as you will about the numpy library.
This library not only provides convenience routines for matrix, vector and higher-order tensor operations,
but also allows you to leverage high-performance vectorized operations if you’re careful about restructuring
your data/code in a vectorizable form. This assignment will briefly expose you to some of these nuances;
those of you with MATLAB experience will find this coding paradigm familiar, but also a little painful, since
myFavouritePrimes = [11, 3, 7, 5, 2]
# In ECSE 343, learning how to sort a list is NOT a core learning objective
myFavouritePrimes.sort() # 100% OK to use this!
print(myFavouritePrimes) # Output: [2, 3, 5, 7, 11]
import matplotlib.pyplot as plt # for plotting
import numpy as np # all of numpy, at least for this assignment
Python 3.7 has a built-in convenience breakpoint() function which will break code execution
into a debugger, where you can inspect variables in the debug REPL and even execute code! This
is a very powerful was to test your code as it runs and to tinker (e.g., inline in the REPL) with the
function calling conventions and input/output behaviour of code.
Be careful, as you can change the execution state (i.e., the debug environment is not isolated
from your scripts execution stack and heap), if you insert REPL code and then continue the
execution of your script from the debugger.
ⓘ
You must not use any additional imports for your solution, other than the ones provided by
us.
Doing so will result in a score of zero (0%) on the assignment.
☒
many foundational conventions are different (e.g., 0- versus 1-indexing, column- vs. row-major default
behaviours, broadcasting conventions, etc.)
1.4 Let’s get started… but let’s also not forget…
With these preliminaries out of the way, we can dive into the assignment tasks. Future assignment
handouts will not include these preliminaries, although they will continue to hold true. Should you
forget, they will remain online in this handout for your reference throughout the semester.
2 Floating point number systems
Any floating point number system can be characterized by four parameters, , where
is the base of the number system,
is its precision (i.e., number of significant digits),
is the lower bound on the exponent , and
is the upper bound on the exponent .
Given an instance of any such system, we can express a real number in its floating point representation
as:
where the base is an integer greater than 1, the exponent is an integer between and (inclusive, i.e.,
) and the digits are integers in the range . The number is usually an
approximation of , unless it happens to fall on one of the (relatively few) numbers that can be perfectly
represented in the floating point system. For any non-zero value, we normally force by adjusting the
exponent so that leading zeros are dropped. As such, the smallest (in magnitue) perfectly representable
non-zero number has a mantissa of .
(𝛽, 𝑡, 𝐿, 𝑈)
𝛽
𝑡
𝐿 𝑒
𝑈 𝑒
𝑥
fl(𝑥)
fl(𝑥) ≡ ± ×
s⏟ign
⎛
⎝
⎜
⎜
⎜
⎜
+ + ⋯ +
𝑑0
𝛽
0
𝑑1
𝛽
1
𝑑𝑡−1
𝛽
𝑡−1

mantissa
⎞
⎠
⎟
⎟
⎟
⎟
𝛽
𝑒
exp
⏟
onent
𝛽 𝑒 𝐿 𝑈
𝐿 ≤ 𝑒 ≤ 𝑈 𝑑𝑖 0 ≤ 𝑑𝑖 ≤ 𝛽 − 1 fl(𝑥)
𝑥
𝑑0 ≠ 0
𝑒
(1.0 ⋯ 0)𝛽
2.1 Fictitious Floating Point Systems
Let’s get a better sense of how far removed representable floating point numbers can get from real
numbers.
Here’s a code snippet with an example visualization of the expected output of this function for a fictitious
floating point system with , , and . Any real numbers that don’t fall exactly on
the stars in the plot below cannot be perfectly represented in this fictitious floating point number system.
The perfectly representable real numbers in a fictitious floating point system example.
# Plot an asterisk at each perfectly representable value along the real line
plt.title(‘Perfectly representable numbers for the $(2,2,-1,2)$ floating point
system’)
tmp = FloatSystem(2, 2, -1, 2)
plt.plot(tmp, zeros(tmp.shape[0]), ‘*’)
plt.yticks([])
plt.show()
Deliverable 1 [20%]
Complete the implementation of the function FloatSystem that returns a 1D numpy.array of all
the perfectly representable numbers (in base 10) for a given floating point number system. You
can safely ignore the NaN and cases, and explicitly add 0 to your representable numbers
output (i.e., without having to treat any special case exponent settings).
☑
±∞
𝛽 = 2 𝑡 = 2 𝐿 = −1 𝑈 = 2
3 Vectorizing Python with numpy
Writing numerical code requires balancing several (sometimes competing) design parameters: correctness
of the code, numerical stability (i.e., in the floating point sense of the word), and the efficiency and
scalability of the code are among these parameters.
Python is undoubtedly a flexible and powerful language, affording numerical software developers with many
tools to tackle their modeling and simulation tasks — however, as an interpreted language, Python’s
performance cannot compete with lower-level optimized code generated from, e.g., a compiler. Luckily,
Python allows for callable modules and libraries that need not be implemented in Python but rather in any
number of higher performance compiled languages. Moreover, Python’s vast ecosystem of specialized
libraries often comprise high-performance compiled backends: in this sense, Python serves just as much as
a high-level “glue” language as it does as a standalone one.
For numerical computation, numpy is one such library that is implemented in highly-optimized machine
code. When used appropriately, numerical code implemented in a manner that leverages numpy’s ability to
efficiently perform data-parallel operations over vectors and higher-order tabular data can be several
orders of magnitude more efficient than its Python-only equivalent.
One could easily teach an entire course on how to write efficient numpy code, and that is not the main goal
of ECSE 343; however, learning to think about numerical operations in vectorized form whenever
appropriate will open up the opportunity for cleaner, more readable and (much) more efficient code.
3.1 Slicing and dicing numpy.arrays
Multi-dimensional arrays are fundamental data structures in numerical computation, and numpy implements
sophisticated indexing schemes that respect specialized broadcasting rules, in addition to treating multidimensional arrays as first-class objects in all the library’s exposed functions.
The next deliverable will give you a brief sense of the power and flexibility of some of numpy’s indexing
notation. It is not meant, by any means, to be comprehensive; instead, the learning goal here is to open the
door to your independent exploration of numpy in order to facilitate implementation tasks, e.g., in future
assignments.
Deliverable 2 [20%]
This is a written answer-only deliverable: answer these questions using Python comments (i.e.,
not code) in your submission .py file. Feel free to experiment with indexing schemes using the
Python REPL or the __main__, in support of your written answers (i.e., do not regurgitate online
☑
Without coming anywhere close to exhaustively enumerating effecient numpy coding practices, to first-order
approximation, the following advice is a good place to start:
avoid for loops by restructuring data (if needed) into (potentially high-dimension) numpy.arrays, and
then performing operations across subsets of the data,
map and reduce strategies, often applied across numpy.array dimensions are a common strategy,
and
leverage numpy’s many built-in vectorized conditional, mathematical, and logical utility functions.
3.2 Avoiding for loops
documentation but, rather, run tiny code snippets to support any understanding you gain with the
support of online documentation.)
Answer the questions below given two numpy.array variables, the first (a) has a shape of (3,)
and the second (b) has a shape of (4,5):
1. What do b[0:3] and b[:, 0:3] do?
2. Why does b[:,:] = b[1,:] work and how would you make b[:,:] = b[:,1] work? Hint:
one of the more elegant solutions relies on using None.
3. Write a one-line Python/numpy expression that returns a numpy.array with shape of (5,)
with elements [a[0], a[1], a[2], a[0], a[1]], but without this explicit parameter list
(i.e., your answer should not be numpy.array([a[0], a[1], a[2], a[0], a[1]]) but
should use an indexing expression on a.) Hint: one of the many valid ways to do this
requires using a Python list comprehension.
4. How are the outputs of a[2], a[[2]] and a[[2,np.newaxis]] different?
. Why is a[a % a.shape[0]] guaranteed to work whereas a[a] may not?
Deliverable 3 [20%]
Perhaps the simplest example of a vectorizable mathematical operation is the computation of the
scalar dot product of two 1D vectors. Complete the implementation of the SlowDotProduct and
FastDotProduct functions, both of which accept two 1D numpy.arrays (you can assume they
have equal shape) and returns the scalar dot product of the two vectors. The SlowDotProduct
routine should not use any numpy utility functions and, instead, rely on Python language math
operators and for loops. The FastDotProduct function should instead use the full functionality
of numpy.
☑
4 Computer arithmetic
Let’s learn about errors due to floating point number representations, using simple numerical algorithms as
examples.
4.1 Catch the NaN
In the IEEE floating point standard, NaNs are “infectious”. You can use the Python REPL to explore how NaNs
behave with different arithmetic, inequality and logical operations. To do so, you can use numpy’s NaN as nan
= np.float64(“nan”) and then tinker with expressions like, e.g., nan < 4 or nan == nan or 3 + nan, etc.
4.2 Floating point errors
There are many types of floating point error your numerical code will be susceptible to, some of which you
can guard against… some of the time. Examples include
catastrophic cancellation resulting in loss in significant digits when, e.g.,
summing numbers with different scales, e.g., , or
In Python, unlike other languages like C/C++/Java, any expression containing an integer and floating
point arithmetic is automatically cast to the highest precision (64-bit) floating point representation to
favour precision over efficiency. While Python itself is a loosely typed language, numpy has explicit
facilities to force underlying numerical type representations, such as numpy.float64 and
numpy.float32 for double- and single-precission floating point values.
Deliverable 4 [20%]
Tracking down a NaN in a numerical algorithm can sometimes be a real pain. You will implement a
routine CatchTheNaN that identifies the source of a single spurious NaN in a very simple numerical
algorithm. The function takes as input a single (potentially very large) 2D matrix with many
NaNs in it. Behind the scenes (i.e., before the function is called), it turns out that was
constructed as the outer product of two 1D column vectors and , i.e., such that .
Both and each have a single NaN (i.e., in only one of their vector elements), and CatchTheNaN
should return an np.array of shape (2,) of the indices of these NaNs, where element [0] is the
index of the NaN in and element [1] is the index of the NaN in . We will grade you on the
correctness and performance of your code, i.e., avoid using brute force search with for loops.
☑
𝐌
𝐌
𝐱 𝐲 𝐌 = 𝐱 𝐲
T
𝐱 𝐲
𝐱 𝐲
10 +
38 10
−24
summing numbers with differences in low-precision bits, e.g.,
,
round-off errors that accumulate, e.g., , and
overflow and underflow, e.g., , to name a few.
Consider the following concrete example: when evaluating the natural logarithm of the sum of
exponentials of a set of values ,
the differences in the magnitude of the exponential terms of the summands may lead to overflow during
accumulation, depending on the value of .
After careful mathematical manipulation, leveraging the properties of logarithms and exponentials, we can
rewrite this equation as
where is the maximum value among the ,
This alternative mathematical formulation — when implemented in a floating point system — will be more
stable to the variations in the scale of the exponential terms.
1.38383 × 10 − 1.38382 ×
38 10
38
∑ 1.0
10
38
𝑖=0
10 ×
24 10
50
𝑁
𝑥𝑖
𝑦 = log (
exp( )) ∑ ,
𝑖=1
𝑁
𝑥𝑖
𝑁
𝑦 = 𝑚 + log (
exp( − 𝑚)) ∑ ,
𝑖=1
𝑁
𝑥𝑖
𝑚 𝑥𝑖
𝑚 = max .
𝑖
𝑥𝑖
Deliverable 5 [20%]
Implement and compare a naive version of this sum (in LogSumExpNaive) and a more numerically
robust version (in LogSumExpRobust). Each function takes a single 1D numpy.array with a shape
of (N,) as input and returns the scalar sum. You can use the __main__ test suite to explore the
differences in the outputs of these functions, particularly for large and/or elements in
different ranges of minimum and maximum magnitude.
Include a brief (i.e., one- to two-sentence long) comment in your solution implementation
of LogSumExpRobust that explains why this reformulation leads to more robust floating
point calculations.
☑
𝑁 𝑥𝑖
5 You’re Done!
Congratulations, you’ve completed the zero
th assignment. Review the submission procedures and
guidelines at the start of the handout before submitting the Python script file with your assignment solution.
formatted by Markdeep 1.13
✒

ECSE 343 Assignment 1: Matrix Factorizations and Solving Linear Systems

1 LU Solver
You will implement a simplified LU linear system solver. This will comprise a bare bones LU decomposition, as
well as the forward and backward substitution algorithms. Later in the assignment, you’ll test your solver on
polynomial regression problems.
1.1 LU Decomposition
Your first task is to decompose a matrix as the product of a lower triangular matrix and an upper
triangular matrix . Your decomposition will not treat pivoting. In this setting, the elements of and can
be expressed algebraically as:
𝐀 𝐋
𝐔 𝐋 𝐔
1.2 Forward and Backward Substitution
Given a lower triangular linear system , forward substitution solves for as:
Given an upper triangular linear system , backward substitution solves for as:
𝑈𝑖𝑗 =
⎧
⎩
⎨
⎪
⎪
𝐴𝑖𝑗
𝐴𝑖𝑗 − ∑
𝑘=0
𝑖−1
𝐿𝑖𝑘𝑈𝑘𝑗
for 𝑖 = 0,
i > 0
(1)
𝐿𝑖𝑗 =
⎧
⎩
⎨
⎪
⎪
⎪
⎪
𝐴𝑖𝑗
𝑈𝑗𝑗
𝐴𝑖𝑗 − ∑
𝑗−1
𝑘=0 𝐿𝑖𝑘𝑈𝑘𝑗
𝑈𝑗𝑗
for 𝑗 = 0,
j > 0
(2)
Deliverable 1 [10%]
Perform a simplified LU decomposition: Use equations (1) and (2) to complete
LU_Decomposition in the base code. The routine takes an numpy.array and returns the
lower and upper triangular matrices and .
☑
(𝑛, 𝑛) 𝐀
𝐋 𝐔
𝐋𝐲 = 𝐛 𝐲
𝑦𝑖 =
⎧
⎩
⎨
⎪
⎪
⎪
⎪
𝑏1
𝐿11
(
− )
1
𝐿𝑖𝑖
𝑏𝑖 ∑
𝑗=1
𝑖−1
𝐿𝑖𝑗𝑦𝑗
for 𝑖 = 1,
otherwise.
(3)
Deliverable 2 [10%]
Perform forward substitution: Use equation (3) to complete ForwardSubstitution in the base
code. The routine takes a lower triangular numpy.array and a numpy.array ; it
returns an numpy.array .
☑
(𝑛, 𝑛) 𝐋 (𝑛,) 𝐛
(𝑛,) 𝐲
𝐔𝐱 = 𝐲 𝐱
𝑥𝑖 =
⎧
⎩
⎨
⎪
⎪
⎪
⎪
𝑦𝑛
𝑈𝑛𝑛
(
− )
1
𝑈𝑖𝑖
𝑦𝑖
𝑗
∑=𝑖+1
𝑛
𝑈𝑖𝑗𝑥𝑗
for 𝑖 = 𝑛,
otherwise.
(4)
You now have all the pieces needed to piece together a linear system solver. A solution to the general system
can be obtained by solving the two simplified systems:
2 Solving Polynomial Regression
Let’s test out your linear solver on a few polynomial regression problems. Here, you’ll formulate polynomial
regression problems as solutions to linear systems of equations, before applying your solver to them.
2.1 Polynomial Regression System
Given data points we want to find the degree polynomial
that best fits the data points and where the coefficient vector fully defines the polynomial.
We can formulate this problem as the solution to a linear system of equations
where the Vandermonde matrix can be formed using the independent variables of our data points, .
Deliverable 3 [10%]
Perform backward substitution: Use equation (4) to complete BackwardSubstitution in the
base code. The routine takes an upper triangular numpy.array and a numpy.array
; it returns an numpy.array .
☑
(𝑛, 𝑛) 𝐔 (𝑛,)
𝐲 (𝑛,) 𝐱
𝐀𝐱 = 𝐋𝐔𝐱 = 𝐛
𝐋𝐲 = 𝐛 and 𝐔𝐱 = 𝐲 . (5)
Deliverable 4 [10%]
Put together a simple solver: Complete the routine LU_solver to solve the linear system
using your simplified LU decomposition, forward and backward substitution.
☑
𝐀𝐱 = 𝐛
𝑚 (𝑡 , ) 𝑖 𝑏𝑖 𝑛 − 1
𝑝𝑛−1(𝑡) = ∑
𝑗=1
𝑛
𝑥𝑗 𝑡
𝑗−1
𝐱 = [𝑥 , . . . , ] 1 𝑥𝑛
𝐀𝐱 = = = 𝐛
⎡
⎣
⎢
⎢
⎢
⎢
⎢
1
1
⋮
1
𝑡1
𝑡2
⋮
𝑡𝑚
⋯
⋯
⋱
⋯
𝑡
𝑛−1
1
𝑡
𝑛−1
2
⋮
𝑡
𝑛−1
𝑚
⎤
⎦
⎥
⎥
⎥
⎥
⎥
⎡
⎣
⎢
⎢
⎢
⎢
𝑥1
𝑥2
⋮
𝑥𝑛
⎤
⎦
⎥
⎥
⎥
⎥
⎡
⎣
⎢
⎢
⎢
⎢
𝑏1
𝑏2
⋮
𝑏𝑚
⎤
⎦
⎥
⎥
⎥
⎥
𝐀 𝑡𝑖
2.2 Fully-constrained and Overdetermined Polynomial Fitting
If the number of data points matches the number of unknowns , we can perfectly solve .
Consequently, the Vandermonde matrix here would be square.
If, on the other hand, we have more data points than degrees of freedom for our polynomial, we have an
overdetermined problem. A perfect fit will not generally exist, but we can solve for the fit that minimizes the
squared residual -norm . The normal equations allow us to express the least-squares solution
using the modified system .
3 Image Reconstruction using Deconvolution
Most real-world problems are sufficiently complex to require more robust solvers than what we developed
above: without pivoting and careful performance-minded vectorization/implementation, your LU solver won’t
likely scale to larger problems.
Setting up a problem and using a solver is just as important as know how to write one. The final set of
deliverables will focus on understanding a more complex problem, formulating its solution as a linear system,
and solving it with an industry-caliber LU solver.
3.1 Image Filtering — a motivating example
Imagine taking a photo with your smartphone only to realize, after the fact, that the photo came out blurry.
Luckily, your phone’s accelerometer was able to register a horizontal motion at the time of capture.
Deliverable 5 [10%]
Given an numpy.array and positive integer , complete the implementation
of GetVandermonde to construct and return an numpy.array Vandermonde matrix for the
degree polynomial.
☑
(𝑚,) 𝐭 = [𝑡1, . . . , 𝑡𝑚] 𝑛
(𝑚, 𝑛) 𝐀
𝑛 − 1
𝑚 𝑛 𝐀𝐱 = 𝐛
𝐀
𝑚
2 ||𝐀𝐱 − 𝐛||
2
2
𝐀 𝐀𝐱 = 𝐛
𝐓 𝐀
𝐓
Deliverable 6 [20%]
Solve the polynomial regression problem: Complete the implementation of PolynomialFit. It
takes as input an numpy.array of the data point and a positive integer to
denote the polynomial degree . You should use your GetVandermonde and LU_solver
routines.
☑
(𝑚, 2) 𝑚 (𝑡𝑖, 𝑏𝑖) 𝑛
(𝑛 − 1)
We could model the forward process that polluted the image as a linear operator. Specifically, by convolving
the discretized (1D) horizontal blur kernel (that the accelerometer registered) along the horizontal axis of the
uncorrupted image, we can arrive at the blurry image.
The blur kernel, which we’ll visualize as a few dots ●●●, has an odd number of elements, each with a
numerical weight. You can imagine centering the kernel atop every pixel in the original uncorrupted image,
weighting each pixel it overlaps with by the corresponding kernel value, and summing accross the weighted
pixel values in order to obtain a corrupted image pixel value (⬣, below). When any part of the kernel falls
outside the bounds of the image, it “loops over” to the opposing end of the image.
For example below, to compute the ⬣ pixel value in the corrupted image on the right, we weight the
intensities of solid pixels from the uncorrupted image on the left by the one overlapping element (of three
elements, in this example) of the blur kernel.
⬚ ⬚ ⬚ ⬚ ⬚
⬚ ⬚ ⬚ ⬚ ⬚
⬚ ⬚ ⬣ ⬚ ⬚
⬚ ⬚ ⬚ ⬚ ⬚
⬚ ⬚ ⬚ ⬚ ⬚
By “sliding” the kernel, and repeating this weighted sum, along each pixel of each row of the uncorrupted
image, we construct the final blurred image:
⬣ ⬚ ⬚ ⬚ ⬚ ⬚ ⬣ ⬚ ⬚ ⬚ ⬚ ⬚ ⬣ ⬚ ⬚
⬚ ⬚ ⬚ ⬚ ⬚ ⬚ ⬚ ⬚ ⬚ ⬚ ⬚ ⬚ ⬚ ⬚ ⬚
⬚ ⬚ ⬚ ⬚ ⬚ , ⬚ ⬚ ⬚ ⬚ ⬚ , ⬚ ⬚ ⬚ ⬚ ⬚ ,
⬚ ⬚ ⬚ ⬚ ⬚ ⬚ ⬚ ⬚ ⬚ ⬚ ⬚ ⬚ ⬚ ⬚ ⬚
⬚ ⬚ ⬚ ⬚ ⬚ ⬚ ⬚ ⬚ ⬚ ⬚ ⬚ ⬚ ⬚ ⬚ ⬚
⬚ ⬚ ⬚ ⬚ ⬣ ⬚ ⬚ ⬚ ⬚ ⬚
⬚ ⬚ ⬚ ⬚ ⬚ ⬚ ⬚ ⬚ ⬚ ⬚
. . . , ⬚ ⬚ ⬚ ⬚ ⬚ , . . . , ⬚ ⬚ ⬚ ⬚ ⬚
⬚ ⬚ ⬚ ⬚ ⬚ ⬚ ⬚ ⬚ ⬚ ⬚
⬚ ⬚ ⬚ ⬚ ⬚ ⬚ ⬚ ⬚ ⬚ ⬣
3.2 Blurring an Image
If we flatten the uncorrupted input image row-wise into a single column vector , where is the
number of image pixels, then we can model this linear corruption process with a matrix that relies only on
the blur kernel elements.
𝑛 × 1 𝐱 𝑛
𝐀
𝑘 × 1
Once you form , you can obtain the linearized corrupted image (of size ) as .
3.3 Deblurring an Image
We’re lucky that, in our setting, the accelorometer was able to provide us with an estimate of the blurring
kernel. Given this, and the forward model of how the uncorrupted image was corrupted with the kernel, we
can solve the inverse problem: given only the corrupted image and the blur kernel, we aim to retrieve the
uncorrupted image. This problem is referred to a non-blind deconvolution; if we were not given the blur
kernel, the (much harder) problem is referred to as blind deconvolution.
Here, we can retriev uncorrupted image given the blur matrix derived from the 1D horizontal blur kernel
and the corrupted (blurred) image by solving for in .
Deliverable 7 [20%]
Construct the blur matrix:Complete the implementation of CreateBlurMatrix. It takes as input
a numpy.array with the 1D horizontal blur kernel elements, as well as two positive integers
denoting the width and height of the images — we assume that the corrupted and uncorrupted
images have the same dimensions.
☑
(𝑘, 1)
𝐀 𝐛 𝑛 × 1 𝐀 𝐱
Deliverable 8 [5%]
Blur an image:Complete the implementation of BlurImage. It takes as input a numpy.array
of the blur matrix computed with CreateBlurMatrix and a width,height numpy.array of the
uncorrupted grayscale image. It should output a width, height numpy.array of the corrupted
grayscale image.
☑
(𝑛, 𝑛)
( )
( )
𝐱 𝐀
𝐛 𝐱 𝐀𝐱 = 𝐛
Deliverable 9 [5%]
Deblur an image:Complete the implementation of DeblurImage. It takes as input a
numpy.array of the blur matrix computed with CreateBlurMatrix and a width,height
numpy.array of the corrupted grayscale image. It should output a width, height numpy.array
of the uncorrupted grayscale image. Note: Your LU solver will not scale to the sizes of images
we will be using. You should use scipy’s LU solving routines, which we have imported for
you.
☑
(𝑛, 𝑛)
( )
( )
The test images we provide you for deliverables seven through nine. Your code will generate the two
missing images, which should match those of the upper row.
4 You’re Done!
Congratulations, you’ve completed the 1
st assignment. Review the submission procedures and guidelines at
the start of the Assignment 0 handout before submitting your Python script file assignment solution.
formatted by Markdeep 1.13
✒

ECSE 343 Assignment 2: Advanced Model Fitting

The learning objectives of this assignment are to
explore and implement a:
maximum likelihood estimator,
maximum a posteriori estimator, and
the expectation-maximization algorithm.
In doing so, you’ll revisit the polynomial fitting problem before moving on to a more sophisticated mixture model
classification problem.
1 Polynomial regression with MLE and MAP [40%]
We begin by exploring a maximum likelihood estimator (MLE) and maximum a posteriori (MAP) estimator on the
familiar polynomial regression problem.
Comparing degree-7 polynomial fits to noisy data and across estimators.
1.1 Polynomial regression
The overdetermined degree- polynomial regression problem — with an explicit corruption/noise model on the
data — seeks an interpolant across data samples that satisfies:
where:
is the independent coordinate of sample , with ,
is the dependent coordinate of sample , with ,
is standard Gaussian noise corrupting the outputs , and
are the polynomial coefficients we are solving for.
Note that one common way to rewrite this model is by “folding in” the deterministic component into the mean of
the Gaussian, as:
1.2 Maximum likelihood estimation (MLE) [20%]
You will implement a MLE for the polynomial parameters that maximize the data’s likelihood function:
where — assuming the samples are drawn i.i.d. — the likelihood function can be expressed using the
normal distribution’s density, as
Taking the log of the likelihood before taking the argmax — which is a valid transformation under argmax, given
log’s monotonicity — yields:
which we optimize for by setting the derivatives w.r.t. the parameters to zero,
before combining these partial derivatives to form a linear system , with
𝑚
𝑛 (𝑥𝑖, 𝑦𝑖)
𝑦𝑖 = ∑ + (𝑖 = 1,…, 𝑛)
𝑗=0
𝑚
𝜃𝑗𝑥
𝑗
𝑖 𝜖𝑖
𝑥𝑖 𝑖 𝐱 = {𝑥0 …𝑥𝑛−1}
𝑦𝑖 𝑖 𝐲 = {𝑦0 …𝑦𝑛−1}
𝜖𝑖 ∼ (𝜇, 𝜎 )
2 𝑦𝑖
𝜽 = {𝜃0, 𝜃1,…, 𝜃𝑚} 𝑝 = 𝑚 + 1
∼ 
(
𝜇 + , )
𝑦 . 𝑖 ∑
𝑗=0
𝑚
𝜃𝑗𝑥
𝑗
𝑖 𝜎
2
𝜽
𝜽MLE := argmax 𝑝(𝐲|𝐱, 𝜽) ,
𝜽
(𝑥𝑖, 𝑦𝑖)
𝑝(𝐲|𝐱, 𝜽) = ∏𝑝( | , 𝜽) = exp − .
𝑖=0
𝑛−1
𝑦𝑖 𝑥𝑖 ∏
𝑖=0
𝑛−1
1
2𝜋𝜎
2 √‾‾‾‾‾
⎛
⎝
⎜
⎜
1
2𝜎
2 [
− 𝜇 −
]
𝑦𝑖 ∑
𝑗=0
𝑚
𝜃𝑗𝑥
𝑗
𝑖
2⎞
⎠
⎟
⎟
log(𝑝(𝐲|𝐱, 𝜽)) = log(
𝑝( | , 𝜽)) ∏ = −𝑛 log( )) −
𝑖=0
𝑛−1
𝑦𝑖 𝑥𝑖 2𝜋𝜎
2 √‾‾‾‾‾
1
2𝜎
2 ∑
𝑖=0
𝑛−1
[
− 𝜇 −
]
𝑦𝑖 ∑
𝑗=0
𝑚
𝜃𝑗𝑥
𝑗
𝑖
2
log(𝑝(𝐲|𝐱, 𝜽)) =
[
− 𝜇 −
]
= 0 ,
∂
∂𝜃𝑘
1
𝜎
2 ∑
𝑖=0
𝑛−1
𝑥
𝑘
𝑖 𝑦𝑖 ∑
𝑗=0
𝑚
𝜃𝑗𝑥
𝑗
𝑖
𝐀 𝜽MLE = 𝐛
Solving this system yields the MLE for .
1.3 Maximum a posteriori (MAP) [20%]
Maximum a posteriori (MAP) estimators consider an additional prior over the parameters and solve for
parameters that maximize the (log) posterior distribution of the parameters, given the data:
In the polynomial regression setting, we will assume that the are drawn i.i.d from a normal distribution
, so their joint prior probability density is given by:
The derivative with respect to each parameter of is
and we combine these partials together to form a linear system that we can solve for the MAP
estimate of the parameters, where
𝐴𝑖,𝑗 = ∑ and = ( − 𝜇) ; 𝑖,𝑗 ∈ [0,𝑀]
𝑠=0
𝑛−1
𝑥
𝑖+𝑗
𝑠 𝑏𝑖 ∑𝑠=0
𝑛−1
𝑥
𝑖
𝑠 𝑦𝑠
𝜽MLE
Deliverable 1 [20%]
Complete the implementation of the function PolRegMLE that takes inputs (equal to ), , ,
and and outputs the system matrix and RHS for the MLE. Notice how the system does not
depend on , as discussed in class.
☑
𝑝 𝑚 + 1 𝜇 𝐱
𝐲 𝐀 𝐛
𝜎
Unlike the first part of A1, deliverables 1, 2 and 3 below won’t solve the linear system, but only form it.
You can explore numpy’s many linear system solving routines (e.g, numpy.linalg.solve) in order to
compute and assess the solutions of your various estimator systems.
ⓘ
𝜽∗
𝜽MAP
𝜽MA := 𝑝(𝜽|𝐲, 𝐱) = 𝑝(𝐲|𝐱, 𝜽) 𝑝(𝜽) . P argmax
𝜽
argmax
𝜽
𝜽
(𝜇 , ) 𝜽 𝜎
2
𝜽
𝑝(𝜽) = 𝑝( ) = exp (
− ( − ) ∏ .
𝑘=0
𝑚
𝜃𝑘 ∏
𝑘=0
𝑚
1
2𝜋𝜎
2
𝜽
‾‾‾‾‾ √
1
2𝜎
2
𝜽
𝜃𝑘 𝜇𝜽)
2
𝜃𝑘
log(𝑝(𝐲|𝐱, 𝜽) 𝑝(𝜽)) ≡ 𝑙(𝜽|𝐱, 𝐲)
𝑙(𝜽|𝐱, 𝐲) =
[
− 𝜇 −
]
− ( − ) ,
∂
∂𝜃𝑘
1
𝜎
2 ∑
𝑖=0
𝑛−1
𝑥
𝑘
𝑖 𝑦𝑖 ∑
𝑗=0
𝑚
𝜃𝑗𝑥
𝑗
𝑖
1
𝜎
2
𝜽
𝜃𝑘 𝜇𝜽
𝐀𝜽MAP = 𝐛
𝐴𝑖,𝑗 = and = ( − 𝜇) + ; 𝑖,𝑗 ∈ [0, 𝑚]
⎧
⎩
⎨
⎪
⎪
⎪
⎪
∑ +
𝑠=0
𝑛−1
𝑥
𝑖+𝑗
𝑠
𝜎
2
𝜎
2
𝜽
∑𝑠=0
𝑛−1
𝑥
𝑖+𝑗
𝑠
if 𝑖 = 𝑗
otherwise
𝑏𝑖 ∑𝑠=0
𝑛−1
𝑥
𝑖
𝑠 𝑦𝑠 𝜇𝜽
𝜎
2
𝜎
2
𝜽
Deliverable 2 [15%]
Complete the implementation of the function PolRegMAP with inputs , , , , , , and
and that outputs the MAP linear system matrix and RHS .
☑
𝑝 = 𝑚 + 1 𝜇 𝜎 𝜇𝜽 𝜎𝜽 𝐱
𝐲 𝐀 𝐛
As discussed in class, the linear least squares solution to the polynomial regression problem is equivalent to the
MLE when the noise is assumed to be normally distributed with zero mean. As such, we should expect that MLE
and LLS yield the same parameters — with the exception of the constant term — when the noise is not of zero
mean.
2 Mixture Models with Expectation-Maximization [60%]
We will now explore the expectation-maximization (EM) algorithm through the lens of a Gaussian mixture model
(GMM) example.
Fitting a GMM to a randomly generated, multi-modal 2D dataset using the EM algorithm.
Consider many -dimensional data samples drawn from a mixture model parameterized by .
Our example will use a mixture of multivariate normal distributions (i.e., high-dimensional, anisotropic
Gaussians), with mixture components, and so our , where
each of the -dimensional Gaussians are parameterized by a -dimensional mean vector , a
covariance matrix (see below), and
the normalized mixture proportions weight the amount of each of the components in the mixture, with
.
We use the standard definition of a -dimensional multivariate normal distribution, with density:
𝜃0
Deliverable 3 [5%]
Revisit LLS by completing the implementation of PolRegLLS with inputs , , and and
that outputs the appropriate fully-constrained linear system matrix and RHS needed to construct
the LLS system.
☑
𝑝 = 𝑚 + 1 𝑛 𝐱 𝐲
𝐀 𝐛
Our evaluation of your GMM code will impose a (sufficiently lenient) time limit on the execution your
algorithm — be mindful to avoid loops and to leverage numpy vectorization whenever suitable. ⓘ
𝑑 𝐱𝑖 𝜽
(𝐱|𝝁,𝚺)
𝐾 𝜽 = {𝝁0,…, 𝝁𝐾−1,𝚺0,…,𝚺𝐾−1, 𝜋0,…, 𝜋𝐾−1}
𝐾 𝑑 𝑑 𝝁𝑖 𝑑 × 𝑑
𝚺𝑖
𝜋𝑐 𝐾
∑ = 1
𝐾−1
𝑐=0 𝜋𝑐
𝑑
(𝐱|𝝁,𝚺) = exp[− (𝐱 − 𝝁 (𝐱 − 𝝁)]
,
1
(2𝜋)
𝑑/2
|𝚺|
1/2
1
2
)
⊺𝚺
−1
and so the marginal distribution of is:
As such, the joint probability of i.i.d. data samples is
And so, given the samples and knowing only the number of mixture components , our goal is to estimate the
mixture’s parameters that best match the data: the means (i.e., centers), covariances (i.e., anisotropic
“shapes”), and weights of the mixture components.
2.1 From MLE to Expectation-Maximization
A natural first step is to derive an MLE for . The likelihood function is simply the joint probability above,
and so the log-likelihood is
We immediately run into difficulties when differentiating with respect to, say, :
as we are unable to isolate the : the non linearity induced by the
term. We face similar problems when differentiating with respect
to other parameters in .
To sidestep this issue, we will introduce the concept of a latent label for each data sample :
this latent label corresponds to the ground truth association of each data sample to its corresponding mixture
component, and so .
By the law of total probability, the marginal conditional probability distribution of given can now be rewritten as
Intuitively, knowing the latent variables should help us in computing the MLE. To do so, we first compute the
vector posterior with elements using Bayes’ theorem:
𝐱𝑖
𝑝(𝐱 |𝜽) = ( | , ) . 𝑖 ∑𝑐=0
𝐾−1
𝜋𝑐 𝐱𝑖 𝝁𝑐 𝚺𝑐
𝑛 ˆ𝐱 = {𝐱0,…, 𝐱𝑛−1}
𝑝(ˆ𝐱|𝜽) = ∏ ( | , ) .
𝑖=0
𝑛−1
∑𝑐=0
𝐾−1
𝜋𝑐 𝐱𝑖 𝝁𝑐 𝚺𝑐
ˆ𝐱 𝐾
𝜽 𝝁𝑐 𝚺𝑐
𝜋𝑐
𝜽 𝑝(ˆ𝐱|𝜽)
log 𝑝( |𝜽) = log(
( | , ))
ˆ𝐱 ∑ .
𝑖=0
𝑛−1
∑𝑐=0
𝐾−1
𝜋𝑐 𝐱𝑖 𝝁𝑐 𝚺𝑐
𝝁𝑘
log 𝑝( |𝜽) = ( − ) = 0 ,
∂
∂𝝁𝑗
ˆ𝐱 ∑
𝑖=0
𝑛−1
⎛
⎝
⎜
⎜
⎜
⎜
𝜋𝑗(𝐱𝑖|𝝁𝑗
,𝚺𝑗)
∑ ( | , )
𝑐=0
𝐾−1
𝜋𝑐 𝐱𝑖 𝝁𝑐 𝚺𝑐
⎞
⎠
⎟
⎟
⎟
⎟
𝐱𝑖 𝝁𝑗 𝚺
−1
𝑗
𝝁𝑘
( ( | , )) / ( ( | , ))
𝜋𝑗 𝐱𝑖 𝝁𝑗 𝚺𝑗 ∑
𝑐=0
𝐾−1
𝜋𝑐 𝐱𝑖 𝝁𝑐 𝚺𝑐
𝜽
𝑧𝑖 ∈ {1,…, 𝐾} 𝐱𝑖
𝜋𝑐 = 𝑝(𝑧𝑖 = 𝑐|𝜽)
𝐱𝑖 𝜽
𝑝(𝐱 |𝜽) = . 𝑖 ∑𝑐=0
𝐾−1
𝑝(𝑧𝑖 = 𝑐|𝜽)

𝜋𝑐
𝑝(𝐱𝑖|𝑧𝑖 = 𝑐, 𝜽)

(𝐱𝑖|𝝁 , ) 𝑐 𝚺𝑐
𝑧𝑖
𝜶 𝛼𝑖,𝑗 = 𝑝(𝑧𝑖 = 𝑘|𝐱𝑖, 𝜽)
We can now rewrite the partial derivative of the log-likelihood with respect to and set it to zero as:
Even though depends on , we can define an iterative process that assumes we have an approximate value
of and solving for as:
where can be interpreted as the effective number of data samples assigned to the mixture
component . We can similarly obtain expressions for the other parameters as:
After computing these (updated) values for the , and , we can recompute new estimates for from the
expression we derived using Bayes’ theorem. We repeat this process until “convergence”.
2.2 Formalizing the EM Algorithm
To summarize:
if we knew the parameters , we could compute the posterior probabilities , and
if we knew the posteriors , we could compute the parameters .
The EM algorithm proceeds as follows:
0. Initialize the parameters and evaluate the log-likelihood with these parameters,
1. E-step: evaluate the posterior probabilities using the current values of the parameters ,
2. M-step: estimate new parameters with the current values of ,
3. Evaluate the log-likelihood with the new parameter estimates — if the log-likelihood has
changed by less than some (typically small) convergence threshold, terminate; otherwise, repeat steps 1
through 3.
𝛼 := 𝑝( = 𝑗| , 𝜽) = = . 𝑖,𝑗 𝑧𝑖 𝐱𝑖
𝑝(𝑧𝑖 = 𝑗|𝜽)𝑝(𝐱𝑖|𝑧𝑖 = 𝑗, 𝜽)
𝑝(𝐱𝑖|𝜽)
𝜋𝑗(𝐱𝑖|𝝁𝑗
,𝚺𝑗)
∑ ( | , )
𝑐=0
𝐾−1
𝜋𝑐 𝐱𝑖 𝝁𝑐 𝚺𝑐
𝝁𝑗
log 𝑝( |𝜽) = ( − ) = 0 .
∂
∂𝝁𝑗
ˆ𝐱 ∑
𝑖=0
𝑛−1
𝛼𝑖,𝑗 𝐱𝑖 𝝁𝑗 𝚺
−1
𝑗
𝛼𝑖,𝑗 𝝁𝑗
𝛼𝑖,𝑗 𝝁𝑗
𝝁𝑗 =
1
𝜔𝑗 ∑
𝑖=0
𝑛−1
𝛼𝑖,𝑗 𝐱𝑖
𝜔𝑗 = ∑
𝑛−1
𝑖=0 𝛼𝑖,𝑗
𝑗
𝚺 = ( − ( − ) and = . 𝑗
1
𝜔𝑗 ∑
𝑖=0
𝑛−1
𝛼𝑖,𝑗 𝐱𝑖 𝝁𝑗)
⊺ 𝐱𝑖 𝝁𝑗 𝜋𝑗
𝜔𝑗
𝑛
𝝁𝑗 𝚺𝑗 𝜋𝑗 𝛼𝑖,𝑗
𝜽 𝜶
𝜶 𝜽
𝜽 log 𝑝(ˆ𝐱|𝜽)
𝛼𝑖,𝑗 𝜽
𝜽 𝛼𝑖,𝑗
log 𝑝(ˆ𝐱|𝜽)
EM is sensitive to the initial parameter values, so special care must be taken in the first step. Given
“valid” initial values, however, the EM algorithm guarantees that the log-likelihood increases after
every iteration.
ⓘ
We provide a basic class structure to encapsulate GMM parameters and some utility methods.
2.3 EM Implementation
Compute normalization of the mixture for all the samples
Compute the mixture normalization factors , given the current parameters , for all the data
samples, as:
You will use the result of this function in the expectation step and log-likelihood computation.
Expectation step
Using the values of , you have to compute the posterior probabilities with elements for the mixture
components, as:
Maximization step
Using the values of you have to update the following quantities:
During EM iterations with GMMs, covariance matrices may become singular. To circumvent this
problem, consider regularizing the covariance matrices, when necessary; note, however, that using too
large a regularization can lead to divergence (e.g., a possible reduction in the log-likelihood estimate
between iterations.)
ⓘ
When completing your implementation, you should only require explicit loops over the mixture
component. All other operations can be vectorized with numpy routines, such as sum( , axis=), dot(
, axis=), argmax( , axis=) and the component-wise array operators +, – , * and /. We similarly
provide you with access to scipy’s multivariate_normal class whose (efficient and vectorized) pdf
method implements the Gaussian density evaluation.
ⓘ
Carefully review the GMM class we provide in the base code.
𝜷 = {𝛽0,…, 𝛽𝑛−1} 𝜽
𝛽 := ( | , ) . 𝑖 ∑𝑐=0
𝐾−1
𝜋𝑐 𝐱𝑖 𝝁𝑐 𝚺𝑐
Deliverable 4 [10%]
Complete the implementation of normalization which takes as input the GMM instance and updates
its property gmm.beta.
☑
𝜷
𝛽𝑖 𝜷 𝜶 𝛼𝑖,𝑗
𝛼𝑖,𝑗 =
𝜋𝑗(𝐱𝑖|𝝁𝑗
,𝚺𝑗)
𝛽𝑖
Deliverable 5 [10%]
Complete the implementation of expectation which takes as input the GMM instance and updates its
property gmm.alpha.
☑
𝜶
𝛼𝑖,𝑗 𝜶
Computing the log-likelihood
We can finally compute the log-likelihood estimate using the current :
formatted by Markdeep 1.13
𝑧𝑖 = argmax = = =
𝑗
𝛼𝑖,𝑗 𝜔𝑗 ∑
𝑖=0
𝑛−1
𝛼𝑖,𝑗 𝜋𝑗
𝜔𝑗
𝑛
𝝁𝑗
1
𝜔𝑗 ∑
𝑖=0
𝑛−1
𝛼𝑖,𝑗 𝐱𝑖
𝚺𝑗 = ( − ( − ) +
1
𝜔𝑗 ∑
𝑖=0
𝑛−1
𝛼𝑖,𝑗 𝐱𝑖 𝝁𝑗)
⊺ 𝐱𝑖 𝝁𝑗 𝚺regularization
Deliverable 6 [20%]
Complete the implementation of maximization which updates the , , , , in the input GMM
instance; the associated member variables are gmm.Z, gmm.weight, gmm.pi, gmm.mu, gmm.cov.
Afterwards, you can update the multivariate normal parameters of the mixture elements in gmm.gauss
with and .
☑
𝑧𝑖 𝜔𝑗 𝜋𝑗 𝝁𝑗 𝚺𝑗
𝝁𝑗 𝚺𝑗
ⓘ Don’t forget to regularize the covariance matrices.
𝜷
log 𝑃(ˆ𝐱|𝜽) = ∑log( )
𝑖=0
𝑛−1
𝛽𝑖
Deliverable 7 [10%]
Complete the implementation of logLikelihood which computes the log-likelihood
returning nothing and updating the GMM instance variable gmm.log_likelihoods.
☑
log 𝑃(ˆ𝐱|𝜽)
✒

ECSE 343 Assignment 3: 2D Deconvolution, Regularization and SVD

1 Revisiting Image Convolution & Deconvolution
Let’s revisit and generalize the image convolution and deconvolution problems we saw in Assignment 1. This
time, instead of only consider a discrete 1D horizontal blur kernel, we will consider more general 2D blurs
(i.e., horizontal and vertical).
To illustrate the forward blurring process we extend our previous diagrams, however the “sweep and stamp”
process still applies. Our diagram below illustrates a few instances of applying a blur kernel (i.e., with 9
elements, which we’ve numbered for clarity: ❶, ❷, , ❾) to a input image, and generating a
output image, e.g.,
⬚ ⬚ ⬚ ⬚ ⬚
❶❷❸ ⬚ ⬚ ⬚ ⬚ ⬚
❹❺❻ ⬚ ⬚ ⬚ ⬚
❼❽❾ ⬚ ⬚ ⬚ ⬚ ⬚
⬚ ⬚ ⬚ ⬚ ⬚
3 × 3
… 5 × 5 5 × 5
You can assume that we will only test square kernels, input/output image resolutions greater than
or equal to the kernel size, and only square images. That being said, writing your code in a manner
that doesn’t assume square kernels can help with debugging/testing.
ꋏ
The kernel is centered — this time horizontally and vertically — on an input pixel at the location of the
corresponding output pixel (illustrated as a circle with black contour and grey fill), and if portions of the
kernel happen to fall outside the bounds of the input image, they wrap around (vertically and horizontally,
now).
By “sliding” the kernel along each pixel of the uncorrupted input image, weighting the underlying
source/input pixels by the overlapping kernel values and summing across these weighted input pixels, we
construct the final blurred image one pixel at a time:
❺❻ ❹ ⬚ ⬚ ⬚ ⬚ ❹❺❻ ⬚ ⬚ ⬚ ⬚ ❹❺❻ ⬚ ⬚ ⬚ ⬚
❽❾ ❼ ⬚ ⬚ ⬚ ⬚ ⬚ ❼❽❾ ⬚ ⬚ ⬚ ⬚ ⬚ ❼❽❾ ⬚ ⬚ ⬚ ⬚ ⬚
⬚ ⬚ ⬚ ⬚ ⬚ , ⬚ ⬚ ⬚ ⬚ ⬚ , ⬚ ⬚ ⬚ ⬚ ⬚ , ⬚ ⬚ ⬚ ⬚ ⬚ ⬚ ⬚ ⬚ ⬚ ⬚ ⬚ ⬚ ⬚ ⬚ ⬚
❷❸ ❶ ⬚ ⬚ ⬚ ⬚ ⬚ ❶❷❸ ⬚ ⬚ ⬚ ⬚ ⬚ ❶❷❸ ⬚ ⬚ ⬚ ⬚ ⬚
❻ ❹❺ ⬚ ⬚ ⬚ ⬚ ❾ ❼❽ ⬚ ⬚ ⬚ ⬚ ⬚
❾ ❼❽ ⬚ ⬚ ⬚ ⬚ ⬚ ⬚ ⬚ ⬚ ⬚ ⬚ . . . , ⬚ ⬚ ⬚ ⬚ ⬚ , . . . , ⬚ ⬚ ⬚ ⬚ ⬚
⬚ ⬚ ⬚ ⬚ ⬚ ❸ ❶❷ ⬚ ⬚ ⬚ ⬚ ⬚
❸ ❶❷ ⬚ ⬚ ⬚ ⬚ ⬚ ❻ ❹❺ ⬚ ⬚ ⬚ ⬚
Recall, as in Assignment 1, while the diagrams above illustrate the blurring process in 2D and for 2D images
and kernels, you will express this linear map as a matrix — and thus — must treat the input (and output)
images as “flattened” 1D vectors. Treating this change of indexing from 2D pixel coordinates to flattened 1D
“image” vector elements is one part of the learning objective for this task.
1.1 Isotropic Gaussian Blur Kernel
One common blur kernel is the discrete isotropic Gaussian kernel. For one-dimensional blurs, the continuous
kernel is defined as
where is the standard deviation of the Gaussian, which defines its fall-off rate, and the continuous kernel is
defined over the entire real line, i.e., for . In practice, one way to take this continuous kernel
and discretize it to have length is to evaluate the kernel at integer lattice values .
For two-dimensional blur kernels, the continuous isotropic Gaussian is
Deliverable 1 [25%]
Construct the blur matrix: Complete the implementation of CreateBlurMatrix. It takes as input
a numpy.array with the 2D blur kernel elements, as well as two positive integers denoting
the width and height of the image the matrix will be applied to — you can assume that the
corrupted and uncorrupted images have the same dimensions, and that all images are square (i.e.,
width = height). The routine returns the blur matrix as a numpy.array with shape (width
height, width height).
☑
(푘, 푘)
×
×
퐺1D (푥; 휎) = 1
√‾2‾휋‾휎
푒
− 푥2
2휎 2 (1)
휎
−∞ ≤ 푥 ≤ ∞
푘 푥 ∈ {−⌊ ⌋, … , ⌊ ⌋} 푘
2
푘
2
and its discretized version is obtained by evaluating on a 2D integer lattice with
and . Note that can be decomposed as the
product of two 1D Gaussians, a fact you can optionally leverage to simplify your implementation.
1.2 Conditioning & Regularization
As in Assignment 1, we provide you with a benchmark image set. Here, in addition to providing the clean and
blurred images, we also provide a version of the blurred image polluted with subtle noise. We will use this
noisy blurred image to explore the conditioning of your blur operator.
In the current scenario, with a fully-constrained blur matrix and blurred image solving for the leastsquares unpolluted image as
amounts to solving the original linear system as . For many blurring kernels, the blur matrix can
become ill-conditioned, meaning the solution can be catastrophically sensitive to minor deviations in the
input .
We can employ regularization strategies to try to mitigate the poor conditioning of the system. One common
strategy is Tikhonov regularization, which aims instead to solve the modified least-squares problem in
where is the Tikhonov regularization parameter that allows you to control the degree of regularization. This
regularized least-squares problem can be expressed as the solution to the following augmented system of
linear equations:
퐺2D (푥, 푦; 휎) = 1
2휋휎2
푒
− 푥2+푦2
2휎 2 (2)
푘 × 푘 (푥, 푦)
푥 ∈ {−⌊푘/2⌋, … , ⌊푘/2⌋} 푦 ∈ {−⌊푘/2⌋, … , ⌊푘/2⌋} 퐺2D
Deliverable 2 [10%]
2D Gaussian Blur Kernel: Complete the implementation of Gaussian2D to generate the 2D
discretized Gaussian blur kernel. It takes as input the (odd) kernel extent (i.e., for a square 2D
kernel of size ) and the standard deviation for the isotropic Gaussian sigma. It outputs a
numpy.array of the 2D Gaussian evaluated at integer lattice locations.
☑
푘 푘
(푘, 푘) 휎
(푘, 푘)
Consider different ways that you can test your Gaussian2D implementation, as well as how it can
be used to test your implementation of CreateBlurMatrix.
ꋏ
퐀 퐛
퐱
min ‖퐛 − 퐀퐱 퐱 ‖2
2 (3)
퐱 = 퐀 퐛 −1
퐛
min (‖퐛 − 퐀퐱 + ‖퐱 ) = , 퐱 ‖2
2 휆2 ‖2
2 min
퐱 [ ] − [ ] 퐱
‖
‖
‖ 퐛
0
퐀
휆퐈
‖
‖
‖
2
2
(4)
휆
(퐀 퐀 + 퐈) 퐱 = 퐛 . T 휆2퐈
T 퐀T (5)
Tikhonov regularization is arguably the most basic regularization scheme but, even then, a judicious choice
of its parameter can yield much better results than a more naïve setting. One generalization of the
Tikhonov formulation,
admits an additional operator — a so-called regularization matrix — that can be designed to specialize the
regularization based on the needs of a specific application. In the case of image deblurring, some choices for
that can outperform vanilla Tikhonov regularization include local derivative and/or Laplacian operators.
2 Solving Linear Systems with Singular Value Decomposition
Let’s explore how singular value decompositions (SVDs) can be used to solve (implicitly regularized) linear
systems. As a warm up, we will implement a simple singular vector solver, before moving on to the
deconvolution problem using numpy’s built-in SVD algorithm. Recall that the SVD of is and has
non-zero singular values.
Deliverable 3 [10%]
Regularized deblurring: Complete the implementation of DeblurImage. It takes as input the
unregularized blur matrix, the (potentially noisy) blurred image, and the Tikhonov regularization
constant lambda_. The function should construct the augmented/regularized linear system of
equations and solve them (using whichever numpy routine you like) to solve for and output the
(unflattened, 2D) deblurred image.
☑
휆
• You may find it useful to test your routines on smaller images.
• You may want to complete the rest of the assignment before returning to complete the next
deliverable.
ꋏ
휆
min (‖퐛 − 퐀퐱 + ‖퐃퐱 ) = , 퐱 ‖2
2 휆2 ‖2
2 min
퐱 [ ] − [ ] 퐱
‖
‖
‖ 퐛
0
퐀
휆퐃
‖
‖
‖
2
2
(6)
퐃
퐃
Deliverable 4 [15%]
Deblurring competition: DeblurImage implements the basic Tikhonov regularization scheme,
however other regularization methods for deblurring exist. Complete the implementation of
DeblurImageCompetition to implement any set of regularization approaches with the goal of
generating your best possible deblurred image result — here, you have full latitude in how you go
about regularizing, augmenting, etc. the linear system in order to improve the deblurring result.
The function signature is the same as DeblurImage from A1 (i.e., the same as DeblurImage
above, but without lambda_; you should hardcode any regularization constants you use within
your DeblurImageCompetition function body). Your grade in this deliverable will be based on
how well your deblurring works compared to your colleagues’, and how well you document your
solution. Deblurring error will be computed as the difference between your deblurred image
and the ground truth deblurred (i.e., source) image.
☑
퐿2
퐀 퐀 = 퐔횺퐕T
푘 ≤ min(푚, 푛)
2.1 Power Iteration
We discussed how the Power Iteration Method could be used to compute the eigenvectors associated to the
largest and smallest eigenvalues in a non-singular matrix . We can extend this approach to compute the left
and right singular vectors associated to the largest and smallest singular values of , by applying the method
to two separate eigenvector settings:
by applying the power iteration method to you can compute the left singular vectors in
associated to the largest ( for ) or smallest ( for ) singular values, and
by applying the power iteration method to you can compute the left singular vectors in
associated to the largest ( for ) or smallest ( for ) singular values.
2.2 Rank-reduced Deconvolution
We can use the SVD to solve linear systems. Consider taking the SVD of your blur matrix as and
then solving for in the unregularized deblurring problem. Clearly if we use the full-rank
reconstruction of as we will run into conditioning-related instabilities when, e.g., the blurred image
is polluted by noise.
Instead of using the full rank SVD reconstruction of , you can implement a rank-reduced reconstruction
that retains only a subset (i.e., a percentage up to 100%) of the system’s “energy”. To do so whilst retaining
of the system’s energy, you can form an effective rank-reduced approximation of that only
retains those singular values above the -percentage threshold — and zeros out the remaining singular
values.
Specifically, your rank-reduced approximation of will include only those singular values that
satisfy this inequality:
Your rank-reduced has and can be used to solve
Here, we assume that (i.e., ), otherwise the degree to which is zero-padded would differ.
Keep in mind that, while implementing Equation as-is will not lead to incorrect results, there are several
avenues for optimizing its performance without sacrificing any accuracy — in fact, some optimizations may
improve the numerics of the solution.
퐀
퐀
퐀퐀T 퐔
퐮1 휎1 퐮푘 휎푘
퐀 퐀T 퐕
퐯1 휎1 퐯푘 휎푘
Deliverable 5 [15%]
Power Iteration: Complete the implementation of PowerMiniSVD to return the left and right
singular vectors associated to the largest singular value, and . The function takes an input
matrix and number of power iterations as input.
☑
퐮1 퐯1
퐀 = 퐔횺퐕T
퐱 = 퐀 퐛 −1
퐀 퐔횺퐕T
퐛
퐀
0 ≤ 푅 ≤ 1 퐀
푅
퐀 {휎1, … , 휎푟}
(∑ / ) ≤ 푅 . 푟
푖=1 휎2
푖 ∑푘
푗=1 휎2
푗 (7)
퐀̂= 퐔횺̂
퐕T 횺 = diag( , … , , 0, … , 0) ̂ 휎1 휎푟
퐱 = 퐀 퐛 = 퐕 퐛 . ̂
−1
횺̂
−1
퐔T (8)
푟 < 푘 푅 < 1 횺̂
8
3 You’re Done!
Congratulations, you’ve completed the 3rd assignment. Review the submission procedures and guidelines at
the start of the Assignment 0 handout before submitting your Python script file assignment solution.
formatted by Markdeep 1.13
Deliverable 6 [25%]
SVD Deblurring: Complete the implementation of DeblurImageSVD which takes as input the SVD
factors of a blur matrix (computed outside the function), the blurred (and potentially noisy) input
image, the amount of energy to retain in the rank-reduced reconstruction of the blur matrix, and
returns the deblurred image as in Equation .
☑
푅
8
Some notes:
• Computing the SVD is expensive. You may want to implement a caching scheme in your test
code using numpy.save and numpy.load
• Do your rank-reduced image deblurring results require additional explicit regularization? You
don’t have to actually answer this in your solution, but it’s worthwhile reflecting on.
ꋏ
✒

ECSE 343 Assignment 4: Discrete Fourier Transform

2 Discrete Fourier Transform
Given an input continuous function sampled at equally-spaced discrete points on an integer lattice
in the primal domain (e.g., time, space, etc.), we denote this discretized signal as
. Note that, in general, can be a complex-valued function (and, so too the sampled values
); however, our examples will only consider real-valued input functions .
The discrete Fourier transform (DFT) of the sampled signal is then given by:
which can be equivalently expressed as a matrix-vector operation (see lecture slides) as
You must not use any additional imports for your solution, other than the ones provided by
us.
Do not edit the imports at the top of the .py file. Doing so will result in an assignment grade of
0%.
☒
It is very difficult to successfully complete this assignment without writing your own tests in
__main__. We strongly recommend against relying exclusively on the example tests we provide.
Think about what kind of input/output behaviours you can write and test against.
ⓘ
𝑓 𝑁
𝑥 ∈ {0,…, 𝑁 − 1}
𝑓(𝑥) ≡ 𝑓[𝑥] 𝑓
𝑓[⋅] 𝑓 : ℝ → ℝ
𝑓[𝜔]
̂
[𝜔] = 𝑓[𝑘] exp( )
𝑓 ,
̂ ∑
𝑘=0
𝑁−1 −2𝜋ı𝑘𝜔
𝑁
(1)
ˆ𝐟 = 𝐌𝑁 𝐟 ,
where is a vector of all the sampled function values, encodes the
linear operator form of the DFT for a length- input, and is a vector of
the sampled DFT of . The (discrete) signal can be interpreted as a discrete sampling of the
continuous Fourier transform of an -periodic summation of — i.e., repeated
periodically.
Note that, regardless of whether the input signal is real- or complex-valued, the DFT always has
complex-valued elements: encodes the frequency amplitudes and their phases.
The inverse discrete Fourier transform (iDFT) allows us to transform back from this (discrete, complexvalued) sampled frequency spectrum to the original (possibly only real-valued) sampled signal in the primal
domain, as
Note the similarities in the mathematical definitions of the DFT and iDFT in Equations and , comprising of
a sign change in the complex exponential and different normalization constants. You will leverage these
similarities to avoid code duplication: your DFT and iDFT implementations will rely on an implementation of a
generalized (i.e., bidirectional) DFT.
Specifically, you will implement this generalized DFT formula for a sampled input ,
where we obtain the DFT by setting and (and and ), and the iDFT with
and (and and ). Your implementations will act on, and return, an entire
vector of sampled primal- (i.e., ) or frequency-domain (i.e., ) values.
As discussed above, the generalized DFT equation can be used to realize both a forward and inverse DFT.
𝐟 = (𝑓[0], 𝑓[1],…, 𝑓[𝑁 − 1]) 𝐌𝑁
𝑁 ˆ𝐟 = (𝑓[0], [1],…, [𝑁 − 1])
̂ 𝑓̂ 𝑓̂
𝑓 𝑓[𝑥]
̂
{𝑓𝑁 (𝑥)} 𝑁 𝑓(𝑥) 𝑓(𝑥)
𝑓[𝑥]
̂
Re(𝑓)
̂ Im(𝑓)
̂
𝑓[𝑥] = [𝑘] exp( )
.
1
𝑁 ∑
𝑘=0
𝑁−1
𝑓̂
2𝜋ı𝑘𝑥
𝑁
(2)
1 2
Υ[⋅]
Ψ[𝑎;s] = Υ[𝑏] exp(
s ) ∑ ,
𝑏=0
𝑁−1
2𝜋ı𝑏𝑎
𝑁
Υ = 𝑓 𝑠 = −1 𝑎 ≡ 𝜔 𝑏 ≡ 𝑥
Υ =
1
𝑁
𝑓̂ 𝑠 = 1 𝑎 ≡ 𝑥 𝑏 ≡ 𝜔
𝐟 ˆ𝐟
Deliverable 1 [17.5%]
Implement the Generalized DFT Routine: Complete the implementation of the generalized DFT
function DFT as specified by and with the default parameter setting of .
☑
Ψ s = −1
Υ
Deliverable 2 [2.5%]
Implement the inverse DFT Routine: Complete the implementation of the inverse DFT function
iDFT that relies on calling your generalized DFT routine, above.
☑
3 2D Discrete Fourier Transform
Let’s now consider a generalization to 2D, where we sample a continuous function on the
(square) 2D integer lattice leading to a discretized 2D signal at .
We can similarly define the sampled 2D DFT of as
which admits a similarly compact matrix-vector expression as
where we can arrange the sampled primal function values in a matrix with elements , with
encoding the linear operator form of the 2D DFT for a square length- sampled input, and
is the matrix of the sampled 2D DFT coefficients with .
Again, similarly to the 1D setting, DFT elements are always complex-valued with encoding
frequency amplitudes and their phases.
The 2D inverse DFT is now
and a similar extension of the generalized DFT routine applies in the 2D setting.
3.1 Leveraging Factorization
After little manipulation, one can observe that the 2D DFT and iDFT in Equations and can be factorized
into several applications of their 1D counterparts.
A careful numpy implementation should not rely on any for loops or related Python operations (e.g.,
loop comprehensions). You may also wish to leverage the fact that exp (ı𝑥) = cos 𝑥 + ısin 𝑥.
𝑓 : ℝ → ℝ
2
𝑓[𝑥, 𝑦] (𝑥, 𝑦) ∈ {0,…, 𝑁 − 1}
2
𝑓[ , ]
̂𝜔𝑥 𝜔𝑦 𝑓[𝑥, 𝑦]
[ , ] = 𝑓[ , ] exp ( )
𝑓 ,
̂𝜔𝑥 𝜔𝑦
𝑘
∑𝑥=0
𝑁−1
𝑘
∑𝑦=0
𝑁−1
𝑘𝑥 𝑘𝑦
−2𝜋ı(𝑘𝑥𝜔𝑥 + 𝑘𝑦𝜔𝑦)
𝑁2
(3)
ˆ𝐅 = 𝐌𝑁 𝐅 , ,𝑁
𝐅 (𝐅)𝑖,𝑗 = 𝑓[𝑖,𝑗]
𝐌𝑁,𝑁 𝑁 × 𝑁 ˆ𝐅
(ˆ𝐅 )𝑖,𝑗 = 𝑓[𝑖,𝑗]
̂
𝑓[⋅, ⋅]
̂ Re(𝑓)
̂
Im(𝑓)
̂
𝑓[𝑥, 𝑦] = [ , ] exp ( )
1
𝑁2
𝑘
∑𝑥=0
𝑁−1
𝑘
∑𝑦=0
𝑁−1
𝑓̂𝑘𝑥 𝑘𝑦
2𝜋ı(𝑘 𝑥 + 𝑦) 𝑥 𝑘𝑦
𝑁2
(4)
3 4
By applying a 1D DFT independently to each of the rows of , before applying the 1D DFT independently to
the columns of the (partially-)transformed , we arrive at the 2D DFT transform of . In matrix form, this
amounts to
You should leverage the factorization to accelerate your implementation.
4 Fourier for Image Convolution & Deconvolution
In Assignments 1 and 3, we explored convolution and deconvolution using linear operators that act on
sampled signals in the primal (i.e., pixel, spatial) domain.
We constructed the convolution operators as a function of the underlying convolution kernel, and we
implemented a wrap-based boundary condition in instances during convolution where kernels only partially
overlapped the input domain.
Manually treating the indexing, which additionally required “flattening” 2D images into 1D vectors, required
care.
Frequency-space representations are particularly well-suited to certain applications, among them
convolution. We will use the DFT to perform image convolution and (regularized) deconvolution, leveraging
these advantages:
1. convolution in the primal domain amounts to multiplication in the frequency domain, and
2. in the frequency domain, we can perform computation directly with native 2D signals (e.g., the image
and kernel), instead of having to “flatten” them for the sake of formulating the forward (and backward)
problems as systems of linear equations.
𝐅
𝐅 𝐅
ˆ𝐅 = 𝐌𝑁 𝐅 = ,
,𝑁 𝐌𝑁 (𝐌𝑁 𝐅 )
T T
(5)
Deliverable 3 [22.5%]
Implement a Generalized 2D DFT Routine: Complete the implementation of the generalized 2D
DFT function DFT2D, which extends the 1D generalized routine to 2D.
☑
Keep in mind that you don’t have to explicitly form these matrices at any point in order to take
advantage of the aforementioned factorization property.
Deliverable 4 [2.5%]
Implement the inverse 2D DFT Routine: Complete the implementation of the inverse 2D DFT
function iDFT2D that relies on calling your generalized DFT2D routine.
☑
4.1 2D Fourier Convolution
The (forward) convolution of a 2D function with a 2D kernel yields a 2D output function
. In the image convolution settings we have explored, is a
discretely-sampled input image, is a (typically square, odd edge-length) discrete kernel,
and is the blurred output image.
Given a matrix of the 2D sampled image pixel values, i.e., with , we obtain its DFT as per
Equation . For the kernel, we can also form a matrix of its sampled values before computing its DFT ;
here, however, we need take into account two additional points:
1. the kernel matrix needs to be 0-padded to match the size of the image matrix — this will be
needed to allow for an element-wise product of their DFTs, later on; and,
2. the “location”/placement of the non-zero kernel elements in the matrix needs to be chosen carefully
so to respect the same “wrap” boundary condition behaviour that we maintained in the primal
domain
1
.
1 The shifting property of convolutions can give you a hint of exactly how to structure and place the non-zero kernel
matrix elements in , before taking its DFT.
Given the appropriately constructed , and its DFT , we can obtain the DFT of the sampled output
image (i.e., with ) as
where is an element-wise product. Taking the real component of the inverse DFT of , we arrive at the
final output blurred image matrix .
𝐼(𝑥, 𝑦) 𝑘(𝑥, 𝑦)
𝐵(𝑥, 𝑦) = 𝐼(𝑥, 𝑦) ⊛ 𝑘(𝑥, 𝑦) 𝐼(𝑥, 𝑦) ≡ 𝐼[𝑥, 𝑦]
𝑘(𝑥, 𝑦) ≡ 𝑘[𝑥, 𝑦]
𝐵(𝑥, 𝑦) ≡ 𝐵[𝑥, 𝑦]
𝐈 (𝐈)𝑖,𝑗 = 𝐼(𝑖,𝑗) ˆ𝐈
5 𝐊 ˆ𝐊
𝐊 𝐈
𝐊
Deliverable 5 [25%]
Fourier Kernel Matrix: Complete the implementation of the FourierKernelMatrix that takes as
input an input blur kernel (with odd ), a 2-tuple of the shape (i.e., dimensions) of the
input/output image (image_shape), and (optionally) a forward 2D DFT function (DFT2D_func)
that you will call in your FourierKernelMatrix solution. The function returns the matrix
corresponding to the DFT of . The non-zero kernel elements in should be placed in the
matrix so to admit the same post-convolution wrap-based blurring behaviour as we saw in
Assignments 1 and 3.
☑
𝑛 × 𝑛 𝑛
ˆ𝐊
𝐊 𝐊
When implementing and debugging FourierKernelMatrix, you don’t necessarily have to rely on
your DFT2D routine (e.g., by passing something other than DFT2D as the DFT2D_func parameter).
ⓘ
𝐊
𝐊 ˆ𝐊 ˆ𝐁
𝐁 (𝐁)𝑖,𝑗 = 𝐵(𝑖,𝑗)
ˆ𝐁 = ˆ𝐈 ⋆ ˆ𝐊 ,
⋆ ˆ𝐁
𝐁
4.2 2D Fourier Deconvolution
In an unregularized setting, we can readily express the deconvolution problem using the DFT as
where is an element-wise division. Taking the real component of the inverse DFT of yields the
unblurred image matrix . Since the conditioning of the DFT and iDFT are 1, the underlying conditioning of
the convolution/deconvolution system remain as (un)stable as they were in the primal domain.
4.2.1 Laplacian Regularization
We provide you with a blurred image polluted by mild noise, as well as an blur kernel (which you will
need to extend appropriately into , as above).
As mentioned briefly in Assignment 3, Laplacian regularization is well-suited to image deblurring problems:
it adds a regularization term that penalizes large second derivatives in image-space, a statistical property of
many natural images.
In the primal domain, the Laplacian regularization kernel is
which you can extend and arrange into a matrix with dimensions that match the image resolution and with
structure that respects the wrapping boundary condition treatment (i.e., similar to the blur kernel’s matrix
structure).
We can solve the regularized least squares deconvolution problem, expressed in the primal domain,
in the frequency domain as
ˆ𝐈 = ˆ𝐁 ⋆ ,
−1 ˆ𝐊 (6)
⋆
−1 ˆ𝐈
𝐈
Unlike Assignments 1 and 3, we will not be providing you with the ground truth deblurred image for
your reference, below. You can code up example tests for forward convolution and nonnoisy/unregularized deblurring to sanity check your DFT/iDFT implementations; your code from A3 can
come in handy, here!
𝑛 × 𝑛
𝐊
3 × 3
𝑙(𝑥, 𝑦) ≡ 𝑙[𝑥, 𝑦] =
⎡
⎣
⎢
⎢
0
1
0
1
−4
1
0
1
0
⎤
⎦
⎥
⎥
(7)
𝐋
𝐊
min ‖𝐵(𝑥, 𝑦) − 𝐼(𝑥, 𝑦) ⊛ 𝑘(𝑥, 𝑦) + 𝜆‖𝐼(𝑥, 𝑦) ⊛ 𝑙(𝑥, 𝑦)
𝐼
‖
2
2 ‖
2
2
(8)
where is the DFT of and the absolute value of a complex value is . Taking the
real component of the inverse DFT of yields our regularized and deblurred image matrix . Note that in
Equation the regularized kernel is already in an inverted form, which is why an element-wise
multiplication ( ) is employed instead of an element-wise division ( , as in Equation ); i.e., with ,
Equations and are equivalent.
5 You’re Done!
Congratulations, you’ve completed the 4
th assignment. Review the submission procedures and guidelines
before submitting your Python script file assignment solution. Recall that any test code you incude in your
mainline (__main__) will not be graded, however it must still respect the collaboration/plagiarism polices.
formatted by Markdeep 1.13
ˆ𝐈 = ˆ𝐁 ⋆ ,
⎛
⎝
⎜
⎜
⎜
ˆ𝐊
+ 𝜆
∣
∣ˆ𝐊∣
∣
2
∣
∣ˆ𝐋∣
∣
2
⎞
⎠
⎟
⎟
⎟

ˆ𝐊
−1
reg
(9)
ˆ𝐋 𝐋 | ⋅ | 𝑎 + 𝑏ı 𝑎
2 + 𝑏 √‾‾‾‾‾‾‾2
ˆ𝐈 𝐈
9 ˆ𝐊
−1
reg
⋆ ⋆
−1 6 𝜆 = 0
6 9
Deliverable 6 [30%]
Laplacian-regularized Fourier Deconvolution Kernel: Complete the implementation of
LaplacianInverseFourierKernel to form and return . The function parameters match
those of FourierKernelMatrix. Choose and hardcode a regularization coefficient value that
yields a qualitatively good deblurred image (i.e., when applied to our blurred_noisy_image test
data using Equation ; see our test code in __main__: part of your grade for this deliverable will
be based on whether we can discern important details from the deblurred image.
☑
ˆ𝐊
−1
reg
𝜆
9
✒