関数の微分値を計算する方法はいくつかあります。
この記事では、二重数を用いて自動で一次微分を計算する方法を紹介します。
二重数とは、\(a + b\varepsilon\) で表される数値です。
ここで \(\varepsilon\) は、二乗すると 0 になる元を表します。
例えば、２つの二重数の掛け算は次のように計算されます。
$$
\begin{split}
(a + b\varepsilon) \times (c + d\varepsilon)
&= ac + ad\varepsilon + bc\varepsilon + bd{\varepsilon}^2 \\
&= ac + (ad + bd)\varepsilon
\end{split}
$$
二重数の原理を利用することで、複雑な数式の微分値を直観的に計算することができるようになります。
二重数のプログラム例
構造体と四則演算の定義
構造体Dualを定義します。
メンバの名称は難しいところですが、安直にaとbにしておきます。
#include <cmath>
#include <iostream>
struct Dual {
double a;
double b;
Dual() : a(0.0), b(0.0) {}
Dual(double a) : a(a), b(0.0) {}
Dual(double a, double b) : a(a), b(b) {}
Dual operator+() const { return *this; }
Dual operator-() const { return Dual(-a, -b); }
};
Dual operator+(Dual lhs, Dual rhs) {
return Dual{lhs.a + rhs.a, lhs.b + rhs.b};
}
Dual operator-(Dual lhs, Dual rhs) {
return lhs + (-rhs);
}
Dual operator*(Dual lhs, Dual rhs) {
return Dual{lhs.a * rhs.a, lhs.a * rhs.b + lhs.b * rhs.a};
}
Dual operator/(Dual lhs, Dual rhs) {
return Dual{lhs.a / rhs.a, (lhs.b * rhs.a - lhs.a * rhs.b) / (rhs.a * rhs.a)};
}
乗算の式は前述のとおり \( ac + (ad + bd)\varepsilon \) です。
除算については、分母の \(\varepsilon\) を消去するようにして、次のように計算できます。
$$
\begin{split}
\frac{a + b\varepsilon}{c + d\varepsilon}
&= \frac{(a + b\varepsilon)(c – d\varepsilon)}{(c + d\varepsilon)(c – d\varepsilon)} \\
&= \frac{ac – ad\varepsilon + bc\varepsilon -bd\varepsilon}{c^2 – d^2{\varepsilon}^2} \\
&= \frac{ac + (bc – ad)\varepsilon}{c^2} \\
&= \frac{a}{c} + \frac{bc – ad}{c^2}\varepsilon
\end{split}
$$
その他、加算代入演算子 += なども定義する方が利便性が高いですが、この記事では省略します。
初等関数の定義
二重数に対する初等関数もいくつか定義します。
初等関数については、テイラー展開の原理によって \(f(a + b\varepsilon) = f(a) + bf'(a)\varepsilon\) として計算することができます。
Dual exp(Dual dual) {
double a = std::exp(dual.a);
return Dual{a, dual.b * a};
}
Dual log(Dual dual) {
return Dual{std::log(dual.a), dual.b / dual.a};
}
Dual sqrt(Dual dual) {
double a = std::sqrt(dual.a);
return Dual{a, -0.5 * dual.b / a};
}
Dual sin(Dual dual) {
return Dual{std::sin(dual.a), dual.b * std::cos(dual.a)};
}
Dual cos(Dual dual) {
return Dual{std::cos(dual.a), -dual.b * std::sin(dual.a)};
}上記に定義した初等関数については、次のような数式となります。
$$
\begin{split}
e^{a + b\varepsilon} &= e^a + be^a\varepsilon \\
\ln(a + b\varepsilon) &= \ln{a} + \frac{b}{a}\varepsilon \\
\sqrt{a + b\varepsilon} &= \sqrt{a} – \frac{b}{2\sqrt{a}}\varepsilon \\
\sin{a + b\varepsilon} &= \sin{a} + b\varepsilon\cos{a} \\
\cos{a + b\varepsilon} &= \cos{a} – b\varepsilon\sin{a} \\
\end{split}
$$
初等関数は他にもありますが、ここでは省略します。
標準出力の定義
利便性のため、Dual構造体に対する標準出力も定義しておきます。
std::ostream &operator<<(std::ostream &os, Dual dual) {
os << "(" << dual.a << ", " << dual.b << ")";
return os;
}自動微分の計算例
ここまで準備すると、自動微分を計算することができます。
以下のコードは、関数 \(f(x) = x\sin{x}\) に対して \(f'(2))\) を計算するためのサンプルプログラムです。
int main() {
Dual x{2.0, 1.0};
Dual y = x * sin(x);
std::cout << y << std::endl;
// (1.81859, 0.0770038)
return 0;
}まず、微分したい変数 \(x\) を、\(b=1\) である二重数として定義します。
\(a\) には \(x\) に代入したい値 \(2\) を指定します。
つまり、ここでは \(x = 2 + \varepsilon\) と定義します。
そのあとは、二重数 \(x\) を用いて \(x\sin{x}\) をそのまま計算するだけです。
計算結果には \(f(2) + f'(2)\varepsilon\) の二重数が格納されることになります。
実行結果としては、\(f(2) = 1.81859, f'(2) = 0.0770038\) となります。
後者の数値が求めたい微分値となります。
なお、前者の方に出力されているように、ついでに関数 \(f(2)\) そのものの計算結果も取得されることとなります。
微分式を求めて計算する場合
検算のため、微分式を直接計算するプログラムを確認します。
関数 \(f(x) = x\sin{x}\) を実直に微分すると \(f'(x) = \sin{x} + x\cos{x}\) となります。
これに 2 を代入すると、次のように 0.0770038 が得られます。
int main() {
double x = 2.0;
std::cout << std::sin(x) + x * std::cos(x) << std::endl;
// 0.0770038
return 0;
}自動微分の結果と一致することが確認できました。
自動微分の応用例
自動微分の応用例として、ニュートン法への利用があります。
以下の記事で解説します。
コメント