多次元点位置検索データ構造のGPUを用いた高速化とトラジェクトリパターン照合への応用

DEIM Forum 2016 G3-3

多次元点位置検索データ構造の GPU を用いた高速化と

トラジェクトリパターン照合への応用

山本

雅大

†

笹川

裕人

†

有村

博紀

†

†

北海道大学大学院情報科学研究科

〒 060-0814 北海道札幌市北区北 14 条西 9 丁目

E-mail:

†{

m yamamoto,sasakawa,arim

}

@ist.hokudai.ac.jp

あらまし 本論文では、自動車の GPS データ等のような軌跡データ (トラジェクトリデータ) に対するトラジェクト

リパターン照合法と、関連する点位置検索データ構造と呼ばれるデータ構造の GPU を用いた高速化について考察

する。はじめに、トラジェクトリパターン照合問題を考察し、山本らによって提案された、文字列照合アルゴリズ

ム Shift-And 法を拡張したトラジェクトリパターン照合アルゴリズム FastShift-AndTrajMatch を紹介する。このアル

ゴリズムでは、クエリ点に所属する面を答える点位置検索クエリのコストが非常に高く、実応用時には応用になると

考えられている。そこで本稿では、 FastShift-AndTrajMatch を、GPU を用いた並列化によって高速化することを考

える。基本的なアイディアは、クエリをブロックに分割することによる処理の並列化と、GPU の特性を考慮したアル

ゴリズムの改善である。

キーワード 並列データ構造, GPU, ストリーム検索, 時空間データ

1.

は じ め に

GPSやスマートフォンなどの位置センサーの普及により、自 動車や、人、野生動物などの移動データに代表される、トラジェ クトリデータの種類と量が増え続けている。また、これらの応 用として、これに関連する多次元点位置検索データ構造と呼ば れるデータ構造は、大規模データ解析や、時空間情報検索、コ ンピューターグラフィックスなどの分野で、様々な応用を持つ。 例えば、HONDAは、カーナビを搭載した車の移動データを元 に、現在通行が可能な道路の情報を示した、通行実績情報マッ プを公開している[6]。これらの技術において、位置センサーや ポイントクラウドなどのセンサー技術の発展により、入力点集 合データのサイズは年々増大しており、高速な位置検索データ 構造が求められている。 そこで本論文では、上記のようなトラジェクトリデータに対 する高速なパターン照合法と、点位置検索データ構造のGPU を用いた高速化について考察する。本論文で取り組む問題は以 下の通りである。 1. 1 トラジェクトリパターン照合問題 第1の問題として、3.節では、2次元空間のL∞距離に関す るトラジェクトリパターン照合問題について考察する。 はじめに、よく知られた文字列照合法であるShift-And法を 取り上げ、この手法に基づく素朴なトラジェクトリパターン照 合法を与える。文字列上のShift-And法は、パターンの前処理 によりマスクを構築して、高速化を実現していた。それに対し て、Shift-And法を用いた素朴なトラジェクトリ照合では、入 力が与えられないとパターンを前処理することができない。そ のため、この素朴なアルゴリズムでは、実行時に入力点が与え られるたびにマスクを生成し、照合を行う。ワード長wの計算 機上で、この素朴なアルゴリズムは、長さmの軌跡パターンと 長さnの軌跡テキストに対して、O(m)領域を用いて、O(mn) 実行時間で軌跡パターン照合を行う。 次に、山本らは、この手法を改良し、点位置検索データ構造 とビット並列計算を組み合わせて、より高速なトラジェクトリ パターン照合アルゴリズムFastShift-AndTrajMatchを提案し た[11]。このアルゴリズムは、点位置検索データ構造を用いる ことで、次のように動作する。 はじめに、前処理において、軌跡パターンP を受け取り、 点位置検索データ構造を構築する。実行時には、軌跡テキ ストT 上で、点位置検索データ構造を用いて、各入力点を 含む面の集合を求め、対応するNFAの状態遷移を行う。こ のアルゴリズムは、ワード長wの計算機上で、アルゴリズ ムFastShift-AndTrajMatchは、長さmの軌跡パターンと長さ nの軌跡テキストに対して、O(m log m + m⌈m/w⌉)前処理時 間とO(m⌈m/w⌉)領域を用いて、O(n(log m +⌈m/w⌉))実行 時間で軌跡パターン照合を行う。この計算時間は、素朴なトラ ジェクトリパターン照合法よりも高速である。 1. 2 多次元点位置検索データ構造のGPU上の並列実装 第2の問題として、4.節では、3.節で用いている点位置検索 データ構造が、巨大なサイズの軌跡パターンに対して、現実的 な時間で処理を終えることが出来ない問題について考察する。 2次元の点位置検索データ構造として、3.節では、高々m個 の正方形の配置による平面分割に対して、O(m2)語の領域と O(m2log m)前処理時間を用いて、1個の点位置検索クエリに O(log m)時間で答えるデータ構造を用いている。これはCPU の逐次計算を用いることを前提としているが、巨大なmに対 する効率の良い実現は難しい。そこで、この手法をGPUを用 いた並列計算によって高速化することを考える。

本章で用いるGPU(Graphic Processing Unit)は、数100∼ 数1000の多数のコアを持つ並列計算デバイスである。GPUは

元々画像処理用に開発されたが、現在では一般的な計算の高速化 にも広く用いられており、その並列化技術はGPGPU(General Purpose computing on GPU, GPU汎用計算)と呼ばれている。

提案手法は、与えられた入力点集合を定数cのサイズに分割 し、それぞれをマクロブロックとして、GPU上の独立したコ ア集団で並列に解くことを考える。各マクロブロック内の計算 については、次の4つのアプローチを提案した。 • データ構造A：マイクロ線形走査 • データ構造B：分岐レス二分探索 • データ構造C：ハイブリッドアプローチ • データ構造D：数値演算を用いたビット内並列比較 各アプローチに共通する考えは、(1)データを小さなブロッ クに分割し、個々のブロックをGPUのコアで並列に実行する ことと、(2)その際にGPUの特性を考慮して、分岐処理を減 らすことで高速化を行うことである。本論文では、上記の4つ のアプローチについて議論する。上記のアプローチの実装と計 算機実験による性能評価は、今後の課題である。 1. 3 本論文の構成 本論文の構成を述べる。はじめに、2章では、以降の章に必 要な定義や概念を導入した後、3.節では、2次元空間のL∞距 離に関するトラジェクトリパターン照合問題について考察する。 4.節では、多次元点位置検索データ構造をGPU上で並列実装 するための効率良い手法について考察する。最後に、5.で本論 文の結論と、今後の課題を述べる。 1. 4 関 連 研 究 GPUにおける並列化では、その特性から、プログラムに条件 分岐が少ない方がより並列化の効率が上がることが知られてい る。橋口ら[12]は、GPUを用いた論理シミュレーションの並 列高速化において、分岐演算数の削減の効果について議論して いる。また、坂本ら[13]は、超解像度画像処理におけるGPU 上での二分探索を議論しており、この部分の条件分岐が性能低 下を引き起こす可能性を指摘している。

2.

準

備

本節では、文献[2]にしたがって、計算幾何学の基本的な概 念と記法を導入する。その後、本稿で考察する多次元点位置検 索データ構造を導入する。以下に含まれない、アルゴリズム一 般の用語に関しては教科書[3]を参照されたい。文字列照合に 関しては教科書[9]を参照されたい。 2. 1 空間と距離 はじめに、正整数d >_{= 1}を次元と呼ぶ。集合Rd _{= { p =} (x1, . . . , xd)| xi∈ R, 1 <_{= i <= d }}をd次元平面と呼び、各要素 p = (x1, . . . , xd)∈ Rdを点という。ここで、任意の1 <= i <= d に対して、点pの第i座標をp.i = xiと書く。特に、d = 2の とき、第1座標と第2座標を、x座標とy座標と呼ぶ。Rd上 の距離関数dist :R2_{× R}2 _{→ Rを、任意の点x, y, zに対して} 次の条件を満たすものとして定義する： (i)非負性：dist(x, y) >_{= 0;}

(ii)対称性：dist(x, y) = dist(y, x);

(iii)三角不等式：dist(x, y) + dist(y, z) >_{= dist(x, z);}

(iv)非退化性：dist(x, y) = 0⇔ x = y; p, q∈ Rd_{を任意の点とし、点pとqの間のL} ∞距離を次のよ うに定義する： dist∞(p, q) := max 1<_=i< =d {|p.i − q.i|}. また、任意の正整数α = 1, 2, . . .に対して、点pとqの間の Lα距離を次のように定義する： distα(p, q) :=1 d { _d ∑ i=1 |p.i − q.i|α }1/α . 以後、考えているLα距離が明らかならば、これを単に距離と よび、距離関数をdist(p, q)と書くことがある。 次に、任意のα = 1, 2, . . .または∞に対して、与えられた 点p∈ R2を中心とする半径r > 0のLα超球を Shapeα(p, r) :={ q ∈ Rd| distα(p, q) <_{= r }⊂=R}d と定義する。ここで、2次元の場合に、もしα = 2ならば、 Shape1(p, r)はpを中心とする半径rの円盤である。同様に、 α = 1ならば、Shape2(p, r)はpを中心とする1辺 √ 2rの正 方形を原点を中心に45度回転させて得られる図形(ダイアモン ド形)であり、α =∞ならば、Shape∞(p, r)はpを中心とす る1辺2rの正方形∏d_i=1[xi− r, xi+ r]⊂_=Rdである。 2. 2 軌跡データと軌跡パターン 文献[4]にしたがって、軌跡データを定義する。 ［定義1］ (R2上の)軌跡テキスト(テキストともいう)、及び は軌跡パターン(パターンともいう)、長さnとmの点列 T = t1, ..., tn∈ (R2)∗, P = p1, ..., pm∈ (R2)∗ である。ここに、各j >_{= 1}に対して、tj, pjは点である。 2. 3 近似軌跡パターン照合 ［定義2］ 正数をrとする。軌跡パターンPが軌跡テキストT に最大距離rで近似照合するとは、ある位置1 <_{= k <= n}が存 在して、すべての1 <_{= i <= m}に対して、 d_∞(pi, tk+i−1) <= r が成立することをいう。このとき、位置kをPのT における 出現位置という。 ［定義3］ 近似軌跡パターン照合問題とは、軌跡テキストT と 軌跡パターンP に対し、tj...tj+mとp1...pmがそれぞれの点 で近似照合するようなTの開始位置jを求める問題をいう。 2. 4 r-配 置 任意のα = 1, 2, . . .または∞と非負実数r >_{= 0}を考える。 入力点集合とは、m >_{= 1}個の点の集合P ={p1. . . pm}であ り、Pが定めるLα超球集合とは、Pの点を中心とする半径r

P1 P3 P2 P4 P5 P6 f1 f2 f3 f4 f5 f10 f9 f6 f7 f8 f11 f12 f13 図 1 L∞距離に対する r-配置 (r-arrangement)A∞(P, r)と面 fk のLα超球のなす集合

Shapeα(P, r) :={Shapeα(p1)、. . .、Shapeα(pm)}⊂_=R

である。ここで、RはRd_{上の閉領域の族である。}

Shapeα(P, r)と任意の点q∈ Rdに対して、それを含む超球 の中心点の集合、すなわち、qのPによるr-被覆集合(cover set)を、

Coverα(q; P, r) := { 1 <_{= i <= m | q ∈ Shape}α(pi, r)} (1)

= { 1 <_{= i <= m | dist}α(p, q) <_{= r }} (2) と定める。このとき、 全空間Rd_{は、 超球集合Shapeα(P, r)に} 含まれる超球の表面によって、いくつかの連続した領域に区切 られる。この区切られ方、すなわち、表面の組み合わせをPに よるr-配置(r-arrangement)と呼び、Aα(P, r)で表す。また、 連続した領域のそれぞれを面(face)と呼び、この配置における 面すべての集合をF aceα(P, r)で表す。任意の面fに対して、 それに含まれる任意の点は、同じr-被覆集合Coverα(q; P, r) をもつことに注意されたい。 二次元の場合には、α =∞の時、各超球は正方形であり、こ れらの正方形の各面で区切られた多角形領域が面となる。一方、 α = 2の時、各超球は円盤であり、対応する円で区切られた領 域が面である。一般に、超球の表面は曲線(または直線)であ り、曲線と曲線の交点を頂点(vertex)と呼び、頂点と頂点を結 ぶ面の境界線を辺(edge)と呼ぶ。以後、αが明確ならば、省略 することがある。 2. 5 多次元点位置検索データ構造 非負実数r >_{= 0}に対して、入力点集合P ={p1. . . pm}で定ま るr-配置Aα(P, r)を考える。任意の点q∈ Rdの配置Aα(P, r) における出現位置とは、点qを含む面f ∈ F aceα(P, r)のこ とを表す。ここで、点qの点位置検出問題を以下のように定義 する。 ［定義4］（点qの点位置検出問題） あ るr-配 置 Aα(P, r) が 存 在 す る 時 、任 意 の 点 q に 対 し て 、q の 出 現 位 置 f ∈ F aceα(P, r) s.t. q∈ fを求めよ。 図2に点qの点位置検出問題の例を示す。 最後に、本稿で考察する多次元点位置検索データ構造を以下 のように定義する。 ［定義5］（多次元点位置検索データ構造） 入力として、非負実 数r >_{= 0}、入力点集合P ={p1. . . pm}を受け取り、次の点位 置検出クエリP Lをサポートするようなデータ構造PL • P L(q): 点qの点位置検出問題を解き、その解fを出力 p1 p2 p3 p4 p5 f1 f2 f3 f4 f7 f5 f6 f8 f11 _f12 f10 f9 f13 q 図 2 点位置検出問題の例。点 q は面 f4に出現している する。 以後、このデータ構造PLを、多次元点位置検索データ構 造と呼ぶ。

2. 6 GPUとGPU汎用計算(GPGPU)

GPU(Graphic Processing Unit)とは、複数のStream

Multi-processor(SM)を備えた、コンピュータ上での画像処理に特化 したハードウェアである。それぞれのSMには16∼32個の SIMD演算コアが内蔵されており、各演算コアは、SMに搭載 された低容量・高速な共有メモリと、GPUボード上に搭載さ れた、大容量・低速のグローバルメモリを用いて並列に計算を 行うことで、高速に画像処理演算を行う。 このような多数の演算コアと、メモリ階層を利用した高速な データアクセスによる高い並列演算能力を、画像処理以外の汎用 的な計算に用いる動きが増えている。これはGPGPU(General Purpose computing on GPU, GPU汎用計算)と呼ばれ、近年 注目を集めている技術である[7]。

3.

空間パターン索引とビット並列計算を用いた

アルゴリズム

次に3.節では、2次元空間のL_∞距離に関するトラジェクト リパターン照合問題について考察する。はじめに、よく知られ た文字列照合法であるShift-And法を取り上げ、この手法に基 づく素朴なトラジェクトリパターン照合法を与える。次に、こ の手法を改良し、点位置検索データ構造とビット並列計算を組 み合わせて、より高速なトラジェクトリパターン照合アルゴリ ズムFastShift-AndTrajMatchを提案する。 3. 1 Shift-And法を用いた素朴なアルゴリズム はじめに、よく知られた文字列照合法であるShift-And法を 取り上げ、この手法に基づく素朴なトラジェクトリパターン照 合法を与える。

Shift-Andアルゴリズムは、Baeza-YatesとGonnet [1]及び

WuとManber [10]によって提案された、ビット並列手法に

基づく文字列照合アルゴリズムである。パターンを受理す る非決定性有限オートマトン(N F A)にテキストを1文字ず つ入力し、その各状態、及び状態遷移をビットで表すことに よってO(m + n)時間で照合を行う。Algorithm 1にShift-And

法を示す。このShift-And法のトラジェクトリ照合への拡張

は、直感的には照合の条件を変えるだけで可能である。

Algo-rithm 2にShift-And法に基づく近似軌跡パターン照合アルゴ リズムNaiveShift-AndTrajMatchを示す。

Algorithm 1 Shift-And法に基づく文字列照合アルゴリズ ムShift-AndStringMatch 1: Preprocessing 2: P← p1p2...pm ▷ piは P 中の i 番目の文字 3: T ← t1t2...tn ▷ tjは T 中の j 番目の文字 4: for c∈ Σ do 5: B[c]← 0m 6: end for 7: for i∈ 1...m do 8: B[pi]← B[pi]| 0m−i10i−1 9: end for 10: Searching 11: D← 0m 12: for pos∈ 1...n do 13: D← ((D << 1) | 0m−11) & B[tpos] 14: if D & 10m−1 |= 0m_then

15: report an occurrence at pos− m + 1

16: end if 17: end for Algorithm 2 Shift-And法に基づく素朴なトラジェクトリ照 合アルゴリズムNaiveShift-AndTrajMatch 1: Preprocessing 2: P← p1p2...pm ▷ piは P 中の i 番目の文字 3: T ← t1t2...tn ▷ tjは T 中の j 番目の文字 4: Searching 5: D← 0m _{▷状態集合マスク} 6: for pos∈ 1...n do ▷現在の入力点に対して文字マスクを作る 7: B← 0m _{▷文字マスク} 8: for i∈ 1...m do 9: if d_∞(pi, tpos) <= r then 10: B← B | 0m−i₁₀i−1 _{▷ i番目のビットを 1 にする} 11: end if 12: end for 13: D← ((D << 1) | 0m−11) & B 14: if D & 10m−1 _{|= 0}m_then

15: report an occurrence at pos− m + 1

16: end if 17: end for このアルゴリズムの問題点を述べる。Shift-And法では、ビッ トマスク表は|Σ|個の文字のそれぞれに対して前処理で1回だ け作る。しかし、近似軌跡パターン照合においては、前処理の 段階では、可能な入力点の候補は無限個存在するため、あらか じめビットマスクパターンを準備することができない。そのた め、NaiveShift-AndTrajMatchアルゴリズムでは各入力点に対 して毎回ビットマスクパターンを作っている。したがって、実 行時間はO(mn)時間になる。これは素朴な文字列照合法の実 行時間O(mn)と同じであり、効率が悪い。そこで、次節では、 前処理でビットマスクパターンを1回だけ作る方法を提案する。 3. 2 点位置検索データ構造 本小節では、提案アルゴリズムで用いる、2次元のL_∞距離 に対する点位置検索データ構造の実現方法について述べる。こ れは、入力点に対してビットマスクを検索するために用いられ る。はじめに、長さMの軌跡Pが含む点すべてに対して、それ らを中心とする半径rの正方形の全体D(P ) = {D0, . . . , DM} を考える。ここに、各i = 0, . . . , M に対して、Di∈ D(P )は Pのi番目の点P (i)を中心とする半径rの正方形を表す。 次に、これらの正方形全体の平面上へのL∞距離に対する r-配置(r-arrangement)A(P )を考える(図1参照)。この時、次 の補題が知られている[2]。 ［補題1］ Pを2次元平面上のm個の点からなる軌跡パター ンとする。A(P)をPで決まる配置とする。このとき、次が成 立する。 （1） A(P)の頂点数はO(m2)個である。 （2） A(P)の辺の数はO(m2_{)個である。} （3） A(P)の面の数はO(m2_{)個である。} 配置A(P )の中で、面の全体をF(P )で表す。これらの面を f1, . . . , fK(K >_{= 1)}と番号づける。表記として、面fk∈ F(P ) が、D(P )中の正方形Di ∈ D(P )に含まれているときに、 fk∈ Diと書くことにする。 アルゴリズムでは、面の集合F (P )を、次の演算をサポート する点位置検索(point location)データ構造PLに格納する： • 点位置検索演算P ointLocation(p)：任意の点p∈ R2を 入力として、それを含む面fk∈ F (P )を返す。そのような面が なければ、N ullを返す。 • 集合取得演算BoxList(fk)：任意の面fkを入力として、 それを作る正方形Diの添字iの集合をリストとして返す。す なわち、 BoxList(fk) ={ 1 <_{= i <= m | f}k∈ Di, Di∈ D(P ) } そのような正方形が1つもなければ、∅を返す。 2次元平面上のL_∞距離に関する点位置検索データ構造につ いては、次のようなデータ構造が容易に構築できる。 ［補題2］（山本他[11]） 二次元平面上のm個の正方形の配置 に対して，それをO(⌈m_w2⌉)語の領域とO(m log m +⌈m_w2⌉)構 築時間を用いて，点位置検索演算をO(log m +⌈m_w⌉)時間で実 現する索引構造が存在する．

証明: 基本的なアイデアは，m個の正方形の配置A(P )をx軸 とy軸に射影して得られる区間の集合Segx(P )とSegy(P )

に 分 解 し ，各 面 を そ れ ら の 区 間 の 直 積 の 直 和 で 表 す 事 で あ る ．索 引 構 築 演 算 で は ，全 正 方 形 の 両 端 点 のx座 標 を ， x1<· · · < xh(h <_{= 2m)}のようにO(m log m)時間でソートし，

区間の集合Segx(P ) = { Ix = [xi, xi+1)| i = 1, . . . , h }を計

算する．ここに，x0=−∞, xh+1=∞である．次に，各区間

sx∈ Segx(P )について，sxを含む正方形の集合SegListx(sx)

を，長さmのビットマスク表現Bx[sx]∈ {0, 1}mで記録し，

ビットマスク配列Bx を構築する．同様にy軸上の区間の集

合をSegy(P )と，ビットマスク配列Byを構築する．このと

き領域は，Segx(P ), Segx(P )にO(m)語を，Bx, Byの保持に

O(⌈m2

w ⌉)語をそれぞれ用いる．点位置検索演算P ointLocation

は，入力点のx座標とy座標を含む区間sxとsyを，それぞ

れ，Segx(P )とSegy(P )上の二分探索で求める．得られた2

Algorithm 3 Shift-And法に基づく提案のトラジェクトリパ ターン照合アルゴリズムFastShift-AndTrajMatch 1: Preprocessing 2: P← p1p2...pm ▷ piは P 中の i 番目の文字 3: T ← t1t2...tn ▷ tjは T 中の j 番目の文字 4: for all fkdo 5: B[fk]← 0m ▷領域 fkに対するビットマスクの初期化 /*空間パターン索引の作成 */ 6: list← BoxList(fk) ▷正方形の添字をリストとして受け取る 7: while list.length > 0 do 8: i← list.top ▷ リストの先頭の添字を取り出し、i に代入 する 9: erase(list.top) ▷リストの先頭要素を削除する 10: B[fk]← B[fk]| 0m−i10i−1▷正方形 Diに対応するビッ トを 1 にする 11: end while 12: end for 13: Searching 14: D← 0m _{▷状態集合マスク} 15: for pos∈ 1...n do

16: fk← P ointLocation(tpos) ▷点 tposを含む面 fkを見つける

17: D← ((D << 1) | 0m−11) & B[fk] ▷状態遷移を行う

18: if D & 10m−1 |= 0m_then

19: report an occurrence at pos− m + 1

20: end if 21: end for りO(log m +⌈m_w⌉)時間で実行する．以上より示された． □ 他 の 距 離 に 関 す る 点 位 置 検 索 デ ー タ 構 造 に つ い て は Ma-tousek [8]を参照されたい。 3. 3 提案アルゴリズム 次に、3. 1節の素朴なアルゴリズムを改良し、前節の点位置検 索データ構造とビット並列計算を組み合わせて、より高速なトラ ジェクトリパターン照合アルゴリズムFastShift-AndTrajMatch を提案する。Algorithm 3に提案アルゴリズムを示す。このア ルゴリズムは、前節の素朴なアルゴリズムと同様に、Shift-And 法に基づくが、マスク表の構築法が大きく異なっている。この アルゴリズムは、点位置検索データ構造を用いることで、前処 理でマスク表を構築することが出来る。 提案アルゴリズムは次のように動作する。はじめに前処理時 において、与えられたトラジェクトリパターンを、配置A(P ) から得られる面の集合F (P )と、Pの各点が構成する正方形を 辺ラベルとする非決定性有限オートマトン(NFA)に変換する。 続いて、F (P )とNFAからShift-And法と同様のマスク表を作 成する。 実行時には、Shit-AndMatchアルゴリズムと同様にして、こ のマスク表を用いて入力点を順に読み込みながら、NFAの対応 する状態遷移を行い、パターン照合を効率よく実行する。この 時、点位置検索データ構造を用いて、各入力点を含むボックス の集合に対応するマスク表を検索することで、高速化を行う。 このアルゴリズムの計算量は以下の様に与えられる[11]。 ［定理1］ 長さmの軌跡パターンP と長さnの軌跡テキスト 表 1 並列点位置検索データ構造の前処理アルゴリズム 入力: パターン点集合 P ={p1, . . . , pm}⊂_=Rd,正実数 r > 0. タスク: Pの各点を中心とするボックス (正方形) をデータ構造に格納 （ 1 ） ホスト上で、マスタープログラムは、点集合 P をホストから GPUに送信し、GPU のグローバルメモリ上におく。この点集合は、c 個ずつのブロック P1, . . . , P_⌈m/c⌉に分割される。 （ 2 ） ホスト上で、マスタープログラムは、c 個のコア C1, . . . , Cc を GPU 上に起動する。 （ 3 ） GPU 上では、1   <_{= i   <= c に対して、各コア C}iはそれぞ れ並列に、自分の担当するブロック Piを、グローバルメモリから自分 の属するコアグループの共有メモリに読み込む。 （ 4 ） GPU 上 で 、各 コ ア Ci は ブ ロック Pi に 含 ま れ る 各 点 {pi,1, . . . , pi,⌈m/c⌉} に対して、それらを中心とする 1 辺が 2r のボッ クス (正方形) Bi,1, . . . , Bi,⌈m/c⌉を計算する。 （ 5 ） GPU 上 で 、各 コ ア Ci は 、各 次 元 1 <= j <= d に 対 し て 、こ れ ら の ボック ス の 辺 に よ る 、次 元 j の 座 標 軸 の 領 域 分 割

Ai,j = { pi,k.j| pi,k ∈ Pi, 1 <= k <= ⌈m/c⌉ }⊂_{=R を求めて、}1次 元の順序辞書 (predecessor dictionary) Aj= Ai,jに格納する。ここ

に、最小値と最大値を除く、Ajの異なるキーの数を 1 <_{= n(i, j) <= 2m} とおくと、昇順にソートした A のキーは x0<· · · < xn(i,j)+1とかけ る。ここで、各ランク 0 <_{= k <= n(i, j) に対して、辞書 A}jには、キー xkに対する値として、m-ビットマスク Mki,j∈ {0, 1}mを格納する。 ここに、マスク M_ki,jは、値 xkの左に隣接する区間 Ik= [xk−1, xk] を第 j 座標の辺が含むような P のボックスの番号の集合を表すビット マップである。 （ 6 ） GPU 上で、各コアはクエリを待って、待機する。 T、正数rに対して、アルゴリズムFastShift-AndTrajMatchは O(m log m + m⌈m/w⌉)前処理時間とO(m⌈m/w⌉)領域を用い て、O(n(log m +⌈m/w⌉))実行時間で、トラジェクトリパター ン照合問題を解く。ここでwはワード長である。

4.

多次元点位置検索データ構造の

GPU

を用いた並列化

本節では、2次元空間のL∞距離に関して、GPU上で高速 に動作する多次元点位置検索データ構造のアルゴリズムを考察 する。GPUでは、共有メモリを共有する複数のコアがグルー プになっており、各コアは自分のローカルメモリを所有すると 仮定する。さらに、ホストとすべてのコアがグローバルメモリ にアクセスできると仮定する。以降では4つの異なるアプロー チからアルゴリズムを考察する。 準備として、実数配列Aに対して、1次元の順序辞書 (pre-decessor dictionary)とは、実数の集合を格納するデータ構造 であり、挿入(insert)と、削除(delete)、前者(predecessor)、 後者(successor)の4つのクエリを実現するものである。ここ に、実数qに対する前者クエリは、q以下の最大のA中の要素 を返し、qに対する後者クエリは、q以上の最小のA中の要素 を返す。ただし、これらのクエリが常に回答を返すために、辞 書は元のデータに含まれない最小値minと最大値maxを必ず 含むと仮定する。例えば、C++言語のSTLライブラリのmap は順序辞書である。

表 2 並列点位置検索データ構造の実行時処理アルゴリズム 入力: クエリ点 q∈ Rd_. 出力: クエリ点 q を含むボックスの番号の集合を表す m-ビットマス ク M ∈ {0, 1}m_. （ 1 ） ホスト上で、マスタープログラムはクエリ点 q を受け取ると、 GPU上の c 個のコア C1, . . . , Ccに q を、グローバルメモリを用いて 一斉送信する。 （ 2 ） GPU 上で、各 1 <_{= i <= ⌈m/c⌉ に対して、コア i はグローバ} ルメモリから q を、共有メモリ経由で局所メモリに並列に受け取る。 （ 3 ） GPU 上で、各コア i は、各次元 1 <_{= j <= d に対して順に、} データを格納した順序辞書 Ai,j上で、クエリの第 j 座標 q.j をキーと する後者クエリを実行する。これにより、第 j 座標が q.j を含むボッ クス名 1 <_{= k <= m の集合を表す m-ビットマスク M}i,jが返される。 （ 4 ） GPU 上で、各コア i は、すべての次元に対するビットマスク の論理積 Mi_{= M}i,1_{∧ · · · ∧ M}i,d_{∈ {0, 1}}m_{を求めて、共有メモリ} とグローバルメモリを経由して、並列にホストへ送る。 （ 5 ） ホ ス ト 上 で 、マ ス タ ー プ ロ グ ラ ム は 、GPU 上 の c 個 の コア C1, . . . , C_⌈m/c⌉から、グローバルメモリ経由でビットマスク M1_{, . . . , M}⌈m/c⌉_{を受け取る。最終的な解として、論理和 M}1_{∨ . . . ∨} M⌈m/c⌉∈ {0, 1}mを返す。 4. 1 GPU上の並列点位置検索データ構造のアルゴリズム ここでは、3. 2節の補題2で与えた点位置検索データ構造を、 並列化したデータ構造を与える。この並列データ構造の高レベ ルのアルゴリズムは、4つのアプローチすべてで共通している。 表1と表2に、並列点位置検索データ構造の前処理と実行時 処理の並列プログラムを示す。ここに、パラメータとして、正 整数cとdが与えられると仮定する。cはGPUが持つレジス タの語長wの定数倍、すなわちc = ℓw(ℓは定数)と仮定する。 上のアルゴリズムを実装するには、各コアの上で1次元の順 序辞書の後者クエリを解く方法を与えれば良い。これは、表1 のステップ5と表2のステップ3に対応する。表2のステップ 4と5におけるビットマスクの論理積と論理和を用いた集約演 算は、平衡二分木スキームによる並列化を用いると、総ブロッ ク数s =⌈m/c⌉に対してt = O(log s) = O(log m)並列ステッ プで実行できる。仮定からビットマスクの長さはワード幅wの 定数倍となるので、論理積と論理和それぞれは１演算当たり定 数時間で実行できる。 以下に、表1のステップ5と表2のステップ3に対応するブ ロック内の順序辞書検索の部分問題を複数のコアを用いて高速 に解くためのアプローチを4つ示す。各アプローチは、ブロッ ク内の小問題の解き方に設計上の相違点がある。これ以後、サ イズcのブロックをマクロブロック(大ブロック)と呼ぶ。アプ ローチによっては、各マクロブロックをさらにサイズdのマイ クロブロック(小ブロック)に再分割することがある。 4. 2 データ構造A：マイクロ線形走査 このアプローチでは、ある小さな定数h(1 <_{= h <= c)}をあら かじめ選び、前処理において、マクロブロック内のc個の入力 点をさらに、h個ずつからなる⌈c/h⌉個のマイクロブロックに 分割する。次に、各i = 1, . . . ,⌈c/h⌉に対して、第i番目のマ イクロブロックが含むh個の入力点pi1, . . . , pihによって定めら Algorithm 4マイクロ線形走査 1: array A← ai 1, . . . , aih−1 2: procedure Search(x) 3: for i in 0..h− 1 do

4: if x exists in A[i] then return i;

5: end for

6: return−1

7: end procedure

Algorithm 5逐次版の二分探索

1: procedure Search(x: real, A = A[1] · · · A[c]: array of reals)

2: A[0] = min; A[c + 1] = max, 3: L← 0; R ← c; 4: while L < R do ▷不変条件：x∈ [L, R) 5: M = (L + R)/2; 6: if x < A[M ] then R = M ; 7: else L = M ; 8: end while 9: return M ; 10: end procedure れる区間ai1, . . . , a i h−1を、長さhの配列Aiに格納する。実行 時には、1次元のクエリに対して、配列上を線形探索すること で、1次元の点位置検索問題を解く。その後、それらの解を集 約し、目的のクエリに対する解を得る。 このアプローチの簡単な実装をAlgorithm 4に示す。この コードの通り、実装は極めて単純である。また、このアプロー チの性能は、マクロブロックのサイズhの選び方に大きく依存 すると考えられる。 4. 3 データ構造B：分岐レス二分探索 このアプローチでは、マクロブロック内のc個の入力点はそ れ以上分割しない。前処理では入力点によって定められる区 間a1, . . . , ac−1をソートして長さcの配列Aに格納し、1つの ソート済みブロックとして扱う。実行時には、与えられたクエ リ値xに対して、このソート済み配列上で二分探索を行い、1 次元の点位置検索問題を解く。 このアプローチを単純に実装した場合の問題点は、分岐演算 を使用することである。アルゴリズム5に通常の二分探索のア ルゴリズムを示す。このコードから、4行目から8行目に条件 文が存在し、分岐が生じていることがわかる。 ほとんどのGPUにおいて、このような条件分岐は、if節と else節の両方を同時に実行し、条件が成立しなかった方の結果 を破棄することで実装されている。したがって、1つの物理コ アしか使えなければ、条件分岐の実行に2倍の時間がかかる場 合があると予想される[12], [13]。この問題を解決するために、 以下では、分岐演算を用いずに二分探索をGPU上で実現する 方法を与える。 基本的なアイディアは、二分探索の計算における各種の分岐 演算を、ビット並列計算と配列アクセスを用いて除去すること である。実装の都合上、入力配列Aの要素およびクエリの値x は、w− 1ビットの符号なし正整数と仮定する。これは、入力

Algorithm 6分岐レス比較演算 1: procedure Compare(x, y ∈ {0, 1}w_{: 符号なし w ビット整数)} 記法： w ビットレジスタ Z を、MSB と LSB をそれぞれ左と右 として、Z = bw−1· · · b0のように表記する。 出力： w ビット整数 V . もし x < y ならば V = 0, x >_{= y ならば} V = 1. 2: X = (X| 10w−1_{); Y = (Y | 0}w_{); ▷ Xと Y の最上位ビット} は、1 と 0. 3: Z = X− Y ▷ Zの最上位ビットは、X < Y のとき 0、 X >_{= Y のとき 1.} 4: return V = (Z & 10w−1)≫ (w − 1); ▷ 最上位ビットを取 り出す 5: end procedure Algorithm 7分岐レス二分探索

1: procedure Search(x: real, A = A[1] · · · A[c]: array of reals)

2: A[0] = min; A[c + 1] = max, 3: L← 0; R ← c + 1;

4: for level = 0, . . . , log c + 1 do ▷不変条件：x∈ [L, R)

5: M =⌊(L + R)/2⌋;

6: E[0] = L; E[1] = M ; E[2] = R;

7: of f set = Compare(x, A[M ]); ▷ Test if x < A[M ] 8: LI = 0 + of f set; 9: RI = 1 + of f set; 10: L = E[LI]; R = E[RI]; 11: end for 12: return M ; 13: end procedure の実数を、定められた解像度と原点で変換することで行われる。 Algorithm 6に、分岐レス比較演算のアルゴリズムを示す。 このアルゴリズムは、レジスタ上の減算とビット演算を用いて、 比較演算“x− y”の結果の正負を、それぞれ、整数1または0 で返す。レジスタのビット幅をwとおく。このアルゴリズムで は、入力の２つの数をw− 1ビット非負整数XとY ∈ {0, 1}w で表す。さらに、Xの最上位ビットに1を立て、Y の最上位ビッ トは0とする。このとき、減算Z = X− Y を行うと、X < Y のときにはキャリーの桁下がりにより、差Zの最上位ビットは 1から0になる。逆にX >_{= Y} のときには差Zの最上位ビット は1のままである。したがって、Zの最上位ビットを取り出せ ば、比較演算の結果を取り出すことができる。以上の計算は、 一つの減算と定数個のビット演算で実行できる。 Algorithm 7に、分岐レス二分探索のアルゴリズムを示す。こ のアルゴリズムは、上記の分岐レス比較演算をサブルーチンとし て用いて、c個の整数値を昇順に格納した配列A = A[1]· · · A[c] 上の二分探索を模倣し、与えられた整数値x以上の最小の要素 を見つける。 このアルゴリズムは、ループの不変条件としてx∈ [L, R)が 常に成立するように、探索対象の区間を表す変数LとRを管 理する。配列には両端に仮想的な最小値A[0] = minと最大値 A[c + 1] = maxを加えているので、初期値L = 0, R = c + 1 では仮定は常に成立する。さらに、ループの最初で不変条件 x∈ [L, R)が成立すると仮定すると、７行目から10行目の コード

of f set = Compare(x, A[M ]); LI = 0 + of f set; RI = 1 + of f set; L = E[LI]; R = E[RI]; を実行後に、もしx < A[M ]ならばL := LとR := Mが実行 される。逆にx >_{= A[M ]}ならばL := MとR := Rが実行され る。したがって、常に不変条件x∈ [L, R)が満たされることが わかる。 さらに、このアルゴリズムは停止条件についても分岐命令、 つまり元の逐次型二分探索のwhile文の条件分岐をもたず、定 数回のループで置き換えられている。これは、L + 1 < Rが成 立する間は、ループの実行により必ずR− (L + 1)が真に減少 することと、一度条件L + 1 = Rが成立すると、それ以後の ループの実行においてもL + 1 = Rが成立し続けることが保証 されることからわかる。したがって、ループのlog c + 1回の実 行後には、必ずL + 1 = Rが成立し、Lがクエリ値xの前者 クエリの解となっている。後者クエリの解は、L = xのときL であり、それ以外の時はL + 1である。 以上より、次の補題が示せる。 ［補題3］ アプローチBに関する上のアルゴリズムは、マクロ ブロックにおいて、c個のwビットの非負整数値の配列A[1, c] 上で、wビットの非負整数であるクエリ値xに関する二分探 索を、分岐演算を用いずに正しく実行する。このとき、一つの マクロブロックにおいて、使用するメモリ量はO(c)語であり、 計算時間はO(log c)時間である。 4. 4 データ構造C:ハイブリッドアプローチ これは、データ構造Aのマイクロブロックのアイディアと データ構造Bの分岐演算なしの二分探索のアイディアの複合的 アプローチである。このアプローチでは、データ構造Aと同様 に、マクロブロック内のc個の入力点を、h個ずつのマイクロ ブロックにさらに分割する。各マイクロブロック内で、データ 構造Bの二分探索の分岐演算なしの模倣アルゴリズムを用いて 探索を実行する。 このアプローチでは、データ構造Aのアプローチと同様に、 性能向上にはマクロブロックのサイズhの選び方が重要である。 4. 5 データ構造D：数値演算を用いたビット内並列比較

このアプローチは、Fredman及びWillardによるFusion

Tree [5]の手法にに基づくものである。はじめに、先のすべて のデータ構造と同様に、入力値をサイズcのマクロブロックに 分割する。あらかじめマクロブロック内のc個の値は p0<_{= . . . <= p}c−1 のように昇順にソートする。 次に、最小値pmin= p0との差分をとって、 次のように各値 を差分圧縮し、ベクトル q0:= p0− pmin, . . . , 1c−1:= pc−1− pmin.

を、最小値pminと一緒に格納する。 次にこれらの差分値を、左から一つのレジスタに入るだけと り、語長w(の定数倍)の長さをもつマイクロブロックに詰め込 む。各差分値は、これらのうちで最大のビット長δに合わせて、 同じビット数で表す。 実行時には、与えられたクエリ値に対して、各マイクロブ ロックで独立に1次元の点位置検索問題を、次の述べる数値演 算を用いたビット内並列比較を用いて解く。 • 初期値として、一個のマイクロブロックを格納したレジ スタXに、δビットのk個の差分値が格納されていると仮定す る。このとき、δビットフィールドの境界のビットには1を立 てる。 • ク エ リx か ら 、最 小 値pmin を 引 い て 差 分 ク エ リ 値 y := x− pminを求める。 • 差分クエリ値yが、δビットに入らないようなら、最大 の入力差分値より大きいことがわかるので、直ちに解を返す。 • もし入るようなら、yのコピーをk個詰め込んだレジス タY を、1回の乗算命令を用いて作成する。このとき、δビット フィールドの境界のビットには0を立てる。これは、各kビッ トフィールドの最下位ビットに1をたてた定数レジスタと1個 のyを最下位に詰め込んだレジスタの乗算で求まる。 • レ ジ ス タ X か ら レ ジ ス タ Y を 減 算 し 、レ ジ ス タ Z := X− Y を得る。これにより、対応する各δビットフィー ルド同士の減算が行われる。このとき、X とY の境界フィー ルドビットが、それぞれ、1と0であることから、減算結果の 境界フィールドビットは、減算結果が負ならば0となり、正な らば1となる。 • レジスタZの境界フィールドビット値を、適切なマス クをかけて取り出す。単調性から、下位から上位へ見て、境 界フィールドビットは最初は減少することはない(0から1の 向きに変化する)。したがって、定数回のビット演算と1回の

MSB(Most Significant Bit)演算を用いて、この変化点にビッ

ト1をたてたレジスタを解として返す。

5.

お わ り に

本論文では、上記のようなトラジェクトリデータに対する高 速なパターン照合法と、関連する点位置検索データ構造と呼ば れるデータ構造のGPUを用いた高速化について考察した。 3.節では、2次元空間のL∞距離に関するトラジェクトリパ ターン照合問題について考察した。点位置検索データ構造と ビット並列計算を組み合わせて、高速なトラジェクトリパター ン照合アルゴリズム FastShift-AndTrajMatchを提案した。こ のアルゴリズムは、ワード長wの計算機上で、長さmの軌跡 パターンと長さnの軌跡テキストに対して、O(m3⌈m_w⌉)前処 理時間とO(m3⌈m_w⌉)領域を用いて、O(n(log m +⌈m_w⌉))実行 時間で動作する。 4.節では、多次元点位置検索データ構造をGPU上で並列実 装するための効率良い手法について考察した。提案手法は、与 えられた入力点集合を定数cのサイズに分割し、それぞれをマ クロブロックとして、GPU上の独立したコア集団で並列に解 くことを考えた。各マクロブロック内の計算については、4つ のアプローチを提案した。 今後の課題について述べる。3.節のトラジェクトリパターン 照合アルゴリズムについては、その応用を考えることが課題で ある。また、4.節の並列データ構造については、時間の都合に よりGPU上の実装が出来なかった。これに関して、実装と実 験による性能評価が課題である。 文 献

[1] Ricardo Baeza-Yates and Gaston H. Gonnet. A new ap-proach to text searching. Communications of the ACM,

Vol. 35, No. 10, pp. 74–82, October 1992.

[2] Mark de Berg, Otfried Cheong, Marc van Kreveld, and Mark Overmars. Computational Geometry: Algorithms and

Applications. Springer-Verlag, Santa Clara, CA, USA, 3rd

ed. edition, 2008.

[3] Thomas H. Cormen, Clifford Stein, Ronald L. Rivest, and Charles E. Leiserson. Introduction to Algorithms. McGraw-Hill Higher Education, 2nd edition, 2001.

[4] Thomas Eiter and Heikki Mannila. Computing discrete Fr´echet distance. Technical Report Technical Report CD-TR 94/64, Christian Doppler Laboratory for Expert Sys-tems, TU Vienna, Austria, 1994.

[5] Michael L. Fredman and Dan E. Willard. Surpassing the information theoretic bound with fusion trees. Journal of

Computer and System Sciences, Vol. 47, No. 3, pp. 424 –

436, 1993.

[6] HONDA. 通行実績情報マップ. http://www.honda.co.jp/ internavi/service/disastermap/.

[7] David Kirk. Nvidia cuda software and gpu parallel com-puting architecture. In Proceedings of the 6th International

Symposium on Memory Management, ISMM ’07, pp. 103–

104, New York, NY, USA, October 2007. ACM.

[8] Jiˇr´ı Matouˇsek. Geometric range searching. ACM

Comput-ing Surveys, Vol. 26, No. 4, pp. 422–461, December 1994.

[9] Gonzalo Navarro and Mathieu Raffinot. Flexible Pattern

Matching in Strings: Practical On-line Search Algorithms for Texts and Biological Sequences. Cambridge University

Press, New York, NY, USA, 2002.

[10] Sun Wu and Udi Manber. Fast text searching: Allowing errors. Communications of the ACM, Vol. 35, No. 10, pp. 83–91, October 1992. [11] 山本雅大, 栗田和宏, 笹川裕人, 有村博紀. D-032 2 次元軌跡デー タに対する高速なパターン照合アルゴリズム (d 分野:データベー ス, 一般論文). 情報科学技術フォーラム講演論文集, Vol. 13, No. 2, pp. 153–154, 8月 2014. [12] 橋口拓哉, 青野寛之, 豊永昌彦, 村岡道明. GP-GPU を用いた高 速並列論理シミュレータ. 情報処理学会研究報告. SLDM, [シス テム LSI 設計技術], Vol. 2015, No. 12, pp. 1–6, 1 月 2015. [13] 坂本博之, 佐々木仁, 雫譲, 黒木修隆, 廣瀬哲也, 沼昌宏. B-004

ウェーブレット変換に基づく学習型超解像の GPU による高速 化手法 (並列処理,b 分野:ソフトウェア). 情報科学技術フォーラ ム講演論文集, Vol. 11, No. 1, pp. 173–174, 9 月 2012.