電気学会論文誌 C（電子・情報・システム部門誌）
IEEJ Transactions on Electronics, Information and Systems
Vol.133 No.3 pp.583–595
DOI: 10.1541/ieejeiss.133.583
特集解説
GPGPU を用いた並列アントコロニー最適化法による
2 次割当て問題の高速解法
非会員
筒井
茂義∗a)
Parallel ACO Algorithm on GPU for Fast Solution of QAPs
Shigeyoshi Tsutsui∗a) , Non-member
（2012 年 12 月 8 日受付）
Recently, GPGPU (General Purpose computation on Graphics Processing Units) has become popular with great success, especially in scientific fields such as fluid dynamics, image processing, and visualization using particle methods.
In this paper, we discuss how the GPGPU is used in implementations of parallel ant colony optimization (ACO) for
fast solution of quadratic assignment problems (QAPs). As for the ACO, we use the cunning ant system (cAS) which
is one of the most promising ACO algorithms.
キーワード：進化計算，graphics grocessing units (GPUs)，general purpose computation on GPUs (GPGPU)，アントコロニー
最適化法 (ACO)，タブーサーチ，2 次割当て問題 (QAP)
Keywords: evolutionary computation, graphic grocessing units (GPUs), general purpose computation on GPUs (GPGPU), ant
colony optimization (ACO), tabu search, quadratic assignment problem (QAP)
1.
処理の並列性が高く，並列計算に適した計算手法である。し
まえがき
たがって，GPU に進化計算を高速に並列実行させる研究は
比較的早く行われ，2005 年ごろまでさかのぼる (13) 。2008 年
近年，パソコン (PC) の画面表示機能を担当する GPU
（Graphics Processing Unit）がプログラミング可能になり，
になって，Computational Intelligence on Consumer Games
画像処理だけでなく一般的な計算プログラムも実行可能と
している。そのピーク演算性能は 3TFLOPS（テラ・フロッ
and Graphics Hardware (CIGPU) と呼ばれる研究会が W.B.
Langdon らによって発足された (15) 。CIGPU は，進化計算の
メジャーな学会である IEEE 主催の CEC（Congress of Evolutionary Computation）と ACM 主催の GECCO（Genetic
and Evolutionary Computation Conference）で特別セッショ
プス）を超える。このため GPU は価格性能比が非常によく，
ンまたはワークショップとして毎年開催され，活発な活動
様々な科学技術計算へ GPU を適用する GPGPU (General
を行なっている。また，多くの国内外の学会や学術雑誌で
Purpose computation on GPUs) に関する研究がここ数年盛
んになってきた。このように，PC に最新の GPU を付加す
も GPU による並列進化計算の論文発表が盛んになってき
ることで，一昔前のスーパコンピュータに匹敵する計算機
GPU では，同じ命令からなる処理が「スレッド」として
多数生成され，それらが SIMD (Single Instruction, Multiple Data) 風に並列に実行される。しかし，この並列実行で
は MIMD（Multiple Instruction, Multiple Data）として並
列実行されるマルチコア CPU のように複数のスレッドを
なった。最新の代表的な GPU である NVIDIA 社の GTX
680 は 1536 個の小さいプロセッサを搭載し，数万の「ス
レッド」が同時実行可能な並列計算環境を低コストで提供
た (28) 。
環境が構築できる。
さて，遺伝的アルゴリズム（Genetic Algorithm, GA）に
代表される進化計算（Evolutionary Computation, EC）は，
個体 (Individual) レベル，集団 (Population) レベルなどでの
独立したプログラムとして柔軟に実行することができない
a) Correspondence to: Shigeyoshi Tsutsui.
hannan-u.ac.jp
∗
（〈4・1・2〉項参照）。このため，並列計算に向いている進化
E-mail: tsutsui@
計算といえども GPU の超並列性を効率的に実現するには，
阪南大学経営情報学部
〒 580-8502 大阪府松原市天美東 5-4-33
個別の問題に応じて進化計算のモデル構築や実装法にそれ
ぞれ工夫が必要となる。
Dept. of Management Information, Hannan University
5-4-33, Amamihigashi, Matsubara-shi, Osaka 580-8502, Japan
c 2013 The Institute of Electrical Engineers of Japan.
本稿では，組合せ最適化問題に有効な進化計算の手法で
583
GPGPU を用いた並列アントコロニー最適化法による 2 次割当て問題の高速解法（筒井茂義）
Table 1. Coding in evolutionary computation.
あるアントコロニー最適化法（Ant Colony Optimization，
ACO）に関する筆者の研究
(26) (30)
を紹介し，その ACO を用
いて，組合せ最適化問題の中で最も困難な問題の一つであ
る 2 次割当て問題 (Quadratic Assignment Problem，QAP)
を GPGPU により高速に解く研究 (27) (29) について述べ，進化
計算の GPGPU による並列化の研究の一事例を紹介したい。
2.
アントコロニー最適化（ACO）法
〈2・1〉 進化計算における ACO の位置づけ
進化計算
は，複数の個体（解候補）からなる集団を，対象問題（環
境）における各個体の評価値を用いて，よりよい個体を持
つ集団に進化させることにより解の探索を行うというメタ
ヒューリスティックスの総称である。特に数学的に定式化が
Fig. 1. Basic framework of the evolutionary computation.
困難な問題や，組合せ爆発により厳密解を得ることが困難
な問題の解法に有効である。進化計算の源流は，その名が
示すように生物進化にヒントを得たアルゴリズムとして研
究が始まった進化戦略（Evolution Strategy, ES），進化プロ
グラミング（Evolutionary Programming, EP）および遺伝
的アルゴリズム (Genetic Algorithm, GA) にあり，1960 年
代にまでさかのぼる。これらの研究は 1980 年代の後半まで
は互いに独立に研究されてきた。1980 年代の後半からは，
国際会議を通してこれら相互の研究交流が始まった。また，
1990 年代に入り，従来の生物進化をベースとする手法に加
え，群れの集団行動をモデルとする ACO や粒子群最適化
(Particle Swarm Optimization, PSO) (14) ，さらには GA に統
計的手法を融合する分布推定型アルゴリズム（Estimation
of Distribution Algorithm, EDA）(16) (20) などの手法が，集団
Fig. 2. Tour construction based on pheromone density.
ばれるアルゴリズムである。
AS はアリの行動原理に比較的忠実なアルゴリズムであ
るので，以下では，巡回セールスマン問題 (Traveling Salesman Problem, TSP) の解法への AS の適用を例に，ACO の
ベースの探索手法である進化計算の仲間に加わってきた。
進化計算により問題を解く際の解の表現法（コーディン
グ法）には，(1) バイナリ表現 (1011...11)，(2) 実数値列表
概要を説明しよう。各アリを以下ではエージェントと呼ぶ。
現 (2.11, 2.55, ..., 5.55)，(3) 順列表現 (3, 2, 0,...,4, 5) の 3
各エージェントは，各都市に均等もしくはランダムに配置
つが代表的である。バイナリ表現はいろいろな問題に適用
され，そこを出発点として TSP の巡回路を形成する。この
できる汎用的な方法であり，理論的な研究が最も進んでい
とき，各エージェントは，フェロモン濃度に比例して確率
る。実数値列表現は，最適化設計問題など，実数値を扱う問
的に経路を選択する。2 つの都市 i, j 間の経路のフェロモン
題で最もよく使われている。順列表現は，スケジューリン
濃度を τi j としよう。このとき，一度行った都市は訪問しな
グ問題など組合せ最適化問題で多く使われる。各進化計算
いという TSP の規則に従うとする。この経路選択の過程を
と表現法との関係を Table 1 にまとめた。この中で GA は，
餌行動の際の経路生成過程にヒントを得た探索手法であり，
Fig. 2 に示す。
同図 (a) は，都市 1 から出発するエージェント k を示し
ている。都市にいるエージェント k が次に訪問できる都市
の集合は {2, 3, 4, 5, 6} である。そこで，これらの都市 j
( j ∈ {2, 3, 4, 5, 6}) を選択する確率 pk1 j は，(a) に示されてい
るようになる。(b) は，この確率に基づいて都市 2 が選ばれ
た状況である。都市 2 にいるエージェントが次に訪問でき
る都市の集合は {3, 4, 5, 6} であるので，このエージェント
がこれらの都市 j ( j ∈ {3, 4, 5, 6}) を選択する確率 pk2 j は (b)
多くの組合せ最適化問題に適用され，その有効性が報告さ
に示されているようになる。以下，同様に訪問する都市を
れている (9) (10) 。アリはフェロモンを介して相互のコミュニ
順次決定して TSP の巡回路を完成させる。
各種表現法に対応できる汎用性を有している。 ES，PSO，
DE は主として実数値最適化問題に有効である。本稿で取り
上げる ACO は，順序表現問題に有効な手法である。Fig.1
に進化計算のフレームワークを示すが，各手法の大きな相
違は，新しい解候補の生成法（Variation operators）にあり，
基本的なフレームワークは同じである。
〈2・2〉 ACO の仕組み
ACO は，アリの群れによる採
ケーションを行いながら群れで行動し，ある種の秩序を形
各都市に配置された m 個のエージェントが TSP の巡回
成する。ACO は，この秩序形成過程を探索に用いる。ACO
路を完成させる動作を 1 サイクルとする。このときフェロ
の基本モデルは，Dorigo らによる Ant System (AS) (8) と呼
モン濃度は，次式によって更新される。
584
IEEJ Trans. EIS, Vol.133, No.3, 2013
GPGPU を用いた並列アントコロニー最適化法による 2 次割当て問題の高速解法（筒井茂義）
τi j (t + 1) = ρ · τi j (t) +
m
Δτkij . · · · · · · · · · · · · · · · · · (1)
k=1
ここで，ρ は蒸発係数と呼ばれ，(1 − ρ) はサイクル t と t+1
との間にフェロモンが蒸発する割合を示す。また，Δτkij は
エージェント k により経路 (i, j) に新たに排出されるフェロ
モン濃度である。この値は，エージェント k の巡回路 T k の
長さ Ck が短いほど大きな値となるようにするため，次式
のように Ck の逆数とする。
Δτkij
⎧
⎪
⎪
⎨ 1/Ck if agent k uses edge (i, j) in its tour T k
=⎪
⎪
⎩ 0 otherwise .
· · · · · · · · · · · · · · · · · · · · · · · (2)
Fig. 3. c-ant and d-ant in solving TSP.
初期状態ではすべての経路に同じ濃度 τ0 のフェロモンが
存在すると考える。すなわち，τi j (0) = τ0 とする。Fig. 2 で
説明した選択確率は，ACO では式 (3) のように一般化され
ている。
⎧
⎪
⎪
⎪
⎨
pkij (t) = ⎪
⎪
⎪
⎩
[τ (t)]·[ηi j ]β
ij
[τis (t)]·[ηis ]β
if j ∈ Jk (i)
0
otherwise.
s ∈ Jk (i)
· · · · · · · · · (3)
ただし，Jk (i) は，都市 i にいるエージェント k が訪問でき
る都市（まだ訪問していない都市）のリストであり，β は
フェロモン濃度 τi j と ηi j との重みの度合いを制御するパラ
メータである。ηi j は式 (3) に基づく選択確率にバイアスを
与えるものであり，TSP の場合には経路の選択でより近く
Fig. 4. Colony model of cAS.
の都市を選ぶことが好ましいというヒューリスティックか
ら，ηi j = 1/di j（di j は都市 i, j 間の距離）が使われる。
AS は，アリの群れの行動を比較的忠実にモデル化したア
ルゴリズムであるが，選択圧が弱いため性能的には高々30
都市程度の TSP を解く能力しか有していない。その後，多
くの改良型 ACO アルゴリズムが提案されてきた。その代表
例が Ant Colony System (ACS) (7) と Max-Min Ant System
(MMAS) (22) である。この中でも，MMAS は現在も多くの
ナーアント (donor ant; d-ant) と呼ぶ。
Fig. 3 に c-ant と d-ant との関係を TSP の場合を例にし
て示す。この例では，c-ant は，d-ant の巡回経路 7 → 0 →
1 → 2 → 3 をそのまま拝借している。残りの都市 4, 5 およ
び 6 に関しては，フェロモンの濃度に基づいて経路を選択
する。
応用で用いられている。
コロニーモデル：cAS では，Fig. 4 に示すように，m 個
筆者は先に新しい高性能 ACO アルゴリズムとして，カ
のユニットからなるコロニーモデルを考える。各ユニット
ニングアントシステム（Cunning Ant System, cAS）を提案
(k = 1, 2, · · · , m) は一つのエージェント ant∗k,t から構成され，
経路が保存されている。フェロモン濃度 τi j はコロニー全体
の記憶であり m 個のユニットで共有される。経路形成に関
しては，各ユニットの現在のエージェント ant∗k,t が d-antk,t
となり，c-antk,t+1 が生成される。次に，この c-antk,t+1 が
d-antk,t と比較され，よい方の経路をもつものがそのユニッ
トの次のサイクルのエージェント ant∗k,t+1 となる。サイクル
t における各ユニット k (k = 1, 2, · · · , m) のエージェントに
よるフェロモン放出は，エージェント ant∗k,t によって行わ
れる。フェロモン濃度の更新は，AS と同様である（式 (1)，
(2) 参照）。
cAS では c-ant が経路を生成する際，フェロモン濃度をも
とに生成する部分と d-ant から借用する部分の割合の決定が
重要である。c-ant が τi j (t) を基に生成する部分解のシーケ
ンス長を l s とし，d-ant から借用する部分解のシーケンス長
した (26) (30) 。本稿では，ACO アルゴリズムにこの cAS を用
いるので，次節で cAS について述べる。
〈2・3〉 カニングアントシステム（ cAS）
cAS は，カ
ニングアント（c-ant）と呼ぶエージェントと新しいコロニー
モデルを導入している点が今までの ACO アルゴリズムと
大きく異なっている。
カニングアント（c-ant）
：今まで提案されてきた ACO ア
ルゴリズムでは，巡回路生成は全てフェロモン濃度に基づ
いて行われている。cAS はこれらとは異なり，巡回路を生
成する際，他のエージェントの経路を部分的に借用する。残
りの部分の経路形成には，従来と同様フェロモン濃度を利
用する。cAS では，このようなエージェントを抜け目のな
いアントという意味でカニングアント (cunning ant; c-ant)
と呼んだ。カニングアントに拝借を許すエージェントをド
585
IEEJ Trans. EIS, Vol.133, No.3, 2013
GPGPU を用いた並列アントコロニー最適化法による 2 次割当て問題の高速解法（筒井茂義）
Table 2. Results of cAS. Bestavg are average best solutions over 25 runs and Error indicates average excess (%)
of Bestavg from optimal in 25 runs.
Fig. 5. Algorithm description of cAS.
を lc とすると，lc = n − l s である。ただし，n は問題サイズ
（TSP の場合は都市数）である。ここで，l s の平均値 E(l s )
が E(l s ) = n × γ と決定できる制御パラメータ γ (0 < γ < 1)
Table 3. Results of cAS with LK. Iavg and T avg are average iterations and time in second to find optimal in successful runs. Error indicates average excess (%) from
optimal in 25 runs. PC is Opteron 275 (2.4GHz) processors.
を導入する。cAS の研究では，l s に確率的変動を与えるた
めに，l s の決定には，c-ant の生成ごとに式 (4) に示すよう
な確率密度関数を用いて決定する（詳しくは文献 (26) (30)
を参照していただきたい）
。Fig. 5 に cAS のアルゴリズムの
全体を記述する。
⎧
⎪
⎪
⎪
⎪
⎨
L(l s ) = ⎪
⎪
⎪
⎪
⎩
1−2γ
1−γ
γ
ls
nγ 1 − n
2γ−1
γ
l s 1−γ
n(1−γ) n
for γ ∈ (0, 0.5]
for γ ∈ (0.5, 1].
· · · · · · (4)
以上，cAS の特徴である c-ant とコロニーモデルについ
て述べたが，この方式による効果は以下のようにまとめる
ことができる。
• 新しい解 i の生成において，ユニット i に保存されて
∗
いる解（anti,t
）の n × (1 − γ) の部分（部分解という）
を利用するので，既存解を有効に活用できる。この度
ble 3 には，大きな問題における結果を MMAS と比較した
合いはパラメータ γ で調整できる。
ものを示した。進化計算により一般に大きな問題を解くに
• cAS における選択操作は，Fig. 4 のコロニーモデルで
は，ベースとするアルゴリズムにローカルサーチを組み合
示したように，d-ant と c-ant とのサイズ 2 のトーナメ
わせるのが一般的である。これは ACO においても同様で
ント選択となっている。この両者のは c-ant のカニング
ある。Table 3 では，ACO に対するローカルサーチとして
Concorde TSP solver (Concorde) (5) を用いる。Concorde は，
TSP のローカルサーチとして最も強力なものの一つである
Lin-Kernighan 法（以下，LK 法）(17) の基本アルゴリズムを繰
返して適用するものであり，Chained LK あるいは Iterated
LK (17) (2) と呼ばれる方法で実装されている。なお，この実験
ではパラメータチューニングにより，m = 5 とし，MMAS
の ρ 値は 0.8，cAS の ρ 値および γ 値は，それぞれ，0.5,
0.4 としている。
Table 3 において，Error の結果に着目すると，cAS およ
び MMAS のいずれのアルゴリズムとも非常に小さい値で
あり，Chained LK 単独に比較して，これらのアルゴリズ
ムと LK 法を結合する方法の効果は明白となっている。最
適解発見回数#OPT の結果に注目すると，MMAS は att532
を除くすべての問題において#OPT=25 が得られていない。
これに対して，cAS はすべての問題に対して#OPT=25 が
行動により類似したもの同士である。すなわち，d-ant
と c-ant とは，解構造が平均して少なくとも n × (1 − γ)
が一致している。このように類似しているもの同士に
よるトーナメント選択は，GA において集団の多様性
を維持する手法として有用な Crowding 法 (6) (18) と同様
の効果を有し，コロニーの多様性を維持する効果があ
る。両者の類似の度合いもパラメータ γ で調整される。
〈2・4〉 cAS の性能
TSP を用いて cAS と有力な ACO
アルゴリズムとされている MMAS および ACS との性能を
比較した結果を Table 2 に示す。実験条件は文献 (22) およ
び (11) と同じとする。最大経路生成回数は，k ×n×10000 と
し，対称 TSP (STSP) のとき，k = 1，非対称 TSP (ATSP)
のとき，k = 2 としている。また，蒸発係数は ρ = 0.98，
エージェント数は m = n である。cAS のパラメータ γ 値
は 0.4 である。結果の数値はいずれも 25 回の実行の平均を
示している。同表において Bestavg はベスト解の平均値であ
り，Error は，最適解からの誤差（%）を示す。cAS は d198
得られている。また，最適解を得るのに必要とした時間の
平均 (T avg ) も比較対象のアルゴリズムの中で最小となって
いる。
を除くすべての問題で最も良い結果を示していることがわ
〈2・3〉節で述べたように，cAS では，パラメータ γ の値
かる。
を適切に設定することが重要となる。Fig. 6 および Fig. 7
Table 2 は比較的小さい問題での比較結果であるが， Ta586
IEEJ Trans. EIS, Vol.133, No.3, 2013
GPGPU を用いた並列アントコロニー最適化法による 2 次割当て問題の高速解法（筒井茂義）
Fig. 6. Variation of Error for various values of γ.
Fig. 8. A simple example of QAP.
Fig. 7. Variations of Error for various values of γ on
fl3795 and rl5934.
に TSP における γ の影響を示す。問題にも依存するが
γ ∈ [0.2, 0.5) が好ましい値であることが，これまでの研
究によって得られている。
3.
Fig. 9. c-ant and d-ant in solving QAP. n = 10, γ = 0.4,
l s =4, and lc = n − ls = 6 are assumed.
2 次割当て問題への ACO の適用
ACO のヒントとなったアリの行動モデルからも分かるよ
うに ACO と TSP とは，経路探索という意味で直接的な関
係があるので，前章までは，TSP を用いて ACO の解説を
試みた。事実 ACO の改良の研究には TSP が用いられるこ
とが多い。しかし ACO は，多くの NP 困難な組合せ最適
2 次割当て問題は，最適割当て問題を一般化した形式になっ
ているので，多くの割当て問題にも適用できる。
QAP の簡単な例を Fig. 8 に示す。割当て状況は順列 φ で
表現される。Fig. 8 の φ = {1, 0, 3, 2} は，部門 0 を場所 1 に，
部門 1 を場所 0 に，部門 2 を場所 3 に，部門 3 を場所 2 に
それぞれ割当てることを示す。同図で f (φ) の値は式 (5) か
ら， f (φ) = 1524 である。QAP はこの f (φ) が最小となる φ
化問題に応用され有用な結果が報告されている (9) (10) 。本章
では，本稿で対象としている 2 次割当て問題（QAP）への
cAS の適用法について述べたい。
〈3・1〉 2 次割当て問題（QAP）
2 次割当て問題
（QAP）は，n 個からなる部門を n 個の場所に，式 (5) で
を求める問題である。
〈3・2〉 QAP における cAS の構成 (31)
〈2・3〉節では
TSP をベースに cAS を述べたが，ここでは，QAP におけ
る cAS について，TSP との相違点を述べる。TSP では，フェ
ロモン濃度 τi j は，各都市間の経路 (i, j) の選択における好ま
定義される関数値が最小になるように割当てを決定する組
合せ最適化問題である。
f (φ) =
n−1
n−1 ai j bφ(i)φ( j) · · · · · · · · · · · · · · · · · · · · · · · · (5)
しさの度合いとして定義され，巡回路生成は連続したシーケ
i=0 j=0
ンスとして式 (3) に基づいて行われた。これに対して，QAP
ここで，A = (ai j ) および B = (bi j ) はそれぞれ n × n のマ
のフェロモン濃度 τi j は，部門 i を場所 j に割当てる好まし
トリックスであり，φ は {0, 1, · · · , n − 1} の順列である。マ
さの度合いと定義される。また，ηi j として適切なヒューリ
トリックス A と B は，それぞれ，場所 i, j 間の距離，部門
i, j 間の流量（物流や人的交流の強さ）を表している。QAP
は式 (5) からもわかるように，評価関数が距離と流量との
積になっているため，TSP に比べてはるかに解くことが困
スティックがないので式 (3) において，ηi j = 1 とされる。
Fig. 9 は QAP における c-ant と d-ant との関係の例示で
ある。TSP では，Fig. 3 に示したように c-ant は d-ant の順
回路から連続した部分順回路を借用した。QAP では φ の各
ノードの隣接関係と f (φ) とには直接的な関係がないが，φ
の各位置におけるノードがどのように配置されるかが f (φ)
に関係する。したがって，c-ant が d-ant から借用する方法
は，順列表現における各位置のノード値である。Fig. 9 の
難であり，アルゴリズムの検証のベンチマーク問題にもよ
く用いられる。また，実際の応用問題も多くある。例えば，
大きなビルにおける部門の最適配置問題や，グローバル企
業における事業所立地の最適配置問題などである。その他，
587
IEEJ Trans. EIS, Vol.133, No.3, 2013
GPGPU を用いた並列アントコロニー最適化法による 2 次割当て問題の高速解法（筒井茂義）
Fig. 11. A neighbor φ of φ.
Δ(φ, r, s) = arr (bφ(s)φ(s) − bφ(r)φ(r) )+
ars (bφ(s)φ(r) − bφ(r)φ(s) )+
a sr (bφ(r)φ(s) − bφ(s)φ(r) )+
a ss (bφ(r)φ(r) − bφ(s)φ(s) )+
⎛
⎜⎜⎜ akr (bφ(k)φ(s) − bφ(k)φ(r) )+
⎜⎜
n−1
⎜⎜⎜⎜ aks (bφ(k)φ(r) − bφ(k)φ(s) )+
⎜⎜⎜
k=0,kr,s ⎜
⎜⎜⎜ ark (bφ(s)φ(k) − bφ(r)φ(k) )+
⎝
a sk (bφ(r)φ(k) − bφ(s)φ(k) )
Fig. 10. Pseudo code of Tabu Search.
例では，c-ant は部門 0，2 および 4 のノード値（場所）を
d-ant から借用している。なお，借用するノード数 lc が与
⎞ · · · · · · · · · · (6)
⎟⎟⎟
⎟⎟⎟
⎟⎟⎟
⎟⎟⎟
⎟⎟⎟
⎟⎠
ここで，もし Δ(φ, r, s) の移動コストを記憶する n × n
えられた時，借用するノードの場所はランダムに決定する。
残りの部門 1 および 3 のノード値（場所）は TSP の場合と
個の記憶領域を用いると，φ から u 番目の要素と v 番目
同様フェロモン濃度 τi j に基づいて決定される。
の要素を交換して得られる近傍解への移動コスト計算は，
チとして 2-opt 法がよく知られている。一方，タブーサー
{u, v} ∩ {r, s} = ∅ が満たされる {u, v} に対しては次式のよう
になり，その計算量は O(1) と高速に行える (23) 。本稿の研究
チ (Tabu Search, TS) (12) は，それ自体が組合せ最適化問題の
では，この方式を用いる。
〈3・3〉 タブサーチの結合
QAP におけるローカルサー
解法に適用される強力なメタヒューリスティックスの一つ
Δ(φ , u, v) = Δ(φ, u, v)+
(aru − arv + a sv − a su )×
(bφ (s)φ (u) − bφ (s)φ (v) + bφ (r)φ (v) − bφ (r)φ (u) )+ · · (7)
(aur − avr + avs − aus )×
(bφ (u)φ (s) − bφ (v)φ (s) + bφ (v)φ (r) − bφ (u)φ (r) )
であるが，進化計算と組合せてローカルサーチとしてもよ
く用いられる (22) 。
本稿では，TS を ACO と組合せ，そのローカルサーチに
用いる。QAP における TS は 2-opt 法によく似たアルゴリ
ズムであるが，2-opt 法では常に最良近傍に移動させるの
に対して，TS では各繰返しにおいて現在解の最良近傍解
4.
の値が現在解の値より悪くてもその最良近傍解に移動させ
GPGPU による ACO の並列化方式
〈4・1〉 GPGPU の概要
良近傍になる場合が当然起こりうる。このような場合，現
ACO の GPGPU による並列
化方式について述べる前に，GPGPU の概要を簡単に述べ
在解が移動元に戻ってしまうことになる。TS はその名のと
ておく。
る。この場合，移動先の近傍探索において移動元の解が最
〈4・1・1〉 GPU アーキテクチャ
おり，そのような近傍への移動をある期間（禁断期間）禁
GPGPU に用いられ
止することを特徴とする。移動が禁止される近傍（禁断近
る GPU は NVIDIA 社の GPU である。同社のアーキテク
傍）は「禁断リスト (Tabu list)」と呼ばれるデータ構造で
チャはほぼ 3 年単位に世代交代が行われており，2012 年に
管理される。
なって Kepler アーキテクチャーが発表された。2012 年第
4 四半期に発表された Tesla シリーズ（GPGPU 専用の機種
で画像表示インターフェースを持っていない）の K20 は，
倍精度浮動小数点演算能力が 1.17 TFLOPS，単精度浮動小
数点演算能力が 3.52 TFLOPS のピーク性能を有している。
しかし K20 の価格は約 30 万円とやや高価である。画像処
理インターフェースを有する GeForce GTX シリーズは，倍
Fig. 10 に本稿で用いた TS の擬似コードを示す。禁断近
傍でもある基準（アスピレーション基準）を満たせば移動
を許す場合がある。ここでは標準的な基準に倣い，最良解
が得られた場合に移動を許すことにした。また，Tillard の
Ro-TS (24) にならって，禁断期間はランダム要素を取入れ，
T S list size × r3 とした。ただし，r は [0,1) の一様乱数であ
り，T S list size はパラメータである。
Fig. 10 からわかるように，TS を QAP に適用する場合，
現在の解のすべての近傍 N(φ) への移動による式 (5) の変化
量（以下，移動コスト）を計算しなければならない。近傍
精度浮動小数点演算の機能は低いが，進化計算用には十分
な高性能を有しているので，これを用いると数万円の価格
で高性能な GPGPU 環境を構築することができる。Kepler
アーキテクチャに基づく GTX 680 は，単精度浮動小数点
としては，Fig. 11 に示すように φ の 2 つの位置の値を交換
演算能力が 3.09 TFLOPS のピーク性能を有している。価
したものとする。
格は約 5 万円である。
φ を φ の r 番目の要素と s 番目の要素を交換して得られ
た近傍解とすると，移動コスト Δ(φ, r, s) = f (φ ) − f (φ) の
計算量は，以下のように O(n) となる。
キテクチャに基づく GTX 480 である。以下では，そのアー
本稿で用いる GPU は，Kepler の 1 世代前の Fermi アー
キテクチャーを Fig. 12 を用いて説明する。GPU 内におけ
588
IEEJ Trans. EIS, Vol.133, No.3, 2013
GPGPU を用いた並列アントコロニー最適化法による 2 次割当て問題の高速解法（筒井茂義）
Fig. 13. CUDA programing model.
Fig. 12. An example of GPU archtecture (GTX 480).
Table 4. Specifications of typical GeForce GTX GPUs.
では，NVIDIA 社が GPU 計算向けに提供している C 言語
〈4・1・3〉
節参照）と呼ばれ，32 個の CUDA core が一つのグ
を拡張した統合開発環境 CUDA（Compute Unified Device
Architecture）が最もよく利用されている。Fig. 13 に CUDA
のプログラミングモデルを示す。CUDA プログラムでは，
スレッドはグリッド（grid）レベルとブロック（block）レ
ベルの 2 階層の並列構成をとる。ブロックは，スレッドの
集合であり，1 次元， 2 次元，または 3 次元構成をとるこ
とができる。一方，グリッドはブロックの集合であり，1 次
元または 2 次元構成をとることができる。
各スレッドは，カーネル関数（kernel function）に記述さ
ループとしてストリーミングマルチプロセッサ（Streaming
れた同じコードを実行する。スレッドのスケジューリング
Multiprocessor, SM）を構成している。GTX 480 は全体で
15 個の SM を有している。SM 内の各 CUDA core は，共
有メモリ（shared memory）を介してデータを共有すること
ができる。各 SM には，64KB の高速オンチップメモリを
持っているが，これを，共有メモリ 48KB と L1 キャッシュ
16KB として利用するか，あるいは，共有メモリ 16KB と
L1 キャッシュ48KB として利用するかが選べる（本並列化
では L1 キャッシュメモリを増やす後者の構成を選択）。一
方，SM 間のデータ共有は，VRAM を介して行われる。各
SM で共通に使われる L1 キャッシュのサイズは，768KB で
ある。GPU 計算では VRAM が GPU のメインメモリとな
る。なお，単精度浮動少数点のピーク性能は 1.345 TFLOPS
である（Table 4 参照）。
〈4・1・2〉 SIMT
プログラムはスレッドとして実
は，ハードウェアにより自動的に行われる。カーネル関数
るプロセッサは，CUDA core（CUDA の名称については
は，通常のデータ引数の他，グリッドとブロックの定義を
引数として持つ。カーネル関数がコールされると，グリッ
ドとブロックの定義にしたがってスレッドが生成され，そ
れらのスレッドが一斉に実行を開始する。
SM のリソースである共有メモリやレジスタ（ローカル
変数はレジスタに割当てられる）は，一つの SM に同時に
割当てられるブロック間で分割して使われる。したがって，
一つの SM に割当てられるブロック数はこれらの制限によっ
て決まり，同時に実行されるスレッド数もこのブロックの
割当て状況に依存する。
〈4・2〉 CUDA アーキテクチャに基づく並列 ACO の全
体構成
QAP を解くための CUDA アーキテクチャに基
づく並列化 ACO の全体構成を，Fig. 14 に示す。今回の実装
行される。ここで注意しなければならないことは，各スレッ
に当たっては，ACO の各ステップの機能は，全て GPU で
ドは，
「ウォープ (warp) 」と呼ばれる 32 スレッドを単位と
実行されるカーネル関数としてコード化した（Fig. 14 には
して，各 SM で SIMD (Single Instruction Multiple Data) 風
生する。GPU におけるこのスレッドの実行法を，NVIDIA
主要な 3 つのカーネル関数，Construct solutions(...)，
Apply TS(...) および Update Pheromone Density(...)
を示している）。CPU は ACO の各カーネル関数を順次呼
出して ACO の繰返し制御を行うのみである。これにより，
CPU と GPU との間では，探索の進行状況および最終解の
では SIMT (Single Instraction, Multiple Threads）と呼んで
データのみを転送すればよい。
に実行されることである。したがって，同一ウォープ内のス
レッドの計算が，分岐などにより相互に大きく異なった処理
を行う場合には，処理待ちに伴うアイドリングタイムが発
いる (21) 。SIMT では，同一ウォープ内では，できるだけ分
したがって，この方式では，CPU と GPU との間の通信
岐処理がないようにプログラムを作成し，アイドリングタ
オーバヘッドは小さい。Table 5 は，上で述べた 3 つのカー
イムを少なくすることが非常に重要となる。
ネル関数のグリッドとブロック（Fig. 13 参照）の構成を示
〈4・1・3〉 CUDA プログラミングモデル
したものである。たとえば，c-ant の生成を行うカーネル関
GPU 計算
589
IEEJ Trans. EIS, Vol.133, No.3, 2013
GPGPU を用いた並列アントコロニー最適化法による 2 次割当て問題の高速解法（筒井茂義）
τ
Fig. 15. Indexing of moves (n = 8).
Table 7. Neighborhood sizes for various problem sizes.
Fig. 14. ACO with TS on a GPU.
Table 5. Three main kernel functions. n is the problem
size, m is the number of agents. As for T total , see Fig 17
in Section〈4・4〉
C とすると C = n(n − 1)/2 である（n は問題サイズ）。TS の
τ
計算時間がアルゴリズム全体の大部分を占める理由は，こ
ρ
れら C 個の移動コストの計算に多くの時間を要するからで
ある。したがって，この処理をブロック内で効率よく並列
実行させることが重要となる。並列化を行う単純な方法は，
各近傍に対して Fig. 15 のように番号を付与し，ブロック
Table 6. Distribution of computation time of cAS in
solving QAP with sequential runs on a CPU.
内でこの番号に対応する C 個のスレッドを並列に実行させ
る方法である。同図において，n = 8，C = 28 であり，現
在解 φ が φ の 2 と 4 の要素を入替えて得られたものと仮
定している。このとき，Δ(φ , u, v) の計算量が O(n) である
ものを太字で示している。その他のものは，計算量が O(1)
である。
しかし，この方法では，次のような問題が発生する。C
個の近傍の内，式 (6) により CO(n) の計算量で移動コスト
の計算ができる近傍数は，CO(n) = 2n − 3 であり，式 (7)
数 Construct solutions(...) では m 個のブロックが並
により O(1) の計算量で移動コストを計算できる近傍数は，
列で実行される。ブロック内では 1 個のスレッドしか使わ
cAS を CPU によりシーケンシャル実行した場合の，各ス
CO(1) = (n − 2)(n − 3)/2 である。Table 7 は，各 n に対す
る C ，CO(1) および CO(n) を示している。これからわかるよ
うに，CO(n) は C の 10%以下であり，問題サイズが大きく
なるにしたがってその割合は小さくなる。各スレッドは 32
スレッドを単位（ウォープ）として SIMD 風に実行される
が（〈4・1・2〉項参照），このようにわずかな割合である O(n)
の計算量のスレッドが，大部分を占める O(1) の計算量のス
テップの処理に要する時間分布を Table 6 に示す。Table 6
レッドと同じウォープ内で混在すると，そのウォープに属す
れていない。また，TS を実行する Apply TS(...) では，
m 個のブロックによって m 個のエージェントの TS が並列
化され，さらに各ブロック内で T total 個（〈4・4〉節参照）の
スレッドによって移動コストの計算が並列実行される。
〈4・3〉 GPU による移動コスト計算の並列化の課題
からわかるように，TS の時間が全体の処理時間の 99.9%以
るスレッドの処理時間は Fig. 16 に示すように O(1)，O(n)
上を占めている。この処理時間の内訳から，高速処理のた
のスレッドとも，アイドリング時間が生じ，処理時間が長
めには TS の効率的な実装が重要であることがわかる。した
くなる。すなわち，式 (7) を用いて移動コスト計算を O(1)
がって以下では，TS の部分（Fig. 14 におけるカーネル関数
Apply TS(...)）を中心にその並列化について述べよう。
カーネル関数 Apply TS(...) では，各ブロック内で TS
の繰返しごとに，現在解 φ の近傍の移動コスト Δ(φ , u, v)
の計算をしなければならない。φ の近傍 N(φ ) のサイズを
で行えるというメリットを生かせない。
もう一つの問題は，一つのブロックにおける最大スレッ
ド数の制限である。Fermi アーキテクチャでは 1 つのブロッ
クにおける最大スレッド数が 1024 である（Table 4 参照）。
Table 7 からも明らかなように，例えば問題サイズ n = 50
590
IEEJ Trans. EIS, Vol.133, No.3, 2013
GPGPU を用いた並列アントコロニー最適化法による 2 次割当て問題の高速解法（筒井茂義）
Fig. 16. Simple assignment of calculations of move
costs to threads in a block. Due to idling time in SIMT,
the parallel computation by thread in a block is inefficient.
Fig. 17. The thread structure in a block in computing
move costs for TS (MATA).
でも近傍サイズ C は 1225 であり，既にこの問題で最大ス
レッドの制限を越えてしまう。
〈4・4〉 移動コスト計算の効率的並列計算法
〈4・3〉節
で述べた問題に対して以下のような工夫を行い実装した。
ドの構成を示す。この方式により，各スレッドのアイドリ
• CO(1) の処理を行なうスレッドと CO(n) の処理を行うス
ング時間を大幅に減少させることができ，式 (7) を用いて
レッドとが，別のウォープとなるようにスレッド番号
移動コスト計算を O(1) で行えるメリットを生かすことがで
を割当てる。
きる。
• CO(n) の計算を行うスレッドには一つの近傍計算のみを
5.
割当てるのに対して，軽い処理である CO(1) の計算を
行なうスレッドには，NO(1)（NO(1) > 1）個の近傍計算
実験結果
〈5・1〉 実験条件
を割当てる。
本実験で用いた計算機は，インテ
示したように番号 k（k = 0, 1, 2, · · · , C − 1) を付与する。
ル社の Core i7 965 (3.2 GHz) プロセッサと NVIDIA 社の
GTX480（Table 4 参照）を 1 個搭載した PC である。OS
は Windows XP Professional で CUDA プログラムのコンパ
イルには，Microsoft Visual Studio 2008 Professional Edition（最適化オプションは/O2）および CUDA SDK 4.0 を
あるスレッド番号 t = k/NO(1) には，k が t × NO(1) から
用いた。
これを実現する方法には各種の方法が考えられるが，こ
こでは処理オーバヘッドを伴わない以下のような方法を用
いる。まず，交換ペア（近傍）(u, v) に対して，Fig. 15 に
t × NO(1) + NO(1) − 1 までの近傍の計算を割当てる（O(1) の
処理の数は (n − 2)(n − 3)/2 であるが，実装の容易さから C
個として，k が O(n) となる近傍計算の場合には何も計算し
ないこととする）。したがって，O(1) の近傍計算が割当て
られるスレッドの総数は T O(1) = C/NO(1) となる。
O(n) の近傍計算は，ウォープのサイズが 32 であるこ
とから，32 単位の区切れとなるスレッド番号 T O(n) start =
T O(1) /32 × 32 から始まる T O(n) = 2n のスレッドに，O(n)
の処理を割当てる（O(n) の個数は 2n − 3 であるがここで
は先の場合と同様，実装の容易性から 2n 個のスレッドを使
い，不要な組合せの 3 個のスレッドは何も実行しない）。
テスト問題には QAPLIB ベンチマークライブラリ (3) の問
題を用いる。QAPLIB のテスト問題は，(i) ランダム生成問
題，(ii) グリッド距離ベースランダム問題，(iii) 実問題，(iv)
実問題に似せて生成した問題，の 4 つのクラスに分けられ
る (22) 。このうち，本実験では，Taillard によるクラス (i) お
よび (iv) の問題を用いる。Taillard のクラス (i) の問題には
tai*a シリーズ，クラス (iv) の問題には tai*b シリーズがあ
るが，ここでは問題サイズが 40 から 150 までの問題，すな
わち，tai40a，tai50a，tai60a，tai80a，tai100a（以上クラス
(i) の問題）および tai50b，tai60b，tai80b，tai100b，tai150b
（以上，クラス (iv) の問題）の 10 個の問題を用いる。クラ
ス (i) のランダム生成問題は極めて困難な問題として知ら
れ，またクラス (iv) の実問題風に生成された問題は，クラ
ス (i) に比べてやや容易である。
ACO により生成された一つの解に対する TS の繰返し
以上から，ブロック内において使用される総スレッド数
は T total = T O(n)
start+2n
となる。このスレッド割当て法を，
本稿では MATA (Move-Cost Adjusted Thread Assignment)
と呼ぶ。Fig. 17 に，MATA による近傍コスト計算のスレッ
591
IEEJ Trans. EIS, Vol.133, No.3, 2013
GPGPU を用いた並列アントコロニー最適化法による 2 次割当て問題の高速解法（筒井茂義）
Table 8. Parameter values (γ is a control paremeter of
cAS in Section〈2・3〉.
6.
複数 GPU モデルへの拡張 (27)
5 章の実験に用いた PC は 4 個の PCIe x16 スロットを有
している。そこで本章では，3 個の GTX 480 を新たに追
加し，4 章で述べた並列 ACO モデルを 4 個の GPU に実
装し，更なる高速化を行った結果について述べよう。複数
の GPU を並列化する方法として，進化計算でよく用いら
れる島モデル（island model）とマスター／スレーブモデル
（master/slave model）を用いる (1) (4) 。
Table 9. Revised results with MATA.
〈6・1〉 複数 GPU モデル
〈6・1・1〉 島モデル（island model）
4 章で述べた
ACO の GPU への実装では，1 個の GPU に 1 つのコロニー
からなる ACO アルゴリズムを並列実装している。ここで
述べる島モデルでは，各 4 個の GPU にそれぞれコロニー
を構築し，4 つの ACO アルゴリズムを並列に実行させる。
この並列実行においては，各コロニーのエージェントをコ
ロニー間で交換する。
回数を IT T S
Fig. 14 に示したように，各エージェントは GPU 上の
VRAM に保存されている。CUDA 4.x 以降の環境では，CPU
を介すことなく “cudaMemcpyPeer(...)” 関数を用いてこ
れらのデータの GPU 間での交換を直接高速に行うことが
できるが，今回は実装の容易性を優先して CPU を介して
行い，どのエージェントをどの GPU に転送するかなどの
判断は CPU で行う方式とした（Fig. 18 参照）。なお，CPU
において 4 個の GPU の制御には OpenMP API を用いた。
10 参照），ACO の最大サイクル数
MAX ×IT T S MAX
は，アルゴリズム中の TS の総繰返し数となる。この実験で
は，IT T OT AL = m × n × 3200 に固定し，アルゴリズムの終
了条件を，TS の総繰返し数が IT T OT AL に達する，または，
最適解が得られたときとした。Table 8 にその他の代表的な
パラメータを示す。これらのパラメータ値は，文献 (29) で
IT ACO
MAX （Fig.
MAX ，とすると IT T OT AL = m×IT ACO
島モデルを構成するには，各種のトポロジーが提案され
ているが，ここでは，以下の 4 つのトポロジーを考える。
(1) 独立実行型島モデル（Island model with independent
runs, IM INDP) このモデルでは 4 個の GPU の各 ACO ア
用いたものをさらにチューニングしている（括弧内の値は
ルゴリズムが互いに干渉することなく実行される。終了条
文献 (29) で用いた値）。実験回数は 25 とした。
件が満たされれば，直ちにすべての ACO アルゴリズムは
〈5・2〉 結
停止する。
果
Table 9 に実験結果を示す。この実
験では，各問題に対して，(1)MATA を用いた GPU 計算，
(2)MATA を用いない，すなわち，ウォープ内で O(1) と O(n)
の混在を許す GPU 計算（以下，non-MATA）および (3) シ
ングルスレッドで行った CPU 計算の 3 種類の実験を行い，
実行時間および最適解からの誤差 (%) の 25 回の実行の平
均を同表に示した。なお，non-MATA では，一つのブロッ
ク内のスレッド数は C/NO(1) としている。
まず，MATA と non-MATA の T avg 実行時間の比を見る
と，tai100a では 5.4 倍，tai100b では 5.5 倍，MATA の方が
高速である。全体の平均で見ると 6.4 倍高速となっている。
これは〈4・4〉節で述べた MATA が有効であることを示して
いる。このように，GPU 計算では問題の性質に応じて並列
(3) リング型島モデル（Island model with ring connected, IM RING) このモデルでは，各 GPU g (g =
0, 1, 2, 3) のコロニーの最良エージェントのコピーを隣の
GPU (g + 1) Mod 4 のコロニーに送る。受取ったコロニー
化の工夫が必要になる。本稿では言及しなかったが，メモ
では，受取ったエージェントがそのコロニー中の最悪エー
リアクセスにおける遅延を最小化するような工夫も必要に
ジェントよりもよい場合に限って，受取ったエージェント
なる。次に GPGPU と CPU 計算の T avg の比を MATA で見
と最悪エージェントとを入替える。
ると，GPU 計算は GPGPU に対して 16.9∼30.5 倍高速に
(4) 多数移動型リングモデル（Island model with elitist
and massive ring connected, IM ELMR) このモデルでは，
まず (2) の IM ELIT と同様，最良エージェント（global best
(2) エリートエージェント配分型島モデル（Island model
with elitist, IM ELIT) このモデルでは，ACO の Iinterval サ
イクルごとに各 GPU のコロニーの最良エージェントのコ
ピーを CPU に送る。CPU はこれらのエージェントから最
良のエージェント（global best agent）を決定し，各 GPU
に送る。受取ったコロニーでは，受取ったエージェントがそ
のコロニー中の最悪エージェントよりもよい場合に限って，
受取ったエージェントと最悪エージェントとを入替える。
なることを示しており，平均すると 23.8 倍の高速化が得ら
れた。
592
IEEJ Trans. EIS, Vol.133, No.3, 2013
GPGPU を用いた並列アントコロニー最適化法による 2 次割当て問題の高速解法（筒井茂義）
Table 10. Results of the island models with 4 GPUs.
Fig. 18. Island model with 4 GPUs.
スピードアップ値を T avg,1 /T avg,4 で測定する。ただし，1 個
の GPU のときの T avg,1 は 4 章で述べた実装方式，すなわち
ACO アルゴリズムの実行はすべて GPU で行われているモ
デルで測定した値とする。
QAP の各問題に対する許容解は，クラス (i) の問題に対
しては，許容解を最適解からの誤差が 0.5%以下の解とし，
クラス (iv) の問題に対しては，tai150b を除いて許容解を
Fig. 19. Master/slave model with 4 GPUs.
最適解とする。tai150b に関しては，最適解からの誤差を
0.2%以下の解とする。島モデルにおけるパラメータ Iinterval
の値は 1 のときが最もよい結果を示したのでこれを用いる。
〈6・2・2〉 島モデルの結果
島モデルの結果を Table 10
に示す。4 個の ACO アルゴリズムを独立に実行する島モデ
ル IM INDP は，最も単純なモデルであるので，ここでの分
析では他の島モデルのベンチマークとして用いる。IM INDP
によるスピードアップ値は 1.2 (tai60b)∼2.4 (tai80a) である。
tai40a を除くすべての問題で，他の島モデルは IM INDP
よりもよいスピードアップ値を示している。Table 10 には，
許容解が得られるまでの ACO のサイクル数の平均（Average IT ACO ）を示しているが，tai40a ではこの値が 1.7 と小
さく，各 GPU の ACO アルゴリズム間のエージェントの交
換の効果が現れる前に終了条件に達している。これが tai40a
agent）を各 GPU のコロニーに分配する。次に各コロニー
からランダムに m × drate 個のエージェントを選び，(3) の
IM RING と同様そのコピーを隣の GPU のコロニーに送る。
受取ったコロニーでは，それらのエージェントは同じくラ
ンダムに選んだ m × drate 個のエージェントとトーナメント
形式で比較され，置き換えが行われる。drate の値としては，
予備実験により 0.5 とした。
〈6・1・2〉 マ ス タ ー ／ ス レ ー ブ モ デ ル（master/slave
model） 〈4・3〉項の Table 6 で示したように本アルゴリ
ズムでは，ローカルサーチとして用いる TS における移動コ
ストの計算が全体の計算時間の大きな割合を占めている。こ
こでのマスター／スレーブモデルは，Fig. 19 に示すように，
ACO アルゴリズムの実行は CPU が行い，4 個の各 GPU は
m/4 個のエージェントに対する TS を担当する。
CPU において ACO アルゴリズムが m 個の解を生成
したとき，まず m/4 個の解を CPU から GPU に転送し，
“Apply TS(...)” カーネル関数を呼出す。このカーネル関
数の実行終了後に結果を GPU から CPU に送り返す。なお，
ここで，m は 4 で割切れなければならないが，割切れない
ときは m = m/4 × 4 と修正した値を用いる。
〈6・2〉 複数 GPU モデルの実験
〈6・2・1〉 実験条件
計算機環境は，
〈5・1〉
節で述べた
ものと同じである（ただし，GTX 480 を 4 個搭載）。この
実験は 4 個の GPU を用いることによる高速化の効果を知
ることが目的であるので，
〈5・1〉
節のアルゴリズム終了条件
と異なる。各 QAP の問題に対して許容解を設定し，その許
容解が得られる時間を測定する。1 個の GPU および 4 個の
GPU を用いて実行した場合の，許容解が得られるまでの実
行時間の 25 回の平均値を，それぞれ T avg,1 ，T avg,4 で表し，
でいずれのモデルでもスピードアップ値が小さくなってい
る原因である。
IM RING と IM ELIT のスピードアップ値は，ほぼ同じ
値を示している。tai80b と tai150b においては，スピード
アップ値が 4 倍を超えるというスーパリニアが観測されて
いる。IM ELIT と IM INDP との比較に関して，t テストの
結果も示したが，この結果からもわかるように，エージェ
ント交換の効果は特にクラス (iv) の問題で明確に示されて
いる。このアルゴリズムでは，クラス (i) の問題の TS の長
さ（IT T S ）は，クラス (iv) の TS の長さよりも大きくとって
いるので（Table 8 参照），クラス (i) の ACO のサイクル数
（IT ACO ）の値はクラス (iv) に比較して小さい（Table 10 参
照）。このためエージェント交換の効果が，クラス (i) の問
題ではクラス (iv) の問題よりも小さくなったと考えられる。
4 つの島モデルの中では，多くのエージェントを交換す
る IM ELMR が最も優れたスピードアップ値を示している。
また，IM ELMR と IM ELIT との t テストの結果からもわ
593
IEEJ Trans. EIS, Vol.133, No.3, 2013
GPGPU を用いた並列アントコロニー最適化法による 2 次割当て問題の高速解法（筒井茂義）
Table 11. Estimation of Speedup.
かるように，ここでも IM ELMR の有利さは，クラス (i) の
Fig. 20. Results of the master/slave model with 4 GPUs.
問題よりもクラス (iv) の問題で顕著に現れている。
いずれの島モデルにおいてもスピードアップ値は，問題ご
とに大きく異なっている。IM INDP でスピードアップ値が
大きい問題では，他の島モデルでもスピードアップ値が大き
くなっている。そこで以下では，IM INDP を用いて，この
ような異なった値が観測される原因についての考察を試み
る。1 個の GPU を用いたときの計算時間の分布を考え，そ
の確率密度関数および確率分布関数をそれぞれ，f (t)，F(t)
とする。ここで， p 個の GPU からなる島モデル IM INDP
Fig. 21. Computation times for various number of
agents.
を考え，その実行時間の確率分布関数を G(t, p) で表すと，
G(t, p) は次式のように得られる。
G(p, t) = 1 − (1 − F(t)) p , · · · · · · · · · · · · · · · · · · · · · · · (8)
で実行したときの，10 回の平均実行時間を示したものであ
p 個の GPU からなる IM INDP の平均実行時間 T avg,p は，
上式の G(p, t) から，
∞
T avg,p =
t · G (p, t)dt · · · · · · · · · · · · · · · · · · · · · · · (9)
る。ここで，エージェント数 m を 1 から 150 まで変化さ
せた。実行時間は m=1 のときの実行時間で正規化してい
る。TS を行うカーネル関数（Apply TS(...)）では，1 つ
のエージェントの処理に 1 つのブロックを割当てている。
0
本稿の研究で用いている GPU（GTX 480）の SM 数は 15
となる。
このように， p 個の GPU を用いることによる島モデ
であるので，m の増加に対して計算時間は，15 間隔で増加
ル IM INDP によるスピードアップ値は，Speedup(p) =
T avg,1 /T avg,p として，F(T )（あるいは f (t)）から求められ
る。各 QAP の問題によりそれらの F(T )（ f (t)）は異なると
していることがわかる。しかし，その増加時間の量は，こ
考えられるので，スピードアップ値は問題に依存して異なっ
よび 46≤ m ≤60 の m の各区間の実行時間の差は小さい。
の 2 つの問題で大きく異なっている。
tai60a では，1 ≤ m ≤ 15, 16 ≤ m ≤30, 31≤ m ≤45 お
たものになる。Table 11 に，いくつかの f (t) の関数型を仮
〈4・1・3〉項で述べたように，SM のハードウェアリソースに
定した場合に得られる Speedup(p) と Speedup(4) の値を示
余裕があるときには，1 つの SM に複数のブロックが割当
す。実際の実行時間の推定を行って更に詳しい解析が必要
てられる。tai60a では，このような現象が起こっている。す
であるが，各問題によりスピードアップ値が異なったもの
なわち，15 を超えるブロックを割当てても十分ハードウェ
になることは，この解析から理解できる。すなわち，T avg,p
アリソース上の余裕がある。
tai60a は問題サイズが 60（n = 60）であり，総エージェ
は 1 個の GPU による各問題の処理時間の確率分布（F(t)
ント数は m = 60 である (m = n,〈5・1〉節および Table 8 参
あるいは f (t)）に依存して決まる。
アルゴ
照)。したがって，Fig. 19 のマスター／スレーブモデルで各
リズムの全体の実行時間に占める TS の計算量が大きいた
〈6・2・3〉 マスター／スレーブモデルの結果
め（Table 6 参照），マスター／スレーブモデルは，理想上
GPU に割当てられるエージェント数は，m/4 = 15 と少な
い。すなわち，この実験の設定では，GPU のハードウェア
のスピードアップ値である 4 を示すと期待した。しかし，
リソースを十分使い切っていない。このため，予想したス
Fig. 20 に示すように，小さい問題サイズ（tai40a, tai50a,
tai60a, tai50b, tai60b) では，スピードアップ値は 4 よりかな
り小さい値を示している。一方，大きな問題サイズ（tai80a,
tai100a, tai80b, tai100b tai150b）では，理想上の値である
4 に近いスピードアップ値が得られている。
この原因について考えよう。Fig. 21 は，〈6・1・2〉項で述
べたマスター／スレーブモデルにおいて GPU を 1 個のみ
用い，tai60a と tai150b において ACO をサイクル数 10 ま
ピードアップ値が得られていない。
一方，Fig. 21 において tai150b (n = 150) では，m の増
加に対して計算時間は，15 間隔で比例的に増大している。
これは，tai150b では，ハードウェアリソースの制約から，
1 つの SM には 1 個のブロックしか同時に割当てられてい
ないということを示している。tai150b では m = 150 であ
るので，Fig. 20 のマスター／スレーブモデルにおいて，各
GPU に割当てられるエージェント数は 150/4 = 37 であ
594
IEEJ Trans. EIS, Vol.133, No.3, 2013
GPGPU を用いた並列アントコロニー最適化法による 2 次割当て問題の高速解法（筒井茂義）
る（この場合，全エージェント数は m = 37 × 4 = 148 と修
systems”, Ph.D. Dissertation, University of Michigan (1975)
（ 7 ） M. Dorigo and L. M. Gambardella: “Ant colony system: A cooperative
learning approach to the traveling salesman problem”, IEEE Trans. on EC,
Vol.1, No.1, pp.53–66 (1997)
（ 8 ） M. Dorigo, V. Maniezzo, and A. Colorni: “The ant system: Optimization by
a colony of cooperating agents”, IEEE Trans. on SMC-Part B, Vol.26, No.1,
pp.29–41 (1996)
（ 9 ） M. Dorigo and T. Stützle: “Ant Colony Optimization”, MIT Press, Massachusetts (2004)
（10） M. Dorigo and T. Stützle: “Ant colony optimization: Overvew and recent
advances”, Handbook of Metaheuristics, 2nd Edition, pp.227–263 (2010)
（11） L. M. Gambardella and M. Dorigo: “Solving symmetric and asymmetric
tsps by ant colonies”, In Proc. the IEEE Int. Conf. on Evolutionary Computation (ICEC’96), pp.622–627 (1996)
（12） F. Glover and M. Laguna: “Tabu Search”, Kluwer, Boston (1997)
（13） Simon Harding (2012) http://gpgpgpu.com/
（14） J. Kennedy and R. C. Eberhart: “Particle swarm optimization”, In Proc. of
the IEEE Int. Conf. on Neural Networks, pp.1942–1948 (1995)
（15） W. B. Langdon (2012) http://www0.cs.ucl.ac.uk/staff/ucacbbl/cigpu/
（16） P. Larranaga and J. A. Lozano: “Estimation of distribution algorithms”,
Kluwer Academic Publishers (2002)
（17） S. Lin and B. W. Kernighan: “An effective heuristic algorithm for the traveling salesman problem”, Operations Research, Vol.21, pp.498–516 (1973)
（18） S.W. Mahfoud: “A comparison of parallel and sequential niching methods”,
In International Conference on Genetic Algorithms, pp.136–145 Morgan
Kaufmann (1995)
（19） O. Maitre, F. Krüger, S. Querry, N. Lachiche, and P. Collet: “EASEA: specification and execution of evolutionary algorithms on GPGPU”, Soft Comput., Vol.16, No.2, pp.261–279 (2012)
（20） H. Mühlenbein and G. Paaß: “From recombination of genes to the estimation of distribution I. Binary parameters”, Parallel Problem Solving from
Nature (PPSN IV), pp.178–187 (1996)
（21） NVIDIA (2010) developer.download.nvidia.com/compute/cuda/3 2 prod/
toolkit/docs/CUDA C Programming Guide.pdf
（22） T. Stützle and H. Hoos: “Max-Min Ant System”, Future Generation Computer Systems, Vol.16, No.9, pp.889–914 (2000)
（23） É. D. Taillard: “Comparison of iterative searches for the quadratic assignment problem”, Location Science, Vol.3, No.2, pp.87–105 (1995)
（24） É. D. Taillard: taboo search tabou qap code (2004) mistic.heigvd.ch/taillard/codes.dir/tabou qap.cpp
（25） TOP500 SUPERCOMPUTER SITES (2012) http://www.top500.org/
（26） S. Tsutsui. cAS: “Ant colony optimization with cunning ants”, Parallel Problem Solving from Nature, pp.162–171 (2006)
（27） S. Tsutsui: ACO on multiple GPUs with CUDA for faster solution of
QAPs”, In Parallel Problem Solving from Nature Part II- PPSN XII, pp.174–
184 (2012)
（28） S. Tsutsui and P. Collet: “Evolutionary computation on graphics processing
units (Natural Computing Series)”, Springer (2013) in press
（29） S. Tsutsui and N. Fujimoto: “ACO with tabu search on a GPU for solving
QAPs using move-cost adjusted thread assignment”, In Genetic and Evolutionary Computation Conference, pp.1547–1554, ACM (2011)
（30） 筒井茂義・cAS：
「カニングアントを用いた ACO の提案」，人工知能
学会論文誌，Vol.22, No.1, pp.29–36 (2007)
（31） 筒井茂義：
「対称型マルチプロセッシングを用いた並列化 ACO によ
る 2 次割当て問題の解法とその評価」，人工知能学会論文誌，Vol.24,
No.1, pp.46–57 (2009)
正される）。各 GPU は，37 個のエージェントの処理（TS）
をハードウェアリソースを十分使い切って実行しているた
め，4 に近いスピードアップ値が得られている。
7.
む す び
本稿では，筆者が提案した ACO アルゴリズムの一つで
あるカニングアントシステム（cAS）を紹介した後，cAS
を用いて組合せ最適化問題の中でも最も困難な問題の一つ
である 2 次割当て問題 (QAP) を取上げ，アントコロニー最
適化 (ACO) 法とタブーサーチ（TS）とを組合せ，GPGPU
により高速に解く方法を紹介した。
GPGPU はエネルギー効率からも CPU 計算よりも有利で
あるといわれている。一般の科学技術計算では，東京工業
大学の TUBAME2 に代表されるように大規模なスーパーコ
ンピュータが GPU を組合わせて開発されている。世界の高
速コンピュータシステムの上位 500 位までを定期的にラン
ク付けしている TOP500 の 2012 年の第 1 位は，クレイ社が
開発した米国オークリッジ国立研究所のスパコン「TITAN」
となっている (25)（ちなみに，2011 年に 1 位であった京は第
3 位に後退）。TITAN は，米 AMD 社の Opteron 約 30 万個
に加え，NVIDIA 社の Kepler アーキテクチャに基づく最新
の GPU「K20x」を約 26 万個使用している。
進化計算の応用には，スケジューリング問題など実時間
で高速に解を得ることが必要な問題が多くある。本稿では，
単一の GPU および 4 個の GPU を用いる並列進化計算の実
装並びにその結果について述べたが，クラスター計算機を
含む各種計算機からなるプラットホームにおいて，CPU ＋
GPU のヘテロジニアス・コンピューティングによる進化計
算のアルゴリズムを GPGPU により実行する研究などが今
後の興味ある課題である。海外では，既にこのような研究
の例が見られる（例えば EASEA (19) ）。GPU の世代交代は，
ほぼ 3 年ごとに進みエネルギー効率の改良を含めて大きな
進歩が続いている。今後は，プログラミング環境を含めて
より汎用的な開発環境が出現すると思われ，大規模な並列
進化計算への展開と実問題への応用が期待される。
文
献
（ 1 ） Enrique Alba: “Parallel Metaheuristics: A New Class of Algorithms”, John
Wiley and Sons (2005)
（ 2 ） D. Applegate, W. Cook, and A. Rohe: “Chained lin-kernighan for large traveling salesman problems”, INFORMS J. on Computing, Vol.15, pp.82–92
(2003)
（ 3 ） R. E. Burkard, E. Çela, S. E. Karisch, and F. Rendl: “QAPLIB - a quadratic
assignment problem library (2009) www.seas.upenn.edu/qaplib
（ 4 ） E. Cantú-Paz: “Efficient and Accurate Parallel Genetic Algorithms”, Kuwer
Academic Publishers (2000)
（ 5 ） Concorde TSP Solver (2005) http://www.tsp.gatech.edu/concorde.html
（ 6 ） K. A. De Jong: “An Analysis of the behavior of a class of genetic adaptive
筒
井
茂
義 （非会員） 1944 年生。1967 年大阪市立大学工学
部電気工学科卒業。1969 年同大学大学院工学研究
科修士課程修了。同年（株）日立製作所入社。1987
年より阪南大学に勤務。現在経営情報学部 経営情
報学科，兼大学院 企業情報研究科 教授。2000 年
4 月–2001 年 3 月イリノイ大学 IlliGAL 客員研究
員。工学博士。主に進化計算の研究に従事。1996
年度人工知能学会論文賞受賞。
595
IEEJ Trans. EIS, Vol.133, No.3, 2013
© Copyright 2024 ExpyDoc
