とりあえず最急降下法

うちの学類は数理計画とかはあるのに、学類で数値計算っぽいのがあんまりないので、基本と思われる最急降下法というのをやってみた。

• http://www.neurosci.aist.go.jp/~kurita/lecture/neuro/node2.html

もちろんRでやるよ!!

問題1

微分したやつ計算してもいいんだけど、後のことも考えてRで微分とかやっておく。こんな感じでできる。

trig.exp <- expression((a-1)^2)
D.sc <- D(trig.exp, "a")

で、evalして、funtionを適当にしてあげると微分された後の関数が返ってくる。

function(a){eval(D.sc)}

で、最急降下法とか言うので最適化問題を解いてやる。とりあえずなんでこれで出たりするのかは置いておく。100回回したら収束したけど、出力ははしょっておく。

> a <- -10
> #初期値
> alpha <- 0.1
> #なんじゃこりゃ
> 
> for(i in seq(100)){
+   a <- a - alpha * (function(a){eval(D.sc)})(a)
+   cat(paste(i,"is",a),fill=T)  
+ }
1 is -7.8
2 is -6.04
3 is -4.632
4 is -3.5056
5 is -2.60448
6 is -1.883584
7 is -1.3068672
8 is -0.84549376
9 is -0.476395008
10 is -0.181116006400000
11 is 0.0551071948800002
12 is 0.24408575590400
13 is 0.3952686047232
14 is 0.51621488377856
15 is 0.612971907022848
16 is 0.690377525618278
17 is 0.752302020494623
18 is 0.801841616395698
19 is 0.841473293116559
20 is 0.873178634493247

応用問題1

今度は関数が変わる。

f <- function(a){(a-1)^2*(a+1)^2}
plot(f,-2,2)
#プラマイ1がもちろんよさげ

で、今度も同様にやる。今回はパラメータの値と初期値に大分依存していることが実験から分かってきた。うまくいった場合のみ掲載。

trig.exp <- expression((a-1)^2*(a+1)^2)
D.sc <- D(trig.exp, "a")
plot(function(a){eval(D.sc)},-2,2)

a <- 0.75
alpha <- 0.1

for(i in seq(100)){
  a <- a - alpha * (function(a){eval(D.sc)})(a)
  cat(paste(i,"is",a),fill=T)  
}

a <- -0.75
alpha <- 0.1

for(i in seq(100)){
  a <- a - alpha * (function(a){eval(D.sc)})(a)
  cat(paste(i,"is",a),fill=T)  
}