Homework assignments for Math 221
Homework #1 due Tuesday, Sept 8:
From Chapter 1: 1,2,3,4,5,10,11
Answers to Homework #1
Homework #2 due Tuesday, Sept 15:
From Chapter 1: 7, 13, 14, 15
From Chapter 1: 16 (parts 1, 2, 4, 5, 6, 8, 10)
From Chapter 1: 20
Answers to Homework #2
Homework #3 due Tuesday, Sept 22:
From Chapter 2: 2, 7, 8, 11, 13 (parts 1 and 2), 17, 18 (part 1)
(Extra credit): Devise a sequential algorithm for solving TX=B for X, where T is n x n and triangular, and B is n x n, which minimizes communication in the following sense: The algorithm should do the same number of arithmetic operations as the naive algorithm (substitution to solve for column j of X from column j of B), but do Omega(n^3/sqrt(M)) loads and stores of data between slow memory and a cache of size M. Hint: Use a divide-and-conquer algorithm analogous to the one in lecture for matrix multiplication.
Answers to Homework #3
Homework #4 due Tuesday, Sept 29:
Question 1: 2.19. Hint for part 1: To show A is nonsingular, we need to show that if the column vector x is nonzero, so is x'*A. Suppose x(i) is the largest component of x in absolute value. What can you say about the i-th component of x'*A?
Question 2: Devise the fastest algorithm you can to compute the sum of all the entries in the even rows and odd columns of inv(A). In other words, give an algorithm for computing
sum_{even i and odd j} (inv(A))(i,j)
in c*n^3 + O(n^2) flops for the smallest constant c that you can.
Question 3, Part 1: Implement, in Matlab, the recursive LU algorithm RLU described in lecture, including partial pivoting. Test it by comparing the output to the output of Matlab's built-in lu() routine on a number of randomly generated matrices, as well as the matrices used in class to cause large pivot growth. Think carefully about how to compare the outputs: how close do they have to be to say the answers are "the same"? They might differ because of round-off (so that you'll have to compare the answers to within some suitably chosen error bounds). They might also differ because of "ties" in choosing the largest pivot. For example, test your solver vs Matlab on a matrix whose first column has all equal entries (and the other entries are all random).
Question 3, Part 2: Show that RLU does 2n^3/3 + O(n^2) flops (floating point operations). Hint: Write down a recurrence for the #flops, using explicit expressions for the number of flops in the calls in RLU to matmul and solving T*X=B (for T triangular and B rectangular). Then hypothesize an expression for #flops for calling RLU on a rectangular matrix, and plug it in and confirm that it satisfies your recurrence (up to O(n^2)).
Question 3, Part 3: Next, assume that both the calls in RLU to matmul and to solve T*X=B are implemented so as to minimize the number of words moved between slow memory and fast memory of size M, i.e. move only
O(max( #flops/sqrt(M), #inputs+#outputs )) words
where #flops is the number of floating point operations (which depends on the dimensions of the input matrices), and #inputs + #outputs is the total size of the input and output matrices (a lower bound for small matrices). Given this assumption, show that RLU moves only the minimum number of words, i.e. O(max( n^3/sqrt(M), n^2 )). Hint: Again, use a recurrence for #words_moved.
Question 3, Extra credit: The recursive LU routine is designed to minimize data movement, without having to choose a block size that depends on the (possibly unknown) cache size. But it needs to do both matmul and solve T*X=B for T triangular and B rectangular. These can be done using built-in Matlab functions, but for extra credit you should implement your own recursive versions of these as well.
Note: Because of Matlab's overheads, this will *not* be faster than the built-in Matab lu() routine (which calls the LAPACK code sgetrf), but illustrates how speed might be attained).
Answers to Homework #4
(updated Oct 9, 7:56pm, fixed notation in answer to Question 3, part 2)
Homework #5 due Tuesday, Oct 6
Question 1, Part 1: Consider an n by m matrix A in Compressed Sparse Row (CSR) data format: if A has nnz nonzero entries then
val(1:nnz) contains the nonzero entries packed row by row, from left to right in each row
colindex(1:nnz) contains the column index of each nonzero entry, i.e. val(i) is in column colindex(i) of A
rowbegin(1:n+1) points to where each row begins: the nonzeros of row j are contained in val( rowbegin(j) : rowbegin(j+1)-1 )
Write, in Matlab, explicit loops to compute the matrix-vector product y = A*x, where x and y are (dense) vectors of length m and n. (In other words, don't use any of Matlab's built in function for sparse matrices, though you may use notation like sum(a(i:j).*b(i:j)) or (a(i:j))'*b(i:j) for the dot product of parts of two column vectors.) Test it by running it on a few (small) sparse matrices, and comparing to Matlab's built-in routine for y=A*x. Make sure to test a matrix A with one or more zero rows. (Matlab has a built in sparse-matrix-vector-multiplication routine, that is called when doing A*x in Matlab on a matrix A that has been declared "sparse." But the point of this question is to understand CSR format, which is very commonly used.)
Question 1, Part 2: Now suppose that A is stored in Compressed Sparse Column (CSC) data format; this means storing A^T in CSR. Again write in Matlab explicit loops to compute y = A*x.
Question 1, Part 3: Now suppose that A is symmetric, and that only the upper triangle of A is stored in CSR to save space (triu(A,0) in Matlab notation). Again write in Matlab explicit loops to compute y = A*x.
Question 2: In class we defined a weighted undirected graph G=(V,E,W) where V is the set of vertices, E the set of edges and W is the set of edge weights. Now we define a directed graph, where each edge (u,v) is considered to be directed from u to v (think of an arrow pointing from vertex u to vextex v). Since (u,v) and (v,u) are now different, we see that there is a natural correspondence between a general, possibly unsymmetric square matrix A and a graph G=(V,E,W) with
one vertex for each row/column
an edge (i,j) if A(i,j) is nonzero, with weight A(i,j)
Recall also that a path of length k from vertex v(1) to vertex v(k+1) in a graph is a set of k edges that connect end to end
(v(1),v(2)), (v(2),v(3)) ,..., (v(k),v(k+1))
where, for a directed graph, the directions of the edges have to match as shown.
Now consider two nxn matrices and their associated graphs on the same set of vertices:
A_r and G_r = (V, E_r, W_r) (the "red" edges)
A_b and G_b = (V, E_b, W_b) (the "blue" edges)
with all edges weights equal to 1; i.e. the nonzeros entries of A_r and A_b are all ones.
Question 2, Part 1: Show that the number of paths consisting of one red edge followed by one blue edge connecting vertex i to vertex j is the (i,j) entry of A_r * A_b. In other words, multiplying matrices counts paths.
Question 2, Part 2: Consider just G_r. Show that the number of paths of length k from vertex i to vertex j in G_r is the (i,j) entry of G_r^k.
Answers to Homework #5
SpMV_test.m
SpMV_CSR.m
SpMV_CSC.m
Homework #6 due Tuesday, Oct 13
Chapter 3, Question 3 (just parts 1, 2 and 3)
Chapter 3, Question 4 (just the normal equations part)
Chapter 3, Questions 9, 12, 13, 15
Answers to Homework #6
Homework #7 due Tuesday, Oct 20
Chapter 3, Question 7
Chapter 3, Question 17 (just part 1)
Chapter 4, Questions 1,2,3,4,5
Answers to Homework #7
QR_Ruv.m
Givens.m
Homework #8 due Thursday, Oct 29 (note delayed due date)
Chapter 4, Question 6, 7, 8, 11, 13
Chapter 4, Question 16. You do not have to write a working Matlab program (though you can for extra credit). Instead, just describe how to find the intersection(s) of a surface S and a curve C by setting it up as an eigenvalue problem. Write down the eigenproblems explicitly for the surface/curve pairs #1 and #6 on page 193. Suggest what might go wrong if either S and C do not intersect, or if C lies entirely in S.
Answers to Homework #8
Extra Credit Answers to Homework #8
Homework #9 due Tuesday, Nov 10
Chapter 5, Question 1, 2
Chapter 5, Question 4 (just part 1)
Chapter 5, Question 16, 18, 28
Extra Credit: Chapter 5, Question 10