Jenkins–Traub algorithm explained
The Jenkins–Traub algorithm for polynomial zeros is a fast globally convergent iterative 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".[1]
This article describes the complex variant. Given a polynomial P,with complex coefficients it computes approximations to the n zeros
\alpha1,\alpha2,...,\alphan
of
P(
z), one at a time in roughly increasing order of magnitude. After each root is computed, its linear factor is removed from the polynomial. Using this
deflation guarantees that each root is computed only once and that all roots are found.
The real variant follows the same pattern, but computes two roots at a time, either two real roots or a pair of conjugate complex roots. By avoiding complex arithmetic, the real variant can be faster (by a factor of 4) than the complex variant. The Jenkins–Traub algorithm has stimulated considerable research on theory and software for methods of this type.
Overview
The Jenkins–Traub algorithm calculates all of the roots of a polynomial with complex coefficients. The algorithm starts by checking the polynomial for the occurrence of very large or very small roots. If necessary, the coefficients are rescaled by a rescaling of the variable. In the algorithm, proper roots are found one by one and generally in increasing size. After each root is found, the polynomial is deflated by dividing off the corresponding linear factor. Indeed, the factorization of the polynomial into the linear factor and the remaining deflated polynomial is already a result of the root-finding procedure. The root-finding procedure has three stages that correspond to different variants of the inverse power iteration. See Jenkins and Traub.[2] A description can also be found in Ralston and Rabinowitz[3] p. 383.The algorithm is similar in spirit to the two-stage algorithm studied by Traub.[4]
Root-finding procedure
Starting with the current polynomial P(X) of degree n, the aim is to compute the smallest root
of
P(x). The polynomial can then be split into a linear factor and the remaining polynomial factor
- for stage 2, if s is close enough to
:
H^(X) = P_1(X) +O\left(\left|\frac\right|^M \cdot \left|\frac\right|^\right) and
s-\frac = \alpha_1+O\left(\ldots\cdot|\alpha_1-s|\right).
H^(X) =P_1(X) +O\left(\prod_^ \left|\frac\right| \right) and s_= s_\lambda-\frac =\alpha_1+O\left(\prod_^ \left|\frac\right| \cdot \frac
\right) giving rise to a higher than quadratic convergence order of
, where
is the
golden ratio.
Interpretation as inverse power iteration
All stages of the Jenkins–Traub complex algorithm may be represented as the linear algebra problem of determining the eigenvalues of a special matrix. This matrix is the coordinate representation of a linear map in the n-dimensional space of polynomials of degree n - 1 or less. The principal idea of this map is to interpret the factorizationP(X)=(X-\alpha_1)\cdot P_1(X)with a root
and
the remaining factor of degree
n - 1 as the eigenvector equation for the multiplication with the variable
X, followed by remainder computation with divisor
P(
X),
M_X(H) = (X\cdot H(X)) \bmod P(X)\,.This maps polynomials of degree at most
n - 1 to polynomials of degree at most
n - 1. The eigenvalues of this map are the roots of
P(
X), since the eigenvector equation reads
0 = (M_X-\alpha\cdot id)(H)=((X-\alpha)\cdot H) \bmod P\,,which implies that
, that is,
is a linear factor of
P(
X). In the monomial basis the linear map
is represented by a
companion matrix of the polynomial
P, as
M_X(H) = \sum_^H_mX^-H_\left(X^n+\sum_^a_mX^m\right) = \sum_^(H_-a_H_)X^m-a_0H_\,,the resulting transformation matrix is
A=\begin0 & 0 & \dots & 0 & -a_0 \\1 & 0 & \dots & 0 & -a_1 \\0 & 1 & \dots & 0 & -a_2 \\\vdots & \vdots & \ddots & \vdots & \vdots \\0 & 0 & \dots & 1 & -a_\end\,.To this matrix the
inverse power iteration is applied in the three variants of no shift, constant shift and generalized Rayleigh shift in the three stages of the algorithm. It is more efficient to perform the linear algebra operations in polynomial arithmetic and not by matrix operations, however, the properties of the inverse power iteration remain the same.
Real coefficients
The Jenkins–Traub algorithm described earlier works for polynomials with complex coefficients. The same authors also created a three-stage algorithm for polynomials with real coefficients. See Jenkins and Traub A Three-Stage Algorithm for Real Polynomials Using Quadratic Iteration.[5] The algorithm finds either a linear or quadratic factor working completely in real arithmetic. If the complex and real algorithms are applied to the same real polynomial, the real algorithm is about four times as fast. The real algorithm always converges and the rate of convergence is greater than second order.
A connection with the shifted QR algorithm
There is a surprising connection with the shifted QR algorithm for computing matrix eigenvalues. See Dekker and Traub The shifted QR algorithm for Hermitian matrices.[6] Again the shifts may be viewed as Newton-Raphson iteration on a sequence of rational functions converging to a first degree polynomial.
Software and testing
The software for the Jenkins–Traub algorithm was published as Jenkins and Traub Algorithm 419: Zeros of a Complex Polynomial.[7] The software for the real algorithm was published as Jenkins Algorithm 493: Zeros of a Real Polynomial.[8]
The methods have been extensively tested by many people. As predicted they enjoy faster than quadratic convergence for all distributions of zeros.
However, there are polynomials which can cause loss of precision[9] as illustrated by the following example. The polynomial has all its zeros lying on two half-circles of different radii. 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,[10] p. 64. The original polynomial was of degree 60 and suffered severe deflation instability.
External links
Notes and References
- 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.
- Jenkins, M. A. and Traub, J. F. (1970), A Three-Stage Variables-Shift Iteration for Polynomial Zeros and Its Relation to Generalized Rayleigh Iteration, Numer. Math. 14, 252–263.
- Ralston, A. and Rabinowitz, P. (1978), A First Course in Numerical Analysis, 2nd ed., McGraw-Hill, New York.
- Traub, J. F. (1966), A Class of Globally Convergent Iteration Functions for the Solution of Polynomial Equations, Math. Comp., 20(93), 113–138.
- Jenkins, M. A. and Traub, J. F. (1970), A Three-Stage Algorithm for Real Polynomials Using Quadratic Iteration, SIAM J. Numer. Anal., 7(4), 545–566.
- Dekker, T. J. and Traub, J. F. (1971), The shifted QR algorithm for Hermitian matrices, Lin. Algebra Appl., 4(2), 137–154.
- Jenkins, M. A. and Traub, J. F. (1972), Algorithm 419: Zeros of a Complex Polynomial, Comm. ACM, 15, 97–99.
- Jenkins, M. A. (1975), Algorithm 493: Zeros of a Real Polynomial, ACM TOMS, 1, 178–189.
- Web site: 8 August 2005. William Kahan Oral history interview by Thomas Haigh. 2021-12-03. The History of Numerical Analysis and Scientific Computing. Philadelphia, PA.
- Wilkinson, J. H. (1963), Rounding Errors in Algebraic Processes, Prentice Hall, Englewood Cliffs, N.J.