Core-chasing algorithms for the eigenvalue problem

This monograph is about a class of methods for solving matrix eigenvalue problems. Of course the methods are also useful for computing related objects such as eigenvectors and invariant subspaces. We will introduce new algorithms along the way, but we are also advocating a new way of viewing and implementing existing algorithms, notably Francis's implicitly-shifted QR algorithm [36]. Our first message is that if we want to compute the eigenvalues of a matrix A, it is often advantageous to store A in QR-decomposed form. That is, we write A = QR, where Q is unitary and R is upper triangular, and we store Q and R instead of A. This may appear to be an inefficient approach but, as we shall see, it often is not. Most matrices that arise in applications have some special structures, and these often imply special structures for the factors Q and R. For example, if A is upper Hessenberg, then Q is also upper Hessenberg, and it follows from this that Q can be stored very compactly. As another example, suppose A is unitary. Then Q = A, and R is the identity matrix, so we don't have to store R at all. Every matrix can be transformed to upper Hessenberg form by a unitary similarity transformation. We will study this and related transformations in detail in Chapter 6, but for the early chapters of the book we will simply take the transformation for granted; we will assume that A is already in Hessenberg form. Thus we consider an arbitrary upper Hessenberg matrix in QR-decomposed form and show how to compute its eigenvalues. Our method proceeds by a sequence of similarity transformations that drive the matrix toward upper triangular form. Once the matrix is triangular, the eigenvalues can be read from the main diagonal. In fact our method is just a new implementation of Francis's implicitly-shifted QR algorithm. The storage space requirement is O(n^2) because we must store the upper-triangular matrix R, and the flop count is O(n^3). Once we know how to handle general matrices, we consider how the procedure can be simplified in special cases. The easiest is the unitary case, where R = I. This results in an immediate reduction in the storage requirement to O(n) and a corresponding reduction of the computational cost to O(n^2) flops. A similar case is that of a companion matrix, which is both upper Hessenberg and unitary-plus-rank-one. This results in an R that is unitary-plus-rank-one. Once we have figured out how to store R using only O(n) storage, we again get an algorithm that runs in O(n 2 ) flops. The unitary case is old [42], but our companion algorithm is new [7]. A structure that arises frequently in eigenvalue problems is symmetry: A = A^T . This very important structure does not translate into any obvious structure for the factors Q and R, so it does not fit into our framework in an obvious way. In Section 4.4 we show that the symmetric problem can be solved using our methodology: we turn it into a unitary problem by a Cayley transform [10]. We did not seriously expect this approach would be faster than all of the many other existing methods for the symmetric eigenvalue problem [73, § 7.2], but we were pleased to find that it is not conspicuously slow. Our solution to the symmetric problem serves as a stepping stone to the solution of the symmetric-plus-rank-one problem, which includes comrade and colleague matrices [14] as important special cases. If a polynomial p is presented as a linear combination of Chebyshev or Legendre polynomials, for example, the coefficients can be placed into a comrade matrix with eigenvalues equal to the zeros of p. A Cayley transform turns the symmetric-plus-rank-one problem into a unitary-plus-rank-one problem, which we can solve by our fast companion solver. Our O(n^2) algorithm gives us a fast way