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 A_k (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
   • the final matrix A_k
   • the approximate eigenvalues of A returned as a vector (this is just the diagonal of A_k), sorted into increasing order; you can use Matlab's sort, which sorts by real part in the case of complex numbers
   • the norm of the lower triangular part of A_k (which is useful since we know that the smaller it is is, the more accurate the eigenvalues are)
   • the number of iterations that were needed
   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.

   • a randomly generated complex matrix, generated by A=randn(n) + i*randn(n) (where i is the imaginary unit; make sure you have not redefined it to be something else). Does the final A_k have small lower triangular entries? If not, either there is a bug in your program, or your iteration limit is too small.
   • a randomly generated complex Hermitian matrix, generated by A=randn(n) + i*randn(n); A = A + A'. What is the structure of the final A_k, neglecting small entries, and why?
   • a randomly generated real symmetric matrix, generated by A=randn(n); A = A + A'. What is the structure of the final A_k, neglecting small entries, and why?
   • a randomly generated real (but not symmetric) matrix, generated by A=randn(n). What is the structure of the final A_k, neglecting small entries, and why?
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 A_k? 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).