2 次錐計画法アルゴリズム

2 次錐計画法の定義

2 次錐計画問題は、次の形式をとります。

$$\underset{x}{\min }{f}^{T}x$$

以下の制約に従います。

$$\begin{array}{c}\Vert A_{\text{sc}}(i)\cdot x-b_{\text{sc}}(i)\Vert \le d_{\text{sc}}^T(i)\cdot x-\gamma (i)\\ A\cdot x\le b\\ \text{Aeq}\cdot x=\text{beq}\\ \text{lb}\le x\le \text{ub}.\end{array}$$

f、x、b、beq、lb、ub はベクトル、A と Aeq は行列です。i ごとに、行列 A_sc(i)、ベクトルの b_sc(i) と d_sc(i)、およびスカラー γ(i) が、secondordercone を使用して作成された 2 次錐制約内に存在します。

つまり、この問題には、線形目的関数と線形制約だけでなく、$$\Vert A_{\text{sc}}(i)\cdot x-b_{\text{sc}}(i)\Vert \le d_{\text{sc}}^T(i)\cdot x-\gamma (i)$$ という形式の 2 次錐制約のセットも含まれています。

coneprog アルゴリズム

coneprog ソルバーは、Andersen、Roos、および Terlaky [1]で説明されているアルゴリズムを使用します。この手法は、linprog 内点法アルゴリズムと同様の内点法アルゴリズムです。

標準形式

アルゴリズムは、問題を "標準形式" で配置することから始まります。アルゴリズムが非負のスラック変数を追加するため、問題の形式は次のようになります。

$$\underset{x}{\min }{f}^{T}x$$

以下の制約に従います。

$$\begin{array}{l}A\cdot x=b\\ x\in K.\end{array}$$

ソルバーは、線形係数ベクトル f と線形制約行列 A のサイズをスラック変数を考慮するように拡大します。

領域 K は、"ローレンツ錐" 式 1と非負象限のクロス積です。各凸錐

$$\Vert A_{\text{sc}}(i)\cdot x-b_{\text{sc}}(i)\Vert \le d_{\text{sc}}^T(i)\cdot x-\gamma (i)$$

をローレンツ錐 式 1 に変換するには、変数 t₁、t₂、…、t_n+1 の列ベクトルを作成します。

$$\begin{array}{c}{t}_1=d^{T}x-\gamma \\ t_{2:(n+1)}=A_{\text{sc}}x-b_{\text{sc}}.\end{array}$$

ここで、各錐 i の変数 n の数は、A_sc(i) 内の行の数です。その定義により、変数ベクトル t は次の不等式を満たします。

$$\Vert t_{2:(n+1)}\Vert \le t_1.$$

(1)

式 1 は、(n+1) 個の変数でのローレンツ錐の定義です。変数 t が、凸領域 K 内の変数 x の代わりに問題内に表示されます。

内部的に、アルゴリズムは、錐制約の再編成で "回転したローレンツ錐" も使用しますが、このトピックでは取り上げません。詳細については、Andersen、Roos、および Terlaky [1]を参照してください。

スラック変数を追加すると、アルゴリズムは、必要に応じて、変数を無効にして、以下のような適切な制約を追加します。

• 境界が 1 つしかない変数の下限は 0 である。

• 境界が 2 つある変数の下限は 0 で、スラック変数を使用しているため、上限はない。

• 境界がない変数はスラック変数を制約付き変数として含むローレンツ錐に配置される。このスラック変数は、他の式、目的、または制約には含まれません。

双対問題

双対錐は以下のように表せます。

$$K_{*}=\left\{s:s^{T}x\ge 0\forall x\in K\right\}.$$

"双対" 問題は以下のように表せます。

$$\underset{y}{\max }{b}^{T}y$$

条件

$$A^{T}y+s=f$$

一部の

$$s\in K_{*}.$$

双対最適解は、双対制約を満たし、双対目的を最大化する点 (y,s) です。

同次自己双対形式

潜在的に実行不可能なまたは非有界の問題を処理するために、アルゴリズムは、2 つの変数の τ と κ を追加して、問題を同次 (0 に等しい) かつ自己双対として定式化します。

$$\begin{array}{c}Ax-b\tau =0\\ A^{T}y+s-f\tau =0\\ -f^{T}x+b^{T}y-\kappa =0\end{array}$$

(2)

および次の制約

$$(x;\tau )\in \overline{K},(s;\kappa )\in {\overline{K}}_{*}.$$

(3)

ここで、$$\overline{K}$$ は、(x;τ) の空間である非負の実数直線に接している錐 K です。同様に、$${\overline{K}}_{*}$$ は、(s;κ) の空間である非負の実数直線に接している錐 $$K_{*}$$ です。この定式化では、次の補助定理が、τ が実行可能解のスケーリングで、κ が実行不可能問題の指標であることを示しています。

補助定理 ([1] 補助定理 2.1)

(x, τ, y, s, κ) が 式 2 の実行可能解と 式 3 の制約であるとします。

• x^Ts + τκ = 0.

• τ > 0 の場合は、(x, y, s)/τ が標準形式の 2 次錐問題の主双対最適解になります。

• κ > 0 の場合は、以下の少なくとも 1 つの狭義の不等式が成り立ちます。

  b^Ty > 0

  f^Tx < 0.

  1 番目の不等式が成り立つ場合は、標準形式の主 2 次錐問題が実行不可能になります。2 番目の不等式が成り立つ場合は、標準形式の双対 2 次錐問題が実行不可能になります。

要約すると、実行可能問題の場合は、変数 τ によって、解が元の標準形式問題と同次自己双対問題の間でスケールされます。実行不可能問題の場合は、最後の反復 (x、y、s、τ、κ) が、元の標準形式問題が実行不可能であることを証明します。

開始点

反復の開始点は実行可能点です。

• 非負の変数の場合は x = 1、ローレンツ錐の 1 番目の変数の場合は 1、それ以外の場合は 0。

• y = 0.

• 錐の場合は s = (1,0,…,0)、非負の変数の場合は 1。

• τ = 1.

• κ = 1.

中央パス

アルゴリズムは、γ が 1 から 0 まで減少する以下の方程式のパラメーター化された解法である "中央パス" を辿ろうとします。

$$\begin{array}{c}Ax-b\tau =\gamma \left(A{x}_0-b{\tau }_0\right)\\ A^{T}y+s-c\tau =\gamma \left(A^T{y}_0+s_0-f{\tau }_0\right)\\ -f^{T}x+b^{T}y-\kappa =\gamma \left(-f^T{x}_0+b^T{y}_0-{\kappa }_0\right)\\ XSe=\gamma {\mu }_{0}e\\ \tau \kappa =\gamma {\mu }_0.\end{array}$$

(4)

• 0 の添字が付いた変数は、変数の開始点を示します。

• 変数の X と S は、それぞれ、x ベクトルと s ベクトルで構成された "矢じり" 行列です。ベクトル x = [x₁,x₂,…,x_n] の場合は、矢じり行列 X が次のように定義されます。

  $$X=\text{mat}(x)=\left[\begin{array}{cc}{x}_1& x_{2:n}^T\\ x_{2:n}& x_{1}I\end{array}\right].$$

  定義上、X は対称です。

• 変数 e は、x₁ ローレンツ錐座標に対応する錐座標が 1 のベクトルです。

• 変数 μ₀ は次のように定義されます。

  $${\mu }_0=\frac{x_0^T{s}_0+{\tau }_0{\kappa }_0}{k+1},$$

  ここで、k は、x₀ 内の非ゼロ要素の数です。

中央パスは、開始点で始まり、同次自己双対問題の最適解で終わります。

Andersen、Roos、および Terlaky [1] は、補助定理 3.1 で、相補性条件 x^Ts = 0 (ここで、x と s はローレンツ錐 L の積内にある) が次の条件と等しいことを示しています。

$$X_i{S}_i{e}_i=S_i{X}_i{e}_i=0$$

この条件はすべての錐 i に適用されます。ここで、X_i = mat(x_i), x_i はローレンツ錐 i, S_i = mat(s_i) に関連付けられた変数で、e_i は該当する次元の単位ベクトル [1,0,0,…,0] です。この議論では、中央パスが終了点で相補性条件を満たすことを示します。

探索方向

パラメーター γ が 1 から 0 まで減少する間に中央パスに近い点を取得するために、アルゴリズムはニュートン法を使用します。検出する変数には、(x, τ, y, s, κ) というラベルが付けられています。d_x は、x 変数の探索方向を表しているとします (その他も同様)。次に、ニュートン ステップが 式 4 から導出された次の連立一次方程式を解きます。

$$\begin{array}{c}A{d}_x-b{d}_{\tau }=(\gamma -1)\left(A{x}_0-b{\tau }_0\right)\\ A^T{d}_y+d_s-f{d}_{\tau }=(\gamma -1)\left(A^T{y}_0+s_0-f{\tau }_0\right)\\ -f^T{d}_x+b^T{d}_y-d_{\kappa }=(\gamma -1)\left(-f^T{x}_0+b^T{y}_0-\kappa \right)\\ X_0{d}_s+S_0{d}_x=-X_0{S}_{0}e+\gamma {\mu }_{0}e\\ {\tau }_0{d}_{\kappa }+{\kappa }_0{d}_{t}au=-{\tau }_0{\kappa }_0+\gamma {\mu }_0.\end{array}$$

アルゴリズムは、d の方向に 1 ステップ進むことによってその次の点を取得します。

$$\left[\begin{array}{c}{x}_1\\ {\tau }_1\\ y_1\\ \begin{array}{l}{s}_1\\ {\kappa }_1\end{array}\end{array}\right]=\left[\begin{array}{c}{x}_0\\ {\tau }_0\\ y_0\\ \begin{array}{l}{s}_0\\ {\kappa }_0\end{array}\end{array}\right]+\alpha \left[\begin{array}{c}{d}_x\\ d_{\tau }\\ d_y\\ \begin{array}{l}{d}_s\\ d_{\kappa }\end{array}\end{array}\right]$$

特定のステップ $$\alpha \in [0,1]$$ の場合。

数値安定性と加速収束のために、アルゴリズムは、Nesterov と Todd の提案 [8] に従ってステップをスケールします。また、アルゴリズムは、Mehrotra の予測子修正子法のバリアント [7] に従ってステップを修正します。(詳細については、Andersen、Roos、および Terlaky [1]を参照してください)。

ステップ ソルバーの種類

ここまでの説明は、値に 'augmented' が指定された LinearSolver オプションに関するものです。ソルバーの値を変更することで、問題のタイプに応じてステップ計算を変更できます。

• 'auto' (既定) — coneprog がステップ ソルバーを選択する。

  • スパースな問題の場合は 'prodchol' をステップ ソルバーとする。

  • それ以外の場合は 'augmented' をステップ ソルバーとする。

• 'normal' — ソルバーは、スパースな問題に適した 'augmented' ステップのバリアントを使用する。詳細については、Andersen、Roos、および Terlaky [1]を参照してください。

• 'schur' — ソルバーは、密な列がほとんどないスパースな問題を扱うために修正されたシューア補行列法を使用する。この手法は大きな錐にも適しています。詳細については、Andersen [2]を参照してください。

• 'prodchol' — ソルバーは、Goldfarb と Scheinberg ([4]と[5]) が説明する、密な列がほとんどないスパースな問題を扱うための手法を使用する。この手法は大きな錐にも適しています。

反復表示と停止条件

反復 k ごとに、アルゴリズムは、次の 3 つの相対収束尺度を計算します。

• 主実行不可能性

  $${\text{Infeas}}_{\text{Primal}}^k=\frac{\Vert A{x}_k-b{\tau }_k\Vert }{\max \left(1,\Vert A{x}_0-b{\tau }_0\Vert \right)}.$$

• 双対実行不可能性

  $${\text{Infeas}}_{\text{Dual}}^k=\frac{\Vert A^T{y}_k+s_k-f{\tau }_k\Vert }{\max \left(1,\Vert A^T{y}_0+s_0-f{\tau }_0\Vert \right)}.$$

• ギャップ実行不可能性

  $${\text{Infeas}}_{\text{Gap}}^k=\frac{\left|-f^T{x}_k+b^T{y}_k-{\kappa }_k\right|}{\max \left(1,\left|-f^T{x}_0+b^T{y}_0-{\kappa }_0\right|\right)}.$$

この 3 つの統計情報は、コマンド ラインで反復表示を指定することによって表示できます。

options = optimoptions('coneprog','Display','iter');

問題が実行可能で、ソルバーが収束すれば、この 3 つすべてが 0 に近づくはずです。実行可能問題では、変数 κ_k が 0 に近づき、変数 τ_k が正定数に近づきます。

停止条件の 1 つは、ギャップ実行不可能性に関係しています。停止条件は、次の最適性の尺度が最適性の許容誤差を下回ったときです。

$${\text{Optimality}}^k=\frac{\left|f^T{x}_k-b^T{y}_k\right|}{{\tau }_k+\left|b^T{y}_k\right|}=\frac{\left|f^T{x}_k/{\tau }_k-b^T{y}_k/{\tau }_k\right|}{1+\left|b^T{y}_k/{\tau }_k\right|}.$$

この統計尺度は、目的値の精度です。

ソルバーも停止し、問題が次の条件下で実行不可能であると宣言します。3 つの相対実行不可能性尺度は、c = ConstraintTolerance 未満で、次のようになります。

$${\tau }_k\le c\max \left(1,{\kappa }_k\right).$$

b^Ty_k > 0 の場合は、ソルバーが主問題が実行不可能であることを宣言します。f^Tx_k < 0 の場合は、ソルバーが双対問題が実行不可能であることを宣言します。

アルゴリズムは、以下の場合にも停止します。

$${\mu }_k\le c{\mu }_0$$

および以下となります。

$${\tau }_k\le c\max \left(1,{\kappa }_k\right).$$

この場合は、coneprog が、問題が数値的に不安定であると報告します (終了フラグ -10)。

残りの停止条件は、少なくとも 1 つの実行不可能性尺度が ConstraintTolerance を上回っており、計算されたステップ サイズが小さすぎる場合に発生します。この場合は、coneprog が、探索方向が小さくなりすぎて、これ以上進むことができないと報告します (終了フラグ -7)。

参照

[1] Andersen, E. D., C. Roos, and T. Terlaky. On implementing a primal-dual interior-point method for conic quadratic optimization. Math. Program., Ser. B 95, pp. 249–277 (2003). https://doi.org/10.1007/s10107-002-0349-3

[2] Andersen, K. D. A modified schur-complement method for handling dense columns in interior-point methods for linear programming. ACM Transactions on Mathematical Software (TOMS), 22(3):348–356, 1996.

[3] Ben-Tal, Aharon, and Arkadi Nemirovski. Convex Optimization in Engineering: Modeling, Analysis, Algorithms. (1998).

[4] Goldfarb, D. and K. Scheinberg. A product-form cholesky factorization method for handling dense columns in interior point methods for linear programming. Mathematical Programming, 99(1):1–34, 2004.

[5] Goldfarb, D. and K. Scheinberg. Product-form cholesky factorization in interior point methods for second-order cone programming. Mathematical Programming, 103(1):153–179, 2005.

[6] Luo, Zhi-Quan, Jos F. Sturm, and Shuzhong Zhang. Duality and Self-Duality for Conic Convex Programming. (1996).

[7] Mehrotra, Sanjay. “On the Implementation of a Primal-Dual Interior Point Method.” SIAM Journal on Optimization 2, no. 4 (November 1992): 575–601. https://doi.org/10.1137/0802028.

[8] Nesterov, Yu. E., and M. J. Todd. “Self-Scaled Barriers and Interior-Point Methods for Convex Programming.” Mathematics of Operations Research 22, no. 1 (February 1997): 1–42. https://doi.org/10.1287/moor.22.1.1.

参考

coneprog | secondordercone

関連するトピック

• 二次計画法と錐計画法