双循环线性模型?

use*_*124 1 r

当n增加时,R中双循环的速度非常慢.有没有办法提高for循环的速度?

    set.seed(1)
    n=1000

    y=rnorm(n)
    x1=rnorm(n)
    x2=rnorm(n)

    lm.ft=function(y,x1,x2)
      lm.fit(cbind(1,x1.bar,x2.bar), y)$coef

    res=array(,dim=c(1,3,n,n))
    for(i in 1:n)
      for(j in 1:n){
       x1.bar=x1-x1[i]
       x2.bar=x2-x2[j]
       res[,,i,j]=lm.ft(y,x1.bar,x2.bar)
      }
Run Code Online (Sandbox Code Playgroud)

Jor*_*eys 7

只是为了给你一个完整的答案:除了代码中的一些奇怪之处(如使用x1.barx2.bar内部lm.ft而不是x1x2),我的分析说:你到底想要实现什么?

如果我在你自己的代码上运行它:

Rprof("profile1.out")
for(i in 1:n)
  for(j in 1:n){
    x1.bar=x1-x1[i]
    x2.bar=x2-x2[j]
    res[,,i,j]=lm.ft(y,x1.bar,x2.bar)
  }
Rprof(NULL)
summaryRprof("profile1.out")
Run Code Online (Sandbox Code Playgroud)

我得到以下有趣的图片:

> summaryRprof("profile1.out")
$by.self
                self.time self.pct total.time total.pct
".Call"              0.96    22.86       0.96     22.86
"lm.fit"             0.92    21.90       4.08     97.14
...
"cbind"              0.22     5.24       0.22      5.24
...

$by.total
                total.time total.pct self.time self.pct
"lm.ft"               4.12     98.10      0.04     0.95
"lm.fit"              4.08     97.14      0.92    21.90
...
"cbind"               0.22      5.24      0.22     5.24
...
Run Code Online (Sandbox Code Playgroud)

98%的时间你只是适合模型.循环并不慢,你试图适应100万个模型的事实让你等待.你真的要重新思考你的问题.

如果这真的是你想要做的,那么优化你的功能将涉及摆脱lm.fit的开销和矢量化减法.节省约50%.

lm.ft=function(y,x1,x2)
  .Call(stats:::C_Cdqrls, cbind(1,x1,x2), y, tol=1e-7)$coef 

x1.bar <- outer(x1,x1,`-`)
x2.bar <- outer(x2,x2,`-`)
for(i in 1:n)
  for(j in 1:n){
    res[,,i,j]=lm.ft(y,x1.bar[,i],x2.bar[,j])
  }
Run Code Online (Sandbox Code Playgroud)