nAGライブラリソルバーで多項式の根を高速に求める

Published 11/05/2021 By Edvin Hopkins

多項式のゼロを見つけることは、数学における長年の問題であり、金融、物理、工学、制御理論、信号処理...などに応用されています。このような古くて古典的な問題は、今までに完全に解決されているに違いないと考えたくなりますが、実はそうではありません。

実際に、単純な二次方程式より難しい問題を調べずとも、これを見る事ができます。以下は、標準的な二次式 $$x=\frac{-b\pm \sqrt{b^2-4ac}}{2a}.$$ を実装したPythonコードの一例です。

from cmath import sqrt
def quad_solve(a,b,c):
    return (-b+sqrt(b**2-4*a*c))/(2*a), (-b-sqrt(b**2-4*a*c))/(2*a)

これを使って方程式 $x^2+{10}^{9}x+1=0$ を解いてみましょう。

>>> quad_solve(1,1e9,1)
(0j, (-1000000000+0j))

上記では、根の1つがゼロですが、二次関数を簡単に調べると、これは明らかにナンセンスです。この結果となってしまう理由は、x 係数が他の係数に比べて非常に大きい場合、式の中で桁落ちが発生してしまうからです。

nAG Library for Python のルーチン quadratic_real と quadratic_complex（それぞれ実数および複素数の二次方程式の根を求める）は、このような落とし穴を避けるように細心の注意を払って作成されています。以下のコード例では同じ二次方程式を解くために quadratic_real を呼び出しています。

>>> from naginterfaces.library import zeros
>>> zeros.quadratic_real(1,1e9,1)
QuadraticRealReturnData(z_small=(-1e-09+0j), z_large=(-1000000000+0j))

上記では、ゼロに近い根がより賢明に評価されていることがわかります。

図1: 係数が+/-1で次数が18以下のすべての多項式の根

もちろん、高次の多項式の場合はさらに厄介で、5次以上では根の閉形式が存在しないため、反復法を用いなければなりません。1つの方法は、根を見つけ（ニュートン法などの標準的な反復法を使用）、その根を割り算して（これをデフレといいます）、プロセスを繰り返すことです。しかし、Wilkinson's polynomial（nAGの初期の支援者であるJames Hardy Wilkinsonにちなんで命名）は、このアプローチが失敗する理由を示しています。多項式のたった一つの係数の小さな乱れが、根の値や性質を劇的に変えてしまうのです。つまり、この問題は悪条件です。より洗練されたアプローチが必要です。

この問題に最適なアルゴリズムの1つは、ラゲール法です。これは、単純な根に対して3次収束し、一般に非常に安定しています。 nAG Library for Pthonでは、それぞれ poly_complex（複素係数の場合）およびpoly_real（実係数の場合）として実装されています。 ただし、ラゲール法でさえ、十分に大きな問題や条件の悪い問題には苦戦します。

最近、Penn State Behrend大学の数学助教授であるThomas R. Cameron氏は、ラゲール法の修正版を開発しました。この修正版では、暗黙のデフレ戦略を用いて手法の性能と精度を向上させています。 nAGではThomasさんとのコラボレーションに成功し、このアルゴリズムはpoly_complex_fpml（複素係数を持つ多項式用）およびpoly_real_fpml（実係数を持つ多項式用）として、nAGライブラリのMark27.1で利用できるようになりました。

図2のグラフは、一般的な多項式において、新しいルーチンが古いルーチンをいかに凌駕しているかを示しています（実際、より大きな次数では、古いルーチンが全く収束しないこともありました）。新しいルーチンは古いルーチンの2倍の速さであり、相対誤差は小さいままですが、古いルーチンでは相対誤差が多項式の次数に応じて直線的に増加しています。

図2: [-1000,1000]でランダムに生成された係数を持つ多項式に対するnAGルーチンの比較

当初我々がこの新しいアルゴリズムについて調べていた時、いくつかの特に悪条件の問題（Wilkinson多項式など）では、新しいルーチンが古いルーチンよりも精度の低い結果をもたらすことがわかりました。 しかしその後の開発で、計算速度を多少犠牲にしますが、より良い結果を得るための追加のPolish引数を導入しました。これによりユーザーはより良い結果を得るために、Polish方法の検討（模索）ができるようになりました。その結果、誤差も桁違いに減少し、新しいルーチンが古いルーチンを上回る事になりました。 したがって、悪条件の問題については、Polishを使用することをお勧めします。 （ルーチンは条件数も返すことに注意してください。）

図 1 は、poly_real_fpml で polish=1 を指定して計算したもので、、係数が +/-1 で次数が 18 以下のすべての多項式の複素半面上の根を示しています。このルーチンでは、これらの根を計算することに問題はありませんでしたが、同等なを オリジナルルーチンであるpoly_realで行おうとしたところ、多くの問題が発生しました。

このようなプロットを自分で生成してみたい場合は、以下のPythonコードをお試しください。

import numpy as np
import matplotlib.pyplot as plt
from naginterfaces.library import zeros
nmax = 18
roots = np.empty(0)
for n in range(1,nmax+1):
    a = np.ones(n+2)
    a[n] = -1
    j = 0
    while j=0.0,z[0])))
        while a[j]==-1:
            j=j+1
        if j<=n:
            a[j]=-1
            a[0:j] = 1
            j=0
plt.scatter(roots.real,roots.imag,s=0.05)

新しい多項式求根ルーチンは、nAGライブラリのMark 27.1で、poly_complex_fpml（複素係数を持つ多項式の場合）およびpoly_real_fpml（実係数を持つ多項式の場合）として利用できます。 古いルーチン（poly_complexおよびpoly_real）は非推奨になりました。 本記事ではPythonインターフェースの例を示しましたが、同様な事は、Fortranインターフェース（c02aafおよびc02abf）およびCインターフェース（c02aacおよびc02abc）を使用しても行う事ができます。

最後に本記事を書くにあたり貴重な情報を提供してくれたLawrence MulhollandとJoe Davisonに感謝します。