R/EI.grad.R

In DiceOptim: Kriging-Based Optimization for Computer Experiments

##' Analytical gradient of the Expected Improvement criterion
##' 
##' Computes the gradient of the Expected Improvement at the current location.
##' The current minimum of the observations can be replaced by an arbitrary
##' value (plugin), which is usefull in particular in noisy frameworks.
##' 
##' 
##' @param x a vector representing the input for which one wishes to calculate
##' \code{\link{EI}}.
##' @param model an object of class \code{\link[DiceKriging]{km}}.
##' @param plugin optional scalar: if provided, it replaces the minimum of the
##' current observations,
##' @param type Kriging type: "SK" or "UK"
##' @param minimization logical specifying if EI is used in minimiziation or in
##' maximization,
##' @param envir an optional environment specifying where to get intermediate
##' values calculated in \code{\link{EI}}.
##' @param proxy an optional Boolean, if TRUE EI is replaced by the kriging mean (to minimize)
##' @return The gradient of the expected improvement criterion with respect to
##' x.  Returns 0 at design points (where the gradient does not exist).
##' @author David Ginsbourger 
##' 
##' Olivier Roustant 
##' 
##' Victor Picheny 
##' @seealso \code{\link{EI}}
##' @references
##' 
##' D. Ginsbourger (2009), \emph{Multiples metamodeles pour l'approximation et
##' l'optimisation de fonctions numeriques multivariables}, Ph.D. thesis, Ecole
##' Nationale Superieure des Mines de Saint-Etienne, 2009.
##' 
##' J. Mockus (1988), \emph{Bayesian Approach to Global Optimization}. Kluwer
##' academic publishers.
##' 
##' T.J. Santner, B.J. Williams, and W.J. Notz (2003), \emph{The design and
##' analysis of computer experiments}, Springer.
##' 
##' M. Schonlau (1997), \emph{Computer experiments and global optimization},
##' Ph.D. thesis, University of Waterloo.
##' @keywords models optimize
##' @examples
##' 
##' set.seed(123)
##' # a 9-points factorial design, and the corresponding response
##' d <- 2; n <- 9
##' design.fact <- expand.grid(seq(0,1,length=3), seq(0,1,length=3)) 
##' names(design.fact)<-c("x1", "x2")
##' design.fact <- data.frame(design.fact) 
##' names(design.fact)<-c("x1", "x2")
##' response.branin <- apply(design.fact, 1, branin)
##' response.branin <- data.frame(response.branin) 
##' names(response.branin) <- "y" 
##' 
##' # model identification
##' fitted.model1 <- km(~1, design=design.fact, response=response.branin, 
##' covtype="gauss", control=list(pop.size=50,trace=FALSE), parinit=c(0.5, 0.5))
##' 
##' # graphics
##' n.grid <- 9  # Increase to 50 for a nicer picture
##' x.grid <- y.grid <- seq(0,1,length=n.grid)
##' design.grid <- expand.grid(x.grid, y.grid)
##' #response.grid <- apply(design.grid, 1, branin)
##' EI.grid <- apply(design.grid, 1, EI,fitted.model1)
##' #EI.grid <- apply(design.grid, 1, EI.plot,fitted.model1, gr=TRUE)
##' 
##' z.grid <- matrix(EI.grid, n.grid, n.grid)
##' 
##' contour(x.grid,y.grid,z.grid,20)
##' title("Expected Improvement for the Branin function known at 9 points")
##' points(design.fact[,1], design.fact[,2], pch=17, col="blue")
##' 
##' # graphics
##' n.gridx <- 5  # increase to 15 for nicer picture
##' n.gridy <- 5  # increase to 15 for nicer picture
##' x.grid2 <- seq(0,1,length=n.gridx) 
##' y.grid2 <- seq(0,1,length=n.gridy) 
##' design.grid2 <- expand.grid(x.grid2, y.grid2)
##' 
##' EI.envir <- new.env()	
##' 	environment(EI) <- environment(EI.grad) <- EI.envir 
##'
##' for(i in seq(1, nrow(design.grid2)) )
##' {
##' 	x <- design.grid2[i,]
##' 	ei <- EI(x, model=fitted.model1, envir=EI.envir)
##' 	eigrad <- EI.grad(x , model=fitted.model1, envir=EI.envir)
##' 	if(!(is.null(ei)))
##' 	{
##' 	suppressWarnings(arrows(x$Var1,x$Var2,
##' 	x$Var1 + eigrad[1]*2.2*10e-5, x$Var2 + eigrad[2]*2.2*10e-5, 
##' 	length = 0.04, code=2, col="orange", lwd=2))
##' 	}
##' }
##' 
##' @export EI.grad
EI.grad <- function(x, model, plugin=NULL, type="UK", minimization = TRUE, envir=NULL, proxy=FALSE){ 

  ########################################################################################
  if (is.null(plugin)){ 
    if (minimization) {
      plugin <- min(model@y)
    } else {
      plugin <- -max(model@y)
    }
  }
  m <- plugin
  
  ########################################################################################
  # Convert x in proper format(s)
  d <- length(x)
  if (d != model@d){ stop("x does not have the right size") }
  newdata.num <- as.numeric(x)
  newdata <- data.frame(t(newdata.num))
  colnames(newdata) = colnames(model@X)
  ########################################################################################
  # Get quantities related to the model
  T <- model@T
  X <- model@X
  z <- model@z
  u <- model@M
  covStruct <- model@covariance
  
  # Get quantities related to the prediction
  if (is.null(envir))
  {  
    predx <- predict(object=model, newdata=newdata, type=type, checkNames = FALSE,se.compute=TRUE,cov.compute=FALSE)
    kriging.mean <- predx$mean
    if(!minimization) kriging.mean <- -kriging.mean
    kriging.sd <- predx$sd
    v <- predx$Tinv.c
    c <- predx$c
    
    xcr <- (m - kriging.mean)/kriging.sd
    xcr.prob <- pnorm(xcr)
    xcr.dens <- dnorm(xcr)    
  } else
  {  # If uploaded through "envir", no prediction computation is necessary 
    toget <- matrix(c("xcr", "xcr.prob", "xcr.dens", "kriging.sd", "c", "Tinv.c"), 1, 6)
    apply(toget, 2, get, envir=envir)
    xcr        <- envir$xcr
    xcr.prob   <- envir$xcr.prob
    xcr.dens   <- envir$xcr.dens
    kriging.sd <- envir$kriging.sd
    c          <- envir$c
    v          <- envir$Tinv.c
  }

  F.newdata <- model.matrix(model@trend.formula, data=newdata)

  ########################################################################################
  # Pursue calculation only if standard deviation is non-zero
  if ( kriging.sd/sqrt(model@covariance@sd2) < 1e-06) 
  { ei.grad <- rep(0,d)
  } else 
  { # Compute derivatives of the covariance and trend functions
    dc <- covVector.dx(x=newdata.num, X=X, object=covStruct, c=c)  
    f.deltax <- trend.deltax(x=newdata.num, model=model)
  
    # Compute gradients of the kriging mean and variance
    W <- backsolve(t(T), dc, upper.tri=FALSE)
    kriging.mean.grad <- t(W)%*%z + t(model@trend.coef%*%f.deltax)
    if (!minimization) kriging.mean.grad <- -kriging.mean.grad
    
    if (proxy) {
      ei.grad <- - kriging.mean.grad
    } else {
      
      if (type=="UK")
      { tuuinv <- solve(t(u)%*%u)
      kriging.sd2.grad <-  t( -2*t(v)%*%W +
                                2*(F.newdata - t(v)%*%u )%*% tuuinv %*%
                                (f.deltax - t(t(W)%*%u) ))
      } else
      { kriging.sd2.grad <-  t( -2*t(v)%*%W) }
      
      kriging.sd.grad <- kriging.sd2.grad / (2*kriging.sd)
      
      # Compute gradient of EI
      ei.grad <- - kriging.mean.grad * xcr.prob + kriging.sd.grad * xcr.dens
    }
  }

  ########################################################################################  
  return(ei.grad)
}