Root-finding procedure

Starting with the current polynomial P(X) of degree n, the aim is to compute the smallest root ${\displaystyle \alpha }$ of P(x). The polynomial can then be split into a linear factor and the remaining polynomial factor

$${\displaystyle P(X)=(X-\alpha ){\bar {H}}(X)}$$

Other root-finding methods strive primarily to improve the root and thus the first factor. The main idea of the Jenkins-Traub method is to incrementally improve the second factor.

To that end, a sequence of so-called H polynomials is constructed. These polynomials are all of degree n − 1 and are supposed to converge to the factor ${\displaystyle {\bar {H}}(X)}$ 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 ${\displaystyle {\bar {H}}}$ polynomials that keeps the coefficients in a numerically sensible range. The construction of the H polynomials ${\displaystyle \left(H^{(\lambda )}(z)\right)_{\lambda =0,1,2,\dots }}$ is guided by a sequence of complex numbers ${\displaystyle (s_{\lambda })_{\lambda =0,1,2,\dots }}$ 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

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

and

$${\displaystyle (X-s_{\lambda })\cdot H^{(\lambda +1)}(X)\equiv H^{(\lambda )}(X){\pmod {P(X)}}\ .}$$

A direct solution to this implicit equation is

$${\displaystyle H^{(\lambda +1)}(X)={\frac {1}{X-s_{\lambda }}}\cdot \left(H^{(\lambda )}(X)-{\frac {H^{(\lambda )}(s_{\lambda })}{P(s_{\lambda })}}P(X)\right)\,,}$$

where the polynomial division is exact.

Algorithmically, one would use long division by the linear factor as in the Horner scheme or Ruffini rule to evaluate the polynomials at ${\displaystyle s_{\lambda }}$ 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

$${\displaystyle \left.{\begin{aligned}P(X)&=p(X)\cdot (X-s_{\lambda })+P(s_{\lambda })\\H^{(\lambda )}(X)&=h(X)\cdot (X-s_{\lambda })+H^{(\lambda )}(s_{\lambda })\\\end{aligned}}\right\}\implies H^{(\lambda +1)}(z)=h(z)-{\frac {H^{(\lambda )}(s_{\lambda })}{P(s_{\lambda })}}p(z).}$$

Since the highest degree coefficient is obtained from P(X), the leading coefficient of ${\displaystyle H^{(\lambda +1)}(X)}$ is ${\displaystyle -{\tfrac {H^{(\lambda )}(s_{\lambda })}{P(s_{\lambda })}}}$. If this is divided out the normalized H polynomial is

$${\displaystyle {\begin{aligned}{\bar {H}}^{(\lambda +1)}(X)&={\frac {1}{X-s_{\lambda }}}\cdot \left(P(X)-{\frac {P(s_{\lambda })}{H^{(\lambda )}(s_{\lambda })}}H^{(\lambda )}(X)\right)\\[1em]&={\frac {1}{X-s_{\lambda }}}\cdot \left(P(X)-{\frac {P(s_{\lambda })}{{\bar {H}}^{(\lambda )}(s_{\lambda })}}{\bar {H}}^{(\lambda )}(X)\right)\,.\end{aligned}}}$$

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

$${\displaystyle R^{n}+|a_{n-1}|\,R^{n-1}+\dots +|a_{1}|\,R=|a_{0}|\,.}$$

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 ${\displaystyle s=R\cdot \exp(i\,\phi _{\text{random}})}$ on the circle of this radius. The sequence of polynomials ${\displaystyle H^{(\lambda +1)}(z)}$, ${\displaystyle \lambda =M,M+1,\dots ,L-1}$, is generated with the fixed shift value ${\displaystyle s_{\lambda }=s}$. 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

$${\displaystyle t_{\lambda }=s-{\frac {P(s)}{{\bar {H}}^{(\lambda )}(s)}}}$$

is traced. The second stage is terminated as successful if the conditions

$${\displaystyle |t_{\lambda +1}-t_{\lambda }|<{\tfrac {1}{2}}\,|t_{\lambda }|}$$

and

$${\displaystyle |t_{\lambda }-t_{\lambda -1}|<{\tfrac {1}{2}}\,|t_{\lambda -1}|}$$

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 ${\displaystyle H^{(\lambda +1)}(X)}$ polynomials are now generated using the variable shifts ${\displaystyle s_{\lambda },\quad \lambda =L,L+1,\dots }$ which are generated by

$${\displaystyle s_{L}=t_{L}=s-{\frac {P(s)}{{\bar {H}}^{(L)}(s)}}}$$

being the last root estimate of the second stage and

$${\displaystyle s_{\lambda +1}=s_{\lambda }-{\frac {P(s_{\lambda })}{{\bar {H}}^{(\lambda +1)}(s_{\lambda })}},\quad \lambda =L,L+1,\dots ,}$$

where ${\displaystyle {\bar {H}}^{(\lambda +1)}(z)}$ is the normalized H polynomial, that is ${\displaystyle H^{(\lambda )}(z)}$ divided by its leading coefficient.

If the step size in stage three does not fall fast enough to zero, then stage two is restarted using a different random point. If this does not succeed after a small number of restarts, the number of steps in stage two is doubled.

Analysis of the H polynomials

Let ${\displaystyle \alpha _{1},\dots ,\alpha _{n}}$ be the roots of P(X). The so-called Lagrange factors of P(X) are the cofactors of these roots,

$${\displaystyle P_{m}(X)={\frac {P(X)-P(\alpha _{m})}{X-\alpha _{m}}}.}$$

If all roots are different, then the Lagrange factors form a basis of the space of polynomials of degree at most n − 1. By analysis of the recursion procedure one finds that the H polynomials have the coordinate representation

$${\displaystyle H^{(\lambda )}(X)=\sum _{m=1}^{n}\left[\prod _{\kappa =0}^{\lambda -1}(\alpha _{m}-s_{\kappa })\right]^{-1}\,P_{m}(X)\ .}$$

Each Lagrange factor has leading coefficient 1, so that the leading coefficient of the H polynomials is the sum of the coefficients. The normalized H polynomials are thus

$${\displaystyle {\bar {H}}^{(\lambda )}(X)={\frac {\sum _{m=1}^{n}\left[\prod _{\kappa =0}^{\lambda -1}(\alpha _{m}-s_{\kappa })\right]^{-1}\,P_{m}(X)}{\sum _{m=1}^{n}\left[\prod _{\kappa =0}^{\lambda -1}(\alpha _{m}-s_{\kappa })\right]^{-1}}}={\frac {P_{1}(X)+\sum _{m=2}^{n}\left[\prod _{\kappa =0}^{\lambda -1}{\frac {\alpha _{1}-s_{\kappa }}{\alpha _{m}-s_{\kappa }}}\right]\,P_{m}(X)}{1+\sum _{m=1}^{n}\left[\prod _{\kappa =0}^{\lambda -1}{\frac {\alpha _{1}-s_{\kappa }}{\alpha _{m}-s_{\kappa }}}\right]}}\ .}$$

Convergence orders

If the condition ${\displaystyle |\alpha _{1}-s_{\kappa }|<\min {}_{m=2,3,\dots ,n}|\alpha _{m}-s_{\kappa }|}$ holds for almost all iterates, the normalized H polynomials will converge at least geometrically towards ${\displaystyle P_{1}(X)}$.

Under the condition that

$${\displaystyle |\alpha _{1}|<|\alpha _{2}|=\min {}_{m=2,3,\dots ,n}|\alpha _{m}|}$$

one gets the asymptotic estimates for

• stage 1:
  $${\displaystyle H^{(\lambda )}(X)=P_{1}(X)+O\left(\left|{\frac {\alpha _{1}}{\alpha _{2}}}\right|^{\lambda }\right).}$$

  
• for stage 2, if s is close enough to ${\displaystyle \alpha _{1}}$:
  $${\displaystyle H^{(\lambda )}(X)=P_{1}(X)+O\left(\left|{\frac {\alpha _{1}}{\alpha _{2}}}\right|^{M}\cdot \left|{\frac {\alpha _{1}-s}{\alpha _{2}-s}}\right|^{\lambda -M}\right)}$$

  
  and
  $${\displaystyle s-{\frac {P(s)}{{\bar {H}}^{(\lambda )}(s)}}=\alpha _{1}+O\left(\ldots \cdot |\alpha _{1}-s|\right).}$$

  
• and for stage 3:
  $${\displaystyle H^{(\lambda )}(X)=P_{1}(X)+O\left(\prod _{\kappa =0}^{\lambda -1}\left|{\frac {\alpha _{1}-s_{\kappa }}{\alpha _{2}-s_{\kappa }}}\right|\right)}$$

  
  and
  $${\displaystyle s_{\lambda +1}=s_{\lambda }-{\frac {P(s)}{{\bar {H}}^{(\lambda +1)}(s_{\lambda })}}=\alpha _{1}+O\left(\prod _{\kappa =0}^{\lambda -1}\left|{\frac {\alpha _{1}-s_{\kappa }}{\alpha _{2}-s_{\kappa }}}\right|\cdot {\frac {|\alpha _{1}-s_{\lambda }|^{2}}{|\alpha _{2}-s_{\lambda }|}}\right)}$$

  
  giving rise to a higher than quadratic convergence order of ${\displaystyle \phi ^{2}=1+\phi \approx 2.61}$, where ${\displaystyle \phi ={\tfrac {1}{2}}(1+{\sqrt {5}})}$ 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 factorization

$${\displaystyle P(X)=(X-\alpha _{1})\cdot P_{1}(X)}$$

with a root ${\displaystyle \alpha _{1}\in \mathbb {C} }$ and ${\displaystyle P_{1}(X)=P(X)/(X-\alpha _{1})}$ 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),

$${\displaystyle 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

$${\displaystyle 0=(M_{X}-\alpha \cdot id)(H)=((X-\alpha )\cdot H){\bmod {P}}\,,}$$

which implies that ${\displaystyle (X-\alpha )\cdot H)=C\cdot P(X)}$, that is, ${\displaystyle (X-\alpha )}$ is a linear factor of P(X). In the monomial basis the linear map ${\displaystyle M_{X}}$ is represented by a companion matrix of the polynomial P, as

$${\displaystyle M_{X}(H)=\sum _{m=0}^{n-1}H_{m}X^{m+1}-H_{n-1}\left(X^{n}+\sum _{m=0}^{n-1}a_{m}X^{m}\right)=\sum _{m=1}^{n-1}(H_{m-1}-a_{m}H_{n-1})X^{m}-a_{0}H_{n-1}\,,}$$

the resulting transformation matrix is

$${\displaystyle A={\begin{pmatrix}0&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_{n-1}\end{pmatrix}}\,.}$$

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.