JP / EN

Wilkinson's Polynomial

数値解析の「怪物」とホモトピー法による全解探索の奇跡

$$W(x) = \prod_{i=1}^{20} (x - i) = (x-1)(x-2)\cdots(x-20) = 0$$

Wilkinson's Polynomial
Roots: $x = 1, 2, 3, \dots, 20$

1. 序論：最も無害な顔をした怪物

ウィルキンソンの多項式は、1959年にジェームズ・H・ウィルキンソンによって提示された、数値解析における「悪条件問題（Ill-conditioned problem）」の代表格です。 解は $x = 1$ から $20$ までの整数であり、一見すると中学生でも解けるほど単純に思えます。

しかし、この式を展開して係数を求め、それを計算機（浮動小数点演算）に解かせようとした途端に牙を剥きます。 展開後の $x^{19}$ の係数である $-210$ に対して、わずか $2^{-23}$ という単精度浮動小数点の1ビットにも満たないノイズを加えただけで、実数だった解の半分以上が複素数に化けてしまうのです。

この「丸め誤差に対する圧倒的な脆弱性」を持つ多項式に対し、あえてホモトピー法を用いて挑むとどうなるのか。3つの異なるアプローチで検証を行いました。

2. ホモトピー法による挑戦

今回は、M.A.L.のコア技術である「SPICE指向型解析法（SPICE-oriented approach）」を用いてホモトピー法を実現し、解の探索を行いました。

① ニュートン法（標準ソルバー）

結果：解ける（ただし局所的）

初期値（例えば $x_0 = 0.1$）を与えて通常のニュートン法を走らせた場合、巨大な勾配に背中を押される形で、最も近い解である $x = 1$ にストンと落ちて収束しました。ただし、見つかるのはこの1つの解のみです。

② 不動点ホモトピー法（Fixed-Point Homotopy）

結果：敗北（Time step too small）

全解探索を目指して以下のホモトピー方程式を組みました。
$$H(x,t) = t \cdot W(x) + (1-t)(x - x_0) = 0$$

結果は、過渡解析ソルバーが破綻してエラー停止。原因は「スケールの不一致」です。
$W(x)$ の局所極大値は約 $10^{13}$ に達しますが、ペナルティ項 $(x-x_0)$ は高々数十のオーダーです。この絶望的なダイナミックレンジの差によりヤコビアン行列が破綻し、シミュレータはタイムステップを極限まで縮小した挙句に計算を放棄しました。

③ ニュートン・ホモトピー法（Newton Homotopy）

結果：勝利（20個の全解を完走）

次に、方程式を以下のように組み直しました。
$$H(x,t) = W(x) - (1-t)W(x_0) = 0$$

すると、エラーを吐くことなくシミュレーションが最後まで完走し、信じられないほど美しい軌跡を描いてすべての解を見つけ出しました。

3. 奇跡の「自己スケーリング（Self-Scaling）」

なぜニュートン・ホモトピー法だけが、倍精度浮動小数点（double）の限界を突破できたのでしょうか？
それは、方程式の構造が引き起こした「自己スケーリング」という魔法のおかげです。

方程式を $t$ について解くと、パスの軌跡は次のように表せます。
$$1 - t = \frac{W(x)}{W(x_0)}$$ ここで、初期値 $x_0 = 0.1$ とした場合の分母 $W(x_0)$ は、約 $2.4 \times 10^{18}$ というとてつもない巨大な値になります。

方程式 $W(x)$ が解と解の間で約 $10^{13}$ まで暴れ狂ったとしても、それを圧倒的なスケールを持つ定数 $10^{18}$ で割る形になるため、$t$ の変動幅は高々 $10^{-5}$（0.00001）のオーダーに完璧に押し込められます。

「巨大な値に小さな値を足す」という歪みがなくなり、方程式自らがダンパーとして機能することで、弧長法（Arc Length）ソルバーがパスを見失わずに進むことができたのです。

4. 解析結果：驚異の全解探索

以下が、SPICE指向型解析法を使って実現したニュートン・ホモトピー法のシミュレーション結果（$t - x$ 平面の軌跡）です。

Graph: Wilkinson Polynomial Newton Homotopy Result 1

図1: 全体図。$t=1.0$ のラインに張り付きながら、$x$ が右方向へスライドしているように見える。

Graph: Wilkinson Polynomial Newton Homotopy Result 2 (Zoomed)

図2: $t=1.0$ 付近の拡大図。微小な振動を繰り返しながら、$x=1$ から $20$ までのすべての解を串刺しにしている。

Dr.WataWata Insight：

拡大図（図2）を見ると、縦軸 $t$ の変動が 0.999990 から 1.000010 の間に見事に収まっています。理論通りの $10^{-5}$ のスケーリングが実証されています。

「悪条件問題」として恐れられるウィルキンソンの多項式も、アルゴリズムの選択（ニュートン・ホモトピーへの切り替え）一つで、ここまで美しく調教することができます。
数学的な性質とソルバーの数値的な挙動（ヤコビアンのスケール）を同時に見極めること。これこそが、数理工学における真の醍醐味です。

5. 解析環境

  • 解析エンジン： Mac SPICE3 (SPICE-oriented approach)
  • PC： MacBook
  • OS： macOS Monterey 12.7.6
  • CPU： 1.2GHz デュアルコア Intel Core m5
  • メモリ： 8GB