On ill-conditioned eigenvalues, multiple roots of polynomials, and their accurate computations

Zhonggang Zeng^{∗}

September 11, 1998

Abstract

Algebraic eigenvalues with associated left and right eigenvectors (nearly) orthogonal, and polynomial roots that are multiple, have been known to be sensitive to perturbations in numerical computation and thereby ill-conditioned. By constructing an extended polynomial system, it is proved that, under certain conditions, traditionally defined ill-conditioned eigenvalues and polynomial roots can be remarkably stable and insensitive to perturbations, if they are computed as numerical multiple eigenvalues/roots. Furthermore, using the extended polynomial system which is easy to solve, those numerical multiple eigenvalues can be computed with high forward accuracy, while backward errors are estimated as by-products through auxiliary variables of the extended system.

Keywords: eigenvalue, ill-condition, polynomial, multiple root

AMS subject classifications: 65F15, 65F35, 65H05, 65H17

1 Introduction

It is known in classical numerical linear algebra theories that an algebraic eigenvalue, whose associated left and right eigenvectors are nearly orthogonal, is sensitive to perturbations and thereby ill-conditioned [5, 23]. For ill-conditioned eigenvalues, existing numerical methods generally give answers with large (forward) errors. Similar condition theories and computational difficulties exist for multiple roots of polynomials [3, 22]. In fact, if a matrix A has ill-conditioned eigenvalues, it is a close neighbor of a matrix with multiple eigenvalues [15, 24]. Namely,

∗Department of Mathematics, Northeastern Illinois University, Chicago, IL 60625-4699, email:

zzeng@neiu.edu. Research at MSRI, Berkeley is supported in part by Northeastern Illinois University COR Grant 1-30906.

det(A − λI) = (λ − λ_{∗})^{m}p(λ) + ε(λ), m > 1,

with small perturbation ε(λ) from the polynomial (λ−λ_{∗})^{m}p(λ), which has an m-fold eigenvalue λ_{∗}. The behavior of ill-conditioned eigenvalues can be illustrated by the matrix

A = (a_{ij}) = ( \begin{matrix}1 & & & 0 \\ 1 & 1 \\ & ... & ... \\ 0 & & 1 & 1 \end{matrix} )_{100×100} ,  (1)

A tiny perturbation ε = 10^{−10} at the entry a_{1,100} results in a perturbed characteristic polynomial (λ − 1)^{100} − 10^{−10} and an (forward) error 10^{−0.1} ≈ 0.79 from the 100-fold eigenvalue λ = 1. Although the computational result 1.794328234724 can hardly be considered an approximation to the target eigenvalue 1.0, it is an exact solution to a nearby problem. In other words, for a given matrix A, one can only expect a numerical algorithm to calculate exact eigenvalues of another matrix A + E with an inevitable tiny noise matrix E from inherent round-off errors. If the perturbation E is small enough, the computational results and the algorithm are of high backward accuracy. The large forward error, or the actual distance between the exact eigenvalue and the calculated one, is attributed to the sensitivity of the problem instead of the algorithm. “If the answer is highly sensitive to perturbations, you have probably asked the wrong question” [19], meaning that the eigenproblem is from a poor model with a wrong matrix to begin with.

Polynomial root finding is one of the oldest mathematical problems with a huge volume of literature. More than 3000 articles and books on this topic were listed in [12] and a related website^{1}. Many of them are dealing with multiple roots, such as [9, 13, 14, 17, 18, 25]. There has been an “attainable accuracy” barrier for numerical multiple-root solvers [9, 25]: if the polynomial f(x) = (x − x_{∗})^{m}p(x) and k is the number of correct digits attainable in evaluating f(x), then the multiple root x_{∗} can only be calculated to the accuracy of k/m correct digits. With IEEE standard double precision, one can evaluate function values up to 16 digits, making multiple roots of 16 plus multiplicities not worth trying since not even a single digit from the result can be trusted anyway.

1As this paper in writing, the website is

http://www.elsevier.nl/homepage/sac/cam/mcnamee/index.html

Base on these condition theories, every multiple root of a polynomial, and every eigenvalue of a matrix corresponding to a non-linear elementary divisor are automatically ill-conditioned. It seems that we are giving up too many problems for being “wrong”. In this paper, we shall argue that, λ = 1 is an ill-conditioned eigenvalue of matrix A in (1) only because it is regarded and computed as a simple eigenvalue. It is a remarkably well-conditioned numerical multiple eigenvalue that can be, and should be, computed to high accuracy.

In later sections, we shall prove that ill-conditioned eigenvalues or polynomial multiple roots are not sensitive to perturbations under certain conditions. Furthermore, they can be accurately calculated via solving an extended polynomial system, with auxiliary variables estimating backward errors of the eigenvalues. Our proposed method is tested with defective and derogatory algebraic eigenvalues, i.e. multiple polynomial roots, with multiplicity in hundreds and accurately locates multiple roots to 12 digits and beyond.

The extended polynomial system is also applicable in determining the Jordan canonical structure of a matrix. One of the most difficult issues in numerical computation of Jordan canonical form is to decide whether a cluster of computed eigenvalues is actually a multiple eigenvalue [1, 6, 10]. We provide a practical way for this purpose.

2 Multiple roots of a polynomial

From common definitions, a point x_{∗} ∈ C is an m-fold root (or zero) of the polynomial f if

f(x_{∗}) = f′(x_{∗}) = ⋯ = f^{(m−1)}(x_{∗}) = 0, and f^{(m)}(x_{∗}) ≠ 0

x_{∗} is a simple root if m = 1 and multiple root if m ≥ 2. The integer m is the multiplicity of x_{∗}. Namely, if x_{∗} is an m-fold root of f, we can write f(x) = (x − x_{∗})^{m}g(x − x_{∗}), where g is a lower ordered polynomial with nonzero constant term g(0). Theoretically [20], when m ≥ 2, Newton’s iteration

x_{k+1} = x_{k} − \frac{f(x_{k})}{f′(x_{k})} k = 0, 1, ⋯  (2)

converges linearly and the modified Newton’s iteration

x_{k+1} = x_{k} − m\frac{f(x_{k})}{f′(x_{k})} k = 0, 1, ⋯  (3)

converges quadratically to the m-fold root x_{∗} of f. In numerical computation, it is almost a certainty that neither iteration (2) nor (3) converges near x_{∗} unless the multiplicity m is small or f(x) is explicitly evaluated in the factored form (x − x_{∗})^{m}(b_{0} + b_{1}x + b_{2}x^{2} + ⋯). QR iteration [4, 5, 23] on the companion matrix would also show similar devastating forward errors. To see why Newton’s iterations above fail, let’s consider the floating point evaluations of f(x) = (x − x_{∗})^{m}g(x − x_{∗}) and its derivative at x near x_{∗}

fl[f(x)] = (x − x_{∗})^{m}g̃(x − x_{∗}) + (x − x_{∗})^{m−1}ε_{m−1}(x) + ⋯ + (x − x_{∗})ε_{1}(x) + ε_{0}(x)
fl[f′(x)] = (x − x_{∗})^{m−1}ĝ(x − x_{∗}) + (x − x_{∗})^{m−2}σ_{m−2}(x) + ⋯ + (x − x_{∗})σ_{1}(x) + σ_{0}(x)

where ε_{i}(x), σ_{i}(x), i = 0, 1, ⋯, are unpredictable tiny numbers corresponding to x, as results of rounding errors in evaluations. When x is close to x_{∗}, the computed Newton correction

\frac{fl[f(x)]}{fl[f′(x)]} ≈ \frac{ε_{0}(x)}{σ_{0}(x)}

behaves like a random number generator and has little meaningful relationship to original data.

There are many other numerical methods intended for multiple roots of polynomials [9, 13, 14, 17, 18, 20, 25]. They all suffer similar limitations and apparently can only applied to finding polynomial roots of low multiplicity.

Newton’s iterations (2) and (3) and other multiple-root finders fail in numerical computation of multiple roots because they are essentially applied to compute simple roots, where numerically nonzero f′ near x_{∗} is expected in the denominator. To calculate multiple roots of a polynomial, the fact that f^{(j)}(x_{∗}) = 0, j = 0, 1, ⋯, m−1 must be taken into consideration. Based on this observation, we construct an extended polynomial operator F : C^{m} ⟶ C^{m} associated with the polynomial f : C ⟶ C, such that an m-fold root of f is a simple zero of F. Here, from common definitions, a point z^{∗} = (z_{1}^{∗}, ⋯, z_{n}^{∗}) is a simple zero of F : C^{n} ⟶ C^{n} if F (z^{∗}) = 0 and the

Jacobian DF (z^{∗}) ≡ (\frac{∂f_{i}}{∂z_{j}} (z^{∗})) is a nonsingular n × n matrix.

Theorem 1 A complex number x_{∗} ∈ C is an m-fold root of a polynomial f if and only if the point (0, 0, ⋯, 0, x_{∗})^{⊤} ∈ C^{m} is a simple zero of the extended polynomial operator F : C^{m} ⟶ C^{m} associated with f with auxiliary variables t_{1}, ⋯, t_{m−1}, defined as

F ( t_{1} ; ⋮ ; t_{m−1} ; x ) = ( \frac{f(x)}{0!} − t_{1} \frac{f^{(m)}(x)}{m!} ; \frac{f′(x)}{1!} − t_{2} \frac{f^{(m)}(x)}{m!} ; ⋮ ; \frac{f^{(m−2)}(x)}{(m−2)!} − t_{m−1} \frac{f^{(m)}(x)}{m!} ; f^{(m−1)}(x) ) .  (4)

Proof: Since F ((0, ⋯, 0, x)^{⊤}) = (\frac{f(x)}{0!}, \frac{f′(x)}{1!}, ⋯, \frac{f^{(m−2)}(x)}{(m − 2)!}, f^{(m−1}(x))^{⊤}, the

point (0, 0, ⋯, 0, x_{∗})^{⊤} is a zero of F if and only if

f(x_{∗}) = f′(x_{∗}) = ⋯ = f^{(m−1)}(x_{∗}) = 0.

Moreover, the Jacobian

DF ((0, ⋯, 0, x_{∗})^{⊤}) = diag (\frac{f^{(m)}(x_{∗})}{m!}, ⋯, \frac{f^{(m)}(x_{∗})}{m!}, f^{(m)}(x_{∗}))

is nonsingular if and only if f^{(m)}(x_{∗}) ≠ 0, where diag(⋯) denotes the diagonal matrix with listed entries. QED.

Corollary 1 Let x_{∗} be an m-fold root of f ∈ C^{m+1}(Ω). Then there is a neighborhood N of x_{∗} in Ω such that Newton’s iteration

z_{k+1} = z_{k} − [DF (z_{k})]^{−1}F (z_{k}), k = 0, 1, ⋯,

converges to (0, ⋯, 0, x_{∗})^{⊤} from every z_{0} = (0, 0, ⋯, 0, x_{0})^{⊤} ∈ C^{m−1} × N.

Proof Standard theory of Newton’s iteration. QED.

3 Conditioning of polynomial roots

As mentioned above, a polynomial with an m-fold root x_{∗} can be written as (x − x_{∗})^{m}g(x − x_{∗}) with g(0) ≠ 0. In numerical computation, we expect the backward accuracy. That is, whenever we calculate an approximate root of f, we are evaluating an exact root of f̂(x) = f(x) + h(x) with polynomial h(x) “small” in some sense. If f has an m-fold root x_{∗}, we can further assume that

h(x) = (x − x_{∗})^{m−1}ε_{m−1} + ⋯ + (x − x_{∗})ε_{1} + ε_{0}

without loss of generality because higher order terms of (x−x_{∗}) can be absorbed into (x − x_{∗})^{m}g(x − x_{∗}) while maintaining the multiplicity of x_{∗} as a root of f. In this section, we consider the polynomial

f(x) = (x − x_{∗})^{m}g(x − x_{∗}) + (x − x_{∗})^{m−1}ε_{m−1} + ⋯ + (x − x_{∗})ε_{1} + ε_{0}  (5)

where g is a polynomial with nonzero constant term g(0). We are interested in cases where

|g(0)| ≫ max_{0≤i≤m−1} {|ε_{i}|}

which implies that f(x) is “near” the polynomial (x − x_{∗})^{m}g(x − x_{∗}) that has an m-fold root x_{∗}.

Lemma 1 Let {a_{1}, a_{2}, a_{3}, ⋯, a_{n}} ⊂ C with a_{1} ≠ 0. Then there is a δ > 0 such that for every θ with |θ| < δ|a_{1}|, the polynomial

f(x) = θ + a_{1}x + a_{2}x^{2} + a_{3}x^{2} + ⋯ + a_{n}x^{n}

has a root σ = −\frac{θ}{a_{1}} + O(|θ|^{2})

Proof Let â_{k} = max_{|x|≤1} \{\frac{|f^{(k)}(x)|}{k!}\}, k = 2, 3, ⋯, n. Notice that â_{k}’s are indepen-

dent of θ. Let δ = min \{\frac{1}{2}, (\frac{4}{|a_{1}|} ∑_{k=2}^{n} kâ_{k})^{−1}\}, and σ_{0} = −\frac{θ}{a_{1}} where |θ| < δ|a_{1}|.

From the Taylor expansion of f(x) at σ_{0}, we have

f(x) = a_{00} + a_{01}(x − σ_{0}) + ⋯ + a_{0n}(x − σ_{0})^{n} where a_{0k} = \frac{f^{(k)}(σ_{0})}{k!}, k = 0, 1, ⋯, n

and thus |a_{0k}| < â_{k}, k = 2, 3, ⋯, n, Particularly, because σ_{0} ∈ (−δ, δ), we have

inequalities |σ_{0}| ∑_{k=2}^{n} â_{k} ≤ |σ_{0}| ∑_{k=2}^{n} kâ_{k} < \frac{1}{4}|a_{1}|, implying

|a_{00}| = |(θ + a_{1}σ_{0}) + ∑_{j=2}^{n} a_{j}σ_{0}^{j}| ≤ |σ_{0}^{2}| |∑_{j=2}^{n} a_{j}σ_{0}^{j−2}| ≤ |σ_{0}| (|σ_{0}| ∑_{k=2}^{n} â_{k}) < \frac{1}{4}|σ_{0}||a_{1}|

|a_{01}| = |a_{1} + ∑_{j=2}^{n} ja_{j}σ_{0}^{j−1}| ≥ |a_{1}| − |σ_{0}| ∑_{k=2}^{n} kâ_{k} > \frac{3}{4}|a_{1}|.

3 CONDITIONING OF POLYNOMIAL ROOTS 7

Let σ_{1} = −\frac{a_{00}}{a_{01}} , we have |σ_{1}| < \frac{\frac{1}{4}|σ_{0}||a_{1}|}{\frac{3}{4}|a_{1}|} = \frac{1}{3}|σ_{0}| < δ. Generally, with mathemat-

ical induction, we can write

f(x) = a_{k0} + a_{k1} (x − ∑_{j=0}^{k} σ_{j}) + ⋯ + a_{kn} (x − ∑_{j=0}^{k} σ_{j})^{n} with σ_{k} = −\frac{a_{k−1,0}}{a_{k−1,1}}

where a_{kj} = \frac{1}{j!}f^{(j)} (∑_{l=0}^{k} σ_{l}), j = 0, 1, ⋯, n. Particularly,

|a_{k0}| = |a_{k−1,0} + a_{k−1,1}σ_{k} + ∑_{j=2}^{n} a_{k−1,j}σ_{k}^{j}| ≤ |σ_{k}| (|σ_{k}| ∑_{j=2}^{n} â_{j}) < \frac{1}{4} (\frac{1}{3})^{k} |σ_{k}||a_{1}|
|a_{k1}| = |a_{k−1,1} + ∑_{j=2}^{n} ja_{k−1,j}σ_{k}^{j−1}| ≥ |a_{k−1,1}| − |σ_{k}| ∑_{j=2}^{n} jâ_{j} > |a_{k−1,1}| − \frac{1}{4} (\frac{1}{3})^{k} |a_{1}|
> [1 − \frac{1}{4} ∑_{j=0}^{k} (\frac{1}{3})^{j}] |a_{1}| > \frac{5}{8}|a_{1}|

and

σ_{k+1} = −\frac{a_{k0}}{a_{k1}} , with |σ_{k+1}| < \frac{\frac{1}{4} (\frac{1}{3})^{k} |σ_{k}||a_{1}|}{\frac{5}{8}|a_{1}|} < \frac{1}{3}|σ_{k}|, k = 1, 2, · · · .

Thus σ = ∑_{j=0}^{∞} σ_{j} converges and the inequality |a_{k0}| < \frac{1}{4} (\frac{1}{3})^{k} |σ_{k}||a_{1}| yields

f(σ) = lim_{k→∞} f (∑_{j=0}^{k} σ_{j}) = lim_{k→∞} a_{k0} = 0

Also, σ = σ_{0} + ∑_{j=1}^{∞} σ_{j} = −\frac{θ}{a_{1}} + ∑_{j=1}^{∞} σ_{j}, while

|∑_{j=1}^{∞} σ_{j}| < \frac{3}{2}|σ_{1}| = \frac{3}{2} |\frac{a_{00}}{a_{01}} | < \frac{3}{2} \frac{|σ_{0}^{2}| ∑_{j=2}^{n} jâ_{j}}{\frac{3}{4}|a_{1}|} = O(|θ|^{2}).

QED.

The following theorem establishes the sensitivity of a multiple root of a polynomial as a simple zero of the extended polynomial operator (4).

Theorem 2 Let p(x) be a polynomial with an m − fold root x_{∗}, i.e, p(x) = (x − x_{∗})^{m}g(x − x_{∗}) with g(0) ≠ 0. Let

h(x) = ε_{0} + (x − x_{∗})ε_{1} + (x − x_{∗})^{2}ε_{2} + ⋯

be a polynomial with ε_{j} ∈ C, j = 0, 1, 2, · · ·. Then there is a δ > 0 such that if |ε_{j}| < δ for j = 0, 1, 2, ⋯, the extended polynomial operator (4) associated with the perturbed polynomial f(x) = p(x) + h(x) has a simple zero

( t_{1}^{∗} ; t_{2}^{∗} ; ⋮ ; t_{m−1}^{∗} ; x^{∗} ) = ( 0 ; 0 ; ⋮ ; 0 ; x_{∗} ) + \frac{1}{g(0)} ( ε_{0} ; ε_{1} ; ⋮ ; ε_{m−2} ; \frac{ε_{m−1}}{m} ) + O(ε^{2})

where ε = max_{0≤j≤m−1} {|ε_{j}|}.

Proof Consider

f^{(m−1)}(x) = \frac{(m − 1)!}{0!}ε_{m−1} + \frac{m!}{1!} \frac{g(0)}{0!}(x − x_{∗}) + \frac{(m + 1)!}{2!} \frac{g′(0)}{1!}(x − x_{∗})^{2} + · · · .

From Lemma 1, f^{(m−1)}(x) has a root at

x − x_{∗} = −\frac{(m − 1)!ε_{m−1}}{m!g(0)} + O(ε^{2})

i.e., x^{∗} = x_{∗} − \frac{ε_{m−1}}{mg(0)} + O(ε^{2}). From (4),

t_{j}^{∗} = \frac{m!f^{(j−1)}(x^{∗})}{(j − 1)!f^{(m)}(x^{∗})} = \frac{ε_{j−1} + O(ε^{2})}{g(0) + O(ε)} = \frac{ε_{j−1}}{g(0)} + O(ε^{2}), j = 1, 2, ⋯, m − 1.

The point (t_{1}^{∗}, ⋯, t_{m−1}^{∗}, x^{∗})^{⊤} is a simple zero of the operator F in (4) since the Jacobian

DF ( t_{1}^{∗} ; t_{2}^{∗} ; ⋮ ; t_{m−1}^{∗} ; x^{∗} ) = ( \begin{matrix}\frac{f^{(m)}(x^{∗})}{m!} & & & & \frac{f′(x^{∗})}{0!} − t_{1}^{∗} \frac{f^{(m+1)}(x^{∗})}{m!} \\ & \frac{f^{(m)}(x^{∗})}{m!} & & & \frac{f′′(x^{∗})}{1!} − t_{2}^{∗} \frac{f^{(m+1)}(x^{∗})}{m!} \\ & & ... & & ⋮ \\ & & & \frac{f^{(m)}(x^{∗})}{m!} & \frac{f^{(m−1)}(x^{∗})}{(m−2)!} − t_{1}^{∗} \frac{f^{(m+1)}(x^{∗})}{m!} \\ & & & & f^{(m)}(x^{∗}) \end{matrix} )

with f^{(m)}(x^{∗}) = m![g(0) + ε_{m}] + O(ε^{2}) ≠ 0 for small enough δ. QED.

Let f(x) be a polynomial and x_{∗} one of its roots. The condition of x_{∗} as a

root of f(x) has traditionally been determined [23] by \frac{1}{f′(x_{∗})}. If f(x) is subject to

perturbation h(x) = ε_{0} + (x − x_{∗})ε_{1} + (x − x_{∗})^{2}ε_{2} + ⋯, the corresponding perturbed root x_{∗} + Δx_{∗} of

f(x) + h(x) = ε_{0} + (x − x_{∗}) (f′(x_{∗}) + ε_{1}) + (x − x_{∗})^{2} (\frac{f′′(x_{∗})}{2!} + ε_{2}) + ⋯

satisfies |Δx_{∗}| = \frac{1}{f′(x_{∗})}|ε_{0}| + O(ε^{2}), with ε = max_{i}\{|ε_{i}|\}, according to Lemma 1.

Therefore, if x_{∗} is a multiple root of f, it is automatically ill-conditioned. Moreover, the (exact) simple roots x_{j}^{∗}, j = 1, 2, ⋯, m, of f(x)+h(x) corresponding to the m-fold root x_{∗} of f are [22, pp39-40]

x_{j}^{∗} = x_{∗} + [(−\frac{ε_{0}}{g(0)})^{\frac{1}{m}} ]_{j} + O (ε^{\frac{2}{m}} ) , j = 1, 2, ⋯, m

where [(·)^{1/m}]_{j} denotes the j − th interpretation of the (·)^{1/m}. This error estimate further confirm that computing simple roots of a polynomial around a multiple root is an ill-conditioned problem.

Theorem 2 shows that, however, if we consider an m-fold root of the polynomial, along with m − 1 auxiliary variables t_{1} = t_{2} = ⋯ = t_{m−1} = 0, as a simple root of the extended polynomial operator (4), a multiple root of a polynomial may not be sensitive to perturbation at all. Actually, from Theorem 2, a small perturbation h(x) upon the the polynomial (x − x_{∗})^{m}g(x − x_{∗}) results in a forward error proportional to the magnitude of the perturbation. Therefore the m-fold root x_{∗} is well-conditioned. Moreover, not only x_{∗} can be computed within an error proportional to the magnitude of the perturbation, the actual perturbation, or the backward error, is well measured by the auxiliary variables t_{1}, ⋯, t_{m−1}, via calculating the simple zero of the extended polynomial system (4).

More specifically, if g(0) = \frac{f^{(m)}(x_{∗})}{m!} is not small, the m-fold root x_{∗} is numerically stable and the forward error of the computed root x^{∗} is proportional to

\frac{ε_{m−1}}{m} = \frac{h^{(m−1)}(x_{∗})}{m!}, each auxiliary variable t_{i} has a deviation from 0 proportional

to ε_{i}, i = 1, ⋯, m − 1, all with ratio \frac{1}{g(0)} = \frac{m!}{f^{(m)}(x_{∗})}. This ratio effectively describes the sensitivity of x_{∗} as an m-fold root of f(x). Therefore we can consider x_{∗} well-conditioned if

|\frac{f^{(m)}(x_{∗})}{m!}| ≫ |\frac{f^{(k)}(x_{∗})}{k!}| , k = 1, 2, ⋯, m − 1.

Definition 1 A complex number x_{∗} is a well conditioned numerical m-fold root of a polynomial f if

|\frac{f^{(m)}(x_{∗})}{m!}| ≫ ∥(\frac{f(x_{∗})}{0!}, \frac{f′(x_{∗})}{1!}, ⋯, \frac{f^{(m−1)}(x_{∗})}{(m − 1)!})^{⊤}∥_{∞} .

Otherwise, it is an ill-conditioned numerical m-fold root of f.

A traditionally defined ill-conditioned root is an ill-conditioned simple root by Definition 1. Every numerical multiple root is automatically sensitive to perturbations if it is subject to numerical computation as a simple one, while not necessarily illconditioned in the sense of Theorem 2 and Definition 1.

Theorem 2 establishes the forward error of the approximate solution x^{∗} compared to the numerical multiple root x_{∗} in (5) and the meaning of the auxiliary variables t_{j}’s. We shall claim further that, if (t_{1}^{∗}, ⋯, t_{m−1}^{∗}, x^{∗}) is a simple zero of the extended polynomial operator F in (4) with tiny t_{j}^{∗}’s, then x^{∗} is an exact m-fold root of a polynomial close to f, and we know exactly how far it is from f toward that polynomial. In other words, the computed solution x^{∗} is itself a numerical multiple root of f and is backward accurate.

Theorem 3 Let f(x) be a polynomial with associated polynomial operator F defined as in (4). Let (t_{1}^{∗}, ⋯, t_{m−1}^{∗}, x^{∗})^{⊤} ∈ C^{m} be a simple zero of F. then there is a polynomial ĝ(x − x_{∗}) with ĝ(0) ≠ 0 such that

f(x) = (x − x^{∗})^{m}ĝ(x − x^{∗})  (6)
+ [ĝ(0)t_{m−1}^{∗}] (x − x^{∗})^{m−2} + ⋯ + [ĝ(0)t_{2}^{∗}] (x − x^{∗}) + ĝ(0)t_{1}^{∗}.

Proof From Taylor’s expansion of f(x) at x^{∗},

f(x) = ∑_{k=0}^{n} c_{k}(x − x^{∗})^{k} = (x − x^{∗})^{m} (c_{m} + ∑_{j=1}^{n−m} c_{m+j}(x − x^{∗})^{j})
+c_{m−1}(x − x^{∗})^{m−1} + ⋯ + c_{1}(x − x^{∗}) + c_{0},
with c_{k} = \frac{f^{(k)}(x^{∗})}{k!},  k = 0, 1, ⋯, n

Since (t_{1}^{∗}, ⋯, t_{m−1}^{∗}, x^{∗}) is a simple zero of F in (4), we have

ĝ(x − x_{∗}) = c_{m} + ∑_{j=1}^{n−m} c_{m+j}(x − x^{∗})^{j} with ĝ(0) = c_{m} = \frac{f^{(m)}(x^{∗})}{m!} ≠ 0

c_{m−1} = \frac{f^{(m−1)}(x^{∗})}{(m − 1)!} = 0,

c_{j} = \frac{f^{(j)}(x^{∗})}{j!} = t_{j}^{∗} \frac{f^{(m)}(x^{∗})}{m!} = ĝ(0)t_{j}^{∗}, j = 0, 1, ⋯, m − 2

QED.

The polynomial f in (5) with “small” h is, in its written form, near the polynomial (x − x_{∗})^{m}g(x − x_{∗}) and thereby has a numerical m-fold root x_{∗}. In fact, f is near a family of polynomials (x − x̂)^{m}gx̂(x − x̂) since we can always write

f(x) = (x − x̂)^{m}gx̂(x − x̂) + σ_{m−1}(x̂)(x − x̂)^{m−1} + ⋯ + σ_{1}(x̂)(x − x̂) + σ_{0}(x̂)

for every x̂ ∈ C, with σ_{j}(x̂) = \frac{f^{(j)}(x̂)}{j!}. Finding a simple zero of F in (4) is picking a particular x̂ = x^{∗} such that σ_{m−1}(x̂) = 0. Theorem 2 guarantees that this x̂ = x^{∗} is indeed close to x_{∗} we are looking for, while Theorem 3 ensures that the backward error

f(x) − (x − x^{∗})ĝ(x − x^{∗}) = [ĝ(0)t_{m−1}^{∗}] (x − x^{∗})^{m−2} + ⋯ + [ĝ(0)t_{2}^{∗}] (x − x^{∗}) + ĝ(0)t_{1}^{∗}

is known to be tiny as by-product of solving F (z) = 0.

4 Conditioning of algebraic eigenvalues

Let A be an n × n matrix. Its eigenvalues are roots of the characteristic polynomial det(A − λI) of A. For an eigenvalue λ, its associated (right) eigenvector x and the left eigenvector y are vectors in C^{n} such that

(A − λI)x = 0, and y^{H}(A − λI) = 0

respectively, where y^{H} is the conjugate transpose of y. Assuming further that x and y are normalized, the condition number of the eigenvalue λ is defined as [5, 23]

cond(λ, A) = \frac{1}{|y^{H}x|}, ∥x∥_{2} = ∥y∥_{2} = 1.

When cond(λ, A) ≫ 1, the eigenvalue λ has been called ill-conditioned, and A is close to a matrix Â with λ as an exact multiple eigenvalue [15, 24]. To be more specific, Wilkinson [24] proved that for ε = |y^{H}x|, there is a matrix E such that

∥E∥_{2} ≤ \frac{ε∥A∥_{2}}{1 − ε^{2}}

and λ is an exact multiple eigenvalue of A + E. On the other hand, if A has an multiple eigenvalue with ascent (i.e., the size of the largest Jordan block associated with the multiple eigenvalue) l, then for ∥ΔA∥_{2} small enough, given any eigenvalue λ′ of A + ΔA, there is an eigenvalue λ such that [2, 3]

|λ′ − λ| ≤ 2(∥V ∥_{2}∥V^{−1}∥_{2}∥X_{∗}∥_{2})^{\frac{1}{l}} ∥ΔA∥2^{\frac{1}{l}}

where V is a Jordan basis of A and X_{∗} is a basis of the left invariant subspace. For very ill-conditioned eigenvalues, traditional algorithms like Francis’ QR [4, 23] generally give answers that make no sense in terms of forward accuracy. We have been satisfied with the explanation that the calculated eigenvalues, no matter how ridiculously inaccurate, are all simple ones of a neighboring matrix. However, based on our analysis, those “ill-conditioned” eigenvalues may be very much well conditioned in a deeper sense. Given the fact that we have to evaluate the eigenvalues of a neighboring matrix, it is apparently more meaningful to evaluate the eigenvalues of that neighboring matrix with exact and well-conditioned multiple eigenvalues.

Through an orthogonal transformation P, the matrix A is similar to an upper-Hessenberg matrix

P^{H}AP = B = (b_{jk}) = ( \begin{matrix}∗ & ∗ & ∗ & … & … & ∗ \\ b_{21} & ∗ & ∗ & … & … & ∗ \\ & b_{32} & ∗ & … & … & ∗ \\ & & ... & ... & … & ⋮ \\ & 0 & & b_{n−1,n−2} & ∗ & ∗ \\ & & & & b_{n,n−1} & ∗ \end{matrix} ) , P^{H}P = I.

Without loss of generality, we may assume that B is irreducible, namely, the subdiagonal entries b_{j,j−1} ≠ 0, j = 2, 3, ⋯, n. In fact, if anyone of those subdiagonal entries is zero or tiny, the eigenproblem of A can be reduced to those of the two submatrices. For irreducible upper-Hessenberg matrices, Hyman’s method [7, 23], which evaluates the characteristic polynomial of A and its derivatives, can be described as follows [11]: for every λ ∈ C, let \frac{f^{(0)}(λ)}{0!} ≡ f(λ) ∈ C and x_{0} = (x_{01}, ⋯, x_{0n})^{⊤} ∈ C^{n} satisfy

( \begin{matrix}−1 ; 0 ; ⋮ ; 0 & B − λI \\ 0 & 0 ⋯ 0 1 \end{matrix} ) ( \frac{f^{(0)}(λ)}{0!} ; x_{01} ; ⋮ ; x_{0n} ) = ( 0 ⋮ 0 1 ) .  (7)

Then

f(λ) ≡ \frac{det[B − λI]}{∏_{k=1}^{n−1}(−b_{k+1,k})} = \frac{det[A − λI]}{∏_{k=1}^{n−1}(−b_{k+1,k})}, for λ ∈ C,

Differentiating both sides of (7) with respect to λ, we have that \frac{f′(λ)}{1!} and x_{1} = (x_{11}, ⋯, x_{1n})^{⊤} satisfy

( \begin{matrix}−1 ; 0 ; ⋮ ; 0 & B − λI \\ 0 & 0 ⋯ 0 1 \end{matrix} ) ( \frac{f′(λ)}{1!} ; x_{11} ; ⋮ ; x_{1n} ) = ( x_{01} ; ⋮ ; x_{0n} ; 0 ) .  (8)

Recursively, after \frac{f^{(k−1)}(λ)}{(k − 1)!} and x_{k−1} = (x_{k−1,1}, ⋯, x_{k−1,n})^{⊤} are calculated, the next

\frac{f^{(k)}(λ)}{k!} and x_{k} = (x_{k1}, ⋯, x_{kn})^{⊤} satisfy

( \begin{matrix}−1 ; 0 ; ⋮ ; 0 & B − λI \\ 0 & 0 ⋯ 0 1 \end{matrix} ) ( \frac{f^{(k)}(λ)}{k!} ; x_{k1} ; ⋮ ; x_{kn} ) = ( x_{k−1,1} ; ⋮ ; x_{k−1,n} ; 0 ) .  (9)

The equation (7) can be written as

(B − λI)x_{0} = f(λ)e_{1}, x_{0n} = 1

with e_{1} = (1, 0, ⋯, 0)^{⊤}. If λ is an eigenvalue of B (or A) and y the associated left eigenvector, then

0 = y^{H}(B − λI)x_{0} = (y^{H}e_{1}) f(λ) = \overline{y}_{1}f(λ), where y_{1} = e_{1}^{H}y since B is irreducible upper-Hessenberg, it is easy to verify that y_{1} ≠ 0 and therefore f(λ) = 0. Thus x_{0} is an eigenvector of B associated with λ. Similarly, equation (8) can be written as

(B − λI)x_{1} = x_{0} + f′(λ)e_{1}, x_{1n} = 0,

that yields

0 = y^{H}(B − λI)x_{1} = y^{H}x_{0} + \overline{y}_{1}f′(λ).

Therefore the condition number of λ

\frac{|y^{H}x_{0}|}{∥y∥_{2}∥x_{0}∥_{2}} = \frac{|\overline{y}_{1}|}{∥y∥_{2}∥x_{0}∥_{2}} \frac{|f′(λ)|}{1!}

is a multiple of the first derivative of the characteristic polynomial det(A − λI). Generally, (9) can be written as

(B − λI)x_{k} = x_{k−1} + \frac{f^{(k)}(λ)}{k!}e_{1}

and thus

\frac{|y^{H}x_{k−1}|}{∥y∥_{2}∥x_{0}∥_{2}} = \frac{|\overline{y}_{1}|}{∥y∥_{2}∥x_{0}∥_{2}} \frac{|f^{(k)}(λ)|}{k!} k = 1, 2, ⋯, n

Based on our analysis in §3, A small |y^{H}x_{0}| does not necessarily mean the eigenvalue is sensitive to perturbations. The numerical m-fold eigenvalue λ of A is well conditioned if

|y^{H}x_{m−1}| ≫ |y^{H}x_{k−1}| , k = 1, 2, ⋯, m − 1.

5 Numerical experiments

In this section, let J_{k}(λ) denote the k × k Jordan block with eigenvalue λ,

J_{k}(λ) = ( \begin{matrix}λ & 0 & ⋯ & ⋯ & 0 \\ 1 & λ & 0 & & ⋮ \\ 0 & 1 & ... & ... & ⋮ \\ ⋮ & ... & ... & ... & 0 \\ 0 & ⋯ & 0 & 1 & λ \end{matrix} )_{k×k}

and R_{m} the m × m random matrix with entries uniformly distributed on [−1, 1]. All computations in this section are conducted with Matlab. The condition numbers of eigenvalues are calculated with the Matlab tool eigsens.m in [8].

5.1 Example 1: A small matrix

We first test the Newton’s iteration on the extended operator (4) with a 10 × 10 matrix A = P ∗diag (J_{8}(2), R_{2}) P^{⊤}, where P is an orthogonal matrix. The eigenvalues calculated by QR and their condition numbers are

                eigenvalues by QR                          Condition numbers
                0.10174440783530 + 0.33490385222369i                   1.89
                0.10174440783530 - 0.33490385222369i                   1.89
                2.01130354558921 + 0.00468662615444i       3045676685097.
                2.01130354558921 - 0.00468662615444i       3046164799003.
                2.00467135675756 + 0.01130329121656i       3055064826680.
                2.00467135675756 - 0.01130329121656i       3055164380992.
                1.99531399688166 + 0.01128863222830i       3067021306586.
                1.99531399688166 - 0.01128863222830i       3067131909543.
                1.98871110077157 + 0.00467199242653i       3074701046616.
                1.98871110077157 - 0.00467199242653i       3074634066591.

The bottom eight eigenvalues have condition numbers about 3 × 10^{12}. They “exploded” from the 8-fold eigenvalue λ = 2 when calculated as simple ones. The errors are as big as 0.011. However, the Newton’s iteration on the extended polynomial operator (4) has no trouble locating the multiple eigenvalue accurately up to 14 digits in 5 iterations. The progress of the iteration is listed below with meaningful digits in λ_{i}’s and the auxiliary variables

\begin{aligned}λ_{0} & = & 2.1 \\ λ_{1} & = & 2.03 \\ λ_{2} & = & 2.003 \\ λ_{3} & = & 2.00003 \\ λ_{4} & = & 2.000000005 \\ λ^{∗} = λ_{5} & = & 2.000000000000000 \end{aligned}  ( \begin{matrix}t_{1}^{∗} \\ t_{2}^{∗} \\ t_{3}^{∗} \\ t_{4}^{∗} \\ t_{5}^{∗} \\ t_{6}^{∗} \\ t_{7}^{∗} \end{matrix} ) = 10^{−15} ( \begin{matrix}0.3831232126338 \\ 0.5108309501784 \\ 0.5746848189507 \\ −0.3831232126338 \\ −0.1277077375446 \\ −0.7343194908815 \\ −3.1926934386151 \end{matrix} ) .

The next test is to perturb matrix A to

B = A + 10^{−6}E

where E is a random matrix with entries uniformly distributed on [−1, 1]. Then “approximate” eigenvalues of B calculated by QR are

                eigenvalues of B by QR                     condition numbers
                0.10174382731821 + 0.33490379864001i           1.887
                0.10174382731821 - 0.33490379864001i           1.887
                2.15380211678568 + 0.06628105515802i       36054.
                2.15380211678568 - 0.06628105515802i       36054.
                2.05686339846006 + 0.15411717201049i       40732.
                2.05686339846006 - 0.15411717201049i       40732.
                1.93336955871027 + 0.14380540950961i       50109.
                1.93336955871027 - 0.14380540950961i       50109.
                1.85596567090505 + 0.05701606186128i       57472.
                1.85596567090505 - 0.05701606186128i       57472.

That is, the perturbation of magnitude 10^{−6} on the matrix causes QR’s results deviate as large as 0.15. However, by computing them as an 8-fold eigenvalue, the iteration and the auxiliary variables are:

\begin{aligned}λ_{0} & = & 2.1 \\ λ_{1} & = & 2.03 \\ λ_{2} & = & 2.003 \\ λ_{3} & = & 2.00003 \\ λ_{4} & = & 1.9999999 \\ λ_{5} & = & 1.999999870 \\ λ^{∗} = λ_{6} & = & 1.99999986978124 \end{aligned}  ( \begin{matrix}t_{1}^{∗} \\ t_{2}^{∗} \\ t_{3}^{∗} \\ t_{4}^{∗} \\ t_{5}^{∗} \\ t_{6}^{∗} \\ t_{7}^{∗} \end{matrix} ) = 10^{−6} ( \begin{matrix}0.4562499978658 \\ 1.4372014231023 \\ 1.7139977872110 \\ −0.4530018548542 \\ −2.1762844787344 \\ 0.1969728713613 \\ 3.3943752282809 \end{matrix} ) .

The numerical multiple eigenvalue has a forward error 1.3 × 10^{−7} as the result of the perturbation E. More importantly, the auxiliary variables t_{i}^{∗}’s indicate that the eigenvalue λ^{∗} = 1.9999927407264 as an 8-fold eigenvalue has backward error in the magnitude of 3.4 × 10^{−6}, which is a close overestimate of the actual magnitude of the backward error.

5.2 Example 2: Defective eigenvalues with large multiplicities

Let A = diag(J_{50}(2.0), J_{30}(−3.0), R_{20}). That is, A has a 50-fold defective eigenvalue λ_{1} = 2.0, a 30-fold defective eigenvalue λ_{2} = −3.0, and 20 well-conditioned simple eigenvalues. The eigenvalue distribution calculated by QR is plotted in Fig 1. The two multiple eigenvalues explode under QR into eighty “ill-conditioned” eigenvalues and form two circles with a radius about 0.5.

We first scramble the matrix to P AP^{⊤} with an orthogonal matrix P. Thus λ_{1} and λ_{2} are in fact very stable numerical multiple eigenvalues. The iterations and magnitudes of auxiliary variables are

5 NUMERICAL EXPERIMENTS 17

\begin{aligned}λ_{1}^{(0)} & = & 2.001 \\ λ_{1}^{(1)} & = & 2.0002 \\ λ_{1}^{(2)} & = & 2.00001 \\ λ_{1}^{(3)} & = & 2.00000003 \\ λ_{1}^{∗} = λ_{1}^{(4)} & = & 2.0000000000004 \\ ∥(t_{1,1}^{∗}, ⋯, t_{1,49}^{∗})^{⊤}∥_{∞} & = & 7.6 × 10^{−12} \end{aligned} \begin{aligned}λ_{2}^{(0)} & = & −2.999 \\ λ_{2}^{(1)} & = & −3.0004 \\ λ_{2}^{(2)} & = & −3.00003 \\ λ_{2}^{(3)} & = & −3.0000002 \\ λ_{2}^{∗} = λ_{2}^{(4)} & = & −2.999999999997 \\ ∥(t_{2,1}^{∗}, ⋯, t_{2,29}^{∗})^{⊤}∥_{∞} & = & 2.7 × 10^{−11} \end{aligned}

Figure 1: The eigenvalue distribution of a 100 × 100 matrix calculated by QR

Both eigenvalues are calculated with 13-digit forward accuracy and backward error in the range of 10^{−11}. We then scramble the matrix A to BAB^{−1} using a random matrix B = R_{100} with condition number cond(B) = 201.8. Roughly speaking, both the forward and the backward errors are magnified by the condition number.:

\begin{aligned}λ_{1}^{(0)} & = & 2.0001 \\ λ_{1}^{(1)} & = & 2.000003 \\ λ_{1}^{(2)} & = & 2.000000003 \\ λ_{1}^{∗} = λ_{1}^{(3)} & = & 2.00000000005306 \\ ∥(t_{1,1}^{∗}, ⋯, t_{1,49}^{∗})^{⊤}∥_{∞} & = & 4.0 × 10^{−9} \end{aligned} \begin{aligned}λ_{2}^{(0)} & = & −2.999 \\ λ_{2}^{(1)} & = & −3.0004 \\ λ_{2}^{(2)} & = & −3.00004 \\ λ_{2}^{(3)} & = & −3.0000003 \\ λ_{2}^{(4)} & = & −2.9999999995 \\ λ_{2}^{∗} = λ_{2}^{(5)} & = & −2.999999999553 \\ ∥(t_{2,1}^{∗}, ⋯, t_{2,29}^{∗})^{⊤}∥_{∞} & = & 2.9 × 10^{−8} \end{aligned} .

5.3 Example 3: A derogatory eigenvalue of large multiplicity

Let A = diag(J_{10}(2.0), J_{20}(2.0), J_{30}(2.0), J_{40}(2.0), R_{20}). That is, A has a defective and derogatory multiple eigenvalue λ_{1} = 2.0. The matrix A is then scrambled to BAB^{−1} using a random 120 × 120 matrix B = R_{120}. The 120 eigenvalues calculated by QR are plotted in Fig. 2. The 100 “ill-conditioned” QR-eigenvalues are plotted in Fig. 3.

Figure 2: The eigenvalue distribution of a 120 × 120 matrix calculated by QR

Figure 3: The 100 QR-eigenvalues exploded from the 100-fold eigenvalue λ = 2.0

The iteration for the numerical multiple eigenvalue and the magnitude of auxiliary variables are:

\begin{aligned}λ_{1}^{(0)} & = & 2.0001 \\ λ_{1}^{(1)} & = & 2.000004 \\ λ_{1}^{(2)} & = & 2.000000005 \\ λ_{1}^{∗} = λ_{1}^{(3)} & = & 2.0000000000003 \end{aligned}  ∥( t_{1}^{∗} ; ⋮ ; t_{99}^{∗} )∥_{∞} = 2.0 × 10^{−10}.

We can further confirm this accuracy by calculating the the singular values of the matrix BAB^{−1} − 2.0000000000003I_{120}. The smallest 6 singular values are

· · · · · ·
0.03243152758926348
0.01804602808483129
0.00000000000000434
0.00000000000000125
0.00000000000000051
0.00000000000000021.

The singular values clearly indicate that λ_{1} = 2.0000000000003 is a derogatory eigenvalue corresponding to 4 Jordan blocks.

5.4 Example 4: Computing Jordan a canonical form

The following matrix is a classical example for testing computations of Jordan canonical form [10]

A = ( \begin{matrix}1 & 1 & 1 & −2 & 1 & −1 & 2 & −2 & 4 & −3 \\ −1 & 2 & 3 & −4 & 2 & −2 & 4 & −4 & 8 & −6 \\ −1 & 0 & 5 & −5 & 3 & −3 & 6 & −6 & 12 & −9 \\ −1 & 0 & 3 & −4 & 4 & −4 & 8 & −8 & 16 & −12 \\ −1 & 0 & 3 & −6 & 5 & −4 & 10 & −10 & 20 & −15 \\ −1 & 0 & 3 & −6 & 2 & −2 & 12 & −12 & 24 & −18 \\ −1 & 0 & 3 & −6 & 2 & −5 & 15 & −13 & 28 & −21 \\ −1 & 0 & 3 & −6 & 2 & −5 & 12 & −11 & 32 & −24 \\ −1 & 0 & 3 & −6 & 2 & −5 & 12 & −14 & 37 & −26 \\ −1 & 0 & 3 & −6 & 2 & −5 & 12 & −14 & 36 & −25 \end{matrix} ) .  (10)

The Hessenberg reduction results in a numerically reducible upper-Hessenberg matrix with diagonal blocks H_{1} of 2 × 2, H_{2} of 4 × 4 and H_{3} of 4 × 4. Their approximate eigenvalues calculated by QR and condition numbers are,

   Eigenvalues by QR                       Condition numbers
  1.99999996667999                                 5001999.
  2.00000003332001                                 5001999.
  3.00000000000000 + 0.00000001937343i            13650592.
  3.00000000000000 - 0.00000001937343i            13650592.
  2.00000000000000 + 0.00000004969745i             1058412.
  2.00000000000000 - 0.00000004969745i             1058412.
  3.00000010568030                                 1161073.
  2.99999989431971                                 1161073.
  2.00000000000000                                       1.47
  1.00000000000000                                       1.26

Checking the derivatives f′ and f′′ of the characteristic polynomial f at all eight QR-eigenvalues above with large condition numbers, we observed a gap of 10^{7} in magnitudes, indicating they are all numerical double eigenvalues of the corresponding diagonal blocks. Using those QR-eigenvalues as starting iterates with multiplicity 2, we recalculate the eigenvalues with their backward errors, grouped as eigenvalues of H_{1}, H_{2} and H_{3} respectively:

Refined eigenvalues                      Backward errors
(repeated according to multiplicity)
2.00000000000000                           9.99e-016
2.00000000000000                           9.99e-016
3.00000000000000 + 0.00000000000000i       1.44e-015
3.00000000000000 - 0.00000000000000i       1.44e-015
2.00000000000000 + 0.00000000000000i       7.47e-015
2.00000000000000 - 0.00000000000000i       7.47e-015
3.00000000000001                           2.61e-014
3.00000000000001                           2.61e-014
2.00000000000000                           3.12e-015
1.00000000000000                           6.88e-015.

Therefore, it is clear that the matrix A in (10) has a numerical 5-fold eigenvalue 2.0, a numerical 4-fold eigenvalue 3.0, and a simple eigenvalue 1.0. The Jordan structure of A can be determined by the singular values of (A − λI)^{k}’s:

Singular values of A − 2I · · · · · · 0.237920273240 0.096366416213 0.000000000000 0.000000000000  Singular values of (A − 2I)^{2} · · · · · · 0.401358847554 0.067192074052 0.000000000000 0.000000000000 0.000000000000 0.000000000000  Singular values of A − 3I · · · · · · 0.441640607758 0.153907600163 0.000000000000 0.000000000000  Singular values of (A − 3I)^{2} · · · · · · 0.581708725435 0.072422782669 0.000000000000 0.000000000000 0.000000000000 0.000000000000 .

Thus, eigenvalues 2 and 3 are both from two Jordan blocks of sizes at least 2. Therefore, A has the Jordan canonical form

J = diag (J_{3}(2.0), J_{2}(2.0), J_{2}(3.0), J_{2}(3.0), J_{1}(1.0))

Even if the matrix A is perturbed with a random matrix E with ∥E∥_{2} ≈ 10^{−8}, the above process still determines the numerical Jordan canonical structure easily.

5.5 Example 5: Ill-conditioned numerical multiple eigenvalues

The test matrix is known as the Frank matrix [8]. Its entries are

A = (a_{ij}) where a_{ij} = \{ \begin{matrix}n + 1 − max\{i, j\}, & j ≥ i − 1 \\ 0, & j < i − 1 \end{matrix} .

Eigenvalues of the Frank matrix with small magnitudes are ill-conditioned. For example, for a 12 × 12 Frank matrix A, the eigenvalues and corresponding condition numbers are

QR-eigenvalues           Condition numbers
32.22889150157214               3.29
20.19898864587707               4.98
12.31107740086851               3.14
6.96153308556711               1.71
3.51185594858070               6.92
1.55398870913212             216.
0.64350531891624           14467.
0.28474972647719          560310.

                      0.14364644045883         6701452.
                      0.08122798609205        26645187.
                      0.04950694629965        38773986.
                      0.03102829015834        18284203..

The smallest five or six eigenvalues can be considered ill-conditioned in traditional sense. We apply Newton’s iteration on the extended polynomial operator (4) with initial iterate λ = 0 and various attempted multiplicity m. The following table lists obtained numerical multiple eigenvalue corresponding to different attempted multiplicity. The backward error suggested that none of the result is backward accurate:

               attempted         numerical multiple         backward
              multiplicity          eigenvalue              errors
                 m = 2           0.03861651128815           0.000076
                 m = 3           0.05031683035298           0.0052
                 m = 4           0.06988614253113           0.0029
                 m = 5           0.10596312822852           0.016
                 m = 6           0.17866718710829           0.097.

It is clear that those eigenvalues are inherently ill-conditioned, either under traditional definition or as numerical multiple eigenvalues.

6 Concluding remarks

In this paper we discussed the computations of numerical multiple eigenvalues used to be considered ill-conditioned. They can be calculated accurately in forward sense, provided that the multiplicity is known and a good initial iterate is given. For eigenvalues with low multiplicity, eigenvalues calculated with any traditional method such as QR can be used as initial iterates. The multiplicity can be obtained by observing

the magnitude gaps between adjacent derivatives \frac{f(λ)}{0!}, \frac{f′(λ)}{1!}, \frac{f′′(λ)}{2!} · · ·. There are also other methods that can be used to estimate the multiplicity of a multiple root [11, 20]. If the multiplicity is high, eigenvalues calculated by QR could be too inaccurate to be used as initial iterates. The multiplicity could be difficult to determine. On the other hand, if a group of QR-eigenvalues exploded from a numerical multiple eigenvalues, computational experience shows that the mean of the eigenvalue group is usually very good approximate of the numerical multiple eigenvalue and should be used as the initial iterate. Actually, this strategy is often used in numerical computation of Jordan canonical forms [10].

The concept of numerical multiple eigenvalue can be traced back to earlier literatures such as [10, 16, 21], where the stability of numerical multiple eigenvalues has been suggested. It has long been known that the mean of perturbed eigenvalues (or polynomial roots) corresponding to the multiple one is stable [3, 22]. The analysis and numerical results in this paper further confirm that, among the traditionally defined ill-conditioned eigenvalues, some are indeed well-conditioned and can be computed with remarkable forward accuracy. Moreover, as a by-product of the computation, the backward errors can be estimated via the auxiliary variables. If the data matrix warrants such stability and accuracy, there seems little reason to settle on the results with huge forward errors.

Acknowledgments The author wishes to thank Jim Demmel for a number of fruitful discussions on the subject and related issues.

References

[1] T. Beelen and P. Van Dooren, Computational aspects of the Jordan canonical form, in Reliable Numerical Computation, M.G. Cox and S. Hammerling, eds., Clarendon Press, Oxford, 1990, pp57-72.

[2] F. Chatelin, Ill conditioned eigenproblems, in Large Scale Eigenvalue Problems, J. Cullum and R. A. Willoughby, eds, Elsevier Science Publishers B. V. North-Holland, Amsterdam, 1986.

[3] F. Chaitin-Chatelin and V. Frayssé, Lectures on Finite Precision Computations, SIAM, Philadelphia, 1996.

[4] J. G. F. Francis, The QR transformation, Parts I and II, Computer J. Vol. 4 (1961, 1962), pp 265–271, pp 332–345.

[5] G. H. Golub and C. F. Van Loan, Matrix Computations, 3rd ed., The Johns Hopkins University Press, Baltimore and London, 1996.

[6] G. H. Golub and J. H. Wilkinson, Ill-conditioned eigensystems and the computation of the Jordan canonical form, SIAM Review, Vol. 18 (1976), pp578-619.

[7] M. A. Hyman, Eigenvalues and eigenvectors of general matrices, Twelfth National Meeting A.C.M., Houston, Texas (1957).

[8] N. J. Higham, The test matrix toolbox for Matlab, Numerical Analysis Report No. 237, Department of Mathematics, University of Manchester (1993).

[9] M. Igarashi and T. Ypma, Relationships between order and efficiency of a class of methods for multiple zeros of polynomials, J. Comput. Appl. Math., Vol 60 (1995), pp101–113.

[10] B. Kågström and A. Ruhe, An algorithm for numerical computation of the Jordan normal form of a complex matrix, ACM Trans. Math. Software, Vol. 6 (1980), pp 398-419.

[11] T. Y. Li and Z. Zeng, The homotopy continuation algorithm for the real nonsymmetric eigenproblem, further development and implementation, to appear: SIAM J. Sci. Comput.

[12] J. McNamee, A bibliography on roots of polynomials, J. Comput. Appl. Math., Vol. 47 (1993), pp 391–394.

[13] T. Miyakoda, Iterative methods for multiple zeros of a polynomial by clustering, J. Comput. Appl. Math., Vol. 28 (1989), pp 315–326.

[14] H. M. Möller and H. J. Stetter, Multivariate polynomial equations with multiple zeros solved by matrix eigenproblems, Numer. Math., Vol. 70 (1995), pp 311-329.

[15] A. Ruhe, Properties of a matrix with a very ill-conditioned eigenproblem, Numer. Math., Vol. 15 (1970), pp 57–60.

[16] A. Ruhe, An algorithm for numerical determination of the structure of a general matrix, BIT, Vol. 10 (1970), pp196-216.

[17] J. A. Stolan, An improved Šiljak’s algorithm for solving polynomial equations converges qadratically to multiple zeros, J. Comput. Appl. Math., Vol 64 (1995), pp 247–268.

[18] M. V. Straeten and H. Van de Vel, Multiple root-finding methods, J. Comput. Appl. Math., Vol. 40 (1992), pp 105-114.

[19] L. N. Trefethen, Maxims about numerical mathematics, computers, science and life, SIAM News, Vol. 31, No. 1, p 4, 1998.

[20] C. W. Ueberhuber, Numerical Computation 2, Springer-Verlag, Berlin, Heidelberg, New York, 1997.

[21] J.M. Varah,Computing invariant subspaces of a general matrix when the eigensystem is poorly conditioned, Math. Comp., Vol. 24 (1970), pp137-149.

[22] J. H. Wilkinson, Rounding Errors in Algebraic Processes, Prentice-Hall, Englewood Cliffs, N.J. , 1963.

[23] J. H. Wilkinson, The Algebraic Eigenvalue Problem, Oxford University Press, New York, 1965.

[24] J. H. Wilkinson, Note on matrices with a very ill-conditioned eigenproblem, Numer. Math. Vol. 19, pp 176–178, 1972.

[25] T. J. Ypma, Finding a multiple zero by transformations and Newton-like methods, SIAM Review, Vol. 25 (1983), pp 365–378.

Department of Mathematics, Northeastern Illinois University, Chicago, IL 60626

   E-mail address: zzeng@neiu.edu,   http://www.neiu.edu/∼zzeng