Алгоритм Дженкинса – Трауба - Jenkins–Traub algorithm

В Алгоритм Дженкинса – Трауба для полиномиальных нулей - это быстрый глобально сходящийся итерационный метод, опубликованный в 1970 г. Майкл А. Дженкинс и Джозеф Ф. Трауб. Они предложили два варианта: один для общих многочленов с комплексными коэффициентами, широко известный как алгоритм «CPOLY», и более сложный вариант для частного случая многочленов с действительными коэффициентами, широко известный как алгоритм «RPOLY». Последнее является «практически стандартом для полиномиальных кортежей черного ящика».^[1]

В этой статье описан сложный вариант. Учитывая многочлен п,

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

с комплексными коэффициентами вычисляет приближения к п нули ${\displaystyle \alpha _{1},\alpha _{2},\dots ,\alpha _{n}}$ из п(z) по одному в примерно возрастающем порядке. После вычисления каждого корня его линейный множитель удаляется из многочлена. Используя это дефляция гарантирует, что каждый корень вычисляется только один раз и что все корни найдены.

Реальный вариант следует тому же шаблону, но одновременно вычисляет два корня, либо два действительных корня, либо пару сопряженных комплексных корней. Избегая сложной арифметики, реальный вариант может быть быстрее (в 4 раза), чем сложный вариант. Алгоритм Дженкинса – Трауба стимулировал значительные исследования теории и программного обеспечения для методов этого типа.

Обзор

Алгоритм Дженкинса – Трауба вычисляет все корни многочлен с комплексными коэффициентами. Алгоритм начинается с проверки полинома на наличие очень больших или очень маленьких корней. При необходимости коэффициенты масштабируются путем изменения масштаба переменной. В алгоритме правильные корни находятся один за другим и, как правило, в увеличивающемся размере. После того, как каждый корень найден, полином сдувается путем деления соответствующего линейного множителя. Действительно, факторизация полинома в линейный множитель и оставшийся спущенный полином уже является результатом процедуры поиска корня. Процедура поиска корней состоит из трех этапов, которые соответствуют разным вариантам обратная степенная итерация. См. Дженкинс и Трауб.^[2]Описание также можно найти в Ральстоне и Рабинович^[3] п. 383. Алгоритм аналогичен по духу двухэтапному алгоритму, исследованному Траубом.^[4]

Процедура поиска корней

Начиная с текущего полинома п(Икс) степени п, наименьший корень из Р (х) вычисляется. С этой целью была проведена серия так называемых ЧАС многочленов построено. Все эти многочлены имеют степень п - 1 и должны сходиться к коэффициенту п(Икс), содержащую все остальные корни. Последовательность ЧАС полиномы встречаются в двух вариантах: ненормализованный вариант, который позволяет легко получить теоретическое понимание, и нормализованный вариант ${\displaystyle {\bar {H}}}$ полиномы, которые удерживают коэффициенты в числовом значении.

Строительство ЧАС многочлены ${\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)}$ и ${\displaystyle (X-s_{\lambda })\cdot H^{(\lambda +1)}(X)\equiv H^{(\lambda )}(X){\pmod {P(X)}}\ .}$

Прямым решением этого неявного уравнения является

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

где полиномиальное деление точное.

Алгоритмически можно использовать, например, Схема Хорнера или же Правило Руффини вычислить полиномы в ${\displaystyle s_{\lambda }}$ и одновременно получить частные. С полученными частными п(Икс) и час(Икс) в качестве промежуточных результатов следующие ЧАС полином получается как

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

Поскольку коэффициент высшей степени получается из Р (Х), старший коэффициент ${\displaystyle H^{(\lambda +1)}(X)}$ является ${\displaystyle -{\tfrac {H^{(\lambda )}(s_{\lambda })}{P(s_{\lambda })}}}$. Если это разделить нормализованное ЧАС многочлен

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

Этап первый: бессменный процесс

За ${\displaystyle \lambda =0,1,\dots ,M-1}$ набор ${\displaystyle s_{\lambda }=0}$. Обычно M = 5 выбирается для многочленов умеренных степеней до п = 50. Этот этап не является необходимым только из теоретических соображений, но полезен на практике. Это подчеркивается в ЧАС полиномы - сомножитель (линейного множителя) наименьшего корня.

Этап второй: процесс фиксированной смены

Сдвиг для этого этапа определяется как некоторая точка, близкая к наименьшему корню полинома. Он квазислучайно расположен на окружности с внутренним радиусом корня, который, в свою очередь, оценивается как положительное решение уравнения

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

Поскольку левая часть является выпуклой функцией и монотонно возрастает от нуля до бесконечности, это уравнение легко решить, например, с помощью Метод Ньютона.

Теперь выберите ${\displaystyle s=R\cdot \exp(i\,\phi _{\text{random}})}$ на окружности этого радиуса. Последовательность многочленов ${\displaystyle H^{(\lambda +1)}(z)}$, ${\displaystyle \lambda =M,M+1,\dots ,L-1}$, генерируется с фиксированным значением сдвига ${\displaystyle s_{\lambda }=s}$. Во время этой итерации текущее приближение для корня

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

прослеживается. Второй этап завершается успешно, если условия

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

одновременно встречаются. Если после некоторого количества итераций успеха не было, пробуется другая случайная точка на круге. Обычно используется 9 итераций для полиномов средней степени со стратегией удвоения в случае множественных отказов.

Этап третий: процесс переменной смены

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

является последней оценкой корня второго этапа и

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

куда ${\displaystyle {\bar {H}}^{(\lambda +1)}(z)}$ нормализованный ЧАС полином, то есть ${\displaystyle H^{(\lambda )}(z)}$ делится на его старший коэффициент.

Если размер шага на третьем этапе недостаточно быстро падает до нуля, второй этап запускается заново с другой случайной точкой. Если это не удается после небольшого числа перезапусков, количество шагов на втором этапе удваивается.

Конвергенция

Можно показать, что при условии L выбирается достаточно большим, s_λ всегда сходится к корню п.

Алгоритм сходится для любого распределения корней, но может не найти все корни многочлена. Кроме того, сходимость немного быстрее, чем квадратичная сходимость итерации Ньютона – Рафсона, однако он использует как минимум в два раза больше операций на шаг.

В чем сила алгоритма?

Сравните с Итерация Ньютона – Рафсона

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

Итерация использует данный п и ${\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 })}}}$

в точности итерация Ньютона – Рафсона, выполняемая на определенных рациональные функции. Точнее, Ньютон – Рафсон выполняется на последовательности рациональных функций

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

За ${\displaystyle \lambda }$ достаточно большой,

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

как можно ближе к полиному первой степени

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

куда ${\displaystyle \alpha _{1}}$ один из нулей ${\displaystyle P}$. Несмотря на то, что стадия 3 представляет собой итерацию Ньютона – Рафсона, дифференцирование не выполняется.

Анализ ЧАС многочлены

Позволять ${\displaystyle \alpha _{1},\dots ,\alpha _{n}}$ быть корнями п(Икс). Так называемые факторы Лагранжа Р (Х) являются кофакторами этих корней,

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

Если все корни разные, то множители Лагранжа образуют базис пространства многочленов степени не выше п - 1. Анализируя процедуру рекурсии, обнаруживаем, что ЧАС многочлены имеют координатное представление

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

Каждый фактор Лагранжа имеет старший коэффициент 1, так что старший коэффициент полиномов H является суммой коэффициентов. Таким образом, нормированные полиномы H имеют вид

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

Приказы конвергенции

Если условие ${\displaystyle |\alpha _{1}-s_{\kappa }|<\min {}_{m=2,3,\dots ,n}|\alpha _{m}-s_{\kappa }|}$ выполняется почти для всех итераций, нормализованные полиномы H будут сходиться по крайней мере геометрически к ${\displaystyle P_{1}(X)}$.

При условии, что

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

получаем асимптотические оценки для

• этап 1:

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

• для этапа 2, если s достаточно близко к ${\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)}$

  и

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

• и для этапа 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)}$

  и

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

приводя к более высокому порядку сходимости, чем квадратичный ${\displaystyle \phi ^{2}=1+\phi \approx 2.61}$, куда ${\displaystyle \phi ={\tfrac {1}{2}}(1+{\sqrt {5}})}$ это Золотое сечение.

Интерпретация как итерация обратной мощности

Все этапы комплексного алгоритма Дженкинса – Трауба могут быть представлены как задача линейной алгебры определения собственных значений специальной матрицы. Эта матрица является координатным представлением линейной карты в п-мерное пространство многочленов степени п - 1 или меньше. Основная идея этой карты - интерпретировать факторизацию

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

с корнем ${\displaystyle \alpha _{1}\in \mathbb {C} }$ и ${\displaystyle P_{1}(X)=P(X)/(X-\alpha _{1})}$ оставшийся фактор степени п - 1 как уравнение на собственный вектор для умножения с переменной Иксс последующим вычислением остатка с делителем п(Икс),

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

Это отображает многочлены степени не выше п - от 1 до многочленов степени не выше п - 1. Собственные значения этого отображения являются корнями п(Икс), поскольку уравнение собственного вектора имеет вид

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

откуда следует, что ${\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}\,,}$

результирующая матрица коэффициентов

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

К этой матрице обратная степенная итерация применяется в трех вариантах: без сдвига, постоянного сдвига и обобщенного рэлеевского сдвига на трех этапах алгоритма. Более эффективно выполнять операции линейной алгебры в полиномиальной арифметике, а не с помощью матричных операций, однако свойства обратной степенной итерации остаются прежними.

Реальные коэффициенты

Описанный ранее алгоритм Дженкинса – Трауба работает для многочленов с комплексными коэффициентами. Те же авторы создали трехэтапный алгоритм для многочленов с действительными коэффициентами. Увидеть Дженкинса и Трауба Трехэтапный алгоритм для вещественных многочленов с использованием квадратичной итерации.^[5] Алгоритм находит линейный или квадратичный фактор, полностью работающий в реальной арифметике. Если комплексный и реальный алгоритмы применяются к одному и тому же действительному многочлену, реальный алгоритм будет примерно в четыре раза быстрее. Реальный алгоритм всегда сходится, и скорость сходимости больше второго порядка.

Связь со сдвинутым QR-алгоритмом

Есть удивительная связь со сдвинутым QR-алгоритмом для вычисления собственных значений матрицы. См. Деккер и Трауб Сдвинутый QR-алгоритм для эрмитовых матриц.^[6] Опять же, сдвиги можно рассматривать как итерацию Ньютона-Рафсона для последовательности рациональных функций, сходящихся к полиному первой степени.

Программное обеспечение и тестирование

Программное обеспечение для алгоритма Дженкинса – Трауба было опубликовано как Jenkins and Traub Алгоритм 419: нули комплексного многочлена.^[7] Программное обеспечение для настоящего алгоритма было опубликовано как Jenkins. Алгоритм 493: нули действительного многочлена.^[8]

Методы были протестированы многими людьми. Как и предполагалось, они пользуются более быстрой, чем квадратичная сходимость для всех распределений нулей.

Однако есть многочлены, которые могут вызвать потерю точности, как показано в следующем примере. Все нули полинома лежат на двух полукругах разного радиуса. Уилкинсон рекомендует, чтобы для стабильного дефляции сначала вычислялись меньшие нули. Сдвиги второй ступени выбираются так, чтобы сначала находились нули на меньшем полукруге. Известно, что после дефляции многочлен с нулями на полукруге плохо обусловлен, если степень велика; см. Уилкинсон,^[9] п. 64. Исходный полином был 60-й степени и страдал от сильной дефляционной нестабильности.

Рекомендации

1. Press, W.H., Teukolsky, S.A., Vetterling, W.T. и Flannery, B.P. (2007), Численные рецепты: Искусство научных вычислений, 3-е изд., Cambridge University Press, стр. 470.
2. Дженкинс, М. А. и Трауб, Дж. Ф. (1970), Трехэтапная итерация сдвига переменных для полиномиальных нулей и ее связь с обобщенной итерацией Рэлея, Нумер. Математика. 14, 252–263.
3. Ральстон А. и Рабиновиц П. (1978), Первый курс численного анализа, 2-е изд., Макгроу-Хилл, Нью-Йорк.
4. Трауб, Дж. Ф. (1966), Класс глобально сходящихся итерационных функций для решения полиномиальных уравнений, Математика. Comp., 20 (93), 113–138.
5. Дженкинс, М. А. и Трауб, Дж. Ф. (1970), Трехэтапный алгоритм для вещественных многочленов с использованием квадратичной итерации, SIAM J. Numer. Анализ., 7 (4), 545–566.
6. Деккер, Т. Дж. И Трауб, Дж. Ф. (1971), Сдвинутый QR-алгоритм для эрмитовых матриц, Лин. Алгебра и приложения, 4 (2), 137–154.
7. Дженкинс, М. А. и Трауб, Дж. Ф. (1972), Алгоритм 419: нули комплексного многочлена, Comm. ACM, 15, 97–99.
8. Дженкинс, М.А. (1975), Алгоритм 493: нули действительного многочлена, АСМ ТОМС, 1, 178–189.
9. Уилкинсон, Дж. Х. (1963), Ошибки округления в алгебраических процессах, Prentice Hall, Englewood Cliffs, N.J.

внешняя ссылка

• Бесплатное загружаемое приложение для Windows, использующее метод Дженкинса – Трауба для полиномов с действительными и комплексными коэффициентами.
• RPoly ++ Оптимизированная для SSE реализация C ++ алгоритма RPOLY.