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.
Background.
As with most eigenvalue algorithms for Hermitian matrices, divide-and-conquer begins with a reduction to tridiagonal form. For an formula_1 matrix, the standard method for this, via Householder reflections, takes formula_2 floating point operations, or formula_3 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 formula_5 are simply those of formula_6 and formula_7, 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 formula_1 real symmetric tridiagonal matrix formula_5. Although the algorithm can be modified for Hermitian matrices, we do not give the details here.
Divide.
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 formula_6 we will call formula_11, and then formula_7 is formula_13. Note that the remark about formula_5 being almost block diagonal is true regardless of how formula_15 is chosen (i.e., there are many ways to so decompose the matrix). However, it makes sense, from an efficiency standpoint, to choose formula_16. 
We write formula_5 as a block diagonal matrix, plus a rank-1 correction:
The only difference between formula_6 and formula_19 is that the lower right entry formula_20 in formula_19 has been replaced with formula_22 and similarly, in formula_23 the top left entry formula_24 has been replaced with formula_25.
The remainder of the divide step is to solve for the eigenvalues (and if desired the eigenvectors) of formula_19 and formula_23, that is to find the diagonalizations formula_28 and formula_29. 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.
Conquer.
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 formula_30, where formula_31 is the last row of formula_32 and formula_33 is the first row of formula_34. 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 formula_36, where formula_37 is diagonal with distinct entries and formula_38 is any vector with nonzero entries.
The case of a zero entry is simple, since if wi is zero, (formula_39,di) is an eigenpair (formula_39 is in the standard basis) of formula_36 since
formula_42.
If formula_43 is an eigenvalue, we have:
where formula_45 is the corresponding eigenvector. Now
Keep in mind that formula_49 is a nonzero scalar. Neither formula_38 nor formula_45 are zero. If formula_49 were to be zero, formula_45 would be an eigenvector of formula_37 by formula_44. If that were the case, formula_45 would contain only one nonzero position since formula_37 is distinct diagonal and thus the inner product formula_49 can not be zero after all. Therefore, we have:
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 formula_61 flops (for an formula_62-degree rational function), making the cost of the iterative part of this algorithm formula_63.
Analysis.
As is common for divide and conquer algorithms, we will use the master theorem for divide-and-conquer recurrences to analyze the running time. Remember that above we stated we choose formula_16. We can write the recurrence relation:
In the notation of the Master theorem, formula_66 and thus formula_67. Clearly, formula_68, so we have
Remember that above we pointed out that reducing a Hermitian matrix to tridiagonal form takes formula_2 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 formula_63 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 formula_3, but the second part of the algorithm takes formula_73 as well. For the QR algorithm with a reasonable target precision, this is formula_74, whereas for divide-and-conquer it is formula_75. The reason for this improvement is that in divide-and-conquer, the formula_73 part of the algorithm (multiplying formula_77 matrices) is separate from the iteration, whereas in QR, this must occur in every iterative step. Adding the formula_3 flops for the reduction, the total improvement is from formula_79 to formula_80 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 formula_77 and the vectors formula_82 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.
The divide-and-conquer algorithm is readily parallelized, and linear algebra computing packages such as LAPACK contain high-quality parallel implementations.