SuperDC: Stable superfast divide-and-conquer eigenvalue decomposition

For dense Hermitian matrices with small off-diagonal (numerical) ranks and in a hierarchically semiseparable form, we give a stable divide-and-conquer eigendecomposition method with nearly linear complexity (called SuperDC) that significantly improves an earlier basic algorithm in [Vogel, Xia, et al., SIAM J. Sci. Comput., 38 (2016)]. We incorporate a sequence of key stability techniques and provide many improvements in the algorithm design. Various stability risks in the original algorithm are analyzed, including potential exponential norm growth, cancellations, loss of accuracy with clustered eigenvalues or intermediate eigenvalues, etc. In the dividing stage, we give a new structured low-rank update strategy with balancing that eliminates the exponential norm growth and also minimizes the ranks of low-rank updates. In the conquering stage with low-rank updated eigenvalue solution, the original algorithm directly uses the regular fast multipole method (FMM) to accelerate function evaluations, which has the risks of cancellation, division by zero, and slow convergence. Here, we design a triangular FMM to avoid cancellation. Furthermore, when there are clustered intermediate eigenvalues or when updates to existing eigenvalues are very small, we design a novel local shifting strategy to integrate FMM accelerations into the solution of shifted secular equations so as to achieve both the efficiency and the reliability. We also provide several improvements or clarifications on some structures and techniques that are missing or unclear in the previous work. The resulting SuperDC eigensolver has significantly better stability while keeping the nearly linear complexity for finding the entire eigenvalue decomposition. In a set of comprehensive tests, SuperDC shows dramatically lower runtime and storage than the Matlab eig function. The stability benefits are also confirmed with both analysis and numerical comparisons.