滑动窗口线性回归的最快实现？R语言代码优化与性能对比

问题与原始实现 #

计算滑动窗口线性回归的最高效方式是什么？以下是需要提速的示例代码：

df <- as.data.frame(cbind(1:100,cumsum(rnorm(100))))
colnames(df) <- c("time","val")

n <- 10
roll_lm <- rep(NA,nrow(df))
for(i in n:nrow(df)){
  roll_df <- df[(i-(n-1)):i,]
  mod <- lm(val ~ time , roll_df)
  roll_lm[i] <- tail(predict(mod, roll_df),1)
} 

plot(df)
lines(roll_lm, col=4,lwd=2)

优化实现方案 #

以下是几种不同的优化实现函数：

fu1（基础循环调用lm） #

fu1 <- function(df,n=10){
  roll_lm <- rep(NA,nrow(df))
  for(i in n:nrow(df)){
    roll_df <- df[(i-(n-1)):i,]
    mod <- lm(val ~ time , roll_df)
    roll_lm[i] <- tail(predict(mod, roll_df),1)
  } 
  return(roll_lm)
}

fu2（apply+lm） #

fu2 <- function(df,n=10){
  roll2 <- apply(embed(df$val, n),  # 创建滑动窗口矩阵
               1,                  # 按行应用函数
               function(x) head(predict(lm(x ~ seq_along(x))), 1))
  roll2 <- c(rep(NA, 9), roll2)
  return(roll2)
}

fu3（手动计算回归系数） #

fu3 <- function(df,n=10){
  roll_lm <- rep(NA_real_, nrow(df))
  for (i in n:nrow(df)) {
    x  <- df$time[(i-(n-1)):i]
    y  <- df$val[(i-(n-1)):i]
    b1 <- cov(x, y) / var(x, x) 
    b0 <- mean(y) - b1*mean(x)
    roll_lm[i] <- b0 + df$time[i] * b1
  } 
  return(roll_lm)
}

fu4（RcppRoll滚动统计） #

fu4 <- function(df,n=10){
  library(RcppRoll)
  roll_lm <- rep(NA,nrow(df))
  x <- df$time
  y <- df$val
  for(i in n:nrow(df)){
    roll_x <- x[(i-(n-1)):i]
    roll_y <- y[(i-(n-1)):i]
    roll_mean_x <- roll_mean(roll_x, n)
    roll_mean_y <- roll_mean(roll_y, n)
    roll_cov_xy <- roll_sum((roll_x - roll_mean_x) * (roll_y - roll_mean_y), n) / (n - 1)
    roll_var_x <- roll_var(roll_x, n)
    beta_hat <- roll_cov_xy / roll_var_x
    alpha_hat <- roll_mean_y - beta_hat * roll_mean_x
    roll_lm[i] <- tail(alpha_hat + beta_hat * roll_x, 1)
  }
  return(roll_lm)
}

fu5（rollRegres包实现） #

fu5 <- function(df,n=10){
  library(rollRegres)
  fit <- roll_regres(val ~ time, data = df, width = n)
  return(fit$coefs[,1] + fit$coefs[,2] * df$time)
}

性能对比 #

使用microbenchmark对各方案进行性能测试，结果如下：

microbenchmark::microbenchmark(fu1(df),fu2(df),fu3(df),fu4(df),fu5(df))

测试输出：

Unit: milliseconds
    expr        min        lq      mean     median         uq       max neval
 fu1(df) 315.878570 338.35775 379.62677 359.243318 390.317795 831.82362   100
 fu2(df) 233.707004 245.90377 273.49740 256.383181 288.005876 697.46965   100
 fu3(df)  13.957729  15.37773  18.01751  16.442809  18.554429  67.16131   100
 fu4(df)  25.305796  27.50037  32.42579  29.063492  32.757900 100.82235   100
 fu5(df)   1.782347   2.02923   2.79950   2.179755   3.030447  16.40774   100

结论 #

从测试结果可以清晰看到：

• fu5（基于rollRegres包）的性能最优，平均耗时仅约2.8毫秒，是所有方案中最快的；
• fu3通过手动计算回归系数，跳过lm函数的冗余流程，性能远优于直接调用lm的fu1和fu2；
• fu1和fu2由于每次循环都调用完整的lm模型，包含大量不必要的计算步骤，性能最差，耗时是fu5的上百倍；
• fu4虽然使用了RcppRoll的高效滚动统计函数，但仍在循环内逐窗口计算，效率不如专门的滚动回归包rollRegres。

因此，推荐使用rollRegres包实现滑动窗口线性回归，这是当前最高效的方案。

内容的提问来源于stack exchange，提问作者mr.T