読者です 読者をやめる 読者になる 読者になる

2変数関数での最急降下法をRでやってみた

最適化理論 R

先生が授業でmathematicaを使ってデモを見せてくれたりしたんだけど、「え、それRでできるよ?」と言いたいがためだけに、実装してみるなどしました。元の関数をplotして、最急降下法で、解が改善されて行く様子をplotする簡単なお仕事です。persp関数はあんまり好きになれないなぁ。outer関数を使わないといけないってのがかわいくない。
f:id:syou6162:20161009204452p:plain
100(x_1-x_2)^2+(x_1-1)^2を無制約条件で最小化しています*1。特徴とか、苦労したところは

  • 数値微分ではなく、数式微分している*2
  • 3次元の関数のplotとラインのプロットを同時にやるところに苦労した
    • 3次元で二次元のようなlinesやpointsを使うにはtrans3dを一回かましてから使わないといけない、というところにはまった

というところか。あと変数が増えてきたら、今のコードではだめなのでベクトル形式にしないと崩壊するんだけど、面倒だから放置(ry。まあ、簡単なところは自分でとりあえずやって、仕組みがある程度分かったら最適化とかしてあるライブラリを使えばいいんじゃないかな。

library(scatterplot3d)

f <- function(x){x1 <- x[1];x2 <- x[2];return(100*(x1-x2)^2+(x1-1)^2)}
trig.exp <- expression(100*(x1-x2)^2+(x1-1)^2)
D.sc1 <- D(trig.exp, "x1")
D.sc2 <- D(trig.exp, "x2")
x_bar1 <- -10
x_bar2 <- -20
epsilon <- 0.3
tau <- 0.9
beta <- 1
n <- 2000
x_vec <- matrix(0,n+1,2,byrow=TRUE)
for(i in seq(n)){
  nabla_f <- c((function(x1,x2){eval(D.sc1)})(x_bar1,x_bar2),
               (function(x1,x2){eval(D.sc2)})(x_bar1,x_bar2))
  while(f(c(x_bar1,x_bar2) - beta * nabla_f) >
        f(c(x_bar1,x_bar2)) - epsilon * beta * sum(c((function(x1,x2){eval(D.sc1)})(x_bar1,x_bar2),(function(x1,x2){eval(D.sc2)})(x_bar1,x_bar2))^2)){
    beta <- tau * beta
  }
  a <- beta
  x_bar1_old <- x_bar1
  x_bar1 <- x_bar1 - a * (function(x1,x2){eval(D.sc1)})(x_bar1,x_bar2)
  x_bar2 <- x_bar2 - a * (function(x1,x2){eval(D.sc2)})(x_bar1_old,x_bar2)
  x_vec[i,] <- c(x_bar1,x_bar2)
  cat(paste(x_bar1),"\t")
  cat(paste(x_bar2),fill=T)
}

x <- seq(-1,1,length=50)
y <- x
z <- outer(x,y,function(x,y){mapply(function(x,y){f(c(x,y))},x,y)})
persp(x,y,z,theta=150,phi=10,expand=0.5) -> res
x <- x_vec[,1]
y <- x_vec[,2]
z <- mapply(function(x,y){f(c(x,y))},x,y)
lines(trans3d(x,y,z, pmat = res),col="red",lwd=5)

*1:授業でもこの関数でやっていたので

*2:Rがやってくれてます