最適化問題

 
#Objective Function 
funZ<-function(x,y) {x^4-16*x^2-5*x+y^4-16*y^2-5*y}
xmin=-5;xmax=5;n_len=100;
xs=ys<-seq(xmin, xmax,length.out=n_len)
z<-outer(xs,ys,funZ)
contour(x=xs,y=ys,z=z,nlevels=30,drawlabels=F,col="darkgreen",
 main=expression(z=x^4-16*x^2-5*x+y^4-16*y^2-5*y))

勾配法の実装例

 
# initialization
x=x.hist=runif(1,min=-5,max=5) # initial value of x
y=y.hist=runif(1,min=-5,max=5) # initial value of y
lambda=0.01                    # step size
zeta=1e-8                      # tol (stopping criterion)
t=0                            # time
g=100                          # maximum gradient
# end initialization

### note ###
# obj fun: z=x^4-16*x^2-5*x+y^4-16*y^2-5*y
# part. deriv. x: 4*x^3-32*x-5
# part. deriv. y: 4*y^3-32*y-5
#############

# main 
while (g>zeta) {
  t=t+1
  g.x=4*x^3-32*x-5
  g.y=4*y^3-32*y-5
  x=x-lambda*g.x
  y=y-lambda*g.y
  g=max(abs(g.x),abs(g.y))
  x.hist=c(x.hist,x)
  y.hist=c(y.hist,y)
}

# plotting results
par(mfrow=c(1,3))
par(mar=c(4.0,4.1,4.1,1.1))
contour(x=xs,y=ys,z=z,nlevels=30,drawlabels=F,col="darkgreen",
 main=expression(z=x^4-16*x^2-5*x+y^4-16*y^2-5*y))
lines(x.hist,y.hist,type='o',pch=20,col='red')
plot(x.hist,type='l',xlab="time",ylab="value of x",lwd=2,
 main="how value of y changes over time")
plot(y.hist,type='l',xlab="time",ylab="value of y",lwd=2,
 main="how value of y changes over time")

単純確率的最適化法の実装例

 
# initialization
x=runif(1,min=-5,max=5)             # initial value of x
y=runif(1,min=-5,max=5)             # initial value of y
f.xy=x^4-16*x^2-5*x+y^4-16*y^2-5*y  # initial value of obj. fun.
w=0.1                               # search width
t=0                                 # time
max.t=1000                          # max iter (stopping criterion)
stop.crit=F                         # stop criterion
x.hist=rep(0,max.t)                 # history of x
y.hist=rep(0,max.t)                 # history of x
x.hist[1]=x
y.hist[1]=y
# end initialization

# main
while (stop.crit==F) {
  t=t+1
  x.temp=x+rnorm(1,mean=0,sd=w)
  y.temp=y+rnorm(1,mean=0,sd=w)
  f.temp=funZ(x.temp,y.temp)
  if (f.temp < f.xy) {
    x=x.temp
    y=y.temp
    f.xy=f.temp
  }
  x.hist[t]=x
  y.hist[t]=y
  if (t==max.t){stop.crit=T}
}

# plotting results
par(mfrow=c(1,3))
par(mar=c(4.0,4.1,4.1,1.1))
contour(x=xs,y=ys,z=z,nlevels=30,drawlabels=F,col="darkgreen",
 main=expression(z=x^4-16*x^2-5*x+y^4-16*y^2-5*y))
lines(x.hist,y.hist,type='o',pch=20,col='red')
plot(x.hist,type='l',xlab="time",ylab="value of x",lwd=2,
 main="how value of y changes over time")
plot(y.hist,type='l',xlab="time",ylab="value of y",lwd=2,
 main="how value of y changes over time")

焼きなまし法の実装例

# initialization
x=runif(1,min=-5,max=5)             # initial value of x
y=runif(1,min=-5,max=5)             # initial value of y
f.xy=x^4-16*x^2-5*x+y^4-16*y^2-5*y  # initial value of obj. fun.
t=0                                 # time
max.t=1000                          # max iter (stopping criterion)
stop.crit=F                         # stop criterion
x.hist=rep(0,max.t)                 # history of x
y.hist=rep(0,max.t)                 # history of x
x.hist[1]=x
y.hist[1]=y
tau=1
c=1
g=0.999
# end initialization
 
# main
while (stop.crit==F) {
  t=t+1
  x.temp=x+rnorm(1,mean=0,sd=tau)
  y.temp=y+rnorm(1,mean=0,sd=tau)
  f.temp=funZ(x.temp,y.temp)
  deltaE=f.temp-f.xy
  p=1/(1+exp(deltaE/(c*tau)))
  if (runif(1) < p) {
    x=x.temp
    y=y.temp
    f.xy=f.temp
  }
  x.hist[t]=x
  y.hist[t]=y
  tau=tau*g
  if (t==max.t){stop.crit=T}
}
 
# plotting results
par(mfrow=c(1,3))
par(mar=c(4.0,4.1,4.1,1.1))
contour(x=xs,y=ys,z=z,nlevels=30,drawlabels=F,col="darkgreen",
 main=expression(z=x^4-16*x^2-5*x+y^4-16*y^2-5*y))
lines(x.hist,y.hist,type='o',pch=20,col='red')
plot(x.hist,type='l',xlab="time",ylab="value of x",lwd=2,
 main="how value of y changes over time")
plot(y.hist,type='l',xlab="time",ylab="value of y",lwd=2,
 main="how value of y changes over time")

３つの手法の比較

 

GD<-function(x.init,y.init,n.iter,lambda){
x=x.init;y=y.init;z.hist=rep(0,n.iter)
# main 
for (i_iter in 2:n.iter) {
  g.x=4*x^3-32*x-5
  g.y=4*y^3-32*y-5
  x=x-lambda*g.x
  y=y-lambda*g.y
  g=max(abs(g.x),abs(g.y))
  x.hist=c(x.hist,x)
  y.hist=c(y.hist,y)
  z.hist[i_iter]=funZ(x,y)
}
return(z.hist)
}

stochOpt<-function(x.init,y.init,n.iter,w) {
x=x.init;y=y.init;f.xy=funZ(x,y)
z.hist=rep(0,n.iter);z.hist[1]=f.xy
# main
for (i_iter in 2:n.iter) {
  x.temp=x+rnorm(1,mean=0,sd=w)
  y.temp=y+rnorm(1,mean=0,sd=w)
  f.temp=funZ(x.temp,y.temp)
  if (f.temp < f.xy) {
    x=x.temp;y=y.temp;f.xy=f.temp
  }
  z.hist[i_iter]=f.xy
}
return(z.hist)
}
 
simAnn<-function(x.init,y.init,n.iter,tau,c,g) {
x=x.init;y=y.init;f.xy=funZ(x,y)
z.hist=rep(0,n.iter);z.hist[1]=f.xy
# main
for (i_iter in 2:n.iter) {
  x.temp=x+rnorm(1,mean=0,sd=tau)
  y.temp=y+rnorm(1,mean=0,sd=tau)
  f.temp=funZ(x.temp,y.temp)
  deltaE=f.temp-f.xy
  p=1/(1+exp(deltaE/(c*tau)))
  if (runif(1) < p) {
    x=x.temp;y=y.temp;f.xy=f.temp
  }
  z.hist[i_iter]=f.xy
  tau=tau*g
}
return(z.hist)
}

n.rep=500
n.iter=500
xs=runif(n.rep,-3,3);ys=runif(n.rep,-3,3)
GD.hist=SA.hist=SO.hist=matrix(0,n.iter,n.rep)
for (i.rep in 1:n.rep) {
  GD.hist[,i.rep]=GD(xs[i.rep],ys[i.rep],n.iter,0.025)
  SO.hist[,i.rep]=stochOpt(xs[i.rep],ys[i.rep],n.iter,2)
  SA.hist[,i.rep]=simAnn(xs[i.rep],ys[i.rep],n.iter,3,1,0.999)
}

plot(rowMeans(GD.hist),type='l',ylim=c(-160,0),lwd=3)
lines(rowMeans(SO.hist),col='red',lwd=3)
lines(rowMeans(SA.hist),col='blue',lwd=3)