Analysis of the Cholesky method with iterative refinement for solving the symmetric definite generalized eigenproblem

A standard method for solving the symmetric definite generalized eigenvalue problem Ax = λBx, where A is symmetric and B is symmetric positive definite, is to compute a Cholesky factorization B = LLT (optionally with complete pivoting) and solve the equivalent standard symmetric eigenvalue problem Cy = λy, where C = L-1AL-T. Provided that a stable eigensolver is used, standard error analysis says that the computed eigenvalues are exact for A + ΔA and B + ΔB with max(∥ΔA∥2/∥A∥2, ∥ΔB∥2/∥B∥2) bounded by a multiple of κ2(B)u, where u is the unit roundoff. We take the Jacobi method as the eigensolver and give a detailed error analysis that yields backward error bounds potentially much smaller than κ2(B)u. To show the practical utility of our bounds we describe a vibration problem from structural engineering in which B is ill conditioned yet the error bounds are small. We show how, in cases of instability, iterative refinement based on Newton's method can be used to produce eigenpairs with small backward errors. Our analysis and experiments also give insight into the popular Cholesky-QR method, in which the QR method is used as the eigensolver. We argue that it is desirable to augment current implementations of this method with pivoting in the Cholesky factorization.