嫌lm( )慢？试试lm.fit( )

youngtf

（如有错误，烦请各位指点，实话说对这结论我有些没把握）

前一段时间要对一组数据做分析，其中的一步是对样本量约900的一组样本的数据做一元回归，变量的数目比较多，估计是8亿多不到9亿。简单地试验估计下，在R中用lm()函数来获得一元回归的结果的话，这个过程大概要20-30天，虽说比较大的数据处理起来就是要比较久，不过想想一次性跑一个月，也实在是有些沮丧。

加快计算的方法也考虑过不少，一方面是并行计算，不过我缺乏相关的经验，对于并行计算所需要的硬件和软件完全是没有了解。虽然学习很快乐，不过毕竟不能立竿见影。另一方面是修改代码，让计算时间尽可能缩短一点。在后面这个过程中，我找到了一个似乎可行的办法，就是使用lm.fit( )来代替lm( )函数。

在R的基础文档中，对lm.fit( )的介绍并不多，可能是因为如文档中所说，一开始这个函数就不建议用户直接使用的原因。然而在glm( )和glm.fit( )的介绍中有提到，在设计矩阵已经确定的情况下，glm.fit( )可能会是比glm( )更有效率的选择。

文档里说的很清楚（我没理解错的话），lm.fit( )是lm( )中计算的核心部分，除非是熟悉的用户否则不要使用。我战战兢兢地试用了下，结果如下：

<br />
	> a = rnorm(1000)<br />
	> b = matrix(rnorm(5000000),1000,5000)</p>
<p>	> func1 = function(y,x){<br />
	    beta = vector("numeric", NCOL(x))<br />
	    for (i in 1:NCOL(x)){<br />
          beta[i] = lm(y~x[,i])$coefficients[2]<br />
	    }<br />
	    beta<br />
	  }<br />
	> func2 = function(y,x){<br />
	    beta = vector("numeric", NCOL(x))<br />
	    for (i in 1:NCOL(x)){<br />
          beta[i] = lm.fit(y = y, x = cbind(rep(1,NROW(x)),x[,i]))$coefficients[2]<br />
	    }<br />
	    beta<br />
	  }</p>
<p>	> system.time(func1(a,b))<br />
	   user  system elapsed<br />
	  11.810  0.368  12.198<br />
	> system.time(func2(a,b))<br />
	   user  system elapsed<br />
	  1.244   0.200   1.446<br />
	> sum(which(func1(a,b) != func2(a,b)))<br />
	  [1] 0<br />

如上文所说，使用lm.fit( )时要自己构造设计矩阵，所以在自变量前面要加上截距对应的一列1。最后一句的结果可以看出，这两个函数计算出来的回归系数是一样的。而计算的速度的话,差距应该算是非常明显了，至少我当时是吓了一跳。

nan.xiao

和此文第二节说的情况差不多 http://cos.name/2009/12/improve-r-computation-efficiency/

R 中也有简单好用的并行方案，比如 foreach, 可以拿来试试。

推荐下 RcppArmadillo::fastLmPure, 速度还可以再提高一倍。

<br />
require(RcppArmadillo)<br />
set.seed(1001)<br />
y = rnorm(1000)<br />
x = matrix(rnorm(10000000), 1000, 10000)<br />
system.time(for (i in 1:ncol(x)) lm(y ~ x[ , i])$coefficients[2])<br />
# 用户   系统   流逝<br />
# 37.891  0.033 38.032<br />
system.time(for (i in 1:ncol(x)) lm.fit(y = y, x = cbind(rep(1,NROW(x)),x[,i]))$coefficients[2])<br />
# 用户  系统  流逝<br />
# 4.206 0.113 4.333<br />
system.time(for (i in 1:ncol(x)) fastLmPure(y = y, X = cbind(rep(1, nrow(x)), x[ , i]))$coefficients[2])<br />
# 用户  系统  流逝<br />
# 2.180 0.000 2.186<br />

youngtf

回复第2楼的肖楠：逛COS不仔细，漏了这么好的文章。非常感谢！

读了才想起来一件事，顺带一提算是补充，apply是我最早试用想加快速度的，结果几乎没效果。

nan.xiao

*apply 在这里不太适用，即使提高也有限。

测试了下 foreach + doMC 并行, 效果不佳。可能是我的核不够多，列数不足。

再往下就是底层矩阵运算子程序的速度了。这里可以看到一个评测结果。把默认的 BLAS 换成 MKL 估计可以节约一定的时间。RcppArmadillo 和原生 R 都可以利用 MKL.

youngtf

回复第4楼的肖楠：QR分解是不是可以用qr( )函数？想要更快还是要乖乖学其他效率更高的语言是吧。。。

nan.xiao

回复第5楼的 youngtf：可以。只是不知道速度怎么样。

我觉得差不多了。现在这样单机两天也能算完了。想要更快就并行吧。

youngtf

回复第2楼的肖楠：才看到2楼的推荐，多谢！长了不少见识。

micro@

1.舍不了内存套不着速度；

2.拒绝一切免费赠品。

<br />
set.seed(1001)<br />
y = rnorm(1000)<br />
x = matrix(rnorm(10000000), 1000, 10000)</p>
<p>system.time(for (i in 1:ncol(x)) lm.fit(y = y, x = cbind(rep(1,NROW(x)),x[,i]))$coefficients[2])<br />
#   user  system elapsed<br />
#  3.661   0.000   3.662</p>
<p>system.time(cor(y,x)*sd(y)/apply(x,2L,sd))<br />
#   user  system elapsed<br />
#  1.173   0.011   1.184</p>
<p>system.time({mx=colSums(x); (y%*%x-mx*mean(y))/(colSums(x*x)-mx*mx/length(y))})<br />
#   user  system elapsed<br />
#  0.173   0.012   0.186</p>
<p>

micro@

回复第6楼的肖楠：回复第5楼的 youngtf：QR分解应该只在分解结果可以重用的情况下比较有优势，比如例中x,y位置互换的时候，只需要一次QR。

youngtf

回复第8楼的 micro@：果然把帖子发到这来是对的。。太有趣了，多谢多谢。。

youngtf

回复第8楼的 micro@：这三种方法计算出来的结果略有不同是吧。。

micro@

回复第11楼的 youngtf：吓唬我。。。

<br />
> set.seed(1001)<br />
> y = rnorm(1000)<br />
> x = matrix(rnorm(10000000), 1000, 10000)</p>
<p>> b0=numeric(ncol(x)); for (i in 1:ncol(x)) lm.fit(y = y, x = cbind(rep(1,NROW(x)),x[,i]))$coefficients[2]->b0[i]</p>
<p>> cor(y,x)*sd(y)/apply(x,2L,sd)->b1</p>
<p>> local({mx=colSums(x); (y%*%x-mx*mean(y))/(colSums(x*x)-mx*mx/length(y))})->b2</p>
<p>> range(b0-b1)<br />
[1] -3.191891e-16  3.885781e-16<br />
> range(b0-b2)<br />
[1] -3.469447e-16  4.163336e-16</p>
<p>

If these are the differences you meant, search similar posts in this forum.

youngtf

回复第12楼的 micro@：多谢，我说的就是这个，找到其他帖子里你们的讨论了。

bjt

回复第8楼的 micro@：估摸着你的机器上有编译好的并行BLAS（MKL）加速，这矩阵运算是必然并行的

frank230

我很想知道什么数据变量维度有8亿之多！

micro@

回复第14楼的刘思喆：

Why do you think so? Now I'm only using element-wise matrix operations:

<br />
> system.time({mx=colSums(x); (colSums(y*x)-mx*mean(y))/(colSums(x*x)-mx*mx/length(y))})<br />
   user  system elapsed<br />
  0.259   0.060   0.319<br />

In this version, I cannot see where BLAS is used. Actually lm.fit uses QR decomposition, which usually depends on LAPACK/BLAS.

<br />
$ R CMD config BLAS_LIBS<br />
-L/home/***/lib64/R/lib -lRblas<br />
$ R CMD config LAPACK_LIBS<br />
-L/home/***/lib64/R/lib -lRlapack<br />

These seem pointing to the default LAPACK/BLAS of R on this system.

In the last version of the code above, it seems to me that R_binary and real_binary haven't even used openmp. But colSums does use openmp. So, let's check the benefit of openmp:

<br />
(default.max=.Internal(setMaxNumMathThreads(4L))<br />
#[1] 1<br />
(default=.Internal(setNumMathThreads(1L))<br />
#[1] 1</p>
<p>set.seed(1001)<br />
y = rnorm(1000)<br />
x = matrix(rnorm(10000000), 1000, 10000)<br />
system.time(for (i in 1:ncol(x)) lm.fit(y = y, x = cbind(rep(1,NROW(x)),x[,i]))$coefficients[2])<br />
#   user  system elapsed<br />
#  3.643   0.090   3.748</p>
<p>system.time(cor(y,x)*sd(y)/apply(x,2L,sd))<br />
#   user  system elapsed<br />
#  1.448   0.044   1.497</p>
<p>system.time(replicate(10L, {mx=colSums(x); (y%*%x-mx*mean(y))/(colSums(x*x)-mx*mx/length(y))}))/10L<br />
#   user  system elapsed<br />
# 0.2331  0.0281  0.2612</p>
<p>system.time(replicate(10L, {mx=colSums(x); (colSums(y*x)-mx*mean(y))/(colSums(x*x)-mx*mx/length(y))}))/10L<br />
#   user  system elapsed<br />
# 0.3136  0.0548  0.3684</p>
<p>.Internal(setNumMathThreads(4L))<br />
#[1] 1</p>
<p>system.time(for (i in 1:ncol(x)) lm.fit(y = y, x = cbind(rep(1,NROW(x)),x[,i]))$coefficients[2])<br />
#   user  system elapsed<br />
#  3.670   0.000   3.671</p>
<p>system.time(cor(y,x)*sd(y)/apply(x,2L,sd))<br />
#   user  system elapsed<br />
#  1.262   0.028   1.290</p>
<p>system.time(replicate(10L, {mx=colSums(x); (y%*%x-mx*mean(y))/(colSums(x*x)-mx*mx/length(y))}))/10L<br />
#   user  system elapsed<br />
# 0.5100  0.0264  0.1997</p>
<p>system.time(replicate(10L, {mx=colSums(x); (colSums(y*x)-mx*mean(y))/(colSums(x*x)-mx*mx/length(y))}))/10L<br />
#   user  system elapsed<br />
# 0.8629  0.0567  0.2889</p>
<p>.Internal(setNumMathThreads(default))<br />
#[1] 4<br />
.Internal(setMaxNumMathThreads(default.max))<br />
#[1] 4<br />

This does not look very impressive...

youngtf

回复第15楼的 frank230：基因互作

micro@

回复第17楼的 youngtf：For GXG testing, normally the two main effects should be included into the null model. So my code for simple linear regression may not work directly. However, the same idea still applies. It's just a bit of more work.

youngtf

回复第18楼的 micro@：恩，正如您所说。谢谢提醒！