BIB-VERSION:: CS-TR-v2.0
ID:: CORNELLCS//TR93-1360
ENTRY:: 1993-10-14
ORGANIZATION:: Cornell University, Computer Science Department
LANGUAGE:: English
TITLE:: Pseudozeros of Polynomials and Pseudospectra of Companion Matrices
AUTHOR:: Toh, Kim-Chuan 
AUTHOR:: Trefethen, Lloyd N.
DATE::  June 1993
PAGES:: 27
ABSTRACT::
It is well known that the zeros of a polynomial $p$ are equal to the 
eigenvalues of the associated companion matrix $A$. In this paper, we take a 
geometric view of the conditioning of these two problems and of the stability 
of algorithms for polynomial zerofinding. The $\epsilon$-pseudozero set 
$Z_{\epsilon}(p)$ is the set of zeros of all polynomials $\hat{p}$ obtained by 
coefficientwise perturbations of $p$ of size $\leq \epsilon$; this is a 
subset of the complex plane considered earlier by Mosier, and is bounded by a 
certain generalized lemniscate. The $\epsilon$-pseudospectrum 
$\Lambda_{\epsilon}(A)$ is another subset of C defined as the set of 
eigenvalues of matrices $\hat{A}=A+E$ with $||E|| \leq\epsilon$; it is 
bounded by a level curve of the resolvent of $A$. We find that if $A$ is 
first balanced in the usual EISPACK sense, then $Z_{\epsilon||p||}(p)$ and 
$\Lambda_{\epsilon||A||}(A)$ are usually quite close to one another. It 
follows that the Matlab ROOTS algorithm of balancing the companion matrix, 
then computing its eigenvalues, is a stable algorithm for polynomial 
zerofinding. Experimental comparisons with the Jenkins-Traub (IMSL) and
Madsen-Reid (Harwell) Fortran codes confirm that these three algorithms have 
roughly similar stability properties.

Key words: polynomial zeros, comapnion matrix, pseudospectrum.
END:: CORNELLCS//TR93-1360
BODY::
Pseudozeros of Polynomials and
Pseudospectra of Companion Matrices
Kim-Chuan Toh*
Lloyd N. Trefethen**
TR 93-1360
June1993
Department of Computer Science
Cornell University
Ithaca, NY 14853-7501
*Center for Applied Mathematics, Cornell University, Ithaca, NY 14853
(kc?cs.cornell.edu) Supported by NSF Grant DM59116110.
**Department of Computer Science, Cornell University, Ithaca, NY 14853
(LNT?cs.comell.edu) Supported by NSF Grant DM59116110.
Pseudozeros of Polynomials and Pseudospectra
of Companion Matrices
Kim-Chuan Toli *
Lloyd N. Trefethen
Abstract
It is well known that the zeros of a polynomial p are equal to the eigenvalues
of the associated companion matrix A. In this paper we take a geometric view
of the conditioning of these two problems and of the stability of algorithms for
polynomial zerofinding. The ?-pseudozero set Zc(p) is the set of zeros of all
polynomials p obtained by coefficientwise perturbations of p of size <e this is
a subset of the complex plane considered earlier by Mosier, and is bounded by a
certain generalized lemniscate. The &psendospectritm Af(A) is another subset
of C defined as the set of eigenvalues of matrices A = A + E with II Eli < cit is
bounded by a level curve of the resolvent of A. We find that if A is first balanced
in the usual EISPACK sense, then and A?jjAjj(A) are usualiy quite close
to one another. It follows that the Matlab ROOTS algorithm of balancing
the companion matrix, then computing its eigenvalues, is a stable algorithm
for polynomial zerofinding. Experimental comparisons with the Jenkins-Traub
(IMSL) and Madsen-Reid (Harweli) Fortran codes confirm that these three
algorithms have roughly similar stability properties.
Key words: polynomial zeros, companion matrix; pseudospectrum.
AMS(MOS) Subject classification. 65F35.
*Center for Applied Mathematics, Cornell University, Ithaca, NY 14853 (kcecs cornell.edu)
Supported by NSF Grant DMS-9116110.
t Department of Computer Science, Cornell University, Ithaca, NY 14853 (LNTccs. cornell.edu)
Supported by NSF Grant DMS-9116110.
1 Introduction
Zeros of polynomials and eigenvalues of nonsymmetric matrices are well-known ex-
amples of problems whose answers may be highly sensitive to perturbations. The
sensitivity of these two problems was made famous among numerical analysts by
Wilkinson in the early 1960s, and contributed to his development of the notions of
stability and conditioning that are now standard in this field [19, 20]. And, of course,
the two problems are related, for the zeros of a polynomial are the same as the
eigenvalues of the associated companion matrix.
Despite the classical nature of the subject, the relationship between these two
problems has received less study than one might suppose. Polynomial zerofinding
has been something of a backwater in numerical analysis, and it is probably fair to
say that although all numerical analysts know that one can find zeros via companion
matrices in principle, most assume that it isn't a good idea to do so. It seems to
be not well known, though it is pointed out in the original paper on the subject [6],
that the central step of the Jenkins-Traub algorithm for polynomial zerofinding is
equivalent to a generalized Rayleigh quotient iteration for finding an eigenvalue of
the companion matrix (for details see the Appendix).
The approach we take in this paper is a geometric one. For a monic polynomial
p(z), let Z(p) denote the zero set (= set of zeros) of p(z), and, for any ? > 0, define
the &pseudozero set of p(z) by
Zc(p) = ? C : z ? Z(p) for some pJ,
(1)
where p ranges over polynomials whose coefficients are those of p modified by per-
turbations of size < e (for details see Section 2). The relevance of such sets to the
conditioning of the zerofinding problem has been discussed by Mosier [13]. Similarly,
for a matrix A, let A(A) denote the spectrum of A, and define the &pseudospectrurn
of A by
Af(A) = C C : z E A(A + E) for some E with IlEli < E?. (2)
Godunov [3], and others
Matrix pseudospectra have been studied by Trefethen [18]?
going back at least to H.J. Landau in 1975 [9].
In this paper we report numerical experiments that show that Z?11p1j(p) and A?jjAj;(A)
are generally quite close to one another when A is a companion matrix of p that has
been "balanced" in the usual EISPACK sense proposed originally by Parlett and
Reinsch [15]. It follows that the polynomial zerofinding problem and the balanced
matrix eigenvalue problem are comparable in conditioning and therefore that finding
zeros via eigenvalues of companion matrices, the method used by the Matlab ROOTS
command, is a stable algorithm. To test these conclusions we compare ROOTS with
the Jenkins-Traub (IMSL) code CPOLY [8] and the Madsen-Reid (llarwell) code
PA16 [10] for a variety of polynomials. Our experiments suggest that all three codes
are reliable, with the highest accuracy typically achieved by PAl6 and the next-best
accuracy by ROOTS.
2
The significance of pseudozero sets and pseudospectra is not just a matter of
rounding errors and numerical stability. In any applied mathematical problem that
apparently depends on polynomial zeros, it is likely that what really matters is
whether p(z) is very small, not necessarily exactly zero. Similarly, in a matrix eigen-
value problem, what really matters may be whether II(zI --H A)-1 is very large, not
necessarily exactly infinity. Thus the investigation of pseudozero sets and pseudospec-
tra is a natural extension of the zerofinding and eigenvalue problems themselves, not
just a tool for analyzing numerical algorithms.
2 Pseudozero sets of polynomials
We begin this section by introducing some notation:
IF) : The set of monic polynomials of degree n with complex coefficients. For each
p ? IF), we will express p as
p(z)=Zciz?, c?=1.
i=O
We will also denote the vector of coefficients ..... . , c??1 )T by p when there is
no danger of confusion.
Z(p) : The set of zeros of p.
: The reciprocal polynomial of p, i.e., p*(z) = z?p(z--H1).
D The set of n x n diagonal matrices with diagonal entries in C. For each D ?
V, we denote its diagonal vector ....... , dn?i)T by d; d-1 denotes the vector
(d0-1			d-l1)T
lixild : For each nonsingular D c V, the norm IIxild is defined over C? by
lixild = IIDxII2, x e
: For each z ? C, z denotes the vector (1, z,.. . , z??1 )T.
..... . , e? : the standard basis in C?
For a fixed D ? V, we can assign a metric on IF) by
lip--Hplid =
fl--Hi			1/2
? IdiI2I?'--Hc??
which measures the perturbations in the coefficients of p relative to the weights given
by d. Note also that the norm liXild induces a matrix norm on C??? (the space of all
n x n matrices over C), which satisfies the formula lAid = IlDAD-1Il2.
3
For each p ? II), we define the t-pseudozero set of p(z) by
Z?(p;d) = E C: ZE Z(p) for some pE II) with IIP?PIId <
(3)
These sets quantify the conditioning of the zerofinding problem and arise naturally
in analyzing the stability of zerofinding algorithms. For a zerofinding algorithm
to be stable, the computed zeros of p should lie in a region Zc?(p; d) for an ap-
propriate choice of d, where C = 0(Ilpild) and u is the machine precision. For
example, the choice d = IIPII2P?1 corresponds to coefficientwise perturbations and
d = ??.,..., 1)T corresponds to normwise perturbations.
The following proposition gives an algebraic characterization of Z?(p; d) in terms of
the level curves of a certain function. This allows us to determine Z?(p; d) numerically.
Proposition 2.1
Z?(p;d) = ? C: Ip(z)I <
IIZIId?1
Mosier [13] introduced these sets and proved this proposition in the oo-norm. llere we
extend the result to the 2-norm; in fact, it is valid in the p-norm for any 1 ? p < 00.
Proof. If z E Z?(p; d), then ? E IF) with p(z) = 0 and lip --H plid < e. By the H5lder
inequality,
lp(z)l = Ip(z)--Hp(z)I = IZ(c?--Hcj)z?I
< lip --H pild iiziid--H1 < Eiiziid--H1.
To show the converse, consider z E C such that p(z) < tiiziid-1. Let 0 = arg(z).
Consider the degree n --H 1 polynomial r defined by
r(w)=?ft1
k
r?wk, where rk =1zk1 e?ikOd??2
Then r(z) = iirii?iizii?-i, and thepdynomiJp E IF) defined byp(w) = p(w)?P%iZz) r(w)
satisfies p(z) =0 and lip--Hpild = iIPjzZ?)iIiiriid ? ?. Thus z E Z?(p;d). El
Note that if d = e1?1 =(1,00,... , 00)T, then Z?(p; d) becomes the region bounded
by the usual &lernniscate of p [5],
L?(p)=fzE C: ip(z)i<?J.
More generally we may call the boundary of Z?(p; d) a generalized lemniscate.
Example 2.1. Consider the monic polynomialp of degree 10 with zeros 1,2,... , 10.
Figure 1 shows the eiipii?-pseudozero sets of p for two different choices of d correspond-
ing to coefficientwise (d = iiPii2P?1) and normwise (d = 7n(1,..., 1)T) perturbations.
The ideas described above apply to arbitrary polynomials and finite perturbations.
For the limiting case of infinitesimal perturbations, a single condition number suffices
4
40			4
ol
.40			.4
0
40			80
(a) normwise
0			5			10
(b) coefficientwise
Figure 1: ?IIpII?-pseudozero sets (e = 1o-?, 10-6) corresponding to normwise and Co-
efficientwise perturbations for the degree-tO monic polynomial with zeros 1,2,.. . , 10.
Note the different scales in (a) and (b).
to describe what may happen to each simple root. Let the condition number of the
root ? of p be defined by
= iip?lpiWia?o su?Y lip?)Id?IIpIId
Then from Prop. 2.1 we can readily derive the following formula:
Proposition 2.2
=			il?lid-i			(4)
ilPlid1y(?)1
(Thus for infinitesimal t, the component of Zciipii?(P; d) about ? is the disk of radius
p; d)). This result is essentially equivalent to formulas obtained by Gautschi in
a sequence of several papers on polynomial condition numbers [2]. We define the
condition number of the zerofinding problem for p to be
= m?x ?(?,p;d) .			(5)
3 Pseudospectra of companion matrices
For each monic polynomial p = ?flj=0 CjZ? E IP, the companion matrix associated with
p is the n x n matrix
o			--H?
1			--Hc1
Ap =
1			Cn?1
5
(6)
The characteristic polynomial of Ap is p itself, and therefore the set of eigenvalues of
Ap coincides with the set of zeros of p. We now introduce concepts for companion
matrices analogous to those just given for polynomials.
For each p E II), we define the &pseudospectrum of Ap by
A?(A?; d) =			Ei C : z E A(A) for some A E Mn with IlA --H Aplid ? eJ.			(7)
The size of A?(A?; d) is related to the conditioning of the companion matrix eigenvalue
problem. The appearance of the norm Ii lid in (7) corresponds to the consideration
of diagonally weighted perturbations in the entries of A?:
ilA --H Apild = iiD(A --H Ap)D?1ll2.
(8)
The main reason we include D in our formulation is that the process of balancing
a matrix, to be discussed in the next section, involves a diagonal similarity trans-
formation; thus this formulation allows us to treat the balanced and unbalanced
cases together. For an eigenvalue algorithm applied to the companion matrix to be
stable, the computed eigenvalues of Ap should lie in a region Ac?(A?; d) for some
C = O(iiAplId).
The following proposition gives an algebraic characterization of A?(A?; d) in terms
of the level curves of the norm of the resolvent of Ap.
Proposition 3.1
A?(A;d) =			E C : li(zI--H A)?1iid >
This result holds for any n x n matrix and any matrix norm induced by a vector norm;
see [1]. For completeness, we include a proof, omitting the subscript d for simplicity.
Proof. It is easily shown that A?(Ap; d) is a subset of the right-hand side. To show
the converse, consider z Ei C such that ii(zI--H A)-1ii = ?0?1 for some ? ?. We can
find u C C? such that
ii(zI --H A)?1uii = ?-? and hull = 1.
Let v = to(zI --H A)?1u; then liv ii = 1. By a standard corollary of the Hahn-Banach
theorem, there exists a linear functional f ? (C?)* (the dual space of C?) such that
f(v) = 1 and ilfil = I. Let a = (f(e1),. . . ,f(e?))*; then f(x) = a*x Vx Ei C? (the
superscript * means conjugate transpose). Consider the matrix A = A + ?0ua*. A
satisfies
Av = Av + ?ua*v = Av + ?u = zv,
showing that z is an eigenvalue of A. Now iiA --H All = licoua*li = ?iifii = E0 < ?, 50
z E A?(A?; d), and this completes the proof.			El
6
x 10'
0			5			10
Figure 2: eiiApii?-pseudospectra (e 1o--H?, 10-6) of the unbalanced and balanced
companion matrices corresponding to the same polynomial as in Figure 1. Note the
different scales in (a) and (b), and the approximate agreement of Figs. 1(b) and 2(b).
0
0			2 xli
Example 3.1. Consider the companion matrices corresponding to the same poly-
nomial p as in Example 2.1. Figure 2 compares e??Ap???-pseudospectra of the unbal-
anced and balanced (see next section) companion matrices, which turn out to be far
larger in the unbalanced case. By contrast, note the approximate agreement of the
curves in Figure 2(b) with those in Figure 1(b).
We define the condition number of the simple eigenvalue A of a matrix B by
(a)			unbalanced			(b) balanced
(9)
where the superscript * denotes conjugate transpose and x and y are left and right
eigenvectors of B, respectively, corresponding to the eigenvalue A. (Just as (4) is
a consequence of Prop. 2.1, (9) can be derived as a consequence of Prop. 3.1.) For
the special case of the companion matrix Ap, the left and right eigenvectors can be
written in closed form:
IA--HAl
?(A,B;d)= IIBJBiWi?o S??Y lIB --H Bild/IlBild
A formula for this quantity is (see [4]):
?(A, B; d) = IBlid IIXIId?1IlYlid
Ix*yl
7
x			=			(1,A,...,An--Hl)T,			(10)
y			=			(bo,bi....,bn?i)T,			(11)
where bo,..., bn?i (dependent on A) are the coefficients of the monic polynomial
(p(z) --Hp(A))/(z--HA) = Z?i=%1biZ? The expression for the condition number then
reduces to the following formula:
Proposition 3.2
?(A, A?; d) = IlApild Ilt(Al)yll(dAllilld?l
where ?A) = (bo,.. .,bn?i)T
(12)
(Thus for infinitesimal e, the component of A?ii4ii?(Ap; d) about A is the disk of
radius etc(A, A?; d)). If the eigenvalues of Ap are simple (equivalently, the zeros of p
are simple), we define the condition number of the problem of finding the eigenvalues
of Ap with respect to perturbations of the form given in (8) to be
= m?x tc(?,A?;d).			(13)
4 The importance of balancing
The conditioning of the eigenvalue problem for a companion matrix Ap may be
changed enormously by a diagonal similarity transformation. The best one could
hope for would be to find such a transformation that makes the eigenvalue problem
no worse conditioned than the underlying polynomial zerofinding problem. Our ex-
periments show that to up to a factor of 10 or so, the similarity transformation known
as balancing tends to achieve just this. In this section we shall give some details, com-
paring balancing with other diagonal similarity transformations and measuring the
results by the scalar condition number K of (13). The plots of the next section carry
these observations further, revealing that balancing tends to achieve an approximate
match not only of condition numbers but also of pseudospectra and pseudozero sets.
Recall from (8) that problems of the condition of Ap can be formulated in two
equivalent ways: we can either leave the 2-norm fixed and change Ap to DApD?1,
where D is a diagonal matrix, or leave Ap fixed and change the 2-norm toll. lid, where
d = diag(D). Using the latter formulation, we consider four possibilities:
1. Pure companion matrix. One possibility is to consider the matrix (8) with
no diagonal similarity transformation. Equivalently, the diagonal similarity trausfor-
mationisdefinedbyct=e= #n(l,...,l)T.
2. Balancing. Balancing the matrix Ap corresponds to finding a diagonal matrix
T ? V such that TApT?1 has the 2-norm of its ith row and ith column approximately
equal for each i = 1,... ,n (we denote ??diag(T-1)??2 diag(T) by t). This idea was
introduced in [15] and is a standard option in EISPACR [17] and also the default
procedure in Matlab [11].
3. Scaling. For each scalar ? > 0, the associated scaled monic polynomial of
p ? IF) is
1			Cj
Pa(Z)=anp(OZ)=Z,			.z.
s=0 an?t
8
The diagonal similarity transformation that corresponds to this scaling operation is
defined by the diagonal matrix Da = diag(d(a)), where d(a) --H II(a" ... ,o1)II2(Q?n,... ,
4. Coefficientwise. If p has all nonzero coefficients, there is a natural diagonal
similarity transformation associated with it, viz., the coefflcientwise diagonal similar-
ity transformation, which is represented by the diagonal matrix C = diag(c), where
C = IIPII2P?1.
We have found empirically that the balancing operation tends to achieve the best
conditioned eigenvalue problem for Ap among these four choices of diagonal similarity
transformations. Specifically, consider the ratios of the eigenvalue condition numbers
for the three other problems to that of the balanced companion matrix, i.e., the ratios
K(Ap; t)
(14)
for d = c, d(a), and e, with t corresponding to the balancing operation as indicated
above. In the case of polynomial scaling, a is chosen to be the optimal value in the
sense that a(A?; d(a)) is minimized. We have found that a(A?; e) and a(A?; d(a)) are
usually much greater than unity, implying that these two eigenvalue problems are
much worse conditioned than the eigenvalue problem associated with the balanced
matrix. The third ratio, a(A?; c), is of order 1 in many cases, but coefflcientwise
transformations have the defect that they are not well defined when some of the
coefficients of p are zero.
Having found that balancing operation achieves approximately the best condi-
tioned eigenvalue problem for Ap, we then compare this matrix condition number
with the condition number of the coefflcientwise perturbed zerofinding problem for p.
That is, we consider the ratio
tc(A?;t)
Our experiments show that this ratio tends to be of order 1.
Table 1 gives empirical evidence to support these statements. It is a tabulation
of the four ratios discussed above for a variety of degree-20 monic polynomials:
(1) "Wilkinson polynomial": zeros 1,2,3,... , 20.
(2) the monic polynomial with zeros equally spaced in the interval [--H2.1,1.9].
(3) p(z) = ?2kO?o zk/k!.
(4) the Bernoulli polynomial of degree 20.
(5) the monic polynomial z20 + z19 + z18 + ... + 1 (zeros are the 21st roots of unity
except 1).
(6) the monic polynomial with zeros 2?1o, 2-?, ....... ,
9
(7) the Chebyshev polynomial of degree 20.
(8) the monic polynomial with zeros equally spaced on a sine curve, viz., (2r/19(k+
0.5)) + i sin(2?/19(k + 0.5)), k = --H1o,--H9,--H8,... `9.
coefficientwise			scaled			(a)			unbalanced tc(A?; t)/?(p; c)
(1)			3.4 x			100			1.9 x			108			(9.0)			2.3 x			10 1			1.4 x			10'
(2)			2.0 x			10			4.2 x			10?			(0.9)			4.5 x			10?			3.6 x			100
(3)			1.4 x			100			1.7 x			10			(10.0)			1.1 x			10?			5.2 x			100
(4)			*			4.6 x			10			(1.6)			2.9 x			106			8.7 x			100
(5)			1.1 x			100			9.4 X			10-'			(1.1)			1.1 x			10			7.6 x			10
(6)			1.7x			101			1.2 x			1024			(0.8)			1.4 x			10 ?			2.1 x			101
(7)			*			6.4 x			10?			(0.6)			6.5 x			10?			4.4 x			100
(8)			*			3.4 x			10?			(1.7)			9.6 x			10			2.6 x			100
* coefficientwise transformation not defined because of zero coefficient.
Table 1. Comparison of the conditioning of the balanced cornpanion matrix elgenvalue
problem with that of the eigenvalue problems of various matrices diagonally similar
to Ap.
The scatter plot of ?(p; c) versus ?(A?; t) in Figure 3 provides further evidence
that the condition of the balanced companion matrix eigenvalue problem is typically
comparable to that of the coefficientwise perturbed zerofinding problem. For this
plot, we considered a random sample of one hundred degree-lO monic polynomials
with coefficients of the form
a, x 10?? + ja2 x 10??			(15)
where a? and e? (i = 1,2) are drawn from the uniform distributions on the intervals
[--H1, 1] and [--H10, 10], respectively. The idea of considering polynomials with random
coefficients of this form to test the quality of zerofinding algorithms was proposed by
Jenkins and Traub in 1974 [7].
We should note that the conditions of the balanced companion matrix eigenvalue
problem and the coefficientwise perturbed zerofinding problem are closer for some
polynomials than others. When the zeros of the polynomials are of widely varying
magnitudes, the match tends to be worse; this accounts for the deviation of some
points from the dashed line in Figure 3.
5 Numerical experiments
We now come to a sequence of more detailed experiments, Figures 4--H11. For the
same eight polynomials described in the last section, each figure presents first of all a
10
10
cond. no. of
balanced
eigenvalue
problem
10			,6
x ???
1? ?66
io5
1010
10			10			10
cond. no. of zerofinding problem
100
10
1020
graphical comparison of two pseudozero sets Z?ii?ii?(p; c) (c = IIpII2p?1) of the polyno-
mial with the corresponding pseudospectra A?IIApIIt (A?; t) of the associated balanced
companion matrix. A reasonably close agreement is observed in all cases.
This agreement of pseudozero sets and pseudospectra suggests that it ought to be
possible to compute zeros of polynomials stably via eigenvalues of balanced companion
matrices. To test this prediction, we have compared three zerofinding methods:
1. J-T. Our first zerofinder is the Jenkins-Traub program CPOLY, available from
ACM TOMS via Netlib and also in the IMSL library [6]. The innermost step of this
algorithm is equivalent to a certain modified generalized Rayleigh quotient iteration
for companion matrices (see Appendix A). A three stage procedure is used, each stage
being characterized by the type of shift involved.
2. M-R. Our second zerofinder is the Madsen-Reid code PAi6 from the Harwell
Library [10], which is a Newton-based method coupled with line search. The algorithm
is designed to find zeros beginning with the zero of smallest magnitude of a given p.
Stage 1 of the method uses Newton formula to find a search direction for the real
function p(z) I. Once the iterates are close enough to a zero of p(z), the algorithm
enters stage 2 and uses the standard Newton iteration.
3. ROOTS. Finally, we consider ROOTS, the Matlab zerofinding code based on
constructing the balanced companion matrix and finding its eigenvalues by standard
methods. (Matlab actually uses the QZ algorithm, but our experiments suggest that
the QR algorithm gives similar results.)
In our experiments, we first find the "exact" roots of p by computing the eigenval-
ues of Ap in quadruple precision ( 34 digits) via a standard EISPACK computational
Figure 3: Comparison of the condition number ?(A?; t) of the balanced companion
matrix eigenvalue problem with the condition number ?(p; c) of the coefflcientwise
perturbed zerofinding problem for 100 degree-lO polynomials with random coefficients
of the form (15).
ii
path for eigenvalues of a general complex matrix, (CBAL, CORTH, COMQR) [17].
The rest of our computations are then carried out in double precision ( 16 digits).
For each of the eight polynomials, Figures 4-11 list the maximum absolute error of
the roots computed by the above algorithms together with the condition numbers of
the coefflcientwise perturbed zerofinding problem for p and the balanced companion
matrix eigenvalue problem for Ap.
Figures 4--H11, and our other numerical experiments, indicate that PA16 is typically
the most accurate of these zerofinding algorithms, followed by ROOTS. Further em-
pirical evidence is presented in Figures 12--H14, which show scatter plots corresponding
to the same sample of one hundred random degree-lO monic polynomials as in Fig.
3. A few remarks are in order concerning these scatter plots:
First, the scatter plots indicate that ROOTS and PA16 are stable. The Jenkins--H
Traub code CPOLY appears somewhat unstable in a few cases.
Second, the error in the zeros computed by ROOTS is typically of the order of mag-
nitude of the condition number of the coefficientwise perturbed zerofinding problem,
not the balanced companion matrix eigenvalue problem. With reference to Figure 3,
we thus observe that ROOTS finds the zeros somewhat more accurately than would
be predicted by the condition number of the balanced companion matrix eigenvalue
problem. We have carried out numerical experiments to explain this phenomenon and
find that it is apparently a consequence of the sparsity of the balanced companion
matrices. The condition number (9) is of course defined for a general square matrix,
not taking sparsity or other structure into account.
Finally, the reader will note the appearance of some points well below the dashed
diagonal line for J-T and M-R, but not ROOTS. These correspond to examples with
zeros that are nearly multiple. Evidently the performance of ROOTS matches lin-
ear perturbation theory in these cases, whereas J-T and M-R actually do somewhat
better. Thus in these cases ROOTS could be said to be mildly unstable.
12
pseudozero sets
-1
0
0
0
0			10			20
error(J-T)
error(M-R)
error(ROOTS)
pseudospeotra
0
0			10			20
=			1.36e--H02
=			1.53e--HO2
=			3.58e--HO3
Figure 4: Wilkinson polynomial (? = 1o?14, 10-13).
pseudozero sets
0
-2			0			2
u?(p;c)			=			2.76e--HOl
=			4.OOe+OO
pseudospectra
0
-2			0			2
error(J-T)			=			5.31e --H 13			u t(p; c)			=			6?57e --H 12
error(MR)			=			2.46e--H 13			u			=			2.34e--H 11
error(ROOTS)			=			6.lle --H 13
Figure 5: Polynomial with zeros equally spaced in [--H2.1, 1.9] (t = 1o-?, 10-2).
13
pseudozero sets
10
0
pseudospeotra
10
0k'
o+10
0
0
10
E
.5
error(J-T)
error(M-R)
error(ROOTS)
0			10
=			1.98e--H11			utc(p;c)			=			3.16e--H11
=			1.OOe--H12			=			1.66e--HlO
=			9.OOe --H 12
Figure 6: Partial sums of the Taylor series for ez (? = io-?, 1o-?).
pseudozero sets
0
pseudospectra
0'o+			0			?o+
-2			0
error(J-T)
error(M-R)
error(ROOTS)
3
-2			0			3
=			4.IOe--H12
=			1.17e--H12
=			1.38e--H12
Figure 7: Bernoulli polynomial (? = io-?, to--H3).
14
u			=			1.lOe--H 11
u'c(A?;t)			=			9.62e--H11
pseudozero sets
2
15
=			5.12e--H13			utc(p;c)			=			6.49e--H12
=			2.84e--H13			=			1.36e--HlO
=			9.66e--H 13
Figure 9: Polynomial with zeros 2?1o, ....... , 2? (e = 1o-?, 10--H2).
0			200			400
ol			0
-2
-2
pseudospectra
-2
-2
0			2			0			2
error(J-T)			=			4.92e--H 13			utc(p;c)			=			4.22e--H 16
error(M-R)			=			2.48e--H 16			u			=			3.19e--H 15
error(ROOTS)			=			1.29e--H 15
Figure 8: Polynomial with zeros equal to the 21st roots of unity except 1 (e =
1o?1.2, 10--H1).
pseudozero sets
200			200
0
-200			-200
0
error(J-T)
error(M.R)
error(ROOTS)
200			400
pseudospectra
pseudozero sets
Th			Th-
-2
-2
4
A
-4
0			2
pseudozero sets
pseudospecira
16
error(J-T)			=			3.7Oe--H 13
error(M-R?)			=			5.63e--H 13
error(ROOTS)			=			5.66e--H 13
Figure 11: Polynomial with zeros equally spaced on the section of the sine curve from
--H7r to 7r			= 1o--H?, 10?2).
0			4
u?(p;c)			=			5.O5e--H12
=			1?9e--H11
error(J-T)
error(M-R)
error(ROOTS)
= 1.98e--H1O
= 8.59e--H 12
= 4.89?--H 12
Figure 10: Ghebyshev polynomial T20(x) (? = io-4, 1o-?).
pseudospectra
.4
0			4
-4
A
2
-2
0			2
utc(p;c)			=			5.55e--H11
uK(Ap;t)			=			2.45e--HlO
2
10
10
10
10
max.
error 10
10
, x
x6-xX
x 6x
6x6Th6
xl
10
10			10			10
machine epsilon times
cond. no. of zerofinding problem
Figure 12: Comparison of the accuracy of the Matlab companion matrix algorithm
ROOTS with the conditioning of the coefflcientwise perturbed zerofinding problem
(100 random problems).
10
10
10
max.
error 10
10
10
x			6,x
x			6? ?
x			?			?x6xX?
?			6
xx			x
I?xx???x?x xx x
x			?
x?			X
6/6Ti
, x
10			10
machine epsilon times
cond. no. of zerofinding problem
Figure 13: Same as Fig. 12 but for the Jenkins-Traub code CPOLY.
17
10
10
10
max,
error			10?
10			4xXx
x			xX X			x
x ?			?
10			1&			10			10
machine epsilon times
cond. no. of zerofinding problem
Figure 14: Same as Fig. 12 but for the Madsen?I?eid code PAl 6.
6 Asymptotic agreement of pseudozero sets and
pseudospectra
In this section we establish various mathematical relationships between pseudozero
sets and pseudospectra of the associated companion matrices, showing in particular
that in the limits e H 0 (infinitesimal perturbations) and ? H 0O (zeros at the origin),
they become identical to one another.
Let us begin by noting that the e-pseudozero set of p, Z?(p; d), is contained in
the e'-pseudospectrum of Ap, A?i(A?; d), if e' = e dn?i I?? This is true because we
may view Z?(p; d) as the set of eigenvalues obtained by structured perturbations of
Ap that preserve the companion matrix form:
Z?(p; d) = C : z ? A(Ap) for some p ? IF) with lAp --H Aplid ? t'1.
For a similar reason, the lemniscatic region L??-i (p) is contained in Z?(p; d), for we
may view this region as the set of zeros obtained by perturbing only the constant
term of p:
= Ei C : z Ei Z(p+ ? for some scalar ? E C with l?dol< e).
In summary, we have the following proposition:
Proposition 6.1 The lemniscatic regions, pseudozero sets and pseudospectra satisfy
L??-i(p) C Z?(p;d) C A???1-1(A?;d)
18
provided d0, dn?1 $ 0.
We showed in Section 2 that the &pseudozero set Z?(p; d) is the region bounded
by the e-level curve of a certain function defined over the complex plane (denoted
by ?(. ; d) later in this section). A similar result (Proposition 6.2) holds for the
e-pseudospectrum A?(Ap; d) when e is sufficiently small (we denote the associated
function, defined below, by ?( ; d)). From these simple characterizations, we are able
to give analytical relationships between:
1.
2.
3.
The pseudozero sets and the psendospectra (limit e 0). Corollary 6.1 shows
that A?(Ap; d(1)) and Z?(p; d(2)) agree with one another in the limit e H
where
CI = e (IIpII?(?)?(Ap; d?1?))/(IIApII??i??(p; d(2))).
(16)
The pseudospectra for different diagonal similarity transformations (limit e H
0). Corollary 6.2 shows that A?(A?; d(1)) and A?i(A?; d(2)) agree with one another
in the limit e H 0? where
CI = e (IIApII???)tc(Ap; d(1)))/(IIApIId(1) tc(A?; d(2))).
(17)
The lemniscatic regions, pseudozero sets and pseudospectra (limit a H oo). For
the special case of polynomial scaling, Corollary 6.3 shows that the lemniscatic
region L?Qn-1 (p), the pseudozero set Z?n-i (p; d(a)) and the pseudospectrum
A?(Ap; d(?)) agree with one another in the limit a H 00.
Lemma 6.1 For each z E ? --H
(zi --H Ap)?1 =
1			B(z)			(18)
p(z)b(z)zT
(19)
zp)lz)h(z?l)iT+ H(z),
where b(z) = (b0(z),... , bn?i(z))T is the vector of coefficients of the polynomial
--H p(z))/(w --H z) = Zt?wo1 bi(z)w? h(z?1) --H (hn?1(z?1),... , ho(z?l))T is the vec-
tor of coefficients of the polynomial (p*(w) --H p*(z--H1))/(w --H z--H1) = Zs%%1 hj(z?1)w?,
and
0 1			... z??2
.1
0
Proof. We can write zi --H Ap as
zi--H Ap = T(z) +penT?
H(z) = z?1			?--H1
1--Hn
z			... z?1
where T(z) is the bidiagonal matrix with z along its main diagonal and --H1 along its
subdiagonal. Application of the Sherman-Morrison formula gives (18), and further
algebraic manipulations yield (19).
19
For each D E D, let the functions r(. ; d), `k( ; d), and ?(. ; d) be defined by
IIb(z)IId			if IzI< 1
r(z;d) =
Iz?1I IIh(z?1)IId if Izi> 1,
_ Ip(41			_ _____
--H IIzIId--H1'			?(z;d)--H r(z;d)
(20)
(21)
Note that r(z;d) is continuous at each z E Z(p) because IIb(z)IId =Iz?1 I IIh(z?1)IId
if p(z) = 0. Also, from Proposition 2.1, Z?(p; d) is the region bounded by the e-level
curve of ?( ; d). Furthermore,
tc(z,p; d) = IIpIIdIIzIId--H1/Ip'(z) I,
?(z, A?; d) = IIApII?(r(z; d)IIzIId--Hl )/ip'(z) I,
(22)
(23)
for all z ?
Proposition 6.2 For each D ?
II(ZI4zA;l?lIId = 1+0			as ?(z;d)?oo.			(24)
Proof. From Lemma 6.1, for z ? 1,
II(zI --H Ap)?1 lid --H IIt(z)IIdIIzIId--Hl < lIB(z)lId ? rn1,
Ip(z)I
wherem1=n?max?Idj/d?I : O<i<k--H1 1<k<n--H1J. ForIzI>1,wehave
IIh(z?1)IIdIIiIId--Hl < II1I(z)IId < ?2,
II(zI --H Ap)?1IId --H			zp(z)
wherem2=n.maxf Idi/dkl : k+1 <i<n--H1 O<k<n--H2J. Combiningthese
estimates gives
which implies (24).
II(zI--HAp)?1IId _ max(mi,m2),
--H1 ?
El
Example 6.1 We may paraphrase Proposition 6.2 as follows: In the limit e 0,
A?(Ap; d) agrees with the region bounded by the e-level curve of ?( ; d). Figure 15
illustrates this agreement pictorially. Consider the monic polynomial p with zeros
1,2,... , 10 and diagonal similarity transformations d?1? = e and d?2) --H c. The figure
20
4
or			0 0
I;
oF
-4			-4
0			5			10
(a) normwise (? --H 10--H10			1o--H?)
0			5			10
(b) coefficientwise (c = io--H? 1o?6)
Figure 15: ?-pseudospectra of Ap (solid) and ?-level curves of?(., d(?)) (dashed) corre-
sponding to normwise and coefflcientwise diagonal similarity transformations for the
degree-lO monic polynomial with zeros 1,2,... , 10. The dashed curves are not visible
because they match the solid ones to plotting accuracy.
shows the &levd curves of (zi--H Ap)?'IId(i) and ?(z;d(?)), i = 1,2. To the resolution
of this figure, these two sets of curves are indistinguishable.
We note that in analyzing the conditioning of the eigenvalue problem of Ap, one
is usually only interested in a region where ?(z; d)  1. In this case, Proposition 6.2
shows that ?(z; d) provides an accurate estimate of II(zI --H Ap)?1 lid. Note that only
0(n) flops are needed to calculate ?(z; d) by llorner's rule as compared with 0(n3)
flops needed to calculate II(zI --H Ap)?1iid by the SVD.
Corollary 6.1 For any D1,D2 ED and z0 E Z(p),
II(zI--HAp)?1iid(1) H IIPIId(2)tc(z0,A?;d(1))
IIApIId(1) tc(z0,p;d(2))			as z H Z0			(25)
Proof. The left-hand side is equal to
II(zJ--HAp)?1IId(1) r(z;d(1))IIzIIl/d(1)
?(z; d(1))			IIZIIl/d(2)
Since ?(z; d) oo as z H Z0, by Proposition 6.2, the first fraction tends to 1 as
z H Z0 The second tends to
r(zo;d(1))iizoiil/?l)
IIZoIIl/d(2)
because r(z; d) is continuous at z = z0. Substituting (22) and (23) into this limit
gives (25).
21
10			20
Figure 16: e-pseudospectra (solid) of the balanced companion matrix of Ap and e'-
pseudozero sets (dashed) of p with respect to coefficientwise perturbations for the
monic polynomial with zeros 1,2,... , 20 (e = 10?12, 10-11).
Example 6.2 We may paraphrase Corollary 6.1 as follows: In the limit e H
the components of A?(Ap; d(1)) and Z?i(p; d(2)) (e' as in (16)) containing z0 agree with
one another. Figure 16 illustrates this agreement for the "Wilkinson polynomial"
p with zeros 1, 2,... , 20. The figure compares the pseudospectra of the balanced
matrix of Ap and the pseudozero sets of p with respect to coefflcientwise perturbation.
Choosing z0 = 15, the figure plots the boundaries of A?(A?; t) and Z?i(p; c) (e' =
e (IIpIIcK(zo,Ap;t))/(IIApIIttc(zo,p;c))) for e = 10?12, 10?11 The agreement is close.
Corollary 6.2 For any D1, D2 c D and z0 ?
II(zI--HAp)?1IId(1) __ IIApIId(2)t:(z0,A?;d(1))
II(zI --H Ap)?1IId(2)			IIApIId(1) K(z0, A?; d(2)) as z			z0.
(26)
Example 6.3 Corollary 6.2 can be paraphrased as follows: In the limit e H 0, the
components of A?(A?; d(1)) and A??(A?; d(2)) (e' as in (17)) containing zo agree with
one another. For illustration, we consider the "Wilkinson polynomial" again as in
Example 6.2. For z0 = 15, Figure 17 plots the boundaries of A?(A?; t) and A?1(A?; c)
for e = 10?12,... , 10?8. Again the agreement is close.
Corollary 6.3 For each scalara >0, let d(a) = (?n,...,?1)2(??n,...,??1) The
following limits hold for each z ?
--HAp)?1IIdip(z) = 1 + 0 as a H 00, (27)
an--Hi
iI'(z;d(?))
Ip(z)I =1+0			as a?+oo.
(28)
22
-10
0			10			20
Figure 17: e-pseudospectra (solid) of the balanced matrix of Ap and e'-pseudospectra
(dashed) of the matrix similar to Ap via the coefficientwise diagonal similarity trans-
formation for the same polynomial as in Figure 16 (c = 1o?12,... ,
Proof. It can be shown that
as a H 00.
Substituting this expression into Proposition 6.2 establishes (27), and (28) follows
from the estimate
IIZiI1/d(Q) = a?			as a H 00.
El
Example 6.4 We may paraphrase Corollary 6.3 as follows: In the limit a H
00, L?n-i (p), Z?n-i (p; d(?)) and A?(A?; d(?)) agree with one another. Figure 18
illustrates this agreement for the monic polynomial p with zeros 1,2, 3,4,5. The figure
shows L?Qn-l(P) , Z??n-i(p;d(?)), and A?(A?;d(?)) for e = 1o--H?,1o--H?, with a = 20.
The figure shows that the agreement is close, especially between the pseudozero sets
and the lemniscatic regions.
7 Conclusions
In this paper we have argued that the pseudozero sets of polynomials tend to match
approximately the pseudospectra of the associated balanced companion matrices.
This correspondence gives a geometric explanation of why it is that eigenvalue algo-
rithms applied to companion matrices appear to compute polynomial zeros stably, an
observation that we have confirmed by experiments comparing the companion matrix
23
-2
0-
-3
c			2			6
Figure 18: ?-pseudospectra of Ap (solid), (?cfl?1 )-pseudozero sets of p (dashed) and
(?a?--H1 )-lemniscatic regions of p (dotted) corresponding to perturhations from the
scaling operation on p (? = io-?, 1O??; a = 20).
code ROOTS (Matlab) with the more traditional zerofinding codes PAi6 (Madsen-
Reid) and CPOLY (Jenkins-Traub).
Although our results become exact in certain limits (Section 6), for most poly-
nomials they are inexact and empirical. We do not claim that they apply to all
polynomials without exception. They are solid enough, however, that a reasonable
degree of confidence in zerofinding via companion matrices seems justified.
If accuracy is not an argument againist the use of ROOTS, what about speed? We
have not discussed the basic fact that whereas matrix eigenvalue algorithms require
0(N3) work, polynomial zerofinding codes such as PAi6 and CPOLY come closer
to 0(N2). For large N this difference will certainly become significant. However,
polynomial zerofinding is a rather unusual problem of scientific computing: one is
not interested in large N! It is very hard to conceive of a genuine application where
it would be a good idea to compute zeros of high-degree polynomials specified by
their coefficients. Thus it is our view that asymptotic complexity is not of decisive
importance for zerofinding algorithms.
Appendix
The generalized Rayleigh quotient iteration (GRQI) associated with an N x N
matrix A with initial vectors g(O), h(0) and scalar so generates two sequences of vectors
?g???Ji>i1 fh(?)Ji>i and a sequence of scalars fsiji>i through the following process:
for? = 1,2..
(A --H 5j1)Tg(i+l) --H 9(i)
(A --H siI)h(?+1) --H
24
Sj+1 = Sj + (i+1)Th(i)
9?+1)Th(t+1)
It has been shown that this iteration converges cubically under mild conditions, with
the sequence tsi?i>i converging to an eigenvalue A of A and ?g???Ji>i and ?h(?)Ji>i
converging respectively to left and right eigenvectors of A associated with A [14]. For
general matrices A, the iteration is carried out numerically by solving two linear sys-
tems, for g(i) and h(?), in each step. llowever if (A --H zI)?1 is known analytically or
certain information concerning the left and right eigenvectors of A is known, this pro-
cess may be simplified. In this Appendix, we show that if A is the companion matrix
Ap associated with a monic polynomial p, then a certain modification of GRQI leads
to the Jenkins-Traub variable-shift iteration for polynomial zeros [6]. In other words,
the Jenkins-Traub iteration can be interpreted as a scheme for taking advantage of
companion matrix structure in a Rayleigh quotient iteration so that the work per step
is reduced from 0(N3) to 0(N). This observation is originally due to Jenkins and
Traub themselves (section 5 of [6]). The following derivation is the reverse of theirs
in that the modified GRQI is derived from the variable-shift iteration rather than the
other way around.
If A is an eigenvalue of Ap, then the vector A = (1,A,... , An--H1)T is a left eigenvec-
tor of Ap associated with A. Therefore, instead of having an independent sequence
(g(?)1j>1 for the left eigenvectors, we may modify the GRQI process for Ap by taking
g(i) to be the vector
9(i) =			= (1, ?i--H1,			, stfl?i1)T
The generalized Rayleigh quotient iteration then reduces to the following:
for? = 1,2,			(29)
(A --H siI)h(?+') --H h???,			(30)
Sj+1 = Sj + Th(i)			(31)
Further simplification is possible because (zi --H Ap)?1 is known analytically:
0 1			... z??2
1			.
(zi --H Ap)?1 =			b(z)zT --H			.
p(z)			.			1
0
where b(z) = (bo(z),... , bn?i(z))T is the vector of coefficients of the polynomial
--H p(z))/(w --H z) = Z2?.-1 b.(z)w?. Applying this formula in (30), we get
--H Hp((??b?(si) --H T1???(si),			(32)
where Hi(z) is the polynomial Hi(z) = ??k=%1 h(;)zk and ?(?)(sj) is the vector of co-
efficients of the polynomial (Hi(z) --H Hi(si))/(z --H si). Note that we have used the
identity 8?Th(i) --H Z?k=o1 hk(?)stk = Hj(sj) in arriving at the above expression.
25
Define Hi+i(z) to be the polynomial H?+i(z) = zTh(i+l). Substituting (32) into
this expression gives
Hi+1(z) =
Hi(si) zTb(s.) --H zTii(i)(si)
p(sj)
--H1
z --H
Hi(z) Hi(Si)p(z)
p(sj)
(33)
(34)
Let Hi+i(z) be the normalized polynomialofH?+1(z), i.e., Hi+i(z)_= (p(s?)/H?(s?)) Hi+i(z).
Since sth(?) --H H?(sj) and 5?Th(?+1) --H Hj+1(sj) = (H?(s?)/p(s?)) Hi+i(z), (31) is equiv-
alent to
5j+1 = 5j + p(sj)			(35)
J'i+1 (5j)
The combination of (34) and (35) is exactly the Jenkins-Traub variable-shift iteration.
References
[1]
[2]
J.L.M. van Dorsselaer, J.F.B.M. Kraaijevanger and M.N. Spijker, Linear stability
analysis in the numerical solution of initial value problems, Acta Numerica 1993
(1993), 199--H237.
W. Gautschi, Questions of numerical condition related to polynomials, in MAA
Studies in Numerical Analysis, 24, G.H. Golub, ed., Math. Assoc. of America,
Washington, DC, 1984.
[3] S.K. Godunov, O.P. Kiriljuk and V.1. Kostin, Spectral portraits of matrices,
Preprint #3, Inst. of Math., Acad. Sci. USSR, Novosibirsk, 1990.
[4] C.ll. Golub and C.F. Van Loan, Matrix Computations, 2nd ed., Johns Hopkins
University Press, Baltimore, MD, 1989.
[5] E. Hille, Analytic Function Theory, Vol. 2, 2nd ed., Chelsea Publishing Company,
New York, 1987
[6]
M.A. Jenkins and J.F. Traub, A three-stage variable-shift iteration for polynomial
zeros and its relation to generalized Rayleigh iteration, Numer. Math. 14 (1970),
252--H263.
[7] M.A. Jenkins and J.F. Traub, Principles for testing polynomial zerofinding pro-
grams, Tech. Report, Department of Computer Science, Carnegie-Mellon Uni-
versity, March 1974.
[8] M.A. Jenkins and J.F. Traub, Algorithm 419--HZeros of a complex polynomial,
Comm. ACM 15,1972, 97--H99.
26
[9] H.J. Landau, On SzegJ's eigenvalue distribution theory and non-Hermitian ker-
nets, J. d'Analyse Math. -28 (1975), 335--H357.
[10]
K. Madsen and J. Reid, Fortran subroutines for finding polynomial zeros, Report
HL.75/1172(C.13), Computer Science and Systems Divisions, A.E.R.E. Harwell,
Oxford, February 1975.
[11] The MathWorks, Inc., MATLAB User's Guide, The MathWorks Inc., 1992.
[12] C.B. Moler, ROOTS--HOf polynomials, that is, MathWorks Newsletter, V. 5, No.
1, (1991), 8--H9.
[13] R.G. Mosier, Root neighborhoods of a polynomial, Math. Comp. 47 (1986), 265--H
273.
[14] A.M. Ostrowski, On the convergence of the Rayleigh quotient iteration for the
computation of characteristic roots and vectors, IV. Generalized Rayleigh quo-
tient for nonlinear elementary divisors, Arch. Rat. Mech. Anal. 3 (1959), 341--H
347.
[15] B.N. Parlett and C. Reinsch, Balancing a matrix for calculation of eigenvalues
and eigenvectors, Numer. Math. 13 (1969), 293--H304.
[16] G. Peters and J.ll. Wilkinson, Practical problems arising in the solution of poly-
nomial equations, J. Inst. Maths. Applics. 8 (1971), 16--H35.
[17] B.T. Smith, et al., Matrix Eigensystem Routines--HEISPACK Guide, Lecture
Notes in Computer Science, V. 6, 2nd ed., Springer-Verlag, 1976.
[18] L.N. Trefethen, Pseudospectra of matrices, in D.F. Griffiths and C.A. Watson,
eds., Numerical Analysis 1991, Longman, 1992.
[19] J.H. Wilkinson, Rounding Errors in Algebraic Processes, Prentice-Hall, Engle-
wood Cliffs, NJ, 1963.
[20] J.H. Wilkinson, The Algebraic Eigenvalue Problem, Clarendon Press, Oxford,
1965.
27
