1#### 2#### h o o k e j e e v e s . R Hooke-Jeeves Minimization Algorithm 3#### 4 5## From: John C Nash <nashjc@uottawa.ca> 6## To: Martin Maechler <maechler@stat.math.ethz.ch>, Hans Werner Borchers 7## <hwborchers@googlemail.com> 8## Subject: Re: Rmpfr for optimization? Minor success. 9## Date: Tue, 5 Jun 2012 12:37:19 -0400 10 11## Changing to hjk routine was a bit easier to deal with. I found main changes 12## were to wrap output with as.numeric() to allow cat() to function. 13 14## I'll get no prizes for tidy code, but it is running an n=2 Chebyquad 15## minimization, and seems to be working on an n=6. This may be a good way to 16## encourage the sale of cpu power. 17 18## Best, JN 19 20hjkMpfr <- function(par, fn, control = list(), ...) { 21 ## Following fails when par is mpfr number JN120605 22 ## if (!is.numeric(par)) 23 ## stop("Argument 'par' must be a numeric vector.", call. = FALSE) 24 n <- length(par) 25 if (n == 1) 26 stop("For univariate functions use some different method.", call. = FALSE) 27 28 ##-- Control list handling ---------- 29 cntrl <- list(tol = 1.e-06, 30 maxfeval = Inf, # set to Inf if no limit wanted 31 maximize = FALSE, # set to TRUE for maximization 32 target = Inf, # set to Inf for no restriction 33 info = FALSE) # for printing interim information 34 nmsCo <- match.arg(names(control), choices = names(cntrl), several.ok = TRUE) 35 if (!is.null(names(control))) cntrl[nmsCo] <- control 36 37 tol <- cntrl$tol; 38 maxfeval <- cntrl$maxfeval 39 maximize <- cntrl$maximize 40 target <- cntrl$target 41 info <- cntrl$info 42 43 scale <- if (maximize) -1 else 1 44 fun <- match.fun(fn) 45 f <- function(x) scale * fun(x, ...) 46 47 ##-- Setting steps and stepsize ----- 48 nsteps <- floor(log2(1/tol)) # number of steps 49 steps <- 2^c(-(0:(nsteps-1))) # decreasing step size 50 dir <- diag(1, n, n) # orthogonal directions 51 52 x <- par # start point 53 fx <- f(x) # smallest value so far 54 fcount <- 1 # counts number of function calls 55 56 if (info) cat(sprintf("step nofc %-12s | %20s\n", 57 "fmin", "xpar")) 58 59 ##-- Start the main loop ------------ 60 ns <- 0 61 while (ns < nsteps && fcount < maxfeval && abs(fx) < target) { 62 ns <- ns + 1 63 hjs <- .hjsearch(x, f, steps[ns], dir, fcount, maxfeval, target) 64 x <- hjs$x 65 fx <- hjs$fx 66 ## found <- hjs$found 67 fcount <- fcount + hjs$finc 68 69 if (info) 70 cat(sprintf("%4d %5d %-12.7g | %-20.15g %-20.15g%s\n", 71 ns, fcount, as.numeric(fx/scale), 72 as.numeric(x[1]), as.numeric(x[2]), 73 if(n > 2)" ....")) 74 } 75 76 conv <- 77 if (fcount > maxfeval) { 78 warning("Function evaluation limit exceeded -- may not converge.") 79 FALSE 80 } else if (abs(fx) > target) { 81 warning("Function exceeds min/max value -- may not converge.") 82 FALSE 83 } else 84 TRUE 85 fx <- fx / scale # undo scaling 86 list(par = x, value = fx, 87 convergence = conv, feval = fcount, niter = ns) 88} 89 90## Search with a single scale ----------------------------- 91.hjsearch <- function(xb, f, h, dir, fcount, maxfeval, target) { 92 xc <- x <- xb 93 finc <- 0 94 hje <- .hjexplore(xb, xc, f, h, dir) 95 x <- hje$x 96 fx <- hje$fx 97 found <- hje$found 98 finc <- finc + hje$numf 99 100 ## Pattern move 101 while (found) { 102 d <- x-xb 103 xb <- x 104 xc <- x+d 105 fb <- fx 106 hje <- .hjexplore(xb, xc, f, h, dir, fb) 107 x <- hje$x 108 fx <- hje$fx 109 found <- hje$found 110 finc <- finc + hje$numf 111 112 if (!found) { # pattern move failed 113 hje <- .hjexplore(xb, xb, f, h, dir, fb) 114 x <- hje$x 115 fx <- hje$fx 116 found <- hje$found 117 finc <- finc + hje$numf 118 } 119 if (fcount + finc > maxfeval || abs(fx) > target) break 120 } 121 122 list(x = x, fx = fx, found=found, finc=finc) 123} 124 125## Exploratory move --------------------------------------- 126.hjexplore <- function(xb, xc, f, h, dir, fbold) { 127 n <- length(xb) 128 x <- xb 129 130 if (missing(fbold)) { 131 fb <- f(x) 132 numf <- 1 133 } else { 134 fb <- fbold 135 numf <- 0 136 } 137 138 fx <- fb 139 xt <- xc 140 found <- FALSE # do we find a better point ? 141 dirh <- h * dir 142 fbold <- fx 143 for (k in sample.int(n, n)) { # resample orthogonal directions 144 p <- xt + (d. <- dirh[, k]) 145 fp <- f(p) 146 numf <- numf + 1 147 148 if (fp >= fb) { 149 p <- xt - d. 150 fp <- f(p) 151 numf <- numf + 1 152 } 153 if (fp < fb) { 154 found <- TRUE 155 xt <- p 156 fb <- fp 157 } 158 } 159 if(found) { 160 x <- xt 161 fx <- fb 162 } 163 list(x = x, fx = fx, found=found, numf = numf) 164} 165