最適化問題
#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")
3つの手法の比較
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)