前記事でインストールしたOptimLibを使って計算してみます．とりあえず，最急降下法とニュートン法を試したサンプルコード(←のページの一番下)のロジスティック回帰¹に少し手を加えて計算してみます．手を加えた内容は乱数シードの導入だけですが，リンクの設定に苦労しました．armadillo関連のエラーが多かったのですが，documentationと検索により，頑張りました．

コード解説

まずは，include部分（不必要なものもあるかもしれません）．

test.cpp

#define ARMA_USE_CXX11
#define ARMA_USE_EXTERN_CXX11_RNG

#include <iostream>
#include <complex>
#include <random>
#include <fstream>
#include <ctime>
#include <sys/time.h>
using namespace std;
#include "./armadillo_bits/compiler_setup.hpp"
#include "./armadillo_bits/typedef_elem.hpp"
#include "./armadillo_bits/arma_rng_cxx11.hpp"
#include "./armadillo_bits/arma_rng.hpp"
#include "optim.hpp"

ロジスティック関数(sigm)の定義とOptimLibの最適化関数(ここではgdとnewton)にインプットするための構造体の定義です．sigmの引数のXはここでは(4000,5)のmatrixです．そのまま，行ごとに足して，(4000,1)のmatrixを返すようです．

test.cpp

// sigmoid function
inline
arma::mat sigm(const arma::mat& X)
{
    return 1.0 / (1.0 + arma::exp(-X));
}

// log-likelihood function data
struct ll_data_t
{
    arma::vec Y;
    arma::mat X;
};

目的関数（対数尤度）の定義です．//A:は僕が追加したコメントです．

test.cpp

// log-likelihood function with hessian
double ll_fn_whess(const arma::vec& vals_inp, arma::vec* grad_out, arma::mat* hess_out, void* opt_data)
{
    ll_data_t* objfn_data = reinterpret_cast<ll_data_t*>(opt_data);

    arma::vec Y = objfn_data->Y; //A: 構造体を元に戻す
    arma::mat X = objfn_data->X;

    arma::vec mu = sigm(X*vals_inp);

    const double norm_term = static_cast<double>(Y.n_elem);

    const double obj_val = - arma::accu( Y%arma::log(mu) + (1.0-Y)%arma::log(1.0-mu) ) / norm_term; //A: 目的関数の算出（対数尤度の平均値）

    if (grad_out) //A: grad_outは目的関数の微分値．最適解の探索で利用．
    {
        *grad_out = X.t() * (mu - Y) / norm_term;
    }

    if (hess_out) //A: hess_outはヘシアン．最適解の探索で利用．
    {
        arma::mat S = arma::diagmat( mu%(1.0-mu) );
        *hess_out = X.t() * S * X / norm_term;
    }

    return obj_val; //A: 目的関数を返す．
}

// log-likelihood function for Adam //A: gdで利用．ヘシアンがいらない．
double ll_fn(const arma::vec& vals_inp, arma::vec* grad_out, void* opt_data)
{
    return ll_fn_whess(vals_inp,grad_out,nullptr,opt_data);
}

main関数の設定．まずは，冒頭のデータ生成と真値の設定です．元のデータは乱数シード固定でしたので，乱数シードを固定あるいはランダムに設定できるよう，追加しました．//CA:は僕の加筆コードです．

test.cpp

int main()
{
    int n_dim = 5;     // dimension of parameter vector
    int n_samp = 4000; // sample length

    //arma::arma_rng::seed_type va = 2; //CA: 乱数シードを自分で指定・記録する場合はこの4行を使う． 
    //arma::arma_rng tt; //CA: class arma_rngのオブジェクトのせんげん 
    //tt.set_seed(va); //CA: シードの設定
    //cout << "va " << va << endl; //CA: シードを書き出す
    arma::arma_rng::set_seed_random(); //CA: ランダムの乱数シードの設定(自分で指定する場合はコメントアウト)
    arma::mat X = arma::randn(n_samp,n_dim);//A: 正規分布から付与
    cout << "X[0][0-2] " << X(0,0) << " " << X(0,1) << " " << X(0,2) << endl; //CA:  試しに書き出し．行・列の指定の仕方はarmaの方法です．
    arma::vec theta_0 = 1.0 + 3.0*arma::randu(n_dim,1); //A: 真のパラメータ：theta_0 

    arma::vec Y(n_samp); //A: Initialize a vector
    for (int i=0; i < n_samp; i++) //A: 0/1の分類
    {
        Y(i) = ( arma::as_scalar(arma::randu(1)) < mu(i) ) ? 1.0 : 0.0; //A: create a true (disrete) result 
    }

インプットデータであるXと0/1分類結果であるYを構造体としてまとめる（最適化関数に渡すため）．

test.cpp

    // fn data and initial values
    ll_data_t opt_data;
    opt_data.Y = std::move(Y);
    opt_data.X = std::move(X); //A: optim::gd等の最適化関数にopt_dataを渡す．

    arma::vec x = arma::ones(n_dim,1) + 1.0; // initial values //A: 求めたいパラメータの初期値．また，最適化計算後はこのxに計算結果が上書きされる．

最適化計算とその結果の出力です．setting.gd_methodの参照はこちらです．

test.cpp

    // run Adam-based optim
    optim::algo_settings_t settings;

    settings.gd_method = 6; //A: 値の更新手法の設定
    settings.gd_settings.step_size = 0.1; //A: ステップサイズα 

    std::chrono::time_point<std::chrono::system_clock> start = std::chrono::system_clock::now();
 
    bool success = optim::gd(x,ll_fn,&opt_data,settings); //A: 最急降下法による計算

    std::chrono::time_point<std::chrono::system_clock> end = std::chrono::system_clock::now();
    std::chrono::duration<double> elapsed_seconds = end-start;

    if (success) { //A: 最適化できたかどうか
        std::cout << "Adam: logit_reg test completed successfully.\n"
                  << "elapsed time: " << elapsed_seconds.count() << "s\n";
    } else {
        std::cout << "Adam: logit_reg test completed unsuccessfully." << std::endl;
    }
 
    arma::cout << "\nAdam: true values vs estimates:\n" << arma::join_rows(theta_0,x) << arma::endl; //A: 真値と推定結果の比較．

最後に，ニュートン法もgdと同様です．(ニュートン法は初期パラメータによっては求解できませんでした)

test.cpp

    // run Newton-based optim

    x = arma::ones(n_dim,1) + 1.0; // initial values //A: reset

    start = std::chrono::system_clock::now();
 
    success = optim::newton(x,ll_fn_whess,&opt_data);

    end = std::chrono::system_clock::now();
    elapsed_seconds = end-start;

    //
 
    if (success) {
        std::cout << "newton: logit_reg test completed successfully.\n"
                  << "elapsed time: " << elapsed_seconds.count() << "s\n";
    } else {
        std::cout << "newton: logit_reg test completed unsuccessfully." << std::endl;
    }
 
    arma::cout << "\nnewton: true values vs estimates:\n" << arma::join_rows(theta_0,x) << arma::endl;
 
    return 0;
}

最後にコンパイル．前記事の通り，ヘッダーオンリーversionを使うので，optim/header_only_versionのフォルダ・ファイルを入れたフォルダの中にtest.cppを保存して，
g++ -std=c++14 test.cpp -lopenblas -llapack -larmadillo
でコンパイルしました．

乱数設定しているの結果は書きませんが，もちろんきちんと求解出来ていました．

感想（メリットデメリット）

(+) 構造体でデータを引き受け，目的関数を算出する関数を定義することで，複数インプット・複数パラメータの最適化問題が求解できる．
(+) 連続量のパラメータの探索アルゴリズムはライブラリで与えられているので，再設計不要．しかも，ほぼ同じ形で色々な手法が試せる．
(+) パラメータ探索範囲の制約もどの最適化手法でも設定できそう．→(10/24追記) 制約を与えても，その範囲に収まってくれないこともある．なぜ？
(-) ただし，求解方法によって必要となる微分や二階微分(ヘシアン)は自分で与える必要がでてきます（ll_fn_whess関数のgrad_outとhess_out）．
(+) とはいっても，Differential Evolution, Particle Swarm Optimization, Nelder-Meadは微分を使わないので，気にならない．
(-) armadilloの構文に多少慣れる必要がありそう．

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

追記（最適化関数のパラメータ）

ロジスティック回帰の　参考１　参考２ ↩