Algorytm Jenkinsa-Trauba

Algorytm Jenkinsa-Trauba dla wielomianów zerowych to szybka globalnie zbieżna iteracyjna wielomianowa metoda wyszukiwania pierwiastków, opublikowana w 1970 roku przez Michaela A. Jenkinsa i Josepha F. Trauba . Podali dwa warianty, jeden dla ogólnych wielomianów o zespolonych współczynnikach, powszechnie znany jako algorytm „CPOLY”, oraz bardziej skomplikowany wariant dla szczególnego przypadku wielomianów o rzeczywistych współczynnikach, powszechnie znany jako algorytm „RPOLY”. To ostatnie jest „praktycznie standardem w wielomianowych wyszukiwarkach pierwiastków z czarną skrzynką”.

W tym artykule opisano złożony wariant. Biorąc pod uwagę wielomian P ,

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

ze złożonymi współczynnikami oblicza przybliżenia do n zer ${\ Displaystyle \ alfa _ {1}, \ alfa _ {2}, \ kropek, \ alfa _ {n}$ } P ( z ), po jednym na raz w mniej więcej rosnącym rzędzie wielkości. Po obliczeniu każdego pierwiastka jego czynnik liniowy jest usuwany z wielomianu. Użycie tej deflacji gwarantuje, że każdy pierwiastek zostanie obliczony tylko raz i że wszystkie pierwiastki zostaną znalezione.

Wariant rzeczywisty jest zgodny z tym samym wzorcem, ale jednocześnie oblicza dwa pierwiastki, albo dwa pierwiastki rzeczywiste, albo parę sprzężonych pierwiastków zespolonych. Unikając złożonej arytmetyki, wariant rzeczywisty może być szybszy (czterokrotnie) niż wariant złożony. Algorytm Jenkinsa-Trauba pobudził znaczne badania nad teorią i oprogramowaniem dla metod tego typu.

Przegląd

Algorytm Jenkinsa-Trauba oblicza wszystkie pierwiastki wielomianu ze złożonymi współczynnikami. Algorytm rozpoczyna się od sprawdzenia wielomianu pod kątem występowania bardzo dużych lub bardzo małych pierwiastków. W razie potrzeby współczynniki są przeskalowane poprzez przeskalowanie zmiennej. W algorytmie odpowiednie pierwiastki są znajdowane jeden po drugim i generalnie w coraz większym rozmiarze. Po znalezieniu każdego pierwiastka wielomian jest deflowany przez podzielenie odpowiedniego współczynnika liniowego. Rzeczywiście, faktoryzacja wielomianu na czynnik liniowy i pozostały wielomian deflowany jest już wynikiem procedury wyszukiwania pierwiastków. Procedura wyszukiwania korzeni składa się z trzech etapów, które odpowiadają różnym wariantom odwrotna iteracja mocy . Patrz Jenkins i Traub . Opis można znaleźć także u Ralstona i Rabinowitza, str. 383. Algorytm jest podobny w duchu do algorytmu dwuetapowego badanego przez Trauba.

Procedura wyszukiwania korzeni

Rozpoczynając od bieżącego wielomianu P ( X ) stopnia n , obliczany jest najmniejszy pierwiastek z P(x) . W tym celu konstruuje się ciąg tak zwanych wielomianów H. Wszystkie te wielomiany są stopnia n - 1 i mają zbiegać się do czynnika P ( X ) zawierającego wszystkie pozostałe pierwiastki. Sekwencja H występuje w dwóch wariantach, nieznormalizowanym wariancie, który umożliwia łatwy wgląd teoretyczny i znormalizowanym wariancie ${\ displaystyle {\ bar {H}}}$ wielomiany, które utrzymują współczynniki w zakresie liczbowym.

Konstrukcja wielomianów H. ${\ Displaystyle \ lewo (H ^ {(\ lambda)} (z) \ prawej) _ {\ lambda = 0, 1,2 \ kropki}}$ zależy od ciągu liczb zespolonych ${\ Displaystyle (s _ {\ lambda}) _ {\ lambda = 0,1,2, \ kropki }}$ zwane przesunięciami. Same te przesunięcia zależą, przynajmniej na trzecim etapie, od poprzednich wielomianów H. Wielomiany H są zdefiniowane jako rozwiązanie niejawnej rekurencji

${\ Displaystyle H ^ {(0)} (z) = P ^ {\ pierwsza} (z)}$ i ${\ Displaystyle (X-s _ {\ lambda}) \ cdot H ^ {(\ lambda +1)} (X) \ równoważnik H ^ {(\ lambda)} (X) {\ pmod {P (X)}} \ .}$

Bezpośrednim rozwiązaniem tego ukrytego równania jest

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

gdzie podział wielomianowy jest dokładny.

Algorytmicznie można by użyć na przykład $Hornera lub reguły Ruffiniego do oszacowania wielomianów w punkcie$ i uzyskania ilorazów w samym czasie. Z wynikowych ilorazów p ( X ) i h ( X ) jako wyników pośrednich otrzymujemy kolejny wielomian H jako

${\ Displaystyle \ lewo. {\ rozpocząć {wyrównane} 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{wyrównane}}\right\} \implikuje H^{(\lambda +1)}(z)=h(z)-{\frac {H^{(\lambda )}(s_{\lambda })}{P(s_{\lambda }) }}p(z).}$

Ponieważ współczynnik najwyższego stopnia uzyskuje się z P ( X ) , wiodącym współczynnikiem ${\ Displaystyle H ^ {(\ lambda + 1)} (X)}$ jest ${\ Displaystyle - {\ tfrac {H ^ {(\ lambda)} (s _ {\ lambda})} {P (s_ {\ lambda})}}}$ . Jeśli to jest podzielone, znormalizowany H wynosi

${\ Displaystyle {\ rozpocząć {wyrównane} {\ bar {H}} ^ {(\ lambda + 1)} (X) & = {\ Frac {1} {X-s _ {\ lambda}}} \ cdot \ lewo (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{wyrównane}} }$

Etap pierwszy: proces bez zmian

dla zestawu s ${\ Displaystyle \ lambda = 0,1, \ kropki, M-1}$$\ Displaystyle s_ {\ lambda} = 0}$ λ . Zwykle M=5 wybiera się dla wielomianów o umiarkowanych stopniach do n = 50. Ten etap nie jest konieczny z samych rozważań teoretycznych, ale jest przydatny w praktyce. Podkreśla w H kofaktor (czynnika liniowego) najmniejszego pierwiastka.

Etap drugi: proces na stałe

Przesunięcie dla tego etapu jest określane jako punkt bliski najmniejszemu pierwiastkowi wielomianu. Znajduje się quasi-losowo na okręgu o wewnętrznym promieniu pierwiastka, który z kolei jest szacowany jako dodatnie rozwiązanie równania

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

Ponieważ lewa strona jest funkcją wypukłą i rośnie monotonicznie od zera do nieskończoności, równanie to jest łatwe do rozwiązania, na przykład metodą Newtona .

Teraz wybierz ${\ Displaystyle s = R \ cdot \ exp (i \, \ phi _ {\ tekst {losowo}})}$ na okręgu o tym promieniu. Ciąg wielomianów ${\ Displaystyle H ^ {(\ lambda + 1)} (z)}$ , ${\ Displaystyle \ lambda = M, M + 1, \ kropki, L-1}$ , jest generowany ze stałą wartością przesunięcia ${\ Displaystyle s _ {\ lambda} = s}$ . Podczas tej iteracji bieżące przybliżenie pierwiastka

${\ Displaystyle t _ {\ lambda} = s - {\ Frac {P (s)} {{\ bar {H}} ^ {(\ lambda )}(S)}}}$

jest śledzony. Drugi etap kończy się pomyślnie, jeśli warunki

${\ Displaystyle | t _ {\ lambda +1} -t _ {\ lambda} | <{\ tfrac {1} {2}} \, | t_ {\ lambda} |}$ i ${\ Displaystyle | t _ {\ lambda} -t _ {\ lambda -1} | <{\ tfrac {1} {2}} \, | t_ {\ lambda -1} |}$

są jednocześnie spełnione. Jeśli po pewnej liczbie iteracji nie było powodzenia, próbowany jest inny losowy punkt na okręgu. Zazwyczaj stosuje się liczbę 9 iteracji dla wielomianów średniego stopnia, ze strategią podwojenia w przypadku wielu niepowodzeń.

Etap trzeci: proces zmiennej zmiany

H. ${\ Displaystyle H ^ {(\ lambda + 1)} (X)} są teraz generowane$ użyciu zmiennych przesunięć ${\ Displaystyle s_ {\lambda },\quad \lambda =L,L+1,\dots }$ które są generowane przez

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

będąc ostatnim oszacowaniem pierwiastkowym drugiego etapu i

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

gdzie ${\ Displaystyle {\ bar {H}} ^ {(\ lambda + 1)} (z)}$ jest znormalizowanym wielomianem H , czyli ${\ Displaystyle H ^ {( \lambda )}(z)}$ podzielone przez wiodący współczynnik.

Jeśli wielkość kroku w etapie trzecim nie spada wystarczająco szybko do zera, etap drugi jest uruchamiany ponownie przy użyciu innego losowego punktu. Jeśli to się nie powiedzie po niewielkiej liczbie ponownych uruchomień, liczba kroków w etapie drugim jest podwojona.

Konwergencja

Można pokazać, że pod warunkiem, że L jest wybrane wystarczająco duże, s _λ zawsze zbiega się do pierwiastka z P .

Algorytm jest zbieżny dla dowolnego rozkładu pierwiastków, ale może nie znaleźć wszystkich pierwiastków wielomianu. Co więcej, zbieżność jest nieco szybsza niż zbieżność kwadratowa iteracji Newtona-Raphsona, jednak wykorzystuje co najmniej dwa razy więcej operacji na krok.

Co daje algorytmowi jego moc?

Porównaj z iteracją Newtona-Raphsona

${\ Displaystyle z_ {i + 1} = z_ {i} - {\ Frac {P (z_ {i})} {P ^ {\ pierwsza} (z_ {i})}}.}$

Iteracja wykorzystuje podane P i ${\ displaystyle \ scriptstyle P ^ {\ prime}}$ . W przeciwieństwie do trzeciego etapu Jenkinsa-Trauba

${\ 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})}}}$

jest dokładnie iteracją Newtona-Raphsona wykonaną na pewnych funkcjach wymiernych . Dokładniej, Newton-Raphson jest wykonywany na sekwencji funkcji wymiernych

${\ Displaystyle W ^ {\ lambda} (z) = {\ Frac {P (z)} {H ^ {\ lambda} (z)}}.}$

Dla wystarczająco dużego, ${\ displaystyle \ lambda}$

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

jest jak najbardziej zbliżony do wielomianu pierwszego stopnia

${\ Displaystyle Z- \ alfa _ {1}, \,}$

gdzie jest jednym z zer ${1}}$$alpha$ . Mimo że etap 3 jest dokładnie iteracją Newtona-Raphsona, różnicowanie nie jest wykonywane.

Analiza wielomianów H

Niech ${\ Displaystyle \ alfa _ {1}, \ kropki, \ alfa _ {n}}$ będą pierwiastkami P ( X ). Tak zwane czynniki Lagrange'a P(X) są kofaktorami tych pierwiastków,

${\ Displaystyle P_ {m} (X) = {\ Frac {P (X) -P (\ alfa _ {m})} {X- \ alfa _ {m}}}.}$

Jeśli wszystkie pierwiastki są różne, to czynniki Lagrange'a tworzą podstawę przestrzeni wielomianów stopnia co najwyżej n - 1. Analizując procedurę rekurencji, stwierdza się, że wielomiany H mają reprezentację współrzędnych

${\ Displaystyle H ^ {(\ lambda)} (X) = \ suma _ {m = 1} ^ {n} \ lewo [\ prod _ {\ kappa = 0} ^ {\ lambda -1} (\ alfa _ {m}-s_{\kappa})\right]^{-1}\,P_{m}(X)\ .}$

Każdy czynnik Lagrange'a ma wiodący współczynnik 1, tak więc wiodący współczynnik wielomianów H jest sumą współczynników. Znormalizowane wielomiany H są zatem

${\ Displaystyle {\ bar {H}} ^ {(\ lambda)} (X) = {\ Frac {\ suma _ {m = 1} ^ {n} \ lewo [\ 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 }}}\ Prawidłowy]}}\ .}$

Zamówienia konwergencji

Jeżeli warunek ${\ Displaystyle |\ alfa _ {1} -s _ {\ kappa} | <\ min {} _ {m = 2,3, \ kropki, n} | \ alfa _ {m} -s_ {\ kappa} |}$ zachodzi dla prawie wszystkich iteracji, znormalizowane wielomiany H będą zbiegać się przynajmniej geometrycznie w kierunku ${\ Displaystyle P_ {1} (X)}$ .

Pod warunkiem, że

${\ Displaystyle | \ alfa _ {1} | <| \ alfa _ {2} | = \ min {} _ {m = 2,3, \ kropki, n} | \ alfa _ {m} |}$

otrzymuje się asymptotyczne oszacowania dla

• etap 1:

  ${\ Displaystyle H ^ {(\ lambda)} (X) = P_ {1} (X) + O \ lewo (\ lewo | {\ Frac {\ alfa _ {1}} {\ alfa _ {2}}} \right|^{\lambda }\right).}$

• dla etapu 2, jeśli s jest wystarczająco blisko ${\ Displaystyle \ alpha _ {1}}$ :

  ${\ Displaystyle H ^ {(\ lambda)} (X) = P_ {1} (X) + O \ lewo (\ lewo | {\ Frac {\ alfa _ {1}} {\ alfa _ {2}}} \right|^{M}\cdot \left|{\frac {\alpha _{1}-s}{\alpha _{2}-s}}\right|^{\lambda -M}\right)}$

  i

  ${\ Displaystyle s - {\ Frac {P (s)} {{\ bar {H}} ^ {(\ lambda)} (s)}} = \ alfa _ {1} + O \ lewo (\ ldots \ cdot |\alpha _{1}-s|\right).}$

• i dla etapu 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)}$

  i

  ${\ 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}|}}\prawo)}$

dając wyższy niż kwadratowy rząd zbieżności , gdzie ${\ Displaystyle \ phi ^ {2} = 1 + \ phi \ około 2,61}$${\ Displaystyle \phi ={\tfrac {1}{2}}(1+{\sqrt {5}})}$ to złoty podział .

Interpretacja jako odwrotna iteracja potęgowa

Wszystkie etapy złożonego algorytmu Jenkinsa-Trauba można przedstawić jako problem algebry liniowej wyznaczania wartości własnych specjalnej macierzy. Ta macierz jest reprezentacją współrzędnych mapy liniowej w n -wymiarowej przestrzeni wielomianów stopnia n -1 lub mniejszego. Główną ideą tej mapy jest interpretacja rozkładu na czynniki

${\ Displaystyle P (X) = (X- \ alfa _ {1}) \ cdot P_ {1} (X)}$

z pierwiastkiem ${\ Displaystyle \ alpha _ {1} \ in \ mathbb {C}}$ i ${\ Displaystyle P_ {1} ( X)=P(X)/(X-\alpha _{1})}$ pozostały współczynnik stopnia n − 1 jako równanie wektora własnego dla mnożenia przez zmienną X , po którym następuje obliczenie reszty z dzielnikiem P ( X ),

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

To odwzorowuje wielomiany stopnia co najwyżej n - 1 na wielomiany stopnia co najwyżej n - 1. Wartości własne tej mapy są pierwiastkami P ( X ), ponieważ równanie wektora własnego brzmi

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

co oznacza, że $\ Displaystyle (X- \ alfa) \ cdot H) = C \ cdot P (X)}$ , czyli ${ \displaystyle (X-\alpha )}$ jest współczynnikiem liniowym P ( X ). W bazie jednomianowej mapa liniowa jest reprezentowana przez towarzyszącą macierz wielomianu ${\ displaystyle M_ {X}}$ P. , jak

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

wynikowa macierz współczynników to

${\ Displaystyle A = {\ rozpocząć {pmatrix} 0&0&\ kropki & 0&-P_ {0}\\1&0&\ kropki & 0&-P_ {1}\\0&1&\ kropki & 0&-P_ {2}\\\ vdots &\ vdots & \ddots &\vdots &\vdots \\0&0&\dots &1&-P_{n-1}\end{pmatrix}}\,.}$

Do tej macierzy stosuje się odwrotną iterację potęgową w trzech wariantach bez przesunięcia, stałego przesunięcia i uogólnionego przesunięcia Rayleigha w trzech etapach algorytmu. Bardziej wydajne jest wykonywanie operacji algebry liniowej w arytmetyce wielomianowej, a nie przez operacje macierzowe, jednak właściwości iteracji potęgi odwrotnej pozostają takie same.

Rzeczywiste współczynniki

Opisany wcześniej algorytm Jenkinsa-Trauba działa dla wielomianów o zespolonych współczynnikach. Ci sami autorzy stworzyli również trójstopniowy algorytm dla wielomianów o rzeczywistych współczynnikach. Zobacz: Jenkins i Traub Trzyetapowy algorytm dla wielomianów rzeczywistych z wykorzystaniem iteracji kwadratowej . Algorytm znajduje czynnik liniowy lub kwadratowy działający całkowicie w rzeczywistej arytmetyce. Jeśli złożone i rzeczywiste algorytmy zostaną zastosowane do tego samego rzeczywistego wielomianu, rzeczywisty algorytm jest około cztery razy szybszy. Prawdziwy algorytm zawsze jest zbieżny, a stopień zbieżności jest większy niż drugiego rzędu.

Połączenie z przesuniętym algorytmem QR

Istnieje zaskakujący związek z przesuniętym algorytmem QR do obliczania wartości własnych macierzy. Patrz Dekker i Traub Przesunięty algorytm QR dla macierzy hermitowskich . Ponownie przesunięcia można postrzegać jako iterację Newtona-Raphsona na sekwencji funkcji wymiernych zbiegających się do wielomianu pierwszego stopnia.

Oprogramowanie i testowanie

Oprogramowanie dla algorytmu Jenkinsa-Trauba zostało opublikowane jako Jenkins and Traub Algorithm 419: Zeros of a Complex Polynomial . Oprogramowanie dla algorytmu rzeczywistego zostało opublikowane jako Jenkins Algorithm 493: Zeros of a Real Polynomial .

Metody zostały gruntownie przetestowane przez wiele osób. Zgodnie z przewidywaniami cieszą się zbieżnością szybszą niż kwadratowa dla wszystkich rozkładów zer.

Istnieją jednak wielomiany, które mogą powodować utratę precyzji, co ilustruje poniższy przykład. Wielomian ma wszystkie swoje zera leżące na dwóch półokręgach o różnych promieniach. Wilkinson zaleca, aby dla stabilnej deflacji pożądane było, aby najpierw obliczyć mniejsze zera. Przesunięcia drugiego etapu są dobierane tak, aby zera na mniejszym półokręgu znajdowały się jako pierwsze. Wiadomo, że po deflacji wielomian z zerami na półkolu jest źle uwarunkowany, jeśli stopień jest duży; patrz Wilkinson, s. 64. Pierwotny wielomian miał stopień 60 i miał poważną niestabilność deflacyjną.

Linki zewnętrzne

• Bezpłatna aplikacja Windows do pobrania wykorzystująca metodę Jenkinsa – Trauba dla wielomianów o rzeczywistych i zespolonych współczynnikach
• RPoly++ Zoptymalizowana pod kątem SSE implementacja algorytmu RPOLY w języku C++.