最急降下法と共役勾配法
前期の授業で習ったまま放置してたのでちょっと復習。
最適化問題を解く際に利用される手法。
数式を書くには体力が足りないので…
今回はケプラーの第3法則をデータとして利用しました(授業で)
D:太陽との距離
R:公転周期
K:定数
みたいな法則があります。
今、この式が
となっていて、DとRだけデータがあるものとする。
D R 水星 36.00 88.0 金星 67.25 224.7 地球 93.00 365.3 火星 141.75 687.0 木星 483.80 4332.1
このとき、上記式に対数変換を施して、
という自乗誤差(目的関数)の式を考える。この目的関数を最小化するようなを求めることで、最初に出てきたケプラーの第三法則を満たすようなパラメータを求めることができる(はず)である。
この時の手法として、前回のの結果を改善するようにを更新していくのが反復改善法になる。
今回はその中でも2つの方法を考える。
1.最急降下法
基本的な流れは
・初期パラメータを決める
・勾配ベクトルを計算
・勾配ベクトルが小さければ終了
・探索方向計算:を計算
・探索幅計算:を計算
・パラメータ更新
・反復数更新
このとの決め方が手法によって異なる
①単純な最急降下法
,
=定数
②慣性項付き最急降下法
=定数
=定数
この2つは定数があるためこの値によって結果が左右される
※以後、
である。
case1:単純最急降下法,
case2:単純最急降下法,
case3:慣性項付き最急降下法,
そして、定数を使わずに最急降下法を行うのが
③最急降下法
,
case4:最急降下法
いい感じですね。ただし、これが非線形になると一気にダメになってしまいます(画像ないですすいません)
さらに非線形にも対応できるのが
④共役勾配法
…式大変だなぁ…
プログラム
#include <iostream> #include <fstream> #include <cstdlib> #include <cmath> #include <iomanip> #include <vector> #include <algorithm> #include<sys/times.h> struct tms tms_buffer; double ST, UT, st, ut, tt; int hour, minute; double second; long ticks; #define INIT_TIME \ ticks = sysconf(_SC_CLK_TCK);\ times (&tms_buffer);\ ST = (double) tms_buffer.tms_stime ;\ UT = (double) tms_buffer.tms_utime #define FORMAT_TIME_AND_PRINT_TIME \ times (&tms_buffer);\ st = ((double) tms_buffer.tms_stime - ST)/ ticks;\ ut = ((double) tms_buffer.tms_utime - UT)/ ticks;\ tt = st + ut;\ hour = (int) (tt / 3600.0);\ minute = (int) (tt / 60.0) - hour * 60;\ second = (tt) - 60 * (int) (tt / 60.0);\ printf("(%dh %dm %fs)\n", hour, minute, second); using namespace std; int Dm, Dn; double *VecY, **MatX, *VecG,*VecTheta,**MatT; std::vector <vector <double> > MatTheta; void readValue(char *fn1){ ifstream fin; int i, j, l, f; double v,s; fin.open(fn1); if(!fin){ cerr << "ERROR: Failed to open file" << endl; exit(1); } fin >> Dm >> Dn; VecY = new double [Dm]; MatX = new double *[Dm]; for(i = 0; i < Dm; i++){ MatX[i] = new double[2]; fin >> v >> s; VecY[i] = log(s); MatX[i][0] = 1.0; MatX[i][1] = log(v); } fin.close(); } void initValue(){ int i, j, k; double v; VecG = new double [Dm]; MatTheta.resize(Dn); MatTheta[0].push_back(-1.3);//初期ベクトル MatTheta[1].push_back(1.3); } void trans() { int i,j; MatT = new double *[2]; for(i=0;i<2;i++){ MatT[i] = new double[Dm]; } for(i=0;i<Dm;i++){ for(j=0;j<2;j++){ MatT[j][i] = MatX[i][j]; } } } void calGradient(){ int i, j, k; double **XtX,*XtY,*XtXTheta,sum,Ramba,sum2,GtG,*GtXtX,*Delta,sum3,Beta; XtX = new double *[2]; XtY = new double [2]; XtXTheta = new double [2]; GtXtX = new double [2]; Delta = new double [2]; for(i=0;i<2;i++){ XtX[i] = new double [2]; for(j=0;j<2;j++){ for(k=0,sum=0.0,sum2=0.0; k<Dm; k++){ sum += MatT[j][k]*VecY[k]; sum2 += MatT[i][k]*MatX[k][j]; } XtX[i][j] = sum2; XtY[j] = sum; } } for(i = 0; i < 100; i++){ for(j = 0; j < 2; j++){ for(k=0,sum=0.0; k<2; k++){ sum += XtX[j][k]*MatTheta[k][MatTheta[k].size()-1]; } XtXTheta[j] = sum; } for(k=0; k<2; k++){ VecG[k] = XtXTheta[k] - XtY[k]; } for(k=0,Ramba=0.0; k<2; k++){ Ramba += powf(VecG[k],2); } Ramba = sqrt(Ramba); if(Ramba < 0.00001) break; if(i==0){ for(j=0;j<2;j++) Delta[j] = (-1)*VecG[j]; }else{ for(k=0,sum=0.0,sum3=0.0;k<2;k++){ for(j=0,sum2=0.0;j<2;j++){ sum2 += Delta[j]*XtX[j][k]; } sum += GtXtX[k]*Delta[k]; sum3 += sum2*Delta[k]; } Beta = sum/sum3; for(k=0;k<2;k++) Delta[k] = (-1)*VecG[k] + Beta*Delta[j]; } for(k=0,sum=0.0,sum3=0.0;k<2;k++){ for(j=0,sum2=0.0;j<2;j++){ sum2 += Delta[j]*XtX[j][k]; } sum += VecG[k]*Delta[k]; sum3 += sum2*Delta[k]; } Ramba = (-1)*sum/sum3; //printf("%f\n",Ramba); for(j=0;j<2;j++){ sum = MatTheta[j][i] + Ramba*Delta[j]; printf("%f ",sum); MatTheta[j].push_back(sum); } printf("\n"); } delete [] XtX; delete [] XtY; delete [] XtXTheta; delete [] GtXtX; delete [] Delta; } int main(int argc, char **argv){ int i, j, k; double v; srand(time(0)); INIT_TIME; readValue(argv[1]); cout << "read OK\n"; initValue(); cout << "init OK\n"; trans(); calGradient(); cout << "cal OK\n"; FORMAT_TIME_AND_PRINT_TIME; delete [] VecY; delete [] VecG; delete [] VecTheta; delete [] MatX; delete [] MatT; }
結果
-1.243547 1.571232 -0.897498 1.499207 -0.897374 1.499800 -0.896618 1.499642 -0.896618 1.499644
プログラムはないですが非線形でも大丈夫です
次回はより高速(?)な反復改善法であるニュートン法、さらに、別の視点からパラメータ推定を行うEMアルゴリズムへと進む予定。
※これは手法とMatlabが使えるようになるのが目的なので、面倒な時はMatlabのコード載せます。
まだまだ勉強不足です。