最適化法の比較

認知情報科学基礎実習 課題03
提出期限:2013.01.31 23:59

勾配法、単純確率的最適化法、焼き鈍し法を比較してみましょう。
目的関数:x^4-16*x^2-5*x+y^4-16*y^2-5*y
初期値:一様分布(-0.5~0.5)に従う乱数
繰り返し数:1000回
各最適化法の結果の平均を可視化して(下の図のようなもの)傾向を考察してみてください。
各最適化は下のものを参考にしてください。
各最適化法のパラメターを変えてみてましょう。

3optims

 #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;
x<-seq(xmin, xmax,length.out=n_len);y<-x;z<-outer(x,y,funZ)
 #gradient descent 
minGD<-function(initXY=c(0,0), maxItr=1000, stepSize=0.01){
  x=initXY[1];y=initXY[2];z=funZ(x,y); 
  opHist=matrix(0,nrow=maxItr,ncol=3)
  opHist[1,]=c(x,y,z)
  for (i_loop in 2:maxItr) {
    x=x-stepSize*(4*x^3-32*x-5);
    y=y-stepSize*(4*y^3-32*x-5);
    z=funZ(x,y);
    opHist[counter,]=c(x,y,z);
  }
  return(opHist[1:counter,])
}

res<-minGD(c(0.5,0.5),1000,0.01)
contour(x,y,z,nlevels=30,drawlabel=F)
lines(res[,1:2],col='cyan',type='o',pch=19)
#Simple Stochastic Optimization
minSTOCH<-function(initXY=c(0,0),maxItr=1000,width=0.1) {
  x=initXY[1];y=initXY[2];z=funZ(x,y); 
  opHist=matrix(0,nrow=maxItr,ncol=3)
  opHist[1,]=c(x,y,z)
  for (i_loop in 2:maxItr) {
    xTemp=x+rnorm(1,mean=0,sd=width);
    yTemp=y+rnorm(1,mean=0,sd=width);
    zTemp=funZ(xTemp,yTemp);
    if (z > zTemp) {
  	x=xTemp;y=yTemp;z=zTemp;
    }
    opHist[i_loop,]=c(x,y,z);
  }
  return(opHist)
}

opHist<-minSTOCH(c(0,0),1000,0.1)
contour(x,y,z,nlevels=30,drawlabel=F)
lines(opHist[,1:2],type='o',lwd=2,col='green',pch=20)
#Simulated Annealing
minSA<-function(initXY=c(0,0),maxItr=1000,C=1,eta=0.99,width=10) {
  x=initXY[1];y=initXY[2];z=funZ(x,y); 
  opHist=matrix(0,nrow=maxItr,ncol=3)
  opHist[1,]=c(x,y,z)
  for (i_loop in 2:maxItr) {
    xTemp=x+rnorm(1,mean=0,sd=width)
    yTemp=y+rnorm(1,mean=0,sd=width)
    zTemp=funZ(xTemp, yTemp);
    delta=zTemp-z;
    prob=1/(1+exp(delta/(C*width)))
    if (runif(1) < prob) {
      x=xTemp;y=yTemp;z=zTemp;
     } 
    opHist[i_loop,]=c(x,y,z);
    width=width*eta
  }
  return(opHist)
}

res<-minSA(c(0,0),1000,1,0.99,10)
contour(x,y,z,nlevels=30,drawlabel=F)
lines(res[,1:2],type='o',lwd=2,col='green',pch=20)