Divide-and-conquer eigenvalue algorithm

As with most eigenvalue algorithms for Hermitian matrices, divide-and-conquer begins with a reduction to tridiagonal form. For an ${\displaystyle m\times m}$ matrix, the standard method for this, via Householder reflections, takes ${\displaystyle {\frac {4}{3}}m^{3}}$ floating point operations, or ${\displaystyle {\frac {8}{3}}m^{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

${\displaystyle T={\begin{bmatrix}T_{1}&0\\0&T_{2}\end{bmatrix}}.}$

The eigenvalues and eigenvectors of ${\displaystyle T}$ are simply those of ${\displaystyle T_{1}}$ and ${\displaystyle T_{2}}$, 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 ${\displaystyle m\times m}$ real symmetric tridiagonal matrix ${\displaystyle T}$. The algorithm can be modified for Hermitian matrices.

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 ${\displaystyle T_{1}}$ we will call ${\displaystyle n\times n}$, and then ${\displaystyle T_{2}}$ is ${\displaystyle (m-n)\times (m-n)}$. ${\displaystyle T}$ is almost block diagonal regardless of how ${\displaystyle n}$ is chosen. For efficiency we typically choose ${\displaystyle n\approx m/2}$.

We write ${\displaystyle T}$ as a block diagonal matrix, plus a rank-1 correction:

The only difference between ${\displaystyle T_{1}}$ and ${\displaystyle {\hat {T}}_{1}}$ is that the lower right entry ${\displaystyle t_{nn}}$ in ${\displaystyle {\hat {T}}_{1}}$ has been replaced with ${\displaystyle t_{nn}-\beta }$ and similarly, in ${\displaystyle {\hat {T}}_{2}}$ the top left entry ${\displaystyle t_{n+1,n+1}}$ has been replaced with ${\displaystyle t_{n+1,n+1}-\beta }$.

The remainder of the divide step is to solve for the eigenvalues (and if desired the eigenvectors) of ${\displaystyle {\hat {T}}_{1}}$ and ${\displaystyle {\hat {T}}_{2}}$, that is to find the diagonalizations ${\displaystyle {\hat {T}}_{1}=Q_{1}D_{1}Q_{1}^{T}}$ and ${\displaystyle {\hat {T}}_{2}=Q_{2}D_{2}Q_{2}^{T}}$. This can be accomplished with recursive calls to the divide-and-conquer algorithm, although practical implementations often switch to the implicitly shifted QR algorithm for small enough submatrices.^[1]

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 ${\displaystyle z^{T}=(q_{1}^{T},q_{2}^{T})}$, where ${\displaystyle q_{1}^{T}}$ is the last row of ${\displaystyle Q_{1}}$ and ${\displaystyle q_{2}^{T}}$ is the first row of ${\displaystyle Q_{2}}$. It is now elementary to show that

${\displaystyle T={\begin{bmatrix}Q_{1}&\\&Q_{2}\end{bmatrix}}\left({\begin{bmatrix}D_{1}&\\&D_{2}\end{bmatrix}}+\beta zz^{T}\right){\begin{bmatrix}Q_{1}^{T}&\\&Q_{2}^{T}\end{bmatrix}}}$

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 ${\displaystyle D+ww^{T}}$, where ${\displaystyle D}$ is diagonal with distinct entries and ${\displaystyle w}$ is any vector with nonzero entries. In this case ${\displaystyle w={\sqrt {|\beta |}}\cdot z}$.

The case of a zero entry is simple, since if w_i is zero, (${\displaystyle e_{i}}$,d_i) is an eigenpair (${\displaystyle e_{i}}$ is in the standard basis) of ${\displaystyle D+ww^{T}}$ since ${\displaystyle (D+ww^{T})e_{i}=De_{i}=d_{i}e_{i}}$.

If ${\displaystyle \lambda }$ is an eigenvalue, we have:

${\displaystyle (D+ww^{T})q=\lambda q}$

where ${\displaystyle q}$ is the corresponding eigenvector. Now

${\displaystyle (D-\lambda I)q+w(w^{T}q)=0}$

${\displaystyle q+(D-\lambda I)^{-1}w(w^{T}q)=0}$

${\displaystyle w^{T}q+w^{T}(D-\lambda I)^{-1}w(w^{T}q)=0}$

Keep in mind that ${\displaystyle w^{T}q}$ is a nonzero scalar. Neither ${\displaystyle w}$ nor ${\displaystyle q}$ are zero. If ${\displaystyle w^{T}q}$ were to be zero, ${\displaystyle q}$ would be an eigenvector of ${\displaystyle D}$ by ${\displaystyle (D+ww^{T})q=\lambda q}$. If that were the case, ${\displaystyle q}$ would contain only one nonzero position since ${\displaystyle D}$ is distinct diagonal and thus the inner product ${\displaystyle w^{T}q}$ can not be zero after all. Therefore, we have:

${\displaystyle 1+w^{T}(D-\lambda I)^{-1}w=0}$

or written as a scalar equation,

${\displaystyle 1+\sum _{j=1}^{m}{\frac {w_{j}^{2}}{d_{j}-\lambda }}=0.}$

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.

Solving the nonlinear secular equation can be done using an iterative technique, such as the Newton–Raphson method. However, each root can be found in O(1) iterations, each of which requires ${\displaystyle \Theta (m)}$ flops (for an ${\displaystyle m}$-degree rational function), making the cost of the iterative part of this algorithm ${\displaystyle \Theta (m^{2})}$. The fast multipole method has also been employed to solve the secular equation in ${\displaystyle \Theta (m\log(m))}$ operations.^[2][1]