[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)
```

# 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)
}

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

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 ]

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.