Jenkins–Traub algorithm: Difference between revisions

Revision as of 09:04, 27 July 2009

Revision as of 15:05, 27 July 2009

Line 75:

For the inverse of the derivative of this map one calculates, using <math>G(X)=X-s</math>

:<math>

(\nabla F(G,\bar H))^{-1}(Y)

=\left(

\frac{Y(s)}{\bar H(s)},\;

\frac1{X-s}\left(Y-\tfrac{Y(s)}{\bar H(s)}\bar H\right)

\right)\,.

</math>.

In the third stage, the full Newton update is applied, that is

:<math>\begin{align}

(G^{k+1},\; \bar H^{k+1})

&=(G^k,\; \bar H^k)-(\nabla F(GG_k,H\bar H_k))^{-1}(G^k\cdot \bar H^k-P)\\

&=\left(

X-s_k+\frac{P(ss_k)}{\bar H^k(ss_k)}, \;

\frac1{X-s_k}\left(P-\tfrac{P(ss_k)}{\bar H^k(ss_k)}\bar H^k\right)

\right)\,.

\end{align}</math>.

From this formula one recovers <math>s_{k+1}=s_k-\tfrac{P(ss_k)}{\bar H^k(ss_k)}</math> and

that <math>\bar H^{k+1}</math> is again a normed polynomial, since ''P'' has leading coefficient ''1'' and <math>\bar H^k</math> has one degree less, not influencing the leading coefficient of the quotient.

In the first and second stage, the shifts <math>s_k</math> and consequently the first polynomial <math>G^k</math> are left constant, that is without update. The ensuing method is best understood as an inverse power iteration.

=== As inverse power iteration ===

This describes the behaviour of the Jenkins-Traub algorithm in its first two stages.

Polynomial multiplication by ''X'' followed by remainder computation with divisor ''P(X)'',

:<math>M_X(H)=(X\cdot H(X) \mod P)\,,</math>

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

:<math>0=(M_X-\xi\cdot id)(H)=((X-\xi)\cdot H) \mod P)\,,</math>

which implies that <math>(X-\xi)</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 polynom ''P'',

:<math> M_X(H)=\sum_{m=1}^{n-1}(H_{m-1}-P_{m}H_{n-1})X^m-P_0H_{n-1}</math>

has as coefficient matrix

:<math>A=\begin{pmatrix}

0 & 0 & \dots & 0 & -P_0 \\

1 & 0 & \dots & 0 & -P_1 \\

0 & 1 & \dots & 0 & -P_2 \\

\vdots & \vdots & \vdots & \vdots & \vdots \\

0 & 0 & \dots & 1 & -P_{n-1} \\

\end{pmatrix}\,.</math>

Although it is more efficient to perform the linear operations in polynomial arithmetic and not by matrix operations, the properties of the inverse power iteration remain the same.

By using the [[inverse power iteration]]

:<math>H^{k+1}=\xi_k\cdot M_X{}^{-1}(H_k)</math>

one approximates the smallest eigenvalue, that is the smallest root of ''P(X)''. The convergence is of linear order with a factor that is the quotient of the smallest (in absolute value) and next-to-smallest eigenvalue.

Using the shifted variant, that is computing the smallest eigenvalue of <math>M_X-s\cdot id</math> by

:<math>H^{k+1}=(\xi_k-s)\cdot (M_X-s\cdot id)^{-1}(H_k)\,,</math>

one may speed up the convergence. If the shift ''s'' is close to the smallest eigenvalue, then the smallest eigenvalue of the shifted matrix is close to zero, which may substantially decrease the convergence factor.

For the inverse power iteration one has to solve the equation

:<math>H^k(X)=(M_X-s\cdot id)(H^{k+1})=((X-s)\cdot H^{k+1}(X) \bmod P(X))</math>

for <math>H^{k+1}</math>. If the smallest eigenvalue is simple, then the polynomials <math>H^k</math> will after a while be almost multiples of their predecessors, <math>H^{k+1}\approx \alpha_kH^k</math>. The approximation of the eigenvalue is then the reciprocal value, <math>(\xi_k-s)=\tfrac1{\alpha_k}</math>.

To solve the polynomial identity representing the inverse power iteration, one has to find a factor <math>Q(X)</math> such that

:<math>H^k(X)=((X-s)\cdot H^{k+1}(X) +Q(X)\cdot P(X))\,.</math>

Evaluating at <math>X=s</math> one arrives at <math>H^k(s)=Q(s)P(s)</math>, and by degree computation one concludes that ''Q(X)=Q(s)'' is a constant. Then the equation is solved by exact polynomial division giving

:<math>H^{k+1}(X)

=\frac1{X-s}\,(H^k(X)-Q(s)\,P(X))

=\frac1{X-s}\,\left(H^k(X)-\tfrac{H^k(s)}{P(s)}\,P(X)\right)\,,

</math>

which is the unnormalized iteration for the polynomial sequence <math>(H^k(s))_{k\in\N}</math>. Comparing the leading coefficients (LC) of this sequence one finds that

:<math>LC(H^{k+1})=-\frac{H^k(s)}{P(s)}=-LC(H^k)\,\frac{\bar H^k(s)}{P(s)}\,,</math>

such that

:<math>-\frac{P(s)}{\bar H^k(s)}</math>

is a current approximation of the smallest eigenvalue of <math>(M_X-s\cdot id)</math> resp.

:<math>s-\tfrac{P(s)}{\bar H^k(s)}</math>

is an improved approximation of the smallest root of ''P(X)''.

==Real coefficients==