ドキュメンテーション

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

二次計画法のアルゴリズム

二次計画法の定義

二次計画法は線形制約に従って、二次関数を最小化するベクトル x を見つける問題です。

minx12xTHx+cTx

条件 A·x ≤ b, Aeq·x = beq, l ≤ x ≤ u。

interior-point-convex quadprog アルゴリズム

interior-point-convex アルゴリズムは以下の手順を実行します。

解決の前処理と後処理

このアルゴリズムは、まず重複を取り除き制約を単純化することで問題を単純化しようと試みます。解決の前処理のステップで実行されるタスクは次のとおりです。

  • 範囲の上限と下限が等しい変数があるかどうかをチェックする。見つかった場合はその実行可能性をチェックし、修正してからそれらの変数を削除する。

  • 含まれている変数が 1 つだけの線形不等式制約があるかどうかチェックする。見つかった場合はその実行可能性をチェックし、線形制約を範囲に変更する。

  • 含まれている変数が 1 つだけの線形等式制約があるかどうかチェックする。見つかった場合はその実行可能性をチェックし、修正してからそれらの変数を削除する。

  • ゼロの行を含む線形制約行列があるかどうかチェックする。見つかった場合はその実行可能性をチェックし、それらの行を削除する。

  • 範囲制約と線形制約に矛盾がないかチェックする。

  • 目的関数内の線形項としてのみ含まれ、どの線形制約にも含まれない変数があるかどうかチェックする。見つかった場合はその実行可能性と有界性をチェックし、それらの変数を修正して適切な範囲制約を設定する。

  • 線形不等式制約があれば、スラック変数を追加して線形等式制約に変更する。

アルゴリズムは、実行不可能または非有界の問題を検出した場合、停止して適切な終了メッセージを表示します。

アルゴリズムが単一の実行可能点を得るような場合、これは解を表しています。

解決の前処理のステップで実行不可能または非有界の問題が検出されなかった場合、アルゴリズムは必要に応じて他のステップに進みます。最後に、アルゴリズムは元の問題を復元し、解決の前処理による変更を元に戻します。この最終ステップが解決の後処理ステップです。

詳細は、Gould および Toint [63] を参照してください。

初期点の生成

アルゴリズムの初期点 x0 は次のように決定されます。

  1. x0ones(n,1) に初期化します。n は H 内の行数です。

  2. 上限 ub と下限 lb の両方をもつ要素の場合、x0 の要素が厳密に範囲内にないときには、その要素は (ub + lb)/2 に設定されます。

  3. 要素の範囲制約が 1 つだけの場合は、必要に応じて要素に変更を加え、厳密に範囲内に収まるようにします。

予測子修正子

fmincon内点法アルゴリズムと同様に、interior-point-convex アルゴリズムも カルーシュ・キューン・タッカー (KKT) 条件が満たされる点を見つけようとします。二次計画法の定義で説明している二次計画問題の場合、これらの条件は次のようになります。

Hx+cAeqTyA¯Tz=0A¯xb¯s=0Aeqxbeq=0sizi=0, i=1,2,...,ms0z0.

ここで

  • A¯ は、線形不等式として書かれた範囲を含む拡張された線形不等式行列です。b¯ は対応する線形不等式ベクトルで、範囲を含みます。

  • s は不等式制約を等式に変換するスラックのベクトルです。s の長さは m ですが、これは線形不等式と範囲の数です。

  • z は s に対応するラグランジュ乗数のベクトルです。

  • y は等式制約に関連付けられたラグランジュ乗数のベクトルです。

アルゴリズムは最初にニュートン・ラフソン式からステップを予測し、その後、修正子ステップを計算します。修正子は非線形制約 sizi = 0 をより適切に適用しようとします。

予測子ステップに関する定義を以下に示します。

  • rd は双対残差です。

    rd=Hx+cAeqTyA¯Tz.

  • req は主等式制約残差です。

    req=Aeqxbeq.

  • rineq は主不等式制約残差で、範囲とスラックを含みます。

    rineq=A¯xb¯s.

  • rsz は相補性の残差です。

    rsz = Sz.

    S はスラック項の対角行列で、z はラグランジュ乗数の列行列です。

  • rc は平均相補性です。

    rc=sTzm.

ニュートン ステップでは、x、s、y、および z における変化が次の式で計算されます。

(H0AeqTA¯TAeq000A¯I000Z0S)(ΔxΔsΔyΔz)=(rdreqrineqrsz).

ただし、s および z に対する正値性制約のため、完全なニュートン ステップが実行不可能な場合があります。したがって、quadprog は正値性を維持するために必要であればステップを短縮します。

さらに、内部での "中央" の位置を維持するために、アルゴリズムは sizi = 0 を解く代わりに正のパラメーター σ をとり、次の式を解こうと試みます。

sizi = σrc.

quadprog はニュートン ステップの方程式内の rsz を rsz + ΔsΔz – σrc1 で置き換えます (1 は 1 のベクトルです)。quadprog は予測子ステップでの計算用にシステムを対称で数値的により安定したものにするために、ニュートン方程式を並べ替えます。

詳細については、Mehrotra [47] を参照してください。

複数の訂正

quadprog は訂正されたニュートン ステップを計算した後さらに計算を実行して、現在のステップを長くするとともに、後続のステップをより良いものにする準備を行うことができます。これらの複数の訂正計算により、性能とロバスト性を両方とも向上させられます。詳細については、Gondzio [62] を参照してください。

総相対誤差

quadprog は各反復で "メリット関数" φ を計算します。メリット関数は実行可能性の尺度であり、"合計相対誤差" とも呼ばれます。メリット関数が大きくなりすぎると quadprog は停止します。この場合、quadprog により問題が実行不可能であると宣言されます。

メリット関数は問題の KKT 条件に関連しています。予測子修正子を参照してください。次の定義を使用します。

ρ=max(1,H,A¯,Aeq,c,b¯,beq)req=Aeqxbeqrineq=A¯xb¯+srd=Hx+c+AeqTλeq+A¯Tλ¯ineqg=xTHx+fTxb¯Tλ¯ineqbeqTλeq.

A¯b¯ の表記は、線形不等式の係数に、範囲を表す項が補足されていることを意味します。同様に λ¯ineq の表記は、範囲制約を含む線形不等式制約のラグランジュ乗数を表します。これは予測子修正子では z と呼ばれ、λeq は y と呼ばれていました。

メリット関数 φ は次のようになります。

1ρ(max(req,rineq,rd)+g).

trust-region-reflective quadprog アルゴリズム

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

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

mins{q(s), sN}.(9-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Δ},(9-2)

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

1Δ1s=0.

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

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

Hs2=g,(9-3)

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

s2THs2<0.(9-4)

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

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

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

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

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

  4. Δ を調節します。

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

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

部分空間の信頼領域法が探索方向を決定するのに使われます。ただし、ステップを非線形最小化の場合と同じく (場合によっては) 1 つの反射ステップに制限する代わりに、区分的な Reflective ライン探索が各反復で行われます。ライン探索の詳細は [45] を参照してください。

前処理付き共役勾配法

線形方程式系 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 が使用されます (非線形最小化に対する信頼領域法を参照)。

線形等式制約

線形制約は、制約なしの最小化問題に対して記述された状況を複雑にします。しかし、前述の基本的な考え方は、明確かつ効率的な形で踏襲されます。Optimization Toolbox のソルバーの信頼領域法は、厳密に実行可能な反復処理を行います。

一般的な線形等号制約付き最小化問題は、次のように表現されます。

min{f(x)  such that  Ax=b},(9-5)

ここで A は m 行 n 列 (m ≤ n) の行列です。いくつかの Optimization Toolbox のソルバーは A を前処理し、AT の LU 分解をベースにした手法 (参考文献 [46]) を使って、厳密な意味で線形従属の部分を除去します。ここで A はランク m とします。

式 9-5 を解くためのメソッドは制約なしアプローチと 2 つの点で異なります。まず初期実行可能点 x0 は、スパースの最小二乗ステップを使って Ax0 = b となるように計算されます。第 2 に、近似退化ニュートン ステップ (または A のヌル空間における負の曲率方向) を計算するために、PCG アルゴリズムは RPCG (退化前処理付き共役勾配法) に置き換えられます ([46] を参照)。キーになる線形代数ステップは、次の型のシステムを解くことです。

[CA˜TA˜0][st]=[r0],(9-6)

ここで、A˜ は A を近似 (ランクを失わないという条件で A の小さな非ゼロがゼロに設定される) し、C は H へのスパース対称正定近似、すなわち C = H です。詳細は [46] を参照してください。

ボックス制約

ボックス制約問題は、次のように表現されます。

min{f(x)  such that  lxu},(9-7)

ここで、l は下限を表すベクトル、u は上限を表すベクトルです。l の要素のいくつか (またはすべて) は –∞ にすることができ、u の要素のいくつか (またはすべて) は ∞ にすることができます。この方法は厳密な意味で実行可能点の列を生成します。ロバストな収束挙動を達成しながら、可能領域を維持するために 2 つの手法が使われます。1 つ目の手法ではスケーリングされた変更ニュートン ステップが (2 次元の部分空間 S を定義するために) 制約のないニュートン ステップと置き換わります。2 つ目の手法では反射がステップ サイズを増加させるために使われます。

スケーリングされた変更ニュートン ステップは 式 9-7 式のキューン・タッカーの必要条件をチェックして生成されます。

(D(x))2g=0,(9-8)

ここで、

D(x)=diag(|vk|1/2),

ベクトル v(x) は 1 ≤ i ≤ n の範囲で次のように定義されます。

  • gi < 0 および ui < ∞ の場合、vi = xi – ui とします。

  • gi ≥ 0 および li > –∞ の場合、vi = xi – li とします。

  • gi < 0 および ui = ∞ の場合、vi = –1 とします。

  • gi ≥ 0 および li = –∞ の場合、vi = 1 とします。

非線形システム 式 9-8 には微分不可能な点があります。vi = 0 のとき微分不可能になります。厳密に実行可能性を維持することにより、このような点を避けることができます。すなわち、l < x < u の制約を適用します。

式 9-8 により与えられる線型方程式系に対する、スケーリングされた変更ニュートン ステップ sk は、

M^DsN=g^(9-9)

の第 k 反復における線形システムの解として定義されています。ここで

g^=D1g=diag(|v|1/2)g,(9-10)

および以下となります。

M^=D1HD1+diag(g)Jv.(9-11)

ここで Jv は |v| のヤコビアンの役割をします。対角行列 Jv の個々の対角要素は 0、–1、1 のいずれかに等しくなります。l と u のすべての要素が有限ならば、Jv = diag(sign(g)) になります。gi = 0 の点で、vi は微分可能にならない可能性があります。そのような点では Jiiv=0 が定義されます。このような微分不可能性は、vi がどの値を取るかは重要でないため、問題とはなりません。さらに、|vi| はこの点でもまだ不連続となりますが、関数 |vi|·gi は連続です。

次に、反射がステップ サイズを増加させるために使われます。(単一) 反射ステップは、次のように定義されます。範囲制約に交差するステップ p を与え、p に交差する最初の範囲制約を与えます。すなわち、これを i 番目の範囲制約 (i 番目の上限または i 番目の下限) であると仮定します。i 番目の要素を除いて、反射ステップは pR = p になります。ここで pRi = –pi です。

active-set quadprog アルゴリズム

問題 quadprog の処理を再度扱います。

minx12xTHx+cTx(9-12)

よって A·x ≤ b, Aeq·x = beq, and l ≤ x ≤ u となります。m は線形制約の合計数で、A と Aeq の行数の和です。

quadprog active-set アルゴリズムは有効制約法 (射影法とも呼ばれます) ですが、[18][17] に記述されている Gill 等の方法に似ています。この方法は修正されて、LP (線形計画法)、QP (二次計画法) のどちらにも使われています。

解を求めるステップは、2 つの部分からなっています。最初の段階は (解が存在する場合) 可能領域点を計算します。2 番目の段階は解に収束する可能領域点の列を生成します。

有効制約法の反復

この手法では有効制約法行列 Sk が、解を計算する点の有効な制約法 (すなわち制約境界上にあること) の推定値になるように維持されます。特に有効制約法のSk は Aeq の行と A の行の部分集合からなります。Sk は各反復 k で更新され、探索方向 dk の基底を作成するために使われます。等式制約は常に有効制約法 Sk 内にあります。探索方向 dk が計算され、任意の有効制約法の境界上に存在したまま目的関数が最小化されます。dk の部分可能空間は、基底 Zk から作成されます。この基底の列は有効制約法 Sk の推定値に直交します (すなわち、SkZk = 0 です)。これにより、Zk の列結合の線形和から形成される探索方向は、有効制約法の境界上にあることが保証されます。

行列 Zk は行列 SkT を QR 分解した行列の最後の m – l 列から形成されます。ここで l はアクティブな制約の数であり、l < m です。つまり、Zk は次のように与えられます。

Zk=Q[:,l+1:m],(9-13)

ここで、

QTSkT=[R0].

Zk が見つかると、dk における目的関数を最小にする探索方向 dk が求められます。ここで dk は有効制約法のヌル空間です。つまり、dk は Zk の列の線形結合で、あるベクトル p に対しdk = Zkp となります。

そして dk を置き換えて、p の関数として二次の目的関数を表示した場合、結果は以下になります。

q(p)=12pTZkTHZkp+cTZkp.(9-14)

p で微分すると、次の結果が得られます。

q(p)=ZkTHZkp+ZkTc.(9-15)

∇q(p) は、 Zk が定義する部分空間の射影勾配なので、二次関数の射影勾配になります。ZkTHZk は射影ヘッシアンと呼ばれます。ヘッセ行列 H は、正定であると仮定しているので、Zk が定義する部分空間での関数 q(p) の最小値は ∇q(p) = 0 で生じます。これは p が次の連立方程式系の場合に成り立ちます。

ZkTHZkp=ZkTc.(9-16)

次の手順は以下になります。

xk+1=xk+αdk,  where dk=Zkp.(9-17)

各反復で目的関数の二次式の性質により、ステップの長さ α の選択には 2 つの方法があります。dk に沿った単位長さのステップは、Sk のヌル空間に制限された関数の最小値に向かう厳密なステップになります。制約に違反せずに、このようなステップをとることができれば、これは QP の解 (式 9-12) になります。そうでない場合、最近傍の制約に向かう dk に沿ったステップが 1 単位より小さくなり、新しい制約は次の反復で有効制約法に含まれます。任意方向 dk の制約境界までの距離は次式で与えられます。

α=mini{1,...,m}{(Aixkbi)Aidk},(9-18)

これは有効制約法内にない制約に対して定義され、方向 dk は制約の境界に向かっています。すなわち、Aidk>0, i=1,...,m

ラグランジュ乗数 λk は次の線形方程式の特異点にならないように計算されます。

SkTλk=c.(9-19)

λk のすべての要素が正の場合、xk は QP の最適解です (式 9-12を参照)。しかし λk に負の要素があり、等式制約に対応していない場合、その要素は有効制約法から削除され、新しい反復が行われます。

初期化

アルゴリズムは実行可能な初期点を必要とします。初期点が実行可能でないならば、可能領域内の点は次の線形計画法を解くことにより求めることができます。

minγ, xnγ  such thatAix=bi,      i=1,...,me (the rows of Aeq)Aixγbi,  i=me+1,...,m  (the rows of A).(9-20)

表記 Ai は行列 A の i 番目の行です。式 9-20 式での実行可能点 (存在する場合) は等式制約を満たす値に x を設定して求めることができます。これは、等式制約の組から作成される過決定または、劣決定のどちらかの線形方程式を解くことにより得られます。この問題の解が存在する場合、スラック変数 γ はこの点で最大不等式制約になります。

LP 問題の前の QP アルゴリズムは、各反復で探索方向 d を最急降下方向に設定することにより変更できます。ここで、gk は目的関数の勾配 (線形目的関数の係数に等しい) です。

d=ZkZkTgk.(9-21)

LP 法を使って実行可能な点が求められたならば、メインの QP 部分に移行します。探索方向 dk は、線形連立方程式を解いて求まる d1 を使って初期化されます。

Hd1=gk,(9-22)

ここで gk は現在の反復 xk での目的関数の勾配です (すなわち Hxk + c です)。

この情報は役に立ちましたか?