Go と Python で解き比べ

TL;DR

Python ではじめる数理最適化（第2版） の 第２章 P 17 の問題を Go と Python で解き比べました。

問題設定

ある工場では製品 $p$ と $q$ を製造しています。製品 $p$ と $q$ を製造するには原料 $m$ と $n$ が必要で、次のことがわかっています。

• 製品 $p$ を $1$ kg 製造するには、原料 $m$ が $1$ kg、原料 $n$ が $2$ kg 必要である。
• 製品 $q$ を $1$ kg 製造するには、原料 $m$ が $3$ kg、原料 $n$ が $1$ kg 必要である。
• 原料 $m$ の在庫は $30$ kg、原料 $n$ の在庫は $40$ kg である。
• 単位量あたりの利得は、製品 $p$ は $1$ 万円、製品 $q$ は $2$ 万円である。

この条件で利得を最大にするためには、製品 $p$ と $q$ をそれぞれ何 kg 製造すればよいでしょうか。

PuLP で解いてみる

Python の数理最適化ライブラリである PuLP を使用してこの問題を解いてみましょう。

import pulp

problem = pulp.LpProblem("LP", pulp.LpMaximize)

x = pulp.LpVariable('x', cat='Continuous')
y = pulp.LpVariable('y', cat='Continuous')

problem += 1 * x + 3 * y <= 30
problem += 2 * x + 1 * y <= 40
problem += x >= 0
problem += y >= 0
problem += x + 2 * y

status = problem.solve()

print('Status:', pulp.LpStatus[status])
print('x= ', x.value(), 'y= ', y.value(), 'obj= ', problem.objective.value())

実行結果は以下のようになります。

Status: Optimal
x=  18.0 y=  4.0 obj=  26.0

製品 $p$ を 18.0 kg、製品 $q$ を 4.0 kg 製造することで、利得が最大値 26.0 万円になることがわかります。

Go で解いてみる

Go の数値計算ライブラリである gonum (https://gonum.org) を使用して、同じ問題を解いてみます。今回は optimize/convex/lp パッケージの Simplex Algorithm を利用します。

package main

import (
	"fmt"
	"log"

	"gonum.org/v1/gonum/mat"
	"gonum.org/v1/gonum/optimize/convex/lp"
)

func main() {
	// 最大化問題 minimize c^T x に変換するため、目的関数係数に -1 を掛ける
	c := []float64{-1, -2, 0, 0} // 製品p, 製品q, スラック変数s1, スラック変数s2 の係数

	// 制約条件 A * x = b の行列 A
	// 1*x + 3*y + 1*s1 + 0*s2 = 30
	// 2*x + 1*y + 0*s1 + 1*s2 = 40
	A := mat.NewDense(2, 4, []float64{
		1, 3, 1, 0,
		2, 1, 0, 1,
	})

	// 制約条件の右辺 b
	b := []float64{30, 40}

	// Simplex 関数を実行
	// 引数: 目的関数係数c, 制約行列A, 右辺b, 許容誤差tol, 初期基底(nilで自動探索)
	z, s, err := lp.Simplex(c, A, b, 0, nil)
	if err != nil {
		log.Fatal(err)
	}

	// 結果出力
	// zは最小化された目的関数の値なので、-1を掛けて元の最大化の値に戻す
	// s[0]がx、s[1]がyに対応
	fmt.Printf("x= %v y= %v obj= %v\n", s[0], s[1], z*-1)
}

実行結果は以下のようになります。

x= 18 y= 4 obj= 26

gonum の Simplex 関数は最小化問題を解くため、目的関数の値は符号が反転しますが、PuLP と同じ解 (x=18, y=4) を得ることができました。

Simplex Algorithm について

シンプレックス法について理解を深めるために、この問題を通してアルゴリズムのステップを手計算で追ってみましょう。

gonum のドキュメント によると、Simplex 関数の実装は Gilbert Strang の "Linear Algebra and Applications." (1976) の第8章に基づいています。また、UC Math 352 の講義動画 (lectures 11-13 of UC Math 352) も参考になります。この動画では、シンプレックス法の標準形、基底解、被約費用、投入変数と離脱変数の選択、非有界LP、縮退といった概念が解説されています。

問題文から立式する

まず、この問題を数理最適化問題として定式化します。製品 $p$ の生産量を $x$ (kg)、製品 $q$ の生産量を $y$ (kg) とすると、目的関数（利得の最大化）と制約条件は以下のようになります。

$$
\begin{align}
\displaylines{
\text{maximize} \quad 1 \cdot x + 2 \cdot y \newline
\text{s.t.} \quad x \geq 0 \newline
\quad y \geq 0 \newline
\quad 1 \cdot x + 3 \cdot y \leq 30 \quad \text{ (原料 m の制約)} \newline
\quad 2 \cdot x + 1 \cdot y \leq 40 \quad \text{ (原料 n の制約)} \newline
}
\end{align}
$$

• maximize : 目的関数
• s.t. : 制約条件 (subject to)

この問題をグラフで描画すると、実行可能領域は以下のようになります。シンプレックス法は、この実行可能領域の頂点をたどりながら最適な解を探します。

標準形にする

シンプレックス法を適用するためには、問題を標準形（minimize $c^T x$ s.t. $Ax = b, x \geq 0$）に変換する必要があります。

まず、目的関数を最大化から最小化に変更するため、目的関数の符号を反転させます。利得 $z = x + 2y$ の最大化は、$(-z) = -x - 2y$ の最小化と等価です。

次に、不等式制約を等式制約に変換するために、非負のスラック変数 (slack variables) $s_1, s_2$ を導入します。

• 原料 m の制約: $x + 3y \leq 30 \quad \rightarrow \quad x + 3y + s_1 = 30, \quad s_1 \geq 0$
• 原料 n の制約: $2x + y \leq 40 \quad \rightarrow \quad 2x + y + s_2 = 40, \quad s_2 \geq 0$

これにより、問題は以下の標準形に変換されます。

$$\displaylines{
\text{minimize} \quad -x - 2y + 0 \cdot s_1 + 0 \cdot s_2 \newline
\text{s.t.} \quad x + 3y + 1 \cdot s_1 + 0 \cdot s_2 = 30 \newline
\quad 2x + y + 0 \cdot s_1 + 1 \cdot s_2 = 40 \newline
\quad x, y, s_1, s_2 \geq 0
}
$$ベクトルと行列を用いて表すと、以下のようになります。ここで、$\boldsymbol{s} = \begin{pmatrix} x & y & s_1 & s_2 \end{pmatrix}^T$ です。

• 目的関数:
  $\text{minimize} \quad \boldsymbol{c}^T \boldsymbol{s}$
  ただし、$\boldsymbol{c}^T = \begin{pmatrix} -1 & -2 & 0 & 0 \end{pmatrix}$

• 制約条件:
  $\boldsymbol{A} \boldsymbol{s} = \boldsymbol{b}, \quad \boldsymbol{s} \geq \boldsymbol{0}$
  ただし、

$$
\boldsymbol{A} = \begin{pmatrix}
1 & 3 & 1 & 0 \newline
2 & 1 & 0 & 1
\end{pmatrix}, \quad \boldsymbol{b} = \begin{pmatrix}
30 \newline
40
\end{pmatrix}
$$制約行列 $\boldsymbol{A}$ の行数 $m=2$、列数 $n=4$ です。

シンプレックス法のステップ (手計算)

標準形に変換した問題を、シンプレックス法のステップに従って手計算で解いてみましょう。

シンプレックス法は、基底変数と非基底変数を定義し、基底変数が非負である実行可能基底解（実行可能領域の頂点）を一つずつ改善していく手法です。

初期化

最初は、スラック変数 $s_1, s_2$ を基底変数、元の変数 $x, y$ を非基底変数とします。非基底変数は $0$ に固定します。

$$\boldsymbol{s_N} =
\begin{pmatrix}
x \newline
y \newline
\end{pmatrix} = \begin{pmatrix}
0 \newline
0
\end{pmatrix}, \quad
\boldsymbol{s_B} =
\begin{pmatrix}
s_1 \newline
s_2 \newline
\end{pmatrix}
$$制約条件の式に代入すると、基底変数の値が求まります。

$$
\displaylines{
1 \cdot (0) + 3 \cdot (0) + s_1 = 30 \quad \implies s_1 = 30 \newline
2 \cdot (0) + 1 \cdot (0) + s_2 = 40 \quad \implies s_2 = 40
}
$$初期基底解は $(x, y, s_1, s_2) = (0, 0, 30, 40)$ となり、これはすべての変数が非負であるため実行可能です。このときの目的関数値（最小化目的）は $-1 \cdot 0 - 2 \cdot 0 + 0 \cdot 30 + 0 \cdot 40 = 0$ です。元の利得としては $0$ 万円です。

基底の更新

現在の基底解が最適かどうかを判定し、もし最適でなければ目的関数値を改善する方向へ基底変数を入れ替えます。このプロセスを繰り返します。

1回目の反復

現在の基底変数は $s_1, s_2$、非基底変数は $x, y$ です。

目的関数 $z = -x - 2y$ （最小化）を、非基底変数を用いて表します。現在の基底では、これはそのまま $z = -1 \cdot x - 2 \cdot y + 0 \cdot s_1 + 0 \cdot s_2$ です。

非基底変数 $x, y$ の係数（被約費用）は $-1, -2$ です。最小化問題の場合、被約費用が負の非基底変数を増加させると目的関数値を減らすことができます。最も負の値を持つのは $y$ の係数 $(-2)$ なので、$y$ を投入変数として選択します。

次に、$y$ を増加させたときに、基底変数である $s_1, s_2$ のどちらが先にゼロになるか、すなわち基底から離脱するかを決定します。制約条件の式を変形して、基底変数 $s_1, s_2$ を非基底変数 $x, y$ で表すことを考えます。

現在の基底行列 $\boldsymbol{B} = \begin{pmatrix} 1 & 0 \newline 0 & 1 \end{pmatrix}$、非基底行列 $\boldsymbol{N} = \begin{pmatrix} 1 & 3 \newline 2 & 1 \end{pmatrix}$、基底変数係数 $\boldsymbol{c_B} = \begin{pmatrix} 0 \newline 0 \end{pmatrix}$、非基底変数係数 $\boldsymbol{c_N} = \begin{pmatrix} -1 \newline -2 \end{pmatrix}$ です（最小化目的なので係数は負）。

制約条件 $\boldsymbol{A} \boldsymbol{s} = \boldsymbol{b}$ は $\boldsymbol{N}\boldsymbol{s_N} + \boldsymbol{B}\boldsymbol{s_B} = \boldsymbol{b}$ と書けます。よって

$$
\boldsymbol{s_B} = \boldsymbol{B}^{-1} \boldsymbol{b} - \boldsymbol{B}^{-1} \boldsymbol{N}\boldsymbol{s_N}
$$

現在の基底では $\boldsymbol{B} = \begin{pmatrix} 1 & 0 \newline 0 & 1 \end{pmatrix}$ なので $\boldsymbol{B}^{-1} = \begin{pmatrix} 1 & 0 \newline 0 & 1 \end{pmatrix}$ です。

$$
\displaylines{
\boldsymbol{s_B} =
\begin{pmatrix}
s_1 \newline
s_2
\end{pmatrix} = \newline
\begin{pmatrix}
1 & 0 \newline
0 & 1
\end{pmatrix}
\begin{pmatrix}
30 \newline
40
\end{pmatrix} -
\begin{pmatrix}
1 & 0 \newline
0 & 1
\end{pmatrix}
\begin{pmatrix}
1 & 3 \newline
2 & 1
\end{pmatrix}
\begin{pmatrix}
x \newline
y
\end{pmatrix} = \newline
\begin{pmatrix}
30 \newline
40
\end{pmatrix} -
\begin{pmatrix}
1 & 3 \newline
2 & 1
\end{pmatrix}
\begin{pmatrix}
x \newline
y
\end{pmatrix}
}
$$

$$\displaylines{
s_1 = 30 - x - 3y \newline
s_2 = 40 - 2x - y
}
$$非基底変数 $x$ を $0$ に固定し、$y = \theta$ とした場合、

$$
\displaylines{
s_1 = 30 - 3\theta \newline
s_2 = 40 - \theta
}
$$基底変数は非負でなければならないので、$s_1 \geq 0$ かつ $s_2 \geq 0$ です。

$$
\displaylines{
30 - 3\theta \geq 0 \implies 3\theta \leq 30 \implies \theta \leq 10 \newline
40 - \theta \geq 0 \implies \theta \leq 40
}
$$

両方の条件を満たす最大の $\theta$ は $\min(10, 40) = 10$ です。これは、非基底変数 $y$ を最大 10 まで増やすことができることを意味します。このとき、$s_1$ がゼロになります。したがって、離脱変数は $s_1$ となります。

ここで、基底変数と非基底変数を入れ替えます。新しい基底変数は $y, s_2$、非基底変数は $x, s_1$ です。

ピボット操作

制約条件の式から、新しい基底変数（$y, s_2$）を新しい非基底変数（$x, s_1$）で表すように変形します。離脱変数 $s_1$ に対応する式（1番目の制約式）を、投入変数 $y$ について解きます。

$$
\displaylines{
x + 3y + s_1 = 30 \newline
\implies 3y = 30 - x - s_1 \newline
\implies y = 10 - \frac{1}{3}x - \frac{1}{3}s_1
}
$$

この $y$ を2番目の制約式に代入します。

$$
\displaylines{
2x + (10 - \frac{1}{3}x - \frac{1}{3}s_1) + s_2 = 40 \newline
2x + 10 - \frac{1}{3}x - \frac{1}{3}s_1 + s_2 = 40 \newline
\frac{5}{3}x - \frac{1}{3}s_1 + s_2 = 30
}
$$

新しい基底変数と非基底変数の関係式が得られました。

$$\displaylines{
y = 10 - \frac{1}{3}x - \frac{1}{3}s_1 \newline
s_2 = 30 - \frac{5}{3}x + \frac{1}{3}s_1
}
$$行列で表すと

$$
\begin{pmatrix}
y \newline
s_2
\end{pmatrix} =
\begin{pmatrix}
10 \newline
30
\end{pmatrix} -
\begin{pmatrix}
1/3 & 1/3 \newline
5/3 & -1/3
\end{pmatrix}
\cdot
\begin{pmatrix}
x \newline
s_1
\end{pmatrix}
$$次に、目的関数 $z = -x - 2y$ を新しい非基底変数で表します。$y = 10 - \frac{1}{3}x - \frac{1}{3}s_1$ を代入します。

$$
\displaylines{
z = -x - 2(10 - \frac{1}{3}x - \frac{1}{3}s_1) \newline
z = -x - 20 + \frac{2}{3}x + \frac{2}{3}s_1 \newline
z = -20 - \frac{1}{3}x + \frac{2}{3}s_1
}
$$

新しい目的関数は $z + \frac{1}{3}x - \frac{2}{3}s_1 = -20$ とも書けます。新しい非基底変数 $x, s_1$ に対する被約費用は $-\frac{1}{3}, \frac{2}{3}$ です。

新しい非基底変数 $x = 0, s_1 = 0$ とおくと、新しい基底解は $(y, s_2) = (10, 30)$ となります。これは実行可能です（すべての変数が非負）。このときの目的関数値（最小化目的）は $z = -20$ です。元の利得としては $20$ 万円です。初期解の $0$ 万円から改善されました。

2回目の反復

現在の基底変数は $y, s_2$、非基底変数は $x, s_1$ です。
目的関数を非基底変数で表した式は $z = -20 - \frac{1}{3}x + \frac{2}{3}s_1$ です。
非基底変数 $x, s_1$ の係数（被約費用）は $-\frac{1}{3}, \frac{2}{3}$ です。負の被約費用を持つ非基底変数は $x$ のみです $(-\frac{1}{3}) < 0$ 。したがって、$x$ を投入変数として選択します。

次に、$x$ を増加させたときに、基底変数である $y, s_2$ のどちらが先にゼロになるかを決定します。

新しい基底変数と非基底変数の関係式

$$\displaylines{
y = 10 - \frac{1}{3}x - \frac{1}{3}s_1 \newline
s_2 = 30 - \frac{5}{3}x + \frac{1}{3}s_1
}
$$

非基底変数 $s_1$ を 0 に固定し、$x = \theta$ とした場合、

$$
\displaylines{
y = 10 - \frac{1}{3}\theta \newline
s_2 = 30 - \frac{5}{3}\theta
}
$$

基底変数は非負でなければならないので、$y \geq 0$ かつ $s_2 \geq 0$ です。

$$
\displaylines{
10 - \frac{1}{3}\theta \geq 0 \implies \frac{1}{3}\theta \leq 10 \implies \theta \leq 30 \newline
30 - \frac{5}{3}\theta \geq 0 \implies \frac{5}{3}\theta \leq 30 \implies \theta \leq 30 \cdot \frac{3}{5} = 18
}
$$

両方の条件を満たす最大の $\theta$ は $\min(30, 18) = 18$ です。これは、非基底変数 $x$ を最大 18 まで増やすことができることを意味します。このとき、$s_2$ がゼロになります。したがって、離脱変数は $s_2$ となります。

ここで、基底変数と非基底変数を入れ替えます。新しい基底変数は $y, x$、非基底変数は $s_1, s_2$ です。

ピボット操作

制約条件の式から、新しい基底変数（$y, x$）を新しい非基底変数（$s_1, s_2$）で表すように変形します。離脱変数 $s_2$ に対応する式（2番目の変形された制約式）を、投入変数 $x$ について解きます。

$$
\displaylines{
\frac{5}{3}x - \frac{1}{3}s_1 + s_2 = 30 \newline
\implies \frac{5}{3}x = 30 + \frac{1}{3}s_1 - s_2 \newline
\implies x = \frac{3}{5}(30 + \frac{1}{3}s_1 - s_2) \newline
= 18 + \frac{1}{5}s_1 - \frac{3}{5}s_2
}
$$

この $x$ を、1番目の変形された制約式 $y = 10 - \frac{1}{3}x - \frac{1}{3}s_1$ に代入します。

$$
\displaylines{
y = 10 - \frac{1}{3}(18 + \frac{1}{5}s_1 - \frac{3}{5}s_2) - \frac{1}{3}s_1 \newline
y = 10 - 6 - \frac{1}{15}s_1 + \frac{1}{5}s_2 - \frac{1}{3}s_1 \newline
y = 4 - \frac{1}{15}s_1 - \frac{5}{15}s_1 + \frac{1}{5}s_2 \newline
y = 4 - \frac{6}{15}s_1 + \frac{1}{5}s_2 \newline
y = 4 - \frac{2}{5}s_1 + \frac{1}{5}s_2 \newline
}
$$

新しい基底変数と非基底変数の関係式が得られました。

$$\displaylines{
y = 4 - \frac{2}{5}s_1 + \frac{1}{5}s_2 \newline
x = 18 + \frac{1}{5}s_1 - \frac{3}{5}s_2
}
$$行列で表すと

$$
\begin{pmatrix}
y \newline
x
\end{pmatrix} =
\begin{pmatrix}
4 \newline
18
\end{pmatrix} -
\begin{pmatrix}
2/5 & -1/5 \newline
-1/5 & 3/5
\end{pmatrix}
\cdot
\begin{pmatrix}
s_1 \newline
s_2
\end{pmatrix}
$$

次に、目的関数 $z = -20 - \frac{1}{3}x + \frac{2}{3}s_1$ を新しい非基底変数で表します。$x = 18 + \frac{1}{5}s_1 - \frac{3}{5}s_2$ を代入します。

$$
\displaylines{
z = -20 - \frac{1}{3}(18 + \frac{1}{5}s_1 - \frac{3}{5}s_2) + \frac{2}{3}s_1 \newline
z = -20 - 6 - \frac{1}{15}s_1 + \frac{1}{5}s_2 + \frac{2}{3}s_1 \newline
z = -26 - \frac{1}{15}s_1 + \frac{3}{15}s_2 + \frac{10}{15}s_1 \newline
z = -26 + \frac{9}{15}s_1 + \frac{3}{15}s_2 \newline
z = -26 + \frac{3}{5}s_1 + \frac{1}{5}s_2 \newline
}
$$

新しい目的関数は $z - \frac{3}{5}s_1 - \frac{1}{5}s_2 = -26$ とも書けます。新しい非基底変数 $s_1, s_2$ に対する被約費用は $\frac{3}{5}, \frac{1}{5}$ です。

最適性の判定

新しい非基底変数 $s_1, s_2$ に対する被約費用は $\frac{3}{5}, \frac{1}{5}$ であり、どちらも非負です。これは、これ以上目的関数値を改善できないことを意味します（最小化目的なので、負の被約費用がなければ最適）。したがって、現在の基底解が最適解となります。

現在の非基底変数 $s_1 = 0, s_2 = 0$ とおくと、基底変数の値は $(y, x) = (4, 18)$ となります。
最適解は $(x, y, s_1, s_2) = (18, 4, 0, 0)$ です。
このときの目的関数値（最小化目的）は $z = -26$ です。元の利得としては $z = -(-26) = 26$ 万円となります。

手計算の結果、製品 $p$ を $18$ kg、製品 $q$ を $4$ kg 製造するときに利得が最大となり、そのときの最大利得は $26$ 万円であることがわかりました。これは PuLP および gonum の結果と一致します。

gonum の Simplex 関数による実装

gonum の lp.Simplex 関数はこの標準形問題を解きますが、その内部処理は、先ほど手計算で追ったシンプレックス法のステップを自動で行っています。コードの内部で行われている主なステップを、手計算のステップと対比させながら見ていきましょう。

1. 入力の検証 (verifyInputs 関数)

コードの冒頭で、verifyInputs 関数が呼び出され、入力された c、A、b の整合性や、制約行列 A が問題を適切に定義しているか（例: 非ゼロの右辺 b に対応するゼロ行がないなど）を確認します。これは手計算の前に問題が解ける形になっているかを確認するステップに相当します。

verifyInputs 関数定義

func verifyInputs(initialBasic []int, c []float64, A mat.Matrix, b []float64) error {
	m, n := A.Dims()
	if m > n {
		panic("lp: more equality constraints than variables")
	}
	if len(c) != n {
		panic("lp: c vector incorrect length")
	}
	if len(b) != m {
		panic("lp: b vector incorrect length")
	}
	if len(c) != n {
		panic("lp: c vector incorrect length")
	}
	if len(initialBasic) != 0 && len(initialBasic) != m {
		panic("lp: initialBasic incorrect length")
	}

	// Do some sanity checks so that ab does not become singular during the
	// simplex solution. If the ZeroRow checks are removed then the code for
	// finding a set of linearly independent columns must be improved.

	// Check that if a row of A only has zero elements that corresponding
	// element in b is zero, otherwise the problem is infeasible.
	// Otherwise return ErrZeroRow.
	for i := 0; i < m; i++ {
		isZero := true
		for j := 0; j < n; j++ {
			if A.At(i, j) != 0 {
				isZero = false
				break
			}
		}
		if isZero && b[i] != 0 {
			// Infeasible
			return ErrInfeasible
		} else if isZero {
			return ErrZeroRow
		}
	}
	// Check that if a column only has zero elements that the respective C vector
	// is positive (otherwise unbounded). Otherwise return ErrZeroColumn.
	for j := 0; j < n; j++ {
		isZero := true
		for i := 0; i < m; i++ {
			if A.At(i, j) != 0 {
				isZero = false
				break
			}
		}
		if isZero && c[j] < 0 {
			return ErrUnbounded
		} else if isZero {
			return ErrZeroColumn
		}
	}
	return nil
}

2. 初期実行可能基底解の探索

シンプレックス法を開始するには、まず実行可能領域のどこか一つの頂点、すなわち初期実行可能基底解を見つける必要があります。

• lp.Simplex 関数は、initialBasic 引数で初期基底のインデックスを指定することもできますが、nil を指定すると自動的に探索します (findInitialBasic 関数)。この自動探索は、手計算で最初に $(x, y, s_1, s_2) = (0, 0, 30, 40)$ という初期基底解を見つけたステップに相当します。

findInitialBasic 関数定義

// findInitialBasic finds an initial basic solution, and returns the basic
// indices, ab, and xb.
func findInitialBasic(A mat.Matrix, b []float64) ([]int, *mat.Dense, []float64, error) {
	m, n := A.Dims()
	basicIdxs := findLinearlyIndependent(A)
	if len(basicIdxs) != m {
		return nil, nil, nil, ErrSingular
	}

	// It may be that this linearly independent basis is also a feasible set. If
	// so, the Phase I problem can be avoided.
	ab := mat.NewDense(m, len(basicIdxs), nil)
	extractColumns(ab, A, basicIdxs)
	xb := make([]float64, m)
	err := initializeFromBasic(xb, ab, b)
	if err == nil {
		return basicIdxs, ab, xb, nil
	}

	// This set was not feasible. Instead the "Phase I" problem must be solved
	// to find an initial feasible set of basis.
	//
	// Method: Construct an LP whose optimal solution is a feasible solution
	// to the original LP.
	// 1) Introduce an artificial variable x_{n+1}.
	// 2) Let x_j be the most negative element of x_b (largest constraint violation).
	// 3) Add the artificial variable to A with:
	//      a_{n+1} = b - \sum_{i in basicIdxs} a_i + a_j
	//    swap j with n+1 in the basicIdxs.
	// 4) Define a new LP:
	//   minimize  x_{n+1}
	//   subject to [A A_{n+1}][x_1 ... x_{n+1}] = b
	//          x, x_{n+1} >= 0
	// 5) Solve this LP. If x_{n+1} != 0, then the problem is infeasible, otherwise
	// the found basis can be used as an initial basis for phase II.
	//
	// The extra column in Step 3 is defined such that the vector of 1s is an
	// initial feasible solution.

	// Find the largest constraint violator.
	// Compute a_{n+1} = b - \sum{i in basicIdxs}a_i + a_j. j is in basicIDx, so
	// instead just subtract the basicIdx columns that are not minIDx.
	minIdx := floats.MinIdx(xb)
	aX1 := make([]float64, m)
	copy(aX1, b)
	col := make([]float64, m)
	for i, v := range basicIdxs {
		if i == minIdx {
			continue
		}
		mat.Col(col, v, A)
		floats.Sub(aX1, col)
	}

	// Construct the new LP.
	// aNew = [A, a_{n+1}]
	// bNew = b
	// cNew = 1 for x_{n+1}
	aNew := mat.NewDense(m, n+1, nil)
	aNew.Copy(A)
	aNew.SetCol(n, aX1)
	basicIdxs[minIdx] = n // swap minIdx with n in the basic set.
	c := make([]float64, n+1)
	c[n] = 1

	// Solve the Phase I linear program.
	_, xOpt, newBasic, err := simplex(basicIdxs, c, aNew, b, 1e-10)
	if err != nil {
		return nil, nil, nil, fmt.Errorf("lp: error finding feasible basis: %s", err)
	}

	// The original LP is infeasible if the added variable has non-zero value
	// in the optimal solution to the Phase I problem.
	if math.Abs(xOpt[n]) > phaseIZeroTol {
		return nil, nil, nil, ErrInfeasible
	}

	// The basis found in Phase I is a feasible solution to the original LP if
	// the added variable is not in the basis.
	addedIdx := -1
	for i, v := range newBasic {
		if v == n {
			addedIdx = i
		}
		xb[i] = xOpt[v]
	}
	if addedIdx == -1 {
		extractColumns(ab, A, newBasic)
		return newBasic, ab, xb, nil
	}

	// The value of the added variable is in the basis, but it has a zero value.
	// See if exchanging another variable into the basic set finds a feasible
	// solution.
	basicMap := make(map[int]struct{})
	for _, v := range newBasic {
		basicMap[v] = struct{}{}
	}
	var set bool
	for i := range xOpt {
		if _, inBasic := basicMap[i]; inBasic {
			continue
		}
		newBasic[addedIdx] = i
		if set {
			mat.Col(col, i, A)
			ab.SetCol(addedIdx, col)
		} else {
			extractColumns(ab, A, newBasic)
			set = true
		}
		err := initializeFromBasic(xb, ab, b)
		if err == nil {
			return newBasic, ab, xb, nil
		}
	}
	return nil, nil, nil, ErrInfeasible
}

• より複雑な問題で自明な初期基底解がない場合は、「フェーズI」と呼ばれる補助問題を解いて初期実行可能基底解を見つけます。これは findInitialBasic 関数の中で実装されています。

フェーズ I 問題の箇所

	// It may be that this linearly independent basis is also a feasible set. If
	// so, the Phase I problem can be avoided.
	ab := mat.NewDense(m, len(basicIdxs), nil)
	extractColumns(ab, A, basicIdxs)
	xb := make([]float64, m)
	err := initializeFromBasic(xb, ab, b)
	if err == nil {
		return basicIdxs, ab, xb, nil
	}

	// This set was not feasible. Instead the "Phase I" problem must be solved
	// to find an initial feasible set of basis.
	//
	// Method: Construct an LP whose optimal solution is a feasible solution
	// to the original LP.
	// 1) Introduce an artificial variable x_{n+1}.
	// 2) Let x_j be the most negative element of x_b (largest constraint violation).
	// 3) Add the artificial variable to A with:
	//      a_{n+1} = b - \sum_{i in basicIdxs} a_i + a_j
	//    swap j with n+1 in the basicIdxs.
	// 4) Define a new LP:
	//   minimize  x_{n+1}
	//   subject to [A A_{n+1}][x_1 ... x_{n+1}] = b
	//          x, x_{n+1} >= 0
	// 5) Solve this LP. If x_{n+1} != 0, then the problem is infeasible, otherwise
	// the found basis can be used as an initial basis for phase II.
	//
	// The extra column in Step 3 is defined such that the vector of 1s is an
	// initial feasible solution.

• 初期化が成功すると、basicIdxs に基底変数の列インデックスが、ab に対応する A の列からなる基底行列が格納されます。xb に基底変数の値が格納されます。
• nonBasicIdx には非基底変数のインデックスが格納されます。

フェーズ I 問題とは

解説

なぜフェーズIが必要か？

シンプレックス法を開始するためには、まず「実行可能基底解（feasible basic solution）」が必要です。これは、すべての変数（元の変数とスラック変数）が非負であり、かつ連立等式制約を満たすような基底解、つまり実行可能領域の頂点の一つを見つけることです。

例題のように、制約条件がすべて「$\leq$」かつ右辺が非負の場合、スラック変数を基底変数、元の変数を非基底変数（ゼロに固定）とすることで、簡単に初期実行可能基底解（この例では $(x,y,s_
1 ,s_2)=(0,0,30,40)$）を見つけることができます。しかし、制約条件に「$\geq$」や等式が含まれる場合、あるいは右辺が負の値を持つ場合、このような自明な初期実行可能基底解が存在しないことがあります。

このような場合に、シンプレックス法を開始するための初期実行可能基底解を見つけるために用いられるのが「フェーズI」です。

フェーズI問題の考え方

フェーズIの目的は、元の線形計画問題の実行可能解を見つけること です。そのために、補助的な線形計画問題を作成します。この補助問題は、必ず自明な初期実行可能基底解を持つように構築され、その最適解が元の問題の実行可能解と密接に関連するように設計されています。

補助問題は、元の問題の制約条件を緩和するために「人工変数 (artificial variables)」を導入し、これらの人工変数の合計を最小化するという目的関数を持ちます。もし、この補助問題の最適解においてすべての人工変数の値がゼロになるならば、それは元の問題の実行可能解に対応します。逆に、補助問題の最適解でいずれかの人工変数が正の値を持つ場合、それは元の問題には実行可能解が存在しない、すなわち実行不可能であることを意味します。

gonum の findInitialBasic 関数におけるフェーズIの実装

findInitialBasic 関数は、与えられた A と b に対して初期実行可能基底解を見つけようとします。

• 1. 線形独立な列の探索:
  • まず、行列 A から m 個の線形独立な列（基底となりうる列）を探します (findLinearlyIndependent 関数を使用）。
• 1. 簡易な初期基底の確認:
  • 見つかった線形独立な列を基底行列 ab とし、ab * xb = b を解いて基底変数の値 xb を求めます。もし、この xb のすべての要素が非負であれば（initializeFromBasic 関数で確認）、それはそのまま実行可能基底解となるため、フェーズIはスキップされます。
• 1. フェーズIの開始:
  • 上記の簡易な方法で実行可能基底解が見つからなかった場合、フェーズIに進みます。
  • 人工変数の導入:
    • 元の問題の変数に加えて、x_{n+1} という人工変数（コード内では n 番目の列として扱われます）を導入します。
  • 補助問題の構築:
    • 元の制約条件 A * x = b に対して、人工変数を含む新しい制約条件を作成します。gonum の実装では、少し特殊な方法で人工変数の列 a_{n+1} を定義しています。これは、元の問題で最も制約違反が大きい（最も値が負である）基底変数 x_j を特定し、その列 a_j を使って $a_{n+1} = b - \sum_{i \in basicIdxs, i \neq j} a_i$ となるように計算しています。そして、この人工変数 x_{n+1} を、元の基底変数の中で最も値が負だった x_j と入れ替えて新しい基底とします。これにより、人工変数を含む拡張された問題に対して、自明な初期実行可能基底解（基底変数の値が1であるベクトル）が得られるように工夫されています。
  • フェーズIの目的関数:
    • 人工変数 x_{n+1} を最小化するという目的関数を設定します。新しい目的関数ベクトル c は、人工変数に対応する要素のみが 1、それ以外は 0 となります。
  • フェーズI問題の解決:
    • 構築した補助問題（拡張された A と新しい c、そして元の b を持つ問題）を、再帰的に simplex 関数を呼び出して解きます。
• 1. フェーズIの結果判定:
  • 補助問題の最適解 xOpt が得られたら、人工変数 xOpt[n] の値を調べます。
  • math.Abs(xOpt[n]) > phaseIZeroTol の場合（人工変数の値がゼロでない場合）：
    • 元の問題には実行可能解が存在しないため、ErrInfeasible を返します。
  • math.Abs(xOpt[n]) <= phaseIZeroTol の場合（人工変数の値がほぼゼロである場合）：
    • 補助問題の最適解は元の問題の実行可能解に対応します。このとき得られた基底 (newBasic) が、フェーズII（元の目的関数を最小化する段階）を開始するための初期実行可能基底となります。
• 1. 人工変数が基底に残っている場合の処理:
  • もし人工変数が基底変数として残っているがその値がゼロである場合、これは縮退した状況です。コードは、基底に残っている人工変数を、元の問題の非基底変数の中から実行可能解を崩さないものと交換しようと試みます。これが成功すれば、人工変数を含まない実行可能基底が得られます。

3. シンプレックス法の主反復 (simplex 関数のループ内) "フェーズ2" 問題

初期実行可能基底解が見つかると、シンプレックス法の主アルゴリズムのループに入ります。ここでは、現在の基底解から目的関数値を改善する方向へ、実行可能領域の辺に沿って隣の頂点へ移動する操作を繰り返します。これは、手計算で1回目の反復、2回目の反復と進めたステップに相当します。

• 被約費用の計算:
  • 現在の基底における非基底変数の被約費用 ("resuced costs" r) を計算します。被約費用は、対応する非基底変数を1単位だけ増加させたときに目的関数値がどれだけ変化するかを示します（最小化問題では、被約費用が負であれば目的関数は減少します）。
  • 手計算で目的関数を非基底変数で表し、その係数を見た部分に相当します。

r を計算している箇所

		// Compute reduced costs -- r = cn - anᵀ ab¯ᵀ cb
		var tmp mat.VecDense
		err = tmp.SolveVec(ab.T(), cbVec)
		if err != nil {
			break
		}
		data := make([]float64, n-m)
		tmp2 := mat.NewVecDense(n-m, data)
		tmp2.MulVec(an.T(), &tmp)
		floats.SubTo(r, cn, data)

• 最適性の判定:
  • 計算されたすべての非基底変数の被約費用が非負（最小化問題の場合）であれば、現在の基底解が最適解です。これ以上目的関数値を改善することはできないため、ループを終了します。
  • コードでは、最小の被約費用が -tol 以上であるかで判定します。
    • 許容誤差（制約条件違反を受け入れるための許容値（"tolerance" tol ））
  • 手計算で被約費用がすべて非負になった時点で計算を終了した判断に相当します。

最適性の判定箇所

		// Replace the most negative element in the simplex. If there are no
		// negative entries then the optimal solution has been found.
		minIdx := floats.MinIdx(r)
		if r[minIdx] >= -tol {
			break
		}

• 投入変数の選択:
  • もし負の被約費用が存在する場合、最も負の値を持つ被約費用に対応する非基底変数を選択し、基底に投入する変数とします。
  • 手計算で目的関数係数が負の非基底変数の中から最も負の係数を持つものを選んだステップに相当します。
• 移動距離の計算と離脱変数の選択 (computeMove 関数):
  • 投入変数をどれだけ増加させられるかを計算します。これは、投入変数を増やしていく際に、どの基底変数が最初にゼロになるか、という限界によって決まります。

呼び出し箇所

		// Compute the moving distance.
		err = computeMove(move, minIdx, A, ab, xb, nonBasicIdx)
		if err != nil {
			if err == ErrUnbounded {
				return math.Inf(-1), nil, nil, ErrUnbounded
			}
			break
		}

conputeMove 関数定義

// computeMove computes how far can be moved replacing each index. The results
// are stored into move.
func computeMove(move []float64, minIdx int, A mat.Matrix, ab *mat.Dense, xb []float64, nonBasicIdx []int) error {
	// Find ae.
	col := mat.Col(nil, nonBasicIdx[minIdx], A)
	aCol := mat.NewVecDense(len(col), col)

	// d = - Ab^-1 Ae
	nb, _ := ab.Dims()
	d := make([]float64, nb)
	dVec := mat.NewVecDense(nb, d)
	err := dVec.SolveVec(ab, aCol)
	if err != nil {
		return ErrLinSolve
	}
	floats.Scale(-1, d)

	for i, v := range d {
		if math.Abs(v) < dRoundTol {
			d[i] = 0
		}
	}

	// If no di < 0, then problem is unbounded.
	if floats.Min(d) >= 0 {
		return ErrUnbounded
	}

	// move = bhat_i / - d_i, assuming d is negative.
	bHat := xb // ab^-1 b
	for i, v := range d {
		if v >= 0 {
			move[i] = math.Inf(1)
		} else {
			move[i] = bHat[i] / math.Abs(v)
		}
	}
	return nil
}

• • computeMove 関数がこの計算を行い、各基底変数について許容される投入変数の増加量を move に格納します。
  • この増加量が最小となる基底変数が、基底から離脱する変数となります。もし投入変数を無限に増やせる場合（computeMove が ErrUnbounded を返す場合）、問題は非有界です。
  • 手計算で $\theta$ を計算し、離脱変数を選択したステップに相当します。
• 縮退の処理 (replaceBland 関数):
  • 計算された移動距離の最小値がゼロである場合、縮退が発生しています。縮退時に対策を講じないと、アルゴリズムが同じ基底を繰り返し計算する「サイクリング」に陥る可能性があります。
  • replaceBland 関数は、Bland の規則と呼ばれるルールに従って投入変数と離脱変数を選択し、サイクリングを回避しようとします。
  • 手計算では通常意識しない部分ですが、実装上重要な考慮事項です。

呼び出し箇所

		// Replace the basic index along the tightest constraint.
		replace := floats.MinIdx(move)
		if move[replace] <= 0 {
			replace, minIdx, err = replaceBland(A, ab, xb, basicIdxs, nonBasicIdx, r, move)
			if err != nil {
				if err == ErrUnbounded {
					return math.Inf(-1), nil, nil, ErrUnbounded
				}
				break
			}
		}

`replaceBland` 関数定義

// replaceBland uses the Bland rule to find the indices to swap if the minimum
// move is 0. The indices to be swapped are replace and minIdx (following the
// nomenclature in the main routine).
func replaceBland(A mat.Matrix, ab *mat.Dense, xb []float64, basicIdxs, nonBasicIdx []int, r, move []float64) (replace, minIdx int, err error) {
	m, _ := A.Dims()
	// Use the traditional bland rule, except don't replace a constraint which
	// causes the new ab to be singular.
	for i, v := range r {
		if v > -blandNegTol {
			continue
		}
		minIdx = i
		err = computeMove(move, minIdx, A, ab, xb, nonBasicIdx)
		if err != nil {
			// Either unbounded or something went wrong.
			return -1, -1, err
		}
		replace = floats.MinIdx(move)
		if math.Abs(move[replace]) > blandZeroTol {
			// Large enough that it shouldn't be a problem
			return replace, minIdx, nil
		}
		// Find a zero index where replacement is non-singular.
		biCopy := make([]int, len(basicIdxs))
		for replace, v := range move {
			if v > blandZeroTol {
				continue
			}
			copy(biCopy, basicIdxs)
			biCopy[replace] = nonBasicIdx[minIdx]
			abTmp := mat.NewDense(m, len(biCopy), nil)
			extractColumns(abTmp, A, biCopy)
			// If the condition number is reasonable, use this index.
			if mat.Cond(abTmp, 1) < 1e16 {
				return replace, minIdx, nil
			}
		}
	}
	return -1, -1, ErrBland
}

• ピボット操作 (基底の更新):
  • 選択された投入変数と離脱変数のインデックスを、basicIdxs と nonBasicIdx の間で交換し、対応する目的関数係数も交換します。
  • また、基底行列 ab と非基底行列 an も新しい基底に合わせて更新します。
  • これは手計算で連立方程式を解き直して基底変数を非基底変数で表し直したステップに相当します。

計算箇所

		// Replace the constrained basicIdx with the newIdx.
		basicIdxs[replace], nonBasicIdx[minIdx] = nonBasicIdx[minIdx], basicIdxs[replace]
		cb[replace], cn[minIdx] = cn[minIdx], cb[replace]
		tmpCol1 := mat.Col(nil, replace, ab)
		tmpCol2 := mat.Col(nil, minIdx, an)
		ab.SetCol(replace, tmpCol2)
		an.SetCol(minIdx, tmpCol1)

• 新しい基底解の計算:
  • 更新された基底行列 ab を用いて、新しい基底変数 xb の値を ab * xb = b を解くことで再計算します。
  • 手計算で新しい非基底変数をゼロとおいて基底変数の値を求めたステップに相当します。

計算箇所

		// Compute the new xb.
		xbVec := mat.NewVecDense(len(xb), xb)
		err = xbVec.SolveVec(ab, bVec)
		if err != nil {
			break
		}

これらのステップを繰り返し、最適解に到達するまで基底解を改善していきます。

4. 最終結果

ループが終了したら、現在の基底解が最適解です。

• 最適な目的関数値は、現在の基底変数に対応する目的関数係数 cb と基底変数の値 xb の内積として計算されます。
• 最終的な解ベクトル xopt は、元の変数の数と同じ次元で作成され、基底変数に対応する位置に xb の値が格納され、非基底変数の位置はゼロとなります。

計算箇所

	// Found the optimum successfully or died trying. The basic variables get
	// their values, and the non-basic variables are all zero.
	opt := floats.Dot(cb, xb)
	xopt := make([]float64, n)
	for i, v := range basicIdxs {
		xopt[v] = xb[i]
	}
	return opt, xopt, basicIdxs, err

参考文献

Python ではじめる数理最適化（第2版）
しっかり学ぶ数理最適化　モデルからアルゴリズムまで
gonum/optimize/convex/lp documentation
lectures 11-13 of UC Math 352