[R] finding a faster way to run lm on rows of predictor matrix

[cypark at princeton.edu]

Hi, everyone.
I need to run lm with the same response vector but with varying predictor vectors. (i.e. 1 response vector on each individual 6,000 predictor vectors)
After looking through the R archive, I found roughly 3 methods that has been suggested.
Unfortunately, I need to run this task multiple times(~ 5,000 times) and would like to find a faster way than the existing methods.
All three methods I have bellow run on my machine timed with system.time 13~20 seconds.

The opposite direction of 6,000 response vectors and 1 predictor vectors, that is supported with lm runs really fast ~0.5 seconds. 
They are pretty much performing the same number of lm fits, so I was wondering if there was a faster way, before I try to code this up in c++.

thanks!!

## sample data ###
regress.y = rnorm(150)
predictors = matrix(rnorm(6000*150), ncol=150, nrow=6000)

## method 1 ##                                                                                                                                                                                                    
data.residuals = t(apply(predictors, 1, function(x)( lm(regress.y ~ -1 + as.vector(x))$residuals)))

user  system elapsed 
 15.076   0.048  15.128

## method 2 ##                                                                                                                                                             
data.residuals = matrix(rep(0, nrow(predictors) * ncol(predictors)), nrow=nrow(predictors), ncol=ncol(predictors) )

for( i in 1:nrow(predictors)){
    pred = as.vector(predictors[i,])
    data.residuals[i, ] = lm(regress.y ~ -1 + pred )$residuals
}

 user  system elapsed 
 13.841   0.012  13.854

## method 3 ##
library(nlme)                                                                                                                                                                                                   

all.data <- data.frame( y=rep(regress.y, nrow(predictors)), x=c(t(predictors)), g=gl(nrow(predictors), ncol(predictors)) )
all.fits <- lmList( y ~ x | g, data=all.data)
data.residuals = matrix( residuals(all.fits), nrow=nrow(predictors), ncol=ncol(predictors))

user  system elapsed 
 36.407   0.484  36.892 


## the opposite direction, supported by lm ### 
lm(t(predictors) ~ -1 + regress.y)$residuals

 user  system elapsed 
 0.500   0.120   0.613