MKLを使った非線形関数の最小化 in Fortran

Intel Math Kernel Library(MKL)と言えば、行列の積は固有値問題などを解くBLASとLAPACKが有名だと思います。
このMKLには、非線形最小二乗問題を解くためのルーチンがありますので、それを紹介します。
アルゴリズムはTrust-Region algorithmです。

バージョン

• gfortran: gcc version 8.3.0 (Homebrew GCC 8.3.0)
• MKL: 2018.1.126

非線形最小二乗問題

非線形最小二乗問題(Nonlinear Least Squares Problem)とは、

\min_{x} ||F(x)||^2 = \min_{x} ||y -f(x)||^2

です。ここで、yはm次元の実数のベクトル、xはn次元の実数のベクトルとなります。m>nとします。
上の式は、m種類の関数F_i(x)にベクトルxを入れた時に、どのF_i(x)も最小となるようなxを求める問題となります。

例えば、

F_i(x) = \Bigl{|}\frac{1}{10i }- \left[ 10i - \sum_{k=1}^{n} \frac{x_k}{10i - x_k} \right]^{-1} \Bigl{|}

みたいな関数が最小となるようなxを求めるというわけです。

MKL

MKLには非線形最小二乗問題のサブルーチンとしてstrnlsp_solveとdtrnlsp_solveがあります。sは単精度実数、dは倍精度実数です。
このコードの説明は、
http://scc.qibebt.cas.cn/docs/compiler/intel/2011/mkl/mkl_manual/osr/osr_NLSPnoC.htm
にあります。
サンプルのコードがMKLをインストールしたディレクトリの
examples\solver\source
にあるので、そちらを参考にして書いてみました。
注意点としては、mkl_rci.fiというインクルードファイルが必要なので、
コンパイルは例えば、

gfortran test.f90 -lmkl_intel_lp64  -lmkl_sequential -lmkl_core  -lpthread -lm -ldl -I/opt/intel/compilers_and_libraries_2018.1.126/mac/mkl/include

のようになります。

コード

Fortran90のコードは以下の通りです。

subroutine test()
    implicit none
    include "mkl_rci.fi"
    external testfunc2

    integer,parameter::N=40
    integer,parameter::M=40
    real(8)::X (N)
    real(8)::fvec (M)
    real(8)::fjac (M,N)
    integer(8) :: handle
    integer::res
    real(8)::eps(6)
    real(8)::jac_eps
    integer::iter1,iter2
    real(8)::rs
    integer::i,j
    integer::RCI_Request
    integer::successful
    !** NUMBER OF STOP-CRITERION
    integer::st_cr
    real(8)::r1,r2
    integer::iter
    integer::icount

    !** PRECISIONS FOR JACOBI CALCULATION
    jac_eps = 1.d-8
    !** SET MAXIMUM NUMBER OF ITERATIONS
    iter1 = 1000
!** SET MAXIMUM NUMBER OF ITERATIONS OF CALCULATION OF TRIAL-STEP
    iter2 = 100
    do i = 1, 6
        eps (i) = 1.d-8
    end do
    !** SET INITIAL STEP BOUND
    rs = 100.d0

    !** SET THE INITIAL GUESS
    call random_number(x)

!** SET INITIAL VALUES
    do i = 1, M
        fvec (i) = 0.d0
        do j = 1, N
            fjac (i, j) = 0.d0
        end do
    end do

    res = dtrnlsp_init(handle, N, M, x, eps, iter1, iter2, rs)


    icount = 0
    RCI_Request = 0
    successful  = 0
    do while (successful == 0)        
        icount = icount + 1

        res = dtrnlsp_solve(handle, fvec, fjac, RCI_Request)
        !write(*,*) icount,RCI_Request
        select case (RCI_Request)
        case (-1, -2, -3, -4, -5, -6)
            !**   STOP RCI CYCLE
            successful = 1
        case (1)
!**   RECALCULATE FUNCTION VALUE
!**     M               IN:     DIMENSION OF FUNCTION VALUE
!**     N               IN:     NUMBER OF FUNCTION VARIABLES
!**     X               IN:     SOLUTION VECTOR
!**     FVEC            OUT:    FUNCTION VALUE F(X)
            call testfunc2(M, N, X, fvec)
        case (2)
!**   COMPUTE JACOBI MATRIX
!**     EXTENDED_POWELL IN:     EXTERNAL OBJECTIVE FUNCTION
!**     N               IN:     NUMBER OF FUNCTION VARIABLES
!**     M               IN:     DIMENSION OF FUNCTION VALUE
!**     FJAC            OUT:    JACOBI MATRIX
!**     X               IN:     SOLUTION VECTOR
!**     JAC_EPS         IN:     JACOBI CALCULATION PRECISION
            !res = djacobi(fcn, n, m, fjac, x, jac_eps)
            res = djacobi(testfunc2, N, M, fjac, X, jac_eps)
        end select
    end do

!** GET SOLUTION STATUSES
!**   HANDLE            IN: TR SOLVER HANDLE
!**   ITER              OUT: NUMBER OF ITERATIONS
!**   ST_CR             OUT: NUMBER OF STOP CRITERION
!**   R1                OUT: INITIAL RESIDUALS
!**   R2                OUT: FINAL RESIDUALS
    res = dtrnlsp_get(handle, iter, st_cr, r1, r2)
    write(*,*) "iter = ",iter
    write(*,*) "x = ",x
    call testfunc(M, N, x, fvec)
    write(*,*) "fvec = ",fvec
    write(*,*) "initial res = ",r1,"final res = ",r2
    res = dtrnlsp_delete(handle)


end subroutine


    subroutine testfunc2(M, N, X, F) !
        implicit none
        integer,intent(in):: M, N
        real(8),intent(in)::X (N)
        real(8),intent(out)::F (M)
        integer i,k

        do i = 1, N
            F(i) = i*10 
            do k=1,M

                F(i) = F(i)- x(k)/(i*10 - x(k))
            end do
            F(i) = abs(1d0/(i*10)-1d0/F(i))
        end do

    end subroutine testfunc2


program main
    implicit none
    call test()
end program

まず初期化をするためにdtrnlsp_initを呼び、その後dtrnlsp_solveをループの中で呼びます。
その後ヤコビアン行列を計算したりサブルーチンを呼んだりします。
djacobiというのはある関数のヤコビアン行列を計算する関数ですが、これもMKLに入っています。

例では、xの次元が40、yの次元が40です。xは乱数で初期化しました。

出力結果は

 iter =           19
 x =   -2.3692537834812321E-003  -9.7268400230563646E-003   3.6446073380145961E-003   1.8391493884423330E-003  -1.0072348936971448E-002   2.1956223501030927E-003   9.0863286208358041E-003  -5.7405536248615061E-003   1.5273677639816479E-003  -1.2091518617019284E-002  -6.0245493944464353E-003   9.0758965850220360E-003   3.7251262242309164E-003  -5.1543024403054912E-003  -3.0718820339978559E-003   7.2074732937779089E-003  -1.2071924622387297E-003   3.2189402835896023E-003   5.9119356813252880E-003  -4.8407902190080706E-003  -5.5791253155157406E-004   6.2848742255312362E-004  -3.9097023300285845E-003   4.1594757345492991E-003   8.2795445826901893E-003  -7.3663899301664847E-003  -3.3012662977024991E-003   2.5150282805543323E-003   7.0307552593837424E-003  -1.8395506788844776E-003   8.9668254124308042E-003  -6.0515507635582151E-003   2.3738226737637994E-003   2.5224512540703726E-003   3.1335234997537198E-003  -3.6931877135098906E-003  -3.8847665528000822E-003  -1.6007260203330085E-003   2.4325727651127111E-003   2.9183453564981714E-003
 fvec =   -9.9637654014044869E-002   4.0371267057609682E-003   2.8954611752227165E-004   5.6006015792153113E-005   1.1883874564059481E-002   3.3153916595763923E-002   2.5526564392639773E-004   5.9338403019988252E-005 -0.11938781840621120       -3.3765623700654897E-002   1.7994418182744251E-009   1.8018748979507199E-004  -4.7817898178824002E-002  -2.2985337277781244E-002   9.7903431267138078E-007   3.8348122511597614E-005   3.0982210373657296E-002   2.4043826056568018E-002   7.4044838875389439E-005   4.1751655220864369E-005   5.7269616939796620E-003  -1.8043230674946339E-002   7.1366880639429659E-005   7.0372542972590443E-005  -6.5384354718974663E-002  -1.3005630054145744E-002   5.8347802786892168E-007   1.0508137415049348E-004  -1.1364751529461033E-002   3.3582110041174791E-002   3.9097949770774650E-004   5.4121348366604376E-004   2.7598335214467525E-002   1.5264990335617498E-002   1.4021997296745383E-005   1.1639906104047464E-004  -1.9892026756130166E-002  -1.0862205359441180E-003   4.1807494908320876E-005   1.4635758343867683E-004
 initial res =    2.5260985425864267E-002 final res =    1.0609263349493416E-008

となります。

Module化

もしこの関数をModuleにするのであれば、

module minimize
    contains
    subroutine find_minimum(calc_fvec,m,n,x0,x,residual)
        implicit none
        include "mkl_rci.fi"
        interface
            subroutine calc_fvec(m,n,x,f)
                integer:: m
                integer::n
                real(8)::x(*) !n-dimensional vector
                real(8)::f(*) !m-dimensional vector
            end subroutine calc_fvec
        end interface


        integer,intent(in)::m,n
        real(8),intent(in)::x0(:)
        real(8),intent(out)::x(:)
        real(8),intent(out)::residual
        real(8),allocatable::fvec(:),fjac(:,:)
        integer(8) :: handle
        integer::res
        real(8)::eps(6)
        real(8)::jac_eps
        integer::iter1,iter2
        real(8)::rs
        integer::i,j
        integer::RCI_Request
        integer::successful
        !** NUMBER OF STOP-CRITERION
        integer::st_cr
        real(8)::r1,r2
        integer::iter
        integer::icount

        jac_eps = 1.d-8
        !** PRECISIONS FOR JACOBI CALCULATION
        jac_eps = 1.d-8
        !** SET MAXIMUM NUMBER OF ITERATIONS
        iter1 = 1000
        !** SET MAXIMUM NUMBER OF ITERATIONS OF CALCULATION OF TRIAL-STEP
        iter2 = 100
        do i = 1, 6
            eps (i) = 1.d-8
        end do
        !** SET INITIAL STEP BOUND
        rs = 100.d0     

        allocate(fvec(1:m),fjac(1:m,1:n))
        x = x0  
        !** SET INITIAL VALUES
        do i = 1, M
            fvec (i) = 0.d0
            do j = 1, N
                fjac (i, j) = 0.d0
            end do
        end do
        res = dtrnlsp_init(handle, N, M, x, eps, iter1, iter2, rs)
        icount = 0
        RCI_Request = 0
        successful  = 0
        do while (successful == 0)        
            icount = icount + 1

            res = dtrnlsp_solve(handle, fvec, fjac, RCI_Request)
            !write(*,*) icount,RCI_Request
            select case (RCI_Request)
            case (-1, -2, -3, -4, -5, -6)
                !**   STOP RCI CYCLE
                successful = 1
            case (1)
!**   RECALCULATE FUNCTION VALUE
!**     M               IN:     DIMENSION OF FUNCTION VALUE
!**     N               IN:     NUMBER OF FUNCTION VARIABLES
!**     X               IN:     SOLUTION VECTOR
!**     FVEC            OUT:    FUNCTION VALUE F(X)                
                call calc_fvec(M, N, X, fvec)
            case (2)
!**   COMPUTE JACOBI MATRIX
!**     EXTENDED_POWELL IN:     EXTERNAL OBJECTIVE FUNCTION
!**     N               IN:     NUMBER OF FUNCTION VARIABLES
!**     M               IN:     DIMENSION OF FUNCTION VALUE
!**     FJAC            OUT:    JACOBI MATRIX
!**     X               IN:     SOLUTION VECTOR
!**     JAC_EPS         IN:     JACOBI CALCULATION PRECISION
            !res = djacobi(fcn, n, m, fjac, x, jac_eps)
                res = djacobi(calc_fvec, N, M, fjac, X, jac_eps)
            end select
        end do           

!** GET SOLUTION STATUSES
!**   HANDLE            IN: TR SOLVER HANDLE
!**   ITER              OUT: NUMBER OF ITERATIONS
!**   ST_CR             OUT: NUMBER OF STOP CRITERION
!**   R1                OUT: INITIAL RESIDUALS
!**   R2                OUT: FINAL RESIDUALS
        res = dtrnlsp_get(handle, iter, st_cr, r1, r2)
        res = dtrnlsp_delete(handle)        
        residual = r2
        write(*,*) "iter = ",iter



        return
    end subroutine




end module minimize


    subroutine testfunc2(M, N, X, F) !
        implicit none
        integer,intent(in):: M, N
        real(8),intent(in)::X (N)
        real(8),intent(out)::F (M)
        integer i,k

        do i = 1, N
            F(i) = i*10 
            do k=1,M

                F(i) = F(i)- x(k)/(i*10 - x(k))
            end do
            F(i) = abs(1d0/(i*10)-1d0/F(i))
        end do

    end subroutine testfunc2

program main
    use minimize
    implicit none
    external testfunc2
    real(8)::residual
    integer::n,m
    real(8),allocatable::x0(:),x(:)
    n = 40
    m = 40
    allocate(x0(n),x(n))
    call random_number(x0)
    call find_minimum(testfunc2,m,n,x0,x,residual)
    write(*,*) "residual = ",residual
end program

とすれば、引数にサブルーチンを取ることができます。