Matrix Algorithms

The front matter includes the title page, copyright page, TOC, list of algorithms, and preface

In 1965 James Hardy Wilkinson published his classic work on numerical linear algebra, The Algebraic Eigenvalue Problem. His purpose was to survey the computational state of the art, but a significant part of the book was devoted to the mathematical foundations, some of which had been developed by Wilkinson himself. This happy mixture of the theoretical and computational has continued to characterize the area—computational advances go hand in hand with deeper mathematical understanding.

The purpose of this chapter is to provide the mathematical background for the first half of this volume. The first section consists of an exposition of the classic algebraic theory of eigensystems. The predominant theme of this section is the systematic reduction of matrices to simpler forms by similarity transformations. One of these forms — the Schur decomposition — is the goal of the QR algorithm, the workhorse of dense eigenvalue computations.

The second and third sections are devoted to the analysis of eigensystems. The second section concerns norms and the their relation to the spectrum of a matrix. Here we will pay special attention to the asymptotic behavior of powers of matrices, which are important in the analysis of many matrix algorithms. The third section is devoted to perturbation theory — how eigensystems behave when the original matrix is perturbed. This is a large subject, and we will only be able to survey the basics.

The majority of eigenvalue problems involve real matrices, and their reality often results in computational savings. However, the theory of eigensystems is best treated in terms of complex matrices. Consequently, throughout this chapter:

Unless otherwise stated A is a complex matrix of order n.

1. The Algebra of Eigensystems

In this section we will develop the classical theory of eigenvalues, eigenvectors, and reduction by similarity transformations, with an emphasis on the reduction of a matrix by similarity transformations to a simpler form. We begin with basic definitions and their consequences.

All general-purpose eigenvalue algorithms are necessarily iterative. This statement is a consequence of Abel's famous proof that there is no algebraic formula for the roots of a general polynomial of degree greater than four. Specifically, to any polynomial p(x) = xⁿ − a_n−1xⁿ⁻¹ ⋯ − a₁x − a₀ there corresponds a companion matrix Cp = ( an−1 an−2 ⋯ a1 a0 10⋯0001⋯00⋮⋮⋮⋮00⋯0000⋯10) whose characteristic polynomial is p. Thus if we had a noniterative general-purpose eigenvalue algorithm, we could apply it to the companion matrix of an arbitrary polynomial to find its zeros in a finite number of steps — in effect producing a formula for the roots.

Thus the history of eigenvalue computations is the story of iterative methods for approximating eigenvalues and eigenvectors. An important theme in this story concerns methods based on powers of a matrix. The climax of that theme is the QR algorithm — the great success story of modern matrix computations. This algorithm, with proper attention to detail, can compute a Schur form of a general matrix in O(n³) operations. It can be adapted to reduce a real matrix in real arithmetic to a real variant of the Schur form. Other forms of the algorithm compute the eigenvalues of and eigenvectors of a Hermitian matrix, the singular value decomposition of a general matrix, and the generalized Schur form of a matrix pencil. And the underlying theory of the method continues to suggest new algorithms.

The purpose of this chapter is to describe the QR algorithm and its variants. The algorithm is closely related to two algorithms that are important in their own right — the power method and the inverse power method. We will therefore begin with an exposition of those methods.

We have seen that symmetric (or Hermitian) matrices have special mathematical properties. In particular, their eigenvalues are real and their eigenvectors can be chosen to form an orthonormal basis for ℝⁿ. They also have special computational properties. Because of their symmetry they can be stored in about half the memory required for a general matrix. Moreover, we can frequently use the symmetry of the problem to reduce the operation count in certain algorithms. For example, the Cholesky algorithm for triangularizing a symmetric matrix requires half the operations that are required by its nonsymmetric counterpart, Gaussian elimination (see Algorithm I:2.2.1).

The first part of this chapter is devoted to efficient algorithms for the symmetric eigenvalue problem. In the first section we will show how to adapt the QR algorithm to solve the symmetric eigenvalue problem. The procedure begins as usual with a reduction to a simpler form — in this case tridiagonal form — followed by a QR iteration on the reduced form. We will treat this reduction in §1.1 and the tridiagonal QR iteration itself in §1.2.

In §2 we will consider some alternative methods for the symmetric eigenvalue problem: a method for updating the symmetric eigendecomposition, a divide-and- conquer algorithm, and methods for finding eigenvalues and eigenvectors of band matrices.

We next turn to the singular value decomposition. It is included in this chapter because the squares of the singular values of a real matrix X are the eigenvalues of the symmetric matrix A = X^TX and the eigenvectors of A are the left singular vectors of X. Thus, in principle, the algorithms of the first section of this chapter can be used to solve the singular value problem. In practice, this approach has numerical difficulties (see §3.2), and it is better to solve the singular value problem on its own terms. Fortunately, the QR algorithm for symmetric matrices has an analogue for the singular value problem.

Finally, in a brief coda, we treat the symmetric positive definite generalized eigenvalue problem Ax = λBx, where A is symmetric and B is positive definite. This problem has appealing analogies to the ordinary symmetric eigenproblem, but it is mathematically and numerically less tractable.

The first half of this book has been concerned with dense matrices whose elements can be stored in the main memory of a computer. We now shift focus to problems that are so large that storage limitations become a problem. This problem has two aspects.

First, if the order of the matrix in question is sufficiently large, we will not be able to store all its elements. This limitation is not as burdensome as it might appear. Many large matrices arising in applications are sparse; that is, their zero elements greatly outnumber their nonzero elements. With appropriate data structures we can get away with storing only the latter. Unfortunately, these data structures restrict the operations we can perform on the matrix, and, in consequence, most algorithms for large sparse eigenvalue problems confine themselves to matrix-vector multiplications or perhaps the computation of an LU factorization.

Second, the eigenvectors of a sparse matrix are in general not sparse. It follows that if the matrix is large enough we cannot afford to store its entire eigensystem. Instead we must be satisfied with computing some subset of its eigenpairs. The space spanned by a set of eigenvectors is an example of an eigenspace, and it turns out that algorithms for large matrices are best explained in terms of eigenspaces. Consequently, this chapter is devoted to the basic facts about eigenspaces and their approximations.

We begin with the basic algebraic theory of eigenspaces, and in §2 we present the analytic theory — residual analysis and perturbation bounds. Because eigenspaces, like eigenvectors, cannot be computed directly, most algorithms for large problems proceed by generating a sequence of subspaces that contain increasingly accurate approximations to the desired eigenvectors. In §3 we treat one of the most popular ways of generating subspaces — the Krylov sequence — and in §4 we analyze the Rayleigh-Ritz method for extracting individual eigenvectors from a subspace.

Although all of the results of this chapter apply to matrices of any size, our chief concern is with large eigenproblems. Therefore:

Unless otherwise specified A will denote a matrix of order n. It will often help to think of n as being very large.

In this chapter and the next we will treat algorithms for large eigenproblems. There are two ways of presenting this subject. The first is by the kind of problem — the symmetric eigenproblem, the generalized eigenproblem, etc. This approach has obvious advantages for the person who is interested in only one kind of problem. It has the disadvantage that, at the very least, a thumbnail sketch of each algorithm must accompany each problem type. The alternative is to partition the subject by algorithms, which is what we will do here.

In the present chapter we will be concerned with methods based on Krylov sequences—specifically, Arnoldi's method and the Lanczos method. In the first subsection we will introduce certain decompositions associated with Krylov subspaces. The decompositions form the basis for the the Arnoldi and Lanczos algorithms to be treated in §2 and §3.

1. Krylov Decompositions

For a fixed matrix A, not every subspace of dimension k is a Krylov subspace of A. This means that any basis of a Krylov subspace must satisfy certain constraining relations, which we will call Krylov decompositions. It turns out that a knowledge of these decompositions enables us to compute and manipulate Krylov subspaces. Consequently, this section is devoted to describing the relations that must hold for a basis to span a Krylov subspace.

The most widely used Krylov decompositions are those that result from orthogonalizing a Krylov sequence. They are called Arnoldi decompositions for non-Hermitian matrices and Lanczos decompositions for Hermitian matrices. We will treat these special decompositions in the first two subsections. Unfortunately, these decompositions are not closed under certain natural operations, and in §1.3 we introduce more general decompositions that are closed — decompositions that will be used in our computational algorithms.

Although the most widely used algorithms for large eigenproblems are based on Krylov sequences, there are other approaches. In this chapter we will consider two alternatives: subspace iteration and methods based on Newton's method — particularly the Jacobi-Davidson method.

Subspace iteration is a block generalization of the power method. It is an older method and is potentially slow. But it is very reliable; its implementation is comparatively simple; and it is easy to use. Moreover, combined with shift-and-invert enhancement or Chebyshev acceleration, it sometimes wins the race.

We next turn to methods based on Newton's method for solving systems of nonlinear equations. Since Newton's method requires the solution of a linear system, these methods appear to be mostly suitable for dense problems. But with some care they can be adapted to large sparse eigenprobems. When we combine Newton's method with Rayleigh-Ritz refinement, we get the Jacobi-Davidson algorithm, with which this chapter and volume conclude.

1. Subspace Iteration

In motivating the power method with starting vector U₀ (§1.1, Chapter 2) we assumed that A was diagonalizable and expanded Ak u0 = γ1 λ 1 k x1 + γ2 λ 2 k x2 + γ3 λ 3 k x3 +⋯ 1.1 in terms of the eigenvectors X_J of A. If | λ1 | >| λ2 |≥| λ3 |≥⋯≥| λn | and γ₁ ≠ 0, then the first term in the expansion dominates, and A^ku′₀, suitably scaled, converges to x₁.

When | λ₂| is near | λ₁ |, the converges of the power method will be slow.

The reader of this volume is assumed to know the basics of linear algebra and analysis as well as the behavior of matrix algorithms on digital computers. This background has been presented in full detail in Volume I of this work. Here we confine ourselves to a review of some of the more important notation and facts.

1. Scalars, vectors, and matrices

1. Scalars are real or complex numbers. The set of real numbers will be denoted by ℝ the set of complex numbers by ℂ. All scalars are denoted by lowercase Latin or Greek letters.

2. An n -vector is a collection of n scalars arranged in a column. The set of real n-vectors is denoted by ℝⁿ; the set of complex n-vectors by ℂⁿ. Vectors are denoted by lowercase Latin letters.

3. The scalars that form a vector are called the components of the vector. If a vector is denoted by a letter, its components will be denoted by a corresponding lowercase Latin or Greek letter. Thus the vector a may be written a= ( a1 a2 ⋮ aπ ) =( α1 α2 ⋮ απ ). Figure A. 1 gives the Greek alphabet and its corresponding Latin letters (some of the correspondences are quite artificial).

4. An m×n matrix is a rectangular array of mn scalars having m rows and n columns. The set of real m×n matrices is denoted by ℝ^m×n; the set of complex m×n matrices by ℂ^m×n. Matrices are denoted by uppercase Latin and Greek letters.

5. The scalars in a matrix are called the elements of the matrix.

The back matter includes references and index