Surrogate-based optimization of noisy black-box function

Description

Optimize (minimize) a noisy black-box function (i.e., a function which may not be differentiable, and may not execute deterministically). A b* tgp model is used as a surrogate for adaptive sampling via improvement (and other) statistics. Note that this function is intended as a skeleton to be tailored as required for a particular application

Usage

1
2
3
4
5
optim.step.tgp(f, rect, model = btgp, prev = NULL, X = NULL,
     Z = NULL, NN = 20 * length(rect), improv = c(1, 5), cands = c("lhs", "tdopt"),
     method = c("L-BFGS-B", "Nelder-Mead", "BFGS", "CG", "SANN",  "optimize"), ...)
optim.ptgpf(start, rect, tgp.obj,
     method=c("L-BFGS-B", "Nelder-Mead", "BFGS", "CG", "SANN", "optimize"))

Arguments

f

A function to be optimized, having only one free argument

rect

matrix indicating the domain of the argument of f over which an optimal should be searched; must have ncol(rect) = 2 and nrow agreeing with the argument of f indicating the dimension of the data. For 1-d data, a vector of length 2 is allowed

model

The b* regression model used as a surrogate for optimization; see btgp, and others, for more detail

prev

The output from a previous call to optim.step.tgp; this should be a list with entries as described the “Value” section below

X

data.frame, matrix, or vector of current inputs X, to be augmented

Z

Vector of current output responses Z of length equal to the leading dimension (rows) of X, i.e., length(Z) == nrow(X), to be augmented

NN

Number of candidate locations (XX) at which to sample from the improvement statistic

improv

Indicates the improv argument provided to a b* model function for sampling from the improvement statistic at the NN candidate locations (XX); see btgp, and others, for more detail

cands

The type of candidates (XX) at which samples from the improvement statistics are gathered. The default setting of "lhs" is recommended. However, a sequential treed D-optimal design can be used with "tdopt" for a more global exploration; see tgp.design for more details

method

A method from optim, or "optimize" which uses optimize as appropriate (when the input-space is 1-d)

...

Further arguments to the b* model function

start

A starting value for optimization of the MAP predictive (kriging) surface of a "tgp"-class object. A good starting value is the X or XX location found to be a minimum in the mean predictive surface contained in "tgp"-class object

tgp.obj

A "tgp"-class object that is the output of one of the b* functions: blm, btlm bgp, bgpllm, btgp, or btgpllm, as can be used by predict.tgp for optimizing on the MAP predictive (surrogate) kriging surface

Details

optim.step.tgp executes one step in a search for the global optimum (minimum) of a noisy function (Z~f(X)) in a bounded rectangle (rect). The procedure essentially fits a tgp model and samples from the posterior distribution of improvement statistics at NN+1 candidates locations. NN of the candidates come from cands, i.e., "lhs" or "tdopt", plus one which is the location of the minima found in a previous run via prev by using optim (with a particular method or optimize instead) on the MAP model predictive surface using the "tgp"-class object contained therein. The improv[2] with the the highest expected improvement are recommended for adding into the design on output.

optim.ptgpf is the subroutine used by optim.step.tgp to find optimize on the MAP (surrogate) predictive surface for the "tgp"-class object contained in prev.

Please see vignette("tgp2") for a detailed illustration

Value

The list return has the following components.

X

A matrix with nrow(rect) columns whose rows contain recommendations for input locations to add into the design

progress

A one-row data.frame indicating the the X-location and objective value of the current best guess of the solution to the (kriging) surrogate optimization along with the maximum values of the improvement statistic

obj

the "tgp"-class object output from the model function

Note

The ellipses (...) argument is used differently here, as compared to optim, and optimize. It allows further arguments to be passed to the b* model function, whereas for optim it would describe further (static) arguments to the function f to be optimized. If static arguments need to be set for f, then we recommend setting defaults via the formals of f

Author(s)

Robert B. Gramacy, rbgramacy@chicagobooth.edu, and Matt Taddy, taddy@chicagobooth.edu

References

Matthew Taddy, Herbert K.H. Lee, Genetha A. Gray, and Joshua D. Griffin. (2009) Bayesian guided pattern search for robust local optimization. Technometrics, to appear.

http://bobby.gramacy.com/r_packages/tgp

See Also

btgp, etc., optim, optimize, tgp.design, predict.tgp, dopt.gp

Examples

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
## optimize the simple exponential function 
f <- function(x) { exp2d.Z(x)$Z }

## create the initial design with D-optimal candidates
rect <- rbind(c(-2,6), c(-2,6))
Xcand <- lhs(500, rect)
X <- dopt.gp(50, X=NULL, Xcand)$XX
Z <- f(X)

## do 10 rounds of adaptive sampling
out <- progress <- NULL
for(i in 1:10) {

  ## get recommendations for the next point to sample
  out <- optim.step.tgp(f, X=X, Z=Z, rect=rect, prev=out)

  ## add in the inputs, and newly sampled outputs
  X <- rbind(X, out$X)
  Z <- c(Z, f(out$X))

  ## keep track of progress and best optimum
  progress <- rbind(progress, out$progress)
  print(progress[i,])
}

## plot the progress so far
par(mfrow=c(2,2))
plot(out$obj, layout="surf")
plot(out$obj, layout="as", as="improv")
matplot(progress[,1:nrow(rect)], main="optim results",
        xlab="rounds", ylab="x[,1:2]", type="l", lwd=2)
plot(log(progress$improv), type="l", main="max log improv",
     xlab="rounds", ylab="max log(improv)")