.-- --

スポンサーサイト

上記の広告は1ヶ月以上更新のないブログに表示されています。
新しい記事を書く事で広告が消せます。
スポンサー広告 comment(-) trackback(-)
.17 2013

Eigenを使って非線形最小2乗法 ~Levenberg-Marquardt法

こんにちわ。Uedaです。

今回は非線形最小二乗法の問題でよく使われるLevenberg-Marquardt法(以下LM法)をEigenを使って解いてみます。
(LM法の理解は筆者は間違っているかもしれません。以下の語尾には「だと思われます」が付きます

参考サイト
http://d.hatena.ne.jp/wildpie/20120712

非線形問題の解き方についてはニュートン法やガウス・ニュートン法があります。
ニュートン法 ~ wiki
http://ja.wikipedia.org/wiki/ニュートン法

ガウス・ニュートン法参考サイト
http://homepage3.nifty.com/skomo/f6/hp6_9.htm

これらの方法は解のにおいのする方向へ探索していき、一歩ずつ解に近づいていくという方法です。(この解のにおいのする方向を調べるのに微分使ってます。

これらの方法は解への収束速度はとても速いのですが、(初期値の取り方を誤ったりすると)誤った解に収束してしまう恐れの大きい方法です。Levenberg-Marquardt法は少し速度を犠牲にする代わりに、この危険性を下げる方法です。
LM法も解のにおいのする方向へ探索していき、一歩ずつ解に近づいていくという方法ですが、その歩幅を工夫しています。LM法の歩幅は最初は小さく、解の方向でありそうだとわかれば次第に歩幅を大きくしていくという具合です。(ガウス・ニュートン法は最初から大股な感じです。

参考サイトをもとにEigenのLM法を試してみました。
<問題> z = a / x + b / y + c の非線形の式のa, b, cを、ノイズを乗せたM個のデータxi, yi, zi(0 < i < M)から推定する。
・正解a, b, cはそれぞれ
a = 1.0 ,b = 2.0, c = 3.0
・初期データ_a, _b, _cはそれぞれ
_a = a + noise2, _b = b + noise2, _b = b + noise
・-0.5 < xi < 0.5, -0.5 < yi < 0.5, zi = z + noise1
・-2.5 < noise1 < 2.5
・-50.0 < noise2 < 50.0

コード



#import "Eigen/Core"
#import "Eigen/Dense"
#import "unsupported/Eigen/NonLinearOptimization"
#import "unsupported/Eigen/NumericalDiff"

・・・

const int DATA_COUNT = 100;
const double ANS_A = 1.0;
const double ANS_B = 2.0;
const double ANS_C = 3.0;
const double NOIZE1 = 5.0;
const double NOIZE2 = 100.0;

double *xa = new double[DATA_COUNT];
double *ya = new double[DATA_COUNT];
double *za = new double[DATA_COUNT];
for (int i = 0; i < DATA_COUNT; i++)
{
double x = 0.5 - [self random];
double y = 0.5 - [self random];
double noize = (0.5 - [self random]) * NOIZE1;
double z = ANS_A/x + ANS_B/y + ANS_C + noize;
xa[i] = x;
ya[i] = y;
za[i] = z;
}

VectorXd p(3); // beta1とbeta2の初期値(適当)
double a = ANS_A + (0.5 - [self random]) * NOIZE2;
double b = ANS_B + (0.5 - [self random]) * NOIZE2;
double c = ANS_C + (0.5 - [self random]) * NOIZE2;
p(0) = a;
p(1) = b;
p(2) = c;

misra1a_functor functor(3, DATA_COUNT, xa, ya, za);

NumericalDiff<misra1a_functor> numDiff(functor);
LevenbergMarquardt<NumericalDiff<misra1a_functor>> lm(numDiff); // 注意 "<"が全角
int info = lm.minimize(p);

NSLog(@"正解 %f, %f, %f", ANS_A, ANS_B, ANS_C);
NSLog(@"推定前 %f, %f, %f", a, b, c);
NSLog(@"推定後 %f, %f, %f, %d", p[0], p[1], p[2], info);
}

template
struct Functor
{
typedef _Scalar Scalar;
enum {
InputsAtCompileTime = NX,
ValuesAtCompileTime = NY
};
typedef Matrix InputType;
typedef Matrix ValueType;
typedef Matrix JacobianType;
};

struct misra1a_functor : Functor
{
misra1a_functor(int inputs, int values, double *x, double *y, double *z)
: inputs_(inputs), values_(values), x(x), y(y), z(z) {}

double *x;
double *y;
double *z;
int operator()(const VectorXd& b, VectorXd& fvec) const
{
for (int i = 0; i < values_; ++i) {
fvec[i] = z[i] - b[0]/x[i] - b[1]/y[i] - b[2];
}
return 0;
}

const int inputs_;
const int values_;
int inputs() const { return inputs_; }
int values() const { return values_; }
};



結果

2013-05-17 17:28:55.085 TestEigen[2043:c07] 正解 1.000000, 2.000000, 3.000000
2013-05-17 17:28:55.086 TestEigen[2043:c07] 推定前 44.168610, -46.338116, 0.831853
2013-05-17 17:28:55.086 TestEigen[2043:c07] 推定後 0.996553, 1.992320, 2.897285, 1

むちゃくちゃな初期値(結果の推定前)から、見事、正解に近づくことができました!
次回コードレビューします。
関連記事
スポンサーサイト

Comment

Post comment

  • comment
  • secret
  • 管理者にだけ表示を許可する

Trackback

trackbackURL:http://appteam.blog114.fc2.com/tb.php/263-8c1e4587

ブログ内検索

関連リンク

製品情報

最新記事

カテゴリ

プロフィール

neoxneo



NEXT-SYSTEM iOS Developers Blog


  • UTO:
    カナダ版iPhone4Sは、マナーモードでシャッター音がならない…


  • Ehara:
    ...


  • Hayate:
    ...


  • Tasaki:
    Developer登録完了...したのはいいけど


  • Ueda:
    ...



リンク

上記広告は1ヶ月以上更新のないブログに表示されています。新しい記事を書くことで広告を消せます。