DEIM Forum 2016 G3-4

三次元空間における RMSD 距離に基づく近似点集合マッチングの高速化
佐々木耀一†

渋谷哲朗††

伊藤公人†††

有村博紀†† 北海道大学大学院情報科学研究科 〒 060–0814 札幌市北区北 14 条西 9 丁目†† 東京大学医科学研究所ヒトゲノム解析センター 〒 108–8639 東京都港区白金台 4-6-1††† 北海道大学人獣共通感染症リサーチセンター 〒 001–0020 北海道札幌市北区北 20 条西 10 丁目
E-mail: †{ysasaki,arim}@ist.hokudai.ac.jp, ††tshibuya@hgc.jp, †††itok@czc.hokudai.ac.jp
あらまし 本稿では近似三次元点集合マッチング問題を考察する．この問題は，入力としてデータ点集合 T とパターン点集合 P が与えられて，T の部分点集合で，最小二乗和平均平方根距離 (MinRMSD) のもとで P と類似しているものを出力する問題である．この問題に対して，先に提案した部分点集合列挙と MinRMSD の下界関数に基づく枝刈りを用いたアルゴリズムを改良し，特に小さな半径とサイズをもつパターンに対して，高速に点集合マッチングを行うアルゴリズムを提案する．このアルゴリズムは，前処理として，パターンの半径 ℓ と，パターンの点の数 k と，距離閾値 r からある距離 L を計算し，実行時には各データ点を中心として，半径 L の球内の点のみに探索を制限することでパターン照合を高速化する．元のアルゴリズムはパターンサイズの指数の最悪時間をもつが，このアルゴリズム
は ℓ と, k, r が定数のとき，空間ポアソン過程で生成されたデータ点集合に対して，O(n logk n) 時間で点集合マッチを行う．さらに，実際的な高速化技法を導入する．
キーワード

パターン照合，分枝限定法，列挙アルゴリズム，RMSD, 最小二乗和平均平方根距離，探索空間制約
RMSD スコアの下界関数を用いた分枝限定法を用いた点集合

1. は じ め に
1. 1 背マッチングアルゴリズムを考案した．2 節では，この BasicMatch

景アルゴリズムを説明する．

高分子やタンパク質の三次元構造においての検索技術は，分本稿では，特に小さな半径とサイズをもつパターンに対して，子生物学において近年，重要な役割を担っている．例えば，あ高速に点集合マッチングを行うアルゴリズム LPSM を 3 節でるタンパク質の構造が与えられ，その性質が知られていないと提案する．このアルゴリズムは，前処理として，パターンの半きに，構造が似ているものは性質も似ているという理論から，

径 ℓ と，パターンの点の数 k と，距離閾値 r からある距離 L を既に構造と性質が知られているタンパク質に対して，構造的に計算し，実行時には各データ点を中心として，半径 L の球内の類似しているかを調べることにより，性質を予想することがで点のみに探索を制限することでパターン照合を高速化する．元きる [1]．また，近年の急速なデータベースの増大に伴い，よりのアルゴリズムはパターンサイズの指数の最悪時間をもつが，早い時間での検索技術が必要とされている．

このアルゴリズムは ℓ, k, r が定数のとき，空間ポアソン過程で1. 2 本研究の目的生成されたデータ点集合に対して，O(n logk n) 時間で点集合

本論文では，データを点集合の入力とし，与えられたパターマッチングを行う．

ンの点集合と構造が似ている部分点集合をデータの点集合の1. 3 関 連 研 究中から見つけ出すという近似点集合マッチング問題を考えてい本論文の関連研究として，Schwartz と Sharir [3] はデータく．この問題は，データとして n 個の点を要素として持つ点集ベースとパターンを共に点列の入力とし，RMSD 指標による合 T と，パターンとして k 個の点からなる点集合 P ，非負実マッチングを考察し，データベースサイズを N としたときに数 r > 0 を入力とし，パターンに対する類似度が r 以下となるO(N log N ) 時間で実行するアルゴリズムを提案している．ま

データ中のサイズ k の部分点集合を抽出する問題である．

た，Shibuya [2] は，同じマッチング問題を，データベースサイナイーブなアルゴリズムは，データベースの点集合から可能ズに対して線形時間で実行するアルゴリズムを提案している．

な k 個の点を抽出し，その後，パターンとの類似度を求め，閾

その他に，Akutsu [4] らは，ハッシュを用いた手法を提案して

値 r より小さければ解となる類似部分点集合として，出力すおり，O(N ) 期待計算時間を実現している．また，Shibuya [1]る．この動作を繰り返すことにより，データベースの全ての部は，幾何接尾辞木 (geometric suﬃx tree) と呼ばれる，三次元分点集合に対して近似マッチングを見つけることができる．計構造に対する索引データ構造を示している．これは，文字列に算量は k · n Pk =k · n · (n − 1) · · · (n − (k + 1)) より O(knk ) 時

対する索引構造である接尾辞木 (suﬃx tree) [6], [7] を基に，タk間となり，データの内容に関わらず，常に O(kn ) 時間を要す

ンパク質のような複雑な三次元構造に対して適応させたデータる．作業領域は O(k) 領域である。

構造である．このデータ構造は，RMSD 指標による効率よいこの問題を解決するために，文献 [9] [10] で，我々は最小近似部分構造マッチングを実現し，データベース中に存在するタンパク質構造同士の頻出する類似部分構造の発見にも用いることができる．

移動するような v が RM SDR,v (S, T ) を最小化することが知ら

また，三次元点集合マッチング問題に対して，分枝限定法を用いた実際的に高速なアルゴリズムが提案されている．[9] [10]2. 準回転行列 R が固定されていれば，R · T の重心を S の中心にれている [2]．つまり R̂ が与えられれば v̂ は O(n) 時間で計算可能である．
二つの点列 S, T の重心を原点に合わせるように移動した後，備

RM SDR,v (S, T ) を最小化する回転行列 R を考える．これは，vu nu∑RM SDR (S, T ) = t|si − (R · ti )|2

2. 1 基本的定義
本稿では，三次元空間 R を考える．長さ n の点列（point′

3sequence）S は，S = (s1 , . . . , sn ) ∈ (R3 )n である．ここに，
各 i = 1, . . . , n に対して，si = (xi , yi , zi ) ∈ R3 は i 番目の点の三次元座標を表すベクトルである．ここに，|S| = n
で S の 長 さを 表 す．任 意 の 1 <
= i <
= j <
= n に対して，P [i . . . j] = (si , si+1 , . . . , sj ) で S の添え字 i から j までの部分点列を表す．点列 S の長さ 1 <=i<= n の前綴り（prefix）とは，最初の i 個からなる列 S[1..i] である．
三 次 元 の 回 転 行 列 R に よって 点 列 S を 回 転 し た 構 造 は

R · S = (Rs1 , . . . , Rsn ) で与えられる．回転行列 R と移動ベク
トル v は三次元空間の点に対する変換 f = fR,v (S) = R · S + vを与える．ベクトル v の転置ベクトルを v t で表し，行列 A の転置行列を At で表す．trace(A) は行列 A のトレースを与え，|v| はベクトル v のノルムを与える．

を最小化する R̂ を求めることと同じである．ここで，J =∑nt
i=1 si · ti と置くと，これは明らかに O(n) 時間で計算可能で∑′ttある．従って RM SD (S, T ) は ni=1 (si si +ti ti )−2·trace(R·J)∑と表すことが出来る．また，UV を J の特異値分解 (SVD)とすると R = V U T のとき最大となる [5]．三次元座標においては J が 3 × 3 行列であり，3 × 3 に対する SVD は O(33 ) 時間で計算可能であることから R̂ は J から定数時間で求められる．
以上より，v̂ ，R̂ がそれぞれ O(n) 時間で求められることから，M inRM SD は定義式より O(n) 時間で求められることは明らかである．

2. 4 二つの点集合間の最小平方平均二乗和距離
本小節では，点列を拡張し，本稿で問題としている，二つの2. 2 二つの点列間の最小平方平均二乗和距離
初めに，Shibuya [1] で三次元点列間の類似度として用いられている二つの点列間の最小平方平均二乗和距離を定義する．
三次元座標で表された長さ n の二つの点列を S = (s1 , . . . , sn )と T = (t1 , . . . , tn ) とおく．任意の回転行列 R と移動ベクトルv が与えられたとき，S と T の間の R と v で決まる変換のもとでの平方平均二乗和距離を式
vunu1 ∑RM SDR,v (S, T ) = t|si − (R · ti + v)|2n i=1

(3)

i=1

点集合間の類似度として用いられる最小平方平均二乗和距離を定義する．P = P [1..n] と Q = Q[1..n] を，要素数 n の三次元空間の点集合とする．ただし，これらの点集合は，点列として表されているものとし，要素である点の順番は任意に固定するとする．
順列 π : {1, . . . , n} → {1, . . . , n} に対して，π による長さ n の点列 P = P [1..n] の並べ替えを π(P ) = (P [π(1)], . . . , P [π(n)])(1)と定める．
点列として表されたサイズ n の集合 P と Q 間の最小平方平均二乗和距離を，すべての順列 π と，回転行列 R，移動ベクトによって定義する．

ル v に関する最小の RMSD 値

このとき，二つの長さ n の点列 S と T 間の最小平方平均二M inRM SD(P, Q) = min min RM SDR,v (π(P ), Q)乗和距離（最小 RMSD スコアともいう）を，

πR,v(4)

M inRM SD(S, T )と定義する．すなわち，M inRM SD(P, Q) は，全ての並べ替= min RM SDR,v (S, T )えに関して点列 P と Q の最小 RMSD スコアをとった最小値

vunu1 ∑= min t|si − (R · ti + v)|2R,vn i=1である．

R,v(2)と定義する．上式で，R と v は全ての回転行列と移動ベク

2. 5 幾何学的距離に関する近似点集合マッチング問題
本論文で考察する最小 RMSD スコアに関する点集合マッチ
ング問題を次のように定義する．三次元空間を考える．

トルの上を動く変数である．すなわち，最小 RMSD スコア

データ点集合とは，n 個の点からなる点集合 T = {t1 , . . . , tn }M inRM SD(S, T ) は，全ての回転行列 R と移動ベクトル v にで あ る ．パ タ ー ン 点 集 合と は ，k 個 の 点 か ら な る 点 集 合2. 3 点列間の最小 RMSD スコアの計算方法

P = {p1 , . . . , pk } である．T と P の要素を，それぞれ，データ点およびパターン点と呼ぶ．ただし，以下では常に k <=nで

ここで，Shibuya [1] に基づいて，二つの点列間の最小 RMSDあり，データ点集合とパターン点集合は，点列として表される

対する RM SDR,v (S, T ) の値の最小値である．

スコアの線形時間アルゴリズムを与える．点列 S と T に対しものと仮定する．て，RM SDR,v (S, T ) を最小とする回転行列，移動ベクトルを非負実数 r > 0 に対して，パターン点集合 P がデータ点集それぞれ R̂，v̂ とする．最小 RMSD スコアは，以下に示す方

合 T に最小 RMSD スコア r でマッチするとは，k 個のデー法で，O(n) 時間で計算可能である．

タ点からなるある点集合 Q⊂=T が存在して，

M inRM SD(P, Q) <=r

Algorithm 1 下界関数による枝刈りを用いた点集合マッチンが成立することである．このとき，T の部分点集合 Q を，T における P の出現位置という．

グの基本アルゴリズム BasicMatch
1: procedure BasicMatch(P, T, r)

入力： パターン P [1..k]，テキスト T [1..k]，正実数 r.

［定義 2.1］ 最小 RMSD スコアに関する点集合マッチング探

出力： T 中の最小 RMSD スコア r に関する P の全ての異なる出索問題は，データ点集合 T と，パターン点集合 P ，非負実数現位置．

BasicFind(ε, 0, P, T, |Q|, |T |) を呼び出す;

r > 0 を受け取り，T 上の最小 RMSD スコア r でのパターン

2:

P の出現位置 Q⊂=T を一つ見つける問題である．

3: end procedure

上記の問題の変種として，P が T に最小 RMSD スコア r でマッチするかどうかを YES または NO で答えるものを点集合マッチングの決定問題と呼ぶ．また，全ての出現位置 Q⊂=T を

4:
5: procedure BasicFind(Q[1..i], i, P, T, k, n)
6:
7:

もれなく，かつ重複なしに見つけるものを点集合マッチングの8:列挙問題と呼ぶ．ただし，点列として異なる出現位置 Q は互い9:に異なるとみなす．

10:

if i = k then

マッチング位置の集合 Q を出力する;

end ifelse
for j = 1, . . . , n do

本稿では，とくに断らなければ，点集合マッチングの列挙版

11:

問題を考察する．決定版と探索版の問題は，列挙版のアルゴリ

12:

x = T [j];13:

if x ∈ Q then14:

continue

ズムで解くことができる．

2. 6 基本アルゴリズム▷ すべての点を選び終わった

if M inRM SD(P [1..k], Q[1..k]) <= r then

15:end if

本小節では，我々が考案した三次元点集合マッチングのアル16:

R := Q の末尾に要素 x を加えて得られる点列;

ゴリズム [9] [10] を紹介する．このアルゴリズムは，最小 RMSD

17:

if LB M inRM SDi+1,k (R, P [1..i + 1]) > r then

スコアの下界関数を用いた探索の枝刈りによる効率化を行う．

18:

まず初めに，素朴なアルゴリズムとして，再帰的な計算によ19:

り解となる候補集合を探索し，探索木の葉，すなわち，候補集合のサイズがパターンサイズと同じ k に達したときに初めて，最小 RMSD スコアによる解の検査を行う手順が考えられる．

continueend ifBasicFind(R, i + 1, P, T, k, n) を再帰的に呼び出す;

20:
21:
22:

end forend if

23: end procedure

しかし，探索のすべての枝を葉までたどるので効率が悪い．
これを改良するために，前綴りに関する最小 RMSD スコアの下界関数（lower-bound function of minimum RMSD forprefices）を導入し，これを用いた分枝限定法（branch-and-可能性をもっている点で実際には利点がある．

3. 提 案 手 法bound method） を取り入れている．
［定義 2.2］（MinRMSD の下界関数） 任意の正整数 1 <=i<=kに対して，下界関数を次のように定める：

本節では，探索空間制約を用いた最小 RMSD スコアに関する３次元点集合マッチの高速化を行う．このアルゴリズムはパターン点集合のサイズが小さく，空間ポアソン過程で生成されLB M inRM SDi,k (P [i], Q[i])( )1/2= ki
M inRM SD(P [i], Q[i]).

たデータ点集合に対し高速に計算できる．

(5)

3. 1 提案アルゴリズムの概要
Algorithm 2 に，提案アルゴリズム LPSM (Locality ofAlgorithm 1 に，下界関数による分枝限定を用いた点集合Point Set Matching) を示す．

マッチングアルゴリズムを示す．このアルゴリズムは，再帰にLPSM アルゴリズムは次のように動く．初めに，データ点集基づいて，すべての候補集合を列挙し，最小 RMSD スコアが合 T を全て内部に含むような立方体を求める．次に，図 1 のよ閾値 r 以下の解を正しく出力する．17 行目から 19 行目で，下うにその立方体を M × M × M のセルに分割し，各セルにデー界関数を用いた枝刈りを行う．この枝刈りは健全であり，いかタ点を一つずつ格納する．その後，パターン点集合 P より，十なる正解を見落とすこともないことが言える．このアルゴリズ分大きな正実数 L を求める．L の求め方は後の 3. 2 節で説明

ムについて次が示されている．

する．各データ点に対し，図 2 に示すように，その点を中心に［定理 2.1］ 点集合マッチングアルゴリズム Algorithm 1 は，辺の長さが 2L の立方体を考え，この立方体を包含する複数個与えられたデータ点集合 T とパターン点集合 P に対して，最のセルを求める．最後に，その複数個のセルに格納されている

小 RMSD スコアが閾値 r 以下となる P のすべての出現位置を点のみに対して，つまり局所的な範囲に制限してパターン点集正しく見つける．

合 P との点集合マッチングを行う．マッチングには，任意の既アルゴリズムの使用領域は O(k) 領域である．計算時間に関存の点集合マッチングアルゴリズムを用いることができる．

しては，総計算時間が O(knk ) 時間であること以外はわからなこのアルゴリズムの健全性は明らかである．もしアルゴリズ

い．ただし，素朴なアルゴリズムが入力に無関係に常に Θ(knk )ムが最小 RMSD スコア r 以下のマッチングを見つければ，その時間を要するのに対して，枝刈りアルゴリズムは早期終了の

Algorithm 2 部分空間制限を用いた点集合マッチングのアルzL

ゴリズム LPSM (Locality of Point Set Match)L

1: procedure LPSM(P, T, r)

入力： n 個のデータ点集合 T と，k 個のパターン点集合 P ，正実数 r > 0．
出力： T 中の最小 RMSD スコア r に関する P の全ての異なる出L

現位置．
2:

データ点集合 T を全て内部に含むような直方体 A を求める．

3:

M × M × M のセルに分割し，各セルにデータ点を格納する．

4:

パ タ ー ン 点 集 合 P に 対 し て パ ラ メ ー タ k, ℓ, r か ら L =y

L

L(k, ℓ, r) > 0 を求める．
5:6:(xi  ,  yi ,  zi )for i = 1, . . . , n do
データ点の部分集合として，選択したデータ点 ti を中心とLし，辺の長さが 2L の立方体 C(ti , 2L) を求める．
7:

C(ti , 2L) を含むセルを，分割された直方体 A の中から取り出す．

8:L

図 2 ある点を始点にした際に探索しなければいけない最小空間取り出したセル内の点のみに対して，任意の RMSD に関する点集合マッチアルゴリズムを用いて，データ点集合 Tc におけるパターン P の RMSD に関する点集合マッチの出現 Q ⊂ Tc を見つけ出し，全て出力する．

まれる点の数を N (B) で与える．
以下で，前小節で述べた L について定理を述べる．

end for［定理 3.1］（マッチングにおける局所性） 中心が c ∈ Rd ，半10: end procedure

径が ℓ > 0，サイズ k の任意のパターン点集合 P と，実数9:r > 0 を与える．また，回転 R と平行移動 v を行う幾何変換 fを与える．このとき，データ点集合 T 内の，任意のサイズ k のz<部分点集合 Q⊂=T に対して，もし M inRM SD(f (P ), Q) = r が成立するならば，Q の全ての点は，あるデータ点 c′ = f (c) ∈ Tを中心とする，半径L = L(k, r, ℓ) := 2(k1/2 r + ℓ)

M(7)の超球に含まれる．
（証明） P = (pi )ki=1 ，Q = (qi )ki=1 とおく. 任意の i にy

対 し て ，||f (c) − qi || <= L が成り立つことを示せば良い．

L2 距 離 に お け る 三 角 不 等 式 を 考 え る と ，||f (c) − qi || <=||f (c) − f (pi )|| + ||f (pi ) − qi || (*) が言える．初めに，仮定
より ||c − pi || <= radius(P ) = ℓ なので，||f (c) − f (pi )|| <=ℓがMM示せる．次に仮定より r >= RM SD(f (P ), Q) = L2 (f (P ) − Q)

xである．k 個の要素をもつベクトル X に対して L2 (X) >=図 1 M × M × M のセルに分割したデータ点集合を包含する直方体

れは正しい解である．一方で，完全性，すなわち解の網羅性は，正実数 L の選び方に依存する．探索範囲を狭める際に，正解を含むために十分大きな L を取ることが重要であることがわかる． 3. 2 節で L の設定法を説明する．
まず初めに，任意のデータ点 c ∈ Rd を中心に持つ，任意の半径 ℓ > 0 の超球を}
Bc,ℓ = q ∈ Rd | ||q − c|| <=ℓ1/2に掛けると, ||f (pi ) − qi || <= k r が導かれる．上記の議論と1/2(*) 式より，||f (c) − qi || <= ℓ + k r が言える．

同様に，f (c) のもっとも近いデータ点 tc ∈ T に対して1/2||tc − f (c)|| <= ε となることが示せる．ただし，ε := ℓ + k rである．

3. 2 点集合マッチングにおける局所性{

1/2(1/k)1/2 L∞ (X) なので L2 (f (P ) − Q) >= (1/k) L∞ (f (P ) −
1/2
1/2>
Q) >
= (1/k) ||f (pi ) − qi ||(∗∗). k= 0 より (**) 式の両辺

従って，三角不等式を再度適用することにより，||tc −f (pi )|| <=1/2||tc − f (c)|| + ||f (c) − qi || <= ε + (ℓ + k r) = L を得ることができる．
（証明終わり）(6)

3. 3 空間ポアソン過程における解析
空間ポアソン過程 [8] とは密度 λ の均一分布モデルである．と定める．
ここで，点集合 P の半径を radius(P ) で与え，P 内の全ての点を包含する球の半径の最小のものとする．また，球 B に含空間ポアソン過程は以下の性質を満たす．

(S1) 任意の球 B に対して，次が成り立つ．P r(N (B) = k) =((λ|B|)k /k!)e−λ|B|(S2) 独立な A1 , . . . , Am に対して, N (A1 ), . . . , N (Am ) も独立である．
理論的解析については，次の定理を示すことができる．証明は [10] を参照されたい．
［定理 3.2］ データ点集合が密度 λ の空間ポアソン過程で生成されるとき，パラメータ θ = (k, ℓ, r): パターンのサイズ k > 0，半径 ℓ > 0，閾値 r > 0 を定数として固定する．このとき，任意の δ > 0 に対して以下が言える:
任意の半径かつ n 点のデータ点集合 T と，ℓ の半径かつ k 点のパターン点集合に対して，t ∈ T を中心に，半径 L = L(k, r, ℓ) 内に含まれる点を集めた部分点集
合 Tt に対してマッチングを行うことにより，1 − δ の確率で O(n logk n) 時間， O(log n) 領域で行うことができる．

4. 結論

本稿では，最小 RMSD スコアに関する３次元空間での点集合マッチングを考察し，空間制約を用いた高速化アルゴリズムについて議論した．
今後の課題として，提案アルゴリズムの実装を行い具体的なデータに対して実験を行いたい．また応用として，本稿の手法に基づいた，分子データベースや，3D モデルデータベース，時空間系列データにおける近似検索への適用は興味深い課題である．
文

献[1] Tetsuo Shibuya, ”Geometric Suﬃx Tree : Indexing Protein
3-D Structures,” Jounal of the ACM, Vol.57, No.3, Article
15, 2010.
[2] Tetsuo Shibuya, ”Searching Protein 3-D Structures in Linear Time,” RECOMB 2009, LNCS 5541, 115, 2009.
[3] Jacob T. Schwartz, Micha Sharir, ”Identification of partiallyobscured objects in two and three dimensions by matching
noisy characteristic curves,” Intl. J. of Robotics Res. 6, 2944,
1987.
[4] Tatsuya Akutsu, Kentaro Onizuka, Masato Ishikawa, ”NewHashing Techniques and Their Application to a Protein
Strucure Database System,” Proc. Hawaii Int. Conf. System Sciences 5, 197-206, 1995.
[5] Gilbert Strang, ”Introduction to Linear Algebra, 3rd Edition,” Wellesley-Cambridge Press, 2003.
[6] Dan Gusfield, ”Algorithms on Strings, Trees and Sequences:
Computer Science and Computational Biology,” Cambridge
University Press, 1997.
[7] Esko Ukkonen, ”On-line construction of suﬃx-trees,” Algorithmica, Vol. 13, No. 3, 249-260, 1995.
[8] Mark Pinsky and Samuel Karlin, ”An introduction tostochastic modeling,” Academic press, 2010.
[9] 佐々木 耀一，渋谷 哲朗，伊藤 公人，有村 博紀，”三次元空間における効率良い近似点集合マッチングと分子パターン照合への応用,” 第 42 回バイオ情報学研究会, 2015.
[10] Yoichi Sasaki, Tetsuo Shibuya, Kimihito Ito, Hiroki
Arimura, ”Eﬃcient Approximate 3-Dimensional PointSet Matching Using Root-Mean-Square Deviation Score,”
SISAP 2015, 191-203, 2015.