2011-10-08

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

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

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

今回はケプラーの第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つの方法を考える。

１．最急降下法

基本的な流れは
・初期パラメータ ${\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$

case2:単純最急降下法， $\lambda=0.005$

case3:慣性項付き最急降下法， $\lambda=0.005,\alpha=0.9$

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

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のコード載せます。

まだまだ勉強不足です。

KIWAM_KEN_DIARY

最急降下法と共役勾配法