R/qEI.R

In DiceOptim: Kriging-Based Optimization for Computer Experiments

##' Analytical expression of the multipoint expected improvement (qEI)
##' criterion
##' 
##' Computes the multipoint expected improvement criterion.
##' 
##' 
##' @param x a matrix representing the set of input points (one row corresponds
##' to one point) where to evaluate the qEI criterion,
##' @param model an object of class \code{km},
##' @param plugin optional scalar: if provided, it replaces the minimum of the
##' current observations,
##' @param type "SK" or "UK" (by default), depending whether uncertainty
##' related to trend estimation has to be taken into account,
##' @param minimization logical specifying if EI is used in minimiziation or in
##' maximization,
##' @param fastCompute if TRUE, a fast approximation method based on a
##' semi-analytic formula is used (see [Marmin 2014] for details),
##' @param eps the value of \emph{epsilon} of the fast computation trick.
##' Relevant only if \code{fastComputation} is TRUE,
##' @param envir an optional environment specifying where to get intermediate
##' values calculated in \code{\link{qEI}}.
##' @return The multipoint Expected Improvement, defined as \deqn{qEI(X_{new})
##' := E[ ( min(Y(X)) - min(Y( X_{new} )) )_{+} | Y(X)=y(X)],}
##' 
##' where \eqn{X} is the current design of experiments, \eqn{ X_{new} } is a
##' new candidate design, and \eqn{Y} is a random process assumed to have
##' generated the objective function \eqn{y}.
##' @author Sebastien Marmin 
##' 
##' Clement Chevalier 
##' 
##' David Ginsbourger 
##' @seealso \code{\link{EI}}
##' @references
##' 
##' C. Chevalier and D. Ginsbourger (2014) Learning and Intelligent
##' Optimization - 7th International Conference, Lion 7, Catania, Italy,
##' January 7-11, 2013, Revised Selected Papers, chapter Fast computation of
##' the multipoint Expected Improvement with applications in batch selection,
##' pages 59-69, Springer.
##' 
##' D. Ginsbourger, R. Le Riche, L. Carraro (2007), A Multipoint Criterion for
##' Deterministic Parallel Global Optimization based on Kriging. The
##' International Conference on Non Convex Programming, 2007.
##' 
##' S. Marmin. Developpements pour l'evaluation et la maximisation du critere
##' d'amelioration esperee multipoint en optimisation globale (2014). Master's
##' thesis, Mines Saint-Etienne (France) and University of Bern (Switzerland).
##' 
##' D. Ginsbourger, R. Le Riche, and L. Carraro. Kriging is well-suited to
##' parallelize optimization (2010), 
##' In Lim Meng Hiot, Yew Soon Ong, Yoel
##' Tenne, and Chi-Keong Goh, editors, \emph{Computational Intelligence in
##' Expensive Optimization Problems}, Adaptation Learning and Optimization,
##' pages 131-162. Springer Berlin Heidelberg.
##' 
##' J. Mockus (1988), \emph{Bayesian Approach to Global Optimization}. Kluwer
##' academic publishers.
##' 
##' M. Schonlau (1997), \emph{Computer experiments and global optimization},
##' Ph.D. thesis, University of Waterloo.
##' @keywords models parallel optimization
##' 
##' @examples
##' 
##' \donttest{
##' set.seed(007)
##' 
##' # Monte-Carlo validation
##' 
##' # a 4-d, 81-points grid design, and the corresponding response
##' d <- 4; n <- 3^d
##' design <- do.call(expand.grid,rep(list(seq(0,1,length=3)),d))
##' names(design) <- paste("x",1:d,sep="")
##' y <- data.frame(apply(design, 1, hartman4))
##' names(y) <- "y"
##' 
##' # learning
##' model <- km(~1, design=design, response=y, control=list(trace=FALSE))
##' 
##' # pick up 10 points sampled from the 1-point expected improvement
##' q <- 10
##' X <- sampleFromEI(model,n=q)
##' # simulation of the minimum of the kriging random vector at X
##' t1 <- proc.time()
##' newdata <- as.data.frame(X)
##' colnames(newdata) <- colnames(model@@X)
##' 
##' krig  <- predict(object=model, newdata=newdata,type="UK",se.compute=TRUE, cov.compute=TRUE)
##' mk <- krig$mean
##' Sigma.q <- krig$cov
##' mychol <- chol(Sigma.q)
##' nsim <- 300000
##' white.noise <- rnorm(n=nsim*q)
##' minYsim <- apply(crossprod(mychol,matrix(white.noise,nrow=q)) + mk,2,min)
##' # simulation of the improvement (minimization)
##' qImprovement <- (min(model@@y)-minYsim)*((min(model@@y)-minYsim) > 0)
##' 
##' # empirical expectation of the improvement and confident interval (95%)
##' eiMC <- mean(qImprovement)
##' confInterv <- c(eiMC - 1.96*sd(qImprovement)*1/sqrt(nsim),eiMC + 1.96*sd(qImprovement)*1/sqrt(nsim))
##' 
##' # MC estimation of the qEI
##' print(eiMC) 
##' t2 <- proc.time()
##' # qEI with analytical formula
##' qEI(X,model,fastCompute= FALSE)
##' t3 <- proc.time()
##' # qEI with fast computation trick
##' qEI(X,model)
##' t4 <- proc.time()
##' t2-t1  # Time of MC computation
##' t3-t2  # Time of normal computation
##' t4-t3  # Time of fast computation
##' }
##' 
##' @importFrom mnormt dmnorm pmnorm
##' @export qEI
qEI <- function(x,model ,plugin = NULL, type = "UK", minimization = TRUE, fastCompute = TRUE, eps = 10^(-5), envir = NULL) {
  d <- model@d
  if (!is.matrix(x)) {
    x <- matrix(x,ncol = d)
  }
  xb <- rbind(model@X,x)
  nExp <- model@n
  x <- unique(round(xb,digits = 8))
  if ((nExp+1)>length(x[,1])) {return (0)}
  x <- matrix(x[(nExp+1):length(x[,1]),],ncol=d)
  if(is.null(plugin) &&  minimization) plugin <- min(model@y)
  if(is.null(plugin) && !minimization) plugin <- max(model@y)

  if (length(x[,1]) == 1) {
    return(EI(x=x,model=model,plugin=plugin,type=type,envir=envir))
  }

  krig  <- predict(object=model, newdata=x,type=type,se.compute=FALSE, cov.compute=TRUE,checkNames = FALSE)

  sigma <- krig$cov
  mu <- krig$mean
  
  q <- length(mu)
  
  pk <- first_term <- second_term <- rep(0,times=q)
  if (!fastCompute) {
    symetric_term <- matrix(0,q,q)
    non_symetric_term <- matrix(0,q,q)
  }
  for(k in 1:q){
    #covariance matrix of the vector Z^(k)
    Sigma_k <- covZk( sigma=sigma, index=k )
    
    #mean of the vector Z^(k)
    mu_k <- mu - mu[k] 
    mu_k[k] <- -mu[k]
    if(minimization) mu_k <- -mu_k
    
    b_k <- rep(0,times=q)
    b_k[k] <- -plugin
    if(minimization) b_k <- -b_k
    pk[k] <- pmnorm(x= b_k - mu_k  ,varcov=Sigma_k,maxpts=q*200)[1]
    
    first_term[k] <- (mu[k] - plugin)*pk[k]
    if(minimization) first_term[k] <- -first_term[k]
    
    if (fastCompute) {
	    second_term[k] <- 1/eps*(pmnorm(x= b_k + eps*Sigma_k[,k]- mu_k  ,varcov=Sigma_k,maxpts=q*200)[1] - pk[k])
    } else {
        for(i in 1:q){
        	non_symetric_term[k,i] <- Sigma_k[i,k]
	       if (i>=k) {
          mik <- mu_k[i]
	         sigma_ii_k <- Sigma_k[i,i]
	            bik <- b_k[i]
	            phi_ik <- dnorm(x=bik, mean=mik , sd= sqrt(sigma_ii_k))
	  
	            #need c.i^(k) and Sigma.i^(k)
	            cik <- get_cik(b = b_k , m = mu_k , sigma = Sigma_k , i=i  )
	            sigmaik <- get_sigmaik( sigma = Sigma_k , i=i )
	            Phi_ik <- pmnorm(x = cik  ,varcov=sigmaik,maxpts=(q-1)*200)[1]
	            symetric_term[k,i] <- phi_ik*Phi_ik
	          }
          }
        }
  }
  if (!fastCompute) {
    symetric_term <- symetric_term + t(symetric_term)
    diag(symetric_term) <- 0.5 * diag(symetric_term) 
    second_term <- sum(symetric_term*non_symetric_term)
  }  
  if (!is.null(envir)) {
    assign("pk",pk,envir=envir)
    if (fastCompute == FALSE) {
      assign("symetric_term",symetric_term,envir = envir)
    }
    assign("kriging.mean", mu, envir = envir)
    assign("kriging.cov", sigma, envir = envir)
    assign("Tinv.c", krig$Tinv.c, envir = envir)
    assign("c", krig$c, envir = envir)
  }
  somFinal <- sum(first_term,second_term)
  return( somFinal )
}



get_sigmaik <- function( sigma , i ){
  
  result <- sigma
  q <- nrow(sigma)
  
  for(u in 1:q){
    for(v in 1:q){
      if((u!=i) && (v!=i)){
        result[u,v] <- sigma[u,v] - sigma[u,i]*sigma[v,i]/sigma[i,i] 
      }else{
        result[u,v] <- 0
      }
    }
  }

  result <- result[-i,-i]

  return(result)
}


get_cik <- function( b , m , sigma , i  ){
  
  sigmai <- sigma[i,] / sigma[i,i]
  
  result <- (b - m) - ( b[i] - m[i] ) * sigmai
  result <- result[-i]
  
  return(result)
  
}



covZk <- function(sigma,index){
  
  result <- sigma
  q <- nrow(sigma)
  
  result[index,index] <- sigma[index,index]
  for(i in 1:q){
    if(i!=index){
      result[index,i] <- result[i,index] <- sigma[index,index] - sigma[index,i]
    }
  }
  for(i in 1:q){
    for(j in i:q){
      if((i!=index)&&(j!=index)){
        result[i,j] <- result[j,i] <- sigma[i,j] + sigma[index,index] - sigma[index,i] - sigma[index,j]
      }
    }
  }
  
  return(result)
  
}