[Rd] An example of very slow computation

[John C Nash]

This message is about a curious difference in timing between two ways of computing the
same function. One uses expm, so is expected to be a bit slower, but "a bit" turned out to
be a factor of >1000. The code is below. We would be grateful if anyone can point out any
egregious bad practice in our code, or enlighten us on why one approach is so much slower
than the other. The problem arose in an activity to develop guidelines for nonlinear
modeling in ecology (at NCEAS, Santa Barbara, with worldwide participants), and we will be
trying to include suggestions of how to prepare problems like this for efficient and
effective solution. The code for nlogL was the "original" from the worker who supplied the
problem.

Best,

John Nash

--------------------------------------------------------------------------------------

cat("mineral-timing.R == benchmark MIN functions in R\n")
#  J C Nash July 31, 2011

require("microbenchmark")
require("numDeriv")
library(Matrix)
library(optimx)
# dat<-read.table('min.dat', skip=3, header=FALSE)
# t<-dat[,1]
t <- c(0.77,  1.69,  2.69,  3.67,  4.69,  5.71,  7.94,  9.67, 11.77, 17.77,
 23.77, 32.77, 40.73, 47.75, 54.90, 62.81, 72.88, 98.77, 125.92, 160.19,
 191.15, 223.78, 287.70, 340.01, 340.95, 342.01)

# y<-dat[,2] # ?? tidy up
 y<- c(1.396, 3.784, 5.948, 7.717, 9.077, 10.100, 11.263, 11.856, 12.251, 12.699,
12.869, 13.048, 13.222, 13.347, 13.507, 13.628, 13.804, 14.087, 14.185, 14.351,
14.458, 14.756, 15.262, 15.703, 15.703, 15.703)


ones<-rep(1,length(t))
theta<-c(-2,-2,-2,-2)


nlogL<-function(theta){
  k<-exp(theta[1:3])
  sigma<-exp(theta[4])
  A<-rbind(
  c(-k[1], k[2]),
  c( k[1], -(k[2]+k[3]))
  )
  x0<-c(0,100)
  sol<-function(t)100-sum(expm(A*t)%*%x0)
  pred<-sapply(dat[,1],sol)
  -sum(dnorm(dat[,2],mean=pred,sd=sigma, log=TRUE))
}

getpred<-function(theta, t){
  k<-exp(theta[1:3])
  sigma<-exp(theta[4])
  A<-rbind(
  c(-k[1], k[2]),
  c( k[1], -(k[2]+k[3]))
  )
  x0<-c(0,100)
  sol<-function(tt)100-sum(expm(A*tt)%*%x0)
  pred<-sapply(t,sol)
}

Mpred <- function(theta) {
# WARNING: assumes t global
kvec<-exp(theta[1:3])
k1<-kvec[1]
k2<-kvec[2]
k3<-kvec[3]
#   MIN problem terbuthylazene disappearance
    z<-k1+k2+k3
    y<-z*z-4*k1*k3
    l1<-0.5*(-z+sqrt(y))
    l2<-0.5*(-z-sqrt(y))
    val<-100*(1-((k1+k2+l2)*exp(l2*t)-(k1+k2+l1)*exp(l1*t))/(l2-l1))
} # val should be a vector if t is a vector

negll <- function(theta){
# non expm version JN 110731
  pred<-Mpred(theta)
  sigma<-exp(theta[4])
  -sum(dnorm(dat[,2],mean=pred,sd=sigma, log=TRUE))
}

theta<-rep(-2,4)
fand<-nlogL(theta)
fsim<-negll(theta)
cat("Check fn vals: expm =",fand,"   simple=",fsim,"  diff=",fand-fsim,"\n")

cat("time the function in expm form\n")
tnlogL<-microbenchmark(nlogL(theta), times=100L)
tnlogL

cat("time the function in simpler form\n")
tnegll<-microbenchmark(negll(theta), times=100L)
tnegll

ftimes<-data.frame(texpm=tnlogL$time, tsimp=tnegll$time)
# ftimes


boxplot(log(ftimes))
title("Log times in nanoseconds for matrix exponential and simple MIN fn")