通过RasterStack对象的单个单元格进行速度优化循环

Question

我正在使用两个RasterStack对象,每个对象由十个层组成,代表单个时间步.

# Mock data
pred.rst.stck <- do.call("stack", lapply(seq(10), function(i) {
  pred.rst <- raster(nrows = 15, ncols = 15, xmn= 0, xmx = 10, ymn = 0, ymx = 10)
  pred.rst[] <- rnorm(225, 50, 10)
  return(pred.rst)
})
resp.rst.stck <- do.call("stack", lapply(seq(10), function(i) {
  resp.rst <- raster(nrows = 10, ncols = 10, xmn = 0, xmx = 10, ymn = 0, ymx = 10)
  resp.rst[] <- rnorm(100, 50, 10)
  return(resp.rst)
})

pred.rst.stck用作预测变量resp.rst.stck集和响应变量集.对于预测器RasterStack的每个单个单元格,我想在响应RasterStack的每个单元格上拟合线性模型,提取每个拟合模型的相应R平方并将它们相加.简而言之,这是迄今为止使用R parallel包的最快方法:

# Parallelization
library(parallel)
n.cores <- detectCores()
clstr <- makePSOCKcluster(n.cores)

# Extract cell values from RasterStack objects
pred.vals <- getValues(pred)
resp.vals <- getValues(resp)
clusterExport(clstr, c("pred.vals", "resp.vals"))

# Loop through all predictor cells
rsq.sums <- parLapply(clstr, seq(nrow(pred.vals)), function(i) {
  # For each predictor cell, loop through all response cells, 
  # fit linear model, extract and sum up single R-squared 
  do.call("sum", lapply(seq(nrow(resp.vals)), function(j) {
    summary(lm(resp.vals[j, ] ~ pred.vals[i, ]))$r.squared
  }))
})

虽然parLapply与普通相比表现要好得多lapply,但我想知道是否有一种优雅的方法来加速整个过程.有什么建议？

干杯,
弗洛里安

Answer 1

你可以尝试一些技巧.我不太了解你创建线性模型的方式,但是r.squared你从线性模型中计算得相当于Pearson相关系数的平方(cor在R中有默认参数),这比线性计算得快得多模型.

使用您的数据比较这两个函数:

# Finding r.squared using a lm
f1 <- function(n){summary(lm(resp.vals[n,] ~pred.vals[n,]))$r.squared}
# Finding r.squared using Pearson's
f2 <- function(n){ cor(resp.vals[n,],pred.vals[n,])^2}

# Both give the same result
f1(1)
[1] 0.0009032986
f2(1)
[1] 0.0009032986

就单个调用这些函数的时间而言:

require(microbenchmark)
microbenchmark( f1(1) , f2(1) )
Unit: microseconds
  expr      min        lq   median        uq      max neval
 f1(1) 2009.328 2037.0680 2071.045 2124.9225 4102.079   100
 f2(2)   73.075   77.7365   84.690   89.7215 5485.368   100

所以你应该能够通过使用cor代替来减少25倍的处理时间lm.

一个快速的系统时间比较交换原始功能使用cov()^2显示了这种情况:

#Using cov()
   user  system elapsed 
  0.013   0.002   1.085

#Using lm()
   user  system elapsed 
  0.159   0.028  26.334