MATLABクローンによる大域的最適化（3）−Octaveはここまでできる一
久野 営人…1州111州Il…‖＝‖＝＝＝＝＝＝‖‖＝‖‖＝＝‖＝‖＝‖‖＝＝＝＝‖‖‖‖州‖‖‖‖＝＝＝‖＝‖‖＝川‖＝‖‖州Il…皿＝＝州＝州‖川11州11州川州Il…冊＝＝ll…………ll州州＝1111111州＝州＝冊＝＝州＝1111111州＝州側＝＝子問題P（C烏）は元の問題（4）と本質的に同一の凹最小

9．矩形分枝限定法化問題であり，直接には解けない．そこで通常，以下のように再帰的な方法が用いられている：

Octave版改訂単体法simpIex（）で料理するのは，アルゴリズム2．

分離可能凹最小化問題［6，7］

ステップ0．ゑ：＝1，Cl：＝C，￡：＝（1），Ⅹ＊：＝0，Z＊：最小化／（Ⅹ）＝∑失1亮（み）＋∑デ＝打1C晶（4）

条 件 Ax≦b，Ⅹ≧0＝／（Ⅹ＊）とおく．

に決まった．前回，仮定したように制約条件は線形計

ステップ1［子問題選択］．￡＝8ならば終了（Ⅹ＊が画問題（1）と同じ，A∈R研×乃でb∈Rmは非負，簡単

最適解）．そうでなければ，￡から適当な添字オ烏．のために実行可能領域を取り出してC：＝C砧とおく．
ステップ2［限定操作］．子問題P（C）の下界値乏を

β＝（Ⅹ∈R乃】Ax≦b，Ⅹ≧0）計算する．／（Ⅹ烏）＜z＊を満たす実行可能解Ⅹ烏が得は有界で内点があるものとする．したがって，〟J＝maX（みlx∈β），ノ＝1，‥リA（5）

はすべて正の有限値であり，各ノ＝1，…，♪に対して￡
は［0，〟J］を含む開区間んで非線形，凹関数である．

られれば，Ⅹ＊：＝Ⅹ烏，Z＊：＝／（Ⅹりと更新する．2≧z＊ならば，Cを破棄してステップ1に戻る．

ステップ3［分枝操作］．Cを2つの矩形C2烏，料理のまえにその段取りとなる矩形分枝限定法（rec−C2烏＋1に分割し，（2ゑ，2々＋1）を￡に加える．ゑ：＝tangularbranchTand−boundalgorithm）について説

カ＋1としてステップ1に戻る．

明しておこう．拙稿［4］と重なる部分もあるが，おそ整数計画法の知識があれば，すぐにこのアルゴリズ

らくは忘却の彼方であろうから省略せずに解説しよう．

ムが混合整数計画問題

9．1整数計画法と同じ？

最小化 cTx

実行可能領域β上では，目的関数／に対して非線

条 件 Ax≦b，Ⅹ≧0，Ⅹク∈（0，1）p形に働く変数Ⅹp＝（れ，‥．，∬♪）Tが（5）で定まる矩形

C＝［0，仇］×…×［0，㍑♪］を解くための分枝限定法と同様のものであることに気づくだろう．分杖操作で（0，1）♪を分割すれば，そのの外にはみ出ることはない．したがって，問題（4）に条要素はアルゴリズムの反復とともに着実に減少し，有件Ⅹァ∈Cを追加して限回の反復後に1点となって（6）の子問題は簡単な線形計画問題に帰着する．ところが，凹最小化問題（4）では

小
p（C）l

b，Ⅹ≧。，Xp∈C

Cをいくら分割しても1点にはならず，P（C）が簡単を解いても（4）と同じ解が得られるはずだ．矩形分枝限な問題になることはない．その一方，（6）では（0，1）♪

定法は，int・で集合の内部を表すことにすると，が1点になるまで分割しないと実行可能解すら見つけ

C＝Uk。￡Ch，intChnintCL＝8，k≠l，ることは難しいが，以下で見るように（4）では限定操作が満たされるようにCを複数の矩形C烏（々∈￡）に分の下界値乏を計算する過程で，副産物として反復の割し，（4）の子問題P（C員）を次々に解いていく方法でたびに実行可能解を得ることができる．

ある．矩形Cを小さなC烏に差し替えたところで，
くの たかひと筑波大学大学院システム情報工学研究科〒305−8573つくば市天王台1−1−1
826（36）■

9．2 一般的な限定操作
矩形Cは々回めの反復で

C＝［sl，ム］×…×【ざ♪，′♪］となっているものとしよう．目的関数の各亮は区間© 日本オペレーションズ・リサーチ学会. 無断複写・複製・転載を禁ず.

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

戸（C）よりも弱い下界値しか得られない．しかし，矩［5J，わ］で凹関数であり，前回にも示した通り，
ム［（1一人）sJ＋砧］≧（1一人）ム（5J）＋Å以わ）（7）

形Cがどのように分割されようと戸（C）の実行可能が任意のス∈［0，1］で満たされる．このパラメータバ

領域は常に不変であり，目的関数をいくらか修正すれを使えば［5ノ，わ］上の点はJJ＝（1一人）5J十勅のようにば前の反復で解いたP（C）の最通解Ⅹト1を初期基底書けるので，これを使って（7）から人を消せば解としてsimpJex（）をすぐさま実行することができる．

血）≧禦…）十鮎）

反復のたびに線形計画問題（1）の利潤ベクトルcに関する感度分析を行うのと同じである．また，Simplex（）となる．この不等式の右辺をが生成する戸（C）の最通解Ⅹカは目標の問題（4）の実行可能解でもあるので，暫定値z＊に対して′（Ⅹ々）＜z＊∬ノーか鮎）め（∬ノ）＝を満たせば，暫定解Ⅹ＊の更新にも用し−ることができと定義すると，ムとアフィン関数めとの間にる．

仇（み）≦ム（∬J），∀み∈［ざノ，わ］（8）

利潤ベクトルC以外のデータが共通の線形計画問題の関係が成り立ち，JJ∈（sノ，わ）のときには仇（ェノ）＝を繰り返し解くことになるので，前回作ったsimplex

ム（み）となる．さらに，￡は非線形なので（）のように，すべての定数，変数が関数内でローカル

仇（JJ）＞ム（み），∀JJ∈ム＼［5J，ん］（9）を示すこともできる．不等式（8）から直ちにわかるのは，

少）1C㌦J子問題P（C）で各ムを仇に置き換えた問題
戸（C）lを解くと，その最適値としてP（C）の下界値が手に入ることだ．問題戸（C）は，防がアフィン関数なので線形計画問題である．早速，関数simp［ex（）を使ってみに定義されていてはいささか都合が悪い．そこで，行

列Aや右辺ベクトルb，基底B，非基底Nなど，
simplex（）が呼び出されるたびに初期化されては困る
ものを外部変数として別のところで定義することにし，
Cを引数ccとしてsimplex（）に渡せば，それに対する最適値valだけを返すように書き換えておこう：
functionval＝Simplex（cc）％％StepO％％

たくなるが，厄介なのは追加した条件Ⅹp∈Cの存在
globalZ8；g10balRF；g10baIk；

である．この条件は単なる有界条件

globa（m；globaIn；g10baIB；g10baIN；

SJ≦み≦ん，ノ＝1，…，久なので，有界単体法（bounding simplex algorithm）［1］を使えば非負条件を扱うのと変わらない手間で処理できるが，残念ながら関数simpIex（）にその機能は

g10baIx；g10baJAA；g10balBi；g10ba‡b；

OPt＝0；

Whileopt＝＝0

持たせていない．そのため，Simplex（）で処理するときには2♪個の不等式条件をAx≦bの中に組み込ん％％Stepl％％y＝CC（B）＊Bi；rC＝CC（N）一y＊AA（：，N）；

で処理せぎるをえない．これはいかにも馬鹿げている．

9．3 関数simplex（）で行う限定操作

endwhile

下界値の計算で問題戸（C）を解く方法を最初に提案

endfunction

したのはFalkとSoland［2］で1969年の大昔だが，例えば，「globalx」は解を収める外部変数xの使用

その5年後にSolandは戸（C）に有界条件Ⅹp∈Cが宣言であり，globaL宣言さえあればsimplex（）以外か

必ずしも必要のないことを示している［5］．つまり，らもxにアクセスできる．このわずかな修正で関数

戸（C）を簡略して線形計画問題（1）とまったく同じ構造

simplex（）は，戸（C）の再最適化を行って子問題P（C）のの下界値乏を生成するサブルーチンとしてアルゴリ

ズム2のプログラムから呼び出すことが可能になる．

小
戸（C）l

b，Ⅹ≧。

9．4 也l分割による分枝操作

を解けばよいのである．この間題P（C）の最通解Ⅹ烏節9．1で触れたように，離散最適化の場合と違ってをsimplex（）で求めれば，子問題P（C）の下界値は乏凹最小化問題（4）に対する分枝限定法では，束巨形Cを＝g（Ⅹ烏）で与えられる．問題戸（C）の制約条件を緩和

いくら分割しても一般にP（C）が自明な最適解をもつしたわけだから至極あたりまえだが，代償としてことはない．そのため，実際に使うときには分枝が無© 日本オペレーションズ・リサーチ学会. 無断複写・複製・転載を禁ず.

2005年12月号（37）82丁限に続くことのないように，ステップ2の打ち切り判

10．1再帰関数による実現

定之≧z＊にe＞0の誤差を許して之≧z＊−eとするこ

凹関数亮（み）＝−qJプに対して区間［も，わ】で下界値とが多い．しかし，e＝0であっても，ステップ2でを与えるアフィン関数めは，む≠わならば戸（C）を解くことに加え，ステップ3でのCの分割に，動（み）＝−q（も＋わ）み＋qsノわ

これから紹介するw分割（w−Subdivision）を用いれと書けるので，（cノ（sノ＋わ），qもげを2×♪行列Gのノば有限回の反複後に必ずアルゴリズムを終了させるこ列に格納する．また，矩形Cも区間［sJ，ん］の端点とができる．Soland［5］の功績は，戸（C）の簡略化な（sムげをノ列にもつ2×♪の行列Cで表すことにする．

どよりむしろ，これを示したところにある．

こうしておけば，ステップ3で［5。，J。］が更新されても，Cノを成分とするベクトルCを使って問題戸（C）の最適基底解Ⅹ々に対して，まず
q∈argmax（f，（xf）−gj（xf）lj＝1，…，b）

G（1，q）＝C（q）＊（C（1．q）＋C（2，q））；（10）

G（2，q）＝C（q）＊（C（1，q）＋C（2，q））；

を選び，∬吉を境に区間［5。，J9］を【ぶ9，ポ］とレ吉，′。］のようにP（C）の目的関数♂を修正できるし，の2つに分割する．そして，新たな矩形C2烏とC2打1

zb＝Sum（G（2，：））を次のように定義する：

C2烏＝［∫9，∬書］×口由【も，ち］C2糾1＝［J吉，J9】×口栖［5ノ，わ］．
）−Simplex（［G（1，：），C（p＋1：m＋n）］）；（11）

によって下界値乏も計算できる．問題（l鋸こらは乃個しかないが，C（∩＋1：m＋∩）はゼロ・ベクトルである．

区間h，〜。】は分割できないかもしれないが，∬綽（∫9，

問題戸（C）の最適解Ⅹhはsimplex（）が算出して外部′。）ならば（軌（9）によってg。（∬吉）≧ム（∬吉）が成り立つ

変数×に格納するので，これを使えばベクトルはずである．添字ヴの選択の仕方から，め（ガ）≧fx＝C（1：P）′．＊x（1：P）．＊x（1：P）；￡（ガ）がすべてのノに対して成り立ち，の第ノ成分が一力（ガ）である．したがって，Ⅹ々にお2＝g（Ⅹ烏）ける（1却の目的関数値は＝∑失1め（ガ）＋∑畏舛1qガ

tz＝−Sum（fxトc（p＋1：n）＊x（p＋1：n）；≧∑失1カ（ガ）＋∑プ＝舛1らガに等しく，また下界値zbが暫定値よりも小さいとき＝／（Ⅹ烏）≧z＊にはα分割の（10）をとなってCはステップ2で破棄されていなければな［v，q］＝maX（G（1，：）′．＊x（1：P）−G（2，：）′一fx）；

らない．この甜分割を各反復で実行すると矩形Cは，で計算すればよい．このときに必要なのはqだけだ

Ⅹ烏が基底解なので，凸多面体βの端点を通って∬1

が，X（q）で区間［5。，J。］を分割，更新し，（11）で定義さ軸から∬♪軸までのいずれかに垂直な平面でしか分割れる新たな矩形2つに以上の手続きを再び施す．

されない．端点の数も座標軸の数も有限であり，C便利な道具をさらにもう一つ紹介しよう．うれしいは有限回しか分割できないことになる．したがって，

ことに，Octaveでは関数の再帰呼び出しが可能であ

アルゴリズム2も有限回の反復後に終了する．結局，る．再帰を使えば，上の手続きの繰り返しがコンパク

誤差ビを使わなくて済み，終了時のⅩ＊に残されるのトな次の関数recursion（）で実現できる：

は凹最小化問題（4）の厳密な大域的最適解である．

functionrecursion（｛，q，W，G，C）

10．Octave版矩形分枝限定法
g10barZ8；grObaJm；gIobaln；g10balp；globalc；

問題（4）の数値的な解決は，その目的関数に含まれる

groba（x；gIobaliz；gJobalix；globarnum；

凹関数上の値が計算できないことには始まらない．

％％Step2％％

ここでは単純にcノ≧0に対し，2次関数＋＋num；

ム（み）＝−qガ，ノ＝1，．‥，AC（「，q）＝W；（l為によって力の値を与えることにするが，後から変更

G（1，q）＝C（q）＊（C（1，q）＋C（2，q））；

することは容易だ．まな板にのるのは凹2次計画問題

G（2，q）＝C（q）＊（C（1，q）＋C（2，q））；
zb＝Sum（G（2，：））一Simprex（［G（1，：），C（p＋1：m＋

最小化 −∑失1CJポー∑左折1qみヨ門監

条 件 Ax≦b，Ⅹ≧0n）］）；である．さあ，料理に取り掛かることにしよう．

fx＝C（1：P）′．＊x（1：P）．＊x（1：P）；© 日本オペレーションズ・リサーチ学会. 無断複写・複製・転載を禁ず.

828（38）

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

tz＝−Sum（fx）−C（p＋1：n）＊x（p＋1：n）；

Z8＝1．Oe−8；RF＝512；k＝0；

if（tz＜iz）［m，n］＝Size（A）；N＝［1：n］；B＝［n＋1：n＋m］；

x＝［zeros（n，1）；rhs］；C＝ZerOS（1，n＋m）；

iz＝tZ；ix＝X（1：∩）；

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

endif

p＝CO］umns（nc）；C＝G＝ZerOS（2，P）；％％Step3％％％％StepO．2％％

if（iz−Zb＞Z8）［v，q］＝maX（G（1，：）′．＊x（1：P）−G（2，‥）′−fx）；

forj＝［1：P］

c（j）＝1．0；C（2，j）＝Simplex（c）；C（j）＝0．0；

w＝X（q）；％％Stepl％％endfo「G（1，：）＝nC．＊C（2，：）；C（1：P）＝nC；C（p＋1：n）＝lc；

recursion（2，q，W，G，C）；reCurSion（1，q，W，G，C）；％％Branch−and−bound％％endif

zb＝−Simplex（［G（1，：），C（p＋1：m＋n）］）；num＝P＋1；

endfunction

fx＝C（1：P）′．＊x（1：P）．＊x（1：P）；

アルゴリズム2のステップ1−3がどの部分に対応するかを示すのは，再帰的に記述したこの関数recur−

iz＝−Sum（fx）−lc＊x（p＋1：n）；ix＝X（1：n）；

sion（）では難しいが，敢えて記せば％で始まる注釈のif（iz−Zb＞Z8）［v，q］＝maX（G（1，：）′．＊x（1：P）−fx）；

ようになる．ステップ1の子問題選択は，再帰呼び出

W＝X（q）；

しで最も簡単に実装できる深さ優先則（depth−firstrecursion（2，q，W，G，C）；reCurSion（1，q，W，G，C）；

rule）を採用した．したがって，子問題の中から常に

endif

最も新しく生成されたものが選択される．変数ixはVa［＝iz；SO］＝ix；nSb＝num；PIV＝k；

暫定解Ⅹ＊，izはその目的関数値z＊＝／（Ⅹ＊）である．
下界値tzとizの比較による打ち切り判定には外部で定義する誤差Z8＝1．Oe−8を許しているが，これは丸め誤差のイタズラを防ぐためであり，理論的には必要ない．また，reCurSion（）の中で戸（C）の目的関数gと矩形Cの情報をそれぞれ収めた行列G，Cが再帰的に呼び出されるrecursion（）に引数として渡されるが，C言語と同じでG，Cそのものではなく，その値のコピーが渡されるに過ぎない．そのため，呼び出されたrecursion（）がコピーを書き換えても，呼出しから帰ったときに元のG，Cにまで変更がおよんでいる心配はない．

endfunction

関数rectBB（）は，問題（13）の行列A，ベクトルb，
それに目的関数の非線形部分の係数（cl，…，C♪）と線形部分の（cp．l，…，Cn）をそれぞれA，rhs，nC，lcとして渡すと，（13）の最適値val，最適解soL，関数simplex（）の呼び出し回数nsb，およびピボット演算の総数pivを返す．注釈のStepO．1で示した部分で関数simplex（）
やrecusion（）と共有する外部変数の初期化を行い，

StepO．2で初期矩形Cl＝［0，ul］×…×［0，uP］を定めている．鉄則に反して，ちゃっかりforループを使っているが，できるだけタイトな矩形を求めるには（5）に相当する♪題の線形計画問題を解かざるをえない．

10．2 関数rectBB（）関数recursion（）は，戸（C）の目的関数や矩形Cに関する情報が与えられたときにステップト3を繰り返
すだけで単独では何もできない．そこで，外部変数に必要な値を与えてrecursion（）を起動させるプログラ

ムが必要になる．それが，次に示す関数rectBB（）

注釈でBranch−and−boundと書いたところから下がアルゴリズム2のステップ1−3になるが，この部分は関数recursion（）と重複する．これは，最初の子問題

P（Cl）をrectBB（）の中で処理しているためだが，P（Cl）の最適解Ⅹ1と／（Ⅹ1）をそれぞれ暫定解ix，暫定値izに代入してから再帰関数recursion（）が呼び出さだ：

function［val，SOl，nSb，Piv］＝reCtBB（A，rhs，nC，lc）％％StepO．1％％
globalZ8；globalRF；globalk；globalm；
g10baln；globalB；globalN；globalx；
globaIAA；globalBi；g10balb；g10balc；
g10balp；g10baliz；globalix；globalnum；

れる．
10．3 性能評価と考察
線形計画問題と違ってプログラムの動作確認に適当な非凸計画問題のテスト問題を捜し出すのはなかなか容易でないが，（1カで定義した凹関数ムに線形項を足して力（J）＝−Cノエプ＋めJノとすれば，問題集［3］にいく
つか例題を見つけることができる．それらを解く場合© 日本オペレーションズ・リサーチ学会. 無断複写・複製・転載を禁ず.

2005年12月号（39）829は，reCtBB（），reCurSion（）で目的関数値の計算に若ans＝4549

干の修正が必要だ．紙幅が少なくなってきたので，動0CtaVe：11〉［tl，ul，Sl］＝CPutime（）；＼

作確認は省略し，（13）のランダムに生成した問題例だけ〉［vaJ，SOr，nSb，Piv］＝reCtBB（A，b，C（1：15），C（16：100））；＼を関数rectBB（）で解いてみることにしよう．
まずは，小手調べに（桝，乃，♪）＝（50，30，10）の問題〉［t2，u2，S2］＝CPutime（）；OCtaVe：12〉t2−tl，nSb，Piv
ans＝27．457

を解いてみると
OCtaVe＝1〉A＝rand（50，30）−rand（50，30）；
OCtaVe＝2〉b＝OneS（50，1）；C＝rand（1，30）；

ans＝2972
ans＝12507
0CtaVe：13〉［tl，ul，Sl］＝CPutime（）；＼OCtaVe：3〉［tl，ul，Sl］＝CPutime（）；＼〉［vaI，SOl，nSb，Piv］＝reCtBB（A，b，C（1：10），C（11：30））；＼〉［val，SOl，nSb，Piv］＝rectBB（A，b，C（1：20），C（21：100））；＼〉［t2，u2，S2］＝CPutime（）；〉［t2，u2，S2］＝cputime（）；OCtaVe：14〉t2−tl，nSb，PIV

OCtaVe：4〉t2−tl，nSb，Piv

ans＝198．45

ans＝0．46800

ans＝291

ans＝39413

ans＝478

ans＝75851になり，0．47秒の間に関数simplex（）が291即平び出

非線形項上の数♪の増加とともにCPU時間も加速的されて計478回のピボット演算を行ったことがわかる．
同じデータで（∽，乃，♪）＝（50，30，20）としても，に増加するが，30年も前のアルゴリズムがベースのプログラムとしては大健闘ではないだろうか．実のとOCtaVe：5〉［tl，ul，Sl］＝CPutime（）；＼ころ，Simplex（）と違って関数rectBB（）もrecurtion〉［val，SOl，nSb，Piv］＝reCtBB（A，b，C（1：20），C（21：30））；＼（）も計算効率より記述の簡潔さを優先してある．便利〉［t2，u2，S2］＝CPutime（）；

な再帰呼び出しも，局所変数の生成・消去を繰り返す

OCtaVe：6〉t2−tl，nSb，PrV

ばかりで無駄が多く，本当に効率のよいプログラムを

ans＝8．3240

書きたいのであれば使うべきではない．アルゴリズム

ans＝9537

2の添字集合￡をスタックとして自前で管理すれば，

ans＝3987

プログラムのCPU時間は劇的に改善するはずだ．

のように8秒強で解けた．注目すべきは，関数sim・Plex（）の呼び出し回数よりもピボット演算の回数の方

11．おわりにがずっと少ない点だ．小規模であるとはいえ，9千も第1回に約束したのでOctaveそのものの高速化にの線形計画問題戸（C）をいちいち端から解いていては，ついても記しておこう．Octaveの心臓とi）いえる2

とても8秒やそこらで（1湖の大域的最適解にたどり着け

つの線形代数ライブラリーBLAS（Basic Linear

そうにない．しかし，単体法の感度分析は各戸（C）を

AlgebraSubprograms）とLAPACK（LinearAlge−

再最適化するのに平均して0．5回に満たないピボット

bra PACKage）をチューン・アップする．まずはイ演算しか必要としなかったのである．

ンターネットで数値計算ライブラリーの宝庫Netlib

もっと大きな問題でも大丈夫そうなので，思い切っ［8］にアクセスし，ATLAS（Automatically Tunedて行列Aのサイズを（研，乃）＝（150，100）にして♪＝LinearAlgebraSoftware）の最新安定版atlas3．6．0

10，15，20の場合を解いてみたところ，以下のようなをダウンロードしよう．圧縮されたソース

結果となった：atlas3．6．0．tgzを解凍するとATLASという名のディ

OCtaVe：7〉A＝rand（150，100）−rand（150，100）；

レクトリが生成され，その下にコンパイルに必要なフOCtaVe：8〉b＝OneS（150，1）；C＝rand（1，100）；

ァイルが展開される．では，ディレクトリATLASに

OCtaVe：9〉［tl，ul，Sl］＝CPutime（）；＼

降りてコンパイル作業を始めよう．

〉［val，SOJ，nSb，Piv］＝rectBB（A，b，C（1：10），C（11：100））；＼〉［t2，u2，S2］＝CPutime（）；OCtaVe：10〉t2−tl，nSb，Piv
ans＝9．0040

最初に計算環境に合わせた設定ファイルを作るため，OSのコマンドラインから‡makeconfigと入力する．ターミナルに出力されるいくつかの質問

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

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

に答えなければならないが，よほど特殊な環境でない‡rpmbuild−−rebuildoctave−2．1．57．src．rpm

かぎり，自動検出されるデフォルトに設定しておけばとやれば，ATLASが自動的に組み込まれたバイナリ

大丈夫．この作業の終りに，例えば‖makeinstaJJー・バッケ，

arch＝Linux＿P4SSE2 のように次に行うべき作業のる．これをインストールすれば作業終了である．

ジoctave−2．1．57．i386．rpmができあがチューン・アップ前のOctaveで1，000変数の連立指示が出力される．この場合，OSとしてはLinux，アーキテクチャとしてPentium4が選択されたこと1次方程式Ax＝bを解いたところ，を意味するが，指示にしたがって

octave：15〉rand（‖seed， ，2005）；b＝rand（1000，1）；‡makeinstallarch＝Linux＿P4SSE2

octave：16〉A＝rand（1000，1000）−rand（1000，1000）；と入力すれば，あとはチューン・アップされたオブジ

octave：17〉t＝tjme（）；A＼b；time（）−tェクト・ファイル＊．0のアーカイブができるのを待

ans＝1．2715てばよい．この工程はかなりの時間がかかるが，
ATLASinsta complete‖が出力されれば終了である．
アーカイフuは，現在のデイレクトリATLASの中のIibの下，この例ではデイレクト1）Linux＿P4SSE2に生成であった結果が，チューン・アップ後にはoctave：1〉rand（ seed＝，2005）；b＝rand（1000，1）；
octave：2〉A＝rand（1000，1000）−rand（1000，1000）；octave：3〉t＝time（）；A＼b；time（）−tans＝0．83288

される．
ディレクトリLinux P4SSE2には7つのファイルがあるが，＊．aファイルの5つがアーカイフて、ある．
これらをOSの共有ライブラリが置かれているデイレクト））／usr／Jibか／usr／JocaJ／Jibなどにコピーすればよいが，アーカイブの1つIiblapack．aと同じ名前のファイルが／usr／他にも存在するので注意しなければならない．新しく生成された方の眈厄pack．aには元もと存在する／usr／lib／liblapack．aを補完するオブジのようになった．1行目の「rand（＝seed，，，2005）」は，データを同一にする乱数の種の指定である．どうです，
ますますOctaveを始めてみたくなったでしょ？
謝辞 この連載で行った研究の一部は日本学術振興会科研曹，基盤研究（C）17560050の補肋を受けています．また，連載の機会を与えてくださった編集委員の成蹟大学池上先生には，原稿に関する貴重なコメン
ト，アドバイスの数々に心から感謝いたします．

ェクト・ファイルが納められており，両者をマージして1つにする必要がある．それには，まず

参考文献［1］Chvまtal，Ⅴ．，LinearP7て智7Bmmin＆Freeman（1983）．

＄mkdirtmp；Cdtmp；arXJLiblapack．a；Cd．．

［2］Falk，］．E．and R．M．Soland，とやってLinux＿P4SSE2の下に一時的なデイレクト1）

Separable nonconvex programming problems

tmpを作り，その中でアーカイブ1iblapack．aを展開する．次に，，Man・dge∽e〃′Scgβ乃Cg15（1969），550−569．

［3］Floudas，A．，et al．，Handbook qf7bst P7Vblemsin＄mv／usr／［ib／liblapack．a．；arr［iblapack．atmp／＊．0のように／usr／（ibの下のIiblapack．aをLinux−P4SSE2

Localand Global（わtimization，Kluwer（1999）．
［4］久野善人， 非凸計画問題≠解けない問題一分枝限定法
による大域的最適化 ，，オペレーションズ・リサーチ44に移動させ，このアーカイブにtmpで展開したオブ（1999），232−236．

ジェクト・ファイル＊．0を追加する．この作業が終わったら，Linux＿P4SSE2の下にある5つのアーカイ［5］soland，R．M．， Optimalfacilitylocationwithcon−

ブ1ibatJas．a．Jibcblas．a，ljbf77blas．a，lib）apack．a，Iibt−

stat］as．aを／usr／Libか／usr／10Cal／［ibにコピMする．

あとは，Octaveを再コンパイルするだけだ．バイナ
リーではなく，ソース・パッケー

An algorithm for

ジの方のoctave−

2．1．57．src．rpmを，例えば／usr／src／redhat／SRPMのCaVeCOStS ，（砂e7tltions Research22（1974），373−382．

［6］Horst，R andH．Tuy，Global（砂timization：Deter−
minわtic App7VaChes，Springer（1993）．

［7］Tuy，H．，ConuexAnab｝Sis and Global（わtimization，Kluwer（1998）．
［8］http：／／www．netlib．org／

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

2005年12月号（41）83l