C++
数値計算

二重数で自動微分する

More than 1 year has passed since last update.

参考文献 [1] に 二重数 (dual number) で自動微分 (automatic differentiation) する方法が載っていたので試してみた。


二重数

二重数は、複素数のように実数を拡張した数で

a + b \epsilon

と記述する数である。ここで、$\epsilon$ が複素数の $i$ にあたるもので、$\epsilon^2 = 0$ という性質を持つものとして定義されている。

加減乗除も複素数と同様に定義されて、C++ で実装すると、


dn.hpp

#pragma once


#include <ostream>

namespace math {

// 二重数 a + b ε
template <typename T = double>
class dual_number {
public:
using this_type = dual_number<T>;

public:
dual_number(T a, T b = T(0)) : a_(a), b_(b) {}

this_type operator - () const {
return this_type(-a_, -b_);
}

this_type& operator += (this_type const& rhs) {
a_ += rhs.a_;
b_ += rhs.b_;
return *this;
}

this_type& operator -= (this_type const& rhs) {
a_ -= rhs.a_;
b_ -= rhs.b_;
return *this;
}

this_type& operator *= (this_type const& rhs) {
b_ = (b_ * rhs.a_ + a_ * rhs.b_);
a_ *= rhs.a_;
return *this;
}

this_type& operator /= (this_type const& rhs) {
b_ = (b_ * rhs.a_ - a_ * rhs.b_) / (rhs.a_ * rhs.a_);
a_ /= rhs.a_;
return *this;
}

this_type operator + (this_type const& rhs) const {
return this_type(*this) += rhs;
}

this_type operator - (this_type const& rhs) const {
return this_type(*this) -= rhs;
}

this_type operator * (this_type const& rhs) const {
return this_type(*this) *= rhs;
}

this_type operator / (this_type const& rhs) const {
return this_type(*this) /= rhs;
}

friend std::ostream& operator << (std::ostream& out, this_type const& rhs) {
return out << '[' << rhs.a_ << ',' << rhs.b_ << ']';
}

private:
T a_;
T b_;
};
}


みたいになる。


自動微分

関数 $f(x)$ のテイラー展開は

f(x + \delta) = f(x) + \left( \frac{d}{dx} f(x) \right) \delta + O(\delta^{2})

である。

ここで $\delta$ を $\epsilon$ とおいて $x + \epsilon$ を二重数とみなすと、 $\epsilon^{n}$ の $n$ が $2$ 以上の項は $0$ となり

f(x + \epsilon) = f(x) + \left( \frac{d}{dx} f(x) \right) \epsilon

が得られる。

これは、関数 $f(x + \epsilon)$ を計算すると、実数部として $f(x)$ が、$\epsilon$ の係数として 微分値 $\frac{d}{dx} f(x)$ が得れることを示している。


実行結果

\begin{align}

f(x) &= 4 x^2 + 3 x + 2 \\
\frac{d}{dx} f(x) &= 8 x + 3
\end{align}

を計算してみる。


main.cpp


#include <iostream>
#include "dn.hpp"

template <typename T>
T f(T x) {
return (T(4) * x + T(3)) * x + T(2);
}

template <typename T>
T df(T x) {
return T(8) * x + T(3);
}

int main() {
// 手計算の結果 (確認用)
std::cout << '[' << f(2.0) << ',' << df(2.0) << ']' << std::endl;

// 自動微分の結果
std::cout << f(math::dual_number<>(2.0, 1.0)) << std::endl;

return 0;
}



出力結果

$ ./a.out

[24,19]
[24,19]

たしかに合ってる。


備考

ところで

f(a + b \epsilon) = f(a) + \left(\frac{d}{dx} f(a) \right) b \epsilon

ということに気がつくと、任意の実数関数 $f(x)$ について二重数の関数を求められる。

実数の関数
二重数の関数

$\sin(x)$
$\sin(a + b \epsilon) = \sin(a) + \cos(a) b \epsilon$

$\cos(x)$
$\cos(a + b \epsilon) = \cos(a) - \sin(a) b \epsilon$

$\exp(x)$
$\exp(a + b \epsilon) = \exp(a) (1 + b \epsilon)$

$\log(x)$
$\log(a + b \epsilon) = \log(a) + \frac{b}{a} \epsilon$

もちろん、これらを組み合わせた複雑な関数の微分値も計算できる。


main.cpp

#include <iostream>

#include <cmath>
#include "dn.hpp"

template <typename T>
inline math::dual_number<T> exp(math::dual_number<T> x) {
T const a = std::exp(x.a());
return math::dual_number<T>(a, a * x.b());
}

// シグモイド関数
template <typename T>
T f(T x) {
using namespace std;
return T(1) / (T(1) + exp(-x));
}

// シグモイド関数の導関数
template <typename T>
T df(T x) {
using namespace std;
T const s = f(x);
return s * (T(1) - s);
}

int main() {
// 手計算の結果 (確認用)
std::cout << '[' << f(2.0) << ',' << df(2.0) << ']' << std::endl;

// 自動微分の結果
std::cout << f(math::dual_number<>(2.0, 1.0)) << std::endl;

return 0;
}



実行結果

$ ./a.out

[0.880797,0.104994]
[0.880797,0.104994]


参考文献


  1. ゲームアプリの数学 Unityで学ぶ基礎からシェーダーまで