Main Content

このページの内容は最新ではありません。最新版の英語を参照するには、ここをクリックします。

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

最小二乗の定義

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

minxF(x)22=minxiFi2(x)

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

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

ソルバーF(x)制約
mldivideC·x – dなし
lsqnonnegC·x – dx ≥ 0
lsqlinC·x – d範囲、線形
lsqnonlin一般 F(x)一般的な線形と非線形
lsqcurvefitF(x, xdata) – ydata一般的な線形と非線形

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

  • lsqlin 内点法

  • lsqlin 有効制約法

  • 信頼領域 Reflective 法 (非線形または線形最小二乗法、範囲制約)

  • レーベンバーグ・マルカート (非線形最小二乗法、範囲制約)

  • 非線形最小二乗ソルバー lsqnonlin および lsqcurvefit 用に変更された fmincon 'interior-point' アルゴリズム (一般的な線形制約と非線形制約)。

  • lsqnonneg に使用されるアルゴリズム

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

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

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

minx12xTHx+cTx

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

Cxd22=(Cxd)T(Cxd)=(xTCTdT)(Cxd)=(xTCTCx)(xTCTddTCx)+dTd=12xT(2CTC)x+(2CTd)Tx+dTd.

これは行列 H を 2CTC、c ベクトルを (–2CTd) に設定すると、quadprog フレームワークに適合します (加法項 dTd は、最小値の位置に影響しません)。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 における最小化 (または近似最小化) によって計算されます。信頼領域の部分問題を次に示します。

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 の重要かつ特別なケースをいくつか扱います。非線形最小二乗法、二次関数、線形最小二乗法を考えてみましょう。しかし、根底に存在するアルゴリズムは、一般的な場合と同じです。これらの特別な場合は後続のセクションで説明します。

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

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

minxifi2(x)=minxF(x)22,(5)

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

minJs+F22,(6)

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

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

JTJs=JTF,

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

大規模な線形最小二乗法

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

f(x)=Cx+d22,

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

CTCx=CTd,

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

部分空間の信頼領域法が探索方向を決定するのに使われます。ただし、ステップを非線形最小化の場合と同じく (場合によっては) 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) を最小化します。

minxf(x)=F(x)22=iFi2(x).(7)

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

minxnt1t2(y(x,t)φ(t))2dt,(8)

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

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

minxnf(x)=i=1m(y(x,ti)φ(ti))2,(9)

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

F(x)=[y(x,t1)φ(t1)y(x,t2)φ(t2)...y(x,tm)φ(tm)].

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

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

G(x)=2J(x)TF(x)H(x)=2J(x)TJ(x)+2Q(x),(10)

ここで、

Q(x)=i=1mFi(x)Di(x).

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

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

mindknJ(xk)dk+F(xk)22.(11)

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

P(xi)={lbiif xi<lbiubiif xi>ubixiotherwise

または

P(xi)=min(max(xi,lbi),ubi).

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

f(x)tol,(14)

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

xP(xf(x))2tolf(x).(15)

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

xP(xf(x))2=f(x)2tolf(x),

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

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

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

制約付き最小二乗法に合わせて変更された fmincon アルゴリズム

線形制約および非線形制約を処理するために、lsqnonlin および lsqcurvefit ソルバーの 'interior-point' アルゴリズムが、変更されたバージョンの fmincon 'interior-point' アルゴリズムを内部的に呼び出します。変更されたバージョンは、fmincon を実行する際に使用される近似ではなく、別のヘッセ近似を使用します。既定では、fmincon は BFGS ヘッセ近似を使用します。このセクションの後半で、最小二乗ソルバーが使用する別のヘッセ近似を確認できます。fmincon 'interior-point' アルゴリズムとそのヘッセ近似の詳細については、fmincon の内点法アルゴリズムを参照してください。

最小二乗のコンテキストでは、目的関数 f(x) はスカラーではなくベクトルです。最小化する二乗和は以下です。

fT(x)f(x).

解決する問題は以下です。

minxf(x)22

条件

c(x) ≤ 0

ceq(x) = 0.

この問題のラグランジュ関数は以下です。

L(x,λE,λI)=f(x)22+λETceq(x)+λITc(x).

ラグランジュ関数の x に対する勾配は以下です。

xL(x,λE,λI)=2FT(x)f(x)+AET(x)λE+AIT(x)λI.

ここで、F は目的関数のヤコビアンです。

Fi,j(x)=fi(x)xj,

また、AE と AI はそれぞれ、等式制約関数 ceq(x) と不等式制約関数 c(x) のヤコビアンです。

この問題のカルーシュ・キューン・タッカー (KKT) 条件は以下です。

2FT(x)f(x)+AET(x)λE+AIT(x)λI=0SλIμe=0ceq(x)=0c(x)+s=0.

ここで、変数 s と λI は非負、S = diag(s) であり、μ は反復回数が増えるとゼロに近付く範囲パラメーターであり、e は 1 のベクトルです。

これらの方程式を解くために、ソフトウェアは反復処理を実行します。反復回数 k を呼び出します。このアルゴリズムは、次の更新方程式を使用して、ステップ Δx、Δs、ΔλE、ΔλI を取ることで方程式を解こうとします。

[H0AETAIT0Sdiag(λi)0SAE000AIS00][ΔxS1ΔsΔλEΔλI]=[2FT(x)f(x)+AET(x)λE+AIT(x)λISλIμeceq(x)c(x)+s].

ここで、H は目的関数のヘッシアンです。BFGS ヘッセ近似を使用すると、ステップ k + 1 におけるヘッセ近似 Bk+1 は以下になります。

Bk+1=BkBkskskTBkskTBksk+ykykTykTsk,(16)

ここで、

sk=xk+1xkyk=xL(xk+1,λk+1)xL(xk,λk+1).

このヘッセ近似の初期値は B0 = I です。

fmincon におけるこのアプローチの詳細については、ヘッセ行列の更新を参照してください。

最小二乗 BFGS ヘッシアンの更新は、fmincon ヘッシアンの更新と比較して以下の違いがあります。最小二乗目的関数のヘッシアンは以下です。

xx2f(x)2=2FT(x)F(x)+2ifi(x)xx2fi(x).

最小二乗ヘッシアンの更新は、初項である 2FTF を個別の量として使用し、2FTF の項を除いた BFGS 式を使用してヘッセ近似を更新します。言い換えれば、BFGS ヘッシアンの更新は、Bk のためのヘッシアンの更新式式 16を使用します。ここで、ラグランジュ関数 (目的関数 + ラグランジュ乗数 × 制約) のヘッセ近似は 2FTF + Bk、量 sk および yk は以下になります。

sk=xk+1xkyk=xL(xk+1,λk+1)xL(xk,λk+1)2FT(xk+1)F(xk+1)(xk+1xk)=(AET(xk+1)AET(xk))λEk+1+(AIT(xk+1)AIT(xk))λIk+1+2FT(xk+1)f(xk+1)2FT(xk)f(xk)2FT(xk+1)F(xk+1)(xk+1xk).(17)

この変更されたヘッセ近似では、同じ問題の fmincon 反復と、lsqnonlin または lsqcurvefit 反復との違いが考慮されます。内部的には、最小二乗ソルバーがこのヘッセ近似式 16式 17を、fmincon ソルバーのカスタム ヘッセ関数として使用します。

参考

| | | |

関連するトピック