[R] linear model benchmarking

[ivo welch]

I cleaned up my old benchmarking code and added checks for missing
data to compare various ways of finding OLS regression coefficients.
I thought I would share this for others.  the long and short of it is
that I would recommend

           ols.crossprod = function (y, x) {
              x <- as.matrix(x)
              ok <- (!is.na(y))&(!is.na(rowSums(x)))
              y <- y[ok]; x <- subset(x, ok)
              x <- cbind( 1, x)

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

for fast and stable coefficients.  (yes, stable using double
precision, even though not as stable as lm().  it works just fine with
X variables that have >99.99% correlation with as few as 100
observations.  if your situation is worse than this, you probably have
an error in your data---or you are looking for the Higgs Boson.)

I added the code below.  feel free to ignore.

/iaw


################################################################
###
### code to test alternatives how fast OLS coefficients can be obtained.
### including tests to exclude missing observations where necessary.
###
### for a more informed article (and person), see Bates, 'Least Squares
### Calculations in R', Rnews 2004-1. the code here does not test his sparse
### matrix examples, or his geMatrix/poMatrix examples.
###
### Basic Results: for the examples that I tried, typical relative time
### factors of the algorithms were about
###
###      lm  lmfit solve crossprod  cholesky  (special-case 2vars)
###      1.0  0.5   0.3    0.15      0.17         0.1
###
### there was no advantage to cholesky, so you may as well use the simpler
### crossprod.
###
### I was also interested in algorithm scaling N and K.  yes, there were
### some changes in the factors across algorithms, but the general pattern
### wasn't too different.  for the cholesky decomposition,
###
###            N=1000   N=10000   N=100000   N=200000
###    K=1       1.0       7        80       160
###    K=10      2.5      26
###    K=50     16
###    K=100    43        70
###    K=200   140
###
### some of this may well be swap/memory access, etc.  roughly speaking, we
### scale ten-times N takes twice as long. 10 and ten-times K takes 25 times
### as long.
###
### of course, ols.crossprod and ols.cholesky are not as stable as ols.lm,
### but they are still amazingly stable, given the default double precision.
### even with a correlation of 0.999999(!) between the two final columns,
### they still produce exactly the same result as ols.lm with 1000
### observations.  frankly, the ill-conditioning worry is overblown with
### most real-world data.  if you really have data THIS bad, you should
### already know it; and you probably just have some measurement errors in
### your observations, and your regression is giving you garbage either way.
###
### if I made the R core decision, I would switch away from lm()'s default
### method, and make it a special option.  my guess is that it is status-quo
### bias that keeps the current method.  or, at least I would say loudly in
### the R docs that for common use, here is a much faster method...
###
################################################################

MC <- 100
N <- 1000
K <- 10
SD <- 1e-3

ols <- list(
            ols.lm = function (y, x) { coef(lm(y ~ x)) },

            ols.lmfit = function (y, x) {
              x <- as.matrix(x)
              ok <- (!is.na(y))&(!is.na(rowSums(x)))
              y <- y[ok]; x <- subset(x, ok)
              x <- as.matrix(cbind( 1, x))
              lm.fit(x, y)$coefficients
            },

            ols.solve = function (y, x) {
              x <- as.matrix(x)
              ok <- (!is.na(y))&(!is.na(rowSums(x)))
              y <- y[ok]; x <- subset(x, ok)
              x <- cbind(1, x)
              xy <- t(x)%*%y
              xxi <- solve(t(x)%*%x)
              b <- as.vector(xxi%*%xy)
              b
            },

            ols.crossprod = function (y, x) {
              x <- as.matrix(x)
              ok <- (!is.na(y))&(!is.na(rowSums(x)))
              y <- y[ok]; x <- subset(x, ok)
              x <- cbind( 1, x)

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

            ols.cholesky = function (y, x) {
              x <- as.matrix(x)
              ok <- (!is.na(y))&(!is.na(rowSums(x)))
              y <- y[ok]; x <- subset(x, ok)
              x <- cbind( 1, x)

              ch <- chol( crossprod(x) )
              backsolve(ch, forwardsolve(ch, crossprod(x,y),
upper=TRUE, trans=TRUE))
            }

            )

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

cat("N=", N, "K=", K, "  (MC=", MC, ")")
if (K>1) {
  sum.cor <- 0
  for (mc in 1:MC) {
    x[,K,mc] <- x[,K-1,mc]+rnorm(N, sd=SD)
    sum.cor <- sum.cor + cor(x[,K,mc], x[,K-1,mc], use="pair")
  }
  options(digit=10)
  cat( " sd=", SD, "The bad corr= 1+", sum.cor/MC-1 )
} else {
  ols$ols.xy.Kis1 <- function(y, x) {
    ok <- (!is.na(y))&(!is.na(x))
    y <- y[ok]; x <- x[ok]
    b <- cov(y,x)/var(x)
    c(mean(y) - b*mean(x), b)
  }
}
cat("\n")

## make sure we do the NaN checks
x[1,1,1] <- NA
y[2,2] <- NA

options(digits=2)
time.factor.first <- first.checksum <- NA
for (ctr in 1:length(ols)) {
  invisible({gc(); gc(); gc()})
  coefsum <- 0.0
  v <- system.time(for (mc in 1:MC) coefsum <- coefsum +
ols[[ctr]](y[,mc], x[,,mc]))
  thiscoefsum <- sum(abs(coefsum))
  thistime <- sum(as.numeric(v)[1:2])/MC
  if (is.na(time.factor.first)) time.factor.first <- thistime
  if (is.na(first.checksum)) first.checksum <- thiscoefsum
  cat(sprintf("%13.13s", names(ols[ctr])), " Checksum = ", thiscoefsum, "(",
      sprintf("%.8f", abs(thiscoefsum - first.checksum)/MC), ")",
      "Time =", sprintf("%.3f", 100*thistime), "msecs.  Factor=",
thistime/time.factor.first, "\n")
}