Glued Matrices and the MRRR Algorithm

Abstract

Abstract. During the last ten years, Dhillon and Parlett devised a new algorithm (multiple relatively robust representations (MRRR)) for computing numerically orthogonal eigenvectors of a symmetric tridiagonal matrix T with O(n 2 ) cost. It has been incorporated into LAPACK version 3.0 as routine stegr.We have discovered that the MRRR algorithm can fail in extreme cases. Sometimes eigenvalues agree to working accuracy and MRRR cannot compute orthogonal eigenvectors for them. In this paper, we describe and analyze these failures and various remedies. 1. Introduction. Starting in the mid-90s, Dhillon and Parlett developed the algorithm of multiple relatively robust representations (MRRR) that computes numerically orthogonal eigenvectors of a symmetric tridiagonal matrix T with O(n 2 ) cost [13,15,16,17,7]. This algorithm was tested on a large and challenging set of matrices and has been incorporated into LAPACK version 3.0 as routine stegr. The algorithm is described, in the detail we need here, in section 2. For a more detailed description see also [8]. An example of MRRR in action is given in section 3.In 2003, one of us (Vömel) came to Berkeley to assist in the modification of stegr to compute a subset of k eigenpairs with O(kn) operations. When testing stegr on more and more challenging matrices, he discovered cases of failure.Investigation of these cases brought to light assumptions made in stegr that hold in exact arithmetic and in the majority of cases in finite precision arithmetic, but that can fail. These assumptions, why they are reasonable, and how they can fail are the subject of section 4.In section 5, we propose and analyze various remedies for the aforementioned shortcomings and show how to incorporate them into MRRR. We also look at the cost of these modifications. We select as most suitable an approach that is based on small random perturbations and introduces artificial roundoff effects. This approach preserves the complexity of the original algorithm.

Experiences with a Lanczos Eigensolver in High-Precision Arithmetic

“…The idea of perturbing problems in order to avoid such singularities is well-established in existing research [4,6,28].…”

Performance and Accuracy of LAPACK's Symmetric Tridiagonal Eigensolvers

“…Moreover, Wilkinson and glued Wilkinson matrices are known to be notoriously difficult for MR due to strongly clustered eigenvalues [13]. Thus we expect to see MR perform poorly on these classes of matrices.…”

“…The results for the latter class are shown in Figure 4.11. The difficulties of MR for these matrices are well understood: since the eigenvalues of glued matrices come in groups of small size but extreme tightness, the representation tree generated by the MR algorithm is very broad and the overhead for the tree generation is considerable, see [13,7]. On top of the difficulties of MR, the fraction deflated in DC is ∈ [59%, 78%], that is DC is extraordinarily efficient and even faster than for practical matrices.…”

“…This class includes distributions that are already used in LAPACK's tester, synthetic distributions, and matrices that are used in [8]. Furthermore, it includes Wilkinson matrices [31,23] and glued Wilkinson matrices (see for example [13]). A detailed description of these matrices is given later, in Table 4.3 of Section 4.3.…”

“…If some eigenvalues of T agree to d digits on average, then the algorithm has to do work proportional to dn 2 . The algorithm uses a random perturbation to ensure with high probability that eigenvalues cannot be too strongly clustered, see [13] for details. MR cannot make use of higher level BLAS.…”