Using lapply to maximize a function along a grid

Question

I rewrote the question to make it reproducible. Suppose I want to maximize the function exp(alpha+eta+gamma) for alpha,eta,gamma along a grid of my own choice. I have done this using for-loops but I want to make use of apply-functions to speed up the procedure. Here's what I have done (eta and gamma is here being held fixed).

eta=0.11
gamma=0.06
alpha=0.5    
alpha_vals=seq(0.1,1,by=0.1)
eta_vals=eta
gamma_vals=gamma


ml_temp=-Inf

lapply(alpha_vals,function(alpha_v){
  lapply(eta_vals,function(eta_v){
    lapply(gamma_vals,function(gamma_v){
      temp=exp(alpha_v+eta_v+gamma_v)
      if (temp >= ml_temp) {
        ml_temp=temp
        mle_matrix=c(alpha_v,eta_v,gamma_v)
      }  
    })  
  })    
})

Outputting mle_matrix I get 0 0 0, so something is clearly not working. Any help is appreciated.

Why are you using `invisible()` at the end of each function? This means you are returning `NULL` as your output at each step. — Andrie, Aug 22 '12 at 08:49
Oh, had some help by someone, and he added those. I thought they had a purpose, but I guess I should remove them. — Stefan Hansen, Aug 22 '12 at 08:50
I removed the `invisible()` but my output is still wrong. In this case I get `mle_matrix = 0 0 0` which clearly is functioning properly. — Stefan Hansen, Aug 22 '12 at 08:52
Can you try to make your code [reproducible](http://stackoverflow.com/questions/5963269/how-to-make-a-great-r-reproducible-example)? At the moment I can't run your code, since `loglike_in` isn't defined. — Andrie, Aug 22 '12 at 08:57
@Andrie: See the edit. I changed the function to be `exp(alpha+eta+gamma)` instead. — Stefan Hansen, Aug 22 '12 at 09:05
If you change each `lapply` to `sapply` you'll get a matrix as the result. Then you have to remember that your assignment of `temp` will have no effect, since it only gets assigned in the scope of the innermost function. — Andrie, Aug 22 '12 at 09:06
That's true, but it's not the matrix i'm interested in. I'm interested in assigning the value of `(alpha,eta,gamma)` to `mle_matrix` that gives the largest value of `exp(alpha+eta+gamma)`. Changing `lapply` to `sapply` still returns `0 0 0` in `mle_matrix`. — Stefan Hansen, Aug 22 '12 at 09:11
Do I assume correctly that ultimately you want to give vectors for `eta_vals` and `gamma_vals` as well? — cbeleites unhappy with SX, Aug 22 '12 at 09:12
It looks like you were trying to use a "closure". Note that you're not setting or accessing global variables here, instead you're defining local variables within the context of the inner-most function. For example, `ml_temp` is only a copy of the outer-most `ml_temp`. Also, you might consider to use `mapply` instead of the nested `lapply`s. — fotNelton, Aug 22 '12 at 09:24

score 5 · Accepted Answer · answered Aug 22 '12 at 09:24

The easiest way is to use expand.grid() and apply()

toTest <- expand.grid(
    alpha = seq(0.1, 1, by = 0.1), 
    eta = seq(0.1, 1, by = 0.1), 
    gamma = seq(0.1, 1, by = 0.1))
ml <- apply(toTest, 1, function(x){
  exp(sum(x))
})
toTest[which.max(ml), ]

cbeleites unhappy with SX · Answer 2 · 2012-08-22T10:36:24.167

I assume that your target function is vectorized. exp is, if possible you should try to write your real target function vectorized, too.
For the moment, I assume that you can hold the whole of the results in memory. If not, you'll have to lapply of one or more of your parameters, but there's probably no need to loop over all of them.

Now here's the suggestion:

search.pts <- expand.grid (alpha_v = alpha_vals, 
                           gamma_v = gamma_vals, 
                           eta_v   = eta_vals)
target.val <- exp (search.pts$alpha_v + search.pts$gamma_v + search.pts$eta_v)
solution <- which.max (target.val)

search.pts [solution,]

or maybe return

list (params = search.pts [solution,], value = target.val [solution])

Here's a benchmark on vectorized calculation vs. apply for each of the three parameter vectors having 100 values:

microbenchmark (
    vectorized = exp (search.pts$alpha_v + search.pts$gamma_v + search.pts$eta_v),
    apply      = apply (search.pts, 1, function (x) exp (sum (x))), 
    times = 10)

## Unit: milliseconds
##         expr        min        lq     median          uq       max
## 1      apply 9569.52277 9687.8617 9940.53103 10140.13430 11413.508
## 2 vectorized   44.37456   45.3286   46.75505    67.38978   314.196

Thank you. Unfortunately it is not vectorized and i'm unable to make it so. I will use the suggestion above, but thanks anyway! — Stefan Hansen, Aug 22 '12 at 09:37
@StefanHansen: I added a benchmark of vectorized vs. apply. If you find your grid search unbearably slow, you may need to reconsider a vectorized version of the target function. — cbeleites unhappy with SX, Aug 22 '12 at 10:37

Using lapply to maximize a function along a grid

2 Answers2