多変数 Newton 法の反復はどのように書けますか？

F: R^n → R^n の零点を求めるとき、Jacobian J(x) = ∂F/∂x を用いて x_{k+1} = x_k - J(x_k)^{-1} F(x_k) と更新します。実装上は逆行列を直接計算せず、線形系 J(x_k) Δ = -F(x_k) を LU 分解 (または QR) で解いて Δ を得て x_{k+1} = x_k + Δ とします。

Broyden 法は何が便利ですか？

Broyden 法は Jacobian の代わりに近似行列 B_k を持ち、各反復で B_k を rank-1 更新します。F の解析的 Jacobian が手元になくても動き、行列の更新と逆行列の更新が Sherman-Morrison で O(n^2) に収まるので、Newton の O(n^3) より反復あたりが軽くなります。収束は超線形 (~ 1.62)、Newton の 2 次より遅いですが、Jacobian が高コストなら総コストで勝ちます。

Sherman-Morrison 式は何を更新しますか？

(A + uv^T)^{-1} = A^{-1} - (A^{-1} u v^T A^{-1}) / (1 + v^T A^{-1} u) という、行列 A に rank-1 修正 uv^T を加えたときの逆行列の公式です。Broyden 法は B_k の rank-1 更新で次の Jacobian 近似を作るので、B_k^{-1} を直接 Sherman-Morrison で更新でき、毎反復 O(n^2) で済みます。

ホモトピー法 (homotopy continuation) は何のため？

初期値 x_0 が解 α から遠いと Newton 法は発散します。ホモトピー法はパラメータ t ∈ [0, 1] で簡単な系 G(x) = 0 から目標系 F(x) = 0 を連続的につなぐ H(x, t) = (1-t) G(x) + t F(x) を作り、t を 0 → 1 と少しずつ動かして解を予測子-修正子で追跡します。初期値依存性を排除し、全実解を取り出すことも可能です。多項式系の全解列挙では Bertini や HomotopyContinuation.jl 等で標準的な手法です。

多変数求根 (N-D Root Finding)

概要

連立非線形方程式 $\mathbf{F}(\mathbf{x}) = \mathbf{0}$ ($\mathbf{F}: \mathbb{R}^n \to \mathbb{R}^n$) の解を求める問題。 1 次元の Newton/secant を直接一般化した手法群と、大域収束を狙うホモトピー法に大別される。

多変数 Newton 法

$\mathbf{F}$ を $\mathbf{x}_k$ の周りで 1 次 Taylor 展開:

$$\mathbf{F}(\mathbf{x}_k + \Delta) \approx \mathbf{F}(\mathbf{x}_k) + J(\mathbf{x}_k) \Delta$$

これを $\mathbf{0}$ にする $\Delta$ で更新:

$$\mathbf{x}_{k+1} = \mathbf{x}_k - J(\mathbf{x}_k)^{-1} \mathbf{F}(\mathbf{x}_k)$$

$J(\mathbf{x}) = \partial \mathbf{F}/\partial \mathbf{x}$ はJacobian 行列 ($n \times n$):

$$J_{ij}(\mathbf{x}) = \frac{\partial F_i}{\partial x_j}(\mathbf{x})$$

収束

1 次元と同じく、単純零点 (Jacobian が正則) かつ滑らかな $\mathbf{F}$ なら二次収束:

$$\|\mathbf{x}_{k+1} - \boldsymbol{\alpha}\| \leq C \|\mathbf{x}_k - \boldsymbol{\alpha}\|^2$$

収束半径外から出発すると発散することがあり、line search やホモトピーで補強する。

sangi の newton_raphson_system(F, J, x0) はユーザーから $\mathbf{F}$ と $J$ をコールバックで受け取る。 Jacobian が解析形で得られない場合は数値微分 (中心差分) で代替できる。

線形系の解法 (LU vs QR)

$\Delta = -J^{-1} \mathbf{F}$ を計算するとき、$J^{-1}$ を陽に作るのは数値的に不安定。代わりに線形系 $J \Delta = -\mathbf{F}$ を解く。

LU 分解

$J = L U$ に分解後、前進代入と後退代入で $\Delta$ を得る。$O(n^3)$ の分解 + $O(n^2)$ の代入。 Jacobian が良条件 (条件数が小さい) なら最速。sangi の newton_raphson_system はデフォルトで部分ピボット付き LU を使う。

QR 分解

$J = QR$ に分解 ($Q$ は直交、$R$ は上三角)。LU より計算量が約 2 倍だが、条件数が大きい (Jacobian が特異に近い) ときに数値的安定性が高い。 Newton 反復の終盤で Jacobian が退化し始める場合は QR が安全。

反復解法

$n$ が大きく $J$ がスパースな場合は、GMRES や BiCGSTAB といった Krylov 部分空間反復法で $J \Delta = -\mathbf{F}$ を解く (Newton-Krylov 法)。 Jacobian を陽に作らず、行列ベクトル積 $J \mathbf{v}$ だけで動くため、メモリ・計算量とも節約できる。

Broyden 準 Newton 法

Jacobian を毎反復で計算する Newton 法は、$\mathbf{F}$ が複雑だと評価コストが支配的になる。 Broyden 法 (1965) はJacobian を近似行列 $B_k$ で置き換え、各反復で rank-1 更新する:

$$B_{k+1} = B_k + \frac{(\mathbf{y}_k - B_k \mathbf{s}_k) \mathbf{s}_k^\top}{\mathbf{s}_k^\top \mathbf{s}_k}$$

ここで $\mathbf{s}_k = \mathbf{x}_{k+1} - \mathbf{x}_k$, $\mathbf{y}_k = \mathbf{F}(\mathbf{x}_{k+1}) - \mathbf{F}(\mathbf{x}_k)$。これは「secant 条件」 $B_{k+1} \mathbf{s}_k = \mathbf{y}_k$ を満たす最小修正である。

収束

収束次数: 超線形 ($\approx 1.62$)

Newton の 2 次より遅いが、Jacobian の評価が不要。 Jacobian 評価コストが $\mathbf{F}$ の数倍以上の場合、総コストで Newton に勝つことが多い。

初期値 $B_0$ は単位行列、または有限差分による 1 回だけの Jacobian 推定を使う。

Sherman-Morrison による逆行列更新

Broyden の rank-1 更新は逆行列も $O(n^2)$ で更新できる。 Sherman-Morrison の式:

$$(A + \mathbf{u} \mathbf{v}^\top)^{-1} = A^{-1} - \frac{A^{-1} \mathbf{u} \mathbf{v}^\top A^{-1}}{1 + \mathbf{v}^\top A^{-1} \mathbf{u}}$$

これを使えば $B_{k+1}^{-1}$ を $B_k^{-1}$ から $O(n^2)$ で得られる。 Newton の毎反復 $O(n^3)$ LU 分解と比較して、$n$ が大きいときに大きな速度差が出る。

注意点: Sherman-Morrison は数値的に LU/QR より不安定なので、何度も rank-1 更新を重ねると条件数が悪化する。 sangi の broyden_method は周期的に Jacobian を再計算する (リスタート) 戦略を採用できる。

line search との組合せ

Newton ステップ $\Delta = -J^{-1} \mathbf{F}$ をそのまま採用すると、収束半径外で発散することがある。 line search は $\mathbf{x}_{k+1} = \mathbf{x}_k + \lambda \Delta$ ($\lambda \in (0, 1]$) として、メリット関数 $\|\mathbf{F}(\mathbf{x}_{k+1})\|^2 / 2$ が十分減少する $\lambda$ を選ぶ。

Armijo 条件 (十分減少条件)

$$\|\mathbf{F}(\mathbf{x}_k + \lambda \Delta)\|^2 \leq (1 - 2 c \lambda) \|\mathbf{F}(\mathbf{x}_k)\|^2$$

$c \in (0, 1/2)$ は十分減少パラメータ (典型値 $10^{-4}$)。 backtracking で $\lambda = 1, 1/2, 1/4, \ldots$ と縮めて条件を満たす最大の $\lambda$ を採用する。

line search を入れた Newton 法は、収束半径外でも発散せず徐々に収束領域に入る大域収束性を持つ (滑らかな $\mathbf{F}$ に対して)。 sangi の現状の newton_raphson_system はオプションとして line search のサポートを計画中。

ホモトピー法 (大域収束)

line search でも収束しない悪条件 (深い極小値が複数ある) では、ホモトピー法 (continuation method) を使う。

アイデア

パラメータ $t \in [0, 1]$ で簡単な系から目標系へ連続的につなぐ:

$$H(\mathbf{x}, t) = (1 - t) G(\mathbf{x}) + t \mathbf{F}(\mathbf{x})$$

$G$ は解が既知の簡単な系 (例: $G(\mathbf{x}) = \mathbf{x} - \mathbf{x}_0$ なら解は $\mathbf{x}_0$)。 $t$ を $0 \to 1$ と少しずつ動かしながら、各 $t$ での解を予測子-修正子 (predictor-corrector) で追跡する:

予測子: $\partial H/\partial t$ から接線方向の Euler ステップで次の $t$ での近似解を予測
修正子: その近似解から固定 $t$ で Newton 反復し、$H(\mathbf{x}, t) = \mathbf{0}$ を解く

$t = 1$ に到達すれば目標系の解。多項式系の全実解列挙では Bertini や HomotopyContinuation.jl のような実装で標準的に使われる。

sangi は現状ホモトピー法を公式 API として提供していないが、多変数 Newton を呼び出す上位ラッパーとして実装可能。

収束半径と初期値選択

多変数 Newton/Broyden は局所収束のみ保証され、収束半径の外から出発すると発散する。 Kantorovich の定理は収束半径の十分条件を与える: 初期点 $\mathbf{x}_0$ で

$$\|J(\mathbf{x}_0)^{-1}\| \cdot \|\mathbf{F}(\mathbf{x}_0)\| \cdot L \leq \frac{1}{2}$$

($L$ は Jacobian の Lipschitz 定数) なら Newton が二次収束する解が $\mathbf{x}_0$ の近傍に存在する。実用上は次のような戦略を組み合わせる:

物理的・問題ドメインの知識: 解の存在範囲が事前に分かっているなら、その中心を $\mathbf{x}_0$ にする
連続化: 関連する easy な問題から徐々にパラメータを動かしていく (ホモトピーの簡単版)
line search / trust region: 発散を抑え、徐々に収束領域に近づける
multistart: 複数の初期値から並列に走らせ、収束したものを採用

初期値が良ければ Newton の 2 次収束は 5〜10 反復で機械精度に達する。反復が 50 回を超えても収束しない場合は、初期値またはアルゴリズム選択を疑うべき。

設計判断

Newton と Broyden を別関数として公開する理由

多変数求根の使い分けは、Jacobian の評価コストに強く依存する:

Jacobian が解析的に書ける + 評価コストが $\mathbf{F}$ の数倍以内 → Newton (2 次収束で速い)
Jacobian が複雑または不明 → Broyden (Jacobian 不要)

sangi は両者を newton_raphson_system と broyden_method として独立に公開し、ユーザーが問題の特性で選べるようにしている。将来的に Newton-Krylov や trust region を追加する場合も同じ命名規約で並べる予定。

1D Newton との一貫性

1D Newton (newton_raphson) と多変数 Newton (newton_raphson_system) は同じ概念の異なる次元への一般化。 sangi は両者でほぼ同じ API 形式 (関数 + 導関数/Jacobian + 初期値) を採用し、コード移行のコストを最小化している。

参考文献

Ortega, J. M. and Rheinboldt, W. C. (1970). Iterative Solution of Nonlinear Equations in Several Variables. Academic Press.
Broyden, C. G. (1965). "A class of methods for solving nonlinear simultaneous equations". Mathematics of Computation, 19(92), 577–593.
Allgower, E. L. and Georg, K. (2003). Introduction to Numerical Continuation Methods. SIAM.
Dennis, J. E. and Schnabel, R. B. (1996). Numerical Methods for Unconstrained Optimization and Nonlinear Equations. SIAM.