自動微分で遊ぼう

@rydotyosh / Ryogo Yoshimura

2016.07.23

自己紹介

CAD/CAM

https://www.cgsys.co.jp/g/products/CAM-TOOL/3dcam_cl.htm

自動微分とは

もくじ

微分

https://en.wikipedia.org/wiki/Derivative

偏微分

http://en.wikipedia.org/wiki/Gradient

利用場面

コンピュータで計算する微分

微分, 導関数

微分の記法

合成関数

合成関数の微分

例 | 合成関数の微分

多変数 | 合成関数の微分

多変数の例 | 合成関数の微分

計算グラフ

微分 | 計算グラフ

多変数 | 計算グラフ

プログラムで書かれた関数

std::vector<double> f(
    const std::vector<double> &x );

制御文 | プログラムで書かれた関数

if ( x > 0 )
  y = x;
else
  y = -x;

再代入 | プログラムで書かれた関数

自動微分

概要 1/2 | ボトムアップ型

概要 2/2 | ボトムアップ型

オーバーロード | ボトムアップ型

計算グラフ | ボトムアップ型

例 1/6 | ボトムアップ型

例 2/6 | ボトムアップ型

例 3/6 | ボトムアップ型

例 4/6 | ボトムアップ型

例 5/6 | ボトムアップ型

例 6/6 | ボトムアップ型

実装例 | ボトムアップ型

#include <iostream>
struct ad { double x, dx; };
ad operator+( const ad &f, const ad &g ) {
  return ad{ f.x + g.x, f.dx + g.dx };
}
ad operator*( const ad &f, const ad &g ) {
  return ad{ f.x * g.x, g.x * f.dx + f.x * g.dx };
}
int main() {
  ad x{ 2, 1 }, a{ 3, 0 }, b{ 4, 0 };
  ad y = ( x + a ) * ( b * x );
  std::cout << y.x << "," <<
               y.dx << std::endl;
}
// --> 40,28

1変数まとめ | ボトムアップ型

多変数 | ボトムアップ型

概要 1/3 | トップダウン型

概要 2/3 | トップダウン型

概要 3/3 | トップダウン型

計算グラフ | トップダウン型

例 1/11 | トップダウン型

例 2/11 | トップダウン型

例 3/11 | トップダウン型

例 4/11 | トップダウン型

例 5/11 | トップダウン型

例 6/11 | トップダウン型

例 7/11 | トップダウン型

例 8/11 | トップダウン型

例 9/11 | トップダウン型

例 10/11 | トップダウン型

例 11/11 | トップダウン型

まとめ | トップダウン型

自動微分ができるC++ライブラリ

Adept の紹介

式テンプレート

使ってみる

#include <iostream>
#include "adept/adept.h"
int main() {
  adept::Stack stack;    // 導関数の情報を格納するオブジェクト
  adept::adouble x = 2;  // 入力変数
  stack.new_recording(); // アルゴリズムの記録を開始
  adept::adouble y =
    ( x + 3 ) * ( 4 * x ); // アルゴリズムを実行
  y.set_gradient( 1.0 );   // 出力変数の勾配を設定
  stack.reverse();         // トップダウン型自動微分を実行
  std::cout << y.value() << "," << // 出力変数の値
               x.get_gradient() << // 導関数の値
               std::endl;
}
// --> 40,28

Boost と組み合わせてみる

※ 思いついた順

Adept x Boost.special_functions

Adept x Boost.special_functions

Adept x Boost.odeint

Adept x Boost.accumulators

Adept x Boost.ublas

Adept x Odeint

天体の動き | Adept x Odeint

初期条件 | Adept x Odeint

std::vector< double > init {
  /*x =*/ 3.0, /*y =*/ 0.0, // 座標
  /*u =*/ 0.3, /*v =*/ 0.2, // 速度
  /*m =*/ 1.0               // 太陽の質量
};

微分方程式 | Adept x Odeint

template< class T >
std::vector< T > simulate( const std::vector< T > &init ) {
  using namespace boost::numeric;
  typedef std::array< T, 4 > state_t;
  state_t x0 { init[ 0 ], init[ 1 ], init[ 2 ], init[ 3 ] };
  T m = init[ 4 ];
  auto system = [&]( const state_t &x,
                     state_t &dxdt, T /*t*/ ) { // 微分方程式を書く
    dxdt[ 0 ] = x[ 2 ]; // u`
    dxdt[ 1 ] = x[ 3 ]; // v`
    T r2 = x[ 0 ] * x[ 0 ] + x[ 1 ] * x[ 1 ];
    T r3 = pow( r2, 3.0 / 2.0 );
    dxdt[ 2 ] = -m * x[ 0 ] / r3; // x`
    dxdt[ 3 ] = -m * x[ 1 ] / r3; // y`
  };
  ...

軌道を求める | Adept x Odeint

template< class T >
std::vector< T > simulate( const std::vector< T > &init ) {
  ...
  auto stepper =
    odeint::controlled_runge_kutta<
      odeint::runge_kutta_dopri5< state_t, T > >();
  std::vector< T > orbit; // 軌道保存用
  auto observer = [&]( const state_t &x, T /*t*/ ) {
    orbit.push_back( x[ 0 ] );
    orbit.push_back( x[ 1 ] );
  };
  T t0 = 0.0, t1 = 15.0, dt = 0.1;
  odeint::integrate_const( // 軌道計算する
    stepper, system, x0, t0, t1, dt, observer );
  return orbit;
}

軌道を求める | Adept x Odeint

int main() {
  std::vector< double > init {
    /*x =*/ 3.0, /*y =*/ 0.0, // 座標
    /*u =*/ 0.3, /*v =*/ 0.2, // 速度
    /*m =*/ 1.0               // 太陽の質量
  };
  std::vector< double > orbit = simulate( init );
  ...

軌道を求める | Adept x Odeint

観測データっぽく | Adept x Odeint

  ...
  std::vector< double > observed = orbit;
  std::mt19937 gen( 0 );
  std::normal_distribution<> d( 0, 0.05 );
  for ( double &x : observed ) {
    x += d( gen );
  }
  ...

観測データっぽく | Adept x Odeint

テキトーな初期条件 | Adept x Odeint

  ...
  std::vector< double > init {
    /*x =*/ observed[0], /*y =*/ observed[1],
    /*u =*/ 0.4, /*v =*/ 0.3,
    /*m =*/ 1.3
  };
  ...

テキトーな初期条件 | Adept x Odeint

誤差 | Adept x Odeint

template< class T >
T error( const std::vector< T > &orbit,
         const std::vector< double > &observed ) {
  size_t n = orbit.size();
  T sum_sqr = 0.0;
  for ( size_t i = 0; i < n; ++i ) {
    T dv = orbit[ i ] - observed[ i ];
    sum_sqr += dv * dv;
  }
  return sum_sqr / n;
}

初期条件で微分 | Adept x Odeint

  ...
  size_t dim = init.size();
  ...
  adept::Stack stack;
  std::vector< adept::adouble > init_( dim ); // 入力変数
  boost::copy( init, init_.begin() );
  stack.new_recording();    // アルゴリズムの記録を開始
  std::vector< adept::adouble > orbit_ = simulate( init_ );
  adept::adouble err_ = error( orbit_, observed ); // 誤差値
  err_.set_gradient( 1.0 ); // 出力変数(誤差値)の勾配を設定
  stack.reverse();          // トップダウン型自動微分を実行
  std::vector< double > grad( dim );
  for ( size_t i = 0; i < dim; ++i  )
    grad[ i ] = init_[ i ].get_gradient(); // 偏導関数の値
  ...

初期条件を変化 | Adept x Odeint

  ...
  adam adam( dim );
  ...
  adam( grad, init );
  ...

Adam法 | Adept x Odeint

struct adam {
  ...
  void operator()( const std::vector< double > &dx,
                   std::vector< double > &x ) {
    double t = static_cast< double >( ++k );
    for ( size_t i = 0; i < dim; ++i ) {
      m1[ i ] = b1 * m1[ i ] + ( 1.0 - b1 ) * dx[ i ];
      m2[ i ] = b2 * m2[ i ] + ( 1.0 - b2 ) * dx[ i ] * dx[ i ];
      double c1 = m1[ i ] / ( 1.0 - pow( b1, t ) );
      double c2 = m2[ i ] / ( 1.0 - pow( b2, t ) );
      x[ i ] -= a * c1 / ( sqrt( c2 ) + e );
    }
  }
};

Adam法 | Adept x Odeint

struct adam {
  size_t dim, k;
  double a = 0.001, b1 = 0.9, b2 = 0.999, e = 1e-8;
  std::vector< double > m1, m2;
  adam( size_t dim ) : dim( dim ), k( 0 ),
    m1( dim ), m2( dim ) {}
  ...
}

初期条件を変化 | Adept x Odeint

元の軌道と比較 | Adept x Odeint

誤差の推移 | Adept x Odeint

感想

参考文献

Q and A