如何避免循环实现矩阵行作为因变量的样条拟合与预测

问题概述 #

需要处理一个最多150万行15列的矩阵，以每一行作为因变量，搭配同一个1×15的预测变量向量拟合样条曲线，随后对x=10000进行预测，最终汇总所有预测结果为向量。当前循环和sapply实现耗时过长，需用非循环的矩阵运算方案大幅提升效率。

示例数据 #

矩阵（3行15列） #

matrix_example <- rbind(
  c(1, 0.9999866, 0.9999833, 0.9999822, 0.9998178, 0.9996189, 0.9994455, 0.007492490, 0.007492490, 0.007492195, 0.007464383, 0.0003291809, 0.0003291808, 0.00002728396, 0.000017999925),
  c(1, 0.9997588, 0.9990516, 0.9990033, 0.9959569, 0.9942259, 0.9920646, 0.063989436, 0.063989428, 0.063980612, 0.063502466, 0.0052701181, 0.0052700809, 0.00079669065, 0.000497011826),
  c(1, 0.9882412, 0.7925734, 0.7920651, 0.7890917, 0.7312206, 0.7283561, 0.424428825, 0.423345436, 0.422478875, 0.409804031, 0.2936134533, 0.2902640241, 0.13727615950, 0.085531730428)
)

预测变量向量 #

vector_example <- c(0.000, 5689.072, 11915.687, 19188.547, 27796.767, 37742.035, 49564.349, 64430.295, 84381.754, 111870.835, 149611.382, 221043.651, 362982.876, 583956.304, 1120546.126)

当前低效实现 #

For循环版本 #

library(splines)
beta = function(matrix, vector){
  predictions = c()
  for (i in 1:nrow(matrix)){
    temp_df = data.frame(y = matrix[i,], x = vector)
    spline = lm(y ~ ns(x, df = 7), data = temp_df)
    predictions = c(predictions, predict(spline, data.frame(x = 10000)))
  }
  return(predictions)
}
system.time(beta(matrix_example, vector_example)) # user 62.89（针对3万行数据）

sapply版本 #

fun = function(i){
  temp_df = data.frame(y = matrix_example[i, ], x = vector_example) 
  predict(lm(y ~ ns(x, df = 7), data = temp_df), data.frame(x = 10000))
}
system.time(sapply(1:nrow(matrix_example), fun)) # user: 53.42（针对3万行数据）

高效优化方案 #

核心思路：所有模型共享同一个预测变量向量，因此样条基函数的设计矩阵固定不变。利用矩阵运算一次性完成所有行的拟合与预测，避免循环中的重复计算开销。

优化代码实现 #

library(splines)

batch_spline_predict <- function(Y, x, x0, df=7) {
  # 1. 构建通用样条设计矩阵（含截距项）
  ns_x <- ns(x, df=df)
  X <- model.matrix(~ ns_x)
  
  # 2. 计算设计矩阵的广义逆，仅需计算一次
  X_crossprod <- crossprod(X)
  X_chol <- chol(X_crossprod)
  X_inv <- chol2inv(X_chol) %*% t(X)
  
  # 3. 构建预测点的样条设计矩阵，需与训练时的节点一致
  ns_x0 <- ns(x0, df=df, knots=attr(ns_x, "knots"), Boundary.knots=attr(ns_x, "Boundary.knots"))
  X0 <- model.matrix(~ ns_x0)
  
  # 4. 批量计算所有行的预测值，利用矩阵乘法完成
  predictions <- Y %*% X_inv %*% t(X0)
  
  # 转换为向量输出
  as.vector(predictions)
}

测试与性能对比 #

# 测试示例数据
system.time(preds <- batch_spline_predict(matrix_example, vector_example, x0=10000))
print(preds)

# 针对150万行数据测试（模拟数据）
set.seed(123)
large_matrix <- matrix(runif(1500000*15), nrow=1500000)
system.time(large_preds <- batch_spline_predict(large_matrix, vector_example, x0=10000))

效率提升原因 #

• 避免重复计算：样条基函数、设计矩阵及其广义逆仅计算一次，无需在循环中重复生成
• 矩阵运算优化：R的矩阵运算基于BLAS/LAPACK底层库，处理百万级数据时远快于逐行循环
• 减少内存开销：无需创建大量临时数据框，直接用原始矩阵运算，内存占用更低

内容的提问来源于stack exchange，提问作者LuBe