The PDF file attached has the assignment details. The .zip file contains an IPYNB file where the assignment needs to be done.
1/29/2020 ECSE 443: Introduction to Numerical Methods in Electrical Engineering www.cim.mcgill.ca/~derek/ecsex43/A1.md.html 1/30 ECSE 443: Introduction to Numerical Methods in Electrical Engineering Assignment 1 — Linear least squares Due: Monday, Febuary 10th, 2020 at 11:59pm EST on myCourses Final weight: 15% Contents 1 Assignment submission process 1.1 Late policy 1.2 Python import statements 2 Forward and backward substitution [20%] 2.1 Forward substitution [10%] 2.2 Backward substitution [10%] 3 Polynomial fitting via linear least squares [50%] 3.1 Learning objectives 3.2 Solving Normal Equations using Cholesky 3.2.1 Cholesky decomposition [10%] 3.2.2 Applying Cholesky to the Normal Equations [5%] 3.3 QR Decomposition 3.3.1 Householder transformations [10%] 3.3.2 QR Decomposition using Householder Reflections [15%] 3.3.3 Solving Linear Least Squares using QR [5%] 3.4 Polynomial fitting [5%] 4 Image Reconstruction using Linear Least Squares [30%] 4.1 Learning objectives https://mycourses2.mcgill.ca/ 1/29/2020 ECSE 443: Introduction to Numerical Methods in Electrical Engineering www.cim.mcgill.ca/~derek/ecsex43/A1.md.html 2/30 4.2 Loading the data 4.3 Regularized Image Reconstruction [30%] 1 Assignment submission process Download and modify the ipynb notebook we provide on myCourses and submit your modified file, renamed according to your student ID, as [YourStudentID].ipynb For example, if your id is 234567890, your submission filename should be 234567890.ipynb. 1.1 Late policy This assignment is to be completed individually. 1.2 Python import statements We provide all import statements required to complete the assignment. 2 Forward and backward substitution [20%] Note that every time you submit a new file, your previous submission will be overwritten. All submissions must be made using myCourses. You can submit as many times as you like, but we will only grade the last submission. ⚠ Late Policy: non-negotiable time-based penalty You will lose 5% per hour past the deadline, up to a maximum of −100%. ☒ import Policy: you must not use any imports other than the ones we provide. Doing so will result in a score of zero (0%) on the assignment. ☒ 1/29/2020 ECSE 443: Introduction to Numerical Methods in Electrical Engineering www.cim.mcgill.ca/~derek/ecsex43/A1.md.html 3/30 We will start by implementing the forward and backward substitution algorithms to solve lower and upper triangular linear systems of equations. These algorithms will also form the foundation for more complex algorithms, later on. 2.1 Forward substitution [10%] Forward substitution is the process of solving a linear system of equations with a lower triangular coefficient matrix . The formula below defines the elements of that are computed using forward substitution: from numpy import array, zeros, tril, triu, absolute Run the following cell in the Jupyter Notebook to include the required modules. Lx = b L xi x =xi ⎧ ⎩ ⎨ ⎪⎪⎪⎪ ⎪⎪⎪⎪ b1 L1,1 ( − )1 Li,i bi ∑ j=1 i−1 Li,jxj if i = 1, otherwise. ☑ Complete the forward substitution algorithm implementation in the cell below: 1/29/2020 ECSE 443: Introduction to Numerical Methods in Electrical Engineering www.cim.mcgill.ca/~derek/ecsex43/A1.md.html 4/30 def ForwardSubstitution(L, b): """ Perform forward substitution of to solve Lx = b, where L is a lower triangular matrix Parameters: L (np.array((n, n))): Left hand side (LHS) lower triangular matrix b (np.array(n)): Right hand side (RHS) vector Returns: x np.array(n): Solution vector """ # Save dimension of U m, n = L.shape # L must be square if m != n: raise Exception("Sorry, it's hip to be square!") # Initialize x x = zeros(n) # INSERT YOUR CODE BELOW (remove the NotImplementedError) raise NotImplementedError() # Outer loop: loop over all rows, from left to right # Set x[i] as b[i] # Inner loop: loop over each column from left to diagonal # Update x[i] # divide by diagonal value to isolate x[i] return x Run the following cell to test your forward substitution solver. 1/29/2020 ECSE 443: Introduction to Numerical Methods in Electrical Engineering www.cim.mcgill.ca/~derek/ecsex43/A1.md.html 5/30 2.2 Backward substitution [10%] Backward substitution is the process of solving a linear system of equations with an upper triangular coefficient matrix . These equation similarly define the solution element of that you will compute using backward substitution: L = tril([[4,12,-16],[12,37,-43],[-16,-43,98]]) b = array([1,2,3]) x = ForwardSubstitution(L, b) print('L = ', L) print('b = ', b) print('x = ', x) print('L.dot(x) = ', L.dot(x)) try: assert((absolute(L.dot(x) - b) < 1e-5).all())="" print('you="" computed="" the="" correct="" result!')="" except:="" raise="" exception("if="" at="" first="" you="" don't="" succeed...")="" ux="b" u="" xi="" x="xi" ⎧="" ⎩="" ⎨="" ⎪⎪⎪⎪="" ⎪⎪⎪⎪="" bn="" un,n="" (="" −="" )1="" ui,i="" bi="" ∑="" j="i+1" n="" ui,jxj="" if i="n," otherwise. ="" ☑="" complete="" the="" backward="" substitution="" algorithm="" implementation="" in="" the="" cell="" below:="" 1/29/2020="" ecse="" 443:="" introduction="" to="" numerical="" methods="" in="" electrical="" engineering="" www.cim.mcgill.ca/~derek/ecsex43/a1.md.html="" 6/30="" def="" backwardsubstitution(u,="" b):="" """="" perform="" back="" substitution="" of="" to="" solve="" ux="b" where="" u="" is="" an="" upper="" triangular="" matrix="" parameters:="" u="" (np.array((n,="" n))):="" lhs="" upper="" triangular="" matrix="" b="" (np.array(n)):="" rhs="" vector="" returns:="" x="" np.array(n):="" solution="" vector="" """="" #="" save="" dimension="" of="" u="" m,="" n="U.shape" #="" u="" must="" be="" square="" if="" m="" !="n:" raise="" exception("sorry,="" it's="" hip="" to="" be="" square")="" #="" initialize="" x="" x="zeros(n)" #="" insert="" your="" code="" below="" (remove="" the="" notimplementederror)="" raise="" notimplementederror()="" #="" outer="" loop:="" loop="" over="" all="" rows="" from="" bottom="" to="" top="" #="" set="" temporary="" variable="" as="" the="" current="" b[i]="" #="" inner="" loop:="" loop="" over="" each="" column="" from="" right="" to="" diagonal.="" #="" update="" temporary="" variable="" from="" the="" already="" known="" x[j]="" and="" u[i,j]="" #divide="" by="" diagonal="" value="" to="" isolate="" x[i]="" return="" x="" run="" the="" following="" cell="" to="" test="" your="" backward="" substitution="" solver.="" 1/29/2020="" ecse="" 443:="" introduction="" to="" numerical="" methods="" in="" electrical="" engineering="" www.cim.mcgill.ca/~derek/ecsex43/a1.md.html="" 7/30="" 3="" polynomial="" fitting="" via="" linear="" least="" squares="" [50%]="" we="" will="" implement="" our="" own="" linear="" least="" squares="" solver="" using="" two="" different="" approaches,="" testing="" them="" on="" a="" polynomial="" fitting="" problem.="" given="" an="" overdetermined="" matrix="" of="" dimension="" ,="" where="" ,="" and="" a="" column="" vector="" ,="" we="" want="" to="" find="" the="" column="" vector="" that="" best="" fits="" the="" linear="" equation="" in="" the="" least-squares="" sense.="" we="" write="" instead="" of="" because="" the="" true="" linear="" system="" is="" not="" generally="" satisfiable.="" specifically,="" our="" fit="" will="" minimize="" the="" squared="" 2-norm="" of="" the="" residual="" vector="" ,="" i.e.:="" an="" example="" least-squares="" fit="" of="" a="" higher-order="" polynomial.="" u="triu([[4.,12.,-16.],[12.,37.,-43.],[-16.,-43.,98.]])" b="array([1.,2.,3.])" x="BackwardSubstitution(U," b)="" print('u=', U) print(' b=', b) print(' x=', x) print(' u.dot(x)=', U.dot(x)) try: assert((absolute(U.dot(x) - b) < 1e-5).all()) print(' you\'re="" on="" a="" roll!')="" except:="" raise="" exception("better="" luck,="" next="" time...")="" a="" m="" ×="" n="" m=""> n m × 1 b n × 1 x Ax ≈ b Ax ≈ b Ax = b x r = Ax − b .argmin x ∥Ax − b∥2 http://www.cim.mcgill.ca/~derek/ecsex43/polyfit.png 1/29/2020 ECSE 443: Introduction to Numerical Methods in Electrical Engineering www.cim.mcgill.ca/~derek/ecsex43/A1.md.html 8/30 3.1 Learning objectives The goals of this assignment component are for you to: 1. implement a Cholesky decomposition and apply it to the linear least squares solution expressed using the normal equations; 2. implement a QR decomposition and apply it to the linear least squares solution expressed using the normal equations; 3. understand when using each of these methods is preferrable to the other; and, 4. to solve a overdetermined polynomial fitting problem using linear least squares. 3.2 Solving Normal Equations using Cholesky 3.2.1 Cholesky decomposition [10%] Suppose that is a square matrix of dimension satisfying the following properties: 1. is symetric, e.g., , and 2. is positive definite, e.g., for all column vectors of dimension such that . Then, admits the unique Cholesky decomposition where is a lower triangular matrix with real and positive diagonal entries. After some algebraic manipulation, we can arrive at: from numpy import array, linspace, identity, copy, dot, zeros, zeros_like, sqrt, outer, polyval, poly1d, absolute from numpy.random import random from numpy.linalg import cond, norm import matplotlib.pyplot as plt Run the following cell to include the required modules. M n × n M M = M⊺ M xM > 0x⊺ x 1 × n ∥x∥ ≠ 0 M M = L ,L⊺ L 1/29/2020 ECSE 443: Introduction to Numerical Methods in Electrical Engineering www.cim.mcgill.ca/~derek/ecsex43/A1.md.html 9/30 We see that we can compute the entry of if we know the entries to the left and above. The Cholesky- Banachiewicz algorithm starts from the upper left corner of the matrix and proceeds to calculate its elements row by row. Access pattern (white) and writing pattern (yellow) for the Cholesky-Banachiewicz algorithm on a matrix. =Li,j ⎧ ⎩ ⎨ ⎪⎪⎪⎪⎪⎪⎪ ⎪⎪⎪⎪⎪⎪⎪ −Mj,j ∑ k=1