ドキュメンテーション

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

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

二次計画法の定義

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

minx12xTHx+cTx

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

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

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

メモ

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

解決の前処理と後処理

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

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

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

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

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

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

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

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

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

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

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

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

初期点の生成

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

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

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

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

  4. 予測子修正子のフル ステップではなく、実現可能性の補正がわずかな予測子ステップ (予測子修正子を参照) を実行します。これにより、予測子修正子のフル ステップのオーバーヘッドを伴うことなく、初期点が "中央パス" に近づきます。中央パスの詳細は、Nocedal と Wright [7] (p. 397) を参照してください。

予測子修正子

スパース内点法凸アルゴリズムとフル内点法凸アルゴリズムは、主に予測子修正子フェーズで異なっています。アルゴリズムは類似していますが、一部の詳細部分が異なっています。アルゴリズムの基本的な説明については、Mehrotra [47]を参照してください。

スパース予測子修正子-  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).(9-1)

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

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

sizi = σrc

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

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

フル予測子修正子-  フル予測子修正子アルゴリズムは、線形制約に範囲を組み合わせないため、範囲に対応する別のスラック変数セットがあります。アルゴリズムは下限をゼロにシフトします。さらに、変数の範囲が 1 つのみの場合、アルゴリズムは上限の不等式の符号を反転して、それをゼロの下限にします。

A¯ は、線形不等式および線形等式の両方を含む拡張された線形行列です。b¯ は、対応する線形等式ベクトルです。A¯ には、不等式制約を等式制約にするスラック変数 s でベクトル x を拡張するための項も含まれています。

A¯x=(Aeq0AI)(x0s),

ここで、x0 は元のベクトル x を意味します。

KKT 条件は以下のとおりです。

Hx+cA¯Tyv+w=0A¯x=b¯x+t=uvixi=0, i=1,2,...,mwiti=0, i=1,2,...,nx,v,w,t0.(9-2)

解 x、式 9-2 へのスラック変数と双対変数を求めるために、アルゴリズムは基本的にニュートン・ラフソン ステップを考慮します。

(HA¯T0IIA¯0000I0I00V00X000W0T)(ΔxΔyΔtΔvΔw)=(Hx+cA¯Tyv+wA¯xb¯uxtVXWT)=(rdrprubrvxrwt),(9-3)

ここで、X、V、W、および T は、ベクトル x、v、w、および t にそれぞれ対応する対角行列です。方程式の最も右辺にある残差ベクトルは以下のとおりです。

  • rd: 双対残差

  • rp: 主残差

  • rub: 上限残差

  • rvx: 下限相補性残差

  • rwt: 上限相補性残差

アルゴリズムは、最初に対称行列形式に変換してから 式 9-3 を解きます。

(DA¯TA¯0)(ΔxΔy)=(Rrp),(9-4)

ここで、

D=H+X1V+T1WR=rdX1rvx+T1rwt+T1Wrub.

D および R の定義に含まれる逆行列は、行列が対角であるため、すべて簡単に計算できます。

式 9-3 から式 9-4 を導くには、式 9-4 の 2 行目が式 9-3 の 2 行目と同じであることに注目します。式 9-4の 1 行目を得るには、式 9-3の最後の 2 行を Δv および Δw について解き、その後 Δt について解きます。

式 9-4 を解くために、アルゴリズムは Altman および Gondzio [1] の基本要素に従います。アルゴリズムは、LDL 分解によって対称システムを解きます。Vanderbei や Carpenter [2] などの作成者によって指摘されているとおり、この分解はピボットを行わず数値的に安定しているため、高速になる場合があります。

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

フル quadprog 予測子修正子アルゴリズムは、linprog'interior-point' アルゴリズムとほとんど同様ですが、二次の項が含まれます。詳細は、予測子修正子を参照してください。

参照

[1] Altman, Anna and J. Gondzio. Regularized symmetric indefinite systems in interior point methods for linear and quadratic optimization. Optimization Methods and Software, 1999. Available for download here.

[2] Vanderbei, R. J. and T. J. Carpenter. Symmetric indefinite systems for interior point methods. Mathematical Programming 58, 1993. pp. 1–32. Available for download here.

停止条件

予測子修正子アルゴリズムは、実行可能 (許容誤差範囲内に収まるようにする制約を満たす) で相対ステップ サイズが小さい点に到達するまで反復します。具体的には、次のように定義します。

ρ=max(1,H,A¯,Aeq,c,b¯,beq)

これらのすべての条件が満たされると、アルゴリズムが停止します。

rp1+rub1ρTolConrdρTolFunrcTolFun,

ここで、

rc=maxi(min(|xivi|,|xi|,|vi|),min(|tiwi|,|ti|,|wi|)).

rc は実質的に、相補性の残差 xv および tw のサイズを測定します。これらの残差はそれぞれ、解でゼロ ベクトルとなるベクトルです。

実行不可能性の検出

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).

このメリット関数が大きくなりすぎると、quadprog は問題が実行不可能であると宣言し、終了フラグ -2 で停止します。

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

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

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

mins{q(s), sN}.(9-5)

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-6)

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

1Δ1s=0.

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

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

Hs2=g,(9-7)

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

s2THs2<0.(9-8)

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

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

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

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

  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-9)

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

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

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

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

ボックス制約

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

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

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

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

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

ここで、

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-12 には微分不可能な点があります。vi = 0 のとき微分不可能になります。厳密に実行可能性を維持することにより、このような点を避けることができます。すなわち、l < x < u の制約を適用します。

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

M^DsN=g^(9-13)

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

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

および以下となります。

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

ここで 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 です。

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