ニュートン法，準ニュートン法

この時期になって就活のメールが届くようになってきました。
正直まだ早いなぁという気持ちですが、そろそろ準備しないとまずいんですかね。M2の皆様のアドバイス大募集中です。（身近にM2の先輩がいないので…割と本気でお願いします）

さて、前回は最急降下法と共役勾配法でしたが、今回はニュートン法・準ニュートン法を忘備録として書いておきます。

ニュートン法では、共役勾配法と同様に目的関数のヘッセ行列を使って探索を行います。

ニュートン法は、線形問題の場合、目的関数のパラメータ数以下の反復で解が収束することが保証されており、式もそれほど難しいものではありません。ただし、初期値を適切に決めないと収束は保証されません。

今回は、前回まで使っていたケプラーの第三法則を非線形問題としたものを考えます。

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

$D = KR^{2/3}$

　　　　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}R_{n}^{\theta}_{2}))^{2}$

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

となります。

ニュートン法の基本的な流れは
・初期パラメータ ${\bf \theta}^{(1)}$ を決める
・勾配ベクトル ${\bf g}$ を計算
・勾配ベクトルが小さければ終了
・探索方向計算： ${\bf \delta}^{(t)}$ を計算
・探索幅は定数： $\lambda=1$
・パラメータ更新 ${\bf \theta}^{(t+1)}={\bf \theta}^{(t)}+\lambda{\bf \theta}^{(t)}$
・反復数更新 $t=t+1$
です。肝は
${\bf \delta}^{(t)}$
を求めるところですが、

ここで、勾配ベクトル
$\frac{\partial E{\bf(\theta)}}{\partial{\bf \theta}} = {\bf g(\theta)}$

を求め、さらにヘッセ行列
${\bf H}({\bf \theta})=\frac{\partial}{\partial{\bf \theta}}{\bf g}({\bf\theta})^{T}$

を求めます。（本当は書くべきですが、これはなかなか大変…）

そして、これらが求まったとするならば、

${\bf \delta}^{(t)}=-{\bf H}({\bf \theta})^{-1}{\bf g(\theta)}$

として、パラメータを更新します。

ここでニュートン法で厄介なのは、 ${\bf H}({\bf \theta})^{-1}$ を求めなければならない部分です。僕がもしやるならば、今のスペックではLU分解をするか、パワー法で固有値と固有ベクトルを求めて…という感じになると思いますが、これはなかなか大変です。

そこで、逆行列を使わず、類似逆行列を用いて反復させるのが準ニュートン法です。

この類似逆行列の算出にはいくつか方法がありますが、今回はBFGS法(Broyden-Flecher-Goldfarb-Shano法)を使いました。

ここで、
$\triangle{\bf \theta}(k)={\bf \theta}^{t+1}-{\bf \theta}^{t}$
$\triangle{\bf g}(k)={\bf g}^{t+1}-{\bf g}^{t}$
${\bf H}(0)=I,k=0$
とすると、

${\bf H}(k+1)=(I-\frac{\triangle{\bf \theta}(k)\triangle{\bf g}(k)^{T}}{\triangle{\bf \theta}(k)^{T}\triangle{\bf g}(k)}){\bf H}(k)(I-\frac{\triangle{\bf g}(k)\triangle{\bf \theta}(k)^{T}}{\triangle{\bf \theta}(k)^{T}\triangle{\bf g}(k)})+\frac{\triangle{\bf \theta}(k)\triangle{\bf \theta}(k)^{T}}{\triangle{\bf \theta}(k)^{T}\triangle{\bf g}(k)}$

となります。

#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, *VecX,*VecTheta,**Hess;

std::vector <vector <double> > MatTheta;
std::vector <vector <double> > Gradient;

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];
  VecX = new double [Dm];
  for(i = 0; i < Dm; i++){
    fin >> v >> s;
    VecY[i] = s;
    VecX[i] = log(v); 
  }
  fin.close();
}

void initValue(){
  int i, j, k;
  double v;
  VecTheta = new double [Dn];
  Hess = new double *[Dn];
  for(i=0;i<Dn;i++){
    Hess[i] = new double[Dn];
  }
  MatTheta.resize(Dn);
  Gradient.resize(Dn);
  MatTheta[0].push_back(0.3);//初期ベクトル
  MatTheta[1].push_back(1.45);
}

void calGradient(){
  int i, j, k, l;
  double sum,Ramba,sum2,*VecE,**MatR2,*delta_Theta,*delta_Grad,TG,**GTt,**TGt,**Idn,**Hess2,**TTt,**R2R2,*VecGrad,**MatR2t;
  delta_Theta = new double [2]; //今回の探索方向-前回の探索方向
  delta_Grad = new double [2]; //今回の勾配ベクトル-前回の勾配ベクトル
  VecE  = new double [Dm]; //
  MatR2 = new double *[2]; //
  MatR2t = new double *[Dm];
  Idn = new double *[2];
  GTt = new double *[2];
  TGt = new double *[2];
  TTt = new double *[2];
  Hess2 = new double *[2];
  R2R2 = new double *[2];
  VecGrad = new double [2];
  for(k=0;k<Dm;k++) MatR2t[k] = new double [2];
  for(k=0;k<Dn;k++){
    MatR2[k] = new double [Dm];
    Idn[k] = new double [2];
    GTt[k] = new double [2];
    TGt[k] = new double [2];
    Hess2[k] = new double [2];
    TTt[k] = new double [2];
    R2R2[k] = new double [2];
    Idn[k][k] = 1;
    for(j=k+1;j<Dn;j++){
      Idn[k][j] = 0;
      Idn[j][k] = 0;
    }
  }
  for(i = 0; i < 100; i++){
    for(j = 0; j < Dm; j++){
      VecE[j] = VecY[j] - MatTheta[0][i]*exp(MatTheta[1][i]*VecX[j]);
      MatR2[0][j] = exp(MatTheta[1][i]*VecX[j]);
      MatR2[1][j] = MatTheta[0][i]*exp(MatTheta[1][i]*VecX[j])*VecX[j];
    } 
    for(j = 0; j < Dm; j++){
      for(k=0; k<2; k++){
	MatR2t[j][k] = MatR2[k][j];
      }
    }
    for(j=0; j<2; j++){
      for(k=0,sum=0.0; k<Dm; k++){
	sum+=MatR2[j][k]*VecE[k];
      }
      Gradient[j].push_back((-1)*sum);
    }
    for(k=0,Ramba=0.0; k<2; k++){
      Ramba += powf(Gradient[k][i],2);
    }
    Ramba = sqrt(Ramba);
    if(Ramba < 0.00001) break;
    if((i+1)%2 == 1){
      for(k=0;k<Dn;k++){
	Hess[k][k] = 1;
	for(j=k+1;j<Dn;j++){
	  Hess[k][j] = 0;
	  Hess[j][k] = 0;
	}
      }
      //
    }else{
      for(j=0;j<Dn;j++){
	delta_Theta[j] = MatTheta[j][i] - MatTheta[j][i-1];
	delta_Grad[j] = Gradient[j][i] - Gradient[j][i-1];
      }
      for(j=0,TG=0.0;j<2;j++) TG += delta_Theta[j]*delta_Grad[j];
      for(j=0;j<2;j++){
	for(k=0;k<2;k++){
	  TGt[j][k] = delta_Theta[j]*delta_Grad[k];
	  GTt[j][k] = delta_Grad[j]*delta_Theta[k];
	}
      }
      for(j=0;j<2;j++){
        for(k=0;k<2;k++){
          TGt[j][k] = Idn[j][k] - TGt[j][k]/TG;
          GTt[j][k] = Idn[j][k] - GTt[j][k]/TG;
        }
      }
      for(j=0;j<2;j++){
	for(l=0;l<2;l++){
	  for(k=0,sum=0.0;k<2;k++){
	    sum += TGt[j][k] * Hess[k][l];
	  }
	  Hess2[j][l] = sum;
	}
      }
      for(j=0;j<2;j++){
        for(l=0;l<2;l++){
          for(k=0,sum=0.0;k<2;k++){
            sum += Hess2[j][k]*GTt[k][l];
          }
          Hess[j][l] = sum;
        }
      }
      for(j=0;j<2;j++){
        for(k=0;k<2;k++){
          TTt[j][k] = (delta_Theta[j]*delta_Theta[k])/TG;
        }
      }
      for(j=0;j<2;j++){
        for(k=0;k<2;k++){
          Hess[j][k] = Hess[j][k] + TTt[j][k];
	  printf("%f ",Hess[j][k]);
        }
	printf("\n");
      }     
      printf("\n");
    }
    for(j=0;j<2;j++){
      for(k=0,sum=0.0;k<2;k++){
	sum += Hess[j][k] * Gradient[k][i];
      }
      VecTheta[j] = sum*(-1);
    }
    //
    for(j=0;j<2;j++){
      for(l=0;l<2;l++){
	for(k=0,sum=0.0;k<Dm;k++){
	  sum += MatR2[j][k] * MatR2t[k][l];
	}
	R2R2[j][l] = sum;
      }
    }
    for(j=0;j<2;j++){
      for(k=0,sum=0.0;k<2;k++){
        sum += VecTheta[k] * R2R2[k][j];
      }
      VecGrad[j] = sum;
    }
    for(j=0,sum=0.0,sum2=0.0;j<2;j++){
      sum += Gradient[j][i]*VecTheta[j];
      sum2 += VecGrad[j]*VecTheta[j];
    }
    Ramba = (-1)*(sum/sum2);
    for(j=0;j<2;j++){
      VecTheta[j] = Ramba*VecTheta[j];
    }
    for(j=0,sum=0.0;j<2;j++) sum += powf(VecTheta[j],2);
    if(sqrt(sum) > 0.2){
      for(j=0;j<2;j++) VecTheta[j] = (0.2/sqrt(sum))*VecTheta[j];
    }
    for(j=0;j<2;j++){
      MatTheta[j].push_back(MatTheta[j][i] + VecTheta[j]);
      printf("%f\n",MatTheta[j][i] + VecTheta[j]);
    }
    printf("\n");
  }
  delete [] VecE;
  delete [] MatR2;
  delete [] delta_Theta;
  delete [] delta_Grad;
  delete [] GTt;
  delete [] TGt;
  delete [] Idn;
  delete [] Hess2;
  delete [] TTt;
  delete [] R2R2;
  delete [] VecGrad;
}

int main(int argc, char **argv){
  srand(time(0));
  INIT_TIME;
  readValue(argv[1]);
  cout << "read OK\n";
  initValue();
  cout << "init OK\n";
  calGradient();
  cout << "cal OK\n";
  FORMAT_TIME_AND_PRINT_TIME;
  delete [] VecY;
  delete [] VecTheta; 
  delete [] VecX;
  delete [] Hess;
}

結果

read OK
init OK
0.357727
1.556252

0.804700 -0.397930
-0.397930 0.196780

0.537005
1.467598

0.533594
1.456421

0.915761 -0.277751
-0.277751 0.084242

0.392323
1.499268

0.394765
1.505147

0.853189 -0.353918
-0.353918 0.146811

0.408218
1.499566

0.408214
1.499555

0.862998 -0.344112
-0.344112 0.137211

0.407436
1.499865

0.407436
1.499866

0.862367 -0.344515
-0.344515 0.137633

0.407437
1.499865

0.407437
1.499865

cal OK
(0h 0m 0.000000s)

 if((i+1)%2 == 1){
      for(k=0;k<Dn;k++){
	Hess[k][k] = 1;
	for(j=k+1;j<Dn;j++){
	  Hess[k][j] = 0;
	  Hess[j][k] = 0;
	}
      }
非線形の問題の場合、ヘッセ行列に初期のころの探索方向を記憶していると、探索がうまくいかないことがあります。

そこで、パラメータ数に応じて、ヘッセ行列をリセットしてあげるのが一般的な解法なようです。

また、探索幅の修正もしていて、なかなかずっしりとした中身になっています。しかもこの問題にしか適用できないので、きついですね…

行列の演算がもっと簡単にできるといいんですけどね。