Menu

     << (newer post) (older post) >>     

自動微分パッケージ ForwardDiff

#numerical-analysis #package


Julia の便利なパッケージの一つがこの ForwardDiff. いや～、ありがたやありがたや．

なんとこいつは、関数 $f$ を与えると自動微分をしてくれる． 自動微分ってのは、関数計算を tree に展開して、その枝に各計算要素の微分係数を重みとして貼り付けることによって、関数計算のわずか数倍程度の計算量で偏微分値関数を生成するという手法 1 で、数値解析のきちんとした本なんかに載っているやつだ．

つまり、$\mbox{grad } f ( = \nabla f)$ やら Hessian やら Jacobian やらを計算してくれるのだ． これらの量は重要で、例えば非線形連立方程式を解くために Newton 法を使おうとすると 2 、その方程式を関数とみなしたときの Jacobian が必要になる． こうした場合、連立の次元数 $N$ は普通に 1000 とか 10000 などという数になるが、Jacobian は $N \times N$ 行列になるので、よほどきれいな行列にならないかぎり、こうした計算を手で行うのは、大変すぎるし、間違えやすい． しかも、境界条件もからんでくると例外処理が発生するので、こうした問題で手計算する場合は大変に気を使うことになる.

こういう計算ほど計算機にやってもらえるとありがたいわけで、そうした意味でこのパッケージの存在はとても助かる．

ちょっと使い方のサンプルを書くと、次のような感じかな．

### パッケージを読み込んで…
using ForwardDiff

### ベクトルをもらって同じ次元のベクトルを返す関数例．周期的境界条件下での u * u_x をイメージしている．
function f(u::Vector)
    v = similar(u)
    N = length(v)
    v[1] = u[1]*(u[2]-u[N])
    v[N] = u[N]*(u[1]-u[N-1])
    for k in 2:N-1
        v[k] = u[k]*(u[k+1]-u[k-1])
    end
    return v
end

### Jacobian を計算するのもコマンド一つで！ 返り値は j(x) という関数(出力は行列)になることに注意．
j = jacobian(f)

# 例えば 5次元ベクトルを試しに用意して…
x = collect(1.0:1.0:5.0)

# 入力が x のときの、$f$ の Jacobian を具体的に出力させる．
j(x)

とすると、この $x$ のときの関数 $f$ の Jacobian が

  5x5 Array{Float64,2}:
   -3.0   1.0   0.0   0.0  -1.0
   -2.0   2.0   2.0   0.0   0.0
    0.0  -3.0   2.0   3.0   0.0
    0.0   0.0  -4.0   2.0   4.0
    5.0   0.0   0.0  -5.0  -3.0

という感じに具体的に値で得ることが出来る．

大変素晴らしいのは、関数の次元数など、具体的な数値を与える「前」に計算が可能な点． 理論的にはそりゃ可能なんだけど、実装がそうなされているというのはほんとうに素晴らしい．ビバ自動微分！

ちょっとだけ残念な制限を書いておくと、現バージョンの ForwardDiff では、関数 $f$ がもてる引数は一つだけで、二つ以上持っていてはいけない． これはまじめに考えると実用上は結構不便な制限で、$f$ が必要とするパラメータは global 変数にするか関数の定義の中に埋め込むかしかなく、プログラムの局所化がしにくくなってしまう． …はずなんだけれども、次のようにして事実上回避可能(グローバル変数にしていることに近いけどね)．ありなんだろかね、これ…

using ForwardDiff

# パラメータ up と、本質的な引数 u をもつ関数 g を考えて、
function g(u,up)
  u .* up
end

# パラメータの実際は vp として、
vp = collect(1.0:1.0:5.0);

# 関数の引数は v に格納して、
v = randn(5);

# 形式的に引数を一つだけ持つ関数 f を定義して、
f(u) = g(u, vp)

# 位置 v での f の Jacobian を計算してみると、
jacobian(f, v)

ちゃんとパラメータの影響を受けた Jacobian が次のように得られる．

  5x5 Array{Float64,2}:
   1.0  0.0  0.0  0.0  0.0
   0.0  2.0  0.0  0.0  0.0
   0.0  0.0  3.0  0.0  0.0
   0.0  0.0  0.0  4.0  0.0
   0.0  0.0  0.0  0.0  5.0

なんか落とし穴がありそうな気がしてるんだけどどうだろう？ ちょっとさらにテストしてみる．

まず、パラメータ $\bm{vp}$ は変えるが、他はきちんと動くように同様に作業する． つまり、パラメータ変更後にあらためて関数 $\bm{f}$ をきちんと再定義しなおしてから、計算して動作が正しく行われるかみてみる．

# パラメータ変更
vp = collect(11.0:1.0:15.0);

# 念の為に f を再定義
f(u) = g(u, vp)

# Jacobian 計算
jacobian(f, v)

ふむ、結果は以下のとおり．正しく動く．まあこれは良い．

  5x5 Array{Float64,2}:
   11.0   0.0   0.0   0.0   0.0
    0.0  12.0   0.0   0.0   0.0
    0.0   0.0  13.0   0.0   0.0
    0.0   0.0   0.0  14.0   0.0
    0.0   0.0   0.0   0.0  15.0

では、関数 $\bm{f}$ の再定義もさぼったらどうなる？ やってみる．

# パラメータをまた変更
vp = collect(21.0:1.0:25.0);

# いきなり Jacobian を計算．
jacobian(f, v)

結果は以下のとおり．やはり正しく動く．う～ん…

  5x5 Array{Float64,2}:
   21.0   0.0   0.0   0.0   0.0
    0.0  22.0   0.0   0.0   0.0
    0.0   0.0  23.0   0.0   0.0
    0.0   0.0   0.0  24.0   0.0
    0.0   0.0   0.0   0.0  25.0

結局のところ、コンパイラが関数評価をどのタイミングでどう行っているかが真面目にわからないといけないんだろうけれども、ともかく遅延評価になってるのは確かなようなので意図したとおりには使えるようだ． 本気で調べるのは面倒だからそのうちいつか、な.

  1. 関数の数値誤差限界の計算と本質的には同じ計算手法である． [return]
  2. もちろん、非線形連立方程式を解く手法は他にもたくさんあるし、それ用のパッケージ NLsolve もある．でもこの NLsolve も「Jacobian がわかるならばより効果的に使えるもの」なので、やっぱこの ForwardDiff はあったほうがいいんだな、これが． [return]
     << (newer post) (older post) >>