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

Question

我编写了一个算法，该算法使用 lm() 拟合线性模型，然后迭代地“更新”响应变量。问题是:在高维场景中，拟合线性模型会产生瓶颈。

另一方面，所需的大部分工作是仅依赖于协变量矩阵 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 会“重复”与第一次 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 创建于 2019-09-06 (v0.3.0)}