本記事は、Chainer Advent Calendar 2017の16日目です．

TL;DR

gy, = chainer.grad([y], [x], enable_double_backprop=True)
ggy, = chainer.grad([gy], [x])

y.backward(enable_double_backprop=True)
gy = x.grad_var
x.cleargrad()  # reset as x.grad = 0
gy.backward()
ggy = x.grad_var

はじめに

最近、Chainerにはdouble backpropagation(微分値の自動微分)の機能が追加されました。それというのも、最近のモデルは1階微分では事足らず、Unrolled GANの実装やWGANの精度改善のために2階微分が必要だったりするからです。そこまで使用頻度は高くないかもしれないですが、Chainerでどのように2階微分を計算すればよいか説明できればと思います。

自動微分とは

自動微分とは、「プログラムで定義された関数を解析し、偏導関数の値を計算するプログラムを導出する技術である。自動微分は複雑なプログラムであっても加減乗除などの基本的な算術演算や基本的な関数（指数関数・対数関数・三角関数など）のような基本的な演算の組み合わせで構成されていることを利用し、これらの演算に対して連鎖律を繰り返し適用することによって実現される。」[1]というものです。

Chainerなどの多くのDeep Learning(DL)フレームワークでは、下記のように数式を計算グラフとして表現し、その情報を用いて自動的に微分を行います。

そもそもどこで自動微分が使われているか

世の中には多くのDLフレームワークが存在しますが、これらのフレームワークが提供する機能は大別すると3つの共通する機能を提供しています。

多次元配列の計算機能
ニューラルネットワーク(NN)の定義機能
NNの最適化機能

自動微分は、「NNの定義機能」と「NNの最適化機能」に関わる話です。ChainerなどのDeep Learingフレームワークを使用していると、数式でモデルを定義し、そのモデルを最適化クラス(Chainerで言うところのTrainer、Optimizerなど)に渡すと自動でモデルパラメータの学習ができてしまいます。それなので、一見すると微分計算なんてどこで必要なのか？と思ってしまうかもしれないですが、実はモデルパラメータの最適化にパラメータによる損失関数の微分値を使用しています。さらにその微分値は、数式の表現できる全てのモデルに対して、微分値を計算できる必要があります。そのため、自動微分という技術が必要になってきます。

もし、自動微分という技術を導入していない場合どうなるかというと、新しくモデルを定義するたびに、そのモデルパラメータによる損失関数の微分値を計算するプログラムを別途定義しなくてはならなくなります。

自動微分、double backpropagationの計算方法

この資料のp15以降がかなり詳しいのでこちらにおまかせします。すみません。。
comparison-of-deep-learning-frameworks-from-viewpoint-of-double-backpropagation

Chainerでのdouble backpropagation

chainer.gradを使う場合は、enable_double_backprop=Trueとすると、微分可能な型Variableとして返却されるので、さらに微分をすることができます。

import chainer
from chainer import Variable
import numpy as np

x = Variable(np.array(2.0))  # x = 2.0
y = x ** 2. + x + 1.0  # y = x^2 + x + 1

gy, = chainer.grad([y], [x], enable_double_backprop=True)  # gy = dy/dx = 2x + 1
print(gy)  # variable(5.0)

ggy, = chainer.grad([gy], [x])  # ggy = dydy/dxdx = 2
print(ggy)  # variable(2.0)

Variable.backwardを使う場合は、enable_double_backprop=Trueとすると、Variable.grad_varで微分可能な型Variableを取得することができ、さらに微分をすることができます。

import chainer
from chainer import Variable
import numpy as np

x = Variable(np.array(2.0))  # x = 2.0
y = x ** 2. + x + 1.0  # y = x^2 + x + 1

y.backward(enable_double_backprop=True)
gy = x.grad_var  # gy = dy/dx = 2x + 1
print(gy)  # variable(5.0)

x.cleargrad()  # reset as x.grad = 0

gy.backward()
ggy = x.grad_var  # ggy = dydy/dxdx = 2
print(ggy)  # variable(2.0)

おまけ

また、最後には戯事ですが、2階微分が必要な応用として金融派生商品のグリークス(価格の微分値)を計算してみました。今回使用しているモデルはブラック–ショールズモデル[3]です。(あとで時間があったら詳しく説明できればと思います。)

import numpy as np
import chainer
from chainer import Variable
import chainer.functions as F
from math import exp


def main():
    stock_price = 100.0
    strike = 100.0
    years = 1.0
    risk_free_rate = 0.01
    volatility = 0.1
    n_paths = 100000
    n_steps = 100

    times = np.linspace(0, years, n_steps, dtype=np.float64)

    r = risk_free_rate
    v = volatility
    k = strike
    T = years
    initial_stocks = Variable(stock_price * np.ones((n_paths,), dtype=np.float64))
    currents = initial_stocks
    for c, n in zip(times[:-1], times[1:]):
        dt = n - c
        zs = np.random.normal(0.0, 1.0, currents.shape)
        currents = currents * F.broadcast(np.exp((r - 0.5 * v * v) * dt + v * dt ** 0.5 * zs))

    call_price = F.average(exp(-r * T) * F.softplus(currents - k, beta=10))  # you should smooth the relu to calculate gamma
    # call_price = F.average(exp(-r * T) * F.relu(currents - k))  
    print('call price: {0}'.format(call_price))
    deltas = F.sum(chainer.grad([call_price], [initial_stocks], enable_double_backprop=True)[0])
    print('call delta: {0}'.format(deltas))
    gammas = chainer.grad([deltas], [initial_stocks], enable_double_backprop=True)
    print('call gamma: {0}'.format(np.sum(gammas)))


if __name__ == '__main__':
    main()