core 設計

linear-program は、変数、線形式、目的関数、制約、および 2 段階シンプレックス法のソルバーを src の 1 つのパッケージにまとめています。このページでは、ソルバーが実装する数学を述べ、停止に用いる判定を導出し、コードがそれらを中心にどう構成されているかを説明します。

設計目標

• 教科書の記法で線形計画を書けるようにすること。名前付きの変数、向きを持つ目的関数、等式制約と不等式制約です。
• 教科書のアルゴリズム、すなわち密なシンプレックス表上の 2 段階シンプレックス法で解き、各ステップ（Lp::pivot、Lp::simplex_iteration、Lp::phase_1、Lp::phase_2）を公開して、アルゴリズムを追跡し教えられるようにすること。
• luna-generic の trait によって係数型をジェネリックに保ち、シンプレックス表を linear-algebra の行列として保持すること。

数学的背景

線形計画と標準形

変数 $x \in \mathbb{R}^n$ 上の線形計画は、線形の等式・不等式のもとで線形の目的関数を最適化します。ソルバーが扱うのは標準形

$$\min_{x \in \mathbb{R}^n} \; c^\top x \quad \text{subject to} \quad A x = b, \quad x \ge 0, \qquad A \in \mathbb{R}^{m \times n},\; b \in \mathbb{R}^m_{\ge 0}.$$

このパッケージで書いた計画はすべて、最適解を変えずにこの形に変換できます。Lp::to_standard は次の書き換えを行います。

$$\begin{aligned} \max\; c^\top x &= -\min\; (-c)^\top x, \\ a^\top x = \beta,\ \beta < 0 \quad &\Longleftrightarrow \quad (-a)^\top x = -\beta, \\ a^\top x \le \beta,\ \beta \ge 0 \quad &\Longleftrightarrow \quad a^\top x + y = \beta,\ y \ge 0, \\ a^\top x \le \beta,\ \beta < 0 \quad &\Longleftrightarrow \quad (-a)^\top x - y = -\beta,\ y \ge 0, \\ a^\top x \ge \beta,\ \beta > 0 \quad &\Longleftrightarrow \quad a^\top x - y = \beta,\ y \ge 0, \\ a^\top x \ge \beta,\ \beta \le 0 \quad &\Longleftrightarrow \quad (-a)^\top x + y = -\beta,\ y \ge 0. \end{aligned}$$

各同値性が成り立つのは、$y$ が不等式のスラックを測るからです。$a^\top x \le \beta$ に対して $y = \beta - a^\top x$ とおくと、これは不等式が成り立つときちょうど非負になります。等式に $-1$ を掛けても何も変わりませんが、不等式に $-1$ を掛けると向きが逆になります。これが $\beta < 0$ と $\beta \le 0$ の場合にスラック変数の符号が反転する理由です。どの場合も新しい右辺は非負になり、これは第 1 段階に必要です。

標準形は、元の計画の変数を含むすべての変数について $x \ge 0$ を仮定します。Variable に保存された境界は使われません。

基底と基底解

$A$ の階数を $m$ とします。基底とは、部分行列 $A_B$ が正則となる $m$ 個の列添字の集合 $B$ で、残りの添字が $N$ をなします。非基底変数をゼロとおくと基底解

$$x_N = 0, \qquad x_B = A_B^{-1} b,$$

が定まり、$x_B \ge 0$ のとき実行可能です。線形計画法の基本定理によれば、計画に最適解があれば最適な基底実行可能解（基底可能解）も存在するため、シンプレックス法は基底だけを訪れればよいのです。¹1 例えば Chvátal, Linear Programming（1983）第 3 章、または Bertsimas と Tsitsiklis, Introduction to Linear Optimization（1997）定理 2.7 を参照してください。

シンプレックス表と被約費用

ソルバーは計画を $(m + 1) \times (n + 1)$ のシンプレックス表として保持します。

$$T = \begin{bmatrix} c^\top & 0 \\ A & b \end{bmatrix}.$$

第 0 行が目的関数行、第 $1, \dots, m$ 行が制約で、最終列が右辺です。$T_{rs} \ne 0$ である位置 $(r, s)$ でのピボットは Gauss–Jordan の 1 ステップです。第 $r$ 行を $T_{rs}$ で割り、第 0 行を含む他のすべての行 $i$ から新しい第 $r$ 行の $T_{is}$ 倍を引きます。その後、第 $s$ 列は単位ベクトル $e_r$ になります。これがまさに Lp::pivot です。

一連のピボットによって基底 $B$ の列が単位ベクトルになったとします。各ピボットは正則行列を左から掛けることであり、第 0 行には制約行の倍数しか加えられていないので、表は次の形をしています。

$$T_B = \begin{bmatrix} c^\top - y^\top A & -y^\top b \\ A_B^{-1} A & A_B^{-1} b \end{bmatrix}$$

ここで $y \in \mathbb{R}^m$ はあるベクトルです。第 0 行の基底列はゼロなので $c_B^\top - y^\top A_B = 0$、すなわち

$$y^\top = c_B^\top A_B^{-1}.$$

$y$ を代入すると正準形のシンプレックス表が得られます。

$$T_B = \begin{bmatrix} \bar c^\top & -z_B \\ \bar A & \bar b \end{bmatrix}, \qquad \bar c^\top = c^\top - c_B^\top A_B^{-1} A, \quad z_B = c_B^\top A_B^{-1} b, \quad \bar A = A_B^{-1} A, \quad \bar b = A_B^{-1} b.$$

$\bar c_j$ が被約費用、$\bar b = x_B$ が基底解で、第 0 行の最後の要素は目的値の符号反転 $-z_B = -c^\top x$ です。これが、Lp::phase_2 が $-z$ を返し、Lp::two_stage が最小化の目的関数に対して再び符号を反転する理由です。

被約費用は、制約に沿って目的関数がどう変化するかを表します。$Ax = b$ を満たす任意の $x$ について基底変数を解くと $x_B = \bar b - A_B^{-1} A_N x_N$ です。したがって

$$\begin{aligned} c^\top x &= c_B^\top x_B + c_N^\top x_N \\ &= c_B^\top \bar b - c_B^\top A_B^{-1} A_N x_N + c_N^\top x_N \\ &= z_B + \bar c_N^\top x_N . \end{aligned}$$

最適性

$\bar c \ge 0$ ならば、その基底実行可能解は最適です。 実行可能な $x$ はすべて $x_N \ge 0$ を満たすので、上の恒等式から $c^\top x = z_B + \bar c_N^\top x_N \ge z_B$ であり、基底解は $z_B$ を達成します。Lp::simplex_iteration は第 0 行に負の要素がなくなると停止します。これがこの判定です（$\bar c_B = 0$ は構成により成り立ちます）。

入る列の選択

ある $\bar c_s < 0$ なら、$x_s$ を $0$ から増やすと目的関数は速さ $|\bar c_s|$ で減少します。ソルバーは最も負の被約費用を選びます。

$$s = \operatorname*{arg\,min}_{j} \; \bar c_j, \qquad \bar c_s < 0,$$

これが Dantzig の規則で、同値の場合は最も左の列を選びます。

比テストと非有界性

$x_s$ を $t \ge 0$ まで増やし、他の非基底変数はゼロのままにします。制約により $x_B(t) = \bar b - t\, \bar a_s$（$\bar a_s$ は $\bar A$ の第 $s$ 列）となり、目的関数は $z_B + t\, \bar c_s$ になります。

$\bar a_s \le 0$ ならば、計画は非有界です。 $x_B(t)$ の各成分は $t$ について非減少なので、すべての $t \ge 0$ で $x(t)$ は実行可能であり、一方 $z_B + t\,\bar c_s \to -\infty$ となります。この場合 Lp::simplex_iteration は Problem is unbounded で abort します。

そうでなければ、実行可能性のために $\bar a_{is} > 0$ のすべての行で $\bar b_i - t\, \bar a_{is} \ge 0$ が必要なので、最大のステップは

$$t^\ast = \min_{i \,:\, \bar a_{is} > 0} \frac{\bar b_i}{\bar a_{is}},$$

で、これが最小比テストです。最小値を達成する行 $r$ が基底から出ます。その基底変数は $t^\ast$ でゼロに達します。$(r, s)$ でピボットすると新しい基底 $B' = B \setminus \{B_r\} \cup \{s\}$ に移り、その基底解は $x(t^\ast) \ge 0$ なので実行可能性は保たれ、目的値は

$$z_{B'} = z_B + t^\ast\, \bar c_s \le z_B .$$

同値の場合、ソルバーは最も上の行を選びます。

停止性と退化

すべてのピボットで $t^\ast > 0$ なら（計画が非退化なら）、目的関数は狭義に減少し、基底は繰り返されず、高々 $\binom{n}{m}$ 回のピボットで停止します。出る行で $\bar b_r = 0$ のときは $t^\ast = 0$ となり、基底は変わっても点と目的値は変わりません。このような退化ピボットが続くと以前の基底に戻り、永遠に繰り返すことがあります。このパッケージが実装する、最小添字で同値を解消する Dantzig の規則は、小さな例で巡回することが知られています。²2 E. M. L. Beale, “Cycling in the dual simplex algorithm”, Naval Research Logistics Quarterly 2（1955）は 3 つの制約を持つ巡回の例を示しています。R. G. Bland, “New finite pivoting rules for the simplex method”, Mathematics of Operations Research 2（1977）は、候補の中から添字が最小の入る変数と出る変数を選べば巡回が防げることを証明しました。 このパッケージには巡回防止の規則がありません。実行は max_iterations 回のピボット（既定は 1000）で停止してメッセージを表示し、そのとき返される表は最適ではありません。

設計上の決定

人工変数による 2 段階法

問題。 シンプレックス法は基底実行可能解から出発しますが、標準形の計画はそれを伴っていません。

選択肢。 (a) ビッグ M 法: 大きな費用 $M$ を持つ人工変数を元の目的関数に加える。(b) 2 段階法: まず補助的な目的関数で実行可能な基底を見つけ、それから最適化する。(c) 実行可能な基底の提供を利用者に求める。

決定。 (b)。ビッグ M 法には、計画に対して十分大きく、かつ浮動小数点で他の係数を埋没させない程度に小さい $M$ の数値が必要です。2 段階法はそのような定数を必要とせず、第 1 段階で実行可能性も判定できます。

第 1 段階。 使える単位列を持たない各制約行 $i$ に、第 $i$ 行の係数が $1$ の人工変数 $a_i \ge 0$ を加えます。これらの行の集合を $R$ とします。第 1 段階は次を解きます。

$$\min\; w = \sum_{i \in R} a_i \quad \text{subject to} \quad A x + \sum_{i \in R} a_i e_i = b, \quad x \ge 0,\ a \ge 0 .$$

初期基底は人工列と既存の単位列からなり、基底解は $x_B = b \ge 0$ です。to_standard が $b$ を非負にしているので、これは実行可能です。

元の計画が実行可能であるのは、第 1 段階の最適値が $w^\ast = 0$ のときに限ります。 $x$ が $Ax = b$ について実行可能なら、$(x, a = 0)$ は第 1 段階で $w = 0$ の実行可能解であり、常に $w \ge 0$ なので $w^\ast = 0$ です。逆に、$w^\ast = 0$ と $a \ge 0$ から $a = 0$ となり、最適な $x$ は $Ax = b$ を満たします。

第 1 段階の表の正準形。 第 1 段階の目的関数の第 0 行は $[\,0 \mid \mathbf 1^\top \mid 0\,]$ から始まりますが、基底の人工列の費用が $1$ なので正準形ではありません。$A_B = I$ として公式 $y^\top = c_B^\top A_B^{-1}$ を適用すると $y = \mathbf 1_R$ となり、正準形の第 0 行は

$$\Bigl[\; -\sum_{i \in R} A_{i\cdot} \;\Big|\; 0 \;\Big|\; -\sum_{i \in R} b_i \;\Bigr],$$

で、第 0 行から各人工行を引くことで得られます。Lp::phase_1 は simplex_iteration を呼ぶ前にちょうどこの引き算を行います。

第 2 段階。 $w^\ast \ne 0$（後述の許容誤差を超える）なら Lp::phase_2 は abort します。計画は実行不可能です。そうでなければ人工列を削除し、元の費用 $c$ を第 0 行に書き込み、各基底列 $j$ についてその行の $c_j$ 倍を引くことで第 0 行を再び正準形にします。これは上と同じ消去です。そして第 1 段階で見つかった実行可能な基底から simplex_iteration を実行します。

単位列を初期基底として再利用する

$b \ge 0$ の <= 行のスラック変数の列は単位ベクトルに等しいので、人工変数なしで基底に入れて開始できます。第 1 段階の表を作る非公開ヘルパーは、ある列で 1 つの要素が $1$、制約行の他の要素がすべて $0$ であるとき単位列と認識し、残りの行にだけ人工変数を加えます。したがって、すべての制約が非負の右辺を持つ <= である計画では人工変数がまったく不要で、第 1 段階は事実上省略されます。その表は元の目的関数を保持するので、第 1 段階ですでにそれを最適化します。

密なシンプレックス表

シンプレックス表は密な @mutable.Matrix です。各ピボットには $(m + 1)(n + 1)$ 回の乗算と減算、すなわち $O(mn)$ かかります。基底を分解して持つ改訂シンプレックス法なら疎な計画で 1 反復あたりのコストは小さくなりますが、密な表は上の導出のすべての量を 1 つの行列の要素として示し、パッケージの教育的な目標にかない、コードも短く保てます。

エラーを返す代わりに abort する

実行不可能・非有界な計画、ゼロピボット、未知の目的関数の向き、未知の関係文字列は、abort でプログラムを終了させます。このパッケージは、構造化エラーを持つ Result を返すという Luna-Flow の規約より前に作られました。呼び出し側はこれらの結果から回復できません。下の境界を参照してください。

数値の許容誤差

シンプレックスの反復は要素をゼロと厳密に比較します。丸めによって生じた $-10^{-17}$ の被約費用は負とみなされ、もう 1 回のピボットを引き起こします。第 2 段階は trait ApproximatelyZero を 2 か所で使います。$w^\ast$ がゼロかどうかの判定と、単位列の認識です。Double のしきい値は $|x| < 10^{-15}$ で、絶対的な境界です。

ピボットにおける丸め誤差は、関係する要素の大きさに対して相対的です。浮動小数点での 1 回の更新 $t_{ij} - t_{is}\, t_{rj}$ は次を満たします。

$$\operatorname{fl}\bigl(t_{ij} - t_{is}\, t_{rj}\bigr) = \bigl(t_{ij} - t_{is}\, t_{rj}(1 + \delta_1)\bigr)(1 + \delta_2), \qquad |\delta_1|, |\delta_2| \le u = 2^{-53} \approx 1.11 \times 10^{-16},$$

したがって、その絶対誤差はおよそ $u\,(|t_{ij}| + 2|t_{is}\, t_{rj}|)$ で抑えられます。大きさがおよそ $M$ の要素に $k$ 回ピボットすると、$w^\ast$ の誤差は $k\,u\,M$ 程度になります。よってしきい値 $10^{-15} \approx 9u$ は、要素が $1$ 付近でピボット回数の少ない、スケールの整った計画に適しています。テストでは $w^\ast = -4.44 \times 10^{-16} = -4u$ のような残差が見られます。係数が数千に及ぶと残差がしきい値を超え、実行可能な計画が実行不可能と報告されることがあります。$b$ の大きさでスケールした相対許容誤差ならこれを避けられますが、実装されていません。

正しさと不変条件

ソルバーはピボットの間で次の不変条件を保ちます。

1. 正準形。 制約行の各基底列は単位ベクトルで、第 0 行はすべての基底列でゼロです。Lp::pivot は構成によりこれを保ち、Lp::phase_1 は人工基底についてこれを確立し、Lp::phase_2 は第 0 行を置き換えた後で再確立します。
2. 主実行可能性。 右辺 $\bar b$ は非負です。to_standard が $b \ge 0$ とし、上で導いたように比テストがそれを保ちます。
3. 第 0 行の目的値。 第 0 行の最後の要素は $-z_B$ です。

これらの不変条件のもとで、停止判定は上で導いたとおり正しいものです。負の被約費用がなければ最適であり、負の被約費用に非正の列が伴えば非有界です。

実装には、特定の場合にこれらの保証を破る既知の欠陥があります。利用者が避けられるようにここに挙げます。コードは変更していません。

• 目的関数行の整列。 Lp は第 0 行を Obj_func::to_vector で作ります。これは保存された係数を変数名の順に取り、ゼロ係数を飛ばします。そのため、目的係数にゼロがある場合や変数が名前順に宣言されていない場合、第 0 行は列の順序と食い違い、ソルバーは並べ替えられた目的関数を最適化します。例えば $x_1 \ge 1$、$x_2 \ge 2$ のもとで $x_2$ を最小化すると、$2$ ではなく $1$ が返ります。
• 基底の認識。 第 2 段階は単位列を基底変数と認識し、各行で最初のそのような列を採ります。2 つの列が同じ単位ベクトルのとき、目的値は正しくても、報告される解が誤った変数を指すことがあります。
• 退化した人工変数。 第 1 段階の後に人工変数が値ゼロで基底に残ると、第 2 段階はその列をピボットで追い出さずに削除し、その行は第 2 段階の間基底列を持ちません。
• 人工変数の添字。 人工変数を持たない行について、Lp::phase_1 は人工列として 0 を受け取り、これを第 0 列と区別できません。その時点で第 0 列の第 1 段階目的係数がたまたま $1$ だと、引くべきでない行が第 0 行から引かれます。

採用しなかった代替案

• ビッグ M 法。 「人工変数による 2 段階法」で述べた数値的な理由で採用しませんでした。
• Bland の規則。 停止は保証されますが、非退化の計画では Dantzig の規則よりはるかに多くのピボットを要することがよくあります。このパッケージは Dantzig の規則を維持し、代わりに反復回数に上限を設けます。
• 改訂シンプレックス法と内点法。 大規模な疎計画ではよりスケールしますが、シンプレックス表を隠してしまいます。このパッケージは意図的に表を公開しています。
• 代入による変数の境界の適用。 境界 $l \le x \le u$ は $x = l + x'$ と追加の制約 $x' \le u - l$ によって標準形に帰着できます。このパッケージは境界を表示のためだけに保持し、そのような制約は利用者に任せます。

対象外

このパッケージは意図的に次のことを行いません。

• エラーを値として返すこと。実行不可能・非有界な計画、ゼロピボット、不正な引数はプログラムを abort させます。
• $x \ge 0$ 以外の変数の境界を適用すること。
• シンプレックス 1 回あたり 1000 反復という上限を超えて巡回を防ぐこと。
• 疎性、ウォームスタート、基底の分解を活用すること。
• 双対値や感度範囲を計算すること。ただし最終基底のベクトル $y = (A_B^{-1})^\top c_B$ は表の中に暗に含まれています。
• 整数計画、混合整数計画、非線形計画を解くこと。

Footnotes

1. 例えば Chvátal, Linear Programming（1983）第 3 章、または Bertsimas と Tsitsiklis, Introduction to Linear Optimization（1997）定理 2.7 を参照してください。 ↩

2. E. M. L. Beale, “Cycling in the dual simplex algorithm”, Naval Research Logistics Quarterly 2（1955）は 3 つの制約を持つ巡回の例を示しています。R. G. Bland, “New finite pivoting rules for the simplex method”, Mathematics of Operations Research 2（1977）は、候補の中から添字が最小の入る変数と出る変数を選べば巡回が防げることを証明しました。 ↩