Jenkins – Traub-algoritm - Jenkins–Traub algorithm

Den Jenkins-Traub algoritm för polynoma nollor är ett snabbt globalt konvergent iterativ polynom rot undersökningsmetod som publicerades i 1970 av Michael A. Jenkins och Joseph F Traub . De gav två varianter, en för allmänna polynomer med komplexa koefficienter, allmänt känd som "CPOLY" -algoritmen, och en mer komplicerad variant för specialfallet av polynomer med verkliga koefficienter, allmänt känd som "RPOLY" -algoritmen. Det senare är "praktiskt taget en standard i black-box polynom root-finders".

Den här artikeln beskriver den komplexa varianten. Med tanke på ett polynom P ,

${\ displaystyle P (z) = \ sum _ {i = 0} ^ {n} a_ {i} z ^ {ni}, \ quad a_ {0} = 1, \ quad a_ {n} \ neq 0}$

med komplexa koefficienter den beräknar approximationer till de n nollställen av P ( z ), en i taget i stort sett ökande storleksordning. När varje rot har beräknats tas dess linjära faktor bort från polynom. Att använda denna deflation garanterar att varje rot beräknas bara en gång och att alla rötter hittas. ${\ displaystyle \ alpha _ {1}, \ alpha _ {2}, \ dots, \ alpha _ {n}}$

Den verkliga varianten följer samma mönster, men beräknar två rötter åt gången, antingen två verkliga rötter eller ett par konjugerade komplexa rötter. Genom att undvika komplex aritmetik kan den verkliga varianten vara snabbare (med en faktor 4) än den komplexa varianten. Jenkins – Traub-algoritmen har stimulerat avsevärd forskning om teori och programvara för metoder av denna typ.

Översikt

Jenkins – Traub-algoritmen beräknar alla rötter till ett polynom med komplexa koefficienter. Algoritmen börjar med att kontrollera polynomet för förekomst av mycket stora eller mycket små rötter. Vid behov omskalas koefficienterna genom en omskalning av variabeln. I algoritmen finns rätt rötter en efter en och i allmänhet i ökande storlek. Efter att varje rot har hittats deflateras polynom genom att dela upp motsvarande linjär faktor. Faktum är att faktoriseringen av polynomet till den linjära faktorn och det återstående deflaterade polynomet redan är ett resultat av rotfyndningsförfarandet. Rotfyndningsförfarandet har tre steg som motsvarar olika varianter av den inversa kraftgeneration . Se Jenkins och Traub . En beskrivning finns också i Ralston och Rabinowitz s. 383. Algoritmen liknar i andan den tvåstegsalgoritm som Traub studerat.

Rotfyndningsförfarande

Från och med det nuvarande polynomet P ( X ) av grad n beräknas den minsta roten av P (x) . För detta ändamål konstrueras en sekvens av så kallade H- polynomer. Dessa polynom är alla av grad n  - 1 och ska konvergera till faktorn P ( X ) som innehåller alla återstående rötter. Sekvensen av H- polynomer förekommer i två varianter, en onormaliserad variant som möjliggör enkla teoretiska insikter och en normaliserad variant av polynomer som håller koefficienterna inom ett numeriskt förnuftigt intervall. ${\ displaystyle {\ bar {H}}}$

Konstruktionen av H- polynomerna beror på en sekvens av komplexa nummer som kallas skift. Dessa skift beror, åtminstone i det tredje steget, av de tidigare H- polynomema. De H polynom definieras som lösningen på den implicita rekursion ${\ displaystyle \ left (H ^ {(\ lambda)} (z) \ right) _ {\ lambda = 0,1,2, \ dots}}$${\ displaystyle (s _ {\ lambda}) _ {\ lambda = 0,1,2, \ dots}}$

${\ displaystyle H ^ {(0)} (z) = P ^ {\ prime} (z)}$ och ${\ displaystyle (X-s _ {\ lambda}) \ cdot H ^ {(\ lambda +1)} (X) \ equiv H ^ {(\ lambda)} (X) {\ pmod {P (X)}} \.}$

En direkt lösning på denna implicita ekvation är

${\ 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) \ höger) \ ,,}$

där polynomuppdelningen är exakt.

Algoritmiskt skulle man till exempel använda Horner-schemat eller Ruffini-regeln för att utvärdera polynomerna vid och erhålla kvoterna samtidigt. Med de resulterande kvotienterna p ( X ) och h ( X ) som mellanresultat erhålls nästa H- polynom som ${\ displaystyle s _ {\ lambda}}$

${\ displaystyle \ left. {\ begin {align} P (X) & = p (X) \ cdot (X-s _ {\ lambda}) + P (s _ {\ lambda}) \\ H ^ {(\ lambda )} (X) & = h (X) \ cdot (X-s _ {\ lambda}) + H ^ {(\ lambda)} (s _ {\ lambda}) \\\ slut {justerad}} \ höger \} \ antyder H ^ {(\ lambda +1)} (z) = h (z) - {\ frac {H ^ {(\ lambda)} (s _ {\ lambda})} {P (s _ {\ lambda}) }} p (z).}$

Eftersom den högsta graden koefficienten erhålls från P (X) , den ledande koefficienten är . Om detta delas ut är det normaliserade H- polynomet ${\ displaystyle H ^ {(\ lambda +1)} (X)}$${\ displaystyle - {\ tfrac {H ^ {(\ lambda)} (s _ {\ lambda})} {P (s _ {\ lambda})}}}$

${\ displaystyle {\ begin {align} {\ 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) \ höger ) \\ [1em] & = {\ frac {1} {X-s _ {\ lambda}}} \ cdot \ left (P (X) - {\ frac {P (s _ {\ lambda})} {{\ stapel {H}} ^ {(\ lambda)} (s _ {\ lambda})}} {\ bar {H}} ^ {(\ lambda)} (X) \ höger) \,. \ slut {justerad}} }$

Steg ett: ingen skiftprocess

För uppsättning . Vanligtvis väljs M = 5 för polynom med måttliga grader upp till n  = 50. Detta steg är inte nödvändigt från enbart teoretiska överväganden, men är användbart i praktiken. Det betonar i H- polynomerna kofaktorn (av den linjära faktorn) för den minsta roten. ${\ displaystyle \ lambda = 0,1, \ prickar, M-1}$${\ displaystyle s _ {\ lambda} = 0}$

Steg två: fast skiftprocess

Skiftet för detta steg bestäms som någon punkt nära polynomets minsta rot. Det är kvasi-slumpmässigt placerat på cirkeln med den inre rotradien, som i sin tur uppskattas som den positiva lösningen av ekvationen

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

Eftersom vänster sida är en konvex funktion och ökar monotont från noll till oändlighet är denna ekvation lätt att lösa, till exempel med Newtons metod .

Välj nu på cirkeln för denna radie. Sekvensen av polynom , , genereras med den fasta skiftvärdet . Under denna iteration, den nuvarande approximationen för roten ${\ displaystyle s = R \ cdot \ exp (i \, \ phi _ {\ text {random}})}$${\ displaystyle H ^ {(\ lambda +1)} (z)}$${\ displaystyle \ lambda = M, M + 1, \ prickar, L-1}$${\ displaystyle s _ {\ lambda} = s}$

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

spåras. Den andra etappen är klar framgångsrikt om villkoren

${\ displaystyle | t _ {\ lambda +1} -t _ {\ lambda} | <{\ tfrac {1} {2}} \, | t _ {\ lambda} |}$ och ${\ displaystyle | t _ {\ lambda} -t _ {\ lambda -1} | <{\ tfrac {1} {2}} \, | t _ {\ lambda -1} |}$

möts samtidigt. Om det inte blev någon framgång efter ett visst antal iterationer prövas en annan slumpmässig punkt på cirkeln. Vanligtvis använder man ett antal nio iterationer för polynom av måttlig grad, med en fördubblingsstrategi för flera fel.

Steg tre: variabel skiftprocess

De genereras nu med de variabla skift som genereras av ${\ displaystyle H ^ {(\ lambda +1)} (X)}$${\ displaystyle s _ {\ lambda}, \ quad \ lambda = L, L + 1, \ dots}$

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

är den sista rotuppskattningen av andra etappen och

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

var är det normaliserade H- polynomet, som divideras med dess ledande koefficient. ${\ displaystyle {\ bar {H}} ^ {(\ lambda +1)} (z)}$${\ displaystyle H ^ {(\ lambda)} (z)}$

Om stegstorleken i steg tre inte faller tillräckligt snabbt till noll startas steg två igen med en annan slumpmässig punkt. Om detta inte lyckas efter ett litet antal omstarter fördubblas antalet steg i steg två.

Konvergens

Det kan visas att, förutsatt L väljs tillräckligt stor, s _λ alltid konvergerar till en rot till P .

Algoritmen konvergerar för någon distribution av rötter, men kan misslyckas med att hitta alla polynomets rötter. Dessutom är konvergensen något snabbare än den kvadratiska konvergensen av Newton – Raphson-iteration, men den använder dock minst dubbelt så många operationer per steg.

Vad ger algoritmen sin kraft?

Jämför med Newton – Raphson-iterationen

${\ displaystyle z_ {i + 1} = z_ {i} - {\ frac {P (z_ {i})} {P ^ {\ prime} (z_ {i})}}.}$

Iterationen använder den givna P och . Däremot den tredje etappen av Jenkins – Traub ${\ displaystyle \ scriptstyle P ^ {\ prime}}$

${\ displaystyle s _ {\ lambda +1} = s _ {\ lambda} - {\ frac {P (s _ {\ lambda})} {{\ bar {H}} ^ {\ lambda +1} (s _ {\ lambda })}} = s _ {\ lambda} - {\ frac {W ^ {\ lambda} (s _ {\ lambda})} {(W ^ {\ lambda}) '(s _ {\ lambda})}}}$

är just en Newton – Raphson-iteration utförd på vissa rationella funktioner . Mer exakt utförs Newton – Raphson på en sekvens av rationella funktioner

${\ displaystyle W ^ {\ lambda} (z) = {\ frac {P (z)} {H ^ {\ lambda} (z)}}.}$

För tillräckligt stora ${\ displaystyle \ lambda}$

${\ displaystyle {\ frac {P (z)} {{\ bar {H}} ^ {\ lambda} (z)}} = W ^ {\ lambda} (z) \, LC (H ^ {\ lambda} )}$

är så nära som önskat till ett första grads polynom

${\ displaystyle z- \ alpha _ {1}, \,}$

var är en av nollorna till . Även om steg 3 exakt är en Newton – Raphson-iteration utförs inte differentiering. ${\ displaystyle \ alpha _ {1}}$${\ displaystyle P}$

Analys av H- polynomema

Låt vara rötterna till P ( X ). De så kallade Lagrange-faktorerna hos P (X) är kofaktorerna för dessa rötter, ${\ displaystyle \ alpha _ {1}, \ dots, \ alpha _ {n}}$

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

Om alla rötter är olika, så bildar Lagrange-faktorerna en grund för rymden hos polynomier av grad högst n  - 1. Genom analys av rekursionsproceduren finner man att H- polynomema har koordinatrepresentationen

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

Varje Lagrange-faktor har ledande koefficient 1, så att den ledande koefficienten för H-polynomema är summan av koefficienterna. De normaliserade H-polynomema är således

${\ displaystyle {\ bar {H}} ^ {(\ lambda)} (X) = {\ frac {\ sum _ {m = 1} ^ {n} \ left [\ prod _ {\ kappa = 0} ^ {\ lambda -1} (\ alpha _ {m} -s _ {\ kappa}) \ höger] ^ {- 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}}} \ höger] \, P_ {m} (X)} {1+ \ sum _ {m = 1} ^ {n} \ vänster [\ prod _ {\ kappa = 0} ^ {\ lambda -1} {\ frac {\ alpha _ {1} -s _ {\ kappa}} {\ alpha _ {m} -s _ {\ kappa}}} \ rätt]}}\ .}$

Konvergensbeställningar

Om tillståndet gäller för nästan allt iterat, kommer de normaliserade H-polynomerna att konvergera åtminstone geometriskt mot . ${\ displaystyle | \ alpha _ {1} -s _ {\ kappa} | <\ min {} _ {m = 2,3, \ dots, n} | \ alpha _ {m} -s _ {\ kappa} |}$${\ displaystyle P_ {1} (X)}$

Under förutsättning att

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

man får de asymptotiska uppskattningarna för

• steg 1:

  ${\ displaystyle H ^ {(\ lambda)} (X) = P_ {1} (X) + O \ left (\ left | {\ frac {\ alpha _ {1}} {\ alpha _ {2}}} \ höger | ^ {\ lambda} \ höger).}$

• för steg 2, om s är tillräckligt nära för att : ${\ displaystyle \ alpha _ {1}}$

  ${\ displaystyle H ^ {(\ lambda)} (X) = P_ {1} (X) + O \ left (\ left | {\ frac {\ alpha _ {1}} {\ alpha _ {2}}} \ höger | ^ {M} \ cdot \ vänster | {\ frac {\ alpha _ {1} -s} {\ alpha _ {2} -s}} \ höger | ^ {\ lambda -M} \ höger)}$

  och

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

• och för steg 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}}} \ höger | \ höger)}$

  och

  ${\ 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}}} \ höger | \ cdot {\ frac {| \ alpha _ {1} -s _ {\ lambda} | ^ {2}} {| \ alpha _ {2 } -s _ {\ lambda} |}} \ höger)}$

ger upphov till en högre än kvadratisk konvergensordning av , var är det gyllene förhållandet . ${\ displaystyle \ phi ^ {2} = 1 + \ phi \ ca 2.61}$${\ displaystyle \ phi = {\ tfrac {1} {2}} (1 + {\ sqrt {5}})}$

Tolkning som omvänd kraft iteration

Alla steg i Jenkins – Traub-komplexalgoritmen kan representeras som det linjära algebra-problemet för att bestämma egenvärdena för en speciell matris. Denna matris är koordinatrepresentationen för en linjär karta i det n -dimensionella utrymmet för polynom med grad n  - 1 eller mindre. Huvudidén med denna karta är att tolka faktoriseringen

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

med en rot och den återstående faktorn av grad n  - 1 som egenvektorekvationen för multiplikationen med variabeln X , följt av återstående beräkning med divisor P ( X ), ${\ displaystyle \ alpha _ {1} \ in \ mathbb {C}}$${\ displaystyle P_ {1} (X) = P (X) / (X- \ alpha _ {1})}$

${\ displaystyle M_ {X} (H) = (X \ cdot H (X)) {\ bmod {P}} (X) \ ,.}$

Detta kartlägger högst grad polynomier n  - 1 till högst polynomier n  - 1. Egenvärdena för denna karta är rötterna till P ( X ), eftersom egenvektorekvationen lyder

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

vilket antyder att det vill säga är en linjär faktor för P ( X ). På monombasis representeras den linjära kartan av en kompletterande matris av polynom P , som ${\ displaystyle (X- \ alpha) \ cdot H) = C \ cdot P (X)}$${\ displaystyle (X- \ alpha)}$${\ displaystyle M_ {X}}$

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

den resulterande koefficientmatrisen är

${\ displaystyle A = {\ begin {pmatrix} 0 & 0 & \ dots & 0 & -P_ {0} \\ 1 & 0 & \ dots & 0 & -P_ {1} \\ 0 & 1 & \ dots & 0 & -P_ {2} \\\ vdots & \ vdots & \ ddots & \ vdots & \ vdots \\ 0 & 0 & \ dots & 1 & -P_ {n-1} \ end {pmatrix}} \ ,.}$

Till denna matris tillämpas den inversa effekt iterationen i de tre varianterna av inget skift, konstant skift och generaliserad Rayleigh skift i algoritmens tre steg. Det är mer effektivt att utföra de linjära algebraoperationerna i polynomiell aritmetik och inte genom matrisoperationer, men egenskaperna hos den inversa kraft iterationen förblir desamma.

Verkliga koefficienter

Den tidigare beskrivna Jenkins – Traub-algoritmen fungerar för polynom med komplexa koefficienter. Samma författare skapade också en trestegsalgoritm för polynom med verkliga koefficienter. Se Jenkins och Traub A Three-Stage Algorithm for Real Polynomials Using Quadratic Iteration . Algoritmen hittar antingen en linjär eller kvadratisk faktor som fungerar helt i verklig aritmetik. Om de komplexa och verkliga algoritmerna tillämpas på samma riktiga polynom är den verkliga algoritmen ungefär fyra gånger så snabb. Den verkliga algoritmen konvergerar alltid och konvergenshastigheten är större än andra ordningen.

En anslutning med den skiftade QR-algoritmen

Det finns en överraskande koppling till den skiftade QR-algoritmen för beräkning av matrisens egenvärden. Se Dekker och Traub Den skiftade QR-algoritmen för Hermitian-matriser . Återigen kan skiften ses som Newton-Raphson-iteration på en sekvens av rationella funktioner som konvergerar till ett första grads polynom.

Programvara och testning

Programvaran för Jenkins – Traub-algoritmen publicerades som Jenkins och Traub Algorithm 419: Zeros of a Complex Polynomial . Programvaran för den verkliga algoritmen publicerades som Jenkins Algorithm 493: Zeros of a Real Polynomial .

Metoderna har testats utförligt av många människor. Som förutsagt njuter de snabbare än kvadratisk konvergens för alla nollfördelningar.

Det finns emellertid polynomer som kan orsaka förlust av precision, såsom illustreras av följande exempel. Polynomet har alla sina nollor som ligger på två halvcirklar med olika radier. Wilkinson rekommenderar att det är önskvärt för stabil deflation att mindre nollor beräknas först. Skift i andra steget väljs så att nollorna på den mindre halvcirkeln hittas först. Efter deflation är polynom med nollor på halvcirkeln känt för att vara dåligt konditionerat om graden är stor; se Wilkinson, s. 64. Det ursprungliga polynomet var av grad 60 och drabbades av allvarlig deflationsinstabilitet.

Referenser

externa länkar

• En gratis nedladdningsbar Windows-applikation med Jenkins – Traub-metoden för polynom med verkliga och komplexa koefficienter
• RPoly ++ En SSE-optimerad C ++ - implementering av RPOLY-algoritmen.