【Intel MKL】 Intel MKL Nonlinear Least Squaresを使う方法

Intel Math Kernel Library (MKL)の中の非線形最小二乗法ルーチンの使い方についてのメモ。ちなみに、自分のマシン環境を前提に書く。具体的には、言語はFortran90、アーキテクチャはIntel(R) 64である。リファレンスはIntel(R) Math Kernel Library for Linux* OS User's Guide January 2010, Doc. N.314774-010US(UG2010)。

• 言語やアーキテクチャによってパスが異なるので注意。

手順は次の通り：

1. MKLのルーチンを呼び出しているソース・コードの頭に、
   program name
   implicit none
   include 'mkl_rci.fi'
   のように、include 'mkl_rci.fi'を付け加える。


2. 次のようにコンパイルする。
   
   (dynamic & threadingの場合)
   ifort filename_1.f90 filename_2.f90 ... filename_x.f90 -lmkl_intel_lp64 -lmkl_intel_thread -lmkl_core -liomp5 -lpthread
   
   (dynamic & sequentialの場合)
   ifort filename_1.f90 filename_2.f90 ... filename_x.f90 -lmkl_intel_lp64 -lmkl_sequential -lmkl_core -lpthread
   
   (static & threadingの場合)
   ifort filename_1.f90 filename_2.f90 ... filename_x.f90 -Wl,--start-group -lmkl_intel_lp64 -lmkl_intel_thread -lmkl_core -Wl,--end-group -liomp5 -lpthread
   
   (static & sequentialの場合)
   ifort filename_1.f90 filename_2.f90 ... filename_x.f90 -Wl,--start-group -lmkl_intel_lp64 -lmkl_sequential -lmkl_core -Wl,--end-group -lpthread
   • ポイントは、staticの場合、cluster components(この場合使用していない)、interface、threading、そしてcomputational librariesの前後を-Wl,--start-groupと-Wl,--end-groupで囲うこと。
   • sequentialの場合、Compiler Support Run-time librariesは関係ないので、オプション-liomp5(compatibility OpenMP run-time libraryで、US2010が推奨)を外している。

Note.1: Intel MKLの構造
Intel MKLは次の4つの階層を持つレイヤー・モデル：

1. Interface layer
2. Threading layer --- threadingがデフォルト
3. Computational layer
4. Compiler Support Run-time libraries

であり、それぞれのレイヤー毎にユーザーが必要なライブラリの選択を行う。なお、linking modelがstaticかdynamicかで選択するファイル名が異なる(アバウトに言えば、staticなら拡張子が.aで、dynamicなら拡張子が.so)。各レイヤーの内容（選択肢）についてはUG2010 Table 3-6、3-7を参照。

注意点は以下の通り（[x]はレイヤーxに関するもの）：

• [1] 2^31-1(約21億)以上の要素を持つ配列を使用する場合、デフォルト・インターフェイスのLP64ではなく、ILP64を指定する必要がある。これは配列のインデックスを与える上での問題。(UG2010, p.3-6)


• [2] MPIを利用する等、特別な理由がない限りthreadingを選択するべき。(UG2010,p.3-5)


• [2] sequentialを選択する場合、Compiler Support Run-time librariesは関係ない。


• [2] sequentialを選択する場合、*sequential.*libraryを選択する。また、個の場合、link lineにPOSIX threads library(pthread)を追加する。


• [2]Optimization Solver routinesのfunction domainはMKLがデフォで使用するスレッディング機能は使われない(UG2010,p.6-1,Using the Intel(R) MKL Parallelism)。そのため、OpenMPのスレッディングの数を指定に関するオプション設定を気にしなくてもOK。


• [3] 非線形最小二乗法はMKL function domainではOptimization (Trust-Region) Solver routinesに対応。この場合、includeファイルmkl_rci.fiを使う。

また、UG2010は、次のように、dynamic linkを推奨している。

You are strongly encouraged to dynamically link in the compatibility OpenMP* run-time library libiomp or legacy OpenMP* run-time library libguide. Link with libiomp and libguide dynamically even if other libraries are linked statically.

Linking to static OpenMP* run-time library is not recommended because it is very easy with complex software to link in more than one copy of the library. This causes performance problems (too many threads) and may cause correctness problems if more than one copy is initialized. (UG2010,p.5-6)

Note.2: 環境設定
.profileに以下の行を追加。

# setting up MKL environment for bash
. absolute_path_to_installed_MKL/tools/environment/mklvarsem64t.sh

.の後は半角スペースが必要。

Note.3: コンパイル&リンク

ifort filename.f90 -L(MKL path) -I(MKL include)
[-I(MKL include)/{32|em64t|{ilp64|lp64}|64/{ilp64|lp64}}]
[-lmkl_blas{95|95_ilp64|95_lp64}]
[-lmkl_lapack{95|95_ilp64|95_lp64}]
[cluster components]
-lmkl_{intel|intel_ilp64|intel_lp64|intel_sp2dp|gf|gf_ilp64|gf_lp64}
-lmkl_{intel_thread|gnu_thread|pgi_thread|sequential}
[-lmkl_lapack] -lmkl_core
{-liomp5|-lguide} [-lpthread] [-lm]

なお、()はパスを表す。ただし、Note.2で環境設定している場合、-L(MKL path)と-I(MKL include)を書く必要はない。

【Fortran】 Intel MKL Nonlinear Least Squares Problem with Bounded Linear Constraint

　非線形連立方程式を解くためにアルゴリズムを探している。

　いくつかあるだろうけど、ここではIntel Compilerに同梱されているMath Kernel Library (MKL)を使う方法を考える。

　MKLの中では、線形制約下の最小二乗法のルーチンを使うことになる。このルーチンはTrust Region Methodを使っている。Trust Region Methodについては、こちらのレジュメがわかりやすい。

（追記）
　この非線形最小二乗のアルゴリズムは使い物にならないかもしれない。というのも、内部でJacobianを求めるときに引数として用意すべき外部サブルーチンの形式への規制が強すぎるから。具体的には、外部サブルーチンが sub(M,N,X,F) （Mは連立方程式の数、Nは変数の数、Xは変数、Fは関数の値）の形でなければならず、連立方程式に登場するパラメータを引数として定義することができなくなっている。もちろん、サブルーチンの中でテキストファイルから読み込むというアプローチも考えられるけど、I/O時間が余計に掛かってしまって計算が遅くなると思われる。

（追記２）
　上記の（追記）は間違い。グローバル変数を使えば大丈夫。具体的には、moduleファイルを使う。

遺伝的アルゴリズム

遺伝的アルゴリズムについてわかりやすく解説を行っているページを発見。

同志社大学工学部知識工学科知的システムデザイン研究室の上浦二郎氏のページ。

制約条件付非線形最適化のアルゴリズム

Mathematicaさんの解説の引用です。

制約条件付きの非線形最適化に対する数値アルゴリズムは，大まかに分けると勾配法と直接探索法とに分けられる．勾配法では，第1導関数（勾配）か第2導関数（ヘッシアン）が使われる．この例には逐次二次計画(SQP)法，拡大ラグランジュ法，非線形内点法がある．直接探索法では，導関数情報は使われない．この例としては，Nelder-Mead法，遺伝的アルゴリズム，微分進化法，焼きなまし法がある．直接探索法の方が収束が遅い傾向があるが，関数と制約条件のノイズの存在への耐性は強い．

目的関数のarugmentの数と最適化処理速度向上

　今、コンパクト集合があって、が連続で、(y,z)に関して厳密に凹な関数であるとする。このとき、次の問題を空間の分割(discretization)を使って、各xについて最適化を行う解を求めたいという状況にあるとする：

　

　さらに、何らかの事情で、この問題を次のように分解して解かなければならないとする：

　

「凸計画問題の分解」と同じ議論を展開すれば、カッコの中はyの厳密な凹関数であることがいえる。

　まず、yの最適化についてはBinary Searchを使えば良いことがわかる。
　では、他にできることはあるだろうか。例えば、xとfの何らかの関係性が仮定されているとしたらどうだろうか。もし、zのargumentがなくてfがxの厳密な増加関数であるなら、は厳密な増加関数であることがいえ、この単調性を使って、関数fの評価の回数を減らすことが可能。しかし、今はzというargumentがあるため、このような単調性を保証するには、さらに何らかの仮定が必要になる。一般論として、fのargumentが増えるほど単調性を導出することが難しくなる。したがって、数値計算がそれだけbrute forceなやり方に近づいてしまうことになる。

数値最適化の実践(その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秒ほど掛かった。問題を少し拡張して操作変数やそれに対応してラグランジュ乗数が約２倍に増えても、計算時間やコンパイルの時間はほとんど変わらなかった。上記のようにテンプレートさえ作ってしまえばプログラムを書くのは非常に楽。一番時間が掛かるのは目的関数を記述する部分。