为什么内置的 lm 函数在 R 中这么慢？
Why the built-in lm function is so slow in R?
我一直认为 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
你忽略了那个
solve()只有returns你的参数lm() returns 你是一个（非常丰富的）对象，具有 许多 组件，用于后续分析、推理、绘图、...- 您的
lm() 调用的主要成本是 而不是 投影，而是需要构建模型矩阵的公式 y ~ . 的分辨率.
为了说明 Rcpp，我们编写了函数 fastLm() 的一些变体，做更多 lm() 所做的事情（即比基础 R 的 lm.fit() 多一点）并对其进行了测量。参见例如this benchmark script 这清楚地表明 较小数据集的主要成本在于解析公式和构建模型矩阵 。
简而言之，您通过使用基准测试做了正确的事，但在尝试比较几乎无法比较的东西时，您做的并不完全正确：一个子集有很多更大的任务。
您的基准测试有问题
没有人观察到这一点，真是令人惊讶！
您在 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 因式分解 v.s。 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倍左右。所以这里真的没有必要做benchmark。 （如果不知道，n是数据个数，p是参数个数）
LINPACK QR v.s。 LAPACK QR
通常，lm.fit() 使用（修改）LINPACK QR 分解算法，而不是 LAPACK QR 分解算法。也许你对BLAS/LINPACK/LAPACK不是很熟悉；这些是提供内核科学矩阵计算的 FORTRAN 代码。 LINPACK 调用 level-1 BLAS，而 LAPACK 使用块算法调用 level-3 BLAS。平均而言，LAPACK QR 比 LINPACK QR 快 1.6 倍。 lm.fit() 不使用 LAPACK 版本的关键原因是需要部分列旋转。 LAPACK 版本进行全列旋转，使其成为summary.lm() 更难使用 QR 分解的 R 矩阵因子来生成 F 统计量和 ANOVA 检验。
旋转 v.s。没有旋转
fastLm() 来自 RcppEigen 包使用 LAPACK 非透视 QR 分解。同样，您可能不清楚 QR 分解算法和旋转问题。你只需要知道带旋转的 QR 分解只有 3 级的 50% 份额 BLAS，而没有旋转的 QR 分解有 3 级的 100% 份额 BLAS。在这方面，放弃旋转将加快 QR 分解过程。当然，最终结果是不同的，当模型矩阵秩亏时，没有主元会给出危险的结果。
有一个很好的问题与您从 fastLM 得到的不同结果有关：Why does fastLm() return results when I run a regression with one observation?。 @BenBolker、@DirkEddelbuettel 和我在对 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 the RcppEigen package has a fuller set of decompositions in its fastLm() example. But here it is as Ben so eloquently stated: "this is part of the price you pay for speed.". We give you enough rope to hang yourself. If you want to protect yourself from yourself, stick with lm() or lm.fit(), or create a hybrid 'fast-but-safe' version.
快速稳定版
检查my answer Here。
我一直认为 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
你忽略了那个
solve()只有returns你的参数lm()returns 你是一个（非常丰富的）对象，具有 许多 组件，用于后续分析、推理、绘图、...- 您的
lm()调用的主要成本是 而不是 投影，而是需要构建模型矩阵的公式y ~ .的分辨率.
为了说明 Rcpp，我们编写了函数 fastLm() 的一些变体，做更多 lm() 所做的事情（即比基础 R 的 lm.fit() 多一点）并对其进行了测量。参见例如this benchmark script 这清楚地表明 较小数据集的主要成本在于解析公式和构建模型矩阵 。
简而言之，您通过使用基准测试做了正确的事，但在尝试比较几乎无法比较的东西时，您做的并不完全正确：一个子集有很多更大的任务。
您的基准测试有问题
没有人观察到这一点，真是令人惊讶！
您在 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 因式分解 v.s。 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倍左右。所以这里真的没有必要做benchmark。 （如果不知道，n是数据个数，p是参数个数）
LINPACK QR v.s。 LAPACK QR
通常，lm.fit() 使用（修改）LINPACK QR 分解算法，而不是 LAPACK QR 分解算法。也许你对BLAS/LINPACK/LAPACK不是很熟悉；这些是提供内核科学矩阵计算的 FORTRAN 代码。 LINPACK 调用 level-1 BLAS，而 LAPACK 使用块算法调用 level-3 BLAS。平均而言，LAPACK QR 比 LINPACK QR 快 1.6 倍。 lm.fit() 不使用 LAPACK 版本的关键原因是需要部分列旋转。 LAPACK 版本进行全列旋转，使其成为summary.lm() 更难使用 QR 分解的 R 矩阵因子来生成 F 统计量和 ANOVA 检验。
旋转 v.s。没有旋转
fastLm() 来自 RcppEigen 包使用 LAPACK 非透视 QR 分解。同样，您可能不清楚 QR 分解算法和旋转问题。你只需要知道带旋转的 QR 分解只有 3 级的 50% 份额 BLAS，而没有旋转的 QR 分解有 3 级的 100% 份额 BLAS。在这方面，放弃旋转将加快 QR 分解过程。当然，最终结果是不同的，当模型矩阵秩亏时，没有主元会给出危险的结果。
有一个很好的问题与您从 fastLM 得到的不同结果有关：Why does fastLm() return results when I run a regression with one observation?。 @BenBolker、@DirkEddelbuettel 和我在对 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 the RcppEigen package has a fuller set of decompositions in its
fastLm()example. But here it is as Ben so eloquently stated: "this is part of the price you pay for speed.". We give you enough rope to hang yourself. If you want to protect yourself from yourself, stick withlm()orlm.fit(), or create a hybrid 'fast-but-safe' version.
快速稳定版
检查my answer Here。