是否有替代 for-loop 的方法来拟合以矩阵作为因变量的样条曲线？

Question

我有一个大约有 30 000 行和 15 列的矩阵，以及一个大小为 1x15 的向量。我想拟合一个样条曲线，矩阵中的每一行作为因变量，向量作为预测变量。对于每个样条，我想用单个 x 值进行预测，并将所有预测添加到一个向量中。

有什么方法可以跳过 for 循环来解决这个问题并降低时间复杂度吗？

这里是一些示例数据：

矩阵：

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
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
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

向量：

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

我目前的解决方案是遍历矩阵的行，将每一行和向量附加到一个临时数据框，拟合样条并进行预测。我也尝试过 sapply，时间复杂度的改进非常有限。

for循环解决方案：

library(splines)
beta = function(matrix, vector){
  predictions = c()
  for (i in 1:nrow(matrix)){
    # Temporary data frame
    temp_df = data.frame(y = matrix[i,], x = vector)
    
    # Fit a spline for each observation
    
    spline = lm(y ~ ns(x, df = 7), data = temp_df)
    # Predict a value and add to vector
    predictions = c(predictions, predict(spline, data.frame(x = 10000)))
  }
  return(predictions)
}
system.time(beta(matrix, vector)) # user 62.89

应用解决方案：

fun = function(i){
  predictions = c()
  temp_df = data.frame(y = matrix[i, ], x = vector) 
  predictions = c(predictions, predict(lm(y ~ ns(x, df = 7), data = temp_df), data.frame(x = 10000)))
  return(predictions)
}
system.time(sapply(1:nrow(matrix), fun)) # user: 53.42

样条的类型对我的解决方案并不重要，使用的包也不重要。我试图同时直接为矩阵的所有行拟合样条，但没有成功。

我需要能够将矩阵扩展到大约 1 500 000 x 15，并在短时间内做出几个不同的预测。真的很感激这方面的帮助。

提前致谢！

Answer 1

您可以将您的 30000x15 矩阵转置为 15x30000 并对其进行多元线性模型，然后预测将全部发生在一行中。

mat <- read.table(text="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
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
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")

vec <- 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)
mat <- t(mat)
colnames(mat) <- paste0("y", 1:ncol(mat))
dat <- cbind(as.data.frame(mat), x=vec)

form <- paste0("cbind(", paste(colnames(mat), collapse=", "), ") ~ ns(x, df=7)")

library(splines)
mod <- lm(form, data=dat)
coef(mod)
#>                        y1         y2         y3
#> (Intercept)     1.0320976  1.0303319  1.0356140
#> ns(x, df = 7)1  0.2965524  0.2774531 -0.1627959
#> ns(x, df = 7)2 -0.6004987 -0.5797548 -0.5066507
#> ns(x, df = 7)3 -1.2040751 -1.1148353 -0.6408537
#> ns(x, df = 7)4 -0.9342747 -0.9346508 -0.6687842
#> ns(x, df = 7)5 -1.0382844 -1.0351010 -0.8216596
#> ns(x, df = 7)6 -1.2107326 -1.1990489 -1.1924609
#> ns(x, df = 7)7 -0.9446160 -0.9469103 -0.8121099

predict(mod, newdata=data.frame(x=10000))
#>          y1        y2        y3
#> 1 0.9415758 0.9442323 0.8442096

^{创建于 2023-04-01 与 reprex v2.0.2}

在我的机器（MacBook Pro、M1 Max）上，for 循环大约需要 27 秒，而在 30k x 15 值矩阵上的多元线性模型大约需要 15 秒。

另一种选择是只进行矩阵乘法来生成系数，然后根据系数生成预测。

matfun <- function(matrix, vector){
  require(Matrix)
  n <- ns(vector, df=7)
  X <- model.matrix(~ns(vector, df=7)) 
  z <- solve(t(X) %*% X) %*% t(X)
  b <- crossprod(t(z), t(matrix))
  predX <- c(1, ns(10000, df=7, knots=attr(n, "knots"), Boundary.knots = attr(n, "Boundary.knots")))
  preds <- t(b) %*% predX
  return(preds)
}

对我来说，

matfun(mat, vet)

对于 30k x 15 矩阵需要 0.007 秒 - 比其他任何一种方法都快。所有三个解决方案都产生相互关联为 1 的预测。如果你想要的只是预测，矩阵函数是最快的。