バックトラックを用いた最急降下法とニュートン法の実装


大学で連続最適化の勉強をすると降下法の代表例として最急降下法とニュートン法を学ぶと思います。また、これらの方法に、直線探索法の代表例であるバックトラック法を組み込むことが効率的であるということも学ぶと思います。私も先日、これらの方法について学んだところなので、実際にC++でプログラムを実装することにしました。

##降下法とは？
連続最適化は大抵の場合、目的関数の最小化を目標とします。例えば、
>$\mathrm{minimize} \ \ f(x)$
>$\mathrm{subject \ to}\ x\in{X}$

と定式化することが多いです。この時$X$を制約集合といい、最適化問題に於いては、この制約集合に包含される点の中から最適解を探すのですが、$X=\mathbb{R}^n$である場合には実数を成分とした全てのn次元ベクトルを制約集合に含むので、制約集合には実質的に制限はありません。今回の実装は対象とする最適化問題がこの場合であることを前提とします。

降下法は、このような場合に関数の値を減少させる点列$\\{x^0,x^1,x^2,...\\}$を生成する反復法で、$i<j$ならば$f(x^i)＞f(x^j)$となるように点列を生成します。具体的な手続きは以下のようなものです。

1. $x^0$(初期点),$\epsilon>0$を決めて$k:=0$とする.
2. $x^k$が$||\nabla{f}(x^k)||<\epsilon$(終了条件)を満たしているときは$x^k$を解として終了する.
3. $(d^k,\nabla{f}(x^k))<0$が満たされるように$d^k$(降下方向)を決める.
4. ステップサイズ$t_k$を決める.
5. $x^{k+1}←x^k+t_kd^k$とし$k←k+1$して終了条件の判定のステップに戻る.

この時、最急降下法は降下方向を$d^k:=-\nabla{f}(x^k)$とし、ニュートン法では$d^k:=-\nabla^2f(x^k)^{-1}\nabla{f}(x^k)$としています。

また、バックトラック法はステップサイズを決めるのに理論的にも実用的にも優れている方法です。具体的な手続きは以下のようなものです。

1. $\xi\in(0,1),\rho\in(0,1),\bar{t}>0$(初期ステップサイズ)を選び、$\bar{t}=t$とする.
2. $f(x^k+t_kd^k)\leq{f}(x^k)+\xi{t}(d^k,\nabla{f}(x^k))$(アルミホのルール)が満たされるまで$t←\rho{t}$とすることを繰り返す.
3. $t_k=t$として終了する.

##実際に実装してみた。
まずは、バックトラック法を用いた最急降下法を実装してみます。

```cpp
#include<iostream>
#include<iomanip>
#include<vector>
#include<math.h>
using namespace std;

#define D 10 //初期点の次数を定める

double f(vector<double> x);
vector<double> grad_f(vector<double> x);

int main(){
    vector<double> x(D,1.0); //初期点は全てが要素１のベクトルとしている
    double epsilon=pow(10,-6);
    //最急降下法
    //反復は最大10000回としている
    for(double k=0;k<=10000;k++){ //反復は最大10000回
        const auto &z=grad_f(a);
      
        //終了条件
        double check=0;
        for(int i=0; i<D; i++) check=max(check, z[i]);
        double epsilon=pow(10,-6);
        if(check<=epsilon){ 
            double ans=f(a); 
            for(int i=0; i<D; i++) cout << a[i] << ' '; //最適解出力
            cout << endl;
            cout << ans << endl; //最適値出力
            break;
        }
　　　　　
        vector<double> d(D,0);
        //降下方向設定
        for(int i=0; i<D; i++) d[i]=-z[i];

        //バックトラック法
        vector<double> tmp(D);
        double xi=pow(10,-4);
        double rho=0.5;
        double t=1.0;
        while(true){ 
            for(int i=0; i<D; i++) tmp[i]=a[i]+t*d[i];
            double kep1=f(tmp);
            double kep2=f(a);
            double kep3=xi*t*(d[0]*z[0]+d[1]*z[1]);
            if(kep1<=kep2+kep3) break; //アルミホのルール
            t*=rho;
        }

        //点列の更新
        for(int i=0; i<D; i++) a[i]=tmp[i];
    }

double f(vector<double> x){
    double y;
    //目的関数で計算
    return y;
}

vector<double> f(vector<double> x){
    vector<double> y(D);
    //目的関数の勾配を計算
    return y;
}
```
以上がバックトラック法を用いた最急降下法のサンプルコードです。目的関数やその勾配については各々の最適化問題によって異なってくるので、問題ごとに適宜書き換えることにします。バックトラック法を用いたニュートン法は、今のコードを幾分か変更することで実装できます。

まず、先のコードに「目的関数のヘッセ行列$\nabla^2f(x)$を求める関数」と、「ヘッセ行列の逆行列求める関数」を追加します。具体的には、今回のサンプルコードの記法にのっとれば、まずmain関数の前に

```cpp
vector<vector<double> > hessian(vector<double> x);
vector<vector<double> > inverse(vector<vector<double> > h);
```
と関数のプロトタイプ宣言をしておいて、main関数の後に「目的関数のヘッセ行列$\nabla^2f(x)$を求める関数」と、「ヘッセ行列の逆行列求める関数」を追記して

```cpp
vector<vector<double> > hessian(vector<double> x){
    vector<vector<double> >　y(D, vector<double>(D));
    //目的関数のヘッセ行列を計算
    return y;
}

vector<vector<double> > inverse(vector<vector<double> > h){
    vector<vector<double> >　y(D, vector<double>(D));
    //ヘッセ行列の逆行列を計算
    return y;
}
```
とします。目的関数のヘッセ行列について問題ごとに異なるので変わるので適宜書き換えることにします。（個人的にはヘッセ行列の逆行列を求めるところで実装に詰まりそうだったので、別の機会に逆行列の求め方についても記事を書きます。）次に、サンプルコードの「降下方向の設定」の箇所を以下のように書き換えます。

```cpp
//降下方向設定
const auto h=hessian(x);
const auto h_=inverse(h);
for(int i=0; i<D; i++){
    for(int j=0; j<D; j++){
        d[i]+=(h_[i][j]*z[j]);
    }
    d[i]=-d[i];
}
```
以上２点の変更で、バックトラック法を用いたニュートン法も実装できます。

###最急降下法とニュートン法どっちがいい？
最後に、最急降下法とニュートン法について両者の実力が発揮できる使い所をを付け加えておきます。$x^*\in\mathbb{R}^n$(最適解)について次数nが小さいときはニュートン法の方が最急降下法よりも圧倒的に早く最適解に収束します。ところが、ヘッセ行列の逆行列を求めるには次数nに対して$\mathcal{O}(n^3)$の計算量が必要になるので、最適解の次数が大きい場合はニュートン法では逆行列の計算負荷が大きくなります。そうした場合は、目的関数の勾配までしか必要としないため次数が大きくなっても計算負荷が増えにくい最急降下法を使うことが無難です。

双方の計算手法については、各々の最適化問題の目的関数の特徴を見て使い分けるのが望ましいでしょう。