作業中のメモ

よく「計算機」を使って作業をする.知らなかったことを中心にまとめるつもり.

レーベンバーグ・マーカート法と準ニュートン法(C 言語)

どうも筆者です.

最近,観測データを既知の関数で近似する必要が生じた.そのための近似方法を以前示したが,あれは,そのままプログラムを組むと,次元が大きくなると解きにくくなってしまうという問題がある. そのため,何か良い方法がないかと調べていたところ,「レーベンバーグ・マーカート法」というものがあることが分かった.

レーベンバーグ・マーカート法は,非線形関数の 2 乗和が最小になるものを求める方法として使われるらしい.これを使えば,誤差の 2 乗和が最小となるという点で,関数フィッティングが行える.

これを使うには,「LAPACK」と「BLAS」が必要になる.なので,参考サイトにあるように,インストールを試みたが,「levmar」が上手くインストール出来なかった.しかし,筆者の PC には,Intel コンパイラがインストールされている.このときに,MKL も一緒にインストールされるので,これを用いれば実行出来るのではないかと考えた.Makefile があったので,Intel 用の「Makefile.icc」を書き換えることにした.現段階では,バージョンが 2.6 であり,Makefile.icc の中身は,1 段目のようになっている.このまま使えば動くはずだったが,上手く動かなかった(f2c をインストールし忘れていたのが原因だと思う)ので,Makefile.icc を 2 段目のように書き換えた.

これで,実行することが出来た.ここで,[IntelVersion] は,各々がインストールした Intel Compiler のバージョンが記述されたディレクトリ名が入る.ここは,各自調べて欲しい.

また,関数の最小値も求めたかったので,色々調べた.すると,「liblbfgs」というパッケージが存在することが分かった.これは,Fortran で実装されていた「準ニュートン法」を C 言語に直したものらしい.よくやったなぁ...作者に感謝しつつ,手順に従ってインストールすることにした.しかし,これも何故かインストールが上手くいかなかったので,仕方なく,毎回コンパイルすることにした.

修正した Makefile とソースコートは以下のようになった.

また,ファイルの階層は,

  • Makefile [libsrc] [src]
    • [libsrc] 「levmar-2.6」に含まれているファイルを全てコピーし,Makefile を書き換えたもの.
    • [src] MT.h arithmetic_ansi.h lbfgs.c lbfgs.h main.c

という風になっている.分かりづらくて申し訳ない.

今回は,

 { \displaystyle 
  y = -2x^{2} + 6x + \dfrac{1}{2}
}

の最大値を求めることにした.真の関数値に ボックス・ミューラー法により,標準正規分布に従う乱数を加えたものを測定値とした.そして,

  1. レーベンバーグ・マーカート法により,測定データを 2 次関数で近似する.そのときのパラメータを p[0],p[1],p[2] とした.
  2. 近似したパラメータを用いて,準ニュートン法により,最大値をとる  x を求めた.

という流れで処理を行った.最初は,levmar で実装されている「dlevmar_dif」の引数「adata」の使い方が分からなかったが,調べたところ,キャストして渡せばよいらしいので,プログラムにあるように渡した.また,これと似た関数に「dlevmar_def」というものがある.これは,ヤコビ行列を計算する関数を自分で与える場合に利用する.今回は,面倒だったので,近似したものを用いた(すなわち,dlevmar_dif を利用した).

また,lbfgs では,float と double の使い分けが出来るようにしてあるため,double ではなく,「lbfgsfloatval_t」というものが利用されている.個人的には,double でしか計算しないので,ここは,double と書いてもかまわない(ヘッダーで typedef により宣言されているだけなので).

dlevmar_dif には,ユーザー側が与えるデータとして,測定データを与えているのに対し,lbfgs では,推定した係数データを与えていることに注意.また,最小値を求める関数であるため,最大値を求める場合は,マイナスをかけておく必要がある.そのため,evaluate 関数では,最後に,g[0] と fx の符号を反転させている.

これを実行すると,p[0] = -2.002300,p[1] = 5.944318,p[2] = 0.611619 と係数を推定でき,最大値を取る  x は, \hat{x} = 1.484373 と求められる.真の関数では, x^{*} = \dfrac{3}{2} のとき最大となるので,正しく推定できていることが確認できた.

参考文献

Levenberg-Marquardt in C/C++

http://kivantium.hateblo.jp/entry/20140408/p1

http://www.sat.t.u-tokyo.ac.jp/~omi/random_variables_generation.html#Prepare_MT