先生が授業でmathematicaを使ってデモを見せてくれたりしたんだけど、「え、それRでできるよ?」と言いたいがためだけに、実装してみるなどしました。元の関数をplotして、最急降下法で、解が改善されて行く様子をplotする簡単なお仕事です。persp関数はあんまり好きになれないなぁ。outer関数を使わないといけないってのがかわいくない。
を無制約条件で最小化しています*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)
