r - 高效的多维动态时间扭曲实现

Question

以下文献解释了如何计算两个时间序列的多维动态时间扭曲:

 library(dtw)
 x<- cbind(1:10,1)
 y<- cbind(11:15,2)
 cxdist <-dist(x,y,method="euclidean")
 dtw(cxdist)$distance

事实上，它首先计算交叉距离矩阵，然后将其用作 dtw 函数的输入。

我想在具有相当大图像的图像分类中使用多维动态时间扭曲。图像值存储在数据框中，如下所示:

 inDf <- data.frame(matrix(rnorm(60), ncol = 6))
 colnames(inDf) <- c('var1t1','var2t1','var1t2','var2t2','var1t3','var2t3')

在此示例中，有两个变量(var1 和 var2)被观察了 3 次。

问题是如何在计算强度方面尽可能高效地得到dtw距离矩阵？

这里有一些想法:- 迭代输入图像矩阵的每个值，将向量 reshape 为矩阵，以便能够计算交叉距离，然后计算 dtw 距离并将其存储在专用矩阵中。这无疑是计算最密集的解决方案

Answer 1

在处理密集计算时，考虑 Rcpp 包总是有意义的。如果你想更快地得到欧式距离的距离矩阵，可以实现相应的Rcpp函数:

library(Rcpp)
library(inline)

# Rcpp function for euclidean distance
fastdist <- cxxfunction(signature(x="matrix", y="matrix"), plugin="Rcpp",
body='
  Rcpp::NumericMatrix dx(x);
  Rcpp::NumericMatrix dy(y);

  const int N = dx.nrow();
  const int M = dy.nrow();

  Rcpp::NumericMatrix res(N, M);

  for(int i=0; i<N; i++){
    for(int j=0; j<M; j++){
      res(i,j) = sqrt(sum((dx(i,_)-dy(j,_))*(dx(i,_)-dy(j,_))));
    }
  }

  return res;
')

它使用Rcpp语法sugar以使代码更加紧凑和可读。然而，有时最好使用包装函数来检查类型、强制转换等。这不是必需的 - 您可以直接调用 fastdist 。但是，无论如何，包装器可以如下所示:

# Wrapper R function
fast.dist <- function(x, y){
  stopifnot(class(x) %in% c("data.frame","matrix") &
            class(y) %in% c("data.frame","matrix") &
            ncol(x)==ncol(y))

  fastdist(as.matrix(x), as.matrix(y))
}

现在我们可以转向文献示例。

library(dtw)

# EXAMPLE 1
x<- cbind(1:10,1)
y<- cbind(11:15,2)
# Check results
all.equal(fast.dist(x,y), dist(x,y,method="euclidean"), check.attributes=F)
# [1] "target is matrix, current is crossdist"
all.equal(fast.dist(x,y), matrix(dist(x,y,method="euclidean"), ncol=nrow(y)))
# [1] TRUE

请注意，dist 返回类crossdist 的结果。因此，为了进行比较，应将其强制转换为矩阵。

现在你的主要问题 - 我们首先生成数据:

# EXAMPLE 2
set.seed(1234)
N <- 100
inDf <- data.frame(matrix(rnorm(6*N), ncol = 6))
colnames(inDf) <- c('var1t1','var2t1','var1t2','var2t2','var1t3','var2t3')

# Extracting variables
var1 <- inDf[,c("var1t1","var1t2","var1t3")]
var2 <- inDf[,c("var2t1","var2t2","var2t3")]

我不完全确定你的数据结构，但无论如何你总是可以根据你的需要准备变量。

比较和基准测试:

library(rbenchmark)

all.equal(fast.dist(var1,var2), matrix(dist(var1,var2), ncol=N))
# [1] TRUE
benchmark(fast.dist(var1,var2), dist(var1,var2), order="relative")[,1:4]
#                    test replications elapsed relative
# 1 fast.dist(var1, var2)          100   0.081    1.000
# 2      dist(var1, var2)          100   0.246    3.037

在这种情况下，

fast.dist 大约比 dist 快 3 倍。然而，当 N 不断增长时，相对加速将会下降。

另请注意，正如评论中提到的，dtw 可以自行计算距离矩阵。尽管如此，预先计算距离矩阵会更有效。请参阅下面的快速测试:

cxdist <- fast.dist(var1,var2)
benchmark(dtw(cxdist)$distance, dtw(var1,var2)$distance, order="relative")[,1:4]
#                       test replications elapsed relative
# 1     dtw(cxdist)$distance          100   0.476    1.000
# 2 dtw(var1, var2)$distance          100   0.736    1.546

另外，如果您只对 $distance 感兴趣，您可以将 distance.only=T 传递给 dtw() - 它给出一些加速。