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. は じ め に

長さ n の軌跡テキストに対して、O(m) 領域を用いて、O(mn)実行時間で軌跡パターン照合を行う。

GPS やスマートフォンなどの位置センサーの普及により、自次に、山本らは、この手法を改良し、点位置検索データ構造

動車や、人、野生動物などの移動データに代表される、トラジェとビット並列計算を組み合わせて、より高速なトラジェクトリ

クトリデータの種類と量が増え続けている。また、これらの応パターン照合アルゴリズム FastShift-AndTrajMatch を提案し用として、これに関連する多次元点位置検索データ構造と呼ばた [11]。このアルゴリズムは、点位置検索データ構造を用いるれるデータ構造は、大規模データ解析や、時空間情報検索、コ

ことで、次のように動作する。

ンピューターグラフィックスなどの分野で、様々な応用を持つ。

はじめに、前処理において、軌跡パターン P を受け取り、

例えば、HONDA は、カーナビを搭載した車の移動データを元

点位置検索データ構造を構築する。実行時には、軌跡テキに、現在通行が可能な道路の情報を示した、通行実績情報マッスト T 上で、点位置検索データ構造を用いて、各入力点をプを公開している [6]。これらの技術において、位置センサーや含む面の集合を求め、対応する NFA の状態遷移を行う。こポイントクラウドなどのセンサー技術の発展により、入力点集のアルゴリズムは、ワード長 w の計算機上で、アルゴリズ

合データのサイズは年々増大しており、高速な位置検索データ

ム FastShift-AndTrajMatch は、長さ m の軌跡パターンと長さ構造が求められている。

n の軌跡テキストに対して、O(m log m + m⌈m/w⌉) 前処理時

そこで本論文では、上記のようなトラジェクトリデータに対間と O(m⌈m/w⌉) 領域を用いて、O(n(log m + ⌈m/w⌉)) 実行

する高速なパターン照合法と、点位置検索データ構造の GPU

時間で軌跡パターン照合を行う。この計算時間は、素朴なトラを用いた高速化について考察する。本論文で取り組む問題は以ジェクトリパターン照合法よりも高速である。

下の通りである。

1. 1 トラジェクトリパターン照合問題
第 1 の問題として、3. 節では、2 次元空間の L∞ 距離に関するトラジェクトリパターン照合問題について考察する。

1. 2 多次元点位置検索データ構造の GPU 上の並列実装
第 2 の問題として、4. 節では、3. 節で用いている点位置検索
データ構造が、巨大なサイズの軌跡パターンに対して、現実的な時間で処理を終えることが出来ない問題について考察する。

はじめに、よく知られた文字列照合法である Shift-And 法を

2 次元の点位置検索データ構造として、3. 節では、高々m 個

取り上げ、この手法に基づく素朴なトラジェクトリパターン照の正方形の配置による平面分割に対して、O(m2 ) 語の領域と合法を与える。文字列上の Shift-And 法は、パターンの前処理

O(m2 log m) 前処理時間を用いて、1 個の点位置検索クエリにによりマスクを構築して、高速化を実現していた。それに対しO(log m) 時間で答えるデータ構造を用いている。これは CPUて、 Shift-And 法を用いた素朴なトラジェクトリ照合では、入の逐次計算を用いることを前提としているが、巨大な m に対力が与えられないとパターンを前処理することができない。そする効率の良い実現は難しい。そこで、この手法を GPU を用のため、この素朴なアルゴリズムでは、実行時に入力点が与えいた並列計算によって高速化することを考える。

られるたびにマスクを生成し、照合を行う。ワード長 w の計算

本章で用いる GPU(Graphic Processing Unit) は、数 100〜

機上で、この素朴なアルゴリズムは、長さ m の軌跡パターンと数 1000 の多数のコアを持つ並列計算デバイスである。GPU は

元々画像処理用に開発されたが、現在では一般的な計算の高速化にも広く用いられており、その並列化技術は GPGPU(General

Purpose computing on GPU, GPU 汎用計算) と呼ばれている。
提案手法は、与えられた入力点集合を定数 c のサイズに分割
し、それぞれをマクロブロックとして、GPU 上の独立したコア集団で並列に解くことを考える。各マクロブロック内の計算については、次の 4 つのアプローチを提案した。

•

データ構造 A：マイクロ線形走査•

データ構造 B：分岐レス二分探索•

データ構造 C：ハイブリッドアプローチ•

データ構造 D：数値演算を用いたビット内並列比較

各アプローチに共通する考えは、(1) データを小さなブロッ(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 {|p.i − q.i|}.
1<
= i<
=d

また、任意の正整数 α = 1, 2, . . . に対して、点 p と q の間のLα 距離を次のように定義する：

クに分割し、個々のブロックを GPU のコアで並列に実行することと、(2) その際に GPU の特性を考慮して、分岐処理を減

distα (p, q) :=

1
d

{ d
∑

}1/α
|p.i − q.i|α.

らすことで高速化を行うことである。本論文では、上記の 4 つ

i=1のアプローチについて議論する。上記のアプローチの実装と計以後、考えている Lα 距離が明らかならば、これを単に距離と算機実験による性能評価は、今後の課題である。

よび、距離関数を dist(p, q) と書くことがある。

1. 3 本論文の構成
本論文の構成を述べる。はじめに、2 章では、以降の章に必要な定義や概念を導入した後、3. 節では、2 次元空間の L∞ 距離に関するトラジェクトリパターン照合問題について考察する。

次に、任意の α = 1, 2, . . . または ∞ に対して、与えられた
点 p ∈ R2 を中心とする半径 r > 0 の Lα 超球を
dShapeα (p, r) := { q ∈ Rd | distα (p, q) <= r }⊂=R

4. 節では、多次元点位置検索データ構造を GPU 上で並列実装と定義する。ここで、2 次元の場合に、もし α = 2 ならば、するための効率良い手法について考察する。最後に、5. で本論

Shape1 (p, r) は p を中心とする半径 r の円盤である。同様に、√
α = 1 ならば、Shape2 (p, r) は p を中心とする 1 辺 2r の正文の結論と、今後の課題を述べる。

1. 4 関 連 研 究

方形を原点を中心に 45 度回転させて得られる図形 (ダイアモン

GPU における並列化では、その特性から、プログラムに条件ド形) であり、α = ∞ ならば、Shape∞ (p, r) は p を中心とす∏dる 1 辺 2r の正方形 di=1 [xi − r, xi + r]⊂=R である。
2. 2 軌跡データと軌跡パターン

分岐が少ない方がより並列化の効率が上がることが知られている。橋口ら [12] は、GPU を用いた論理シミュレーションの並列高速化において、分岐演算数の削減の効果について議論して

文献 [4] にしたがって、軌跡データを定義する。

いる。また、坂本ら [13] は、超解像度画像処理における GPU［定義 1］ (R2 上の) 軌跡テキスト (テキストともいう)、及び上での二分探索を議論しており、この部分の条件分岐が性能低は軌跡パターン (パターンともいう)、長さ n と m の点列

下を引き起こす可能性を指摘している。

2. 準

備

本節では、文献 [2] にしたがって、計算幾何学の基本的な概念と記法を導入する。その後、本稿で考察する多次元点位置検索データ構造を導入する。以下に含まれない、アルゴリズム一般の用語に関しては教科書 [3] を参照されたい。文字列照合に関しては教科書 [9] を参照されたい。

T = t1 , ..., tn ∈ (R2 )∗ ,P = p1 , ..., pm ∈ (R2 )∗である。ここに、各 j >= 1 に対して、tj , pj は点である。

2. 3 近似軌跡パターン照合［定義 2］ 正数を r とする。軌跡パターン P が軌跡テキスト Tに最大距離 r で近似照合するとは、ある位置 1 <=k<= n が存在して、すべての 1 <= m に対して、=i<

2. 1 空間と距離

d∞ (pi , tk+i−1 ) <=rd
はじめに、正整数 d >= 1 を次元と呼ぶ。集合 R = { p =(x1 , . . . , xd ) | xi ∈ R, 1 <= d } を d 次元平面と呼び、各要素=i<dp = (x1 , . . . , xd ) ∈ R を点という。ここで、任意の 1 <=d=i<が成立することをいう。このとき、位置 k を P の T におけるに対して、点 p の第 i 座標を p.i = xi と書く。特に、d = 2 の［定義 3］ 近似軌跡パターン照合問題とは、軌跡テキスト T ととき、第 1 座標と第 2 座標を、x 座標と y 座標と呼ぶ。Rd 上

軌跡パターン P に対し、tj ...tj+m と p1 ...pm がそれぞれの点の距離関数 dist : R2 × R2 → R を、任意の点 x, y, z に対してで近似照合するような T の開始位置 j を求める問題をいう。

次の条件を満たすものとして定義する：

出現位置という。

2. 4 r-配置

任意の α = 1, 2, . . . または ∞ と非負実数 r >= 0 を考える。
入力点集合とは、m >= 1 個の点の集合 P = {p1 . . . pm } であり、P が定める Lα 超球集合とは、P の点を中心とする半径 r

f7

P1

f2

f9

f11

P6
P5

f4
f10

f13
f11

f2

f12

p1

P3

f1

f6

p2

f8
f6

p4

f3

P4

f3
P2

q

f5

f8

f4

p3

f12
f10

p5

f9

f5

f13

f7
f1

図 1 L∞ 距離に対する r-配置 (r-arrangement)A∞ (P, r) と面 fk
図 2 点位置検出問題の例。点 q は面 f4 に出現しているの Lα 超球のなす集合

Shapeα (P, r) := {Shapeα (p1 )、. . . 、Shapeα (pm )}⊂=Rである。ここで、R は R 上の閉領域の族である。
dする。
以後、このデータ構造 PL を、多次元点位置検索データ構
造と呼ぶ。

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

2. 6 GPU と GPU 汎用計算 (GPGPU)GPU(Graphic Processing Unit) とは、複数の Stream Multiprocessor(SM) を備えた、コンピュータ上での画像処理に特化set) を、Coverα (q; P, r) := { 1 <=i<= m | q ∈ Shapeα (pi , r) } (1)
= {1<
=i<
= m | distα (p, q) <
=r}

(2)

したハードウェアである。それぞれの SM には 16〜32 個のSIMD 演算コアが内蔵されており、各演算コアは、SM に搭載された低容量・高速な共有メモリと、GPU ボード上に搭載さと定める。このとき、 全空間 R は、 超球集合 Shapeα (P, r) にれた、大容量・低速のグローバルメモリを用いて並列に計算を含まれる超球の表面によって、いくつかの連続した領域に区切行うことで、高速に画像処理演算を行う。

d

られる。この区切られ方、すなわち、表面の組み合わせを P にこのような多数の演算コアと、メモリ階層を利用した高速な

よる r-配置 (r-arrangement) と呼び、Aα (P, r) で表す。また、

データアクセスによる高い並列演算能力を、画像処理以外の汎用

連続した領域のそれぞれを面 (face) と呼び、この配置における的な計算に用いる動きが増えている。これは GPGPU(General

面すべての集合を F aceα (P, r) で表す。任意の面 f に対して、

Purpose computing on GPU, GPU 汎用計算) と呼ばれ、近年それに含まれる任意の点は、同じ r-被覆集合 Coverα (q; P, r)

注目を集めている技術である [7]。

をもつことに注意されたい。
二次元の場合には、α = ∞ の時、各超球は正方形であり、こ3. 空間パターン索引とビット並列計算を用いた
アルゴリズム

れらの正方形の各面で区切られた多角形領域が面となる。一方、

α = 2 の時、各超球は円盤であり、対応する円で区切られた領

次に 3. 節では、2 次元空間の L∞ 距離に関するトラジェクト

域が面である。一般に、超球の表面は曲線 (または直線) であ

リパターン照合問題について考察する。はじめに、よく知られり、曲線と曲線の交点を頂点 (vertex) と呼び、頂点と頂点を結た文字列照合法である Shift-And 法を取り上げ、この手法に基ぶ面の境界線を辺 (edge) と呼ぶ。以後、α が明確ならば、省略

づく素朴なトラジェクトリパターン照合法を与える。次に、こすることがある。

の手法を改良し、点位置検索データ構造とビット並列計算を組

2. 5 多次元点位置検索データ構造
非負実数 r >= 0 に対して、入力点集合 P = {p1 . . . pm } で定まる r-配置 Aα (P, r) を考える。任意の点 q ∈ Rd の配置 Aα (P, r)

み合わせて、より高速なトラジェクトリパターン照合アルゴリ
ズム FastShift-AndTrajMatch を提案する。

3. 1 Shift-And 法を用いた素朴なアルゴリズム

における出現位置とは、点 q を含む面 f ∈ F aceα (P, r) のこ

はじめに、よく知られた文字列照合法である Shift-And 法をとを表す。ここで、点 q の点位置検出問題を以下のように定義

取り上げ、この手法に基づく素朴なトラジェクトリパターン照する。

合法を与える。

［定義 4］（点 q の点位置検出問題） あ る r-配 置 Aα (P, r) がShift-And アルゴリズムは、Baeza-Yates と Gonnet [1] 及び存 在 す る 時 、任 意 の 点 q に 対 し て 、q の 出 現 位 置 f ∈Wu と Manber [10] によって提案された、ビット並列手法に

F aceα (P, r) s.t. q ∈ f を求めよ。

基づく文字列照合アルゴリズムである。パターンを受理す

図 2 に点 q の点位置検出問題の例を示す。

る非決定性有限オートマトン (N F A) にテキストを 1 文字ず

最後に、本稿で考察する多次元点位置検索データ構造を以下つ入力し、その各状態、及び状態遷移をビットで表すことにのように定義する。

よって O(m + n) 時間で照合を行う。Algorithm 1 に Shift-And［定義 5］（多次元点位置検索データ構造） 入力として、非負実法を示す。この Shift-And 法のトラジェクトリ照合への拡張

数r>= 0、入力点集合 P = {p1 . . . pm } を受け取り、次の点位は、直感的には照合の条件を変えるだけで可能である。Algo-

置検出クエリ P L をサポートするようなデータ構造 PL

rithm 2 に Shift-And 法に基づく近似軌跡パターン照合アルゴ•

P L(q): 点 q の点位置検出問題を解き、その解 f を出力

リズム NaiveShift-AndTrajMatch を示す。

Algorithm 1 Shift-And 法に基づく文字列照合アルゴリズ

る。はじめに、長さ M の軌跡 P が含む点すべてに対して、それ

ム Shift-AndStringMatchらを中心とする半径 r の正方形の全体 D(P ) = {D0 , . . . , DM }を考える。ここに、各 i = 0, . . . , M に対して、Di ∈ D(P ) は

1: Preprocessing
2: P ← p1 p2 ...pm▷ pi は P 中の i 番目の文字

3: T ← t1 t2 ...tn▷ tj は T 中の j 番目の文字

4: for c ∈ Σ do
5:

B[c] ← 0mの補題が知られている [2]。
［補題 1］ P を 2 次元平面上の m 個の点からなる軌跡パター

7: for i ∈ 1...m doンとする。A(P) を P で決まる配置とする。このとき、次が成B[pi ] ← B[pi ] | 0m−i 10i−1

立する。

9: end for（ 1 ） A(P) の頂点数は O(m2 ) 個である。

10: Searching
11: D ← 0m

（ 2 ） A(P) の辺の数は O(m2 ) 個である。

12: for pos ∈ 1...n do
13:D ← ((D << 1) | 0m−1 1) & B[tpos ]14:

if D & 10m−1 =| 0m then15:16:

次に、これらの正方形全体の平面上への L∞ 距離に対する r配置 (r-arrangement)A(P ) を考える (図 1 参照)。この時、次6: end for8:

P の i 番目の点 P (i) を中心とする半径 r の正方形を表す。

（ 3 ） A(P) の面の数は O(m2 ) 個である。
配置 A(P ) の中で、面の全体を F (P ) で表す。これらの面を

f1 , . . . , fK (K >= 1) と番号づける。表記として、面 fk ∈ F (P )report an occurrence at pos − m + 1

が、D(P ) 中の正方形 Di ∈ D(P ) に含まれているときに、

end iffk ∈ Di と書くことにする。

17: end for

アルゴリズムでは、面の集合 F (P ) を、次の演算をサポート

Algorithm 2 Shift-And 法に基づく素朴なトラジェクトリ照1: Preprocessing
2: P ← p1 p2 ...pm▷ pi は P 中の i 番目の文字

3: T ← t1 t2 ...tn▷ tj は T 中の j 番目の文字

4: Searching
5: D ← 0m▷ 状態集合マスク

6: for pos ∈ 1...n do
7:

B ← 0m8:

for i ∈ 1...m do▷ 現在の入力点に対して文字マスクを作る

9:

if d∞ (pi , tpos ) <= r then10:

B ← B | 0m−i 10i−1

11:

する点位置検索 (point location) データ構造 PL に格納する：•

合アルゴリズム NaiveShift-AndTrajMatch▷ 文字マスク▷ i 番目のビットを 1 にする

end if点位置検索演算 P ointLocation(p)：任意の点 p ∈ R2 を入力として、それを含む面 fk ∈ F (P ) を返す。そのような面がなければ、N ull を返す。

•

集合取得演算 BoxList(fk )：任意の面 fk を入力として、それを作る正方形 Di の添字 i の集合をリストとして返す。すなわち、

BoxList(fk ) = { 1 <=i<= m | fk ∈ Di , Di ∈ D(P ) }そのような正方形が 1 つもなければ、∅ を返す。

2 次元平面上の L∞ 距離に関する点位置検索データ構造については、次のようなデータ構造が容易に構築できる。

12:end for13:

D ← ((D << 1) | 0m−1 1) & B［補題 2］（山本他 [11]） 二次元平面上の m 個の正方形の配置

14:if D & 10m−1 =| 0m thenに対して，それを O(⌈ mw ⌉) 語の領域と O(m log m + ⌈ mw ⌉) 構

15:
16:

report an occurrence at pos − m + 1
end if

17: end for2

2

築時間を用いて，点位置検索演算を O(log m + ⌈ m⌉) 時間で実w
現する索引構造が存在する．
証明: 基本的なアイデアは，m 個の正方形の配置 A(P ) を x 軸と y 軸に射影して得られる区間の集合 Segx (P ) と Segy (P )このアルゴリズムの問題点を述べる。Shift-And 法では、ビッに 分 解 し ，各 面 を そ れ ら の 区 間 の 直 積 の 直 和 で 表 す 事 で

トマスク表は |Σ| 個の文字のそれぞれに対して前処理で 1 回だ

あ る ．索 引 構 築 演 算 で は ，全 正 方 形 の 両 端 点 の x 座 標 を ，け作る。しかし、近似軌跡パターン照合においては、前処理のx1 < · · · < xh (h <= 2m) のように O(m log m) 時間でソートし，段階では、可能な入力点の候補は無限個存在するため、あらか

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

じめビットマスクパターンを準備することができない。そのた算する．ここに，x0 = −∞, xh+1 = ∞ である．次に，各区間

め、NaiveShift-AndTrajMatch アルゴリズムでは各入力点に対

sx ∈ Segx (P ) について，sx を含む正方形の集合 SegListx (sx )して毎回ビットマスクパターンを作っている。したがって、実を，長さ m のビットマスク表現 Bx [sx ] ∈ {0, 1}m で記録し，

行時間は O(mn) 時間になる。これは素朴な文字列照合法の実ビットマスク配列 Bx を構築する．同様に y 軸上の区間の集行時間 O(mn) と同じであり、効率が悪い。そこで、次節では、

合を Segy (P ) と，ビットマスク配列 By を構築する．このと前処理でビットマスクパターンを 1 回だけ作る方法を提案する。

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

3. 2 点位置検索データ構造O(⌈ mw ⌉) 語をそれぞれ用いる．点位置検索演算 P ointLocation

本小節では、提案アルゴリズムで用いる、2 次元の L∞ 距離は，入力点の x 座標と y 座標を含む区間 sx と sy を，それぞに対する点位置検索データ構造の実現方法について述べる。これ，Segx (P ) と Segy (P ) 上の二分探索で求める．得られた 2

れは、入力点に対してビットマスクを検索するために用いられつのビットマスク Bx [sx ], By [sy ] のビット AND をとる事によAlgorithm 3 Shift-And 法に基づく提案のトラジェクトリパ

d
入力: パターン点集合 P = {p1 , . . . , pm }⊂=R , 正実数 r > 0.

ターン照合アルゴリズム FastShift-AndTrajMatch

タスク: P の各点を中心とするボックス (正方形) をデータ構造に格納

1: Preprocessing
2: P ← p1 p2 ...pm
3: T ← t1 t2 ...tn▷ pi は P 中の i 番目の文字（ 1 ） ホスト上で、マスタープログラムは、点集合 P をホストから▷ tj は T 中の j 番目の文字

GPU に送信し、GPU のグローバルメモリ上におく。この点集合は、c個ずつのブロック P1 , . . . , P⌈m/c⌉ に分割される。

4: for all fk do
5:

B[fk

] ← 0m▷ 領域 fk に対するビットマスクの初期化

6:list ← BoxList(fk ) ▷ 正方形の添字をリストとして受け取る

7:

while list.length > 0 do
i ← list.top▷ リストの先頭の添字を取り出し、i に代入<<（ 3 ） GPU 上では、1=i= c に対して、各コア Ci はそれぞれ並列に、自分の担当するブロック Pi を、グローバルメモリから自分の属するコアグループの共有メモリに読み込む。
（ 4 ） GPU 上 で 、各 コ ア Ci は ブ ロック Pi に 含 ま れ る 各 点する▷ リストの先頭要素を削除する

9:

erase(list.top)10:

B[fk ] ← B[fk ] | 0m−i 10i−1 ▷ 正方形 Di に対応するビッ

クス (正方形) Bi,1 , . . . , Bi,⌈m/c⌉ を計算する。

▷ 状態集合マスク

end while

12: end for
13: Searching
14: D ← 0m

{pi,1 , . . . , pi,⌈m/c⌉ } に対して、それらを中心とする 1 辺が 2r のボッ（ 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

トを 1 にする
11:（ 2 ） ホスト上で、マスタープログラムは、c 個のコア C1 , . . . , Cc
を GPU 上に起動する。

/* 空間パターン索引の作成 */8:

表 1 並列点位置検索データ構造の前処理アルゴリズム

16:

fk ← P ointLocation(tpos ) ▷ 点 tpos を含む面 fk を見つけるとおくと、昇順にソートした A のキーは x0 < · · · < xn(i,j)+1 とかける。ここで、各ランク 0 <=k<= n(i, j) に対して、辞書 Aj には、キー17:

D ← ((D << 1) | 0m−1 1) & B[fk ]xk に対する値として、m-ビットマスク Mki,j ∈ {0, 1}m を格納する。

18:

if D & 10m−1 =| 0m then

15: for pos ∈ 1...n doi,j

ここに、マスク Mk

report an occurrence at pos − m + 1

19:
20:

▷ 状態遷移を行う

マップである。

end if（ 6 ） GPU 上で、各コアはクエリを待って、待機する。

21: end for

り O(log m + ⌈ m⌉) 時間で実行する．以上より示された．
w

は、値 xk の左に隣接する区間 Ik = [xk−1 , xk ]を第 j 座標の辺が含むような P のボックスの番号の集合を表すビット□

T 、正数 r に対して、アルゴリズム FastShift-AndTrajMatch は

他の 距離 に関 する 点位置検 索デー タ構造については Ma-

O(m log m + m⌈m/w⌉) 前処理時間と O(m⌈m/w⌉) 領域を用い

tousek [8] を参照されたい。

て、O(n(log m + ⌈m/w⌉)) 実行時間で、トラジェクトリパター

3. 3 提案アルゴリズム
次に、3. 1 節の素朴なアルゴリズムを改良し、前節の点位置検索データ構造とビット並列計算を組み合わせて、より高速なトラ
ジェクトリパターン照合アルゴリズム FastShift-AndTrajMatch

ン照合問題を解く。ここで w はワード長である。

4. 多次元点位置検索データ構造のGPU を用いた並列化を提案する。Algorithm 3 に提案アルゴリズムを示す。このア

本節では、2 次元空間の L∞ 距離に関して、GPU 上で高速

ルゴリズムは、前節の素朴なアルゴリズムと同様に、 Shift-Andに動作する多次元点位置検索データ構造のアルゴリズムを考察法に基づくが、マスク表の構築法が大きく異なっている。このする。GPU では、共有メモリを共有する複数のコアがグルー

アルゴリズムは、点位置検索データ構造を用いることで、前処プになっており、各コアは自分のローカルメモリを所有すると理でマスク表を構築することが出来る。

仮定する。さらに、ホストとすべてのコアがグローバルメモリ

提案アルゴリズムは次のように動作する。はじめに前処理時において、与えられたトラジェクトリパターンを、配置 A(P )にアクセスできると仮定する。以降では 4 つの異なるアプロー
チからアルゴリズムを考察する。

から得られる面の集合 F (P ) と、P の各点が構成する正方形を準備として、実数配列 A に対して、1 次元の順序辞書 (pre-辺ラベルとする非決定性有限オートマトン (NFA) に変換する。

decessor dictionary) とは、実数の集合を格納するデータ構造

続いて、F (P ) と NFA から Shift-And 法と同様のマスク表を作であり、挿入 (insert) と、削除 (delete)、前者 (predecessor)、

成する。

後者 (successor) の 4 つのクエリを実現するものである。ここ実行時には、Shit-AndMatch アルゴリズムと同様にして、こに、実数 q に対する前者クエリは、q 以下の最大の A 中の要素のマスク表を用いて入力点を順に読み込みながら、NFA の対応を返し、q に対する後者クエリは、q 以上の最小の A 中の要素する状態遷移を行い、パターン照合を効率よく実行する。このを返す。ただし、これらのクエリが常に回答を返すために、辞時、点位置検索データ構造を用いて、各入力点を含むボックス

書は元のデータに含まれない最小値 min と最大値 max を必ずの集合に対応するマスク表を検索することで、高速化を行う。

含むと仮定する。例えば、C++言語の STL ライブラリの map

このアルゴリズムの計算量は以下の様に与えられる [11]。
［定理 1］ 長さ m の軌跡パターン P と長さ n の軌跡テキスト

は順序辞書である。

表 2 並列点位置検索データ構造の実行時処理アルゴリズム

Algorithm 4 マイクロ線形走査

入力: クエリ点 q ∈ Rd .

1: array A ← ai1 , . . . , aih−1

出力: クエリ点 q を含むボックスの番号の集合を表す m-ビットマス

2: procedure Search(x)ク M ∈ {0, 1}m .

3:（ 1 ） ホスト上で、マスタープログラムはクエリ点 q を受け取ると、

4:

GPU 上の c 個のコア C1 , . . . , Cc に q を、グローバルメモリを用いて

5:

end for

一斉送信する。

6:

return −1（ 2 ） GPU 上で、各 1 <=i<= ⌈m/c⌉ に対して、コア i はグローバ

7: end procedure

ルメモリから q を、共有メモリ経由で局所メモリに並列に受け取る。
（ 3 ） GPU 上で、各コア i は、各次元 1 <=j <= d に対して順に、データを格納した順序辞書 Ai,j 上で、クエリの第 j 座標 q.j をキーとする後者クエリを実行する。これにより、第 j 座標が q.j を含むボッi,j が返される。
クス名 1 <=k<
= m の集合を表す m-ビットマスク M（ 4 ） GPU 上で、各コア i は、すべての次元に対するビットマスクの論理積 M i = M i,1 ∧ · · · ∧ M i,d ∈ {0, 1}m を求めて、共有メモリとグローバルメモリを経由して、並列にホストへ送る。
（ 5 ） ホ ス ト 上 で 、マ ス タ ー プ ロ グ ラ ム は 、GPU 上 の c 個 の
コア C1 , . . . , C⌈m/c⌉ から、グローバルメモリ経由でビットマスク

M 1 , . . . , 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 <= c) をあら=h<
かじめ選び、前処理において、マクロブロック内の c 個の入力点をさらに、h 個ずつからなる ⌈c/h⌉ 個のマイクロブロックに分割する。次に、各 i = 1, . . . , ⌈c/h⌉ に対して、第 i 番目のマ
イクロブロックが含む h 個の入力点 pi1 , . . . , pih によって定めら

for i in 0..h − 1 do
if x exists in A[i] then return i;

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

5:

M = (L + R)/2;

6:

if x < A[M ] then R = M ;7:▷ 不変条件：x ∈ [L, R)

else L = M ;8:

end while

9:

return M ;10: end procedureれる区間 ai1 , . . . , aih−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 分岐レス比較演算x ∈ [L, R) が成立すると仮定すると、７行目から 10 行目の1: procedure Compare(x, y ∈ {0, 1}w : 符号なし w ビット整数)

コード

記法： w ビットレジスタ Z を、MSB と LSB をそれぞれ左と右2:として、Z = bw−1 · · · b0 のように表記する。

of f set = Compare(x, A[M ]);出力： w ビット整数 V . もし x < y ならば V = 0, x >= y ならば

LI = 0 + of f set;V = 1.

RI = 1 + of f set;

X = (X | 10w−1 ); Y = (Y | 0w ); ▷ X と Y の最上位ビット

L = E[LI]; R = E[RI];は、1 と 0.
3:

Z =X −Y▷ Z の最上位ビットは、X < Y のとき 0、4:

return V

を実行後に、もし x < A[M ] ならば L := L と R := M が実行される。逆に x >= A[M ] ならば L := M と R := R が実行されX>= Y のとき 1.
= (Z & 10w−1 ) ≫ (w − 1);▷ 最上位ビットを取り出する。したがって、常に不変条件 x ∈ [L, R) が満たされることがわかる。

5: end procedure

さらに、このアルゴリズムは停止条件についても分岐命令、
つまり元の逐次型二分探索の while 文の条件分岐をもたず、定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

5:

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

6:

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

of f set = Compare(x, A[M ]);8:

LI = 0 + of f set;

9:

RI = 1 + of f set;10:

L = E[LI]; R = E[RI];11:

end for12:

return M ;数回のループで置き換えられている。これは、L + 1 < R が成立する間は、ループの実行により必ず R − (L + 1) が真に減少することと、一度条件 L + 1 = R が成立すると、それ以後の▷ 不変条件：x ∈ [L, R)

ループの実行においても L + 1 = R が成立し続けることが保証されることからわかる。したがって、ループの log c + 1 回の実行後には、必ず L + 1 = R が成立し、L がクエリ値 x の前者▷ Test if x < A[M ]13: end procedure

クエリの解となっている。後者クエリの解は、L = x のとき Lであり、それ以外の時は L + 1 である。
以上より、次の補題が示せる。
［補題 3］ アプローチ B に関する上のアルゴリズムは、マクロブロックにおいて、c 個の w ビットの非負整数値の配列 A[1, c]
上で、w ビットの非負整数であるクエリ値 x に関する二分探索を、分岐演算を用いずに正しく実行する。このとき、一つのの実数を、定められた解像度と原点で変換することで行われる。

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 以上の最小の要素を見つける。

マクロブロックにおいて、使用するメモリ量は 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 <= pc−1= ... <

このアルゴリズムは、ループの不変条件として x ∈ [L, R) が常に成立するように、探索対象の区間を表す変数 L と R を管理する。配列には両端に仮想的な最小値 A[0] = min と最大値

A[c + 1] = max を加えているので、初期値 L = 0, R = c + 1
では仮定は常に成立する。さらに、ループの最初で不変条件のように昇順にソートする。
次に、最小値 pmin = p0 との差分をとって、 次のように各値を差分圧縮し、ベクトルq0 := p0 − pmin , . . . , 1c−1 := pc−1 − pmin .

を、最小値 pmin と一緒に格納する。
次にこれらの差分値を、左から一つのレジスタに入るだけとくことを考えた。各マクロブロック内の計算については、4 つのアプローチを提案した。

り、語長 w(の定数倍) の長さをもつマイクロブロックに詰め込

今後の課題について述べる。3. 節のトラジェクトリパターン

む。各差分値は、これらのうちで最大のビット長 δ に合わせて、照合アルゴリズムについては、その応用を考えることが課題で同じビット数で表す。

ある。また、4. 節の並列データ構造については、時間の都合に実行時には、与えられたクエリ値に対して、各マイクロブ
ロックで独立に 1 次元の点位置検索問題を、次の述べる数値演より GPU 上の実装が出来なかった。これに関して、実装と実験による性能評価が課題である。

算を用いたビット内並列比較を用いて解く。

•

初期値として、一個のマイクロブロックを格納したレジ

スタ 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⌉) 領域を用いて、O(n(log m + ⌈ m⌉)) 実行w
w
時間で動作する。

4. 節では、多次元点位置検索データ構造を GPU 上で並列実装するための効率良い手法について考察した。提案手法は、与えられた入力点集合を定数 c のサイズに分割し、それぞれをマ
クロブロックとして、GPU 上の独立したコア集団で並列に解文

献[1] Ricardo Baeza-Yates and Gaston H. Gonnet. A new approach 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, andMark Overmars. Computational Geometry: Algorithms andApplications. Springer-Verlag, Santa Clara, CA, USA, 3rded. edition, 2008.
[3] Thomas H. Cormen, Cliﬀord Stein, Ronald L. Rivest, andCharles E. Leiserson. Introduction to Algorithms. McGrawHill Higher Education, 2nd edition, 2001.
[4] Thomas Eiter and Heikki Mannila. Computing discreteFréchet distance. Technical Report Technical Report CDTR 94/64, Christian Doppler Laboratory for Expert Systems, TU Vienna, Austria, 1994.
[5] Michael L. Fredman and Dan E. Willard. Surpassing theinformation theoretic bound with fusion trees. Journal ofComputer 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 computing architecture. In Proceedings of the 6th International
Symposium on Memory Management, ISMM ’07, pp. 103–
104, New York, NY, USA, October 2007. ACM.
[8] Jiřı́ Matoušek. Geometric range searching. ACM Computing Surveys, Vol. 26, No. 4, pp. 422–461, December 1994.
[9] Gonzalo Navarro and Mathieu Raﬃnot. 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.