Calculating residual DTW distance for a subset of the alignment (Dynamic Time Warping)

Question

My question is: "How do I find the cost for each step taken using the dtw R package?" My work is below, so if you see an easier solution please point it out!

DTW is very effective at aligning my data. Often the supplied time series contain multiple relevant features which I later consider independently.

The normalized residual distance between the full time series is computed and returned automatically. I am interested in the normalized residual distance of only a subset of the curves. In trying to manually calculate this "subset-distance" from the cost matrix to avoid the computation required to re-align the subset regions.

library(dtw)

mat = structure(c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,             0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.01, 0.01, 0.02, 0.03, 0.04, 0.06,             0.09, 0.11, 0.13, 0.16, 0.18, 0.2, 0.22, 0.24, 0.24, 0.22, 0.22,             0.22, 0.22, 0.21, 0.2, 0.19, 0.2, 0.23, 0.29, 0.34, 0.41, 0.51,             0.62, 0.73, 0.82, 0.9, 0.95, 1, 1, 1, 0.92, 0.92, 0.89, 0.89,             0.84, 0.79, 0.7, 0.53, 0.37, 0.23, 0.17, 0.13, 0.11, 0.09, 0.08,             0.07, 0.07, 0.07, 0.07, 0.07, 0.07, 0.08, 0.08, 0.08, 0.09, 0.1,             0.13, 0.15, 0.19, 0.22, 0.27, 0.29, 0.34, 0.35, 0.36, 0.35, 0.38,             0.37, 0.37, 0.32, 0.3, 0.26, 0.24, 0.21, 0.19, 0.17, 0.15, 0.14,             0.12, 0.1, 0.09, 0.09, 0.08, 0.08, 0.07, 0.07, 0.07, 0.07, 0.06,             0.06, 0.06, 0.05, 0.05, 0.05, 0.05, 0.04, 0.04, 0.04, 0.04, 0.04,             0.04, 0.04, 0.04, 0.04, 0.04, 0.04, 0.04, 0.03, 0.04, 0.04, 0.04,             0.03, 0.03, 0.03, 0.04, 0.04, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,             0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.01, 0.01, 0.01, 0.02,             0.02, 0.03, 0.04, 0.05, 0.06, 0.07, 0.08, 0.1, 0.12, 0.12, 0.13,             0.14, 0.15, 0.17, 0.19, 0.2, 0.21, 0.22, 0.24, 0.23, 0.24, 0.26,             0.3, 0.32, 0.33, 0.35, 0.39, 0.44, 0.49, 0.55, 0.61, 0.67, 0.71,             0.76, 0.83, 0.9, 0.97, 1, 0.99, 0.86, 0.68, 0.5, 0.41, 0.33,             0.28, 0.23, 0.2, 0.17, 0.15, 0.13, 0.12, 0.1, 0.1, 0.1, 0.11,             0.11, 0.11, 0.11, 0.11, 0.11, 0.11, 0.13, 0.15, 0.17, 0.18, 0.2,             0.21, 0.24, 0.25, 0.28, 0.29, 0.32, 0.35, 0.36, 0.34, 0.32, 0.3,             0.3, 0.28, 0.26, 0.23, 0.22, 0.19, 0.17, 0.15, 0.14, 0.12, 0.1,             0.09, 0.09, 0.08, 0.08, 0.07, 0.07, 0.07, 0.06, 0.06, 0.05, 0.05,             0.05, 0.05, 0.05, 0.05, 0.04, 0.04, 0.04, 0.04, 0.04, 0.04, 0.04,             0.04, 0.04, 0.04, 0.04, 0.04, 0.04, 0.04, 0.04, 0.04, 0.04, 0.04), .Dim = c(149L, 2L))

.

tw = dtw(mat[,1], mat[,2], keep.internals = T)

EDIT The below calculations are wrong.

tw$distance
# 2.19
tw$normalizedDistance
# 0.007

# Manually Calculated
sum(tw$costMatrix[ cbind(tw$index1, tw$index2) ])
# 219
sum(tw$costMatrix[ cbind(tw$index1, tw$index2) ]) / sum(tw$M, tw$N)
# 0.74

Above are values for the entire curve, I am looking for the distance of individual regions. In an attempt to get this I try to recalculate the distance from the localCostMatrix. For the case of the step pattern having equal weights it works out.

tw = dtw(mat[,1], mat[,2], keep.internals = T, step.pattern = asymmetricP05)

tw$distance
# 1.42
tw$normalizedDistance
# 0.009
sum(tw$costMatrix[ cbind(tw$index1, tw$index2) ])
# 146
sum(tw$costMatrix[ cbind(tw$index1, tw$index2) ]) / sum(tw$M, tw$N)
# 0.49

The asymetric case does not work. I think my missing link lies in the formula for distance given in the reference below. It depends on the distance of each aligned point but also a step-weight m_phi. How can I retrieve m_phi?

http://cran.r-project.org/web/packages/dtw/vignettes/dtw.pdf http://irit.fr/page-perso/Julien.Pinquier/Docs/TP_MABS/res/dtw-sakoe-chiba78.pdf

plot(tw, type="threeway")

Answer 1

I asked the question above. I had misunderstood tw$costMatrix above. Rather than a pairwise cost it is the total accrued cost as we traversed the path. (tw$costMatrix is $d_{phi}$ in the equation below, and thus its a sum of $d$ at every k along the path $(X,Y)$)

$d$ below is tw$localCostMatrix and is the pairwise distance between every point. In the default case above $d$ is the Euclidean distance.

Target values:

tw = dtw(mat[,1], mat[,2], keep.internals = T, step.pattern = asymmetricP05)

tw$distance
# 1.418333

tw$normalizedDistance
# 0.009519016

Correct calculations:

d.phi = tw$costMatrix[ cbind(tw$index1, tw$index2) ]

tw$distance is the total accrued distance and thus, the last value along the path.

d.phi[length(d.phi)]
# 1.418333

tw$normalizedDistance is the total accrued distance independent of path length asymmetricP05 normalizes to the length of the query.

d.phi[length(d.phi)] / tw$M
# 0.009519016

With this we can find the pre-warping distance and post-warping distance of subsets of the alignment.

i = 0:80
i = 80:120 # Indices we want distance for

d.phi.sub = tw$costMatrix[ cbind(tw$index1[i], tw$index2[i]) ]

(d.phi.sub[length(d.phi.sub)] - d.phi.sub[1])
# i.1 = 0.4383333
# i.2 = 0.79

(d.phi.sub[length(d.phi.sub)] - d.phi.sub[1]) / length(i)
# i.1 = 0.005411523
# i.2 = 0.01926829

Outstanding questions: Finding $m_{phi}$ was unnecessary for my goal because the costMatrix accounts for it but it would still be interesting. tw$costMatrix is the sum along the path - this can't be true for any arbitrary path, so how do values that fall away from the optimal warp path get filled? How can we retrieve the path that generated them?

Thanks everyone.