MATLABクローンによる大域的最適化（2）−Octaveで作る改訂単体法一
久野 誉人……川…＝‖‖‖‖‖‖‖‖‖‖‖＝‖＝‖＝‖＝＝＝‖‖‖‖‖‖＝‖‖‖‖‖‖‖‖‖‖＝＝‖＝＝‖‖＝＝‖‖＝＝‖‖＝‖‖‖‖‖‖‖‖‖‖＝州‖ll……l州l川………＝‖‖‖＝‖‖‖‖‖‖‖‖‖‖‖‖‖＝‖‖‖‖‖‖‖‖‖‖‖‖‖‖‖＝‖‖‖＝‖＝‖‖‖‖＝‖‖‖＝‖＝＝‖‖＝‖‖‖＝‖酬

B＝Ⅰ，N＝A，Cβ＝0丁，C〃＝CT

5．前回の宿題と置くだけで，（1）の実行可能な辞書

前回に出した宿題はすんなりできただろうか？

Ⅹβ＝B￣1b−B−1Nx〃（2）

宿題1．線形計画問題

z＝Cβ8￣1b＋（c〃一CβB−1N）Ⅹ〃が得られる．したがってフェイズⅠを飛ばして，ただ最大化 cTx（1）

条 件 Ax≦b，Ⅹ≧0

ちに改訂単体法のフェイズⅠⅠを実行することが可能だ．

を解くための改訂単体法をOctaveかMATLABを釈迦に説法になるかもしれないが，念のためにフェイ使ってプログラム化しなさい．ただし，A∈RmX〝，C

ズⅠⅠのアルゴリズム（例えば，文献［1，3，6］を参照）∈R乃とし，b∈Rmの成分はすべて非負とする． tを書いておこう：アルゴリズム1．
実は，これと同じ内容の宿題を毎年4月，筆者の研究室の学生たちにアルゴリズムを単体法に限定しない

ステップ0．xβ：＝b，Ⅹ〟：＝0とおく．
ステップ1．連立1次方程式yB＝Cβを解いてyをで出題している．7月の夏休み第一週，学生たちの作求め，被約費用F〃：＝C〃−yNを計算する．百〃≦ってきたプログラムを1台のコンピュータに集め，0ならば終了（Ⅹβ，Ⅹ〃が最適解）．そうでなければ，

Octave上でランダムに生成した様々な大きさの問題

百〃から値が正となる第5成分を選ぶ．

（1）を解いて計算時間を競い合うのだが，優勝者はそのステップ2．：Nの第s列をaとし，Bd＝aを解いて夜の打ち上げで他の参加者から奮ってもらえる約束だ．

dを求める．d≦0ならば終了（問題は非有界）．そ

参加資格は学生に限っていないので，もちろん筆者もうでなければ，5：＝Ⅹβ−′d≧0を満たす最大の′参加するのだが，残念ながら一度も優勝したことがなを求め，古から値ゼロの第γ成分を選ぶ．

い．これは筆者が手を抜いているわけではなく，

ステップ3．xβ：＝石，Ⅹ〟の第ざ成分の値を′とおOctaveやMATLABを使うとき，Bの第γ列とNの第ざ列aを入れ換えるピボ（a）教科書に書かれたプログラムの高速化に係わる常ット演算を行う．これに合わせてcβ，C〃，Ⅹβ，Ⅹ〃

識やノウハウが必ずしも通じない，の成分を入れ換えたのち，ステップ1にもどる．■（b）プログラミングに慣れない学生がてこずるはずの行列演算も簡単に，しかも高速に実行できる

この記述自体，ステップ1から3までのループで構成ことが大きな原因である，と自分を慰めることにしてされており，「ループは可能な限り使うな」といういる．そんなわけで，これから説明する宿題の解答はOctaveによる効率的な計算の鉄則に反するが，残念

必ずしもベストなものといえないし，もっと賢明なやながら単体法に限らず，解を反復して改善する最適化

り方があることも否定しない．

アルゴリズムには少なくとも一つのループが必要で，すでに述べたとおり，スラック変数を基底変数Ⅹβ，元の変数を非基底変数ⅩⅣとし，単位行列をⅠ∈R∽×mで表すことにしてこれを使わないことにはプログラムにならない．

6．ニつの連立1次方程式
改訂単体法の各反復で最も手間のかかるのが，ステ

くの たかひと筑波大学大学院システム情報工学研究科

ッ701のyIi＝Cβとステップ2のBd＝aの二つの連〒305−8573つくば市天王台1一二卜1

立1次方程式の求解である．これらを毎回まともに解2005年11月号© 日本オペレーションズ・リサーチ学会. 無断複写・複製・転載を禁ず.

（23）75丁いていたのでは非効率極まりないアルゴリズムとなる装しようとすると，例の鉄則の壁に真向からぶち当たところだが，BTRAN（backwardtransformation），ってしまう．例えばFTRANならば，それまでに生FTRAN（forwardtransformation）と呼ばれる巧妙成されたEl，E2，…，E々のイータ列とその列番号をそな計算テクニックのお陰で，改訂単体法は誕生かられぞれ格納する行列DとベクトルRを用意して，50年以上を経た今も現役として活用されている．

d＝a；

forj＝［1：k］

6．1基底行列のイ一夕分解

E＝eye（m）；E（：，R（j））＝D（：，J）；

ステップ3でBのγ列とNのぶ列aを入れ換えて

dt＝E＼d；d＝dt；

できる新しい基底行列虫は，Bd＝aが成り立つこと

endforからerを第γ基本列ベクトルとして虫＝B＋（a−Ber）eF＝B【Ⅰ＋（d−er）e7］のように書くことができる．ここでE＝Ⅰ＋（d−er）e7と表すことにすると，Eは単位行列の第γ列だけをdに置き換えたイ一夕行列（Etamatrix）で，各行の非ゼロ成分は高々二つにすぎない．この方法で基底行列の更新を続ければ，初期基底行列は単位行列Ⅰなので々回めの反復後にはのようなループの入ったプログラムを書かざるをえな
い．この3行めで組込み関数eye（）の生成するm次単位行列Eは，R（j）列が行列Dのj列に置き換えられてイータ行列EJとなり，次の行で方程式EJdJ＝dJ＿1が解かれる仕組みである．簡潔ではあるものの，このループを単体法の反復ごとに繰り返すプログラム

がOctaveによって効率よく働いてくれるはずのない
ことは，前回の観察から容易に想像がつくだろう．

B烏＝EIE2・‥E鳥

6．2 0ctaveではどうする？

のように基底行列B烏が々個のイータ行列El，…，E貞の積にイ一夕分解（Etafactorization）される．したがって，次の反復のステップ1でyB烏＝Cβを解くとならば，いっそのことd＝B＼aとやってもOctave

はFTRANのプログラムよりも高速に処理してくれ
るかもしれない．しかし，それでは毎回，基底行列B

きには

y々E鳥＝Cβ，y烏一1Eト1＝y々，…，ylEl＝y2，の順（BTRAN）にyk，yh＿1…，ylを計算し，またス
テップ2でB鳥d＝aを解くときには

Eldl＝a，E2d2＝dl，…，Ekdk＝dh＿1，の順（FTRAN）にdl，d2，．．．，dkを計算することで，それぞれの解y＝yl，d＝d烏を求めることができる．
係数行列はイータ行列であり，分解された各方程式を解くには0（弼）の手間しかかからない．したがって々が大きくなければ，yI主烏＝Cβ，B烏d＝aをガウスの消去法で0（椚3）の手間をかけて愚直に計算するよりもずっと効率がよい．また，反復回数々が大きくなりすぎたときにはB烏を置換行列P烏と三角行列し，Uhに再分解（refactorization）し，再びの逆行列B−1を一から計算するに等しく，「改訂」のありがたみがない．そこでBではなくB−1の方を保持し，これをSherman−Morrison−Woodburyのラン
ク・ワン更新（rank oneupdate）（例えば，文献［4］を参月別で次の基底行列虫の逆行列虫￣1に更新することにしよう．

イータ行列Eを介してBと虫の間には虫＝BEの関係があることは見たとおりだが，両辺の逆行列をとると虫￣1＝E￣1B￣1になり，Eの逆行列さえ求まれば
B￣1を虫￣1に更新できることが分かる．計算すれば容易に確かめられるが，E￣1も単位行列とはγ列だけが異なるイータ行列で，drをdの第γ成分としてE￣1＝Ⅰ−（1／あ）（d−er）e7と書くことができる．したがって

B如上＝P烏L烏U烏E糾1Eた＋2…E烏＋′と更新を続ければよい．行列L烏，U烏はともに擁個のイータ行列の積として表せるので，再分解からの反復回数Jが大きくなければ，yB烏＋イ＝CβもB如〜d＝aもやはりBTRAN，FTRANによって速やかに解をとなり，Octaveで処理しやすいようにdの第r成分
drを−1に置き換えたaを使って，さらに変形すると虫￣1＝（Ⅰ−ere7）B￣1−（1仏）aeFB−1

求められることになる．

改訂単体法で効率の鍵を握るBTRAN，FTRAN
は，CやFortranなどの普通の言語を使って実装したときには期待どおりにすばらしい効果を発揮してくれる．ところが，OctaveやMATLABでこれらを実
丁58（24）虫￣1＝B￣1−（1仏）（d−er）e7B￣1になる．つまり，B￣1の第γ行をすべてゼロにした行

列から，列ベクトル（1／dr）aとB￣1の第γ行との積を差し引けば虫￣1が得られるのである．この更新には
0（∽2）の基本演算しかかからないので，y＝B￣1cβと© 日本オペレーションズ・リサーチ学会. 無断複写・複製・転載を禁ず.

オペレーションズ・リサーチ

d＝B￣1aの計算を合わせても二つの方程式は0（∽2）

1．Oe−8もこのステップで準備しておこう．

で処理できることになる．具体的には，B￣1をBiで表してステップ1．ランク・ワン更新した基底行列の逆行

dr＝d（r）；d（r）＝−1．0；d／＝dr；

列Biが使えるはずなのでyB＝CBの解はy＝CC（B）＊

br＝Bi（（，：）；8i（（，：）＝0．0；Bト＝d＊dr；

Biで求め，被約費用百ⅣをとすればBiは鱈−1に更新される．

rc＝CC（N）−y＊AA（：，N）；

ところで，イータ分解に再分解が必要となる第一の理由はyB＝Cβ，Bd＝aの求解で反復に比例する手間で計算する．同じことをforループで行えば
forj＝［1：n］

rc（j）＝CC（N（j））−y＊AA（：，N（j））；

がかかるためだが，他に0（椚）で増加するイータ行
列El，…，E烏の記憶量，累積する丸め誤差の二つを解

endfor

消するためでもある．80年代の名著［1］によると，大のように書けるが，もちろんこれを使うべきではない．

規模問題に対して当時，20反復に1回程度の高頻度

ただ，ここで注意したいのはrc（j）がx（j）でなく，×で再分解が行われていたようだ．しかし，現在はパソ（N（j））の被約費用である点だ．ピボット列の選択は，

コンでも当時のワークステーションよりずっと大きな

ポピュラーな最大係数規則を用いることにしよう：

主記憶を持っているし，そもそもランク・ワン更新で［z，S］＝maX（rc）；

保持しなければならないのは常に同じサイズ∽×別

ここでzにべクトルrcの最大成分の値が代入されるの行列B￣1だけである．しかも，yB＝Cβ，Bd＝aにが，それがZ8以下ならば現在のX（1：n）が最適解で対する求解の手間は反復の多寡にかかわらず0（∽2）

ある．そうでなければ，値がzのs番めの成分に対応

で，これも変わらない．したがって，気にしなければする非基底変数x（N（s））が基底に入る候補となる．

ならないのは丸め誤差のみで，再分解の頻度は20年前よりもはるかに少なくて大丈夫だろう．また再分解

ステップ2．繁雑なのでns＝N（s）とおくことにして，といっても，ここで欲しいのはLとUではなくて

Bd＝aの解をd＝Bi＊AA（：，nS）で求める．有界性のB−1そのものであるが，Octaveの組込み関数inv（）を判定には，まず使えばBi＝inv（B）のように造作もなく求められる．

D＝find（d＞Z8）；

によって値がZ8よりも大きなdの成分の行番号をD

7．Octave版改訂単体法に格納する．そして組込み関数isempty（）を使って改訂単体法のプログラム化で最も厄介な課題がかた

isempty（D）が1か0，つまり Is D empty？ をチづいたので，残りの部分のプログラム化をアルゴリズ

ェックする．もしも yes ならば，dの成分がすべて

ム1のステップを追いながら，Octaveの便利な道具

Z8以下であり，問題は非有界である．そうでなけれを使って解説していくことにしよう．

ば，［t，（］＝min（x（B（D））．／d（D））；｛＝D（｛）；

ステップ0．入力を研×犯行列A，研次列ベクトル

によって基底から掃き出す変数x（B（「））を特定する．

b，乃次行ベクトルCの三つとして，
N＝［1：n］；B＝［n＋1：n＋m］；

ステップ3．ピボット（「，S）が定まったので，あと

X＝［zeros（n，1）；b］；cc＝［c，ZerOS（1，m）］；

はこれを中心にピボット演算を行うだけだ．その前に，

Bi＝eye（m）；AA＝［A，Bi］；

ステップ2で求めたdとtを使って解を更新しよう：のように下ごしらえする．行列AAは辞書（2）のNとx（B）−＝t＊d；x（ns）＝t；

Bを並べたもので，本当はAA＝［A，eye（m）］とすべ

ピボット演算は

きだが，この段階ではまだ基底行列もその道行列Bi

N（s）＝B（（）；B（（）＝nS；

も単位行列であるのでeye（m）の呼び出しを省いていで完了する．このあと，節6．2に説明した方法で基底る．これ以降，BはAA（：，B），NはAA（：，N）で，行列の逆行列Biを更新すればよい．

cBとcN，XB，ⅩNもそれぞれcc（B），CC（N），X（B），X

7．1関数simplex（）（N）で参照することができる．また，丸め誤差を考慮

以上のレシピをもとに改訂単体法をOctaveの関数

してゼロの代わりに用いる10￣8程度の小さな値Z8＝
2005年11月号にまとめよう：© 日本オペレーションズ・リサーチ学会. 無断複写・複製・転載を禁ず.

（25）丁59

function［val，SOl，k］＝SimpIex（A，b，C）％％StepO％％
Z8＝1．Oe−8；RF＝512；k＝0；のだが512反復に1度の再分解を行うことにしてある．
ステップ3のピボット演算ののちにkを一つ増やし，剰余を返す組込み関数rem（）を使ってkがRFで割り［m，n］＝Size（A）；N＝［1：n］；B＝［n＋1：n＋m］；

切れたときにだけ再分解が行われる．

x＝［zeros（n，1）；b］；cc＝［c，ZerOS（1，m）］；
Bi＝eye（m）；AA＝［A，Bi］；

どうです，ルーフ0はステップ1から3を繰り返す
wh‖eループーつで済んだでしょう．それでは，関数
simpIex（）がコンピュータ上で正しく動いてくれるか

OPt＝0；

どうか確かめてみることにしよう．

Whileopt＝＝0

7．2 動作確認と性能評価％％Stepl％％

とりあえず，線形計画法の適当な教科書から答の分

y＝CC（B）＊Bi；rc＝CC（N）−y＊AA（：，N）；［z，S］＝maX（rc）；

かっている例題を選んでsimplex（）に解かせてみよう．

jf z＜＝Z8

例えば，文献［1］の2章にある例題：
3

一

＋

＋

裁
3

一

J

2

ifisempty（D）

裁∬

D＝find（d＞Z8）；

＋

2 J ∫

∬

一2

ns＝N（s）；d＝Bi＊AA（：，nS）；＋＋

最大化 条 件％％Step2％％

＋

5

2

＋

∫

else

J

5

OPt＝1；vaJ＝y＊b；sol＝X（1：n）；

3J3
J3≦3
3J3≦2
（3）

2」r3≦4∬3≦2

れ，∬2こr3≧0．

OPt＝−1；vaI＝lnf；sol＝X（1：n）；

else

最適解はれ＝32／29，∬2＝8／29，∬3＝30／29で最適値％％Step3％％［t，（］＝min（x（B（D））．／d（D））；｛＝D（｛）；
x（B）一＝t＊d；x（ns）＝t；は10である．ここで間違っても巡回を起こす問題例
は選ばないこと．読者はすでに気づいていると思うが，Simptex（）には何の巡回対策も施していないので終了

N（s）＝B（（）；B（（）＝nS；k＋＋；

しない可能性がある．Octaveを起動して，

if rem（k，RF）％％Rank oneupdate％％

octave：1〉A＝［131；一103；2 −12；23 −1］；

dr＝d（｛）；d（｛）＝−1．0；d／＝dr；octave：2〉b＝［3；2；4；2］；c＝［553］；

br＝Bi（「，：）；‡Bi（「，：）＝0．0；とデータを打ち込み，これをsimplex（）に渡してOCtaVe：3〉［var，SOl，k］＝Simprex（A，b，C）

Bi−＝d＊br；else

Val＝10％％Refactorization％％

SOl＝Bi＝hv（AA（：，B））；

1．10345x（B）＝Bi＊b；x（N）＝0．0；

0．27586

endif

1．03448

endif

k＝3

のように瞬時に答が返ってくればひとまず安心だ．通

endi†

endwhiJe

常は，エラーメッセージがいくつも現れて，それを一

endfunction

つずつデバグすることになるが，メッセージはかなり

この関数simplex（）に問題（1）の係数行列Aと右辺ベ
クトルb，費用ベクトルcを引数として与えると，（1）が有界のときには最適値val，最適解soし それにプ

詳細で，例えばサイズの合わない行列の積をとっていたりするとerror：OPeratOr＊：nOnCOnformant arguments（oplislX

ログラムの終了までに要したピボット回数kが返さ

3，OP2is4×4）

れる．問題（1）が非有界ならばval＝lnfとなるが，error：eValuatjngbinaryoperator ＊ near］inelO，COlumn

Octaveでrnfは＋∞の意味である．記号％から右は，

13

LaTeXと同じように無視されるので注釈と思ってよい．定数RFは再分解の頻度を示し，特に根拠はない

丁印（26）のように報告されるのですぐに修正できる．

© 日本オペレーションズ・リサーチ学会. 無断複写・複製・転載を禁ず.

オペレーションズ・リサーチ

バグ取りが終わって（3）が無事に解けるようになった

k＝1356

ら，今度はもう少し大きな問題でsimpJex（）の動作確

0CtaVe：12〉A＝rand（500，300）rrand（500，300）；

認をしたい．何せ（3）は手で解くことを前提にしたオモ

octave：13〉b＝OneS（500，1）；c＝rand（1，300）；チャ問題で，実際わずかk＝3回の反復しかかかっておらずテストにならない．そうなると，適当なテスト
問題を見つけてくる必要があるが，自分で作った問題

octave：14〉［tl，ul，Sl］＝CPutime（）；＼〉［val，SOl．k］＝Simplex（A，b，C）；＼〉［t2，u2，S2］＝CPutime（）；OCtaVe：15〉t2−tl，u2−ul，kを使っても（1）の双対問題：

ans＝55．652

最小化 bTy

ans＝46′207

条 件 ATy≧c，y≧0，

k＝2613の基底解が−rCで与えられることを思い出せば，線ランダムな問題を各サイズたった1題解いただけで形計画問題の最適性条件：（i）主実行可能性，（ii）双対実は何とも評価しがたいが，とにかく関数simplex（）を行可能性，仙双対性（あるいは相補スラック性）から使って300変数500制約式の問題が50秒前後で解け

simp（ex（）の出力の正否をチェックすることができる．

た．商用コードに比べれば確かに見劣るかもしれない

その場合，SimpLex（）の出力に被約費用rcも加えてが，FortranやCを使ってこれだけの性能を引き出function［vaL，SOl，k，rC］＝Simplex（A，b，C）

すにはとても40行足らずのプログラムでは無理だろとしておくとよいだろう．

つ．

さて，動作確認で障害が発見されなければ次は関数

だいぶオタクな内容になってしまったが，今回の冒

simplex（）の性能を評価する段階だ．本格的な数値実頭，5節に述べた（a），（b）が宿題1を通して実感できた

験は読者に任せることにして，ここでは問題をランダのではないだろうか．Octaveによるプログラミング

ムに生成し，いくつか解いてみるだけに留めよう：

はFortranやCとは根本的に異質で，同じ感覚で行

octave：4〉A＝rand（200，100）−rand（200，100）；

うと効率のよいプログラムにならない．その一方，

octave：5〉b＝OneS（200，1）；c＝rand（1，100）；

FotranやCによるプログラミングに慣れていなくて

octave：6〉［tl，ul，Sl］＝CPutime（）；＼〉［val，SOl，k］＝Simplex（A，b，C）；＼〉［t2，u2，S2］＝CPutime（）；OCtaVe：7〉t2−tl，u2−ul，kも，線形代数に関する知識さえあれば，それがOctaveのプログラミングでは直に役立つのである．

8．大域的最適化に向けて

ans＝0．88600

さて，本稿の目標は改訂単体法のプログラム作成で

ansニ0．77600

k＝309

はなく，それを道具に難しい非凸計画問題を解決する

計算時間を計測するため，time（）に代ってcputime（）

ことである．宿題1はこれくらいにして先に進もう．

を用いているが，二つめの出力（ul，u2）がOctave

難しいといっても，使えるのは関数simplex（）だけ

自身の作業時間で，三つめ（sl，S2）はOSが代行しなので解ける問題も自ずと限られてくる．大域的最適た作業の時間，両者の和が最初の出力（tl，t2）で化問題としては初級の分離可能凹最小化間庵（sepa−

ある．したがって∪2−ul＝0．77600秒が，この100

rableconcaveminimizationproblem）［2，7］変数200制約式の問題を解くのにOctaveがsimpLex

最小化／（Ⅹ）＝∑失1ノ；（JJ）＋∑プ＝打1CJJJ（）で賛した正味の計算時間になる．もう2題ほど解い

条 件 Ax≦b，Ⅹ≧0

（4）

てみよう：

に的を絞り，その大域的最適解を求める矩形分枝限定

octave：8〉A＝rand（300，200）−rand（300，200）；法（rectangular branch−and−bound algorithm）を

octave：9〉b＝OneS（300，1）；c＝rand（1，200）；

次回にOctaveでプログラムすることにしよう．問題

octave：10〉［tl，ul，Sl］＝CPutime（）：＼（4）の制約条件はすべて線形計画問題（1）と同一であるも〉［vaL，SOl，k］＝Simplex（A，b，C）；＼〉［t2，u2，S2］＝CPutime（）；OCtaVe：11〉t2−tl，u2−ul，k
ans＝9．9850
ans＝8．6190

2005年11月号のとし，さらに簡単のために実行可能領域
β＝（Ⅹ∈R乃IAx≦b，Ⅹ≧0）は有界で内点があることを仮定しよう．そうすると，
祝ノ＝maX（ェJl∬∈β），ノ＝1，…，久© 日本オペレーションズ・リサーチ学会. 無断複写・複製・転載を禁ず.

（27）丁61

図1問題（4）の2次元の例
はすべて正の有限値になるが，各ノ＝1，…，カに対して難しさの本質はもちろん凹関数ムにあり，その数

ムは［0，勘］を含む開区間ムで非線形，凹関数である♪が少なければ問題（4）は効率よく解けそうに思える．

ものとする．つまり，どのムに対してもも，わ∈ムな

この直感は半ば正しいが，計算の複雑さからいえば♪

らば＝1の場合も（4）はNP困難である［5］．図1は，乃＝2，￡［（1一人）む＋勅］≧（1一人）ム（む）＋人力（わ）が任意のス∈［0，1］で成り立つ．

♪＝1のときの関数′の例／（れ，∬2）＝−（∬1−5）2−10」r2

問題（4）の目的関数′は凹関数￡と線形関数∑左折1qみの和であるが，これも凹関数になる：

を等高線とともにOctaveで（厳密にいえば，Octave
から呼び出されたGnuplotで）プロットし，後からフリーのお絵書きソフトTgif［8］を使って凸多角形β′［（1一人）s＋加］を描き加えたものだ．したがって問題（4）のイメージを＝∑失1ム［（1一人）む＋耽】＋示しているが，こんなに小さな問題例にさえ点A，B∑左折1q［（1一人）む＋勅］≧（1一人）∑失1［亮（み）十∑プ＝舛1q5J］十のような複数の局所最適解が存在する．高次元の問題人∑某1［ム（わ）＋∑プ＝舛1CJん］ともなれば，その難しさは推して知るべしである．悲＝（1一人）／（s）＋バ／（t），∀ス∈［0，1］．

観的な締めくくりとなったが，次回はこの印象の払拭に努めることにしよう．

この不等式から

f［（1−））s＋）t］≧min（f（s），f（t）），∀）∈【0，1］，

参考文献が成り立つことも分かるが，これはsとtを結ぶ線分上において関数／がs，tのいずれかで最小値をとる［1］Chvatal，Ⅴ．，LinearP7WYamming，Freeman（1983）．

ことを意味している．したがって，S，tともβから［2］Horst，R．andH．Tuy，Globa10timization：

選んだ場合を頭の中に描けば，′は凸多面体βのいずれかの端点で少なくとも局所的に最小化されるはずである．大域的最適解は局所最適解の中にあるので，
このことは（4）を解くうえで非常に重要な手がかりである．白状すると，今のところ手がかりはそれがすべてといって過言でない．大域的な最適性を保証するには，

minねticAj坤7mChes，Springer（1993）．
［3］今野 浩，「線形計画法」，銅斗技連（1987）．

［4］Martin，R．K．，La7ge Scde Linear and（砂timization，Kluwer（1999）．

［5］pardalos，P．M．andS．A．Vavasis，

Quadraticpro・

grammlngwithonenegativeeigenvalueisNP−hard，♪〝用αJ〆Cわ∂αJ（初め期加粛州1（1991），15−22．

すべての局所最適解の値を比較するより手だてがない［6］田村明久，村松正和，「最適化法」，共立出版（2002）．

のである．ところが，局所最適解の候補であるβの［7］Tuy，H．，ConvexAna如isand Global物timizaton，

端点の数は次元乃の指数関数であり，それを単純に

Kluwer（1998）．
［8］http：／／bourbon．usc．edu：8001／tgif／

列挙していては数十変数の問題であっても現実的な計算時間のうちに大域的最適解へはたどり着けない．

丁62（28）© 日本オペレーションズ・リサーチ学会. 無断複写・複製・転載を禁ず.

オペレーションズ・リサーチ