[Rd] An example of very slow computation

From: John C Nash <nashjc_at_uottawa.ca>
Date: Wed, 17 Aug 2011 04:27:21 -0400


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")



R-devel_at_r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel Received on Wed 17 Aug 2011 - 08:30:25 GMT

This quarter's messages: by month, or sorted: [ by date ] [ by thread ] [ by subject ] [ by author ]

All messages

Archive maintained by Robert King, hosted by the discipline of statistics at the University of Newcastle, Australia.
Archive generated by hypermail 2.2.0, at Wed 17 Aug 2011 - 18:40:19 GMT.

Mailing list information is available at https://stat.ethz.ch/mailman/listinfo/r-devel. Please read the posting guide before posting to the list.

list of date sections of archive