モーダルを閉じる工作HardwareHub ロゴ画像

工作HardwareHubは、ロボット工作や電子工作に関する情報やモノが行き交うコミュニティサイトです。さらに詳しく

利用規約プライバシーポリシー に同意したうえでログインしてください。

関数の微分をバックプロパゲーションで求める (自動微分)

モーダルを閉じる

ステッカーを選択してください

お支払い手続きへ
モーダルを閉じる

お支払い内容をご確認ください

購入商品
」ステッカーの表示権
メッセージ
料金
(税込)
決済方法
GooglePayマーク
決済プラットフォーム
確認事項

利用規約をご確認のうえお支払いください

※カード情報はGoogleアカウント内に保存されます。本サイトやStripeには保存されません

※記事の執筆者は購入者のユーザー名を知ることができます

※購入後のキャンセルはできません

作成日作成日
2020/07/04
最終更新最終更新
2022/11/16
記事区分記事区分
一般公開

自動微分というアルゴリズムによって関数の微分値を求める例を記載します。本ページでは、一階微分を対象としており、高階微分は考えません。また簡単のため、関数の出力は一つのテンソルであり、入力となる一つ以上のテンソルおよび出力となる一つのテンソルの階数は零である例を考えます。

更に、関数は、微分可能な関数からなる合成関数であることを仮定します。これは自動微分の応用先の一つである、ディープラーニングにおけるニューラルネットワークの性質の一つです: Chainer Tutorial / ニューラルネットワークの基礎 / 誤差逆伝播法(バックプロパゲーション)

PyTorch による例

PyTorch はテンソル (多次元配列) を扱うことができるライブラリで、ディープラーニングなどに用いられます。スカラは零階のテンソル、ベクトルは一階のテンソル、行列は二階のテンソルです。

以下のようにインストールできます

python3 -m pip install torch torchvision

NVIDIA GPU が利用できる環境であれば以下の値が True になります。

python3 -m IPython

In [1]: import torch
   ...: torch.cuda.is_available()
Out[1]: False

PyTorch には自動微分が実装されており、以下のように利用できます。入力となる二つのテンソル xxyy の階数は零です。出力となる一つのテンソル zz の階数も零です。

import torch

x = torch.tensor(1.0, requires_grad=True)
y = torch.tensor(1.0, requires_grad=True)

z = x ** 2 + y ** 2

z.backward()

print(x.grad, y.grad)  # tensor(2.) tensor(2.)

数式微分を考えても x2+y2x^2 + y^2xx および yy に関する偏微分はどちらも 2 となるため、上記結果は正しいことが分かります。必要に応じて数値微分で自動微分の結果を検証することもできます。

自動微分のアルゴリズム

『ゼロから作るDeep Learning ❸』のソースコードはこちらで公開されており、内容を確認すると、自動微分のアルゴリズムは以下のようになっていることが分かります。

main.py

#!/usr/bin/python
# -*- coding: utf-8 -*-

from autograd.variable import Variable

import numpy as np

def Main():

    # 零階のテンソルを二つ入力します。
    x = Variable(np.array(1.0))
    y = Variable(np.array(1.0))

    # 微分可能な関数からなる合成関数です。
    # すべて一つのテンソルを出力しています。
    for z in (Sphere(x, y),
              Matyas(x, y),
              Goldstein(x, y)):

        # 前のループで格納された、微分結果を忘れます。
        x.ClearGrad()
        y.ClearGrad()

        # バックプロパゲーションによって微分します。
        z.Backward()

        # 結果を確認します。
        print(x.GetGrad(), y.GetGrad())

def Sphere(x, y):
    z = x ** 2 + y ** 2
    return z

def Matyas(x, y):
    z = 0.26 * (x ** 2 + y ** 2) - 0.48 * x * y
    return z

def Goldstein(x, y):
    a = 1 + (x + y + 1) ** 2 * (19 - 14 * x + 3 * x ** 2 - 14 * y + 6 * x * y + 3 * y ** 2)
    b = 30 + (2 * x - 3 * y) ** 2 * (18 - 32 * x + 12 * x ** 2 + 48 * y - 36 * x * y + 27 * y ** 2)
    z = a * b
    return z

if __name__ == '__main__':
    Main()

autograd/__init__.py

Python パッケージのために必要な空ファイルです。

# -*- coding: utf-8 -*-

autograd/variable.py

__repr____len__ は Python の特殊メソッドです。同様に __add____mul__ も特殊メソッドで、これらは Python における演算子オーバーロードのために用いられます。

# -*- coding: utf-8 -*-

import numpy as np

class Variable(object):

    def __init__(self, data, name=None):
        if data is not None:
            if not isinstance(data, np.ndarray):
                raise TypeError('{} is not supported'.format(type(data)))
        self.__name = name
        self.__data = data  # テンソル
        self.__grad = None  # テンソル
        self.__creator = None  # この変数を出力した Function
        self.__generation = 0  # この変数の世代

    def __repr__(self):
        if self.__data is None:
            return 'Variable(None)'
        p = str(self.__data).replace('\n', '\n' + ' ' * 9)
        return 'Variable(' + p + ')'

    def __len__(self):
        return len(self.__data)

    def __neg__(self):
        return neg(self)

    def __add__(self, other):
        return add(self, other)

    def __radd__(self, other):
        return add(self, other)

    def __mul__(self, other):
        return mul(self, other)

    def __rmul__(self, other):
        return mul(self, other)

    def __sub__(self, other):
        return sub(self, other)

    def __rsub__(self, other):
        return sub(other, self)  # 引数の順番を逆に設定

    def __truediv__(self, other):
        return div(self, other)

    def __rtruediv__(self, other):
        return div(other, self)  # 引数の順番を逆に設定

    def __pow__(self, other):
        return pow(self, other)

    # NumPy に実装されている演算子オーバーロードよりも
    # Variable クラスに実装されている演算子オーバーロードを
    # 優先させるための設定。
    # https://numpy.org/doc/stable/reference/arrays.classes.html#numpy.class.__array_priority__
    __array_priority__ = 200

    def GetName(self):
        return self.__name

    def GetData(self):
        return self.__data

    def GetCreator(self):
        return self.__creator

    def SetCreator(self, func):
        self.__creator = func
        # 変数の「世代」は、変数を作り出した関数の「世代」に 1 を加えたものとしています。
        self.__generation = func.GetGeneration() + 1

    def GetGeneration(self):
        return self.__generation

    def GetGrad(self):
        return self.__grad

    def SetGrad(self, grad):
        self.__grad = grad

    def ClearGrad(self):
        self.__grad = None

    def Backward(self, retainGrad=False):

        if self.__grad is None:
            # バックプロパゲーションの起点となる変数の微分値を 1 に設定。
            self.__grad = np.ones_like(self.__data)

        funcs = []
        seenSet = set()

        def AddFunc(f):
            if f in seenSet:
                return
            funcs.append(f)
            seenSet.add(f)
            funcs.sort(key=lambda x: x.GetGeneration())

        AddFunc(self.__creator)

        while len(funcs) > 0:

            # 同じ変数を異なる複数の関数に入力した場合への対応。
            # e.g.
            # --> x --> Square(x) : a --> Add(a, b) : c
            #       --> Exp(x)    : b -/
            # 世代の大きい関数から順に処理します。
            f = funcs.pop()
            gys = [output().GetGrad() for output in f.GetOutputs()]
            gxs = f.Backward(*gys)

            if not isinstance(gxs, tuple):
                gxs = gxs,

            for x, gx in zip(f.GetInputs(), gxs):
                if x.GetGrad() is None:
                    x.SetGrad(gx)
                else:
                    # 同じ変数を同じ関数に複数回数入力した場合への対応
                    # e.g. Add(x, x)
                    x.SetGrad(x.GetGrad() + gx)

                if x.GetCreator() is not None:
                    AddFunc(x.GetCreator())

            if retainGrad:
                return

            # バックプロパゲーションの処理が完了した関数から出力された変数の微分値は、
            # 今後の処理で不要なため、オプションが指定されていない場合はメモリ解放のために忘れます。
            for y in f.GetOutputs():
                y().SetGrad(None)

from .function import neg
from .function import add
from .function import mul
from .function import sub
from .function import div
from .function import pow

autograd/function.py

弱参照 weakref によって、循環参照によるメモリリークを回避しています。

# -*- coding: utf-8 -*-

import numpy as np
import weakref

from .variable import Variable

class Function(object):

    def __init__(self):
        self.__inputs = None
        self.__outputs = None
        self.__generation = None

    def __call__(self, *inputs):
        inputs = [AsVariable(x) for x in inputs]
        xs = [x.GetData() for x in inputs]
        ys = self.Forward(*xs)
        if not isinstance(ys, tuple):
            ys = ys,
        outputs = [Variable(AsArray(y)) for y in ys]

        # 関数の世代を、入力変数の世代の最大値、として設定します。
        self.__generation = max([x.GetGeneration() for x in inputs])
        for output in outputs:
            output.SetCreator(self)

        self.__inputs = inputs
        self.__outputs = [weakref.ref(output) for output in outputs]  # 弱参照

        return outputs if len(outputs) > 1 else outputs[0]

    def GetInputs(self):
        return self.__inputs

    def GetOutputs(self):
        return self.__outputs

    def GetGeneration(self):
        return self.__generation

    def Forward(self, x):
        raise NotImplementedError()

    def Backward(self, gy):
        raise NotImplementedError()

def AsArray(x):
    if np.isscalar(x):
        return np.array(x)
    return x

def AsVariable(obj):
    if isinstance(obj, Variable):
        return obj
    return Variable(obj)

class Neg(Function):

    def Forward(self, x):
        return -x

    def Backward(self, gy):
        return -gy

class Add(Function):

    def Forward(self, x0, x1):
        y = x0 + x1
        return y

    def Backward(self, gy):
        return gy, gy

class Mul(Function):

    def Forward(self, x0, x1):
        y = x0 * x1
        return y

    def Backward(self, gy):
        x0 = self.GetInputs()[0].GetData()
        x1 = self.GetInputs()[1].GetData()
        return gy * x1, gy * x0

class Sub(Function):

    def Forward(self, x0, x1):
        y = x0 - x1
        return y

    def Backward(self, gy):
        return gy, -gy

def Div(Function):

    def Forward(self, x0, x1):
        y = x0 / x1
        return y

    def Backward(self, gy):
        x0 = self.GetInputs()[0].GetData()
        x1 = self.GetInputs()[1].GetData()
        gx0 = gy / x1
        gx1 = gy * (-x0 / x1 ** 2)
        return gx0, gx1

class Pow(Function):

    def __init__(self, c):
        self.__c = c

    def Forward(self, x):
        y = x ** self.__c
        return y

    def Backward(self, gy):
        x = self.GetInputs()[0].GetData()
        c = self.__c
        gx = c * x ** (c - 1) * gy
        return gx

class Exp(Function):

    def Forward(self, x):
        return np.exp(x)

    def Backward(self, gy):
        x = self.GetInputs()[0].GetData()
        gx = np.exp(x) * gy
        return gx

def neg(x):
    return Neg()(x)

def add(x0, x1):
    x0 = AsArray(x0)
    x1 = AsArray(x1)
    return Add()(x0, x1)

def mul(x0, x1):
    x0 = AsArray(x0)
    x1 = AsArray(x1)
    return Mul()(x0, x1)

def sub(x0, x1):
    x0 = AsArray(x0)
    x1 = AsArray(x1)
    return Sub()(x0, x1)

def div(x0, x1):
    x0 = AsArray(x0)
    x1 = AsArray(x1)
    return Div()(x0, x1)

def pow(x, c):
    return Pow(c)(x)

def exp(x):
    return Exp()(x)

実行例

$ python3 main.py 
2.0 2.0
0.040000000000000036 0.040000000000000036
-5376.0 8064.0

確かに PyTorch の出力結果と一致することが分かります。

ヤコビ行列について

(y1,...,ym)(y_1, ..., y_m) それぞれが (x1,...,xl)(x_1, ..., x_l) の関数であり、yjy_jxix_i で偏微分可能である場合について、以下の行列をヤコビ行列とよびます。ヤコビ行列の行列式をヤコビアンとよびます。

Jxy=(y1x1...y1xl.........ymx1...ymxl)J_{xy} = \left( \begin{array}{rr} \frac{\partial y_1}{\partial x_1} & ... & \frac{\partial y_1}{\partial x_l} \\ ... & ... & ... \\ \frac{\partial y_m}{\partial x_1} & ... & \frac{\partial y_m}{\partial x_l} \\ \end{array} \right)

ヤコビ行列は例えば関数の一次近似を扱う場合に利用できます。具体例として、ロボットアームのマニピュレータの位置姿勢 x1,x2,x3,ϕ1,ϕ2,ϕ3x_1, x_2, x_3, \phi_1, \phi_2, \phi_3 は関節値 θ1,...,θl\theta_1, ..., \theta_l の関数です。マニピュレータの位置姿勢はヤコビ行列を用いて以下のように近似できます。

(x1(θ)x2(θ)x3(θ)ϕ1(θ)ϕ2(θ)ϕ3(θ))=y(θ)y(θ0)+J(θ0)(θθ0)\left( \begin{array}{c} x_1(\vec{\theta}) \\ x_2(\vec{\theta}) \\ x_3(\vec{\theta}) \\ \phi_1(\vec{\theta}) \\ \phi_2(\vec{\theta}) \\ \phi_3(\vec{\theta}) \end{array} \right) = \vec{y}(\vec{\theta}) \simeq \vec{y}(\vec{\theta_0}) + J(\vec{\theta_0}) (\vec{\theta} - \vec{\theta_0})

連鎖律 (Chain Rule) について

さらに (z1,...,zn)(z_1, ..., z_n) それぞれが (y1,...,ym)(y_1, ..., y_m) の関数であり、zkz_kyjy_j で偏微分可能である場合について、ヤコビ行列を JyzJ_{yz} とします。

Jyz=(z1y1...z1ym.........zny1...znym)J_{yz} = \left( \begin{array}{rr} \frac{\partial z_1}{\partial y_1} & ... & \frac{\partial z_1}{\partial y_m} \\ ... & ... & ... \\ \frac{\partial z_n}{\partial y_1} & ... & \frac{\partial z_n}{\partial y_m} \\ \end{array} \right)

このとき、(z1,...,zn)(z_1, ..., z_n)(x1,...,xl)(x_1, ..., x_l) を用いた一次近似は以下のようになります。

\begin{eqnarray} \vec{z}(\vec{x}) &\simeq& \vec{z}(\vec{x_0}) + J_{xz}(\vec{x_0}) (\vec{x} - \vec{x_0}) \\ &=& \vec{z}(\vec{x_0}) + J_{yz}(\vec{y_0}) J_{xy}(\vec{x_0}) (\vec{x} - \vec{x_0}) \\ \end{eqnarray}

Jxz=JyzJxyJ_{xz} = J_{yz} J_{xy} の部分は連鎖律 (Chain Rule) とよばれています。

zkxi=j=1mzkyjyjxi\frac{\partial z_k}{\partial x_i} = \sum_{j=1}^{m} \frac{\partial z_k}{\partial y_j} \frac{\partial y_j}{\partial x_i}

さきほどの自動微分のサンプルプログラムにおけるバックプロパゲーションのアルゴリズムでは、この Chain Rule が利用されています。例えば、乗算 class MulBackward 関数は以下のようになっています。

class Mul(Function):

    def Forward(self, x0, x1):
        y = x0 * x1
        return y

    def Backward(self, gy):
        x0 = self.GetInputs()[0].GetData()
        x1 = self.GetInputs()[1].GetData()
        return gy * x1, gy * x0

一般性を失うことなく kk が 1 の場合を考えます。例えば以下のような合成関数 z1(x1,x2)z_1(x_1, x_2)y2(x1,x2)y_2(x_1, x_2)x1x_1x2x_2 の乗算である場合、

\begin{eqnarray} z_1(x_1, x_2) &=& y_1(x_1, x_2) + y_2(x_1, x_2) \\ &=& y_1(x_1, x_2) + x_1 x_2 \end{eqnarray}

連鎖律は以下のようになります。

\begin{eqnarray} \frac{\partial z_1}{\partial x_1} &=& \frac{\partial z_1}{\partial y_1} \frac{\partial y_1}{\partial x_1} + \frac{\partial z_1}{\partial y_2} \frac{\partial y_2}{\partial x_1} \\ &=& \frac{\partial z_1}{\partial y_1} \frac{\partial y_1}{\partial x_1} + \frac{\partial z_1}{\partial y_2} x_2 \end{eqnarray}
\begin{eqnarray} \frac{\partial z_1}{\partial x_2} &=& \frac{\partial z_1}{\partial y_1} \frac{\partial y_1}{\partial x_2} + \frac{\partial z_1}{\partial y_2} \frac{\partial y_2}{\partial x_2} \\ &=& \frac{\partial z_1}{\partial y_1} \frac{\partial y_1}{\partial x_2} + \frac{\partial z_1}{\partial y_2} x_1 \end{eqnarray}

z1y2x2\frac{\partial z_1}{\partial y_2} x_2 および z1y2x1\frac{\partial z_1}{\partial y_2} x_1 の部分が def Backward の実装内容です。

Likeボタン(off)0
詳細設定を開く/閉じる
アカウント プロフィール画像

Gitの便利な機能であるInteractive Rebaseを勉強中

記事の執筆者にステッカーを贈る

有益な情報に対するお礼として、またはコメント欄における質問への返答に対するお礼として、 記事の読者は、執筆者に有料のステッカーを贈ることができます。

>>さらに詳しくステッカーを贈る
ステッカーを贈る コンセプト画像

Feedbacks

Feedbacks コンセプト画像

    ログインするとコメントを投稿できます。

    ログインする

    関連記事