数値最適化の実践(その2)

　インテルのコンパイラー（プロフェッショナル版）に同梱されているMKLを使って「数値最適化の実践(その1)」の問題を解くことを考える。

　Karush-Kuhn-Tucker(KKT)条件の一部から以下の式が得られる：





　一方で、非負制約条件として、

がある。

　したがって、非線形最小二乗法の問題は

となる。これまでの議論から、解は一意である。

　以下のコードはインテルのFortranコンパイラのサンプルを参考に（というかほぼ同じように）書いたもの：

program NLP
implicit none
include 'mkl_rci.fi'
!-------------------------------------------------
! DECLARATIONS
!-------------------------------------------------
! 1. DTRNLSPBC_INIT
integer, parameter :: N=4, M=4
double precision X(N), LW(N), UP(N), eps(6), rs
integer iter1, iter2

! 2. DTRNLSPBC_SOLVE
double precision fvec(M), fjac(M,N)
integer RCI_Request, SUCCESSFUL

! 3. DTRNLSPBC_GET
integer iter, st_cr
double precision r1, r2

! 4. DJACOB
double precision jac_eps

! 5. COMMONS
integer*8 handle

!-------------------------------------------------
! VALUES FOR PARAMETERS
!-------------------------------------------------
! maximul iteration
iter1=1000
iter2=100

! bounds
LW = 0.D0
UP(1)=800D0
UP(2)=1.D0
UP(3)=1000D0
UP(4)=1000D0
rs=100.D0

! criterions
eps=1.D-5
jac_eps=1.D-8

!-------------------------------------------------
! INITIALIZATION
!-------------------------------------------------
! initial guess
X(1) = 0.5D0
X(2) = 0.5D0
X(3) = 10D0
X(4) = 0.D0

! objective function
fvec=0.D0

! Jacobian matrix
fjac=0.D0

!-------------------------------------------------
! EXECUTION
!-------------------------------------------------
! 1. DTRNLSPBC_INIT
if (dtrnlspbc_init(handle,N,M,X,LW,UP,eps,iter1,iter2,rs) /= TR_SUCCESS) then
print *, '| ERROR IN DTRNLSPBC_INIT'
call MKL_FREE_BUFFERS
stop 1
endif

RCI_Request=0
SUCCESSFUL=0

do while (SUCCESSFUL == 0)
if (dtrnlspbc_solve(handle,fvec,fjac,RCI_Request)/=TR_SUCCESS) then
print *, '| ERROR IN DTRNLSPBC_SOLVE'
call MKL_FREE_BUFFERS
stop 1
endif

select case (RCI_Request)
case (-1,-2,-3,-4,-5,-6)
SUCCESSFUL=1
case (1) ! recompute function value
call object(M,N,X,fvec)
case (2) ! compute the Jacobian matrix
if (djacobi(object,N,M,fjac,X,jac_eps)/=TR_SUCCESS) then
print *, '| ERROR IN DJACOBI'
call MKL_FREE_BUFFERS
stop 1
endif
endselect

end do

! 2. DTRNLSPBC_GET
if (dtrnlspbc_get(handle,iter,st_cr,r1,r2)/=TR_SUCCESS) then
print *, '| ERROR IN DTRNLSPBC_GET'
call MKL_FREE_BUFFERS
stop 1
endif

! 3. DTRNLSPBC_DELETE
if (dtrnlspbc_delete(handle)/=TR_SUCCESS) then
print *, '| ERROR IN DTRNLSPBC_DELETE'
call MKL_FREE_BUFFERS
stop 1
endif

! FINAL CRITERION
if (r2<1.D-1) then
print *, '| OPTIMIZATION ............PASS'
stop 0
else
print *, '| OPTIMIZATION ............FAILED'
stop 1
endif

contains

! OBJECTIVE FUNCTION
subroutine object(M,N,X,F)
integer M, N
double precision X(*), F(*)
double precision, parameter :: alpha = 0.64D0, w = 5.D0, y = 0.5D0

! X=(c, ell, lambda, mu)
F(1)=X(1)*((1-alpha)*(X(2)/X(1))**alpha-X(3))
F(2)=X(2)*(alpha*(X(1)/X(2))**(1-alpha)-w*X(3)-X(4))
F(3)=X(3)*(y+w*(1-X(2))-X(1))
F(4)=X(4)*(1-X(2))

end subroutine

end program

　パラメータの値は
α=0.64, w=5, y=1/2
に設定してある。

　計算の結果は
c=1.97999999984257,
ell=0.704000000162305.
一方で、解析的な解は
c=1.98,
ell=0.704.
無視できる誤差で、ほぼ正確に求められている。

　計算時間は約0.047秒。ただし、これはアブソリュート・コード(exeファイル)実行時間。一方で、コンパイルの時間は6秒ほど掛かった。問題を少し拡張して操作変数やそれに対応してラグランジュ乗数が約２倍に増えても、計算時間やコンパイルの時間はほとんど変わらなかった。上記のようにテンプレートさえ作ってしまえばプログラムを書くのは非常に楽。一番時間が掛かるのは目的関数を記述する部分。