勾配降下法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")
plot(cntlist,errorlist,type="l")
はい。収束しました。との違いが半端ないですね_φ(ﾟ▽ﾟ*)♪
次は、面白いです。
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")
plot(cntlist,errorlist,type="l")
収束しないですね。
