ドキュメンテーション

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

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

線形計画法の定義

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

minxfTx

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

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

linprog 内点法アルゴリズム

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

  1. 解決の前処理。つまり、問題を単純化して標準的な形式に変換します。

  2. 初期点の生成。初期点の選択は、内点法アルゴリズムを効率的に解くために特に重要であり、このステップには時間がかかることがあります。

  3. KKT 方程式を解くための予測子修正子の反復。このステップには通常、最も時間がかかります。

解決の前処理

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

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

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

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

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

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

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

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

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

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

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

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

初期点の生成

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

  1. x0ones(n,1) に初期化します。ここで、n は目的関数ベクトル f の要素数です。

  2. すべての有界の成分の下限が 0 になるように変換します。成分 i に有限の上限 u(i) がある場合は、x0(i) = u/2 になります。

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

  4. x0 が中心パスに近くなるようにするには、予測子修正子のステップを 1 回実行して、結果の位置およびスラック変数を変更して範囲内に十分に収まるようにします。中心パスの詳細は、Nocedal と Wright [7] (p. 397) を参照してください。

予測子修正子

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

minxfTx subject to {A¯x=b¯x+t=ux,t0.

  • ここでは、すべての変数に少なくとも 1 つの有限範囲があると仮定します。この仮定は、必要に応じて成分のシフトおよび符号反転を行うと、すべての x 成分の下限が 0 になることを意味します。

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

    A¯x=(Aeq0AI)(x0s),

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

  • t は、上限を等式に変換するスラックのベクトルです。

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

  • y: 線形等式に関連付けられたラグランジュ乗数

  • v: 下限に関連付けられたラグランジュ乗数 (正値性制約)

  • w: 上限に関連付けられたラグランジュ乗数

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

L=fTxyT(A¯xb¯)vTxwT(uxt).

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

fA¯Tyv+w=0A¯x=b¯x+t=uvixi=0witi=0(x,v,w,t)0.

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

(0A¯T0IIA¯0000I0I00V00X000W0T)(ΔxΔyΔtΔvΔw)=(fA¯Tyv+wA¯xb¯uxtVXWT)=(rdrprubrvxrwt),(1)

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

  • rd: 双対残差

  • rp: 主残差

  • rub: 上限残差

  • rvx: 下限相補性残差

  • rwt: 上限相補性残差

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

Primal infeasibility=rp1+rub1Dual infeasibility=rd.

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

(DA¯TA¯0)(ΔxΔy)=(Rrp),(2)

ここで、

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

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

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

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

式 2の 2 行目は

A¯Δx=rp

また、1 行目は

DΔx+A¯TΔy=R.

次を代入し

Δx=D1A¯TΔy+D1R

次の結果が得られます。

A¯D1A¯TΔy=A¯D1Rrp.(3)

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

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

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

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

停止条件

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

ρ=max(1,A¯,f,b¯).

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

rp1+rub1ρTolConrdρTolFunrcTolFun,

ここで、

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

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

レガシ内点線形計画法

はじめに

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

メイン アルゴリズム

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

minxfTx such that {Ax=b0xu.(4)

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

minxfTx such that {Ax=bx+s=ux0, s0.(5)

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

maxbTyuTw  such that  {ATyw+z=fz0, w0,(6)

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

F(x,y,z,s,w)=(Axbx+suATyw+zfxizisiwi)=0,                 x0, z0, s0, w0,(7)

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

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

xTz + sTw

"双対性ギャップ" で (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 です。まず、いわゆる "予測" 方向を計算します。

Δvp=(FT(v))1F(v),

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

Δvc=(FT(v))1F(v+Δvp)μe^,

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

v+=v+α(Δvp+Δvc),

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

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

[x+;z+;s+;w+] > 0

上の予測子/修正子の方向の解法の中で、アルゴリズムは A·AT のコレスキー因子を変更して、(スパースの) 直接因子分解を計算します。A が密な列をもつ場合は、代わりに Sherman-Morrison 式を使います。その解が適切でない (残差が大きすぎる) 場合、拡大系形式の LDL 分解の実行によって解を検出します。(MATLAB® ldl 関数リファレンス ページの 例 4 — D の構造体 (MATLAB) を参照)。

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

max(rbmax(1,b),rfmax(1,f),rumax(1,u),|fTxbTy+uTw|max(1,|fTx|,|bTyuTw|))tol,

ここで、

rb=Axbrf=ATyw+zfru={x}+su

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

fTxbTy+uTw

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

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

前処理

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

minxfTx such that {Ax=b0xu.

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

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

x+={xif x>00if x0.

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

Primal infeasibility=((lbx)+)2+((xub)+)2+((Axb)+)2+|Aeqxbeq|2.

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

maxbTyuTw  such that  {ATyw+z=fz0, w0.

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

Dual infeasibility=ATy+zwf2.

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

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

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

 第 1 段階のアウトライン

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

 第 2 段階のアウトライン

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

基底変数と非基底変数

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

minxfTx such that {Ax=b,lbxub.

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

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

参照

[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] 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 http://www.maths.ed.ac.uk/~gondzio/software/correctors.ps.

[5] 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.

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

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