Motivation for computing pseudospectra of sparse matrices

The matrices in many eigenvalue problems are too large to allow direct computation of their full spectra, and two of the iterative tools available for computing a part of the spectrum are ARPACK and its MATLAB counterpart, eigs. For non-symmetric matrices, the mathematical basis of these packages is the Arnoldi iteration with implicit restarting, which works by compressing the matrix to an `interesting' Hessenberg matrix, one which contains information about the eigenvalues and eigenvectors of interest.

For some matrices, non-normality (non-orthogonality of the eigenvectors) may be physically important. In extreme cases, non-normality combined with the practical limits of machine precision can lead to difficulties in accurately finding the eigenvalues. Perhaps the more common and more important situation is when the non-normality is pronounced enough to limit the physical significance of eigenvalues for applications, without rendering them uncomputable. In applications, users need to know if they are in such a situation. The prevailing practice in large-scale eigenvalue calculations is that users get no information of this kind.

There is a familiar tool available for learning more about the cases in which non-normality may be important: pseudospectra. One such case is described in our second example, which shows some of the pseudospectra of the `Grcar matrix' of dimension 400, the exact spectrum, and converged eigenvalue estimates (Ritz values) returned by a run of ARPACK for this matrix. This is an extreme example where the non-normality is so pronounced that even with the convergence tolerance set to its lowest possible value, machine epsilon, the eigenvalue estimates are far from the true spectrum. From the Ritz values alone, one might not realise that anything was amiss. Once the pseudospectra are plotted too, it is obvious.

Computing the pseudospectra of a matrix of dimension N by the obvious method is an expensive task, requiring an O(N³) singular value decomposition at each point in a grid. For a reasonably fine mesh, this leads to an O(N³) algorithm with the constant of the order of thousands. Recent developments in algorithms for computing pseudospectra have improved the constant, and the asymptotic complexity for large sparse matrices, but these are still fairly costly techniques. The GUI can cheaply compute an approximation to the pseudospectra of a large matrix in a region near the interesting eigenvalues. Our method uses the upper Hessenberg matrix constructed after successive iterations of the implicitly restarted Arnoldi algorithm, as implemented in ARPACK. Among other things, this means that after performing an eigenvalue computation with eigs, a user can quickly obtain a graphical check to indicate whether the Ritz values returned are likely to be physically meaningful. Our vision is that every ARPACK or eigs user ought to plot pseudospectra estimates routinely after their eigenvalue computations as a cheap `sanity check'.

For more discussion see the

Pseudospectra GUI home page.