Printer Friendly

Numerical condition of polynomials in different forms.

Abstract. The zeros of high-degree polynomials are notoriously sensitive to changes in the coefficients, causing problems for available zero-finding software. In this paper, we study how this sensitivity depends on the polynomial representation. We first extend the algebraic characterization of polynomial pseudozero sets from the power basis to general bases. We show that for a polynomial, the numerical conditions of its values and zeros are closely related and can be visualized simultaneously by its pseudozero sets. Comparing the pseudozero sets on a set of testing polynomials in the power, Taylor, Chebyshev, and Bernstein bases reveals that appropriate representation of polynomials gives rise to locally well-conditioned zeros, which then leads to an Iterative Refinement Algorithm that combines symbolic formulation with numeric processing to reduce computational errors of polynomial zeros located in the region of interest.

Key words. polynomial equation, polynomial basis, pseudozero set, iterative refinement algorithm.

AMS subject classifications. 65F35.

1. Introduction. The problem of solving a polynomial equation

p(z) = [c.sub.0] + [c.sub.1]z + [c.sub.2][z.sub.2] + ... + [c.sub.n][z.sup.n] = 0 (1.1)

has deeply influenced the development of mathematics since antiquity and continues to play a major role in many areas of science and engineering. Recent sources of interest include computer algebra, computational algebraic geometry, and computer aided geometric design. In these applications, one usually needs to solve (1.1) for large n, typically well above 100 and sometimes on the order of several thousands [7, 9, 11]. The zeros of high-degree polynomials are notoriously sensitive to changes in coefficients, causing problems for available software. However, the dependency of this sensitivity on the choice of polynomial basis has received little attention. Deeper quantitative understanding of this sensitivity dependency would be helpful for improving computational techniques and design of new algorithms.

Condition numbers are commonly used in the study of the sensitivity of polynomials. They are a kind of derivative, enabling us to estimate the magnitude of the changes of the values and zeros of a given polynomial that result from changes in the coefficients. See the classical work of Wilkinson [15]. Some authors have investigated the effect of polynomial bases on the condition numbers. For instance, Gautschi compared the condition numbers of polynomials in the power, orthogonal, and Lagrangian bases [6]. Farouki and Rajan derived the condition numbers for polynomials in the Bernstein basis [4, 5]. In recent years, a geometric approach has been used, which enables one to compute and view the regions in the complex plane in which the zeros may vary as the coefficients are perturbed [3, 10, 13]. With this approach, one can display detailed geometric information not only for simple polynomial zeros, but also for multiple zeros and clusters of zeros. The mathematical insight provided by this approach leads to well-conditioned representations of the problems as well as enhanced algorithms.

This paper takes the geometric approach. In section 2, we first extend the algebraic characterization of the polynomial pseudozero set represented in the power basis, introduced by Mosier [10], Toh and Trefethen [13], to general polynomial bases. We then show that for a polynomial, the numerical conditions of its values and zeros are closely related, and can be visualized simultaneously by its pseudozero sets. In section 3, we display and compare the pseudozero sets for a set of carefully selected polynomials in the power, Taylor, Chebyshev, and Bernstein bases. In section 4, the insight gained from the comparison leads to an Iterative Refinement Algorithm, which is shown to generate local optimal approximations to polynomial zeros under described practical situations. Finally, we present and discuss results of experiments with the algorithm, and point out issues for further investigation.

2. Condition of polynomial evaluation and polynomial zero. Let [P.sub.n] denote the set of polynomials with complex coefficients and degree at most n, and let

B(z) = ([p.sub.0](z), [p.sub.1](z),..., [p.sub.n](z))

be a basis of [P.sub.n]. For each p(z) [member of] [P.sub.n], there are unique [c.sub.k] [member of] C such that

p(z) = [n.summation over (k=0)] [c.sub.k][p.sub.k](z) [P.sub.k] (z). (2.1)

We will also use p to denote the vector of coefficients [([c.sub.0], [c.sub.1],..., [c.sub.n]).sup.T], when there is no danger of confusion. As usual, [[parallel]p[parallel].sub.[infinity]] := [max.sub.0[less than or equal to]k[less than or equal to]n] {[c.sub.k}. For any vector d = [([d.sub.0], [d.sub.1],...,[d.sub.n]).sup.T], [[parallel]x[parallel].sub.d] defines a pseudo-norm on [P.sub.n] defined by

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

Of course, if all [d.sub.k] [not equal to] 0, [[parallel]x[parallel].sub.d] is a d-weighted [infinity]-norm. In particular, [[parallel]p - [??][parallel].sub.d] = [max.sub.0[less than or equal to]k[less than or equal to]n] [absolute value of [d.sub.k]([c.sub.k] - [[??].sub.k])] measures the perturbations in the coefficients of p(z) relative to the weights given by d. We define the vector [d.sup.-1] = [([d.sup.-1.sub.0],...,[d.sup.-1.sub.n]).sup.T], where

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

and we let

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

Below, p is generally fixed and the [c.sub.k]'s are relative to a basis B. We put d := [[parallel]p[parallel].sub.[infinity]] [p.sup.-1] = [[parallel]p[parallel].sub.[infinity]] [([c.sup.-1.sub.0],...,[c.sup.1.sub.n]).sup.T], which corresponds to coefficient-wise perturbations to p. Under this d, we have

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

We define the [epsilon]-neighborhood and the [epsilon]-pseudozero set of p(z) as

[N.sub.[epsilon]](p,B) = {[??] [member of] [P.sub.n] : [[parallel]p - [??][parallel].sub.d] [less than or equal to] [epsilon]},

and

[Z.sub.[epsilon]](p,B) = {z [member of] C : [??](z) = 0 for some [??] [member of] [N.sub.[epsilon]](p,B)}. (2.2)

They are respectively, the set of all polynomials [??] and their zeros obtained by coefficient-wise perturbations of p of size [less than or equal to] [epsilon].

LEMMA 2.1. For any [??] [member of] C, define the polynomials

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

and

[p.sub.[??]](z) = p(z) - p([??])/r([??]) r(z).

Then, we have,

1. [p.sub.[??]]([??]) = 0,

2. [p.sub.[??]](z) [member of] [N.sub.[epsilon]](p,B), with [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], and

3. [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

Proof. Using [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] and [[parallel]r[parallel].sub.d] = 1, follow the proof given by Mosier (see [10, p. 266]).

An algebraic characterization of the [epsilon]-pseudozero set [Z.sub.[epsilon]](p,B) under the power basis was first introduced by Mosier [10], then by Toh and Trefethen [13]. Here we extend their results to general bases B(z).

PROPOSITION 2.2. Let p(z) [member of] [P.sub.n]. If p(z) has the coefficients p = [([c.sub.0], [c.sub.1],...,[c.sub.n]).sup.T] with respect to the basis B(z) = ([p.sub.0](z), [p.sub.1](z),..., [p.sub.n](z)), then

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (2.3)

Proof. Straightforward application of Lemma 2.1.

The quantity [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], which appeared in the [epsilon]-pseudozero set of p(z), also appears in the condition numbers for the polynomial evaluation at z.

LEMMA 2.3. Let p(z) = [[summation].sup.n.sub.k=0] [c.sub.k][p.sub.k](z). Then the absolute and the relative condition numbers for the evaluation of p(z) at [z.sub.0] are respectively

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (2.4)

Proof. Define the map M : [C.sup.n+1] [right arrow] C, which associates with each vector p = [([c.sub.0], [c.sub.1],..., [c.sub.n]).sup.T] [member of] [C.sup.n+1] the value p(z) = [[summation].sup.n.sub.k=0] [c.sub.k][p.sub.k](z). The absolute condition number of M at p is defined as (see [14, p. 90])

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

Since

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

thus [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. The upper bound [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] can be attained for any given p = [([c.sub.0], [c.sub.1],..., [c.sub.n]).sup.T] [member of] [C.sup.n+1] and [z.sub.0] [member of] C by choosing [??] = [([[??].sub.0], [[??].sub.1],..., [[??].sub.n]).sup.T] with [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], because [[parallel][??] - p[parallel].sub.d] = [delta] and [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII].

The relative condition number is

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

The pseudozero set [Z.sub.[epsilon]](p,B) quantifies the conditioning of the zeros of the polynomial p. Using a stable zero-finding algorithm, the computed zeros of p should lie in a region [Z.sub.Cu](p,B), where C = O([[parallel]p[parallel].sub.d]) and u is the machine precision. (O(f) denotes any function whose absolute value is at most c[absolute value of f] for a positive constant c, see [1, p. 7]). If [epsilon] is chosen of the same order as u, i.e., [epsilon] = O(u), then the computed zeros of p should fall into a region that almost coincides with [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (p,B). On the other hand, the condition number [kappa](z, p,B) provides bounds for possible perturbations of the computed value of p(z). When p is evaluated at z by a stable method with the machine precision [epsilon], the relative computing error in evaluating p(z) should be no more than [epsilon][kappa](z, p,B). Although the polynomial evaluations and polynomial zeros are two distinct problems, Proposition 2.2 and Lemma 2.3 reveal an interesting relationship between their numerical conditions:

COROLLARY 2.4.

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (2.5)

That is, when the coefficients of a polynomial p(z) are perturbed, [kappa] and 1/[kappa] provide information on how the values and the zeros, respectively, of p may change; or in other words, the numerical conditions of a polynomial, for both its values and zeros, can be revealed simultaneously by the magnitudes of either [kappa](z, p,B) or 1/[kappa](z, p,B). Since [kappa](z, p,B) = [infinity] at the zeros of p, it is appropriate to use the quantity 1/[kappa](z, p,B), whose [epsilon]-level set is exactly the boundary of [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]-pseudozero set of p(z).

3. Condition of polynomial basis. The conditions of polynomials are basis-dependent, a fact that has been noticed and investigated by several researchers using the condition numbers. For instance, Gautschi introduced the condition numbers of polynomials in the power, orthogonal, and Lagrangian bases [6]. He compared the condition numbers of the zeros of the polynomials (i) and (vii) given below in the power and the Chebyshev bases, and concluded that "it is quite possible that equations that are ill-conditioned in power form become well-conditioned when expanded as orthogonal polynomials, and vice versa." Farouki and Rajan [4] derived the condition numbers for polynomials in the Bernstein basis, and proved that the condition numbers for both simple and multiple polynomial zeros on a real interval ([alpha], [beta]) are always smaller in the Bernstein basis than in the Taylor basis about any s [not member of] ([alpha], [beta]).

Our discussions in the previous section yield a geometric, closed-form representation (eqs. (2.4) and (2.5)) for the conditions of polynomials with respect to any basis B(z) of [P.sub.n]. To be more specific, the influence of a basis on the conditioning of the polynomial (2.1) is determined by the quantity

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

whose [epsilon]-level sets are boundaries of the pseudozero sets [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (p,B), which can be efficiently computed and displayed in the complex plane.

In this section, we examine four bases of [P.sub.n]:

* the power basis: [B.sup.(0)](z) = (1, z, [z.sup.2],...,[z.sup.n]),

* the Taylor basis about s: [B.sup.(s)](z) = (1, z - s, [(z - s).sup.2],..., [(z _ s).sup.n]) (note that the power basis [B.sup.(0)](z) is a special case of the Taylor basis),

* the Chebyshev basis: [B.sup.c](z) = ([T.sub.0](z), [T.sub.1](z),..., [T.sub.n](z)), where [T.sub.k](z) = cos k(arccos z), -1 [less than or equal to] Re(z) [less than or equal to] 1, and

* the Bernstein basis: [B.sup.b](z) = ([b.sup.n.sub.0] (z), [b.sup.n.sub.1] (z),..., [b.sup.n.sub.n](z)), with [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

The power basis is by far the best understood and most commonly used basis for polynomials. The Taylor basis is simply a shift or a variable exchange of the power basis. The Chebyshev basis is an orthogonal basis and has been used in approximation theory, algebra, and number theory, as well as numerical analysis and computation. The Bernstein basis was originally constructed by probabilistic reasoning for use in approximating continuous functions. One of its applications is in parameterizing Bernstein-Bezier curves, which have enjoyed considerable popularity in computer aided design applications due to their elegant geometric properties and the simple recursive algorithms available for processing them.

For each of the bases listed above, we explore the following degree-20 polynomials, most of which are chosen or modified from [13]:

(i) scaled "Wilkinson polynomial": [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

(ii) the Chebyshev polynomial of degree 20: p(z) = cos 20(arccos z), -1 [less than or equal to] Re(z) [less than or equal to] 1.

(iii) the monic polynomial with zeros equally spaced on the curve z = x + i sin([pi]x), -1 [less than or equal to] x [less than or equal to] 1, namely, 2(k + 0.5)/19+ i sin 2[pi](k + 0.5)/19, k = -10,-9,...,9.

(iv) p(z) = [[summation].sup.20.sub.k=0] [(10z).sup.k]/k!.

(v) the monic polynomial with zeros: 10/11 - [2.sup.-k], k = 1, 2,...,20.

(vi) p(z) = [(z - 10/11).sup.20].

(vii) the monic polynomial with zeros: [2.sup.-19], [2.sup.-18],..., 1.

(viii) [p.sub.1](z) = 1 + z + [z.sup.2] + ... + [z.sup.19] + [z.sup.20], [p.sub.2](z) = [(1 + z + [z.sup.2] + ... + [z.sup.10]).sup.2], and [p.sub.3](z) = [(1 + z + [z.sup.2] + ... + [z.sup.5]).sup.4].

Figures 6.1 - 6.9 display the pseudozero sets [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (p,B) whose boundaries are labeled by [log.sub.10] [epsilon]. The experiments were implemented using MATLAB and its Symbolic Math Toolbox. We first input a polynomial p(z) in the power basis. We then converted the coefficients to the coefficients in three other bases symbolically or using 128 decimal precision. Next, we numerically evaluated the function 1/[kappa](z, p,B) at 100 x 100 grid points in a rectangular region that contains all the zeros of p. We finally drew the level sets of -[log.sub.10] [kappa](z, p,B) using the MATLAB function contour.

All the figures confirm that the pseudozero sets of polynomials are basis-dependent, consequently the numerical conditioning of polynomial evaluation and polynomial zero are basis-dependent.

The first phenomenon observed in our experiments is that the components of the pseudozero sets in the power basis [B.sup.(0)] that contain polynomial zeros with small magnitudes are smaller than those containing zeros with larger magnitudes. In particular, the component that contains a root in a neighborhood of z = 0 is significantly smaller than the rest of the components of [Z.sub.[epsilon]](p,[B.sup.(0)]). This observation is verified by eqs. (2.4)-(2.5), since [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] is a monotonic decreasing function of [absolute value of z]:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

Consequently, the component of [Z.sub.[epsilon]](p,[B.sup.(s)]) in the Taylor basis about s that contains the polynomial zeros in the neighborhood of s would be smaller than the other components of [Z.sub.[epsilon]](p,[B.sup.(s)]). The Bernstein basis [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], behaves analogously around z = [alpha] and z = [beta].

Now let us consider polynomial zeros that are far from the origin z = 0. Figures 6.1 - 6.4 show that at well-isolated simple zeros, either real or complex, the pseudozero sets in any of the three bases, power, Chebyshev, and Bernstein, have areas within an order of magnitude of one another. In general, the Bernstein basis has slightly smaller pseudozero sets than the other two. On the other hand, figures 6.5 - 6.7 show that the conditions of clustered or multiple polynomial zeros are generally worse than those of simple zeros for each of the four bases. Between these bases, the Bernstein basis has the smallest pseudozero sets, followed by the Chebyshev basis, and the power basis is the worst (except the case in which a cluster or a multiple zero is located around z = 0, such as the polynomial (vii) for which the power basis behaves better than the Chebyshev basis; also see [6]).

Figures 6.8 and 6.9 reveal that, as polynomial zeros become closer or multiplicities of polynomial zeros become larger, the conditioning of the polynomials in the power basis, and of course in the Taylor basis, deteriorates rapidly, yet the pseudozero sets of the polynomials in the Bernstein basis are relatively unchanged. This phenomenon is also observed in figures 6.5 and 6.6, in which the disparities between the conditioning of the power, Chebyshev, and Bernstein bases become more striking as the cluster of polynomial zeros [[xi].sub.k] = 10/11 - [2.sup.-k], k = 1,..., 20, merges into a single 20-fold zero [xi] = 10=11.

4. Accuracy refinement of polynomial zeros. The analysis and experiments given in the previous sections demonstrate that polynomial bases that have smaller [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] in the region of interest give rise to locally better conditioned representations, which suggests that conditioning of polynomials can be improved through careful basis transformation. Since the zero-finding procedure used for the power basis can be applied directly to the Taylor basis representation by a variable substitution, an immediate application of this idea is to use the Taylor basis [B.sup.(s)](z) with s chosen in the region of interest to improve local numerical condition of polynomials that are initially represented in the power basis. Another application is to refine approximated polynomial zeros to optimal accuracy using a sequence of Taylor bases [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII].

Let p(z) [member of] [P.sub.n] be represented in the power basis [B.sup.(0)](z). With no loss of generality, we assume p(0) [not equal to] 0. (Otherwise p(z) = [z.sup.m]g(z) with g(0) [not equal to] 0, the same discussions are applied to g(z).) For any two Taylor bases [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII](z) and [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII](z), assuming O([absolute value of p([s.sub.2])]) = O([absolute value of p([s.sub.1])]) as [s.sub.i] [right arrow] [xi] with p(0) = 0, if [absolute value of z - [s.sub.2]] << [absolute value of z - [s.sub.1]], then

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

which implies that when [absolute value of [xi] - [s.sub.2]] << [absolute value of [xi] - [s.sub.1]] with p([xi]) = 0, the conditioning of the zero [xi] under the basis [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII](z) is likely better than using the basis [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII](z). Proposition 4.1 formalizes this argument.

Let p(z) [member of] [P.sub.n]. Suppose p(z) has [n.sub.0] well-separated zeros (here, a cluster of zeros is considered as a multiple zero). For [epsilon] > 0 small enough, the set [Z.sub.[epsilon]](p,B) has [n.sub.0] disjoint connected components, each of which contains a zero of p(z) (see [10]). We use [Z.sup.([xi]).sub.[epsilon]] (p,B) to denote the connected component of [Z.sub.[epsilon]](p,B) that contains the zero [xi].

PROPOSITION 4.1. Let [xi] be an m-fold zero of p(z), and [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] C : [absolute value of z - [xi]] [less than or equal to] [[delta].sub.0]}, where [[delta].sub.0] is chosen such that [Z.sub.1] is a small neighborhood of [xi]. Define

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

where [s.sub.k-1] [member of] [Z.sub.k-1], and [[delta].sub.k] = max{[absolute value of z - [??]] : [for all]z([??] [member of] [Z.sub.k])}: Then,

[absolute value of p(z)] [less than or equal to] O([[epsilon].sup.k]) for any z [member of] [Z.sub.k], and [[delta].sub.]k [less than or equal to] O([[epsilon].sup.k/m]). (4.1)

Proof. We prove the proposition by induction on k. Let T denote the function

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

(note that d = [[parallel]p[parallel].sub.[infinity]] [p.sup.-1] and p(z) = [[summation].sup.n.sub.i=0] [p.sup.(i)](s)/i! [(z - s).sup.i]). For k = 1, [for all]z [member of] [Z.sub.1], [absolute value of p(z)] [less than or equal to] [epsilon]T (z, 0) = O([epsilon]), because T(z, 0) is bounded on [Z.sub.1]. Then, we have,

[absolute value of z - [xi]] = [[absolute value of p(z)/[p.sup.(m)]([xi])/m! + O(z-[xi])].sup.1/m] [less than or equal to] O([[absolute value of p(z)].sup.1/m]) [less than or equal to] O([[epsilon].sup.1/m]),

leading to [absolute value of z - [??]| [less than or equal to] [absolute value of z - [xi]] + [absolute value of [??] - [xi]] [less than or equal to] O([[epsilon].sup.1/m]) for any z and [??] [member of] [Z.sub.1]. Thus [[delta].sub.1] [less than or equal to] O([[epsiln].sup.1/m]).

Now, suppose [absolute value of p(z)] [less than or equal to] O([[epsilon].sup.k]), [for all]z [member of] [Z.sub.k] and [[delta].sub.k] [less than or equal to] O([[epsilon].sup.k/m]) for [member of] [less than or equal to] k < m. Then, for [s.sub.k] [member of] [Z.sub.k],

[absolute value of p([s.sub.k])] [less than or equal to] O([[epsilon].sup.k]),

and

[absolute value of [p.sup.(i)]([s.sub.k])] = [absolute value of [p.sup.(m)/(m - i)! [([s.sub.k] - [xi]).sup.m-i] + O[(([s.sub.k] - [xi]).sup.m-i+1])] [less than or equal to] ([[delta].sup.m-i.sub.k]), i = 1,2...m-1.

Thus, [for all]z [member of] [Z.sub.k+1],

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

which gives

[absolute value of p(z)] [less than or equal to] [epsilon]T (z, sk) [less than or equal to] O([[epsilon].sup.k+1]).

Using the same discussion as in the case k = 1, we obtain [absolute value of z - [??]] [less than or equal to] O([[epsilon].sup.(k+1)/m]), for any z and [??] [member of] [Z.sub.k+1], i.e., [[delta].sub.k+1] [less than or equal to] O([[epsilon].sup.(k+1)/m]).

Proposition 4.1 yields a constructive way to refine the accuracy of approximations to polynomial zeros. Based on this, we designed the following Iterative Refinement Algorithm.

Algorithm: [[??].sub.new] = REFI(p, [[??].sub.old], tol)

Input: p: polynomial with complex coefficients and degree n. [[??].sub.old]: an approximation to [xi], where [xi] is an m-fold zero of p. tol: error tolerance.

Output: [[??].sub.new]: [[absolute value of [[??].sub.new] - [xi]] < [absolute value of [[??].sub.new] x tol, or a new approximation to [xi].

Set [[??].sub.1] = [[??].sub.old].

For k = 1,...,m - 1

(i) Form p(z) = [[summation].sup.n.sub.i=0] [a.sub.i][(z - [[??].sub.k]).sup.i] = f(y), where y = z - [[??].sub.k].

(ii) Compute [[??].sub.f], the root of f(y) with the smallest magnitude by a stable zero-finding numerical method.

(iii) Set [[??].sub.k+1] = [[??].sub.k] + [[??].sub.f].

(iv) If [absolute value of [[??].sub.f]] < [absolute value of [[??].sub.k+1]] x tol, stop.

Set [[??].sub.new] = [[??].sub.k+1].

The algorithm generates a set of approximations {[[??].sub.k]} to [xi], an m-fold zero of p(z). Assuming that [epsilon] is of the order of the machine precision, and that the input [[??].sub.old] [member of] [Z.sub.1] = [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (p,[B.sup.(0)]) is well separated from other connected components of [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], then each [[??].sub.k] is expected to fall into [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], since [[??].sub.k] is computed by a stable zero-finding numerical method for p(z) represented in the basis [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII](z). Although it is impossible to rigorously establish that each [[??].sub.k] also lies inside [Z.sub.k-1] as required by Proposition 4.1, we can reasonably conjecture that when _ is small enough,

[[delta].sub.k] = O([[epsilon].sup.k/m]) << O([[epsilon].sup.(k-1)/m]) = [[delta].sub.k-1], k = 2, 3,...,

i.e., [Z.sub.k] is much smaller than [Z.sub.k-1] if they have similar geometric shapes. Since [xi] is a common interior point of [Z.sub.k-1] and [Z.sub.k], it is likely that [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], and consequently it is highly possible that [[??].sub.k] [member of] [Z.sub.k]. According to Proposition 4.1, {[[??].sub.k]} then converges to [xi] linearly at the rate [[epsilon].sup.1/m], and it would take at most m - 1 refinements to get an optimal approximation within the given machine precision.

The Taylor coefficients p = [([a.sub.0], [a.sub.1],..., [a.sub.n]).sup.T] in Step (i) can be computed by using from about 9n [log.sub.2] n to 18n [log.sub.2] n arithmetic operations [1]. Since basis conversions from the power basis are usually numerically unstable [5], the coefficients in the new bases should be computed either symbolically or numerically with high digit precision.

The algorithm can be used interactively. For example, after a set of desirable zeros of the polynomial p(z) is computed by a zero-finding routine, the pseudozero set can be plotted around the computed zeros as shown in figures 6.1 - 6.9. If a suspect ill-conditioned zero is viewed, say [??], we may implement REFI with [[??].sub.old] = [??] ([[??].sub.old] can also be input by clicking the mouse on the computer screen at the center of the suspect component of the pseudozero set). The zeros can be refined either one at a time, or simultaneously on several subsets. For the former, after each zero is refined, the polynomial is deflated and REFI is then applied to the deflated polynomial (care should be taken to ensure deflation stability, see [8, 12, 15] for details); while for the latter, since the refinement of one zero is independent of the others, REFI can be applied concurrently to several subsets of the computed zeros on parallel computers or on clusters of workstations.

5. Numerical experiments. We implemented REFI for the polynomials given below on a SUN Ultra 170E workstation. Our programs were written in the MATLAB language, which uses double precision numeric arithmetic. We first computed numerical zeros {[[xi].sup.(0).sub.i]} in the power basis by the MATLAB roots function. We then implemented REFI with [[??].sub.old] chosen in the desired region. The Taylor coefficients were computed symbolically using the MATLAB Symbolic Math Toolbox. There was no noticeable time delay for using symbolic formulations in our implementations.

Example 5.1. The square of "Wilkinson polynomial": p(z) = [[product].sup.20.sub.k=1][(z - k).sup.2].

We first computed all 40 numerical zeros in the power basis. We then plotted them (marked by green o) together with [Z.sub.[epsilon]](p,[B.sup.(0)]) (marked by green curves) in figure 6.10 (left), in which [epsilon] = eps[[parallel]p[parallel].sub.d], eps = 2.2 x [10.sup.16] is the MATLAB machine precision. As shown, 34 numerical zeros fell into the same connected component of [Z.sub.[epsilon]](p,[B.sup.(0)]), indicating that they could be ill-conditioned. Therefore, we implemented one iteration of REFI with [[??].sub.old] = 14, because 14 is located roughly in the center of this suspect component. The roots of f(y) were also computed using the MATLAB roots function, which generated 40 new numerical zeros of p(z) (marked by blue x). Next, we plotted [Z.sub.[epsilon]](p,[B.sup.(14)]) (marked by blue curves) in the same figure. Among a total of 80 numerical zeros, those that lie inside the set [Z.sub.[epsilon]](p,[B.sup.(0)]) [intersection] [Z.sub.[epsilon]](p,[B.sup.(14)]) are better approximations. All 40 x's and 9 o's fell into this set. Among them, we picked 9 o's (because they are closer to s = 0 than those x's to s = 14 in the same region) and 31 x's that are closer to s = 14 as new approximations to the zeros of p(z). figure 6.10 (right) shows the approximation errors before REFI (marked by green o), after one iteration of REFI (marked by blue x), and the new approximations (marked by dot-dash line).

Example 5.2. The monic polynomial with the n-fold zero [xi] = 10/11: p(z) = [(z - 10/11).sup.n].

We first computed n zeros {[[xi].sup.(0).sub.i]} in the power basis [B.sup.(0)](z). We then arbitrarily picked an [[xi].sup.(0).sub.i] as [[??].sub.old] and implemented REFI with tol = [10.sup.-15]. The approximations at the k-th iteration are denoted by {[[xi].sup.(k).sub.i]}. Table 5.1 gives, when the multiplicity of the roots increases, the number of refinements (k), the accuracy improvement made by REFI, and the computed [epsilon] that appears in the rate of convergence. In the 5-th column, we list [[delta].sub.0] = [[epsilon].sup.1/m.sub.0], where [[epsilon].sub.0] = [2.sup.-1022] is the underflow threshold for IEEE double precision (see [2, p. 11]), m is the multiplicity of the zero [xi] = 10/11 (here, m = n). It is interesting to notice that the approximations achieved by REFI (the 4-th column) almost coincide with the m-th root of the machine underflow threshold, [[delta].sub.0]. This observation not only confirms the inequalities (4.1) proved by Proposition 4.1, but also further indicates that any numerical approximation [??] to an m-fold zero of p(z), [??], is subject to the constraint

O([absolute value of [??] - [xi]]) [greater than or equal to] [[delta].sub.0] = (the machine underflow threshold).sup.1/m],

i.e., [absolute value of [??] - [xi]].sup.m] = O([absolute value of p([??])]) > 0 has to be machine representable. The numerical results suggest that REFI is able to improve computed polynomial zeros to the optimal numerical approximations under the given machine precision.

For n = 40, figure 6.11 (left) shows {[[xi].sup.(0).sub.i]} (green o's), {[[xi].sup.(k).sub.i]} (x's, each color represents an iteration of REFI) obtained by the first few refinements (k = 1, 2, ...), figure 6.11 (right) shows {[[xi].sup.(k).sub.i]} (x's), obtained from the last few refinements (k = ...,34, 35) and the exact multiple zero [xi] (red +). Figure 6.12 shows [min.sub.i] [absolute value of [[xi].sup.(k).sub.i] - [xi]], the error reduction of REFI for n = 40. It confirms the linear convergence of the algorithm with the rate [[epsilon].sup.1=m].

We attempted to refine all zeros computed in the power basis, [[xi].sup.(0).sub.i], i = 1,...,n, for a set of test polynomials. For a polynomial p(z) satisfying

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (5.1)

where [epsilon] = eps[[parallel]p[parallel].sub.d], and {[[xi].sub.i]} are exact zeros of p(z), REFI is guaranteed to produce the optimal approximations to each [[xi].sub.i] (in the meaning given by Example 5.2) by choosing [[??].sub.old] = [[xi].sup.(0).sub.i]. Table 5.2 gives numerical results of REFI on

Example 5.3. The scaled "Wilkinson polynomial": p(z) = [[product].sup.20.sub.k=1] (z - k/20).

For polynomials that fail to meet condition (5.1), i.e., when there exists [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] that contains more than one distinct zero of p(z), REFI is able to generate at least one optimal approximation to a zero [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. Conceptually, the rest of the zeros in [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] can be approximated to their optimal approximations by removing the accepted approximation of [[xi].sub.j] in the manner similar to the approaches of [8, 12, 15], and then applying REFI to the deflated polynomial. However, we encounted numerical difficulties in deflation stability for high-degree polynomials with multiple zeros, e.g., p(z) in Example 5.1. Since polynomial deflation is a challenging research topic of its own, and is beyond the consideration of this paper, we would leave it for future investigation.

6. Concluding remarks.

* REFI is independent of the zero-finding algorithm by which the original zero approximations {[[xi].sup.(0).sub.i]} are computed, and can be adapted to any numerical scheme that is used for computing zeros of f(y) at Step (ii) to reduce the approximation errors in {[[xi].sup.(0).sub.i]}.

* REFI is a local algorithm mixing symbolic and numerical computations. For low-degree polynomials, computational speed is not a concern. The high-degree polynomials encountered in practice are usually sparse. Moreover, only a small subset of their zero set is generally of interest (see [7]). In this situation, REFI could be implemented at reasonable cost by exploring sparse data structure at Step (i) and by adopting more efficient numerical local zero-finding algorithm at Step (ii). Readers may refer to [11] for available local zero-finding methods.

[FIGURE 6.1 OMITTED]

[FIGURE 6.2 OMITTED]

[FIGURE 6.3 OMITTED]

[FIGURE 6.4 OMITTED]

[FIGURE 6.5 OMITTED]

[FIGURE 6.6 OMITTED]

[FIGURE 6.7 OMITTED]

[FIGURE 6.8 OMITTED]

[FIGURE 6.9 OMITTED]

[FIGURE 6.10 OMITTED]

[FIGURE 6.11 OMITTED]

[FIGURE 6.12 OMITTED]

Acknowledgements. I would like to thank Prof. James Madden of Louisiana State University for interesting discussions and careful reading of the manuscript.

REFERENCES

[1] D. BINI AND V. PAN, Polynomial and Matrix Computations, Vol. 1: Fundamental Algorithms, Birkhauser, Boston, Cambridge, MA, 1994.

[2] J. W. DEMMEL, Applied Numerical Linear Algebra, SIAM, 1997.

[3] A. EDELMAN AND H. MURAKAMI, Polynomial roots from companion matrix eigenvalues, Math. Comp., 64 (1995), pp. 763-776.

[4] R. FAROUKI AND V. RAJAN, On the numerical condition of polynomials in Bernstein form, Comput. Aided Geom. Design, 4 (1987), pp. 191-216.

[5] --, Algorithms for polynomials in Bernstein form, Comput. Aided Geom. Design, 5 (1988), pp. 1-26.

[6] W. GAUTSCHI, Questions of numerical condition related to polynomials, in Studies in Numerical Analysis, G. H. Golub, ed., MAA Studies in Math., 24 (1984), pp. 140-177.

[7] L. GONZALEZ-VEGA, The needs of industry for polynomial system solving, FRISCO Consortium (LTR 21.024), 1997.

[8] M. A. JENKINS AND J. TRAUB, Algorithm 419: Zeros of a complex polynomial, Comm. ACM, 15 (1972), pp. 97-99.

[9] J. MCNAMEE, A bibliography on roots of polynomials, J. Comput. Appl. Math., 47 (1993), pp. 391-394.

[10] R. MOSIER, Root neighborhoods of a polynomial, Math. Comp., 47 (1986), pp. 265-273.

[11] V. Y. PAN, Solving a polynomial equation: some history and recent progress, SIAM Rev., 39 (1997), pp. 187-220.

[12] B. PARLETT, Laguerre's method applied to the matrix eigenvalue problem, Math. Comp., 18 (1964), pp. 464-485.

[13] K. C. TOH AND L. N. TREFETHEN, Pseudozeros of polynomials and pseudospectra of companion matrices, Numer. Math., 68 (1994), pp. 403-425.

[14] L. N. TREFETHEN AND D. BAU, Numerical Linear Algebra, SIAM, 1997.

[15] J. WILKINSON, Rounding Errors in Algebraic Processes, Prentice-Hall, Englewood Cliffs. NJ, 1963.

* Received February 13, 2000. Accepted for publication January 26, 2001. Recommended by G. S. Ammar.

HONG ZHANG ([dagger])

([dagger]) Department of Computer Science, Illinois Institute of Technology, Chicago, IL 60616, U.S.A. Email: hzhang@mcs.anl.gov
TABLE 5.1

Error reduction for the zeros of p(z) = [(z - 10/11).sup.n]

 No. of
 Refi. Error before REFI Error after REFI

 [[min.sub.i] | [[min.sub.i] |
 [xi].sub.i.sup.(0)] [xi].sub.i.sup.(k)]
n k - [xi]| - [xi]|

10 12 4.45070e-02 1.1111e-16
20 27 2.90500e-01 1.7870e-15
30 32 2.46660e-01 1.5771e-11
40 35 7.37440e-01 8.0332e-09
50 37 9.26390e-01 3.3792e-07

 Computed

n [[delta].sub.0] [epsilon]

10 1.7169e-31 1.5859e-13
20 4.1435e-16 1.0656e-11
30 5.5579e-11 1.0890e-10
40 2.0356e-08 3.5192e-10
50 7.0299e-07 8.0302e-10

TABLE 5.2

Error reduction for the zeros of p(z) = [PI].sub.k=1.sup.20])
(z - k/20).

 Error before REFI Error after REFI No. of
 Refi.
[[xi].sub.i] |[[xi].sub.i.sup.(0)] |[[xi].sub.i.sup.(k)]
 - [[xi].sub.i] - [[xi].sub.i] k

1/20 1.6453e-13 0 1
2/10 2.2462e-12 0 1
3/20 2.7104e-10 0 1
4/20 1.5041e-09 0 1
5/20 1.3142e-07 0 2
6/20 1.9093e-06 0 2
7/20 1.5435e-05 0 2
8/20 8.5185e-05 0 2
9/20 3.5601e-04 0 2
10/20 1.1012e-03 0 2
11/20 2.9596e-03 0 2
12/20 5.2658e-03 0 2
13/20 9.2826e-03 0 2
14/20 1.0416e-02 0 2
15/20 8.8301e-03 0 2
16/20 6.5150e-03 0 2
17/20 2.7471e-03 0 1
18/20 9.8058e-04 0 1
19/20 1.9272e-04 1.1102e-16 1
20/20 1.8238e-05 0 1
COPYRIGHT 2001 Institute of Computational Mathematics
No portion of this article can be reproduced without the express written permission from the copyright holder.
Copyright 2001 Gale, Cengage Learning. All rights reserved.

Article Details
Printer friendly Cite/link Email Feedback
Author:Zhang, Hong
Publication:Electronic Transactions on Numerical Analysis
Date:Jan 1, 2001
Words:6403
Previous Article:Error analysis of QR algorithms for computing Lyapunov exponents.
Next Article:On parallel two-stage methods for Hermitian positive definite matrices with applications to preconditioning.


Related Articles
A perturbation approach for approximate inertial manifolds.
Special volume on orthogonal polynomials, approximation theory, and harmonic analysis.
Polynomial eigenvalue problems with Hamiltonian structure.
Stability and sensitivity of Darboux transformation without parameter.
Some theoretical results derived from polynomial numerical hulls of Jordan blocks.
On recurrence relations for radial wave functions for the N-th dimensional oscillators and hydrogenlike atoms: analytical and nuberical study.
Bivariate interpolation at Xu points: results, extensions and applications.
On the estimation of the q-numerical range of monic matrix polynomials.
Improved initialization of the accelerated and robust QR-like polynomial root-finding.

Terms of use | Copyright © 2017 Farlex, Inc. | Feedback | For webmasters