R: Optimization of homogenously noisy functions based on Kriging

noisy.optimizer {DiceOptim}

R Documentation

Optimization of homogenously noisy functions based on Kriging

Description

Sequential optimization of kriging-based criterion conditional on noisy observations, with model update after each evaluation. Eight criteria are proposed to choose the next observation: random search, sequential parameter optimization (SPO), reinterpolation, Expected Improvement (EI) with plugin, Expected Quantile Improvement (EQI), quantile minimization, Augmented Expected Improvement (AEI) and Approximate Knowledge Gradient (AKG). The criterion optimization is based on the package rgenoud.

Usage

noisy.optimizer(
  optim.crit,
  optim.param = NULL,
  model,
  n.ite,
  noise.var = NULL,
  funnoise,
  lower,
  upper,
  parinit = NULL,
  control = NULL,
  CovReEstimate = TRUE,
  NoiseReEstimate = FALSE,
  nugget.LB = 1e-05,
  estim.model = NULL,
  type = "UK"
)

Arguments

`optim.crit`	String defining the criterion to be optimized at each iteration. Possible values are: "random.search", "SPO", "reinterpolation", "EI.plugin", "EQI", "min.quantile", "AEI", "AKG".
`optim.param`	List of parameters for the chosen criterion. For "EI.plugin": optim.param$plugin.type is a string defining which plugin is to be used. Possible values are "ytilde", "quantile" and "other". If "quantile" is chosen, optim.param$quantile defines the quantile level. If "other" is chosen, optim.param$plugin directly sets the plugin value. For "EQI": optim.param$quantile defines the quantile level. If not provided, default value is 0.9. For "min.quantile": optim.param$quantile defines the quantile level. If not provided, default value is 0.1. For "AEI": optim.param$quantile defines the quantile level to choose the current best point. If not provided, default value is 0.75.
`model`	a Kriging model of "km" class
`n.ite`	Number of iterations
`noise.var`	Noise variance (scalar). If noiseReEstimate=TRUE, it is an initial guess for the unknown variance (used in optimization).
`funnoise`	objective (noisy) function
`lower`	vector containing the lower bounds of the variables to be optimized over
`upper`	vector containing the upper bounds of the variables to be optimized over
`parinit`	optional vector of initial values for the variables to be optimized over
`control`	optional list of control parameters for optimization. One can control `"pop.size"` (default : [N=32^dim for dim<6 and N=32dim otherwise]]), `"max.generations"` (N), `"wait.generations"` (2) and `"BFGSburnin"` (0) of function `"genoud"` (see `genoud`). Numbers into brackets are the default values
`CovReEstimate`	optional boolean specfiying if the covariance parameters should be re-estimated at every iteration (default value = TRUE)
`NoiseReEstimate`	optional boolean specfiying if the noise variance should be re-estimated at every iteration (default value = FALSE)
`nugget.LB`	optional scalar of minimal value for the estimated noise variance. Default value is 1e-5.
`estim.model`	optional kriging model of "km" class with homogeneous nugget effect (no noise.var). Required if noise variance is reestimated and the initial "model" has heterogenenous noise variances.
`type`	"SK" or "UK" for Kriging with known or estimated trend

Value

A list with components:

`model`	the current (last) kriging model of "km" class
`best.x`	The best design found
`best.y`	The objective function value at best.x
`best.index`	The index of best.x in the design of experiments
`history.x`	The added observations
`history.y`	The added observation values
`history.hyperparam`	The history of the kriging parameters
`estim.model`	If noiseReEstimate=TRUE, the current (last) kriging model of "km" class for estimating the noise variance.
`history.noise.var`	If noiseReEstimate=TRUE, the history of the noise variance estimate.

Author(s)

Victor Picheny

References

V. Picheny and D. Ginsbourger (2013), Noisy kriging-based optimization methods: A unified implementation within the DiceOptim package, Computational Statistics & Data Analysis

Examples


##########################################################################
### EXAMPLE 1: 3 OPTIMIZATION STEPS USING EQI WITH KNOWN NOISE         ###
### AND KNOWN COVARIANCE PARAMETERS FOR THE BRANIN FUNCTION            ###
##########################################################################

set.seed(10)
library(DiceDesign)
# Set test problem parameters
doe.size <- 9
dim <- 2
test.function <- get("branin2")
lower <- rep(0,1,dim)
upper <- rep(1,1,dim)
noise.var <- 0.1

# Build noisy simulator
funnoise <- function(x)
{     f.new <- test.function(x) + sqrt(noise.var)*rnorm(n=1)
      return(f.new)}

# Generate DOE and response
doe <- as.data.frame(lhsDesign(doe.size, dim)$design)
y.tilde <- funnoise(doe)

# Create kriging model
model <- km(y~1, design=doe, response=data.frame(y=y.tilde),
     covtype="gauss", noise.var=rep(noise.var,1,doe.size), 
     lower=rep(.1,dim), upper=rep(1,dim), control=list(trace=FALSE))

# Optimisation with noisy.optimizer (n.ite can be increased)
n.ite <- 2 
optim.param <- list()
optim.param$quantile <- .9
optim.result <- noisy.optimizer(optim.crit="EQI", optim.param=optim.param, model=model,
		n.ite=n.ite, noise.var=noise.var, funnoise=funnoise, lower=lower, upper=upper,
		NoiseReEstimate=FALSE, CovReEstimate=FALSE)

new.model <- optim.result$model
best.x    <- optim.result$best.x
new.doe   <- optim.result$history.x

## Not run: 
##### DRAW RESULTS #####
# Compute actual function on a grid
n.grid <- 12
x.grid <- y.grid <- seq(0,1,length=n.grid)
design.grid <- expand.grid(x.grid, y.grid)
names(design.grid) <- c("V1","V2")
nt <- nrow(design.grid)
func.grid <- rep(0,1,nt)

for (i in 1:nt)
{ func.grid[i] <- test.function(x=design.grid[i,])}

# Compute initial and final kriging on a grid
pred <- predict(object=model, newdata=design.grid, type="UK", checkNames = FALSE)
mk.grid1 <- pred$m
sk.grid1 <- pred$sd

pred <- predict(object=new.model, newdata=design.grid, type="UK", checkNames = FALSE)
mk.grid2 <- pred$m
sk.grid2 <- pred$sd

# Plot initial kriging mean
z.grid <- matrix(mk.grid1, n.grid, n.grid)
filled.contour(x.grid,y.grid, z.grid, nlevels=50, color = topo.colors,
plot.axes = {title("Initial kriging mean");
points(model@X[,1],model@X[,2],pch=17,col="black"); 
axis(1); axis(2)})

# Plot initial kriging variance
z.grid <- matrix(sk.grid1^2, n.grid, n.grid)
filled.contour(x.grid,y.grid, z.grid, nlevels=50, color = topo.colors,
plot.axes = {title("Initial kriging variance");
points(model@X[,1],model@X[,2],pch=17,col="black"); 
axis(1); axis(2)})

# Plot final kriging mean
z.grid <- matrix(mk.grid2, n.grid, n.grid)
filled.contour(x.grid,y.grid, z.grid, nlevels=50, color = topo.colors,
plot.axes = {title("Final kriging mean");
points(new.model@X[,1],new.model@X[,2],pch=17,col="black"); 
axis(1); axis(2)})

# Plot final kriging variance
z.grid <- matrix(sk.grid2^2, n.grid, n.grid)
filled.contour(x.grid,y.grid, z.grid, nlevels=50, color = topo.colors,
plot.axes = {title("Final kriging variance");
points(new.model@X[,1],new.model@X[,2],pch=17,col="black"); 
axis(1); axis(2)})

# Plot actual function and observations
z.grid <- matrix(func.grid, n.grid, n.grid)
tit <- "Actual function - Black: initial points; red: added points"
filled.contour(x.grid,y.grid, z.grid, nlevels=50, color = topo.colors,
plot.axes = {title(tit);points(model@X[,1],model@X[,2],pch=17,col="black"); 
points(new.doe[1,],new.doe[2,],pch=15,col="red");
axis(1); axis(2)})

## End(Not run)

##########################################################################
### EXAMPLE 2: 3 OPTIMIZATION STEPS USING EQI WITH UNKNOWN NOISE       ###
### AND UNKNOWN COVARIANCE PARAMETERS FOR THE BRANIN FUNCTION          ###
##########################################################################
# Same initial model and parameters as for example 1
n.ite <- 2 # May be changed to a larger value 
res <- noisy.optimizer(optim.crit="min.quantile", 
optim.param=list(type="quantile",quantile=0.01),
model=model, n.ite=n.ite, noise.var=noise.var, funnoise=funnoise, 
lower=lower, upper=upper, 
control=list(print.level=0),CovReEstimate=TRUE, NoiseReEstimate=TRUE)

# Plot actual function and observations
plot(model@X[,1], model@X[,2], pch=17,xlim=c(0,1),ylim=c(0,1))
points(res$history.x[1,], res$history.x[2,], col="blue")

# Restart: requires the output estim.model of the previous run
# to deal with potential repetitions
res2 <- noisy.optimizer(optim.crit="min.quantile", 
optim.param=list(type="quantile",quantile=0.01), 
model=res$model, n.ite=n.ite, noise.var=noise.var, funnoise=funnoise, 
lower=lower, upper=upper, estim.model=res$estim.model,
control=list(print.level=0),CovReEstimate=TRUE, NoiseReEstimate=TRUE)

# Plot new observations
points(res2$history.x[1,], res2$history.x[2,], col="red")

[Package DiceOptim version 2.1.1 Index]