線形計画法
∗
半場 滋
†1
目的
本実験の目的は, プログラムを作成することを通じ て, システム工学の重要な一分野である線形計画法に 関する理解を深めることである。2
原理
2.1 線形計画問題
線形計画問題とは, 与えられた線形な等式および不 等式制約のもとで, 線形目的関数を最大化あるいは最 小化する問題である。まず例を挙げる。 例 1 パソコンショップの店主である A さんは, ある顧 客から「PC/AT 互換機のシステムを組んで欲しい」と いう依頼を受けた。顧客の要望は, メモリ 100MB 以上, ディスク 5GB 以上, 予 算はディスクとメモリを合わせて 10 万円以 下, 予算の範囲でメモリもディスクも可能な 限りたくさん載せたい, というものであった。 話を簡単にするため, (非現実的な仮定だが) メモリ およびディスクの大きさは任意の正の実数値を取るも のとする。メモリの価格は 1MB あたり 100 円, ディス クの価格は 1GB あたり 2,500 円とする。また, 顧客の 依頼したパソコンは, メモリは最大 800MB まで, ディ スクは無限に増設できるものとする。A さんの店では, メモリは 1MB あたり 10 円, ディスクは 1GB あたり 200円の利益があるものとする。A さんとしては, 顧 客の予算を目一杯使って, なるべく多くの利益をあげ たい。さて, どのような意思決定をするのが良いだろ うか? メモリの大きさをx MB, ディスクの大きさを y GB としよう。すると, A さんの利益 (目的関数) は, 10x + 200y (1) ∗この文書は琉球大学工学部電気電子工学科で開講されている 4 年次学生向けの専門実験 (電力工学実験, 電子工学実験, システム工 学実験) 向けに執筆されたテキストに若干手を加えたものである。 内容を改変しない限り自由に再配布して頂いて構わない。 †琉球大学工学部電気電子工学科 ([email protected]) となる。よって, A さんが解くべき問題は, 「10x + 200y(目的関数) を最大化すること」となる。 仮にメモリの大きさx とディスクの大きさ y を勝手 に設定できるなら, A さんはいくらでもお金を儲けら れることになる。だが実際にはいくつか制約条件があ る。それは, • メモリの大きさは正 • メモリの大きさは 100 MB 以上 (顧客の要望) • メモリの大きさは 800 MB 以下 (計算機の仕様) • ディスクの大きさは正 • ディスクの大きさは 5 GB 以上 (顧客の要望) • メモリとディスクの合計金額は 10 万円以下 (顧客 の懐具合) というものである。これらを書き直すと, x ≥ 0 x ≥ 100 x ≤ 800 y ≥ 0 y ≥ 5 100x + 2500y ≤ 100000 (2) となる。 図 1 は, 制約条件を満たす領域の形と目的関数のよ うすをあらわしたものである。制約条件を満たす領域 は図中の網のかかった台形の内部である。目的関数は 図中に引かれた傾き−1 の一点鎖線上で一定値を取る。 この直線が図中の右上に移動するほど目的関数の値は 大きくなる。この図から, 台形の右上の頂点の部分で 目的関数が最大値を取ることがわかる。 例 1 は, 結局 変数 (x, y), 目的関数 10x + 200y および制約 条件 (2) が与えられたとき, 与えられた制約 条件のもとで目的関数を最大化するようなx, y を求めよx y 0 400 800 メモリ (MB) HDD(GB) 20 40 x = 100 x = 800 y = 5 Aさんの利益 6,000円 Aさんの利益 4,000円 Aさんの利益 2,000円 最適解 制約条件を満たす領域 予算に関する制約: 100x + 2500y ≤ 100000 図 1 制約条件を満たす領域と目的関数の挙動 という問題であった。これを一般化したのが線形計画 問題である。 線形計画問題 (一般形) n 個の変数 x1, . . . , xnに対 し,m 個の等式あるいは不等式の制約条件 a11x1+. . . + a1nxn ≥ b1 . . . . ak1x1+. . . + aknxn = bk . . . . al1x1+. . . + alnxn ≤ bl . . . . (3) が与えられ, かついくつかの変数に符号に関する制約 条件xi1≥ 0, . . . , xir ≥ 0, xj1 ≤ 0, . . ., xjs ≤ 0 が課さ れているものとする。このとき, 制約条件 (3) を満たす 範囲内で, 目的関数 c1x1+. . . + cnxn (4) を最小化 (あるいは最大化) せよ。 ここに, 最小化 (最大化) とは, 目的関数を最小 (最大) にするx1からxn の値とその最小値 (最大値) を求め ることを言う。 線形計画問題の例として, 工場において一定の原料 やエネルギー消費量の制約のもとで利益を最大化する ための生産計画を立てる問題などが挙げられる。実用 的な線形計画問題では, 変数および制約条件の数は極 めて大きくなる。 線形計画問題の解法としてよく知られているものに, シンプレックス法と内点法がある。シンプレックス法 は 1947 年に G. B. Dantzig によって提案された線形計 画問題の解法である。この方法は, 幾何学的には, 凸多 面体の形をした制約領域の隣接する頂点をたどりなが ら解を改善する方法である。シンプレックス法は, 提 案されて以降 40 年ほどの間, 最も効率的な線形計画 問題の解法とされてきた。これに対し, 1984 年に N. Karmarkarによって提案された内点法は, 理論的な収 束特性がシンプレックス法より良いばかりでなく, 極 めて大規模な問題に対しては実用上の計算効率もシン プレックス法に勝るため, 提案されて以来, 爆発的に研 究が進んできた。(ただし, 規模の小さい問題ではシン プレックス法の方が圧倒的に速い)。 本実験ではシンプレックス法と内点法の双方を取り 扱う。 以下ではシンプレックス法と内点法の基本について 説明してゆくのだが, 説明がかなり長いので, 読者の負 担を減らすために, 読まなくても実験に支障がない部 分には印がつけてある。記号 が欄の左端に記載され ているところの記号 が欄の右端に記載されていると ころあいだは読み飛ばして構わない。具体的には, ペー ジ 5 からページ 7 までとページ 10 からページ 14 の該 当部分である。
2.2 標準形と基準形
線形計画問題の一般形は, 異なる向きの不等号や制 約条件が付いた変数と制約条件のない変数が混在して いるために, 数学的に取り扱いにくい。そこで, ふつう は, 線形計画問題の一般形を特定の標準的な形に変換 してから解く。 標準的な形としてよく用いられるのは, 標準形と呼 ばれるものと, 基準形と呼ばれるものである。以下で これらについて説明する。 以下では, 特に断らない限り,A は m 行 n 列の実行 列,b は m 次元のベクトル, c, x は n 次元のベクトル とする。また, A が m 行 n 列の実行列であることを A ∈ Rm×nという記法であらわし,x が n 次元の実ベ クトルであるということをx ∈ Rnという記法であら わすi)。さらに, 行列A の第 i 番目の列ベクトルを ai などと書き, ベクトルb, c などの第 i 要素をそれぞれ bi, ciなどと書く。また, x と w がともに n 次のベク トルであるとき, x ≥ w という記法によって, ベクト ルx の各成分がベクトル w の対応する成分以上であ る, すなわちxi ≥ wi (i = 1, . . ., n) となっているとい うことをあらわす。零ベクトルあるいは零行列をあら わすのには記号 0 を使う。次元を明示する必要がある i)記号Êは実数体をあらわす。ときには 0n, 0m×nなどのように書くii)。 線形計画問題の標準形とは, 以下のような形式の問 題である。 標準形: min{z = cTx : Ax = b, x ≥ 0} この式の読み方は, 「制約条件Ax = b, x ≥ 0 のも とで目的関数z = cTx を最小化せよ」である。なお, minは minimize(最小化する) の略である。 続いて, 線形計画問題の一般形を標準形に変換する 手順について説明する。このためには, 以下に述べる ような一連の手順を実行すればよい。 1. 不等式制約ai,1x1+· · · + ai,nxn ≥ bi: 新たな変 数yi ∈ R を導入し, 制約条件を ai,1x1+· · · + ai,nxn− yi ≥ bi,yi ≥ 0 と書き直す。目的関数は 変更しない。 2. 等式制約ai,1x1+· · · + ai,nxn=bi: 変更しない。 3. 不等式制約ai,1x1+· · · + ai,nxn ≤ bi: 新たな変 数yi ∈ R を導入し, 制約条件を ai,1x1+· · · + ai,nxn+yi ≥ bi,yi ≥ 0 と書き直す。目的関数は 変更しない。 4. 符号制約なしの変数xi: 新たに 2 個の正の変数 xi−, xi+ を導入し, xi = xi+ − xi−, xi+ ≥ 0, xi− ≥ 0 と書き直す。このとき, 各制約条件にあ
らわれる項ap,ixiはap,ixi+− ap,ixi−で置き換え られる。また, 目的関数にあらわれるcixiという 項もcixi+− cixi−で置き換えられる。 5. 負でない変数xi≥ 0: 変更しない。 6. 正でない変数xi ≤ 0: xi =−xi と変数を置き直 す。このとき, 制約条件にあらわれる項ap,ixiは −ap,ixiで置き換えられる。また, 目的関数にあら われるcixiという項も−cixiで置き換えられる。 7. cTx を最大化する問題 (最大化問題): c =−c と 定義して, 新たな目的関数 (c)Tx を最小化する。 上記の手続きのうち, 1, 3, 4 では変数の数がそれぞれ 1個ずつ増えることを注意しておく。 これに対し, 線形計画問題の基準形とは, 次のような 問題である。 基準形: min{z = cTx : Ax ≥ b, x ≥ 0} ii)零ベクトルと零行列に同じ記号が用いられるので, 混乱しない よう注意すること。 標準形の場合と同様の手順によって, 線形計画問題の 一般形を基準形に変換することができる。この場合は, 制約条件や変数の数は増加する。基準形への変換が標 準形への変換と異なる点は, 上記とは異なる方法を使 うと, 等式制約条件や制約条件なしの変数を除去する 際に, 制約条件や変数の数を減らすこともできる, とい うことである。この方法はやや繁雑なのでここでは述 べない。興味のある読者は文献 [5] を参照すること。
2.3 シンプレックス法
本節では, 線形計画問題を解くためのシンプレック ス法というアルゴリズムについて簡単に説明する。こ の部分の説明はおもに文献 [3],[4] による。 標準形で記述された線形計画問題: (STD) minz = cTx : Ax = b, x ≥ 0 (5) が与えられているものとする。先に注意したように, A ∈ Rm×n,b ∈ Rm,c ∈ Rnである。また,n ≥ m で, 行列A はフルランクであると仮定しておく。また, 説 明の便宜上, (STD) の制約条件を満たす領域は空でな いと仮定する。 シンプレックス法は, まず制約条件を満たす (最適と は限らない) 基底解と呼ばれる解をひとつ見付け, 続い て最適解が得られるまで基底解を改善していく方法で ある。標準形の問題が一定の条件を満たすときは, シ ンプレックス法を用いると有限回の計算で最適解が得 られることが保証される。 ではシンプレックス法の説明に入ろう。 行列A の列ベクトル a1, . . . , anの中から適当なm 個の線形独立なベクトルを選び出し, これらを並べて B = [ai1, . . . , aim] (6) という行列を作る。また, 行列A の列ベクトルから行 列B の列ベクトルを除いた残りを並べて, N = [aj1, . . . , ajn−m] (7) という行列を作る。行列B, N はそれぞれ基底行列, 非基底行列と呼ばれる。さらに, 行列B に対応する変 数の組xi1, . . . , ximを集めて xB = [xi1, . . . , xim]T (8) なるベクトルを作る。xBを基底変数ベクトルと呼ぶ。 続いて, 行列N に対応する変数の組 xj1, . . . , xjn−mを 集めて xN = [xj1, . . . , xjn−m]T (9)なるベクトルを作る。xNを非基底変数ベクトルと呼ぶ。 ここで, 目的関数に対応するベクトルc を基底行列 に対応した形に分割し, cB= [ci1, . . . , cim], cN = [cj1, . . . , cjn−m] (10) とおくと, 目的関数z = cTx は z = cT BxB+cTNxN (11) と書き直される。 さて, 変数ベクトルx を xB =B−1b, xN = 0 (12) となるように選んでみよう。等式制約条件Ax = b は BxB+NxN =b (13) と書き直せるから, (12) から決まるx は制約条件 Ax = b を満たしている。このようにして (12) から定められ た解x を基底解と呼ぶ。ただし, この段階では, 符号制 約x ≥ 0 が満たされているかどうかは明らかでない。 ところで, B−1b ≥ 0 という不等式が成り立つ場合 には, xB ≥ 0, xN ≥ 0 より, 基底解 x は符号制約 x ≥ 0 を満たす。よって, この場合には, 基底解 x は 問題 (STD) のひとつの解になっていることがわかる。 (12)から定まる基底解x が問題 (STD) の解となると き, この解を実行可能基底解と呼び, 対応する行列B を実行可能基底行列と呼ぶiii)。 以下では, まず実行可能基底行列がひとつ得られて いる場合に最適解を構成するためのアルゴリズムにつ いて説明し, 続いて実行可能基底行列が未知の場合へ の拡張について述べる。 2.3.1 実行可能基底行列が得られている場合 実行可能基底行列B がすでに得られ, 制約条件 Ax = b が (13) の形に書き直されているものとする。ここで, ¯ b = B−1b, A¯N =B−1N (14) と定義し, (13) を xB= ¯b − ¯ANxN (15) のように書き直す。(15) を (11) に代入し, zB=cTBb, ¯c¯ TN =cTN − cTBA¯N (16) iii)ここで言う「解」とは制約条件をすべて満たしているもののこ とである。与えられた最適化問題の解 (制約条件をすべて満たし, か つ目的関数を最小化するもの) とは異なるので, 混乱しないよう注意 すること。 とおくと, z = zB+ ¯cTNxN (17) なる表現が得られる。 ここで, 行列 ¯AN の第i 列ベクトルを ¯aiと書くこと にする。すなわち, ¯ AN = [¯a1, . . . , ¯an−m] (18) と書く。(15) と (17) をまとめると, 問題 (SP) と等価な (BF) min{z = zB+ ¯cTNxN : xB = ¯b − ¯ANxN, xB ≥ 0, xN ≥ 0} (19) という問題が得られる。問題 (BF) を実行可能基底行 列B に対する基底形式表現と呼ぶ。 さて, 基底解 (12) は (BF) の制約条件を満たしてい るので, 最適とは限らない解である。そこで, 基底解 (12)を変更して目的関数z をより小さくすることを考 える。このために, (15) と (17) を見比べると, 以下の ことがわかる。 • (16) のベクトル ¯cN の要素の中から値が最小のも のを探す。ここで, 添字s を持つ要素 ¯cN,sが最小 値であったものとしよう。¯cN,s ≥ 0 であれば, 非 基底変数をどのように変更しても目的関数は減少 しないから, 現在のz がすでに最小値である。逆 に ¯cN,s < 0 のときには, (17) で ¯cN,sに対応する 非基底変数xjsの値を零から増加させればz は減 少するから, 現在の解を改善できる可能性がある。 • 非基底変数 xjsは (19) の制約条件を崩さない範囲 で選ばれなければならない。xjs の値を零からθ に変更したとき, xB = ¯b − ¯asθ となる。ベクト ル ¯b の各成分を ¯b1, . . . , ¯bm,ベクトル ¯asの各成分 を ¯as;1, . . ., ¯as;mと書く。xB≥ 0 という制約から, ¯
as,iが正のときには, θ を ¯bi/¯as;i以下に取る必要
がある。これに対し, ¯as,i = 0のときは,θ をいく ら大きく取っても制約条件はつねに満たされる。 よって, 添字集合{1, . . ., m} から ¯as;i> 0 となる ものを抜き出し,
¯
θ = min{¯bi/¯as;i: ¯as;i> 0, 1 ≤ i ≤ m} とすると, この ¯θ が θ の取りうる上限となる。一 方,as<= 0 のときには, θ をいくらでも大きく取
ることができるため, z の値をいくらでも小さく
以上のような考え方からシンプレックス法のアルゴリ ズムが得られるiv)。 シンプレックス法 (初期化) 実行可能基底行列B が与えられているもの とする。B を使って基底形式表現 (19) を求める。 (ステップ 1) {¯cN,1, . . . , ¯cN,n−m} の中で最小のものを ¯ cN,sとする。条件を満たすものが複数あるときに は, どれを選んでもよい。¯cN,s≥ 0 であれば終了。 この場合には, x∗ B=B−1b, x∗N = 0 (20) が最適解である。¯cN,s< 0 であればステップ 2 に 進む。 (ステップ 2) 非基底行列N の第 s 列 ajsを取り, 線 形方程式 B¯as=ajs (21) を解いてベクトル ¯asを求める。¯as ≤ 0 なら終 了 (無限解が存在する)。そうでないときは, ¯θ = min{¯bi/¯as;i : ¯as;i > 0, 1 ≤ i ≤ m} とする。ある
r に対して ¯θ = ¯br/¯as;rとなっているので, このよ うな添字r を選ぶ (条件を満たす添字が複数ある 場合は, どれを選んでもよい)。 (ステップ 3) ステップ 1 とステップ 2 で得られた添字 r と s を用いて基底行列を B = [ai1, . . ., air−1, air, air+1, . . . , aim] (22) から B= [a i1, . . . , air−1, ajs, air+1, . . ., aim] (23) に変更し, 新しい基底変数と非基底変数から対応 する ¯cNを求めてステップ 1 に戻る。 なお, (14) の第 1 式の ¯b を求めるときには, 逆行列は 使わず, 線形方程式 B¯b = b (24) iv)ここで挙げたアルゴリズムは「改訂シンプレックス法」と呼ば れるものである。ところで, 「改訂」シンプレックス法が存在する 以上, ただのシンプレックス法も存在する。シンプレックス法と改 訂シンプレックスの違いは, シンプレックス法が基底形式表現その ものを繰り返しの各ステップで直接行列の基本変形によって更新し てゆくのに対し, 改訂シンプレックス法は各ステップで基底行列 (に 対応する添字) を更新してから対応する線形方程式を解く, という点 にある。高速な線形方程式のソルバが利用可能である場合には, ふ つうは改訂シンプレックス法の方が速い。 を解く。一般に, 逆行列の数値計算は計算量が多くか つ誤差の影響を受けやすいので, 数値解析の分野では 逆行列の計算は可能な限り回避される。 また, 各ステップでベクトル ¯cN を求めるためには cT BB−1N を計算する必要があるが, これを求めるため には • 線形方程式 BTc N0 =cBを解いてベクトルcN0 を求める • cN0に行列N の転置を左から掛ける という手順で計算をおこなう。 (cTBB−1N)T =NT(BT)−1cB だから, このようにしても正しい結果が得られること がわかる。 ¯AN =B−1N のすべての成分を計算して も正しい結果は得られるが, これは計算量の観点から は無駄である。 2.3.2 もとの問題の最適解の抽出 シンプレックス法終了時に得られている実行可能基 底行列をB∗とすると, 線形方程式 B∗x∗ B=b (25) を解くことにより最適解の一部x∗Bが得られる。残り の変数x∗N は零である。x∗Bとx∗N からもとの問題の 最適解x∗を構成するためには, どの変数が基底変数に 入っていたかを確認する必要がある。 たとえば, 基底変数が x1, x3, x5 であり, x∗B = [10, 20, 30] となっている場合には, もとの問題の最適 解は [10, 0, 20, 0, 30, 0, ...] である (非基底変数はすべて零である)。 例 2 シンプレックス法の挙動を調べるため, 以下のよ うな 2 変数の簡単な線形計画問題を考える。 min{x1+x2: 1≤ x1≤ 2, 1 ≤ x2≤ 2, x1≥ 0, x2≥ 0} (26) 問題 (26) の制約条件を満たす領域を図 2 に示す。 目的関数はx1+x2であるが, これは傾き−1 の直線 上で一定値となる。そして, この直線が左下にずれる ほど目的関数の値は小さくなる。よって, この問題の 最適解はx1= 1, x2= 1であることがわかる。
1 2 3 0 1 2 3 x1 x2 図 2 (26) の制約条件を満たす領域 では, この問題を標準形に書き直そう。不等式 1≤ x1 は, 新しい変数x3を導入して, x1− x3= 1, x3≥ 0 と書ける。不等式x1 ≤ 2 は, 新しい変数 x4を導入す ることで, x1+x4= 2, x4≥ 0 のようになる。x2に関する不等式も同様にして x2− x5= 1, x5≥ 0 x2+x6= 2, x6≥ 0 という等式に変わる。さて, x = [x1, x2, x3, x4, x5, x6]T とおくと, この問題は min{cTx : Ax = b, x ≥ 0}, A = ⎡ ⎢ ⎢ ⎢ ⎢ ⎣ 1 0 −1 0 0 0 1 0 0 1 0 0 0 1 0 0 −1 0 0 1 0 0 0 1 ⎤ ⎥ ⎥ ⎥ ⎥ ⎦, b = ⎡ ⎢ ⎢ ⎢ ⎢ ⎣ 1 2 1 2 ⎤ ⎥ ⎥ ⎥ ⎥ ⎦, c = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ 1 1 0 0 0 0 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ (27) という線形計画問題に書き直される。これは標準形で ある。 次に, (27) の実行可能基底行列と実行可能基底解を 求めてみよう。基底変数を決めるときには 1 から 6 ま での添字の中から 4 個の添字を選ぶことになるので, 可 能な組み合わせの数は 6 4 = 15通りであるv)。一 v) n r =r!(n−r)!n! である。 方, 実行可能基底解は図 2 の制約条件を満たす領域 (正 方形) の頂点であり, 頂点の数は 4 個しかないから, 上 記の 15 種類の添字の選び方のうち実行可能なものは 4 個しかないと推察される。この推察は正しく, 実際に 計算してみると, 特定の 4 通りの添字の選び方に対し てのみ実行可能基底解が得られ, それ以外の選び方で は,B が正則でなくなるか, B−1b の成分の中に負のも のがあらわれるということがわかる。表 1 に実行可能 な基底変数ベクトルと, 実行可能基底解に対応するx1 とx2の値を示す。表 1 の右の欄の値から, これらが図 表 1 実行可能な基底変数ベクトルと実行可能基底解 xB (x1, x2)の値 [x1, x2, x3, x5]T (2, 2) [x1, x2, x3, x6]T (2, 1) [x1, x2, x4, x5]T (1, 2) [x1, x2, x5, x6]T (1, 1) 2の 4 個の頂点に対応していることがわかる。 続いて, 基底変数ベクトルをxB = [x1, x2, x3, x5]T とした場合のシンプレックス法の動作を見よう。シンプ レックス法は制約条件を満たす領域の頂点をたどって 最適解を探す方法なのだから, 点 (x1, x2) = (2, 2) から 出発した場合, 2 回の繰り返し計算で最適点 (x1, x2) = (1, 1) に到達できるはずである。 ループ第 1 回目xB = [x1, x2, x3, x5]T, xN = [x4, x6]T から出発する (点 (x1, x2) = (2, 2) から始め ることに対応)。これに対応する基底行列は B = [a1, a2, a3, a5] = ⎡ ⎢ ⎢ ⎢ ⎢ ⎣ 1 0 −1 0 1 0 0 0 0 1 0 −1 0 1 0 0 ⎤ ⎥ ⎥ ⎥ ⎥ ⎦ であり, 非基底行列は N = [a4, a6] = ⎡ ⎢ ⎢ ⎢ ⎢ ⎣ 0 0 1 0 0 0 0 1 ⎤ ⎥ ⎥ ⎥ ⎥ ⎦ である。また, cB= [1, 1, 0, 0]T, cN = [0, 0]T となる。ここから,zB= 4, ¯ b = ⎡ ⎢ ⎢ ⎢ ⎢ ⎣ 2 2 1 1 ⎤ ⎥ ⎥ ⎥ ⎥ ⎦, ¯AN = ⎡ ⎢ ⎢ ⎢ ⎢ ⎣ 1 0 0 1 1 0 0 1 ⎤ ⎥ ⎥ ⎥ ⎥ ⎦, ¯cN = −1 −1
が得られる。 (ステップ 1-1) ベクトルcN の第 1 要素, 第 2 要素はと もに−1 で負なので, このうち番号が若い s = 1 を選択 して次に進む。 (ステップ 1-2) ¯a1 = [1, 0, 1, 0]Tで, このうち第 1 成分 と第 3 成分が正である。¯b1/¯a1;1= 2 ¯b3/¯a1;3= 1であ り, ¯b3/¯a1;3の方が小さいので, 添字r = 3 を選択して 次に進む。 (ステップ 1-3) 基底行列は B = [ai1, ai2, ⇓ ai3, ai4] から B= [ai1, ai2,a⇓j1, ai4] に変わる 。ここで xB = [x1, x2, x3, x5]T, xB = [x4, x6]T であったことを思い出すと, i1 = 1, i2 = 2, i3= 3,i4= 5,j1= 4,j2= 6であることがわかる。ゆ えに, 新しい基底行列は B = [a1, a2, a4, a5] と な る 。対 応 す る 基 底 変 数 ベ ク ト ル は xB = [x1, x2, x4, x5]T である。 ループ第 2 回目 ループ第 1 回目の終了時の状態を継 承して,xB= [x1, x2, x4, x5]T, xN = [x3, x6]T から出 発する。表 1 から, これが点 (x1, x2) = (1, 2) に対応す ることがわかる。先と同様の計算により,zB = 3, ¯ b = ⎡ ⎢ ⎢ ⎢ ⎢ ⎣ 1 2 1 1 ⎤ ⎥ ⎥ ⎥ ⎥ ⎦, ¯AN = ⎡ ⎢ ⎢ ⎢ ⎢ ⎣ −1 0 0 1 1 0 0 1 ⎤ ⎥ ⎥ ⎥ ⎥ ⎦, ¯cN = 1 −1 が得られる。 (ステップ 2-1) ベクトルcN の第 2 要素のみ負なので, s = 2 を選択して次に進む。 (ステップ 2-2) ¯a1 = [0, 1, 0, 1]Tで, このうち第 2 成分 と第 4 成分が正である。¯b2/¯a2;2= 2, ¯b4/¯a2;2= 1であ り, ¯b4/¯a2;4の方が小さいので, 添字r = 4 を選択して 次に進む。 (ステップ 2-3) 基底行列は B = [ai1, ai2, ai3, ⇓ ai4] から B= [a i1, ai2, ai3, ⇓ aj2] に変わる 。ここで xB = [x1, x2, x4, x5]T, xB = [x3, x6]T であったことを思い出すと, i1 = 1, i2 = 2, i3= 4,i4= 5,j1= 3,j2= 6であることがわかる。ゆ えに, 新しい基底行列は B= [a 1, a2, a4, a6] と な る 。対 応 す る 基 底 変 数 ベ ク ト ル は xB = [x1, x2, x4, x6]T である。 ループ第 3 回目 ループ第 2 回目の終了時の状態を継 承して,xB = [x1, x2, x4, x6]T,xN = [x3, x5]T から出 発する。表 1 から, これが点 (x1, x2) = (1, 1) に対応す ることがわかる。この点で目的関数は最小となるから, アルゴリズムはこの段階で終了するはずである。これ を確認しよう。先と同様の計算により,zB= 2, ¯ b = ⎡ ⎢ ⎢ ⎢ ⎢ ⎣ 1 1 1 1 ⎤ ⎥ ⎥ ⎥ ⎥ ⎦, ¯AN = ⎡ ⎢ ⎢ ⎢ ⎢ ⎣ −1 0 0 −1 1 0 0 1 ⎤ ⎥ ⎥ ⎥ ⎥ ⎦, ¯cN = 0 0 (28) が得られる。 (ステップ 3-1) ベクトルcNのすべての要素は零である。 よって, すでに最適解が得られているので, 計算を終了 する。 もとの問題の解の復元 計算終了時の基底変数ベクト ルがxB = [x1, x2, x4, x6]T であったことと最適解が xB = ¯b, xN = 0とあらわされることを思い出すと, (28)から, x1 = 1, x2 = 1, x3 = 0, x4 = 1, x5 = 0, x6= 1が最適解となる。よって,x1= 1,x2= 1が目 的関数を最小とする点であることがわかる。 以上の説明から, シンプレックス法のアルゴリズム が, 点 (x1, x2) = (2, 2) から出発し, (x1, x2) = (1, 2) と いう頂点を経由して最適点 (x1, x2) = (1, 1) に到達し ていることがわかる。 2.3.3 ビッグ M 法 第 2.3.1 節で述べたアルゴリズムを利用して与えら れた線形計画問題を解くためには実行可能基底行列B が必要であった。ところで, 問題の規模が大きいときに は実行可能基底行列をひとつ求めること自体も容易で ないことが多い。そのような状況に対処するためには, 第 2.3.1 節のアルゴリズムを実行可能基底行列が未知 の場合にも適用できるように拡張しておく必要がある。
このような拡張としてよく知られているものに, 2 段 階法とビッグ M 法がある。本節では, これらのうちビッ グ M 法について述べる。 (STD) にビッグ M 法が適用可能であるためには (STD)の制約条件の右辺にあらわれるベクトルb の 各成分が非負でなければならない。そこで,b の第 i 成 分biが負で, 対応する行列A の行ベクトルが αTi であ るときには,biを−biで置き換え,αTi を−αTi で置き 換える。これにより, もとの制約条件と等価でb ≥ 0 を満たす制約条件が得られる。 次に, 人為変数と呼ばれる新たな変数 ξ = [ξ1, . . . , ξm]T (29) を導入し, c2= [cT, M · 1Tm]T, A2= [A, Im] (30) と定義して (ただし 1mはm 次の要素がすべて 1 のベ クトル, M は十分大きい正の定数, Imはm 次の単位 行列), (5) から次のような問題を作る。 min cT 2 x ξ :A2 x ξ =b, x ≥ 0, ξ ≥ 0 (31) もとの問題はb ≥ 0 となるように書き換えられていた から,x = 0, ξ = b は (31) の制約条件を満たす。よっ て, これは実行可能基底解である。ゆえに, 問題 (31) に おいて基底変数ベクトルとしてξ を取ると, 実行可能 基底行列が得られるから, ここから出発してシンプレッ クス法を実行し, 最適解を求めることが可能である。 シンプレックス法の終了条件が満たされ, 繰り返し 計算が終了したものとしよう。このとき問題となるの は, もとの問題の最適解をどのようにして復元するか, ということである。 パラメータM が十分大きく, かつ計算終了時に有 限の最小値が得られた場合を考える。この最小値に対 応するx と ξ をそれぞれ x∗とξ∗とする。このとき, ξ∗= 0であればx = x∗が問題 (5) の最適解であるこ とが言え,ξ∗= 0 であれば問題 (5) は解を持たないと いうことが言える。証明は繁雑なので略す。ただし, パ ラメータM をどの程度取れば十分であるかについて の基準はない。よって, 計算終了時にξ∗= 0 となった ときには,M を大きく取り直してシンプレックス法を 再度実行する必要がある。 なお, 計算終了時の結果からもとの問題の最適解x∗ を構成する手順は第 2.3.2 節と同一である。
2.4 内点法
本節では内点法について述べる。 内点法は大きく分けると主内点法と主-双対法の 2 種 類に分類されるのだが, 本節では主-双対法と呼ばれるア ルゴリズムのうち比較的簡単なもの 2 種類を文献 [5, 6] にしたがって紹介する。 まず, いくつか記号の定義をしておく。ベクトルのノ ルムは, x =ni=1x2i によって定義される。ベク トル同士の Hadamard 積を, 以下によって定義する。 xs = [x1s1, . . . , xnsn]T また, 記法の簡単のため, ベクトルx に対し, 1/x = [1/x1, . . . , 1/xn]T, x2 = [x2 1, . . . , x2n]T などと書く。 2.4.1 内点法のアルゴリズム 基準形で記述された最適化問題 (P) minz = cTx : Ax ≥ b, x ≥ 0 (32) を考える。(32) を主問題と呼ぶ。(32) に対し, 双対問 題と呼ばれる問題が, (D) max w = bTy : ATy ≤ c, y ≥ 0 (33) で定義される。ここに, max は maximize の略であり, 目的関数を最大化する問題をあらわす。なお,x ∈ Rn, y ∈ Rmであるvi)。 主-双対法は, 簡単に言うと, (P) と (D) をまとめて 解く方法である。以下では, まず繰り返し計算のため の初期値を定め, 続いて (P) と (D) を同時に解くため に新たな問題 (SP) を定義し, さいごに繰り返し計算の アルゴリズムを述べる。 まず, x と y の適当な初期値 x0 > 0, y0 > 0 を決 める。これは正であればどのように取っても構わない。 次に,p0= 1/x0,t0= 1/y0とおく。パラメータ ¯b, ¯c, β を ¯ b = t0+b − Ax0 ¯ c = p0− c + ATy0 β = 1 − bTy0+cTx0 (34) vi)主問題が標準形 (STD) で書き表されているときには, 対応する 双対問題は (D) max w = T : T +=,≥ 0 という形になる。標準形では, 主問題と双対問題は対称な形になら ない。また, 標準形に対する双対問題では制約条件が付かない変数 があらわれる。とはいっても, これらの相異はどちらかというと形 式的なものである。内点法は標準形の問題にも適用可能である。ア ルゴリズム自体は標準形の場合も基準形の場合もほとんど同じであ る。なお, 歴史的経緯から, 多くの内点法の解説書では標準形が採用 されている。と定義する。続いて, 行列M およびベクトル q を M = ⎡ ⎢ ⎢ ⎢ ⎢ ⎣ 0m×m A −b ¯b −AT 0 n×n c ¯c bT −cT 0 β −¯bT −¯cT −β 0 ⎤ ⎥ ⎥ ⎥ ⎥ ⎦ q = ⎡ ⎢ ⎢ ⎢ ⎢ ⎣ 0m 0n 0 n + m + 2 ⎤ ⎥ ⎥ ⎥ ⎥ ⎦ (35) と定義し, 次の問題を考える。 (SP) minqTξ : M ξ + q ≥ 0, ξ ≥ 0 (36) ここに, ξ = [yT, xT, κ, θ]T は n + m + 2 次元のベク トルである。問題 (SP) のように, 制約条件に対応する 行列M が歪対称行列vii) であり, 非負のベクトルq が 制約条件と目的関数の双方にあらわれる問題を, 歪対 称問題と呼ぶ。 与えられたξ に対し, (SP) の制約条件であらわれた 式Mξ + q が取る値を s(ξ) とおく。すなわち, s(ξ) = Mξ + q (37) と定義する。また, 行列S, Ξ を S = diag[s1, . . . , sm+n+2], Ξ = diag[ξ1, . . . , ξm+n+2] (38) と定義する。ただし, diag[x,. . . , xk] は対角要素が x1, . . . , xk でそれ以外の要素が零の正方行列をあら わす。 以上の準備のもとで, 2 種類の内点法のアルゴリズム について述べる。最初のアルゴリズムは, Dikin ステッ プ法と呼ばれるものである。 内点法 (Dikin ステップ法) (初期化) 精度パラメータε を適当に決める。また, ξ := [(y0)T, (x0)T, 1, 1]T,s := s(ξ) とする。さらに, ス テップ幅α を決める。ただし 0 < α ≤ 1 である。 (ループ) 以下のプログラムを実行する。 while ξTs ≥ ε do begin ξ := ξ + α∆ξ s := s(ξ) end vii) が正方行列で, T= − となるとき, を歪対称行列 という。 ここに, ∆ξ は (S + ΞM)∆ξ = −ξ 2s2 ξs (39) を解くことで得られる。 このアルゴリズムによって最適解が得られるための 十分条件は, ステップ幅α が不等式 0< α ≤ 1 2√n + m + 2 (40) を満たすことであるviii)。とくに, α = 1/(2√n + m + 2) としたときには, 最適解を得るために必要な繰り返し 計算の回数は高々 2(n + m + 1) logn + m + 2 ε (41) 回であることが示せる [5]。なお, (41) において, 小数 点以下は切り上げて考えるものとする。 問題 (36) を初期値ξ0のもとで解き,qTξ∗= 0なる 最適解ξ∗= [y∗, x∗, κ∗, θ∗]T が得られたとしよう。こ のときκ∗= 0 であれば, もとの問題の最適解は ¯ x =xκ∗∗, ¯y = yκ∗∗ (42) となる。κ∗= 0のときには, 最適解は存在しない。す なわち, 制約条件を満たす領域が存在しないか, 目的関 数の値は有限にならない。これについては, 第 2.4.2 節 においてもう少し詳しく説明する。 続いて, ショートステップパス追跡法と呼ばれるア ルゴリズムについて述べる。これは, 次のようなアル ゴリズムである。 内点法 (ショートステップパス追跡法) (初期化) 精度パラメータε を適当に決める。また, ξ := ξ0,s := s(ξ) とする。さらに, σ = 1 − 0.4/√n + m + 2 (43) とおく。 viii)n + m + 2 はベクトル の次元である。平方根の中で数 2 が加 算されていることに特別な意味があるわけではない。この点に関す る誤解を避けるためには, ベクトルの次元を dimと書き, α に 関する制約条件 (40) を 0 < α ≤ 1 2√dim と書き直した方が良いであろう。
(ループ) 以下のプログラムを実行する。 while ξTs ≥ ε do begin µ := sTξ/(n + m + 2) ξ := ξ + α∆ξ s := s(ξ) end ここに, ∆ξ は線形方程式 (S + ΞM)∆ξ = −sξ + σµ1 (44) の解である。各々のステップで (44) を解き直す必 要があることに注意せよ。 Dikinステップ法は最急降下法の一種である。これ に対して, ショートステップパス追跡法は Newton 法 の一種である。一般にショートステップパス追跡法は Dikinステップ法よりは高速である。Dikin ステップ法, ショートステップパス追跡法のいずれも, 理論的な特 性は良いが, 実際の問題に適用すると実行が極めて遅 いという特徴がある。より実用的なアルゴリズムとし ては, ロングステップパス追跡法, Mehrotra の予測子・ 修正子法などがある. 以下では, まず第 2.4.2 節において, 本節で利用され た行列M などの意味を説明する。続いて, 第 2.4.3 節 において Dikin ステップ法の意味を説明する。第 2.4.4 節では, Newton 法に基づく内点法の一般的なアルゴリ ズムについて説明する。これらは実験とは直接関係な いので, 興味のない読者は読み飛ばしてよい。 2.4.2 線形計画問題 (SP) の意味 本節では, 第 2.4 節の行列M などの意味についてい くつかの証明を省いた上で述べる。この節は実験その ものには関係ないので, 興味のない読者は読み飛ばし てよい。 主問題 (P) とその双対問題 (D) が与えられているも のとする。主問題の制約条件 (Ax ≥ b, x ≥ 0) を満た す点の集合をP と書き, 双対問題の制約条件 (ATy ≤ c, y ≥ 0) を満たす点の集合を D と書く。P や D が空 集合となることもある。このような場合には, 主問題 や双対問題には解が存在しない。 さて, 適当に取られたx ∈ P と y ∈ D に対し, bTy ≤ (xTAT)y = xT(ATy) ≤ xTc より, bTy ≤ cTx と いう不等式が成り立つことがわかる。すなわち, 双対 問題の目的関数値は必ず主問題の目的関数値以下とな る。このことから, あるx∗ ∈ P と y∗ ∈ D が等式 bTy∗ =cTx∗ を満たせば, x∗は主問題の最適解であ り,y∗は双対問題の最適解であるということがわかる。 実は, この逆も成り立つ。すなわち, 主問題と双対問 題がそれぞれ有限の最適解x∗とy∗を持つとき, 等式 bTy∗=cTx∗が成り立つということが言える。後半の 結果を双対定理という。双対定理の証明はやや繁雑な のでここでは述べない。 さて, 最適解のための等式bTy = cTx に不等号を追 加した不等式 bTy ≥ cTx (45) を考えよう。不等式bTy ≤ cTx はつねに成り立つの だから, これは実質的にもとの等式と変わらない。(45) の右辺を左辺に移項し, 主問題の制約条件, 双対問題の 制約条件の両辺に−1 を乗じたものと合わせて整理す ると, Ax − b ≥ 0 −ATy + c ≥ 0 bTy − cTx ≥ 0 という不等式が得られる。これを行列の形でまとめると, ⎡ ⎢ ⎣ 0 A −b −AT 0 c bT −cT 0 ⎤ ⎥ ⎦ ⎡ ⎢ ⎣ y x 1 ⎤ ⎥ ⎦ ≥ ⎡ ⎢ ⎣ 0 0 0 ⎤ ⎥ ⎦ , x ≥ 0, y ≥ 0 (46) という形になる。この段階ではまだ (46) が解を持つか どうかは明らかでないことに注意しよう。 ここで, (46) の左辺のベクトルの最下段にあらわれ た数 1 を変数κ で置き変えた不等式 ⎡ ⎢ ⎣ 0 A −b −AT 0 c bT −cT 0 ⎤ ⎥ ⎦ ⎡ ⎢ ⎣ y x κ ⎤ ⎥ ⎦ ≥ ⎡ ⎢ ⎣ 0 0 0 ⎤ ⎥ ⎦ , ⎡ ⎢ ⎣ y x κ ⎤ ⎥ ⎦ ≥ ⎡ ⎢ ⎣ 0 0 0 ⎤ ⎥ ⎦ (47) を考える。(47) を満たす (yT, xT, κ) が求められてい て, かつκ > 0 であれば, (yT/κ, xT/κ, 1) は (46) を満 たす。すなわち, 新しい変数κ を導入することはもとの 問題 (46) の解を求めることの障害とはならない。では, 新しい変数κ を導入することには何か利点はあるのだ ろうか。まず (47) が自明な解x = 0, y = 0, κ = 0 を 持つことに注意する。ところで, 内点法を用いて問題 (47)の解を求めると, 問題 (46) に解があるときにはκ は零以外の値に収束し, 問題 (46) に解がないときには κ は零に収束することが示せる (証明は略す)。すなわ ち,κ を (46) の解の存在判定に使うことができる。
さて, (47) の左端の行列は先に内点法のアルゴリズム で使った (35) と似てはいるが, 同一ではない。(35) の 行列M の右端および下段の要素は (47) にはないもの である。これはどのような経緯で出て来たのだろうか? 実は, 線形計画問題を内点法で解くためには初期値 として制約領域の内点 (制約条件から等号を取り除い た不等式を満たす点) が必要になるのだが, 内点を 1 個 定めることはそれほど容易でない。そのため, 代替策 として, 新たに変数を 1 個追加して, (47) から内点が既 知の拡大された問題に変換する。(35) は (47) にそのよ うな手法を適用した結果なのである。以下では, この 方法について見てゆこう。 初期値x0> 0, y0> 0 を適当に決め (各成分が 0 で なければ何でもよい),p0 = 1/x0,t0 = 1/y0とおく。 パラメータ ¯b, ¯c, β を (34) によって定義する。また, 行 列M およびベクトル q を (35) によって決め, (36) に よって定められた問題 (SP) を考える。すると, 初期値 ξ0= ⎡ ⎢ ⎢ ⎢ ⎢ ⎣ y0 x0 1 1 ⎤ ⎥ ⎥ ⎥ ⎥ ⎦ (48) が問題 (36) の内点となるのである。これを見ておこう。 ξ0> 0 は定義より明らかだから, Mξ0+q > 0 となる ことを見ればよい。これは, (34) を使えば, Mξ0+q = ⎡ ⎢ ⎢ ⎢ ⎢ ⎣ Ax0− b + ¯b −ATy0+c + ¯c bTy0− cTx0+β −¯bTy0− ¯cTx0− β + (n + m + 1) ⎤ ⎥ ⎥ ⎥ ⎥ ⎦ = ⎡ ⎢ ⎢ ⎢ ⎢ ⎣ t0 p0 1 1 ⎤ ⎥ ⎥ ⎥ ⎥ ⎦ (49) となることからわかる。 2.4.3 Dikinステップ法の意味 続いて, 繰り返し計算で用いられる探索ベクトル∆ξ を求める式 (39) の意味を簡単に説明しておく。 s(ξ) = Mξ + q > 0, ξ > 0 を満たす ξ が与えられ ているものとする。記法の簡単のために,s(ξ) を s と 略記する。さて,M は歪対象行列だから, ξTMξ = 0 となるix)。ゆえに, 目的関数の値qTξ は qTξ = ξTs とも書けるx)。さて,ξ を ξ + ∆ξ に変更し, ∆s = s(ξ + ∆ξ) − s(ξ) とおくと, 目的関数の値は (ξ + ∆ξ)T(s + ∆s) − ξTs = ξT∆s + ∆ξTs だけ変動するxi)。そこで, 単純に考えると, ξT∆s + ∆ξTs が最も減少する方向に ∆ξ を取れば良さそうに 思える。しかし, これは必ずしも望ましくない。解が制 約条件を満たす領域の境界に近付きすぎると計算効率 が落ちることがあるのだが, 上述の ∆ξ の取り方では この問題が考慮されていないからである。そこで, 新 たに ∆ξξ+∆s s ≤ 1 (50) という制約条件を課し, (50) を満たす領域 (この領域を Dikin楕円体と呼ぶ) の中からξT∆s + ∆ξTs を最小 化する ∆ξ を求めたものが (39) である。(50) の制約条 件の形から, ξ あるいは s の成分のうち値が 0 に近い ものについては, 対応する ∆ξ あるいは ∆s の成分の 値はあまり大きく取れず, 結果としてξ あるいは s が 0に近付くことが抑制される。この方法は,ξ や s の大 きさに応じて重みを付けた, 最急降下法の一種と考え られる。 ここで, ごく簡単な歪対称問題について, 制約領域の 各点で Dikin 楕円体がどのような形になるかを見てお こう。ただし, ここで取り扱うのは 2 変数の人工的な 歪対称問題であり, 主問題 (P) と双対問題 (D) から決 まる歪対称問題ではない。これは, 主問題 (P) および 双対問題 (D) が決まる歪対称問題は最低でも 4 変数と なり, 2 次元のグラフでは Dikin 楕円体をうまく書き表 せないからである。 例 3 2 変数の歪対称問題 min ⎧ ⎨ ⎩ 1 0 T ξ1 ξ2 : 0 −1 1 0 ξ1 ξ2 + 1 0 ≥ 0 0 , ξ1 ξ2 ≥ 0 0 (51) ix)T がスカラーであることと T = − から, T = ( T ) T=T T = − TT , よって 2T= 0。 x)T = 0 となることに注意すれば, T = T( +) = T がいえる。 xi)(∆)T∆= (∆)T∆= 0
を考える。 (51)の制約条件は, 0 ≤ ξ1, 0≤ ξ2 ≤ 1 と書き直さ れる。目的関数はξ1なので,ξ1= 0を満たす任意の点 で目的関数は最小となる。 図 3 に (51) の制約条件を満たす領域の各点における Dikin楕円体の形を示す。制約条件を満たす領域は図 0.2 0.4 0.6 0.8 1 1.2 0.2 0.4 0.6 0.8 1 0 ξ1 ξ2 図 3 2 変数の歪対称問題に対する Dikin 楕円体 中の網のかかった領域である。この図から, 境界に近 付くにつれ Dikin 楕円体が潰れてゆくことがわかる。 2.4.4 Newton法に基づく内点法の一般的なアルゴ リズム 本節では, Newton 法に基づく内点法のアルゴリズム について簡単に述べる。 第 2.4.2 節の結果から, 主問題と双対問題を解くこと は, 問題 (SP) minqTξ : Mξ + q ≥ 0, ξ ≥ 0 を解くことに帰着されることがわかっている。ここに, M は歪対称行列であり, q は各成分が非負のベクトル である。そこで, 本節では (SP) を出発点とする。 2.4.4.1 強相補性条件 行列M が歪対称であること から, 目的関数値はqTξ = sTξ を満たすxii)。 さて, 問題 (SP) は自明な最適解ξ = 0 を持つ。しか し, 第 2.4.2 で説明したように, 我々はξ のある成分が 収束する値が 0 であるか否かに応じて主問題 (P) と双 対問題 (D) の解の存在性の判定をするので, 自明な解 は我々の望むものではない。ところで,ξ が最適であり xii)先と同様に,=+とする (sTx = 0), かつ ξ の要素の中の零でないものがあった としよう。ξ ≥ 0, s ≥ 0 であることを思いだすと, ξ と s の成分は, • ξi= 0ならsi> 0, • ξi> 0 なら si= 0 という条件を満たすことがわかる。この条件を強相補 性条件という。強相補性条件は, より簡単に, s + ξ > 0 (52) とも書ける。以上の議論により, 我々は問題 (SP) の強 相補性条件を満たす最適解を (もしあれば) 見付けなけ ればならないことがわかる。 2.4.4.2 中心パス ここで, 中心パスという概念を導 入しておく。1 をすべての要素が 1 のベクトルとし,µ を正の定数として, sξ = µ1 (53) なる等式を考える。ここでは証明を省くが, (SP) の制 約条件を満たす領域が内点を持つ (開集合を含む) とき, (53)と (SP) の制約条件を連立させた非線形連立方程式 s = Mξ + q, sξ = µ1 (54) をξ について解くと, 与えられた µ に対して制約領域 の内点にあるξ(µ) が一意的に定まることが言える。さ らに, 次の事実が成り立つ。 • ξ(µ) は µ の連続関数である。 • µ → 0 とすると, ξ(µ) は強相補性条件を満たす点 に収束する。 これらの事実から, 軌道ξ(µ)(これを中心パスと呼ぶ) を追跡することにより強相補性条件を満たす解が求め られることがわかる。中心パスは, 制約条件を満たす 領域の内部を通って最適解に到達する, 幾何学的に性 質の良い軌道である。 2.4.4.3 アフィンスケーリング方向と中心化方向 以 上の議論を踏まえた上で, 2 種類の異なる探索ベクトル の決め方について述べよう。 まず最初に, (SP) を変数x と s に関する非線形関数 F (ξ, s) = Mξ − s + q sξ (55)
の零点を求める問題とみなす。点ξ と s = Mξ+q が与 えられているとき,ξ =ξ + ∆ξaおよびs=s + ∆sa としてF (ξ, s) = 0を ∆ξa と ∆saについて Tailor 展開によって近似的に解くと, ∆ξaと ∆saは方程式 M −I S Ξ ∆ξa ∆sa = 0 −sξ (56) の解となる。なお, 適当な正の定数α を取り, ξ =ξ + α∆ξa およびs =s + α∆saとすると, M(ξ) +q = s + αM∆ξa =s となる。よって, 任意のα に対し, sとξは等式s= Mξ+q を満たす。(56) を解くことによって得られる 探索ベクトル (∆ξTa, ∆sTa)T をアフィンスケーリング 方向と呼ぶ。 アフィンスケーリング方向は目的関数を減少させる には都合が良いが, 制約条件を見たす領域の境界に解 を誘導する可能性があるという点では不都合である。 ところで, 解を内点に誘導するには, 中心パス上の点を 利用するのが良い。そこで, パラメータµ を適当に決 め, それ対応する中心パス上の点を求める問題を考え る。パラメータµ の選び方にはいろいろなものが考え られるが, たとえば目的関数値を現在の値の 1/ dim ξ にすることを目指す, すなわちµ = (1/ dim ξ)sTξ と するという方法が考えられる。この問題も, 先と同様 に, 関数 G(ξ, s) = Mξ − s + q sξ − µ1 (57) の零点を求める問題に帰着される。ξ=ξ +∆ξc,s= s + ∆scとし, (57) を Tailor 展開によって近似的に解 くと, ∆ξcと ∆scは方程式 M −I S Ξ ∆ξc ∆sc = 0 −sξ + µ1 (58) の解となる。(58) を解くことによって得られる探索ベ クトル (∆ξTc, ∆sTc)T を中心化方向と呼ぶ。 これらの探索ベクトルを用いた解法は, Newton 法に よる非線形方程式の求解の一種となる。 2.4.4.4 アルゴリズムの一般形 アフィンスケーリ ング方向は目的関数を減少させる方向, 中心化方向は 中心パスに向かう方向であった。Newton 法型の内点 法のアルゴリズムでは, アフィンスケーリング方向と 中心化方向の線形結合を使って解を更新してゆく。 σ を 0 以上 1 以下の定数とし, 探索ベクトルを ∆ξ ∆s = (1− σ) ∆ξa ∆sa +σ ∆ξc ∆sc (59) とすると, 探索ベクトル (∆ξT, ∆sT)T は M −I S Ξ ∆ξ ∆s = 0 −sξ + σµ1 (60) の解となる。 なお, (60) は ∆ξ と ∆s に関する線形方程式である が, (60) の上半分は ∆s について解けて ∆s = M∆ξ となるので, 実際にはこれを (60) の下半分に代入した S + ΞM∆ξ = −sξ + σµ1 (61) を ∆ξ について解けばよい. (60)の解 (∆ξT, ∆sT)T を探索ベクトルとして利用 する内点法のアルゴリズムの一般形は, 以下に述べる ようなものである。 内点法 (Newton 法による解法) (初期化) 精度パラメータε を適当に決める。また, ξ := ξ0,s := s(ξ) とする。 (ループ) 以下のプログラムを実行する。 while ξTs ≥ ε do begin σ を決める µ := sTξ/(n) とする (61)を解いて ∆ξ, ∆s を決める α を決める ξ := ξ + α∆ξ s := s + α∆s end 第 2.4.1 節で述べたショートステップパス追跡法は, 上 記のアルゴリズムにおいてα = 1, σ = 1 − 0.4/√dimξ と固定した特別な場合である。大規模な問題では, こ のようにσ を決めると, 探索ベクトルが中心化方向と ほぼ一致してしまうため, ショートステップパス追跡 法は極めて低速になる。 ループ内の各ステップでパラメータα や σ の値を切 り換えると, より効率の良いアルゴリズムが構成しう ることが知られている。興味のある読者はたとえば文 献 [2, 5, 6] を参照してほしい。 例 4 例 3 の歪対称問題を用いて, アフィンスケーリン グ方向と中心化方向の各点における変動を見よう。まず 最初に, 例 3 の問題 (51) において, 中心パスはξ2= 1/2
という直線であり, 強相補解は (ξ1, ξ2) = (0, 1/2) とい う点である。図 4 に各点におけるアフィンスケーリング 方向と中心化方向を示す。この例ではアフィンスケー リング方向は中心パスと平行になる。それゆえ, アフィ ンスケーリング方向のみを用いて解を更新すると強相 補解に到達できない。また,µ = (1/2)sTξ の場合の中 心化方向は中心パスと垂直となり, これを用いて解を 更新すると目的関数は減少しない。 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 ξ 1 1 ξ2 中心パス (a)アフィンスケーリング方向 ξ1 ξ2 中心パス 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 (b)中心化方向 図 4 各点における探索ベクトルの方向
