特集「最適化技術に基づく統計的推論」［原著論文］

統計数理（2013）
第 61 巻 第 1 号 47–78
c 2013 統計数理研究所


制約整数計画ソルバ SCIP の並列化
品野 勇治

1,4

2

1

1

・Tobias Achterberg ・Timo Berthold ・Stefan Heinz ・1

3

Thorsten Koch ・Stefan Vigerske ・Michael Winkler

1（受付 2012 年 7 月 2 日；改訂 9 月 6 日；採択 10 月 12 日）

要

旨

制約整数計画（CIP: Constraint Integer Programs）は，制約プログラミング（CP: Constraint
Programming），混合整数計画（MIP: Mixed Integer Programming），充足可能性問題（SAT: Satisfability Problem）の研究分野におけるモデリング技術と解法を統合している．その結果，制約整数計画は，広いクラスの最適化問題を扱うことができる．SCIP（Solving Constraint Integer
Programs）は，CIP を解くソルバとして実装され，Zuse Institute Berlin（ZIB）の研究者を中心として継続的に拡張が続けられている．本論文では，著者らによって開発された SCIP に対する
2 種類の並列化拡張を紹介する．一つは，複数計算ノード間で大規模に並列動作する ParaSCIPである．もう一つは，複数コアと共有メモリを持つ 1 台の計算機上で（スレッド）並列で動作す
る FiberSCIP である．ParaSCIP は，HLRN II スーパーコンピュータ上で，一つのインスタン
スを解くために最大 7,168 コアを利用した動作実績がある．また，統計数理研究所の Fujitsu
PRIMERGY RX200S5 上でも，最大 512 コアを利用した動作実績がある．統計数理研究所のFujitsu PRIMERGY RX200S5 上では，これまでに最適解が得られていなかった MIPLIB2010のインスタンスである dg012142 に最適解を与えた．
キーワード： 混合整数計画，制約整数計画，SCIP，並列化．
1. はじめにSCIP（Solving Constraint Integer Programs, http://scip.zib.de）は，制約整数計画（CIP: Constraint Integer Program）を解くために開発されたソフトウェア・フレームワークである（Achterberg, 2007, 2009）．CIP は，混合整数計画（MIP: Mixed Integer Programming，本論文では，混合整数線形計画: Mixed Integer Linear Programming を単に MIP と記述する）を完全に包含する最適化問題のクラスである（CIP の詳細は第 2 節で紹介する）．つまり，SCIP は MIP ソルバ
としても機能し，H. D. Mittelmann の WEB ページ（http://plato.asu.edu/bench.html）に示されているベンチマーク結果を見ると，ソースコードが公開されている MIP ソルバの中では，現在，最高性能を示すソルバである．本節では，まず MIP ソルバ開発の現状を紹介する．なお，
本論文では一般性を失うことなく，すべて最小化問題を対象として説明する．
ここでは，特定の問題専用ではなく，あらゆる MIP に対応できる汎用ソルバが対象である．過1
2
3
4

Zuse Institute Berlin, Takustr. 7, D-14195 Berlin-Dahlem, Germany
ILOG, IBM Deutschland GmbH, Ober-Eschbacher Str. 109, 61352 Bad Homburg v.d.H., Germany
GAMS Software GmbH, P.O. Box 40 59, 50216 Frechen, Germany
JST CREST：〒 102–0076 東京都千代田区五番町 7 K’s 五番町

48

統計数理

第 61 巻

第 1 号 2013

去 10 年間の MIP ソルバの性能向上は著しく，商用の MIP ソルバ，Gurobi（http://www.gurobi.
com/），CPLEX（http://www-01.ibm.com/software/integration/optimization/cplex-optimizer/），
Xpress（http://www.fico.com/en/Products/DMTools/xpress-overview/Pages/Xpress-Optimizer.
aspx）が高性能であることは広く知られるようになってきている．Koch et al.（2011）には，過去
10 年間の MIP ソルバの性能向上を示す結果が定量的に示されている．MIP ソルバは基本的には分枝限定法の実装である．緩和問題として線形計画（LP: Linear Programming）を解く分枝限
定法は LP ベースの分枝限定法とよばれる．現在の MIP ソルバは， LP を解いた結果から得られる情報に基づく洗練されたアルゴリズムを実行する．
ソフトウェア的には，LP ベースの分枝限定法内で構成要素となる，前処理，分枝戦略などの機能は十分に整理されてきている．各機能の概略は第 2 節で紹介する．構成要素となる各機能において，どのようなタイミングで，どのアルゴリズムを，どの程度の頻度，どの程度の論
理的な時間動作させるかは，大量の数値実験結果に基づいて決定されている．決定された内容
は，実行時に自動調整されるようにプログラム中に組み込まれることもあるが，一般的にはパラメタ値として実行時に指定できるように実装される．その理由は，適切なパラメタ値は，問
題のインスタンスに強く依存するためである．これまでに解けなかったインスタンスを解けるようにするため，それぞれの機能に新しいアルゴリズムが継続的に追加されている．新たに追加されたアルゴリズムが，あらゆるインスタンスに有効ということは稀であり，一般的には，そのアルゴリズムを動作させるべきかどうかを判定する手間を要することになり，これまでに解けていたインスタンスを解く時間はむしろ長くなる可能性もある．そこで，高性能な汎用ソルバを実現するためには，過去に高速に解けたインスタンスを含む，大量のインスタンスを用いた数値実験により，パラメタの再調整が必要となる．
一般的に，構成要素となる各機能に対して実装されるアルゴリズムのそれぞれがパラメタを持ち，SCIP を例にするとパラメタ数は 1000 を超える．この値は，商用ソルバよりも多いが，商用ソルバにも公開されていないパラメタが多数あるはずで，その実情は変わらないはずである（ユーザに対しては，全てのパラメタの設定には困難を伴うため，メタ・パラメタ設定が用意されるケースが多い．SCIP では，ヒューリスティックやカットなどに対して，「短時間の適用（fast）」，「できる限り適用（aggressive）」，「適用しない（none）」などを指定できる）．つまり，現在のソルバは，各機能毎に多くのパラメタを持ったヒューリスティックの塊が，大量のインス
タンスによる数値実験結果に基づいて高度にチューニングされた大規模ソフトウェアである．
実際，SCIP の現在のコード量は 450,000 行を超える． CPLEX は，コメント量は SCIP と比較して少ないにもかかわらず，コード量は 600,000 行を超える（Achterberg, 2011）．
ソルバの性能を限られたインスタンス数で効率的に評価する上で，適当なインスタンスを収集するということは極めて重要である．ランダムに生成されたデータに基づくインスタンス
は，極端に易しいか極端に困難なインスタンスとなることが多いため，現実問題，あるいは，具体的な組合せ最適化問題から生成されたインスタンスを商用ソルバ開発会社は持っており，
データ自体が資産であるとも言える．研究目的としても，MIPLIB（Bixby et al., 1992, 1998;
Achterberg et al., 2006; Koch et al., 2011）は，1991 年に公開されて以来，最新版の MIPLIB2010（http://miplib.zib.de）まで継続的に更新されてきた．このインスタンス・セットをまとめるという研究は，単なるデータの収集ではなく，ソルバ開発者の視点で何らかの意味を持つインスタンスが慎重に選択されライブラリに加えられてきている．その結果，極めて限られた数のイ
ンスタンス・セットにおいて，ソルバの性能が概観できるものとなっており，共通のインスタ
ンス・セットによってソルバの性能が評価されてきたことが，現在の高性能ソルバ開発につながっている．
本論文で扱う SCIP の並列化は，基本的には分枝限定法における木探索の並列化である．しか制約整数計画ソルバ SCIP の並列化

図 1.

49

Parallel search tree generated by ParaSCIP and FiberSCIP.

し，SCIP により解かれる問題のクラスは MIP を含んでおり，前述のような大規模ソフトウェ
アである MIP ソルバの性能を維持し，その各機能の性能を継承して動作しなければ，大規模な並列探索が実現されたとしても逐次処理ソルバの性能すら達成することは困難である．実際に，一般的な MIP ソルバの並列化に関する研究は過去にも行われてきている（Bixby et al., 1999;
Bussieck et al., 2009; Chen and Ferris, 1999; Eckstein, 1997; Phillips et al., 2006; Linderoth, 1998;
Ralphs et al., 2003; Shinano et al., 2003, 2007, 2008; Xu et al., 2009）が，スーパーコンピュータを利用しても，商用の逐次ソルバの性能を超えることはほとんどなかった．著者らの知る限り，
商用ソルバで解けないインスタンスを商用ソルバより先に解いたケースは，ほぼ手動で部分問
題を GRID 上で解いた Bussieck et al.（2009）のみである．しかし，商用ソルバ CPLEX の探索
木データを保存するファイルを解析して部分問題を生成し，分散計算させるもので，CPLEX
が木構造データの出力を中止した時点で利用不可となった． SCIP に限らず，MIP ソルバの並列化の困難さは Ralphs（2006）に整理されている．
2005 年以降，単独 CPU コアの性能は大きく向上しなくなり，複数コアを持つ CPU がデス
クトップ環境でも用いられるようになるにつれ，商用ソルバの並列化が活発になる．並列処理を適用できる環境がどこにでも存在するという環境の変化は大きく，商業的には環境の変化への追随が必然であったはずである．最も劇的に，絶妙のタイミングで開発されたのは Gurobi
である．それまでは，単独の CPU コアの性能向上が直接ソルバの性能向上に寄与したが，複数コアを利用する場合には，プログラム構造の変更を伴う．各社，既存のコードを拡張して並列ソルバを実現していたが，Gurobi は設計の当初から並列実行を前提に設計された．その基本的なコンセプトは，“Tree of Trees Concept” と呼ばれ，探索木を部分木の集まりとみなし，部
分木は既に高性能な MIP ソルバ自身に解かせるというコンセプトである．これは，著者らが試みてきた手法（Shinano et al., 2008）をソルバ内部に取り込んだコンセプトとみなすことができる．本論文で紹介する SCIP の並列化も同様の手法であり，並列探索木の様子は図 1 のようになる．図 1 において，影のついた各領域が部分木であり，部分木のルートとなる部分問題が独立して動作する MIP ソルバに転送され解かれる．
SCIP は，新しいアイデアをソルバ内部の機能毎に plug-in として柔軟に組み込める構造を持っている．実際，標準で提供されている MIP ソルバとしての機能ですら，plug-in として実現されているものが多く，構造的には簡単に置き換え可能な作りとなっている．しかし，内部に並列処理を組み込むことは当初の設計には含まれていなかった．スレッドセーフな設計にはなっていたが，実装はスレッドセーフでない部分も存在し，本研究を通して，そのような部分

統計数理

50

第 61 巻

第 1 号 2013

が検出され修正された．そのため，全ての商用ソルバが実現している共有メモリを利用した並列版さえ存在しなかった．本論文では，著者らが開発した SCIP の並列化拡張を紹介する．この並列化により，SCIP は共有メモリを利用した並列ソルバだけでなく，大規模分散メモリ環境上で動作するソルバとなり，最適解が知られていなかった MIPLIB のインスタンスのいくつかに初めて最適解を与えることができた（Shinano et al., 2012）．
次節以降の構成は次の通りである．第 2 節では CIP と SCIP を紹介する．第 3 節では SCIP の並列化を実現しているソフトウェア・フレームワーク UG（Ubiquity Generator）を紹介し， UG
上に開発された ParaSCIP と FiberSCIP を紹介する．また，UG の機能として実現される並列処理について説明する．第 4 節では数値実験の結果を示す．数値実験は，MIP を対象とし，共有メモリ上での実行結果と，大規模分散メモリ環境上での実行結果を紹介する．第 5 節では今後の課題等について述べる．
2. 制約整数計画（CIP）と SCIP

制約整数計画は次のように定義される最適化問題である．
定義 1.（制約整数計画（CIP: Constraint Integer Program））制約整数計画 CIP = (C, I, c) は，次の問題を対象とする:(CIP) c = min{c x | C(x), x ∈n, xj ∈

 for all j ∈ I}.

ここで，C = {C1 , . . . , Cm } は制約式 Ci : n → {0, 1}，i = 1, . . . , m からなる有限集合， I ⊆ N ={1, . . . , n} は変数の添字集合の部分集合，c ∈ n は目的関数の係数ベクトルである．加えて，
CIP は次の条件を満たさなければならない:(2.1)∀x̂I ∈ ∃(A , b ) : {x ∈I

ここで，C := N \ I ，A ∈





k×C

C

C

，b ∈

k

| C(x̂I , xC )} = {xC ∈であり，k ∈C

| A xC ≤ b }.

≥0 である．

条件（2.1）が要請しているのは，全ての整数変数の値が固定された後に残る部分問題は，常に LP になることである（注：SCIP が MINLP をサポートするように拡張された際に CIP の定義は拡張され，整数変数の値が固定された後に残る部分問題に対して最適解を求める手段があれば，LP である必要は無くなっている）．よって，整数変数だけが参照されるような制約式に対しては，二次式はもちろん，あらゆる非線形制約が許される． MIP は，全ての制約が線形，
つまり，C(x) = {x ∈ n | Ax ≤ b}，A ∈ m×n ，b ∈ m である，CIP の特殊ケースであり以下で定義される:(MIP) c = min{c x | x ∈ n , xj ∈ for all j ∈ I}.



この MIP の最適値に対する下界値は，整数条件を取り除いた緩和問題である次の LP の最適値により与えられる:(LP) c = min{c x | x ∈ n , Ax ≤ b}.
SCIP（Solving Constraint Integer Programs）は，CIP を解くために開発されたソフトウェア・
フレームワークである．CIP は問題記述能力の極めて高い最適化問題であり，現在の SCIP では，CIP として記述できる全ての最適化問題が解けるわけではない．しかし，SCIP は必要になった新たな制約式のクラスを扱うことができるように，必要に応じて plug-in を追加することで拡張が可能なソフトウェア・フレームワークである．その中心的な働きをするのは，制約式ハンドラ（constraint handler）であり，制約式 C(x) のクラスに対する plug-in を追加すること
で，扱える制約式のクラスが増える．この plug-in は，扱うクラスの制約式の意味を解釈し，そ制約整数計画ソルバ SCIP の並列化

51のクラスに属する制約式に対する操作を実装する．MIP の場合は，線形制約式ハンドラ（linear
constraint handler）が線形制約式を扱う．制約式ハンドラの主な機能は，解が与えられた際に，解いているインスタンス中で，その制約式ハンドラが扱う制約式のクラスに属する全ての制約が満足されているかどうかをチェックすることである．与えられたインスタンスに含まれる全ての制約式のクラスに対する制約式ハンドラが用意されれば，全ての整数変数が固定されるとLP になるため，この機能だけで列挙解法が構成できる．しかし，その解法は全列挙と同じであり，極めて遅い．そこで，SCIP の LP ベースの分枝限定法のフレームワーク側から必要なタイミングで，コールバック・ルーチンが呼び出される．例えば，LP 緩和問題を解いた後，制約式ハンドラのコールバック関数が呼び出される．ここでは，緩和問題の解を調べ，各変数の定義域の縮小や，カットとなる制約式を生成し追加することなどが実現できる．
SCIP により実現される LP ベースの分枝限定法の主たる処理の流れを図 2 に示す．構成要素となる主たる機能は次の通りである．
• Presolvers: 問題のサイズを縮小する前処理• Separators: 下界値を強化するカット生成• Propagators: 変数の上界値・下界値の強化（制約式ハンドラとは別途行われる部分）• Branching Rules: 分枝変数の選択方法の指定• Node Selectors: 分枝限定木における葉ノードの中で次に実行するものを指定• Primal Heuristics: より良い上界値を与えるためのヒューリスティック解法• Relaxation Handlers: 下界値を与える緩和問題の解法

これらの各機能単位に，plug-in として既に多くのアルゴリズムの実装が用意されており，独自に追加することも可能である．各構成要素となる機能は，図 2 に示される，SCIP が実行する
LP ベースの分枝限定法の中の事前に決められた箇所において，事前に決められた順番で呼び出される．ユーザは自由に plug-in を追加できるが，デバッグモードにおいて，提供されている各関数は呼び出し可能なタイミングであるかどうかのチェックが行なわれるため，ユーザが事前に決められた順序を無視したプログラムを書くことは禁止される．一つの構成要素となる機能に対して複数の plug-in が存在し，その呼び出し順はパラメタで制御できる．

図 2.

Flowchart of SCIP.

52

統計数理

第 61 巻

第 1 号 2013

上記の機能以外にも，特に SCIP において特徴的な機能として，Conflict Analysis や Pricer
がある．Conflict Analysis は，SAT ソルバで利用されてきた，探索過程で実行不可能となった部分問題から学習する仕組みを MIP へ拡張した機能であり，plug-in として拡張可能である．
Pricer は，列生成法の実装を可能とするもので，列生成時の変数の追加を可能とする． CIP，
SCIP に関するより詳細な内容に関しては，Achterberg（2007）を参照されたい．ここに述べていない plug-in の機能や plug-in の開発方法に関しては，SCIP の WEB ページ（http://scip.zib.de/）のドキュメントを参照されたい．
3. UG（Ubiquity Generator framework）を利用した並列ソルバ ParaSCIP と FiberSCIP
SCIP の並列化は，Tree of Trees Concept を特定のソルバに依存せず実現する単一のソフト
ウェア・フレームワーク Ubiquity Generator（UG）framework 上で実現している．本論文では，与えられた部分問題を MIP として解く，独立して動作するプログラム単位を単に Solver とよぶ．Solver は，一つのプロセッサ・コアに対して複数動作させることも可能であるが，通常は，プロセッサ・コアあたり，1 Solver を動作させる．実行時の Solver の動作形態は，通信がどのように実現されるかに応じて，スレッド，または，プロセスとして動作する．つまり，本論文で扱う並列ソルバの実体は，一つの問題を解くために，複数の Solver が同時に協調して動作するものである．
ここでは，まず，UG を紹介し，UG 上に開発された ParaSCIP と FiberSCIP を紹介する．UG
は，C++ で記述されたソフトウェア・フレームワークである．並列化対象となる逐次，あるいは，共有メモリ上でスレッド並列に動作する分枝限定法による最適化ソルバ（SCIP に限らない）を，分散メモリ並列処理環境で動作させることを意図して開発された．基本的には，並列化対象ソルバ自身を多数並列に動作させ，それらの間の動的負荷分散を実現する．通信処理関係の関数に関しては，UG 内部で再定義されているので，通信を扱うクラスの継承クラスを開発することで，他の動作環境への移植を容易にしている．同様に，ソルバに関しても，特定のソルバに依存しない作りとなっている．主たる設計思想は，最新の高性能最適化ソルバの性能を維持した並列ソルバの開発を容易にすることである．
一般に，分枝限定法を実行する並列ソルバは，次の三つの実行フェーズを持つ（Ralphs, 2006）．

（1） Ramp-Up フェーズ：計算開始から，全ての Solver が同時に動作可能となる十分な部分問題群が生成されるまでの期間（2） Primary フェーズ：動作している Solver 内に十分な部分問題群が存在して，全てのSolver において並列探索が実現されている期間（3） Ramp-Down フェーズ：計算終了が近づいてきたため，全ての Solver が十分な部分問題を生成できず，常に，少なくとも一つの Solver は，処理する部分問題が得られない状
態（遊休状態）が生じている期間
UG の主たる機能は，動的負荷分散であるが，その他に，複数の Ramp-Up 手法を実現する機能，および，チェックポイント・再スタート機能を，フレームワークの機能として提供する．
前述のように，通信手段は選択可能な構造となっているので，現在は，共有メモリ上でのマル
チ・スレッド並列による比較的小規模の並列処理を実現するために Pthreads ライブラリ，分散
メモリ環境上で大規模並列処理を実現するために MPI（Message Passing Interface）ライブラリが SCIP ソルバの並列動作に利用可能である．UG を利用して開発されたソルバの命名規則は次の通りである:ug[並列化対象ソルバ名，通信に利用するライブラリ名]．

制約整数計画ソルバ SCIP の並列化

図 3.

53

Initialization step.

この命名規則は並列ソルバの特徴を直接的に表現している．この命名規則を使うと， ParaSCIP
は ug[SCIP, MPI] であり，FiberSCIP は ug[SCIP, Pthreads] である．UG を利用することにより，論理的には全く同じ並列化手法を，異なる動作環境で実行することができる．このような実装
により，論理的な挙動と動作環境の影響を，より明確に詳細に調べる手段が提供できる．
3.1 UG が提供する並列処理
ここでは，UG が提供する機能により，どのような並列処理が実現されるのかについて，計算の開始から終了までの流れの中で説明する．ここで説明する並列処理は，プログラムの実行時形態としてプロセスとスレッドの違いはあるが，論理的には，ParaSCIP と FiberSCIP の両方で共通である．プログラム実行時には，ParaSCIP の場合 2 種類の MPI プロセス，FiberSCIP の場
合 2 種類のスレッドがシステム中に存在する．一つは，既に述べた部分問題を解く Solver であり，もう一つは，Solver 間の負荷分散を司り，プログラム全体を制御する LoadCoordinator
である．
まず，初期化処理について説明する．プログラムが起動すると，LoadCoordinator が問題のデータを読み込む．読み込んだ問題データは，LoadCoordinator 内で，一度，前処理（presolve）される．現在の MIP ソルバの前処理は極めて強力で，固定可能な変数の値の固定，冗長な制約式の除去，変数の上・下界値の強化等を行い，問題のサイズが劇的に小さくなることが多い．本論文では，前処理前の問題を元インスタンス（Original instance），前処理後の問題を
前処理後インスタンス（Presolved instance） とよぶ．図 3 に示すように，LoadCoordinator
は，前処理後インスタンスを全 Solver へ転送する．各 Solver は，前処理後インスタンスをSolver 内部では元インスタンスとして扱う．つまり，部分問題は，Solver 内部で，再度，前処理される．この前処理後インスタンスの転送が，プログラム実行中で行われる最も大量のデータ転送になるが，それは Solver が初期化される，この 1 回のみである．各 Solver は，
前処理後インスタンスを Solver 内部に保持し，部分問題の表現として基本的には前処理後イ
ンスタンスと部分問題における変数の上・下界値変化の差分だけを保持する．つまり，部分問
題は，計算過程では常に前処理後インスタンスを対象として解かれ，最終的に求まった解は，
LoadCoordinator 内に保存されている前処理時の情報を利用して，元インスタンスに対する解へ変換され出力される．
初期化が終わると， LoadCoordinator は分枝限定木のルートとなる部分問題を作成する．
このルートとなる部分問題は前処理後インスタンスと同じなので，UG が利用する管理情報を除く部分問題としての情報（前処理後インスタンスとの差分）は空である．この部分問題の情報を含む転送されるデータを ParaNode とよぶ．まず，ルートとなる ParaNode が Solver に転送され計算を開始する．図 1 において，ルートノードが転送されたノードとなっている点に注意されたい．LoadCoordinator 内部では，転送された部分問題が解き終わるまで ParaNode
は Solver の状態を管理するデータ構造内に保管される．各 Solver は，ParaNode を受け取

54

統計数理

第 61 巻

第 1 号 2013

ると，Solver 内部に保持している 前処理後インスタンスと ParaNode の情報により，完全な部分問題を生成して，それを解き始める．
UG は，2 種類の Ramp-Up 手法を提供する．
Normal Ramp-Up. ParaNode を受け取って，遊休状態から部分問題を解いている状態（ 活性化状態とよぶ）となった Solver は，分枝毎に，生成される二つの部分問題の一方をParaNode へ変換して LoadCoordinator へ転送し，もう一方は，Solver 内部で継続して解く．LoadCoordinator は，その内部にノード・プール を持ち，ParaNode を受け取ると，未処理の ParaNode をその下界値の昇順に保持する．ただし，遊休状態の Solver が存在する限
り，ただちに，ParaNode は遊休状態の Solver に割り当てられる．また，遊休状態の Solverが存在しない場合でも，ノード・プールに保存されている「良好な」ParaNode の数が，実行時パラメタで指定された値 p を下回る間は， ParaNode の収集を継続する．ここで，「良好な」
ParaNode とは，その部分問題の下界値と，LoadCoordinator が管理する全ての Solver の計算過程で得た下界値の最小値（これを GlobalDualBound とよぶ）との差が，実行時パラメ
タで指定された相対誤差範囲内にあるものを指す．良好な ParaNode の数が p を超えると，全ての Solver に Ramp-Up 終了の通知を転送する．各 Solver は，その通知を受け取ると，
ParaNode の作成と転送を中止する．
Racing Ramp-Up. LoadCoordinator は，ルートとなる ParaNode を全ての Solver に転送する．それと同時に，Solver 数の異なるパラメタセットを生成し，それぞれを各 Solverに転送する．そして，各 Solver は受け取ったパラメタセットにより，独立にルートとなる問
題を解き，異なる分枝限定木を生成する．各 Solver における下界値や生成されている部分問
題の数，あるいは，制限時間などの情報（実行時パラメタで指定される）により，勝者となる
Solver が LoadCoordinator により決定される．勝者が決まれば，LoadCoordinator は当該 Solver へ勝者となったことを伝える．勝者となった Solver は，生成されている未処理の部分問題を ParaNode に変換して，LoadCoordinator からの送信停止指示を受け取るまで転送する．一方，LoadCoordinator は，勝者以外の Solver に対しては，計算の停止を要求し，その要求を受け取った Solver は実行中の計算を破棄し，新たな部分問題が転送されるのを待つ．本論文では，全ての Solver が独立にルート問題を解いている段階を，Racing 段階とよぶ．通常，動作している Solver の全てに部分問題を与えることが可能になるまで，Racing
段階が継続するようにパラメタを設定して動作させるが，実行環境と扱うインスタンスによっては，そのような設定が現実的でない場合がある．例えば，7,000 Solver を動作させているが，勝者が 7,000 の未処理の部分問題を生成するのを待つのに数日を要する場合もある．そのような場合には， Racing 段階は，十分な部分問題が無い状態で終了し，その後の処理は， Normal
Ramp-Up と同じになる．Racing Ramp-Up は，大規模な並列処理を適用した際に，遊休状態にある Solver を利用して，パラメタの学習，あるいは，チューニングを動的に行うことを意図して設計された．

動的負荷分散方式は，基本的には Shinano et al.（2008）で設計された方式を踏襲している．
MIP を解く上での負荷分散は，最適値に対する上界値と下界値に強く依存する．未処理の部分
問題数だけによる単なる負荷分散では，本来解く必要の無い部分問題の計算に時間を費やすことになる．上界値は，暫定解（それまでに見つかっている実行可能解）よりも小さな目的関数値を持つ実行可能解が見つかることで更新される．ある Solver において，暫定解が更新されると，その解は LoadCoordinator に転送され，その目的関数値，つまり，上界値は，できる限り迅速に全ての Solver へ伝達され，各 Solver では受け取った上界値により計算不要となる制約整数計画ソルバ SCIP の並列化

55

部分問題が削除される．
各 Solver は，定期的に，その状態を LoadCoordinator へ転送する．Solver の状態を示す情報は，その Solver 内の下界値，解いた部分問題数，未処理の部分問題数などであり，
このデータ量は極めて少ない．伝達の頻度は，実行時パラメタで指定する．この情報により，
LoadCoordinator は全ての Solver の状態を把握する．Racing 段階においても，この情報は通知され，この情報に基づき勝者の決定が行われる ．下界値は，分枝限定木の各ノードにおいて，LP を繰り返し解いた結果から求められる．各 Solver は，この情報を転送すると同時に，LoadCoordinator から，LoadCoordinator 中の ノード・プールに存在する ParaNodeが持つ下界値の最小値（これを BestDualBound とよぶ）を受け取る．BestDualBound には，Solver で解いている部分問題に関する情報は含まれていない．
LoadCoordinator がノード・プールに ParaNode を持つ限り，遊休状態の Solver が生
じると ParaNode が即座に転送される．よって，同時に複数の Solver が遊休状態となった場合にも，常に Solver へ部分問題を割り当てるためには，LoadCoordinator は，常にある程度の ParaNode を保持する必要がある．そこで，少なくとも p 個の「良好な」ParaNode を
LoadCoordinator 中のノード・プールに保持するため，Shinano et al.（2008）と同様に， 収集「良好な」とは，その部分木における下界値
モード（collecting mode）を導入した．ここでも，（NodeDualBound）と，全ての Solver で生成されている探索木を考慮した下界値（GlobalDualBound）の差が小さいという意味である．
もし，LoadCoordinator が収集モードではなく，次のような条件:(3.1)

NodeDualBound − GlobalDualBound< Thresholdmax{|GlobalDualBound|, 1.0}を満たす「良好な」ParaNode を p 個未満しか ノード・プールに持たない場合，LoadCoordinator は収集モードとなり，条件（3.1）を満たす部分問題を持つ Solver へ 収集モードになるよう要求を出す．この要求を受け取った Solver は収集モードとなる．現在のソルバは 一般的に探索規則を柔軟に変更できるので，その Solver 内で実行中の分枝限定木の探索規則を保存
し，探索規則を下界値優先探索へ切り替え，LoadCoordinator から収集モードを停止するよ
うに要求が来るまで，分枝毎に 1 個の部分問題を ParaNode へ変換して LoadCoordinatorへ転送する．
LoadCoordinator は，「良好な」ParaNode の数が mp · p（mp も実行時パラメタ）以上になると，収集モード停止要求を収集モード中の Solver へ送る．Solver は，収集モード停止要求を受け取ると，転送を中止し元の探索規則に切り替える．
分枝限定法では，良い上界値が求まれば限定されるはずの部分問題が限定されずに計算されることがある．このような不要な計算は，大規模な並列処理を実現した場合に助長されることが多い．理由の一つは上界値の伝達の遅れであり，もう一つは，ある Solver が極めて悪い下界値の部分問題しか持たないという状態が容易に生じるためである． UG では，この悪い下界値の部分問題しか持たない状態は， LoadCoordinator のノード・プールにある部分問題の下界値 BestDualBound と，Solver 内で計算中の部分問題の下界値を比較することで判断できる．BestDualBound の値が，Solver で計算中の部分問題の下界値に比べて極めて小さい場
合には，LoadCoordinator 中の ParaNode の計算へ移行した方が良い．そこで，Solver は，
新たな ParaNode を LoadCoordinator へ要求し，新たな ParaNode を受け取ると，計算中の分枝限定木を破棄し，受け取った ParaNode の計算を開始する．LoadCoordinator は，
複数の Solver から要求を受ける可能性があるので，常に ParaNode 要求のあった Solver へ
ParaNode を送れるわけではない．新たな ParaNode を受け取れなかった Solver は，そのま
ま計算を続行する．一方，LoadCoordinator が新たな ParaNode を送る場合には，Solver

56

統計数理

第 61 巻

第 1 号 2013

で解かれていた ParaNode は，ノード・プールへ戻され，Solver で計算されていた部分問題の計算は遅らされる．
計算終了フェーズは，LoadCoordinator が，全ての Solver が遊休状態であり，ノード・プールに ParaNode が無いという状態を検出すると始まり，全ての Solver に統計情報の転送を要求する．それらが集まると，解と統計情報を出力して計算が終了する．
3.2 チェックポイントと再スタート
分枝限定法は，一般的に計算終了時間の予測が困難である．スーパーコンピュータでは，各ジョブに実行時間の上限が設定されているのが一般的である．そのため，何らかのチェックポ
イントと再スタートの機能を利用することが必然となる．本論文が扱う問題は，解くことが困
難な問題であり，巨大な分枝限定木が生成されることを仮定している．それらを定期的にチェックポイントファイルへ全て書き出すと，時間を要する上に計算機システムに対しても大きな負荷を与えることになる．そこで，UG は専用のチェックポイントと再スタート機構を持つ．
UG の LoadCoordinator は，分枝限定木全体は保持せず，Solver へ転送された部分問
題だけを保持し，その数は，未処理の部分問題数全体と比較すると極めて少ない．そこで，
LoadCoordinator 中の部分問題だけをチェックポイントで保存することが考えられる．しかし，ある Solver で生成された部分問題が他の Solver へ転送されるため，これらの部分問題
間に祖先・子孫関係が存在する場合がある．再スタート時に，それらを読み込むと，祖先とその子孫となる部分問題がチェックポイントファイル中に存在する場合，子孫から生成される部
分木が冗長に解かれることとなる．従って，LoadCoordinator 中では，部分問題間の祖先・子孫の関係を常に追跡し，祖先が LoadCoordinator 中に存在しない部分問題だけをチェッ
クポイント時に保存する．この部分問題の情報は，具体的には ParaNode の情報であり，変数の上界値・下界値の変更だけを保持しているのでデータ量は極めて少ない．その他には，解法が動作している状態における統計情報である．この統計情報は，各 Solver が送受信した部
分問題数，各 Solver で生成された分枝木のノード数などである．インスタンスデータそのものは，チェックポイント時には保存されない．
再スタート時には，インスタンスデータは，再度読み込まれ，その後，チェックポイントファ
イルから部分問題群が読み込まれ計算が再スタートする．このような再スタートは，大量の計算結果を破棄し，再スタート時にはチェックポイントの状態を回復するために多大な時間を要することになる．このような自明な欠点はあるが，一方でチェックポイントファイルに保存された部分問題の下界値は，その部分問題が生成された時点よりも正確な値に更新されている．
再スタート時には，その更新された情報により部分問題は並べ替えられ，その順に Solver へ送られて解かれるため，チェックポイント時点と全く同じ状態は生成せず，より最適解が得られる可能性の高い方向への探索が進むという解法上の利点がある．定量的な評価は極めて困難であるが，再スタートを繰り返すことで，商用ソルバでは発見することが困難な実行可能解がこれまでに何度も見つかっているため（具体的には，2009 年当時オープンであった MIPLIB2003のインスタンス 6 問中，5 問において最良解を更新している． 1 問は既知の最良解が最適解であった．
），真に最適解を求めることが困難な問題ではかなり有効に機能しているのではないかと思われる．
3.3 ParaSCIP と FiberSCIP の違いについて
並列化は UG フレームワークによって実現されるため，ParaSCIP と FiberSCIP の最も大きな違いは，Solver の実行形態だけであり，論理的には全ての機能が共通となる．しかし，実装
上は全てのデータが共有メモリ内にある FiberSCIP と，共有メモリを持たない環境での動作を実現している ParaSCIP では，SCIP の拡張に対する柔軟性に対して大きな差が生じる．

制約整数計画ソルバ SCIP の並列化

57

FiberSCIP では，共有メモリに存在する SCIP の環境に常にアクセスでき，SCIP の API（Application Program Interface）を直接利用して，並列ソルバを実現している．SCIP の機能が拡張した場合には，その拡張に関連する API は必然的に実装されるので，SCIP の拡張と同時にFiberSCIP の機能も拡張する．一方，ParaSCIP の場合は，SCIP の環境から必要なデータを取り出し，UG フレームワーク内で利用するデータへ変換後，転送し，転送先では逆に UG フレー
ムワーク内のデータから SCIP の環境への設定が必要となる． SCIP の拡張に対する柔軟性を最も損ねるのは，前処理後インスタンスデータを転送する部分である．SCIP は，制約式ハンドラ
が扱うデータ形式を決め，そのクラスに属する制約式のデータを管理している．新たな制約式
ハンドラが追加され，サポートされる最適化問題のクラスが拡張された場合，現在の仕組みを維持して ParaSCIP が拡張されたクラスの最適化問題を解けるようにするには，新たなプログラ
ムコードの追加を余儀なくされる．一方，FiberSCIP は，共有メモリ内に存在する SCIP の環境へのポインタのみを Solver 側へ渡し，SCIP が提供する環境をコピーする関数群を利用して，前処理後インスタンスの転送を実現している．従って，SCIP がサポートする最適化問題のクラスを拡張した場合には，自動的に FiberSCIP がサポートする最適化問題のクラスも拡張する．
Solver 間で転送するデータの種類を追加し，Solver 間で共有する情報により並列ソルバの性能向上を図る場合，プログラム開発とデバッグは，通常のデバッガが利用できるので FiberSCIPの方が格段に容易である．そして，デバッグを終えた後は，わずかなコードの追加により大規
模分散メモリ並列環境で ParaSCIP として動作させることが可能である． FiberSCIP は，ソル
バとしての性能だけが重要なのではなく，大規模分散メモリ並列環境上で動作させるための
ParaSCIP の性能改善を容易に行える開発環境を与えている，という点が重要なのである．
4. 数値実験結果

まず，MIPLIB2010 の Benchmark テストセット（http://miplib.zib.de/miplib2010-benchmark.
php）を利用した FiberSCIP による数値実験の結果を示す．その後，ParaSCIP による大規模計算の代表的な結果をいくつか紹介する．
4.1 FiberSCIP の数値実験結果
数値実験に利用した MIPLIB2010 の Benchmark テストセット（87 インスタンス）の特徴を
表 1 に示す．Benchmark テストセットは，MIP ソルバの性能を 100 未満のインスタンス数で概観することを意図して慎重に選ばれている．表 1 において，Vars. は変数の数，Cons. は制約式の数，NZs は制約行列の非ゼロ要素の数，Bin は 0-1 変数の数，Int は整数変数の数，Contは連続変数の数である．表 1 には，元インスタンスと 1 回の前処理後インスタンスの情報が示されている．現在の MIP ソルバの前処理は問題の規模をかなり縮小することがわかる．最後のカラムの Vars.*Cons. は定常的に必要とするメモリ容量の見積もり値として利用するために示した．
数値実験には，ZIB の Alibaba クラスタの同じタイプの 40 の計算ノードを利用した．各ノードは，2.5 GHz の Quad-Core Xeon E5420 CPU を 2 個搭載する PowerEdgeTM 2950 でメモリ
容量は 16GB である．数値実験に利用した SCIP のバージョンは 2.1.1 であり，LP ソルバとして SoPlex 1.6.0.1 を利用した．FiberSCIP の全ての数値実験は，2 時間の計算時間制限を設けて実施した．Racing Ramp-Up における Racing 段階の停止条件は，2 時間の制限時間の 0.5％ にあたる 36 秒とした．ただし，300 個の未処理の部分問題を勝者 Solver が保持しない場合は，制限時間の 25％ にあたる最大 1800 秒まで Racing 段階を継続する．加えて，実行可能解が見つかっていない状態の場合，Racing 段階は制限時間まで継続する．これは，FiberSCIP では実行不可能性の検出に有効な情報を，他の Solver へ転送していないので，転送した部分問題が表 1.

Benchmark set for MIP (MIPLIB2010 Benchmark set).

58
統計数理
第 61 巻
第 1 号 2013

制約整数計画ソルバ SCIP の並列化

59

表 2.

Sequential executions of SCIP and FiberSCIP.

60
統計数理
第 61 巻
第 1 号 2013

制約整数計画ソルバ SCIP の並列化

61

62

統計数理

第 61 巻

第 1 号 2013

巨大な分枝限定木を生成する可能性が高いためである．本数値実験では，制限時間を考慮してこのような設定を行った．SCIP のパラメタ設定は，基本的にはデフォルトである．ただし，利用メモリ量は 1 Solver スレッド当たり 2GB で収まるように設定しており，数値実験を通して，ディスクアクセスを伴うページフォルトは生じなかった．
まず，逐次処理の SCIP との性能の違い，および，マルチスレッドで Solver を動作させた際の挙動を明確にするための予備実験を実施し，並列化の効果を測定するための基準となる逐
次計算結果を準備する．表 2 は，逐次処理の SCIP と，Solver を 1 スレッドだけ動作させた
FiberSCIP による実験結果である．LoadCoordinator スレッドは存在し，Solver の状態を伝える通信は行うが，Solver スレッドは一つなので逐次処理である．LoadCoordinator による前処理を禁止し，元インスタンスをそのまま Solver スレッドへ転送した場合（Presolve once）と，通常の動作である，LoadCoordinator 側と Solver 側で計 2 回の前処理が実行された場
合の結果が示されている．各設定における実験結果として，解かれた部分問題数と計算時間が示されている．計算時間に “>” が示されているのは制限時間を超えたことを示す．また，“ ! ”が付加されている結果は，プログラムが異常終了したことを示す（制限時間を超えた場合と異常終了の場合は，解かれた部分問題数が正しく出力できないことがあり，正しく出力できない
場合にも 1 と表示されている）．最後の Solution Status は，SCIP および二つの FiberSCIP での計算結果の最悪の計算の終了状態と，正常終了した場合には解の状態を示す．この解の状態は，実行不可能問題の場合には実行不可能であることを出力したら “ok”．最適解をもつ問題の場合に，計算が中断したとき，実行可能解をみつけていれば “ok”，計算が終了し最適値が必要な精度で一致していれば “ok” を示している．MIPLIB2010 には，解のチェッカーが付属しており，生成された解の実行可能性等がチェックできる．本論文に示されている FiberSCIP の全ての数値実験結果は，何らかの解チェッカーにより解が必要な精度で求められているかが常にチェッ
クされている．表の下部には，“解けたインスタンス数（solved）/ 制限時間を超えたインスタン
ス数（timelimit）/ 異常終了したインスタンス数（abort）” が示されている．また，各カラムに対して，幾何平均の値が示されている．この幾何平均の計算には，計算対象のうち 3 ケース全てにおいて解けたインスタンスの結果だけを用い，その解けたインスタンス数は幾何平均の行の最左端に示されている．解けたインスタンス数等の情報の次の行は，3 ケース全てを計算対象としており，全てのケースで解けたインスタンス数は 56 である．それに続く 2 行に示されているのは，左側から 2 つのケースを計算対象とした場合と，右側から 2 つのケースを計算対象とした場合である．最下行の幾何平均は，時間制限に達した，または，異常終了した場合を全て 7200 秒として計算した計算時間の幾何平均であり， H. D. Mittlemann によるベンチマークで採用されている評価方法である．本論文における FiberSCIP による数値実験結果の表は，以降同じ表記方法であり差分を適宜補足する．
SCIP と前処理を 1 回だけ行った FiberSCIP の結果は極めて近かった．前処理を 1 回だけ行う
場合には，SCIP の環境へ読み込まれたデータが，単に Solver 側の環境にコピーされて計算が始まる．計算時間は LoadCoordinator との通信のオーバーヘッドを含むため動作環境の状態によって若干異なるが，ほとんどの結果においてほぼ同じ値を示し，通信のオーバーヘッド
が無視できることがわかる．一方，解かれる部分問題数に関しては SCIP と同じになることが自然であるが，実際には若干の違いを生じた．これは，SCIP のパラメタによりメモリ使用量を制限しているため，探索が深さ優先に切り替わるタイミングの違いから生じていると考えられる．本来の FiberSCIP の実装である 2 回の前処理が実行される場合には，結果は SCIP と比較して大きく異なった．表 2 の下部に示されるように，幾何平均で 47 秒程度遅くなっている．各
インスタンスによる相違を明確にするため，SCIP による計算時間を横軸に取り，対応するインスタンスの FiberSCIP の加速率（SCIP による計算時間/FiberSCIP による計算時間）を縦軸に取っ制約整数計画ソルバ SCIP の並列化

63

図 4. Relative computing time ratio of FiberSCIP compared to SCIP.

たグラフを図 4 に示す．図 4 より，1 回の前処理のケースが，SCIP とほとんど同じ振る舞いをするのに対して，2 回の前処理を行った結果の振る舞いはインスタンスによって大きく異なり，
2 回の前処理のオーバーヘッドはあるが，それが必ずしも減速に繋がるわけではなく，加速するケースもかなりあることがわかる． 2 回の前処理を行うということは， 1 回の前処理とはソ
ルバのアルゴリズムが変わったことを意味し，また，2 回の前処理を FiberSCIP のデフォルト
設定としている．よって，2 回の前処理を行った結果を，並列化の効果を調べるための基準となる逐次処理の結果として採用した．
次に，予備実験として，制限時間内 Racing 段階だけを実行し，Solver 間での上界値の通信を行わない設定で， Solver スレッド数を変化させた FiberSCIP の数値実験を行った．つまり，Solver 数分の異なるパラメタ設定で，それぞれ Solver を実行し，最も速く終了した Solverの実行時間を計算終了の時間とする．表 3 は，Solver スレッド数を，2, 4，6, 8 と増やして行った数値実験結果である．表 2 の値に加えて，加速率と，基準としている 1 Solver スレッドの結果と比較した，解けたインスタンス数の差が示されている． Solver 数を一つ増やすと一つ
多くのパラメタ設定を試す Solver が追加され，上界値の通信はないので，論理的には Solver
数を増やすほど速くなる（上界値の通信を行うと一般的には加速するが，受け取った上界値がSolver に設定されることで，より良い実行可能解の発見が遅れることがあり， Solver 数を増やすほど速くなるとは論理的にも言えない）．しかし実際には，Solver 数が増えると，各Solver がメモリへアクセスする速度は，キャッシュミスやメモリアクセス時のバスの競合などにより低下し，1 Solver 当たりの性能が低下する．実際，表 3 が示すように，最も良い性能を示しているのは，4 Solver スレッドでの結果であり，8 コアのコンピュータを使用しディス
クアクセスが生じないように SCIP のパラメタを設定しているにも関わらず，6 Solver スレッ
ド では既に性能低下しており，8 Solver スレッドでの計算時間の幾何平均値は逐次処理よりも遅く，制限時間内に解けるインスタンス数も逐次処理より少ない．
1 Solver 当たりの性能低下を調べるため， 表 3 のインスタンス名に（i）と（i*）の印をつけた
結果を詳細に調べた．
（i）は，1, 2, 4, 6 Solver スレッド数で同じ数の部分問題を解いたインス
タンスであり，勝者 Solver が同じであると考えられる．一方，（i*）は，6 Solver スレッドだ
けが異なる部分問題数であり，かなり性能の近い競合する複数の Solver の存在を示唆する．
競合する複数の Solver が存在するために 8 Solver による結果だけが異なるケースがあるか表 3.

Only racing stage runs without upper bounds communications.

64
統計数理
第 61 巻
第 1 号 2013

制約整数計画ソルバ SCIP の並列化

65

統計数理

66

図 5.

第 61 巻

第 1 号 2013

Relative computing time ratio of FiberSCIP for each number of Solver threads
used.

もしれないが，その理由で異なるのかどうかの識別ができないため対象から外した．
（i）および（i*）のインスタンスに対する結果について，横軸に定常的なメモリ容量の見積もりとしての（変数の数×制約式の数）を取り，縦軸に加速率をプロットしたグラフが図 5 である．現在のコンピュータアーキテクチャの内部構造は，簡単なモデルで説明できるほど単純ではない．しかし，簡易に入手できる情報だけから予測したいという要求はある．そこで，x を定常的なメモリ利用量として f (x) = 1.0 − c log(x) の割合で 1 Solver スレッド当たりの性能は低下するという仮定で，スレッド数ごとにあてはめを行い，その結果も図 5 に示した．メモリを全く使用しない
ような計算においては性能低下は起こらないという仮定と，メモリ不足などで異常終了せず，計算が実行される限り漸近的にある値に近づくという仮定のもと f (x) を設定した．かなりの上下変動は存在するが，全体としての挙動は，ある程度 f (x) と合致している．係数 c の値は，同時に動作させる Solver スレッドの数が増えれば大きくなる傾向にあり，この値は利用する
コンピュータ環境に依存すると考えられる． FiberSCIP では，スレッド数を増やすと，基本的には 1 Solver スレッド当たりの速度低下を生じるため，実装されているコア数を全て動作させる場合には，並列探索の効果が，このような速度低下を上回る必要がある．
次の予備実験として，先の予備実験と同様であるが，Solver 間で見つかった上界値の通信を行い，それにより限定操作を働かせた場合の効果について調べた．上界値の通信を行わない
場合，計算終了時間が極めて近い Solver スレッドが存在するときには，非決定的な振る舞いはあるが，その例を除くと，非決定的な振る舞いはない．よって，先の予備実験は，各設定で1 回の試行結果を示した．一方，Solver 間で上界値の通信を行うと，それが伝達されるタイミ
ングにより，同じ動作環境，同じパラメタ設定でも，結果が大きく異なる非決定的な振る舞いを示す．従って，この予備実験は，5 回の試行による平均値により評価する．数値実験は先の予備実験で最も良い結果を得た 4 Solver スレッドと，実装されているコアを全て使うことを意図した 8 Solver スレッドにより行った．4 Solver スレッドの結果を表 4，8 Solver スレッドの結果を表 5 に示す．ここでは，上界値の通信による効果を調べたいので，比較対象として
上界値の通信を行わない結果を再度示し，続いて 5 回の試行の結果が示されている．実際，同
じパラメタ設定，同じ Solver スレッド数でも結果が大きく異なるインスタンスがある．左端のカラムには，加速率と 5 回の各試行において解けたインスタンス数の平均が示されており，

制約整数計画ソルバ SCIP の並列化

67

4 Solver スレッドでは，平均で 5％程度加速し，制限時間内に解けるインスタンス数も平均 2
問程度増えている．一方，8 Solver スレッドでは，加速率は，4 Solver スレッドの場合より大きく，16％程度加速し，制限時間内に解けるインスタンス数は平均 1 問程度増えている．並列探索では，早期に良い実行可能解を発見することが多いため，早期に発見された上界値による限定操作の効果と考えられる．ただし，逐次処理と比較すると，8 Solver スレッドでの性
能は表 3 では 15％低下しているので，上界値通信による Racing 段階のみの実行性能は逐次処理の結果と同程度となる．
これまでの予備実験の結果を踏まえて，Racing Ramp-Up と Normal Ramp-Up の性能比較と，2 回の前処理を実行することの効果について調べた．全てのインスタンスの結果を掲載すると多大なページ数を必要とするので，解けたインスタンス数と加速率をまとめた表 6 を示す．
表 6 の全ての結果は，2 回の前処理を実行する逐次処理の結果と比較している．n Solver ス
レッドによる加速率は，Speedup(n) =FiberSCIP 逐次処理の計算時間の幾何平均
FiberSCIP n Solver スレッドによる計算時間の幾何平均として求めているが，幾何平均の対象は，各パラメタ設定において全ての試行で解けたインス
タンスの計算時間であり，異なるパラメタ設定における結果間で比較すると，解けたインスタ
ンスは異なり，解けたインスタンス数も異なる．左端 Settings のカラムに，基準となる逐次処理の結果を含む全ての試行において解けたインスタンス数が示されている．この数のインスタ
ンスが，加速率の計算に用いられた．それと同時に，基準となる逐次処理を除く全ての試行で解けたインスタンス数が（ ）内に示されている．基準となる逐次処理において解けたインスタ
ンス数は 58 である．この値を考慮すると，逐次処理では解けたが並列処理を適用すると解けなかったインスタンスが存在する一方，（ ）内の値は多くのケースで 58 を超えており，逐次処
理では解けないが，並列処理では全ての試行において定常的に解けるインスタンスもあったことがわかる．
表 6 の結果をもとに考察を行う．まず，2 回の前処理を実行することの効果であるが，LoadCoordinator において前処理を行わない場合（no presolve）に，計算時間の観点から性能を落としていることがわかる．この場合，加速率の計算には，1 回だけ前処理を行う逐次処理の結果を基準として利用した方が好ましいが，ここでは他と共通して比較するために，2 回の前処理による結果と比較している．前処理の回数にかかわらず，逐次処理を含む数値実験では，計算対象となるインスタンス数は少なくなる傾向がある．そこで，逐次処理を除いて解けたイン
スタンス数を示す（ ）内の値を見ると，5 回の全試行により解けたインスタンス数という観点においても性能が落ちていることがわかる．この結果は，LoadCoordinator において前処理を行わない場合，各 Solver が定常的に持つインスタンスデータの量が増えることによって，Solver 当たりの性能も低下するので，予想された通りである．したがって，FiberSCIP においては，前処理後インスタンスを転送するというデフォルト設定は妥当であると考えられる．
次に，Normal Ramp-Up と Racing Ramp-Up の結果について考察する．表 6 から，4 Solver
スレッドに対しては，加速率や解けたインスタンス数の観点から Normal Ramp-Up の方が優位と考えられる．しかし，加速率計算の対象としているインスタンス数は， Racing Ramp-Upの方が多い．これは，安定して常に解けるインスタンス数が多いことを意味する．この傾向は，逐次処理の結果を除いた場合でも同様である．4 Solver しか利用しない場合には，Racing 段
階による動的なパラメタ調整の優位性が出ず，最初から並列探索を行った方が加速しているとみなせる．それでも，Racing Ramp-Up は，Racing 段階の結果により動的に挙動を変えること
で，より多くのインスタンスを安定的に解いているという意味では優位性がある．8 Solver ス
レッドにおいては，Racing 段階でのパラメタ調整の効果が顕著になり，制限時間内に解けるイ表 4.

Eﬀects of upper bounds communication (4 Solver Threads).

68
統計数理
第 61 巻
第 1 号 2013

制約整数計画ソルバ SCIP の並列化

69

表 5.

Eﬀects of upper bounds communication (8 Solver Threads).

70
統計数理
第 61 巻
第 1 号 2013

制約整数計画ソルバ SCIP の並列化

71

72

統計数理
表 6.

第 61 巻

第 1 号 2013

Summary of speedups.

ンスタンス数に関しては，Racing Ramp-Up は Normal Ramp-Up と比較して相対的に優位性を増している．加速率は，Normal Ramp-Up の方が良いが，Racing Ramp-Up の方がやはり安定的に多くのインスタンスを解くといえる．全ての試行で解かれたインスタンス数の平均という観点では， Normal Ramp-Up が 4 Solver スレッドから 8 Solver スレッドに増えるにつれ少なくなったのに対して，Racing Ramp-Up の方は増えている．どちらをデフォルト設定にするかの判断は難しいが，Racing Ramp-Up の方が安定的に解けている点と，並列化の規模が増えた際の効果が有利に働いていることから，Racing Ramp-Up をデフォルト設定としている．
最後に，本論文で行った数値実験とは実行環境，実験の条件，評価方法などの点において違いはあるが，FiberSCIP の並列 MIP ソルバとしての性能は H.D.Mittelmann により評価され
WEB ページに公開されている（http://plato.asu.edu/bench.html，ここでは，ug[SCIP/*] ソル
バとして示されている）．最新の並列 MIP ソルバとしての性能に関しては WEB ページを参照されたい．また，本論文においては，ページ数の制限もあるため MIP に対する結果のみを制約整数計画ソルバ SCIP の並列化

73

示した．FiberSCIP は，本原稿執筆時点でも，混合整数非線形計画問題（MINLP: Mixed Integer
Nonlinear Programming）を扱うことができる． SCIP の MINLP ソルバとしての性能は，商用ソルバと比較しても遜色はないため， FiberSCIP はその性能をより強化することが期待されている．
4.2 ParaSCIP の数値実験結果
ParaSCIP による最初の成果は， スーパーコンピュータ HLRN II（https://www.hlrn.de/home/view）を最大 2,048 コア利用し，最適解が得られていなかった MIPLIB2003（http://miplib.zib.de/
miplib2003/index.php）の 2 つのインスタンス ds，stp3d に最適解を与えたことである．ds は，
656 制約式，67,732 の 0-1 変数をもつ問題である．stp3d は，159,488 制約式，204,880 の 0-1 変数をもつ問題であるが，手動で 9 回前処理を行い，88,388 制約式，123,637 の 0-1 変数をもつ問題へ変換したものを解いている．いずれも大規模な問題ではあるが，特に stp3d は緩和問題も巨大で，その計算に時間を要する．計算結果は 10 回以上の再スタートにより求められた．詳細に関しては Shinano et al.（2012）を参照されたい．ここでは，Shinano et al.（2012）に書かれている内容は除き，その後に行った数値実験について述べる．論文としては，結果の新規性に欠けるかもしれないが，研究分野の実情を知る上では有意義ではないかと考える．
まず，当時の状況として，いずれのインスタンスに対しても，ZIB 内では長時間の CPLEX の実行により，暫定解は保持しており，いずれの暫定解も最適解，あるいは，最適解にかなり近い解であると考えられていた．stp3d に関しては，暫定解が最適解であることの確認計算に終わったが，一方，ds に関しては，暫定解の目的関数値を 20％近く更新した．そこで，CPLEX 12.0に得られた解を与えて，CPLEX による計算を，ZIB が所有する 8 Quad-Core AMD Opteron
8384 プロセッサ， クロック 2.7 GHz を 4 個（32 コア）搭載し，512 GB メモリを持つ Sun Galaxy
4600 により行った．約一ヶ月程度の計算の後，CPLEX は異常停止し，その間，ParaSCIP により求まった解を更新することは無かった．そこで，再スタートを行わないで，1 ジョブにより
HLRN II 上で ds インスタンスを ParaSCIP を用いて解くことを試みた．
HLRN II のジョブキューは，1 ジョブ当たり最大利用可能コア数は 2,048 で，最長実行時間
は 12 時間が標準である．この標準仕様のジョブキューでは，再スタートを行わない実行により解くことは無理であると考えられたので，1 ジョブの上限利用可能コア数 4,096，最長実行時
間 48 時間のジョブキューが用意され，解ける見込みがある場合には制限時間の手動による延長が可能となり，1 ジョブで解くことに必要なハードウェアが用意された．一方，ParaSCIP には，当時，Racing Ramp-Up の実装はなかった．Racing Ramp-Up は，1 ジョブによる計算を実現するために開発された．Racing Ramp-Up は，単純なアイデアの実装ではあるが，4,000 以上の異なるパラメタセットに対して安定して動作させるという実装上の困難さがあった．また，動的負荷分散に関しても，1 ジョブでの計算を実現するためには，収集モードとなる Solver
数を制限するとともに，下界値が上昇した場合には，より迅速に収集モードとなる Solver を切替えるなどの工夫を要した．さらに，実行時に 1 計算ノードでもメモリ容量を超えると，その計算ノードで動作していた Solver は停止し，1 Solver でも停止すると， MPI の仕様により全プロセスが停止するので，1 ジョブでの実行はメモリ使用量に注意する必要があった．しかし，本研究で使用メモリ容量を抑えるパラメタ設定をすることは，探索規則を早期に深さ優
先へ切り替えることを意味し，これは無駄に多くの部分問題を解くことになる可能性が高くなる．よって，可能な限りメモリ容量制限まで利用することが好ましい．この種のトレードオフがあるため，パラメタ設定は容易ではない．
表 7 に，HLRN II 上で実現した 1 ジョブでインスタンスを解いた際の実行結果を示す．stp3d
に関しては，最適解を与えて計算を開始している．したがって，最適性の保証計算である．各

統計数理

74

表 7.

表 8.

第 61 巻

第 1 号 2013

Single job computations.

Checkpointing and restarted jobs computation.

インスタンスに対して，計算時間，解かれた部分問題数，転送された部分問題数，転送された部分問題数が解かれた部分問題数に占める割合，動作させた全ての Solver における遊休時間の割合の平均値および最大値が示されている．stp3d に関しては 4,096 コアを利用した Normal
Ramp-Up の結果，および，7,168 コアを利用した Racing Ramp-Up の結果も示している．いずれの結果においても，Racing Ramp-Up 後の負荷分散が良好であることがわかる． ds では，解かれた部分問題数に対して 1.4％の部分問題の転送により負荷分散が実現しており，全 Solverの遊休時間の割合は平均 4.3％，最大でも 10％である．7,168 コアによる stp3d においても，解かれた部分問題の 5.3％にあたる部分問題の転送により負荷分散が実現しており，全 Solver の
遊休時間の割合は平均 6.8％，最大でも 8％である．加えて，4,096 コア，7,168 コアの両 Racing
Ramp-Up の結果において解かれた部分問題数は，4,096 コアによる Normal Ramp-Up の結果において解かれた部分問題数より少ない．stp3d では，最適解を与えて計算を開始しているので，Normal Ramp-Up の方が最初から並列に解く分，速く解けるという予測もできる．このような予測も，インスタンスに依存して変わる．stp3d の場合は，緩和問題の規模も大きいために，Normal Ramp-Up にも時間を要する．このことは，Normal Ramp-Up における遊休時間の割合の大きさからもわかる．遊休は，主に Ramp-Up 時に生じている．このような場合には，むしろ 1 Solver である程度の部分問題を生成してから分配する Racing Ramp-Up が有効となる．
最適解を与えて開始した場合ですら，分枝変数の選択順序により解かれる部分問題の数は変わる．Racing Ramp-Up の勝者の生成した有望な分枝限定木を継続して解くメカニズムは，最適解を与えて計算を開始した場合においても有効であった．4,096 コアを 1.75 倍スケールアップ
した 7,168 コアにより stp3d の計算は 1.34 倍速くなっている．この規模のスケールアップとしては十分な加速が得られたと考えられる．
次に，チェックポイント・再スタートによる計算結果を， 1 ジョブによる結果と比較する．
Shinano et al.（2012）では ds インスタンスが 16 回再スタートによる 17 ジョブで解かれ，計算時間が約 86 時間と報告されている．この時間は，チェックポイントまでの時間を足し合わせて求められている．つまり，チェックポイント以降無駄にした時間は除外されている．本論文では，実際に利用した計算資源の総量を 1 コアで計算した場合の計算時間に変換して，1 ジョブでの計算結果と比較する．表 8 に，各ジョブと利用したコア数および実行時間（Time[h]）を示し，1コアに換算した総計算時間を示す．実際に，計算資源の利用時間は，チェックポイント後の無駄にした時間を含めると約 96 時間であり，1 コア利用の場合に換算すると全体で 181,248 時間になる．一方，4,096 コア利用して 1 ジョブで解いた結果を，1 コア利用の場合に換算すると，

制約整数計画ソルバ SCIP の並列化

75

4, 096 × 76 = 311, 296 時間である．チェックポイント・再スタートでは，再スタート時に，前のジョブで行った計算の多くを破棄し，再スタート時には状態の回復に時間を要する．この回復の計算と，チェックポイント後に無駄にした計算を含めてさえ，チェックポイント・再スター
トの方が少ない計算資源量で問題を解いたことになる．これは，1 ジョブでの実行が，早期に深さ優先探索に切り替わるメモリ不足に関して保守的な設定になっているのに対して，チェックポイント・再スタートを利用した場合には，前の結果に応じて実行時パラメタを調整しており，メモリ不足に関して保守的な設定をしていないことに起因する．実際，表 8 の実験では常に時間制限 12 時間のジョブとして実行していたが，各ジョブの計算時間は大きく異なっている．これは，多くの場合，メモリ不足による異常終了で計算が中断しているためである．極めて困難なインスタンスの場合には，最適解が得られやすい方向へ計算資源を集中させることが期待できるチェックポイント・再スタートの方が，1 ジョブでの実行と比較した場合，大量の計算を破棄しても有効に機能する場合があることがわかる．一方，大規模環境において 1 ジョブで問題を解くということに対するチャレンジによって，負荷分散方式を大幅に改善する結果に繋がったと確信している．大量の計算資源を必要とするので，それを確認するような計算実
験の実施は現実的ではない．しかし，負荷分散がある程度有効に機能しているかどうかを確認するために，大規模環境において 1 ジョブで解くことを試みることは極めて有効である．分枝
限定木全体を集中管理しない場合，1000 を超える Solver 数になると，多くの Solver が無駄な部分木の探索に多くの時間を費やすことを回避するためには実装上の細かな工夫が必要となる．そして，そのような工夫が有効に機能しないとオープンインスタンスを 1 ジョブで解くことは困難である．
MIPLIB2010 が公開され，最適解が得られていない多くのインスタンス（オープン・インスタンス）がライブラリに加わった． ParaSCIP は，それらに最適解を与えることが期待されており，
Koch et al.（2011）が掲載される時点で，4 つのインスタンスには最適解を与えていた．その後，商用ソルバのバージョンが上がる毎に，多くのオープン・インスタンスが商用ソルバにより解かれている（http://miplib.zib.de/ の News に最新情報が示されている）．ParaSCIP は，統計数理研究所のスパコン Distributed-memory supercomputer Fujitsu PRIMERGY RX200S5 上でも動作しており，LP ソルバとして CPLEX 12.3 を利用して，MIPLIB2010 の dg012142（6,310 制
約，640 の 0-1 変数，1,440 の連続変数を持つ MIP）を解いた．このインスタンスは，256 コア
を利用して，約 42 時間の 1 ジョブによって解かれた．解かれた部分問題数は 941,693,415 で，最適値は 2,300,867 であり，MIPLIB2010 のページにも掲載されている．
5. おわりに

本論文では，現在も活発に開発が続けられている，制約整数計画ソルバ SCIP の並列化と，その性能を示す数値実験結果を紹介した．SCIP は商用のソルバと比較すると MIP ソルバとしての性能差はあるが，拡張性のある極めて質の高いソースコードが公開されている．したがって研究目的として利用する，あるいは，特別な用途のための専用ソルバを開発するという用途には適している．FiberSCIP は，SCIP が扱える全てのクラスの最適化問題を並列化できる状態になっており，MIP ソルバとしての性能は商用ソルバに劣っても，MINLP に対する並列ソルバとしての性能は，FiberSCIP が安定して動作する状態になればかなり期待できる．一方，ParaSCIP
は，現在 MIP しか扱えない状態であるため，SCIP の扱える全てのクラスの最適化問題に対して動作することが期待されており，それを実現することは今後の課題である．また，並列ソルバ開発の効率化には，性能が著しく悪くても，決定性並列で動作する並列ソルバ，つまり，同じ動作環境，同じパラメタ設定に対しては，同じ分枝限定木を生成するソルバの開発が必要で

統計数理

76

第 61 巻

第 1 号 2013

あり，これも今後の課題である．
最適化ソルバの性能を示す結果を見る際には，どのような評価により比較されているのかという点に十分注意する必要がある．本論文の実験結果で示したように，分枝限定法をベースとする並列ソルバの性能評価は，どのような評価用インスタンスセットと評価値を用いるかによって大きく変わる．制限時間内に解けるインスタンスを解く速度を重視するのか，制限時間内に解けるインスタンス数を重視するのか，あるいは，H. D. Mittelmann による評価のように総合的な評価指標を重視するのか，評価方法も様々であり確立されているとは言えない．FiberSCIP，
ParaSCIP は，ベータ版ではあるが SCIP Optimization Suite に同胞され，ダウンロードして利用することが可能である．まだ，並列ソルバの研究としては道具が揃ってきた段階であるため，アルゴリズム改善等には多くの課題がある．例えば，並列動作する Solver 数の規模に応じて最も効率的に動作させるには，どのような機構が有効であるかということや，インスタンスとアーキテクチャの適性に関しても調べる必要がある．スーパーコンピュータで解くことを考えると，省電力で動作することも重要な課題である．しかし，現状では，再スタートによる実験結果が示すように，途中で多くの計算結果を破棄することが良い場合もあり，アルゴリズム改
善を優先する必要があると思われる．加えて，ソルバ性能に対する評価方法も決定的なものはない．ベータ版ではあるものの， FiberSCIP と ParaSCIP が公開されたことを契機に，研究を加速することが期待されている．
謝辞

品野勇治と Thorsten Koch は，本研究に対して Google Research Grant のサポートを受けている．JST CREST : ポストペタスケールシステムにおける超大規模グラフ最適化基盤研究代
表の藤澤克樹教授には，開発・デバッグ時に研究室の計算資源を一部使わせて頂いた．ご厚意に感謝致します． HLRN II の大量の計算資源を長期間利用させて下さり，かつ，HPC の技術サポートをしてくださった ZIB の HPC 部門のスタッフの皆様に感謝致します．統計数理研究
所のスーパーコンピュータを利用させて下さった，統計数理研究所の伊藤聡教授に感謝致します．また，技術的なサポートをして下さった統計科学技術センターの蛭田智則様に感謝致します．本論文の原稿を読んでコメントを下さった兵庫県立大学の藤江哲也教授に感謝致します．

参 考 文 献Achterberg, T.（2007）. Constraint Integer Programming, Ph.D. Thesis, Technische Universität Berlin.
Achterberg, T.（2009）. SCIP: Solving constraint integer programs, Mathematical Programming Computation, 1（1），1–41.
Achterberg, T.（2011）. Comparing SCIP to CPLEX, INFORMS 2011 Annual Meeting.
Achterberg, T., Koch, T.and Martin, A.（2006）. MIPLIB 2003, Operations Research Letters, 34（4），361–372.
Bixby, R. E., Boyd, E. A.and Indovina, R. R.（1992）. MIPLIB: A test set of mixed integer programming problems, SIAM News, 25, p.16.
Bixby, R.E., Ceria, S., McZeal, C.and Savelsbergh, M.（1998）. An updated mixed integer programming library: MIPLIB 3.0, Optima, 58, 12–15.
Bixby, R. E., Cook, W., Cox, A. and Lee, E. K.（1999）. Computational experience with parallel mixed
integer programming in a distributed environment, Annals of Operations Research, 90, 19–43.
Bussieck, M. R., Ferris, M. C. and Meeraus, A.（2009）. Grid enabled optimization with GAMS,（3），349–362.
INFORMS Journal on Computing, 21

制約整数計画ソルバ SCIP の並列化

77Chen, Q. and Ferris, M. C.（1999）. FATCOP: A fault tolerant condor-PVM mixed integer program
solver, Mathematical Programming Technical Report 99–05, Madison.
Eckstein, J.（1997）. Distributed versus centralized storage and control for parallel branch and bound:
Mixed integer programming on the CM-5, Computational Optimization and Applications,（2），199–220.
7
Koch, T., Achterberg, T., Andersen, E., Bastert, O., Berthold, T., Bixby, R. E., Danna, E., Gamrath,G., Gleixner, A. M., Heinz, S., Lodi, A., Mittelmann, H., Ralphs, T. K., Salvagnin, D., Steﬀy,D. E. and Wolter, K.（2011）. MIPLIB 2010, Mathematical Programming Computation, 3,
103–163.
Linderoth, J. T.（1998）. Topics in parallel integer optimization, Ph.D. Thesis, Georgia Institute ofTechnology.
Phillips, C., Eckstein, J. and Hart, W.（2006）. Massively parallel mixed-integer programming: Algorithms and applications, Parallel Processing for Scientific Computing（eds. M. Heroux, P.
Raghavan and H. Simon），323–340.
Ralphs, T. K.（2006）. Parallel Combinatorial Optimization, Chapter 3, 53–101, Wiley, Hoboken, New
Jersey.
Ralphs, T. K., Ládanyi, L. and Saltzman, M. J.（2003）. Parallel branch, cut and price for large-scale
discrete optimization, Mathematical Programming, 98, 253–280.
Shinano, Y. and Fujie, T.（2007）. ParaLEX: A parallel extension for the CPLEX mixed integer optimizer, Recent Advances in Parallel Virtual Machine and Message Passing Interface, Lecture
Notes in Computer Science, Vol. 4757, 97–106, Springer, Berlin, Heidelberg.
Shinano, Y., Fujie, T. and Kounoike, Y.（2003）. Eﬀectiveness of parallelizing the ILOG-CPLEX mixed
integer optimizer in the PUBB2 framework, Euro-Par 2003 Parallel Processing, Lecture Notes
in Computer Science, Vol. 2790, 451–460, Springer, Berlin, Heidelberg.
Shinano, Y., Achterberg, T. and Fujie, T.（2008）. A dynamic load balancing mechanism for new
ParaLEX, 14th IEEE International Conference on Parallel and Distributed Systems (ICPADS’08), 455–462.
Shinano, Y., Achterberg, T., Berthold, T., Heinz, S. and Koch, T.（2012）. ParaSCIP—A parallel
extension of SCIP, Competence in High Performance Computing 2010 （eds. Bischof, C.,
Hegering, H. G., Nagel, W. E., Wittum, G.），135–148, Springer, Berlin, Heidelberg.
Xu, Y., Ralphs, T. K., Ladányi, L. and Saltzman, M. J.（2009）. Computational experience with asoftware framework for parallel integer programming, INFORMS Journal on Computing, 21,
383–397.

78

Proceedings of the Institute of Statistical Mathematics Vol. 61, No. 1, 47–78(2013)

Parallelizing the Constraint Integer Programming Solver SCIP
Yuji Shinano1,4 , Tobias Achterberg2 , Timo Berthold1 , Stefan Heinz1 ,Thorsten Koch1 , Stefan Vigerske3 and Michael Winkler1
1 Zuse Institute Berlin
2 ILOG, IBM Deutschland GmbH
3 GAMS Software GmbH
4 JST CREST

The paradigm of constraint integer programming (CIP) combines modeling and solving techniques from the ﬁelds of constraint programming (CP), mixed-integer programming (MIP) and satisfability problem (SAT). This paradigm allows us to address a wide
range of optimization problems. SCIP is an implementation of the idea of CIP and is
now being continuously extended by a group of researchers centered at Zuse Institute
Berlin (ZIB). This paper introduces two parallel extensions of SCIP. One is ParaSCIP,
which is intended to run on a large scale distributed memory computing environment,
and the other is FiberSCIP, intended to run on a shared memory computing environment.
ParaSCIP has been run successfully on the HLRN II supercomputer utilizing up to 7,168
cores to solve a single diﬃcult MIP. It has also been tested on an ISM supercomputer(Fujitsu PRIMERGY RX200S5 using up to 512 cores). The previously unsolved instancedg012142 from MIPLIB2010 was solved by using the ISM supercomputer.

Key words: Mixed integer programming, constraint integer programming, SCIP, parallelization.