CSCI 2033: Elementary Computational Linear Algebra (Fall 2020) Assignment 2 (100 points) Due date: November 13th, 2020 11:59pm In this assignment, you will implement a Matlab function to decompose a matrix into lower and upper triangular matrices (L and U), i.e., A = LU, as outlined in Section 2. In addition, you must complete either Section 3,4 (Matlab implementation of LU with partial pivoting, and inverse computation) or Section 5 (written questions). If you choose to complete both, your highest grade out of the two will be used to compute your final grade for this assignment. For each function that you are asked to implement, you will need to complete the corresponding .m file with the same name that is already provided to you in the zip file. Failing to do so may result in points lost. In the end, you will zip up all your complete .m files and upload the zip file to “HW2 Programming” on Gradescope. For the written portion, please scan or compile your solution into a PDF and submit via “HW2 Written” on Gradescope. Note that you do not need to submit anything to “HW2 Written” if you wish to complete Section 3 and 4 of the assignment instead, but you have to submit your completed code for Section 2 since that part is mandatory. Note on plagiarism A submission with any indication of plagiarism will be directly reported to University. Copying others’ solutions or letting another person copy your solutions will be penalized equally. Protect your code! 1 Submission Guidelines Your zip file should contain the following .m files on Gradescope (you only need to submit the first two files if you are doing the written portion): • replacement_lu.m • my_lu.m • interchange_lup.m • my_lup.m • find_inverse.m You may use Matlab’s high-level array manipulation syntax, e.g. size(A), A(i,:), abs, and [A,B], but the built-in linear algebra functions such as inv, lu, rank, lin_solve, det, rref, and A\b are not allowed. Please contact the TAs for further questions about specific function(s). 1 2 LU Decomposition (50 points) In this part, you will write functions that perform LU decomposition without pivoting. We will deal with pivoting in the next part of the assignment. First, you will implement a simple modular function that will help you with LU decomposition. Notation: for subsequent sections, Ai indicates the i th row of A, and Aij indicates the (i, j) entry of A. Specification: function [U_new, L_new] = replacement_lu(U, i, j, s, L) Input: two square matrices U and L, two integers i and j, and a scalar s. Output: • Unew: updated U by subtracting s times row j from row i, i.e. performing the row replacement operation Ui ← Ui − sUj . • Lnew: updated L by filling in its (i, j) entry with s, i.e., Lij ← s. Warning! If you are adapting your code from Assignment 1, be careful that we are now subtracting s times row j from row i instead of adding it. Now, you will implement of the main portion of LU decomposition. function [L, U] = my_lu(A) Input: an n× n square matrix A. Output: • L: an n× n lower triangular matrix where the diagonal entries are all one, e.g., 1 0 0∗ 1 0 ∗ ∗ 1 where ∗ is a potentially nonzero element. • U: an n× n upper triangular matrix. To get full credit, your function should handle the following case: • Early termination: Due to round-off error, there may exist free variables whose pivots are extremely small but not precisely zero. You should terminate your LU decomposition if the absolute value of a pivot is less than 10−12. The process of LU decomposition uses Gaussian elimination that transforms A to an upper triangular matrix U while recording the pivot multipliers in a lower triangular matrix L. 1. Initialize L to the identity matrix, and U to A. You can use Matlab’s built-in function eye(n). 2. At the ith step, for each row k below the ith row, 2 (a) Record the pivot multiplier to L at (k, i), i.e., Lk,i = Uk,i/Ui,i as shown in Figure 1. Note that this fills in the ith column of L. (b) Reduce the kth row of U using the pivot multiplier, i.e., Uk = Uk − (Uk,i/Ui,i)Ui where Ui is the ith row. U = i k i L = i k i 1 0 0 0 1 1 0/ki iiU U ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗∗∗ 0 0 0 11U iiU ∗ ∗0 0 0 00 0 Figure 1: LU decomposition at the ith step. We provide pseudocode for LU decomposition in Algorithm 1. Algorithm 1 LU decomposition of A 1: n← the number of columns of A 2: L← In where In is n× n identity matrix. 3: U ← A 4: for 1 ≤ i ≤ n− 1 do 5: if Uii is nearly zero* then 6: return L, U . Early termination 7: end if 8: for i + 1 ≤ k ≤ n do 9: p← Uki/Uii 10: [U , L] = replacement_lu(U , k, i, p, L) 11: end for 12: end for 13: return L, U *Consider a number to be nearly zero if its absolute value is smaller than 10−12. Note: When terminating early, L is not unique, i.e., the order of rows correspond- ing to zero rows in U can be different. As long as the resulting decomposition satisfies A = LU, the solution L and U are valid. Test cases: • [L,U] = my_lu([4,-2,2;-2,5,3;2,3,9]) L = 1 0 0−0.5 1 0 0.5 1 1 , U = 4 −2 20 4 4 0 0 4 3 • Early termination: [L,U] = my_lu([4,-2,2,1;-2,5,3,3;4,-2,2,1;-2,5,3,3]) L = 1 0 0 0 −0.5 1 0 0 1 0 1 0 −0.5 1 0 1 , U = 4 −2 2 1 0 4 4 3.5 0 0 0 0 0 0 0 0 . Without early termination, it will return L and U with NaN (Not-a-Number) entries. • You can test your algorithm by constructing a random n×n matrix A, e.g., A=rand(10,10), and testing whether multiplying L and U reconstructs A. You also need to check whether L is a lower triangular matrix with diagonal entries equal to 1, and U is an upper triangular matrix. 4 3 LU Decomposition with Partial Pivoting (20 points) Based on your my_lu, you will write numerically stable LU decomposition with partial pivoting. At the ith step of LU decomposition (ith pivot column), you will find the row that has the largest absolute value in the pivot column (say row j), and swap the ith and jth rows of U as usual. Simultaneously, you will swap the partial entries of the ith and jth rows of L, and record the row exchange in a permutation matrix P. For further details, please see http://www.math.kent.edu/~reichel/courses/intr.num.comp.1/fall09/lecture9/lecture4.pdf First of all, similar to above, you will implement a simple modular function that will help you with LU decomposition with partial pivoting. Specification: function [U_new, L_new, P_new] = interchange_lup(U, i, j, L, P) Input: three same size square matrices U, L, and P, and two integers i and j. Output: • Unew: updated U by swapping rows i and j, i.e. Ui ↔ Uj . • Lnew: updated L by swapping only the first (i− 1) entries of rows i and j, i.e., Li,1 ↔ Lj,1, . . . , Li,(i−1) ↔ Lj,(i−1) as shown in Figure 2. • Pnew: updated P by swapping rows i and j, i.e. Pi ↔ Pj . i j i 1 0 0 0 1 0 0 1 0 0 ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ i j i 1 0 0 0 1 0 0 1 0 0 ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ Partial row exchange n n i-1 i-1 Figure 2: Partial row exchange in L. Now, you will implement of the main portion of LU decomposition with partial pivoting. function [L, U, P] = my_lup(A) Input: an n× n square matrix A. Output: • L: an n× n lower triangular matrix where the diagonal entries are all 1. • U: an n× n upper triangular matrix. 5 • P: an n× n permutation matrix. The process of LU decomposition with partial pivoting needs to compute an additional row permutation matrix P. 1. Initialize L and P to the identity matrix, and U to A. You can use Matlab’s built-in function eye(n). 2. At the ith step, (a) Similar to Assignment 1, perform partial pivoting in U. (b) Record the row exchange to the permutation matrix P by exchanging its corresponding rows. (c) Exchange the corresponding row entries in L to reflect the row ex- change in U. Note that this exchange only involves with the row entries that have been recorded, i.e., not entire row exchange. (d) For each row k below the ith row, record the pivot multiplier to L and replace the row in U using the pivot multiplier, like in the previous part of this assignment. We provide pseudocode for LUP decomposition in Algorithm 2. Algorithm 2 LUP decomposition of A 1: n← the number of columns of A 2: L← In where In is n× n identity matrix. 3: P ← In . Permutation initialization 4: U ← A 5: for 1 ≤ i ≤ n− 1 do 6: j ← the row index of the largest absolute value among rows i to n in the ith column of U 7: [U , L, P ] = interchange_lup(U , i, j, L, P ) 8: if Uii is nearly zero* then 9: return L, U , and P . Early termination 10: end if 11: for i + 1 ≤ k ≤ n do 12: p← Uki/Uii 13: [U , L] = replacement_lu(U , k, i, p, L) 14: end for 15: end for 16: return L, U , P The red lines specify the major modification for partial pivoting. *As before, consider a number to be nearly zero if its absolute value is smaller than 10−12. 6 Note: When terminating early, L and P are not unique, i.e., the order of rows corresponding to zero rows in U can be different. As long as the resulting decomposition satisfies PA = LU, the solution L and P are valid. Test cases: • Partial pivoting: [L,U,P] = my_lup([-2,5,3;2,3,9;4,-2,2]) L = 1 0 00.5 1 0 −0.5 1 1 , U = 4 −2 20 4 8 0 0 −4 , and P = 0 0 10 1 0 1 0 0 • Early termination: [L,U,P] = my_lup([4,-2,2;-2,5,3;4+5E-13,-2+2E-13,2+7E-13]) L = 1 0 0−0.5 1 0 1 0 1 , U = 4 −2 20 4 4 0 0 0 , and P = 0 0 10 1 0 1 0 0 • You can test your algorithm by comparing its results with Matlab’s built-in function, [l, u, p] = lu(A), on a randomly generated n × n matrix A, e.g., A=rand(10,10). 7 4 Finding Inverse via Identity Matrix (30 points) In this section, you will implement the algorithm taught in class, which finds the inverse of a matrix via reduction of the identity matrix (i.e., [A|I]⇒ [I|A−1]). Using any other methods will result in zero credit. You are not allowed to use functions inv and det in any circumstances. In your function, you should return the inverse of a matrix A if there is one, and return false if there is not. Specification: function output = find_inverse(A) Input: a matrix A of size m x n, where m is not necessarily equal to n. Output: return A−1 when A has an inverse, and return false otherwise. The key part of this implementation is to reduce A using RREF and apply the corresponding operations on I. Below, we provide pseudocode for how to do so by slightly modifying the pseudocode from HW1. Note that this pseudocode does not include the part of checking whether A has an inverse or not. Also, make sure when you compare values, allow a 10−12 tolerance because of the Matlab approximation issues. Algorithm 3 Reducing I via RREF of A 1: initialize I to an identity matrix 2: initialize pivot row k = 1, pivot column l = 1 3: while 1 ≤ k ≤ m and 1 ≤ l ≤ n do 4: among rows k to m in A, find row p with the biggest* entry in column l 5: Rk ↔ Rp for R in A 6: Rk ↔ Rp for R in I 7: if Akl is zero** then 8: l← l + 1 9: else 10: Rk ← (1/Akl)Rk for R in I 11: Rk ← (1/Akl)Rk for R in A 12: for i = 1, . . . , k − 1, k + 1, . . . ,m do 13: Ri ← Ri − (Ail/Akl)Rk for R in I 14: Ri ← Ri − (Ail/Akl)Rk for R in A 15: end for 16: k ← k + 1, l← l + 1 17: end if 18: end while *Here “biggest” means having the largest absolute value (abs in Matlab). **Consider a number to be zero if its absolute value is smaller than 10−12. 8 Test cases: find_inverse([1 2; 1 3]) should return [3 -2; -1 1], i.e. A = [ 1 2 1 3 ] → output = [ 3 −2 −1 1 ] find_inverse([1 4; 2 8]) should return false, i.e. A = [ 1 4 2 8 ] → output = false find_inverse([1 ; 2 ; 3]) should return false, i.e. A = 12 3 → output = false find_inverse([1 0 0 2 0 0 3 0 0]) should return false, i.e. A = [ 1 0 0 2 0 0 3 0 0 ]→ output = false find_inverse([3,0,0,1; 0,3,0,1; 0,0,3,1; 1,1,1,1]) should return false, i.e. A = 3 0 0 1 0 3 0 1 0 0 3 1 1 1 1 1 → output = false find_inverse([3,0,0,1; 0,3,0,1; 0,0,3,1; 1,1,1,0]) should return the following: A = 3 0 0 1 0 3 0 1 0 0 3 1 1 1 1 0 → output = 2/9 −1/9 −1/9 1/3 −1/9 2/9 −1/9 1/3 −1/9 −1/9 2/9 1/3 1/3 1/3 1/3 −1 9 5 Written Problems (50 points) 1. For each of the following statements, determine whether it is guaranteed to be true. If yes, provide an explanation. If not, provide a counterexample. (a) Given two invertible matrices A and B, their sum must also be invertible. (b) Given two square matrices A and B, if AB is invertible, then A and B must be both invertible. 2. Consider the matrix A = a a aa a 2 a 1 0 . Answer the following questions: i) For what values of a is A not invertible? Show your steps. ii) For each of these values, why is A not invertible? Show all steps for full credits. 10 3. Let A = 1 2 31 2 5 3 7 10 . Find its LU decomposition with partial pivoting. Show all steps for full credits. 4. Let A = 3 3 33 7 4 3 11 7 . Find the elementary matrices E1,E2,E3 such that E3E2E1A = U. Then, using only the elementary matrices you found, calculate L. Do not compute L based on your reduction from A to U. Note that L and U here refer to the LU decomposition of A. Show all steps. 11
欢迎咨询51作业君
