数理計画の理論と実装

非線形計画法アルゴリズムの実装と応用
田辺 隆八
州＝m‖‖…ll………………ll…川………州川…州川Il……ll………………＝…………………………lll………………州＝111……………tl……川州‖州…ll………ll……………l州Il……川仙‖ll………l州Il…

「ピボッティング」と呼ばれる操作によって入れ換え
1．はじめにながら解を探索する．有効制約法による求解はactiveな制約のよい見積もりが得られる場合には高最適化問題：
最小化：／（∬）

速化されるため，有効制約法をSQPの求解に用いる

制 約：g（∬）＝0，J≧0 ∬∈R乃，g（∬）∈R∽ （1）場合には，あるSQPの反復で解くQPに対するにおいて，目的関数／（∬），制約式g（∬）が一般の非線activeな制約を，次の反復で解くQPのactiveな制形関数であるような問題を非線形計画（NLP）問題

約の見積もりとし，全体の求解の高速化を試みる．こと呼ぶ．本稿ではこの求解アルゴリズムの実装とそのの事情は混合整数線形計画問題における分杖限定法の応用分野について具体的に解説する．まずアルゴリズ

解法に単体法が利用されるのと類似している．

ムと実装を概括し，線形計画法（LP）の実装技術が有効制約法は線形計画（LP）問題：

NLPに継承されていることを示す．次に内点法の実最小化：cfェ

装においてもっとも重要な一次方程式解法部分の実装

制 約：A∬＝占，J≧0に着目，LPの場合との類似点および相違点について述べる．続いてNLPの応用に共通する話題として，（3）に対する単体法の拡張であり，単体法の実装技術である：

モデリングおよび自動微分について解説する．最後に初期基底の作成（crash）

現実の問題においてNLPが利用される分野を具体的実行可能解の生成（phase−I）に解言見する．

基底の分解・更新

プライシング

2．アルゴリズムと実装などのノウハウをそのまま活かすことができる．

NLP求解アルゴリズムで現在主流なものとしてはすなわち，単体法⇒有効制約法⇒SQP

逐次二次最適化法（SequentialQuadratic Program−ming：SQP）と内点法の二つが知られている．SQPという形で実装技術の継承が行われているといえる．の反復は現在の反復点のまわりで目的関数・制約式を実際の数値実験においても単体法と有効制約法のパフそれぞれ二次・一次関数で近似した二次計画（QP）

ォーマンス上の差は比較的小さい．大規模QPを解く

問題：

ことの多い金融分野の例を見てみよう．
次はポートフォリオ選択問題においてリスク尺度最小化：∬亡（わ＋cf∬

制 約：A∬＝∂，∬≧0（2）（ポートフォリオの評佃基準）を分散で行った場合とを解くことに帰着される．したがってSQPの各反復絶対偏差［2］で行った場合の実行速度の比較である．

では性質の似通ったQPを連続して解くことになるが，前者は目的関数が二次関数となるのでQPに，後者はその場合の解法としては厳密解法である有効制約法が絶対値となるのでLPになるが，変数5の絶対値の最用いられることが多い．有効制約法は解において小化を定式化するにあたってactiveな（等式として満たされる）制約の集合をたなべ たかひと
㈱数理システム〒160−0022新宿区新宿2−4−3
2003年11月号

最小化l5l⇒最小化sl＋s2，
ただしざ＝∫1＋∫2，Sl，ざ2≧0（4）という書き換えを行っているので変数の数は異なる．
表1の結果によると実行速度はピボソティングの数© 日本オペレーションズ・リサーチ学会. 無断複写・複製・転載を禁ず.

（41）833

表1ポートフォリオ問題（LP，QPの実行結果）

変数 制約 ピボット 計算時間

5000

5060

10000 10060

621110．25秒

620．91秒

621．86秒

012
l1060

2 97 61000

626262

銘柄数

LP＋単体法
ピボット 計算時間

轍QP＋有効制約法

453332
303

0．49秒
2．09秒
3．60秒

（利用マシンPentiuml．5GHz＋1Gバイトメモリ ソフトウエア：NUOPT5．2．1，以下の計算例についてすべて同一）が少ないQPの方がむしろ高速という結果になっていR椚×りG∈R乃×乃はラグランジる．これは類似した性質を持つLP・QP問題ならば，−Jfz）のヘッセ行列あるいはその近似，行列β∈QPであるがために特段実行速度が低下することはな

R乃×乃は主・双対変数の値から構成される通常要素値いということを示している．このことから，例えばプが正の対角行列である．

ラント運転計画や．電力潮流計算［12］などにおいて，

ュ関数（上…／（∬）−g亡〝

表2は様々な数理計画問題に対する内点法の計算例

LPモデルをQP化して精度を向上するというアプロである．QP，NLPでも規模が小さければ計算時間はーチが現実的であることがわかる．

短いことから計算時間と種別の間に目立った相関はな

もう一つの流れである内点法によるNLP解法についこと，計算時間は変数や制約式の数，すなわち式（6）いても同様の実装技術の継承が見られる．のサイズによることがわかる．

内点法は最適任の必要条件であるKKT条件：

次の項では式（6）の解法に着目し，LPのケースにおいて培われた実装技術がどのようにQP，NLPの解g（J）＝0，法に継承されているかを見る．

∇／−A亡〝一之＝0，（5）

ズ2お＝0，

3．内点法の一次方程式解法∬≧0，Z≧0の後半の2式をまず，式（6）の構造の特殊性に着目した効率化が考え

ズ2お＝仏られる．LPの場合にはCが零であることからこの行

J＞0，Z＞0

列の左上部分は対角行列となる．その場合には一次方のように変形し（〃はパラメータ），非線形方程式系

程式の左辺をと見て解く手法である．主変数J（双対変数z）のみ⇒トAβ−1A亡］を解く主（双対）内点法，両方を同時に解く主双対内点法，大域的収束性の保障のためにメリット関数を導と変換することができる［3］．この変形は特にAの行入して直線探索・信頼領域を用いる方法，局所的収束の数（制約式の数）が少ない場合の行列サイズの削減性の向上のためにメリット関数の非単調ステップを許による高速化や，行列が正定値であることを利用したす方法，KKT条件のより高次の情報を利用する方法

行列分解手法が使えるなどのメリットがある．さらになど，様々なバリエーションが知られており，その適LPの場合にはAが変化しないことを利用してデー用範囲（NLP，QP，LP）も異なる［13］．ただ，いタ構造の効率化が可能である．反面，Aが非零要素ずれも非線形方程式系を解くというアプローチにおいを多く含む列（dense column）を持つとき，この変て共通している（例えばLPでもKKT条件は非線形

形は行列の疎行列性を破壊する（変換後の行列が密行方程式系となる）ため，次の構造の対称不定値行列を列になる）結果を招くので，式（6）の行・列を並べ替え左辺とする一次方程式解法に大部分の計算時間が消費てから式（7）の変形を部分的に行う方法であるされる．具体的には1万変数を超えるサイズのQP問AugumentedSystemApproach［4］が提案された．

題で一次方程式解法は全求解時間の約75％以上の時

β50 −A三0 上）d 一月昌間を占めるという実験例［14］がある．

一Aβ −Ad［℃￣ごf］O（6）（8）

ここで，Aは制約式のヤコピ行列（A…∇g（∬）∈
834（42）（7）［ヱ￣ご亡］© 日本オペレーションズ・リサーチ学会. 無断複写・複製・転載を禁ず.

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

表2 内点法の実行例
種別

変数

制約式

反復

計算時間

資源計画

LP

76308

28129

4437秒

ポートフォリオI

QP

6273

7874

483秒

ポートフォリオⅠI

NLP

123

267

200．5秒

原料購入計画

NLP

4166

2172

20103秒問題名

ここで，AdはAの各列のうち，非零要素の多いものの対称性を保つことができる点に特徴がある．

内点法ではその反復（20〜50回前後）のそれぞれのみから構成された行列，βdはβのうちでAdに対応する要素を集めた対角行列である．A5，βざはそれにおいて式（9）を解く必要が生じるが，ピボット選択をぞれ残りの列に対応する．

常に行うと計算量が増大するので，数値的条件の悪化

この方法はC≠0であるQP，NLPの場合に拡張が見られないかぎりピボット選択順序を再利用して計可能である．もっとも単純なケースとして凸な非線形算効率を向上させる方法が考えられる．式（9）に示され計画問題でCが対角行列であるケースでは，G＋βている内点法の各反復によって解かれる行列は単体法が対角要素で各要素が正なので式（7）あるいは式（8）の変と違い，常に構造（非零要素の場所）が同一で，数値形をそのまま使うことができる．Cが一般の構造を的な性質も急激には変化しないのでこの方策はおおむ持つ場合にも，Gのうち，0または対角行列となり，ね成功し，実際に分解が破綻してピボット選択のやりかつ対応するAの列が密でない部分をGざ，それ以外なおしが生じるのは内点法の全反復中，おおむね1回の部分をCβとするなら式（8）と同様に以下で済むことが実験的に確かめられている．

Gぶ＋βざ 0

このように内点法の大部分の反復ではピボット順を−A三固定して行列を解くことを繰り返している場合には，O Gβ＋βd −A良
一Aざ−Adさらなる高速化の工夫が可能である．具体的には，行O列の構造とピボット順に依存した前処理を行い，行列⇒［竺芸d＿。ざ（5，−1。ぶ ］（9）の分解演算手順の減少や処理系（CPU，メモリ）のなる変形が可能となる．実際のNLPには線形制約に最適化機能の利用促進を図る工夫が重要となる．前者

しか現れていない変数（Cの対応する行，列が0），の一つが，オーダリング（列，行の再番号付け）でああるいは別の変数とのクロスタームが存在しない変数り，後者がSupernodal法として知られている．これ（Gの対応する行，列の非対角要素が0）も多く現れらのノウハウは構造解析，流体解析，半導体のデバイるので，式（9）の変形によって行列の構造や性質がほとス・プロセスシミュレーションの文脈で蓄積されたもんどLPと同一となる場合も多い．すなわち，この方のが多い［10］．特にオーダリングの効果は劇的で，古法を取ることにより，LPと性質が似通った問題は，くから知られるminimum−degreeアルゴリズム［20］LPと同様のパフォーマンスの享受が可能となる．たを用いるのみでも一次方程式の解法時間が数十分の一だ，一般には変形の結果現れた式（9）の右辺の正志値性に短縮できることで疎行列の分解には必須の技術となが保証できないことに注意する必要がある．LPの場っている．オーダリングについてはいくつかのヒュー合にはAがrow−fullrankであるなどの適当な仮定

1）ステイタスが提案されているが，最近はgraph−の下で式（8）の右辺行列の分解時に零のピボットが現れpartitioningを用いた方法によって大規模問題に対すない（分解が破綻しない）ことを保証できるが［4］，るパフォーマンスの向上が報告されている［16］．

その保証が成り立たないNLPの場合には常にピボッSupernodal法は疎行列の演算を密行列演算に帰着ト選択付きの例えばBunch−Parlett分解［15］を行うさせるものである．そのため，密行列演算を処理系の必要がある．Bunch−Parlett分解は（2×2）ピボット最適化機能を活かして高速に実行するライブラリ（例を用いることによって，ピボット選択を行っても行列えばBLASカーネル［17］）の利用が前提となる．以

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

（43）835

表3 Supernodal法の効果

制約 SupernOdal利用

変数

問題種別

SupernOdalなし

LP（大規模ポートフォリオ）

12230

6072

75秒

150秒

NLP（非線形ネットワーク）

112769

44757

2102秒

8828秒

前，そのような実装は処理系のベンダによってのみ開発され，有料で配布されていたため汎用ソフトへの組み込みは難しかったが，近年には高品質で安価なBLASカーネルがWindowsを含めた多くの環境においてフリーで提供されるようになり，マシンの高速化と相まって，WindowsPCなどの一般的な環境でも高パフォーマンスな計算が可能となっている．例えば単純なコードで1200×1700の対称密行列の掛け算を行うのに，CPUがPentiumIV，1．9GHzのWindowsマシンでは22秒を所要するが，BLASカーネルの実図1計算グラフの実際例

装の一つ，ATLAS［6］を利用すると，所要時間は1．7秒に高速化される．これは1Gaopsを超えるパフォーマンスである．密行列ライブラリの性能向上によっノード間の依存関係をグラフの有向アークで表現したて，あえて疎行列性を無視して帰着できる密行列演算ものである．図1は次の二つの式を定義した計算グラ量を増やすことにより，パフォーマンス向上が観測さフの概念図である．

れる例もある．表3は大規模なLP，NLP問題につい／＝（J＋〝）・y・ZてSupernodal法の効果を実測したものである．

g＝eXp（∬＋y）

Supernodal法の効果の度合いには差があり，問題がこの形式を利用するとNLP問題全体は一つの計算大規模で疎性が低い（密行列に近い）ほど有効に機能グラフとして表現される．NLPの複雑さは問題を記する．大規模LPのベンチマーク問題から得られた結述する計算グラフのノード数として計測することがで果によるとSupernodal法の導入によって平均してきる．計算グラフのサイズによる複雑さの計測におい30％程度の速度向上が得られることが実証されている．ては同一の関数の反復構造を除外することができることに注意されたい．例えば上の例で（∬＋〝）に対応す

4．モデリングるコンポーネントが／，g両方の表現に利用されていNLPの計算機上の作成や入出力について考えてみる．

る．式の次数が限られているLP，QPの場合は問題制約式や目的関数（以下総称して「関数」）からそ全体を制約式の係数行列やヘッセ行列という定数行列の依存しているノードを辿っていくと最後には変数にで記述することができるので問題の取り扱いは定数行行きつく．逆に変数からその寄与しているノードを辿列の作成や入出力に帰着する（例えば古くから知られるとその変数に依存している関数に行きつく．自動微るMPS形式［21］はそのことを利用したLPの表現で分算法はこのような計算グラフの横断を行うことによある）．対して一般のNLPの場合には，式の次数やって任意の関数，任意の変数に対する徴係数を求める

関数形に制限がないため行列のみでは問題を表現する手法である．自動微分算法のうち，関数側からノードことができない．このような場合に行列に代わる式のを辿る方法をtop−down算法，逆をbottom−up算法

表現として計算グラフと呼ばれる新しい形式が知られと呼び，どちらが効率的かは変数と制約式の数の大小ている．計算グラフは式（中間変数）の計算の各ステに依存する［18］．一階微係数を求める場合，関数の数ップを初等的な演算（四則演算や初等関数）に分割，よりも変数の数が多い場合には前者が，逆の場合には各ステップによって計算される値をグラフのノード，後者が有利である．

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

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

計算グラフの作成は式の記述に等価である．大規模数の徴係数の計算に多くの時間を所要している．直線で複雑な計算グラフを手動で生成するのは一般に困難探索法は一階微係数のみを必要とするのに対して信頼なので，モデリング言語などの問題入力支援ツールが領域法は二階微係数をあわせて必要とする．二階微係

NLPでは特に重要となる．例えばあるプラント最適数の計算には多くの時間を所要するので，1回あたり化問題（NLP）では変数の数440，制約式の数294との反復では直線探索法が高速になっている．しかし，比較的小規模であるのにも関わらずモデル記述にはモ

信頼領域法は二階徴係数を用いているので，直線探索デリング言語の記述約10000行を必要とし，そこから法に比べて解法自体の収束が速く，全体の計算時間でノード数62574個の計算グラフが生成される．例えばは信頼領域法が有利となっている．表4の最後の行も汎用のモデリングツールであるSIMPLE［7］はモデル類似の格付け推移確率行列の期間構造の推定問題［19］に現れる繰り返し構造を簡潔に記述する機能を持ち，の例（計算グラフノード数65316個）であるが，類似LP，QPの記述に対しては行列を，NLPの記述に対の状況が発生している．最初に述べたプラント最適化しては計算グラフを生成，自動微分アルゴリズムによの問題と同様に上記の2例でも信頼領域法を適用したって非線形関数の二階までの微係数を提供する．

場合，自動微分演算に全体の90％以上を所要しておNLPでは，自動微分算法の時間が実際の計算に影り，総じて変数や制約式の数が比較的少ない複雑な非響する場合も多い．計算グラフが複雑な問題において線形最適化の高速化には自動微分の高速化が重要といは自動微分による徴係数計算の手間が無視できないこえる．と，さらに，LPやQPでは制約式の係数行列（ヤコモデリング言語や自動微分の実行効率は実装の方法ピ行列に相当）や目的関数のヘッセ行列は不変であるに大きく依存する．上記期間構造の推定問題においてのに対して，非線形関数の場合にはそれらは変数の値は，自動微分法の実装を変数が行列の形をしているとに依存して変化するので，アルゴリズムが変数値を更いう構造を活かしてチューニングしたところ，計算効新するたびに徴係数の計算が必要となることが理由と率は約40倍になったという報告［9］がある．一般に最して挙げられる．例えば先のプラント最適化問題では適化問題の記述には添字付けを用いた繰り返し構造が最適化全体に必要な時間130秒のうち，微係数の計算現れることから，計算グラフや自動微分算法そのものに全体の114秒（88％）の時間を消費している．にもこの構造を取り入れることによって，計算効率を微係数の計算手法は計算アルゴリズムの選択にも影向上させることが考えられる［11］．

響を及ぼす．次は信用リスク管理に必要な推移確率行

5．NLPの具体例

列推定問題で，与えられた非対称行列Oyg。rの52乗NLP（QP）を用いた定式化が有効な例としては例

根となる行列で特定の条件を満たすものを求める問題

えば次のようなものが挙げられる．

である［8］．
変 数：0∈R18×18

1．マルコフ過程による信用リスクモテル

最小化：l庵仰−0521】

2．債券ポートフォリオ最適化

制 約：∑0む＝1，0む≧0（10）ノ

表4の最初の行は上記の問題を二つのアルゴリズム（直線探索法・信頼領域法）で求解した場合の計算時

3．プラント解析モデル
4．化学反応系を含む生産計画モデル
5．電力システムにおける最適潮流計算間である．この間題の場合にも先のプラント最適化問

6．実験計画法の応用

題と同様に，変数・制約式の数に比較して計算グラフ

NLPモデルの非線形性はモデルの本質的構造に由のノード数が46901個と規模が大きなことから目的開来の場合もあるが，精度的な要求によるものも多い．

表4 格付け推移確率行列推定問題の結果
直線探索法

信頼領域法

問題

計算時間

反復

計算時間

52乗根

2125秒

781616秒

期間構造

875秒

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

（45）83丁モデルの本質的な構造によるものとしては，例えば

for theimplementation of Karmarkar s algorithm・1．，2．が挙げられ，1．はマルコフ過程を記述する確率

ORSAJ．onComput．，1（2）：84−106，1989．

推移行列が行列積の繰り返しとなることから，2．は金［4］Ⅰ．MarosandCs．Meszaros，Theroleoftheaugment−

利計算にlog／exp関数が現れることから非線形性がedsystemininteriorpointmethods，TechnicalReport
TR／06／95，BrunelUniversity，Department of Math−現れる．

一方，3．〜6．は後者の例で，全体はよく知られたLPの構造を持つが精度を追及する過程で非線形性の導入が必要となる．3．はプラントを構成する基幹部品の状態を記述する物理量の記述から，4．は化学反応（例えば混合・燃焼）の前後における物質，例えば原油などの化石燃料の体積や性状の変化から，5．は電力

ematicsandStatistics，London，1995．

［5］E．D．Andersen，）．Gondzio，Cs．Meszaros，andX．Xu，Implementation of interior point methods for large
scalelinearprograms，InT．Terlaky，editor，Interiorpoint methods of mathematicallprogramming，189−
252，KluwerAcademicPublishers，1996．

［6］R．C．Whaley，A．Petitet，andJ．J．Dongarra，

に対する燃料コスト関数から，6．は仮定するフィッテ

Automated Empirical Optimization of Software andィングの式の次数により，それぞれ非線形性が持ち込

the ATLAS project，Technicalreport，University of

まれる．このようなモデルは一般に大規模だが，よくTennessee，Knoxville，TN，Department of Computer

知られたLPの構造に非線形関数を「はめ込んだ」形

Science，Univ．ofTN，Knoxvi11e，TN37996，2000．となることが多い．このような問題を解くにあたって，
LP⇒NLPへの実装技術の継承は特に重要である．
NLPを実際に応用する場合に注意する点として，一般のNLPには解の一意性が保障できず，局所解が存在する場合があるため，通常の数理計画法ソフトウ
エアでは大域的な最適解が得られる保証がないという［7］山下浩，田辺隆人，逸見宣博，数理科学のためのモデリング言語SIMPLE，研究報告「数理モデル化と問題解決」アブストラクト，No．004−005（http：／／www．ipsj．or．
jp／members／SIGNotes／Jpn／23／1995／004／articleOO5・html）．

［8］楯岡成雄，青沼君明，中川秀敏，クレジット・リスクモ
デル，きんざい，2001．

ことが挙げられる．これは多くの数理計画法ソフトウ［9］青沼君明，田辺隆人，An Estimation for The Term

エアが，最適性に関して一次の必要条件のみしか考慮

StructureofYieldSpread，日本オペレーションズ・リしていないということに起因する．そのようなケース

サーチ学会，金融工学研究部会，2001．には実務上の要件から制約式を追加して実行可能領域を狭めて最適化を実行する，あるいは初期値を変化させて複数回の求解を繰り返すことになる．
モデルに含まれている関数の性質が既知の状態では問題の凸緩和と分枝限定法を利用することによって大域的最適解を求めることが可能な場合がある．ただし，その場合にはNLPの手続きをツール化して別の手続きから繰り返し呼ぶ手続きが必須である．その場合にはNLPの実装のソフトウエアとしてのインタフェースも重要になる．

［10］田辺隆人，対称行列の分解の高速化テクニカルレポ
ート，㈱数理システム，2001．
［11］田辺隆人，添字付け記法による自動微分高速化テク
ニカルレポート，㈱数理システム，2002．
［12］高橋，石原，石岡，小野，古塩，須藤，中村，「火力・揚水・水系水力の協調を考慮した需給運用計画機能の開発」，電気学会研究会資料，PE−02−130，2002．
［13］小島政和，土谷隆，水野眞治，矢部博，「内点法」，朝倉書店，2001．

［14］田辺隆人，「非線形最適化のアルゴリズムとソフトウエア」，最適化とアルゴリズム研究部会（SOA）発表資料，
1999．

参考文献

［15］A．M．ErismanandJ．K．Reid，Direct Methodsfor［1］H．Yabe，H．Yamashita，andT．Tanabe，Agloballysparsematrices，0ⅩfordUniversityPress，NewYork，and superlinearly convergent primal−dualinterior
pointtrustregionmethodforlargescaleconstrained

1986．

［16］E．Rothberg and B．Hendrickson，Sparse Matrix

optimization，TechnicalReport，MathematicalSys−

OrderingMethodsforInteriorPointLinearProgram−temsInc．，Tokyo，Japan，July1997（revisedJuly1998）・

ming，Manuscript，1996（http：〝www−unix．mcs．anl・gOV／［2］今野浩，理財工学Ⅰ，日科技連，1995．

otc／InteriorPoint／abstracts／Rothberg−Hendrickson．

［3］Ⅰ．Adler，N．Karmarkar，M．G．C．Resende，and G・html）．
Veiga，Data structuresandprogrammlngteChniques［17］J．J．Dongarra，J．DuCroz，I．S．Duff，andS・H
838（46）

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

marling，．Aなorithm 679：A set d Leve13 蝕ic
Linear Aなebm SubpY瑠mmS，ACM Trans．Math．

FinancialMarkets，2001．
［20］A．GeorgeandJ．W．H．Liu，771e euOlution〆theminimum d留ree aなorithm，SIAM Review，31：1−19，

Soft．，16，pp．18−28，1990．
［18］伊理正夫，久保田光一，高速自動微分法（1），（2）応用数

1989．

［21］Murtagh，Advanced Linear Programming，理，1，17−35，153−163，1991．

［19］K．AonumaandT．Tanabe，AnEstimationModel

McGraw−Hill，1981．

for Term Structure of Yield Spread ，Asia−Pacific2003年11月号© 日本オペレーションズ・リサーチ学会. 無断複写・複製・転載を禁ず.

（47）839