r - 计算时间序列矩阵的近似熵

Question

引入了近似熵来量化时间序列中波动的规律性和不可预测性。

功能

approx_entropy(ts, edim = 2, r = 0.2*sd(ts), elag = 1)

来自包裹 pracma ，计算时间序列的近似熵 ts .

我有一个时间序列矩阵(每行一个系列) mat我会估计每个人的近似熵，将结果存储在一个向量中。例如:

library(pracma)

N<-nrow(mat)
r<-matrix(0, nrow = N, ncol = 1)
for (i in 1:N){
     r[i]<-approx_entropy(mat[i,], edim = 2, r = 0.2*sd(mat[i,]), elag = 1)
}

但是，如果 N大这个代码可能太慢了。建议提速吗？谢谢!

Answer 1

我还要说并行化，因为应用函数显然没有带来任何优化。

我尝试了 approx_entropy() 函数:

申请

lapply

ParApply

foreach(来自@Mankind_008)

data.table 和 ParApply 的组合

ParApply 似乎是 比其他 2 个并行函数稍微高效 。

因为我没有得到与@Mankind_008 相同的时间，所以我用 microbenchmark 检查了它们。
这些是 10 次运行的结果:

Unit: seconds
expr      min       lq     mean   median       uq      max neval cld
forloop     4.067308 4.073604 4.117732 4.097188 4.141059 4.244261    10   b
apply       4.054737 4.092990 4.147449 4.139112 4.188664 4.246629    10   b
lapply      4.060242 4.068953 4.229806 4.105213 4.198261 4.873245    10   b
par         2.384788 2.397440 2.646881 2.456174 2.558573 4.134668    10   a 
parApply    2.289028 2.300088 2.371244 2.347408 2.369721 2.675570    10   a 
DT_parApply 2.294298 2.322774 2.387722 2.354507 2.466575 2.515141    10   a

完整代码:

library(pracma)
library(foreach)
library(parallel)
library(doParallel)


# dummy random time series data
ts <- rnorm(56)
mat <- matrix(rep(ts,100), nrow = 100, ncol = 100)
r <- matrix(0, nrow = nrow(mat), ncol = 1)      


## For Loop
for (i in 1:nrow(mat)){     
  r[i]<-approx_entropy(mat[i,], edim = 2, r = 0.2*sd(mat[i,]), elag = 1)
}

## Apply
r1 = apply(mat, 1, FUN = function(x) approx_entropy(x, edim = 2, r = 0.2*sd(x), elag = 1))

## Lapply
r2 = lapply(1:nrow(mat), FUN = function(x) approx_entropy(mat[x,], edim = 2, r = 0.2*sd(mat[x,]), elag = 1))

## ParApply
cl <- makeCluster(getOption("cl.cores", 3))
r3 = parApply(cl = cl, mat, 1, FUN = function(x) {
  library(pracma); 
  approx_entropy(x, edim = 2, r = 0.2*sd(x), elag = 1)
})
stopCluster(cl)

## Foreach
registerDoParallel(cl = 3, cores = 2)  
r4 <- foreach(i = 1:nrow(mat), .combine = rbind)  %dopar%  
  pracma::approx_entropy(mat[i,], edim = 2, r = 0.2*sd(mat[i,]), elag = 1)
stopImplicitCluster()  

## Data.table
library(data.table)
mDT = as.data.table(mat)
cl <- makeCluster(getOption("cl.cores", 3))
r5 = parApply(cl = cl, mDT, 1, FUN = function(x) {
  library(pracma); 
  approx_entropy(x, edim = 2, r = 0.2*sd(x), elag = 1)
})
stopCluster(cl)


## All equal Tests
all.equal(as.numeric(r), r1)
all.equal(r1, as.numeric(do.call(rbind, r2)))
all.equal(r1, r3)
all.equal(r1, as.numeric(r4))
all.equal(r1, r5)


## Benchmark
library(microbenchmark)
mc <- microbenchmark(times=10,
  forloop = {
    for (i in 1:nrow(mat)){
      r[i]<-approx_entropy(mat[i,], edim = 2, r = 0.2*sd(mat[i,]), elag = 1)
    }
  },
  apply = {
    r1 = apply(mat, 1, FUN = function(x) approx_entropy(x, edim = 2, r = 0.2*sd(x), elag = 1))
  },
  lapply = {
    r1 = lapply(1:nrow(mat), FUN = function(x) approx_entropy(mat[x,], edim = 2, r = 0.2*sd(mat[x,]), elag = 1))
  },
  par = {
    registerDoParallel(cl = 3, cores = 2)  
    r_par <- foreach(i = 1:nrow(mat), .combine = rbind)  %dopar%  
      pracma::approx_entropy(mat[i,], edim = 2, r = 0.2*sd(mat[i,]), elag = 1)
    stopImplicitCluster()
  }, 
  parApply = {
    cl <- makeCluster(getOption("cl.cores", 3))
    r3 = parApply(cl = cl, mat, 1, FUN = function(x) {
      library(pracma); 
      approx_entropy(x, edim = 2, r = 0.2*sd(x), elag = 1)
    })
    stopCluster(cl)
  },
  DT_parApply = {
    mDT = as.data.table(mat)
    cl <- makeCluster(getOption("cl.cores", 3))
    r5 = parApply(cl = cl, mDT, 1, FUN = function(x) {
      library(pracma); 
      approx_entropy(x, edim = 2, r = 0.2*sd(x), elag = 1)
    })
    stopCluster(cl)
  }
)

## Results
mc
Unit: seconds
expr      min       lq     mean   median       uq      max neval cld
forloop 4.067308 4.073604 4.117732 4.097188 4.141059 4.244261    10   b
apply 4.054737 4.092990 4.147449 4.139112 4.188664 4.246629    10   b
lapply 4.060242 4.068953 4.229806 4.105213 4.198261 4.873245    10   b
par 2.384788 2.397440 2.646881 2.456174 2.558573 4.134668    10  a 
parApply 2.289028 2.300088 2.371244 2.347408 2.369721 2.675570    10  a 
DT_parApply 2.294298 2.322774 2.387722 2.354507 2.466575 2.515141    10  a 

## Time-Boxplot
plot(mc)

内核数量 也会影响速度，更多并不总是更快，因为在某些时候，发送给所有工作人员的开销会消耗一些获得的性能。我用 2 到 7 个内核对 ParApply 函数进行了基准测试，在我的机器上，用 3/4 个内核运行该函数似乎是最好的选择，尽管偏差不大。

mc
Unit: seconds
expr      min       lq     mean   median       uq      max neval  cld
parApply_2 2.670257 2.688115 2.699522 2.694527 2.714293 2.740149    10   c 
parApply_3 2.312629 2.366021 2.411022 2.399599 2.464568 2.535220    10 a   
parApply_4 2.358165 2.405190 2.444848 2.433657 2.485083 2.568679    10 a   
parApply_5 2.504144 2.523215 2.546810 2.536405 2.558630 2.646244    10  b  
parApply_6 2.687758 2.725502 2.761400 2.747263 2.766318 2.969402    10   c 
parApply_7 2.906236 2.912945 2.948692 2.919704 2.988599 3.053362    10    d

完整代码:

## Benchmark N-Cores
library(microbenchmark)
mc <- microbenchmark(times=10,
                     parApply_2 = {
                       cl <- makeCluster(getOption("cl.cores", 2))
                       r3 = parApply(cl = cl, mat, 1, FUN = function(x) {
                         library(pracma); 
                         approx_entropy(x, edim = 2, r = 0.2*sd(x), elag = 1)
                       })
                       stopCluster(cl)
                     },
                     parApply_3 = {
                       cl <- makeCluster(getOption("cl.cores", 3))
                       r3 = parApply(cl = cl, mat, 1, FUN = function(x) {
                         library(pracma); 
                         approx_entropy(x, edim = 2, r = 0.2*sd(x), elag = 1)
                       })
                       stopCluster(cl)
                     },
                     parApply_4 = {
                       cl <- makeCluster(getOption("cl.cores", 4))
                       r3 = parApply(cl = cl, mat, 1, FUN = function(x) {
                         library(pracma); 
                         approx_entropy(x, edim = 2, r = 0.2*sd(x), elag = 1)
                       })
                       stopCluster(cl)
                     },
                     parApply_5 = {
                       cl <- makeCluster(getOption("cl.cores", 5))
                       r3 = parApply(cl = cl, mat, 1, FUN = function(x) {
                         library(pracma); 
                         approx_entropy(x, edim = 2, r = 0.2*sd(x), elag = 1)
                       })
                       stopCluster(cl)
                     },
                     parApply_6 = {
                       cl <- makeCluster(getOption("cl.cores", 6))
                       r3 = parApply(cl = cl, mat, 1, FUN = function(x) {
                         library(pracma); 
                         approx_entropy(x, edim = 2, r = 0.2*sd(x), elag = 1)
                       })
                       stopCluster(cl)
                     },
                     parApply_7 = {
                       cl <- makeCluster(getOption("cl.cores", 7))
                       r3 = parApply(cl = cl, mat, 1, FUN = function(x) {
                         library(pracma); 
                         approx_entropy(x, edim = 2, r = 0.2*sd(x), elag = 1)
                       })
                       stopCluster(cl)
                     }
)

## Results
mc
Unit: seconds
expr      min       lq     mean   median       uq      max neval  cld
parApply_2 2.670257 2.688115 2.699522 2.694527 2.714293 2.740149    10   c 
parApply_3 2.312629 2.366021 2.411022 2.399599 2.464568 2.535220    10 a   
parApply_4 2.358165 2.405190 2.444848 2.433657 2.485083 2.568679    10 a   
parApply_5 2.504144 2.523215 2.546810 2.536405 2.558630 2.646244    10  b  
parApply_6 2.687758 2.725502 2.761400 2.747263 2.766318 2.969402    10   c 
parApply_7 2.906236 2.912945 2.948692 2.919704 2.988599 3.053362    10    d

## Plot Results
plot(mc)

随着 矩阵变大 ，使​​用 ParApply 和 data.table 似乎比使用矩阵更快。以下示例使用了一个包含 500*500 个元素的矩阵，导致这些时间(仅用于 2 次运行):

Unit: seconds
expr      min       lq     mean   median       uq      max neval cld
ParApply 191.5861 191.5861 192.6157 192.6157 193.6453 193.6453     2   a
DT_ParAp 135.0570 135.0570 163.4055 163.4055 191.7541 191.7541     2   a

最小值要低得多，尽管最大值几乎相同，该箱线图中也很好地说明了这一点:


完整代码:

# dummy random time series data
ts <- rnorm(500)
# mat <- matrix(rep(ts,100), nrow = 100, ncol = 100)
mat = matrix(rep(ts,500), nrow = 500, ncol = 500, byrow = T)
r <- matrix(0, nrow = nrow(mat), ncol = 1)      

## Benchmark
library(microbenchmark)
mc <- microbenchmark(times=2,
                     ParApply = {
                       cl <- makeCluster(getOption("cl.cores", 3))
                       r3 = parApply(cl = cl, mat, 1, FUN = function(x) {
                         library(pracma); 
                         approx_entropy(x, edim = 2, r = 0.2*sd(x), elag = 1)
                       })
                       stopCluster(cl)
                     },
                     DT_ParAp = {
                       mDT = as.data.table(mat)
                       cl <- makeCluster(getOption("cl.cores", 3))
                       r5 = parApply(cl = cl, mDT, 1, FUN = function(x) {
                         library(pracma); 
                         approx_entropy(x, edim = 2, r = 0.2*sd(x), elag = 1)
                       })
                       stopCluster(cl)
                     }
)

## Results
mc
Unit: seconds
expr      min       lq     mean   median       uq      max neval cld
ParApply 191.5861 191.5861 192.6157 192.6157 193.6453 193.6453     2   a
DT_ParAp 135.0570 135.0570 163.4055 163.4055 191.7541 191.7541     2   a

## Plot
plot(mc)