Jenkins–Traub algorithm: Difference between revisions

Revision as of 15:12, 3 August 2009

Revision as of 16:25, 3 August 2009

Line 73:

=== 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

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

:<math>P(X)=(X-\alpha_1)\cdot P_1(X)</math>

Polynomialwith a root <math>\alpha_1\in\C</math> and <math>P_1(X)=P(X)/(X-\alpha_1)</math> the remaining factor of degree ''n-1'' as the eigenvector equation for the multiplication bywith the variable ''X'', followed by remainder computation with divisor ''P(X)'',

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

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

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

which implies that <math>(X-\xialpha)\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 polynompolynomial ''P'', as

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

hasthe asresulting coefficient matrix is

:<math>A=\begin{pmatrix}

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

Line 89:

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

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

AlthoughTo itthis matrix the [[inverse power iteration]] is applied in three variants 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.

By''Stage using1:'' Here the [[inverse power iteration]] without shift

:<math>H^{k+1}=M_X{}^{-1}(H^k)</math> resp. <math>\xi_kbar H^{k+1}=q_k\cdot M_X{}^{-1}(H_k\bar H^k)</math>,

oneis used. In the second variant <math>\bar H^k</math> is a normalized multiple of <math>H^k</math>, with the normalization fixing the norm of the vector or some coordinate of it to be 1. This variant 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. This only works if the there is only one eigenvalue with a minimal distance to zero.

Using''Stage 2:'' This step uses the shiftedinverse power iteration in the variant with a constant shift ''s''. In consequnce, thatit is computing the smallest eigenvalue of <math>M_X-s\cdot id</math> by

:<math>H^{k+1}=(\xi_kM_X-s\cdot id)^{-1}(H^k)\,,</math> resp. <math>\bar H^{k+1}=q_k\cdot (M_X-s\cdot id)^{-1}(H_k)\,,bar H^k)</math>,

one may speed up the convergence. If the shift ''s'' is close to the smallest eigenvalue <math>\alpha_1</math>, then the smallest eigenvalue of the shifted matrix is <math>\alpha_1-s</math> which is close to zero,. whichThe mayconvergence substantiallyhas decreasean theorder convergenceof factor.

:<math>\alpha_1-s-q_k=O\left(\left(\tfrac{|\alpha_1-s|}{|\alpha_2-s|}\right)^k\right)</math>

which is the faster the smaller <math>|\alpha_1-s|<\!<|\alpha_2-s|=\min{}_{j=2,\dots,n}|\alpha_j-s|</math> is.

''Stage 3:'' Here the shift is updated every step using

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

to obtain an update to <math>s_k</math> from the growth factor between <math>H^{k+1}</math> and <math>H_1^{k+1}</math>,

:<math>s-s_{k+1}=s_k+\tfrac{PL(sH^{k+1})}{\bar HL(H_1^{k(s_k+1})}</math>

with some linear functional ''L''. If <math>\delta_k</math> denotes the distance of the normalized iteration vector <math>\bar H^k</math> from the eigenspace of <math>\alpha_1</math>, then one gets the for the convergence speed the recursions

:<math>\alpha_1-s_{k+1}=O(\delta_k\,|\alpha_1-s_k|^2)</math> and <math>\delta_{k+1}=O(\delta_k\cdot |\alpha_1-s_k|)</math>

resulting in an asymptotic convergence order of <math>\phi^2=1+\phi\approx 2.61</math>, where <math>\phi=\tfrac{1+\sqrt{5}}2</math> is the [[golden ratio]].

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-ss_k)\cdot H^{k+1}(X) +Q(X)\cdot P(X))\,.</math>

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

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

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

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

</math>

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

:<math>LC(H^{k+1})=-Q(s_k)

=-\frac{H^k(s_k)}{P(s_k)}

=-LC(H^k)\,\tfrac{\bar H^k(s_k)}{P(s_k)}\,,

</math>

such that the iteration for the normalized ''H'' polynomials reads as

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

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

=\frac1{X-s_k}\,\left(P(X)-\tfrac{P(s_k)}{\bar H^k(s_k)}\,\bar H^k(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_k)}</math>

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

In the additional evaluation of the third stage one only needs the fraction of the leading coefficients,

Chosing instead

:<math>s-\tfrac{PLC(s)}{\bar HH_1^{k+1}(s_kX))}</math>

=-\tfrac{H^{k+1}(s_k)}{P(s_k)}

amounts to a second inverse power iteration with shift <math>s_k</math>

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

from where the update formula

:<math>s_{k+1}

=s_k+\tfrac{LC(H^{k+1}(X))}{LC(H_1^{k+1}(X))}

=s_k-\tfrac{P(s_k)}{\bar H^{k+1}(s_k)}

</math> results.

==Real coefficients==