はじめに
「最適化問題を解きたいな」と思い,最初に思いついたのが最急降下法でしたので,解いてみました.
最急降下法とは
関数$f(x)$の導関数$f'(x)$を用いて最小値もしくは最大値を求める手法(極値となることもある)で,勾配法の一つ.
ある$x$から計算を始めて,変数$x$は以下のように更新する.ここで,$\alpha$は定数とする.
$$
x_{next} = x - \alpha \cdot f'(x)
$$
$\alpha$の値によって$x$の変化速度が変わる.注意しなくてはならないのが,$\alpha$が小さすぎると変化が遅くなりすぎ,大きすぎると発散してしまうことである.
使用する関数や定数について
実際に作成したプログラムが以下である関数および定数に特に意味はない.
ここで,解く関数は,
$$
f(x) = ax^2 + bx + c
$$
とし,各定数は$a=1,b=2,c=3$とした.
初期条件は$x=-1.5$と,終了条件は$-0.00001<f'(x)<0.00001$とした.
導関数を算出するとき$\Delta x=0.0001$とし,中心差分を用いた.
手計算によると,
$$
f(x) = ax^2 + bx + c = (x+1)^2+2
$$
となり,$x=-1$の時に最小値$f(x)=2$をとる.
プログラム
上述した更新式の定数$\alpha=0.0001$とした.
#include <iostream>
#include <fstream>
#include <cmath>
using namespace std;
constexpr double ALPHA = 0.0001;
constexpr double INI_X = -1.5;
constexpr double TH = 0.00001;
constexpr double DX = 0.0001;
constexpr double A = 1.00;
constexpr double B = 2.00;
constexpr double C = 3.00;
inline double func(const double x){
return A*x*x + B*x + C;
}
inline double diff(const double x){
return ( func(x+DX)-func(x-DX) )/(2.0*DX);
}
int main(){
double x = INI_X;
bool loopFlag = true;
ofstream fout ("gradiendDescent.txt");
fout << x << " " << func(x) << endl;
double difference;
while(loopFlag){
difference = diff(x);
x -= ALPHA*difference;
if(-TH<difference && difference<TH) loopFlag = false;
fout << x << " " << func(x) << endl;
}
cout << "x : " << x << endl;
cout << "extremum : " << func(x) << endl;
return 0;
}
そして,出力結果と出力ファイルを描画したものが以下である.
x : -1
extremum : 2
解けたようだ.
収束しない場合について
振動
上のプログラムの定数は収束するような値を用いたのだが,定数の設定によっては収束せずに振動する事がある.
振動するように少々強引に定数を以下にして,cout
で出力して変数の状態を確認すると,
constexpr double ALPHA = 1.0;
constexpr double DX = 0.1;
// 出力確認をループ内に入れる
cout << "now : " << x << ", dif : " << difference << endl;
now : -1.5, dif : 1
now : -0.5, dif : -1
now : -1.5, dif : 1
now : -0.5, dif : -1
now : -1.5, dif : 1
now : -0.5, dif : -1
now : -1.5, dif : 1
now : -0.5, dif : -1
now : -1.5, dif : 1
now : -0.5, dif : -1
now : -1.5, dif : 1
これは同じ値を繰り返している.
更新された$x$が関数の軸を対象とする位置となる場合,次の更新式も関数の軸を対象とする位置(元の位置)に戻ってくるため,振動する.
発散
収束でも振動でもなく,発散する場合もあり,同様に確かめた.
constexpr double ALPHA = 10.0;
constexpr double DX = 0.1;
// 出力確認をループ内に入れる
cout << "now : " << x << ", dif : " << difference << endl;
now : 8.5, dif : -1
now : -181.5, dif : 19
now : 3428.5, dif : -361
now : -65161.5, dif : 6859
now : 1.23805e+06, dif : -130321
now : -2.35229e+07, dif : 2.4761e+06
now : 4.46936e+08, dif : -4.70459e+07
now : -8.49178e+09, dif : 8.93872e+08
now : 1.61344e+11, dif : -1.69836e+10
now : -3.06575e+12, dif : 3.2271e+11
now : 5.82986e+13, dif : -6.13643e+12
now : -1.12368e+15, dif : 1.18197e+14
now : 2.70238e+16, dif : -2.81475e+15
now : 2.70238e+16, dif : 0
x : 2.70238e+16
extremum : 7.30287e+32
プログラム自体は終了しているのだが,この値は間違い.
更新された$x$が軸より遠ざかってしまった場合は,次の更新式においても,より遠ざかってしまうため,発散してしまう.
おわりに
結構場当たり的にやりましたが,収束・振動・発散を確認できたので楽しかったです.
記憶を辿りながら作ったので,何か間違いなどあれば教えていただけると幸いです.