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

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