MA332 lecture notes, Fall 1998
Week 8 Thursday

Eigenvalues and eigenvectors continued

Faires & Burden Chapter 9

QR iteration
------------

Don't read section 9.5; it will make your head explode

Here's the basic idea:

Any matrix can be factored A = QR, where Q is orthogonal and R is
upper triangular.

Forget for now HOW we find Q and R, just assume we have.

Then if we perform the following iteration:

     factorize (Aold) -> Q, R
     Anew = RQ

It can be proved (not by me) that the off-diagonal elements get
smaller after every iteration.

If you start with a symmetric matrix you will eventually get
a diagonal matrix, with the eigenvalues on the diagonal.

If you start with a non-symmetric matrix you will eventually get an
upper-triangular matrix, again with the eigenvalues on the diagonal.

Furthermore, it can be proved (this time by me) that all the As
are similar (they have the same eigenvalues, and eigenvectors that
are a matrix multiplication away).

Similarity
----------

Two matrixes are similar if there is a matrix S such
that A = S-1 B S

If Ax = lx, then S-1 B S x = lx

multiply both sides by S,  and BSx = lSx, which implies
that Sx is an eigenvector of B and l is an eigenvalue.

Find an eigenvector of B, multiply it by S-1, and you have
an eigenvector of A.

Form a matrix with the eigenvectors of B.  Multiply it by
S-1 and you have the eigenvectors of A.

RQ iterations are similar
-------------------------

    Anew = RQ

    QR = Aold (multiply by Qt)  Qt Q R = Qt Aold

    Qt Q = I (because Q is orthogonal) R = Qt Aold

    Anew = Qt Aold Q

    Which implies that each iteration is a similarity transformation.


Performance issues
------------------

As a practical matter, QR converges very slowly for general
matrices.  Fortunately, there are two things we can do to
speed it up

1)  get rid of as many off-diagonal elements as possible using
    direct methods

2)  shift the eigenvalues so that the ratio of adjacent eigenvalues
    is maximized (similar to inverse shifted power method)

There are lots of ways to do (1).  For non-symmetric matrices,
you can use a special form of Gaussian elimination to get into
Hessenberg form (one subdiagonal).

For symmetric matrices, we can use Householder's method
to transform A into a similar matrix T that is tridiagonal.

Householder
-----------

A householder transformation is one with the form

  P = I - 2 w wt

  A' = P-1 A P = P A P

(since P is symmetric and orthogonal, P-1 = P)

Clearly A and A' are similar.

The goal is to find a Householder transformation that clobbers
the off diagonal elements in one row/column of A.

We start with k = 1 (to clobber the first row), and continue
until k = n-2.

Once again, the mathematical notation can be enormously simplified
by the introduction of some intermediate variables:

1) t is a copy of the kth row of A0

2) zero out the first k elements of t

3) m is the k+1th element of t

4) alpha = -sign(m) * ||t||     (that's a two norm)

5) r = sqrt(alpha/2 * (alpha - m))

6) subtract alpha from the k+1th element of t

7) w = t divided through by 2r

8) P = I - 2 * w * wt

9) A1 = P A0 P

Repeat for successive values of k.

We have to keep track of the various Ps so we can multiply them together.

   A2 = P1 A1 P1
   A3 = P2 A2 P2 = P2 P1 A1 P1 P2
   ...
   A = P T P
   P = P1 P2 ...Pn-2 

Compare this notation with the instructions on the top of page
388.

Here are some Maple hints.

1)  sign(m) does what you expect

2)  you can use the row command to select a row from a matrix,
    but make sure you load linalg before you use it

3)  use the norm command from the linalg package

4)  use evalm to perform matrix opertions, even things like
    dividing a vector by a scalar:  w:=evalm(t/2/r);

    (don't forget that division and mutiplication are evaluated from
     left to right, hence  t/ 2*r does not do what you expect)

5) you can create an identity matrix i:=array(1..4, 1..4, identity);
   but you can't modify it!

6) transpose(w) does what you expect

Example from book, page 388


Shifting
--------

The convergence of the off-diagonals to zero goes like

    (li / lj) ^ s         where s is the number of iterations
			  li and lj are eigenvalues

To eliminate a particular off-diagonal, we can estimate li with z
and then find the eigenvalues of

    A - z i

Now convergence is like  (li - z) / (lj - z), which is small
if z is closer to li than to lj.

Once the targeted element is zero (and it will go faster than the
others), you can reduce the matrix and repeat the process.

How do we get an initial estimate of li?

1) use prior estimate or application-specific information

2) use a diagonal element

3) form a 2x2 matrix using the ith and jth rows and columns,
   and find the eigenvalues of det(M - l I) (use the smaller)


Summary
-------

If you can factor A into QR, you can find eigenvalues and
eigenvectors of A by iterating  Anew = R Q

This goes much faster if you get rid of as many off-diagonals
as possible (for example, with Householder).

It also goes faster if you use an estimated eigenvalue to
shift the matrix.

We haven't talked yet about HOW to get the QR factorization.

For now, just use Maple:

    R := QRdecomp (A, Q='Q');

QRdecomp is in the linalg package.  Read the documentation!