# Re: [Rd] An example of very slow computation

Date: Wed, 17 Aug 2011 21:24:03 +0000

```	eA <- eigen(A)
T <- eA\$vectors
r <- eA\$values
x0 <- c(0,100)
Tx0 <- T %*% x0

sol <- function(t) 100 - sum(T %*% diag(exp(r*t)) %*% Tx0)
pred <- sapply(dat[,1], sol)
-sum(dnorm(dat[,2], mean=pred, sd=sigma, log=TRUE))
```
}
This is much faster than using expm(A*t), but much slower than "by hand" calculations since we have to repeatedly do the diagonalization. An obvious advantage of this fuunction is that it applies to *any* linear system of ODEs for which the eigenvalues are real (and negative).

Ravi.

Assistant Professor,
Division of Geriatric Medicine and Gerontology School of Medicine Johns Hopkins University

Ph. (410) 502-2619

Yes, the culprit is the evaluation of expm(A*t). This is a lazy way of solving the system of ODEs, where you save analytic effort, but you pay for it dearly in terms of computational effort!

Even in this lazy approach, I am sure there ought to be faster ways to evaluating exponent of a matrix than that in "Matrix" package.

Ravi.

Assistant Professor,
Division of Geriatric Medicine and Gerontology School of Medicine Johns Hopkins University

Ph. (410) 502-2619

-----Original Message-----
From: r-devel-bounces_at_r-project.org [mailto:r-devel-bounces_at_r-project.org] On Behalf Of cberry_at_tajo.ucsd.edu Sent: Wednesday, August 17, 2011 1:08 PM To: r-devel_at_stat.math.ethz.ch
Subject: Re: [Rd] An example of very slow computation

John C Nash <nashjc_at_uottawa.ca> writes:

> 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.

Whew!

The time difference between 'expm( diag(10)+1 )' and 'expm( as( diag(10)+1, "dgeMatrix" ))' is a factor of 10 on my box.

Dunno 'bout the other factor of 100.

Chuck

> 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)
> A<-rbind(
> c(-k, k),
> c( k, -(k+k))
> )
> 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)
> A<-rbind(
> c(-k, k),
> c( k, -(k+k))
> )
> 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
> k2<-kvec
> k3<-kvec
> # 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)
> -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")
>
```--
Charles C. Berry                         cberry_at_tajo.ucsd.edu

______________________________________________
R-devel_at_r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel

______________________________________________
R-devel_at_r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel

______________________________________________
R-devel_at_r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel
```
Received on Wed 17 Aug 2011 - 21:27:40 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 Thu 18 Aug 2011 - 23:40:20 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.