These R files contains an analytical implementation of the criterion “IECI” proposed for Optimization under unknown constraints (2010) by Robert Gramacy and Herbert Lee.

This implementation is based on DiceKriging R package (CRAN version):

if (!is.element("DiceKriging",installed.packages())) install.packages("DiceKriging")

# More robust (and expensive) EI optimization
if (!is.element("rgenoud",installed.packages())) install.packages("rgenoud")

# Also load DiceView for easy plot
if (!is.element("DiceView",installed.packages())) install.packages("DiceView")

# IECI will need randtoolbox, some KrigInv routines
if (!is.element("randtoolbox",installed.packages())) install.packages("randtoolbox")
if (!is.element("KrigInv",installed.packages())) install.packages("KrigInv")

Objective function

fun=function(x){-(1-1/2*(sin(12*x)/(1+x)+2*cos(7*x)*x^5+0.7))}

plot(fun,type='l',ylim=c(-1.2,0.2),ylab="y")

plot of chunk unnamed-chunk-2

X=data.frame(X=matrix(c(.0,.33,.737,1),ncol=1))
y=fun(X)

Expected Improvement (EI)

par(mfrow=c(2,1))

library(DiceKriging)
set.seed(123)
kmi <- km(design=X,response=y,control=list(trace=FALSE),optim.method = 'gen')

library(DiceView)
## Loading required package: DiceEval
## Loading required package: rgl
## Warning in rgl.init(initValue, onlyNULL): RGL: unable to open X11
## display
## Warning: 'rgl_init' failed, running with rgl.useNULL = TRUE
sectionview.km(kmi,col_surf='blue',col_points='blue',xlim=c(0,1),ylim=c(-1.1,0),title = "",Xname = "x",yname="Y(x)")
plot(fun,type='l',add=T)
abline(h=min(y),col='blue',lty=2)
text(labels="m",x=0.05,y=min(y)+0.05,col='blue')
legend(bg = rgb(1,1,1,.4),0.3,0,legend = c("(unknown) objective function","conditionnal gaussian process","observations"),col=c('black','blue','blue'),lty = c(1,2,0),pch = c(-1,-1,20))

source("https://raw.githubusercontent.com/IRSN/IECI/master/EI.R")
.xx=seq(f=0,t=1,l=501)
plot(.xx,EI(.xx,kmi),col='blue',type='l',xlab="x",ylab="EI(x)")

plot of chunk unnamed-chunk-3

Expected Conditional Improvement (ECI) by Monte Carlo estimation (! slow !)

par(mfrow=c(2,1))

xn=0.4

sectionview.km(kmi,col_surf='blue',col_points='blue',xlim=c(xn-0.2,xn+0.2),ylim=c(-1.1,-.7),title = "",Xname = "x")
plot(fun,type='l',add=T)
abline(h=min(y),col='blue',lty=2)
text(labels="m",x=0.35,y=min(y)-0.05,col='blue')

abline(v=xn,col='red')
text(labels="xn",x=xn+0.01,y=-0.75,col='red')

Yn=predict.km(kmi,newdata=xn,type="UK",checkNames = F)
.yy = seq(f=Yn$lower,t=Yn$upper,l=11)
#lines(x=dnorm(.xxx,yn$mean,yn$sd)/100+xn,y=.xxx,col='red')
points(x=rep(xn,length(.yy)),y=.yy,col=rgb(1,0,0,0.2+0.8*dnorm(.yy,Yn$mean,Yn$sd)*sqrt(2*pi)*Yn$sd),pch=20)
## Error in rgb(1, 0, 0, 0.2 + 0.8 * dnorm(.yy, Yn$mean, Yn$sd) * sqrt(2 * : alpha level 1, not in [0,1]
points(x=xn,y=.yy[8],col=rgb(0,1,0),pch=1)
points(x=xn,y=.yy[4],col=rgb(0,1,0),pch=1)

plot(.xx,EI(.xx,kmi),col='blue',type='l',xlab="x",ylab="ECI(xn,x)",xlim=c(xn-0.2,xn+0.2),ylim=c(0,0.05))
text(labels="EI(x)",x=.4,y=0.03,col='blue')

for (yn in .yy) {
  kmii <- km(design=rbind(X,xn),response=rbind(y,yn),control = list(trace=FALSE))
  lines(.xx,5*EI(.xx,kmii),col=rgb(1,0,0,0.2+0.8*dnorm(yn,Yn$mean,Yn$sd)*sqrt(2*pi)*Yn$sd),lty=2)
}
## Error in rgb(1, 0, 0, 0.2 + 0.8 * dnorm(yn, Yn$mean, Yn$sd) * sqrt(2 * : alpha level 1, not in [0,1]
yn = .yy[8]
kmii <- km(design=rbind(X,xn),response=rbind(y,yn),control = list(trace=FALSE))
lines(.xx,5*EI(.xx,kmii),col=rgb(0,1,0),lty=2)
yn = .yy[4]
kmii <- km(design=rbind(X,xn),response=rbind(y,yn),control = list(trace=FALSE))
lines(.xx,5*EI(.xx,kmii),col=rgb(0,1,0),lty=6)

source("https://raw.githubusercontent.com/IRSN/IECI/master/IECI_mc.R",echo = FALSE)
## Loading required package: rngWELL
## This is randtoolbox. For overview, type 'help("randtoolbox")'.
lines(.xx,Vectorize(function(.x)5*ECI_mc_vec.o2(.x,xn,kmi))(.xx),col='red',type='l')
## Error in predict_update_km(newXvalue = y0, newXmean = my0, newXvar = sdy0^2, : could not find function "predict_update_km"
text(labels="ECI(xn,x) = E[ EI(x) | yn~Y(xn)] (x50)",x=.3,y=0.045,col='red')

plot of chunk unnamed-chunk-4

Expected Conditional Improvement (ECI) by exact calculation (! fast !)

par(mfrow=c(2,1))

sectionview.km(kmi,col_surf='blue',col_points='blue',xlim=c(xn-0.2,xn+0.2),ylim=c(-1.1,-.7),title = "",Xname = "x")
plot(fun,type='l',add=T)
abline(h=min(y),col='blue',lty=2)
text(labels="m",x=0.35,y=min(y)-0.05,col='blue')


abline(v=xn,col='red')
text(labels="xn",x=xn+0.01,y=-0.75,col='red')

Yn=predict.km(kmi,newdata=xn,type="UK",checkNames = F)
.yy = seq(f=Yn$lower,t=Yn$upper,l=11)
#lines(x=dnorm(.xxx,yn$mean,yn$sd)/100+xn,y=.xxx,col='red')

plot(.xx,EI(.xx,kmi),col='blue',type='l',xlab="x",ylab="ECI(xn,x)",xlim=c(xn-0.2,xn+0.2),ylim=c(0,0.05))
text(labels="EI(x)",x=.4,y=0.03,col='blue')

source("https://raw.githubusercontent.com/IRSN/IECI/master/IECI.R",echo = FALSE)
## Warning in file(filename, "r", encoding = encoding): cannot open file
## 'int_Phi.R': No such file or directory
## Error in file(filename, "r", encoding = encoding): cannot open the connection
lines(.xx,5*ECI(matrix(.xx,ncol=1),xn,kmi),col='red',type='l')
## Error in ECI(matrix(.xx, ncol = 1), xn, kmi): could not find function "ECI"
text(labels="ECI(xn,x) = E[ EI(x) | yn~Y(xn)] (x50)",x=.3,y=0.045,col='red')

plot of chunk unnamed-chunk-5

And finally, Integrated Expected Conditional Improvement (IECI) calculated at the point xn:

> IECI_mc(xn,kmi)
[1] 0.002014966
> IECI(xn,kmi)
[1] 0.002115793

Overall criterions values (EI and IECI):

par(mfrow=c(1,1))

sectionview.km(kmi,ylim=c(-1,0))
abline(v=X$X,lty=2)

n=function(x){(x-min(x))/(max(x)-min(x))-1}

xx=seq(f=0,t=1,l=1000)
lines(xx,(fun(xx)),type='l')

lines(xx,n(IECI(x0 = xx,model = kmi,lower=0,upper=1)),type='l',col='blue')
## Error in IECI(x0 = xx, model = kmi, lower = 0, upper = 1): could not find function "IECI"
text(labels="IECI(x)",x=0.6,y=0.9-1,col='blue')

lines(xx,n(EI(xx,model=kmi)),col='red')
text(labels="EI(x)",x=0.6,y=0.1-1,col='red')

plot of chunk unnamed-chunk-6

Analytical development of IECI:

Computable form of ECI

ECI is viewed in two main parts : updated realizations where the new conditional point $y_n$ is below current minimum $m$, and updated realizations where the new conditional point $y_n$ is over current minimum $m$ (which then does not change this current minimum for EI calculation):

  • $y_n>m$ part:

Change variable: $z := \frac {y_n-\mu(x_n)}{\sigma(x_n)}$

  • $y_n=\sigma(x_n) z + \mu(x_n)$
  • $dz = dy_n/\sigma(x_n)$
  • $y_n \in [m,+\infty[ \sim z \in [m’ = (m-\mu(x_n))/\sigma(x_n), + \infty[$

So,

Considering that ([@KrigingUpdate]),

Finally,

Given,

  • $y_n<m$ part:

Change variable: $z := \frac {y_n-\mu(x_n)}{\sigma(x_n)}$

  • $y_n=\sigma(x_n) z + \mu(x_n)$
  • $dz = dy_n/\sigma(x_n)$
  • $y_n \in [m,+\infty[ \sim z \in [m’ = (m-\mu(x_n))/\sigma(x_n), + \infty[$

So,

Finally,

Given,

These integrals need the following expressions to become tractable:

  • $\int_{x}^{y} (a+b z) \Phi(a+b z) e^{-\frac{z^2}{2}} dz$
  • $\int_{x}^{y} \phi(a+b z) e^{-\frac{z^2}{2}} dz$

We have (see below):

And

It should be noticed that this exact calculation of ECI is vectorizable, which means that synchronized computations may be performed in an efficient manner.

Detailed calculation of integrals in $ECI^y_n$ and $ECI^m_n$

$\int_{x}^{y} \phi(a+b z) e^{-\frac{z^2}{2}} dz$

$\int_{x}^{y} (a+b z) \Phi(a+b z) e^{-\frac{z^2}{2}} dz$

Changing variable: $t:=a+b z$,

  • Then,

  • And,

We use the bi-normal cumulative density approximation [@Genz1992]:

Changing variable $t:=v-(a+b u)$,

and thus,