📂알고리즘

동적 시간 워핑

개요

두 시계열 $X = \left\{ x_{1}, \cdots , x_{n} \right\}$, $Y = \left\{ y_{1}, \cdots , y_{m} \right\}$ 이 주어져 있다고 하자. 여기서 두 시계열은 길이가 같을 필요가 없어 $n \ne m$ 이어도 상관 없다. 시계열을 이루는 성분들은 공통된 거리공간에 속하는 원소로써, $d (x, y)$ 와 같이 두 성분 사이의 거리를 정의할 수 있다고 가정한다.

이 두 시계열의 유사도^similarity를 측정 혹은 두 시계열을 정렬^alignment하는 알고리즘을 동적 시간 워핑^{Dynamic Time Warping, DTW}이라고 한다.

알고리즘 ¹

근본적으로 동적 프로그래밍이기 때문에 니들맨-분쉬 알고리즘 같은 알고리즘과 매우 유사하다. 각 지점을 한번씩만 비교하는 행렬을 미리 작성한 후 그 성분을 거꾸로 타고가는 방식으로 최적의 경로를 찾는다.

전역비용행렬

이렇게 얻어지는 행렬 $D$ 의 우측하단 $D_{nm}$ 에서 좌측상단 $D_{11}$ 까지 이어지는 패스를 인덱스의 시퀀스 $P$ 와 같이 나타낸다고 할 때, 다음과 같이 $P$ 의 성분을 따라가며 얻어지는 $D$ 의 합이 작을수록 두 시계열 $X$ 와 $Y$ 는 유사하다고 할 수 있다. $$ C_{D} (P) = \sum_{(i, j) \in P} D_{ij} $$ 이 $C_{D}$ 를 비용 함수^{cost function}라 한다.

최적워핑경로

$C_{D} (P)$ 를 최소화하는 경로를 특별히 최적워핑경로^{optimal warping path} $P^{\ast}$ 라 한다. 우리는 이 최적워핑경로에 관심이 있다.

기본적으로 비용이 작아지는 방향으로 한칸 내지 대각선으로 이동할 뿐이다. 행렬의 한쪽 벽에 닿으면 그 떄부터는 수직방향으로만 이동하고, 결과적으로 $P_{1} = (1, 1)$, $P_{|P|} = (n, m)$ 을 만족한다. 이를 경계조건^{boundary condition}이라 한다. 또한, 한칸씩 이동하거나 대각선으로 이동하는 것만 허용하기 때문에 인덱스를 두 칸 이상 건너뛰는 경우는 없으며, 이를 연속조건^{continuity condition}이라 한다. 그와 동시에 $i_{1} \le \cdots \le i_{|P|}$, $j_{1} \le \cdots \le j_{|P|}$ 를 만족해서 인덱스가 증가하는 방향으로만 이동하게 되는데, 이를 단조조건^{monotonicity condition}이라 한다.

설명

tx = 0:0.5:9
ty = -0.5:0.4:8.5
x_ = sin.(tx)
y_ = sin.(ty) .+ 0.01

예로써 위와 같은 두 개의 사인파가 주어져있다고 하자. 주파수와 진폭이 묘하게 다르지만, 우리는 이들이 근본적으로 유사하다는 직관이 서있다. 그러나 이들을 단순히 지점마다의 오차로 접근하면 전혀 닮지 않은 두 시계열이 되고, 애초에 배열의 길이도 첫번째는 19, 두번째는 23이라서 단순하게 비교하는것부터가 어렵다.

이에 DTW를 사용해서 양쪽 인덱스 간의 매핑을 수행하고 그 최적경로를 찾아보면 위와 같다.

julia> z = warpingpath(D);

julia> stack(z)
2×24 Matrix{Int64}:
 1  1  2  2  3  4  4  5  6   7   8   9  10  10  11  12  13  14  14  15  16  17  18  19
 1  2  3  4  5  6  7  8  9  10  11  12  13  14  15  16  17  18  19  20  21  22  23  23

새로운 매핑에 따라 두 시계열을 재배열하면 위와 같이 서로 닮은 부분이 맞춰진다. 파란 실선이 계단모양으로 꺾이는 건 중복되는 인덱스가 있기 때문이다. 이런 부분들은 선형 보간으로 잇는 등의 후보정을 해줄 수 있다.

전체 코드

using Plots

function dtw(X, Y, d = (a, b) -> abs(a - b))
    n = length(X)
    m = length(Y)

    D = zeros(n + 1, m + 1)
    D[2:end, 1] .= Inf
    D[1, 2:end] .= Inf

    for i in 1:n
        for j in 1:m
            cost = d(X[i], Y[j])
            D[i + 1, j + 1] = cost + min(D[i, j + 1], D[i + 1, j], D[i, j])
        end
    end

    return D[2:end, 2:end]
end

function warpingpath(D)
    n, m = size(D)
    i, j = n, m
    P = [[i, j]]

    while i > 1 || j > 1
        if i == 1
            j -= 1
        elseif j == 1
            i -= 1
        else
            min_cost, idx = findmin([D[i - 1, j], D[i, j - 1], D[i - 1, j - 1]])
            if idx == 1
                i -= 1
            elseif idx == 2
                j -= 1
            else
                i -= 1
                j -= 1
            end
        end
        push!(P, [i, j])
    end

    return reverse(P)
end
cost(D, P) = sum(D[i, j] for (i, j) in P)

tx = 0:0.5:9
ty = -0.5:0.4:8.5
x_ = sin.(tx)
y_ = sin.(ty) .+ 0.01
D = dtw(x_, y_)

plot(tx, x_)
plot!(ty .+ 0.5, y_)

z = warpingpath(D)
cost(D, z)
stack(z)

heatmap(D, yflip = true, xmirror = true, framestyle = :box, colorbar = false, size = [460, 380], xticks = [1, length(y_)], yticks = [1, length(x_)], xlims = (0.5, length(y_)+0.5), ylims = (0.5, length(x_)+0.5))
plot!(last.(z), first.(z), lw = 2, c = :white, legend = :none, xlims = (0.5, length(y_)+0.5), ylims = (0.5, length(x_)+0.5))

plot(x_[first.(z)])
plot!(y_[last.(z)])