最小二乗 (モデル当てはめ) アルゴリズム

最小二乗の定義

一般的に最小二乗は二乗和の関数を極小にするベクトル x を見つける問題です。場合によってはいくつかの制約をもちます。

$\min_{x} {‖ F (x) ‖}_{2}^{2} = \min_{x} \sum_{i} F_{i}^{2} (x)$

A·x ≤ b, Aeq·x = beq, lb ≤ x ≤ ubが条件になります。

いくつかの Optimization Toolbox™ のソルバーは、さまざまなタイプの F(x) と制約に利用できます。

ソルバー	F(x)	制約
`mldivide`	C·x – d	なし
`lsqnonneg`	C·x – d	x ≥ 0
`lsqlin`	C·x – d	範囲、線形
`lsqnonlin`	一般 F(x)	範囲
`lsqcurvefit`	F(x, xdata) – ydata	範囲

mldivideで使用されるアルゴリズムに加えて、Optimization Toolbox のソルバーには 5 つの最小二乗法アルゴリズムがあります。

lsqlin内点法
lsqlin有効制約法
信頼領域 Reflective 法 (非線形または線形最小二乗法)
レーベンバーグ・マルカート (非線形最小二乗法)
lsqnonnegに使用されるアルゴリズム

lsqlin有効制約法を除くすべてのアルゴリズムは大規模です。大規模アルゴリズムと中規模アルゴリズムを参照してください。非線形最小二乗法の一般的な調査結果については、Dennis[8]を参照してください。レーベンバーグ・マルカート法の詳細については、Moré[28]を参照してください。

線形最小二乗法:内点法または有効制約法

lsqlin'interior-point'アルゴリズムはinterior-point-convex quadprog アルゴリズムを使用し、lsqlin'active-set'アルゴリズムはアクティブセットquadprogアルゴリズムを使用します。quadprogの問題の定義は、次の二次関数を最小化することです。

$\min_{x} \frac{1}{2} x^{T} H x + c^{T} x$

ここで線形制約と範囲制約が条件になります。関数lsqlinは、線形制約と範囲制約に従いながら、ベクトル Cx – d の 2 乗ノルムを最小化します。言い換えれば、lsqlinは次を最小化します。

$\begin{matrix} {‖ C x - d ‖}_{2}^{2} = {(C x - d)}^{T} (C x - d) \\ = (x^{T} C^{T} - d^{T}) (C x - d) \\ = (x^{T} C^{T} C x) - (x^{T} C^{T} d - d^{T} C x) + d^{T} d \\ = \frac{1}{2} x^{T} (2 C^{T} C) x + {(- 2 C^{T} d)}^{T} x + d^{T} d . \end{matrix}$

これは行列 H を 2C^TC、c ベクトルを(–2C^Td)に設定すると、quadprogフレームワークに適合します (加法項 d^Td は、最小値の位置に影響しません)。lsqlin問題をこのように定式化した後、quadprogは解を計算します。

メモ

quadprog'interior-point-convex'アルゴリズムには 2 つのコードパスがあります。一方はヘッセ行列 H が通常 (フル) の double 行列の場合に使用され、もう一方は H がスパース行列の場合に使用されます。スパースデータ型の詳細については、スパース行列を参照してください。一般に、H をsparseとして指定すると、アルゴリズムは、比較的非ゼロの項が少ない大規模な問題の場合により高速になります。同様に、H をfullとして指定すると、アルゴリズムは、小規模または比較的密な問題の場合により高速になります。

信頼領域 Reflective 法の最小二乗

信頼領域 Reflective 法の最小二乗アルゴリズム

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

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

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

(1)

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 Δ},

(2)

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

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

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

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

H \cdot s_{2} = - g,

(3)

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

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

(4)

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

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

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

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

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

大規模な非線形最小二乗法

f(x) の重要で特別なケースは、非線形最小二乗問題です。

\min_{x} \sum_{i} f_{i}^{2} (x) = \min_{x} {‖ F (x) ‖}_{2}^{2},

(5)

ここで F(x) は f_i(x) に等しい F(x) の要素 i を使ったベクトル値関数です。この問題を解くために使われる基本的な方法は非線形最小化に対する信頼領域法に記述されている一般的な場合と同じです。しかし、非線形最小二乗問題の構造は、効率性を強調したものです。特に、近似的なガウス・ニュートンの方向、すなわち以下の解 s は

\min {‖ J s + F ‖}_{2}^{2},

(6)

2 次元部分空間 S を定義するのに使います (ここで、J は F(x) のヤコビアンです)。要素関数 f_i(x) の 2 次導関数は使われません。

各反復で前処理付き共役勾配法は次の正規方程式を近似的に解くために使われます。

$J^{T} J s = - J^{T} F,$

ここで、正規方程式は、明示的な形をしていません。

大規模な線形最小二乗法

この場合、解く関数 f(x) は以下になります。

$f (x) = {‖ C x + d ‖}_{2}^{2},$

場合によっては線形制約に従います。アルゴリズムは、ある範囲内で局所的な解に収束する厳密に実行可能な反復を行います。各反復は大規模線形システム (n次元: n は x の長さ) の近似解を含みます。反復行列は行列 C の構造をもっています。特に前処理付き共役勾配法は、正規方程式を、近似的に解くために使われます。

$C^{T} C x = - C^{T} d,$

ここで、正規方程式は、明示的な形をしていません。

部分空間の信頼領域法が探索方向を決定するのに使われます。ただし、ステップを非線形最小化の場合と同じく (場合によっては) 1 つの反射ステップに制限する代わりに、区分的な Reflective ライン探索が二次型の場合と同じく各反復で行われます。ライン探索の詳細は[45]を参照してください。一般に、線形システムは、解の点で 1 次の最適性条件をもつニュートンアプローチを表しています。すなわち、非常に強い局所的な収束率をもつことになります。

ヤコビ乗算関数-lsqlinは、行列 C を陽的に使用せずに線形制約付き最小二乗問題を解くことができますが、代わりにユーザーが与えたヤコビ乗算関数jmfun

W = jmfun(Jinfo,Y,flag)

を使用します。この関数は行列 Y に対して以下の乗算を計算しなければなりません。

flag == 0の場合、W = C'*(C*Y)です。
flag > 0の場合、W = C*Yです。
flag < 0の場合、W = C'*Yです。

これは C が大規模な場合は有用ですが、C を陽的に作成せずにjmfunを記述できるような構造が必要になります。例については、線形最小二乗付きヤコビ乗算関数を参照してください。

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

最小二乗問題は、二乗和である関数 f(x) を最小化します。

\min_{x} f (x) = {‖ F (x) ‖}_{2}^{2} = \sum_{i} F_{i}^{2} (x) .

(7)

この種の問題は、特に、非線形パラメーター推定のようにモデル関数をデータに当てはめる実用的応用でよく見られます。このような問題は、目的が出力 y(x,t) をベクトル x とスカラー t の連続モデル軌跡 φ(t) に乗せる制御システムでも見られます。この問題は、次のように表現されます。

\min_{x \in ℜ^{n}} \int_{t_{1}}^{t_{2}} {(y (x, t) - φ (t))}^{2} d t,

(8)

ここで y(x,t) と φ(t) はスカラー関数です。

近似値を求めるために積分を離散化します。

\min_{x \in ℜ^{n}} f (x) = \sum_{i = 1}^{m} {(y (x, t_{i}) - φ (t_{i}))}^{2},

(9)

ここで、t_iは等間隔です。この問題では、ベクトル F(x) が次のようになります。

$F (x) = [\begin{matrix} y (x, t_{1}) - φ (t_{1}) \\ y (x, t_{2}) - φ (t_{2}) \\ ... \\ y (x, t_{m}) - φ (t_{m}) \end{matrix}] .$

この種の問題では、残差∥F(x)∥が最適条件で小さくなる傾向があります。これは、現実的に到達可能な目標軌道を設定することが一般的だからです。制約なしの最適化の基礎で説明されているように、制約なしの最小化のための一般的な手法を用いて式 7の関数を最小化できますが、問題の特性によっては解法手順を繰り返す効率が改善される場合があります。式 7の勾配とヘッシアン行列は特別な構造をもちます。

F(x) の m 行 n 列のヤコビ行列を J(x) とし、f(x) の勾配ベクトルを G(x) とし、f(x) のヘッシアン行列を H(x) とし、各 F_i(x) のヘッシアン行列を D_i(x) とすると、次のようになります。

\begin{array}{l} G (x) = 2 J (^{x) T} F (x) \\ H (x) = 2 J (^{x) T} J (x) + 2 Q (x), \end{array}

(10)

ここで、

$Q (x) = \sum_{i = 1}^{m} F_{i} (x) \cdot D_{i} (x) .$

行列 Q(x) には次のような特性があります。x_kが解に近づくと残差∥F(x)∥が 0 に近づく傾向があるのと同時に、Q(x) も 0 に近づく傾向があります。そのため、∥F(x)∥が解で小さくなる場合は、ガウス・ニュートンの方向を最適化手順の基礎として使用する方法が効果的です。

ガウス・ニュートン法では、線形最小二乗問題の解である探索方向 d_kが大反復 k で得られます。

\min_{d_{k} \in ℜ^{n}} {‖ J (x_{k}) d_{k} + F (x_{k}) ‖}_{2}^{2} .

(11)

この方法で導かれる方向は、項Q(x) = 0のときのニュートン方向と同じです。アルゴリズムは、探索方向 d_kを直線探索戦略の一部として使用して、関数 f(x) が各反復で減少することを確保できます。

ガウス・ニュートン法は、2 次の項 Q(x) が無視できない場合に、問題が発生することがあります。レーベンバーグ・マルカート法がこの問題を解決します。

レーベンバーグ・マルカート法 ([25]と[27]を参照) は、次の線形方程式の解である探索方向を使用します。

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

(12)

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

(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}),

(13)

ここで、スカラー λ_kは d_kの大きさと方向の両方を制御し、diag(A)はAの対角項の行列を意味します。式 12を選択するにはオプションScaleProblemを'none'に設定します。式 13を選択するにはScaleProblemを'Jacobian'に設定します。

パラメーター λ₀の初期値は、ユーザーがInitDampingオプションで設定します。状況によっては、このオプションの既定値0.01が適していないこともあります。レーベンバーグ・マルカート法アルゴリズムでの初期の進行状況が良くない場合は、InitDampingを1e2などの既定値以外の値に設定してみてください。

λ_kがゼロのとき、方向 d_kはガウス・ニュートン法の方向と一致します。λ_kが無限大方向になると,d_kはゼロベクトル方向へ向かい、最急降下方向になります。その結果、一部の十分に大きい λ_kに対して、項F(x_k+ d_k) < F(x_k)が成り立ちます。したがって、アルゴリズムがガウス・ニュートン法の効率を制限する 2 次の項に遭遇した場合でも減少するように項 λ_kを制御できます。ステップが成功した (つまり関数値が低くなった) 場合、アルゴリズムは λ_k+1= λ_k/10 と設定します。ステップが失敗した場合、アルゴリズムは λ_k+1= λ_k*10 と設定します。

内部的に、レーベンバーグ・マルカート法アルゴリズムは、関数の許容誤差の1e-4倍となる最適性の許容誤差 (停止条件) を使用します。

そのため、レーベンバーグ・マルカート法はガウス・ニュートンの方向と最急降下方向の間で交差する探索方向を使用します。

レーベンバーグ・マルカート法のもう 1 つのメリットは、ヤコビアン J がランク落ちしている場合です。この場合は、ガウス・ニュートン法で数値問題が発生する可能性があります。これは、式 11での最小化問題が不良設定されるためです。対照的に、レーベンバーグ・マルカート法は、各反復でフルランクをもつため、このような問題が発生しません。

以下の図は、ローゼンブロックの関数の最小化 (最小二乗形式の、難しいことで有名な最小化問題) におけるレーベンバーグ・マルカート法の反復を示しています。

Rosenbrock 関数上のレーベンバーグ・マルカート法

反復点を生成するスクリプトを含めたこの図の詳細については、バナナ関数の最小化を参照してください。

レーベンバーグ・マルカート法の範囲制約

問題に範囲制約が含まれている場合は、lsqcurvefitとlsqnonlinがレーベンバーグ・マルカート法の反復を変更します。提案された反復点 x が範囲外に存在する場合は、アルゴリズムがステップを最も近い実行可能点に投影します。つまり、実行不可能点を実行可能領域に投影する"射影"演算子として定義された P を使用して、アルゴリズムが提案点 x を P(x) に変更します。定義により、射影演算子 P は、以下に従って、成分 x_iごとに個別に作用します。

$P (x_{i}) = {\begin{array}{l} {lb}_{i} & if x_{i} < {lb}_{i} \\ {ub}_{i} & if x_{i} > {ub}_{i} \\ x_{i} & otherwise \end{array}$

または

$P (x_{i}) = \min (\max (x_{i}, {lb}_{i}), {ub}_{i}) .$

アルゴリズムが 1 次の最適性の尺度の停止条件を変更します。x を提案された反復点とします。制約なしの場合は、停止条件が次のようになります。

{‖ \nabla f (x) ‖}_{\infty} \leq tol,

(14)

ここで、tol は、1e-4*FunctionToleranceである最適性の許容誤差値です。有界の場合は、停止条件が次のようになります。

{‖ x - P (x - \nabla f (x)) ‖}_{\infty}^{2} \leq tol {‖ \nabla f (x) ‖}_{\infty} .

(15)

この条件を理解するために、まず、x が実行可能領域の内部に存在する場合は、演算子 P が機能しないことに注意してください。そのため、停止条件は次のようになります。

${‖ x - P (x - \nabla f (x)) ‖}_{\infty}^{2} = {‖ \nabla f (x) ‖}_{\infty}^{2} \leq tol {‖ \nabla f (x) ‖}_{\infty},$

これは、元の制約なしの停止条件 ${‖ \nabla f (x) ‖}_{\infty} \leq tol$ と同じです。境界制約が有効な場合、つまり、x – ∇f(x)が実行可能でない場合は、アルゴリズムが停止すべき点で、境界上の特定の点での勾配が境界に対して垂直になります。したがって、次の図に示すように、点 x は、最急降下ステップの投影であるP(x – ∇f(x))と一致します。

レーベンバーグ・マルカート法の停止条件

Sketch of x minus the projection of x minus gradient of f(x)

参考

quadprog|lsqlin|lsqcurvefit|lsqnonlin|lsqnonneg