Megiddo-Dyerのアルゴリズム ― 2変数線形計画問題の高速解法

まえがき

線形計画問題の解法としては，単体法あるいはそれをベースとした2段階法が最も有名，かつ変数が数百個程度の問題ならば十分実用的です．これについては別記事に記しました．

変数の数がもっと多いならば，内点法（主双対内点法）が有利になります．

逆に変数の数が少ない問題に対しては，それに特化した効率の良い解法もあります．Megiddo-Dyerのアルゴリズムはそのようなものの一つです．

• Nimrod Megiddo, Linear-Time Algorithms for Linear Programming in $R^{3}$ and Related Problems, SIAM Journal on Computing, 12(4):759--776 (1983).

• Martin E. Dyer, Linear Time Algorithms for Two- and Three-Variable Linear Programs, SIAM Journal on Computing, 13(1):31--45 (1984).

単体法も内点法も，不等式制約条件はスラック変数によって等式制約条件に変形して扱うので，必然的に(制約条件の数)$<$(変数の数)となります．Megiddo-Dyerのアルゴリズムではそのような考え方はせず，不等式制約条件を陽に扱います．制約条件の数は多いが変数の数は少ない，という問題に向きます． 本記事では特に，変数の数が2個である場合に絞ってこれを説明します．なお，元論文で用いられていない記号・表記を本記事で独自に用いている箇所があります．

準備

変数の数が2個である線形計画問題は，一般的に次のように書けます．

\displaylines{
(x_{1}^{*},x_{2}^{*})=\mathop{\mathrm{arg~min}}_{(x_{1},x_{2})}\left\{
c_{1}x_{1}+c_{2}x_{2}
\left|(x_{1},x_{2})\in\mathcal{D}
\right.
\right\}
\\
\mathcal{D}=\left\{
(x_{1},x_{2})\left|a_{k1}x_{1}+a_{k2}x_{2}\leq b_{k}\quad\mbox{for}~\forall k=1,\cdots,N
\right.\right\}
}

ただし，$(x_{1}^{*},x_{2}^{*})$は最適解，$N$は不等式制約条件の数です．$\mathcal{D}$は実行可能領域と呼ばれ，$(x_{1},x_{2})$が$\mathcal{D}$に含まれるときに「$(x_{1},x_{2})$は実行可能である」と言います．
$c_{1}x_{1}+c_{2}x_{2}$は評価関数（ほか目的関数，損失関数など）と呼ばれます．
$c_{1}$と$c_{2}$が同時に$0$になることは（問題が意味を持たなくなるため）無いので，$c_{2}\neq 0$と仮定しても一般性を失いません（$c_{2}=0$の場合には$x_{1}$と$x_{2}$を入れ替えれば良いです）．

ここで

\displaylines{
x=x_{1}
\\
y=c_{1}x_{1}+c_{2}x_{2}
}

とおけば，不等式$a_{k1}x_{1}+a_{k2}x_{2}\leq b_{k}$は次のように変形できます．

\displaylines{
\frac{a_{k2}}{c_{2}}>0~\mbox{のとき}\quad y\leq m_{k} x + y_{0k} \\
\frac{a_{k2}}{c_{2}}<0~\mbox{のとき}\quad y\geq m_{k} x + y_{0k} \\
\frac{a_{k2}}{c_{2}}=0~\mbox{かつ}~a_{k1}<0~\mbox{のとき}\quad x\geq \bar{x}_{k} \\
\frac{a_{k2}}{c_{2}}=0~\mbox{かつ}~a_{k1}>0~\mbox{のとき}\quad x\leq \bar{x}_{k}
}

ただし，

\displaylines{
m_{k}=c_{1}-\frac{a_{k1}c_{2}}{a_{k2}} \\
y_{0k}=\frac{b_{k}c_{2}}{a_{k2}} \\
\bar{x}_{k}=\frac{b_{k}}{a_{k1}}
}

とそれぞれおきました．
これに基づけば，実行可能領域を次の4つの部分領域に分けて考えることができます．

\displaylines{
\mathcal{D}=\left\{(x_{1},x_{2})\left|
(x_{1},y_{2})=\left(x,\frac{y-c_{1}x}{c_{2}}\right),~
(x,y)\in\mathcal{D}_{\mathrm{U}}\cap\mathcal{D}_{\mathrm{L}}\cap\mathcal{D}_{\mathrm{l}}\cap\mathcal{D}_{\mathrm{r}}
\right.\right\}
\\
\mathcal{D}_{\mathrm{U}}=\left\{ (x,y) \left| y\leq m_{k}x+y_{0k},~\forall k\in\mathcal{I}_{\mathrm{U}} \right.\right\}
,\qquad
\mathcal{I}_{\mathrm{U}}=\left\{ k \left| \frac{a_{k2}}{c_{2}}>0 \right.\right\}
\\
\mathcal{D}_{\mathrm{L}}=\left\{ (x,y) \left| y\geq m_{k}x+y_{0k},~\forall k\in\mathcal{I}_{\mathrm{L}} \right.\right\}
,\qquad
\mathcal{I}_{\mathrm{L}}=\left\{ k \left| \frac{a_{k2}}{c_{2}}<0 \right.\right\}
\\
\mathcal{D}_{\mathrm{l}}=\left\{ (x,y) \left| x\geq \bar{x}_{k},~\forall k\in\mathcal{I}_{\mathrm{l}} \right.\right\}
,\qquad
\mathcal{I}_{\mathrm{l}}=\left\{ k \left| \frac{a_{k2}}{c_{2}}=0, a_{k1}<0 \right.\right\}
\\
\mathcal{D}_{\mathrm{r}}=\left\{ (x,y) \left| x\leq \bar{x}_{k},~\forall k\in\mathcal{I}_{\mathrm{r}} \right.\right\}
,\qquad
\mathcal{I}_{\mathrm{r}}=\left\{ k \left| \frac{a_{k2}}{c_{2}}=0, a_{k1}>0 \right.\right\}
}

なお，$a_{k1}$と$a_{k2}$が同時に$0$になることは（制約条件が意味を持たなくなるので）無いと仮定して構いません．

$\mathcal{D}_{\mathrm{U}}$ および $\mathcal{D}_{\mathrm{L}}$は，それぞれ実行可能領域の上側境界および下側境界を決めています．もし$\mathcal{I}_{\mathrm{L}}=\emptyset$ならば，実行可能領域が下側境界を持たないので$y$を無限に小さくできることになります．このような状況は「線形計画問題が無限解を持つ」と言い，異常ケースですので即座に終了します．

一方，$\mathcal{D}_{\mathrm{l}}$および$\mathcal{D}_{\mathrm{r}}$は，それぞれ$x$の左側境界および右側境界を決めており，

\displaylines{
x_{\mathrm{min}}=\begin{cases}
-\infty & (\mathcal{I}_{\mathrm{l}}=\emptyset)
\\
\mathop{\mathrm{max}}\left\{\bar{x}_{k}\left| k\in\mathcal{I}_{\mathrm{l}}\right.\right\} & (\mathcal{I}_{\mathrm{l}}\neq\emptyset)
\end{cases} \\
x_{\mathrm{max}}=\begin{cases}
\infty & (\mathcal{I}_{\mathrm{r}}=\emptyset)
\\
\mathop{\mathrm{min}}\left\{\bar{x}_{k}\left| k\in\mathcal{I}_{\mathrm{r}}\right.\right\} & (\mathcal{I}_{\mathrm{r}}\neq\emptyset)
\end{cases}
}

とおけば，解きたい問題は次のように変形できます．

(x^{*},y^{*})=\mathop{\mathrm{arg~min}}_{(x,y)}\left\{
y \left|
(x,y)\in\mathcal{D}_{\mathrm{U}}\cap\mathcal{D}_{\mathrm{L}},~x_{\mathrm{min}}\leq x\leq x_{\mathrm{max}}
\right.\right\}

以降は，この問題を考えることにします．
なお，$(x_{1}^{*},x_{2}^{*})$は$(x^{*},y^{*})$から次のように求められます．

(x_{1}^{*},x_{2}^{*})=\left(x^{*},\frac{y^{*}-c_{1}x^{*}}{c_{2}}\right)

上図中央の薄紫色に塗られた多角形領域が実行可能領域であり，その上で$y$が最小となる一番下の頂点$(x^{*},y^{*})$が最適解です．上記の問題を解くことは，すなわちこの点を見つけること，ととらえれば良いです．

制約条件境界の直線群から最適解を直接求める

前節のように実行可能領域を4つの部分領域に分けて考えたとき，制約条件の境界となる直線群は$\bar{\mathcal{D}}_{\mathrm{U}}=\left\{y=m_{k}x+y_{0k}\left|\forall k\in\mathcal{I}_{\mathrm{U}}\right.\right\}$，$\bar{\mathcal{D}}_{\mathrm{L}}=\left\{y=m_{k}x+y_{0k}\left|\forall k\in\mathcal{I}_{\mathrm{L}}\right.\right\}$，$x=x_{\mathrm{min}}$，および$x=x_{\mathrm{max}}$と表せます．問題の性質上，最適解は

i) $\bar{\mathcal{D}}_{\mathrm{U}}$の1要素と$\bar{\mathcal{D}}_{\mathrm{L}}$の1要素の交点
ii) $\bar{\mathcal{D}}_{\mathrm{L}}$の2要素の交点
iii) $\bar{\mathcal{D}}_{\mathrm{L}}$の1要素と$x=x_{\mathrm{min}}$の交点
iv) $\bar{\mathcal{D}}_{\mathrm{L}}$の1要素と$x=x_{\mathrm{max}}$の交点

の中で，実行可能領域にあり，かつ$y$座標が最小のものとなります．

$\bar{\mathcal{D}}_{\mathrm{U}}$の要素数を$N_{\mathrm{U}}$，$\bar{\mathcal{D}}_{\mathrm{L}}$の要素数を$N_{\mathrm{L}}$とそれぞれおくと，交点i)は$N_{\mathrm{U}}N_{\mathrm{L}}$個，ii)は${}_{N_{\mathrm{L}}}C_{2}=N_{\mathrm{L}}(N_{\mathrm{L}}-1)/2$個，iii)は高々$N_{\mathrm{L}}$個，iv)は高々$N_{\mathrm{L}}$個あります．よって最適解の候補となる交点の数は，最大で$\displaystyle\left(N_{\mathrm{U}}+\frac{N_{\mathrm{L}}+3}{2}\right)N_{\mathrm{L}}$個となります．たとえば$N_{\mathrm{U}}=N_{\mathrm{L}}=2$ならば$9$個です．この程度ならば，全て調べ上げて直接的に最適解を見つけるのは大した計算ではありません．一方で$N_{\mathrm{U}}$や$N_{\mathrm{L}}$が大きい時は，交点をしらみつぶしに調べるのは大変なコストとなるので，調べるべき交点を効率よく絞り込んでいく工夫が必要です．

最適解の存在する可能性のある範囲を絞り込む

実行可能領域の上側境界$y=f_{\mathrm{U}}(x)$および下側境界$y=f_{\mathrm{L}}(x)$は，次のように表せます．

\displaylines{
f_{\mathrm{U}}(x)&=\mathrm{min}\left\{m_{k}x+y_{0k} \left|\,k\in \mathcal{I}_{\mathrm{U}}\right.\right\}
\\
f_{\mathrm{L}}(x)&=\mathrm{max}\left\{m_{k}x+y_{0k} \left|\,k\in \mathcal{I}_{\mathrm{L}}\right.\right\}
}

$\mathcal{I}_{\mathrm{U}}$と$\mathcal{I}_{\mathrm{L}}$は既知なので，与えられた$x$に対して$f_{\mathrm{U}}(x)$と$f_{\mathrm{L}}(x)$はどちらも$O(N)$で求まることに注意しましょう．

$y=f_{\mathrm{U}}(x)$は上に凸な折れ線関数（上図の赤太線），$y=f_{\mathrm{L}}(x)$は下に凸な折れ線関数（同青太線）で，実行可能領域はこれらに挟まれた凸領域に含まれます．すなわち実行可能領域が存在するならば，$x$が$-\infty$から$\infty$へと増加するにつれて，$f_{\mathrm{L}}(x)-f_{\mathrm{U}}(x)$は無限大から減少していき，どこかで正から負に変わり，その後に減少から増加に転じ，また負から正に戻って，また無限大にまで増加する，という変化を辿ります．このとき，上側境界$y=f_{\mathrm{U}}(x)$の勾配は必ず単調減少，下側境界$y=f_{\mathrm{L}}(x)$の勾配は必ず単調増加となります．この性質に着目することで，次のように最適解が存在し得る範囲を徐々に絞り込んでいくことができます．

まず，ある（$x_{\mathrm{min}}\leq \tilde{x}\leq x_{\mathrm{max}}$を満たす）$x=\tilde{x}$に対し$f_{\mathrm{L}}(\tilde{x})>f_{\mathrm{U}}(\tilde{x})$であった（つまり直線$x=\tilde{x}$が実行可能領域を通過しない）としましょう．
$f_{\mathrm{L}}(x)-f_{\mathrm{U}}(x)$が$x=\tilde{x}$において局所的に減少関数（つまり$f'_{\mathrm{L}}(\tilde{x})-f'_{\mathrm{U}}(\tilde{x})<0$）ならば$\tilde{x}$の右側に，逆に局所的に増加関数（つまり$f'\_{\mathrm{L}}(\tilde{x})-f'\_{\mathrm{U}}(\tilde{x})>0$）ならば$\tilde{x}$の左側に実行可能領域があるはずです．
前述の通り，$f_{\mathrm{U}}(\tilde{x})$を与える$k\in\mathcal{I}_{\mathrm{U}}$は$O(N)$で見つかるのですが，折れ線関数なので複数存在し得ます．$f_{\mathrm{L}}(\tilde{x})$を与える$k\in\mathcal{I}_{\mathrm{L}}$も同様です．そしてその$k$の数だけ傾き$m_{k}$も存在します．

$(\tilde{x},f_{\mathrm{U}}(\tilde{x}))$を通る（複数の）直線の傾きの最小値および最大値を$m_{\mathrm{Umin}}$および$m_{\mathrm{Umax}}$，$(\tilde{x},f_{\mathrm{L}}(\tilde{x}))$を通る（複数の）直線の傾きの最小値および最大値を$m_{\mathrm{Lmin}}$および$m_{\mathrm{Lmax}}$，とそれぞれおきましょう．$f_{\mathrm{L}}(x)-f_{\mathrm{U}}(x)$の変化率は，$x=\tilde{x}$の前後で$m_{\mathrm{Lmin}}-m_{\mathrm{Umax}}$から$m_{\mathrm{Lmax}}-m_{\mathrm{Umin}}$に変わることになります．したがって，

i) $m_{\mathrm{Lmax}}<m_{\mathrm{Umin}}$ならば$f_{\mathrm{L}}(x)-f_{\mathrm{U}}(x)$は$x=\tilde{x}$から先で減少関数なので，実行可能領域は$\tilde{x}$よりも右側にある
ii) $m_{\mathrm{Lmin}}>m_{\mathrm{Umax}}$ならば$f_{\mathrm{L}}(x)-f_{\mathrm{U}}(x)$は$x=\tilde{x}$より前で増加関数なので，実行可能領域は$\tilde{x}$よりも左側にある
iii) どちらでも無い（すなわち$m_{\mathrm{Lmax}}\geq m_{\mathrm{Umin}}$かつ$m_{\mathrm{Lmin}}\leq m_{\mathrm{Umax}}$）ならば，$f_{\mathrm{L}}(x)-f_{\mathrm{U}}(x)$は$x=\tilde{x}$において最小値をとるので，実行可能領域は存在しない．

が言えます．

次に，（$x_{\mathrm{min}}\leq \tilde{x}\leq x_{\mathrm{max}}$を満たす）$\tilde{x}$に対し$f_{\mathrm{L}}(\tilde{x})<f_{\mathrm{U}}(\tilde{x})$であった（つまり直線$x=\tilde{x}$が実行可能領域を通過する）としましょう．$\tilde{x}$の右側か左側どちらかに最適解が存在するわけですが，

i) $m_{\mathrm{Lmax}}<0$ならば$f_{\mathrm{L}}(x)$は$x=\tilde{x}$から先で減少するので最適解は$\tilde{x}$よりも右側
ii) $m_{\mathrm{Lmin}}>0$ならば$f_{\mathrm{L}}(x)$は$x=\tilde{x}$から先で増加するので最適解は$\tilde{x}$よりも左側
iii) どちらでも無い，すなわち$m_{\mathrm{Lmin}}\leq 0\leq m_{\mathrm{Lmax}}$ならば，$\tilde{x}$が最適解

が言えます．

$f_{\mathrm{L}}(\tilde{x})-f_{\mathrm{U}}(\tilde{x})=0$であった場合は特殊な状況で，

i) $m_{\mathrm{Lmax}}<0$かつ$m_{\mathrm{Lmax}}\leq m_{\mathrm{Umin}}$ならば$\tilde{x}$は$y=f_{\mathrm{U}}(x)$と$y=f_{\mathrm{L}}(x)$が交わる左側の点の$x$座標，したがって最適解は$\tilde{x}$よりも右側
ii) $m_{\mathrm{Lmin}}>0$かつ$m_{\mathrm{Lmin}}\geq m_{\mathrm{Umax}}$ならば$\tilde{x}$は$y=f_{\mathrm{U}}(x)$と$y=f_{\mathrm{L}}(x)$が交わる右側の点の$x$座標，したがって最適解は$\tilde{x}$よりも左側
iii) どちらでも無いならば，$\tilde{x}$は$y=f_{\mathrm{U}}(x)$と$y=f_{\mathrm{L}}(x)$の接点，すなわち唯一$f_{\mathrm{L}}(\tilde{x})=f_{\mathrm{U}}(\tilde{x})$を満たす点の$x$座標であり，最適解．

が言えます．

以上のように，$x_{\mathrm{min}}\leq\tilde{x}\leq x_{\mathrm{max}}$を満たす任意の$\tilde{x}$に対し，最適解がどちらに存在するか / $\tilde{x}$自体が最適解か / 最適解が存在しないかを$O(N)$で判別できることが分かります．
後者二つのケースならば，ただちに計算を終了すれば良いです．
残り一つのケースでは，$\tilde{x}$より右側に最適解があると分かれば$x_{\mathrm{min}}$を$\tilde{x}$に，左側にあると分かれば$x_{\mathrm{max}}$を$\tilde{x}$に，それぞれ置き換えられます．

冗長な不等式制約条件の削除

前節の方法で，最適解が存在し得る$x$の範囲を絞り込みながら，同時に明らかに最適解に関係しない（冗長な）不等式制約条件を取り除いていく，というのがMegiddo-Dyerのアルゴリズムの基本的な考え方です．具体的には次のようにします．

$\mathcal{I}_{\mathrm{U}}$の要素をランダムに（頭から順番に，でも構いません）二つずつペアリングしていきます．$\mathcal{I}_{\mathrm{U}}$の要素数が奇数ならば一つ余りますが，それはそのままにしておきます．$\mathcal{I}_{\mathrm{L}}$の要素も同様にペアリングします．

ペアとなった2直線$y=m_{i}x+y_{0i}$および$y=m_{j}x+y_{0j}$が平行，すなわち$m_{i}=m_{j}$であるならば，冗長な制約条件をただちに削除することができます．たとえば$i,j\in\mathcal{I}_{\mathrm{U}}$であるペアの場合，$y_{0i}\geq y_{0j}$ならば任意の$x$に対し$y\leq m_{i}x+y_{0i}\Rightarrow y\leq m_{j}x+y_{0j}$なので，$i$が冗長な制約条件です．同様に，$i,j\in\mathcal{I}_{\mathrm{L}}$であるペアの場合，$y_{0i}\leq y_{0j}$ならば$i$が冗長な制約条件です．

$m_{i}\neq m_{j}$である場合は，2直線の交点の座標$(x_{ij},y_{ij})$が次のように求まります．

(x_{ij},y_{ij})=\left(-\frac{y_{0i}-y_{0j}}{m_{i}-m_{j}}, \frac{m_{i}y_{0j}-m_{j}y_{0i}}{m_{i}-m_{j}}\right)

たとえば，ペアリングされた直線群の交点が上図の黒点のように求まったとしましょう．運が良ければこの中に最適解がありますが，今のところそれは定かではありません．範囲を絞り込むために，$x$座標の値が中央値$x_{\mathrm{m}}$である交点をとりあえず選びます．

上図の場合，$f_{\mathrm{U}}(x_{\mathrm{m}})>f_{\mathrm{L}}(x_{\mathrm{m}})$です．$m_{\mathrm{Umax}}=m_{\mathrm{Umin}}=f'_{\mathrm{U}}(x_{\mathrm{m}})$と$m_{\mathrm{Lmax}}=m_{\mathrm{Lmin}}=f'_{\mathrm{L}}(x_{\mathrm{m}})$もすぐに求まり，$m_{\mathrm{Lmax}}<m_{\mathrm{Umin}}$であるので，前節の事柄を適用することにより，最適解が$x_{\mathrm{m}}$よりも右側にあると分かります．そこで，$x_{\mathrm{min}}\leftarrow x_{\mathrm{m}}$と更新します．

さて，$i,j\in\mathcal{I}_{\mathrm{U}}$であるペアの交点$(x_{ij},y_{ij})$について，$x_{ij}<x_{\mathrm{min}}$である場合，最適解は右にありますから$x=x_{i,j}$において$f_{\mathrm{L}}(x)-f_{\mathrm{U}}(x)$は減少関数であるはずです．したがって$m_{i}>m_{j}$ならば$i$が，$m_{i}<m_{j}$ならば$j$が冗長な制約条件です．
$i,j\in\mathcal{I}_{\mathrm{L}}$であるペアの交点$(x_{ij},y_{ij})$についてはこれが逆転して，$m_{i}>m_{j}$ならば$j$が，$m_{i}<m_{j}$ならば$i$が冗長な制約条件です．

$x_{\mathrm{max}}\leftarrow x_{\mathrm{m}}$と更新されたときも全く同じ考え方に基づけば，$i,j\in\mathcal{I}_{\mathrm{U}}$であるペアの交点$(x_{ij},y_{ij})$について，$x_{ij}>x_{\mathrm{max}}$である場合，最適解は左にありますから，$x=x_{i,j}$において$f_{\mathrm{L}}(x)-f_{\mathrm{U}}(x)$は増加関数であるはずです．したがって$m_{i}>m_{j}$ならば$j$が，$m_{i}<m_{j}$ならば$i$が冗長な制約条件です．
$i,j\in\mathcal{I}_{\mathrm{L}}$であるペアの交点$(x_{ij},y_{ij})$については，$m_{i}>m_{j}$ならば$i$が，$m_{i}<m_{j}$ならば$j$が冗長な制約条件です．

冗長な制約条件を削除し，残った制約条件を再度ペアリングして，同じことを反復します．

こうして冗長な制約条件を取り除きながら，途中で最適解（$m_{\mathrm{Lmin}}\leq 0\leq m_{\mathrm{Lmax}}$となる交点）が見つかればそこで終了し，見つからないまま制約条件の数が十分少なくなったら，前々節で述べた方法で最適解を直接求めます．

アルゴリズム

以上の内容をアルゴリズムとしてまとめます．

まずは全体のフローです．

1. 評価関数の係数$c_{1}$，$c_{2}$，制約条件（不等式の係数）群${(a_{k1},a_{k2},b_{k})}$（$k=1,\cdots,N$）が与えられる．
2. 初期化処理．$(x_{1},x_{2})$を$(x,y)$に座標変換したときの$\mathcal{D}_{\mathrm{U}}$，$\mathcal{D}_{\mathrm{L}}$，$x_{\mathrm{max}}$，$x_{\mathrm{min}}$を得る．不適切な条件設定がなされていたら終了．
3. $|\mathcal{D}_{\mathrm{U}}|+|\mathcal{D}_{\mathrm{L}}|\leq 4$となるまで次を繰り返す．
4. $x_{\mathrm{min}}$，$x_{\mathrm{max}}$を更新．最適解が存在しないと分かったら終了．最適解$(x^{*},y^{*})$が求まったら5へ．
5. 冗長な制約条件を削除．
6. $x_{\mathrm{min}}$，$x_{\mathrm{max}}$，$\mathcal{D}_{\mathrm{U}}$および$\mathcal{D}_{\mathrm{L}}$から直接最適解$(x^{*},y^{*})$を求める．
7. $(x^{*},y^{*})$を$(x_{1}^{*},x_{2}^{*})$に変換．
8. 終了．

「初期化処理」は，与えられた$(a_{1k},a_{2k},b_{k})$（$k=1,\cdots,N$）を$(m_{k},y_{k0})$に変換し，$\mathcal{D}_{\mathrm{U}}$，$\mathcal{D}_{\mathrm{L}}$に振り分ける処理です．$c_{1}$，$c_{2}$も含め，不適切な値が与えられていないかのチェックもここで行います．

1. $c_{2}=0$なら終了（$c_{1}=0$なら異常ケース，$c_{1}\neq 0$なら$x_{1}$と$x_{2}$の入れ替えを提案すべきケース）．
2. $x_{\mathrm{min}}\leftarrow -\infty$，$x_{\mathrm{max}}\leftarrow \infty$，$\mathcal{D}_{\mathrm{U}}\leftarrow\emptyset$，$\mathcal{D}_{\mathrm{L}}\leftarrow\emptyset$とする．
3. $k=1,\cdots,N$について
   1. $a_{k2}=0$なら
      1. $a_{k1}=0$ならスキップ（無意味な制約条件）．
      2. $\bar{x}_{k}\leftarrow b_{k}/a_{k1}$
      3. $a_{k1}<0$かつ$\bar{x}_{k}>x_{\mathrm{min}}$なら$x_{\mathrm{min}}\leftarrow \bar{x}_{k}$
      4. $a_{k1}>0$かつ$\bar{x}_{k}<x_{\mathrm{max}}$なら$x_{\mathrm{max}}\leftarrow \bar{x}_{k}$
   2. $m_{k}\leftarrow c_{1}-a_{k1}c_{2}/a_{k2}$
   3. $y_{0k}\leftarrow b_{k}c_{2}/a_{k2}$
   4. $c_{2}/a_{k2}>0$なら$\mathcal{D}_{\mathrm{U}}\leftarrow\mathcal{D}_{\mathrm{U}}\cup\left\{(m_{k},y_{0k})\right\}$
   5. $c_{2}/a_{k2}<0$なら$\mathcal{D}_{\mathrm{L}}\leftarrow\mathcal{D}_{\mathrm{L}}\cup\left\{(m_{k},y_{0k})\right\}$

「$x_{\mathrm{min}}$，$x_{\mathrm{max}}$の更新」は，上述した内容そのままです．

1. $\mathcal{D}_{\mathrm{U}}$要素ペアリスト$\mathcal{P}_{\mathrm{U}}$および$\mathcal{D}_{\mathrm{L}}$要素ペアリスト$\mathcal{P}_{\mathrm{L}}$を作成．
2. ペアリスト交点$x$座標の中央値$x_{\mathrm{m}}$を求める．
3. $x_{\mathrm{m}}$に対する$f_{\mathrm{L}}(x_{\mathrm{m}})$の値$y_{\mathrm{L}}$，傾き最小値$m_{\mathrm{Lmin}}$，傾き最大値$m_{\mathrm{Lmax}}$を求める．
4. $x_{\mathrm{m}}$に対する$f_{\mathrm{U}}(x_{\mathrm{m}})$の値$y_{\mathrm{U}}$，傾き最小値$m_{\mathrm{Umin}}$，傾き最大値$m_{\mathrm{Umax}}$を求める．
5. $y_{\mathrm{L}}>y_{\mathrm{U}}$かつ
   1. $m_{\mathrm{Lmax}}<m_{\mathrm{Umin}}$ならば，$x_{\mathrm{min}}\leftarrow x_{\mathrm{m}}$
   2. $m_{\mathrm{Lmin}}>m_{\mathrm{Umax}}$ならば，$x_{\mathrm{max}}\leftarrow x_{\mathrm{m}}$
   3. どちらでもない（$m_{\mathrm{Lmax}}\geq m_{\mathrm{Umin}}$かつ$m_{\mathrm{Lmin}}\leq m_{\mathrm{Umax}}$）ならば解なし，終了．
6. $y_{\mathrm{L}}<y_{\mathrm{U}}$かつ
   1. $m_{\mathrm{Lmax}}<0$ならば$x_{\mathrm{min}}\leftarrow x_{\mathrm{m}}$
   2. $m_{\mathrm{Lmin}}>0$ならば$x_{\mathrm{max}}\leftarrow x_{\mathrm{m}}$
   3. どちらでもないならば，$(x^{*},y^{*})\leftarrow (x_{\mathrm{m}},y_{\mathrm{L}})$（$x_{\mathrm{m}}$が最適解）として終了．
7. $y_{\mathrm{L}}=y_{\mathrm{U}}$かつ
   1. $m_{\mathrm{Lmax}}<0$かつ$m_{\mathrm{Lmax}}\leq m_{\mathrm{Umin}}$ならば，$x_{\mathrm{min}}\leftarrow x_{\mathrm{m}}$
   2. $m_{\mathrm{Lmin}}>0$かつ$m_{\mathrm{Lmin}}\geq m_{\mathrm{Umax}}$ならば，$x_{\mathrm{max}}\leftarrow x_{\mathrm{m}}$
   3. どちらでもないならば，$(x^{*},y^{*})\leftarrow (x_{\mathrm{m}},y_{\mathrm{L}})$（$x_{\mathrm{m}}$が最適解）として終了．

「$\mathcal{D}_{\mathrm{U}}$要素ペアリスト$\mathcal{P}_{\mathrm{U}}$作成」は，次のようにします．

1. $\mathcal{P}_{\mathrm{U}}\leftarrow\emptyset$，$i\leftarrow 1$とする．
2. $j\leftarrow i+1$
3. $j>|\mathcal{D}_{\mathrm{U}}|$なら終了
4. $m_{i}=m_{j}$なら，
   1. $y_{0i}\geq y_{0j}$なら$\mathcal{D}_{\mathrm{U}}\leftarrow\mathcal{D}_{\mathrm{U}}-\left\{(m_{i},y_{0i})\right\}$，$i\leftarrow j$．2に戻る．
   2. （$y_{0i}<y_{0j}$なので）$\mathcal{D}_{\mathrm{U}}\leftarrow\mathcal{D}_{\mathrm{U}}-\left\{(m_{j},y_{0j})\right\}$，$j\leftarrow j+1$．3に戻る．
5. （$m_{i}\neq m_{j}$なので）$\mathcal{P}_{\mathrm{U}}\leftarrow\mathcal{P}_{\mathrm{U}}\cup\left\{((m_{i},y_{0i}), (m_{j},y_{0j}))\right\}$
6. $i\leftarrow j+1$
7. 2に戻る．

「$\mathcal{D}_{\mathrm{U}}$要素ペアリスト$\mathcal{P}_{\mathrm{L}}$作成」も同様の手続きとなります．4(a)(b)で冗長な制約条件を判別するための$y_{0i}$と$y_{0j}$の比較演算子が入れ替わることに注意して下さい．

そして「冗長な不等式制約条件の削除」ですが，$\mathcal{D}_{\mathrm{U}}$に限れば次のアルゴリズムで行えます．

1. $\mathcal{P}_{\mathrm{U}}$の全ての要素$((m_{i},y_{0i}), (m_{j},y_{0j}))$について
   1. $x_{i,j}<x_{\mathrm{min}}$かつ
      1. $m_{i}<m_{j}$ならば$\mathcal{D}_{\mathrm{U}}\leftarrow\mathcal{D}_{\mathrm{U}}-\left\{(m_{j},y_{0j})\right\}$
      2. $m_{i}>m_{j}$ならば$\mathcal{D}_{\mathrm{U}}\leftarrow\mathcal{D}_{\mathrm{U}}-\left\{(m_{i},y_{0i})\right\}$
   2. $x_{i,j}>x_{\mathrm{max}}$かつ
      1. $m_{i}<m_{j}$ならば$\mathcal{D}_{\mathrm{U}}\leftarrow\mathcal{D}_{\mathrm{U}}-\left\{(m_{i},y_{0i})\right\}$
      2. $m_{i}>m_{j}$ならば$\mathcal{D}_{\mathrm{U}}\leftarrow\mathcal{D}_{\mathrm{U}}-\left\{(m_{j},y_{0j})\right\}$

$\mathcal{D}_{\mathrm{L}}$にも同様のアルゴリズムを適用すれば良いですが，前節で述べた通り，これも判別の比較演算子が入れ替わります．

最後に，$x_{\mathrm{min}}$，$x_{\mathrm{max}}$，$\bar{\mathcal{D}}_{\mathrm{U}}$および$\bar{\mathcal{D}}_{\mathrm{L}}$から直接最適解$(x^{*},y^{*})$を求めるアルゴリズムです．

1. $y^{*}\leftarrow\infty$とする．
2. $\mathcal{D}_{\mathrm{U}}$の要素$(m_{k},y_{0k})$（$k=1,\cdots,|\mathcal{D}_{\mathrm{U}}|$）について
   1. $\mathcal{D}_{\mathrm{L}}$の要素$(m_{l},y_{0l})$（$l=1,\cdots,|\mathcal{D}_{\mathrm{L}}|$）について
      1. 交点$(x_{k,l},y_{k,l})$を求める．
      2. $(x_{k,l},y_{k,l})$が実行可能かつ$y_{k,l}<y^{*}$ならば$(x^{*},y^{*})\leftarrow (x_{k,l},y_{k,l})$
3. $\mathcal{D}_{\mathrm{L}}$の要素$(m_{k},y_{0k})$（$k=1,\cdots,|\mathcal{D}_{\mathrm{L}}|-1$）について
   1. $\mathcal{D}_{\mathrm{L}}$の要素$(m_{l},y_{0l})$（$l=k+1,\cdots,|\mathcal{D}_{\mathrm{L}}|$）について
      1. 交点$(x_{k,l},y_{k,l})$を求める．
      2. $(x_{k,l},y_{k,l})$が実行可能かつ$y_{k,l}<y^{*}$ならば$(x^{*},y^{*})\leftarrow (x_{k,l},y_{k,l})$
4. $x_{\mathrm{min}}$のときの値$y_{l}$を求める．$(x_{\mathrm{min}},y_{l})$が実行可能かつ$y_{l}<y^{*}$ならば$(x^{*},y^{*})\leftarrow (x_{\mathrm{min}},y_{l})$
5. $x_{\mathrm{max}}$のときの値$y_{l}$を求める．$(x_{\mathrm{max}},y_{l})$が実行可能かつ$y_{l}<y^{*}$ならば$(x^{*},y^{*})\leftarrow (x_{\mathrm{max}},y_{l})$

このアルゴリズムの計算量を$c(N)$とおきましょう．1回の反復計算で，全ペアの交点のうち$x$座標中央値を境に最適解と反対側にあるもの（全体の2分の1）の，交点を作る制約条件2本のうち1本が取り除かれるので，制約条件の数はつごう4分の3に減ることになります．$x$座標中央値を求めるのも，最適解と反対側にある交点を見つけるのも，計算量は$N$に比例すると考えて良いので，$c(N)$は次式を満たします．

c(N)\leq \alpha N+c(3N/4)

ただし，$\alpha$は定数です．このことより，$c(N)=O(N)$であると分かります．

実装上の注意

元論文に書かれている内容は以上なのですが，書かれていないこととして，実装する上で工夫できることや注意すべきことをまとめておきます．

コードの重複をなるべく減らす工夫

Megiddo-Dyerのアルゴリズムには，「処理としてはほとんど同じなのに，対象が$\mathcal{D}_{\mathrm{U}}$か$\mathcal{D}_{\mathrm{L}}$かによって，ある値の大小関係に基づいて条件分岐する際に不等号の向きが反対になるもの」が幾つかあります．

たとえば，与えられた$x$に対し$f_{\mathrm{U}}(x)$と$f_{\mathrm{L}}(x)$を求める処理ですが，前者は$\mathcal{D}_{\mathrm{U}}$に含まれる制約条件の最小値を求めるのに対し，後者は$\mathcal{D}_{\mathrm{L}}$に含まれる制約条件の最大値を求めます．

上の例の他には，$x$に対する$f_{\mathrm{U}}(x)$と$f_{\mathrm{L}}(x)$の勾配の最小値と最大値を求める処理，平行な境界直線を持つ制約条件の一方を削除する処理，非平行な境界直線のペアのうち冗長なものを判別する処理があります．これらを，できるだけ共通部品を用いて実装しようとするのは自然な発想でしょう．

そこで，二つの値$a$，$b$に対し真偽値$a>b$を返す関数と$a<b$を返す関数を用意し，これらの関数自体を引数にとるような一段階抽象化された関数を用意することで，コード重複を最小限に抑えることが可能になります．このような手法は多くのプログラミング言語で使えます．

制約条件ペアリストの交点計算と中央値の探索

「$x_{\mathrm{min}}$，$x_{\mathrm{max}}$の更新」の説明は，制約条件ペアリストを作成してからそれぞれのペアの交点を計算する，というように読めますが，交点はペアが平行かどうかの判別に続けて求めてしまうのが効率良いですし，他の処理でも必要になるので，ペア情報$\left\{((m_{i},y_{0i}), (m_{j},y_{0j}))\right\}$と一緒に保存しておくのが良いです．

論文ではさらに，「ペアリスト交点$x$座標の中央値$x_{\mathrm{m}}$を求める」こととなっています．ペアリストは$\mathcal{P}_{\mathrm{U}}$と$\mathcal{P}_{\mathrm{L}}$に分かれていますので，$x_{\mathrm{m}}$を求めるためにはこれらをマージしなければなりません．
一方，「冗長な不等式制約条件の削除」においては$\mathcal{P}_{\mathrm{U}}$と$\mathcal{P}_{\mathrm{L}}$を区別して処理しなければならないので，マージしたリストを再度元の二つに分離する必要があります．
この手続きはちょっと煩雑です．
加えて言えば，$\mathcal{P}_{\mathrm{L}}$のペア交点は最適解となる可能性がありますが，$\mathcal{P}_{\mathrm{U}}$のペア交点は決して最適解になりません．

そこで，基本的に$x$座標中央値を探索する対象は$\mathcal{P}_{\mathrm{L}}$の交点群のみとし，制約の削除が進んで$\mathcal{P}_{\mathrm{L}}=\emptyset$となった場合に探索対象を$\mathcal{P}_{\mathrm{U}}$の交点群に切り替える，という方法にアレンジします．これならば，$\mathcal{P}_{\mathrm{U}}$と$\mathcal{P}_{\mathrm{L}}$をマージ・再分離する手間はなくなりますし，選ばれた中央値が最適解となる可能性も上がります．

制約条件ペアリストの効率的更新

「$\mathcal{D}_{\mathrm{U}}$要素ペアリスト$\mathcal{P}_{\mathrm{U}}$作成」では，削除されず残った全ての制約条件群から毎回ペアリングしています．これをもう少し効率良くするために，「ペアリングされていない制約条件群」の集合を作ることにします．
最初は全ての制約条件をこれに登録します．また，ある（冗長な）制約条件が削除されたら，もう一方の制約条件をこの集合に戻してやります．ペアリング処理は，この集合から要素を二つずつ，残りの数が1以下になるまで取り出して行うようにすれば良いでしょう．

無限解の存在判別

実行可能領域が下側境界を持たないとき，線形計画問題は無限解を持つのでした．一方，下側境界があったとしても無限解が存在するケースがあります．元論文で示されたアルゴリズムにはこれらのケースを判別する手続きがありませんので，本記事で独自に設けることにします．

下側境界が存在し，かつそれが下限を持つならば，無限解が存在しないと言えます．したがって，次の条件のうち少なくとも一つが満たされるならば，$x=-\infty$は無限解ではありません．

• $x_{\mathrm{min}}>-\infty$である．
• $f_{\mathrm{L}}(x)$の傾き最小値を$m_{\mathrm{Lmin}}$としたとき，$m_{\mathrm{Lmin}}<0$である．
• 同じく$f_{\mathrm{U}}(x)$の傾き最大値を$m_{\mathrm{Umax}}$としたとき，$m_{\mathrm{Umax}}>m_{\mathrm{Lmin}}$である．

同様に，少なくとも次のうち一つが満たされるならば，$x=\infty$は無限解ではありません．

• $x_{\mathrm{max}}<\infty$である．
• $f_{\mathrm{L}}(x)$の傾き最大値を$m_{\mathrm{Lmax}}$としたとき，$m_{\mathrm{Lmax}}>0$である．
• 同じく$f_{\mathrm{U}}(x)$の傾き最小値を$m_{\mathrm{Umin}}$としたとき，$m_{\mathrm{Umin}}<m_{\mathrm{Lmax}}$である．

よって，初期化処理において制約条件を登録する際に上記をチェックし，最後までどの条件も満たされなければ無限解が存在する，と判断すれば良いです．

退化

実行可能領域に$f_{\mathrm{L}}(x)$の傾きが$0$となる区間があるならば，その区間に無数の最適解が存在することになります．このような状況を，「線形計画問題が退化している」と言います．

退化していても，その区間が右か左に他の制約条件境界と交点を持つならば，アルゴリズムの性質上その交点が最適解として選ばれるので，大きな問題はありません．その区間が$x=-\infty$から$x=\infty$まで広がっている，すなわち問題がただひとつの制約条件$y\geq y_{0}$（$y_{0}$は任意の定数）しか持たないならば，アルゴリズムは正常に動作しません．
このようなケースも，最適化計算前に除外すべきです（判別方法は自明です）．

実行可能領域の左端と右端の整合性チェック

「上側境界と下側境界が交わらないため実行可能領域が存在しない」ケースは，「$x_{\mathrm{min}}$，$x_{\mathrm{max}}$の更新」の中で検出されます．一方，「左側境界と右側境界が交わらないため実行可能領域が存在しない」ケースについては言及されていません．反復計算が開始すればこのようなことは起こり得ないので，これの検出は初期化処理の中で行えば良いです．
つまり，ある$k$において$a_{k2}=0$ならば$x_{\mathrm{min}}$または$x_{\mathrm{max}}$が更新されるわけですが，その結果$x_{\mathrm{max}}<x_{\mathrm{min}}$となったら解なし，と判定できます．

ZMでの実装

前節で述べた注意点も全て反映したMegiddo-Dyerのアルゴリズムを，ZMに実装しています．

ヘッダはzm_opt_lp_megiddodyer.h，ソースはzm_opt_lp_megiddodyer.cをご参照下さい．