A Block Algorithm for Matrix 1-Norm Estimation, with an Application to 1-Norm Pseudospectra.
The matrix 1-norm estimation algorithm used in LAPACK and various other software libraries and packages has proved to be a valuable tool. However, it has the limitations that it offers the user no control over the accuracy and reliability of the estimate and that it is based on level 2 BLAS operations. A block generalization of the 1-norm power method underlying the estimator is derived here and developed into a practical algorithm applicable to both real and complex matrices. The algorithm works with n  t matrices, where t is a parameter. For t=1 the original algorithm is recovered, but with two improvements (one for real matrices and one for complex matrices). The accuracy and reliability of the estimates generally increase with t and the computational kernels are level 3 BLAS operations for t > 1. The last t-1 columns of the starting matrix are randomly chosen, giving the algorithm a statistical flavor. As a by-product of our investigations we identify a matrix for which the 1-norm power method takes the maximum number of iterations. As an application of the new estimator we show how it can be used to efficiently approximate 1-norm pseudospectra.
Introduction
. Research in matrix condition number estimation began in the
1970s with the problem of cheaply estimating the condition number
and an approximate null vector of a square matrix A given some factorization of
it. The earliest algorithm is one of Gragg and Stewart [10]. It was improved by
Cline, Moler, Stewart and Wilkinson [4], leading to the 1-norm condition estimation
algorithm used in LINPACK [8] and later included in Matlab (function rcond).
During the 1980s, attention was drawn to various componentwise condition numbers
and it was recognized that most condition estimation problems can be reduced
to the estimation of kAk when matrix-vector products Ax and A T x can be cheaply
computed [2], [16, Sec. 14.1]. Hager [12] derived an algorithm for the 1-norm that is
a special case of the more general p-norm power method proposed by Boyd [3] and
later investigated by Tao [18]. Hager's algorithm was modified by Higham [14] and
incorporated in LAPACK (routine xLACON) [1] and Matlab (function condest).
The LINPACK and LAPACK estimators both produce estimates that in practice
are almost always within a factor 10 and 3, respectively, of the quantities they are
estimating [13], [14], [15]. This has been entirely adequate for applications where only
an order of magnitude estimate is required, such as the evaluation of error bounds.
However, in some applications an estimate with one or more correct digits is required
(see, for example, the pseudospectra application described in section 4), and for these
the LINPACK and LAPACK estimators have the drawback that they offer the user
no way to control or improve the accuracy of the estimate. Here, "accuracy" refers
to average case behaviour. Also of interest for a norm estimator is its worst case
behaviour, that is, its "reliability".
This work was supported by Engineering and Physical Sciences Research Council grants
GR/L76532 and GR/L94314.
y Department of Mathematics, University of Manchester, Manchester, M13 9PL, England
(higham@ma.man.ac.uk, http://www.ma.man.ac.uk/~higham/).
z Department of Mathematics, University of Manchester, Manchester, M13 9PL, England
(ftisseur@ma.man.ac.uk, http://www.ma.man.ac.uk/~ftisseur/).
N. J. HIGHAM AND F. TISSEUR

Table
Empirical probabilities that minf e OE s one A of the form
inv(randn(100)) and N(0; 1) vectors x j .
Since estimating kAk cheaply appears inevitably to admit the possibility of arbitrarily
poor estimates (although proving so is an open problem [6]), one might look
for an approach for which probabilistic statements can be made about the accuracy of
the estimate. The definition of a subordinate matrix norm
suggests the estimate
ae
oe
where s is a parameter and the x j are independently chosen random vectors. In
the case of the 2-norm (kxk and an appropriate distribution of the x j ,
explicit bounds are available on the probability of such estimates being within a given
factor of kAk [7]. As is done in [11] for certain estimates of the Frobenius norm we can
scale our estimates e OE s / ' s OE s , where the constant ' s is chosen so that the expected
value of e
OE s is kAk (note that e
OE s can therefore be greater or less than kAk). For
the 1-norm (kxk
which is our interest here, we investigate this approach
empirically. For a fixed matrix A of the form, in Matlab notation, inv(randn(100)),

Table

1.1 shows the observed probabilities that minf e
various ff and s, based on 1000 separate evaluations of the OE s with vectors x j from
the normal N(0; 1) distribution, and where ' s is determined empirically so that the
mean of the e
OE s is
The table shows that even with only 35% of the estimates
were within a factor 0.9 of the true norm. The statistical sampling technique is
clearly too crude to be useful for obtaining estimates with correct digits. One way to
exploit the information contained in the vectors Ax j is to regard them as first iterates
from the 1-norm power method with starting vectors x j and to continue to iterate.
These considerations motivate the block generalization of the 1-norm power method
that we present in this paper. Our block power method works with a matrix with t
columns instead of a vector. To give a feel for how our new estimator compares with
the sampling technique, we applied the estimator (Algorithm 2.4) to 1000 random
matrices of the form inv(randn(100)). The results are shown in Table 1.2; the
estimates for are obtained at approximately the same cost as the estimates
e
OE s for respectively. The superiority of the new estimator is clear.
In section 2 we derive the block 1-norm power method and develop it into a
practical algorithm for both real and complex matrices. In section 3 we present
numerical experiments that give insight into the behaviour of the algorithm. An
application involving complex matrices is given in section 4, where we describe how
the algorithm can be used to approximate 1-norm pseudospectra. Conclusions are
presented in section 5.

Table
Empirical probabilities that est - for est from Algorithm 2.4 and A of the form
inv(randn(100)).
Finally, we note that although our work is specific to the 1-norm, the 1-norm
can be estimated by applying our algorithm to A   , since
2. Block 1-Norm Power Method. The 1-norm power method is a special case
of Boyd's p-norm power method [3] and was derived independently by Hager [12]. For
a real matrix A we denote by sign(A) the matrix with (i; according
as a ij - 0 or a ij ! 0. The jth column of the identity matrix is denoted by e j .
Algorithm 2.1 (1-norm power method). Given A 2 R n\Thetan this algorithm computes
and x such that
repeat
quit
(smallest such
Algorithm 2.1 was modified by Higham [14, Alg. 4.1] (see also [15], [16, Alg. 14.4])
to improve its reliability and efficiency. The modifications that improve the reliability
are, first, to force at least two iterations and, second, to take as the final estimate the
maximum of that produced by the algorithm and
The vector b is a heuristic choice intended to "pick out" any large elements of A in
those cases where such elements fail to be revealed during the course of the algorithm.
Efficiency is improved by terminating the algorithm after computing - if it is the same
as the previous -, since it can be shown that convergence would otherwise be declared
after the subsequent computation of z.
To obtain a more accurate and reliable estimate than that provided by Algorithm
2.1 we could run the algorithm t times in succession on t different starting
vectors. This idea was suggested in [12], with each starting vector being the mean
of the unit vectors e j not already visited and with the algorithm being prohibited
from visiting unit vectors previously visited. Note that the estimates obtained this
way are nondecreasing in t. This approach has two weaknesses: it allows limited
communication of information between the t different iterations and the highest level
computational kernel remains matrix-vector multiplication. We therefore develop a
block algorithm that works with an n \Theta t matrix as a whole instead of t separate
4 N. J. HIGHAM AND F. TISSEUR
n-vectors. The block approach offers the potential of better estimates, through providing
more information on which to base decisions, and it allows the use of level
3 BLAS operations, thus promising greater efficiency. The following algorithm estimates
not just the 1-norm of A, but, as a by-product, the 1-norms of the t columns
of A having largest 1-norms.
Algorithm 2.2 (block 1-norm power method). Given A 2 R n\Thetan and a positive
integer t, this algorithm computes vectors g and ind with
t such that g j is a lower bound for the 1-norm of the column of A of jth
largest 1-norm.
Choose starting matrix X 2 R n\Thetat with columns of unit 1-norm.
repeat
Sort g so that
ind best
re-order ind correspondingly.
Like the basic 1-norm power method Algorithm 2.2 has the attractive
property that it generates increasing sequences of estimates. Denote with a superscript
"(k)" quantities from the kth iteration of the loop in Algorithm 2.2 and let a j denote
the jth column of A.
Lemma 2.3. The sorted vectors g (k) and h (k) satisfy
and
2:
Proof. First, we have
1-j-t
r
But
r
y
r
r
Thus
r
r k1 kx (k)
r k1 - max
Furthermore, if h (k)
r s (k)
To prove (2.3), assume without loss of generality that x (k\Gamma1)
t. Then
has the form
Z
ff t+1 ff t+1
where each ff i in row i is, in general, different, and ff i - ka n. The algorithm
chooses the ind values corresponding to the t rows of Z (k) with largest 1-norm, and
the g (k+1)
on the next stage are at least as large as these 1-norms. Since row j
of Z (k) has 1-norm ka
follows.
Algorithm 2.2 has three possible sources of inefficiency. First, the columns of S are
vectors of \Sigma1s and so a pair of columns s i and s j may be parallel
case z and in computing Z a matrix-vector multiplication is redundant. There
can be up to t\Gamma1 redundant matrix-vector products per iteration, this maximum being
achieved on the second iteration with S the matrix of ones when A has nonnegative
elements. The second possible inefficiency is that a column of S may be parallel to one
from the previous iteration, in which case the formation of
redundant computation. We choose to detect parallel columns and replace them by
random vectors randf\Gamma1; 1g not already in S or the previous
denotes a vector with elements from the uniform distribution on the set f\Gamma1; 1g.
The detection is done by forming inner products between columns and looking for
elements of magnitude n. The total cost of these computations is O(nt 2 ) flops, which
is negligible compared with the 2n 2 t flops required for each matrix product in the
algorithm, since t - n in practice.
Our strategy could be extended to check for parallel columns between the current
S and all previous S; we return to this possibility in section 3. If all the columns
of S are parallel to columns from the previous iteration then it is easy to see that
Algorithm 2.2 is about to converge; we therefore immediately terminate the iteration
without computing Z.
Finally, on step k we can have x (k)
so that the computation of Ax (k)
then repeats an earlier computation. Repeated
vectors e j are easily avoided by keeping track of the indices of all the previously used
e j and selecting ind ind t from among the indices not previously used. If all
the indices are repeats, then again we prematurely terminate the iteration, saving a
matrix product.
Note that the first and third inefficiencies are possibly only for t ? 1, while the
second can happen for the original 1-norm power method.
Our strategy of detecting redundant computations and replacing them by ones
that provide new information has three benefits. First, it can reduce the amount of
computation, through premature detection of convergence. Second, it can lead to better
estimates. For the columns of S this depends on the random replacement vectors
generated, but for the e j the improvement is deterministic up to future replacements
6 N. J. HIGHAM AND F. TISSEUR
of columns of S. The third benefit is that the dimensions of the matrix multiplications
remain constant at each iteration, as opposed to varying if we simply skip redundant
computations; this helps us to make efficient use of the computing resources.
The next algorithm incorporates these modifications. The algorithm is forced to
take at least 2 (and at most itmax) iterations so that it computes at least t columns
of A; it also explicitly identifies the approximate maximizing vector that achieves the
norm estimate.
Algorithm 2.4 (practical block 1-norm estimator). Given A 2 R n\Thetan and positive
integers t and itmax - 2, this algorithm computes a scalar est and vectors v and
w such that est -
Choose starting matrix X 2 R n\Thetat with columns of unit 1-norm.
recording indices of used unit vectors e j .
est
est
est ? est old or
ind best
est - est old , est = est old , goto (6), end
est
(2) If every column of S is parallel to a column of S old , goto (6), end
Ensure that no column of S is parallel to another column of S
or to a column of S old by replacing columns of S by randf\Gamma1; 1g.
best , goto (6), end
re-order ind correspondingly.
(5) If ind(1: t) is contained in ind \Gamma hist, goto (6), end
Replace ind(1: t) by the first t indices in ind(1: n) that are
not in ind \Gamma hist.
hist ind(1: t)]
best
Note that this algorithm does not explicitly compute lower bounds for the 1-norms
of all t largest columns of A. If this information is required (and we are not aware of
any applications in which it is needed) then it can be obtained by keeping track of
the largest
Inequality (2.2), now expressed as est (k) - h (k)
est (k+1) , is still valid, except
that h (k)
est (k+1) is possible on the last iteration if the original ind (k)
1 is a repeat
(this event is handled by the test (1)). However, (2.3) is no longer true, because of
the avoidance of repeated indices.
How do Algorithms 2.4 and 2.2 compare? Ignoring the itmax test, Algorithm 2.4
terminates in at most n=t+1 iterations, since t vertices e j are visited on each iteration
after the first and no vertex can be visited more than once. The same cannot be said
of Algorithm 2.2 because of the possibility of repeated vertices. Algorithm 2.4 can
produce a smaller estimate than Algorithm 2.2: when a redundant computation is
avoided, the new information computed can lead to an apparently more promising
vertex (based on the relative sizes of the h i ) replacing one that actually corresponds
to a larger column of A. However it is more likely, when Algorithms 2.4 and Algorithm
2.2 produce different results, that Algorithm 2.4 produces a better estimate, as
in the following example.
Example 2.5. For a certain starting
obtained the following results. For Algorithm 2.2:
1: (0, 8.10e-001) (0, 6.13e-001)
2: (10, 1.77e+000) (4, 1.76e+000)
Underestimation ratio: 1.99e-001
The first column denotes the iteration number. The kth row gives the sorted g (k)
t, each one preceded by the corresponding index ind i . (Since X 1 does not have
columns e j , the ind i for the first iteration are shown as zero.) For Algorithm 2.4:
1: (0, 8.10e-001) (0, 6.13e-001)
2: (10, 1.77e+000) (4, 1.76e+000)
parallel column between S and S-old
3: (2, 8.87e+000) (5, 4.59e+000)
Exact estimate!
Algorithm 2.2 converges after 2 iterations and produces an estimate too small by a
factor 5. However, on the second iteration Algorithm 2.4 detects a column of S parallel
to one of S old and replaces it. The new column produces a different Z matrix and
causes the convergence test (4) to be failed. The extra iteration visits the (unique)
column of maximum 1-norm, so an exact estimate is obtained.
When 2.4 differs from the modified version of Algorithm 2.1
used in LAPACK 3.0 [14, Alg. 4.1], [16, Alg. 14.4] in two ways. First, Algorithm 2.4
does not use the "extra estimate" (2.1). Second, the LAPACK algorithm checks for
but not for - does Algorithm 2.4. This is an oversight.
We recommend that LAPACK's xLACON be modified to include the extra test. This
change will not affect the estimates produced but will sometimes reduce the number
of iterations.
We now explain our choice of starting matrix. We take the first column of X to
be the vector of 1s, which is the starting vector used in Algorithm 2.1. This has the
advantage that for a matrix with nonnegative elements the algorithm converges with
an exact estimate on the second iteration, and such matrices arise in applications,
for example as a stochastic matrix or as the inverse of an M-matrix. The remaining
columns are chosen as randf\Gamma1; 1g, with a check for and correction of parallel columns,
exactly as for S in the body of the algorithm. We choose random vectors because it
is difficult to argue for any particular fixed vectors and because randomness lessens
the importance of counterexamples (see the comments in the next section).
Next, we consider complex matrices, which arise in the pseudospectrum application
of section 4. Everything in this section remains valid for complex matrices
provided that sign(A) is redefined as the matrix (a
transposes are replaced by conjugate transposes. The matrix S is now complex with
elements of unit modulus and we are much less likely to find parallel columns of S
8 N. J. HIGHAM AND F. TISSEUR
from one iteration to the next or within the current S. Therefore for complex matrices
we omit the tests for parallel columns. However, we take the same, real, starting
matrix. There is one further question in the complex case. In the analogue of Algorithm
2.1 for complex matrices in [14, Alg. 5.1] z is defined as z = Re(A   -), based
on subgradient considerations. In our block algorithm should we take
The former can be justified from heuristic considerations and preserves
more information about A. We return to this question in the next section.
The motivation for Algorithm 2.4 is to enable more accurate and reliable estimates
to be obtained than are provided by the 1-norm power method. The question arises
of how the accuracy and reliability of the estimates varies with t. Little can be
said theoretically because, unlike the approach of [12] mentioned at the start of this
section, the estimates are not monotonic in t. If we run Algorithm 2.4 for t 1 and for
using a common set of t 1 starting vectors, we can obtain a smaller estimate
for t 2 than for t 1 , because a less promising choice of unit vector e j can turn out to
be better than a more promising choice made with more available information. Non-monotonicity
is unlikely, however, and we argue that it is a price worth paying for
the other advantages that accrue. In the next section we investigate the behaviour of
Algorithm 2.4 empirically.
3. Numerical Experiments. Our aim in this section is to answer the following
questions about Algorithm 2.4, bearing in mind that for the algorithm is an
implementation of the well understood 1-norm power method.
1. How does the accuracy and reliability of the norm estimates vary with t?
2. How good are the norm estimates in general?
3. How does the number of iterations behave for t ? 1?
Note that we are not searching for counterexamples, as was done for previous condition
estimators [4], [5], [14]. We know that for a fixed starting matrix and any t - n there
must be families of matrices whose norm is underestimated by an arbitrarily large
factor, since the algorithm samples the behaviour of the n \Theta n matrix A on fewer than
vectors. But since the algorithm uses a random starting matrix for t ? 1, each
counterexample will be valid only for particular starting matrices.
All our tests have been performed with Matlab.
Our first group of tests deals with random real matrices. Amongst the matrices
A we used were:
1. A from the normal N(0; 1) distribution (denoted randn) and its inverse, orthogonal
QR factor, upper triangular part, and inverse of the upper triangular
part.
2. A from the uniform distribution on the set f\Gamma1; 0; 1g (denoted rand(-1,0,1-)),
A \Gamma1 from the uniform distribution on the interval [0; 1].
3. A and A \Gamma1 of the form U \Sigma V T where U and V are random orthogonal matrices
and with the singular values oe i distributed exponentially,
arithmetically or with all except the smallest equal to unity, and with 2-norm
condition number ranging from 1 to 10 16 .
Note that we omit, for example, matrices from the uniform [0; 1] and uniform f\Gamma1; 1g
distributions, because for such matrices Algorithm 2.4 is easily seen to produce the
exact norm for all t.
We chose n and t in the range
For each test matrix we recorded a variety of statistics including the underestimation
ratio averaged and minimized over each type of A for fixed n and
t, the relative error jest and the number of iterations. We declared an
estimate exact if the relative error was no larger than 10 \Gamma14 (the unit roundoff is of
For a given matrix A we first generated a starting
columns, where t max is the largest value of t to be used, and then ran Algorithm 2.4
using starting t). In this way we could see the effect of
increasing t. In particular, we checked what percentage of the estimates for a given t
were at least as large as the estimates for all smaller t; we denote this "improve%".
We set Algorithm 2.4.
First, we give some general comments on the results.
1. Increasing t usually gave larger average and minimum underestimation ratios,
though there were exceptions. The quantity improve% was 100 about half
the time and never less than 78.
2. The number of iterations averaged between 2 and 3 throughout, with maxima
ranging from 2 to 5 depending on the type of matrix. Thus increasing t from
1 has little effect on the number of iterations-an important fact that could
not be predicted from the theory. Although there are specially constructed
examples for which the 1-norm power method requires many iterations (one
is described below), it is rare for the limit of 5 iterations to come into effect.
3. Throughout the tests we also computed the extra estimate (2.1) used by
the LAPACK norm estimator [14]. As expected, in none of our tests (with
random matrices) was the extra estimate larger than the estimate provided
by Algorithm 2.4 with
In

Tables

3.1 and 3.2 we show detailed results for two particular types of random
matrix from among those described above, with
rand(-1,0,1-). The columns headed "Products" show the average and maximum
total number of matrix products AX and A T S. In each case 5000 matrices were used.
For the matrices inv(randn) taking significantly improves the worst-case and
average estimates and the proportion of exact estimates over
estimate is exact almost 98 percent of the time. For the matrices rand(-1,0,1-)
the improvements as t increases are less dramatic but still useful; notice that exactly
four matrix products were required in every case. As well as recording the
number of products, we checked how convergence was achieved. For the matrices
rand(-1,0,1-) convergence was always achieved at the test (4) in Algorithm 2.4,
while for inv(randn) convergence was declared at tests (4) and (2) in approximately
96 and 4 percent of the cases, respectively (with just a few instances of convergence
at (5) for t - 2). The last two columns of Table 3.1 show the average and maximum
number of parallel columns of S that were detected. A small number of repeated e j
vectors were detected and replaced (the largest average was 0.03, occurring for
and the maximum number of 7 occurred for columns or repeated
were detected for the matrices in Table 3.2.
The strategy of replacing parallel columns of S and repeated e j vectors has little
effect on the overall performance of Algorithm 2.4 in our tests with random matrices.
Since particular examples can be found where it is beneficial (see Example 2.5) and
the cost is negligible, we feel its use is worthwhile. However, we see no advantage
to extending the strategy to compare the columns of S with those of all previous S
matrices.
Higham [15] gives a tridiagonal matrix for which the 1-norm power method (Algo-
rithm 2.1) requires n iterations to converge. We have constructed a matrix for which
N. J. HIGHAM AND F. TISSEUR

Table
Results for 5000 matrices inv(randn) of dimension 100.
Underest. ratio Products Parallel cols.
t min average % exact average max improve% average max
9 0.893 1.000 99.88 4.0 4 100.00 1.22 15

Table
Results for 5000 matrices rand(f-1,0,1g) of dimension 100.
Underest. ratio Products
t min average % exact average max improve%
9 0.775 0.951 28.56 4 4 90.74
the maximum iterations is required. It is the inverse of a bidiagonal matrix:
An
1 .
. ff3
1 .
. \Gammaff3
(The minus sign in front of the matrix is necessary!) It is straightforward to show
that when Algorithm 2.1 is applied to An (ff) it produces, for

Table
Results for matrix A 100 repetitions.
Underest. ratio Products
t min average % exact average max improve%
6 1.000 1.000 100.00 4.6 11   100.00
7 1.000 1.000 100.00 4.3 11   100.00
8 1.000 1.000 100.00 4.2 8 100.00
9 1.000 1.000 100.00 4.1 8 100.00
Thus every column of An (ff) is computed, in order from first to last, and the exact
norm is obtained. If the algorithm is terminated after p iterations then it produces
the estimate (1 \Gamma ff
behaves in exactly the same way.
We applied Algorithm 2.4 1000 times to A 100 as in all
our tests). The results are shown in Table 3.3; in the 6th column "11   " denotes that
convergence was declared because the iteration limit was reached (the percentage of
such occurrences varied from 100% for to 0.1% for 7). The underestimation
ratio for agrees with the theory and is unacceptably small (note that there is
no randomness, and hence only one estimate, for 1). But for all t - 2 the average
norm estimates are satisfactory. The extra estimate (2.1) has the value 0.561; thus it
significantly improves the estimate for but is worse than the average estimates
for all greater t.
For complex matrices we have tried both at (3) in
Algorithm 2.4. Tables 3.4 and 3.5 compare the two choices for 5000 100 \Theta 100 random
complex matrices of the form inv(rand+i*rand), where rand is a matrix from the
uniform distribution on the interval [0; 1]. In this test, larger underestimation ratios
are obtained for the percentage of exact estimates is higher, and the
statistics on the number of matrix products are slightly better. In other tests we have
found the complex choice of Z always to perform at least as well, overall, as the real
choice (see, for example, the riffle shuffle example in section 4). We therefore keep Z
complex in Algorithm 2.4. In version 3.0 of LAPACK the norm estimator has been
modified from that in version 2.0 to keep the vector z complex.
Finally, how does Algorithm 2.4 compare with the suggestion of Hager mentioned
at the beginning of section 2 of running the 1-norm power method t times in suc-
cession? In tests with random matrices we have found Hager's approach to produce
surprisingly good norm estimates, but they are generally inferior to those from Algorithm
2.4. Since Hager's approach is based entirely on level 2 BLAS operations
Algorithm 2.4 is clearly to be preferred.
4. Computing 1-Norm Pseudospectra. In this section we apply the complex
version of Algorithm 2.4 to the computation of 1-norm pseudospectra. For ffl - 0 and
N. J. HIGHAM AND F. TISSEUR

Table
Results for 5000 matrices inv(rand+i*rand) of dimension 100, with
Underest. ratio Products
t min average % exact average max improve%
9 0.859 1.000 99.86 4.0 4 100.00

Table
Results for 5000 matrices inv(rand+i*rand) of dimension 100, with
Underest. ratio Products
t min average % exact average max improve%
9 0.819 0.999 98.12 4.0 6 99.74
any subordinate matrix norm the ffl-pseudospectrum of A 2 C n\Thetan is defined by [22]
z is an eigenvalue of A +E for some E withkEk - ffl g
or, equivalently, in terms of the resolvent (zI
Most published work on pseudospectra has dealt with the 2-norm and the utility of 2-
norm pseudospectra in revealing the effects of non-normality is well appreciated [19],
The 2-norm and any other p-norm of an n \Theta n matrix differ by a factor at most
n. For small n, pseudospectra therefore do not vary much between different p-norms.
However, J'onsson and Trefethen have shown [17] that in Markov chain applications
the choice of norm for pseudospectra can be crucial. In Markov chains representing
a random walk on an m-dimensional hypercube and a riffle shuffle of m cards, the
transition matrices are of dimension exponential in m and factorial in m, respectively.
It is known that when measured in an appropriate way, these random processes converge
to a steady state not gradually but suddenly after a certain number of steps.
As these processes involve powers of matrices of possibly huge dimension, and as the
1-norm is the natural norm for probability 1 , 1-norm pseudospectra are an important
tool for explaining the transient behaviour [17].
One of the most useful graphical representations of pseudospectra is a plot of level
curves of the resolvent. We therefore consider the standard approach of evaluating
on an equally spaced grid of points z in some region of interest in the
complex plane and sending the results to a contour plotter. A variety of methods for
carrying out these computations for the 2-norm are surveyed by Trefethen [21], but
most of the ideas employed are not directly applicable to the 1-norm.
Explicitly forming (zI at each grid point is computationally expensive.
A more efficient approach is to factorize P LU at each point z by LU
factorization with partial pivoting and to use Algorithm 2.4 to estimate k(zI \Gamma
the matrix multiplications in the algorithm become triangular solves with multiple
right-hand sides. This approach can take advantage of sparsity in A. However, the
method still requires O(n 3 ) operations per grid point when A is full.
We consider instead a more efficient approach that is applicable when we can
compute a Schur factorization [9, Chap. 7]
where Q is unitary and T is upper triangular. Given this factorization we have
so that forming a matrix product with (zI its conjugate transpose reduces
to solving a multiple right-hand side triangular system and multiplying by Q and
. Given this initial decomposition the cost per grid point of estimating the resolvent
norm using Algorithm 2.4 is just O(n 2 t) flops-a substantial saving over the
first approach, and of the same order of magnitude as the cost of standard methods
for computing 2-norm pseudospectra (provided t is small). In place of the Schur
decomposition we could use a Hessenberg decomposition, computed either by Gauss
transformations or Householder transformations [9, Sec. 7.4]; these decompositions
are less expensive but the cost per grid point is a larger multiple of n 2 flops because
of the need to factorize a Hessenberg matrix.
In more detail our approach is as follows.
Algorithm 4.1 (1-norm pseudospectra estimation). Given A 2 C n\Thetan and a
positive integer t, this algorithm estimates k(zI on a specified grid of points
in the complex plane, using Algorithm 2.4 with parameter t.
Compute the Schur factorization
for each grid point z
Apply the complex version of Algorithm 2.4 to (zI
with parameter t, using the representation (4.2).
Our experience is that Algorithm 4.1 frequently leads to visually acceptable contour
plots even for 1. As the following example shows, however, a larger value of t
may be needed. We take an example from Markov chains: the Gilbert-Shannon-Reeds
model of a riffle shuffle on a deck of n cards. The transition matrix P is of dimension
n!. Remarkably, as J'onsson and Trefethen explain [17], the dimension of the matrix
1 Because of the use of row vectors in the Markov chain literature, it is actually the 1-norm that
is relevant. By using k1 we can continue to work with the 1-norm.
14 N. J. HIGHAM AND F. TISSEUR
can be reduced to n by certain transformations that preserve the 1-norms of powers
and of the resolvent. For this experiment we took and, as in [17], we computed
pseudospectra of the "decay matrix" working with the reduced
form of A.

Figure

4.1 shows approximations to the 1-norm pseudospectra computed on a
100 \Theta 100 grid, with 2.4. Contours are plotted for
the dashed line marks the unit circle, and the eigenvalues are
plotted as dots. We did not exploit the fact that, since A is real, the pseudospectra
are symmetric about the real axis. The contour plot for clearly incorrect in
the outer contour, while yields an improvement and the plot for
to visual accuracy. Table 4.1 summarizes the key statistics from these computations,
showing that on average the norm estimates had about t correct significant digits for
3. The figure and the table confirm that it is better to keep Z complex in
Algorithm 2.4.
We give a further example, in which A is a spectral discretization of an integral
operator of Landau [21, Sec. 21]; like the operator, A is complex and symmetric. We
took dimension Fresnel number 8. As noted in [21] for the 2-norm,
a fine grid is needed to resolve the details for this example; we used a 200 \Theta 200
grid.

Figure

4.2 shows the computed pseudospectra for
summarizes the statistics for these values and for 2. Contours are plotted for
the dashed line marks the unit circle, and the eigenvalues
are plotted as dots. In Table 4.2, "11   " denotes that convergence was declared because
the iteration limit was reached (the percentage of such occurrences was 0.12% for
and 0.025% for one of the contour lines misses an eigenvalue in
the north-west corner of the plot. For the plot differs from the exact one only
by some tiny oscillations in two outer contours. While Algorithm 2.4 performs well
even for small t, as measured by the underestimation ratio, quite accurate values of
the 1-norm of the resolvent are needed in this example in order to produce smooth
contours.
5. Conclusions. We have derived a new matrix 1-norm estimation algorithm,
Algorithm 2.4, with a number of key features. Most importantly, the algorithm has a
parameter t that can be used to control the accuracy and reliability of the estimate
(which is actually a lower bound). While there is no guarantee that increasing t increases
the estimate (leaving aside the fact that the starting matrix is partly random),
the estimate typically does increase with t, leading quickly to one or more correct significant
digits in the estimate. A crucial property of the algorithm is that the number
of iterations and matrix products required for convergence is essentially independent
of t (for random matrices about 2 iterations are required on average, corresponding
to 4 products of n \Theta n and n \Theta t matrices). The algorithm avoids redundant computations
and keeps constant the size of the matrix multiplications. In future work we
intend to investigate how the choice of t affects the efficiency of the algorithm in a
high performance computing environment.
Unlike the statistically-based norm estimation techniques in [7], [11] which currently
apply only to real matrices, our algorithm handles both real and complex
matrices.
Since our algorithm uses a partly random starting matrix for t - 2, it is natural
to ask whether bounds, valid for all A, can be obtained on the probability of the
estimate being within a certain factor of kAk 1 . We feel that the very features of the
algorithm that make it so effective make it difficult or impossible to derive useful
bounds of this type.
For our algorithm is very similar to the estimator in LAPACK, the differences
being that our algorithm omits the extra estimate (2.1) and for real matrices
we test for parallel rather than simply repeated sign vectors (which improves the ef-
ficiency). Unlike in the estimator in LAPACK 2.0, for complex matrices we do not
take the real part of the z vector (which improves both the quality of the estimates
and the efficiency), and this change has been incorporated into LAPACK 3.0.
The new algorithm makes an attractive replacement for the existing LAPACK
estimator. The value would be the natural default choice (with the extra
estimate (2.1) included for extra reliability and backward compatibility) and a user
willing to pay more for more accurate and reliable 1-norm estimates would have the
option of choosing a larger t.

Acknowledgements

. We thank Nick Trefethen for providing M-files that compute
the riffle shuffle and Landau matrices used in section 4.
N. J. HIGHAM AND F. TISSEUR

Table
Results for riffle shuffle example.
Underest. ratio Products
t min average % exact average max
Fig. 4.1. 1-norm pseudospectra for the riffle shuffle. Clockwise from top left:

Table
Results for Landau matrix example.
Underest. ratio Products
t min average % exact average max
Fig. 4.2. 1-norm pseudospectra for the Landau matrix. Clockwise from top left:
N. J. HIGHAM AND F. TISSEUR



--R


Solving sparse linear systems with sparse backward error.
The power method for
An estimate for the condition number of a matrix.
A set of counter-examples to three condition number estimators
Open problems in numerical linear algebra.
Estimating extremal eigenvalues and condition numbers of matrices.
Matrix Computations.
A stable variant of the secant method for solving nonlinear equations.

Condition estimates.
A survey of condition number estimation for triangular matrices.
FORTRAN codes for estimating the one-norm of a real or complex matrix
Experience with a matrix norm estimator.
Accuracy and Stability of Numerical Algorithms.
A numerical analyst looks at the
Convergence of a subgradient method for computing the bound norm of matrices.
Pseudospectra of matrices.
Pseudospectra of linear operators.
Computation of pseudospectra.
Spectra and Pseudospectra: The Behavior of Non-Normal Matrices and Operators
--TR

--CTR
J. R. Cash , F. Mazzia , N. Sumarti , D. Trigiante, The role of conditioning in mesh selection algorithms for first order systems of linear two point boundary value problems, Journal of Computational and Applied Mathematics, v.185 n.2, p.212-224, 15 January 2006
J. R. Cash , F. Mazzia, A new mesh selection algorithm, based on conditioning, for two-point boundary value codes, Journal of Computational and Applied Mathematics, v.184 n.2, p.362-381, 15 December 2005
LAPACK;matrix condition number;p-norm power method;level 3 BLAS;1-norm pseudospectrum;condition number estimation;matrix 1-norm;matrix norm estimation