論文
DBSJ Letters Vol.4, No.1――――――――――――――――――――――――――――――――――――

タイムワーピングに基づく時系列
データの類似検索：次元縮小による効率化
Efficient Search of Similar Time Series
under Time Warping with Dimensionality
Reduction

大桃 諭♥
古瀬 一隆♦陳 漢雄♦
大保 信夫♠

Satoshi OOMOMO Hanxiong CHEN
Kazutaka FURUSE Nobuo OHBO
時系列データの類似度にユークリッド距離が使用されていたが，時間軸方向のずれに柔軟性を持たせたタイムワーピン
グ距離が近年注目をあびている．しかし，タイムワーピング
距離の計算時間は非常に長く，しかも距離公理の三角不等式
が成り立たないため，伝統的な索引技法の適用は困難である．
そこで本研究では，時系列データに対する新しい次元縮小法と下界関数を提案する．そして，この下界関数をフィルタリ
ングに使用したときの検索効率が，関連研究で提案されている下界関数によるものよりも優れていることを実験によって実証する．
Recently, similarity between time series data measured
by time warping (TW) distance which allows out of phasein the time axis has got a lot of attention in recent years.
However, since TW violates the triangle inequality,
traditional indexing techniques is powerless. To overcome
this problem, in this paper we propose a new lower boundfunction based on dimensionality reduction. This enables
efficient indexing hence filtering in similarity search oftime series data. We show by experimental results that
our lower bound function outperforms other ones proposedin related works.

1. はじめに
時系列データとは，時刻iの時の値をviとすると，V= v1,
v2…,vNで表される実数のシーケンスである．時系列データの例としては，株価や為替レート，心電図，気温・湿度変化などがあり，ここに挙げた以外にも多く存在する．このように科学，医療，経済，工学などの様々な分野で時系列データが扱われている．さらに近年では，扱うデータ量も非常に多くなってきている．このような現状において，時系列データの効率的な類似検索の必要性が高まっている．例えば，最近の株価の動向と類似したものを過去のデータから見つけ出して，今後の動向を予測するといったことが挙げられる．

類似検索を行うためには，2つの時系列間の類似度を定義する必要がある．現在，時系列の類似度として最もよく使用されているものはユークリッド距離である．しかしユークリ
ッド距離には，2つの時系列の長さが等しくなければならな
い，時間軸方向のわずかなずれに対しても大きな影響を受けてしまうなどの問題がある．そこで，ユークリッド距離に代わるものとしてタイムワーピング距離が使用されている．
類似検索に要する計算コストは，次の2つの要因に大きく依存する．1つ目は時系列の長さである．ユークリッド距離は時系列の長さに比例したCPUコストで済むが，タイムワー
ピング距離は長さの2乗に比例したCPUコストを要する．この問題に対処するために使用される手法の1つが次元縮小である．離散フーリエ変換（DFT），離散ウェーブレット変換（DWT），
Piecewise Aggregate Approximation（PAA），Adaptive
Piecewise Constant Approximation（APCA）などが使用されている．2つ目はデータベースに含まれる時系列の数である．
時系列の数が増えることによって，CPUコストだけではなくI/Oコストも大きく増加する．そこで，R-treeなどのインデックスを使用する手法が提案されている．
しかしこれらの手法の中には，false dismissalを生じてしまうものがある．false dismissalとは，本来検索結果として返されるべき解が返されないことを意味する．このfalse dismissalを生じさせないで計算コストを小さくする
手法として使用されているものが，下界関数を使用したフィ
ルタリングである[2]~[4]．
本研究では，タイムワーピング距離を使用して類似検索を行う場合における効率的な検索手法を提案する．提案手法では，与えられた時系列に対してその時系列を包含する領域ボックス列というものを作成する．そして，この領域ボックス列を使用して下界関数を定義する．この下界関数をフィルタリングに使用することで，効率的な検索を可能にする．

2. タイムワーピング
この章では，タイムワーピング距離の計算方法について簡単に述べる([1])．
長さn,mの2つの時系列Q=q1,…,qi,…,qn，C=c1,…,cj,…,cmに対して，まずは，以下の式で定義される2つの点qi,cj間の距離d(qi,cj)を(i,j)要素の値とするn×m行列を作成する．

d (qi , c j ) = (qi − c j )2(1)

次に，ワーピングパスW=w1,w2,…,wk,…,wKを求める．ワーピ
ングパスは，以下の3つの条件を満たす行列の要素の順列で表現される．ワーピングパスの例を図1に示す。

♥

学生会員 筑波大学大学院システム情報工学研究科
oomomo@dblab.is.tsukuba.ac.jp♦
正会員 筑波大学大学院システム情報工学研究科{chx, furuse}@cs.tsukuba.ac.jp♠
筑波大学大学院システム情報工学研究科
ohbo@cs.tsukuba.ac.jp

1

図 1
Fig.1

ワーピングパスの例

An example of warping path

日本データベース学会 Letters Vol.4, No.1

論文
DBSJ Letters Vol.4, No.1――――――――――――――――――――――――――――――――――――[境界条件] w1=(1,1),wK=(n,m)とする．
[連続性] wk=(a,b),wk-1=(a',b')とすると，a-a'≦1
かつb-b'≦1となる．
[単調性] wk=(a,b),wk-1=(a',b')とすると，a-a'≧0
かつb-b'≧0となる．
上記の3つの条件の下で得られるワーピングパスは非常に多く存在するが，それぞれのワーピングパスに対して(2)式で示すワーピングコストと呼ばれるものが存在する．

C (W ) =

∑

K
k =1

wk

ワーピングパスは非常に多く存在するので，全てのワーピン
グパスに対してワーピングコストを計算して最小値を求めるのは現実的ではない．そこで，以下の再帰関数を使用してタイムワーピング距離を計算する．

γ (0,0) = 0, γ (i,0) = γ (0, j ) = ∞

(4)

(5)

この再帰関数を使用すると，式(3)は以下のように表される．

D tw (Q, C) =γ (N , M )(6)

これは動的プログラミングによってO(NM)で計算可能である．

3. 関連研究2 つの時系列 Q, C に対して本来の距離を D(Q,C)とすると，以下の式を満たす D lb (Q, C) が下界関数である．

D lb (Q, C) ≤ D ( Q , C )この下界関数をフィルタリングに使用することで，効率的な検索を行うことができる．例えば，問い合わせ Q との距離が e 以下の時系列を検索する場合，データベースから取り出した時系列 C に対して下界距離を計算し，e より小さい場合にのみ，本来の距離 D(Q,C)を計算すればよい．
また，下界関数には性能の善し悪しがある．1 つは下界距離自身の計算コストである．最低でも本来の距離の計算コス
トよりも小さくなければ，下界関数をフィルタリングに使用することは意味を成さない．計算コストが小さくなるほどフィルタリングに要する時間が短くなるので，結果として検索時間も短くなる．もう 1 つは下界関数の精度，すなわち下界
距離と本来の距離との差である．極端な例であるが，0 も当
然下界距離になれる。しかし，e がどんなに小さな正の値であっても，結局は D(Q,C)を計算することになってしまい，フィルタリング効果がまったくない．下界距離が D(Q,C)に近い
ほどフィルタリング効果が大きくなり，検索時間が短くなる．
Yi らは[5]において以下のような下界関数を定義している．
第 1 章で述べたように，タイムワーピング距離の計算では一方の時系列の点は他方の時系列のいずれかの点と必ず対応させなくてはならない．このことを利用して，一方の時系列の点の中で他方の時系列の最大の点より大きいもの，もしくは最小の点より小さいものを見つけて，それらの点の値の差を用いて下界関数を定義している．図で表すと，図 2 の点線で表される距離の和が下界距離になる．この手法は CPU コス
トが O(N+M)で済むが，データベースから全ての時系列を取り

2

図 2

Yi による下界関数（左）と Kim による下界関数（右）
Fig.2

(2)

この式から得られるワーピングコストの中で最小のものがタイムワーピング距離になる．
(3)
D tw (Q, C) = min{ C (W )}⎧ γ ( i − 1, j − 1) ⎫⎪⎪γ ( i , j ) = d ( q i , c j ) + min ⎨ γ ( i − 1, j ) ⎬⎪ γ ( i , j − 1)⎪⎩⎭

出す必要があるので I/O コストは減少しない．

Lower bound by Yi (left) and Kim (right),
respectively

Kim らも[4]においてそれぞれの時系列から 4 次元の特徴
ベクトルを抽出する．特徴ベクトルは時系列の最初の値
First(C)，最後の値 Last(C)，最大値 Max(C)，最小値 Min(C)から成る．そして 2 つの特徴ベクトルが与えられたとき，それぞれの要素の値の差の中で最大のものが下界距離になる．
図 2 では最後の値の差 Last(Q)-Last(C)が下界距離となる．
この特徴ベクトルを 4 次元空間上に配置して R-tree などのインデックスを構築することで，計算コストを減少させている．この手法はデータベースが大きい場合には有効であるが，下界関数自身の精度はよくない．
Keogh らは[2][3]において時系列を固定長な領域ボックスに区切り MBR を用いる下界関数を提案した．この方法は下界関数の精度をあげた一方，事前に領域ボックスの数を決めなければならないため規則性のない時系列データには無駄が生じる．

4. 提案手法
我々は時系列データにあわせた可変長な領域ボックスを提案し，それを用いてKeoghの下界関数([3])を改良する．つ
ぎにまず，ある時系列Q=q1, …, qNが与えられたときに，その時系列を包含する領域ボックス列RQ=rQ1, …, rQNを作成する方法を述べる．

4.1

領域ボックス列の作成

まずは，時系列Qを重なりのない複数の部分時系列Q=<q1, …,
qh>,<qh+1, …, qk>,…, <qk+1, …, qN>,に分割する．そして，それぞれの部分時系列に対して 1 つの領域ボックスを作成す
る．1 つの領域ボックスは上限U，下限L，サイズSの 3 つの要素から成り，部分時系列<qh+1, …, qk>に対して作成された領域ボックスをrQiとすると，以下のように定義される．

rQi.U=max{qh+1, …, qk}, rQi.L=min{qh+1, …, qk}, rQi.S=k-h,
領域ボックス列の作成において中心となる処理は部分時系
列への分割である．そして，この分割方法が検索効率に大きな影響を与える．影響を与える 1 つの要因は領域ボックスの数 N' で あ り ， も う 1 つ は 領 域 ボ ッ ク ス の 面 積N '∑ (r i .U - r i .L ) r i .S である．下界距離の計算方法についi =1ては次の章で詳しく説明するが，数が少なくなるほど下界距離の計算コストは小さくなり，面積が小さくなるほど下界距
離によるフィルタリング効果は上がるが，数と面積はトレードオフの関係にある．そこで，本研究では図 3 のアルゴリズ
ムを用いて領域ボックス列を生成する．
1.
2.

Input: 許容振幅率ρ∈(0,1]τ=ρ・(max(Q)-min(Q))

3.

i=1, n=1

4.

find minimum j(i≦j<N) such thatmax {q i , L , q j }− min {q i , L , q j } ≤ τ

日本データベース学会 Letters Vol.4, No.1

論文
DBSJ Letters Vol.4, No.1――――――――――――――――――――――――――――――――――――
max {q , L , q }− min {q , L , q } > τ⎧ ( ri Q . L − r jC .U ) 2 , if ri Q . L − r jC .U
j +1

i

5.
6.

i

j +1

if j dose not exist then j=N
r i q .U = max{q , L , q }, r . L = min {q , L , q }
q

i

j

i

i

j

ri q . S = j − i + 1

if j==N return RQ, otherwise goto 4
図 3

領域ボックス列作成アルゴリズム

Fig.3 Region box series creation algorithm
このアルゴリズムに従って，τ=８として領域ボックス列を作成する様子を図４に示す．点 A から点 B までの部分時系
列の最大値と最小値の差は 7 であり，τ 以下の値となっているが，点 B の次の点である点 C まで進むと差は 10 になり，
閾値より大きくなる．よって，点 A から点 B までが 1 つの領域ボックスとなり，点 C が次の領域ボックスの開始点になる．

⎪⎪
d ' ( ri Q , r jC ) = ⎨ ( ri C . L − r jQ .U ) 2 ,
⎪
0,
⎩⎪

if

ri C . L − r jQ .U

otherwise⎧γ 'd (i − 1, j − 1) + d ' ( ri Q , r jC ) • r jC .S⎪⎪
γ 'r (i , j ) = min ⎨γ 'c (i − 1, j ) + d ' ( riQ , r jC )⎪Q
CC⎪⎩γ 'r (i , j − 1) + d ' ( ri , r j ) • r j .S
⎧γ ' d (i − 1, j − 1) + d ' ( riQ , r jC ) • riC .S⎪⎪
γ ' c (i, j ) = min ⎨γ ' c (i − 1, j ) + d ' ( riQ , r jC ) • riC .S⎪
QC⎪⎩γ ' r (i, j − 1) + d ' ( ri , r j )⎧γ ' (i , j )γ ' d (i , j ) = max ⎨ r⎩γ ' c (i , j )γ ' d ( 0 , 0 ) = 0, γ ' r ( i , 0 ) = γ ' c ( 0, j ) = ∞この再帰関数を使用して下界関数を以下のように定義する．

D tw_lb (R Q , R C ) = γ ' d ( N ' , M ' )
図 4
Fig.4

領域ボックス列の例

An example of region box series.

閾値 τ の値によって領域ボックスの数や面積は決まるが，時系列の取り得る値の範囲が事前にわかっていない限り，閾
値 τ を直接指定することは困難である．そこで，許容振幅率ρというものをパラメータとし，このρから閾値 τ を計算する（図３の２行目）．ρは時系列の取り得る値の範囲に関係なく，0 から 1 の値で指定することができる．ρが小さいほど領域ボックスの面積は小さくなるが数は多くなる．面
積と数はトレードオフの関係にあり，より効率的な検索を行うためにはユーザがρに適切な値を指定する必要がある．

図 5
Fig.5

4.2

領域ボックス間の距離

Distance between region boxes.

改良下界関数および検索アルゴリズム

2 つの時系列Q=q1,…,qi,…,qN，C=c1,…,cj,…,cMに対して，それぞれの領域ボックス列RQ=rQ1, …, rQNʼ, RC=rC1, …, rCMが作成されたとする．この 2 つの領域ボックス列を使用して，タイムワーピング距離Dtw(Q,C)の下界関数Dtw̲lb(RQ, RC)を定義する．本文は領域ボックスが可変長であるため，本質的には異なるが，表現上の便宜から[3] と同じ形式の計算式を利用する．まずは式(1)で表される 2 点間の距離の代わりとなる，図 5 で示すような 2 つの領域ボックスrQi，rCj間の距離dʼ(rQi，rCj)を定義する．
次に再帰関数(4)，(5)をもとに，下界関数のための新たな再帰関数 γ ' d を定義する．

3

これは動的プログラミングによって O(N'M')で計算可能である．紙面の関係で証明を省略するが、このように定義された関数は以下の式を満たし，よって下界として機能する．

D tw_lb (R Q , R C ) ≤ D tw_lb (Q, C)
この下界関数を用いて類似検索に対するフィルタリングが可能である．類似検索には大きく分けると近傍検索と範囲
検索の 2 種類がある．ここでは k=1 のときの最近傍検索アルゴリズムを説明するが，k> 1 の場合の近傍検索や範囲検索もアルゴリズムを少し変更することで行うことができる．
実際に検索を行う前に，データベース中の全ての時系列に対して領域ボックス列を作成してファイルに格納しておく．
このファイルを圧縮ファイルと呼ぶことにする．問い合わせ
時系列 Q が与えられると，次のアルゴリズムを用いて最近傍検索を行う．
1. Qから領域ボックス列RQを作成する．
2. 圧縮ファイルから領域ボックス列RCを 1 つ取り出して下界距離Dtw̲lb(RQ, RC)を計算する．
3. 下界距離がこれまでの最小距離より小さい場合には，
データベースから時系列を取り出して本来の距離D Dtw(Q,C)
を計算する．これも最小距離より小さい場合には最小距離を更新する．
4. 圧縮ファイルから次の領域ボックス列を取り出して同様の処理を行う．これを全ての時系列に対して行って，最終的に最小距離となった時系列を検索結果として返す．

5. 実験
今回の実験では株価データセットとvt=vt-1+rt(rtは[-1,1]の乱数)から生成されるRandomWalk合成データセットの 2 種類を使用する．合成データセットは長さが 256，512，1024 のRandomWalkデータをそれぞれ 100，1000，10000 個生成して作成することで，合計 9 つの合成データセットを作成した．
株 価 デ ー タ セ ッ ト は S&P 500 (http://biz.swcp.com/
stocks/)から 479 個の銘柄の始値，高値，安値，終値の 4 種
類のデータ 252 日分（データ長が 252，データ数が 479×4）

日本データベース学会 Letters Vol.4, No.1

論文
DBSJ Letters Vol.4, No.1――――――――――――――――――――――――――――――――――――からなる．問い合わせデータは，データセット中のデータと同じ長さのRandomWalkデータを新たに生成して使用する．また，それぞれのデータに対して平均値が 0，標準偏差が 1 になるように正規化を行う．
表 1

Table 1
datasetLen. data
Num. data
Serial
Yiʼs
Kimʼs
RB
(ρ)

フィルタリング率の比較

Comparison on filtering rate
Synthetic dataset
256
512
1024
10000
10000
10000
0
0
0
17.7
25.0
23.9
36.8
26.5
15.9
30.2
32.3
40.1
(0.20)(0.20)(0.20)

表 2

ディスクアクセス数の比較

Table 2Comparison on disk accessdatasetLen. data
Num. data
Serial
Yiʼs
Kimʼs
RB
(ρ)

Synthetic dataset
256
512
1024
10000
10000
10000
20040
40040
80040
20040
40040
80040
19581
37286
76233
11505
16313
26402
(0.13)(0.13)(0.09)
表 3

Table 3

stock
252
1916
0
10.6
24.4
37.3
(0.20)

stock
252
1916
3780
3780
4400
2102
(0.13)

実行時間の比較

Comparison on execution time

datasetSynthetic dataset
stock
Len. data
256
512
1024
252
Num. data
10000
10000
10000
1916
Serial
35.6
136.4
543.9
6.3
Yiʼs
29.6
105.3
419.9
5.7
Kimʼs
22.9
102.7
459.9
4.8
RB
6.26
17.31
49.31
1.16
(ρ)
(0.10)(0.09)(0.06)(0.09)
本研究で提案した領域ボックスを使用した下界関数をフィルタリングに使用する RB 検索手法と 3 つの手法の検索効率を比較する．すなわち，下界関数でのフィルタリングを行わずに，直接タイムワーピング距離を計算する逐次検索手法，
第 3 章で述べた Yi らや Kim らの下界関数によるフィルタリ
ング検索手法である．領域ボックス列を作成するときに指定する許容振幅率ρには，0.01〜0.20 の 0.01 刻みの値を使用する．50 個の異なった問い合わせデータで最近傍検索を行って平均結果を計測し，表１〜３でその一部の結果を示す。
まずフィルタリング率（表１）について比較した．この値が大きいほど実際のタイムワーピング距離の計算が減り，検索効率が良いと判断できる．RB 検索手法では許容振幅率を大きくするとフィルタリング率が小さくなっている．これは，許容振幅率を大きくすると領域ボックスの面積が大きくなるので，下界距離が小さくなってフィルタリング効果が下がるためと考えられる．しかし，許容振幅率を 0.2 まで大きく

4

してもフィルタリング率は 30〜40%であり，Yi や Kim の 10〜40%に比べて高いことがわかる．
つぎにディスクアクセス数（表２）について比較した．

いうまでもなく，この値が小さい検索効率が良いと判断できる．RB 検索手法では，圧縮ファイルへのアクセスとデー
タベースへのアクセスのバランスが取れた箇所でディスク
アクセス数が最小になる．ほかの方法より顕著に減少した．
最後に検索効率を総合的に評価するため問い合わせを与えてから検索結果が得られるまでの実行時間（表３）についても比較した．データの数や長さにもよるが，Yi や Kim では最大でも逐次検索手法の 65%程度までしか実行時間が短くならないが，RB 検索手法では 30%〜10%程度まで実行時間が短くなっている．

6. おわりに本研究では，時系列から可変長領域ボックス列を作成してKeogh の下界関数を新しく定義し，false dismissal を生じることなく効率的な類似検索を行うことを試みた．最近傍検索の実験によって，提案手法が関連研究の手法よりも優れた性能を持っていることが実証された．今後の課題としては，今回は時系列を領域ボックス列に次元縮小して圧縮ファイルを作成しただけなので，R-tree などのインデックスを構築
して下界距離の計算回数をさらに減らすことが考えられる．

［文献］[1] J. B. Kruskall and M. Liberman, ``The SymmetricTime-Warping Problem: From Continuous to Discrete,''Time Warps, String Edits, and Macromolecules: The Theory andPractice of Sequence Comparison, pp.125-161, Addison-Wesley,
1983[2] E. Keogh, K. Chakrabarti, S. Mehrotra, and M. Pazzani,``Locally Adaptive Dimensionality Reduction for Indexing
Large Time Series Databases,'' In Proc. of 2001 ACM SIGMOD,
2001.
[3] E. Keogh, ``Exact Indexing of Dynamic Time Warping'' In
Proc, VLDB, 2002., pp. 406-417.
[4] S.-W. Kim, S. Park, and W. W. Chu, ``An Index-BasedApproach for Similarity Search Supporting Time Warping in
Large Sequence Databases,'' In Proc. of 17th ICDE, 2001.
[5] B.-K. Yi, H. V. Jagadish, and C. Faloutsos, ``Efficient Retrievalof Similar Time Sequences under Time Warping,'' In Proc. of
14th ICDE, 1998.

大桃 諭

Satoshi OOMOMO

2003 年筑波大学情報学類卒業．現在，同大大学院博士課程シ
ステム情報工学研究科在学中．データベースシステムに関する研究に従事．日本データベース学会学生会員．

陳 漢雄

Hanxiong CHEN

1993 年筑波大学大学院博士課程工学研究科了．現在，同大電
子・情報工学系講師．データベースシステムに関する研究に従事．工博．日本データベース学会正会員．

古瀬 一隆

Kazutaka FURUSE

1993 年筑波大学大学院博士課程工学研究科了．現在，同大電
子・情報工学系講師．データベースシステムに関する研究に従事．工博．日本データベース学会正会員．

大保 信夫

Nobuo OHBO

1968 年東京大学理学部卒．1970 年同大大学院修士課程了．
同年同大理学部助手．1980 年筑波大学電子・情報工学系講師．
1995 年同大電子・情報工学系教授．現在に至る．データベー
スシステムに関する研究に従事．理博．

日本データベース学会 Letters Vol.4, No.1