Large-Scale Computation of Pseudospectra Using ARPACK and Eigs.
ARPACK and its {\sc Matlab} counterpart, {\tt eigs}, are software packages that calculate some eigenvalues of a large nonsymmetric matrix by Arnoldi iteration with implicit restarts. We show that at a small additional cost, which diminishes relatively as the matrix dimension increases, good estimates of pseudospectra in addition to eigenvalues can be obtained as a by-product. Thus in large-scale eigenvalue calculations it is feasible to obtain routinely not just eigenvalue approximations, but also information as to whether or not the eigenvalues are likely to be physically significant. Examples are presented for matrices with dimension up to 200,000.
Introduction
. 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 [10, 11] and its Matlab counter-
part, eigs. 1 For nonsymmetric matrices, the mathematical basis of these packages is
the Arnoldi iteration with implicit restarting [11, 23], which works by compressing the
matrix to an "interesting" Hessenberg matrix, one which contains information about
the eigenvalues and eigenvectors of interest. For general information on large-scale
nonsymmetric matrix eigenvalue iterations, see [2, 21, 29, 31].
For some matrices, nonnormality (nonorthogonality of the eigenvectors) may be
physically important [30]. In extreme cases, nonnormality combined with the practical
limits of machine precision can lead to di#culties in accurately finding the eigenvalues.
Perhaps the more common and more important situation is when the nonnormality 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
nonnormality may be important: pseudospectra. Figure 1 shows some of the pseudospectra
of the "Grcar matrix" of dimension 400 [6], the exact spectrum, and converged
eigenvalue estimates (Ritz values) returned by a run of ARPACK (seeking
the eigenvalues of largest modulus) for this matrix. In the original article [23] that
described the algorithmic basis of ARPACK, Sorensen presented some similar plots
of di#culties encountered with the Grcar matrix. This is an extreme example where
the nonnormality 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 realize that anything was
# Received by the editors June 6, 2000; accepted for publication (in revised form) January 15,
2001; published electronically July 10, 2001.
http://www.siam.org/journals/sisc/23-2/37322.html
Oxford University Computing Laboratory, Parks Road, Oxford OX1 3QD, UK (TGW@comlab.
ox.ac.uk, LNT@comlab.ox.ac.uk).
In Matlab version 5, eigs was an M-file adapted from the Fortran ARPACK codes. Starting
with Matlab version 6, the eigs command calls the Fortran ARPACK routines themselves.
G. WRIGHT AND LLOYD N. TREFETHEN
Fig. 1. The #-pseudospectra for of the Grcar matrix (dimension
400), with the actual eigenvalues shown as solid stars and the converged eigenvalue estimates (for
the eigenvalues of largest modulus) returned by ARPACK shown as open circles. The ARPACK
estimates lie between the 10 -16 and 10 -17 pseudospectral contours.
amiss. Once the pseudospectra are plotted too, it is obvious.
Computing the pseudospectra of a matrix of dimension N is traditionally an expensive
task, requiring an O(N 3 ) singular value decomposition at each point in a grid.
For a reasonably fine mesh, this leads to an O(N 3 ) algorithm with the constant of the
order of thousands. Recent developments in algorithms for computing pseudospectra
have improved the constant [28], and the asymptotic complexity for large sparse matrices
[3, 12], but these are still fairly costly techniques. In this paper we show that
for large matrices, we can cheaply compute an approximation to the pseudospectra
in a region near the interesting eigenvalues. Our method uses the upper Hessenberg
matrix constructed after successive iterations of the implicitly restarted Arnoldi al-
gorithm, as implemented in ARPACK. Among other things, this means that after
performing an eigenvalue computation with ARPACK or 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."
Some ideas related to ours have appeared in earlier papers by Nachtigal, Reichel,
and Trefethen [17], Ruhe [19], Sorensen [23], Toh [24], and Toh and Trefethen [25].
For example, Sorensen plotted level curves of filter polynomials and observed that
they sometimes approximated pseudospectra, and Ruhe showed that pseudospectra
could be approximated by a rational Krylov method. What is new here is the explicit
development of a method for approximating pseudospectra based on ARPACK. Of
course, one could also consider the use of di#erent low-dimensional compressions of
a matrix problem such as those constructed by the Jacobi-Davidson algorithm [22].
Preliminary experiments, not reported here, show that this kind of Jacobi-Davidson
approximation of pseudospectra can also be e#ective.
We start by giving an overview of pseudospectra calculations and the implicitly
restarted Arnoldi iteration, followed by the practical details of our implementation
along with a discussion of some of the problems we have had to deal with. After this
we give some examples of the technique in practice. We also mention our Matlab
graphical user interface (GUI), which automates the computation of pseudospectra
LARGE-SCALE COMPUTATION OF PSEUDOSPECTRA 593
after the eigenvalues of a matrix have been computed by eigs in Matlab.
The computations presented in this paper were all performed using eigs in Matlab
version 6 (which essentially is ARPACK) and our GUI rather than the Fortran
ARPACK, although our initial experiments were done with the Fortran code.
2. Pseudospectra. There are several equivalent ways of defining # (A), the
#-pseudospectrum of a matrix A. The two most important (see, e.g., [27]) are perhaps
and
E) for some E with #E#}.
When the norms are taken to be the 2-norm, the definitions are equivalent to
where # min (-) denotes minimum singular value. This provides the basis of many
algorithms for computing pseudospectra. The most familiar technique is to use a grid
over the region of the complex plane of interest and calculate the minimum singular
value of zI - A at each grid point z. These values can then be passed to a contour
plotter to draw the level curves. For the rest of this paper, we consider the 2-norm;
other norms are discussed in [8, 28].
The reason for the cost of computation of pseudospectra is now clear: the amount
of work needed to compute the minimum singular value of a general matrix of dimension
N is O(N 3 ) (see, e.g., [5]). However, several techniques have been developed to
reduce this cost [28]. Here are two important ones.
(I) Project the matrix onto a lower dimensional invariant subspace via, e.g., a
partial Schur factorization (Reddy, Schmid, and Henningson [18]). This works well
when the interesting portion of the complex plane excludes a large fraction of the
eigenvalues of the matrix. In this case, the e#ect of the omitted eigenvalues on the
interesting portion of the pseudospectra is typically small, especially if the undesired
eigenvalues are well conditioned. Projection can significantly reduce the size of the
matrix whose pseudospectra we need to compute, making the singular value computation
dramatically faster. In general, the additional cost of projecting the matrix is
much less than the cost of repeatedly computing the smallest singular value for the
shifted original matrix.
Perform a single matrix reduction to Hessenberg or triangular form before
doing any singular value decompositions (Lui [12]), allowing the singular value calculations
to be done using a more e#cient algorithm.
One way of combining these ideas is to do a complete Schur decomposition of the
and then to reorder the diagonal entries of the triangular matrix
to leave the "wanted" eigenvalues at the top. The reordered factorization can then
be truncated leaving the required partial Schur factorization. We can now find the
singular values of the matrices shifted for each grid point z using either the original
matrix A or the triangular matrix
This allows us to work solely with the triangular matrix T once the O(N 3 ) factorization
has been completed. The minimum singular value of zI - T can be determined
594 THOMAS G. WRIGHT AND LLOYD N. TREFETHEN
in O(N 2 ) operations 2 using the fact that # min (zI - T
This can be calculated using either inverse iteration or inverse Lanczos iteration, which
require solutions to systems of equations with the matrix (zI - T
can be solved in two stages, each using triangular system solves.
By combining these techniques with more subtle refinements we have an algorithm
which is much more e#cient than the straightforward method. It is suggested in [28]
that the speedup obtained is typically a factor of about N/4, assuming the cost of
the Schur decomposition is negligible compared with that of the rest of the algorithm.
This will be the case on a fine grid for a relatively small matrix (N of the order
of a thousand or less), but for larger matrices the Schur decomposition is relatively
expensive, and it destroys any sparsity structure.
3. Arnoldi iteration. The Arnoldi iteration for a matrix A of dimension N
works by projecting A onto successive Krylov subspaces K j of dimension
starting vector v 1 [1, 20]. It builds an orthonormal basis for the Krylov
subspace by the Arnoldi factorization
is an upper Hessenberg matrix of dimension j, the columns of V j form an
orthonormal basis for the Krylov subspace, and f j is orthogonal to the columns of V j .
The residual term f j e # j can be incorporated into the first term V j H j by augmenting
the matrix H j with an extra row, all zeros except for the last entry, which is #f j #,
and including the next basis vector, v #, in the matrix V j . The matrix
H j is now rectangular, of size (j
The matrix
being rectangular, does not have any eigenvalues, but we can
define its pseudospectra in terms of singular values by (2.3). (With
the definition occasionally used that # is an eigenvalue of a rectangular matrix A if
A - #I is rank-deficient, where I is the rectangular matrix of appropriate dimension
with 1 on the diagonal and 0 elsewhere. In general, a rectangular matrix will have
no eigenvalues, but it will have nonempty #-pseudospectra for large enough #.) It can
then be shown [14, 25] that the pseudospectra of successive
are nested.
Theorem 3.1. Let A be an N - N matrix which is unitarily similar to a Hessenberg
matrix H, and let
denote the upper left (j section (in particular
could be created using a restarted Arnoldi iteration). Then for any # 0,
Thus as the iteration progresses, the pseudospectra of the rectangular Hessenberg
matrices better approximate those of A, which gives some justification for the
approximation
Unfortunately, this is only the case for the rectangular
matrices
There do not appear to be any satisfactory theorems to justify
a similar approximation for the square matrices H j , and of course for # su#ciently
small the #-pseudospectra of H j must be disjoint from those of A, since they will be
small sets surrounding the eigenvalues of H j , which are in general distinct from those
of A. This is not the case for the rectangular matrix as there will not be points in
the complex plane with infinite resolvent norm unless a Ritz value exactly matches an
can also be used for Hessenberg matrices [12], but for those we do not
have the advantage of projection.
LARGE-SCALE COMPUTATION OF PSEUDOSPECTRA 595
eigenvalue of the original matrix. That is, # (
typically empty for su#ciently
small #.
Although the property (3.2) is encouraging, theorems guaranteeing rapid convergence
in all cases cannot be expected. The quality of the approximate pseudospectra
depends on the information in the Krylov subspace, which in turn depends on the
starting vector v 1 . Any guarantee of rapid convergence could at best be probabilistic.
3.1. Implicitly restarted Arnoldi. In its basic form, the Arnoldi process may
require close to N iterations before the subspace contains good information about the
eigenvalues of interest. However, the information contained within the Hessenberg
matrix is very dependent on the starting vector relatively small
components of the eigenvectors corresponding to the eigenvalues which are not re-
quired, convergence may be quicker and the subspace size need not grow large. To
avoid the size of the subspace growing too large, practical implementations of the
Arnoldi iteration restart when the subspace size j reaches a certain threshold [20]. A
new starting vector -
v 1 is chosen which has smaller components in the directions of
eigenvectors corresponding to unwanted eigenvalues, and the process is begun again.
Implicit restarting [11, 23] is based upon the same idea, except that subspace is
only implicitly compressed to a single starting vector - v 1 . What is explicitly formed
is an Arnoldi factorization of size k based on this new starting vector, where k is the
number of desired eigenvalues, and this Arnoldi factorization is obtained by carrying
out implicitly shifted steps of the QR algorithm, with shifts possibly corresponding
to unwanted eigenvalue estimates. The computation now proceeds in an
accordion-like manner, expanding the subspace to its maximum size p, then compressing
to a smaller subspace. 3 This is computationally more e#cient than simple
restarting because the subspace is already of size k when the iteration restarts, and
in addition, the process is numerically stable due to the use of orthogonal transformations
in performing the restarting. This technique has made the Arnoldi iteration
competitive for finding exterior eigenvalues of a wide range of nonsymmetric matrices.
3.2. Arnoldi for pseudospectra. In a 1996 paper, Toh and Trefethen [25]
demonstrated that the Hessenberg matrix created during the Arnoldi process can
sometimes provide a good approximation to the pseudospectra of the original matrix.
They provided results for both the square matrix H j and the rectangular matrix
We choose to build our method around the rectangular Hessenberg matrices
though this makes the pseudospectral computation harder than if we worked with the
square matrix. The advantage of this is that we retain the properties of Theorem 3.1,
and the following in particular:
For every # 0, the approximate #-pseudospectrum generated by our ARPACK
algorithm is a subset of the #-pseudospectrum of the original matrix,
This is completely di#erent from the familiar situation with Ritz values, which
are, after all, the points in the 0-pseudospectrum of a square Hessenberg matrix. Ritz
values need not be contained in the true spectrum. Simply by adding one more row
to consider a rectangular matrix, we have obtained a guaranteed inclusion for every #.
The results presented by Toh and Trefethen focus on trying to approximate the
full pseudospectra of the matrix (i.e., around the entire spectrum) and they do not use
3 Our use of the variable p follows eigs. In ARPACK and in Sorensen's original paper [23], this
would be p
596 THOMAS G. WRIGHT AND LLOYD N. TREFETHEN
any kind of restarting in their implementation of the Arnoldi iteration. While this is
a useful theoretical idea, we think it is of limited practical value for computing highly
accurate pseudospectra since good approximations are often obtained generally only
for large subspace dimensions.
Our work is more local; we want a good approximation to the pseudospectra
in the region around the eigenvalues requested from ARPACK or eigs. By taking
advantage of ARPACK's implicit restarting, we keep the size of the subspace (and
hence
reasonably small, allowing us to compute (local) approximations to the
pseudospectra more rapidly, extending the idea of [25] to a fully practical technique
(for a more restricted problem).
4. Implementation. In deciding to use the rectangular Hessenberg matrix
we have made the post-ARPACK phase of our algorithm more di#cult. While the
simple algorithm of computing the minimum singular value of zI -A at each point has
approximately the same cost for a rectangular matrix as a square one, the speedup
techniques described in section 2 are di#cult to translate into the rectangular case.
The first idea, projection to a lower dimensional invariant subspace, does not
make sense for rectangular matrices because there is no such thing as an invariant
subspace. The second idea, preliminary triangularization using a Schur decompo-
sition, also does not extend to rectangular matrices, for although it is possible to
triangularize the rectangular matrix while keeping the same singular values (by performing
a QR factorization, for example), doing so destroys the vital property of
shift-invariance (see (2.4)).
However, our particular problem has a feature we have not yet considered: the
matrix is Hessenberg. One way to exploit this property is to perform a QR factorization
of the matrix obtained after shifting for each grid point. The upper triangular
matrix R has the same singular values as the shifted matrix, and they are also unchanged
on removing the last row of zeros, which makes the matrix square. We can
now use the inverse Lanczos iteration as in section 2 to find its smallest singular value.
The QR factorization can be done with an O(N 2 ) algorithm (see, e.g., [5, p. 228]),
which makes the overall cost O(N 2 ). Unfortunately, the additional cost of the QR
factorization at each stage makes this algorithm slightly slower for the small matrices
(dimensions 50-150) output from ARPACK than for square matrices of the same size,
but this appears to be the price to be paid for using matrices which have the property
of (3.3).
4.1. Refinements. In some cases we have found that inverse iteration to find
the minimum eigenvalue of (zI - R) # (zI - R) is more e#cient than inverse Lanczos
iteration but only when used with continuation (Lui [12]). Continuation works by
using the vector corresponding to the smallest singular value from the previous grid
point as the starting guess for the next grid point.
This sounds like a good idea; if the two shifted matrices di#er by only a small
shift, their singular values (and singular vectors) will be similar. When it works, it
generally means that only a single iteration is needed to satisfy the convergence crite-
rion. However, as Lui indicates, there is a problem with this approach if the smallest
and second smallest singular values "change places" between two values of z: the iteration
may converge to the second smallest singular value instead of the smallest, since
the starting vector had such a strong component in the direction of the corresponding
singular vector. This leads to the convergence criterion being satisfied for the wrong
singular value (even after several iterations).
LARGE-SCALE COMPUTATION OF PSEUDOSPECTRA 597
Choose max subspace size p-larger p for better pseudospectra.
Choose number of eigenvalues k-larger k for better pseudospectra.
Run ARPACK(A, p, k) to obtain
Define a grid over a region of C enclosing converged Ritz values.
For each grid point z:
Perform reduced QR factorization of shifted matrix: z -
I - #
Get # max (z) from Lanczos iteration on (R # R) -1 , random starting vector.
end.
Start GUI and create contour plot of the # min values.
Allow adjustment of parameters (e.g., grid size, contour levels) in GUI.
Fig. 2. Pseudocode for our algorithm.
In the course of our research, we have found several test matrices which su#er
from this problem, including the so-called Tolosa matrix [4]. Accordingly, because of
our desire to create a robust algorithm, we do not use inverse iteration. In theory it is
also possible to use continuation with inverse Lanczos iteration, but our experiments
indicate that the benefit is small and it again brings a risk of misconvergence.
Our algorithm (the main loop of which is similar to that in [28]) is summarized
in

Figure

2.
5. Practical examples. While one aim of our method is to approximate the
pseudospectra of the original matrix accurately, this is perhaps no more important
than the more basic mission of exhibiting the degree of nonnormality the matrix has,
so that the ARPACK or eigs user gets some idea of whether the Ritz values returned
are likely to be physically meaningful. Even in cases where the approximations of
the sets # (A) are inaccurate, a great deal may still be learned from their qualitative
properties.
In the following examples, ARPACK was asked to look for the eigenvalues of
largest real part except where otherwise indicated. However, the choice of region of
the complex plane to focus on is unimportant for our results and is determined by
which eigenvalues are of interest for the particular problem at hand. The number
of requested eigenvalues k was chosen rather arbitrarily to be large enough so that
the approximate pseudospectra clearly indicate the true behavior in the region of
the complex plane shown, and the maximum subspace size p was chosen to ensure
convergence of ARPACK for the particular choice of k. Experiments show that the
computed pseudospectra are not very sensitive to the choices of k and p, provided
they are large enough, but we have not attempted to optimize these choices.
5.1. Two extremes. Our first example (Figure 3), from Matrix Market [15],
shows a case where the approximation is extremely good. The matrix is the Jacobi
matrix of dimension 800 for the reaction-di#usion Brusselator model from chemical
engineering [7], and one seeks the rightmost eigenvalues. The matrix is not highly non-
normal, and the pseudospectra given by the approximation almost exactly match the
G. WRIGHT AND LLOYD N. TREFETHEN
Pseudospectra Approximate pseudospectra
Fig. 3. Pseudospectra for the matrix rdb800l (left) computed using the standard method, and
pseudospectra of the upper Hessenberg matrix of dimension computed using ARPACK
(right) in about 9% of the computer time (on the fine grid used here). Levels are shown for
. The number of matrix-vector products needed by ARPACK (nv ) is
1,493.
true pseudospectra around the converged Ritz values. This is a case where the pseudospectra
computed after running ARPACK indicate that the eigenvalues returned
are both accurate and physically meaningful, and that no further investigation is nec-
essary. In this computation we used a maximum subspace dimension of
requested eigenvalues.
The second case we consider is one where the matrix has a high degree of non-
normality-the Grcar matrix. As seen in Figure 1, ARPACK can converge to Ritz
values which are eigenvalues of a perturbation of order machine precision of the original
matrix, and the nonnormality of this particular matrix (here of dimension 400)
means that the Ritz values found can lie a long way from the spectrum of the matrix.

Figure

4 shows that the pseudospectra of the Hessenberg matrix (computed using
asking for eigenvalues of largest modulus) in this case are not
good approximations to the pseudospectra of the original one.
This is typical for highly nonnormal matrices-the Hessenberg matrix cannot capture
the full extent of the nonnormality, particularly when more than p eigenvalues
of the original matrix lie within the region of the complex plane in which the pseudospectra
are computed. In other words, the approximation is typically not so good in
areas away from the Ritz values computed, and then only accurately approximates the
pseudospectra of the original matrix when the Ritz values are good approximations
to the eigenvalues. Despite this, a plot like that of Figure 4 will instantly indicate to
the ARPACK user that the matrix at hand is strongly nonnormal and needs further
investigation.
5.2. A moderately nonnormal example. While the above examples show
two extreme cases, many important applications are more middle-of-the-range, where
LARGE-SCALE COMPUTATION OF PSEUDOSPECTRA 599
Pseudospectra ARPACK approximate pseudospectra
Fig. 4. The pseudospectra of the Grcar matrix of dimension 400 (left) computed using the
standard method, and the pseudospectra of the upper Hessenberg matrix of dimension 50 computed
using ARPACK (right) in about 8% of the computer time (on this fine grid). Contours are shown
-0.4
-0.3
-0.2
-0.4
-0.3
-0.2
-0.10.10.3Pseudospectra Approximate pseudospectra
Fig. 5. Pseudospectra for linearized fluid flow through a circular
pipe at Reynolds number 10,000 (streamwise-independent disturbances with azimuthal wave number
machine precision is su#cient to accurately converge the eigenvalues, but pronounced
nonnormality may nevertheless diminish the physical significance of some of them. A
good example of a case in which this is important is the matrix created by linearization
about the laminar solution of the Navier-Stokes equations for fluid flow in an infinite
circular pipe [26]. (Our matrix is obtained by a Petrov-Galerkin spectral discretization
of the Navier-Stokes problem due to Meseguer and Trefethen [16]. The axial
and azimuthal wave numbers are 0 and 1, respectively, and the matrix dimension is
402.) The pseudospectra are shown in Figure 5, and although the eigenvalues all have
negative real part, implying stability of the flow, the pseudospectra protrude far into
the right half-plane. This implies pronounced transient growth of some perturbations
of the velocity field in the pipe, which in the presence of nonlinearities in practice may
lead to transition to turbulence [30]. The approximate pseudospectra also highlight
G. WRIGHT AND LLOYD N. TREFETHEN
-226
-0.50.5Brusellator Wave Model
Crystal
Fig. 6. Left: The pseudospectra for the Brusselator wave model,
nv= 16,906. Right: Pseudospectral contours for a matrix of dimension
10,000 from the Crystal set at Matrix Market,
this behavior. The parameters used here were
5.3. Larger examples. We now consider four larger examples. The first is the
Brusselator wave model from Matrix Market (not to be confused with the very first
example), which models the concentration waves for reaction and transport interaction
of chemical solutions in a tubular reactor [9]. Stable periodic solutions exist for a
parameter when the rightmost eigenvalues of the Jacobian are purely imaginary. For a
matrix of dimension 5,000, using a subspace of dimension 100 and asking ARPACK for
20 eigenvalues, we obtained the eigenvalue estimates and approximate pseudospectra
shown in Figure 6 (left). The departure from normality is evidently mild, and the
conclusion from this computation is that the Ritz values returned by ARPACK are
likely to be accurate and the corresponding eigenvalues physically meaningful.

Figure

6 (right) shows approximate pseudospectra for a matrix of dimension
10,000, taken from the Crystal set at Matrix Market, which arises in a stability analysis
of a crystal growth problem [32]. The eigenvalues of interest are the ones with
largest real part. The fact that we can see the 10 -13 pseudospectrum (when the axis
scale is O(1)) indicates that this matrix is significantly nonnormal, and although the
matrix is too large for us to be able to compute its exact pseudospectra for compar-
ison, this is certainly a case where the nonnormality could be important, making all
but the rightmost few eigenvalues of dubious physical significance in an application.
The ARPACK parameters we used in this case were and the computation
took about one hour on our Sun Ultra 5 workstation. Although we do not
have the true pseudospectra in this case, we would expect that the rightmost portion
should be fairly accurate where there is a good deal of Ritz data and relatively little
nonnormality. We expect that the leftmost portion is less accurate where the e#ect
of the remaining eigenvalues of the matrix unknown to the approximation begins to
become important.
The third example, Figure 7 (left), shows the Airfoil matrix created by performing
transient stability analysis of a Navier-Stokes solver [13], also from Matrix Market. In
this case the matrix appears fairly close to normal, and the picture gives every reason
to believe that the eigenvalues have physical meaning. Using
ARPACK took about 9 hours to converge to the eigenvalues, while we were able to
plot the pseudospectra in about 3 minutes (even on the fine grid used here).
Our final example is a matrix which is bidiagonal plus random sparse entries
LARGE-SCALE COMPUTATION OF PSEUDOSPECTRA 601
-0.4
Airfoil
Random
Fig. 7. Left: The pseudospectra for the Airfoil matrix from
Matrix Market of dimension 23,560, nv= 72,853. Right: The
for our random matrix, with nv= 61,347.
elsewhere, created in Matlab by
The approximation to the pseudospectra of the matrix of dimension 200,000 is shown
in

Figure

7 (right), from which we can conclude that the eigenvalue estimates returned
are probably accurate, but that the eigenvalues toward the left of the plot would likely
be of limited physical significance in a physical application, if there were one, governed
by this matrix. We used a subspace size of 50 and requested 30 eigenvalues from this
example, and the whole computation took about 26 hours.
6. MATLAB GUI. We have created a Matlab GUI to automate the process
of computing pseudospectra, and Figure 8 shows a snapshot after a run of ARPACK.
Initially the pseudospectra are computed on a coarse grid to give a fast indication
of the nonnormality of the matrix, but the GUI allows control over the number of
grid points if a higher quality picture is desired. Other features include the abilities
to change the contour levels shown without recomputing the underlying values, and
to select parts of the complex plane to zoom in for greater detail. The GUI can
also be used as a graphical front end to our other pseudospectra codes for computing
pseudospectra of smaller general matrices. The codes are available on the World Wide
Web from http://web.comlab.ox.ac.uk/oucl/work/nick.trefethen/.
7. Discussion. The examples of section 5 give an indication of the sort of pictures
obtained from our technique. For matrices fairly close to normal, the approximation
is typically a very close match to the exact pseudospectra, but for more
highly nonnormal examples the agreement is not so close. This is mainly due to the
e#ect of eigenvalues which the Arnoldi iteration has not found: their e#ect on the
pseudospectra is typically more pronounced for nonnormal matrices.
The other point to note is that if we use the Arnoldi iteration to look (for ex-
G. WRIGHT AND LLOYD N. TREFETHEN
Fig. 8. A snapshot of the Matlab GUI after computing the pseudospectra of a matrix.
ample) for eigenvalues of largest real part, the rightmost part of the approximate
pseudospectra will be a reasonably good approximation. This can clearly be seen in

Figure

4 where we are looking for eigenvalues of largest modulus: the top parts of the
pseudospectra are fairly good and only deteriorate lower down where the e#ect of the
"unfound" eigenvalues becomes important.
However, as mentioned in the introduction, creating accurate approximations of
pseudospectra was only part of the motivation for this work. Equally important has
been the goal of providing information which can help the user of ARPACK or eigs
decide whether the computed eigenvalues are physically meaningful. For this purpose,
estimating the degree of nonnormality of the matrix is more important than getting
an accurate plot of the exact pseudospectra.
One of the biggest advantages of our technique is that while the time spent on the
computation grows as the dimension of the matrix increases, the time spent
on the pseudospectra computation remains roughly constant. This is because the
pseudospectra computation is based just on the final Hessenberg matrix, of dimension
typically in the low hundreds at most. Figure 9 shows the proportion of time spent
on the pseudospectra part of the computation for the examples we have presented
here. These timings are based on the time to compute the initial output from our
LARGE-SCALE COMPUTATION OF PSEUDOSPECTRA 603
Airfoil
Crystal
Grcar
Pipe Flow
R.-D. Brusselator
Random
Brusellator Wave Model
Dimension, N
Pseudospectra
time
Total
time
Fig. 9. The proportion of the total computation time spent on computing the pseudospectra
for the examples presented in this paper. For large N, the pseudospectra are obtained at very little
additional cost.
R.-D. Brusselator
Grcar
Fig. 10. Lower-quality plots of pseudospectra produced by our GUI on its default coarse
Such plots take just a few seconds. The matrices shown are rdb800l with

Figure

and the Grcar matrix with

Figure

4).
GUI using a coarse grid such as those illustrated in Figure 10, which is our standard
resolution for day-to-day work. For the "publication quality" pseudospectra of the
resolution of the other plots in this paper, the cost is about thirty times higher, but
G. WRIGHT AND LLOYD N. TREFETHEN
this is still much less than the cost of ARPACK for dimensions N in the thousands.

Acknowledgments

. We would like to thank Rich Lehoucq for his advice on
beginning during his visit to Oxford in October-November 1999, Penny
Anderson of The MathWorks, Inc., for her help with the beta version of the new eigs
command, and Mark Embree for his comments on drafts of the paper.



--R

The principle of minimized iteration in the solution of the matrix eigenvalue problem
Vorst, eds., Templates for the Solution of Algebraic Eigenvalue Problems: A Practical Guide
Computing the field of values and pseudospectra using the Lanczos method with continuation
Stability analysis in aeronautical industries
Matrix Computations
Operator Coe
Theory and Applications of Hopf Bifurcation
A block algorithm for matrix 1-norm estimation
Waves in Distributed Chemical Systems: Experiments and Computations
Package, http://www.

Computation of pseudospectra by continuation
Eigenvalue calculation procedure for an Euler- Navier-Stokes solver with applications to flows over airfoils
On hybrid iterative methods for nonsymmetric systems of linear equations
Matrix Market
A Spectral Petrov-Galerkin Formulation for Pipe Flow I: Linear Stability and Transient Growth
A hybrid GMRES algorithm for non-symmetric linear systems
Pseudospectra of the Orr-Sommerfeld operator
Rational Krylov algorithms for nonsymmetric eigenvalue problems.
Variations on Arnoldi's method for computing eigenelements of large unsymmetric matrices
Numerical Methods for Large Eigenvalue Problems
A Jacobi-Davidson iteration method for linear eigenvalue problems
Implicit application of polynomial filters in a k-step Arnoldi method
Matrix Approximation Problems and Nonsymmetric Iterative Methods
Calculation of pseudospectra by the Arnoldi iteration
Spectra and pseudospectra for pipe Poiseuille flow
Pseudospectra of matrices
Computation of Pseudospectra

Hydrodynamic stability without eigenvalues
Computational Methods for Large Eigenvalue Problems
Numerical Computation of the Linear Stability of the Di
--TR

--CTR
S.-H. Lui, A pseudospectral mapping theorem, Mathematics of Computation, v.72 n.244, p.1841-1854, October
Lorenzo Valdettaro , Michel Rieutord , Thierry Braconnier , Valrie Frayss, Convergence and round-off errors in a two-dimensional eigenvalue problem using spectral methods and Arnoldi-Chebyshev algorithm, Journal of Computational and Applied Mathematics, v.205 n.1, p.382-393, August, 2007
Kirk Green , Thomas Wagenknecht, Pseudospectra and delay differential equations, Journal of Computational and Applied Mathematics, v.196 n.2, p.567-578, 15 November 2006
C. Bekas , E. Kokiopoulou , E. Gallopoulos, The design of a distributed MATLAB-based environment for computing pseudospectra, Future Generation Computer Systems, v.21 n.6, p.930-941, June 2005
Computing smallest singular triplets with implicitly restarted Lanczos bidiagonalization, Applied Numerical Mathematics, v.49 n.1, p.39-61, April 2004
implicit restarting;pseudospectra;ARPACK;arnoldi;eigenvalues