very fast OLS regression?

Dear R experts:

I just tried some simple test that told me that hand computing the OLS
coefficients is about 3-10 times as fast as using the built-in lm()
function.   (code included below.)  Most of the time, I do not care,
because I like the convenience, and I presume some of the time goes
into saving a lot of stuff that I may or may not need.  But when I do
want to learn the properties of an estimator whose input contains a
regression, I do care about speed.

What is the recommended fastest way to get regression coefficients in
R?  (Is Gentlemen's weighted-least-squares algorithm implemented in a
low-level C form somewhere?  that one was always lightning fast for
me.)

regards,

/ivo



bybuiltin = function( y, x )   coef(lm( y ~ x -1 ));

byhand = function( y, x ) {
  xy<-t(x)%*%y;
  xxi<- solve(t(x)%*%x)
  b<-as.vector(xxi%*%xy)
  ## I will need these later, too:
  ## res<-y-as.vector(x%*%b)
  ## soa[i]<-b[2]
  ## sigmas[i]<-sd(res)
  b;
}


MC=500;
N=10000;


set.seed(0);
x= matrix( rnorm(N*MC), nrow=N, ncol=MC );
y= matrix( rnorm(N*MC), nrow=N, ncol=MC );

ptm = proc.time()
for (mc in 1:MC) byhand(y[,mc],x[,mc]);
cat("By hand took ", proc.time()-ptm, "\n");

ptm = proc.time()
for (mc in 1:MC) bybuiltin(y[,mc],x[,mc]);
cat("By built-in took ", proc.time()-ptm, "\n");

check the following options:

ols1 <- function (y, x) {
? ?coef(lm(y ~ x - 1))
}

ols2 <- function (y, x) {
? ?xy <- t(x)%*%y
? ?xxi <- solve(t(x)%*%x)
? ?b <- as.vector(xxi%*%xy)
? ?b
}

ols3 <- function (y, x) {
? ?XtX <- crossprod(x)
? ?Xty <- crossprod(x, y)
? ?solve(XtX, Xty)
}

ols4 <- function (y, x) {
? ?lm.fit(x, y)$coefficients
}

byhand(y, x)

ans <- byhand.qr(y, x)
ans

y - c(x %*% ans)

require(MASS)
ans.min <- c(ginv(x) %*% y)
ans.min

Dear R experts:
 
 I just tried some simple test that told me that hand computing the OLS
 coefficients is about 3-10 times as fast as using the built-in lm()
 function.   (code included below.)  Most of the time, I do not care,
 because I like the convenience, and I presume some of the time goes
 into saving a lot of stuff that I may or may not need.  But when I do
 want to learn the properties of an estimator whose input contains a
 regression, I do care about speed.
 
 What is the recommended fastest way to get regression coefficients in
 R?  (Is Gentlemen's weighted-least-squares algorithm implemented in a
 low-level C form somewhere?  that one was always lightning fast for
 me.)
 
 regards,
 
 /ivo
 
 
 
 bybuiltin = function( y, x )   coef(lm( y ~ x -1 ));
 
 byhand = function( y, x ) {
   xy<-t(x)%*%y;
   xxi<- solve(t(x)%*%x)
   b<-as.vector(xxi%*%xy)
   ## I will need these later, too:
   ## res<-y-as.vector(x%*%b)
   ## soa[i]<-b[2]
   ## sigmas[i]<-sd(res)
   b;
 }
 
 
 MC=500;
 N=10000;
 
 
 set.seed(0);
 x= matrix( rnorm(N*MC), nrow=N, ncol=MC );
 y= matrix( rnorm(N*MC), nrow=N, ncol=MC );
 
 ptm = proc.time()
 for (mc in 1:MC) byhand(y[,mc],x[,mc]);
 cat("By hand took ", proc.time()-ptm, "\n");
 
 ptm = proc.time()
 for (mc in 1:MC) bybuiltin(y[,mc],x[,mc]);
 cat("By built-in took ", proc.time()-ptm, "\n");
 

check the following options:

ols1 <- function (y, x) {
? ?coef(lm(y ~ x - 1))
}

ols2 <- function (y, x) {
? ?xy <- t(x)%*%y
? ?xxi <- solve(t(x)%*%x)
? ?b <- as.vector(xxi%*%xy)
? ?b
}

ols3 <- function (y, x) {
? ?XtX <- crossprod(x)
? ?Xty <- crossprod(x, y)
? ?solve(XtX, Xty)
}

ols4 <- function (y, x) {
? ?lm.fit(x, y)$coefficients
}

lm is slow because it has to set up the design matrix (X) each time. See
 ?model.matrix and ?model.matrix.lm  for how to do this once 
separately from
 lm (and then use lm.fit).
 
 I am far from an expert on numerical algebra, but I'm pretty sure your
 comments below are incorrect in the sense that "naive" methods are **not**
 universally better then efficient numerical algebra methods 
using,say, QR
 decompositions. It will depend very strongly on the size and specific 
nature
 of the problem. Big Oh Complexity statements ( O(n) or O(n^2), etc.) 
are,
 after all, asymptotic. There is also typically a tradeoff between
 computational accuracy (especially for ill-conditioned matrices)  and 
speed
 which your remarks seem to neglect.
 
 -- Bert
 
 
 Bert Gunter
 Genentech Nonclinical Biostatistics
 650-467-7374
 
 -----Original Message-----
 From: r-help-bounces at r-project.org [ On
 Behalf Of ivo welch
 Sent: Wednesday, March 25, 2009 2:30 PM
 To: Dimitris Rizopoulos
 Cc: r-help
 Subject: Re: [R] very fast OLS regression?
 
 thanks, dimitris.  I also added Bill Dunlap's "solve(qr(x),y)"
 function as ols5.   here is what I get in terms of speed on a Mac Pro:
 
 ols1 6.779 3.591 10.37 0 0
 ols2 0.515 0.21 0.725 0 0
 ols3 0.576 0.403 0.971 0 0
 ols4 1.143 1.251 2.395 0 0
 ols5 0.683 0.565 1.248 0 0
 
 so the naive matrix operations are fastest.  I would have thought that
 alternatives to the naive stuff I learned in my linear algebra course
 would be quicker.   still, ols3 and ols5 are competitive.  the
 built-in lm() is really problematic.  is ols3 (or perhaps even ols5)
 preferable in terms of accuracy?  I think I can deal with 20% speed
 slow-down (but not with a factor 10 speed slow-down).
 
 regards,
 
 /iaw
 
 
 On Wed, Mar 25, 2009 at 5:11 PM, Dimitris Rizopoulos
 <d.rizopoulos at erasmusmc.nl> wrote:

 > check the following options:
 >
 > ols1 <- function (y, x) {
 > ? ?coef(lm(y ~ x - 1))
 > }
 >
 > ols2 <- function (y, x) {
 > ? ?xy <- t(x)%*%y
 > ? ?xxi <- solve(t(x)%*%x)
 > ? ?b <- as.vector(xxi%*%xy)
 > ? ?b
 > }
 >
 > ols3 <- function (y, x) {
 > ? ?XtX <- crossprod(x)
 > ? ?Xty <- crossprod(x, y)
 > ? ?solve(XtX, Xty)
 > }
 >
 > ols4 <- function (y, x) {
 > ? ?lm.fit(x, y)$coefficients
 > }
 >

Dear R experts:

I just tried some simple test that told me that hand computing the OLS
coefficients is about 3-10 times as fast as using the built-in lm()
function.   (code included below.)  Most of the time, I do not care,
because I like the convenience, and I presume some of the time goes
into saving a lot of stuff that I may or may not need.  But when I do
want to learn the properties of an estimator whose input contains a
regression, I do care about speed.

What is the recommended fastest way to get regression coefficients in
R?  (Is Gentlemen's weighted-least-squares algorithm implemented in a
low-level C form somewhere?  that one was always lightning fast for
me.)

regards,

/ivo



bybuiltin = function( y, x )   coef(lm( y ~ x -1 ));

byhand = function( y, x ) {
  xy<-t(x)%*%y;
  xxi<- solve(t(x)%*%x)
  b<-as.vector(xxi%*%xy)
  ## I will need these later, too:
  ## res<-y-as.vector(x%*%b)
  ## soa[i]<-b[2]
  ## sigmas[i]<-sd(res)
  b;
}


MC=500;
N=10000;


set.seed(0);
x= matrix( rnorm(N*MC), nrow=N, ncol=MC );
y= matrix( rnorm(N*MC), nrow=N, ncol=MC );

ptm = proc.time()
for (mc in 1:MC) byhand(y[,mc],x[,mc]);
cat("By hand took ", proc.time()-ptm, "\n");

ptm = proc.time()
for (mc in 1:MC) bybuiltin(y[,mc],x[,mc]);
cat("By built-in took ", proc.time()-ptm, "\n");

thanks, dimitris.  I also added Bill Dunlap's "solve(qr(x),y)"
function as ols5.   here is what I get in terms of speed on a Mac Pro:

ols1 6.779 3.591 10.37 0 0
ols2 0.515 0.21 0.725 0 0
ols3 0.576 0.403 0.971 0 0
ols4 1.143 1.251 2.395 0 0
ols5 0.683 0.565 1.248 0 0

so the naive matrix operations are fastest.  I would have thought that
alternatives to the naive stuff I learned in my linear algebra course
would be quicker.   still, ols3 and ols5 are competitive.  the
built-in lm() is really problematic.  is ols3 (or perhaps even ols5)
preferable in terms of accuracy?  I think I can deal with 20% speed
slow-down (but not with a factor 10 speed slow-down).

regards,

/iaw


On Wed, Mar 25, 2009 at 5:11 PM, Dimitris Rizopoulos
<d.rizopoulos at erasmusmc.nl> wrote:

check the following options:

ols1 <- function (y, x) {
? ?coef(lm(y ~ x - 1))
}

ols2 <- function (y, x) {
? ?xy <- t(x)%*%y
? ?xxi <- solve(t(x)%*%x)
? ?b <- as.vector(xxi%*%xy)
? ?b
}

ols3 <- function (y, x) {
? ?XtX <- crossprod(x)
? ?Xty <- crossprod(x, y)
? ?solve(XtX, Xty)
}

ols4 <- function (y, x) {
? ?lm.fit(x, y)$coefficients
}

Yes, Bert.  Any least-squares solution that forms X'X and then inverts it is 
not to be recommended.  If X is nearly rank-deficient, then X'X will be more 
strongly so.  The QR decomposition approach in my byhand.qr() function is 
reliable and fast.

check the following options:

ols1 <- function (y, x) {
     coef(lm(y ~ x - 1))
}

ols2 <- function (y, x) {
     xy <- t(x)%*%y
     xxi <- solve(t(x)%*%x)
     b <- as.vector(xxi%*%xy)
     b
}

ols3 <- function (y, x) {
     XtX <- crossprod(x)
     Xty <- crossprod(x, y)
     solve(XtX, Xty)
}

ols4 <- function (y, x) {
     lm.fit(x, y)$coefficients
}

# check timings
MC <- 500
N <- 10000

set.seed(0)
x <- matrix(rnorm(N*MC), N, MC)
y <- matrix(rnorm(N*MC), N, MC)

invisible({gc(); gc(); gc()})
system.time(for (mc in 1:MC) ols1(y[, mc], x[, mc]))

invisible({gc(); gc(); gc()})
system.time(for (mc in 1:MC) ols2(y[, mc], x[, mc]))

invisible({gc(); gc(); gc()})
system.time(for (mc in 1:MC) ols3(y[, mc], x[, mc]))

invisible({gc(); gc(); gc()})
system.time(for (mc in 1:MC) ols4(y[, mc], x[, mc, drop = FALSE]))

On Wed, 2009-03-25 at 16:28 -0400, ivo welch wrote:

Dear R experts:

I just tried some simple test that told me that hand computing the OLS
coefficients is about 3-10 times as fast as using the built-in lm()
function. ? (code included below.) ?Most of the time, I do not care,
because I like the convenience, and I presume some of the time goes
into saving a lot of stuff that I may or may not need. ?But when I do
want to learn the properties of an estimator whose input contains a
regression, I do care about speed.

What is the recommended fastest way to get regression coefficients in
R? ?(Is Gentlemen's weighted-least-squares algorithm implemented in a
low-level C form somewhere? ?that one was always lightning fast for
me.)

No one has yet mentioned Doug Bates' article in R News on this topic,
which compares timings of various methods for least squares
computations.

Douglas Bates. Least squares calculations in R. R News, 4(1):17-20, June
2004.

A Cholesky decomposition solution proved fastest in base R code, with an
even faster version developed using sparse matrices and the Matrix
package.

You can find Doug's article here:

http://cran.r-project.org/doc/Rnews/Rnews_2004-1.pdf

regards,

/ivo



bybuiltin = function( y, x ) ? coef(lm( y ~ x -1 ));

byhand = function( y, x ) {
? xy<-t(x)%*%y;
? xxi<- solve(t(x)%*%x)
? b<-as.vector(xxi%*%xy)
? ## I will need these later, too:
? ## res<-y-as.vector(x%*%b)
? ## soa[i]<-b[2]
? ## sigmas[i]<-sd(res)
? b;
}


MC=500;
N=10000;


set.seed(0);
x= matrix( rnorm(N*MC), nrow=N, ncol=MC );
y= matrix( rnorm(N*MC), nrow=N, ncol=MC );

ptm = proc.time()
for (mc in 1:MC) byhand(y[,mc],x[,mc]);
cat("By hand took ", proc.time()-ptm, "\n");

ptm = proc.time()
for (mc in 1:MC) bybuiltin(y[,mc],x[,mc]);
cat("By built-in took ", proc.time()-ptm, "\n");

--
%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%
?Dr. Gavin Simpson ? ? ? ? ? ? [t] +44 (0)20 7679 0522
?ECRC, UCL Geography, ? ? ? ? ?[f] +44 (0)20 7679 0565
?Pearson Building, ? ? ? ? ? ? [e] gavin.simpsonATNOSPAMucl.ac.uk
?Gower Street, London ? ? ? ? ?[w] http://www.ucl.ac.uk/~ucfagls/
?UK. WC1E 6BT. ? ? ? ? ? ? ? ? [w] http://www.freshwaters.org.uk
%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%