日本ソフトウェア科学会第 32 回大会 (2015 年度) 講演論文集

数値制約ソルバーのスケーラブルな並列化
石井 大輔

美添 一樹

鈴村 豊太郎本研究では，数値制約充足問題の区間計算に基づくソルバーをスケーラブルに並列化するための手法を提案する．提案する並列化手法は各 CPU コアにワーカーを配置し，大域的負荷分散スキーム (Global Load Balancing, GLB)に基づきワーカ同士が協調しながら並列に探索処理を行う．GLB 法では，ランダムなワーカ選択と，ライフライン(超立方格子状のワーカ間ネットワーク) に基づいた 2 段階のワークスティーリングを実施し，branch and prune アルゴリズムが扱う探索空間を分散する．提案手法を並列プログラミング言語 X10 を用いて実装した．TSUBAME2.5の 900 コア上で，ベンチマーク問題を用いた実験結果では，最大 751 倍の速度向上を達成した．

以外の) 並列 CSP ソルバーでは数 100 コア程度での1 はじめに並列化にとどまっている (2 節)．最近，Saraswat ら

数値制約充足問題 (NCSP, 3 節) と NCSP の自動[18] が大域的負荷分散 (GLB) の枠組みを提案した．

求解ツールであるソルバーは，実数領域の諸々の問題

GLB は，不均一な並列タスクのための大域的なワーのモデリングと高信頼な解析のための有用な技術であクロードの分散機能とタスクの終了判定機能を提供

る [20] [6] [4]．NCSP では，算術的な制約と，その充足するものである (4. 1 節)．

解を探索すべき領域を表す矩形 (区間ベクトル) が与本研究では，GLB を用いた並列 NCSP ソルバーをえられ，branch and prune アルゴリズムが，矩形を提案する (4. 2 節)．ソルバーは並列プログラミング言区分けしたり，矛盾する箇所を削減したりしながら，語 X10 とその標準ライブラリが含む GLB ライブラ

解集合を指定精度で包む矩形集合を計算する．しかリにより簡潔に実装され，GLB の単位タスクとしてし，NCSP の求解処理は計算量が大きく，既存の逐IBEX ライブラリが提供する数値制約伝播器の逐次実次ソルバーでは扱えるインスタンスの数が限られて装を利用している．5 節では，提案する並列ソルバーいた．そこで，多数の CPU コア上でスケールするよを TSUBAME2.5 の 900 コア上で実行した実験結果

うに NCSP ソルバーを並列化することができれば有を報告する．

用である [7]．
探索空間を区分け・分散する手法は，一般に CSP

2 関連研究ソルバーを並列化する際の単純かつ効率的なアプロー

Gent らによるサーベイ [8] では，並列 CSP ソルチだと考えられている．しかし最新の (実数ドメイン

バーを，(1) 探索空間の区分け・分散によるもの，(2)異質なワーカ群の協調によるもの (cf. ポートフォリ

Scalable Parallelization of a Numerical Constraint
Solver.
Daisuke Ishii, 東 京 工 業 大 学, Tokyo Institute ofTechnology.
Kazuki Yoshizoe, 東京大学, University of Tokyo.
Toyotaro Suzumura, IBM ワトソン研究所, IBM T.J.
Watson Research Center.

オ)，(3) 制約伝播処理の並列化によるもの，の 3 つに分類している．本研究は (1) のアプローチをとる．
探索空間の分散において難しいのは，部分木を均一になるように分散し，各ワーカをアクティブに保つことである．探索木が不均衡になるとき，探索木の適切な区分け方法を予測するのは難しく，求解中に動的に数値制約充足問題 (numerical constraint satisfac-負荷分散することが必要となる．そのために，ワーク

tion problem, NCSP) を 3 つ組として定義する．各スティーリング手法を用い，あるワーカのワークロー

要素は，変数 v = (v1 , . . . , vn )，矩形として与えらドが枯渇しているとき，当該ワーカが他ワーカへリクれる初期ドメイン v 0 ∈ In ，連立式の形をした制エストを送信し，その対応としてワークロードを得る

約 c(v) := f (v) = 0 ∧ g(v) < 0 (f : Rn → R#f ,ことができる．多くの既存並列ソルバーでは，多少な

g : Rn → R#g , #f, #g ∈ N) となっている．NCSPりとも中央集権的な構成をとっており，スケール度合の解は制約 c(ṽ) を充足する変数への付値 ṽ ∈ v 0 であいが限定されていた (例: 40 コア [19], 64 コア [13] [3],る．集合 {ṽ ∈ v 0 | c(ṽ)} を解集合という.

256 コア [2])．
著者らは，GLB を用いて同一の処理を行うワー

3. 1 Branch and prune アルゴリズム

カを協調させた並列 NCSP ソルバーを開発してきた

Branch and prune アルゴリズム [20] は NCSP

[7] [12]．本研究では，NCSP 実装を Realpaver からを求解するための標準的な方法である．アルゴリズム

IBEX に変更するとともに (3. 1 節参照)，よりスケーは，NCSP と精度 ϵ ∈ R>0 を入力とし，解集合を精ラブルな並列化を目指す．

度 ϵ で近似する矩形集合 S ∈ In を出力する．

最適化問題のための branch and bound アルゴリ

アルゴリズムを図 1 に示す．矩形を格納するキュー

ズムについて，ワークスティーリングに基づき探索木L を用意し，2–11 行目のループにおいてキューが空を分散させ，並列化を行う手法が多数提案されているになるまで以下の計算を行う．3 行目では，キューか[9] [14] [16]．アルゴリズムが NCSP の求解処理と似てら矩形を取り出し，Prune 手続きを適用して矩形のいるものの，既存の研究は離散ドメインを対象として端領域の削減を行う．本研究では Prune の実装としいた．本研究では，実数ドメインを区間近似する求解て，HC4Revise [1] 手続きや区間ニュートン法に基づ処理の並列化に取り組む．

く基本的な実装を用いた．矩形に Prune を適用すると，矩形は (i) 解を含まない，(ii) 精度が良い (1 辺3 数値制約プログラミングの最大幅が ϵ 以下)，(iii) 解集合の内部，(iv) その他，数値制約プログラミングは離散ドメインの制約のいずれかに評価される．(i) の場合，矩形は捨てら

プログラミング [17] の拡張として提案され，区間解れ (4 行目)，(ii) と (iii) の場合，矩形は S に追加さ

析 [15] の技法を援用しながら研究が行われてきた．

れ (6 行目)，(iv) の場合，Branch に渡される (8 行

数値制約は実数 R の部分集合を変数ドメインとす目)．Branch 手続きは矩形をある辺の中点で分割し，る．とくに，区間 [a, b] := {r ∈ R | a ≤ r ≤ b} (aその結果をキュー L に追加する (8 行目)．本研究でと b は浮動小数点数) と矩形 (box, 区間ベクトル)は，辺の選択にラウンドロビン法を用いた．

([a1 , b1 ], . . . , [an , bn ]) を扱う．本稿ではボールド体

NCSP ソルバーの (逐次) 実装として，IBEX†1 や(例: v) で区間と矩形を表す．I で区間の集合を表

Realpaver [10] が開発されている．本研究では，Prune

し，I で n 次元の矩形の集合を表す．数値制約の求の実装として IBEX の一部を用いた．

n解は精度保証付き数値計算と branch and prune ア

数値ドメインの CSP に特徴的な，並列求解に関すルゴリズムに基づき行われる．数値制約ソルバーは，る性質をいくつか挙げることができる．まず，Prune

制約を記述する実数関数の区間拡張を計算し，高信の計算量が大きいため (本研究の実験では平均で約 1

頼に制約を評価する．例として，加算の区間拡張をミリ秒かかる)，求解プロセスのボトルネックとなっ[a1 , b1 ] + [a2 , b2 ] := [⌊a1 + a2 ⌋, ⌈b1 + b2 ⌉] と与えるこてしまう点が挙げられる．つぎに，Prune の適用によとができる (⌊·⌋ と ⌈·⌉ は下向き・上向き丸めを表す)．

り，探索木が不均衡となる点が挙げられる．Prune はBranch and prune アルゴリズムは，制約を (区間評価により) 充足する付値を探索する分割統治法である．

†1 http://www.ibex-lib.org/

3. 2 例

Input: NCSP (v, v0 , c), precision ϵ(v1 , v2 ) 平 面 の 2 つ の 円 盤 の 共 通 部 分 を NCSP

Output: set of boxes S(v, v 0 , c) として表すことができる．ただし，v :=1: L := {v0 }; S := ∅(v1 , . . . , v4 ), v 0

2: while L ̸= ∅ do
3:v := Prunec (PollFirst(L))

4:if v ̸= ∅ then5:

S := S ∪ {v}
L := AddLast(L, Branch(v))

9:end if10:

end ifムで求解した結果を図 2 に示す．

本研究では，探索木を分割し，並列動作するワーカ群へと分散することにより，branch and prune アルゴリズムを並列化する手法を提案する．しかし，効率的な並列化のために，探索木を均等に分散する方法

11: end whileは自明でない．本研究では，並列プログラミング言語

12: return S
図1(v1 − 1)2 + v22 − v4 ) = 0 と4 並列化手法else

8:([−1, 1], [−1, 1], [0, 1], [0, 1]),する．ϵ = 0.01 として branch and prune アルゴリズ

if wid v ≤ ϵ ∨ v is an inner box then6:7:c :=:=(v12 + v22 − v3 ,X10 [5] とその大域的負荷分散ライブラリを用いた手

Branch and prune アルゴリズム

法を提案する．

4. 1 X10 GLB ライブラリ

1 v2

GLB は X10 の標準ライブラリに含まれる大域的負荷分散 (global load balancing) 機能の実装である．
ライフライングラフに基づくワークスティーリングア

0ルゴリズム [18] を実装している．GLB は，AI 分野における探索のような，部分タスクの負荷の予測が難し-10v1

1いタスクの並列化に適している．

GLB の計算は協調的なワーカ群によってなされる．
図2

解集合を近似する矩形集合

各 X10 プレースにワーカが 1 つ配置される．各ワー
カはアクティブかアイドルの 2 状態をとりつつ，同探索の対象である矩形を大きく削減したり，さらには一の手続きにより保持するワークロードを処理する．

矩形全体を除去することもある．そのため，探索木のアクティブ状態: ワークロードを保持しているワー

各節点において Prune を適用して探索空間を削減す

カは，タスクを一定量終えるごとに，他ワーカからのるのが望ましいが，他方，均等な分散が難しくなってリクエストを受け取るキューを確認し，保持するワーしまう．3 つ目に，NCSP の制約の数が変数の数よりクロードの一部をリクエスト元のワーカへ送信する．

も少なく，かつ小さい ϵ が与えられると，解の矩形のアイドル状態: ワークロードが枯渇したワーカは，数が大きくなり，探索木が末広がりの形となる点が挙2 フェーズからなるワークスティーリング処理を行う．

げられる．そのため，各枝を並列求解することの効果

まず第 1 フェーズとして，当該ワーカは，ランダムにが期待できる．最後に，探索木の異なる枝の求解処理

他ワーカを選び，リクエストを送信し，リクエスト先間では通信が不要である点が挙げられる．そのため，からワークロードが分け与えられることを期待する．

各枝の処理には既存の逐次実装をそのまま利用する

リクエスト先からワークロードが得られなかった場合ことができる．

は，第 2 フェーズとして，ライフライン (lifeline) と呼ばれる超立方体状のワーカ間ネットワーク上で自分と隣接するワーカに対してリクエストを送信する．

5. 1 TSUBAME2.5 スーパーコンピュータ

第 2 フェーズのライフライン越しの通信は全体のタTSUBAME2.5 は 東 京 工 業 大 学 の ス ー パ ー コ ンピュータである．†2 TSUBAME2.5 の各ノードは 2スク終了処理も兼ねている．

GLB は X10 で実装されており，複数ベンチマークつの Intel Westmere EP 2.93GHz プロセッサ (計 12タスクについて 16,000 プレース程度までスケールすコア) と 54GB のローカルメモリを備える．実験は，ることが報告されている [21]．

75 ノード (900 コア) 上で，高々900 X10 プレースをGLB にはおもに 4 つのパラメタがある．n ∈ N は，用いて行った．X10 は，バージョン 2.4.3.2 のネイティ

主タスクの個数を表し，ワーカがワークロードの分散

ブコンパイラと MPI バックエンド (Open MPI 1.6.5)処理を行う間隔を指定する．w ∈ N はワークスティーを使用した．

リングの第 1 フェーズで送信するランダムリクエス
トの個数を指定する．l ∈ N と z ∈ N は，ライフライ

5. 2 実験結果

ングラフの 1 辺のノード数と，1 ノードでの分岐数を経済学モデルに関する NCSP [11] の 2 インスタン指定する．w := 0 とすると，第 1 フェーズが無効化

スと，球と平面の交差領域を表した NCSP [6] の 2 イされる．l と z については，まず l を指定し，z を lzンスタンスを求解した．GLB のパラメタは，n := 1，がワーカ数以上の最小値となるようにすればよい．lw :∈ {0, 1}，l := 2，z := ⌈logl P ⌉ と設定した．各を最小値 2 にすると直径が最小の密なライフラインインスタンスの仕様と実験結果を表 1 に示す．列グラフとなる．

“problem”, “size”, “ϵ”, “# sol”, “# br”, “w” は，問題名，変数の数，精度，解の矩形の数，Prune 適用

4. 2 GLB を用いた NCSP ソルバーの実装の数，パラメタ w の値を示している．以降の列は実GLB の TaskQueue イ ン タ ー フェー ス を NCSP

験結果を示している．tj は j プレースを用いた実行に 適 合 さ せ た ク ラ ス QueueImpl の 実 装 を 行った ．

時間 (最短値に下線)，ar900 は，各ワーカが branchQueueImpl は 探 索 す べ き 矩 形 を 格 納 し た キュー L

and prune 処理をしている時間の求解時間全体に占と解の矩形集合 S を保持し，branch and prune アルめる割合の平均値 (900 プレース使用時)，# sb900 はゴリズムをカプセル化する．TaskQueue の各メソッ

負荷分散のために 900 プレース使用時にプレース間ドを以下のように実装した．

で送信された矩形の総数を表している．並列求解処理• process(n) は，branch and prune アルゴリズムの速度向上を図 3 に示す．のメインループの繰り返し n 回分，求解処理を行う．

5. 3 考察• split() は L を均等に分割し，片方を L に保持し，実験において，並列ソルバーは 900 プレース・コ

もう片方を戻り値とし，GLB ライブラリに分散アまでスケールし，751 倍の速度向上を達成した (並させる．

列化効率 0.83)．

• merge(L̃) は矩形の集合 L̃ を (リクエスト先から) 受け取り，L に格納する．

サイズが小さい 2 インスタンス (eco8 と sp2-2 ) において，最高の速度向上は設定 w = 0 により得られ実装は，X10 部分が約 1,000 行，C++部分が約た．これらの求解処理では，w = 0 のときワーカのア2,400 行となっており，https://github.com/dsksh/

クティブ率 (表 1 の ar900 を参照) が高くなった．負icpx10 にて入手可能である．

荷分散の第 1 フェーズを無効化し，最初からライフラインを用いた通信オーバーヘッドの大きい負荷分散

5 評価実験を実施している (# sb900 を参照)．小さいインスタ

文献由来の問題を用い，TSUBAME2.5 上で並列NCSP ソルバーの評価を行った．

†2 http://tsubame.gsic.titech.ac.jp/en

表1

problemEconomics(eco8 )
Economics(eco10 )2D sphere andplane (sp2-2 )4D sphere andplane (sp2-4 )

実験結果

size
8

ϵ
10−8

# sol
8

# br
42 2791010−816

1 967 440

2+20.001

3 042 524

3 312 827

2+40.001

19 269 626

26 859 765

w
0
1
0
10
10
1

t1
37s

t600
0.36s
0.68s
6.9s
7.1s
1.4s
1.7s
16.7s
17.0s

3 090s610s
8 550s

t900
0.31s
0.57s
5.1s
4.8s
1.2s
1.5s
11.7s
11.4s

ar900
40%
17%
75%
84%
84%
66%
91%
94%

# sb900
36 K
21 K
1 440 K
601 K
3 370 K
1 080 K20 300 K
6 290 K

900

Speedup

àæ
à
æ

600

751649æ

515æà

æàæà

à

300

àæ

æ
à

æ
à

æ

æ
à

à

æ

120æà

à
æ
à
æ

à

300

600
Places9000æàæàæà

æ
à0àæ

æàà
æ

æ
à

300

600
Places

図39000æà

300600
Places

9000300600
Places900

求解処理の速度向上(左から順に eco8, eco10, sp2-2, sp2-4 の結果を示す．実線は w = 0，点線は w = 1 の結果に対応する．)

ンスでは求解時間が限られるため，高速な負荷分散を行い，すぐに終了処理へ移ることが有利になっている

6 まとめと今後の課題と考えられる．

本論文では，branch and prune アルゴリズムによ

サイズが大きい 2 インスタンス (eco10 と sp2-4 )る NCSP 求解を GLB を用いて並列化する手法につについては，ランダムなワークスティーリングを有効いて述べた．実験では 900 プレース・コアまで線形にした設定 (w = 1) が無効にした設定よりも効率がに近い速度向上が得られることを確認した．

良かった．w = 1 の場合，負荷分散のための通信回数今後の課題として，数値最適化問題の求解手続きをを節約できるメリットがあり，これが大きいインスタ並列化したり，より大規模な並列化に取り組むことがンスでは有効に働いていると考えられる．

挙げられる．

ワークロード分散処理の間隔を指定する n に関し

謝 辞 本 研 究 は 科 研 費 25880008, 15K15968,ては，最も頻繁な設定 n = 1 とした．branch and25700038, 26280024, 23240005 と，JST ERATO のprune アルゴリズムの各繰り返し毎に分散処理が行わ補助を得て行った．

れることになるが，Prune 1 回の処理は 1 ミリ秒前後かかるため，十分なワークロード分散のためにこの頻

参 考 文 献度が必要となる．

[ 1 ] Benhamou, F., Granvilliers, L., Goualard, F.,
and Puget, J.-F.: Revising Hull and Box Consistency, ICLP, 1999, pp. 230–244.
[ 2 ] Bergman, D., Cire, A. A., Sabharwal, A., Samulowitz, H., Saraswat, V., and Hoeve, W.-J. V.: Parallel Combinatorial Optimization with Decision Diagrams, CPAIOR, LNCS 8451, 2014, pp. 351–367.
[ 3 ] Bordeaux, L., Hamadi, Y., and Samulowitz,
H.: Experiments with Massively Parallel ConstraintSolving, IJCAI, 2006, pp. 443–448.
[ 4 ] Caro, S., Chablat, D., Goldsztejn, A., Ishii, D.,

他に本実験の特徴として，求解にかかる実行時間が短くてもある程度の速度向上が得られた点が挙げられる．たとえば，eco8 と sp2-2 は逐次処理に 37 秒および 610 秒かかるが，900 プレース使用時に 120 倍および 515 倍の速度向上が得られている．

and Jermann, J.: A branch and prune algorithm
for the computation of generalized aspects of parallel robots, Artificial Intelligence, Vol. 211(2014),
pp. 34–50.
[ 5 ] Charles, P., Grothoﬀ, C., and Saraswat,
V.: X10: an object-oriented approach to nonuniform cluster computing, Proc. of OOPSLA,
2005, pp. 519–538.
[ 6 ] Ishii, D., Goldsztejn, A., and Jermann,
J.: Interval-based projection method for underconstrained numerical systems, Constraints Journal, Vol. 17, No. 4(2012), pp. 432–460.
[ 7 ] Ishii, D., Yoshizoe, K., and Suzumura, T.:
Scalable Parallel Numerical CSP Solver, CP,
LNCS 8656, 2014, pp. 398–406.
[ 8 ] Gent, I. P., Jeﬀerson, C., Miguel, I., Moore, N.
C. A., Nightingale, P., Prosser, P., and Unsworth,C.: A Preliminary Review of Literature on Parallel
Constraint Solving, Workshop on Parallel Methods
for Constraint Solving, 2011.
[ 9 ] Grama, A., Gupta, A., Karypis, G., and Kumar,
V.: Introduction to Parallel Computing, Addison
Wesley, 2003.
[10] Granvilliers, L. and Benhamou, F.: Algorithm
852: RealPaver: An Interval Solver using Constraint Satisfaction Techniques, ACM Transactions
on Mathematical Software, Vol. 32, No. 1(2006),
pp. 138–156.
[11] Hentenryck, P. V., Michel, L., and Benhamou,F.: Newton: Constraint Programming over Nonlinear Constraints, Science of Computer Programming, Vol. 30, No. 1-2(1998), pp. 83–118.
[12] Ishii, D., Yoshizoe, K., and Suzumura, T.: Scal-able Parallel Numerical Constraint Solver Using
Global Load Balancing, ACM SIGPLAN Workshop
on X10, 2015, pp. 33–38.
[13] Jaﬀar, J., Santosa, A., Yap, R., and Zhu, K.:
Scalable distributed depth-first search with greedy
work stealing, ICTAI, IEEE, 2004, pp. 98–103.
[14] Lüling, R., Monien, B., Reinefeld, A., andTschöke, S.: Mapping Tree-Structured Combinatorial Optimization Problems onto Parallel Computers, Solving Combinatorial Optimization Problems
in Parallel, Vol. 7141, No. 7141, 1996, pp. 115–144.
[15] Moore, R. E.: Interval Analysis, Prentice-Hall,1966.
[16] Otten, L. and Dechter, R.: Towards Parallel Search for Optimization in Graphical Models,
ISAIM, 2010.
[17] Rossi, F., Beek, P. V., and Walsh, T.: Handbookof Constraint Programming, Foundations of Artificial Intelligence, Vol. 2, Elsevier, 2006.
[18] Saraswat, V., Kambadur, P., Kodali, S., Grove,D., and Krishnamoorthy, S.: Lifeline-based global
load balancing, PPoPP, 2011, pp. 201–212.
[19] Schulte, C.: Parallel search made simple, TRICS(Techniques foR Implementing Constraint programming Systems), 2000, pp. 41–57.
[20] Van Hentenryck, P., McAllester, D., and Kapur,D.: Solving Polynomial Systems Using a Branch
and Prune Approach, SIAM Journal on NumericalAnalysis, Vol. 34, No. 2(1997), pp. 797–827.
[21] Zhang, W., Tardieu, O., Grove, D., Herta, B.,
Kamada, T., and Saraswat, V.: GLB : Lifelinebased Global Load Balancing Library in X10,
PPAA, 2014, pp. 31–40.