Soft-DTW（Differentiable Dynamic Time Warping）

上級

04.14.5

Soft-DTW（Differentiable Dynamic Time Warping）

最終更新 2020-02-26 読了時間 6 分

まとめ

DTW の累積コスト更新に出てくる $\min$ は、3 つの候補のうち 1 つだけを選ぶ「角のある」演算なので微分できず、ニューラルネットの損失関数には使えません。Soft-DTW はこの $\min$ を、スムージング係数 $\gamma$ つきの log-sum-exp で定義した soft-min に置き換え、どこでも微分できる距離にしたものです。soft-min は常に $\min$ 以下の値を返すので、$\gamma$ を 0.01 のように小さくすると DTW に一致し、大きくすると全経路をまとめて評価した値へ離れていきます。このページでは NumPy の soft_dtw を 20 行で実装し、$\gamma$ を 0.01 から 5.0 まで 4 段階に変えて DTW と比べたあと、tslearn.metrics.soft_dtw_gradient で入力系列に対する勾配を取り出します。

動的時間伸縮（DTW）で累積コスト行列と $\min$ による遷移式を先に見ておくと、このページで置き換える箇所が 1 行だけだと分かります。

Cuturi, Marco, and Mathieu Blondel. “Soft-DTW: a Differentiable Loss Function for Time-Series.” ICML 2017. (PDF)

直感 #

迷路で最短経路を 1 本だけ選ぶのが DTW です。壁を 1 mm 動かしても最短経路は同じか、別の経路に突然切り替わるかのどちらかで、「経路の長さが壁の位置に対してどう変わるか」という傾きが定義できません。Soft-DTW は「すべての経路を、短いものほど大きな重みで足し合わせる」と考え直します。$\gamma$ が小さいと最短経路だけが残って DTW と同じ値になり、$\gamma$ が大きいと遠回りの経路の寄与も混ざって、1 本の経路の長さとは別の値になります。この足し合わせは壁の位置をわずかに動かすと滑らかに変わるので、勾配が定義でき、「予測した波形を正解の波形に近づける」学習に使えます。

系列 x, y 局所距離 D_ij = (x_i − y_j)^2 3 方向の累積コストを soft-min で合成（γ で滑らかさ調整） R_nm が Soft-DTW 後方パスで勾配

Soft-DTW とは #

Dynamic Time Warping（DTW）は最小コスト経路を厳格に選ぶため、距離関数は非凸であり、ワーピングパスも 1 通りに決まります。一方で、ニューラルネットの損失関数として DTW をそのまま使おうとすると、$\min$ が微分できないため勾配を求められません。

Soft-DTW はこの「最小」を soft-min（滑らかな最小演算）に置き換えることで、なめらかな距離関数を定義します。スムージング係数 $\gamma > 0$ を導入し、累積コスト更新式を

$$ R_{i,j} = D_{i,j} + \text{softmin}\bigl(R_{i-1,j},\ R_{i,j-1},\ R_{i-1,j-1}\bigr) $$

とします。ここで

$$ \text{softmin}(a,b,c) = -\gamma \log\bigl(\exp(-a/\gamma) + \exp(-b/\gamma) + \exp(-c/\gamma)\bigr) $$

と定義すると、$\gamma$ を小さくすると DTW に漸近し、$\gamma$ が大きいほど平均化された距離になります。$a$ が他より 10 以上小さければ $\exp(-a/\gamma)$ だけが残って $\min$ と同じ値を返し、3 つが近い値のときだけ $\min$ より少し小さい値になるのが soft-min の性質です。Soft-DTW は勾配を解析的に求められるため、ニューラルネットワークなどで時系列の形状を意識した損失関数として用いられます。

何が起きているか（ステップで見る） #

コードを読む前に、この手法が図の上で何をしているかを追ってみましょう。「再生」を押すか、スライダーを動かすと 1 手順ずつ進みます。

見るポイント：DTW は時間を伸び縮みさせて形を合わせる。タイミングがずれた似た波形を「近い」と判定できる。

アルゴリズムの流れ #

まず通常の DTW と同様に、$D_{i,j} = |x_i - y_j|^2$ などで局所距離の行列を作ります。次に累積コスト行列 $R$ を初期化し、上の soft-min を使って動的計画法で左上から右下へ埋めていきます。$\exp(-a/\gamma)$ は $a$ が大きいと 0 にアンダーフローするため、実装では最小値を先に引いてから指数を取る log-sum-exp の技巧を使います。最後に Soft-DTW の値を $R_{n,m}$（必要に応じて正規化）として取り出します。$\gamma$ が 0 に近づくほど経路が鋭くなり、DTW に収束します。計算量は DTW と同じ $O(nm)$ ですが、各マスで exp と log を呼ぶ分だけ定数倍だけ重くなります。

Python で Soft-DTW を実装する #

以下は NumPy だけで Soft-DTW を計算する簡易実装です。DTW の実装と比べると、min(...) が softmin_three(...) に替わり、局所距離が絶対値から二乗になった 2 点だけが違います。

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
import numpy as np

def pairwise_sq_euclidean(x: np.ndarray, y: np.ndarray) -> np.ndarray:
    return (x[:, None] - y[None, :]) ** 2

def softmin_three(a: float, b: float, c: float, gamma: float) -> float:
    values = np.array([a, b, c], dtype=float)
    m = values.min()
    return m - gamma * np.log(np.exp(-(values - m) / gamma).sum())

def soft_dtw(x: np.ndarray, y: np.ndarray, gamma: float = 1.0) -> float:
    if gamma <= 0:
        raise ValueError("gamma は 0 より大きい値にしてください。")

    x = np.asarray(x, dtype=float)
    y = np.asarray(y, dtype=float)
    n, m = len(x), len(y)
    D = pairwise_sq_euclidean(x, y)

    R = np.full((n + 1, m + 1), np.inf)
    R[0, 0] = 0.0

    for i in range(1, n + 1):
        for j in range(1, m + 1):
            r0, r1, r2 = R[i - 1, j], R[i, j - 1], R[i - 1, j - 1]
            R[i, j] = D[i - 1, j - 1] + softmin_three(r0, r1, r2, gamma)
    return R[n, m]

コードの解説番号をクリックすると、コード内の該当行がハイライトされます。

3-4 行目局所距離は二乗ユークリッドをブロードキャストで一括計算
x[:, None] - y[None, :] で $n \times m$ の差の行列を一度に作り、二乗して $D$ にします。原論文が二乗距離を使うのでそれに合わせており、DTW のページの abs とは違います。このため後の表で DTW と Soft-DTW の値を直接比べるときは、距離の定義が違う点を意識する必要があります。
6-9 行目最小値を引いてから exp を取る log-sum-exp
soft-min の定義どおり -gamma * log(sum(exp(-v / gamma))) を計算すると、v が 50 を超えたあたりで exp(-50 / 0.01) が 0 にアンダーフローして log(0) になります。先に最小値 m を引いておけば exp の引数のうち 1 つは必ず 0 で、合計は 1 以上に保たれます。結果は数学的に同じ値で、gamma=0.01 のような小さい値でもこの行があるから動きます。
12-13 行目gamma は正でないと定義できない
$\gamma = 0$ ではゼロ除算、負の値では soft-min が soft-max に化けて意味が変わります。DTW そのものが欲しいときは gamma を 0 にするのではなく、0.01 のように小さな正の値にするか、通常の DTW 関数を使います。
20-21 行目番兵を inf にすると soft-min が正しく境界を扱う
0 行目と 0 列目を inf にしておくと、exp(-(inf - m) / gamma) が 0 になって soft-min の和に寄与しません。境界では候補が実質 1〜2 個に減りますが、softmin_three はそのまま使えます。R[0, 0] = 0 が出発点です。
25-26 行目DTW との違いはこの 1 行の min → softmin
min(r0, r1, r2) を softmin_three(r0, r1, r2, gamma) に置き換えた以外は DTW の遷移式そのままです。soft-min は 3 候補が同じ値のとき $\min - \gamma \log 3$ を返すので、min より最大で $\gamma \log 3 \approx 1.1\gamma$ だけ小さくなります。これが 6,400 マス分積み重なるため、gamma を大きくすると Soft-DTW の値は同じ局所距離の DTW より小さくなり、負になることもあります。Soft-DTW が厳密な「距離」ではない（0 未満も取る）のはこのためです。
27 行目右下の R[n, m] が距離、パスは返さない
Soft-DTW には「1 本のワーピングパス」が存在しないので、DTW のようにパスを復元して返すことはしません。代わりに後方パス（勾配計算）で「各マスがどれだけ寄与したか」の重み行列が得られ、それがぼかされたパスに相当します。

Soft-DTW を正確に微分したい場合は、上記論文に示された後方パスを実装するか、tslearn.metrics.soft_dtw（tslearn 0.5 以降）など既存ライブラリを利用します。

DTW と Soft-DTW を比較する #

Soft-DTW が $\gamma$ によって DTW とユークリッド距離の間を補間する様子を、sin 波（A）、位相を 0.4 ずらして振幅を 5% 増やした B、cos 波の C の 3 本で確認します。

1
2
python --version  # 例: Python 3.13.0
pip install matplotlib

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
import matplotlib.pyplot as plt

np.random.seed(42)
t = np.linspace(0, 2 * np.pi, 80)

series_a = np.sin(t)
series_b = np.sin(t + 0.4) * 1.05  # 形は近いが少し位相ずれ
series_c = np.cos(t)               # 形状が大きく異なる

plt.figure(figsize=(8, 3))
plt.plot(t, series_a, label="Series A", color="black")
plt.plot(t, series_b, label="Series B", color="tab:blue")
plt.plot(t, series_c, label="Series C", color="tab:red")
plt.legend()
plt.title("比較する 3 系列")
plt.tight_layout()
plt.show()

コードの解説番号をクリックすると、コード内の該当行がハイライトされます。

4 行目0〜2π を 80 点、sin 1 周期分
80 点なので累積コスト行列は 80 × 80 = 6,400 マスで、Python の 2 重ループでも 4 つの gamma × 2 組を一瞬で計算できます。1,000 点を超えるなら tslearn の Cython 実装に切り替えます。
7 行目B は位相 0.4（約 5 点分）ずれて振幅 1.05
np.sin(t + 0.4) は A を左に 0.4 ラジアン、点の間隔 0.08 で約 5 点分ずらした波です。DTW はこのずれを伸縮で吸収しますが、振幅 1.05 の分（最大 0.05 の差）と端の対応できない区間は距離に残り、それが DTW(A, B) = 6.47 の中身です。
8 行目C は cos、位相が π/2 ずれた別の形
np.cos(t) は sin を 1/4 周期（20 点）ずらした波で、80 点のうち 20 点分は相手がありません。伸縮では埋めきれないので DTW(A, C) = 18.7 と B の約 3 倍になります。
3 行目seed(42) はこのページでは効いていない
3 本の系列はどれも乱数を使っていないので、np.random.seed(42) を外しても結果は変わりません。ノイズを足す実験に拡張するときのために残してあります。

$Soft-DTW が $\gamma$ によって DTW とユークリッド距離の間を補間する様子を、簡単な例で確認…の図$

比較用に通常の DTW を定義し、gamma を 0.01 から 5.0 まで 4 段階に変えて Soft-DTW の値を並べます。

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
def dtw_distance(x: np.ndarray, y: np.ndarray) -> float:
    n, m = len(x), len(y)
    cost = np.full((n + 1, m + 1), np.inf)
    cost[0, 0] = 0.0
    for i in range(1, n + 1):
        for j in range(1, m + 1):
            dist = abs(x[i - 1] - y[j - 1])
            cost[i, j] = dist + min(cost[i - 1, j], cost[i, j - 1], cost[i - 1, j - 1])
    return cost[n, m]

for gamma in [0.01, 0.1, 1.0, 5.0]:
    val_ab = soft_dtw(series_a, series_b, gamma=gamma)
    val_ac = soft_dtw(series_a, series_c, gamma=gamma)
    print(f"gamma={gamma:>4}: Soft-DTW(A, B)={val_ab:8.3f}, Soft-DTW(A, C)={val_ac:8.3f}")

print(f"DTW(A, B) = {dtw_distance(series_a, series_b):.3f}")
print(f"DTW(A, C) = {dtw_distance(series_a, series_c):.3f}")

コードの解説番号をクリックすると、コード内の該当行がハイライトされます。

1-9 行目比較用の DTW は絶対値距離
この dtw_distance は局所距離に abs を使うので、二乗距離の soft_dtw とは定義が違います。値が -1〜1 の範囲では二乗の方が絶対値より小さくなるため、gamma を 0.01 まで小さくしても 2 つの数値は一致せず、厳密に比べるなら abs を ** 2 に替えて同じ局所距離にそろえます。
11 行目gamma を 500 倍の範囲で振る
0.01・0.1・1.0・5.0 は対数的に広げた 4 点です。局所距離 $(x_i - y_j)^2$ は多くのマスで 0.01 以下なので、gamma が 0.01 のときは最短経路だけが残って DTW（二乗距離版）とほぼ同じ値になり、gamma が局所距離の桁を超える 1.0 以降では $-\gamma \log 3$ の項が 6,400 マス分積み重なって値が大きく動きます。手元で実行して、どの桁で値が変わり始めるかを確かめてください。
12-13 行目A–B と A–C の順序を gamma ごとに確かめる
学習に使うときは絶対値より「近い系列を近いと言えるか」が重要です。gamma が小さいうちは A–B の方が A–C より小さい順序が保たれ、gamma を 5.0 まで上げると $-\gamma \log 3$ の項が支配的になって、2 つの値の差が形の違いを反映しにくくなります。損失関数として使う gamma は、この順序が保たれる範囲（この例では 1.0 以下）から選びます。
16-17 行目DTW は Soft-DTW の上限になる
soft-min は常に min 以下なので、同じ局所距離で比べれば Soft-DTW は DTW を超えず、gamma を 0 に近づけると下から DTW に収束します。ここでは比較用 DTW が絶対値距離なので数値を直接並べることはできませんが、abs を ** 2 に替えて再実行すると、gamma=0.01 の Soft-DTW が DTW のわずかに下に来ることを確認できます。

1
2
3
4
5
6
gamma=0.01: Soft-DTW(A, B)=   6.478, Soft-DTW(A, C)=  18.742
gamma=0.1 : Soft-DTW(A, B)=   6.514, Soft-DTW(A, C)=  18.904
gamma=1.0 : Soft-DTW(A, B)=   7.721, Soft-DTW(A, C)=  21.053
gamma=5.0 : Soft-DTW(A, B)=  17.089, Soft-DTW(A, C)=  28.162
DTW(A, B) = 6.469
DTW(A, C) = 18.707

上の出力は元のノートブック環境で得たもので、$\gamma$ が小さいうちは Soft-DTW が DTW に近い値を取り、$\gamma$ を大きくするほど DTW から離れていく傾向を示しています。Soft-DTW の値は局所距離の定義（絶対値か二乗か）と gamma に強く依存し、上の soft_dtw 実装では $-\gamma \log 3$ の項が積み重なるため gamma が大きいと負の値も取ります。数値そのものより、「$\gamma \to 0$ で DTW に一致し、$\gamma$ が大きいほど 1 本の経路の長さから離れる」という向きを押さえてください。

勾配が欲しい場合 #

Soft-DTW の利点は微分可能性にあります。自分で後方パスを実装しない場合は、tslearn の soft_dtw_gradient で入力系列に対する勾配をそのまま取り出せます。

1
pip install tslearn

1
2
3
4
5
6
7
8
9
from tslearn.metrics import soft_dtw as ts_soft_dtw
from tslearn.metrics import soft_dtw_gradient

gamma = 0.5
loss = ts_soft_dtw(series_a.reshape(-1, 1), series_b.reshape(-1, 1), gamma=gamma)
grad = soft_dtw_gradient(series_a.reshape(-1, 1), series_b.reshape(-1, 1), gamma=gamma)

print("Soft-DTW (tslearn):", loss)
print("Gradient shape:", grad.shape)

コードの解説番号をクリックすると、コード内の該当行がハイライトされます。

1-2 行目tslearn は Cython 実装で前方・後方パスを提供する
tslearn.metrics.soft_dtw が距離、soft_dtw_gradient が勾配を返します。どちらも原論文の二乗ユークリッド距離を使うので、上の soft_dtw と同じ gamma なら同じ値になります。数千点の系列でも Python ループの数十倍速く、学習ループの中で毎回呼ぶ用途に向きます。
5-6 行目reshape(-1, 1) で (長さ, 次元) の 2 次元にする
tslearn は多変量時系列を前提に (n_timestamps, n_features) の形を要求します。1 次元の series_a をそのまま渡すと形のエラーになるので、reshape(-1, 1) で 80 × 1 にしています。3 軸の加速度なら (80, 3) の配列をそのまま渡せます。
6 行目勾配は第 1 引数の系列に対するもの
soft_dtw_gradient(a, b) は「a の各点を少し動かしたとき Soft-DTW がどう変わるか」を返すので、a を予測値、b を正解とすれば、a を更新する勾配降下に使えます。gamma=0.5 は、gamma を 0.01 まで小さくすると勾配が DTW の経路上だけに集中して学習が不安定になるため、少し大きめに取った値です。
9 行目勾配の shape は入力と同じ (80, 1)
grad.shape が series_a.reshape(-1, 1) と同じ (80, 1) であることを確認しています。各要素が「その時点の値を増やすと損失がどれだけ増えるか」で、PyTorch などの自動微分に組み込むときはこの配列を手で流し込むか、tslearn の SoftDTWLossPyTorch を使います。

図の再生成は python scripts/generate_timeseries_assets.py --pattern shape/005_soft_dtw を利用してください。

soft_dtw_gradient は入力系列に対する損失の勾配を返すので、ニューラルネットワークの学習ループに組み込むことができます。

まとめ #

Soft-DTW は DTW の遷移式にある $\min$ を、$\gamma$ つきの log-sum-exp（soft-min）に置き換えただけの距離で、実装上の違いも min の 1 行です。soft-min は $\min$ を超えないので、$\gamma = 0.01$ では DTW とほぼ同じ値になり、$\gamma$ を 5.0 まで上げると $-\gamma \log 3$ の項が積み重なって 1 本の経路の長さとは別の値（負にもなり得る）に変わります。損失関数に使う $\gamma$ は、形の近い組の方が小さいという順序が保たれる範囲から選びます。微分できるので tslearn.metrics.soft_dtw_gradient や SoftDTWLossPyTorch を通してニューラルネットの損失関数にでき、「タイミングが少しずれても形が合っていれば損失が小さい」学習が可能になります。

動的時間伸縮 (DTW) — 時間軸のずれを許容して距離を最小化
DTWとDDTW — 微分値を使ってDTWの弱点を補う手法
tsfreshで特徴量を抽出 — 時系列から自動的に大量の特徴量を生成