最適化問題を解く（その04/X） - 「大人の教養・知識・気付き」を伸ばすブログ

　単体法では、初回の実行可能基底形式が必要だった。しかし一般には、そうした実行可能基底形式が標準形から直接見つかるとは限らない。そこで2段解法が提案された。
　2段解法は、第1段階で実行可能性を判定し実行可能基底形式を求める。第2段階では、第1段階で得られた基底形式を初回の実行可能基底形式として、単体法のアルゴリズムを適用して最適解を求める。2段解法では、第1段階と第2段階のいずれでも単体法が用いられる。

　線形計画問題の標準形

$\begin{aligned} \mathrm{Minimize}\ \ &\boldsymbol{c}^{\mathrm{T}}\boldsymbol{x}\\ \mathrm{subject\ to}\ \ &A\boldsymbol{x}=\boldsymbol{b},\\ &\boldsymbol{x}\geq\boldsymbol{0} \end{aligned}$

が与えられたとき、第1段階では、強制的に実行可能基底形式を作るために人為変数と呼ばれる非負変数 $v_1,\cdots,v_m$ を導入し、標準形の等式制約を以下のように変形する。

$\begin{aligned} \begin{cases} a_{11}x_{1}+\cdots+a_{1n}x_{n}&+v_1& & & &=b_1\geq0\\ a_{21}x_{1}+\cdots+a_{2n}x_{n}& &+v_2& & &=b_2\geq0\\ \vdots\ \ \ \ \ \ \ +\ \ \ \ \ \ \vdots& & &\ddots& &\vdots\\ a_{m1}x_{1}+\cdots+a_{mn}x_n{}& & & &v_m&=b_1\geq0\\ x_i\geq0,& & & & & i=1,\cdots,n,\\ v_i\geq0,& & & & & i=1,\cdots,m \end{cases} \end{aligned}$

ただし右辺 $b_i,$ $i=1,\cdots,m$ はすべて非負であるとする。そうでない場合は両辺に $-1$ を掛けて $b_i\geq0$ とする。ここで

$\begin{aligned} (x_1,\cdots,x_n,v_1,\cdots,v_m)=(0,\cdots,0,b_1,\cdots,b_m) \end{aligned}$

が上の変形式において実行可能基底解になることは明らかである。この問題に単体法を適用して人為変数がすべて零になるような実行可能解 $(\bar{x}_1,$ $\cdots$ $\bar{x}_n,$ $0,$ $\cdots,$ $0)$ が見つかれば、 $(\bar{x}_1,$ $\cdots$ $\bar{x}_n$ が元の線形計画問題の実行可能基底解になる。このためには目的関数 $w$ の代わりに

$\begin{aligned} u=\displaystyle{\sum_{i=1}^{m}v_i} \end{aligned}$

を最小化することを考える。元の問題が実行可能解を持つならば、 $u$ の最小値 $u^{*}$ は零であり、そのとき人為変数 $v_1,$ $\cdots,$ $v_m$ はすべて零になる。他方でもし最小値が $u^{*}\gt0$ を満たすならば、すべての人為変数を零にすることはできないため、元の問題には実行可能解が存在しないことが分かる。ここまでが第1段階である。
　変形式それぞれにおいて $v_i$ について解き、これを $u=\displaystyle{\sum_{i=1}^{m}v_i}$ に代入することで

$\begin{aligned} u=\left(-\displaystyle{\sum_{i=1}^{m}a_{i1}}\right)x_1+\cdots+\left(-\displaystyle{\sum_{i=1}^{m}a_{in}}\right)x_n+\displaystyle{\sum_{i=1}^{m}b_{i}} \end{aligned}$

を得る。ここで

$\begin{aligned} d_j&=-\left(\displaystyle{\sum_{i=1}^{m}a_{ij}}\right),j=1,\cdots,n,\\ u_0&=\displaystyle{\sum_{i=1}^{m}b_i}\geq0 \end{aligned}$

とおいて、変形式に

$\begin{aligned} -u+\displaystyle{\sum_{i=1}^{n}d_i x_i}=-u_0 \end{aligned}$

を追加すれば実行可能基底形式が得られる。また第1段階では $w$ を考慮する必要はないものの、第2段階で考慮すべく、当初から第1段階の基底形式に

$\begin{aligned} -w+\displaystyle{\sum_{i=1}^{n}c_i x_i}=0 \end{aligned}$

を追加しておく。
　線形計画問題の標準形が退化していない場合、各反復で新しい実行可能基底解に移る際に必ず目的関数の値が減少するため、同じ実行可能基底解が選ばれることはない。また実行可能基底解の選び方は高々 ${}_n\!C_m$ 通りしかないため、単体法は有限回の反復で最適解を得る、もしくは解が有界でないと判定して終了することが分かる。

　他方で、退化している場合、目的関数値が減少せず、実行可能基底解の入れ替えを繰り返して同じ実行可能基底解を選ぶことがあり得る。このために単体法が終了しない現象が生じ得、これを巡回という。これはピボットを選択する際に $p,q$ を選ぶ自由度があることが原因である。

　巡回対策を施した単体法のアルゴリズムとして、 $\mathrm{Bland}$ の巡回対策を用いた単体法がある：

	(1)	相対費用係数が $\bar{\boldsymbol{c}}_{N}\geq\boldsymbol{0}$ を満たすならば最適解を得たことになるため終了する。そうでなければ(2)に移る。
	(2)	$i=\displaystyle{\mathrm{arg}\ \min_{\bar{c}_i}i}$ とする。
	(3)	$\bar{\boldsymbol{a}}_q\leq\boldsymbol{0}$ ならば、目的関数が有界でないため終了する。そうでなければ(4)に移る。
	(4)	$\displaystyle{\min_{\bar{a}_{iq}\gt0}\frac{\bar{b}_i}{\bar{a}_{iq}}}=\displaystyle{\frac{\bar{b}_p}{\bar{a}_{pq}}}$ を満たすような添字 $p$ を求める。もし最小値を与えるような $i$ が複数存在する場合、そのうち最小のものを $p$ とする。
	(5)	$\bar{a}_{pq}$ をピボットとする掃き出しを実行して、 $x_p$ の代わりに $x_q$ を基底変換とする新しい実行可能基底形式を構築する。
	(6)	(1)に戻る。

$\mathrm{Bland}$ の巡回対策を用いた単体法は退化しているか否かに関わらず有限回の反復で終了することが知られている。
　単体法はたいてい比較的少ない回数で終了するものの、最適解に到達するまでにすべての端点を巡らなければならない例が存在することが知られている。とはいえそのような事態は殆ど生じ得ない。

2.3.3　内点法

　単体法と対照的に、内点法は $\boldsymbol{x}$ の成分がすべて正であるような点を辿りながら最適解に到達する方法で、1984年に $\mathrm{Karmarkar}$ によって提案されたアルゴリズムを契機に大きく発展したものである。以下では内点法の中でも主双対内点法を考える。
　標準的な形態での線形計画問題

$\begin{aligned} \mathrm{Minimize}\ \ &{}^{t}\boldsymbol{c}\boldsymbol{x}\\ \mathrm{subject\ to}\ \ &A\boldsymbol{x}=\boldsymbol{b},\\ &\boldsymbol{x}\geq\boldsymbol{0} \end{aligned}$

および双対問題

$\begin{aligned} \mathrm{Maximize}\ \ &{}^{t}\boldsymbol{b}\boldsymbol{y}\\ \mathrm{subject\ to}\ \ &{}^{t}A\boldsymbol{y}+\boldsymbol{s}=\boldsymbol{c}\\ &\boldsymbol{s}\geq\boldsymbol{0} \end{aligned}$

の最適性条件は、

$\begin{aligned} A\boldsymbol{x}=\boldsymbol{b},\\ {}^{t}A\boldsymbol{y}+\boldsymbol{s}=\boldsymbol{c},\\ x_js_j=0,j=1,\cdots,n,\\ \boldsymbol{x}\geq\boldsymbol{0},\boldsymbol{s}\geq\boldsymbol{0} \end{aligned}$

で与えられる。内点法では、このうち相補性条件 $x_js_j=0,j=1,\cdots,n$ にパラメータ $\nu\gt0$ を導入して変形した

$\begin{aligned} A\boldsymbol{x}=\boldsymbol{b},\\ {}^{t}A\boldsymbol{y}+\boldsymbol{s}=\boldsymbol{c},\\ x_js_j=\nu,j=1,\cdots,n,\\ \boldsymbol{x}\geq\boldsymbol{0},\boldsymbol{s}\geq\boldsymbol{0} \end{aligned}$

を考える。
　この条件で $\nu\rightarrow0$ としたときの極限が最適性条件に帰着できる。ここで $\nu\gt0$ を固定すると、上記の解は一意に存在することが知られている。この解を $(\boldsymbol{x}(\nu),\boldsymbol{y}(\nu),\boldsymbol{s}(\nu))$ とおくと、 $\nu\gt0$ を変化させたときのそれらの軌跡は実行可能領域の内部を通る滑らかな曲線となることが知られているが、この軌跡は中心曲線と呼ぶ。中心曲線は $\nu\rightarrow0$ で線形計画問題の最適解に収束する。内点法は、この中心曲線を近似的にトレースすることで最適解を得る手法である。

　まず初期点 $(\boldsymbol{x}^{(0)},\boldsymbol{y}^{(0)},\boldsymbol{s}^{(0)})$ を $\boldsymbol{x}^{(0)},\boldsymbol{s}^{(0)}$ の成分がすべて正であるように選択する。そして $k=0,1,2,\cdots$ に対して中心曲線上の点

$\begin{aligned} \begin{bmatrix} \boldsymbol{x}\left(\nu^{(k)}\right)\\ \boldsymbol{y}\left(\nu^{(k)}\right)\\ \boldsymbol{s}\left(\nu^{(k)}\right) \end{bmatrix}=\begin{bmatrix} \boldsymbol{x}^{(k)}\\ \boldsymbol{y}^{(k)}\\ \boldsymbol{s}^{(k)} \end{bmatrix}+ \begin{bmatrix} \Delta\boldsymbol{z}\\ \Delta\boldsymbol{y}\\ \Delta\boldsymbol{s} \end{bmatrix} \end{aligned}$

を考える。ただし $\nu^{(k)}\gt0$ は $\nu^{(k-1)}$ より小さくなるように選択する。これを最適性条件に代入することで

$\begin{aligned} A\Delta\boldsymbol{x}=\boldsymbol{b}-A\boldsymbol{x}^{(k)},\\ {}^{t}A\Delta\boldsymbol{y}+\Delta\boldsymbol{s}=\boldsymbol{c}-{}^{t}A\boldsymbol{y}^{(k)}-\boldsymbol{s}^{(k)},\\ s_j^{(k)}\Delta x_j+x_j^{(k)}\delta s_j=\nu^{(k)}-x_j^{(k)}s_j^{(k)}-\Delta x_j\Delta s_j,j=1,\cdots,n \end{aligned}$

が得られる。ここで近似として最後の条件を無視すると、

$\begin{aligned} s_j^{(k)}\Delta x_j+x_j^{(k)}+x_j^{(k)}\Delta s_j=\nu^{(k)}-x_j^{(k)}s_j^{(k)},j=1,\cdots,n \end{aligned}$

であり、これらをまとめると $(\Delta\boldsymbol{x},\Delta\boldsymbol{y},\Delta\boldsymbol{s})$ を未知数とする連立方程式

$\begin{aligned} \begin{cases} A\Delta\boldsymbol{x}=\boldsymbol{b}-A\boldsymbol{x}^{(k)},\\ {}^{t}A\Delta\boldsymbol{y}+\Delta\boldsymbol{s}=\boldsymbol{c}-{}^{t}A\boldsymbol{y}^{(k)}-\boldsymbol{s}^{(k)},\\ s_j^{(k)}\Delta x_j+x_j^{(k)}+x_j^{(k)}\Delta s_j=\nu^{(k)}-x_j^{(k)}s_j^{(k)},j=1,\cdots,n \end{cases} \end{aligned}$

であり、その解は探索方向と呼ぶ。探索方向を $(\Delta\boldsymbol{x}^{(k)},\Delta\boldsymbol{y}^{(k)},\Delta\boldsymbol{s}^{(k)})$ で表すと、内点法の位置反復では現在の点を

$\begin{aligned} \begin{bmatrix} \boldsymbol{x}^{(k+1)}\\ \boldsymbol{y}^{(k+1)}\\ \boldsymbol{s}^{(k+1)} \end{bmatrix}= \begin{bmatrix} \boldsymbol{x}^{(k)}\\ \boldsymbol{y}^{(k)}\\ \boldsymbol{s}^{(k)} \end{bmatrix}+\alpha_k \begin{bmatrix} \Delta\boldsymbol{x}^{(k)}\\ \Delta\boldsymbol{y}^{(k)}\\ \Delta\boldsymbol{s}^{(k)} \end{bmatrix} \end{aligned}$

と更新する。ここで $\alpha_k\gt0$ はステップ幅と呼ばれ、 $\boldsymbol{x}^{(k+1)},\boldsymbol{s}^{(k+1)}$ の成分がすべて正になる範囲の値に決められる。

2.3.4　内点法の多項式性

　内点法の本質部分に当たる以下の定理を議論する。

定理2.7　内点法の定理　主双対内点法において、 $\gamma=1-\displaystyle{\frac{0.1}{\sqrt{n}}}$ とし、ステップ幅 $\alpha_k\equiv1$ としたものを考える。初期点 $(\boldsymbol{x}_0,\boldsymbol{y}_0,\boldsymbol{s}_0)$ が、 $X_0=\mathrm{diag}(\boldsymbol{x}_0)$ として条件

$\begin{aligned} \|X_0\boldsymbol{s}_0-\rho_0\boldsymbol{1}\|\leq0.2\rho_0 \end{aligned}$

を満たすものとする。このとき主双対内点は、すべての $k=0,1,\cdots$ について、 $X_k=\mathrm{diag}(\boldsymbol{x}_k)$ として条件

$\begin{aligned} (i)&{}^{t}\boldsymbol{x}_{k+1}\boldsymbol{s}_{k+1}=\left(1-\displaystyle{\frac{0.1}{\sqrt{n}}}\right){}^{t}\boldsymbol{x}_k\boldsymbol{s}_k,\\ (ii)&\|X_k\boldsymbol{s}_k-\rho_k\boldsymbol{1}\|\leq0.2\rho_k \end{aligned}$

を満たすような内点実行可能解の点列を生成する。

この(i)は双対ギャップ ${}^{t}\boldsymbol{x}_k\boldsymbol{s}_k$ がデータの次元のみに依存する1次収束をすることを意味する。これがアルゴリズムの多項式性を担う本質的部分である。(ii)は $\boldsymbol{x}_k,$ $\boldsymbol{s}_k$ を要素ごとに掛けたベクトル $X_k\boldsymbol{s}_k$ の各要素が概ね同水準にあることを示し、これが $(\boldsymbol{x}_k,\boldsymbol{y}_k,\boldsymbol{s}_k)$ が中心曲線の近くに存在することを具体的に表現している。
　(i)より、双対ギャップをある水準 $g_1$ から $g_2$ にまで減らすのに反復回数は $10\sqrt{n}\log\left(\displaystyle{\frac{g_2}{g_1}}\right)$ で抑えることができる。また1回の反復に必要な計算量は約 $O(n^2m)$ であり、その主要な部分は探索方向の計算である。

2.3.5　Pythonによる実装例

####################
### 線形計画問題 ###
####################

import cvxpy as cp

x1, x2 = cp.Variable(), cp.Variable() # 決定変数の定義
obj = cp.Maximize(20 * x1 + 60 * x2)  # 目的関数
cons = [5 * x1 + 4 * x2 <= 80,
        2 * x1 + 4 * x2 <= 40,
        2 * x1 + 8 * x2 <= 64,
        x1 >= 0,
        x2 >= 0]                      # 制約
P = cp.Problem(obj, cons)             # 最適化問題の定義
P.solve(verbose = True)               # 最適化の実施
print(x1.value, x2.value)             # 最適解の実施
print(P.value)                        # 最適解における目的関数値


# numpyによるベクトル・行列表記

import cvxpy as cp
import numpy as np

x = cp.Variable(2) # 変数の数を引数に
c = np.array([20.0, 60.0]).T
G = np.array([[5.0, 4.0],
              [2.0, 4.0],
              [2.0, 8.0],
              [-1.0, 0],
              [0, -1.0]])      # 制約
h = [80.0, 40.0, 64.0, 0.0, 0.0]
obj = cp.Maximize(c.T @ x)
cons = [G @ x <= h]
P = cp.Problem(obj, cons)
P.solve(verbose = True)

print(x.value[0])
print(x.value[1])
print(P.value)

参考文献

寒野善博[著]/駒木文保[編](2019)「最適化手法入門」(講談社)
東京大学工学教程編纂委員会[編]、寒野　善博[著]、土谷　隆[著](2014)「東京大学工学教程基礎系数学最適化と変分法」(丸善出版)
矢部博(2006)「工学基礎　最適化とその応用」(数理工学社)
Kochenderfer, Mykel J., Wheeler, Tim A. (2019), "Algorithms for Optimization", The MIT Press (https://algorithmsbook.com/optimization/files/optimization.pdf)

*1:双対性は別に章を立てて扱う予定である。

はじめに

前回

2. 線形計画と凸2次計画

2.3 線形計画問題の解法

2.3.2 2段解法

2.3.3 内点法

2.3.4 内点法の多項式性

2.3.5 Pythonによる実装例

参考文献

2.　線形計画と凸2次計画

2.3　線形計画問題の解法

2.3.2　2段解法

2.3.3　内点法

2.3.4　内点法の多項式性

2.3.5　Pythonによる実装例