線形計画法のアルゴリズム

線形計画法の定義

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

$\min_{x} f^{T} x$

次の 1 つ以上の式が成り立つ必要があります。

A·x ≤ b
Aeq·x = beq
l ≤ x ≤ u.

`linprog` 内点法アルゴリズム

linprog 'interior-point' アルゴリズムは、interior-point-convex quadprog アルゴリズムとよく似ています。また、linprog 'interior-point-legacy' アルゴリズムと多くの機能を共有しています。これらのアルゴリズムに共通する概要は次のとおりです。

解決の前処理。つまり、問題を単純化して標準的な形式に変換します。
初期点の生成。初期点の選択は、内点法アルゴリズムを効率的に解くために特に重要であり、このステップには時間がかかることがあります。
KKT 方程式を解くための予測子修正子の反復。このステップには通常、最も時間がかかります。

解決の前処理

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

範囲の上限と下限が等しい変数があるかどうかをチェックする。見つかった場合はその実行可能性をチェックし、修正してからそれらの変数を削除する。
含まれている変数が 1 つだけの線形不等式制約があるかどうかチェックする。見つかった場合はその実行可能性をチェックした後、線形制約を範囲に変更する。
含まれている変数が 1 つだけの線形等式制約があるかどうかチェックする。見つかった場合はその実行可能性をチェックし、修正してからそれらの変数を削除する。
ゼロの行を含む線形制約行列があるかどうかチェックする。見つかった場合はその実行可能性をチェックした後、それらの行を削除する。
範囲制約と線形制約に矛盾がないか判断する。
目的関数内の線形項としてのみ含まれ、どの線形制約にも含まれない変数があるかどうかチェックする。見つかった場合はその実行可能性と有界性をチェックした後、それらの変数を修正して適切な範囲制約を設定する。
線形不等式制約があれば、スラック変数を追加して線形等式制約に変更する。

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

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

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

簡単にするため、解決の前処理のステップで問題が解決しない場合、アルゴリズムはすべての有限の下限をゼロにシフトします。

初期点の生成

初期点 x0 を設定するために、アルゴリズムは次の処理を行います。

x0 を ones(n,1) に初期化します。ここで、n は目的関数ベクトル f の要素数です。
すべての有界の成分の下限が 0 になるように変換します。成分 i に有限の上限 u(i) がある場合は、x0(i) = u/2 になります。
要素の範囲制約が 1 つだけの場合は、必要に応じて要素に変更を加え、厳密に範囲内に収まるようにします。
x0 が中央パスに近くなるようにするには、予測子修正子のステップを 1 回実行して、結果の位置およびスラック変数を変更して範囲内に十分に収まるようにします。中央パスの詳細については、Nocedal と Wright [8] (p. 397) を参照してください。

予測子修正子

fmincon の内点法アルゴリズムと同様に、interior-point アルゴリズムもカルーシュ・キューン・タッカー (KKT) 条件が満たされる点を見つけようとします。線形計画問題についてこれらの方程式を記述するために、前処理後の線形計画問題の標準的な形式を考えます。

$\min_{x} f^{T} x subject to {\begin{matrix} \bar{A} x = \bar{b} \\ x + t = u \\ x, t \geq 0. \end{matrix}$

ここでは、すべての変数に少なくとも 1 つの有限範囲があると仮定します。この仮定は、必要に応じて成分のシフトおよび符号反転を行うと、すべての x 成分の下限が 0 になることを意味します。
$\bar{A}$ は、線形不等式および線形等式の両方を含む拡張された線形行列です。 $\bar{b}$ は、対応する線形等式ベクトルです。 $\bar{A}$ には、不等式制約を等式制約にするスラック変数 s でベクトル x を拡張するための項も含まれています。
$\bar{A} x = (\begin{matrix} A_{e q} & 0 \\ A & I \end{matrix}) (\begin{matrix} x_{0} \\ s \end{matrix}),$
ここで、x₀ は元のベクトル x を意味します。
t は、上限を等式に変換するスラックのベクトルです。

このシステムのラグランジュ関数には、次のベクトルが含まれます。

y: 線形等式に関連付けられたラグランジュ乗数
v: 下限に関連付けられたラグランジュ乗数 (正値性制約)
w: 上限に関連付けられたラグランジュ乗数

ラグランジュ関数は次のようになります。

$L = f^{T} x - y^{T} (\bar{A} x - \bar{b}) - v^{T} x - w^{T} (u - x - t) .$

そのため、このシステムの KKT 条件は次のようになります。

$\begin{matrix} f - {\bar{A}}^{T} y - v + w = 0 \\ \bar{A} x = \bar{b} \\ x + t = u \\ v_{i} x_{i} = 0 \\ w_{i} t_{i} = 0 \\ (x, v, w, t) \geq 0. \end{matrix}$

linprog アルゴリズムは、返されるラグランジュ乗数に対して、ここでの説明とは異なる符号規則を使用しています。この説明では、ほとんどの文献と同じ符号を使用しています。詳細については、lambda を参照してください。

アルゴリズムは最初にニュートン・ラフソン式からステップを予測し、その後、修正子ステップを計算します。修正子は非線形相補性方程式 s_iz_i = 0 の残差を減らそうとします。ニュートン・ラフソンステップは次のようになります。

(\begin{matrix} 0 & - {\bar{A}}^{T} & 0 & - I & I \\ \bar{A} & 0 & 0 & 0 & 0 \\ - I & 0 & - I & 0 & 0 \\ V & 0 & 0 & X & 0 \\ 0 & 0 & W & 0 & T \end{matrix}) (\begin{matrix} Δ x \\ Δ y \\ Δ t \\ \begin{array}{l} Δ v \\ Δ w \end{array} \end{matrix}) = - (\begin{matrix} f - {\bar{A}}^{T} y - v + w \\ \bar{A} x - \bar{b} \\ u - x - t \\ \begin{array}{l} V X \\ W T \end{array} \end{matrix}) = - (\begin{matrix} r_{d} \\ r_{p} \\ r_{u b} \\ \begin{array}{l} r_{v x} \\ r_{w t} \end{array} \end{matrix}),

(1)

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

r_d: 双対残差
r_p: 主残差
r_ub: 上限残差
r_vx: 下限相補性残差
r_wt: 上限相補性残差

反復表示は以下の量をレポートします。

$\begin{array}{l} Primal infeasibility = {‖ r_{p} ‖}_{1} + {‖ r_{u b} ‖}_{1} \\ Dual infeasibility = {‖ r_{d} ‖}_{\infty} . \end{array}$

式 1を解くには、まず対称行列形式に変換します。

(\begin{matrix} - D & {\bar{A}}^{T} \\ \bar{A} & 0 \end{matrix}) (\begin{matrix} Δ x \\ Δ y \end{matrix}) = - (\begin{matrix} R \\ r_{p} \end{matrix}),

(2)

ここで、

$\begin{array}{l} D = X^{- 1} V + T^{- 1} W \\ R = - r_{d} - X^{- 1} r_{v x} + T^{- 1} r_{w t} + T^{- 1} W r_{u b} . \end{array}$

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

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

式 2 は対称ですが、項 –D があるため正定値ではありません。そのため、コレスキー分解を使用して解くことができません。さらにいくつかのステップを実行して、正定値であり、コレスキー分解によって効率的に解くことができる、別の方程式を導きます。

式 2の 2 行目は

$\bar{A} Δ x = - r_{p}$

また、1 行目は

$- D Δ x + {\bar{A}}^{T} Δ y = - R .$

次を代入し

$Δ x = D^{- 1} {\bar{A}}^{T} Δ y + D^{- 1} R$

次の結果が得られます。

\bar{A} D^{- 1} {\bar{A}}^{T} Δ y = - \bar{A} D^{- 1} R - r_{p} .

(3)

通常、ニュートンステップを求めるための最も効率的な方法は、コレスキー分解を使用して式 3を Δy について解くことです。Δy の乗数である行列は明らかに対称行列であり、退化がない場合には正定値になるため、コレスキー分解の使用が可能です。その後、ニュートンステップを求めるために、結果を元の式に代入して Δx、Δt、Δv、および Δw を求めます。ただし、 $\bar{A}$ が密な列をもつ場合は、代わりに式 2を解く方が効率的なこともあります。linprog 内点法アルゴリズムは、列の密度に基づいて解のアルゴリズムを選択します。

アルゴリズムについての詳細については、Mehrotra [7]を参照してください。

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

予測子修正子アルゴリズムは、二次の項を除き、完全な quadprog の 'interior-point-convex' バージョンとほぼ同じです。詳細については、フル予測子修正子を参照してください。

停止条件

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

$ρ = \max (1, ‖ \bar{A} ‖, ‖ f ‖, ‖ \bar{b} ‖) .$

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

$\begin{matrix} {‖ r_{p} ‖}_{1} + {‖ r_{u b} ‖}_{1} \leq ρ TolCon \\ {‖ r_{d} ‖}_{\infty} \leq ρ TolFun \\ r_{c} \leq TolFun, \end{matrix}$

ここで、

$r_{c} = \max_{i} (\min (| x_{i} v_{i} |, | x_{i} |, | v_{i} |), \min (| t_{i} w_{i} |, | t_{i} |, | w_{i} |)) .$

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

レガシ内点線形計画法

はじめに

レガシ内点法は LIPSOL ([52]) をベースにしています。この方法は、主双対内点法である Mehrotra 予測子修正子アルゴリズム ([47]) のバリアントです。

メインアルゴリズム

アルゴリズムは一連の前処理ステップ (前処理を参照) を適用することから始めます。処理後、問題は次の型をしています。

\min_{x} f^{T} x such that {\begin{matrix} A \cdot x = b \\ 0 \leq x \leq u . \end{matrix}

(4)

上限の制約は制約行列 A に暗黙的に含まれています。主スラック変数 s を加えると、式 4 は以下になります。

\min_{x} f^{T} x such that {\begin{matrix} A \cdot x = b \\ x + s = u \\ x \geq 0, s \geq 0. \end{matrix}

(5)

これは "主問題" と考えられます。x は主変数で s は主スラック変数からなります。"双対" 問題は以下のように表せます。

\max b^{T} y - u^{T} w such that {\begin{matrix} A^{T} \cdot y - w + z = f \\ z \geq 0, w \geq 0, \end{matrix}

(6)

ここで、y と w は双対変数で、 z は双対スラック変数です。この線形計画法に対する最適性、すなわち主式 5および双対式 6は以下になります。

\begin{array}{l} F (x, y, z, s, w) = (\begin{matrix} A \cdot x - b \\ x + s - u \\ A^{T} \cdot y - w + z - f \\ x_{i} z_{i} \\ s_{i} w_{i} \end{matrix}) = 0, \\ x \geq 0, z \geq 0, s \geq 0, w \geq 0, \end{array}

(7)

ここで、x_iz_i と s_iw_i は要素単位の乗算を表します。

二次方程式 x_iz_i = 0 および s_iw_i = 0 は線形計画法に対する "相補性" 条件と呼ばれ、他の (線形) 方程式は "実行可能" 条件と呼ばれます。次の式で計算される量は、

x^Tz + s^Tw

"双対性ギャップ" で (x,z,s,w) ≥ 0 のとき F の相補性部分の残差を測定します。

このアルゴリズムは、"主双対アルゴリズム" です。すなわち、主計画と双対計画が同時に解かれます。これは式 7の中で線形二次システム F(x,y,z,s,w) = 0 に適用されるニュートン法に似た手法と考えられます。一方、反復の x、z、w および s は正に保たれ、そのため内点法の名前をもちます (反復は式 5 の中の不等式制約により表現される厳密な意味の内点領域に存在します)。

アルゴリズムは、Mehrotra が提唱する予測子修正子アルゴリズムのバリアントです。反復 v = [x;y;z;s;w] を考えてみましょう。ここでは [x;z;s;w] > 0 です。まず、いわゆる "予測" 方向を計算します。

$Δ v_{p} = - {(F^{T} (v))}^{- 1} F (v),$

これはニュートン方向になります。そして "修正子" の方向は次のようになります。

$Δ v_{c} = - {(F^{T} (v))}^{- 1} F (v + Δ v_{p}) - μ \hat{e},$

ここで、μ > 0 は "センタリング" パラメーターと呼ばれ、注意深く選択しなければなりません。 $\hat{e}$ は F(v) の二次方程式に対応する 1 をもつ、0 と 1 を要素とするベクトルです。すなわち、すべてが二次である相補性条件のみに摂動が適用され、すべて線形である実行可能条件には適用されません。2 つの方向はステップ長パラメーター α > 0 と結合され、新たな反復 v⁺ を得るために v が更新されます。

$v^{+} = v + α (Δ v_{p} + Δ v_{c}),$

ここで、ステップ長パラメーター α は

v⁺ = [x⁺;y⁺;z⁺;s⁺;w⁺]

が

[x⁺;z⁺;s⁺;w⁺] > 0.

上の予測子/修正子の方向の解法の中で、アルゴリズムは A·A^T のコレスキー因子を変更して、(スパースの) 直接因子分解を計算します。A が密な列をもつ場合は、代わりに Sherman-Morrison 式を使います。その解が適切でない (残差が大きすぎる) 場合、拡大系形式のステップ方程式の LDL 分解を実行して解を求めます。ldl を参照してください。

そしてアルゴリズムは反復が収束するまで反復します。メイン停止条件は標準的なものです。

$\max (\frac{‖ r_{b} ‖}{\max (1, ‖ b ‖)}, \frac{‖ r_{f} ‖}{\max (1, ‖ f ‖)}, \frac{‖ r_{u} ‖}{\max (1, ‖ u ‖)}, \frac{| f^{T} x - b^{T} y + u^{T} w |}{\max (1, | f^{T} x |, | b^{T} y - u^{T} w |)}) \leq t o l,$

ここで、

$\begin{matrix} r_{b} = A x - b \\ r_{f} = A^{T} y - w + z - f \\ r_{u} = {x} + s - u \end{matrix}$

は、それぞれ、主残差、双対残差、実行可能上限 ({x} は上限が有限であるそれらの x を意味する) です。

$f^{T} x - b^{T} y + u^{T} w$

は、主目的値と双対目的値の間の差で tol は許容誤差です。停止条件の和は式 7の最適化条件で相対的な誤差の総和を計測します。

主問題の実行不可能性の尺度は ||r_b|| で、双対問題の実行不可能性の尺度は ||r_f|| です。ここでノルムはユークリッドノルムです。

前処理

範囲の上限と下限が等しい変数があるかどうかをチェックする。見つかった場合はその実行可能性をチェックし、修正してからそれらの変数を削除する。
含まれている変数が 1 つだけの線形不等式制約があるかどうかチェックする。見つかった場合はその実行可能性をチェックした後、線形制約を範囲に変更する。
含まれている変数が 1 つだけの線形等式制約があるかどうかチェックする。見つかった場合はその実行可能性をチェックし、修正してからそれらの変数を削除する。
ゼロの行を含む線形制約行列があるかどうかチェックする。見つかった場合はその実行可能性をチェックした後、それらの行を削除する。
範囲制約と線形制約に矛盾がないか判断する。
目的関数内の線形項としてのみ含まれ、どの線形制約にも含まれない変数があるかどうかチェックする。見つかった場合はその実行可能性と有界性をチェックした後、それらの変数を修正して適切な範囲制約を設定する。
線形不等式制約があれば、スラック変数を追加して線形等式制約に変更する。

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

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

簡単にするため、アルゴリズムはすべての下限をゼロにシフトします。

これらの前処理ステップがアルゴリズムの反復部分をスピードアップすることができるとき、ラグランジュ乗数が必要であれば、アルゴリズムの処理中に計算される乗数は変換された問題に使用され、元のものには使用されないため、前処理ステップは未処理でなければなりません。そのため乗数が必要でない場合、この変換の逆が計算 "されず"、計算時間が節約されるかもしれません。

レガシ双対シンプレックス法アルゴリズム

上位レベルでは、linprog 'dual-simplex-legacy' アルゴリズムは実質的に、シンプレックス法アルゴリズムを "双対問題" に対して実行します。

このアルゴリズムは、まず前処理で説明されている前処理を行います。詳細については、Andersen と Andersen [1]、および Nocedal と Wright [8] の第 13 章を参照してください。この前処理によって、元の線形計画問題を式 4の形に簡約します。

$\min_{x} f^{T} x such that {\begin{matrix} A \cdot x = b \\ 0 \leq x \leq u . \end{matrix}$

A と b は、元の制約行列を変換したバージョンです。これが主問題になります。

主実行可能性は、⁺ に関して次のように定義できます。

$x^{+} = {\begin{array}{l} x & if x > 0 \\ 0 & if x \leq 0. \end{array}$

主問題の実行不可能性の尺度は次のようになります。

$Primal infeasibility = \sqrt{{({(lb - x)}^{+})}^{2} + {({(x - ub)}^{+})}^{2} + {({(A x - b)}^{+})}^{2} + {| Aeq x - beq |}^{2}} .$

式 6 で説明されているように、双対問題とは、次を解くベクトル y と w およびスラック変数ベクトル z を求めることです。

$\max b^{T} y - u^{T} w such that {\begin{matrix} A^{T} \cdot y - w + z = f \\ z \geq 0, w \geq 0. \end{matrix}$

双対問題の実行不可能性の尺度は次のようになります。

$Dual infeasibility = {‖ A^{T} y + z - w - f ‖}_{2} .$

双対問題の任意の有限解は主問題への解を与え、また主問題の任意の有限解が双対問題の解を与えることはよく知られています ([8] を参照)。さらに、主問題と双対問題のどちらかが非有界である場合、もう一方の問題は実行不可能です。そして、主問題または双対問題のどちらかが実行不可能な場合には、もう一方の問題は実行不可能であるか、非有界であるかのどちらかです。したがって、この 2 つの問題は、解が存在する場合には、有限解を求める点で等価です。主問題と双対問題は数学的に等価ではあるものの計算手順は異なるため、双対問題を解くことによって主問題を解く方が良い場合もあります。

退化を緩和するために (Nocedal と Wright [8] の 366 ページを参照)、双対シンプレックス法アルゴリズムは最初に目的関数に摂動を加えます。

双対シンプレックス法アルゴリズムの第 1 段階は、双対実行可能点を見つけることです。このアルゴリズムは、補助線形計画問題を解くことによってこれを求めます。

第 1 段階のアウトライン

第 1 段階では、このアルゴリズムは補助的な区分的線形計画問題を解くことにより、初期の基底実行可能解 (定義については基底変数と非基底変数を参照) を見つけます。補助問題の目的関数は "線形ペナルティ関数" $P = \sum_{j} P_{j} (x_{j})$ です。

ここで、P_j(x_j) は以下のように定義されます。

$P_{j} (x_{j}) = {\begin{array}{l} x_{j} - u_{j} & if x_{j} > u_{j} \\ 0 & if l_{j} \leq x_{j} \leq u_{j} \\ l_{j} - x_{j} & if l_{j} > x_{j} . \end{array}$

P(x) は、点 x が上限と下限の条件をどれだけ超えているかを測ります。補助問題は以下になります。

$\min_{x} \sum_{j} P_{j} subject to {\begin{matrix} A \cdot x \leq b \\ A e q \cdot x = b e q . \end{matrix}$

オリジナルの問題は、補助問題が最小値 0をもつ場合に限り、実行可能な基底点をもちます。

アルゴリズムは、必要に応じてスラック変数と人為変数を追加するというヒューリスティックな方法により補助問題に対する初期点を見つけます。すると、アルゴリズムは、補助問題を解くために、この初期点をシンプレックスアルゴリズムと共に使用します。この解は、メインアルゴリズムの段階 2 に対する初期点です。

第 2 段階では、ソルバーが entering 変数と leaving 変数を繰り返し選択します。アルゴリズムは、Forrest と Goldfarb [3] により提唱された双対最急エッジプライシング (dual steepest-edge pricing) という手法に従って leaving 変数を選択します。また、アルゴリズムは Koberstein [6] により提唱された Harris の比率判定法のバリエーションを使用して entering 変数を選択します。退化を緩和するため、アルゴリズムによって第 2 段階で追加の摂動が加えられる場合もあります。

第 2 段階のアウトライン

ソルバーは、主問題の実行不可能性を低減し、かつ双対問題の実行可能性を維持できるような方法で反復を行います。反復は、低減され摂動された問題について主問題と双対問題の両方が実行可能であるような解が見つかるまで続けられます。加えられた摂動はアルゴリズムによって取り除かれます。摂動のない元の問題において、摂動を加えた問題の解では双対問題が実行不可能な場合、ソルバーは主シンプレックス法または第 1 段階のアルゴリズムを使用して双対問題の実行可能性を復元します。最後に、ソルバーが前処理の手順を元に戻し、解を元の問題に返します。

HiGHS 双対シンプレックス法アルゴリズム

'dual-simplex-highs' アルゴリズムは HiGHS オープンソースソフトウェアに基づいています。上位レベルでは、linprog 'dual-simplex-highs' アルゴリズムは実質的に、双対実行可能性を維持しながら主実行可能性を実現するまで反復するシンプレックス法アルゴリズムを実行します。詳細については、主実行可能性と双対実行可能性を参照してください。アルゴリズムの詳細については、Huangfu と Hall [4]のセクション 2 を参照してください。

このアルゴリズムは、まず前処理で説明されている前処理を行います。背景については、Andersen と Andersen [1]、および Nocedal と Wright [8]の第 13 章を参照してください。通常、この前処理によって、元の線形計画問題がより小さく解きやすい問題に簡約されます。

退化を緩和するために (Nocedal と Wright [8] の 366 ページを参照)、双対シンプレックス法アルゴリズムは最初に目的関数の係数に摂動を加えます。

双対シンプレックス法アルゴリズムの第 1 段階の目的は、双対実行可能点を見つけることです。このアルゴリズムは、第 1 段階の問題に第 2 段階のアルゴリズム (後述) を当てはめることによってこれを求めます。この問題の (摂動を加えた) 目的関数では、変数と制約のコスト係数と範囲が次のように変更されます。上限 (下限) の片側制限付きの変数と制約の範囲は [0, 1] ([–1, 0]) に設定され、固定変数と固定方程式の範囲はゼロに設定され、自由変数と自由制約の範囲は [–1000, 1000] に設定されます。すべての範囲が有限であるため、主変数を適切な範囲に設定するといずれの基底も双対実行可能になります。(非基底) 変数の双対値が元の範囲に関して実行可能でない場合、双対目的関数値に対してマイナスに寄与します。双対目的関数の上限がゼロになることがわかります。第 1 段階の問題を最適性が得られるまで解いたときに目的値がゼロであれば、現在の点は元の範囲に関して双対実行可能です。負の目的値は、摂動を加えた目的関数のコスト係数では元の問題が双対実行不可能であることを意味し、摂動のない元の問題の双対実行不可能性を強く示唆します。このシナリオを評価するには、摂動のない目的関数の係数に戻し、引き続き現在の双対解から第 1 段階の問題を解くよう試みます。

第 1 段階の解がメインアルゴリズムの第 2 段階の初期点になります。

第 2 段階では、ソルバーが leaving 変数と entering 変数を繰り返し選択します。アルゴリズムは、Forrest と Goldfarb [3] により提唱された双対最急エッジプライシング (dual steepest-edge pricing) という手法に従って leaving 変数を選択します。また、アルゴリズムは Koberstein [6] により提唱された Harris の比率判定法のバリエーションを使用して entering 変数を選択します。丸めによる小さな双対実行不可能性は、コスト係数のシフトによってすべて除去されます。

第 2 段階では、第 1 段階で得られた双対実行可能点から始まる双対シンプレックス法アルゴリズムを使用して元の問題を解きます。アルゴリズムは、各反復で最適性条件 (主実行可能性) をテストし、現在の解が最適であれば停止します。現在の解が最適でない場合、アルゴリズムは定義について基底変数と非基底変数を使用して次のことを行います。

主値が実行不可能である基底変数から "leaving 変数" と呼ばれる変数を 1 つ選択します。
非基底変数から "entering 変数" と呼ばれる変数を 1 つ選択し、それを使用して基底の leaving 変数を置き換えます。
現在の主解および双対解と、現在の目的値を更新します。

いわゆる双対比テストでは、双対実行可能性を維持しながら leaving 変数の双対値を可能な限り大きくします。entering 変数は、この制限内で双対値がゼロになる非基底変数です。双対比のデータの計算には 1 つの方程式系の解が必要です。基底が変わった後は、2 つ目の方程式系を解いて主解を更新する必要があります。双対比テストに制限がなければ、摂動を加えた問題は双対非有界 (したがって主実行不可能) になります。

ソルバーは、双対実行可能性を維持しながら主実行不可能性を低減するように反復を行います。反復は、低減して摂動を加えた問題について、主問題と双対問題の両方が実行可能となる解が見つかるまで続けられます。その後、アルゴリズムは目的関数のコストの摂動を除去します。摂動のない元の問題において、摂動を加えた問題の解では双対問題が実行不可能な場合、ソルバーは主シンプレックス法の第 2 段階のアルゴリズムを使用して双対実行可能性を復元します。最後に、ソルバーが前処理の手順を元に戻し、解を元の問題に返します。

基底変数と非基底変数

このセクションでは線形計画問題に対して用語として"基底"、"非基底"、および"実行可能基底解" を定義します。この定義では、問題が次の標準形で与えられることを仮定します。

$\min_{x} f^{T} x such that {\begin{matrix} A \cdot x = b, \\ l b \leq x \leq u b . \end{matrix}$

(A および b は元の問題で不等式を定義する行列やベクトルではないことに注意してください)。A がランク m < n の m 行 n 列の行列で、その列が {a₁, a₂, ..., a_n} であると仮定します。 ${a_{i_{1}}, a_{i_{2}}, ..., a_{i_{m}}}$ がインデックス集合 B = {i₁, i₂, ..., i_m} として A の列空間に対する基底であり、N = {1, 2, ..., n}\B が B の補集合であると仮定します。部分行列 A_B は "基底" と呼ばれ、相補性部分行列 A_N は "非基底" と呼ばれます。"基底変数" のベクトルは x_B であり、"非基底変数" のベクトルは x_N です。段階 2 の各反復でアルゴリズムは現在の基底の 1 列を非基底の列で置き換え、それに応じて変数 x_B と x_N を更新します。

x が A·x = b の解であり、x_N の非基底変数すべてがそれらの下限または上限に等しいとき、x は "基底解" と呼ばれます。さらに x_B の基底変数がそれらの上限と下限を満たし、したがって x が実行可能な点である場合、x は "基底実行可能解" と呼ばれます。

主実行可能性と双対実行可能性

主実行不可能性の尺度は最大絶対制約違反として報告されます。これは、元の問題の変数で示すと次のようになります。

max(0, lb – x, x – ub, abs(Aeqx – beq), Ax – b)

(8)

双対実行可能性を定義するには、まず線形行列 A (基底変数と非基底変数の行列を使用) を 2 つのブロックに分割します。線形独立列を含む基底行列 B と非基底行列 N です。

$\begin{matrix} A x = [\begin{matrix} B & N \end{matrix}] [\begin{matrix} x_{B} \\ x_{N} \end{matrix}] \\ = B x_{B} + N x_{N} = b . \end{matrix}$

すると、方程式 A x = b の変数 x は、必然的に x_B と x_N に分割されます。

$\begin{matrix} A x = [\begin{matrix} B & N \end{matrix}] [\begin{matrix} x_{B} \\ x_{N} \end{matrix}] \\ = B x_{B} + N x_{N} = b . \end{matrix}$

B は線形独立列で構成されるため可逆であり、次のようになります。

$\begin{matrix} x_{B} = B^{- 1} (b - N x_{N}) \\ = B^{- 1} b - B^{- 1} N x_{N} . \end{matrix}$

目的関数 f^Tx を同じように記述できます。

$\begin{matrix} f^{T} x = f_{B}^{T} x_{B} + f_{N}^{T} x_{N} \\ = f_{B}^{T} (B^{- 1} b - B^{- 1} N x_{N}) + f_{N}^{T} x_{N} \\ = f_{B}^{T} B^{- 1} b + (f_{N}^{T} - f_{B}^{T} B^{- 1} N) x_{N} . \end{matrix}$

被約費用は、非基底変数の調整されたコスト係数です。

$c_{N}^{T} = f_{N}^{T} - f_{B}^{T} B^{- 1} N .$

すべての非基底変数の下限でこの被約費用値 c_j (c_N^T の j 番目のエントリ) が 0 以上であり、すべての上限で c_j が 0 以下であれば、現在の解は双対実行可能です。双対シンプレックス法アルゴリズムは、この双対実行可能性を維持しながら、シンプレックス法の反復で得られる解が主実行可能になるように試みます。

参照

[1] Andersen, E. D., and K. D. Andersen. Presolving in linear programming. Math. Programming 71, 1995, pp. 221–245.

[2] Applegate, D. L., R. E. Bixby, V. Chvátal and W. J. Cook, The Traveling Salesman Problem: A Computational Study, Princeton University Press, 2007.

[3] Forrest, J. J., and D. Goldfarb. Steepest-edge simplex algorithms for linear programming. Math. Programming 57, 1992, pp. 341–374.

[4] Huangfu, Q. and J. A. J. Hall. Parallelizing the dual revised simplex method. Math. Prog. Comp. 10, pp. 119–142, 2018. Available at https://link.springer.com/article/10.1007/s12532-017-0130-5.

[5] Gondzio, J. “Multiple centrality corrections in a primal dual method for linear programming.” Computational Optimization and Applications, Volume 6, Number 2, 1996, pp. 137–156. Available at https://www.maths.ed.ac.uk/~gondzio/software/correctors.ps.

[6] Koberstein, A. Progress in the dual simplex algorithm for solving large scale LP problems: techniques for a fast and stable implementation. Computational Optim. and Application 41, 2008, pp. 185–204.

[7] Mehrotra, S. “On the Implementation of a Primal-Dual Interior Point Method.” SIAM Journal on Optimization, Vol. 2, 1992, pp 575–601.

[8] Nocedal, J., and S. J. Wright. Numerical Optimization, Second Edition. Springer Series in Operations Research, Springer-Verlag, 2006.