kivantium活動日記

プログラムを使っていろいろやります

レーベンバーグ・マーカート法(Levenberg-Marquardt Method)を使うには

レーベンバーグ・マーカート法とは

レーベンバーグ・マーカート法は、非線形関数の自乗の和の形で表された関数の極小を求める反復法の一種です。
最急降下法ニュートン法を組み合わせた方法とされ、非線形最小二乗問題(非線形な関数の二乗和の最小値を求める問題)を解く標準的な方法になっています。現在の解が正解から遠い場合は遅いが収束することが保証されている最急降下法と同じように動作し、正解から近い場合はニュートン法を実行するものです。

使うライブラリ

ここでは、レーベンバーグ・マーカート法のANSI C実装であるlevmarを使います。ライセンスはGPLで、公式サイトはここです。

インストール

筆者の環境(Ubuntu 13.10)でのインストール方法を紹介します。
まず、LEVMARの動作に必要なライブラリをインストールします

sudo apt-get install liblapack-dev liblas-dev

次に、LEVMARを配布しているPPAを利用してインストールします。(PPA名はUbuntuのバージョンによって異なるようです)

sudo add-apt-repository ppa:olebole/astro-saucy
sudo apt-get update
sudo apt-get install levmar

使用例

3点を通る1次方程式を求める

3点A(-3.1, 0.2), B(0.3, 0.8), C(1.2, 1.0)を通る直線を求めます。コードは次の通りです

#include <stdio.h>
#include <stdlib.h>
#include <levmar.h>

void func(double *p, double *hx, int m, int n, void *adata);

int main()
{
       double* p;
	   p = (double*)malloc(2*sizeof(double));
	   p[0]=0.0; p[1]=0.0;       //最初の推定値
       double* x = NULL;
       int m = 2;                //求めたい変数の数
       int n = 3;                //方程式の数
       int itmax = 300;          //反復法を行う回数
       double* opts = NULL;      //デフォルト値を使用
       double* info = NULL;      //最小化に関するパラメーター、ここでは使わない
       double* work = NULL;      //メモリを自動確保
       double* covar = NULL;     //必要ないのでNULL
       void* adata = NULL;       //追加のデータを入れる配列。これも使わない

       dlevmar_dif(func, p, x, m, n, itmax, opts, info, work, covar, adata); //pを求める

       printf("p1 = %f\np2 = %f\n", p[0], p[1]);
  
       return 0;
}


void func(double *p, double *hx, int m, int n, void *adata)     //p:parameter vector, hx:objective function vector
{
       double  datax[3] = {-3.1, 0.3, 1.2};     //x座標を集めたもの
       double  datay[3] = {0.18, 0.86, 1.04};   //y座標を集めたもの	

       hx[0] = (-datay[0] + p[0] * datax[0] + p[1]) * (-datay[0] + p[0] * datax[0] + p[1]); //(-y_1+p_0*x_1+p_1)^2を計算
       hx[1] = (-datay[1] + p[0] * datax[1] + p[1]) * (-datay[1] + p[0] * datax[1] + p[1]); 
       hx[2] = (-datay[2] + p[0] * datax[2] + p[1]) * (-datay[2] + p[0] * datax[2] + p[1]); 
       return;
}

main.cppに保存して

g++ main.cpp -llevmar -llapack -lblas

コンパイルし、./a.outで実行すると

p1 = 0.188845
p2 = 0.806997

という出力を得ます。

求めたい直線は
{y=p_{1}x+p_{2}}という形をしています。
求めたいp1, p2は、xにx1〜x3を代入したときにyの値とy1〜y3の誤差が最小になるような値となります。
すなわち、{(p_{1}x_{1}+p_{2}-y_{1})^2+(p_{1}x_{2}+p_{2}-y_{2})^2+(p_{1}x_{3}+p_{2}-y_{3})^2}が最小になるようなp1, p2を求めれば、求めたい直線が分かることになります。
これを実現するコードが上のものです。実際に出力された値を入れてグラフを書いたものが下図になります。

全ての点を直線上に乗せることはできませんが、それっぽい直線が求まったことが分かると思います。
同じような方法で二次関数なども求めることができるはずです。

非線形関数の自乗の和の極小値を求める

ここでは{\sin(x)^2+\cos(\frac{x}{2})^2}の極小値を求めます。
コードは以下の通りです。

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <levmar.h>

void func(double *p, double *hx, int m, int n, void *adata);

int main()
{
       double* p;
       p = (double*)malloc(2*sizeof(double));
       p[0]=2.5;                 //最初の推定値
       double* x = NULL;
       int m = 1;                //求めたい変数の数
       int n = 1;                //方程式の数
       int itmax = 1000;         //反復法を行う回数
       double* opts = NULL;      //デフォルト値を使用
       double* info = NULL;      //最小化に関するパラメーター、ここでは使わない
       double* work = NULL;      //メモリを自動確保
       double* covar = NULL;     //必要ないのでNULL
       void* adata = NULL;       //追加のデータを入れる配列。これも使わない

       dlevmar_dif(func, p, x, m, n, itmax, opts, info, work, covar, adata); //pを求める

       printf("x = %f\n", p[0]);
  
       return 0;
}

void func(double *p, double *hx, int m, int n, void *adata)     //p:parameter vector, hx:objective function vector
{
       hx[0] = sin(p[0]) * sin(p[0]);
       hx[1] = cos(p[0]/2) * cos(p[0]/2);
       return;
}

実行結果は

x = 3.141448

となります。ここで、最初の推定値を0.5にした結果は次のようになります。(p[0]=2.5の行をp[0]=0.5に書き換える)

x = 0.000050

最初の推定値によって異なる結果が出ました。この原因を知るためにy={\sin(x)^2+\cos(\frac{x}{2})^2}のグラフを見てみましょう。

グラフからも分かるように、x=0.5に近い極小値はx=0の時で、x=2.5に近い時は極小値がx=3.14の時になります。
レーベンバーグ・マーカート法で求められるのはあくまで極小値なので、得られる値が真の最小値ではない可能性があることに気をつけなければいけません。