勾配降下法3(最急降下法) - INFINITY -数学とかプログラミングとか-

2つ解がある場合。

x=seq(-10,10,,10000)

func4=function(t){
return((t+3)^2*(t-3)^2)
}

plot(x,func4(x),xlim=c(-10,10),ylim=c(0,100),type="l")


x=seq(-10,10,,10000)

dffunc4=function(t){
return(2*(t+3)*(t-3)^2+(t+3)^2*(t-3))
}

plot(x,dffunc4(x),xlim=c(-10,10),ylim=c(-100,100),type="l")



xlist=c()
cntlist=c()
errorlist=c()

grad4=function(x,alpha,times){
	cnt=1
	for(i in 1:times){
		xbefore=x
		x=x-alpha*dffunc4(x)
		#if(cnt==1||0==cnt%%100){
			xlist<<-c(xlist,x)
			cntlist<<-c(cntlist,cnt)
			errorlist<<-c(errorlist,abs(x-xbefore))
			printf("cnt:%d;x=%f,abs(x-xbefore)=%.13f\n",cnt,x,abs(x-xbefore))
		#}
		if(abs(x-xbefore)<=1.0e-8){
			#printf("cnt:%d;x=%f,abs(x-xbefore)=%.13f\n",cnt,x,abs(x-xbefore))		
			break
		}
		cnt=cnt+1
	}

}

初期値が離れすぎてる場合は、発散しちゃう。

> grad4(10,0.1,100000)
cnt:1;x=-235.700000,abs(x-xbefore)=245.7000000000000
cnt:2;x=3944049.544900,abs(x-xbefore)=3944285.2449000012130
cnt:3;x=-18405525867132454912.000000,abs(x-xbefore)=18405525867136399360.0000000000000
cnt:4;x=1870535459532385043044026668888600002026004600862808044044.000000,abs(x-xbefore)=1870535459532385043044026668888600002026004600862808044044.0000000000000
cnt:5;x=-1963446586185536064848800206820440606442864844482246402440408468802820884668660608868602668468244402846888446686228086846842688044468004844826422608262846082866864022468868.000000,abs(x-xbefore)=1963446586185536064848800206820440606442864844482246402440408468802820884668660608868602668468244402846888446686228086846842688044468004844826422608262846082866864022468868.0000000000000
cnt:6;x=Inf,abs(x-xbefore)=Inf
cnt:7;x=NaN,abs(x-xbefore)=NaN
 以下にエラー if (abs(x - xbefore) <= 1e-08) { : 
   TRUE/FALSE が必要なところが欠損値です

同じく、

> grad4(10,0.01,100000)
cnt:1;x=-14.570000,abs(x-xbefore)=24.5700000000000
cnt:2;x=80.384377,abs(x-xbefore)=94.9543767900000
cnt:3;x=-15286.798903,abs(x-xbefore)=15367.1832797677998
cnt:4;x=107176419039.288223,abs(x-xbefore)=107176434326.0871276855469
cnt:5;x=-36933373827666237250648062244666.000000,abs(x-xbefore)=36933373827666237250648062244666.0000000000000
cnt:6;x=1511395762460033041124802820688646600400622204806208242208224484804484400044080820064224040480.000000,abs(x-xbefore)=1511395762460033041124802820688646600400622204806208242208224484804484400044080820064224040480.0000000000000
cnt:7;x=-103575217854207314670264642682424206026248604088200064668042206024644288482422064608086202884268088042062820886644084404680002868228248246080448840246626888488406880408008820826882882048466820888242064682688648284682800486240420262446686486622684826686486624420624246444486868428.000000,abs(x-xbefore)=103575217854207314670264642682424206026248604088200064668042206024644288482422064608086202884268088042062820886644084404680002868228248246080448840246626888488406880408008820826882882048466820888242064682688648284682800486240420262446686486622684826686486624420624246444486868428.0000000000000
cnt:8;x=Inf,abs(x-xbefore)=Inf
cnt:9;x=NaN,abs(x-xbefore)=NaN
 以下にエラー if (abs(x - xbefore) <= 1e-08) { : 
   TRUE/FALSE が必要なところが欠損値です

alphaを小さくすると良いようです。

> grad4(10,0.001,100000)
cnt:1;x=7.543000,abs(x-xbefore)=2.4570000000000
cnt:2;x=6.602833,abs(x-xbefore)=0.9401672490210
cnt:3;x=6.021302,abs(x-xbefore)=0.5815303431274
cnt:4;x=5.610719,abs(x-xbefore)=0.4105831008715
cnt:5;x=5.299770,abs(x-xbefore)=0.3109492773220
cnt:6;x=5.053554,abs(x-xbefore)=0.2462163858180
cnt:7;x=4.852436,abs(x-xbefore)=0.2011179286358
cnt:8;x=4.684322,abs(x-xbefore)=0.1681141195830
cnt:9;x=4.541265,abs(x-xbefore)=0.1430570732112
cnt:10;x=4.417783,abs(x-xbefore)=0.1234812393203

略

cnt:100;x=3.025785,abs(x-xbefore)=0.0009808316557

略

cnt:200;x=3.000648,abs(x-xbefore)=0.0000242106286

略

cnt:301;x=3.000016,abs(x-xbefore)=0.0000005962107

略

cnt:400;x=3.000000,abs(x-xbefore)=0.0000000158131
cnt:401;x=3.000000,abs(x-xbefore)=0.0000000152438
cnt:402;x=3.000000,abs(x-xbefore)=0.0000000146950
cnt:403;x=3.000000,abs(x-xbefore)=0.0000000141660
cnt:404;x=3.000000,abs(x-xbefore)=0.0000000136560
cnt:405;x=3.000000,abs(x-xbefore)=0.0000000131644
cnt:406;x=3.000000,abs(x-xbefore)=0.0000000126905
cnt:407;x=3.000000,abs(x-xbefore)=0.0000000122336
cnt:408;x=3.000000,abs(x-xbefore)=0.0000000117932
cnt:409;x=3.000000,abs(x-xbefore)=0.0000000113687
cnt:410;x=3.000000,abs(x-xbefore)=0.0000000109594
cnt:411;x=3.000000,abs(x-xbefore)=0.0000000105649
cnt:412;x=3.000000,abs(x-xbefore)=0.0000000101845
cnt:413;x=3.000000,abs(x-xbefore)=0.0000000098179

初期値が離れすぎていると、alphaを小さくしなければならないようですね。

plot(cntlist,xlist,type="l")

plot(cntlist,errorlist,type="l")

予想外の結果：

0からすこしでも3に近ければつまり、初期値が0.00001なら、3に収束すると思っていましたが、違いましたね。

> grad4(0.00001,0.01,100000)
cnt:1;x=-0.269987,abs(x-xbefore)=0.2699972999970
cnt:2;x=-0.610107,abs(x-xbefore)=0.3401193700552
cnt:3;x=-1.026856,abs(x-xbefore)=0.4167488935507
cnt:4;x=-1.509991,abs(x-xbefore)=0.4831355400126
cnt:5;x=-2.016000,abs(x-xbefore)=0.5060086995937
cnt:6;x=-2.462586,abs(x-xbefore)=0.4465861597730
cnt:7;x=-2.767536,abs(x-xbefore)=0.3049503309645
cnt:8;x=-2.919075,abs(x-xbefore)=0.1515389217677
cnt:9;x=-2.975392,abs(x-xbefore)=0.0563170973292
cnt:10;x=-2.992929,abs(x-xbefore)=0.0175363201582
cnt:11;x=-2.998005,abs(x-xbefore)=0.0050763934341
cnt:12;x=-2.999440,abs(x-xbefore)=0.0014351870847
cnt:13;x=-2.999843,abs(x-xbefore)=0.0004029521182
cnt:14;x=-2.999956,abs(x-xbefore)=0.0001129132171
cnt:15;x=-2.999988,abs(x-xbefore)=0.0000316225010
cnt:16;x=-2.999997,abs(x-xbefore)=0.0000088548336
cnt:17;x=-2.999999,abs(x-xbefore)=0.0000024793952
cnt:18;x=-3.000000,abs(x-xbefore)=0.0000006942339
cnt:19;x=-3.000000,abs(x-xbefore)=0.0000001943858
cnt:20;x=-3.000000,abs(x-xbefore)=0.0000000544280
cnt:21;x=-3.000000,abs(x-xbefore)=0.0000000152399
cnt:22;x=-3.000000,abs(x-xbefore)=0.0000000042672

grad4(0.9,0.001,100000)#-3に収束

grad4(2,0.001,100000)#3に収束

前者は-3に収束。後者は3に収束という結果でした。

なんででしょうかね。

それに、

> grad4(1,0.01,100000)
cnt:1;x=1.000000,abs(x-xbefore)=0.0000000000000

も、なぜだか分からない。

式見ればわかるとか？

今日は、この辺で筆(？)を置きましょう壁|)≡ｻｯ!!

では。

おしまい。

すべておしまい。

次は、多次元版をやる必要がありますね。

追記(6/23)：

> grad4(0.1,0.01,100000)
cnt:1;x=-0.142730,abs(x-xbefore)=0.2427300000000
cnt:2;x=-0.450569,abs(x-xbefore)=0.3078387141761
cnt:3;x=-0.833388,abs(x-xbefore)=0.3828190599212
cnt:4;x=-1.290202,abs(x-xbefore)=0.4568141296647
cnt:5;x=-1.794187,abs(x-xbefore)=0.5039849717840
cnt:6;x=-2.278774,abs(x-xbefore)=0.4845868986786
cnt:7;x=-2.653261,abs(x-xbefore)=0.3744874502147
cnt:8;x=-2.868095,abs(x-xbefore)=0.2148342138415
cnt:9;x=-2.957916,abs(x-xbefore)=0.0898204899029
cnt:10;x=-2.987687,abs(x-xbefore)=0.0297714469916
cnt:11;x=-2.996507,abs(x-xbefore)=0.0088196656529
cnt:12;x=-2.999018,abs(x-xbefore)=0.0025112716531
cnt:13;x=-2.999725,abs(x-xbefore)=0.0007065259288
cnt:14;x=-2.999923,abs(x-xbefore)=0.0001980936312
cnt:15;x=-2.999978,abs(x-xbefore)=0.0000554871485
cnt:16;x=-2.999994,abs(x-xbefore)=0.0000155380437
cnt:17;x=-2.999998,abs(x-xbefore)=0.0000043507810
cnt:18;x=-3.000000,abs(x-xbefore)=0.0000012182288
cnt:19;x=-3.000000,abs(x-xbefore)=0.0000003411048
cnt:20;x=-3.000000,abs(x-xbefore)=0.0000000955094
cnt:21;x=-3.000000,abs(x-xbefore)=0.0000000267426
cnt:22;x=-3.000000,abs(x-xbefore)=0.0000000074879

対称的にならないといけないのではないのか？

> grad4(-0.1,0.01,100000)
cnt:1;x=-0.396670,abs(x-xbefore)=0.2966700000000
cnt:2;x=-0.767178,abs(x-xbefore)=0.3705080412404
cnt:3;x=-1.213113,abs(x-xbefore)=0.4459352480744
cnt:4;x=-1.712946,abs(x-xbefore)=0.4998331510877
cnt:5;x=-2.206633,abs(x-xbefore)=0.4936868950518
cnt:6;x=-2.604009,abs(x-xbefore)=0.3973758783528
cnt:7;x=-2.843943,abs(x-xbefore)=0.2399335992785
cnt:8;x=-2.949112,abs(x-xbefore)=0.1051690386810
cnt:9;x=-2.984978,abs(x-xbefore)=0.0358665390428
cnt:10;x=-2.995726,abs(x-xbefore)=0.0107479656826
cnt:11;x=-2.998798,abs(x-xbefore)=0.0030715464587
cnt:12;x=-2.999663,abs(x-xbefore)=0.0008650764164
cnt:13;x=-2.999906,abs(x-xbefore)=0.0002426207819
cnt:14;x=-2.999974,abs(x-xbefore)=0.0000679652193
cnt:15;x=-2.999993,abs(x-xbefore)=0.0000190327251
cnt:16;x=-2.999998,abs(x-xbefore)=0.0000053293562
cnt:17;x=-2.999999,abs(x-xbefore)=0.0000014922349
cnt:18;x=-3.000000,abs(x-xbefore)=0.0000004178270
cnt:19;x=-3.000000,abs(x-xbefore)=0.0000001169916
cnt:20;x=-3.000000,abs(x-xbefore)=0.0000000327577
cnt:21;x=-3.000000,abs(x-xbefore)=0.0000000091721

>

？ $\alpha$ をどう選ぶと発散しなくて済むかという話ですが、ゼミで苦し紛れに

$| \alpha (t_{k-1}) \frac{\partial J(w)}{\partial w}|_{t_{k-1}}| \leq | \alpha (t_{k-1})\frac{\partial J(w)}{\partial w}|_{t_{k}}|$ となるような $\alpha(t_{k})$ を選ばなければ収束しないと答えましたが、これは検証する必要がありますね。