Divide-and-conquer eigenvalue algorithm
Divide-and-conquer eigenvalue algorithms are a class of eigenvalue algorithms for Hermitian or real symmetric matrices that have recently (circa 1990s) become competitive in terms of stability and efficiency with more traditional algorithms such as the QR algorithm. The basic concept behind these algorithms is the divide-and-conquer approach from computer science. An eigenvalue problem is divided into two problems of roughly half the size, each of these are solved recursively, and the eigenvalues of the original problem are computed from the results of these smaller problems.
Here we present the simplest version of a divide-and-conquer algorithm, similar to the one originally proposed by Cuppen in 1981. Many details that lie outside the scope of this article will be omitted; however, without considering these details, the algorithm is not fully stable.
As with most eigenvalue algorithms for Hermitian matrices, divide-and-conquer begins with a reduction to tridiagonal form. For an matrix, the standard method for this, via Householder reflections, takes flops, or if eigenvectors are needed as well. There are other algorithms, such as the Arnoldi iteration, which may do better for certain classes of matrices; we will not consider this further here.
In certain cases, it is possible to deflate an eigenvalue problem into smaller problems. Consider a block diagonal matrix
The eigenvalues and eigenvectors of T are simply those of T1 and T2, and it will almost always be faster to solve these two smaller problems than to solve the original problem all at once. This technique can be used to improve the efficiency of many eigenvalue algorithms, but it has special significance to divide-and-conquer.
For the rest of this article, we will assume the input to the divide-and-conquer algorithm is an real symmetric tridiagonal matrix T. Although the algorithm can be modified for Hermitian matrices, we do not give the details here.
The divide part of the divide-and-conquer algorithm comes from the realization that a tridiagonal matrix is "almost" block diagonal.
The size of submatrix T1 we will call , and then T2 is . Note that the remark about T being almost block diagonal is true regardless of how n is chosen (i.e., there are many ways to so decompose the matrix). However, it makes sense, from an efficiency standpoint, to choose .
We write T as a block diagonal matrix, plus a rank-1 correction:
The only difference between T1 and is that the lower right entry tnn in has been replaced with tnn − β and similarly, in the top left entry tn + 1,n + 1 has been replaced with tn + 1,n + 1 − β.
The remainder of the divide step is to solve for the eigenvalues (and if desired the eigenvectors) of and , that is to find the diagonalizations and . This can be accomplished with recursive calls to the divide-and-conquer algorithm, although practical implementations often switch to the QR algorithm for small enough submatrices.
The conquer part of the algorithm is the unintuitive part. Given the diagonalizations of the submatrices, calculated above, how do we find the diagonalization of the original matrix?
First, define , where is the last row of Q1 and is the first row of Q2. It is now elementary to show that
The remaining task has been reduced to finding the eigenvalues of a diagonal matrix plus a rank-one correction. Before showing how to do this, let us simplify the notation. We are looking for the eigenvalues of the matrix D + wwT, where D is diagonal with distinct entries and w is any vector with nonzero entries.
If wi is zero, (ei,di) is an eigenpair of D + wwT since (D + wwT)ei = Dei = diei.
If λ is an eigenvalue, we have:
- (D + wwT)q = λq
where q is the corresponding eigenvector. Now
- (D − λI)q + w(wTq) = 0
- q + (D − λI) − 1w(wTq) = 0
- wTq + wT(D − λI) − 1w(wTq) = 0
Keep in mind that wTq is a nonzero scalar. Neither w nor q are zero. If wTq were to be zero, q would be an eigenvector of D by (D + wwT)q = λq. If that were the case, q would contain only one nonzero position since D is distinct diagonal and thus the inner product wTq can not be zero after all. Therefore, we have:
- 1 + wT(D − λI) − 1w = 0
or written as a scalar equation,
This equation is known as the secular equation. The problem has therefore been reduced to finding the roots of the rational function defined by the left-hand side of this equation.
All general eigenvalue algorithms must be iterative, and the divide-and-conquer algorithm is no different. Solving the nonlinear secular equation requires an iterative technique, such as the Newton–Raphson method. However, each root can be found in O(1) iterations, each of which requires Θ(m) flops (for an m-degree rational function), making the cost of the iterative part of this algorithm Θ(m2).
In the notation of the Master theorem, a = b = 2 and thus log ba = 1. Clearly, Θ(m2) = Ω(m1), so we have
- T(m) = Θ(m2)
Remember that above we pointed out that reducing a Hermitian matrix to tridiagonal form takes flops. This dwarfs the running time of the divide-and-conquer part, and at this point it is not clear what advantage the divide-and-conquer algorithm offers over the QR algorithm (which also takes Θ(m2) flops for tridiagonal matrices).
The advantage of divide-and-conquer comes when eigenvectors are needed as well. If this is the case, reduction to tridiagonal form takes , but the second part of the algorithm takes Θ(m3) as well. For the QR algorithm with a reasonable target precision, this is , whereas for divide-and-conquer it is . The reason for this improvement is that in divide-and-conquer, the Θ(m3) part of the algorithm (multiplying Q matrices) is separate from the iteration, whereas in QR, this must occur in every iterative step. Adding the flops for the reduction, the total improvement is from to flops.
Practical use of the divide-and-conquer algorithm has shown that in most realistic eigenvalue problems, the algorithm actually does better than this. The reason is that very often the matrices Q and the vectors z tend to be numerically sparse, meaning that they have many entries with values smaller than the floating point precision, allowing for numerical deflation, i.e. breaking the problem into uncoupled subproblems.
Variants and implementation
The algorithm presented here is the simplest version. In many practical implementations, more complicated rank-1 corrections are used to guarantee stability; some variants even use rank-2 corrections.
There exist specialized root-finding techniques for rational functions that may do better than the Newton-Raphson method in terms of both performance and stability. These can be used to improve the iterative part of the divide-and-conquer algorithm.
- Demmel, James W. (1997), Applied Numerical Linear Algebra, Philadelphia, PA: Society for Industrial and Applied Mathematics, ISBN 0-89871-389-7, MR1463942 .
Numerical linear algebra Key concepts Problems Hardware Software
Wikimedia Foundation. 2010.
Look at other dictionaries:
Divide and conquer — (derived from the Latin saying Divide et impera) may refer to: Divide and rule, in politics, sociology and economics, a strategy to gain or maintain power Defeat in detail, in warfare, a tactical maneuver to efficiently deal with a numerous… … Wikipedia
Eigenvalue algorithm — In linear algebra, one of the most important problems is designing efficient and stable algorithms for finding the eigenvalues of a matrix. These eigenvalue algorithms may also find eigenvectors. Contents 1 Characteristic polynomial 2 Power… … Wikipedia
Divide — Articles concerning Divide include: In Geography: Drainage divide or watershed, a ridge of land between two drainage basins Continental divide, a water divide between the drainage of two oceans Continental Divide, the North American divide… … Wikipedia
List of numerical analysis topics — This is a list of numerical analysis topics, by Wikipedia page. Contents 1 General 2 Error 3 Elementary and special functions 4 Numerical linear algebra … Wikipedia
List of mathematics articles (D) — NOTOC D D distribution D module D D Agostino s K squared test D Alembert Euler condition D Alembert operator D Alembert s formula D Alembert s paradox D Alembert s principle Dagger category Dagger compact category Dagger symmetric monoidal… … Wikipedia
Characteristic polynomial — This article is about the characteristic polynomial of a matrix. For the characteristic polynomial of a matroid, see Matroid. For that of a graded poset, see Graded poset. In linear algebra, one associates a polynomial to every square matrix: its … Wikipedia
Singular value decomposition — Visualization of the SVD of a 2 dimensional, real shearing matrix M. First, we see the unit disc in blue together with the two canonical unit vectors. We then see the action of M, which distorts the disk to an ellipse. The SVD decomposes M into… … Wikipedia
Numerical analysis — Babylonian clay tablet BC 7289 (c. 1800–1600 BC) with annotations. The approximation of the square root of 2 is four sexagesimal figures, which is about six decimal figures. 1 + 24/60 + 51/602 + 10/603 = 1.41421296... Numerical analysis is the … Wikipedia
Arrowhead matrix — In the mathematical field of linear algebra, an arrowhead matrix is a square matrix containing zeros in all entries except for the first row, first column, and main diagonal. [cite book |author=Stewart, G. W. |title=Matrix Algorithms… … Wikipedia