Description
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
✒
