[R] Full enumeration, how can this be vectorized

[Uwe Ligges]

Daniel Hoppe wrote:
> 
> Hi all,
> 
> I am currently using the code snippet below to minimize a function which I
> wasn't able to minimize properly with optim (cliffy surface, constraints,
> poles). Obviously this iteration is probably the poorest way of implementing
> such a function and it is expectedly slow. I've already written some
> vectorized functions, but I'm afraid that I'm missing the "big picture" in
> this case, I don't know how to start. Probably I could build a 3 x
> (length(x) * length(y) * length(z)) matrix holding each possible parameter
> combinations. With my current step-size and value range that would be
> slightly more than 2,000,000 triples, so that should not be a problem. Then
> comes the part I'm totally clueless about, that is applying the target
> function to each of the triples in the matrix in an efficient way. From the
> resulting vector I could then easily take the minimum value index and get
> the parameters from the parameter matrix.
> 
> Could someone kindly give me a hint how this could be implemented?


Without looking closer, I think this is a very nice example for a piece
of code that should be implemented in either C or Fortran (on the one
hand it is easy to implement, on the other hand its "S interpreted"
execution takes quite a long time). OK, quite a few possible
improvements of the R code are obvious, but I don't think this is the
best point to start ...

Nevertheless, the code looks like you are searching on a grid. Are you
really sure whether there is no method in optim() which performs better?

Uwe Ligges



> Thanks a lot,
> 
> Daniel
> 
> Code Snippet:
> fullEnumeration <- function(par, stepSize = c(.03, .03, .03))
> {
>    res <- Inf
> 
>    best.x <- -1
>    best.y <- -1
>    best.z <- -1
> 
>    stepx <- stepSize[1]
>    stepy <- stepSize[2]
>    stepz <- stepSize[3]
>    x <- seq(.01, 2, by=stepx)
>    y <- seq(.01, 2, by=stepy)
>    z <- seq(.01, 15, by=stepz)
> 
>    for (i in 1:length(x))
>     for (j in 1:length(y))
>         for (k in 1:length(z))
>         {
>                 # Pass the current best result to the function being minimized to allow
> for early termination
>             newRes <- fun(c(x[i], y[j], z[k]), par, res)
>             if (!is.na(newRes) && newRes < res)
>             {
>                 res = newRes
>                 best.x <- x[i]
>                 best.y <- y[j]
>                 best.z <- z[k]
>             }
> 
>         }
>    return (c(
>         best.x,
>         best.y,
>         best.z))
> 
> }
> 
> fun <- function(x, par, currentBestValue=Inf)
> {
>   x = x[1]
>   y = x[2]
>   z = x[3]
>   TMax <- length(par)
> 
>   result <- 0
>   for (myT in 1:TMax)
>   {
>     result <- result +
>     (
>     ((x * z) / (y-1) * (1 - (z / (z + myT))^(y-1) ))
>     - par[myT]
>     )^2
> 
>     # Allow for short-cut evaluation if running in complete enumeration mode
> 
>     if (is.na(result) || currentBestValue < result)
>     {
>         return (Inf)
>     }
> 
>   }
>   return (result)
> 
> }
> 
> -.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
> r-help mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
> Send "info", "help", or "[un]subscribe"
> (in the "body", not the subject !)  To: r-help-request at stat.math.ethz.ch
> _._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._
-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
r-help mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
Send "info", "help", or "[un]subscribe"
(in the "body", not the subject !)  To: r-help-request at stat.math.ethz.ch
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._