template <class T> class autodif {
public:
T v;
boost::numeric::ublas::vector<T> d;
autodif() {
v = 0.;
d.resize(0);
}
autodif(const T& x) {
v = x;
d.resize(0);
}
};
} // namespace kv
のような感じで、メンバとして値v(スカラー)と微分値d(ベクトル)を持っており、
例えば
一応ヘルパー関数として、
具体的な使用例は、test/test-autodif.ccを見るか、以下を見て欲しい。
#include <kv/autodif.hpp>
kv::autodif<double> f(const kv::autodif<double>& x) {
kv::autodif<double> tmp;
tmp = x * x - 1.;
return cos(tmp + sin(tmp));
}
int main()
{
kv::autodif<double> a1, a2;
double d1, d2;
a1 = kv::autodif<double>::init(1.5);
a2 = f(a1);
kv::autodif<double>::split(a2, d1, d2);
std::cout << d1 << "\n"; // f(1.5)
std::cout << d2 << "\n"; // f'(1.5)
}
これを実行すると、
-0.58768 -3.19266と表示される。 これは、f(x) = cos(x2 -1 + sin(x2 - 1)) に対して、f(1.5)を計算すると同時にf'(1.5)を計算している。
このように、initで初期化し、splitと値を微分値を分離する、という流れで使用する。
kv::autodif<double>::split(f(kv::autodif<double>::init(1.5)), d1, d2);
のように一気に書いても構わない。これでd1にf(1.5)が、d2にf'(1.5)が入る。
もしf(x)=xという関数なら(x.v,x.d(0))=(1.5,1)なので、 init関数はそのように初期化している。 また、splitはd1にa2.vを、d2にa2.d(0)をただコピーしてるだけ。
上の例だとfの定義が煩雑な上に他の型を流し込めなくなっているが、 当然template関数を使うことも出来る。 (autodif-sample2.cc)
#include <kv/autodif.hpp>
template <class T> T f(const T& x) {
T tmp;
tmp = x * x - 1.;
return cos(tmp + sin(tmp));
}
int main()
{
kv::autodif<double> a1, a2;
double d1, d2;
a1 = kv::autodif<double>::init(1.5);
a2 = f(a1);
kv::autodif<double>::split(a2, d1, d2);
std::cout << d1 << "\n"; // f(1.5)
std::cout << d2 << "\n"; // f'(1.5)
d1 = 1.5;
d2 = testfunc2(d1); // double can be also passed
std::cout << d2 << "\n"; // f(1.5)
}
これだとfの定義がすっきりする上に、全く同じfにdoubleを流しこむことも
可能になる。精度保証付き数値計算のプログラムを書くにはこれが出来ると
とても楽になる。
#include <boost/numeric/ublas/vector.hpp>
#include <boost/numeric/ublas/matrix.hpp>
#include <boost/numeric/ublas/io.hpp>
#include <kv/autodif.hpp>
namespace ub = boost::numeric::ublas;
template <class T> T f(const ub::vector<T>& x) {
return atan(x(0) * x(0) + 2. * x(1) * x(1) - 1.);
}
int main()
{
double d1;
kv::autodif<double> a1;
ub::vector<double> v1, v2;
ub::vector< kv::autodif<double> > va1;
v1.resize(2);
v1(0) = 5.; v1(1) = 6.;
va1 = kv::autodif<double>::init(v1);
a1 = f(va1);
kv::autodif<double>::split(a1, d1, v2);
std::cout << d1 << "\n"; // f(5, 6)
std::cout << v2 << "\n"; // gradient
}
initとsplitを使うのは1変数の場合とほぼ同じ。
initの引数とsplitの第3引数がベクトルになるのが違い。
initは、
va1(0).v = 5.; va1(0).d(0) = 1.; va1(0).d(1) = 0,;
va1(1).v = 6.; va1(1).d(0) = 0.; va1(1).d(1) = 1,;
のような動作をする。splitは、
y = a1.v; v2(0) = a1.d(0); v2(1) = a1.d(1)
のような動作をする。
実行すると、
1.56038 [2](0.00108495,0.00260388)と表示される。 それぞれ、
関数値 第1変数に関する微分、第2変数に関する微分である。
#include <boost/numeric/ublas/vector.hpp>
#include <boost/numeric/ublas/matrix.hpp>
#include <boost/numeric/ublas/io.hpp>
#include <kv/autodif.hpp>
namespace ub = boost::numeric::ublas;
template <class T> ub::vector<T> f(const ub::vector<T>& x) {
ub::vector<T> y(2);
y(0) = 2. * x(0) * x(0) * x(1) - 1.;
y(1) = x(0) + 0.5 * x(1) * x(1) - 2.;
return y;
}
int main()
{
ub::vector<double> v1, v2;
ub::vector< kv::autodif<double> > va1, va2;
ub::matrix<double> m;
v1.resize(2);
v1(0) = 5.; v1(1) = 6.;
va1 = kv::autodif<double>::init(v1);
va2 = f(va1);
kv::autodif<double>::split(va2, v2, m);
std::cout << v2 << "\n"; // f(5, 6)
std::cout << m << "\n"; // Jacobian matrix
}
前の例と大体同じであるが、splitの第2引数がベクトル、第3引数が行列になる。
initは、
va1(0).v = 5.; va1(0).d(0) = 1.; va1(0).d(1) = 0,;
va1(1).v = 6.; va1(1).d(0) = 0.; va1(1).d(1) = 1,;
のような動作をする。splitは、
v2(0) = va2(0).v; m(0)(0) = va2(0).d(0); m(0)(1) = va2(0).d(1);
v2(1) = va2(1).v; m(1)(0) = va2(1).d(0); m(1)(1) = va2(1).d(1);
のような動作をする。
実行すると、
[2](299,21) [2,2]((120,50),(1,6))のように表示される。それぞれ関数値(ベクトル)とヤコビ行列である。