Dynamic Time Warping（動的時間伸縮法）

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

デモ

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

手法

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

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

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

定義

まず、波形同士のアライメントを表現する記法を説明します。 波形 $x$ の時刻 $t=w^x$ と波形 $y$ の時刻 $t=w^y$ がアライメントによって対応することを $w = (w^x, w^y)$ とし、波形のアライメント全体を $K$ 個の対応関係の集合 $W$ で表します。

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

1. $w^x_{k-1} - w^x_k \leq 1$ かつ $w^y_{k-1} - w^y_k \leq 1$
2. $w^x_{k-1} - w^x_k \geq 0$ かつ $w^y_{k-1} - w^y_k \geq 0$
3. $w_1 = (1, 1)$ かつ $w_K = (T, T)$

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

アルゴリズム

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

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

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

を意味します。

様々な経路が考えられますが、このうち経路上にある $D$ の要素の総和が最小となるのが最適な $W$ というわけです。 そのときの経路上にある要素の総和が $\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);
  }
}

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

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 に関係する論文が出ているようです³。 実装も簡易でわかりやすいですし、時系列データに「とりあえず」で試してみてもいいかもしれません。