为什么内置的lm功能在R中如此之慢？

Question

我一直认为lmR 中的函数非常快,但正如本例所示,使用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

Answer 1

你忽略了那个

• solve() 只返回你的参数
• lm()返回一个(非常丰富的)具有许多组件的对象,用于后续分析,推理,绘图,......
• 您lm()调用的主要成本不是投影,而是y ~ .需要构建模型矩阵的公式的分辨率.

为了说明Rcpp,我们写了一些函数的变体,它们fastLm()做了更多的事情lm()(即比lm.fit()基数R 更多)并测量它.参见例如此基准脚本,其清楚地表明较小数据集的主要成本在于解析公式和构建模型矩阵.

简而言之,您通过使用基准测试来做正确的事情,但是在尝试比较大多数无法比较的事情时,您所做的并不是那么正确:具有更大任务的子集.

Answer 2

您的基准测试有问题

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

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

在我的旧英特尔Nahalem(2008英特尔酷睿2双核处理器)上,我有:

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 vs LAPACK QR

通常,lm.fit()使用(修改)LINPACKQR分解算法,而不是LAPACKQR分解算法.也许你不是很熟悉BLAS/LINPACK/LAPACK; 这些是FORTRAN代码,提供内核科学矩阵计算.LINPACK调用level-1 BLAS,同时使用块算法LAPACK调用level-3 BLAS.平均而言,LAPACKQR比LINPACKQR 快1.6倍.不使用版本的关键原因是需要部分列旋转.lm.fit()LAPACK LAPACK版本进行全列旋转,使得summary.lm()使用RQR分解的矩阵因子来生成F统计量和ANOVA测试变得更加困难.

枢轴与无枢轴

fastLm()从RcppEigen包使用LAPACK非旋转QR分解.同样,您可能不清楚QR分解算法和旋转问题.您只需要知道具有旋转的QR分解仅具有级别3的50%份额BLAS,而没有旋转的QR分解具有级别3的100%份额BLAS.在这方面,放弃旋转将加速QR分解过程.当然,最终的结果是不同的,当模型矩阵排名不足时,没有旋转会产生危险的结果.

有一个很好的问题与你得到的不同结果有关fastLM:为什么fastLm()当我用一次观察运行回归时返回结果？.@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 RcppEigen包在其fastLm()示例中具有更全面的分解集.但就像Ben如此雄辩地说:"这是你为速度付出代价的一部分." 我们给你足够的绳索来吊死自己.如果你想保护自己免受自己,坚持lm()或lm.fit(),或创建一个混合"快但安全"的版本.

快速稳定的版本

检查我的答案在这里.