Hello, dear friend, you can consult us at any time if you have any questions, add WeChat: daixieit

CS 314 Numerical Methods

Fall 2022

Homework 6

You are permitted to discuss these problems with other students, but you should write down the names of your collaborators, including those who are not current CS 314 students, in your answer to the rst problem in this set. You should write up your own solution, and you are not permitted to either share a written copy of your solution or copy some one else’s written solution. Similarly, you are not permitted to share your code with another student. Please read the collaboration policy on HW and Codes included in Lecture 1, and ask if you have questions on whether something is permissible or not.

If you copy someone else’s work, or let your work be copied, you will get zero points for the entire HW, a loss of one letter grade in the course, and you will be reported to the Dean of Students. I follow the course policies described by Professor Gene Spafford (of Purdue Computer Science department) at the URL http://spaf .cerias .purdue .edu/ ˜spaf/cpolicy .html  If you are not familiar with this policy, please read it!    Please submit your assignments using Gradescope and Brightspace, following the instructions given by the Teaching Assistants.

0. If you collaborated with anyone in answering these questions, please write down your col- laborator’s name. If you did not collaborate with anyone, please indicate that too. If you do not answer this question, your HW will not be graded.

1. 20 points

We consider solving an n × n upper triangular system of equations Ux = y; here we are given the matrix U and the vector y, and we need to compute x. The diagonal elements of the matrix U are not equal to 1 in general (unlike the unit lower triangular matrix L in the Gauss elimination algorithm we learned). If any of the diagonal elements is equal to 0, or has absolute value less than max1oi,jon luij l * ∈ (here ∈ is the unit roundoff and uij  is the (i, j) element of U), then the algorithm should terminate with a failure message.

(a) Write down an algorithm to solve the system Ux = y. You can write your algorithm using pseudo-code as we have done in class, or Matlab code similar to the code for solving lower triangular systems of equations shown in the textbook on pages 140 and 141. Implement your algorithm in Matlab. Generate a random upper triangular matrix of order 15 using the commands triu and rand.  Compute the vector y to be the product of U and the vector xt, the vector of all ones with 15 components. Now solve for xc in the equation U xc = y using your code. Submit your code and the 1-norm of the vector xc - xt.

(b) Compute the number of operations that your algorithm or code takes to solve an n × n upper triangular system of equations. Your answer should report the the exponent e of the highest order term in n and the constant C multiplying it, i.e., Cne .

2.  15 points

Let x be an n-vector. Prove the following inequalities.

(a)  lxl|  < lxl1  < nlxl| .

(b)  lxl|  < lxl2  < ^n lxl| .

(c)  lxl2  < lxl1 .

3. 20 points

Let

A = .

Recall that we can compute a matrix norm by

lAlp  =  max  lA xlp ,

l&lp=1

for any vector p-norm, and also that we have obtained simpler expressions for these matrix norms when p = 1, o and 2. Do the following problems using paper and pencil; you can verify some of your answers using Matlab, but your must show your work.

(a) Compute the 1-norm of the matrix A, and find a vector x such that lA xl1  = lAl1 . (b) Compute the o-norm of the matrix, and nd a vector x such that lA xl|  = lAl| .

(c) Can you generalize?  Given an n × n matrix A write down a vector x such that the o-norm of Ax is the o-norm of A.

(d) Given a matrix A with real-valued elements, the matrix  lAl  is the matrix obtained by replacing each element of A by its absolute value. What is the relationship between the 1-norms of A and lAl?

4.  10 points

Recall that the p-norm condition number of a matrix is the product of the p-norm of a matrix and the p-norm of its inverse.

K(A)p  = lAlp  lA<1lp .

(a) Compute the o-norm condition number of the matrix

A = 0(1)   1(2).

You can compute the matrix A<1 by solving two systems of equations: Ax¹  = e1 , and Ax2  = e2 , where ei  is the i-th coordinate vector with 1 in the i-th component and 0

elsewhere. Then A<1  = (x1 x2 ). Again do this problem using paper and pencil.       (b) Let A be a square, non-singular matrix. Show that lA<1lp  > 1/ lAlp , where p > 1.

5. In this explore we will explore how sparse matrices and graph models are used to solve systems of linear equations Ax  =  b, with symmetric positive definite coefficient (SPD) matrices A.

Let us begin with a small, sparse matrix:

 

 -1    0      3    -1   -1    0   

( 0      0      0    -1   -1    3

By typing A=  sparse(A), you can get a sparse matrix in Matlab, i.e., a matrix whose nonzero elements are the only ones stored.  You can verify that the matrix A is SPD by computing its eigenvalues using the eigs command.  We will set up a linear system of equations to solve using A, by setting xt = àne&(n, 1), and b = A * xt.  Then we can compute the solution xc from xc = A/b, and compare it with the true solution xt.

However, for large, sparse matrices, we do not know if this method will terminate in a reasonable amount of time, say a few hours. So we take the five-fold path that helps us solve a sparse system of equations, to first estimate the cost and storage needed to compute the Cholesky factor.  If this cost is not prohibitive, we can proceed to compute the Cholesky factorization and the solution.

First we will solve this problem using the small matrix A given above, for practice. Then when you are able to solve the problem correctly, generate a larger, sparse matrix called wathen using the Matlab command A  =  gallery(’wathen’,  100,  100).  Note that this matrix has about 30, 000 rows (and columns) and about 470, 000 nonzeros. Report only the results for the wathen matrix.  (doc  gallery will give you more informa- tion about this finite element matrix.  You can generate smaller versions of this matrix to experiment by making the two numerical arguments smaller, say 10.)  Also generate the corresponding vectors xt and b for the linear system of equations.

(a) 5 points The first step on the path is to form a graph model. By typing H  =  graph(A) and then plot(H), you can see what the graph looks like (and how the vertices are numbered for small graphs). If we do this, the graph will look like a hairball due to the size of the graph and the resolution of the figure. Hence, only for this part, form the smaller matrix gallery(’wathen’,  10,  10) and plot its graph. Submit a plot of this smaller graph, and the number of vertices and edges in it. (The graph will look nicer if you remove the loops corresponding to the diagonal elements in the matrix. You can do this using the spdiags command. spdiags(A,  0) is the diagonal of the sparse matrix A. )

(b)  10 points The second step is to compute an ordering. By using nested dissection order- ing dissect and the Approximate Minimum Degree ordering amd in Matlab, you can compute permutation vectors, P 1 and P2 that find orderings that reduce the fill elements in the factorization.

The minimum degree ordering is another way to order vertices for elimination in the elimination game.  Here we choose a vertex v of minimum degree to eliminate first. (The degree of a vertex is the number of neighbors it has, not including the loop, the edge to itself.)  Since we add ll edges necessary to make all of v’s neighbors com- pletely connected (a clique in the graph), this is a greedy strategy to reduce the ll. We can reorder the matrix to get a permuted matrix AP = A(P, P), where P is the permutation vector obtained from the ordering algorithm. The construction above sym- metrically reorders the rows and columns of the matrix A using the permutation vector P. Submit spy plots of the permuted matrices from the two orderings.

(c) 5 points The third step is to compute an elimination tree from the ordering. The com- mand etree will compute the elimination tree of the permuted matrix, and treeplot will plot the tree. Report the heights of the elimination trees obtained from the nested dissection and AMD orderings.

(d) 5 points The fourth step computes the storage needed for the Cholesky factor L (or its adjacency graph, the filled graph) using symbolic factorization. symbfact computes the number of nonzeros in the rows of the upper triangular matrix R (equivalently the nonzeros in the columns of the lower triangular Cholesky factor L, where LT  = R). Report the total number of nonzeros in the Cholesky factor R1 obtained from the nested dissection ordering; do this also for the number of nonzeros in R2 obtained from the AMD ordering. Based on the results you have obtained so far, which ordering is better

for computing the Cholesky factorization, nested dissection or AMD?

(e)  10 points The fth step is to compute the numerical values of the nonzeros in the Cholesky factors R1 and R2 obtained from the two orderings, using the respective permuted matrices. Recall that Matlab’s version of the Cholesky factorization, chol, computes A = RT R, where R is upper triangular.

Hence we can solve Axc = b by rst computing R, and then solving the two triangular systems of equations RT y = b, and then R xc = y. Report the 1-norm of the error vector abs(xc  -  xt).

Remark: In the next HW, we will use our knowledge of sparse systems of equations to compute the Pagerank that ranks web pages obtained from a search of the WWW in order of importance.