Jenkins–Traub algorithm: Difference between revisions

Revision as of 20:49, 10 August 2023

Latest revision as of 22:14, 16 August 2025

(6 intermediate revisions by 6 users not shown)

Line 1:

{{Short description|Root-finding algorithm for polynomials}}

The '''Jenkins–Traub algorithm for polynomial zeros''' is a fast globally convergent iterative [[Root-finding algorithms#Roots of polynomials|polynomial root-finding]] method published in 1970 by [[Michael A. Jenkins]] and [[Joseph F. Traub]]. They gave two variants, one for general polynomials with complex coefficients, commonly known as the "CPOLY" algorithm, and a more complicated variant for the special case of polynomials with real coefficients, commonly known as the "RPOLY" algorithm. The latter is "practically a standard in black-box polynomial root-finders".<ref>Press, W. H., Teukolsky, S. A., Vetterling, W. T. and Flannery, B. P. (2007), Numerical Recipes: The Art of Scientific Computing, 3rd ed., Cambridge University Press, page 470.</ref>

Line 18 ⟶ 19:

To that end, a sequence of so-called ''H'' polynomials is constructed. These polynomials are all of degree ''n''&nbsp;&minus;&nbsp;1 and are supposed to converge to the factor <math>\bar H(X)</math> of ''P''(''X'') containing (the linear factors of) all the remaining roots. The sequence of ''H'' polynomials occurs in two variants, an unnormalized variant that allows easy theoretical insights and a normalized variant of <math>\bar H</math> polynomials that keeps the coefficients in a numerically sensible range.

The construction of the ''H'' polynomials <math>\left(H^{(\lambda)}(z)\right)_{\lambda=0,1,2,\dots}</math> is guided by a sequence of [[complex numbersnumber]]s <math>(s_\lambda)_{\lambda=0,1,2,\dots}</math> called shifts. These shifts themselves depend, at least in the third stage, on the previous ''H'' polynomials. The ''H'' polynomials are defined as the solution to the implicit recursion

<math display="block">

H^{(0)}(z)=P^\prime(z)

Line 183 ⟶ 184:

\frac{|\alpha_1-s_\lambda|^2}{|\alpha_2-s_\lambda|}

\right)

</math> giving rise to a higher than quadratic convergence order of <math>\phi^2=1+\phi\approx 2.61618</math>, where <math>\phi=\tfrac12(1+\sqrt5)</math> is the [[golden ratio]].

=== Interpretation as inverse power iteration ===

Line 192 ⟶ 193:

This maps polynomials of degree at most ''n''&nbsp;&minus;&nbsp;1 to polynomials of degree at most ''n''&nbsp;&minus;&nbsp;1. The eigenvalues of this map are the roots of ''P''(''X''), since the eigenvector equation reads

<math display="block">0 = (M_X-\alpha\cdot id)(H)=((X-\alpha)\cdot H) \bmod P\,,</math>

which implies that <math>(X-\alpha)\cdot H)=C\cdot P(X)</math>, that is, <math>(X-\alpha)</math> is a linear factor of ''P''(''X''). In the monomial basis the linear map <math>M_X</math> is represented by a [[companion matrix]] of the polynomial ''P'', as

<math display="block"> M_X(H) = \sum_{m=0}^{n-1}H_mX^{m+1}-H_{n-1}\left(X^n+\sum_{m=0}^{n-1}a_mX^m\right) = \sum_{m=1}^{n-1}(H_{m-1}-a_{m}H_{n-1})X^m-a_0H_{n-1}\,,</math>

the resulting transformation matrix is

Line 215 ⟶ 216:

The methods have been extensively tested by many people.{{Who|date=December 2021}} As predicted they enjoy faster than quadratic convergence for all distributions of zeros.

However, there are polynomials which can cause loss of precision<ref>{{Cite web|date=8 August 2005|title=William Kahan Oral history interview by Thomas Haigh|url=http://history.siam.org/oralhistories/kahan.htm|url-status=live|access-date=2021-12-03|website=The History of Numerical Analysis and Scientific Computing|publication-place=Philadelphia, PA}}</ref> as illustrated by the following example. The polynomial has all its zeros lying on two half-circles of different radii. [[James H. Wilkinson|Wilkinson]] recommends that it is desirable for stable deflation that smaller zeros be computed first. The second-stage shifts are chosen so that the zeros on the smaller half circle are found first. After deflation the polynomial with the zeros on the half circle is known to be ill-conditioned if the degree is large; see Wilkinson,<ref>Wilkinson, J. H. (1963), Rounding Errors in Algebraic Processes, Prentice Hall, Englewood Cliffs, N.J.</ref> p.&nbsp;64. The original polynomial was of degree 60 and suffered severe deflation instability.

==References==

Line 228 ⟶ 229:

{{DEFAULTSORT:Jenkins-Traub Algorithm}}

[[Category:Numerical analysis]]

[[Category:Root-findingPolynomial factorization algorithms]]