Incremental Dynamic Time Warping

Introduction to Dynamic Time Warping: dtw()

The distance measure DTW calculates the minimal costs of the shortest non-linear alignment of two time series \(Q\) (the query time series) and \(C\) (the candidate time series), conditional to

• The first element of \(Q\) is aligned to the first element of \(C\)
• The last element of \(Q\) is aligned to the last element of \(C\)
• A 1:k alignment is applied (unlike to the Euclidean distance calculation). That is, one observation of \(Q\) can be applied to k-many of \(C\), and vice versa.

Remark: The DTW distance measure is not a distance metric in the conventional sense, since it does not fulfill the triangle inequality.

In the following we show by aid of toy example how to use the functions of IncDTW and what they are doing. We start with the basic functions. Suppose we have two time series and we need to calculate the DTW distance of these:

set.seed(1090)
Q <- sin(1:10)#+rnorm(20)
C <- sin(-2:10)#+rnorm(15)
tmp <- IncDTW::dtw(Q = Q, C = C, return_diffM = TRUE, return_wp = TRUE,
                   return_QC = TRUE, return_cm = TRUE, return_diffp = TRUE)
names(tmp)

##  [1] "gcm"                 "dm"                  "diffp"               "ii"                  "jj"                  "wp"                  "distance"            "normalized_distance" "cm"                  "diffM"               "Q"                   "C"

Another toy example:

set.seed(1090)
Q <- c(1, 1, 2, 3, 2, 0)
C <- c(0, 1, 1, 2, 3, 2, 1)

tmp <- IncDTW::dtw(Q = Q, C = C, return_diffM = TRUE, return_wp = TRUE,
                   return_QC = TRUE, return_cm = TRUE, return_diffp = TRUE)
names(tmp)

##  [1] "gcm"                 "dm"                  "diffp"               "ii"                  "jj"                  "wp"                  "distance"            "normalized_distance" "cm"                  "diffM"               "Q"                   "C"

The matrix of differences diffM simply stores the differences of the two time series \(Q\) and \(C\).

tmp$diffM

##      [,1] [,2] [,3] [,4] [,5] [,6] [,7]
## [1,]    1    0    0   -1   -2   -1    0
## [2,]    1    0    0   -1   -2   -1    0
## [3,]    2    1    1    0   -1    0    1
## [4,]    3    2    2    1    0    1    2
## [5,]    2    1    1    0   -1    0    1
## [6,]    0   -1   -1   -2   -3   -2   -1

The global cost matrix is calculated by walking through the matrix of absolute differences abs(diffM) from top left to bottom right applying the following rule:

cm <- abs(diffM) # cost matrix
gcm <- cm # initialize the global cost matrix
for(j in 1:m){
   for(i in 1:n){
      gcm[ , 1] <- cm[ ,1]
      gcm[1,  ] <- cm[1, ]
      gcm[i, j] <- cm[i, j] + min(c(gcm[i-1, j-1], 
                                    gcm[i  , j-1],
                                    gcm[i-1, j  ]))}}

The direction matrix dm is calculated simultaneously to the global cost matrix and stores for each position where to go next to make the cheapest step.

dm <- matrix(NA, ncol=length(C), nrow = length(Q))
dm[1,  ] <- 3
dm[ , 1] <- 2
dm[1, 1] <- NA

for(j in 1:m){
   for(i in 1:n){
      min_index <- which.min(c(gcm[i-1, j-1], 
                               gcm[i  , j-1], 
                               gcm[i-1, j  ]))
      if( min_index == 1){
         dm[i, j] <- 1
      } else if( min_index == 2){
         dm[i, j] <- 2
      } else if( min_index == 3){
         dm[i, j] <- 3
      } } }

In case the cost matrix, or the matrix of differences is already in storage, and you want to save calculation time, you can use these matrices also as input:

IncDTW::dtw(Q = tmp$diffM, C = "diffM")
IncDTW::dtw(Q = abs(tmp$diffM), C = "cm")

Incremental Calculation of DTW: idtw()

The main intention of idtw() is to compare the DTW distance for different subsets of two time series. That is to compare the DTW(\(Q\), \(C_0\)) with DTW(\(Q\), \(C_1\)) where \(C_1\) could be the extension of \(C_0\) for couple of observations.

Another possible application is a live data stream and you want to update the DTW calculation each time a set of values of \(Q\) are observed. Of course this depends on the application and certainly does not make sense in any case, since a time warp of two time series where the lengths differ too much is questionable. A reasonable range for the ratio n/m depends on the application and purposes of the data analysis.

Assume we already have calculated the DTW of \(C\) and \(Q\), when an update of \(C\) provides new observations:

set.seed(1030)
C0 <- cumsum(rnorm(1000))
Q <- cumsum(rnorm(800))
tmp0 <- IncDTW::dtw(Q = Q, C = C0)
gcm0 <- tmp0$gcm
dm0 <- tmp0$dm

To update the DTW distance for new observations of \(C\) we recycle as much from previous results as possible to keep calculation time low:

C_new <- cumsum(rnorm(10))
C_update <- c(C0, C_new) 

# result from incremental calculation
res_inc <- IncDTW::idtw(Q = Q, C = C_update, newObs = C_new, 
                    gcm = tmp0$gcm, dm = tmp0$dm)

Finally we compare the results from the incremental calculation with the one from scratch:

# result from scratch
res_scratch <- IncDTW::dtw(Q = Q, C = C_update) 
sapply(names(res_inc), function(x){identical(res_inc[[x]], res_scratch[[x]])})

## distance      gcm       dm 
##     TRUE     TRUE     TRUE

Decremental Dynamic Time Warping: dec_dm()

Now that we have understood the incremental approach and tested the calculation time and correctness of results, we want to introduce another concept, which is just the reverse: Instead of adding new observations to the time series x, in this section we introduce an easy and fast way of reducing most current observations, conditional to already existing calculation results including all observations of \(C\).

Of course this is only interesting if you are interested not only in the DTW distance, but also want to know about the warping path. Meanwhile the DTW distance could be directly taken from the global cost matrix from previous calculations (by selecting the value of the last row and last but k column, k is the number of reduced observations of \(C\)), we cannot update or modify the previous warping path but need to calculate it again from scratch. This is required, since adding or neglecting observations causes an updated alignment of the two updated time series, and in general the new alignment could be completely different from the original.

set.seed(1150)
Q <- cos(1:100)
C <- cumsum(rnorm(80))
Ndec <- 4
# the ordinary calculation
result_base <- IncDTW::dtw(Q=Q, C=C) 
gcm0 <- result_base$gcm

To update the DTW distance by decreasing x for some observations (say 4), again we recycle as much from previous results as possible to keep calculation time low:

# the ordinary calculation without the last 4 observations
result_decr1 <- IncDTW::dtw(Q=Q, C=C[1:(length(C) - Ndec)], 
                            return_wp = TRUE)
gcm1 <- result_decr1$gcm

# the decremental step: reduce C for 4 observation
result_decr2 <- IncDTW::dec_dm(result_base$dm, Ndec = Ndec) 

# compare the results: ii, jj, wp and the gcm of 
# result_decr1 (conventional) and result_decr2 (recycling previous results)
comparison_0 <- c(
identical(result_decr1$ii, result_decr2$ii),
identical(result_decr1$jj, result_decr2$jj),
identical(result_decr1$wp, result_decr2$wp))
comparison_0

## [1] TRUE TRUE TRUE

Vector-based Implementation

The matrix based implementation is perfect for a detailed analysis of the alignment of two time series. However for calculating the distances of many time series it has one major drawback, the computation time. If the data analysis task is to cluster or classify a database of time series a distance matrix of pairwise distances is required. Obtaining such a distance matrix is computationally expensive, and so the classic DTW implementation maybe is considered as too slow. The vector-based implementation offers a solution which is much faster than the matrix based implementation since no matrices are allocated. The computation principle is the same as for the matrix-based method, however instead of allocating the two cost matrix, global cost matrix and direction matrix only two vectors are allocated. To achieve the DTW distance between the time series Q and C the calculation of the j-th column of the global cost matrix solely depends on the values of the former column and the respective distances of the time series C and the j-th entry of Q. Since there is no dependence on the second last column the algorithm overwrites this column. This principle is implemented in the functions dtw2vec() and the incremental equivalent idtw2vec().

Q <- rnorm(100)
C <- rnorm(120)
IncDTW::dtw2vec(Q, C, ws = 30)$distance

## [1] 90.13429

IncDTW::dtw(Q, C, ws = 30)$distance

## [1] 90.13429

We can also perform the incremental calculation with the vector based implementation with idtw2vec(), which can distinguish between an initial calculation and the incremental. For this example we simulate time series as univariate random walks.

set.seed(1210)
Q <- cumsum(rnorm(800))
C_initial <- Q[1:750] + rnorm(750)
resV_init <- IncDTW::idtw2vec(Q = Q, newObs = C_initial, gcm_lc = NULL)

C_newObs <- Q[751:800] + rnorm(50)
resV_inc <- IncDTW::idtw2vec(Q = Q, newObs = C_newObs, 
                             gcm_lc = resV_init$gcm_lc_new)

Also we compare the DTW distance of the incremental calculation idtw2vec() with the one from scratch dtw2vec() and see that the distance measures are equal.

C_update <- c(C_initial, C_newObs)
resV_scr <- IncDTW::dtw2vec(Q = Q, C = C_update)
resV_scr$distance - resV_inc$distance

## [1] 0

resV_init <- idtw2vec(Q = Q, newObs = C_initial, gcm_lc = NULL)
resV_inc <-  idtw2vec(Q = Q, newObs = C_newObs, nC = length(C_initial),
                      gcm_lc = resV_init$gcm_lc_new,
                      gcm_lr = resV_init$gcm_lr_new)

Q_newObs <-cumsum(rnorm(10))
resV_inc2 <- idtw2vec(Q = C_update, newObs = Q_newObs, nC = length(Q),
                      gcm_lc = resV_inc$gcm_lr_new,
                      gcm_lr = resV_inc$gcm_lc_new)

res_scratch <- dtw2vec(Q = c(Q, Q_newObs), C = c(C_initial, C_newObs))
res_scratch$distance - resV_inc2$distance

## [1] 0

res_scratch$normalized_distance - resV_inc2$normalized_distance

## [1] 0

Q <- sin(1:50)
C <- cos(1:1000)
tmp <- IncDTW::dtw(Q = Q, C = C)
IncDTW::dtw_partial(tmp, partial_Q = FALSE, partial_C = TRUE)

## $rangeQ
## [1]  1 50
## 
## $rangeC
## [1]  1 61
## 
## $normalized_distance
## [1] 0.08396951

Q <- sin(1:100)
C <- c(rnorm(50), Q)
tmp <- IncDTW::dtw(Q = rev(Q), C = rev(C))
dtw_partial(tmp, partial_Q = FALSE, partial_C = TRUE, reverse = TRUE)

## $rangeQ
## [1]   1 100
## 
## $rangeC
## [1]  51 150
## 
## $normalized_distance
## [1] 0

Q <- sin(seq(1, 50, length.out = 200))
C <- cos(seq(1, 5, length.out = 20))
tmpi <- idtw2vec(Q = Q, newObs = C)
dtw_partial(tmpi, partial_Q = TRUE, partial_C = FALSE)

## $rangeQ
## [1]  1 24
## 
## $rangeC
## [1]  1 20
## 
## $normalized_distance
## [1] 0.08470878

C_new <- cos(seq(5, 10, length.out = 20))
tmp1 <- idtw2vec(Q = Q, newObs = C_new, gcm_lc = tmpi$gcm_lc_new, 
                 gcm_lr = tmpi$gcm_lr_new, nC = length(C))
dtw_partial(tmp1, partial_Q = TRUE, partial_C = FALSE)

## $rangeQ
## [1]  1 44
## 
## $rangeC
## [1]  1 40
## 
## $normalized_distance
## [1] 0.0643759

Analysing a Database of Time Series

There are three more functions in IncDTW that are highly useful when comparing a set of time series:

• dtw_dismat() calculates distance matrices of pairwise DTW distances, so comparing many time series with many. Especially useful for clustering or classification of a set of time series.
• dtw_disvec() calculates vectors of DTW distances, so comparing one time series with many distance matrices of pairwise DTW distances. Especially useful for k-nearest neighbor searches.
• DBA() calculates the DTW Barycenter Average of a set of time series that are non-linearly aligned. Especially useful to define a cluster representative.