R/AEI.grad.R

In DiceOptim: Kriging-Based Optimization for Computer Experiments

##' AEI's Gradient
##' 
##' Analytical gradient of the Augmented Expected Improvement (AEI) criterion.
##' 
##' 
##' @param x the input vector at which one wants to evaluate the criterion
##' @param model a Kriging model of "km" class
##' @param new.noise.var the (scalar) noise variance of the new observation.
##' @param y.min The kriging predictor at the current best point (point with
##' smallest kriging quantile). If not provided, this quantity is evaluated.
##' @param type Kriging type: "SK" or "UK"
##' @param envir environment for inheriting intermediate calculations from AEI
##' @return Gradient of the Augmented Expected Improvement
##' @author Victor Picheny 
##'
##' David Ginsbourger 
##' @examples
##'
##' set.seed(421)
##' 
##' # Set test problem parameters
##' doe.size <- 12
##' dim <- 2
##' test.function <- get("branin2")
##' lower <- rep(0,1,dim)
##' upper <- rep(1,1,dim)
##' noise.var <- 0.2
##' 
##' # Generate DOE and response
##' doe <- as.data.frame(matrix(runif(doe.size*dim),doe.size))
##' y.tilde <- rep(0, 1, doe.size)
##' for (i in 1:doe.size)  {
##' y.tilde[i] <- test.function(doe[i,]) + sqrt(noise.var)*rnorm(n=1)
##' }
##' y.tilde <- as.numeric(y.tilde)
##' 
##' # 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))
##' 
##' # Compute actual function and criterion on a grid
##' n.grid <- 8 # Change to 21 for a nicer picture
##' x.grid <- y.grid <- seq(0,1,length=n.grid)
##' design.grid <- expand.grid(x.grid, y.grid)
##' nt <- nrow(design.grid)
##' 
##' crit.grid <- apply(design.grid, 1, AEI, model=model, new.noise.var=noise.var)
##' crit.grad <- t(apply(design.grid, 1, AEI.grad, model=model, new.noise.var=noise.var))
##' 
##' z.grid <- matrix(crit.grid, n.grid, n.grid)
##' contour(x.grid,y.grid, z.grid, 30)
##' title("AEI and its gradient")
##' points(model@@X[,1],model@@X[,2],pch=17,col="blue")
##'
##' for (i in 1:nt)
##' {
##'  x <- design.grid[i,]
##'  suppressWarnings(arrows(x$Var1,x$Var2, x$Var1+crit.grad[i,1]*.6,x$Var2+crit.grad[i,2]*.6,
##' length=0.04,code=2,col="orange",lwd=2))
##' }
##' 
##' @export AEI.grad
AEI.grad <- function(x, model, new.noise.var=0, y.min=NULL, type = "UK", envir=NULL){

  ########## 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)
  
  # Compute y.min if missing
  if (is.null(y.min))
  {
    pred <- predict(object=model, newdata=model@X, type=type, checkNames = FALSE)
    mk <- pred$mean
    sk <- pred$sd
    qk <- mk + qnorm(0.75)*sk
    y.min <- mk[which.min(qk)]
  }
  
  T <- model@T
  X <- model@X
  z <- model@z
  u   <- model@M
  covStruct <- model@covariance
  
  if (is.null(envir))
  {  predx <- predict.km(object=model, newdata=newdata, type=type, checkNames = FALSE)
     mk  <- predx$mean
     sk  <- predx$sd
     c   <- predx$c
     v   <- predx$Tinv.c
     xcr <- (y.min - mk)/sk
     xcr.prob <- pnorm(xcr)
     xcr.dens <- dnorm(xcr)
  } else
  { toget <- matrix(c("xcr", "xcr.prob", "xcr.dens", "c", "Tinv.c", "mk","sk"),1,7)
    apply(toget, 2, get, envir=envir)
    xcr      <- envir$xcr
    xcr.prob <- envir$xcr.prob
    xcr.dens <- envir$xcr.dens
    c        <- envir$c
    v        <- envir$Tinv.c
    mk       <- envir$mk
    sk       <- envir$sk
  }
  F.newdata <- model.matrix(model@trend.formula, data=newdata)
  
  ###########################################################################
  
  if (sk < sqrt(model@covariance@sd2)/1e6)
  { aei.grad.val <- rep(0,d)
  } else 
  { # AEI
    ei.val <- (y.min - mk) * xcr.prob + sk * xcr.dens
    pen <- (1- sqrt(new.noise.var)/sqrt(new.noise.var + sk^2))
    aei.val <- ei.val * pen
    
    # EI gradient
    # 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)
    mk.grad <- t(W)%*%z + t(model@trend.coef%*%f.deltax)
    
    if (type=="UK")
    { tuuinv <- solve(t(u)%*%u)
      sk2.grad <-  t( -2*t(v)%*%W +
                      2*(F.newdata - t(v)%*%u )%*% tuuinv %*%
                      ( f.deltax - t(t(W)%*%u) ))
    } else
    { sk2.grad <-  t( -2*t(v)%*%W)
    }
    sk.grad <- sk2.grad / (2*sk)
    
    # Compute gradient of EI
    ei.grad <- - mk.grad * xcr.prob + sk.grad * xcr.dens
    
    # Penalization gradient
    pen.grad <- sqrt(new.noise.var)/2*sk2.grad*(new.noise.var + sk^2)^(-3/2)
    
    # AEI gradient
    aei.grad.val <- ei.grad*pen + ei.val*pen.grad
    
#     pen.grad <- - ei.val * (1 - sqrt(new.noise.var))*sk*sk.grad*(new.noise.var + sk^2)^(-3/2)
#     aei.grad.val <- ei.grad*pen + ei.val^2*pen.grad
  }
  return(aei.grad.val)
}