160414駒場

数理工学のすすめ
http://researchmap.jp/sada/
東京大学大学院情報理工学系研究科
数理情報学専攻
定兼邦彦
2016年4月14日
簡潔データ構造
• データを極限（情報理論的下界）まで圧縮
• 無視できるサイズの索引を用い
• 高速な問い合わせ等を実現
2
木構造
• 情報検索等で用いられる基本的なデータ構造
• データを再帰的に分割することで効率よく検索
2
13
1
3
8
7
5
10
6
15
9
3
• ある値 (7) より大きい数と小さい数に分ける
2
13
1
3
5 6
7
8
9
15
10
4
• それぞれをさらに分ける
2
1
3 5
6
7
9
8
13
10
15
5
• 分割の仕方を木で表現する
• 根の値 (7) より小さい値は左の部分木，
大きい値は右の部分木にある
• ある値を探すときは，どちらかの部分木だけ
探せばいい→高速化
根 (root)
7
3
10
2
1
6
5
8
15
9
13
6
かな漢字変換辞書
• 「読み」の辞書順に従って単語を木に格納する
• 同じ文字が何度も出てくるので効率化したい
辞書
あさ  朝
あさひ  朝日
あし  足
あしか  アシカ
あしかせ  足枷
あした 明日
あせ  汗
あそ  阿蘇
いえ  家
いか  以下
いし  石
うし  牛
うしろ  後ろ
うそ  嘘
7
• 接頭辞 (prefix) が同じ部分は1つにまとめる
⇒トライ (trie) と呼ばれる木構造
辞書
あ
い
さ
しせ
足汗阿蘇
朝
ひ
そ
か
た
朝日アシカ明日
せ
え
う
し
かし
家以下石
牛
ろ
後ろ
そ
嘘
あさ  朝
あさひ  朝日
あし  足
あしか  アシカ
あしかせ  足枷
あした 明日
あせ  汗
あそ  阿蘇
いえ  家
いか  以下
いし  石
うし  牛
うしろ  後ろ
うそ  嘘
足枷
8
辞書の探索で必要な木の演算
• ある点 v から，文字 c のラベルがついている
v
枝をたどった点 w を求める Child(v, c)
あ
い
• 点 v の数字を求める Preorder(v)
w
1
あ
2
さ
3
ひ
4
せ
10
9
か
た
6
8
せ
7
11
そ
し
5
う
い
え
12
15
し
かし
13
14
16
ろ
17
そ
18
3朝
4  朝日
5足
6  アシカ
7  足枷
8  明日
9汗
10  阿蘇
12  家
13  以下
14  石
16  牛
17  後ろ
9
18  嘘
う
Preorderの求め方
• 木構造を括弧列 P で表現する
• 木の点 v は P 中の開き括弧の位置で表される
• v のPreorderは，P の先頭から v までの開き括弧
の数と等しい
⇒ P でのrankで求まる
1
2
3
4
9
5
6
11
10
8
12
13
15
14
16
18
17
7
1 234
5 67
8
12
16
9 10 11 13 14 15 17 18
P (((())((())())()())(()()())((())()))
10
Childの求め方
• Child(v, c) を求めるために，次の演算を使う
– First(v): v の一番左の子 w
– Next(w): w の次の子 x
• FirstとNextを使って，v の子 w, x, y を
順番に見ていく
• 枝のラベルが c である子が出てきたら答え
v
あ
w
い
x
う
y
11
Firstの求め方
• v を表す括弧 ( … ) の間に，v の子 w, x, y が
書かれている
• v の一番左の子は v の括弧の直後にある
• つまり First(v) = v + 1
v
w
x
y
vw x y
P (()()())
12
Nextの求め方
• w の次の子 x は，w を表す括弧 ( … ) の右にある
• findclose(P, w) を，位置 w の開き括弧に対応する
閉じ括弧の位置とする
• Next(w) = findclose(P, w) + 1
1
w
3
9
5
4
6
w
2
x
2
10
8
5 67
13
15
14
16
y
18
17
7
1 34
12
11
8
y
x
11
15
9 10 1213 14 1617 18
P (((())((())())()())(()()())((())()))
13
括弧列表現
•
•
•
•
BP (Balanced Parentheses) 表現とも言う
各節点を対応する括弧のペアで表現
n 節点で 2n bits
1
サイズの下限と一致
2
3
6
4
5
7
8
1
2
6
3
BP
P
4
5
7
8
((()()())(()()))
14
括弧列での基本操作
• ノードは( の位置で表す
• findclose(P,i): P[i] の( と対応する )の位置を返す
• enclose(P,i): P[i] の( を囲む括弧の位置を返す
enclose
1
3
2
4
5 6
11
8
7
9
findclose
1 2 3 4 5 6
7
8 9 10
11
P(()((()())())(()())())
10
15
木の巡回操作
•
•
•
•
parent(v) = enclose(P,v)
firstchild(v) = v + 1
sibling(v) = findclose(P,v) + 1
lastchild(v) = findopen(P, findclose(P,v)1)
1
enclose
3
2
4
5 6
11
8
7
9
findclose
1 2 3 4 5 6
7
8 9 10
11
(()((()())())(()())())
10
16
子孫の数 (部分木の大きさ)
• v を根とする部分木の大きさは
subtreesize(v) = (findclose(P,v)v+1)/2
• degree (子の数) は findclose の繰り返しで求まるが
子の数に比例する時間がかかる
1
3
2
4
5 6
11
8
7
9
1 2 3 4 5 6
7
8 9 10
11
P (()((()())())(()())())
10
17
葉に関する操作
• 葉はBP中の () で表される
• i 番目の葉の位置 = select()(P, i)
• 部分木中の葉の数，部分木中の最左葉，最右葉
も求まる
1
3
2
4
5 6
11
8
7
9
10
1 2 3 4 5 6
7
8 9 10
11
P(()((()())())(()())())
3 を根とする部分木
18
ノードの深さ
• 超過配列 E[i] = rank((P,i)  rank)(P,i) とすると
depth(v) = E[v]
• E は格納せず P と rank 計算用索引を格納
1
3
2
4
11
8
7
9
10
1 2 3 4 5 6
7
8 9 10
11
P (()((()())())(()())())
E 1212343432321232321210
5 6
19
findcloseの求め方
• findclose(P, v) = w とする
• v から括弧列を右に見ていき，E[w] = E[v]1 と
なる最初の w が答え
• 1つずつ見ていくと遅い
1
3
2
4
11
8
7
9
10
1 2 3 4 5 6
7
8 9 10
11
P (()((()())())(()())())
E 1212343432321232321210
v
w
5 6
20
区間最大最小木
• 超過配列 E を長さ s のブロックに分割
• 木の葉は1つのブロックに対応し，ブロック中の
最小値と最大値を格納
• 内部節点は l 個の子を持ち，子の最大/最小値を
格納
0/4
1/4
1/3
0/2
m/M 1/2 2/4 3/4 2/3 1/3 2/3 1/2 0/0
E 1212343432321232321210
(()((()())())(()())())
s=l=3
21
区間最大最小木の性質
• 各節点は配列のある区間に対応
• 配列の任意の区間は，節点に対応する O(lh) 個の
区間と2つの葉の部分区間の和集合で表される
(h は木の高さ)
0/4
1/4
1/3
0/2
m/M 1/2 2/4 3/4 2/3 1/3 2/3 1/2 0/0
E 1212343432321232321210
(()((()())())(()())())
s=l=3
22
超過配列の性質
• 全ての i に対し，E[i+1] = E[i]1 またはE[i]+1
• 区間 E[u,v] の最小値が a, 最大値が b のとき，
区間内には a  e  b の全ての整数 e が存在し，
それ以外は存在しない (中間値の定理)
• 長さ l の区間の最大値と最小値の差は l1 以下
–
⇒ 値を少ないビットで表現可能
1
3
2
4
5 6
11
8
7
9
10
1 2 3 4 5 6
7
8 9 10
11
P (()((()())())(()())())
E 1212343432321232321210
23
findclose(P,i)の計算
• 区間 E[i+1,N1] を前述のように分割 (N: 配列長)
• 分割した区間を左から順に見て行き，E[i]1 を
含む区間を探す
• O(lh+s) 時間
0/4
1/4
1/3
0/2
m/M 1/2 2/4 3/4 2/3 1/3 2/3 1/2 0/0
E 1212343432321232321210
(()((()())())(()())())
24
定理: n 節点の順序木は 2n + o(n) ビットで表現でき，
各種基本操作は O(1) 時間で実行できる．
基本操作: i-th child, degree, parent, level-ancestor,
leftmost/rightmost leaf, preorder, postorder, inorder,
depth, height, lca, subtree size, …
実装も容易で，速度も速い
2.3n bits
簡潔ではないデータ構造のサイズの 1/50 以下
25
文字列検索
• 情報検索で必須の処理
– サーチエンジン
– ゲノム情報処理
• データ量が莫大
– Web: >20億ページ, 数テラバイト
– DNA配列: ヒト=28億文字, 総計>150億文字
• 索引を構築し，高速検索
26
基本問題
• 存在問い合わせ see: No sea: Yes
• 頻度問い合わせ sea: 2回 sh: 3回
• 列挙問い合わせ sh: 1, 14, 31
1
14
31
She sells seashells on the seashore.
27
検索方法
文字列 T (長さ n) からパタン P (長さ m)を検索
• データ量小⇒逐次検索 (KMP法, BM法)
– O(n+m) 時間
• データ量大⇒索引を使った検索
– 索引の作成 O(n) 時間
– 検索 O(m+output) 時間 (接尾辞木)
28
代表的な索引
• 転置ファイル
• 接尾辞木
• 接尾辞配列
29
転置ファイル
• 文字列を単語(形態素)に分解
• 単語ごとに出現位置(出現文書)を列挙
• 出現回数も記憶
1
4
8
12
15
19
いるかいないかいないかいるかいるいるいるか
T = いるか|いないか|いないか|いるか|いるいる|いるか
いるか: 3回 (1, 12, 19)
いないか: 2回 (4, 8)
いるいる: 1回 (15)
30
索引の完備化
• 任意の文字列を検索したい
• 部分文字列の数 = nC2 = O(n2)
• 全ての部分文字列を索引に格納
⇒索引サイズ: O(n2) 語
31
接尾辞 (suffix)
• 文字列 T の先頭の何文字を除いたもの (n 種類)
• T の任意の部分文字列は，ある接尾辞の接頭辞
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
いるかいないかいないかいるかいるいるいるか
るかいないかいないかいるかいるいるいるか
かいないかいないかいるかいるいるいるか
いないかいないかいるかいるいるいるか
ないかいないかいるかいるいるいるか
いかいないかいるかいるいるいるか
かいないかいるかいるいるいるか
いないかいるかいるいるいるか
ないかいるかいるいるいるか
いかいるかいるいるいるか
かいるかいるいるいるか
いるかいるいるいるか
るかいるいるいるか
かいるいるいるか
いるいるいるか
るいるいるか
いるいるか
るいるか
いるか
るか
か
=T
32
接尾辞木 [Weiner 73]
• 全ての接尾辞を格納した木
か
な
る
い
い
か
か
な
い
な
る
る
い
か
い
い
6
10
4
8
15
17
19
1
12
21
3
7
14
11
5
9
16
18
20
2
13
いかいないかいるかいるいるいるか
いかいるかいるいるいるか
いないかいないかいるかいるいるいるか
いないかいるかいるいるいるか
いるいるいるか
いるいるか
いるか
いるかいないかいないかいるかいるいるいるか
いるかいるいるいるか
か
かいないかいないかいるかいるいるいるか
かいないかいるかいるいるいるか
かいるいるいるか
かいるかいるいるいるか
ないかいないかいるかいるいるいるか
ないかいるかいるいるいるか
るいるいるか
るいるか
るか
るかいないかいないかいるかいるいるいるか
33
るかいるいるいるか
接尾辞配列 [Manber, Myers 93]
• 接尾辞のポインタを辞書順にソートした配列
SA
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
いるかいないかいないかいるかいるいるいるか
るかいないかいないかいるかいるいるいるか
かいないかいないかいるかいるいるいるか
いないかいないかいるかいるいるいるか
ないかいないかいるかいるいるいるか
いかいないかいるかいるいるいるか
かいないかいるかいるいるいるか
いないかいるかいるいるいるか
ないかいるかいるいるいるか
いかいるかいるいるいるか
かいるかいるいるいるか
いるかいるいるいるか
るかいるいるいるか
かいるいるいるか
いるいるいるか
るいるいるか
いるいるか
るいるか
いるか
るか
か
6
10
4
8
15
17
19
1
12
21
3
7
14
11
5
9
16
18
20
2
13
いかいないかいるかいるいるいるか
いかいるかいるいるいるか
いないかいないかいるかいるいるいるか
いないかいるかいるいるいるか
いるいるいるか
いるいるか
いるか
いるかいないかいないかいるかいるいるいるか
いるかいるいるいるか
か
かいないかいないかいるかいるいるいるか
かいないかいるかいるいるいるか
かいるいるいるか
かいるかいるいるいるか
ないかいないかいるかいるいるいるか
ないかいるかいるいるいるか
るいるいるか
るいるか
るか
るかいないかいないかいるかいるいるいるか
34
るかいるいるいるか
接尾辞配列を用いた検索
• SA の上で二分探索
P = いるか
3回 (1, 12, 19)
SA
6
10
4
8
15
17
19
1
12
21
3
7
14
11
5
9
16
18
20
2
13
いかいないかいるかいるいるいるか
いかいるかいるいるいるか
いないかいないかいるかいるいるいるか
いないかいるかいるいるいるか
いるいるいるか
いるいるか
いるか
いるかいないかいないかいるかいるいるいるか
いるかいるいるいるか
か
かいないかいないかいるかいるいるいるか
かいないかいるかいるいるいるか
かいるいるいるか
かいるかいるいるいるか
ないかいないかいるかいるいるいるか
ないかいるかいるいるいるか
るいるいるか
るいるか
るか
るかいないかいないかいるかいるいるいるか
35
るかいるいるいるか
索引のサイズと検索時間
サイズ
転置ファイル
接尾辞配列
接尾辞木
< n bytes
頻度問い合わせ時間
O(m)
4n bytes + |T|
O(m log n)
>10n bytes + |T|
O(m)
注: 転置ファイルは文書が単語に分かれている場合
36
圧縮接尾辞配列 [Grossi, Vitter 00]
• 接尾辞配列を n log nO(n log ||) bits に圧縮
– 文書自身と合わせて |T|+O(n log ||) bits
– ヒトのDNAで13GB 1GB
• SA[i]の計算時間: O(1)O(log n) (01の定数)
– 頻度問い合せ: O(m log1+ n) 時間
– 列挙問い合せ: O(m log1+ n + occ log n) 時間
37
BW変換
T = acagcagg$
• BW[i] = T[SA[i]1]
• ソートした各接尾辞の1つ
前の文字を並べたもの
• T の並び替えになる
• BWから T を復元可能
• BWは圧縮しやすい
• bzip2で使われている
BW = g$ccaggaa
SA
9
1
3
6
2
5
8
4
7
g$
$ acagcagg
c agcagg$
c agg$
a cagcagg$
g cagg$
g g$
a gcagg$
a gg$
38
逆BW変換
• T は BW から復元できる
• SA も復元できる
g
$
c
c
a
g
g
a
a
$
acagcagg$
agcagg$
agg$
cagcagg$
cagg$
g$
gcagg$
gg$
39
FM-index
BWだけでパタン検索が可能
1
2
3
4
5
6
7
8
9
g
$
c
c
a
g
g
a
a
$
a
a
a
c
c
g
g
g
P=P[1, p] を検索
c = P[p]
l = C[c]+1, r = C[c+1]
while (--p  0) {
c = P[p]
l = C[c]+rankc(BW,l1)+1
r = C[c]+rankc(BW,r)
}
[l,r] は P の辞書順
C
$: 0
a: 1
c: 4
g: 6
40
FM-index
•
•
•
•
索引サイズ: nHk(S) + o(n) bits
パタンの検索: O(|P|) time
接尾辞配列の値の復元: O(log1+ n) time
長さ l の部分文字列の復元: O(l + log1+ n) time
41
de novo アセンブリ
DNA配列(未知)
シーケンサで read
(部分配列) を読む
assemble
元の配列を得る
人のDNAは3G文字
• readの長さは 30 ~ 1500
• 「次世代シーケンサ」では
– 超高速に読める
– readの長さが短い (100文字程度)
– エラーが多いので複数回読む (40~100回)
42
de Bruijn グラフ
• 文字列 T の k 次元 de Bruijn グラフ G(V,E)
– V: T 中の全ての長さ k の部分文字列 (k-mer) に対応
– E: {u  v | u, v  V, u = c1s, v = sc2, i c1sc2 = T[i..i+k]}
(s は (k-1)-mer)
• アルファベットサイズ  T = $$$TACGACGTCGACT$
• |T| = N, |E| = m
G
TAC
G
A
CGA
C
T
CGT
C
A
T
ACT
$
ACG
C
$TA
GAC
$$T
T
$$$
GTC
G
TCG
k=3
A
43
de Bruijn グラフを用いたアセンブリ
• 決定したいDNA配列の長さを n とする
• read の長さの合計を N とする (N  100n)
– 1つの read の長さは100程度
• 全 read 中の k-mer から de Bruijn グラフを作る
– k = 23 ~ 90
A
• オイラーパスを求める
GAC
C
T1 = TACAC
T2 = TACTC
T3 = GACAC
A
ACA
C
CAC $
$GA
A
TAC
T
ACT
C
CTC
$
$$G
C
$TA
G
A
$$T
T
44
$$$
ヒトゲノムでの実験結果
•
•
•
•
•
•
Readの長さ: 100, #reads 1,408,414,537
k = 27
異なる (k+1)-mer の数: 36,248,317,498
d = 3 回以上現れる(k+1)-mer の数: 4,998,165,473
dummy edge の数:
343,459,610
データ構造のサイズ: 2.5GB
– ABySS: 336GB
– Gossamer: 32GB
– Minia: 5.7GB
45
Gossamer [Conway, Bromage 11] のグラフ表現法
– グラフの枝は (k+1)-mer で表せる
節点のラベル (k 文字) と枝ラベル (1文字)
– 全 (k+1)-mer を疎ベクトルで表現し，圧縮
– 全ての枝は
 k 1 
  m1.44  (k  1) log   log m
log 2 

m


bits で表現できる
この表現は冗長
– ラベルの重なりがあるので，(k+1)-mer は非独立
– de Bruijnグラフの性質を使って圧縮する
46
Succinct de Bruijn Graph (SDG)
• (k+1)-merを，根元の節点の k-mer の逆辞書順に
ソートする
• サイズ: m(2+log ) bits
last Node
W
$
1
T
$$$
DNAだと 4m bits
1
C
CGA
A
T = $$$TACGACGTCGACT$
C
G
TAC
G
A
CGA
C
T
CGT
C
A
T
GTC
G
ACT
$
TCG
A
G
ACG
C
$TA
GAC
$$T
T
$$$
T
1
0
1
1
1
0
1
1
1
1
1
C
G
T
GG
A
T
AA
$
C
$TA
GAC
GAC
TAC
GTC
ACG
ACG
TCG
$$T
ACT
CGT
47
SDGの定義
• last[i] = 1 iff Node[i]  Node[i+1]
– Nodeが等しい i の範囲を表す
• W[i]  A  A
(A: T のアルファベット，A: フラグつき)
last

• W[i]  A iff  j < i s.t. W[j] = W[i] 1
1
かつ Node の長さ k1 の接尾辞が 1
0
等しい
1
1
• last[i] = 1 と W[j]  A に1対1対応
1
0
1
1
1
1
1
Node
$$$
CGA
$TA
GAC
GAC
TAC
GTC
ACG
ACG
TCG
$$T
ACT
CGT
48
W
T
C
C
G
T
GG
A
T
AA
$
C
SDGでの演算
• outdegree(v): 節点 v の出次数
• outgoing(v, c): 節点 v からラベルが c の枝を
たどって到達する節点
• indegree(v): 節点 v の入次数
• incoming(v, c): u = cs, v = sc’ を満たす節点 u
• node(v): 節点 v のラベル文字列 Node[v]
• index(s): Node[v] = s となる節点 v
49
fwd と bwd
• last[i] = 1 である Node[i] とグラフの節点は1対1
• Node[u] = cs, W[u] = c’ とする．
Node’[u] = sc’ と定義する．
• W[u]  A である全ての u で，Node’[u] は異なり，
last Node
W
Node[v] = sc’ である v が存在
1
T
$$$
C
• i.e. last[i] = 1 と W[j]  A に1対1対応 11 CGA
C
$TA
0
G
GAC
• fwd(u) = v, bwd(v) = u と定義
1
T
GAC
1
GTAC

• 注: W[i]  A のとき，
1
G
GTC
0
A
ACG
 j < i s.t. Node’[j] = Node[i]
1
T
ACG
1
ATCG
fwd(i) = fwd(j) と定義
1
A
$$T
1
1
ACT
CGT
50
$
C
fwd と bwd の計算
• fwd(u)
–
–
–
–
c = W[u]
r = rankc(W, u)
x = F[c] + r  1
v = select1(last, x)
注: W[u]  A なら対応する
A の文字に変換
• bwd(v)
–
–
–
–
x = rank1(last, v)
c = F1[x]
r = x  F[c] + 1
u = selectc(W, r)
• 共に O(1) 時間
$
A
C
G
T
F
0
1
3
6
8
last
1
1
1
0
1
1
1
0
1
1
1
1
1
Node
$$$
CGA
$TA
GAC
GAC
TAC
GTC
ACG
ACG
TCG
$$T
ACT
CGT
51
W
T
C
C
G
T
GG
A
T
AA
$
C
outdegree(v)
• 節点 v から出る枝ラベルは，W[vd+1..v]に格納
(d は出次数)
その範囲の last は 000…1 となっている
last Node
• d は last でのrank/selectで求まる
1
$$$
1
CGA
• O(1) 時間
1
$TA
G
TAC
G
A
CGA
C
GAC
T
ACT
$
T
CGT
C
GTC
G
TCG
A
ACG
C
$TA
A
$$T
T
$$$
0
1
1
1
0
1
1
1
1
1
GAC
GAC
TAC
GTC
ACG
ACG
TCG
$$T
ACT
CGT
52
W
T
C
C
G
T
GG
A
T
AA
$
C
outgoing(v, c)
• 節点 v から出る枝ラベルの格納されている
範囲で c または c を探す
–
–
–
–
x = selectc(W, rankc(W, v))
y = selectc(W, rankc(W, v))
x, y 共に範囲外なら枝は存在しない
w = fwd(x)
• O(1) 時間
last
1
1
1
0
1
1
1
0
1
1
1
1
1
$
A
C
G
T
Node
$$$
CGA
$TA
GAC
GAC
TAC
GTC
ACG
ACG
TCG
$$T
ACT
CGT
F
0
1
3
7
10
W
T
C
C
G
T
GG
A
T
AA
$
C
53
まとめ
• 簡潔データ構造は，通常のデータ構造を限界
まで圧縮しつつ，高速な問い合わせを実現
• 最近の話題
–
–
–
–
動的 / オンライン (末尾追加)
文字列の文法圧縮
BDD, ZDD
最短路探索
• 詳しい資料
– http://researchmap.jp/sada/succinct/
54
レポート課題
• 次のページを参照してください
http://www.keisu.t.u-tokyo.ac.jp/kmb_seminar/details/frontier_2016mistS.html
• 締切5月25日
55

Download Report