多面体的ホモトピー法のソフトウェアと並列計算

全文

(1)Shonan Institute of Technology. MEMOIRS OF SHONAN INSTITUTE OF TECHNOLOGY. Vol. 42, No. 1, 2008. 多面体的ホモトピー法のソフトウェアと並列計算郡司貴之* The Software of Polyhedral Homotopy Continuation Method and Parallel Computation of its Takayuki GUNJI* A system of multi variables polynomial equations arises in many fields of science and engineering. The polyhedral homotopy continuation method is one of the methods for finding all isolated complex solutions of systems of polynomial equations. First, the method constructs the homotopy systems adding the parameter for the systems of polynomial equations. Next, it finds the solutions of starting systems of them. It traces the homotopy curves which can be defined as one of the solution set of homotopy systems from the solutions of starting system. As a result of tracing, the method gives the solutions of system of polynomial equations. We implemented the software PHoM (a Polyhedral HOmotopy continuation Method for polynomial systems) which can automatically enumerate all isolated solutions by polyhedral homotopy continuation method, including construction of the homotopy systems, tracing the homotopy curve and verification of the obtained solutions. In numerical experiment, we show efficient of the software. A system of polynomial equations of increasing size have an increasing number of the homotopy curves. The software using single CPU has been limited the size of problems. For solving large size problems, we consider parallel computation of the polyhedral homotopy continuation method. We implemented parallelization through all stages of the method as PHoMpara (PARAllel implementation of the Polyhedral HOmotopy continuation Method for polynomial systems). The software PHoMpara uses Ninf-1 and PC cluster as parallel computer environment. Furthermore, we propose improving efficiency of parallelization. Specially, we indicate a weak point of existing parallel calculation of mixed volume and we suggest improving calculation of it. In numerical experiments, we show effect of the proposed parallel techniques. And we present that can enumerate all isolated solutions of the problems that had not been solved.. 計算量が抑えられる別の手法として，ホモトピー法を. 1．はじめに. 用いた手法がある。ホモトピー法を多項式方程式系の全. 多項式方程式系の全ての孤立解を求めることが必要と. ての孤立解の列挙に用いることをはじめに提案したのは. されている。スプラインなどによる曲線の補間，幾何学. Garcia, Zangwill [10] と Drexler [7] である。これらの論文の. の交点問題，ロボットの動力学，化学の分子状態などか. 中でホモトピー法による多項式方程式系の孤立解の全列. ら多くの多項式方程式系が導出される。また，Morgan. 挙が数学的に可能であることが示された。彼らの提案の. の著書 [19], PoSSo [21] やそれを引き継いだ FRISCO [8] でま. 後，効率良く計算できるように改良されたホモトピー法. とめられた問題集 [26]，Verschelde の web サイト [27] など. が各種提案されている。多面体的ホモトピー法 [14, 28] は. で多くの多項式方程式系が紹介されている。孤立解を全. この提案の中の 1 つである。ホモトピー法の理論が研究. て求める手法の 1 つとして Grobner 基底を構築するため. されると同時に，ホモトピー法を用いたソフトウェアの. に Buchberger のアルゴリズム [5] を利用した消去定理を. 開発も進んでいる。とくに郡司らが過去に発表したソフ. もととする手法がある。この手法は計算量が膨大なもの. トウェア PHoM [11] では解の検証などを行うなどの特徴. であり，方程式の本数が少ない問題でしか解くことがで. を持っている。第 2 節でホモトピー法の概要を簡単に紹. きない。 * コンピュータ応用学科助教平成 19 年 11 月 5 日受付. 介し，第 3 節では郡司らが開発した PHoM の概要を紹介する。第 3 節の最後に数値実験によって，PHoM の性能を簡潔に示す。 — 71 —. NII-Electronic Library Service.

(2) Shonan Institute of Technology. 湘南工科大学紀要第 42 巻第 1 号ステップ 1.. さらに方程式系のサイズが大きな問題を解くために複. 条件を満たすホモトピー方程式系 h (x,. 数台の計算機を並列に使うことによって，計算能力を増. s)0 を構築し，s∞ である初期方程式系 h (x , ∞)0. 大させる並列化に注目した。ホモトピー法に関しては過. の解を求める。. 去からさまざまな並列化の手法が提案されている。過去. ステップ 2.. 初期方程式系の解から h (x , s)0 を満た. に発表したソフトウェア PHoMpara [12] はそのようにパー. しつつパラメータ s を徐々に変化させて，s0 に近づけ. ツごとに提案された並列化手法をまとめ，多面体的ホモ. る。. トピー法の全ての段階を通した並列化の実装を行った。. パラメータ s の領域が [∞, 0] とする方法については [15]. とくにその一部分である混合体積の並列計算について. に従ったものである。ステップ 2.において s0 のとき，. は [25] で提案されているが，この論文では [25] の手法の. つまり h (x , 0)0 は与えられた方程式系であり，求めた. 問題点を提起し，改善できる手法を提案する。第 4 節で. い解が得られるというのがホモトピー法である。ステッ. は並列計算に関する議論を行っていく。とくに混合体積. プ 2.は先ほど定義したホモトピー曲線を追跡しているこ. の並列計算についての詳細を述べ，提案した手法の優位. とに対応する。ただし，ホモトピー曲線の追跡を数値的. 性について数値実験で示す。. に行っているとすると，追跡したいホモトピー曲線から. この論文は郡司が発表した博士論文 [13] をまとめ，一. 別のホモトピー曲線へ移り変わる可能性がある。そこで解が正しく追跡しているかどうかの検証が必要となろう。. 部最新の研究動向も取り上げた。. PHoM [11] や PHoMpara [12] はホモトピー方程式系の構築. 2．ホモトピー法の概要. Start System やホモトピー曲線の追跡 CMPSc に加えて解. 孤立解を求めたい多項式方程式 f (x)0 を f (x)(f1(x),. の検証 Verify を加えた構成となっている。条件を満たすものの一つに線形ホモトピー (linear. f2(x), . . . , f n(x))としたとき，多項式 fi (x)を次のように表記. homotopy) が古くから知られている。しかし，RPS 問. する。. ( ) ∑ c (a ) x (i 1, 2,. . . , n). fi x . a. i. 題 [24] では線形ホモトピー法のホモトピー曲線の本数と (1). なる総次数が 262,144 であるのに対して，孤立解の個数. a ∈ i. が 1,024 である。このように，多項式方程式系の孤立解ただしin (i1, 2, . . . , n,), ci (a){0}(i1, 2, . . . , n, a . の個数と総次数に著しく違いが生じる問題が多く存在す. n),. ることが知られている。孤立解の個数とホモトピー曲線. a. x. x1a1 x2a2. . .xnan,. n は. n 次元の非負整数ベクトル全体. の集合である。. の本数の違いはホモトピー曲線が発散する現象としてホ. 孤立解を求めたい方程式系 f (x)0 に対して，簡単に解を求めることができる方程式系 g(x)0 を考え，f と g. モトピー法にあらわれ，計算時間を大幅に増大される原因の一つであった。. を 1 変数のパラメータ s で結んだホモトピー方程式系. それを改善させるためにいくつか提案された手法の 1. h(x, s)0 を構築する。ホモトピー方程式系を構築するさ. つとして， Huber, Sturmfels [14] と Verschelde, Verlinden,. いに，s∞ のときには h(x, ∞)g(x), s0 のときには. Cools [28] が提案したのが多面体的ホモトピー法である。. h(x, 0)f (x)となるようにしなければならない。また全て. 多面体的ホモトピー法を用いることによって全ての孤立. の孤立解を列挙するために，ホモトピー方程式系はいく. 解を列挙できることは総次数以下の多項式方程式系の孤. つかの条件を満たさなければならない [17]。ホモトピー. 立解の上限を与えるベルンシュタインの定理 [2] に基づ. 方程式系がその条件を満たすことにより，{(x, s)h (x,. き，数学的に保証される。したがって，多面体的ホモト. s)0, s ∞, 0)}が曲線となり，分岐なども存在しない. ピー法のホモトピー曲線の本数は線形ホモトピー法のホ. ことが保障される。この曲線をホモトピー曲線と呼ぶ。. モトピー曲線の数以下になることが保証されている。多. ホモトピー曲線を追跡した結果，得られたものが与えら. 面体的ホモトピー法の特徴の 1 つとして，1 つの多項式. れた多項式方程式系の解となる。. 方程式系に対して複数のホモトピー方程式系を構築して，. ホモトピー法はこのホモトピー方程式系とホモトピー. その全てのホモトピー方程式系を用いて全ての孤立解を. 曲線が大きな役割を果たし，次のようなアルゴリズムで. 列挙することが挙げられる。多面体的ホモトピー法につ. 記述できる。. いては次の節で詳細を述べる。. アルゴリズム 2.1（ホモトピー法の概略） — 72 —. NII-Electronic Library Service.

(3) Shonan Institute of Technology. 多面体的ホモトピー法のソフトウェアと並列計算（郡司貴之）れたものである。. 3．多面体的ホモトピー法. 混合体積の計算は組み合わせ問題であり，非常に計算. この節では，多面体的ホモトピー法について，郡司ら. 量が多い問題である。現在，主流となっている計算方法. が開発したソフトウェア PHoM [11] や PHoMpara [12] を例. は [14] で提案されている混合分割 (mixed subdivision) を. にとり簡単な概要を述べる。その前に PHoM 以外の. 作って，これを基に混合体積を計算し，ホモトピー方程. 1CPU で稼働するソフトウェアについて簡単に紹介する。. 式系を構築する方法である。その後，この計算方法につ. 多面体的ホモトピー法のソフトウェアとして代表的な存. いて速く解くための研究がさまざま行われており，代表. 在なものは PHCpack [29] である。他にも T. Y. Li の web. 的なものとして [17, 25, 18] などがある。PHoM や PHoM-. サイト [16] 上に HOM4PS がある。また多面体的ホモト. para では [25] に基づいて計算を行う。また，これらの開. ピー法にとらわれず，広い意味でのホモトピー法のソフ. 発後の 2007 年に [18] が発表され，1CPU で非常に速い計. トウェアとしては HomLab [23] や HOMPACK [31], HOM-. 算が実現された。PHoM でも導入を検討しなければなら. PACK90 [32] など数多く存在する。. ない課題である。また，[18] のアルゴリズムの並列化は. その中で PHoM は解の検証まで含めた多項式方程式系. 実現されておらず，1 つの研究課題である。この混合体. の全ての複素孤立解の速い列挙を行うために，多面体的. 積の計算については，次の節で並列計算に関連して具体. ホモトピー法を実装し，1CPU で稼働するソフトウェア. 的な説明を行う。 s の次数 r jp(a) (a j , j1, 2, . . . , n) が非常に小さくなっ. として公開した。また，PHoMpara はそれを複数の CPU を用いて並列に計算するソフトウェアとして実装した。. たり，非常に大きくなる可能性があり，そのことは今後. PHoM や PHoMpara は以下のような部分で構成されてい. の数値計算に不安定さを与える。そこでバランシン. る。. グ (balancing) [9] と呼ばれる手法を用いることで，大きい値と小さい値の比を抑える。. 1. 混合体積の計算. このホモトピー方程式系 (2) は [17] の 3 つの条件を満た. 2. 次数の削減. すことはすでに述べられている。このホモトピー方程式. 3. 初期方程式の求解 4. ホモトピー曲線の追跡. 系 (2) で構成されるホモトピー曲線 {(x (s), s): s ∞, 0)} は. 5. 解の検証. s 軸方向に U ターンせず，s に対して単調性を持ってい. これらについて述べる前に，PHoM や PHoMpara で用い. る。PHoM や PHoMpara では，この単調性を利用してア. られているホモトピー方程式系について述べたあと，各. ルゴリズムを構築している。また，与えられた方程式系. 論について簡単に紹介をする。. の係数が一般の位置にあるならば，構成されるホモト. 前の節でも述べたように多面体的ホモトピー法では，複数のホモトピー方程式系を構築し，その全てを利用す. ピー曲線の本数と方程式系の孤立解の個数は一致する特徴がある。. る。PHoM や PHoMpara で用いているホモトピー方程式. 修正ホモトピー方程式系の初期方程式系は二項式方程. 系は孤立解を求めたい方程式系が式 (1) で記述されると. 式系となる。この二項式方程式系はユークリッドの互除. き，次の式で表される方程式系である。. 法を用いることで簡単に解ける。詳細については [17] を. (c˜ (a ) exp (ρ (a ) s) (c (a ) c˜ (a )) exp (( ρ (a ) 1) s)) x. ( ) ∑. h jp x , s ≡. 参照。. p j. j. 初期方程式系の解を求めた後，ホモトピー曲線を追跡. a ∈ j j. ( j 1, 2, . . ., n,. p j. j. p 1, 2, . . ., p *. a. 0 (2). ). する。追跡する数値手法として予測子・修正子法 [1] を使う。予測子・修正子法は 2 つのステップから成り立っている。一つは予測ステップと呼ばれるステップである。このステップでは少しホモトピー曲線から離れるけれど. ただし，c˜ j (a) \{0}(a j , j1, 2, . . . , n) は乱数で選ばれ. も，早く s0 を達成するために s を大きく増加させる役. た数とする。このホモトピー方程式系を構築するために. 割を担う。もう一つは修正ステップと呼ばれる。このス. は，f (x)0 の混合体積を計算しなければならない。方程. テップは予測ステップでホモトピー曲線から離れた分，. 式系 (2) の正の整数 p* と非負数 r jp(a) (p1, 2, . . . , p*,. ホモトピー曲線の追跡の正確性をより高めるためにホモ. a j , j1, 2, . . . , n) はその混合体積の計算によって得ら. トピー曲線に戻る役割を担う。この予測ステップと修正. — 73 —. NII-Electronic Library Service.

(4) Shonan Institute of Technology. 湘南工科大学紀要第 42 巻第 1 号ステップを交互に繰り返すことによって，s を早く 0 に. 問題は cyclic-7 問題 [3]，economic-12 問題 [19]，katsura-11. 近づけながらホモトピー曲線をより良く数値的に追う方. 問題 [4]，noon-n 問題 [20]，reimer-6 問題 [26] である。これ. 法が予測子・修正子法である。. らの問題は問題の後ろの数字によって方程式系の変数の. ベクトル x 0 を初期方程式系 h p(x, ∞)0 の解とする。また exp(s )が十分に 0 に近くなる定数 s を ∞ の近似値 0. 0. 個数が変化する問題である。cyclic-7 問題を除く 4 つの問題については PHoM で解を求めることができる方程式系. として用いる。exp(20)108 より PHoM および PHoM-. の変数の個数が PHCpack [29] で報告されている個数より. para では s020.0 を用いている。PHoM や PHoMpara で. 大きい。この数値実験は AMD Opteron 850 (2.4 GHz) の. は (x 0, s0) から，修正ホモトピー方程式系 h p(x, s)0 に. CPU と 8 GB のメモリを積んだ計算機で行った。また計. よって定義されるホモトピー曲線を s0 まで追跡する。. 算の際には同じ乱数の種を与えている。. 予測子・修正子法は数値解法であり近似解を用いて計. それぞれの問題に対する数値実験結果を表 1 で示す。. 算を行っている。したがって，ホモトピー曲線を完全に. 表 1 は CPU 時間と得られた孤立解の個数を示す。PHoM. 追跡しているわけではない。そこで，ホモトピー曲線を. では曲線を何度か再追跡を行い，値が得られることを確. 正しく追跡したかどうかを事後的に確かめる手法として. 認するが，Iter の列は全ての孤立解を得るまでに何回再. 得られた解の検証を行う。. 追跡を繰り返したかを示す。Trace. i (i1, 2, . . . , 5) の列は. 全てのホモトピー曲線からそれぞれ異なる値の解が得. それぞれ i 回目のホモトピー曲線の追跡と解の検証の. られれば，ベルンシュタインの定理によって，他に孤立. CPU 時間を示す。. 解は存在しないことが保証される。よって全てが異なる. economic-12 問題では 1 回目のホモトピー曲線の追跡に. 値の解を得られたならば計算を終了しても良い。しかし. おいて同じ解に行き着くものが生じた。しかし，パラ. 異なるホモトピー曲線から同じ値の解が得られたならば，. メータを調整して再度追跡することによって，全てのホ. その解が重根か，もしくは単に異なるホモトピー曲線の. モトピー曲線から異なる解を列挙することに成功した。. 解に行き着くような，曲線の乗り換えが発生しているの. これは解の検証のプロセスが効果を発揮したことを示し. かどうかを見極める必要がある。したがって，PHoM や. ている。一方，reimer-n 問題は非常に多くの発散するホ. PHoMpara では得られた値の二点間の距離を用いて，解. モトピー曲線を持つことが知られている。発散するホモ. の検証を行い，近い値を出した場合には曲線を再度追跡. トピー曲線の判定でも解の検証のプロセスが効果を発揮. するようにしている。また発散をした場合や，数値的な. する。reimer-6 問題では 5 回の反復を経過することによっ. 不安定の為に計算を終了してしまった場合についても所. て，ホモトピー曲線を孤立解を持つものと，発散するも. 定の回数，曲線を追跡する。再度曲線を追跡する際には. のに自動的に判別される。noon-n 問題については n6, 7,. パラメータを変更し，細かく追跡するように改良を加え. 8 の 3 種類を示した。解の個数が増加するにつれて，計. ている。. 算時間も増加していることが表 1 よりわかる。これは他. この節の終わりとして，いくつかの問題を 1CPU を用. の問題についても同様であり，詳細は [11, 13] を参照して. いて PHoM で解いたときの数値実験を示す。取り上げた. 表1. 欲しい。. CPU 時間（単位は秒）と計算された孤立解の個数 CPU 時間（単位は秒）. 問題. cyclic-7 eco-12 katsura-11 noon-6 noon-7 noon-8 reimer-6. Iter. Total. S.S.. Trace. 1. Trace. 2. Trace. 3. Trace. 4. Trace. 5. #Sol. 1 2 1 1 2 1 5. 177.3 1497.5 4072.5 113.3 620.5 2898.6 5217.8. 0.1 226.1 215.1 0.1 0.1 0.6 0.1. 177.2 1266.0 3857.4 113.2 619.4 2898.0 1385.4. – 5.4 – – 1.0 – 1322.4. – – – – – – 1827.9. – – – – – – 633.7. – – – – – – 48.3. 924 1024 2048 717 2173 6545 576. — 74 —. NII-Electronic Library Service.

(5) Shonan Institute of Technology. 多面体的ホモトピー法のソフトウェアと並列計算（郡司貴之） (3)を満たす (a , b ) と，任意の j に対して. 4．多面体的ホモトピー法の並列計算. g j,d j j , g jdd j が存在して. 前の節の数値実験で示したように，方程式の変数の個. wj (g )g , a b j 0 ( j1, 2, . . . , n). 数が大きくなるとホモトピー曲線の数が増加し，その結. wj (d j )d j, a b j 0 ( j1, 2, . . . , n). j. 果，計算時間が長くなる。そのため，より変数の個数が. . . (4). j. . . 多い問題を解くためには計算能力の増大を計る必要があ. いま問題 (4) の全ての解 (a , b )(a 1, a 2, . . . , a n, b 1, b 2, . . . ,. る。計算能力の増大を図る様々な方法の中で複数のプロ. b n )を見つけることを考える。線形不等式系 (3) が非退化. セッサを同時に利用することで問題を解決する並列化に. であると仮定すると，問題 (4) は任意の解に対してちょうど 2n 本が等式となる。持ち上げ wj (a ) (a j , j1, 2, . . . ,. よって，大きな計算資源を得ることを検討した。郡司らは変数の個数の多い問題を解くことを目標に並. n) を乱数で選ぶことによって，線形不等式系 (3) が非退. 列計算を行うソフトウェア PHoMpara [12] の実装を行っ. 化であることは確率 1 で満たされる。したがって，問. た。同様に複数の計算機を用いて多面体的ホモトピー法. 題 (4) の全ての解はある j に対して 2 本のみ等式が満たさ. の並列計算を行おうとしているグループに Verschelde ら. れる。問題 (4) の全ての解を (a p, b p) ( p1, 2, . . . , p*) とす. のグループ [30] がある。以下では郡司らが開発した. ると，r jp(a)を r jp(a)wj (a)a, a pb jp (a j , j1, 2, . . . ,. PHoMpara [12] の並列計算に関連する部分について述べて. n) と定めたものをホモトピー方程式系 (2) の r jp(a)と定め. いく。. る。また，. まず，並列計算機環境と通信手段などについて述べる。PHoMpara [12] ではマスタワーカ方式によって並列計算を行った。マスタワーカ方式では，複数の計算機を他の計算機の挙動を管理し，計算に必要なデータの管理を. C jp{a j : wj (a)a, a pb jp0} ( j1, 2, . . . , n),  C p(C1p, C2p, . . . , Cnp)(1, 2, . . . , n). . (5). . で定めた集合 C p (p1, 2, . . . , p*) を混合セルと呼ぶ。. 行うマスタと呼ばれる計算機と，マスタからデータを受. それぞれの j ( j1, 2, . . . , n) に対して要素 g j, d j j を. け取り計算を行い，マスタに計算結果を返すワーカと呼. 決めてから，問題 (4) を満たすかを検証することは組み合. ばれる計算機群に分ける。また，計算機群の間の通信に. わせ的な計算量となるため，非常に効率が悪い。そこで. は Ninf-1 [22] を用いた。. 分枝限定法に似た手法を用いて計算量の削減を図る。以. 多面体的ホモトピー法では複数の段階を経て，全ての解を得ていたが，並列化はその全ての段階で行うことが. 下で問題 (4) の全ての解を求めるために用いる解法の概要を示す。まず説明を容易にするために，次のような集合 S(k). できる。この論文ではもっとも主要な混合体積の計算を行う部分に対しての並列計算について述べる。他の部分. (k0, 1, . . . , n) を定義する。 . の詳細は [12, 13] を参照して欲しい。. . 混合体積の計算と既存の並列手法. まず，並列計算について述べる前に，1CPU の際の混. S(k)  . 4.1. 合体積の計算の方法について簡単に述べる。混合体積は. N (N1, N2, . . . , Nn ): Nj j ,. . #Nj 2 ( j1, 2, . . . , k,),. .  Nj f ( jk1, k2 . . . , n) . 混合分割を構築し，その要素である混合セルを得ること. ここで #Nj は Nj の要素数をあらわし，f は空集合をあら. で計算できる。混合分割の数学的な定義は [14] を参照。. わす。S(n)は混合セルとなりうる全ての候補を含む。. PHoM で用いている手法は持ち上げ (lifting) を用いる方法. さらに次を満たす列挙木を構築する。. であり，[25] に基づいて行われている。. (a) ノード(f , f , . . . , f )S(0) を根とする。. 与えられた方程式 (1) の要素 a j ( j1, 2, . . . , n) に対. (b) 葉のノードである N S(n) は子ノードを持たない。 (c) ノード N S(k) の任意の子ノード N S(k1) の最. して wj (a) を R の開区間から選ばれた乱数とし，これを持ち上げ (lifting) と呼ぶ。このように wi を定めた上で，. 初の k 個目までの要素 N1, N2, . . . , Nk は親である N. 次のような線形不等式系を考える。. の要素と同じである。. wj (a)a, a b j 0 (a j , j1, 2, . . . , n). この列挙木を図にしたものが図 1 である。このように定. (3). 義した列挙木を深さ優先探索で調べることによって，問. 線形不等式系 (3) に条件を追加した以下の問題を考える。. 題 (4) の全ての解を列挙できる。. — 75 —. NII-Electronic Library Service.

(6) Shonan Institute of Technology. 湘南工科大学紀要第 42 巻第 1 号表 2 調べるノードの数による区分に対する部分木の数と 1 ノードあたりの稼動時間の平均（単位は秒） katsura-14 #subtrees. #nodes in a subtree 図1. n2 で #1#23 のときの列挙木の例. (1, 1.0 105) (1.0 105, 2.0 105) (2.0 105, 4.0 105) (4.0 105, 8.0 105) (8.0 105, 16.0 105) (16.0 105, 32.0 105) (32.0 105, 64.0 105) (64.0 105, 128.0 105) (128.0 105, ∞) avr. time for a node. ノード N が混合セルとなるための条件は，N S(n) であり，かつ次の不等式系が実行可能なときである。. b ja, a wj (a) (a Nj, j1, 2, . . . , n), b ja, a wj (a) (a , j \Nj, j1, 2, . . . , n).. (6). 不等式系 (6) の実行可能性は k {1, 2, . . . , n} の列挙木の任. 13342 151 179 218 189 194 130 88 29 4.36 10 105. eco-15 #subtree 6214 1024 820 596 427 271 150 37 16 3.64 105. 意のノード NS(k) において，不等式系 (6) を制約として持つ線形計画問題として確かめられる。もしあるノード N S(n) において不等式系 (6) が実行不可能ならば，その. 分木のノードがワーカによって実行可能性を調べるが，. 子 N S(k1) に対応する不等式系は必ず実行不可能で. 部分木の各ノードに対して不等式系 (6) の実行可能性を. ある。したがって，N をルートとする部分木には混合セ. 調べ，実行不可能なときにはそのノードの下の部分木の. ルは絶対に存在しないことがわかる。その結果，不等式. ノードは実行可能性を調べないため，それぞれの部分木. 系 (6) が実行不可能なノード N S(k) の部分木は調べる必. で実行可能性を調べたノードの数は異なる。そして部分. 要がなくなり，高速化が図られる。. 木を探索する前に，実行可能性を調べたノードの数を知. ここまでが 1CPU における混合体積の計算の概要であ. ることはできない。もしそれぞれの部分木で実行可能性. る。混合体積の並列計算を行う際にもっとも容易な方法. を調べたノードの数が非常にばらつくならば，計算時間. は，あるレベル l を定め，集合 S(l ) の中に ml 個のノード. もワーカごとに著しく異なったものになるであろう。. i. i. 表 2 は katsura-14 問題と，economic-15 問題において，. N (i1, 2, . . . ,ml ) があるとする。この N をルートとする部分木を，それぞれ別のワーカに割り当てることで並列計. 部分木の中で実行可能かどうかを調べたノードの数を区. 算を実現する方法である。これは，それぞれの部分木の. 別した部分木の数を示している。例えば，katsura-14 問. 間のノードは交わることはないため，S(l ) のあるノード. 題は調べたノードの数が (1,1.0 105) の範囲にある部分木. i. N (1 i ml ) をそれぞれ異なるワーカに割り当てて計算. が 13,342 個ある。また economic-15 問題は調べたノード. しても，他のノードに関係なく計算できることから並列. の数が (1.0 105, 2.0 105) の範囲にある部分木は 1,024 個あ. 計算が可能となっている。この方法は [25] で述べられて. る。表 2 の最終行はワーカの稼働時間を調べたノードの. いた方法である。. 数で割ったものを示している。1 つのワーカで 1 つの部. 以下の議論では PC クラスタの中でワーカの台数を w*. 分木を計算するのにどれぐらい時間を消費したかを示す. 台とし，レベル l を w* #S(l ) が満たすように選ぶ。また，. 指標として「調べたノードの数」「1 つのノードの実. このような l が常に存在するものと仮定する。. 行可能性を調べることに対する平均計算時間」を用い. 次にワーカの割り当てについて [25] で述べられていた 1. 既存手法について簡単に説明する。まず S(l ) の N から N. w*. をそれぞれワーカ 1 からワーカ w* に割り当てる。そ. れ以外のノード N i (w*1 i ml) について，ワーカが計. る。表 2 を見ると，katsura-14 問題は多くても 1.0 105 個のノードしか調べない 13,342 個の部分木と，128.0 105 個以上のノードを調べる 29 個の部分木がある。この大きな違いはワーカ間の稼動時間の不均衡を引き起こす。そこで，混合セルの並列計算の効率良い実装としてア. 算を停止したとき，残っているノードを停止したワーカ. ルゴリズム 4.1 を提案する。アルゴリズム 4.1 の基本的な. に割り当てる。 4.2. 考え方は以下の 2 点である。. 提案する混合体積の並列計算手法. [25] の既存手法をそれぞれの問題に適用したとき，部. (a) 割り当てられた部分木のノードにおいて不等式. — 76 —. NII-Electronic Library Service.

(7) Shonan Institute of Technology. 多面体的ホモトピー法のソフトウェアと並列計算（郡司貴之）系 (6) の実行可能性を調べることに対して，長い計. （ワーカにおける計算）. 算時間を要するワーカは全ての計算が完了する前. 入力： N と ub ；出力：スタック R と Mw.. に強制的に停止させる。. ステップ 1. カウンタ co←0 を設定する。R を空のスタックとし，N をスタック R に挿入する。. (b) (a) によって停止させられ，実行可能性のチェック. ステップ 2. co←co1 を設定する。スタック R から 1 つのノード N を取り出す。N S(l¯) を満たす ¯l と. が終わっていないノードの計算を空いている他のワーカに再分配させる。. する。R←R {N } を設定する。. 今回の実装方法においてマスタはワーカが計算してい. ステップ 3.. る間，計算時間が長いワーカを知ることができない。そ. もし N N とする線形方程式系 (6) が実. 行可能ならば，ステップ 4 に進む。そうでないな. こで，強制的に停止させるために数字 ub をワーカに送. らば，ステップ 5 に進む。. る。ub を受け取ったワーカは ub 個のノードの実行可能. ステップ 4. もし N が列挙木の葉のノード，すなわ } を設定する。そうち ¯l n ならば，M ←M {N. 性のテストを終了した後で一度仕事を終了する。そのようにして強制的に停止したワーカは実行可能性のチェッ. w. w. マスタはワーカからの仕事の終了の通知を待つだけでよ. でないならば，N の全ての子ノード N 1, N 2, . . . , ¯ N S (l 1) を列挙木から取り出し，それらを. く実装は容易である。マスタとワーカに対するアルゴリ. スタック R に積む。. クが完全に終了していない旨をマスタに通知する。一方，. m ¯l 1. ステップ 5.. ズムを以下に与える。アルゴリズム 4.1（混合セルの並列計算）. もし R0 ならば空のスタック R と Mw を. マスタに送り，終了する。そうでなくもし coub. （マスタにおける計算）. ならば，スタック R と Mw をマスタに送り，終了する。そうでないならば，ステップ 2 に進む。. 入力：初期レベル l ；出力：混合セルの集合 M. ステップ 1. 数字 ub を選び，M←0 を設定する。列挙 ml. アルゴリズム 4.1 の中で数字 ub は再分割のために使われ. 木の中から，S(l ) の全ての要素 N1, N2, . . . , N を取り出す。. る。アルゴリズム 4.1 のワーカにおいて，スタック R を. Q←{N w*1, N w*2, . . . , N ml} と W←{1, 2, . . . , w*} を設定す. 使うことで列挙木に対して深さ優先探索を実行している。図 1 の列挙木を例に用いて，アルゴリズム 4.1 のワー. る。ステップ 2. ∀w W に対してノード N w と ub をワー. カの挙動を説明する。まずワーカは N 0 と ub3 を入力. カ w に送り，N w を根に持つ部分木に含まれる混. として受け取ったと仮定する（図 2 の(i)）。ステップ 2 に. 合セルの計算をワーカ w に要求する。. おいて，co1 となりスタック R から N 0 を取り出す。. ステップ 3. 1 つのワーカが停止するまで待つ。もし. N N 0 における線形方程式系 (6) が実行可能であると仮. ワーカが止まったならば，そのワーカをワーカ w˜. 定する。そのとき，ステップ 4 において，N 1, N 2, N 3 が. とし，スタック R と集合 Mw˜ をワーカ w˜ から受け. スタック R に積まれる（図 2 の (ii)）。次にふたたびステッ. 取る。スタック R の要素を Q に加える。. プ 2 において，co2 となり，スタック R から N 1 を取り. W←W {w˜}と M←Mw˜ を設定する。. 出す。N N 1 における線形方程式系 (6) もまた実行可能. ステップ 4.. もし Q0 ならばステップ 4A に進む。そ. うでないならば，ステップ 4B に進む。ステップ 4A.. がスタック R に積まれる（図 2 の (iii)）。ステップ 2 にお. もし W0 ならば M を出力し終了. する。そうでないならば，ステップ 3 に進む。ステップ 4B.. であると仮定する。ステップ 4 において，N 11, N 12, N 13 いて，co3 となりスタック R から N 11 を取り出す。もし N N 11 における線形方程式系 (6) が実行不可能なら. もし #Ww* ならばステップ 3 に. 進む。そうでないならば，ステップ 5 に進む。ステップ 5.. Q から 1 つのノード N を選び，N S(l). を満たす l とする。{1, 2, . . . , w*}W から w を選ぶ。N と ub をワーカ w に送り根として N を持つ部分木に含まれる混合セルの計算をワーカ w に要求する。Q←Q{N } と W←W{w} を設定す図2. る。ステップ 4 に進む。. ワーカにおけるスタックの挙動例. — 77 —. NII-Electronic Library Service.

(8) Shonan Institute of Technology. 湘南工科大学紀要第 42 巻第 1 号ば，スタック R から N 11 が取り除かれるだけである（図. 式系の構築の稼動時間の増加はホモトピー曲線の追跡に. 2 の (iv)）。ステップ 5 において，coub であると判定さ. 対するそれよりも大きいことがわかる。他方，n が 9 か. れる。このとき，スタック R は {N 12, N 13 ,N 2, N 3} であ. ら 12 まで増加する中で，noon-n 問題におけるホモトピー. る。ワーカはこの R と Mw0 をマスタに送る。アルゴリ. 方程式系の構築に対する稼動時間の消費はホモトピー曲. ズム 4.1 のマスタにおけるステップ 3 の中で，マスタは. 線の追跡に対する稼動時間が増加するのに比べて，それ. スタック R を受け取り，Q に加える。. ほどの増加を示していない。ホモトピー曲線の追跡の最. アルゴリズム 4.1 のワーカにおいて，ワーカはスタッ. 初の反復におけるホモトピー曲線 1 本あたりの平均と最. ク R をマスタに送ると記述している。しかし，列挙木の. 大の計算時間の増加は問題のサイズの増加に比例してい. 任意のレベル l において部分木のノード N S(l ) を順序. ないこともわかる。例えば，noon-10 問題の平均計算時. 付けることができれば，スタック R の先頭要素のみを送. 間は 1.0 秒であるが，noon-11 問題は 2.1 秒，noon-12 問題. るように修正することができる。詳細は [12] を参照して. は 3.3 秒である。. ほしい。 4.3. 表 3 は reimer-n 問題がホモトピー曲線の追跡において，. 数値実験. いくつかの困難があることも示している。reimer-7 問題. この節では，郡司らが実装した PHoMpara [12] を用い. の孤立解の数は 2,880 である。この数は他の問題に比べ. て行った数値実験について述べる。取り上げた問題は. て大きな数字ではない。しかし reimer-n 問題は非常に多. cyclic-n 問題 [3]，economic-n 問題 [19]，katsura-n 問題 [4]，. くの発散するホモトピー曲線を持つ。孤立解は最初の反. noon-n 問題 [20]，reimer-n 問題 [26]，RPS-10 問題 [24] であ. 復か 2 回目の反復で見つけられる。しかし，いくつかの. る。RPS-10 問題以外の 5 つの問題については n が大きく. 発散するホモトピー曲線は 2 回目の反復以降も発散と判. なるにつれて方程式系の変数の個数が増える問題であり，. 定されるまで追跡される。例えば，reimer-7 問題におい. RPS-10 問題は固有の問題である。. て，最大で 6 回繰り返しホモトピー曲線を追跡する。. この数値実験は 2006 年に東京電機大学に居られた藤. cyclic-8 と cyclic-9 は非特異な解と特異な解の両方を持. 澤克樹研究室に設置された SDPA クラスタを使用させて. つことが知られている。一方 cyclic-10 と cyclic-11 は非特. いただいた。SDPA クラスタは 40 台の PC からなり，そ. 異な解しか持たない。cyclic-n 問題の特異な解についての. れぞれの PC は AMD Athlon 2.0 GHz が使用されていた。. 詳細は [27] を参照。これらの非特異な解と特異な解を持. なお，cyclic-n と RPS-10 の実験時には 1 台が使用できな. つ問題に対しても，PHoMpara は解の分類を行えている。. い状況であったので 39 台で実験を行った。. 4.4. 数値実験結果をまとめたのが表 3 である。表 3 の中で. ワーカの台数による比較. economic-14 問題をワーカの台数が 1 台から 37 台まで. “Total”は総稼動時間（単位は秒），“StSy”はホモトピー方. 変えたときの比較を行った数値実験結果を表 4 で示す。. 程式系の構築で消費した稼動時間，“Tr. i”は i 回目のホモ. 他の問題についても同様の数値実験を行っており，[12,. トピー曲線の追跡と解の検証に対する稼動時間，“Tr. all”. 13] を参照してもらいたい。‘ratio’の列はワーカ 1 台に対. は Í 6i1“Tr. i”, “#Sol.”は孤立非特異解の数である。. するそれぞれの稼動時間の比を表す。表 4 は総稼動時間，. 表 3 にある economic-n 問題，noon-n 問題，katsura-n. ホモトピー方程式系の構築に対する稼動時間，ホモト. 問題，RPS-10 問題の数値実験結果から，全てのホモト. ピー曲線の追跡と解の検証に対する稼動時間，そして，. ピー曲線から孤立解が得られており，発散する曲線が存. それぞれの比を含んでいる。時間の単位はいずれも秒で. 在しないことがわかる。ホモトピー曲線を追跡するとき. ある。. に他のホモトピー曲線に乗り換えたなどにより解の検証. economic-14 問題のホモトピー曲線の追跡は 2 台のとき. によって，不正確に追跡されたと判定されたホモトピー. に 2.00 倍，37 台のときに 34.87 倍という効率を示してい. 曲線は Tr. k (k2)の列が示すように繰り返し追跡が行わ. る。37 台のときに 34.87 倍である原因を追求するために，. れる。noon-12 問題において，531,417 個の孤立解を得る. 単一の CPU だけで計算を行っても並列で計算しても変化. ために 49,458 秒を要した。いくつかのホモトピー曲線は. しない時間や，並列計算を行うことによって，単一の. 最大で 4 回繰り返し追跡された。katsura-n 問題におい. CPU だけで計算していた時には生じない時間をまとめた. て，n が 11 から 15 までの間，孤立解の個数は n が 1 増. ものを表 5 で示す。表 5 に掲載したデータを引いて計算. えることに倍増していることがわかる。ホモトピー方程. すると 22649/626.636.14 倍となる。表 5 に掲載した時間. — 78 —. NII-Electronic Library Service.

(9) Shonan Institute of Technology. 多面体的ホモトピー法のソフトウェアと並列計算（郡司貴之）表3. PHoMpara の稼動時間（単位は秒）と得られた孤立解の個数。#cur: ホモトピー曲線の本数 CMPSc. Prob.. Total. time. StSy. time. time. cyclic-8. 573. 3. 570. cyclic-9. 1016. 4. 1012. cyclic-10. 1037. 10. 1027. cyclic-11. 8085. 49. 8056. eco-14. 626. 388. 238. eco-15. 2964. 2300. 664. eco-16. 12051. 10470. 1581. katsura-11. 160. 58. 102. katsura-12. 604. 142. 462. katsura-13. 2137. 531. 1606. katsura-14. 4187. 2231. 1956. katsura-15. 18964. 13638. 5326. noon-9. 314. 20. 294. noon-10. 1797. 27. 1770. noon-11. 10225. 39. 10186. noon-12. 49458. 78. 49380. reimer-5. 18. 1. 17. reimer-6. 110. 3. 107. reimer-7. 4229. 9. 4220. RPS-10. 894. 638. 256. Tr. all. Tr. 1. Tr. 2. Tr. 3. Tr. 4. Tr. 5. Tr. 6. time #cur. ave. max.. time #cur. time #cur. time #cur. time #cur. time #cur. 26 2560 242 11016 813 35940 7883 184756. 0.3 0.8 0.7 2.2 0.8 9.2 1.5 18.2. 25 1408 180 5022 69 1000 86 99. 67 896 78 648 31 60 67 22. 452 896 512 648 114 60 – –. – – – – – – – –. – – – – – – – –. 238 4096 637 8192 1566 16384. 2.1 6.7 2.9 10.4 3.6 19.7. – –. – – 19 6 – –. – – – – – –. – – – – – –. – – – – – –. 102 2048 324 4096 799 8196 1907 16384 5224 32768. 1.9 4.6 3.2 8.2 3.7 10.0 4.5 10.4 6.1 21.5. – – 7 16 9 8 17 34 45 61. – – 17 8 20 7 32 10 57 25. – – 114 2 25 2 – – – –. – – – – 753 2 – – – –. – – – – – – – – – –. 279 19665 1752 59029 9848 177125 46737 531417. 0.5 1.3 1.0 2.6 2.1 5.5 3.3 9.1. 6 36 18 2 94 180 860 333. 9 12 – – 97 66 871 127. – – – – 147 12 912 26. – – – – – – – –. – – – – – – – –. 4 720 31 5040 398 40320. 0.1 0.2 0.2 0.5 0.3 1.5. 4 577 30 4464 399 37448. 4 173 38 1695 524 15512. 5 1 8 4 1363 5299. – – – – 978 608. – – – – 558 8. 256 1024. 9.4 22.1. – –. – –. – –. – –. – –. 8 6 15 8. #Sol.. 1152 5994 34940 184756. 4096 8192 16384. 2048 4096 8192 16384 32768. 19665 59029 177125 531417. 144 576 2880. 1024. — 79 —. NII-Electronic Library Service.

(10) Shonan Institute of Technology. 湘南工科大学紀要第 42 巻第 1 号表 4 economic-14 問題に対するワーカ台数が異なる実験における稼動時間（単位は秒）と比。#wks: ワーカの台数 Total 問題. time economic-14. 1 2 5 10 20 37. 表5 4.0 秒 3.1 秒 38.8 秒. StSy. CMPSc. #wks. 22653 11317 4536 2292 1178 681. ratio. time. 1.00 2.00 4.99 9.88 19.23 33.26. 13620 6804 2727 1383 718 422. ratio. time. 1.00 2.00 4.99 9.84 18.96 32.81. 9033 4513 1809 909 460 259. ratio 1.00 2.00 4.99 9.93 19.63 34.87. economic-14 問題における，ワーカを 37 台使用した際に生じた時間. 並列計算できない部分（線形計画問題を解く部分など）混合体積の計算における通信に要した時間と通信を開始するまでに要した時間の各ワーカにおける平均 (*1) 混合体積の計算における次の問題をマスタから割り当てられたり，他のワーカが計算を終了するまでの待機時間の各ワーカにおける平均 (*2) 曲線の追跡における通信に要した時間と通信を開始するまでに要した時間の各ワーカにおける平均 (*1) 曲線の追跡における次の問題をマスタから割り当てられたり，他のワーカが計算を終了するまでの待機時間の各ワーカにおける平均 (*2). 0.1 秒 8.4 秒. によって台数と同じだけの効果を得ることが難しいこと. 計算は以下の式によって計算した。. がわかる。. ∑. 他の数値実験も含めて多くの数値実験において，ホモ. w* w1 τ w. w * × max ww1 τ w. トピー方程式系の構築に比べて，ホモトピー曲線の追跡. *. (7). において並列化効率がより良いことが示されている。これは 1 本のホモトピー曲線の追跡は他のホモトピー曲線. ただし t w はワーカ w の稼動時間である。したがって，稼. と独立に計算することができるという特性から来るもの. 働率が大きいならば，休眠しているワーカが少ないこと. である。. を表す。たとえば，稼働率が 1 ならば全てのワーカは同. 4.5. 再分割手法の効果の比較. 時に計算を終了し，他のワーカの仕事を待っていたワー. 4.2 節で述べた再分割手法の比較実験を行った。kat-. カは存在しない。一方，1 台だけが計算していて，他の. sura-14 問題，katsura-15 問題，eco-15 問題，eco-16 問題. 全てのワーカが休眠していたならば，稼働率は 1/w* とな. に対してアルゴリズム 4.1 を適用した。今回の数値実. る。. 験では ub2.00 106 と定めた。 katsura-14 問題におい. 表 6 は多くのテスト問題において，再分割手法が 1 に. て調べたノードの総数を数え上げたカウンター co は (4.36. 近い良い稼働率を実現することを意味している。結果と. 10 5 ) (2.00 10 6 )87.2 秒で ub の値 2.00 10 6 になる。. して，katsura-14 問題や eco-15 問題において最大の稼動. economic-15 問題ではその計算結果は 72.8 秒である。. 時間が 10% 減少した。とくに，katsura-15 問題は再分割. 6. ub2.00 10 の値を使うことによって，それぞれのワー. 手法によって，非常に高速に解くことに成功した。再分. カは割り当てられた仕事の実行を 100 秒以下に分割して. 割手法は最大の稼動時間を減少させることによって，大. いる。. きな規模の問題を解くことに効果があることを確認した。. 数値実験結果を表 6 に示す。ワーカごとの稼動時間の平均，最大，最小の時間をそれぞれ平均，最大，最小の. また，最小値と最大値の間も狭くなり，偏りが抑えられていることがわかる。. 列に示した。いずれも単位は秒である。最大の稼動時間. 本論文で採用した再分割手法による動的な負荷分散を. は全体の計算時間に影響することに注意せよ。稼働率の. 行う方法とは別に，最初から細かい問題規模に分割して. — 80 —. NII-Electronic Library Service.

(11) Shonan Institute of Technology. 多面体的ホモトピー法のソフトウェアと並列計算（郡司貴之）表6. 再分割手法を適用しない場合と適用した場合の各ワーカの稼働時間（単位は秒）と稼働率再分割手法を適用しない. 問題. 再分割手法を適用する. 台数平均. 最大. 最小. 稼働率. 平均. 最大. 最小. 稼働率. katsura-14. 10 20 40. 8650 4339 2164. 8804 4596 2462. 8577 4251 2087. 0.982 0.944 0.878. 8684 4338 2183. 8700 4380 2209. 8675 4312 2172. 0.998 0.990 0.988. katsura-15. 40. 13459. 17097. 12269. 0.787. 13461. 13536. 13388. 0.994. economic-15. 10 20 40. 8862 4427 2262. 8919 4506 2446. 8843 4399 2238. 0.993 0.982 0.927. 8888 4432 2253. 8907 4457 2276. 8881 4411 2247. 0.997 0.994 0.990. economic-16. 40. 10324. 10632. 10257. 0.971. 10393. 10428. 10377. 0.996. 表 7 再分割手法と問題規模を変化させたときの各ワーカの稼働時間，稼働率，通信時間（時間の単位はいずれも秒）問題. l. 再分割. katsura-14. 2. しないするしないするしないするしないする. 3. economic-15. 2 3. 部分木数. 平均. 標準偏差. 最大. 最小. 稼働率. 通信時間. 14520 30842 95172 98210. 2393 2416 3078 3121. 103.1 13.2 10.0 10.1. 2693 2463 3122 3173. 2311 2402 3069 3113. 0.888 0.980 0.985 0.983. 5.00 20.54 626.94 666.24. 9555 17478 428318 428519. 2504 2505 18974 18824. 32.9 8.3 8.0 12.9. 2688 2530 19012 18893. 2492 2500 18966 18818. 0.931 0.990 0.998 0.996. 1.91 3.40 16126 16114. から並列計算を行う方法も考えられる。この問題におい. 分割を行うよりも改善される。しかし，通信時間が顕著. て，どちらが望ましいかを検証するために詳細な数値実. に増大し，その結果ワーカの稼動時間が l 2 の時に比べ. 験を行った。この実験におけるワーカの台数は 37 台で. て長くなる。レベル l を増やすことで部分木数が飛躍的. ある。. に増加し，マスタがワーカに問題を配ったり，結果を受. 4.2 節で述べたレベル l を調整することで問題規模の大. け取ったりする時間が増加してしまった。また，通信時. きさを調整することが可能である。そこで katsura-14 と. 間は計算開始時に比べて，計算終了時の方がより時間が. economic-15 に対して，l を 2 と 3 の 2 通り，再分割手法. かかっており，economic-15 の l 3 に至っては最後の方で. を行う，行わないの 2 通りを組み合わせて計 4 通りを. は 1 回の通信に対して約 3 秒の通信時間がかかっていた。. 行った。数値実験結果を表 7 に示す。部分木数は探索し. 今回採用した 2 層によるマスタワーカ方式では最初から. た部分木の個数であり，再分割を行った際には部分木が. 部分木を細かく分割するよりも，探索回数が多い部分木. さらに分割されるため，数字が増える。平均，標準偏. のみを動的に再分割する方が望ましいことがわかる。た. 差，最大，最小は各ワーカにおける稼動時間の平均と標. だし，再分割を行うよりもレベル l を上げる方が稼働率. 準偏差，最も長い間稼動していたワーカの稼動時間，最. は良くなっているため，マスタのボトルネックを解消で. も短く稼動していたワーカの稼動時間である。稼働率は. きれば，あらかじめ部分木を細かく分割した方が早く計. 式 (7) で計算した値である。. 算できる可能性を秘めている。. レベル l を 3 にすることで，稼働率は l を 2 のまま再. 図 3 は katsura-14 における各ワーカの稼動状況を示し. — 81 —. NII-Electronic Library Service.

(12) Shonan Institute of Technology. 湘南工科大学紀要第 42 巻第 1 号. 図3. katsura-14 に対するワーカの稼動状況，上: l2 再分割しない，中: l2 再分割する，下: l3 再分割しない. たもので，上から l 2 で再分割手法を適用しないケー. でも過剰に分割されていることを表している。. ス，l 2 で再分割手法を適用したケース，l 3 で再分割. 5．おわりに. 手法を適用しないケースを同じ時間軸で並べたものである。1 つの長方形が 1 つの部分木の探索を意味する。l 2. この論文では，多項式方程式系の複素数を含めた孤立. の再分割手法を適用しないケースにおいては長い長方形. 解を全て列挙するための手法である多面体的ホモトピー. が存在している。これが再分割手法を適用することに. 法について，PHoM [11] を例に概要を述べた。数値実験. よって，こまめに分割されて計算されている様子が見て. を通して，解の検証の必要性や方程式系の変数の個数が. とれる。再分割をすることなく l を 3 にすることでも同. 増えるにつれて計算時間が増加することが見られた。. 様の効果が得られていることは確認できる。レベル l を. そこで PHoM を元にさらに変数の個数が多い問題を解. 2 のままにして再分割を行う場合と l 3 にする場合では，. くために並列計算について検討した。とくにこの論文で. 稼動時間が異なることも図 3 から確認できる。レベル l. は混合体積の計算においての既存の並列計算手法を説明. を 3 にしたときに図 3 が黒くなって長方形が密集してい. し，その問題点を指摘し，その問題点を改善する手法を. る部分が多く見られるのは，分割が必要でない部分木ま. 提案し，数値実験を通して提案したアルゴリズムの有効. — 82 —. NII-Electronic Library Service.

(13) Shonan Institute of Technology. 多面体的ホモトピー法のソフトウェアと並列計算（郡司貴之）性を検証した。他の部分については説明できなかったが，いままでバラバラに提案されていた多面体的ホモトピー法に関する並列計算技術をまとめて，多面体的ホモトピー法の全ての段階に対して並列計算を行えるようにソフトウェア PHoMpara [12] を実装した。数値実験を通して，PHoMpara の有用性についても簡単ではあるが示した。今後の課題として，[18] で発表された混合体積の計算の実装を利用することで，混合体積の計算の高速化をはかり，さらに変数の個数の多い問題の解決を実現することがひとつの課題である。また変数の個数の多い問題となれば，ホモトピー曲線の追跡の安定性を向上させるための改良や関数値の計算の高速化，マスタワーカ方式の並列計算におけるサブマスタの導入なども検討しなければならない。. 謝辞この論文は [11, 12, 13] などの研究をまとめたものです。 [11, 12, 13] を書くにあたって，さまざまな意見や厳しい指導をしていただきました小島政和先生，Sunyoung Kim 先生，藤澤克樹先生に感謝します。参考文献 [1] E. Allgower and K. Georg, Numerical continuation methods, Springer-Verlag, 1990. [2] D. N. Bernshtein, “The number of roots of a system of equations,” Functional Analysis and Appl. 9 (1975) 183– 185. [3] G. Bjorck and R. Fröberg, “A faster way to count the solutions of inhomogeneous systems of algebraic equations, with applications to cyclic n-roots,” Journal Symbolic Computation 12 (1991) 329–336. [4] W. Boege, R. Gebauer, and H. Kredel, “Some examples for solving systems of algebraic equations by calculating Groebner bases,” J. Symbolic Computation 2 (1986) 83–98. [5] B. Buchberger, “Grobner basis: An algorithmic method in polynomial ideal theory”, Multidimensional System Theory, D. Reidel Publishing Compary, (1985) 184–232. [6] Y. Dai, S. Kim and M. Kojima, “Computing all nonsingular solutions of cyclic-n polynomial using polyhedral homotopy continuation methods,” Journal of Computational and Applied Mathematics, 151 1–2, 83–97 (2003). [7] F. J. Drexler, “Eine methode zur Berechnung sämtlicher Lösungen von Polynomgleichungssystemen,” Numer. Math. 29 (1977) 45–58. [8] FRISCO Project: http://www.nag.co.uk/projects/FRISCO. html. [9] T. Gao, T. Y. Li, J. Verschelde and M. Wu, “Balancing the lifting values to improve the numerical stability of polyhedral homotopy continuation methods,” Applied Mathematics and Computation 114 (2000) 233–247. [10] C. B. Garcia and W. I. Zangwill, “Determining all solutions to certain systems of nonlinear equations,” Mathematics of Operations Research 4 (1979) 1–14. [11] T. Gunji, S. Kim, M. Kojima, A. Takeda, K. Fujisawa and T. Mizutani, “PHoM—a Polyhedral homotopy continuation method for polynomial systems,” Computing. 73 (2004) 55–77. [12] T. Gunji, S. Kim, K. Fujisawa and M. Kojima, “PHoMpara—Parallel Implementation of the Polyhedral Homotopy Continuation Method for Polynomial Systems,” Computing. 77 (2006) 387–411. [13] 郡司貴之“多面体的ホモトピー法の並列実装”東京工業大学博士論文 (2007). [14] B. Huber and B. Sturmfels, “A Polyhedral method for solving sparse polynomial systems,” Mathematics of Computation 64 (1995) 1541–1555. [15] S. Kim and M. Kojima, “Numerical stability of path tracing in polynomial homotopy continuation methods,” Computing, 73, 329–348 (2004). [16] T. Y. Li’s web site: http://www.mth.msu.edu/~li/ [17] T. Y. Li, “Solving polynomial systems by polyhedral homotopies,” Taiwan Journal of Mathematics 3 (1999) 251–279. [18] T. Mizutani and A. Takeda, “DEMiCs: A software package for computing the mixed volume via dynamic enumeration of all mixed cells.” Dept. of Math. and Comp. Sciences, Tokyo Inst. of Tech. B-441, (2002). [19] A. Morgan, “Solving polynomial systems using continuation for engineering and scientific problems,” PrenticeHall, 1987. [20] V. W. Noonberg, A neural network modeled by an adaptive Lotka-Volterra system, SIAM J. Appl. Math. 49 (1989) 1779–1792. [21] PoSSo Project: http://posso.dm.unipi.it/ [22] M. Sato, H. Nakada, S. Sekiguchi, S. Matsuoka, U. Nagashima and H. Takagi, “Ninf: A Network based Information Library for a Global World-Wide Computing Infrastructure,” HPCN’97 (LNCS-1225), pp. 491–502, 1997. [23] A. J. Sommese and C. W. Wampler, “The Numerical Solution of Systems of Polynomials Arising in Engineering and Science,” World Scientific, 2005. [24] H. J. Su, J. M. McCarthy, M. Sosonkina and L. T. Watson, “POLSYS_GLP: A Parallel General Linear Product Homotopy Code for Solving Polynomial Systems of Equations,” To appear in the Algorithms section of ACM Trans. Math. Softw. [25] A. Takeda, M. Kojima, and K. Fujisawa, “Enumeration of all solutions of a combinatorial linear inequality system arising from the polyhedral homotopy continuation. — 83 —. NII-Electronic Library Service.

(14) Shonan Institute of Technology. 湘南工科大学紀要第 42 巻第 1 号. [26] [27] [28]. [29]. method,” J. of Operations Society of Japan 45 (2002) 64–82. C. Traverso, The PoSSo test suite examples. Available at http://www-sop.inria.fr/saga/POL/ J. Verschelde, The database of polynomial systems is in his web site: “http://www.math.uic.edu jan/”. J. Verschelde, P. Verlinden and R. Cools, “Homotopies exploiting Newton polytopes for solving sparse polynomial systems,” SIAM J. Numerical Analysis 31 (1994) 915–930. J. Verschelde, “Algorithm 795: PHCpack: A general-purpose solver for polynomial systems by homotopy continuation,” ACM Trans. Math. Softw. 25 (1999) 251–276.. [30] J. Verschelde and Y. Zhuang, “Parallel Implementation of the Polyhedral Homotopy Method,” Proceedings of the 2006 International Conference on Parallel Processing Workshops (2006) 481–488. [31] L. T. Watson, S. C. Billips, and A. P. Morgan, “Algorithm 652: HOMPACK: a suite for codes for globally convergent homotopy algorithms,” ACM Trans. Math. Softw. 13 (1987) 281–310. [32] L. T. Watson, M. Sosonkina, R. C. Melville, A. P. Morgan, and H. F. Walker, “HOMPACK90: A suite of Fortran 90 codes for globally homotopy algorithms,” ACM Trans. Math. Softw. 23 (1997) 514–549.. — 84 —. NII-Electronic Library Service.

(15)