DEIM Forum 2014 D5-4

GPU を用いた Canopy クラスタリングの高速化
林

史尊†

小澤

佑介†

天笠

俊之††

北川

博之††

† 筑波大学大学院システム情報工学研究科 〒 305–8573 茨城県つくば市天王台一丁目 1 番 1 号†† 筑波大学システム情報系 〒 305–8573 茨城県つくば市天王台一丁目 1 番 1 号
E-mail: †{hayashi,kyusuke}@kde.cs.tsukuba.ac.jp, ††{amagasa,kitagawa}@cs.tsukuba.ac.jp
あらまし

大規模データに対するクラスタリングは計算コストが高く，高速化が求められている．高速化を実現する手法の一つに Canopy クラスタリングがある．Canopy クラスタリングは，凝集法や k-means 法といった一般的なクラスタリングアルゴリズムを高速に実行するための事前処理として行われる．また，クラスタリングアルゴリズムの並列計算による高速化もさかんに研究されており，その中で並列計算に優れた GPU の利用が注目されている．一般的なアルゴリズムは GPU を用いた高速化が研究されているが，Canopy クラスタリングを GPU 上で行う試みは現在までなされていない．本研究では，GPU を用いて Canopy クラスタリングを実行し，高速化の程度を検証した．
キーワード

データクラスタリング, Canopy クラスタリング, GPU

1. は じ め に一方，計算機処理を並列化する手段の一つとして，GPU

近年，科学分野で扱われるデータは非常に膨大なものとなっ（Graphics Processing Unit）コンピューティングがさかんにている．それに合わせて，大規模データの解析に要するコス

なっている [2]．GPU はグラフィックス処理に特化した演算装

トは肥大化しており，効率的にデータを処理する手法が求めら置である．小規模なプロセッサを大量に持つという構造の性質れている．データ・クラスタリング（以降，単に「クラスタリ上，単純な計算の並列処理に優れる．GPU コンピューティン

ング」と呼ぶ）はデータ解析によく用いられる手法で，データ

グは，GPU を利用して，グラフィックス処理以外の計算を行う

集合をデータ間の類似度に基づきグループ化する処理である．

ものである．ただし，GPU の構造は特殊であるため，利用のクラスタリング手法は多く提案されており，データの特徴やグ際にはアルゴリズムを最適化する必要がある．

ループ化の目的によって用いる手法を選ぶ．しかし，クラスタ

GPU コンピューティングにおける凝集法や k-means 法など

リングは一般に計算コストが高いため，大規模データに適用すについてはすでに研究されているが，Canopy クラスタリング

るのは難しく，計算の高速化が求められている．

を GPU 上で行う手法は現在まで提案されていない．本研究でクラスタリングの計算高速化方法の一つとして，Canopy ク

は，GPU を利用して，Canopy クラスタリング並列的に処理

ラスタリング [1] が注目を集めている．Canopy クラスタリン

する．GPU 上では，データ間の距離計算などを並列に行って，グは，凝集法や k-means 法といった，よく知られるクラスタリ

順次 Canopy を生成する．さらに，ナイーブな Canopy クラス

ング手法の前準備として行われる．Canopy と呼ばれるおおま

タリングを GPU 上で実行するだけでなく，より高効率に並列かなクラスタを，簡単な計算によって事前に生成しておくこと処理を行えるよう，セル構造を用いた独自の手法を提案する．

で，後に用いるクラスタリングアルゴリズムでの計算を軽減で独自手法では，データ空間をいくつかのセルに分割し，それをきる．Canopy の生成は，Canopy の中心点を選んで，その中利用して Canopy 生成時の計算範囲を削減することを狙いとす

心から一定の距離内にあるデータを Canopy として登録し，まる．各データが所属するセルとそれに隣接するセルを求めるこた異なる一定の距離内にあるデータを中心点の候補から取り除とで，各データにおける距離計算の対象データは所属・隣接セ

くことを繰り返して実現される．生成された Canopy によってル内のもののみに限定される．これにより，特に大規模なデー

距離の近い（類似度の大きい）データを求めているため，そのタセットに対するクラスタリング計算を高速化できる．独自手

後厳密なクラスタリング計算を行う段階では，すべてのデータ法についても GPU 上での処理を行った．セル構造を用いるこ間の距離を求める必要がなくなる．

とで 1 個のデータに対応する計算範囲が制限されたことにより，また，異なるクラスタリングの高速化方法として並列化が挙並列して複数の Canopy を生成することが可能となる．

げられる．これは，データ集合を複数に分割し，各部分集合でまた，本研究では，ナイーブな Canopy クラスタリング，セル処理を同時実行した後，結果を結合させるというアイディアに構造を用いた Canopy クラスタリングそれぞれについて，CPU基づいている．例えば，クラスタリングの各データ間の距離をおよび GPU 上で実行した際の計算時間を比較する実験を行っ計算する処理は，並列的に行うことで計算時間を削減できる．

た．実験により，2 次元 1,000,000 データを対象にした場合，

ただし，データ分割の方法によって並列処理の効率は大きく変GPU 上の独自手法が CPU 上のナイーブな手法に比べておよ

わるため，適切な方法でデータ分割を行うことが重要となる．

そ 15 倍高速化しているという結果が得られた．

モリは，CPU 側からアクセス可能で，かつどの SM からもア

2. 基本的事項

クセスでき数百 MB〜数 GB と大容量だが，アクセス遅延が大2. 1 クラスタリング

きい（処理実行の単位時間に比べ数十〜数百倍）．GPU での並クラスタリングとは，データ集合をデータ間の類似性に基づ

列処理を行う際は，各メモリの特徴を捉えてデータ分割を行う

きグループ（クラスタ）化する処理であり，規則性の発見やデー必要がある．

タの分類に利用される．具体的には，顔認識，マーケティング

GPU を用いた処理の基本的な流れは次のようになる．

リサーチ，文書のトピック別分類などといった，非常に幅広い（ 1 ） CPU 上で，GPU を使った計算に必要なデータを準備分野で用いられている．クラスタリングは，その結果が階層構する（GPU 上のメモリ確保等）

造を持つか否かで，階層型と非階層型の二種類に分けられる．

（ 2 ） データを GPU に転送する

階層型クラスタリングでは，逐次的にクラスタの結合または（ 3 ） カーネル（GPU 上の処理を定義した関数）を呼び出す

分割を繰り返すことで，階層的クラスタ構造を得る．凝集法が（ 4 ） GPU による処理

代表的手法であり，デンドログラムと呼ばれる階層構造が得ら（ 5 ） 処理後のデータや結果を必要に応じて CPU に転送れる．凝集法は，各クラスタ間の類似度を計算し，類似度が最するも高い 2 個のクラスタを結合することを，クラスタ数が 1 個になお，基本的にカーネル内で動的なメモリ確保は行えない．

なるまで繰り返す（終了条件のクラスタ数を設定することも可

そのため，どれほどのサイズが必要か分からないデータが存在

3能）．凝集法の計算量は，データ数を n とすると O(n ) と大する場合であっても，カーネル実行前に適当なサイズ分のメモ

きい．

リを確保しなければならない．また，CPU–GPU 間のデータ

非階層型クラスタリングでは，クラスタの分割精度を評価する関数を用意し，評価値が最適となる分割を探索する．代表的のやり取りは時間がかかるため，データ転送およびカーネル呼び出しの回数はなるべく抑える必要がある．

な手法に k-means 法がある．k-means 法は，重心をいくつか

GPU コ ン ピュー ティン グ 向 け の 統 合 開 発 環 境 の 一 つ に ，定め，各重心について近いデータを求め，各重心を再計算する

NVIDIA が提供する CUDA（Compute Unified Device Ar-という処理を，すべての重心の位置が変化しなくなるまで繰り

chitecture）がある．CUDA は，GPU が行う処理の指示，主返す．k-means 法の計算量は，データ数を n ，求めるクラスタ

記憶と GPU 用メモリ間のデータ転送などをサポートする．プ

数を k ，反復回数を t とすると O(tkn) となる．

ログラミングレベルでは，基本的に C/C++言語のライブラリ

これらの手法を単純に処理するのは非効率的である．そのたとして利用する．CUDA プログラミングは，GPU アーキテクめ，アルゴリズムの改良や事前処理などが行われる．前者の例

チャと対応する階層構造（スレッディングモデル）を基軸としとしては，F. Murtagh が最近傍グラフを利用し凝集法の計算て行われる．スレッディングモデルは，スレッド，ブロック（ス

量を O(n2 ) に抑えている [3] ．後者の一例には Canopy クラス

レッドブロック），グリッド（カーネルグリッド）からなる．ス

タリング [1] がある．これは，あらかじめ簡単な計算によりおレッド処理は各 SP で実行され，ブロック処理は各 SM で実行

おまかなクラスタを生成しておくことで，後に行う一般的なクされる．ブロックの集合をグリッドという．グリッドは GPU

ラスタリング計算を軽減するという手法である（詳細なアルゴ上で行う一連の処理に当たる．通常，1 ブロックでは数百個のリズムは 3.1 節で述べる）．

スレッドを実行し，1 グリッド内では数百〜数千万のブロック

2. 2 GPU コンピューティングが実行される （注 1） ．実際に使用するスレッド数とブロック数

GPU はグラフィックス処理に特化した演算装置である．一は，処理内容に基づき適切な値をあらかじめ設定しておく必要

般的演算を行う CPU（Central Processing Unit）は逐次的ながある．

処理に優れるのに対し，GPU は簡単な処理の並列実行に優れ論理的なメモリ構造としては，例えば，ブロック処理中にアる．それまで逐次的に行っていた処理を GPU 上で並列的に行クセスできるシェアードメモリや，CPU 側からもアクセス可うことで，単に CPU のみを使う場合に比べ実行アルゴリズム

能なグローバルメモリなどがある．それぞれ，シェアードメモ

全体の計算時間を大きく短縮できると期待される．GPU をグ

リはオンチップメモリに，グローバルメモリはデバイスメモリ

ラフィックス処理とは異なる用途で利用することが注目されてに対応する．CPU から転送されたデータは，基本的にグローいる [2] ．このように，GPU を汎用的計算に活用する技術を

バルメモリに置かれる．また，先述の通り，グローバルメモリ

GPU コンピューティングと呼ぶ．

（デバイスメモリ）はアクセス遅延が大きいため，データの再利

CPU は大きなプロセッサを数個持つが，GPU は小さなプロ用などを行う際にはシェアードメモリを活用する．なお，GPU

セッサを数百個持つ．GPU 内には複数のストリーミングマル

はグローバルメモリやシェアードメモリ以外にも数種類のメモ

チプロセッサ（SM）が存在し，SM はまた複数のストリーミン

リを内蔵するが，ここでは説明を割愛する．

グプロセッサ（SP）から構成される．SM は一つの処理を行う単位となる．

GPU が利用できるメモリもまた，複数の層から構成される．
例えば，各 SM が持つオンチップメモリは，対応する SM から（注 1）：本研究で用いる Fermi アーキテクチャ [5] では，1 ブロック内最大ス

高速アクセス可能だが，数十 KB と容量が小さい．デバイスメ

レッド数は 1024 ，1 グリッド内最大ブロック数は 65535 × 65535 である．

3. ナイーブな Canopy クラスタリング
3. 1 アルゴリズム
Canopy クラスタリングは，McCallum, A. らによって提案された手法であり [1]，凝集法や k-means 法といったクラスタリン
グ手法の前処理として利用される．これまで大きな計算時間を要していた大規模データクラスタリングを，実用的なレベルまで高速化できるため，近年注目されている．例えば，Zhi-Gang

Fan らは，画像データのクラスタリングにおいて，Canopy ク
ラスタリングを用いた高速化を図った [6] ．また，Tuli De らは， 700000 × 1500 という大きなサイズの銀河系外スペクトル
データに対し，Canopy クラスタリングを利用してクラスタリング行った [7] ．
この手法を用いたクラスタリングでは，まず，簡単な距離計算を行ってデータ集合を Canopy と呼ばれるグループに分割す

図 1 Canopy の生成．アルファベットのついた点は各 Canopy の中る（1 個のデータが複数の Canopy に属してもよい）．その後，

心．実線は中心点からの距離 T1 を示す．破線は中心点からの距凝集法や k-means 法のような厳密なクラスタリングを行うが，離 T2 を示す．この例では，まず初めに，すべての点からランダ

ここでの各データ間の距離計算はそれぞれの Canopy の中のみ

ムに 1 個を選択し，選ばれた点（A）を 1 個目の Canopy の中で行われる．クラスタとなり得るデータをあらかじめ Canopy

心としている．その後，距離 T1 内の点を Canopy に含め，距離

としてまとめておくことで，実際のクラスタリング処理におけ

T2 内の点を中心点候補から除外する．これで 1 個の Canopy が生成される．2 個目の Canopy 中心は，残っている中心点候補かる無駄な距離計算を省ける．

ら新たに選出される．選ばれた点（B）を中心に，A の時と同じ

具体的な例を挙げると，基本的な凝集法では，全てのクラス

ように Canopy を生成する．同様の手順で， C，D，E をそれ

タ組合せについて距離計算を行い，最も距離が小さいものを新ぞれ中心とする Canopy が生成される．E の Canopy が生成さたなクラスタとして結合するという処理を繰り返す．特に初めれた時点で，すべての点が中心点候補から除外された，すなわちは，各々のデータが 1 個のクラスタとして扱われるため，距離各点がいずれか 1 個以上の T2 円内に含まれたため，Canopy ク

計算が必要な組合せ数は非常に大きい．10000 個のデータにつラスタリングの処理を完了する．
（図は”Eﬃcient Clustering of

いて凝集法クラスタリングを行う場合であれば，1 個のデータHigh Dimensional Data Sets with Application to Reference
Matching” [1] Figure 1 を引用）につき，他 9999 個のデータとの距離を求めなければならない．
ここで，Canopy クラスタリングをあらかじめ行った場合を考える．あるデータが所属する Canopy が，仮に 1000 個のデー•

複数の Canopy に，同じデータが所属し得る．

タを持っているとする．この状態で凝集法を開始したとき，そ•すべての Canopy の中心点について，互いの距離は T2のデータと距離計算を行うデータは，同じ Canopy に存在すより大きい．

る他 999 個のデータのみでよい．Canopy 外のデータは，元々なお， T1 および T2 の値は，データの性質や距離尺度等をも同じクラスタにはなり得ない（なりにくい）ものとして考えるとに設定する．T1 は各 Canopy が含むデータ数に関わり， T2ためである．このように，最終的に同じクラスタを形成しないは Canopy 生成数に関わる．生成した Canopy の大きさや数データとの距離計算を行わずに済むため，計算回数を抑え処理

は，この後に行われる本格的なクラスタリングの計算結果の妥

時間を削減できる．

当性，および計算時間に影響する（Canopy クラスタリング後

Canopy を生成するアルゴリズムの流れは以下のようになる（図 1 を参照）．

の処理については，本研究では扱わない）．

4. GPU 上での並列処理（ 1 ） クラスタリングを行いたいデータセットを用意し，各データのインデックスを中心点候補とする．

GPU を用いてナイーブな Canopy クラスタリングを実現す（ 2 ） 中心点候補の中から 1 個のデータをランダムに選択る場合，データセットの準備および終了条件判定を除く以下のする．

処理を GPU 上で行うことができる．

（ 3 ） 2 で選んだデータと他のすべてのデータについて，各々（ 2 ） 中心点候補の中から 1 個のデータをランダムに選択の距離を計算する．

する．

（ 4 ） 距離が T1 以下のデータを Canopy に含める．

（ 3 ） 2 で選んだデータと他のすべてのデータについて，各々（ 5 ） 距離が T2 (< T1 ) 以下のデータのインデックスを中心の距離を計算する．

点候補から除外する．

（ 4 ） 距離が T1 以下のデータを Canopy に含める．

（ 6 ） 中心点候補が 0 個になるまで，2〜5 の処理を繰り返す（ 5 ） 距離が T2 (< T1 ) 以下のデータのインデックスを中心生成された Canopy は，以下の性質を持つ．

点候補から除外する．

ここでは主に，3, 4, 5 について説明する．なお，Canopy ク
ラスタリングは他の一般的なクラスタリング手法の事前処理として行われるため，後の処理を GPU 上で行うことを前提として実装を行っている. 具体的には，各 Canopy にどのデータが所属しているかという情報を GPU 上に残すようにした．

4. 1 中心点の選択と距離計算
処理の簡単のため，新しく選択する Canopy の中心は，常に図 3 Canopy の更新．T1 フラグ列（1 段目）と Inclusive Scan 後の中心点候補の先頭要素とする．アルゴリズムの性質上，どのT1 フラグ列（2 段目）を用意する．T1 フラグ列中で値が 1 でデータが中心として選択されたとしても，T2 近傍にあるデー

あるインデックスについて，Inclusive Scan 後の T1 フラグ列

タは必ず中心点候補から削除される．そのため，中心点選択をランダムに行わないことの影響は小さいと考えられる．
続いて，中心点となったデータと各データとの距離計算を，の値を確認する．インデックス値を，Canopy メンバ列中の確認した値の位置に格納する．なお，2 回目の Canopy 生成処理以
降は，既存の Canopy メンバ数 + Inclusive Scan 後の T1 フラ
グ列の値の位置にインデックス値を格納すればよい．

各スレッド上で並列に行う．1 個のスレッドで中心点データと比較対象のデータ 1 個を読込み，その距離を算出する．さらに，同一スレッド上では，算出距離と T1 および T2 値との比較を行い，T1 より距離が小さい場合は T1 フラグとして 1 を，T2
より距離が小さい場合は T2 フラグとして 1 を，それぞれ T1
フラグ列， T2 フラグ列に記録する（1 でない場合は 0 を記録する）．
距離計算および T1 , T2 フラグを求める処理について，簡単な図例を用いて詳しく説明する．図 2 は，8 個のデータについて処理を行おうとするものである．ここで， T1 = 7, T2 = 4と設定し，中心点には data3 が選択されたものとする．各ス
レッドは，まず，スレッド番号に対応するインデックスのデー
タ（data0〜data7）と，中心点となるデータ（center，この例
では data3）を読込む．その後，距離計算を行ったところ，それぞれのスレッドで 3 段目のような結果が得られた．各スレッ
ドは自身の中で求められた距離と T1 , T2 の値を比較して，割当てられたデータが中心点に対しどのような位置にあるかを，

T1 フラグ， T2 フラグとして求める．例えば data0 ならば，中心点から距離が 2 であり，これは T1 および T2 以下であるため，それぞれのフラグに 1 が記録される．

した中心点に対応する Canopy の所属データ（Canopy メンバ）を記録する．このとき，Canopy メンバを更新する場合に必要な情報は，計算対象とする Canopy にどのインデックスのデータが所属しているかである．これは， T1 フラグ列中で 1 の値を持
つインデックスを集めることで得られる．フラグの立ったイン
デックスを抽出するには，Inclusive Scan [9] と呼ばれる処理を行う．Inclusive Scan では，配列 Qin = [q0 , q1 , q2 , ..., qn ] から
配列 Qout = [q0 , q0 + q1 , q0 + q1 + q2 , ..., q0 + q1 + q2 + ... + qn ]が得られる．

Inclusive Scan 後の T1 フラグ列は，T1 フラグが立っている
データの数を先頭から順に記録しているものとみなせる．この特徴を利用することで， T1 フラグ列と Inclusive Scan 後の T1
フラグ列から，T1 フラグの立ったインデックスを抽出した配列が得られる．まず， T1 フラグ列において値が 1 であるインデックスの，Inclusive Scan 後の値を得る．抽出配列では，この値の位置にインデックス値が格納される（図 3 参照．ただし，配列のインデックスは 0 から数えるため，実際に格納される位
置は，Inclusive Scan 後の値 −1 となる）．この抽出配列は計算対象である Canopy のメンバであるため，これを Canopy メ
ンバ列に追加する．

4. 3 中心点候補の更新
続いて， T2 フラグ列を用いて，中心点候補の除去を行う．
その時点で残っている中心点候補それぞれの T2 フラグ列を調べ，値が 1 であるものを除去，あるいは値が 0 であるものを抽出すればよい．これは，以下に示す 3 段階の処理によって実現される．
（ 1 ） 中心点候補に残っているデータのインデックスについてのみ， T2 フラグ列から値を抽出する．また，3 の処理のた
図 2

距離計算の並列化．各スレッドは，スレッド番号に対応するイ
ンデックスのデータと，中心点となるデータを読込む（この例

め，抽出時にフラグの値（ 0, 1 ）を反転させる．
（ 2 ） 抽出 T2 フラグ列の Inclusive Scan を計算する．

では，data3 が中心点として選ばれている）．中心点データと各（ 3 ） その時点での中心点候補，抽出 T2 フラグ列，Inclusive

データで距離計算を行い，その結果によって，T1 , T2 フラグのScan 後の抽出 T2 フラグ列から，中心点候補として残るデータ値を決定する．

のインデックスを求める．

1 の処理は単純に，中心点候補であるデータのインデックス
4. 2 Canopy の更新

それぞれについて T2 フラグ列の値を確認し，抽出 T2 フラグ

距離計算処理によって得られた T1 フラグ列を用いて，選択列に格納する（図 4 参照）．このとき，3 の処理で T2 フラグ図4

中心点候補の更新 (1)．更新前の中心点候補（1 段目）それぞれの T2 フラグの値を，抽出 T2 フラグ列として得る．更新後はT2 フラグの値が 0 であるもの（中心点から距離 T2 以内に存在しないデータ）を残したいため，後の処理に向けて T2 フラグの値を反転させておく．

図 6 セル構造を用いた複数 Canopy の並列生成．この例は，2 次元
データ空間にセルを生成したものである．中心点（緑）を，計算
範囲（赤および橙のセル）が重複しないように選択する．これにより，複数の Canopy を並列して生成することが可能となる．

る．データ空間を間隔 T1 のグリッド状に分割し，分割された
個々の領域をセルと呼ぶ．このとき，すべてのデータは必ずい
図 5 中心点候補の更新 (2)．抽出 T2 フラグ列と Inclusive Scan 後の抽出 T2 フラグ列から，中心点候補が更新後，どの位置に格納されるかが分かる．ここから，T2 フラグが 1 になっている中心点候補を抜き出せる．

ずれか 1 個のセルに含まれる．ここで，ある 1 個のデータがCanopy の中心点に選ばれたとする．この中心点に対する計算
範囲は，中心点が所属するセルと，それに隣接するセルになる．
セルのサイズから，その Canopy に含まれる（中心点から距離が 0 であるものを抽出するため，T2 フラグの値を反転させる．

3 の処理の基本的な流れは，Canopy 更新時の処理と同様である．抽出 T2 フラグ列中で値が 1 であるインデックスについて，Inclusive Scan 後の抽出 T2 フラグ列の値を確認する．中心点候補の同インデックスに格納されている値を，新しい配列の確認した値の位置に格納し，これを更新後の中心点候補とする（図 4 参照）．なお，抽出 T2 フラグ列の末尾の値は，更新後の中心点候補の総数を意味する．この情報は CPU 側に転送され，Canopy クラスタリングの終了条件判定に用いられる（値が 0 であれば Canopy クラスタリングを終了する）．

5. セル構造を用いた Canopy クラスタリング
5. 1 基本アイディアとアルゴリズム
本研究で提案する独自の手法は，複数 Canopy の並列生成を狙いとする．4 節では，ナイーブな Canopy クラスタリングのアルゴリズムと，GPU による並列化の方法について述べた．この方法は，1 個の Canopy 中心点について，データセット中の全データとの距離計算を行い，Canopy メンバと中心点候補の更新を行うという処理を繰り返すものだった．しかしながら，1個の Canopy に含まれるデータは，中心点から距離 T1 以内に存在するもののみである．中心点から明らかに離れているデー
タについても計算を行うのは非効率的である．中心点近傍に存在するデータだけを計算対象にすると，処理に必要な計算回数が減少する．こうして Canopy 1 個あたりの計算範囲を制限することで，複数の Canopy を並列に生成できると考えられる．
計算範囲を決定するために，独自手法ではセル構造を用いT1 以内にある）データは，中心点の所属セルおよび隣接セル
内にしか存在し得ないためである．同様に，中心点と計算範囲が重複しないように，異なる中心点を可能な限り選択する（図6 参照）．選択後，すべての計算範囲内のデータについて，距
離計算，Canopy および中心点候補の更新を行う．
独自手法のアルゴリズムは以下のようになる．
（ 1 ） クラスタリングを行いたいデータセットを用意し，各データのインデックスを中心点候補とする．
（ 2 ） 各データの所属セルを計算し，それぞれの所属セルに対応する隣接セルを求める．
（各所属セルはセル番号を与えて管理する）．具体的には，全データについて以下の処理を行う．

(a) 各次元の値と T1 との商（小数点以下切り捨て）を求め，所属セルの位置情報とする．

(b) 同じ位置情報を持つセルがそれまでに現われていな
ければ，その所属セルに新しいセル番号（直前の番号 +1）を与える．

(c) 新しいセル番号が現われたとき，これまでに現れたセルの位置情報と比較して，すべての次元について差が 1 以下であるセルを，その所属セルの隣接セルとする．
（ 3 ） 中心点候補の中から，計算範囲が重複しないように複数の中心点を選択する．具体的には，以下の処理を行う．

(a) この時点での中心点候補と同一の配列を用意し，これを cp 中心点候補とする．

(b) 選択した中心点を格納する配列を用意し，これを中心点列とする．

(c) cp 中心点候補の先頭を選択して取り出し，中心点列に格納する．

(d) 選択した中心点の所属セル，またその隣接セルに所属するデータを，cp 中心点候補から除去する．

(e) cp 中心点候補が 0 個になるまで，c，d の処理を繰り返す．
（ 4 ） 各中心点に対応する計算範囲内のデータを抽出する．
具体的には，以下の処理を行う．

(a) 各中心点に対応する計算範囲内のデータのインデックスを格納する配列を用意し，これを各中心点対応データメ
ンバ列とする．

(b) 各中心点について，その所属セルおよび隣接セルに

図 7 複数中心点の選択．1 段目は cp 中心点候補，2 段目は各候補の所属セル番号，3 段目は所属セルの隣接セルのセル番号である．

cp 中心点候補には所属・隣接セルに所属しないデータを残すため，これらに残存フラグを立てる．その後は，3.2.3 節の処理と同様に進め，cp 中心点候補を更新する．

所属するデータを，各中心点対応データメンバ列に加える．
（ 5 ） 3 で選んだ複数の中心点と，それぞれの中心点に対応する計算範囲内のデータについて，各々の距離を計算する．
（ 6 ） 距離が T1 以下のデータを Canopy メンバに追加する．
（ 7 ） 距離が T2 (< T1 ) 以下のデータのインデックスを中心点候補から除外する．
（ 8 ） 中心点候補が 0 個になるまで，3〜7 の処理を繰り返す
ナイーブな Canopy クラスタリングと比較して，2〜4 の処理が追加されている．4 の処理後の基本的な流れはナイーブな手法と変わらない．

2 については，処理が逐次的で GPU 上での並列処理に不向きであるため，CPU 上で行うものとする．ここでは，全デー
タについて，各次元の値からそのデータがどのセルに所属しているかを求める．実際には，それぞれのデータに対し，各次元の値とセルのサイズの商（小数点以下切り捨て）をそのデータが所属するセルの座標とする．セルにはセル番号が与えられ，既存しない座標のセルが現われた場合には，そのセルに新しい番号を与える．データとセルは，データ番号およびセル番号を用いて互いの所属情報を持つ．ここで，セルのサイズを一辺T1 とすることにより，あるデータが Canopy 中心点として選択されたとき，その Canopy のメンバとなり得るデータは，その中心点の所属セルとその隣接セルのみであると判定できる．

Canopy 計算時には，各中心点に対して所属・隣接セル内のデー
タを求めることで，距離計算の回数を削減することができる．
あるセル Cn = (cn,1 , cn,2 , ..., cn,d ) において，その隣接するセルの座標は， Cn との各次元の差が ±1 であるという条件を満たす．これを利用することで，隣接セルを以下のように求められる．
（ 1 ） 既存セルを，小さい次元を優先して昇順ソートする（辞書順ソート）．
（ 2 ） k の初期値を 0 として，以下の処理を再帰的に繰り返す．
（ a ） cn,k と各既存セルの k 次元目の座標を比較する．
（ b ） 初めて cn,k − 1 < ci−1 ,k となる位置 i−1 ，初めてcn,k > cj−1 ,k となる位置 j−1 を求める．
（ c ） 初めて cn,k − 1 < ci+1 ,k となる位置 i+1 ，初めてcn,k > cj+1 ,k となる位置 j+1 を求める．
（ d ） [i−1 , j−1 ], (j−1 , i+1 ), [i+1 , j+1 ] の範囲内それぞれで，

k + 1 次元目の座標を比較する．k = d のとき，範囲内にある

cn,d ± 1 と等しい座標を持つセルを，Cn の隣接データとする．
5. 2 GPU 上での並列処理
ここでは，5.1 節で示したアルゴリズム中の，主として 3 と4 の処理について詳しく説明する．なお，2 の処理は CPU 上で行うものとする．これは，2 の処理が逐次的であり，並列化に適していないためである．

5. 2. 1 複数中心点の選択
中心点を，その計算範囲（所属セル＋隣接セル）が重複しな
いように複数選択する．そのためには，新しく中心点を選択す
る度に，所属・隣接セル内のデータを候補から除去していけばよい．GPU では，以下のような処理を行う．

(a) この時点での中心点候補と同一の配列を用意し，これを cp(copied) 中心点候補とする．

(b) 選択した中心点を格納する配列を用意し，これを中心点列とする．

(c) cp 中心点候補の先頭を選択して取り出し，中心点列に格納する．

(d) 選択した中心点の所属セル，またその隣接セルに所属しないデータについて残存フラグを立てる（図 7 参照）．

(e) 残存フラグ列の Inclusive Scan を計算する．
(f) cp 中心点候補，残存フラグ列，Inclusive Scan 後の残存フラグ列を利用して，4.2 節で説明した処理と同様に cp 中心点候補を更新する．

(g) cp 中心点候補が 0 個になるまで，(c)〜(f) の処理を繰り返す．

5. 2. 2 各中心点に対応する計算範囲内のデータ列抽出
この処理は，複数中心点の選択と同一のカーネルで行うことができる．各中心点に対応するデータメンバ列として抽出する
データは，中心点の所属セルとその隣接セル中にあるため，計算内容は複数中心点の選択とほぼ変わらない．
複数中心点の選択時，選択した中心点に対応する計算範囲内のデータは，候補から除去するデータと同一であるため，これを利用して各中心点に対応する計算範囲内のデータを抽出することができる（図 8 参照）．また，Inclusive Scan 後の残存フ
ラグ列の末尾の値は，加えた対応データメンバ列の総数と一致するため，これを対応データ境界列に格納しておく．

図 8 各中心点に対応する計算範囲内のデータ列抽出．この例では，注目する中心点の所属セルのセル番号を 0，その隣接セルのセル番
号を 1，2，3，14 とした場合を考える．これらのセルに含まれるデータは，複数中心点の選択時にその候補から除去するデー
タと同一であるため，これをそのとき注目している中心点に対応するデータメンバとして得る．その後，別の中心点に対して同じ処理を繰り返す場合は，新しく得られるデータ列を対応デー
タメンバ列の末尾に加えていく．
図9

実験に用いたデータセットの例．100000 × 100000 の大きさを持つ 2 次元空間に，10000 個の点を生成したもの．正規分布を

5. 2. 3 Canopy の生成と中心点候補の更新

利用し，クラスタとなるようなむらを作っている．

Canopy の生成の際は，1 個の中心点に対する距離計算を 1個のブロックで行い Canopy メンバを決定する．また，各ブ
ロックで中心点から T2 以下の距離であるデータをグローバル

6. 2 実 験 結 果

メモリ中の T2 フラグ列に記録しておく．すべてのブロックで以下の 4 手法について，データ数を固定し次元数を変化させCanopy メンバを求めた後で， T2 フラグ列を基に中心点候補

る場合と，次元数を固定しデータ数を変化させる場合の 2 通りの更新を行う．

について，計算時間を比較する実験を行った．

•

CPU 上でのナイーブな Canopy クラスタリング•

GPU 上でのナイーブな Canopy クラスタリング

ここでは，ナイーブな Canopy クラスタリング，提案手法で•CPU 上でのセル構造を用いた Canopy クラスタリング

あるセル構造を用いた Canopy クラスタリングについて，CPU•GPU 上でのセル構造を用いた Canopy クラスタリング

および GPU 上で行った計算速度の比較実験の結果と考察を

以下に，それぞれの場合における実験結果を示す．

示す．

6. 2. 1 実験 1：次元数の変化による計算時間の変化

6. 評 価 実 験

6. 1 実 験 環 境

データ数を 100000 個に固定し，次元数を 2〜6 に変化させて，今回の実験で用いるマシンの性能を以下に示す．

4 手法の計算時間を比較した．また， T1 = 100000, T2 = 70000•に固定した．結果のグラフを図 10 に示す．

CPU

R
– Intel ⃝Xeon(R)E5630 (2.53GHz * 4 core)– OS：Ubuntu 12.04 LTS (64bit)– メモリ：4.0GB•

GPU– NVIDIA GeForce GTX 460 (1350MHz * 336 core)– メモリ：767MB
実験に用いたデータセットは人工のものであり，2〜6 次元のベクトルの集合である．データ数は 100〜1000000 個まで変化させる．このデータセットでは， (データ数 × 10)2 の大きさのデータ空間中に，正規分布を利用して各データをある程度偏在するように配置している（一例を図 9 に示す）．

図 10

実験結果（次元数を変化）

今回，データ間の距離尺度としてユークリッド距離を採用した．T1 , T2 の値は実験ごとに固定とする．また， T2 = 0.7T1とする．ただし，これはあくまで本実験のために定めたものであり，実際の T1 , T2 値は，データの性質や求める結果などに応じて定められるべきである．
計算時間の測定については，Canopy クラスタリングの処理部分のみを対象とする．CPU–GPU 間のデータ転送時間もこれに含めるものとする．

いずれの手法においても，次元数に比例して計算時間が増加している．CPU 同士で比較すると，次元数が大きい場合，セル構造を用いた手法が計算高速化に貢献していることが分かる．次元数が大きくなるほどセルの数が増加し，各々の中心点に対応する計算範囲内のデータが減少するためであると考えられる．GPU 上の手法では，両者に大きな差は見られなかった．
これは，提案手法において，CPU 上で行ったセルを生成する

処理が計算時間に含まれているためである．セル生成処理に費所 CPE25149），JSPS 科研費（24240015），「エクサスケール

やされる時間は，例えばデータ数 100,000 のとき全体のおよそ計算技術開拓による先端学際計算科学教育研究拠点の充実」事

10〜25%ほどであった．また，本研究で実装した GPU 化アル

業によるものである．

ゴリズムは，完全には最適化されていないと考えられ，工夫によってさらなる高速化が見込める．CPU と GPU で比較すると，GPU 上での計算によって大きく高速化できていることが見て取れる．

6. 2. 2 実験 2：データ数の変化による計算時間の変化
次元数を 2 に固定し，データ数を 100〜1000000 に変化させて，4 手法の計算時間を比較した．また， T1 = 0.5 ∗ データ数

, T2 = 0.35 ∗ データ数 に固定した．結果のグラフを図 11 に示す．

図 11 実験結果（データ数を変化）

データ数が小さいときに GPU 上の手法の計算時間が大きくなっているのは，CPU–GPU 間のデータ転送のオーバヘッドが大きいためである．また，CPU 上の独自手法もデータ数が小さい場合に計算時間がナイーブなものより大きくなっている
が，これは事前のセル生成処理による計算回数削減の影響が小さいためである．次元数 2 ，データ数 1, 000, 000 のとき，ナイーブな Canopy クラスタリングと GPU 上の独自手法を比較すると，計算時間がおよそ 15 倍速くなっていることが分かる．

7. 結

論

本研究では，大規模データを対象とするクラスタリングを高速化するため，Canopy クラスタリングの GPU を利用した並列計算を試みた．また，セル構造を用いてさらに並列度を高める独自の手法を提案した．データ数，次元数を変化させ実験を行った結果，独自手法がナイーブな Canopy クラスタリングと比較し，次元数 2，データ数 1, 000, 000 であるとき，およそ15 倍高速であることが分かった．
今後解決すべき課題として，GPU のメモリ容量に関する問題がある．GPU のメモリ容量は CPU に比べ小さく，データ
サイズによってはすべてのデータを GPU に転送できない場合がある．この問題の対策として，複数の GPU を用いる研究 [8]を参考に，CPU 側でデータセットを複数に分割し，分割されたサブデータセットごとに GPU へ転送し処理するという方法が考えられる．
謝辞 本研究の一部は，平成 25 年度共同研究（富士通研究文

献[1] A. McCallum, K. Nigam, L. H. Ungar, Eﬃcient Clusteringof High Dimensional Data Sets with Application to Reference Matching, Proceedings of the sixth ACM SIGKDD
international conference on Knowledge discovery and data
mining, pp. 169–178 , 2000.
[2] J. D. Owens, M. Houston, D. Luebke, S. Green, J. E. Stone,
J. C. Phillips, GPU Computing, Proceedings of the IEEE,
volume 96, Number 5, pp. 879–899, 2008.
[3] F. Murtagh, A Survey of Recent Advances in Hierarchical
Clustering Algorithms, The Computer Journal, Vol.26, pp.
354–359, 1983.
[4] NVIDIA, http://www.nvidia.com.
[5] NVIDIA, NVIDIA Partnership Training, 2011.
[6] Zhi-Gang Fan, Yadong Wu, Bo Wu, Maximum NormalizedSpacing for Eﬃcient Visual Clustering, Proceedings of 19th
ACM International Conference on Information and Knowledge Management (CIKM 2010), pp. 409–417, 2010.
[7] Tuli De, Didier Fraix Burnet, Asis Kumar Chattopadhyay,
Clustering large number of extragalactic spectra of galaxies
and quasars through canopies, Communication in Statistics- Theory and Methods (2013), 2013.
[8] Mohiuddin K. Wasif, P. J. Narayanan, Scalable clustering
using multiple GPUs, Proceedings of the 2011 18th International Conference on High Performance Computing, pp.1–
10, 2011.
[9] Yan Shengen, Long Guoping, Zhang Yunquan, StreamScan:
Fast Scan Algorithms for GPUs Without Global Barrier
Synchronization, Proceedings of the 18th ACM SIGPLAN
Symposium on Principles and Practice of Parallel Programming, pp. 229–238, 2013.