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

最急降下法と共役勾配法

前期の授業で習ったまま放置してたのでちょっと復習。

最適化問題を解く際に利用される手法。

数式を書くには体力が足りないので…

今回はケプラーの第3法則をデータとして利用しました(授業で)

D:太陽との距離
R:公転周期
K:定数

D = KR^{2/3}

みたいな法則があります。

今、この式が

D = KR^{\alpha}となっていて、DとRだけデータがあるものとする。

    D    R
水星 36.00  88.0
金星 67.25  224.7
地球 93.00  365.3 
火星 141.75  687.0
木星 483.80  4332.1

このとき、上記式に対数変換を施して、

E(\theta)=\frac{1}{2}\sum^{N}_{i=1}(d_{i}-(\theta_{1}+\theta_{2}r_{i}))^{2}

 {\bf \theta}=(\theta_{1},\theta_{2})^{T}

logD=d, logK = \theta_{1}, \alpha= \theta_{2}, logR = r

という自乗誤差(目的関数)の式を考える。この目的関数を最小化するような\hat{{\bf \theta}を求めることで、最初に出てきたケプラーの第三法則を満たすようなパラメータを求めることができる(はず)である。

この時の手法として、前回の{\bf \theta}^{(t-1)}の結果を改善するように{\bf \theta}^{(t)}を更新していくのが反復改善法になる。

今回はその中でも2つの方法を考える。

1.最急降下法

基本的な流れは
・初期パラメータ{\bf \theta}^{(1)}を決める
・勾配ベクトル{\bf g}を計算
・勾配ベクトルが小さければ終了
・探索方向計算:{\bf \delta}^{(t)}を計算
・探索幅計算:\lambda^{(t)}を計算
・パラメータ更新{\bf \theta}^{(t+1)}={\bf \theta}^{(t)}+\lambda^{(t)}{\bf \theta}^{(t)}
・反復数更新t=t+1

この{\bf \delta}^{(t)}\lambda^{(t)}の決め方が手法によって異なる

①単純な最急降下法
{\bf \delta}^{(t)}=-{\bf g},
\lambda=定数
②慣性項付き最急降下法
{\bf \delta}^{(t)}=-{\bf g}+\frac{\alpha}{\lambda}({\bf \theta}^{(t)}-{\bf \theta}^{(t-1)})
\lambda=定数
\alpha=定数

この2つは定数があるためこの値によって結果が左右される

※以後、
{\bf \theta}^{(1)} = (-1.3,1.3)^{T}である。


case1:単純最急降下法\lambda=0.05
f:id:A_Koide0519:20111008231714p:image

case2:単純最急降下法\lambda=0.005
f:id:A_Koide0519:20111008231715p:image

case3:慣性項付き最急降下法\lambda=0.005,\alpha=0.9
f:id:A_Koide0519:20111008231716p:image

そして、定数を使わずに最急降下法を行うのが
最急降下法
{\bf \delta}^{(t)}=-{\bf g},
\lambda^{(t)}=\frac{{\bf g}^{T}{\bf g}}{{\bf g}^{T}{\bf X}^{T}{\bf X}{\bf g}}

case4:最急降下法
f:id:A_Koide0519:20111008231717p:image
いい感じですね。ただし、これが非線形になると一気にダメになってしまいます(画像ないですすいません)

さらに非線形にも対応できるのが
共役勾配法
…式大変だなぁ…

プログラム

#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

f:id:A_Koide0519:20111008231718p:image

プログラムはないですが非線形でも大丈夫です

f:id:A_Koide0519:20111008231719p:image


次回はより高速(?)な反復改善法であるニュートン法、さらに、別の視点からパラメータ推定を行うEMアルゴリズムへと進む予定。

※これは手法とMatlabが使えるようになるのが目的なので、面倒な時はMatlabのコード載せます。

まだまだ勉強不足です。