Jenkins–Traub algorithm: Difference between revisions

Revision as of 19:43, 28 April 2023

Revision as of 15:02, 10 August 2023

Line 15:

=== Root-finding procedure ===

Starting with the current polynomial ''P''(''X'') of degree ''n'', the aim is to compute the smallest root <math>\alpha</math> of ''P(x)'' is computed. ToThe thatpolynomial end,can athen sequencebe ofsplit so-calledinto ''H''a polynomialslinear is constructed. These polynomials are all of degree ''n''&nbsp;&minus;&nbsp;1factor and arethe supposedremaining to converge to thepolynomial factor of<math ''display="block>P''(''X'')=(X-\alpha)\bar containing all the remaining roots. The sequence of ''H''(X)</math> polynomialsOther occursroot-finding inmethods twostrive variants,primarily anto unnormalizedimprove variantthe thatroot allowsand easythus theoreticalthe insightsfirst andfactor. aThe normalizedmain variantidea of <math>\barthe H</math>Jenkins-Traub polynomialsmethod thatis keepsto theincrementally coefficientsimprove inthe asecond numerically sensible rangefactor.

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> dependsis onguided by a sequence of complex numbers <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 33 ⟶ 34:

where the polynomial division is exact.

Algorithmically, one would use forlong instancedivision by the linear factor as in the [[Horner scheme]] or [[Ruffini rule]] to evaluate the polynomials at <math>s_\lambda</math> and obtain the quotients at the same time. With the resulting quotients ''p''(''X'') and ''h''(''X'') as intermediate results the next ''H'' polynomial is obtained as

<math display="block">

\left.\begin{align}

Line 55 ⟶ 56:

==== Stage one: no-shift process ====

For <math>\lambda = 0,1,\dots, M-1</math> set <math>s_\lambda=0</math>. Usually ''M=5'' is chosen for polynomials of moderate degrees up to ''n''&nbsp;=&nbsp;50. This stage is not necessary from theoretical considerations alone, but is useful in practice. It emphasizes in the ''H'' polynomials the cofactor(s) (of the linear factor) of the smallest root(s).

==== Stage two: fixed-shift process ====

The shift for this stage is determined as some point close to the smallest root of the polynomial. It is quasi-randomly located on the circle with the inner root radius, which in turn is estimated as the positive solution of the equation

<math display="block">

Line 64 ⟶ 67:

Since the left side is a convex function and increases monotonically from zero to infinity, this equation is easy to solve, for instance by [[Newton's method]].

Now choose <math>s=R\cdot \exp(i\,\phi_\text{random})</math> on the circle of this radius. The sequence of polynomials <math>H^{(\lambda+1)}(z)</math>, <math>\lambda=M,M+1,\dots,L-1</math>, is generated with the fixed shift value <math>s_\lambda = s</math>. This creates an asymmetry relative to the previous stage which increases the chance that the ''H'' polynomial moves towards the cofactor of a single root.
During this iteration, the current approximation for the root

<math display="block">t_\lambda=s-\frac{P(s)}{\bar H^{(\lambda)}(s)}</math>

is traced. The second stage is finishedterminated successfullyas successful if the conditions

<math display="block">

|t_{\lambda+1}-t_\lambda|<\tfrac12\,|t_\lambda|

Line 72 ⟶ 77:

|t_\lambda-t_{\lambda-1}|<\tfrac12\,|t_{\lambda-1}|

</math>

are simultaneously met. This limits the relative step size of the iteration, ensuring that the approximation sequence stays in the range of the smaller roots. If there was no success after some number of iterations, a different random point on the circle is tried. Typically one uses a number of 9 iterations for polynomials of moderate degree, with a doubling strategy for the case of multiple failures.

==== Stage three: variable-shift process ====

The <math>H^{(\lambda+1)}(X)</math> polynomials are now generated using the variable shifts <math>s_{\lambda},\quad\lambda=L,L+1,\dots</math> which are generated by

<math display="block">s_L = t_L = s- \frac{P(s)}{\bar H^{(\lambdaL)}(s)}</math>

being the last root estimate of the second stage and

<math display="block">s_{\lambda+1}=s_\lambda- \frac{P(s_\lambda)}{\bar H^{(\lambda+1)}(s_\lambda)}, \quad \lambda=L,L+1,\dots,</math>

Line 86 ⟶ 91:

It can be shown that, provided ''L'' is chosen sufficiently large, ''s''_λ always converges to a root of ''P''.

The algorithm converges for any distribution of roots, but may fail to find all roots of the polynomial. Furthermore, the convergence is slightly faster than the [[Rate of convergence|quadratic convergence]] of the Newton–Raphson iterationmethod, however, it uses at least twiceone-and-half as many operations per step, two polynomial evaluations for Newton vs. three polynomial evaluations in the third stage.

==What gives the algorithm its power?==