方程式を解くためのアルゴリズム

方程式の解法の定義

n 個の一連の非線形関数 F_i(x) があるとします。ここで、n はベクトル x の要素数であり、方程式の求解のゴールはすべての F_i(x) が 0 になるようなベクトル x を見つけることです。

fsolve は、要素の二乗和を最小にすることにより方程式系を解こうとします。二乗和がゼロの場合、方程式系は解かれます。fsolve には次の 3 つのアルゴリズムがあります。

信頼領域法
Trust-region-dogleg 法
レーベンバーグ・マルカート法

すべてのアルゴリズムは大規模です。大規模アルゴリズムと中規模アルゴリズムを参照してください。

関数 fzero は単一の 1 次元方程式を解きます。

関数 mldivide は線形方程式系を解きます。

信頼領域法アルゴリズム

Optimization Toolbox™ のソルバーで使用される多くのメソッドは、"信頼領域法" を基にしています。信頼領域法はシンプルなものですが最適化では重要な概念です。

最適化の信頼領域法のアプローチを理解するために制約なし最小化問題を考え、f(x) を最小化します。ここで関数はベクトル引数を取り、スカラーを出力します。現在の点を n 空間の点 x と想定し、より小さい関数値をもつ点へ移動して最小化を行う場合を考えてみましょう。これを行うために、このアルゴリズムでは、シンプルな関数 q で f を近似します。この関数は、点 x の近傍 N で関数 f の挙動をよく表すものです。この近傍が信頼領域です。ソルバーは、テストステップ s を N における最小化 (または近似最小化) によって計算します。信頼領域の部分問題を次に示します。

$\min_{s} {q (s), s \in N} .$

f(x + s) < f(x) の場合、ソルバーは現在の点を x + s に更新します。そうでない場合は現在の点は変更されず、ソルバーは N (信頼領域) を縮小し、テストステップの計算を繰り返します。

f(x) を最小化するための信頼領域を決める上で重要な問題は、近似 q (現在の点 x で定義) の選択および計算方法、信頼領域 N の選択および変更方法、信頼領域の部分問題を解く精度です。

標準的な信頼領域法 ([48]) では二次近似 q は、x における F のテイラー近似のはじめの 2 項によって決められます。近傍 N は通常、球形または楕円体です。数学的に、信頼領域の部分問題は、次のように表現できます。

\min {\frac{1}{2} s^{T} H s + s^{T} g such that ‖ D s ‖ \leq Δ},

(1)

ここで g は現在の点 x における f の勾配です。H はヘッセ行列 (2 次導関数の対称行列) です。D は対角スケーリング行列であり、Δ は正のスカラーです。‖ . ‖ は 2 ノルムです。式 1を解くために、アルゴリズム ([48] を参照) で H のすべての固有値を計算し、次の永年方程式にニュートン過程を適用できます。

$\frac{1}{Δ} - \frac{1}{‖ s ‖} = 0.$

このようなアルゴリズムにより式 1の正確な解が与えられます。しかし、H の因子分解に比例して時間が必要になります。このため、信頼領域の問題では異なったアプローチが必要となります。式 1に基づく近似とヒューリスティックな方法は文献 ([42] と [50]) で提案されています。Optimization Toolbox のソルバーは、信頼領域法の部分問題を 2 次元の部分空間 S ([39] と [42]) に制限する近似アプローチをサポートしています。ソルバーが部分空間 S を計算した後は、式 1を解く作業が簡単になります。部分空間では問題は 2 次元になるためです。ここでの主な作業は、部分空間を決定することになります。

ソルバーは、前処理付き共役勾配法 (次のセクションで説明) を用いて、2 次元の部分空間 S を決めます。ソルバーは s₁ と s₂ で決まる線形空間として S を定義します。ここで、s₁ は勾配 g の方向にあります。s₂ は近似ニュートン方向、すなわち以下の解

$H \cdot s_{2} = - g,$

または負の曲率の方向です。

$s_{2}^{T} \cdot H \cdot s_{2} < 0.$

このように S を選択する背景の考え方は、最急降下方向または負の曲率方向にグローバルな収束を進めることと、ニュートンステップが存在する場合は、これを介して迅速にローカルな収束を達成することです。

信頼領域法のアプローチを使用した制約なしの最小化のプロセスは以下になります。

2 次元の信頼領域の部分問題の定式化
テストステップ s を決めるため、式 1 を解きます。
f(x + s) < f(x) の場合、x = x + s とします。
Δ を調節します。

ソルバーは収束するまでこの 4 つのステップを繰り返し、信頼領域の大きさ Δ を標準的な規則に基づいて調整します。特に、ソルバーはテストステップを受け入れない場合、f(x + s) ≥ f(x) であれば、信頼領域のサイズを小さくします。詳細については、[46]と[49]を参照してください。

Optimization Toolbox のソルバーは特定の関数 f の重要なケースを扱います。非線形最小二乗法、二次関数、線形最小二乗法を考えてみましょう。しかし、根底に存在するアルゴリズムは、一般的な場合と同じです。

前処理付き共役勾配法

線形方程式系 Hp = –g の大きな対称正定値システムを解く一般的な方法は、前処理付き共役勾配法 (PCG) です。この反復法は、形式 H·v の行列ベクトル積を計算する機能を必要とします。ここで v は任意のベクトルです。対称正定値行列 M は H の "前提条件子" です。すなわち、M = C² です。ここで、C^–1HC^–1 は、条件数の良い行列または分類分けされた固有値をもつ行列です。

最小化の過程で、ヘッセ行列 H が対称と仮定します。しかし、H は、強い最小化子の近傍の中でのみ正定値であることが保証されます。PCG アルゴリズムは、負または 0 の曲率方向が検出された場合に終了します (d^THd ≤ 0)。PCG 出力方向 p は、負の曲率方向またはニュートンシステム Hp = –g への近似解のどちらかです。いずれにせよ、p は信頼領域法のアプローチで使用される 2 次元部分空間を定義するために役立ちます (非線形最小化に対する信頼領域法を参照)。

Trust-region-dogleg 法アルゴリズム

他の方法としては線形方程式系を解いて探索方向を求める方法があります。ニュートン法で次のような探索方向 d_k を求めます。

J(x_k)d_k = –F(x_k)
x_{k + 1} = x_k + d_k,

ここで、J(x_k) は n 行 n 列のヤコビアンです。

$J (x_{k}) = [\begin{matrix} \nabla F_{1} {(x_{k})}^{T} \\ \nabla F_{2} {(x_{k})}^{T} \\ ⋮ \\ \nabla F_{n} {(x_{k})}^{T} \end{matrix}] .$

ニュートン法は問題が生じることがあります。J(x_k) が特異なためにニュートンステップ d_k が定義されない可能性があります。また、正確なニュートンステップ d_k の計算に時間がかかる可能性もあります。さらに、ニュートン法では開始点が解とかけ離れている場合は収束しない可能性があります。

信頼領域法 (非線形最小化に対する信頼領域法で説明) は、J(x_k) が特異な場合にも対応でき、また開始点が解とかけ離れている場合のロバスト性を向上させます。信頼領域法を使用するには、x_k よりも x_{k + 1} の方が良いかどうかを判断するためのメリット関数が必要です。次の選択も可能です。

$\min_{d} f (d) = \frac{1}{2} F {(x_{k} + d)}^{T} F (x_{k} + d) .$

但し、 f(d) の最小値が必ずしも F(x) の根にはなりません。

ニュートンステップ d_k は、

M(x_k + d) = F(x_k) + J(x_k)d,

の根であるため、m(d) の最小値にもなります。ここで

\begin{matrix} \min_{d} m (d) = \frac{1}{2} {‖ M (x_{k} + d) ‖}_{2}^{2} = \frac{1}{2} {‖ F (x_{k}) + J (x_{k}) d ‖}_{2}^{2} \\ = \frac{1}{2} F {(x_{k})}^{T} F (x_{k}) + d^{T} J {(x_{k})}^{T} F (x_{k}) + \frac{1}{2} d^{T} J {(x_{k})}^{T} J (x_{k}) d . \end{matrix}

(2)

メリット関数の選択としては m(d) の方が f(d) よりも優れており、Trust-Region の部分問題は

\min_{d} [\frac{1}{2} F {(x_{k})}^{T} F (x_{k}) + d^{T} J {(x_{k})}^{T} F (x_{k}) + \frac{1}{2} d^{T} J {(x_{k})}^{T} J (x_{k}) d],

(3)

‖D·d‖ ≤ Δ が条件になります。この部分問題は dogleg 法を使用して効率的に解くことができます。

信頼領域法の概要は Conn [4] と Nocedal [31] を参照してください。

Trust-region-dogleg 法の実装

この Trust-region-dogleg 法アルゴリズムの最大の特徴は Powell の dogleg 手順を使用して式 3を最小化する手順 d を計算する点にあります。詳細については、Powell [34]を参照してください。

このアルゴリズムは、コーシーステップ (最急降下方向に沿ったステップ) と f(x) に対応するガウス・ニュートンステップを凸に組み合わせることで、ステップ d を構成します。コーシーステップは、次のように計算されます。

d_C = –αJ(x_k)^TF(x_k),

ここで、α は式 2を最小化します。

ガウス・ニュートンステップは

J(x_k)·d_GN = –F(x_k),

MATLAB^® の mldivide (行列の左除算) 演算子を使って上の式を解くことで計算されます。

アルゴリズムは、次のようにステップ d を選択します。

d = d_C + λ(d_GN – d_C),

ここで λ は ‖d‖ ≤ Δ を満たすような区間 [0,1] の最大値です。これは J_k が (ほぼ) 特異である場合、d はコーシーの方向と一致します。

Trust-region-dogleg 法アルゴリズムが効率的なのは、(ガウス・ニュートンステップ計算のための) 反復の各回で、1 つの線形解しか必要としないからです。また、場合によっては、直線探索を伴うガウス・ニュートン法を使用する方法よりもこのアルゴリズムの方がロバストです。

レーベンバーグ・マルカート法

レーベンバーグ・マルカート法アルゴリズム ([25] および [27]) は、探索方向として次の線形方程式の解を使います。

(J {(x_{k})}^{T} J (x_{k}) + λ_{k} I) d_{k} = - J {(x_{k})}^{T} F (x_{k}),

(4)

またはオプションとして以下の方程式を使います。

(J {(x_{k})}^{T} J (x_{k}) + λ_{k} diag (J {(x_{k})}^{T} J (x_{k}))) d_{k} = - J {(x_{k})}^{T} F (x_{k}),

(5)

ここでスカラー λ_k は d_k の大きさと方向を共にコントロールします。式 4を使用するには、fsolve のオプション ScaleProblem を 'none' に設定します。式 5を使用するには、このオプションを 'jacobian' に設定します。

λ_k がゼロのとき、方向 d_k はガウス・ニュートン法です。λ_k が無限大方向になると、d_k はゼロベクトル方向へ向かい、最急降下方向になります。これは、一部の十分に大きい λ_k に対して、項 F(x_k + d_k) < F(x_k) が成り立つことを意味します。そのため、このアルゴリズムは、ガウス・ニュートン法の効果を制限する二次の項が存在するときでさえ、減少を保証するように項 λ_k をコントロールできます。そのため、レーベンバーグ・マルカート法アルゴリズムはガウス・ニュートンの方向と最急降下の方向の間の探索方向を使います。詳細については、最小二乗に関するドキュメンテーションのレーベンバーグ・マルカート法を参照してください。