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

二次計画法の定義

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

$$\underset{x}{\min }\frac{1}{2}{x}^{T}Hx+c^{T}x$$

(1)

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

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

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

• 解決の前処理と後処理

• 初期点の生成

• 予測子修正子

• 停止条件

• 実行不可能性の検出

メモ

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

解決の前処理と後処理

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

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

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

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

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

• 範囲制約と線形制約に矛盾がないか判断する。

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

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

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

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

アルゴリズムが解決の前処理のステップで実行不可能または非有界の問題を検出しなかった場合、および解決の前処理が解を生成しなかった場合、そのアルゴリズムは次のステップに進みます。停止条件に到達すると、アルゴリズムは元の問題を復元し、解決の前処理による変更を元に戻します。この最終ステップが解決の後処理ステップです。

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

初期点の生成

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

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

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

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

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

予測子修正子

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

このアルゴリズムでは、はじめに線形不等式 Ax <= b を、A と b に -1 を乗算することにより Ax >= b という形式の不等式に変換します。これは解に影響を与えませんが、同じ形式の問題を一部の文献で見つけることができるようになります。

• スパース予測子修正子

• フル予測子修正子

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

$$\begin{array}{c}Hx+c-A_{eq}^{T}y-{\overline{A}}^{T}z=0\\ \overline{A}x-\overline{b}-s=0\\ A_{eq}x-b_{eq}=0\\ s_i{z}_i=0,i=1,2,\mathrm{...},m\\ s\ge 0\\ z\ge 0.\end{array}$$

ここで

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

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

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

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

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

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

• r_d は双対残差です。

  $$r_d=Hx+c-A_{eq}^{T}y-{\overline{A}}^{T}z.$$

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

  $$r_{eq}=A_{eq}x-b_{eq}.$$

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

  $$r_{ineq}=\overline{A}x-\overline{b}-s.$$

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

  r_sz = Sz.

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

• r_c は平均相補性です。

  $$r_c=\frac{s^{T}z}{m}.$$

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

$$\left(\begin{array}{cccc}H& 0& -A_{eq}^T& -{\overline{A}}^T\\ A_{eq}& 0& 0& 0\\ \overline{A}& -I& 0& 0\\ 0& Z& 0& S\end{array}\right)\left(\begin{array}{c}\Delta x\\ \Delta s\\ \Delta y\\ \Delta z\end{array}\right)=-\left(\begin{array}{c}{r}_d\\ r_{eq}\\ r_{ineq}\\ r_{sz}\end{array}\right).$$

(2)

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

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

s_iz_i = σr_c.

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

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

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

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

$$\overline{A}x=\left(\begin{array}{cc}{A}_{eq}& 0\\ A& I\end{array}\right)\left(\begin{array}{c}{x}_0\\ s\end{array}\right),$$

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

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

$$\begin{array}{c}Hx+c-{\overline{A}}^{T}y-v+w=0\\ \overline{A}x=\overline{b}\\ x+t=u\\ v_i{x}_i=0,i=1,2,\mathrm{...},m\\ w_i{t}_i=0,i=1,2,\mathrm{...},n\\ x,v,w,t\ge 0.\end{array}$$

(3)

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

$$\left(\begin{array}{ccccc}H& -{\overline{A}}^T& 0& -I& I\\ \overline{A}& 0& 0& 0& 0\\ -I& 0& -I& 0& 0\\ V& 0& 0& X& 0\\ 0& 0& W& 0& T\end{array}\right)\left(\begin{array}{c}\Delta x\\ \Delta y\\ \Delta t\\ \begin{array}{l}\Delta v\\ \Delta w\end{array}\end{array}\right)=-\left(\begin{array}{c}Hx+c-{\overline{A}}^{T}y-v+w\\ \overline{A}x-\overline{b}\\ u-x-t\\ \begin{array}{l}VX\\ WT\end{array}\end{array}\right)=-\left(\begin{array}{c}{r}_d\\ r_p\\ r_{ub}\\ \begin{array}{l}{r}_{vx}\\ r_{wt}\end{array}\end{array}\right),$$

(4)

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

• r_d: 双対残差

• r_p: 主残差

• r_ub: 上限残差

• r_vx: 下限相補性残差

• r_wt: 上限相補性残差

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

$$\left(\begin{array}{cc}-D& {\overline{A}}^T\\ \overline{A}& 0\end{array}\right)\left(\begin{array}{c}\Delta x\\ \Delta y\end{array}\right)=-\left(\begin{array}{c}R\\ r_p\end{array}\right),$$

(5)

ここで、

$$\begin{array}{l}D=H+X^{-1}V+T^{-1}W\\ R=-r_d-X^{-1}{r}_{vx}+T^{-1}{r}_{wt}+T^{-1}W{r}_{ub}.\end{array}$$

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

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

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

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

フル 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.

停止条件

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

$$\rho =\max \left(1,\Vert H\Vert ,\Vert \overline{A}\Vert ,\Vert A_{eq}\Vert ,\Vert c\Vert ,\Vert \overline{b}\Vert ,\Vert b_{eq}\Vert \right)$$

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

$$\begin{array}{c}{\Vert r_p\Vert }_1+{\Vert r_{ub}\Vert }_1\le \rho \text{TolCon}\\ {\Vert r_d\Vert }_{\infty }\le \rho \text{TolFun}\\ r_c\le \text{TolFun,}\end{array}$$

ここで、

$$r_c=\underset{i}{\max }\left(\min \left(\left|x_i{v}_i\right|,\left|x_i\right|,\left|v_i\right|\right),\min \left(\left|t_i{w}_i\right|,\left|t_i\right|,\left|w_i\right|\right)\right).$$

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

実行不可能性の検出

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

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

$$\begin{array}{c}\rho =\max \left(1,\Vert H\Vert ,\Vert \overline{A}\Vert ,\Vert A_{eq}\Vert ,\Vert c\Vert ,\Vert \overline{b}\Vert ,\Vert b_{eq}\Vert \right)\\ r_{\text{eq}}=A_{\text{eq}}x-b_{\text{eq}}\\ r_{\text{ineq}}=\overline{A}x-\overline{b}-s\\ r_{\text{d}}=Hx+c+A_{\text{eq}}^T{\lambda }_{\text{eq}}+{\overline{A}}^T{\overline{\lambda }}_{\text{ineq}}\\ g=\frac{1}{2}{x}^{T}Hx+c^{T}x-{\overline{b}}^T{\overline{\lambda }}_{\text{ineq}}-b_{\text{eq}}^T{\lambda }_{\text{eq}}.\end{array}$$

$$\overline{A}$$ と $$\overline{b}$$ の表記は、線形不等式の係数に、スパース アルゴリズムの範囲を表す項が補足されていることを意味します。同様に $${\overline{\lambda }}_{\text{ineq}}$$ の表記は、範囲制約を含む線形不等式制約のラグランジュ乗数を表します。これは予測子修正子では z と呼ばれ、$${\lambda }_{\text{eq}}$$ は y と呼ばれていました。

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

$$\frac{1}{\rho }\left(\max \left({\Vert r_{\text{eq}}\Vert }_{\infty },{\Vert r_{\text{ineq}}\Vert }_{\infty },{\Vert r_{\text{d}}\Vert }_{\infty }\right)+g\right).$$

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

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

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

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

$$\underset{s}{\min }\left\{q(s),s\in N\right\}.$$

(6)

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

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

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

$$\min \left\{\frac{1}{2}{s}^{T}Hs+s^{T}g\text{ such that }\Vert Ds\Vert \le \Delta \right\},$$

(7)

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

$$\frac{1}{\Delta }-\frac{1}{\Vert s\Vert }=0.$$

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

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

$$H\cdot s_2=-g,$$

(8)

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

$$s_2^T\cdot H\cdot s_2<0.$$

(9)

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

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

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

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

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 = C² です。ここで、C^–1HC^–1 は、条件数の良い行列または分類分けされた固有値をもつ行列です。

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

線形等式制約

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

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

$$\min \left\{f(x)\text{ such that }Ax=b\right\},$$

(10)

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

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

$$\left[\begin{array}{cc}C& {\tilde{A}}^T\\ \tilde{A}& 0\end{array}\right]\left[\begin{array}{c}s\\ t\end{array}\right]=\left[\begin{array}{c}r\\ 0\end{array}\right],$$

(11)

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

ボックス制約

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

$$\min \left\{f(x)\text{ such that }l\le x\le u\right\},$$

(12)

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

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

$${\left(D(x)\right)}^{-2}g=0,$$

(13)

ここで、

$$D(x)=\text{diag}\left({\left|v_k\right|}^{-1/2}\right),$$

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

• g_i < 0 および u_i < ∞ の場合、v_i = x_i – u_i とします。

• g_i ≥ 0 および l_i > –∞ の場合、v_i = x_i – l_i とします。

• g_i < 0 および u_i = ∞ の場合、v_i = –1 とします。

• g_i ≥ 0 および l_i = –∞ の場合、v_i = 1 とします。

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

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

$$\widehat{M}D{s}^N=-\widehat{g}$$

(14)

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

$$\widehat{g}=D^{-1}g=\text{diag}\left({\left|v\right|}^{1/2}\right)g,$$

(15)

および以下となります。

$$\widehat{M}=D^{-1}H{D}^{-1}+\text{diag}(g)J^v.$$

(16)

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

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

active-set quadprog アルゴリズム

解決の前処理のステップを完了した後、active-set アルゴリズムは 2 つのフェーズに進みます。

• フェーズ 1 — すべての制約について実行可能点を取得します。

• フェーズ 2 — 目的関数を反復的に小さくしながら、アクティブな制約のリストを維持するとともに、各反復における実行可能性を維持します。

active-set 手法 (射影法とも呼ばれる) は、[18]と[17]に記述されている Gill 等による方法に似ています。

解決の前処理のステップ

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

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

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

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

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

• 範囲制約と線形制約に矛盾がないか判断する。

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

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

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

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

アルゴリズムが解決の前処理のステップで実行不可能または非有界の問題を検出しなかった場合、および解決の前処理が解を生成しなかった場合、そのアルゴリズムは次のステップに進みます。停止条件に到達すると、アルゴリズムは元の問題を復元し、解決の前処理による変更を元に戻します。この最終ステップが解決の後処理ステップです。

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

フェーズ 1 のアルゴリズム

フェーズ 1 では、アルゴリズムは目的関数を考慮せずにすべての制約を満たす点 x を見つけようと試みます。quadprog は、提供された初期点 x0 が実行不可能な場合にのみ、フェーズ 1 のアルゴリズムを実行します。

はじめに、アルゴリズムは、X = Aeq\beq など、すべての等式制約について実行可能な点を見つけようと試みます。方程式 Aeq*x = beq に対する解 x が存在しない場合、アルゴリズムは停止します。解 X が存在する場合、範囲および線形不等式を満たすための次のステップに進みます。等式制約がない場合は初期点 X = x0 を設定します。

アルゴリズムは X から開始して M = max(A*X – b, X - ub, lb – X) を計算します。M <= options.ConstraintTolerance のとき、点 X は実行可能となり、フェーズ 1 のアルゴリズムは停止します。

M > options.ConstraintTolerance のとき、アルゴリズムは補助線形計画問題の非負のスラック変数 γ を導入します。

$$\underset{x,\gamma }{\min }\gamma$$

条件

$$\begin{array}{c}Ax-\gamma \le b\\ \text{Aeq}x=\text{beq}\\ \text{lb}-\gamma \le x\le \text{ub}+\gamma \\ \gamma \ge -\rho .\end{array}$$

ここで、ρ は、A および Aeq の最大要素の絶対値で乗算された ConstraintTolerance オプションです。γ = 0 であり等式および不等式を満たす点 x をアルゴリズムが見つけた場合、x はフェーズ 1 の実行可能点となります。γ = 0 で補助線形計画問題に対する解 x がない場合、フェーズ 1 の問題は実行不可能です。

補助線形計画問題を解くために、アルゴリズムは γ₀ = M + 1、x₀ = X と設定し、アクティブ セットを固定変数 (存在する場合) およびすべての等式制約として初期化します。アルゴリズムは、線形計画変数 p が現在の点 x₀ から x のオフセット、つまり x = x₀ + p となるように再定式化します。アルゴリズムは反復を使用して、この線形計画問題を解きます。フェーズ 2 においても、適切に修正されたヘッシアンと共にこの同じ反復を使用して二次計画問題を解きます。

フェーズ 2 のアルゴリズム

変数 d に関して、問題は次のとおりです。

$$\begin{array}{c}\underset{d\in {\Re }^n}{\min }q(d)=\frac{1}{2}{d}^{T}Hd+c^{T}d,\\ A_{i}d=b_i,i=1,\mathrm{...},m_e\\ A_{i}d\le b_i,i=m_e+1,\mathrm{...},m.\end{array}$$

(17)

ここで、A_i は m 行 n 列の行列 A における i 番目の行を指します。

フェーズ 2 の間、アクティブ セット $${\overline{A}}_k$$ は、解の点におけるアクティブな制約 (制約境界上にある制約) の推定値となります。

アルゴリズムは、各反復 k において $${\overline{A}}_k$$ を更新し、探索方向 d_k の基底を作成します。等式制約は常にアクティブ セット $${\overline{A}}_k$$ 内にあります。探索方向 d_k が計算され、任意の有効制約境界上に存在したまま目的関数が最小化されます。アルゴリズムは、d_k の部分可能空間を基底 Z_k から作成します。この基底の列はアクティブ セット $${\overline{A}}_k$$ の推定値に直交します (すなわち、$${\overline{A}}_k{Z}_k=0$$)。したがって、Z_k の列結合の線形和から作成される探索方向は、アクティブな制約条件の境界上にあることが保証されます。

アルゴリズムは、行列 $${\overline{A}}_k^T$$ を QR 分解した行列の終わりの n – l 列から行列 Z_k を作成します。ここで l はアクティブな制約の数であり、l < n です。つまり、Z_k は次のように与えられます。

$$Z_k=Q\left[:,l+1:n\right],$$

(18)

ここで、

$$Q^T{\overline{A}}_k^T=\left[\begin{array}{c}R\\ 0\end{array}\right].$$

アルゴリズムは、Z_k を見つけると、q(d) を最小にする新しい探索方向 d_k を求めます。ここで d_k はアクティブな制約条件のヌル空間です。つまり、d_k は Z_k の列の線形結合で、あるベクトル p に対して $${\widehat{d}}_k=Z_{k}p$$ となります。

d_k を置き換え、二次式を p の関数として見ると次のようになります。

$$q(p)=\frac{1}{2}{p}^T{Z}_k^{T}H{Z}_{k}p+c^T{Z}_{k}p.$$

(19)

この方程式を p で微分すると、次の結果が得られます。

$$\nabla q(p)=Z_k^{T}H{Z}_{k}p+Z_k^{T}c.$$

(20)

∇q(p) は、Z_k が定義する部分空間の射影勾配であるため、二次関数の射影勾配になります。$$Z_k^{T}H{Z}_k$$ は射影ヘッシアンと呼ばれます。射影ヘッセ行列 $$Z_k^{T}H{Z}_k$$ は半正定値であると仮定しているので、Z_k が定義する部分空間における関数 q(p) の最小値は ∇q(p) = 0 で生じます。これは次の線形方程式系の解の場合に成り立ちます。

$$Z_k^{T}H{Z}_{k}p=-Z_k^{T}c.$$

(21)

その後、アルゴリズムは次を作成します。

$$x_{k+1}=x_k+\alpha d_k,$$

ここで、

$$d_k=Z_{k}p.$$

目的関数の二次式の性質により、各反復におけるステップ長 α の選択肢は 2 つしかありません。d_k に沿った単位長さのステップは、$${\overline{A}}_k$$ のヌル空間に制限された関数の最小値に向かう厳密なステップになります。アルゴリズムが制約に違反せずにこのようなステップを取ることが可能な場合、そのステップが二次計画法 (式 18) の解です。そうでない場合、最近傍の制約に向かう d_k に沿ったステップが 1 単位より小さくなり、アルゴリズムは次の反復でアクティブ セットに新しい制約を含めます。任意方向 d_k の制約境界までの距離は次式で与えられます。

$$\alpha =\underset{i\in \left\{1,\mathrm{...},m\right\}}{\min }\left\{\frac{-\left(A_i{x}_k-b_i\right)}{A_i{d}_k}\right\},$$

これはアクティブ セット内にない制約に対して定義され、方向 d_k は制約の境界に向かっています。すなわち、$$A_i{d}_k>0,i=1,\mathrm{...},m$$ です。

アクティブ セットが n 個の独立した制約を含み、最小値の位置を含んでいない場合、アルゴリズムは、次の線形方程式の特異点にならないラグランジュ乗数 λ_k を計算します。

$${\overline{A}}_k^T{\lambda }_k=c+H{x}_k.$$

(22)

λ_k のすべての要素が非負の場合、x_k は二次計画問題式 1の最適解です。しかし、λ_k に負の要素があり、その要素が等式制約に対応していない場合、最小化は未完了となります。アルゴリズムは、最小の負の乗数に対応する要素を削除し、新しい反復を開始します｡

すべての非負のラグランジュ乗数でソルバーが完了するような場合、1 次の最適性の尺度が許容誤差を上回っているか、制約の許容誤差が満たされていません。このような場合、ソルバーは[1]に記述されている再開手続きに従って、より良い解を得ようと試みます。この手続きにおいて、ソルバーは現在のアクティブな制約セットを破棄し、制約をわずかに緩和してから新しいアクティブな制約セットを作成することにより、問題の求解におけるサイクル (同一の状態に繰り返し戻ること) を回避します。必要に応じて、ソルバーはこの再開手続きを何度も実行できます。

メモ

アルゴリズムを早期に停止するために ConstraintTolerance オプションおよび OptimalityTolerance オプションに大きな値を設定することは避けてください。一般に、ソルバーは、潜在的な停止点に到達するまでこれらの値を確認せずに反復します。そして、停止点に到達した場合にのみ、許容誤差が満たされているかどうかを確認します。

active-set アルゴリズムでは、問題が非有界になる場合の検出が困難になることもあります。この問題は、非有界 v の方向が二次の項 v'Hv = 0 である方向となった場合に発生する可能性があります。数値安定性を確保し、コレスキー分解を可能にするために、active-set アルゴリズムは、厳密に凸である小さな項を二次目的関数に追加します。この小さな項により、目的関数は -Inf から外れた状態で有界となります。この場合、active-set アルゴリズムは反復制限に到達し、解が非有界であることを報告しません。言い換えれば、このアルゴリズムは終了フラグ -3 ではなく 0 で停止します。

参照

[1] Gill, P. E., W. Murray, M. A. Saunders, and M. H. Wright. A practical anti-cycling procedure for linearly constrained optimization. Math. Programming 45 (1), August 1989, pp. 437–474.

ウォーム スタート

ウォーム スタート オブジェクトを開始点に使用して quadprog または lsqlin の 'active-set' アルゴリズムを実行すると、ソルバーはフェーズ 1 およびフェーズ 2 のステップの多くをスキップしようとします。ウォーム スタート オブジェクトには制約のアクティブ セットが含まれており、このセットは新しい問題に対して適切かほぼ適切である可能性があります。したがって、ソルバーは反復を回避し、制約をアクティブ セットに追加することができます。また、初期点は新しい問題の解に近い可能性があります。詳細については、optimwarmstart を参照してください。