Pseudospectra¶
The \(p\)-norm \(\epsilon\)-pseudospectrum of a square matrix \(A\) is the set of all shifts \(\xi \in \mathbb{C}\) such that the inverse of \(A - \xi I\) is sufficiently large (greater than \(1/\epsilon\)) in the \(p\)-norm:
with the convention that \(\| (A-\xi I)^{-1} \|_p = +\infty\) if the resolvent, \((A-\xi I)^{-1}\), does not exist. In the most common case, where \(p=2\), the definition of the \(\epsilon\)-pseudospectrum reduces to
Elemental currently supports two methods for computing two-norm pseudospectra: the first is a high-performance improvement of the triangularization followed by inverse iteration approach of [VanLoan1984] and [Lui1997]. In particular, Elemental begins by computing the Schur decomposition of the given matrix, which preserves the (two-norm) \(\epsilon\)-pseudospectrum, up to round-off error, and then simultaneously performs many Implicitly Restarted Arnoldi (IRA) iterations with the inverse normal matrix for each shift in a manner which communicates no more data than a standard triangular solve with many right-hand sides (with vectors deflated immediately after convergence).
The second approach is quite similar and, instead of reducing to triangular form, reduces to Hessenberg form and performs multi-shift triangular solves in a manner similar to [Henry1994] and [BDP1996]. This algorithm is not yet performance-tuned in Elemental, but could perhaps be preferred in situations where a Schur decomposition is infeasible.
Elemental also supports proof-of-concept implementations of high-performance one-norm pseudospectra which are based upon the Hager-Higham algorithm ([Hager1984], [Higham1988], and [Higham1990]) for black-box one-norm estimation, but the blocked variant of [HighamTisseur2006] provides qualitatively superior results.
Please see [TrefethenEmbree2009] or [Trefethen1999] for comprehensive histories of the computation of pseudospectra.
References¶
- BDP1996
Christian H. Bischof, Biswa Nath Datta, and Avijit Purkayastha, A parallel algorithm for the Sylvester Observer Equation, SIAM Journal on Scientific Computing, Vol. 17, No. 3, pp. 686–698, 1996. DOI: http://dx.doi.org/10.1137/S1064827591223276
- Hager1984
William W. Hager, Condition estimates, SIAM Journal on Scientific and Statistical Computing, Vol. 5, No. 2, pp. 311–316, 1984. DOI: http://dx.doi.org/10.1137/0905023
- Henry1994
Greg Henry, The shifted Hessenberg system solve computation, Technical Report, Cornell University, Ithica, NY, 1994. Last accessed from: http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.33.7403
- Higham1988
Nicholas J. Higham, FORTRAN codes for estimating the one-norm of a real or complex matrix, with applications to condition estimation (Algorithm 674), ACM Transactions on Mathematical Software, Vol. 14, Issue 4, pp. 381–396, 1988. DOI: http://dx.doi.org/10.1145/50063.214386
- Higham1990
Nicholas J. Higham, Experience with a matrix norm estimator, SIAM Journal on Scientific and Statistical Computing, Vol. 11, No. 4, pp. 804–809, 1990. DOI: http://dx.doi.org/10.1137/0911047
- HighamTisseur2006
Nicholas J. Higham and Francoise Tisseur, A block algorithm for matrix 1-norm estimation, with an application to 1-norm pseudospectra, SIAM Journal on Matrix Analysis and Applications, Vol. 21, No. 4, pp. 1185–1201, 2006. DOI: http://dx.doi.org/10.1137/S0895479899356080
- Lui1997
S.H. Lui, Computation of pseudospectra by continuation, SIAM Journal on Scientific Computing, Vol. 18, No. 2, pp. 565–573, 1997. DOI: http://dx.doi.org/10.1137/S1064827594276035
- Trefethen1999
Lloyd N. Trefethen, Computation of pseudospectra, Acta Numerica, pp. 247–295, 1999. Last accessed from: http://eprints.maths.ox.ac.uk/1294/1/NA-99-03.pdf
- TrefethenEmbree2009
Lloyd N. Trefethen and Mark Embree, Spectra and pseudospectra: The behavior of nonnormal matrices, Princeton University Press, Princeton, 2009. Official website: http://press.princeton.edu/titles/8113.html
- VanLoan1984
Charles Van Loan, How near is a stable matrix to an unstable matrix?, Technical Report, Cornell University, Ithica, NY. Last accessed from: http://hdl.handle.net/1813/6488
C++ pseudospectra example driver
C++ ChunkedPseudospectra example driver