DEIM Forum 2016 A8-3

GPU を用いたラベル伝搬法によるグラフクラスタリングの高速化
小澤

佑介†

天笠

俊之††

北川 博之††† 筑波大学大学院システム情報工学研究科 〒 305–8573 茨城県つくば市天王台 1–1–1†† 筑波大学システム情報系 〒 305–8573 茨城県つくば市天王台 1–1–1
E-mail: †kyusuke@kde.cs.tsukuba.ac.jp, ††{amagasa,kitagawa}@cs.tsukuba.ac.jp
あらましグラフデータからコミュニティ構造を抽出するための技術として，グラフクラスタリングが注目されている．しかし，近年のグラフデータの大規模化に伴い，クラスタリング処理の高速化が重要な課題となっている．本研
究では， GPU (graphics processing unit) を用いて，グラフクラスタリング手法の一つであるラベル伝搬法の高速化を行なう．GPU は，非常に多くの演算ユニットを搭載しており，高い並列度の処理を得意とする．提案手法では，ラベル伝搬法の各ステップを，GPU 上での高並列処理に適した形式に変換することで高性能化を図る．さらに，実デー
タと合成データを用いた実験により，提案手法の性能を評価する．
キーワード

1. 序

GPU，グラフクラスタリング，ラベル伝搬法

論

グラフデータ分析手法の一種であるグラフクラスタリングが界のグラフデータにおいて，次数分布はべき乗則に従うことが注目を集めている [11]．グラフクラスタリングは，グラフデー

知られており [1]，かつ GPU 上では千以上のスレッドが並列にタから，密に接続した頂点の集合をクラスタとして抽出する技

動作するので，単純に頂点ごとに処理を振り分けるような並列術である．グラフクラスタリングは，Web やソーシャルネット

化では非常に非効率的となってしまうためである．GPU を用ワークから生物学など，多種多様なグラフデータに対する応用

いてラベル伝搬法の高速化を行なった研究 [23] が既に存在するが考えられている [19]．しかし，近年の情報通信技術の進歩にが，負荷分散に関してはあまり考慮されていないという問題点

伴い，グラフデータの大規模化が進んでいるという問題がある．

がある．

例えば，Facebook の月間アクティブユーザ数は，2015 年 3 月

本研究では，既存研究と異なり，GPU 上での負荷分散を考時点で 14.4 億であると報告されている．（注 1） このような巨大慮したグラフクラスタリングの高速化を行なう．まず，ラベル

なネットワークに対してグラフクラスタリングを適用するため伝搬法をベースに，アルゴリズムを GPU 上での高並列処理にには，クラスタリング処理の高速化が非常に重要な課題となる．

適した形式へと変換する．これにより，各スレッドが行なう仕高速なグラフクラスタリングアルゴリズムの一つとして，ラ事量を均一化すると同時に，GPU に適した処理である，配列

ベル伝搬法が存在する [21]．ラベル伝搬法は，グラフの辺数に上の単純な操作のみでクラスタリング処理が可能となる．さら比例する計算量を持ち，他の多くのアルゴリズムより効率的でに，segmented sort や segmented reduce [2] などの，次数のある．また，処理性能だけでなく，クラスタリング精度におい偏りが存在するデータにおいても効率的に並列処理が可能な演ても良い性能を示すという結果が報告されている [25]．さらに，算を活用する．本稿では，実世界のグラフデータと合成データ

局所的な情報のみをもとにクラスタリング処理を行うため，並を用いて評価実験を行なった．その結果，処理性能においては，列化に適しているという利点もある [16]．これは，現在のプロ

既存手法と比較して最大で 117 倍の高速化に達成した．また，セッサのマルチコア化・メニーコア化を考慮すると，処理の高

正解クラスタ情報を用いたクラスタリング精度の評価において速化を行なうにあたって極めて重要な特性であると考えられる．

も，提案する GPU を用いたラベル伝搬法が良い性能を示した．

一方，近年特に注目されている並列処理の技術として，GPU

本稿の構成は以下の通りである．まず，2. 節において，前提(graphics processing unit) コンピューティングがある [18]．

となる知識に関して説明する．その後，本稿で提案する GPU

GPU は，数百から数千の処理ユニットを搭載しており，CPUを用いたラベル伝搬法について 3. 節で述べ，その評価を 4. 節と比較して非常に処理性能が高い．GPU コンピューティング

で行なう．関連研究については 5. 節で説明し，6. 節で結論と今

では，通常は CPU で行われる計算負荷の高い処理を，GPU を後の予定を述べる．

用いて並列処理することで全体の処理の高速化を図る．しかし，

GPU の性能を引き出すためにはいくつかの技術的課題を解決

2. 前 提 知 識する必要がある．特に，GPU 上でグラフデータを処理するに

本節では，まず本研究で扱う問題に関して，2. 1 節で述べる．

あたって，負荷分散が重要な課題となる．これは，多くの実世

2. 2 節では，ラベル伝搬法によるグラフクラスタリングを説明する．最後に，GPU コンピューティングに関して 2. 3 節で述（注 1）：http://investor.fb.com/releasedetail.cfm?ReleaseID=908022

べる．

v0
v1

A
C

D

v2

v3

v4

G

v0

A

v6
v5

B

E

F

(a) 初期状態

v1

A

A

v2

v3

v4

D
D

v0

A

v6
v5

A

v1

D

A

D

v2

v3

D
D
v6

v5

A

(b) ラベルを (v2 , v4 , v6 , v1 , v3 , v0 , v5 ) の順序で更新

v4D(c) 最終状態

図 1: ラベル伝搬法の例：各頂点 vi の側の文字は vi のラベルを示している．

2. 1 問 題 定 義

Algorithm 1: ラベル伝搬法

グ ラ フ を G = (V, E) と 表 す．こ こ で ，V は 頂 点 集 合 ，Input: グラフ G = (V, E)

E ⊂ V × V は辺集合である．また，基本的に，頂点の数を
n = |V |，辺の数を m = |E| と表す．本研究では，簡単化のた

Output: クラスタリング結果 C
1

ラベル割当て L を初期化

めに，重み無し無向グラフを仮定する．ただし，重み付きや有2

repeat向グラフを対象とするような拡張も容易に可能である．本研究

3で扱う問題であるグラフクラスタリングは以下のように定義で

4

きる．

5
6

◃ Li は頂点 vi のラベルを表す

π ← 数列 (1, 2, . . . , n) のランダムな置換
for i ∈ π do!"#!
Li ← arg maxl∈L ! j | (vi , vj ) ∈ E ∧ Lj = l !

end

問題 (グラフクラスタリング). グラフ G = (V, E) が与えられた

7

until 終了条件を満たすとき，以下の条件を満たすクラスタ集合 C = {C1 , C2 , . . . , Cl }8

C ← 同じラベルの頂点を同じクラスタとして抽出を求める：(1) 任意の i についてクラスタ Ci は V の部分集合，!
(2) 任意の i, j (i =| j) に関して Ci ∩Cj = ∅，(3) li=1 Ci = V ，

9return C(4) 同じクラスタ内の頂点は密に接続されており，異なるクラ以上更新が行われない状態となる．そのため，最終的な結果とし

スタ間の頂点は疎に接続されている．

て，二つのクラスタ CA = {v0 , v1 , v2 } と CD = {v3 , v4 , v5 , v6 }2. 2 ラベル伝搬法によるグラフクラスタリング
グラフクラスタリング手法の一つとしてラベル伝搬に基づくが得られる．

2. 3 GPU コンピューティング

アルゴリズムが提案されている [21]．ラベル伝搬法では，あるGPU コンピューティングは，元来グラフィックス処理用プロ頂点が属するクラスタは多くの隣接頂点が属するクラスタと同セッサであった GPU を，汎用的な計算の高速化に用いる技術

じである，という考え方を基にクラスタリングを行う．

のことである [18]．GPU はメニーコアプロセッサとして近年

ラベル伝搬法の擬似コードは Algorithm 1 のように書ける．

注目を集めており，科学計算やデータベース処理 [13]，データ

まず，各頂点に対してユニークなラベルを付与する（1 行目）．

マイニングアルゴリズムの高速化 [24] など，幅広い応用が検そして，各頂点をランダムな順番で調べ，ラベルを更新してい討されてきている．GPU を用いたアプリケーションを開発すく（3–6 行目）．ある頂点のラベルは，隣接頂点のラベルを調る際には，NVIDIA 社製の GPU と，NVIDIA 社が提供する

べ，最も多く出現するラベルへと更新する（5 行目）．出現回

CUDA フレームワーク（注 2） を用いることが主流となっている．

数が同じラベルが存在する場合には，ランダムに更新先のラベ

GPU は複数の Streaming Multiprocessor (SM) から成り，

ルを選択する．ラベル更新のイテレーションを終了条件を満た

各 SM は数十から数百の Scalar Processor (SP) を搭載していすまで継続し，最終的に同じラベルを持つ頂点を同じクラスタる．例えば，NVIDIA Tesla K40 は 15 個の SM を持ち，各 SMとして抽出する（8–9 行目）．終了条件としては，ラベルの更には 192 の SP が含まれている．各 SP は単純な処理のみが可

新が行われなくなる，またはラベル更新のイテレーションを一能な小さなプロセッサであるため，GPU は複雑な条件分岐な

定回数行なう，などが考えられる．

どが含まれる処理を苦手とする．一方，比較的単純な処理を大ラベル伝搬法の例を図 1 に示す．まず，各頂点にユニークなラベルを付与する（図 1(a)）．その後，頂点をランダムな順序量に行うような場合には，多数の SP による並列処理が有効となり，高い性能を発揮できる．

で調べていき，ラベル更新処理を行う．ここでは，更新順序がまた GPU 上には主要なメモリが三種類存在する．一つは，(v2 , v4 , v6 , v1 , v3 , v0 , v5 ) である場合を考える．図 1(a) における

GPU 上で最も容量の大きいデバイスメモリである．例えば，

v2 の隣接頂点のラベルを見ると，A, B, D がそれぞれ 1 回ずつ

NVIDIA Tesla K40 は 12 GB のデバイスメモリを搭載してい出現している．この場合，上述の通り，実際にはランダムに更る．しかし，デバイスメモリは SM の外に存在するため，デバ新ラベルを決定するが，簡単のためにアルファベット順で最小

イスメモリへのアクセスは比較的低速となる．SM 上には，容のラベルに更新することを考える．すなわち，v2 のラベルは A

量は数十 KB と小さいが高速にアクセス可能なオンチップメモ

へと更新される．同様に，v4 のラベルは D へ，v6 のラベルはリがある．このメモリは，L1 キャッシュや SM 内でデータをや

D へと更新されていき，全頂点を更新すると図 1(b) のようになりとりするための共有メモリとして用いられる．また，最も高る．この後，もう一度全頂点を調べると，図 1(c) となり，これ（注 2）：http://www.nvidia.com/object/cuda_home_new.html

速にアクセス可能なメモリとして，各 SM 上にレジスタが存在
する．

0

2

3

4

5

6

7

0 2 4 7 10 13 16 18

id

1 2 0 2 0 1 3 2 4 5 3 5 6 3 4 6 4 5CUDA では上記の階層的なプロセッサやメモリ構成に対応するスレッディングモデルを定義している．まず，GPU 全体に

1

ptr

0

1

2

3

4

5

6

7

8

9

10 11 12 13 14 15 16 17v1 の隣接頂点

対応するものとしてグリッドがあり，グリッドは複数のブロックを含む．ブロックは SM に対応しており，一つのブロックは

図 2: GPU 上のグラフデータ：二つの配列 ptr と id によってある SM に割り当てられ，処理される．一方，一つの SM は複

図 1 中のグラフを表現している．

数のブロックを同時に動作させることが可能となっている．このとき，SM の資源（オンチップメモリなど）は複数ブロック

べる．ラベル伝搬法では，ある頂点のラベルを更新する際，隣で共有して用いられる．一つのブロックは多数のスレッドを持接頂点のラベル（以下，隣接ラベルと呼ぶ）を集計し，最も出

ち，各スレッドの処理は SP によって実行される．

現回数が多い隣接ラベルを求める必要がある．このとき，単純

本節の残りでは，本研究で用いるデータ並列プリミティブにに頂点ごとに GPU のブロックを割り当て並列処理すると，次ついて説明する．データ並列プリミティブは，並列計算におい数の偏りのため，非効率的な並列化となってしまう．そこで，て頻出する処理であり，GPU 向けに最適化された提案・ライ

本研究では，次数の偏りによる影響を受けないような並列化をブラリが存在する [2, 13]．本研究では，これらのプリミティブ

行なうことで，GPU 上で効率的にラベル伝搬法を処理する．

を活用することで，負荷分散を考慮した，高速な GPU 上のラ
ベル伝搬法を実現する．

2. 3. 1 Reduce"#
Reduce は，入力として配列 a0 , a1 , . . . , an−1 と二項演算子 ⊕ を受け取り，スカラ値 a0 ⊕ · · · ⊕ an−1 を返す演算である [13]．例えば，二項演算子として加算 + が与えられた場合，reduce は入力配列中の値の和を返す．
2. 3. 2 Scan

以降では，まず，3. 1 節で GPU 上のデータ表現形式について説明する．その後，隣接ラベルの出現回数を並列に集計する"

方法について 3. 2 節で，ラベル更新処理について 3. 3 節で述べる．

3. 1 GPU 上のデータ表現形式
ラベル伝搬法を処理するためには，以下の二種類のデータを

GPU 上に格納しておく必要がある：(1) 入力グラフ，(2) 各頂#

Scan は，入力として配列 a0 , a1 , . . . , an−1 ，二項演算子 ⊕，"単位元 I を受け取り，配列 I, a0 , a0 ⊕ a1 , . . . , a0 ⊕ a1 ⊕ · · · ⊕#an−2 を返す演算である [12]．Scan は，並列に出力を行なう際のオフセットの計算などによく用いられる．

2. 3. 3 Segmented Sort
Segmented sort は ，入 力 と し て 二 つ の 配 列 A =#"#a0 , a1 , . . . , an−1 と S = s0 , s1 , . . . , sm−1 , sm = n を受け"#取り，A の各部分配列 asi , . . . , asi+1 −1 をそれぞれソートす"る演算である [2]．配列 S の各要素 si は部分配列の先頭から点のラベル情報．GPU 上でこれらのデータを扱うためには，並列処理に適していると同時に，グラフの巨大さを考慮して，省メモリな構造を用いることが望ましい．ラベル情報は，単純に，各頂点のラベルを要素として持つ配列に格納すれば十分である．
なお，本稿では説明の都合上，ラベルを文字で表すが，実際には 0 から n − 1 までの整数値によって表現している．一方，実世界のグラフデータを隣接行列として表した場合，疎行列になると考えられる．そこで，本研究では，疎行列表現形式の一つである CSR (compressed sparse row) [3] を用いてグラフを表現する．

のオフセットを示している．ナイーブな実装として，一つの部

CSR を用いたグラフデータの例を図 2 に示す．一つのグラフ

分配列に一つのブロックを割り当てソートするものが考えられは二つの配列 ptr と id によって表現される．配列 id は各頂点るが，これは負荷分散の観点から見て良い実装とは言えない．

の隣接頂点 ID を格納しており，配列 ptr はある頂点の隣接頂

Modern GPU [2] では，負荷分散を考慮して，各スレッドがほ点 ID が配列 id 中のどこに格納されているかを示している．

ぼ同じ処理量となる実装をしている．

例えば，図 1 中の v1 の隣接頂点が v0 と v2 であるという情報

2. 3. 4 Segmented Reduce

は，配列 id 中の ptr [1] = 2（注 3） から ptr [2] = 4 番目の要素としSegmented reduce は ，入 力 と し て 二 つ の 配 列 A =#"#a0 , a1 , . . . , an−1 , S = s0 , s1 , . . . , sm−1 , sm = n と二項演算"#子 ⊕ を受け取り，A の各部分配列 asi , . . . , asi+1 −1 に対して"#reduce を適用した結果の配列 r0 , r1 , . . . , rm−1 を返す演算でて格納されている．通常，CSR では行列の各成分の値を保持

S の各要素 si は，segmented sort の場合と同様に，部分配列

数を並列に計算する手法について説明する．この処理は，GPUの先頭からのオフセットを示している．Segmented reduce も，上での効率的な並列化のために，以下の 4 ステップに分けて行Modern GPU の一関数として，負荷分散を考慮した高速な実

われる：(1) 頂点 ID から隣接ラベルへの射影，(2) 頂点毎に隣装が公開されている．

接ラベルをソート，(3) 隣接ラベル境界の抽出，(4) 出現回数"ある [2]．ここで，ri = asi ⊕ asi +1 · · · ⊕ asi+1 −1 であり，配列

3. GPU を用いたラベル伝搬法の高速化
本節では，GPU を用いたラベル伝搬法の高速化について述するために配列をもう一つ用いるが，重み無しグラフでは行列中の値が 1 か 0 のみであるため，ここでは省略している．

3. 2 隣接ラベル出現回数の計算
上述のデータ構造をもとに，GPU 上で隣接ラベルの出現回を計算．この 4 ステップへの分割は以下の考え方に基づいてい（注 3）：A[i] は配列 A の i 番目の要素を表す．

0

id

0

L

1

2

3

4

5

6

7

8

9

10 11 12 13 14 15 16 17

1 2 0 2 0 1 3 2 4 5 3 5 6 3 4 6 4 5
1

2

3

4

5

0

1

2

3

4

5

6

2

3

4

5

A =
| B?
→ 1

6

B B A C D C D
0

1

0
7

8

9

10 11 12 13 14 15 16 17

F

8

9

10 11 12 13 14 15 16 17

1

2

B =
| B?
→ 0
3

4

5

6

7

8

9

10 11 12 13 14 15 16 17

1 1 1 1 0 1 1 1 1 1 0 1 1 1 0 1 1 10

図 3: 頂点 ID から隣接ラベルへの射影：配列 L の i 番目の要F′

素は vi のラベルを示している．

+1

S

2

3

0 =
| 1?

4

5

6

7

8

9

10 11 12 13 14 15 16 17

1

4 =
| 4?
2

3

4

5

6

7

8

9

10 11 12 13 14

1 2 3 4 6 7 8 9 10 12 13 14 16 17 18
0

Sptr

1

0 1 2 3 4 4 5 6 7 8 9 9 10 11 12 12 13 140

する同じ隣接ラベルは配列中に連続して格納されている．すな

7

Scan

L′ B A B A B B C A D C C C D C D D D C

る．頂点毎に隣接ラベルがソートされていれば，ある頂点に関

6

L′ A B A B B B C A C D C C D C D D C D

1

2

3

4

5

6

7

0 2 4 6 9 11 13 15

わち，配列中で隣接ラベルが変わる位置（隣接ラベル境界）の

図 4: 隣接ラベル境界の抽出

情報があれば，同じ隣接ラベルが何回出現するかの計算が，特別にデータ構造を用意しなくとも，配列上の簡単な操作のみで可能となる．以降では，各ステップの具体的な処理方法について説明していく．

0

Sptr

0

W

1

2

3

C, D, C, D である場合を考えている．このステップは，単純

6

7

1

2

3

4

5

6

7

8

9

10 11 12 13 14

Segmented reduce

図 3 に示しているように，隣接頂点 ID の配列 id と現在のする．この例では，v0 から v6 のラベルが，それぞれ B, B, A,

5

1 1 1 1 2 1 1 1 1 2 1 1 2 1 1

3. 2. 1 隣接ラベルへの射影
ラベル情報 L をもとに，各頂点の隣接ラベルの配列 L′ を作成

4

0 2 4 6 9 11 13 15

0

1

2

3

4

5

6

Wmax 1 1 2 1 2 2 1
0

I

1

2

3

4

5

6

1 3 4 8 9 12 14

0

1

2

3

4

5

6

L B B B D C D D

図 5: ラベルの更新処理に，1 スレッドが配列 id と L の各要素を調べ，出力を行うことで高速に並列処理可能である．具体的には，i 番目のスレッ
ドは L′ [i] ← L[id[i]] という処理をする．

3. 2. 2 隣接ラベルのソート

次に，前ステップで得られた隣接ラベルの配列 L′ を，頂点毎にソートする．すなわち，頂点 vi に関する部分配列をそれ
ぞれ別々にソート対象として処理する．これは，2. 3. 3 節で説明した，segmented sort を利用することで効率的に実現可能である．具体的には，入力配列 A として隣接ラベルの配列 L′ を，S として ptr を渡して segmented sort を実行することで達成
できる．

3. 2. 3 隣接ラベル境界の抽出
ソートされた隣接ラベル配列をもとに，隣接ラベル境界の抽出を行なう（図 4）．このステップは，以下の三つの処理に分けて行われる：(1) 隣接ラベル境界判定，(2) 隣接ラベル境界
出力オフセットの計算，(3) 隣接ラベル境界抽出．隣接ラベル
境界判定は，隣接ラベルがどの位置で変わるかを示すフラグの
配列 F を作成する．ここでは，配列 L′ の 1 要素に 1 スレッドを割り当て並列に処理する．具体的には，i 番目の要素を担当するスレッドは，L′ 中の右隣りの要素を調べ，配列 F の i 番
目にフラグを格納する．すなわち，L[i] =| L[i + 1] ならば F [i]に 1 を，そうでなければ 0 を格納する．また，ラベルが等しい場合でも，異なる頂点に関する隣接ラベルである場合には，Fに 1 を格納する．例えば，図 4 では，L′ [3] = L′ [4] であるが，F [3] = 1 としている．
次に，フラグ配列 F を用いて，隣接ラベル境界出力オフセッ
トの計算を行なう．そのためには，2. 3. 2 節で説明した，scan
プリミティブを用いる．入力として，配列 F ，加算 +，単位元0 を渡すことで，図 4 の配列 F ′ が出力される．

この配列をもとに，隣接ラベル境界の配列 S を作成する．隣
接ラベル境界判定の場合と同様に，配列 F ′ の 1 要素に 1 スレッ
ドを割り当て並列に処理する．このとき，i 番目の要素を担当するスレッドは，F ′ [i] と F ′ [i + 1] を比較し，異なればラベル境界を出力し，そうでなければ何も出力せずに処理を終える．
出力を行う場合，i 番目のスレッドは，配列 S 中の F ′ [i] 番目に

i + 1 を書き込む．すなわち，S[F ′ [i]] ← i とする．また，同時に，配列 S において隣接ラベル境界がどの頂点に関するものかを示すために，配列 Sptr も計算しておく．この配列は，グラ
フデータにおける ptr と同様の役割を果たしている．Sptr は，

Sptr [i] ← F ′ [ptr [i]] とすることで並列に計算可能である．
3. 2. 4 出現回数の計算

前ステップで計算した配列 S を用い，隣接ラベル出現回数を持つ配列 W を計算する．ここでも，配列 S の 1 要素に 1 ス
レッドを割り当て，並列に処理する．W [i] は，配列 S 中の隣
接する 2 要素の差，すなわち W [i] ← S[i + 1] − S[i] となる．
配列の境界部分においては（例えば，図 4 の配列 S から W [0]を計算する場合），S[i] がそのまま出現回数となる．

3. 3 ラベル更新処理
3. 2 節で計算した結果を用いて，実際にラベルの更新を行なう．そのためには，隣接ラベル出現回数の配列 W から，各頂点毎に，最も出現回数が多い隣接ラベルを求める必要がある．
この処理は，2. 3. 4 節で説明した，segmented reduce を利用することで，高速に実現可能となる．入力配列として W と Sptr ，また二項演算子として二つの値を受け取り大きい値を返す関数

max を渡すことで，各頂点毎に隣接ラベルの最大出現回数を求めることができる．本来，segmented reduce では，reduce 結Processing Time (s)Louvain

1.8
1.6
1.4
1.2
1.0
0.8
0.6
0.4
0.2
0.0

GPU-LP

7

2.5
1.53

CPU-LP

120
6.20

2.10

2.0

100

150

5
1.5

80

4

0.79

60

1.10

40

2
0.5

(a) amazon

0.03

0.08

0

(b) dblp

83.49

30.46

50

20

1

0.0

100

2.78

3

1.0

0.03

200 192.75
103.756(c) youtube

3.32

0.88

0

0

(d) livejournal(e) orkut

図 6: クラスタリング処理時間

表 1: 使用した実データセットnmamazon334 863925 872dblp

317 080

1 049 866

youtube1 134 890

2 987 624

livejournal

3 997 962

34 681 189

orkut

3 072 441

117 185 083

Processing Time (s)データセット名

CPU-LP

8
7
6
5
4
3
2
1
0

GPU-LP

25
20
15
10
5
1

2

果の配列のみを計算していたが，それだけでは実際にどのラベ

3

4

5

6

7

8

9 10

0

1

2

3

4

5

6

7

Iteration

Iteration(a) livejournal(b) orkut

8

9 10

図 7: 各イテレーションにおける処理時間

ルに更新すれば良いかが判別できない．そのため，本研究では，

segmented reduce を拡張し，各 reduce 結果の値が入力配列のどの位置から来たものかを示す配列 I も返すようにした．この

速かつ高いモジュラリティの結果を出力できることで知られて結果を用いることで，最終的に，ラベルの更新処理が可能とな

いる [19]．Louvain 法の実装には，著者らによって公開されてる．具体的には，L[i] ← L′ [S[I[i]]] を計算することで，並列に更新可能である．

4. 評 価 実 験いるものを利用した．（注 4） なお，本実験ではラベル伝搬法の終了条件を，「ラベルが更新された頂点数が 10−5 n 以下になる，またはラベル更新処理のイテレーションを 10 回行なう」と設定した．

本節では，3. 節で提案した，GPU を用いたラベル伝搬法の性

4. 2 実データによる評価

能を評価する．まず，評価実験に用いた環境と比較手法を 4. 1

本節では，実データを用いたグラフクラスタリング処理時

節で説明し，その後，4. 2 節において実データを用いたグラフ

間の評価を行なう．具体的には，CPU 側のメモリ上にグラフ

クラスタリング処理時間の評価を，4. 3 節で合成データを用いデータが準備できている状態から，クラスタリング結果を計たグラフクラスタリング処理時間とクラスタリング精度の評価

算し終えるまでの時間を計測した．提案手法 GPU-LP では，を行なう．

CPU と GPU 間のデータ転送時間も含めている．データセット

4. 1 実 験 設 定

は，Stanford Network Analysis Project [14] で公開されてい実験には，CentOS release 6.5 (Final) が動作するコンピュー

るデータセットの一部を利用した．用いたデータセットを表 1

タを用いた．CPU として Intel Xeon Processor E5-2687W v2にまとめる．これらは，全て実世界のソーシャルネットワーク(3.40 GHz) を，メモリは 128 GB を搭載している．また，GPUやウェブグラフを表している．

は NVIDIA Tesla K40 (15 SM, 2880 SP, 12 GB) を備えている．

図 6 に，五つのデータセットを用いた場合のクラスタリング
処理時間を示す．全体としては，GPU-LP が最も高速で，次提案の GPU によるラベル伝搬法 (GPU-LP ) は，CUDA 7.0に CPU-LP ，Louvain という順の性能になっている．この結果を用いて実装した．比較手法として，ラベル伝搬法の CPU にから，GPU-LP は GPU を用いることで，大幅なクラスタリ

よる逐次処理 (CPU-LP ) と近年注目されているグラフクラス

ング処理の高速化に成功していることがわかる．具体的には，

タリングアルゴリズムの一つである Louvain 法 [4] を用いた．

Louvain と比較して，amazon と orkut と用いた場合には約 55

Louvain 法は，グラフクラスタリング結果の評価指標であるモ
ジュラリティ[17] を，貪欲法により最適化する手法であり，高（注 4）：https://sites.google.com/site/findcommunities/

CPU-LP

GPU-LP

1.0

1.0

0.8

0.8

0.8

0.6

0.6

4
3

0.6
0.4

2

0.2

1
0

1

2

3

4

5

6

7

8

9 10

Number of vertices (×105)(a) クラスタリング処理時間

0.0

ARI

1.0

5

F-measure

6

NMI

Processing Time (s)Louvain

0.4
0.2

1

2

3

4

5

6

7

8

9 10

Number of vertices (×105)

0.0

0.4
0.2

1

2

3

4

5

6

7

8

0.0

9 10

Number of vertices (×105)

1

2

3

(c) F 値(b) NMI

4

5

6

7

8

9 10

Number of vertices (×105)(d) ARI

図 8: 合成データを用いた実験結果

倍の高速化，dblp と youtube を用いた場合には約 70 倍の高速

は 0 から 1 の値を取り，もっとも良い場合，すなわち正解と同化に成功している．livejournal の場合（図 6(d)）には，特に高じクラスタが得られた場合に 1 となる．

速化の度合いが大きく，117 倍の高速化を達成している．

GPU-LP と CPU-LP を比較すると，GPU-LP が 25 倍から 36 倍の処理性能となった．具体的には，三つのデータセッ
ト dblp, youtube, livejournal の場合には，35 倍前後の高速化

F 値は二つのクラスタリング結果間のクラスタのマッチングに基づき，類似度を求める指標である：

F(C, C ′ ) =に達成した．一方，今回使用したデータセット中で最も巨大な

2|Ci ∩ Cj′ |1 %|Ci | max′n C ∈CCj′ ∈C ′ |Ci | + |Cj |iorkut の場合には，25 倍の性能向上と，他のデータセットにおF 値は，NMI と同様に，0 から 1 の値を取り，最良の場合に 1

ける高速化度合いよりは多少劣る結果となった．これは，図 7となる．

に示すように，CPU-LP はイテレーションが進み，ラベル数が少なくなるに連れて処理時間が大きく減少するためである．図から分かるように，livejournal の場合と比較して，orkut の場
合の方が処理時間の減少度合いがより急であるため，CPU-LP

ARI はペアの数え上げによって，クラスタリング結果を評価する指標である：′ARI(C, C ) =が相対的に短い時間で処理を完了した．一方，GPU-LP では，負荷分散の達成のためにラベル数に依存しない実装としている
ため，イテレーション毎の処理時間はほぼ一定であり，概ねグ
ラフの辺数に比例した時間となる．以上のことから，orkut の場合には，25 倍の性能向上に収まったと考えられる．

1
2

$

($ &
i

i,j

Ni·
2

'

&Nij '−M$ &N·j ')

2

+

j

2

−M

ここで，N は NMI の場合と同様の混同行列であり，M ="$ &Ni· ' $ &N·j '# &n'/ 2 である．ARI も，他の二つの指標とij22
同様に，最良の場合に 1 となる．

LFR ベンチマークによる合成データを用いた実験結果を図 8

4. 3 合成データによる評価に示す．データ生成の際，パラメタは平均次数を 20，最大次数

本節では，クラスタリング処理時間だけでなくクラスタリン

を 50 とし，頂点数を 10 万から 100 万まで 10 万刻みで変化さグ結果の精度も評価するために，合成データを用いた実験結果せた．図 8(a) は横軸を頂点数とし，縦軸に処理時間を示していを示す．合成データとしては，グラフクラスタリングアルゴリ

る．また，図 8(b) から図 8(d) は各評価指標の値を示している．

ズムの評価において良く用いられている，LFR (Lancichinetti–

GPU-LP は，Louvain と比較して平均 24 倍の性能，CPU-

Fortunato–Radicchi) ベンチマーク [15] によって生成したもの

LP と比較して平均 10 倍の性能となった（図 8(a)）．実データを用いた．LFR ベンチマークは，正解クラスタ情報を持ったグを用いた場合よりも性能向上率が低い原因は GPU-LP 収束のラフを生成するため，この情報を用いて精度の評価を行う．精遅さである．このデータでは，CPU-LP は 3 回または 4 回の度の評価には以下の三つの指標を用いた [5]：(1) NMI (normal-

イテレーションでラベルの更新を終了するが，GPU-LP は上ized mutial information)，(2) F 値，(3) ARI (adjusted Rand

限の 10 回まで行なっている．これは，ラベルの更新タイミン

index)．

グの違いに起因すると考えられる：CPU-LP では，ある頂点の

NMI は情報理論的考え方に基づき，二つのクラスタリング

更新先ラベルを計算したら，すぐにその頂点のラベルを更新す

結果 C, C ′ 間の類似度を測る指標である：$−2 i,j Nij log(Nij Nt /Ni· N·j )$NMI(C, C ′ ) = $i Ni· log(Ni· /Nt ) +j N·j log(N·j /Nt )る．一方，GPU-LP では，現在のラベルの状態をもとに，全てここで，N は混同行列であり，Nij は正解クラスタ Ci とアル案された半同期的方式 [6, 9] を採用することで解消できると考ゴリズムによって検出されたクラスタ Cj′ に共通して現れる頂えられる．

の頂点のラベルを一斉に計算し，同時に更新するという方式となっている．この方式は，CPU-LP が取っている方式と比較して，収束が遅いことが知られている [16]．この欠点は，近年提点の数である．また，Ni· と N·j はそれぞれ行列 N の i 行の和精度に関しては，Louvain と比較して，ラベル伝搬法に基づと j 列の和であり，Nt は行列 N 全体の和を示している．NMI

く手法（CPU-LP ，GPU-LP ）が非常に良い性能を示している

（図 8(b)–(d)）．CPU-LP と GPU-LP は同等の性能を示してお一方，ラベル伝搬法は処理の性質上，並列化に適した手法でり，また頂点数によらず一定の値を取っている．一方，Louvain

あり，その並列化も検討されてきている．Cordasco ら [6] は，は頂点数を増やすにつれて，精度が劣化していくという現象がラベル伝搬法におけるラベル更新処理を半同期的に行い，並列見られる．これは，Louvain が利用している指標であるモジュ処理する手法を提案した．ラベル更新処理において，各頂点のラリティが持つ欠点の一つの，resolution limit [10]（5. 節）にラベルを逐次的に更新する方式を非同期的，全ての頂点のラベ起因するものと考えられる．以上の結果から，ラベル伝搬法が，ルを同時に更新する方式を同期的と呼ぶ [21]．同期的更新は，処理性能と精度の両方において優れた性能を達成可能なことが並列化が容易であるが，収束に時間がかかるという問題が指摘確認できた．

されている [16]．そこで，Cordasco らは，グラフカラーリン

5. 関 連 研 究

グを用いて入力グラフをいくつかの集合に分割し，各集合毎に同期的処理を行うことで，収束を早めると同時に効率的に並列グラフクラスタリングに関する様々なアルゴリズムや評価指化可能な手法を考案した．しかし，Cordasco らは，実際に並標，並列化手法などが提案されている [25]．本節では，グラフ

列処理環境における処理性能に関する評価は行なっていなかっクラスタリングとその並列処理に関する研究について述べる．

た．Duriakova ら [9] は，Cordasco らの半同期的処理方式を拡

グラフクラスタリング． グラフクラスタリング結果の評価

張・実装し，並列処理環境における処理性能の評価を行なった．

指標としてモジュラリティ[17] が広く用いられており，モジュ

また，GPU を用いてグラフクラスタリング処理を高速化しラリティの最適化に基づく手法も数多く提案されている．モ

た研究もいくつか存在する．Stovall ら [24] は，構造類似度にジュラリティは，グラフクラスタリング結果を，ランダムグラ

基づくグラフクラスタリング手法 SCAN [26] を，GPU によっフとの差を見ることで評価する．具体的には，クラスタ内に辺て処理する GPUSCAN を提案した．ラベル伝搬法の GPU にが多く存在し，クラスタ間の辺が少ない場合に高い値となる．

よる並列化に関しては，Soman ら [23] が提案を行っている．彼

モジュラリティ最適化手法として有名なものの一つに Louvain

らは，マルチコア CPU と GPU の両方に向けた並列ラベル伝

法 [4] がある．Louvain 法は，貪欲的にモジュラリティを最適搬法を開発した．しかし，GPU 上での実装に関する議論はあ化することで高速に近似解を求める．また，Shiokawa ら [22]

まりされておらず，負荷分散等もあまり考慮されていない．一は，Louvain 法を拡張し，逐次的にクラスタの集約を行なうこ方，本研究では，負荷分散を考慮し，かつ GPU 上での処理にとで処理を高速化する手法を提案した．一方，モジュラリティ

適した処理のみで並列化をすることで高速化を達成している．

の問題点として，resolution limit というものが指摘されている [10]．これは，グラフが巨大になるほど，クラスタ構造を正確に検出できなくなるというものである．

6. 結

論

本稿では，GPU を用いたグラフクラスタリングの高速化手

別種のグラフクラスタリングアルゴリズムとして，本研究で法を提案した．具体的には，ラベル伝搬法をベースに，GPU 上もベースとしている，ラベル伝搬に基づく手法がある [16, 21]．

での高並列な処理に適した並列化を行なうことで高速化を図っ

Raghavan ら [21] は，2. 2 節で説明した，最も基本的なラベル

た．特に，segmented sort や segmente reduce などの，次数の伝搬に基づくグラフクラスタリング手法を提案した．Leung

偏りが存在するデータにおいても効率的に並列処理可能な演算

ら [16] は，ラベル伝搬法に，hop attenuation と node prefer-を活用した．実世界のデータと合成データを用いた評価実験に

ence と呼ばれる二つのヒューリスティックを導入することで，より，GPU によるラベル伝搬法が，既存手法と比較して，最大性能向上を図った．Wang ら [25] は，10 種類のグラフクラス

で 117 倍の高速化を達成可能なことを示した．また，クラスタ

タリングアルゴリズムを共通のフレームワークを用いて実装すリング精度の観点においても，良い性能を示すことを確認した．

ることで，公平な比較を行った．結果として，ラベル伝搬法が，今後の課題としては，GPU によるラベル伝搬法のさらなる処理時間と精度の両方の観点において良い性能となった．その高速化が挙げられる．例えば，収束を早めるために半同期的方ため，本研究ではラベル伝搬法をベースに，GPU によるグラ式を導入することで高速化が可能であると考えられる．

フクラスタリングの高速化を行なった．
並列グラフクラスタリング． 巨大なグラフデータを高速に謝辞 本研究の一部は科学研究費補助金基盤研究 B（26280037）お よ び 科 学 研 究 費 補 助 金 特 別 研 究 員 奨 励 費処理するためには，効率的な並列処理が重要である．そのため， （15J02121）によるものである．
様々な並列グラフクラスタリング手法が提案されてきている．
特に，その著名さから，モジュラリティに基づくアルゴリズムの並列化に関する研究が多く見られる．例えば，Djidjev ら [8]は，モジュラリティ最適化問題が最小カット問題へと帰着できるという性質を利用し，既存の並列グラフ分割手法を活用することで処理の高速化を図った．また，Que ら [20] は，8192
ノードからなる大規模分散処理環境において，ハッシュ表を活用したグラフ表現を用いることで，Louvain 法の高速な並列化を行なった．

文

献[1] A.-L. Barabási and R. Albert, “Emergence of scaling in random networks,” Science, vol. 286, pp. 509–512, Oct. 1999.
[2] S. Baxter, “Modern GPU,” http://nvlabs.github.io/
moderngpu/, 2013.
[3] N. Bell and M. Garland, “Implementing sparse matrixvector multiplication on throughput-oriented processors,” inProceedings of the Conference on High Performance Computing Networking, Storage and Analysis (SC), pp. 18:1–
18:11, Nov. 2009.

[4] V. D. Blondel, J.-L. Guillaume, R. Lambiotte, and E. Lefebvre, “Fast unfolding of communities in large networks,”
Journal of Statistical Mechanics: Theory and Experiment,
vol. 2008, no. 10, pp. P10008:1–P10008:12, Oct. 2008.
[5] M. Chen, K. Kuzmin, and B. K. Szymanski, “Community detection via maximization of modularity and its variants,” IEEE Transactions on Computational Social Systems, vol. 1, no. 1, pp. 46–65, Mar. 2014.
[6] G. Cordasco and L. Gargano, “Label propagation algorithm: A semi-synchronous approach,” International Journal of Social Network Mining, vol. 1, no. 1, pp. 3–26, 2012.
[7] L. Danon, A. Daz-Guilera, J. Duch, and A. Arenas, “Comparing community structure identification,” Journal of Statistical Mechanics: Theory and Experiment, vol. 2005,
no. 09, pp. P09008:1–P09008:10, Sep. 2005.
[8] H. N. Djidjev and M. Onus, “Scalable and accurate
graph clustering and community structure detection,” IEEE
Transactions on Parallel and Distributed Systems, vol. 24,
no. 5, pp. 1022–1029, May 2013.
[9] E. Duriakova, N. Hurley, D. Ajwani, and A. Sala, “Analysis of the semi-synchronous approach to large-scale parallel community finding,” in Proceedings of the Second ACM
Conference on Online Social Networks (COSN), pp. 51–62,
Oct. 2014.
[10] S. Fortunato and M. Barthlemy, “Resolution limit in community detection,” Proceedings of the National Academy of
Sciences of the United States of America, vol. 104, no. 1,
pp. 36–41, Dec. 2007.
[11] S. Fortunato, “Community detection in graphs,” Physics
Reports, vol. 486, no. 3–5, pp. 75–174, Feb. 2010.
[12] M. Harris, “Parallel prefix sum (scan) with CUDA,” http://developer.download.nvidia.com/compute/cuda/2_2/sdk/
website/projects/scan/doc/scan.pdf, 2009.
[13] B. He, M. Lu, K. Yang, R. Fang, N. K. Govindaraju, Q. Luo,and P. V. Sander, “Relational query coprocessing on graphics processors,” ACM Transactions on Database Systems,
vol. 34, no. 4, pp. 21:1–21:39, Dec. 2009.
[14] J. Leskovec and A. Krevl, “SNAP Datasets: Stanford large
network dataset collection,” http://snap.stanford.edu/
data, Jun. 2014.
[15] A. Lancichinetti, S. Fortunato, and F. Radicchi, “Benchmark graphs for testing community detection algorithms,”
Physical Review E, vol. 78, no. 4, pp. 046110:1–046110:5,
Oct. 2008.
[16] I. X. Y. Leung, P. Hui, P. Liò, and J. Crowcroft, “Towards
real-time community detection in large networks,” PhysicalReview E, vol. 79, no. 6, pp. 066107:1–066107:10, Jun. 2009.
[17] M. E. J. Newman and M. Girvan, “Finding and evaluating community structure in networks,” Physical Review E,
vol. 69, no. 2, pp. 026113:1–026113:15, Feb. 2004.
[18] J. D. Owens, M. Houston, D. Luebke, S. Green, J. E. Stone,
and J. C. Phillips, “GPU computing,” Proceedings of theIEEE, vol. 96, no. 5, pp. 879–899, May 2008.
[19] S. Papadopoulos, Y. Kompatsiaris, A. Vakali, and P. Spyridonos, “Community detection in social media,” Data Mining and Knowledge Discovery, vol. 24, no. 3, pp. 515–554,
May 2012.
[20] X. Que, F. Checconi, F. Petrini, and J. Gunnels, “Scalable
community detection with the Louvain algorithm,” in Proceedings of the 2015 IEEE International Parallel & Distributed Processing Symposium (IPDPS), pp. 28–37, May
2015.
[21] U. N. Raghavan, R. Albert, and S. Kumara, “Near linear
time algorithm to detect community structures in largescale networks,” Physical Review E, vol. 76, no. 3, p. 036106,
Sep. 2007.
[22] H. Shiokawa, Y. Fujiwara, and M. Onizuka, “Fast algorithm
for modularity-based graph clustering,” in Proceedings of

[23]

[24]

[25]

[26]

the Twenty-Seventh AAAI Conference on Artificial Intelligence (AAAI), pp. 1170–1176, Jul. 2013.
J. Soman and A. Narang, “Fast community detection algorithm with GPUs and multicore architectures,” in Proceedings of the 2011 IEEE International Parallel & Distributed
Processing Symposium (IPDPS), pp. 568–579, May 2011.
T. R. Stovall, S. Kockara, and R. Avci, “GPUSCAN: GPUbased parallel structural clustering algorithm for networks,”
IEEE Transactions on Parallel and Distributed Systems,
vol. 26, no. 12, pp. 3381–3393, Dec. 2015.
M. Wang, C. Wang, J. X. Yu, and J. Zhang, “Community
detection in social networks: An in-depth benchmarking
study with a procedure-oriented framework,” Proceedings of
the VLDB Endowment, vol. 8, no. 10, pp. 998–1009, Jun.
2015.
X. Xu, N. Yuruk, Z. Feng, and T. A. J. Schweiger, “SCAN:
A structural clustering algorithm for networks,” in Proceedings of the 13th ACM SIGKDD International Conferenceon Knowledge Discovery and Data Mining (KDD), pp. 824–
833, Aug. 2007.