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

INFINITY -数学とかプログラミングとか-

統計とプログラムを使って役に立たせたい

TeX用コマンド入力を支援するための辞書をご利用ください。
sanctuary's blogは,適当なことが書いてあります。

勾配降下法2(最急降下法)

記事が長くなったので、新しい記事に書きます。

J(t)=t^2とします。

func1=function(t){
return(t^2)
}

dffunc1=function(t){
return(2*t)
}



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

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

}

> grad1(0.1,0.1,1e+5)
cnt:1;x=0.0800000000,abs(x-xbefore)=0.0200000000000
cnt:2;x=0.0640000000,abs(x-xbefore)=0.0160000000000
cnt:3;x=0.0512000000,abs(x-xbefore)=0.0128000000000
cnt:4;x=0.0409600000,abs(x-xbefore)=0.0102400000000
cnt:5;x=0.0327680000,abs(x-xbefore)=0.0081920000000
cnt:6;x=0.0262144000,abs(x-xbefore)=0.0065536000000
cnt:7;x=0.0209715200,abs(x-xbefore)=0.0052428800000
cnt:8;x=0.0167772160,abs(x-xbefore)=0.0041943040000
cnt:9;x=0.0134217728,abs(x-xbefore)=0.0033554432000
cnt:10;x=0.0107374182,abs(x-xbefore)=0.0026843545600
cnt:11;x=0.0085899346,abs(x-xbefore)=0.0021474836480
cnt:12;x=0.0068719477,abs(x-xbefore)=0.0017179869184
cnt:13;x=0.0054975581,abs(x-xbefore)=0.0013743895347
cnt:14;x=0.0043980465,abs(x-xbefore)=0.0010995116278
cnt:15;x=0.0035184372,abs(x-xbefore)=0.0008796093022
cnt:16;x=0.0028147498,abs(x-xbefore)=0.0007036874418
cnt:17;x=0.0022517998,abs(x-xbefore)=0.0005629499534
cnt:18;x=0.0018014399,abs(x-xbefore)=0.0004503599627
cnt:19;x=0.0014411519,abs(x-xbefore)=0.0003602879702
cnt:20;x=0.0011529215,abs(x-xbefore)=0.0002882303762
cnt:21;x=0.0009223372,abs(x-xbefore)=0.0002305843009
cnt:22;x=0.0007378698,abs(x-xbefore)=0.0001844674407
cnt:23;x=0.0005902958,abs(x-xbefore)=0.0001475739526
cnt:24;x=0.0004722366,abs(x-xbefore)=0.0001180591621
cnt:25;x=0.0003777893,abs(x-xbefore)=0.0000944473297
cnt:26;x=0.0003022315,abs(x-xbefore)=0.0000755578637
cnt:27;x=0.0002417852,abs(x-xbefore)=0.0000604462910
cnt:28;x=0.0001934281,abs(x-xbefore)=0.0000483570328
cnt:29;x=0.0001547425,abs(x-xbefore)=0.0000386856262
cnt:30;x=0.0001237940,abs(x-xbefore)=0.0000309485010
cnt:31;x=0.0000990352,abs(x-xbefore)=0.0000247588008
cnt:32;x=0.0000792282,abs(x-xbefore)=0.0000198070406
cnt:33;x=0.0000633825,abs(x-xbefore)=0.0000158456325
cnt:34;x=0.0000507060,abs(x-xbefore)=0.0000126765060
cnt:35;x=0.0000405648,abs(x-xbefore)=0.0000101412048
cnt:36;x=0.0000324519,abs(x-xbefore)=0.0000081129638
cnt:37;x=0.0000259615,abs(x-xbefore)=0.0000064903711
cnt:38;x=0.0000207692,abs(x-xbefore)=0.0000051922969
cnt:39;x=0.0000166153,abs(x-xbefore)=0.0000041538375
cnt:40;x=0.0000132923,abs(x-xbefore)=0.0000033230700
cnt:41;x=0.0000106338,abs(x-xbefore)=0.0000026584560
cnt:42;x=0.0000085071,abs(x-xbefore)=0.0000021267648
cnt:43;x=0.0000068056,abs(x-xbefore)=0.0000017014118
cnt:44;x=0.0000054445,abs(x-xbefore)=0.0000013611295
cnt:45;x=0.0000043556,abs(x-xbefore)=0.0000010889036
cnt:46;x=0.0000034845,abs(x-xbefore)=0.0000008711229
cnt:47;x=0.0000027876,abs(x-xbefore)=0.0000006968983
cnt:48;x=0.0000022301,abs(x-xbefore)=0.0000005575186
cnt:49;x=0.0000017841,abs(x-xbefore)=0.0000004460149
cnt:50;x=0.0000014272,abs(x-xbefore)=0.0000003568119
cnt:51;x=0.0000011418,abs(x-xbefore)=0.0000002854495
cnt:52;x=0.0000009134,abs(x-xbefore)=0.0000002283596
cnt:53;x=0.0000007308,abs(x-xbefore)=0.0000001826877
cnt:54;x=0.0000005846,abs(x-xbefore)=0.0000001461502
cnt:55;x=0.0000004677,abs(x-xbefore)=0.0000001169201
cnt:56;x=0.0000003741,abs(x-xbefore)=0.0000000935361
cnt:57;x=0.0000002993,abs(x-xbefore)=0.0000000748289
cnt:58;x=0.0000002395,abs(x-xbefore)=0.0000000598631
cnt:59;x=0.0000001916,abs(x-xbefore)=0.0000000478905
cnt:60;x=0.0000001532,abs(x-xbefore)=0.0000000383124
cnt:61;x=0.0000001226,abs(x-xbefore)=0.0000000306499
cnt:62;x=0.0000000981,abs(x-xbefore)=0.0000000245199
cnt:63;x=0.0000000785,abs(x-xbefore)=0.0000000196159
cnt:64;x=0.0000000628,abs(x-xbefore)=0.0000000156928
cnt:65;x=0.0000000502,abs(x-xbefore)=0.0000000125542
cnt:66;x=0.0000000402,abs(x-xbefore)=0.0000000100434
cnt:67;x=0.0000000321,abs(x-xbefore)=0.0000000080347
cnt:68;x=0.0000000257,abs(x-xbefore)=0.0000000064278
cnt:69;x=0.0000000206,abs(x-xbefore)=0.0000000051422
cnt:70;x=0.0000000165,abs(x-xbefore)=0.0000000041138
cnt:71;x=0.0000000132,abs(x-xbefore)=0.0000000032910
cnt:72;x=0.0000000105,abs(x-xbefore)=0.0000000026328
cnt:73;x=0.0000000084,abs(x-xbefore)=0.0000000021062
cnt:74;x=0.0000000067,abs(x-xbefore)=0.0000000016850
cnt:75;x=0.0000000054,abs(x-xbefore)=0.0000000013480
cnt:76;x=0.0000000043,abs(x-xbefore)=0.0000000010784
cnt:77;x=0.0000000035,abs(x-xbefore)=0.0000000008627
cnt:78;x=0.0000000028,abs(x-xbefore)=0.0000000006902
cnt:79;x=0.0000000022,abs(x-xbefore)=0.0000000005521
cnt:80;x=0.0000000018,abs(x-xbefore)=0.0000000004417
cnt:81;x=0.0000000014,abs(x-xbefore)=0.0000000003534
cnt:82;x=0.0000000011,abs(x-xbefore)=0.0000000002827
cnt:83;x=0.0000000009,abs(x-xbefore)=0.0000000002262
cnt:84;x=0.0000000007,abs(x-xbefore)=0.0000000001809
cnt:85;x=0.0000000006,abs(x-xbefore)=0.0000000001447
cnt:86;x=0.0000000005,abs(x-xbefore)=0.0000000001158
cnt:87;x=0.0000000004,abs(x-xbefore)=0.0000000000926
cnt:88;x=0.0000000003,abs(x-xbefore)=0.0000000000741
cnt:89;x=0.0000000002,abs(x-xbefore)=0.0000000000593
cnt:90;x=0.0000000002,abs(x-xbefore)=0.0000000000474
cnt:91;x=0.0000000002,abs(x-xbefore)=0.0000000000379
cnt:92;x=0.0000000001,abs(x-xbefore)=0.0000000000304
cnt:93;x=0.0000000001,abs(x-xbefore)=0.0000000000243
cnt:94;x=0.0000000001,abs(x-xbefore)=0.0000000000194
cnt:95;x=0.0000000001,abs(x-xbefore)=0.0000000000155
cnt:96;x=0.0000000000,abs(x-xbefore)=0.0000000000124
cnt:97;x=0.0000000000,abs(x-xbefore)=0.0000000000099


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

f:id:infinity_th4:20120615004357j:image:w360

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

f:id:infinity_th4:20120615004356j:image:w360

はい。収束しました。t^4との違いが半端ないですね_φ(゚▽゚*)♪


次は、面白いです。
http://ja.wikipedia.org/wiki/%E6%9C%80%E6%80%A5%E9%99%8D%E4%B8%8B%E6%B3%95の問題点の箇所にalpha>1のときに反復ごとに悪い解へと向かう。とあったので、じゃあalpha=1のときどうやねんということでやってみました。

> grad1(1,1,1e+2)
cnt:1;x=-1.0000000000,abs(x-xbefore)=2.0000000000000
cnt:2;x=1.0000000000,abs(x-xbefore)=2.0000000000000
cnt:3;x=-1.0000000000,abs(x-xbefore)=2.0000000000000
cnt:4;x=1.0000000000,abs(x-xbefore)=2.0000000000000
cnt:5;x=-1.0000000000,abs(x-xbefore)=2.0000000000000
cnt:6;x=1.0000000000,abs(x-xbefore)=2.0000000000000
cnt:7;x=-1.0000000000,abs(x-xbefore)=2.0000000000000
cnt:8;x=1.0000000000,abs(x-xbefore)=2.0000000000000
cnt:9;x=-1.0000000000,abs(x-xbefore)=2.0000000000000
cnt:10;x=1.0000000000,abs(x-xbefore)=2.0000000000000

略

cnt:100;x=1.0000000000,abs(x-xbefore)=2.0000000000000

となって、times=1e+5であろうとも同じです。


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

f:id:infinity_th4:20120615005436j:image:w360


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

f:id:infinity_th4:20120615005435j:image:w360


収束しないですね。