codingprolab@gmail.com

$30.00

Order Now
Category: ECSE 343

Description

5/5 - (5 votes)

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

✒

WhatsApp us