如何计算响应矩阵每列的最小但快速的线性回归？

Question

我想在不使用"lm"的情况下计算R中的普通最小二乘(OLS)估计,这有几个原因.首先,考虑到数据大小在我的情况下是一个问题,"lm"还计算了许多我不需要的东西(例如拟合值).其次,我希望能够在R中自己实现OLS,然后再使用其他语言(例如,在C中使用GSL).

您可能知道,模型是:Y = Xb + E; 用E~N(0,sigma ^ 2).如下所述,b是具有2个参数的向量,即平均值(b0)和另一个系数(b1).最后,对于我将要进行的每个线性回归,我想要b1(效应大小)的估计,其标准误差,sigma ^ 2的估计(残差方差)和R ^ 2(确定系数).

这是我的数据.我有N个样本(例如个体,N~ = 100).对于每个样本,我有Y个输出(响应变量,Y~ = 10 ^ 3)和X个点(解释变量,X~ = 10 ^ 6).我想分开处理Y输出,即.我想推出Y线性回归:一个用于输出1,一个用于输出2,等等.此外,我想一个一个地使用解释变量:对于输出1,我想在第1点回归,然后在第2点回归,然后......终于在X点了.(我希望它很清楚......!)

这是我的R代码,通过矩阵代数检查"lm"与计算OLS估算的速度.

首先,我模拟虚拟数据:

nb.samples <-  10  # N
nb.points <- 1000  # X
x <- matrix(data=replicate(nb.samples,sample(x=0:2,size=nb.points, replace=T)),
            nrow=nb.points, ncol=nb.samples, byrow=F,
            dimnames=list(points=paste("p",seq(1,nb.points),sep=""),
              samples=paste("s",seq(1,nb.samples),sep="")))
nb.outputs <- 10  # Y
y <- matrix(data=replicate(nb.outputs,rnorm(nb.samples)),
            nrow=nb.samples, ncol=nb.outputs, byrow=T,
            dimnames=list(samples=paste("s",seq(1,nb.samples),sep=""),
              outputs=paste("out",seq(1,nb.outputs),sep="")))

以下是我自己使用的函数:

GetResFromCustomLinReg <- function(Y, xi){ # both Y and xi are N-dim vectors
  n <- length(Y)
  X <- cbind(rep(1,n), xi)  #
  p <- 1      # nb of explanatory variables, besides the mean
  r <- p + 1  # rank of X: nb of indepdt explanatory variables
  inv.XtX <- solve(t(X) %*% X)
  beta.hat <- inv.XtX %*% t(X) %*% Y
  Y.hat <- X %*% beta.hat
  E.hat <- Y - Y.hat
  E2.hat <- (t(E.hat) %*% E.hat)
  sigma2.hat <- (E2.hat / (n - r))[1,1]
  var.covar.beta.hat <- sigma2.hat * inv.XtX
  se.beta.hat <- t(t(sqrt(diag(var.covar.beta.hat))))
  Y.bar <- mean(Y)
  R2 <- 1 - (E2.hat) / (t(Y-Y.bar) %*% (Y-Y.bar))
  return(c(beta.hat[2], se.beta.hat[2], sigma2.hat, R2))
}

这是我使用内置"lm"的代码:

res.bi.all <- apply(x, 1, function(xi){lm(y ~ xi)})

这是我的自定义OLS代码:

res.cm.all <- apply(x, 1, function(xi){apply(y, 2, GetResFromCustomLinReg, xi)})

当我使用上面给出的值运行此示例时,我得到:

> system.time( res.bi.all <- apply(x, 1, function(xi){lm(y ~ xi)}) )
   user  system elapsed
  2.526   0.000   2.528
> system.time( res.cm.all <- apply(x, 1, function(xi){apply(y, 2, GetResFromCustomLinReg, xi)}) )
   user  system elapsed
  4.561   0.000   4.561

(当然,增加N,X和Y时会变得更糟.)

当然,"lm"具有分别"响应矩阵"(y~xi)的每一列"自动"拟合的良好特性,而我必须使用"应用"Y次(对于每个yi~xi).但这是我的代码更慢的唯一原因吗？你们其中一个人知道如何改善吗？

(对不起这么长的问题,但我真的试图提供一个最小但又全面的例子.)

> sessionInfo()
R version 2.12.2 (2011-02-25)
Platform: x86_64-redhat-linux-gnu (64-bit)

Answer 1

看看CRAN上RcppArmadillo包中的fastLm()功能.在此之前还有类似的RcppGSL - 但你可能想要基于犰狳的解决方案.我在旧的演示文稿中有一些幻灯片(在带有R的HPC上)显示速度增加.fastLm()

还要注意帮助页面中关于更好的"旋转"方法的提示,而不是X'X的直反,这可能与简并模型矩阵有关.