読者です 読者をやめる 読者になる 読者になる

MyEnigma

とあるエンジニアのブログです。#Robotics #Programing #C++ #Python #MATLAB #Vim #Mathematics #Book #Movie #Traveling #Mac #iPhone

GoogleのC++最適化ツールCeres solverを使って複雑な最適化を解いてみた。

はじめての最適化

はじめての最適化

目次

はじめに

先日、上げた記事の続きです。

myenigma.hatenablog.com

 

今回は、もう少しCeresのシステムの詳しい説明と、

複雑な最適化を実際に実装してみたいと思います。

 

Ceresの基本的な部分は先程の記事を参照下さい。

Ceresにおける導関数の指定方法

Ceresは、他の最適化のライブラリと同様に、

評価関数の各変数における導関数を使って、最適化を実施します。

よって、精度が良く安定した最適化を実施するためには、

正しく導関数を与えることが重要になります。

 

大きく分けて2つの方法でCeresに導関数を与えることができます。

一つ目は数値演算的に導関数を与える方法、

二つ目は解析的に導関数を与える方法です。

 

数値演算的に導関数を与える方法

対象とするシステムにおいては、

先日の記事で紹介したようなテンプレートを使ったコスト関数を

設定できないことがあります。

例えば、コスト関数を計算する時に、

外部のライブラリを使っている時などです。

 

そんな時のために、数値演算的に導関数を与える方法があります。

まず初めに下記のように、数値導関数を与える構造体を作ります。

struct NumericDiffCostFunctor {
  bool operator()(const double* const x, double* residual) const {
    residual[0] = 10.0 - x[0];
    return true;
  }
};

先日の記事との違いは、テンプレートを使わなくしたことです。

続いて、上記の構造体を

NumericDiffCostFunctionというCeresのAPIを使うことにより、

Problemクラスに与えます。

CostFunction* cost_function =
  new NumericDiffCostFunction<NumericDiffCostFunctor, ceres::CENTRAL, 1, 1, 1>(
      new NumericDiffCostFunctor)
problem.AddResidualBlock(cost_function, NULL, &x);

先日の記事における、AutoDiffCostFunctionを使った場合と全く同じAPIですが、

導関数を数値的にどうのように計算するかを指定している所が違います。

(ceres::CENTRALの部分)

この導関数の計算の方法に関しては下記のAPIドキュメントが参考になります。

 

ただ、一つの注意点として、

ceresの開発者達は、この数値導関数ではなく、

先日の記事にある自動導関数(Automatic differentiation)の方法を使うことを推奨しています。

なぜなら、C++のテンプレートを使ったほうが、

計算の効率性が高く、数値導関数は、

計算が遅く、精度も悪くなり、収束も遅いことが原理的に示されているからです。

 

解析的に導関数を与える方法

前述の通り、多くの場合、自動導関数を使えば良いですが、

もし、評価関数が閉じた系になっており、

導関数を解析的に計算出来る場合は、

そちらを使った方が効率的になることがあります。

 

このような場合は、

残差項とヤコビ行列の両方をコードに記述する必要があります。

コードとして与える場合は、

CostFunctionクラスか、SizedCostFunctionクラスを継承したクラスを設定する必要があります。

class QuadraticCostFunction : public ceres::SizedCostFunction<1, 1> {
 public:
  virtual ~QuadraticCostFunction() {}
  virtual bool Evaluate(double const* const* parameters,
                        double* residuals,
                        double** jacobians) const {
    const double x = parameters[0][0];
    residuals[0] = 10 - x;

    // Compute the Jacobian if asked for.
    if (jacobians != NULL && jacobians[0] != NULL) {
      jacobians[0][0] = -1;
    }
    return true;
  }
};

CeresはEvaluate関数で定義された内容を元に、

残差項とヤコビ行列を設定します。

 

しかし、上記のヤコビ行列の設定は、

コードの可読性を落とすので、

ceresの開発者としては、

最適化のパフォーマンスに問題が無い場合は、

解析的に導関数を計算出来る時であっても、

AutoDiffCostFunctionを使うことを推奨しているようです。

その他の導関数の求め方

導関数の求め方は、現在でも様々な方法が提案されており、

それらの幾つかの手法はceresにも実装されています。

  • DynamicAutoDiffCostFunction

  • CostFunctionToFunctor

  • NumericDiffFunctor

  • ConditionedCostFunction

詳しい方法はAPIのリファレンスを確認下さい。

 

複雑な最適化の実施

少し複雑な最適化をCeresソルバで解いてみたいと思います。

今回は、パウエル関数の最適化を考えます。

 

パウエル関数は4つの入力を持ち、

f:id:meison_amsl:20160924115040p:plain

下記の式で4次元のベクトルが計算される関数です。

f:id:meison_amsl:20160924115132p:plain

ここでこのパウエル関数の二乗ノルムを最小化する

xをCeresで計算するコードを作ってみたいと思います。

 

ちなみに、パウエル関数は

x=[0,0,0,0]でF(x)=0が極小値になることが

わかっています。

 

まず初めに、

それぞれのベクトルの要素を計算する

評価関数のコードを作成します。

 

struct F1{
  template <typename T>
  bool operator()(
      const T* const x1,  
      const T* const x2,
      T* residual) const{
    residual[0]=x1[0]+T(10.0)*x2[0];
    return true;
    }
};

struct F2{
  template <typename T>
  bool operator()(
      const T* const x3,  
      const T* const x4,
      T* residual) const{
    residual[0]=T(sqrt(5.0))*(x3[0]-x4[0]);
    return true;
    }
};

struct F3{
  template <typename T>
  bool operator()(
      const T* const x2,  
      const T* const x3,
      T* residual) const{
    residual[0]=(x2[0]-T(2.0)*x3[0])*(x2[0]-T(2.0)*x3[0]);
    return true;
    }
};

struct F4{
  template <typename T>
  bool operator()(
      const T* const x1,  
      const T* const x4,
      T* residual) const{
    residual[0]=T(sqrt(10.0))*(x1[0]-x4[0])*(x1[0]-x4[0]);
    return true;
    }
};

テンプレートTを使って評価関数を計算しているので、

根号で計算されたdouble型をTにキャストするのを忘れないことと、

入力の変数は配列のポインタの形なので0番目の要素を指定することを

忘れないようにしましょう。

 

続いて、最適化を実施する関数を書きます。

int main(int argc, char** argv){
  google::InitGoogleLogging(argv[0]);

  //最適化問題を解く変数と初期値の設定
  double x1= 3.0;
  double x2=-1.0;
  double x3= 0.0;
  double x4= 1.0;


  //最適化問題を解く用のオブジェクトの生成
  Problem problem;

  //最適化問題に残差項と変数を設定
  problem.AddResidualBlock(
      new AutoDiffCostFunction<F1,1,1,1>(new F1),NULL,&x1,&x2);
  problem.AddResidualBlock(
      new AutoDiffCostFunction<F2,1,1,1>(new F2),NULL,&x3,&x4);
  problem.AddResidualBlock(
      new AutoDiffCostFunction<F3,1,1,1>(new F3),NULL,&x2,&x3);
  problem.AddResidualBlock(
      new AutoDiffCostFunction<F4,1,1,1>(new F4),NULL,&x1,&x4);


  //最適化の実行
  Solver::Options options;//最適化のオプション設定用構造体
  options.linear_solver_type=ceres::DENSE_QR;
  options.minimizer_progress_to_stdout=true;//最適化の結果を標準出力に表示する。
  Solver::Summary summary;//最適化の結果を格納するよう構造体
  Solve(options,&problem,&summary);//最適化の実行

  //結果の表示
  std::cout<<summary.BriefReport()<<std::endl;
  std::cout<<"x1:"<<x1<<std::endl;
  std::cout<<"x2:"<<x2<<std::endl;
  std::cout<<"x3:"<<x3<<std::endl;
  std::cout<<"x4:"<<x4<<std::endl;

  return 0;
}

コスト関数が複数の場合は、

AddResidualBlockを複数回呼んで、

それぞれ追加すればOKです。

 

ちなみに

AutoDiffCostFunction<F1,1,1,1>

のそれぞれの入力は、

1つめ:評価関数の構造体

2つめ:residualの次元数(F1のresidualは次元が1)

3つめ:1つめの次元数(F1のx1は次元が1)

4つめ:1つめの次元数(F1のx2は次元が1)

です。

詳しくは下記を参照下さい。

 

すると下記のような結果が得られます。

f:id:meison_amsl:20160924124746p:plain

ちゃんとx=[0,0,0,0]でF(x)=0の

極小値が探索できていることがわかります。

 

こちらのコードは下記のGithubリポジトリでも公開しています。

github.com

より深く最適化を学びたい人は

下記の資料が体系的に学べるのでおすすめです。

これなら分かる最適化数学―基礎原理から計算手法まで

これなら分かる最適化数学―基礎原理から計算手法まで

工学基礎 最適化とその応用 (新・工科系の数学)

工学基礎 最適化とその応用 (新・工科系の数学)

はじめての最適化

はじめての最適化

非線形最適化の基礎

非線形最適化の基礎

  

参考資料

myenigma.hatenablog.com

myenigma.hatenablog.com

myenigma.hatenablog.com

myenigma.hatenablog.com

myenigma.hatenablog.com

myenigma.hatenablog.com

myenigma.hatenablog.com

myenigma.hatenablog.com

はじめての最適化

はじめての最適化