The Chebyshev Polynomials of a Matrix.
A Chebyshev polynomial of a square matrix A is a monic polynomial p of specified degree that minimizes |p (A)|2.  The study of such polynomials is motivated by the analysis of Krylov subspace iterations in numerical linear algebra.  An algorithm is presented for computing these polynomials based on reduction to a semidefinite program which is then solved by a primal-dual interior point method.  Examples of Chebyshev polynomials of matrices are presented, and it is noted that if A is far from normal, the lemniscates of these polynomials tend to approximate pseudospectra of A.
Introduction
. Let A be an N \Theta N matrix and n a nonnegative integer. The
degree n Chebyshev polynomial of A is the unique monic polynomial p
n of degree n
such that
(1)
denotes the matrix 2-norm. To be precise, p
n is unique provided that n is
less than or equal to the degree of the minimal polynomial of A; otherwise we have
and the problem ceases to be interesting.
This notion of the polynomial that minimizes kp(A)k seems so simple and natural
that one would expect it to be a standard one. We suspect it may have been considered
before, perhaps decades ago in the literature of approximation theory. Nevertheless,
we have been unable to find any literature on this problem before our 1994 paper with
Greenbaum [7]. In that paper, Chebyshev polynomials of matrices are defined and it
is proved that they exist (obvious by compactness) and that they are unique under
the condition just mentioned (not obvious).
Even if they are not discussed explicitly, Chebyshev polynomials of matrices are
never far away from any discussion of convergence of Krylov subspace iterations in
numerical linear algebra. For these iterations, convergence depends on certain vector
norms kp(A)qk being as small or nearly as small as possible, where q is a starting
vector. Most of the convergence properties of applied interest do not depend too
strongly on q, and thus it is the near-minimality of kp(A)k that is often the heart
of the matter [22]. For finding eigenvalues, the principal iterative method in this
category is the Arnoldi iteration, which becomes the Lanczos iteration if A is real and
symmetric. For solving systems of equations, the analogous methods include GMRES,
biconjugate gradients, CGS, QMR, and Bi-CGSTAB in the general case and conjugate
gradients if A is symmetric positive definite. (For systems of equations, the notion
of a Chebyshev polynomial of A should be normalized differently by the condition
instead of the condition that p is monic. In [7], a Chebyshev polynomial
The work of both authors was supported by NSF grant DMS-9500975CS and DOE grant DE-
FG02-94ER25199.
y Department of Mathematics, National University of Singapore, 10 Kent Ridge Crescent, Singapore
z Oxford University Computing Laboratory, Wolfson Building, Parks Road, Oxford OX1 3QD,
United Kingdom (Nick.Trefethen@comlab.ox.ac.uk).
K.-C. TOH AND L. N. TREFETHEN
of a matrix is called an ideal Arnoldi polynomial, and its analogue with this other
normalization is called an ideal GMRES polynomial.)
The motivation for the term "Chebyshev polynomial of a matrix" is as follows.
All readers will be familiar with the classical Chebyshev polynomials fTn g, which are
times monic polynomials of minimal k \Delta k1-norm on the interval [\Gamma1; 1]. This
notion was generalized by Faber in 1920 to the idea of the Chebyshev polynomials
of S, where S is a compact set in the complex plane C: the monic polynomials of
Now suppose that A is a hermitian or more
generally a normal matrix, having a complete set of orthogonal eigenvectors. Then
by a unitary reduction to diagonal form, it is easily shown that the nth Chebyshev
polynomial of A as defined by (1) is precisely the nth Chebyshev polynomial of S in
this latter sense, where S is the spectrum of A. Such a polynomial can be computed,
for example, by generalizations of the Remez algorithm [15].
Chebyshev polynomials of normal matrices, then, are trivial; the matrix problem
reduces to a scalar problem. But what if A is an arbitrary square matrix, with
non-orthogonal eigenvectors or perhaps no complete set of eigenvectors? This is the
subject of this paper, and our purpose is two-fold.
First, we describe an algorithm for computing Chebyshev polynomials of matrices.
The optimization problem implicit in (1) is far from smooth, and unless the degree
is very small, these problems are quite difficult if approached by general methods of
unconstrained optimization. The algorithm we describe, which we believe is the first
to have been developed for this problem, is based instead on interior point methods for
semidefinite programming. With this algorithm, we can reliably compute Chebyshev
polynomials for matrices of order - 50 in less than a minute on workstations available
in 1996. No parameters are involved that must be tuned. We should mention, however,
that although our algorithm is reasonably fast, it is not fast enough to easily handle
matrix dimensions of the order of a thousand or more.
Second, we present computed examples, the first we know of to have been pub-
lished. A few numerical coefficients are listed for possible comparison by later authors,
but our main aim is to give insight into the behavior of Chebyshev polynomials of
matrices, largely with the aid of pictures. A natural question is, how are the coefficients
of the polynomials affected by the degree and nature of the non-normality of
A? For a partial answer, we plot lemniscates jp
constant of our polynomials
and find that in many cases they approximate pseudospectra of A.
2. Reduction to a semidefinite program. Let fB be a linearly
independent set of matrices in C N \ThetaN . The Chebyshev problem (1) is a special case
of a norm minimization problem involving linear functions of matrices:
min
(2)
For our special case, B
and the numbers x k are the coefficients (actually their negatives) of the Chebyshev
polynomial of A.
It is well known that (2) can be expressed as a semidefinite program [11], [23]. We
shall not show in detail how this is done. One difference between our work and what
CHEBYSHEV POLYNOMIALS OF A MATRIX 3
has been done before is that the existing literature, as far as we are aware, considers
only real matrices.
Theorem 1. The norm minimization problem (2) is equivalent to the following
semidefinite program involving hermitian matrices:
s:t:
I
\GammaiB
means that Z is positive semidefinite.
Proof. Problem (2) is equivalent to the problem of minimizing \Gamma- such that
Using the fact that for any M 2 C N \ThetaN ,
'-
where -max (\Delta) denotes the maximum eigenvalue, equation (5) can be rewritten as
'-
. But this is equivalent to
I
By writing this equation out in full, we get (3).
3. Transformation to a better-conditioned basis. Before we discuss how
the semidefinite program (3) can be solved by interior point methods, we must address
the issue of change of basis in (2), as the numerical stability of these algorithms
depends on the conditioning of the basis fB g. This is an essential point for
the computation of Chebyshev polynomials of matrices. The power basis is usually
highly ill-conditioned, as can be seen by considering the special case of a diagonal
matrix, where we get the phenomenon of ill-conditioning of the basis of monomials
fx k g, familiar in approximation theory. For numerical stability in most cases, the
power basis must be replaced by a better-conditioned alternative.
4 K.-C. TOH AND L. N. TREFETHEN
Suppose f b
Bn g is another linearly independent set of matrices in C N \ThetaN
related linearly to fB by
Bn
Bn
t;
where T is an n \Theta n nonsingular matrix, t is an n-vector, and c is a nonzero scalar.
(The notation here means that
Bn ).) The following theorem describes how (2) is modified by
this change of basis. The proof is straightforward, and we shall omit it.
Theorem 2. The minima
min
and
are the same, and the unique vectors x and - x that achieve them are related by
We are aware of three choices of basis that are particularly attractive for practical
computations.
Scaled power basis. Suppose are given as in (2). A simple way to get
a better conditioned basis is to scale the norm of A to 1. With such a scaling, we
have b
and
Faber polynomial basis. Even the best-scaled power basis is often highly ill-
conditioned. A more powerful idea is to consider the basis fF defined
by the Faber polynomials F associated with some
region\Omega in the complex
plane containing the spectrum of A. The Faber polynomials fFng are the natural
analogues for a general
region\Omega in C of the monomials fz n g for the unit disk or the
Chebyshev polynomials fTn g for [\Gamma1; 1]; see [2]. In most cases, fFn (A)g will be far
better conditioned than any power basis.
For the Faber basis, the matrix T in (2) is upper triangular, with columns containing
the coefficients of F . The scalar c is the positive number cap(\Omega\Gamma n ,
where cap(\Omega\Gamma is the logarithmic capacity of \Omega\Gamma The vector t is the vector of co-efficients
of the expansion of the degree in terms of
Of course, one must choose a
region\Omega for which the associated Faber polynomials
can be obtained either analytically or numerically.
If\Omega is chosen to be an ellipse or
CHEBYSHEV POLYNOMIALS OF A MATRIX 5
an interval, then the Faber polynomials are simply the scaled Chebyshev polynomials
g. More generally,
if\Omega is chosen to be a polygonal domain, the Faber polynomials
can be computed numerically via Schwarz-Christoffel mapping. We have used the
Matlab Schwarz-Christoffel Toolbox for this purpose, due to Driscoll [4].
Orthonormal basis. Finally, our third idea is a more elementary one, but powerful
in practice. One may simply orthonormalize the power basis fI; with
respect to the "trace inner product" hA; to obtain a basis
that is typically well-conditioned even in the 2-norm. This can
be done by a modified Gram-Schmidt procedure similar to that used in the Arnoldi
iteration:
To obtain the matrix T in (2), we note that there is a unique (n
triangular matrix R such that
and the columns of R can be computed from the following recurrence relation (in
Matlab
It is now easy to
see that
again in Matlab notation.
For simplicity, we use the orthonormal basis in the examples reported in this
paper. Although it is more expensive to compute than the other two bases, the
amount of time taken remains small compared to the time required for solving (3).
We note that transformation to a better-conditioned basis does not eliminate any
ill-conditioning that is inherent in the Chebyshev minimization problem itself.
4. Solution by primal-dual interior point method. Assuming a suitable
basis has been chosen, we now turn to the problem of how (3) can be solved by interior
point methods similar to those in linear programming, specifically, by Mehrotra-type
primal-dual predictor-corrector algorithms. Extensive research has been done on both
the algorithms and the theory of SDP. We refer the reader to [1], [9], [10], [11], [12],
[16], [23] and [26] for details.
6 K.-C. TOH AND L. N. TREFETHEN
A general SDP has the form:
s:t:
where C, Z, A k , hermitian matrices and b 2 IR n . The idea
behind an interior point method is to use a suitable barrier function, \Gamma log det(Z) in
the case of SDP, to transform the semidefinite constrained convex problem (D) into
a parametrized family (by -) of equality constrained convex problems whose optimal
solutions (X(-); y(-); Z(-)) satisfy the optimality conditions
where X and Z are hermitian positive definite. The parameter - ? 0 is to be driven
explicitly to zero (as fast as possible), and in the limit - ! 0, an optimal solution of
(6) is obtained.
Mehrotra-type primal-dual predictor-corrector algorithms essentially consist of a
sequence of modified Newton iterations. Usually, one step of Newton's iteration is
applied to (7) for each new -.
It is readily shown that application of Newton's method to (7) gives rise to the
equations
In order to keep \DeltaX hermitian (this is desirable since the fundamental objects in an
SDP are hermitian matrices), equation (9) is usually symmetrized with respect to an
invertible matrix P , whereupon it becomes
where
Different choices of P give rise to different Newton steps. For example,
gives rise to what is known as the Alizadeh-Haeberly-Overton (AHO) direction [1];
rise to the Monteiro direction [10]; and
rise to the Nesterov-Todd (NT) direction [12].
The general algorithmic framework of a Mehrotra-type predictor-corrector method
is as follows.
Algorithm. Given an initial iterate (X positive definite,
CHEBYSHEV POLYNOMIALS OF A MATRIX 7
(Let the current and the next iterate be (X;
1. (Predictor step)
Compute the Newton step (ffiX; ffiy; ffiZ) from (8) and (10) with
2. Determine the real parameter
Here ff and fi are suitable steplengths chosen to ensure that X
are positive definite. Generally, ff and fi have the form
3. (Corrector step)
Compute the Newton step (\DeltaX; \Deltay; \DeltaZ ) from (8) and (10) with the right-hand side
matrix R given by
4. Update (X;
where ff and fi are defined by (12) with ffiX, ffiZ replaced by \DeltaX , \DeltaZ .
We shall not discuss implementation details of the above algorithm-for example,
how to solve efficiently for the search directions (ffiX; ffiy; ffiZ) and (\DeltaX; \Deltay; \DeltaZ ) from
the linear systems of equations (8) and (10). We refer the reader to [16]
for such details. Instead, we just note that the search directions are typically computed
via a Schur complement equation. For such an implementation, each iteration
has a complexity of O(nN 3 which is equal to O(nN 3 ) for our Chebyshev
approximation problem since n ! N . Computations have shown that careful
implementations of the predictor-corrector algorithm that use a Schur complement
equation can typically reduce the duality gap of an SDP to about ffl 2=3
mach for the three
search directions mentioned above, namely, the AHO, Monteiro, and NT directions.
For these three directions, each iteration has a complexity of at most 12nN 3 , and
the number of iterations needed to reduce the duality gap by a factor of
exceeds 20.
In all of our computations, we use the NT direction, for the following reasons.
Although the orders of complexity for computing these three directions are the same,
computing the AHO direction is about twice as expensive as computing the Monteiro
or NT directions. Of the latter two, the NT direction has the virtue of being primal-dual
symmetric. This implies that primal-dual predictor-corrector algorithms based
on the NT direction are likely to be more robust than those based on the Monteiro
direction, in the sense that the problems of stagnation such as taking very small
steplengths are less likely to occur. It has been observed that algorithms based on
the Monteiro direction often encounter such a stagnation problem for a class of SDPs
known as the ETP problems [23].
8 K.-C. TOH AND L. N. TREFETHEN
5. The special case when A is normal. It is worth setting down the form our
algorithm takes in the special case where A is normal, i.e., unitarily diagonalizable.
As we have already mentioned in the Introduction, we may assume in this case that
A is diagonal, so that the Chebyshev problem (1) reduces to the classical Chebyshev
approximation problem on the spectrum  (A) of A, i.e.,
For this special case, the Chebyshev polynomials of A can be computed cheaply by
the predictor-corrector algorithm discussed in the last section, by exploiting the block
diagonal structure present in the associated SDP problem.
As in the general case, we consider the norm minimization problem (2), but the
are now diagonal: each k. Since the
2-norm of a diagonal matrix is the k \Delta k1-norm of its diagonal vector, (2) is equivalent
to the minimax problem
min
where d (l)
k denotes the lth component of the N-vector d k . As before, (13) can be
expressed as an SDP.
Theorem 3. The minimax problem (13) is equivalent to the following SDP involving
block diagonal hermitian matrices:
s:t:
diag
d (l)
\Gammai d (l)
n. The matrices A k consist of N blocks of 2 \Theta 2 matrices on the diagonal.
A proof of the above theorem is similar to that of Theorem 1, based on the observation
that for any complex number a, we have
'-
0 a
- a 0
We omit the details. Also, the process of transformation to a better-conditioned basis
for (14) is exactly the same as for the general case. However, note that (14) cannot
be obtained as a direct consequence of Theorem 1 by specializing the matrices B k to
diagonal matrices.
CHEBYSHEV POLYNOMIALS OF A MATRIX 9
If the initial iterate (X chosen to have the same block diagonal structure
as the matrices A k , then this structure is preserved throughout for (X k ; Z k ). By
exploiting this block diagonal structure, the work for each iteration of the predictor-corrector
algorithm is reduced to O(n 2 N) flops as opposed to O(nN 3 ) for non-normal
matrices. In practice, we can compute the degree-25 Chebyshev polynomial of a
normal matrix of dimension 1000 in Matlab in about 12 minutes on a Sun Ultra
Sparcstation.
It would be interesting to know how this special case of our algorithm for normal
matrices compares with other methods for linear complex Chebyshev approximation,
such as the Remez / semiinfinite programming methods discussed in [15], but we have
not investigated this.
6. Computed examples. We turn now to computed examples of Chebyshev
polynomials of matrices. Our aim is to demonstrate the effectiveness of our algorithm
and to give some insight into the behavior of these polynomials. This is not a subject
we fully understand, but the experimental observations are fascinating.
Most of our experimental results will be presented as plots. To "plot" a polynomial
n , we show its roots in the complex plane and also the boundary of a region that
we call the Chebyshev lemniscate 1 for that polynomial and the given matrix A. This
region is defined by the equation
The Chebyshev lemniscates characterize where in the complex plane the Chebyshev
polynomials of A "live", just as the spectrum or the pseudospectra characterize
(though not precisely, unless A is normal) where in the complex plane A itself "lives".
As a minimum, since kp
n (A)k, we know that the Chebyshev lemniscate
contains the spectrum
ae Ln (A):
In each example we present, the dimension of the matrix A is 48 \Theta 48 or 100 \Theta 100,
though we typically print only its 5 \Theta 5 or 6 \Theta 6 analogue. For each example, we give
plots showing the Chebyshev lemniscates (solid curves) of A, typically of degrees
The zeros of the Chebyshev polynomials are shown as small circles, and
the eigenvalues of A are shown as solid dots.
For comparison with the Chebyshev lemniscate, each of our plots also shows a
dotted curve. This is the boundary of an ffl-pseudospectrum of A. The value of ffl has
been chosen "by hand" to make the match with the Chebyshev lemniscate a good
one. (The ffl-pseudospectrum of A is the set   in the
complex plane; see [14] and [21].)
For all of these examples, the Chebyshev polynomials were computed in Matlab
by the methods described in the previous sections.
Primal-dual predictor-corrector algorithms are highly efficient and robust for solving
SDPs. For the set of examples we present here, it takes an average of 12 iterations
to reduce the duality gap by a factor of (This number is rather insensitive to
the dimension of A; it would be essentially the same for matrices of dimensions 5 \Theta 5
or 200 \Theta 200. This insensitivity to problem size is one of the remarkable features of
Properly speaking, the word lemniscate refers to the boundary of Ln , and Ln itself is a
lemniscatic region, but this expression is cumbersome and we shall avoid it.
K.-C. TOH AND L. N. TREFETHEN
primal-dual interior point methods.) For a 48 \Theta 48 real matrix, each iteration takes
about 5 and 7 seconds for respectively, on a Sun Ultra Sparcstation.
The corresponding numbers for a 48 \Theta 48 complex matrix are about seconds and
seconds.
Here are our examples. Omitted entries are all zero.
Example 1. Diagonal.
where d is a vector whose first entry is 1 and the rest of whose entries are distributed
uniformly in the interval [\Gamma1; 0:8]. Thus the spectrum of A consists of points that
densely fill the interval [\Gamma1; 0:8] and an outlier at z = 1.
Example 2. Bidiagonal.
. 0:2
where the vector d is the same as that in Example 1. The spectrum is the same as in
Example 1.
Example 3. Grcar [21].
Example 4. Ellipse.
Example 5. Bull's head [14].
CHEBYSHEV POLYNOMIALS OF A MATRIX 11
Example 6. Lemniscate1 [14].
Example 7. Lemniscate2 [20].
Example 8. Gauss-Seidel [14].
This is the matrix that corresponds to a Gauss-Seidel iteration applied to the standard
3-point discrete Laplacian on a grid of N points.
Example 9. Beam-Warming [20].
\Gamma1:5 2:0 \Gamma1:5
0:7 \Gamma2:6 2:1C C C C C C A
Example 10. Wilkinson [21].
12 K.-C. TOH AND L. N. TREFETHEN
Example 11. Chebyshev points.
where
Example 12. Random [21].
where by random we mean that the entries of A are independently drawn from the
real normal distribution with mean 0 and variance 1=N .
Example 13. Random triangular [21].
by which we mean that A is the strictly upper triangular part of the random matrix
of Example 12.
Example 14. Convection-diffusion matrix [18].
The matrix A is the projection of the 2N \Theta 2N Chebyshev spectral discretization
matrix of a convection-diffusion operator onto the invariant
subspace associated with the N eigenvalues of maximal real part
In

Table

1, for later authors who may wish to compare the coefficients of some
Chebyshev polynomials of matrices, we list the coefficients of p
8 for the matrices of
Examples 3 and 5. In Table 2, we list kp
8 (A)k for all fourteen examples.
The plots for our fourteen examples are shown in Figures 1-14.
Let us first consider Example 1, the special case where A is diagonal. For any
Chebyshev polynomial of a matrix, we know that the Chebyshev lemniscate must
contain the spectrum (16). In the present case, by the characterization theorems for
the classical complex Chebyshev approximation problem [3, p. 143], we know that
the nth Chebyshev lemniscate must in fact touch the spectrum  (A) at no fewer
than points. This property is evident in Figure 1, where we see that Ln (A)
rather closely, and increasingly so as n increases (see the cover illustration
of [22]). It is interesting also to note how quickly one of the roots of the polynomials
n , which are analogous to the "Ritz values" often taken as eigenvalue estimates in
Arnoldi or Lanczos iterations, converges to the outlier eigenvalue at z = 1. By
one of the roots of p
6 is already very close to the outlier, and the distance between
them decreases geometrically as n increases. In the remainder of the spectrum, on
the other hand, no individual Ritz value is converging rapidly to any one eigenvalue
of A. Rather, it is the Chebyshev lemniscate generated by these Ritz values jointly
that is capturing the spectrum.
CHEBYSHEV POLYNOMIALS OF A MATRIX 13

Table
Computed coefficients of p  for the Grcar and Bull's head matrices (Examples 3 and 5). All
but perhaps the last two digits printed are believed to be correct.
Grcar Bull's head
computed
computed

Table
Norms kp  (A)k for Examples 1-14. All digits printed are believed to be correct, as the estimated
relative accuracies are all less than 10 \Gamma11 .
Example computed kp
1. Diagonal 0:0063675408
2. Bidiagonal 0:0551494047
3. Grcar 1766:3135313
4. Ellipse 7710:2711611
5. Bull's head 1239:4186173
6. Lemniscate1 1:0000000000
7. Lemniscate2 834:73857463
8. Gauss-Seidel 0:0049251285
9. Beam-Warming 7:4348443860
10. Wilkinson 6:2747795054
11. Chebyshev points 46:395131600
12. Random 2:9537221027
13. Random triangular 0:0039633789
14. Convection-diffusion 2623904:6097
One might expect Ln (A) to approximate  (A) even if A is non-normal. But from

Figures

2-14, the reader will see that this does not happen. Nonetheless, though
does not always approximate  (A) very closely, it still gains some information
about A. The plots show that for these examples, to a rather startling degree,
for some ffl - 0, where   ffl (A) is again the ffl-pseudospectrum of A. In particular, the
agreement of the Chebyshev lemniscate of p
n with a pseudospectrum of A is far closer
in most of these examples than the agreement of the roots of p
n with the eigenvalues
14 K.-C. TOH AND L. N. TREFETHEN
Fig. 1. Diagonal. Since A is normal, the Chebyshev lemniscate touches the spectrum at at least
and the roots of p
n lie in the convex hull of the spectrum.
Fig. 2. Bidiagonal-a non-normal analogue of Example 1. The dotted curves are the
pseudospectrum of A.
of A. For example, consider Figure 2, the bidiagonal matrix that is the non-normal
analogue of Example 1 with the same spectrum. Except for the outlier eigenvalue, the
roots of p
bear no resemblance to individual eigenvalues of A. On the other hand,
the Chebyshev lemniscates of these polynomials show a striking resemblance to the
pseudospectrum of A. Clearly the Chebyshev polynomial is approximating
A in a fashion that goes beyond approximation of individual eigenvalues.
CHEBYSHEV POLYNOMIALS OF A MATRIX 15
The other examples illustrate the same effect. In every case, the lemniscate of
the Chebyshev polynomial shows a compelling approximation to the pseudospectrum.
We do not claim that this effect is universal; these examples have been picked for their
pronounced and cleanly structured non-normality. But it is certainly common.
A partial explanation of this phenomenon is as follows. It is well known that a
matrix polynomial p(A) can be expressed as a Cauchy integral
Z
where the integration is over any closed contour or union of contours enclosing the
spectrum of A once in the counterclockwise direction [8]. Taking absolute values gives
the inequality
Z
Now suppose we seek p such that kp(A)k is small. When the degree of p is smaller
than the dimension of A, it is impossible to achieve this in general by putting zeros of
wherever A has eigenvalues, which would make the integral zero. Instead, we must
settle for making jp(z)j small where k(zI large. This immediately suggests
a link between lemniscates of p and pseudospectra of A.
From this kind of reasoning we can derive bounds on kp
(A)k. For example, to
minimize kp(A)k one might seek to minimize kpk   ffl (A) for some ffl that is not too small
(hence jp(z)j is small over the region where k(zI larger
(18) and the minimality of kp
n (A)k we conclude that
where L ffl is the arclength of the boundary of   ffl (A). At this point one runs into the
fact that min p kpk   ffl (A) can be huge if ffl is not small, since the minimum typically
increases geometrically with ffl. Therefore, a compromise must be made on ffl so that
the quantity min p kpk   ffl (A) =ffl on the right-hand side of (19) is as small as possible.
For some matrices A and choices of n and ffl, the estimate just described can be
quite good. It is not always very good, however, and so far, our attempts to make a
more precise link between lemniscates of Chebyshev polynomials and pseudospectra of
the underlying matrix have been unsuccessful except in certain limiting cases n !1
described in [17]. Rather than present partial results that do not bring this matter to
a very satisfactory state, we prefer to leave the explanation of the behavior of Figures
2-14 as an open problem.
7. Conclusions. This paper has made two contributions. The first is a reasonably
fast and apparently robust algorithm for computing the Chebyshev polynomials
of a matrix, based on primal-dual interior point methods in semidefinite program-
ming. The second is an experimental study of these polynomials that indicates that
the associated lemniscates sometimes closely approximate certain pseudospectra of A.
We have said little about applications in iterative numerical linear algebra, though
that was our original motivation. There are many possibilities here that might be explored
now that an algorithm is available. For example, our algorithm may prove
useful in analyzing the convergence of Krylov subspace iterations, or the construction
of preconditioners for such iterations, by means of model problems of moderate
dimension.
K.-C. TOH AND L. N. TREFETHEN
It was mentioned in the Introduction that for applications to iterative solution of
equations rather than eigenvalue calculations, it is appropriate to minimize kp(A)k
with the normalization 1. Plots of lemniscates
for these "ideal GMRES polynomials" can be found in the first author's dissertation
[17]. Because this normalization gives a special status to the origin, these problems
are no longer translation-invariant in the complex plane, and the lemniscates
take special pains to avoid the origin. They also tend to display scallop patterns near
the spectrum or pseudospectra.
Interesting connections can also be made to the notion of a generalized Kreiss
matrix theorem. The usual Kreiss matrix theorem relates the norms kA n k to the
behavior of the pseudospectra of A near the unit disk. Generalizations are obtained by
looking at norms kp(A)k for other polynomials p and the behavior of the pseudospectra
near other regions. These matters are investigated in [19].
We consider the idea of the Chebyshev polynomials of a matrix a natural one,
suggesting many questions to be explored. We hope that more will be learned about
the behavior of these polynomials in the years ahead and that applications in other
areas will be found.



--R


Faber polynomials and the Faber series
Interpolation and Approximation
A MATLAB Toolbox for Schwarz-Christoffel mapping

Matrix Computations
GMRES/CR and Arnoldi/Lanczos as matrix approximation problems
Perturbation Theory for Linear Operators
An infeasible start predictor-corrector method for semi-definite linear programming



Pseudospectra of the convection-diffusion operator
Eigenvalues and pseudo-eigenvalues of Toeplitz matrices
A fast algorithm for linear complex Chebyshev approximations

Matrix Approximation Problems and Nonsymmetric Iterative Methods
Calculation of pseudospectra by the Arnoldi iteration
The Kreiss matrix theorem on a general complex domain
Spectra and Pseudospectra: The Behavior of Non-Normal Matrices and Op- erators
Pseudospectra of matrices


User's guide to SP: software for semidefinite programming
Extremal polynomials associated with a system of curves in the complex plane
On extending some primal-dual interior-point algorithms from linear programming to semidefinite programming
--TR

--CTR
Michel X. Goemans , David Williamson, Approximation algorithms for MAX-3-CUT and other problems via complex semidefinite programming, Proceedings of the thirty-third annual ACM symposium on Theory of computing, p.443-452, July 2001, Hersonissos, Greece
Michel X. Goemans , David P. Williamson, Approximation algorithms for MAX-3-CUT and other problems via complex semidefinite programming, Journal of Computer and System Sciences, v.68 n.2, p.442-470, March 2004
matrix polynomial;krylov subspace iteration;chebyshev polynomial of a matrix;semidefinite programming