haripo.com

Dynamic Time Warping(動的時間伸縮法)

Dynamic Time Warping(DTW, 動的時間伸縮法)を最近知りました。 波形同士のアライメントと差分計算を同時に行い、類似度を求める手法です。 初出は 1960 年代で、音声認識や文字認識などの分野で活発に研究されていた1ようです。 シンプルな手法なので、勉強を兼ねてデモを実装しました。

デモ

枠の中はマウスドラッグで波形が描けるようになっています。 上の枠と下の枠に波形を描くと、波形同士のアライメント(warping path)が表示されます。

手法

2 つの波形データ xx, yy の類似度を計算することを考えます。 波形データは、時刻 t(1tT)t (1 \leq t \leq T) における波形の値 xt,ytx_t, y_t で表されているとします。

一番簡単な類似度の計算は、波形同士の差分を求めることでしょう。

txtyt \sum_t |x_t - y_t|

しかし、この差分は波形がずれているときにうまくいきません。 波形の形に応じてアライメント処理を行い、差分を計算する箇所を変化させる必要があります。

定義

まず、波形同士のアライメントを表現する記法を説明します。 波形 xx の時刻 t=wxt=w^x と波形 yy の時刻 t=wyt=w^y がアライメントによって対応することを w=(wx,wy)w = (w^x, w^y) とし、波形のアライメント全体を KK 個の対応関係の集合 WW で表します。 W={w1,w2,,wK} W = {w_1, w_2, \cdots, w_K }

wk=(wkx,wky)w_k = (w^x_k, w^y_k) とおくと WW は次の条件を満たします。

  1. wk1xwkx1w^x_{k-1} - w^xk \leq 1 かつ wk1ywky1w^y{k-1} - w^y_k \leq 1
  2. wk1xwkx0w^x_{k-1} - w^xk \geq 0 かつ wk1ywky0w^y{k-1} - w^y_k \geq 0
  3. w1=(1,1)w_1 = (1, 1) かつ wK=(T,T)w_K = (T, T)

このとき、Dynamic Time Warping による波形同士の距離は以下のように定義されます。

DTW(x,y)=minKkxwkxywky \mathrm{DTW}(x, y) = \min\sumK^k |x{wk^x} - y{w_k^y}|

アルゴリズム

数式で書くと添字が複雑ですが、アルゴリズムはシンプルです。 まず、波形同士の距離行列 DD を考えます。

D=(x1y1x2y1xTy1x1y2x2y2xTy2x1yTx2yTxTyT) D = \left( \begin{array}{cccc} |x_1 - y_1| & |x_2 - y_1| & \cdots & |x_T - y_1| \ |x_1 - y_2| & |x_2 - y_2| & \cdots & |x_T - y_2| \ \vdots & \vdots & \ddots & \vdots \ |x_1 - y_T| & |x_2 - y_T| & \cdots & |x_T - y_T| \end{array} \right)

WW は距離行列 DD 上の経路として捉えることができます。 D0,0D{0, 0} から DT,TD{T, T} までをつなぐ経路です。 たとえば W={(0,0),(1,1),(1,2),}W={(0, 0), (1, 1), (1, 2), \cdots} とすれば、 これは行列 DD 上の D0,0,D1,1,D1,2D{0, 0}, D{1, 1}, D_{1, 2} \cdots を通る経路とみなせるわけです。 上述した WW の満たす 3 つの条件はそれぞれ経路が

  1. 連続である(経路が途切れない)こと
  2. 単調である(経路が後戻りしない)こと
  3. 経路の始点が (1,1)(1,1) で終点が (T,T)(T,T) であること

を意味します。

様々な経路が考えられますが、このうち経路上にある DD の要素の総和が最小となるのが最適な WW というわけです。 そのときの経路上にある要素の総和が DTW(x,y)\mathrm{DTW}(x, y) です。 動的計画法で解けますね!

実装

擬似コードでの実装は Wikipedia が詳しいです。 ここでは上記デモの実装を紹介します。 デモの実装全体は haripo/DynamicTimeWarping にあります。

var map = [];
var cache = [];

for (var i = 0; i < stream1.length; i++) {
  distance.push([]);
  cache.push([]);
  for (var j = 0; j < stream2.length; j++) {
    distance[i].push(Math.abs(stream1[i] - stream2[j]) + 1);
    cache[i].push(Math.abs(i - j) < 150 ? -1 : Number.MAX_VALUE);
  }
}

距離行列 DD を二次元配列 distance に作って、動的計画法をメモ化するための配列 cache を初期化しています。 波形の遠い箇所がマッチングしてしまうのを防ぐため Math.abs(i - j) < 150 の範囲で探索するようにしています。 このような境界条件をつけることで探索が高速になるほか、実データではロバストになることが分かっています2

var search = null;
search = function(x, y) {
  if (cache[x][y] <= 0) {
    cache[x][y] = distance[x][y] + Math.min(
      (x > 0 ? search(x - 1, y) : Number.MAX_VALUE),
      (y > 0 ? search(x, y - 1) : Number.MAX_VALUE),
      (x > 0 && y > 0 ? search(x - 1, y - 1) : Number.MAX_VALUE)
    );
  }
  return cache[x][y];
};

search(stream1.length - 1, stream2.length - 1);

再起関数を使った普通の動的計画法です。1 行目の var search = null; がないと search 関数中で search 関数が未定義状態となるので必要です。 この時点で cache[stream1.length - 1][stream2.length - 1] が DTW 距離となります。

var x = 0;
var y = 0;
var match = [];
while (x > stream1.length - 1 || y > stream2.length - 1) {
  if (cache[x + 1][y + 1] <= cache[x + 1][y] &&
      cache[x + 1][y + 1] <= cache[x][y + 1]) {
    x += 1;
    y += 1;
  } else if (cache[x][y + 1] <= cache[x + 1][y]) {
    y += 1;
  } else {
    x += 1;
  }
  match.push([x, y])
}

今回のデモでは Warping path を表示したかったので、stream1 と stream2 の各点がどの位置にアライメントされるのかを求めます。 cache を辿って最短ルートで cache[stream1.length - 1][stream2.length - 1] に至るパスを探索しています。

DTW の実装はこれだけです。シンプル。

おわりに

DTW を紹介しました。 わりと古くからある手法ですが、現在も DTW に関係する論文が出ているようです3。 実装も簡易でわかりやすいですし、時系列データに「とりあえず」で試してみてもいいかもしれません。

  1. Senin, Pavel. “Dynamic time warping algorithm review.” Information and Computer Science Department University of Hawaii at Manoa Honolulu, USA 855 (2008): 1-23.
  2. Xi, Xiaopeng, et al. “Fast time series classification using numerosity reduction.” Proceedings of the 23rd international conference on Machine learning. ACM, 2006.
  3. Google Scholar で 2016 年以降の DTW 論文を検索した結果