Homework 4, assigned Oct 19, due Nov 3

  1. Write a function that implements the basic QR iteration for finding eigenvalues (Algorithm 4.7 in the text). This basic algorithm is very simple, as all the work is done by a call to the built-in function qr. All you need to do is add a termination criterion, stopping when the norm of the lower triangular part of Ak (square root of the sum of the squares of the strictly lower triangular entries) is sufficiently small, or a maximum iteration limit is exceeded. The input parameters should be the matrix A, the tolerance and the iteration limit. The output parameters should be According to the discussion in the text, this algorithm should work for matrices with the property that all the eigenvalues have different magnitudes (complex modulus), although depending on the input, it might need a lot of iterations.

    Test the program on the following examples, comparing your results with the built-in Matlab function eig(A). For printing small examples, use n=7 for real matrices and n=3 for complex matrices, so you can easily display the matrix with format short e. Try some runs with bigger n as well (e.g. n=10, n=30), but don't print these matrices.

  2. Modify your program so that it has a different stopping condition when A is real but not symmetric, which can easily be detected before running the iteration. The stopping condition should be that the matrix is close to being in real Schur form, that is block triangular with 1 by 1 and 2 by 2 diagonal blocks (see p.172). This is a little tricky to check and will require some thought. Once you have this working, modify the program to return the eigenvalues as a vector of real and complex numbers, where the complex ones are found by computing the eigenvalues of the 2 by 2 diagonal blocks using the standard quadratic formula. As before, compare with the results obtained by eig.
  3. Write another function to compute an eigenvector x of a given matrix A corresponding to a specified approximate eigenvalue lambda. This is very simple: just take one step of inverse iteration, computing (A-lambda*I)\v, where I = eye(size(A) and v is randomly generated, and then normalize the result to have norm 1. The reason this works is that if lambda is very close to being an eigenvalue of A, then A-lambda*I is very close to being singular, so the result of the operation is close to being a null vector of the nearby exactly singular matrix. (You can compare it with the output of [X,D]=eig(A) if you want, but this can be confusing, since an eigenvector can be normalized in different ways.) Test this on the same examples you used previously, checking the residual norm(A*x-lambda*x); if this is not small, there is a bug in your program or lambda is not close to being an eigenvalue. Which matrix should you pass to the eigenvector function: the original A, or the final Ak? Why?
  4. Modify your eigenvalue program (from Questions 1 and 2) further to include shifts (Algorithm 4.8). This makes it quite a bit more complicated; you have to terminate the iteration one eigenvalue at a time. An input parameter should be added to specify whether or not shifts are required. Since this is even more complicated in the real nonsymmetric case, you can assume the input will not be real nonsymmetric (and if it is, print an informative message and ignore the request to use shifts). You first shift by the bottom right entry (the Rayleigh quotient shift). When the bottom row, except the bottom right entry, is small enough, you declare the bottom right entry to have converged to an eigenvalue and now restrict your attention to the leading (n-1) by (n-1) matrix. Repeat this until all n eigenvalues have converged. If you do this correctly, the number of iterations required should be greatly decreased, and you should see quadratic convergence (doubling of the number of correct digits in an eigenvalue each step) in the general complex case, and cubic convergence (tripling of the number of correct digits each step) in the complex Hermitian or real symmetric case. Include printing of examples with n=3 (complex) and n=7 (real symmetric) examples, using format short e to demonstrate this, showing the convergence. This may require a few pages of output (but only a few).