最新のリリースでは、このページがまだ翻訳されていません。 このページの最新版は英語でご覧になれます。

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

方程式の解法の定義

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

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

  • 信頼領域法

  • Trust-region-dogleg 法

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

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

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

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

信頼領域 fsolve アルゴリズム

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

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

mins{q(s), sN}.(1)

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

f(x) を最小化するための信頼領域を決める上で重要な問題は、近似 q (現在の点 x で定義) の選択および計算方法、信頼領域 N の選択および変更方法、信頼領域の部分問題を解く精度です。この節では制約なしの問題に焦点をあてます。変数上に制約があるために複雑度が増す問題については、後節で取り扱います。

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

min{12sTHs+sTg  such that  DsΔ},(2)

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

1Δ1s=0.

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

2 次元の部分空間 S は、以下に示す前処理を使用した共役勾配法過程を用いて決められます。ソルバーは s1 と s2 で決まる線形空間として S を定義します。ここで、s1 は勾配 g の方向にあります。s2 は近似ニュートン方向、すなわち以下の解

Hs2=g,(3)

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

s2THs2<0.(4)

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

信頼領域法の概念を使用した制約なしの最小化の概要は以下になります。

  1. 2 次元の信頼領域の部分問題の定式化

  2. テスト ステップ s を決めるため、式 2 を解きます。

  3. f(x + s) < f(x) の場合、x = x + s とします。

  4. Δ を調節します。

これら 4 つのステップは、収束するまで繰り返されます。信頼領域の大きさ Δ は標準的な規則に基づいて調整されます。特に、使用するステップが適用されない場合、すなわち f(x + s) ≥ f(x) の場合、内点集合は小さくなります。詳細は [46][49] を参照してください。

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

前処理付き共役勾配法

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

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

Trust-region-dogleg 法

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

J(xk)dk = –F(xk)
xk + 1 = xk + dk,

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

J(xk)=[F1(xk)TF2(xk)TFn(xk)T].

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

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

mindf(d)=12F(xk+d)TF(xk+d).

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

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

M(xk + d) = F(xk) + J(xk)d

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

mindm(d)=12M(xk+d)22=12F(xk)+J(xk)d22=12F(xk)TF(xk)+dTJ(xk)TF(xk)+12dTJ(xk)TJ(xk)d.(5)

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

mind[12F(xk)TF(xk)+dTJ(xk)TF(xk)+12dTJ(xk)TJ(xk)d],(6)

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

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

Trust-region-dogleg 法の実装

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

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

dC = –αJ(xk)TF(xk)

ここで、α は 式 5 を最小化するように選択されます。

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

J(xk)·dGN = –F(xk),

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

ステップ d は次の式が成り立つように選択されます。

d = dC + λ(dGN – dC)

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

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

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

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

(J(xk)TJ(xk)+λkI)dk=J(xk)TF(xk),(7)

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

(J(xk)TJ(xk)+λkdiag(J(xk)TJ(xk)))dk=J(xk)TF(xk),(8)

ここでスカラー λk は dk の大きさと方向を共にコントロールします。式 7 を選択するにはオプション ScaleProblem'none' に設定します。式 8 を選択するには ScaleProblem'Jacobian' に設定します。

λk がゼロのとき、方向 dk はガウス・ニュートン法です。λk が無限大方向になると、dk はゼロ ベクトル方向へ向かい、最急降下方向になります。これは、一部の十分に大きい λk に対して、項 F(xk + dk) <  F(xk) が成り立つことを意味します。そのため λk はガウス・ニュートン法の効果を制限する二次の項が存在するときでさえ、減少を保証するようにコントロールされます。そのため、レーベンバーグ・マルカート法はガウス・ニュートンの方向と最急降下の方向の間の探索方向を使います。

\ アルゴリズム

このアルゴリズムは mldivide についての MATLAB 算術演算子の節で説明されています。

fzero アルゴリズム

fzero はスカラー変数 x のスカラー関数 f の根を見つけようとします。

fzero は f(x) によって符号が変化する初期点のまわりの区間を調べます。初期点ではなく初期区間を与える場合、fzero は f(x) が区間の端点で符号を変化させることを確認します。初期区間は有限である必要があります。すなわち、±Inf は含めることができません。

fzero は f(x) の根を探すために、区間の 2 分割、線形内挿、逆二次内挿を組み合わせて使用します。詳細は、fzero を参照してください。