r - 通过 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 相当于 PIL 逊相关系数的平方 (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 而不是 lm，您应该能够将处理时间缩短 25 倍。

将原始函数替换为使用 cov()^2 的快速系统时间比较表明情况确实如此:

#Using cov()
   user  system elapsed 
  0.013   0.002   1.085

#Using lm()
   user  system elapsed 
  0.159   0.028  26.334