r-是否有替代for循环拟合样条与矩阵作为因变量?

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 with reprex v2.0.2}

在我的机器(MacBook Pro, M1 Max)上，在30k x 15的值矩阵上，for循环大约需要27秒，多元线性模型大约需要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)
}

对我来说，30k x 15矩阵的matfun(mat, vet)花费了0.007秒-比其他任何一种替代方案都快得多。所有三种解决方案产生的预测在1时相互关联。如果你想要的只是预测，矩阵函数是最快的。