minpack.md

minpack で関数フィッティング

Levenberg-Marquardt法 (wiki様) による最小二乗法ソルバ. 今回は fortran で実装してみた.

インストール

準備 : fpmのインストール

インストールには Fortran Package Manager (fpm) を使うので、まずはfpmをインストール. 公式HP (リンク)でバイナリが配布されているので、fpm-0.11.0-linux-x86_64-gcc-12 をダウンロードしてパスの通ったところに置く. ついでに名前も fpm に変えておく.

sudo cp fpm-0.11.0-linux-x86_64-gcc-12 /usr/local/bin/fpm
fpm --version

バージョンが出てくればOK.

minpack のインストール

Github (リンク) からソースコードをダウンロードしてコンパイル.

git clone https://github.com/fortran-lang/minpack.git
cd minpack
fpm build
fpm test
fpm install

＊＊＊という場所に .a と .mod が生成されたと表示されればOK.

自分の場合 ~/.local/lib/ に libminpack.a が、~/.local/include/ に3つの .mod ファイルが生成されたので、/usr/localに移動しておく (任意).

sudo cp ~/.local/lib/libminpack.a /usr/local/lib/
sudo cp ~/.local/include/*.mod /usr/local/include/

完了!

lmdif1 による関数フィッティング

ガウス関数 : $$y = a \exp \left[ - (\frac{x-b}{c})^2 \right]$$ でフィッティングしてみる. 与えるデータは $$x=0-10$$ の範囲を500等分したグリッド上で、$$(a,b,c) = (5.25, 3.14, 1.41)$$ としてあらかじめ用意しておいた (fort.55という名前).

program gauss_fit
use minpack_module, only: wp, enorm, lmdif1
implicit none

integer, parameter :: n = 3, m = 500, lwa = m*n+5*n+m
real(8) :: tol, x(n), fvec(m), wa(lwa)
integer :: iwa(n), info

real(8) :: xi(M), yi(M)
integer :: i

call rd_data

x = [1.d0, 1.d0, 1.d0]  ! initial parameter
tol = 1d-8

call lmdif1(fcn, m, n, x, fvec, tol, info, iwa, wa, lwa)
print*, enorm(m,fvec), info, x

contains
subroutine fcn(m, n, x, fvec, iflag)

    integer, intent(in) :: m
    integer, intent(in) :: n
    real(8), intent(in) :: x(n)
    real(8), intent(out) :: fvec(m)

    do i = 1, M
        fvec(i) = yi(i) - ( x(1)*exp( -((xi(i)-x(2))/x(3))**2 ) )
    end do

end subroutine fcn
subroutine rd_data

    open(101, file="fort.55", status="old")
        do i = 1, M
            read(101,*) xi(i), yi(i)
        end do
    close(101)

end subroutine rd_data

end program gauss_fit

コンパイルと実行

gfortran gauss_fit.f90 -I/usr/local/include -L/usr/local/lib -lminpack
./a.out
   0.0000000000000000                2   5.2500000000000000        3.1400000000000001        1.4099999999999999

左から残差二乗誤差、フィッティングの収束INFO、3つのパラメータが並んでいる.

(a,b,c)を見ると、あらかじめ決めておいたパラメータが求まったことがわかる.

INFOは公式HPを見ると次の通り:

    INFO = 0  Improper input parameters.

    INFO = 1  Both actual and predicted relative reductions in the
              sum of squares are at most FTOL.

    INFO = 2  Relative error between two consecutive iterates is
              at most XTOL.

    INFO = 3  Conditions for INFO = 1 and INFO = 2 both hold.

    INFO = 4  The cosine of the angle between FVEC and any column
              of the Jacobian is at most GTOL in absolute value.

    INFO = 5  Number of calls to FCN has reached or exceeded
              MAXFEV.

    INFO = 6  FTOL is too small.  No further reduction in the sum
              of squares is possible.

    INFO = 7  XTOL is too small.  No further improvement in the
              approximate solution X is possible.

    INFO = 8  GTOL is too small.  FVEC is orthogonal to the
              columns of the Jacobian to machine precision.

   (翻訳)
   INFO = 0 入力パラメータが不適切です。

   INFO = 1 二乗和の実際の相対的減少と予測された相対的減少はどちらも最大で FTOL です。

   INFO = 2 2 つの連続する反復間の相対誤差は最大で XTOL です。

   INFO = 3 INFO = 1 と INFO = 2 の両方の条件が成立します。

   INFO = 4 FVEC とヤコビアンの任意の列の間の角度のコサインの絶対値は最大で GTOL です。

   INFO = 5 FCN の呼び出し回数が MAXFEV に達したか、それを上回りました。

   INFO = 6 FTOL が小さすぎます。二乗和をこれ以上減らすことはできません。

   INFO = 7 XTOL が小さすぎます。近似解 X をこれ以上改善することはできません。

   INFO = 8 GTOL が小さすぎます。FVEC はヤコビアンの列に対して機械精度に対して直交しています。

後は、例えば許容誤差 tol を変えてみたり、パラメータの初期値 (この手の方法では初期値から大きく外れるとNG)をいじるとか、いろいろと微調整は必要そう.

数値計算・プラズマの備忘録

2025年3月25日火曜日