The block version of the classical Gram--Schmidt (BCGS) method is often employed to efficiently compute orthogonal bases for Krylov subspace methods and eigenvalue solvers, but a rigorous proof of its stability behavior has not yet been established. It is shown that the usual implementation of BCGS can lose orthogonality at a rate worse than O(epsilon)kappa(2)(chi), where chi is the input matrix and epsilon is the unit roundoff.
A useful intermediate quantity denoted as the Cholesky residual is given special attention and, along with a block generalization of the Pythagorean theorem, this quantity is used to develop more stable variants of BCGS. These variants are proven to have O(epsilon)kappa(2)(chi) loss of orthogonality with relatively relaxed conditions on the intrablock orthogonalization routine satisfied by the most commonly used algorithms.
A variety of numerical examples illustrate the theoretical bounds.