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

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

В данной статье описан сложный вариант. Учитывая полином P , $$P(z)=\sum _{i=0}^{n}a_{i}z^{n-i},\quad a_{0}=1,\quad a_{n}\neq 0$$ с комплексными коэффициентами он вычисляет приближения к n нулям $\alpha _{1},\alpha _{2},\dots ,\alpha _{n}$ P ) ( z , по одному в примерно возрастающем порядке. После вычисления каждого корня его линейный коэффициент удаляется из многочлена. Использование этого дефлятирования гарантирует, что каждый корень вычисляется только один раз и что все корни найдены.

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

Обзор

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

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

Начиная с текущего многочлена P ( X ) степени n , цель состоит в том, чтобы вычислить наименьший корень $\alpha$ P (x) . Затем полином можно разделить на линейный множитель и оставшийся полиномиальный множитель. $$P(X)=(X-\alpha ){\bar {H}}(X)$$ Другие методы поиска корня направлены в первую очередь на улучшение корня и, следовательно, первого фактора. Основная идея метода Дженкинса-Трауба заключается в постепенном улучшении второго фактора.

Для этого последовательность так называемых H- строится полиномов. Все эти многочлены имеют степень n - 1 и должны сходиться к множителю ${\bar {H}}(X)$ P ) , ( X содержащий (линейные множители) всех остальных корней. Последовательность полиномов H встречается в двух вариантах: ненормализованном варианте, который позволяет легко теоретическое понимание, и нормализованном варианте. ${\bar {H}}$ полиномы, которые удерживают коэффициенты в численно разумном диапазоне. Построение H- полиномов $\left(H^{(\lambda )}(z)\right)_{\lambda =0,1,2,\dots }$ руководствуется последовательностью комплексных чисел $(s_{\lambda })_{\lambda =0,1,2,\dots }$ называются сменами. Сами эти сдвиги зависят, по крайней мере на третьем этапе, от предыдущих H. полиномов Полиномы H определяются как решение неявной рекурсии $$H^{(0)}(z)=P^{\prime }(z)$$ и $$(X-s_{\lambda })\cdot H^{(\lambda +1)}(X)\equiv H^{(\lambda )}(X){\pmod {P(X)}}\ .$$ Прямое решение этого неявного уравнения имеет вид $$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)\,,$$ где полиномиальное деление является точным.

Алгоритмически для оценки при полиномов $s_{\lambda }$ и одновременно получить частное. С полученными факторами p ( X ) и h ( X ) в качестве промежуточных результатов следующий полином H получается как $$\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).$$ Поскольку коэффициент наивысшей степени получается из P(X) , старший коэффициент $H^{(\lambda +1)}(X)$ является $-{\tfrac {H^{(\lambda )}(s_{\lambda })}{P(s_{\lambda })}}$. Если это разделить, нормализованный полином H будет равен $${\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}}$$

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

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

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

Сдвиг на этом этапе определяется как некоторая точка, близкая к наименьшему корню многочлена. Он квазислучайно расположен на окружности с внутренним корневым радиусом, который, в свою очередь, оценивается как положительное решение уравнения $$R^{n}+|a_{n-1}|\,R^{n-1}+\dots +|a_{1}|\,R=|a_{0}|\,.$$ Поскольку левая часть является выпуклой функцией и монотонно возрастает от нуля до бесконечности, это уравнение легко решить, например, методом Ньютона .

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

$$t_{\lambda }=s-{\frac {P(s)}{{\bar {H}}^{(\lambda )}(s)}}$$ прослеживается. Второй этап завершается как успешный, если выполняются условия $$|t_{\lambda +1}-t_{\lambda }|<{\tfrac {1}{2}}\,|t_{\lambda }|$$ и $$|t_{\lambda }-t_{\lambda -1}|<{\tfrac {1}{2}}\,|t_{\lambda -1}|$$ одновременно встречаются. Это ограничивает относительный размер шага итерации, гарантируя, что последовательность аппроксимации останется в диапазоне меньших корней. Если после некоторого количества итераций успех не увенчался успехом, пробуется другая случайная точка круга. Обычно для полиномов средней степени используется девять итераций со стратегией удвоения в случае множественных сбоев.

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

The $H^{(\lambda +1)}(X)$ полиномы теперь генерируются с использованием переменных сдвигов $s_{\lambda },\quad \lambda =L,L+1,\dots$ которые генерируются $$s_{L}=t_{L}=s-{\frac {P(s)}{{\bar {H}}^{(L)}(s)}}$$ являющаяся последней корневой оценкой второго этапа и $$s_{\lambda +1}=s_{\lambda }-{\frac {P(s_{\lambda })}{{\bar {H}}^{(\lambda +1)}(s_{\lambda })}},\quad \lambda =L,L+1,\dots ,$$ где ${\bar {H}}^{(\lambda +1)}(z)$ – нормированный полином H , т.е. $H^{(\lambda )}(z)$ разделенное на его ведущий коэффициент.

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

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

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

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

Что дает алгоритму его мощь?

Сравните с итерацией Ньютона – Рафсона. $$z_{i+1}=z_{i}-{\frac {P(z_{i})}{P^{\prime }(z_{i})}}.$$

Итерация использует заданные P и $\scriptstyle P^{\prime }$. В отличие от третьей стадии Дженкинса-Трауба $$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 })}}$$

это в точности итерация Ньютона-Рафсона, выполняемая над некоторыми рациональными функциями . Точнее, метод Ньютона–Рафсона выполняется над последовательностью рациональных функций. $$W^{\lambda }(z)={\frac {P(z)}{H^{\lambda }(z)}}.$$

Для $\lambda$ достаточно большой, $${\frac {P(z)}{{\bar {H}}^{\lambda }(z)}}=W^{\lambda }(z)\,LC(H^{\lambda })$$ максимально близко к полиному первой степени $$z-\alpha _{1},\,$$ где $\alpha _{1}$ является одним из нулей $P$. Несмотря на то, что этап 3 представляет собой именно итерацию Ньютона-Рафсона, дифференцирование не выполняется.

Анализ H полиномов

Позволять $\alpha _{1},\dots ,\alpha _{n}$ быть корнями P ( X ). так называемые факторы Лагранжа P(X) Кофакторами этих корней являются . $$P_{m}(X)={\frac {P(X)-P(\alpha _{m})}{X-\alpha _{m}}}.$$ Если все корни различны, то факторы Лагранжа образуют базис пространства многочленов степени не выше n - 1. Путем анализа рекурсивной процедуры обнаруживается, что многочлены H имеют координатное представление $$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 имеют вид $${\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]}}\ .$$

Порядки конвергенции

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

При условии, что $$|\alpha _{1}|<|\alpha _{2}|=\min {}_{m=2,3,\dots ,n}|\alpha _{m}|$$ получают асимптотические оценки для

• этап 1: $$H^{(\lambda )}(X)=P_{1}(X)+O\left(\left|{\frac {\alpha _{1}}{\alpha _{2}}}\right|^{\lambda }\right).$$
• для этапа 2, если s достаточно близко к $\alpha _{1}$: $$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)$$ и $$s-{\frac {P(s)}{{\bar {H}}^{(\lambda )}(s)}}=\alpha _{1}+O\left(\ldots \cdot |\alpha _{1}-s|\right).$$
• и для этапа 3: $$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)$$ и $$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)$$ что приводит к более высокому, чем квадратичный, порядку сходимости $\phi ^{2}=1+\phi \approx 2.61$, где $\phi ={\tfrac {1}{2}}(1+{\sqrt {5}})$ это золотое сечение .

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

Все этапы комплексного алгоритма Дженкинса–Трауба можно представить как задачу линейной алгебры об определении собственных значений специальной матрицы. Эта матрица является координатным представлением линейного отображения в n -мерном пространстве многочленов степени n - 1 или меньше. Основная идея этой карты заключается в интерпретации факторизации $$P(X)=(X-\alpha _{1})\cdot P_{1}(X)$$ с корнем $\alpha _{1}\in \mathbb {C}$ и $P_{1}(X)=P(X)/(X-\alpha _{1})$ оставшийся множитель степени n - 1 как уравнение собственного вектора для умножения на переменную X с последующим вычислением остатка с делителем P ( X ), $$M_{X}(H)=(X\cdot H(X)){\bmod {P}}(X)\,.$$ Это отображает полиномы степени не выше n - 1 в полиномы степени не выше n - 1. Собственные значения этого отображения являются корнями P ( X ), поскольку уравнение собственных векторов имеет вид $$0=(M_{X}-\alpha \cdot id)(H)=((X-\alpha )\cdot H){\bmod {P}}\,,$$ что подразумевает, что $(X-\alpha )\cdot H)=C\cdot P(X)$, то есть, $(X-\alpha )$ является линейным фактором P ( X ). В мономиальном базисе линейное отображение $M_{X}$ представляется сопутствующей матрицей полинома P , как $$M_{X}(H)=\sum _{m=0}^{n-1}H_{m}X^{m+1}-H_{n-1}\left(X^{n}+\sum _{m=0}^{n-1}a_{m}X^{m}\right)=\sum _{m=1}^{n-1}(H_{m-1}-a_{m}H_{n-1})X^{m}-a_{0}H_{n-1}\,,$$ результирующая матрица преобразования $$A={\begin{pmatrix}0&0&\dots &0&-a_{0}\\1&0&\dots &0&-a_{1}\\0&1&\dots &0&-a_{2}\\\vdots &\vdots &\ddots &\vdots &\vdots \\0&0&\dots &1&-a_{n-1}\end{pmatrix}}\,.$$ К этой матрице применяется обратная степенная итерация в трех вариантах: без сдвига, с постоянным сдвигом и обобщенным сдвигом Рэлея на трех этапах алгоритма. Операции линейной алгебры эффективнее выполнять в полиномиальной арифметике, а не матричными операциями, однако свойства обратной степенной итерации остаются прежними.

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

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

Связь со смещенным алгоритмом QR

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

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

Программное обеспечение для алгоритма Дженкинса-Трауба было опубликовано как « Алгоритм Дженкинса и Трауба 419: нули комплексного полинома» . Программное обеспечение для реального алгоритма было опубликовано как «Алгоритм Дженкинса 493: нули вещественного полинома» .

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

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

Внешние ссылки

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