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

MyEnigma

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

GoogleのC++最適化ライブラリCeres Solverを使った曲線フィッティングの方法

optimization

最適化の数学 (共立講座 21世紀の数学 13)

最適化の数学 (共立講座 21世紀の数学 13)

目次

はじめに

今回の記事もこれらの記事の続きです。

myenigma.hatenablog.com

myenigma.hatenablog.com

 

今回は、Ceres Solverを使って

曲線フィッティングをしてみようと思います。

 

Ceres Solverによる曲線フィッティング

先日までの記事で扱っていたサンプルでは、

データはなく、関数の最小値などを計算していました。

しかし、元々最小二乗法や非線形最適化が開発された目的は、

データの曲線当てはめでした。

今回は、このデータの曲線当てはめのサンプルコードを作成してみたいと思います。 

 

今回は下記の指数関数に、

f:id:meison_amsl:20160925124328p:plain

m=0,σ=0.2のガウス分布を付加したデータを

準備します。

 

今回は下記のコードを使って、

#!/usr/bin/env python
# -*- coding: utf-8 -*-
# author:Atsushi Sakai

import matplotlib.pyplot as plt
import numpy as np

x=np.arange(0,10,0.1)
y=[np.exp(0.3*tx+0.1)+np.random.normal(0.0,0.2) for tx in x]

plt.plot(x,y,'+r')
plt.grid(True)
plt.show()

f=open('data.csv', 'w')
for (x,y) in zip(x,y):
    f.write(str(x)+","+str(y)+"\n")
f.close()

下記のグラフのような、

データを生成しました。

f:id:meison_amsl:20160925134432p:plain

 

続いて、このデータファイルを読み込ませて、

Ceresを使って曲線当てはめをしてみます。

 

まず初めにコスト関数を設定します。

それぞれのデータに対してコストを計算し、

その総和をコストとすると、下記のようになります。

/**
 * @brief コスト関数
 **/
struct ExponentialResidual{
  ExponentialResidual(double x, double y)
    :x_(x), y_(y){}

  template <typename T>
  bool operator()(const T* const m, const T* const c, T* residual) const {
    residual[0]=T(y_)-exp(m[0]*T(x_)+c[0]);
    return true;
  }

  private:
    //Observations for a sample
    const double x_;
    const double y_;
};

データを使って、コスト関数を計算する場合は、

データを格納するメンバ変数を作り、

コンストラクタでその値を設定するような

コスト関数構造体を作る必要があります。

 

続いて、先ほどのコスト関数構造体を使って、

最適化の設定をします。

データのx,yのデータがそれぞれ格納された配列があるとして、

そのデータをそれぞれループで取得し、

先程のコスト関数に入力したオブジェクトを

AutoDiffCostFunctionでコスト関数オブジェクトに追加していく形です。

//最適化問題を解く変数と初期値の設定
double m = 0.0;//exp(mx+c)
double c = 0.0;

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

for(int i=0; i<x.size(); i++){
   CostFunction* cost_function=new AutoDiffCostFunction<ExponentialResidual, 1, 1, 1>(
          new ExponentialResidual(x[i],y[i]));

   problem.AddResidualBlock(cost_function, NULL, &m, &c);
}

 

あとは以前の記事で紹介した最適化の方法と同じように、

最適化を実施すると下記のような結果が得られます。

f:id:meison_amsl:20161004135829p:plain

上記の結果を見るとわかるように、

m=0.0, c=0.0から初めた最適化が、

m=0.300439, c=0.0966003と、

正解であるm=0.3,c=0.1にかなり近い値で計算出来ていることがわかります。

 

上記の最適化のコードはこちらで公開しています。

github.com

 

Ceres Solverによるロバスト曲線フィッティング

前述の方法は、

入力データがある程度正確であれば問題ありませんが、

もしデータの中にOutlierが含まれる場合(ノイズモデルに従わないデータ)、

下記の図のように、最適化された結果がOutlierに引きずられてしまいます。

f:id:meison_amsl:20161004142832p:plain

 

そんな時に利用すると良いのが、損失関数です。

損失関数をうまく設定することで、outlierに強い最適化を実施することができます。

CeresSolverの場合、

AddResidualBlockの2つ目の引数を、

problem.AddResidualBlock(cost_function, NULL , &m, &c);

下記のように変更することにより、

problem.AddResidualBlock(cost_function, new CauchyLoss(0.5) , &m, &c);

損失関数を設定することができます。

CauchyLossはCeresソルバの損失関数の一つであり、

0.5は損失関数のスケールパラメータです。

ちなみにCauchyLossを使う場合は、

using ceres::CauchyLossを宣言するのを忘れないようにしましょう。

 

上記のように、損失関数を設定すると、

下記のようにOutlierが含まれるデータであっても、

真の値に近い、最適化を実施することができます。

f:id:meison_amsl:20161004145806p:plain

 

上記の最適化のコードは下記で公開しています。

github.com

 

Ceres Solver関連記事

myenigma.hatenablog.com

myenigma.hatenablog.com

myenigma.hatenablog.com

myenigma.hatenablog.com

 

参考資料

myenigma.hatenablog.com

myenigma.hatenablog.com

myenigma.hatenablog.com

myenigma.hatenablog.com

最適化の数学 (共立講座 21世紀の数学 13)

最適化の数学 (共立講座 21世紀の数学 13)