[R] Full enumeration, how can this be vectorized
Uwe Ligges
ligges at statistik.uni-dortmund.de
Mon Nov 25 22:14:39 CET 2002
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
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._
More information about the R-help
mailing list