为啥内置的 lm 函数在 R 中这么慢？

Posted 2023-03-12

【中文标题】为啥内置的 lm 函数在 R 中这么慢？

【英文标题】：Why the built-in lm function is so slow in R?为什么内置的 lm 函数在 R 中这么慢？

【发布时间】：2016-08-02 23:35:14

【问题描述】：

我一直认为 lm 函数在 R 中非常快，但正如这个例子所暗示的，使用 solve 函数计算的封闭解要快得多。

data<-data.frame(y=rnorm(1000),x1=rnorm(1000),x2=rnorm(1000))
X = cbind(1,data$x1,data$x2)

library(microbenchmark)
microbenchmark(
solve(t(X) %*% X, t(X) %*% data$y),
lm(y ~ .,data=data))

如果这个玩具示例是一个坏示例，或者lm 实际上很慢，有人可以解释一下吗？

编辑：正如 Dirk Eddelbuettel 所建议的那样，lm 需要解析公式，比较是不公平的，所以最好使用不需要解析公式的 lm.fit

microbenchmark(
solve(t(X) %*% X, t(X) %*% data$y),
lm.fit(X,data$y))


Unit: microseconds
                           expr     min      lq     mean   median       uq     max neval cld
solve(t(X) %*% X, t(X) %*% data$y)  99.083 108.754 125.1398 118.0305 131.2545 236.060   100  a 
                      lm.fit(X, y) 125.136 136.978 151.4656 143.4915 156.7155 262.114   100   b

【问题讨论】：

请记住，lm() 所做的不仅仅是解决最小二乘问题。 lm()生成的列表有12项。

请显示基准的输出。这样其他人就不需要运行它们，并且如果更高版本的 stats 包更改了这一点，它们会被保留。

您的编辑没有意义。如果你想使用lm.fit()，你需要lm.fit(X,data$y)——在这种情况下它更接近solve，尽管平均值仍然略大。

我可以，但它似乎更慢。还是要我做？

@jogo 是的，最好只更改第一个表达式

【参考方案1】：

你忽略了这一点

solve() 只返回你的参数 lm() 为您返回一个（非常丰富的）对象，其中包含 许多 组件，用于后续分析、推理、绘图…… lm() 调用的主要成本

不是

投影，而是公式 y ~ . 的分辨率，需要从中构建模型矩阵。

为了说明 Rcpp，我们编写了函数 fastLm() 的几个变体，它做了更多 lm() 所做的事情（即比基础 R 中的 lm.fit() 多一点）并对其进行了测量。参见例如this benchmark script 这清楚地表明较小数据集的主要成本是解析公式和构建模型矩阵。

简而言之，您正在通过使用基准测试来做正确的事情，但在尝试比较几乎无法比较的内容时，您做的并不完全正确：具有更大任务的子集。 p>

【讨论】：

出色的答案。实际上，使用lm.fit（我想不需要解决公式），差距就大大缩小了。您认为现在剩下的唯一区别是由于所有其他元素的计算吗？

是的，正如基准测试所证实的那样 :) 如果您喜欢这个答案，请随时接受它（通过点击“勾号”；点击“向上箭头”投票是也欢迎）。

你绝对需要两个使用线性代数到更快的版本到lm？不知道调用lm的C代码是否需要返回所有其他元素。

AFAIK 用于在 lm() 中进行线性代数的 QR 分解是 FORTRAN 代码。

【参考方案2】：

您的基准测试有问题

没有人观察到这一点真是令人惊讶！

您在solve() 中使用了t(X) %*% X。您应该使用crossprod(X)，因为X'X 是一个对称矩阵。 crossprod() 只计算矩阵的一半，同时复制其余部分。 %*% 强制计算所有。所以crossprod() 会快两倍。这解释了为什么在您的基准测试中，solve() 和 lm.fit() 之间的时间大致相同。

在我的旧 Intel Nahalem（2008 Intel Core 2 Duo）上，我有：

X <- matrix(runif(1000*1000),1000)
system.time(t(X)%*%X)
#    user  system elapsed 
#   2.320   0.000   2.079 
system.time(crossprod(X))
#    user  system elapsed 
#    1.22    0.00    0.94 

如果您的计算机速度更快，请尝试改用X <- matrix(runif(2000*2000),2000)。

下面，我将解释所有拟合方法所涉及的计算细节。

---

QR 分解与Cholesky 分解

lm() / lm.fit() 是基于 QR 的，而 solve() 是基于 Cholesky 的。 QR分解的计算成本为2 * n * p^2，而Cholesky方法为n * p^2 + p^3（n * p^2用于计算矩阵叉积，p^3用于Cholesky分解）。所以你可以看到，当n 远大于p 时，Cholesky 方法比 QR 方法快大约 2 倍。所以这里真的没有必要进行基准测试。 （如果你不知道，n 是数据个数，p 是参数个数。）

---

LINPACK QR 对比LAPACK 二维码

通常，lm.fit() 使用（修改后的）LINPACK QR 分解算法，而不是 LAPACK QR 分解算法。可能你对BLAS/LINPACK/LAPACK不是很熟悉；这些是提供内核科学矩阵计算的 FORTRAN 代码。 LINPACK 调用 1 级 BLAS，而LAPACK 使用块算法调用 3 级 BLAS。平均而言，LAPACK QR 比LINPACK QR 快 1.6 倍。 lm.fit() 不使用LAPACK 版本的关键原因是需要部分列旋转。 LAPACK 版本进行全列旋转，使summary.lm() 更难使用R QR 分解矩阵因子产生 F 统计量和ANOVA 检验。

---

旋转 vs.没有旋转

来自RcppEigen 包的fastLm() 使用LAPACK 非枢轴QR 分解。同样，您可能不清楚 QR 分解算法和旋转问题。你只需要知道带旋转的 QR 分解在 level-3 BLAS 中只有 50% 的份额，而没有旋转的 QR 分解在 level-3 BLAS 中的份额为 100%。在这方面，放弃旋转将加快 QR 分解过程。当然，最终结果是不同的，当模型矩阵秩不足时，不进行旋转会产生危险的结果。

有一个很好的问题与您从fastLM 获得的不同结果有关：Why does fastLm() return results when I run a regression with one observation?。 @BenBolker、@DirkEddelbuettel 和我在 cmets 中对 Ben 的回答进行了非常简短的讨论。

---

结论：你想要速度还是数值稳定性？

在数值稳定性方面，有：

LINPACK pivoted QR > LAPACK pivoted QR > pivoted Cholesky > LAPACK non-pivoted QR

在速度方面，有：

LINPACK pivoted QR < LAPACK pivoted QR < pivoted Cholesky < LAPACK non-pivoted QR

正如德克所说，

FWIW RcppEigen 包在其fastLm() 示例中具有更完整的分解集。但在这里，正如 Ben 雄辩地指​​出的那样：“这是你为速度付出的一部分代价。”。我们给你足够的绳子让你上吊。如果您想保护自己免受自己的伤害，请坚持使用lm() 或lm.fit()，或者创建一个混合“快速但安全”的版本。

---

快速稳定的版本

检查my answer Here。