非線形(連立)方程式をOptimLibで数値解析．

OptimLibを使って，非線形方程式を解くためのc++コードです（OptimLibのインストール・最適化はこちらを参照ください）．ここで紹介するのは，微分も使わない力技な解き方です．OptimLibで提供されている関数のうち，微分なしで最適化問題を解くoptim::nm, optim::pso, optim::deが対象となります．

今回は，こちらにある
$f(x) = x^2 - 5x + c - \exp x = 0$
を例として解きます．

optim::nmのインプットは，パラメータ・目的関数・外生変数（データ）・設定となっています．それぞれ，

• para: パラメータ．解を得たい内生変数をarma::vecで指定．(例では$x$)
• objfn: 目的関数．連立方程式の差分の和．（例では$|f(x)|$）
• opt_data: 外生変数やデータ．変更しながら計算したい変数があれば，objfn内ではなく，opt_dataで定義する．（例では$c$）
• settings: optim::nmの数値計算の設定．

をコードでは定義していきます．

コード解説

まずは，単なるヘッダー部分です．armadilloやoptim.hppを読み込みます．

main.cpp

#define ARMA_USE_CXX11
#define ARMA_USE_EXTERN_CXX11_RNG

#include <complex>
#include <iostream>
#include <random>
#include <fstream>
#include <ctime>
#include <sys/time.h>
#include <iomanip>

#include "../header_only_version/armadillo_bits/compiler_setup.hpp"
#include "../header_only_version/armadillo_bits/typedef_elem.hpp"
#include "../header_only_version/armadillo_bits/arma_rng_cxx11.hpp"
#include "../header_only_version/armadillo_bits/arma_rng.hpp"
#include "../header_only_version/optim.hpp"

using namespace std;

#include "function.h"

armadilloが利用している関数の関係上，色々とincludeが必要です．objfnをfunction.h内で定義します．

main.cppの冒頭はこちらです．

main.cpp

    arma::vec para = {4.5};// xの初期値
    var_exo opt_data;
    opt_data.ex0 = 2.0;// A=2.0と設定

    optim::algo_settings_t settings;
    settings.lower_bounds = {-10.0};// xの探索範囲を設定 
    settings.upper_bounds = { 10.0};

objfnとopt_dataの定義はこちらです．

function.h

struct var_exo{ //外生変数は構造体として定義
    double ex0;
};

// 目的関数
double objfn(const arma::vec& para, arma::vec* grad_out, void* opt_data){
    var_exo* objfn_data = reinterpret_cast<var_exo*>(opt_data);
    double icept = objfn_data->ex0;// 可読性向上のため，構造体からdoubleへ．
    
    double out = para(0)*para(0) - 5.0 * para(0) + icept - exp(para(0));//対象式
    cout << out << " " << para(0) <<endl;//目的関数値とパラメータの確認
    return abs(out);
}

目的関数値とパラメータの確認はなくても，もちろんOKです．
optimの関数では最小化計算をしているので，目的関数値は絶対値で与えます．
最小化問題を解いており，局所解回避の保証はないので，必ずf(x)=0となるxを得られるとは限りません．outの値が0になれば解けたことになります．

上記の構造体及びobjfnを使うmainの後半はこちらです．

main.cpp

    bool success;
    success = optim::pso(para, objfn, &opt_data, settings); // optimの戻り値はbool

    if (success){
        cout << endl << "Estimation was completed successfully." << endl;
    } else cout << endl << "Estimation was completed UNsuccessfully." << endl;
    
    cout << std::fixed << std::setprecision(10) << "solution " << para(0) << endl; // 念のため10桁表示．

最適化計算の収束有無はsuccessで確認できます．最後に，xの解を表示します．

結果

g++ -std=c++14 -O3 main.cpp -lopenblas -llapack -larmadillo
./a.out
でコンパイル・実行します．

optim::nmは探索が粗すぎて，うまくいきませんでした．
optim::psoでは，

-118.94 4.78682
-71.3249 4.25048
-110.49 4.71067
-138.969 4.94667
-57.6023 4.01927
-82.4725 4.4049
-71.3807 4.25132
-57.7837 4.02271
～中略～
0 0.168904
0 0.168904
0 0.168904
4.44089e-16 0.168904
0 0.168904

Estimation was completed successfully.
solution 0.1689043579

といった感じで表示されます．
解はx=0.1689043579です．参考とも解は一致しています．
（optim::deでも解けました）

連立方程式(例：f(x,y)=0, g(x,y)=0)にする場合は，パラメータをxyのvecで構成して，objfnの戻り値を|f(x,y)|+|g(x,y)|とすればOKのはずです．

参考・所感

参考：

• 非線形方程式の数値解法
  
• Pythonで方程式を解く方法

所感：
微分を得るのが難しい問題や最適化のパラメータの調整を最小限にして，c++で解くことを目的にOptimLibを導入したので，わりと思惑通りです．

以上です．お粗末さまでした．