a logo for the course

Matrix Computations

David Gleich

Purdue University

Fall 2017

Course number CS-51500

Tuesday and Thursday, 10:30-11:45am

Location Forney B124

Homework 5

Homework 5

Please answer the following questions in complete sentences in a typed, clearly prepared manuscript and submit the solution by the due date on Blackboard (Monday, October 2nd, 2017, around 4-5 am)

Problem 0: Homework checklist

Problem 1: The Cholesky Factorization

Through this problem, assume that the matrices are all symmetric. You may assume this in your code as well.

  1. Implement a Julia function for the Cholesky factorization.

  2. Research the Cholesky factorization and answer the following question, either providing a proof or a reference. Is the Cholesky factorization unique?

  3. Compare your Cholesky factorization code to Julia's chol function on the following matrix.

    n = 10
    Q,R = qr(randn(n,n))
    D = abs.(randn(n))
    A = triu(Q*diagm(D)*Q')
    A = A + triu(A,1)'
  4. Your professor boldy asserted that using the Cholesky factorization was faster than using the LU factorization. Make sure he's still right. Show the performance of Julia's LU and Cholesky factorization for matrices with size that varies from 10 to 500.

  5. Your professor also boldy asserted that we preserve positive definiteness after one step of the Cholesky factorization. Recall that this was: So we set and . Then, we recursively compute the Cholesky factorization of This only works if is positive definite. Prove that this is true. (Hint: you can do it from the definition of positive definiteness: for all .)

  6. Briefly explain why step 5 justifies that a Cholesky factorization exists for any positive definite matrix. (This is really like a one or two sentence answer.)

  7. One of the most useful aspects of the Choleskfy factorization is as a way to check if a matrix is positive definite or negative definite. A matrix is negative definite if for all . Change your Julia implementation to report if a matrix is not positive definite. (Hint: this relates to why a Cholesky factorization always exists for a positive definite matrix in the previous problem.) Use this to test if the matrix for Poisson's equation from Homework 4, Problem 7 is positive or negative definite.

Problem 2: Stability analysis

  1. Suppose that are a sequence of samples from a random variable. The sample variance of is: Compute the condition number of the sample variance.

  2. This is a real problem that I ran into! I was using a program written by someone else to compute a particular statistical estimate. The details don't really matter, but at one step, the code needed to compute the Hurwitz zeta function: where ranged from 1 to 7 and where ranged from 1 to 500. They were using Matlab, which doesn't provide a function for the Hurtwitz zeta function, but does provide a function for the Riemann zeta function: We will use Julia (which also has the zeta function) to continue this analysis. The final result doesn't differ much from Matlab's. Note that the only difference between these two functions is that the first term in the Riemann zeta function is and the first term in the Hurwitz zeta function is . To account for this difference, the code I was using did the following:

     function h=hzeta(s,q)
     z = zeta(s)
     h = z - sum((1:(q-1)).^(-s));

    In Julia,

     using SpecialFunctions
     function hzeta(s,q)
     z = zeta(s)
     h = z - sum(Float64.((1:(q-1))).^(-s))
     return h

    This requires loading the SpecialFunctions package

Was I happy when I realized this? You should use what you've learned about finite precision computations to answer this question. You can evaluate the Hurwitz zeta function to arbitrary precision via Wolfram Alpha. There is an answer that's more correct on this problem, but the most critical part of your answer is how you justify it.

Problem 3: Backwards stability

Only 1 part! Use the properties of floating point arithmetic to show that computing is backwards stable.