是否可以缓存 lm() 矩阵以适应新数据？

Question

我编写了一个算法，它用 lm() 拟合线性模型，然后迭代地 "updates" 响应变量。问题是：在高维场景中，拟合线性模型会产生瓶颈。

另一方面，所需的大部分工作是仅依赖于协变量矩阵 X 的矩阵求逆，即系数由下式给出：solve(t(X) %*% X) %*% X %*% y.

阅读lm()代码，了解到R使用QR分解。

是否可以恢复使用的内部矩阵运算并更快地用新的 y 值拟合新模型？

这是一个最小的例子：

set.seed(1)
X <- matrix(runif(400*150000), nrow = 150000)
y1 <- runif(150000)
y2 <- runif(150000)

mod1 <- lm(y1 ~ X)
mod2 <- lm(y2 ~ X)

理论上，mod2 "repeats" 昂贵的矩阵运算与第一个 lm() 调用中的相同。

我想继续使用 lm() 因为它的高效实施和自动处理不完整秩矩阵的能力。

Answer 1

# Data
set.seed(1)
n = 5
X <- matrix(runif(5*n), nrow = n)
y1 <- runif(n)
y2 <- runif(n)

# lm models
mod1 <- lm(y1 ~ X)
mod2 <- lm(y2 ~ X)

# Obtain QR decomposition of X
q = qr(X)

# Reuse 'q' to obtain fitted values repeatedly
mod1_fv = qr.fitted(q, y1)
mod2_fv = qr.fitted(q, y2)

# Compare fitted values from reusing 'q' to fitted values in 'lm' models
Vectorize(all.equal)(unname(mod1$fitted.values), mod1_fv)
#> [1] TRUE TRUE TRUE TRUE TRUE

Vectorize(all.equal)(unname(mod2$fitted.values), mod2_fv)
#> [1] TRUE TRUE TRUE TRUE TRUE

^{由 reprex package (v0.3.0)}

创建于 2019-09-06

Answer 2

您试过只拟合多元模型吗？我没有检查代码，但在我的系统上，它的速度几乎是单独安装的一半，所以如果它按照你在幕后的建议进行操作，我不会感到惊讶。也就是说，

mods <- lm(cbind(y1, y2) ~ X)

是否可以缓存 lm() 矩阵以适应新数据？

Is it possible to cache `lm()` matrices to fit new data?

r

lm