The Godunov-Inverse Iteration: A Fast and Accurate Solution to the Symmetric Tridiagonal Eigenvalue Problem

Anna M. Matsekh

Institute of Computational Technologies, Siberian Branch of the Russian Academy of Sciences, Lavrentiev Ave. 6, Novosibirsk 630090, Russia.


Abstract

We present a new hybrid algorithm based on Godunov's method for computing eigenvectors of symmetric tridiagonal matrices and Inverse Iteration, which we call the Godunov-Inverse Iteration. Godunov's method [Godunov, S.K., Antonov, A.G., Kiriljuk, O.P., Kostin, V.I., 1993. Guaranteed accuracy in numerical linear algebra. Kluwer Academic Publishers Group, Dordrecht, translated and revised from the 1988 Russian original.] uses two-sided Sturm sequences to compute eigenvectors of a tridiagonal symmetric matrix with guaranteed accuracy in N floating point operations per an eigenvector. The use of two-sided Sturm sequences avoids loss of accuracy associated with rounding errors that is encountered in methods based on conventional Sturm sequences. The algorithm gives provably accurate solutions to the symmetric tridiagonal unreduced eigenvalue problem using a specially designed floating point arithmetic with extended precision and directed rounding. Unfortunately, in IEEE arithmetic, eigenvectors computed according to Godunov's algorithm are susceptible to division by zero and overflow errors, while for closely clustered interior eigenvalues it produces nearly collinear eigenvectors, taking no measures for reorthogonalization. In empirical studies Godunov's method, which is a direct method by it's nature, consistently delivers residuals that are approximately two orders of magnitude larger than those of the eigenvectors computed according to some of the Inverse Iteration implementations for tridiagonal symmetric matrices such as LAPCK implementation of the Inverse Iteration xSTEIN.

To overcome shortcomings of both Godunov's method and Inverse Iteration (often nondeterministic character of starting vectors, need to introduce disturbances into the shift, high worst case complexity) we constructed a new hybrid procedure -- the Godunov-Inverse Iteration. Changing any nonnumeric components of the Godunov eigenvectors to random uniformly distributed numbers, we apply Inverse Iteration to these vectors, which usually achieves desired error bounds in one step. This contrasts with other implementations of the Inverse Iteration algorithm, which require a few more steps to achieve the same accuracy. This is most advantageous in the case of closely clustered eigenvalues when a large fraction of the eigenvectors require reorthogonalization. Godunov-Inverse Iteration is very robust with respect to the choice of the Inverse Iteration shift -- we use right-hand bounds of the eigenvalue intervals computed by the bisection method as extremely accurate shifts in the Godunov-Inverse Iteration. We resort to reorthogonalization within the iteration only in cases of computationally coincident or closely clustered eigenvalues. As a result Godunov-Inverse Iteration produces accurate and robust solutions to the symmetric eigenvalue problem with higher accuracy than Godunov's method and in fewer steps than existing implementations of the Inverse Iteration.