黄金分割法のRcppでの実装例

1.はじめに

表題の通り．Rだと計算に時間がかかる際にC++ライクに書いて実行速度を挙げる対処法がよく利用されます．前回の記事(囲い込みによる非線形最適化入門(私が))の黄金分割法のRcppによる実装を行いました．

2. サンプルプログラム

cpp, goldensection.cpp

#include<Rcpp.h>
using namespace Rcpp;
//[[Rcpp::export]]

List Gsearch(Function f, List interval, double tol = 1e-6, int niter = 1000){
  double upper, lower, gam, ai, bi;
  List output(3), fs(2);
  lower = interval[0];
  upper = interval[1];
  gam = (3-sqrt(5))/2;
  //std::cout << gam << std::endl;
  ai = lower + gam*(upper - lower);
  bi = lower + (1-gam)*(upper - lower);
  output[2] = niter;
  for(int i = 0; i < niter; i++){
    //std::cout << lower << " " <<  ai << " " << bi << " " << upper << std::endl;
    double fai = as<double>(f(ai));
    double fbi = as<double>(f(bi));
    if(fai <= fbi){
      upper = bi;
      bi = ai;
      ai = lower + gam * (upper - lower);
    }else{
      lower = ai;
      ai = bi;
      bi = lower + (1-gam)*(upper - lower);
    }
    if(abs(upper - lower) < tol){
      output[2] = i;
      break;
    }
  }
  output[0] = lower;
  output[1] = upper;
  return output;
}

f(ai) <= f(bi)の評価の際に，そのまま評価すると，結果が安定しないという点に詰まりました．f(ai),f(bi)の返り値の型をdoubleに変換することで解決しました．