トップ サイトマップ
関連情報
ホーム > 製品&サービス > コンサルティングサービス > 自動微分ソリューション > 自動微分ツール:dco/c++ユーザーガイド(トライアル・バージョン)

自動微分ツール:dco/c++ユーザーガイド(トライアル・バージョン)

バージョン3.1、文書バージョン1.1

Uwe Naumann,
The Numerical Algorithms Group Ltd. (NAG),
and
LuFG Informatik 12
Software and Tools for Computational Engineering (STCE)
Department of Computer Science
RWTH Aachen University

The Numerical Algorithms Group Limited, 2015 All rights reserved.


[2016年4月掲載]

序文

dco/c++は、C++の演算子オーバーローディングによる一次微分および高次微分のタンジェントモード(tangent mode)アジョイントモード(adjoint mode)の自動微分(Algorithmic DifferentiationあるいはAutomatic Differentiationとも呼ばれます)の柔軟かつ効果的な実装です。これは、直観的かつパワフルなアプリケーション開発インターフェイスと、C++テンプレートによるキャッシュ最適化された内部表現の組合せにより実現されています。dco/c++は、例えば膨大な量のパラメータ較正や形状最適化等の場面において、様々な数値計算シミュレーションへの適用に成功してきました。

ここで、$n$個の不確定なパラメーター$x\in\mathbb{R}^n$(例えば、マーケットボラティリティ)に依存するスカラー値(例えば、金融商品の価格)を日常的に数値シミュレーションする場合を考えましょう。そのプライシング関数$y=f(x)$計算のための既存のC++プログラムが、あなたのコンピュータ上で一回の計算に1分掛かるとします。さらにこの関数$y$の他に、全$x$に関する$y$の微分を求めるとします。この$n$個の微分を求めるコストは、有限差分近似においては$f$の計算の$O(n)$倍となります。また、これらの精度は、切り捨ておよび浮動小数点演算中のキャンセレーションや丸め誤差などの数値演算上の影響を受けます。さらにもし、例えば$n$=1000とした場合、グラジェントの近似計算は少なくとも1000分(16時間以上)掛かるでしょう。アジョイントモードでのdco/c++は、同等な精度のグラジェント計算を、20分(10分以下のケースが多い)以下に抑えることが期待されます。

ADは、マニュアルで実装することも可能です。その時与えられた任意の目的関数の実装に対して、AD専門家がそのアジョイントモード・バージョンを作成し、ツールベースの実装に対して実行性能が上回ることも可能です。しかしこの作業は冗長で、エラーを生じやすく、デバッグも極めて困難です。さらに重要なことは、維持性や保守性といった現代のソフトウェア工学の基礎的な要求を満たしません。オリジナルソースの変更があれば、アジョイントモードの対応する変更も必要になります。両方のソースコードを矛盾なく維持していくことは極めて困難な作業となるでしょう。

dco/c++は、大規模な実業務のシミュレーションコードに適用されることを想定して設計されています。利用可能なメモリーリソース範囲内で、チェックポイント、陰関数のシンボリックな微分、および既存のアジョイントモード・ソースコードを一つのdco/c++アジョイントコードへ統合することを組み合わせ、ユーザー・アプリケーションのロバストで効果的なアジョイントモードを得ることが出来ます。ユーザーには、dco/c++アジョイントコードとリンクするように最適化された、専門分野に特化したより高いレベルの組み込み関数のライブラリーを構築することが奨励されます。dco/c++アジョイントコードのコールバック・インターフェイスは、GPUのようなアクセラレーターの利用や、既存の高度に最適化された数値計算ライブラリーとの併用を容易にします。NAGライブラリーに含まれる様々な手法は、別のNAG ADライブラリーとのリンクによって組込み関数としてサポートされます。


免責事項:このユーザーマニュアルは幾つかの例題を用いて作成されています。これらの多くは自明なものですが、それ以外は、dco/c++の文法の背後に在るより深い意味解釈が必要となるものです。これらについては本文書の更新版で提供される予定であり、現状内容は概要版です。


本文書はADについて理解の進んだユーザーを対象にしています。ADのコミュニティーにより、広範な文献目録と共に、ポータルサイトwww.autodiff.orgが運営されています。また、本文書は、文献[Naumann, U. (2012). The Art of Differentiating Computer Programs. An Introduction to Algorithmic Differentiation. Number 24 in Software, Environments, and Tools. SIAM.]の記法に沿って記述をしています。


もくじ

1. 機能概要

基本データ型$D$(例えば、単精度や高次精度の浮動小数点、時間間隔、convex/concave緩和、ベクトル、集合等)を用いたコンピュータープログラムとして、以下の多変数ベクトル関数の実装を考える。
(ここで、$f$内の演算は$D$の変数に対して(基本関数(elemental function)のオーバーローディングを通して)完全に定義されていると仮定する)

$$ f : D^{n} \times D^{n'} \rightarrow D^m \times D^{m'} : (\mathbf{y}, \mathbf{y}') = f(\mathbf{x},\mathbf{x}') $$

$n$個の(独立な)アクティブ(active)入力と$n'$個のパッシブ(passive)入力が、$m$個の(独立でない)アクティブ入力と$m'$個のパッシブ出力へ写像される。この実装では、対象とするすべての点について$k$回微分可能と仮定され、以下のヤコビアンとヘッシアンおよび高次の微分テンソル($k$≥2)の存在と有界性が与えられる。

ヤコビアン:

$$ \nabla f = \nabla_\mathbf{x} f(\mathbf{x},\mathbf{x}') \equiv \frac{\partial \mathbf{y}}{\partial \mathbf{x}}(\mathbf{x},\mathbf{x}') \in D^{m \times n}, $$

ヘッシアン:

$$ \nabla^2 f = \nabla^2 f(\mathbf{x},\mathbf{x}') \equiv \frac{\partial^2 \mathbf{y}}{\partial \mathbf{x}^2}(\mathbf{x},\mathbf{x}') \in D^{m \times n \times n}, $$

高次微分テンソル:

$$ \nabla^k f = \nabla^k f(\mathbf{x},\mathbf{x}') \equiv \frac{\partial^k \mathbf{y}}{\partial \mathbf{x}^k}(\mathbf{x},\mathbf{x}') \in D^{m \times \overset{k}{\underset{l=1}{\times}} n}. $$

ここで、

$$ \nabla^k f=\left(\left [\nabla^k f \right ]_i^{j_1,\ldots,j_k} \right)_{i=0,\ldots,m-1}^{j_1,\ldots,j_k=0,\ldots,n-1}, $$

であり、記号$*$はインデックス全体を示しものとする。例えば、$\left [\nabla f \right ]_i^*$はヤコビアンの$i$番目の行を意味し、$\left [\nabla f \right ]_*^j$は$j$番目の列を意味する。

1.1 一次タンジェントモード

1.1.1 ジェネリック型一次スカラータンジェントモード

ジェネリック型一次スカラータンジェントモード(tangent mode)は、任意の基本データ型$D$を用いたジェネリック型一次スカラータンジェントデータ型を規定することにより、ヤコビアンとベクトル$\mathbf{x}^{(1)}$$\in D^n$の積を計算する。

$$ \mathbf{y}^{(1)}=\nabla f \cdot \mathbf{x}^{(1)} \in D^m. $$

1.1.2 ジェネリック型一次ベクトルタンジェントモード

ジェネリック型一次ベクトルタンジェントモードは、任意の基本データ型$D$を用いたジェネリック型一次ベクトルタンジェントデータ型を規定することにより、ヤコビアンと行列$\mathbf{X}^{(1)} \in D^{n \times l}$の積を計算する

$$ Y^{(1)}=\nabla f \cdot X^{(1)} \in D^{m \times l}. $$

1.1.3 式テンプレートの利用

式テンプレートを用いて、個々の割当てレベルにおいて静的に最適化されたグラジェントコードを生成可能である。これはベクトルタンジェントモードで効果を発揮する。

1.2 一次アジョイントモード

1.2.1 ジェネリック型一次スカラーアジョイントモード

ジェネリック型一次スカラーアジョイントモード(adjoint mode)は、任意の基本データ型$D$を用いたジェネリック型一次スカラーアジョイントデータ型を規定することにより、転置ヤコビアンとベクトル$\mathbf{y}_{(1)}$$\in D^m$の積を計算する

$$ \mathbf{x}_{(1)}=\nabla f^T \cdot \mathbf{y}_{(1)} \in D^n. $$

1.2.2 ジェネリック型一次ベクトルアジョイントモード

ジェネリック型一次ベクトルアジョイントモードは、任意の基本データ型$D$を用いたジェネリック型一次ベクトルアジョイントデータ型を規定することにより、転置ヤコビアンと行列$\mathbf{Y}_{(1)} \in D^{m \times l}$の積を計算する

$$ X_{(1)}=\nabla f^T \cdot Y_{(1)} \in D^{n \times l}. $$

1.2.3 グローバル “blob” テープ

指定されたサイズのメモリーが確保され、テープ(tape)内容の保存に用いられる。境界値チェックはしない。

1.2.4 グローバル “chunk” テープ

このテープは物理的なメモリー限界まで可能である。

1.2.5 (スレッド)ローカルテープ

例えば、スレッド化されたアジョイントコード実装用に、“blob”テープあるいは“chunk”テープを複数確保可能である。

1.2.6 式テンプレートの利用

式テンプレートを用いて、個々の割当てレベルにおいて静的に最適化されたグラジェントコードを生成可能。これによりテープのメモリー要求を減少させることが出来る。

1.2.7 テープの繰返し計算

テープは与えられた地点で記録され、繰返し再生可能である。

1.2.8 アジョイント・コールバック・インターフェイス

アジョイント・コールバック・インターフェイスは以下の機能をサポートする。

  1. チェックポイント機能
  2. 記号的アジョイント手法
  3. タンジェントモードとアジョイントモードの組合せ
  4. ローカルな微分テンソルの式テンプレート
  5. 領域指定の高次基本関数群の生成

1.2.9 ユーザー指定ローカルヤコビアン・インターフェイス

他所で生成した偏微分式テンプレートを、テープの再生系列内へ直接挿入することが可能である。

1.3 二次および高次タンジェントモード

1.3.1 ジェネリック型二次スカラータンジェントモード

非微分形の基本データ型$D$を用いたジェネリック型一次スカラータンジェントデータ型のジェネリック型一次スカラータンジェントデータ型のインスタンスは、二次スカラータンジェントデータ型である。これは、ヘッシアン$\nabla^2 f \in D^{m \times n \times n}$の2つの長さ$n$次元領域におけるスカラー射影を以下のように計算する

$$ \mathbf{y}^{(1,2)}=\langle \nabla^2 f,\mathbf{x}^{(1)},\mathbf{x}^{(2)}\rangle \equiv \left(\mathbf{x}^{{(1)}^T} \cdot \left [\nabla^2 f \right ]_i^{*,*} \cdot \mathbf{x}^{(2)} \right)_{i=0,\ldots,m-1} \in D^m. $$

ここで、$\mathbf{x}^{(1)},\mathbf{x}^{(2)} \in D^n$は任意の基本データ型$D$である。$\mathbf{x}^{(1)},\mathbf{x}^{(2)}$を$D$のカルテシアン基底ベクトル全体に渡らせた際の全ヘッシアンの計算コストは$O(n^2) \cdot \text{Cost}(f)$である。

1.3.2 ジェネリック型二次ベクトルタンジェントモード

非微分形の基本データ型$D$を用いたジェネリック型一次ベクトルタンジェントデータ型のジェネリック型一次ベクトルタンジェントデータ型のインスタンスは、二次スカラータンジェントデータ型である。これは、ヘッシアン$\nabla^2 f \in D^{m \times n \times n}$の2つの長さ$n$次元領域におけるベクトル射影を以下のように計算する

$$ Y^{(1,2)}=\langle \nabla^2 f,X^{(1)},X^{(2)}\rangle \equiv \left(X^{{(1)}^T} \cdot \left [\nabla^2 f \right ]_i^{*,*} \cdot X^{(2)} \right)_{i=0,\ldots,m-1} \in D^{m \times l_1 \times l_2}. $$

ここで、$X^{(1)} \in D^{n \times l_1},$ $X^{(2)} \in D^{n \times l_2}$は任意の基本データ型$D$である。$\mathbf{X}^{(1)},\mathbf{X}^{(2)}$を$D$の単位元全体に渡らせた際の$D$上の全ヘッシアンの計算コストは$O(n^2) \cdot \text{Cost}(f)$になる。

1.3.3 他のジェネリック型二次タンジェントモード

非微分形の基本データ型$D$を用いたジェネリック型一次スカラータンジェントデータ型のジェネリック型一次ベクトルタンジェントデータ型のインスタンスは、二次タンジェントデータ型である。同様に、非微分形の基本データ型$D$を用いたジェネリック型一次ベクトルタンジェントデータ型のジェネリック型一次スカラータンジェントデータ型のインスタンスは、二次タンジェントデータ型である。

1.3.4 ジェネリック型三次および高次タンジェントモード

$k$次のタンジェント型を用いたタンジェント型のインスタンスは、$(k+1)次タンジェント型である。

1.4 二次および高次アジョイントモード

1.4.1 ジェネリック型二次スカラーアジョイントモード

非微分形の基本データ型$D$を用いたジェネリック型一次スカラータンジェントデータ型のジェネリック型一次スカラーアジョイントデータ型のインスタンスは、二次スカラーアジョイントデータ型である。これは、長さ$n$の二次元領域におけるヘッシアン$\nabla^2 f \in D^{m \times n \times n}$のスカラー射影を以下のように計算する

$$ \mathbf{x}_{(1)}^{(2)}=\langle \mathbf{y}_{(1)},\nabla^2 f,\mathbf{x}^{(2)}\rangle \equiv \left(\mathbf{y}_{(1)}^T \cdot \left [\nabla^2 f \right ]_*^{j,*} \cdot \mathbf{x}^{(2)} \right)_{j=0,\ldots,n-1} \in D^n. $$

ここで、$\mathbf{y}_{(1)} \in D^m$ and $\mathbf{x}^{(2)} \in D^n$は任意の基本データ型$D$である。$\mathbf{y}_{(1)}$および$\mathbf{x}^{(2)}$を$D$のカルテシアン基底ベクトル全体に渡らせた際の全ヘッシアンの計算コストは$O(m\cdot n) \cdot \text{Cost}(f)$となる。

1.4.2 ジェネリック型二次ベクトルアジョイントモード

非微分形の基本データ型$D$を用いたジェネリック型一次ベクトルタンジェントデータ型のジェネリック型一次ベクトルアジョイントデータ型のインスタンスは、二次ベクトルアジョイントデータ型である。これは、長さ$n$の二次元領域におけるヘッシアン$\nabla^2 f \in D^{m \times n \times n}$のベクトル射影を以下のように計算する

$$ X_{(1)}^{(2)}=\langle Y_{(1)},\nabla^2 f,X^{(2)}\rangle \equiv \left(Y_{(1)}^T \cdot \left [\nabla^2 f \right ]_*^{j,*} \cdot X^{(2)} \right)_{j=0,\ldots,n-1} \in D^{l_1 \times n \times l_2}. $$

ここで、$Y_{(1)} \in D^{m \times l_1}$および$X^{(2)} \in D^{n \times l_2}$は任意の基本データ型$D$である。$Y_{(1)}$および$X^{(2)}$を$D$の単位元全体に渡らせた際の$D$上の全ヘッシアンの計算コストは$O(m\cdot n) \cdot \text{Cost}(f)$となる。

1.4.3 他のジェネリック型二次アジョイントモード

非微分形の基本データ型$D$を用いたジェネリック型一次スカラータンジェントデータ型のジェネリック型一次ベクトルアジョイントデータ型のインスタンスは、二次アジョイントデータ型である。同様に、非微分形の基本データ型$D$を用いたジェネリック型一次ベクトルタンジェントデータ型のジェネリック型一次スカラーアジョイントデータ型のインスタンスは、二次アジョイントデータ型である。

二次元領域のヘッシアンの対称性から、非微分形の基本データ型$D$を用いたジェネリック型一次アジョイントデータ型のジェネリック型一次タンジェントデータ型のインスタンス化により、以下の二次アジョイントデータ型が得られる

$$ \langle \mathbf{y}_{(2)}^{(1)},\nabla^2 f,\mathbf{x}^{(1)}\rangle \equiv \left(\mathbf{y}_{(2)}^{(1)^T} \cdot \left [\nabla^2 f \right ]_*^{j,*} \cdot \mathbf{x}^{(1)} \right)_{j=0,\ldots,n-1} \in D^n. $$

同様に、非微分形の基本データ型$D$を用いたジェネリック型一次アジョイントデータ型のジェネリック型一次アジョイントデータ型のインスタンス化により、以下の二次アジョイントデータ型が得られる

$$ \langle \mathbf{x}_{(1,2)},\mathbf{y}_{(1)},\nabla^2 f\rangle \equiv \left(\mathbf{y}_{(1)}^T \cdot \left [\nabla^2 f \right ]_*^{j,*} \cdot \mathbf{x}_{(1,2)} \right)_{j=0,\ldots,n-1} \in D^n. $$

1.4.4 ジェネリック型三次および高次アジョイントモード

$k$次のアジョイント型を用いたタンジェント型のインスタンスは、$(k+1)$次アジョイント型である。同様に、$k$次のタンジェントあるいはアジョイント型を用いたアジョイント型のインスタンスは、$(k+1)$次アジョイント型である。

2. 一次タンジェントモード

2.1 目的

dco/c++のデータ型gt1s<DCO_BASE_TYPE>::typeにより、一次スカラータンジェントモードが実装される。

多変数関数$f:\mathbb{R}^n \rightarrow \mathbb{R}^m, \mathbf{y}=f(\mathbf{x})$の一次タンジェントモードバージョン:

$$ \left( \begin{array}{c} \mathbf{y} \\ \mathbf{y}^{(1)} \\ \end{array} \right) :=f^{(1)}(\mathbf{x},\mathbf{x}^{(1)}) $$

これは、多変数関数$f$へのdco/c++の適用により生成され、以下の計算を行う

\begin{eqnarray} \mathbf{y} &:=& f(\mathbf{x}) \\ \mathbf{y}^{(1)} &:=& \langle \frac{\partial f}{\partial \mathbf{x}},\mathbf{x}^{(1)} \rangle \equiv \frac{\partial f}{\partial \mathbf{x}} \cdot \mathbf{x}^{(1)}. \tag{2.1} \end{eqnarray}

2.2 例題

2.2.1 例題ソース

ここで、以下の関数$f:\mathbb{R}^4 \rightarrow \mathbb{R}^2$の実装を考える:


 1 template<typename T>
 2 void f(const vector<T>& x, vector<T>& y) {
 3   T v = tan(x[2] * x[3]);
 4   T w = x[1] - v;
 5   y[0] = x[0] * v / w;
 6   y[1] = y[0] * x[1];
 7 }

そのドライバーは、式(2.1)を$\mathbf{x}\hat{=}$xv,$\mathbf{x}^{(1)}\hat{=}$xt,$\mathbf{y}\hat{=}$yv,$\mathbf{y}^{(1)}\hat{=}$ytとして計算する。


 1 #include <vector>
 2 using namespace std;
 3
 4 #include "dco.hpp"
 5 using namespace dco;
 6
 7 typedef double DCO_BASE_TYPE;
 8 typedef gt1s<DCO_BASE_TYPE> DCO_MODE;
 9 typedef DCO_MODE::type DCO_TYPE;
10
11 #include "f.hpp"
12
13 void driver(
14   const vector<double>& xv, const vector<double>& xt,
15   vector<double>& yv, vector<double>& yt
16 ) {
17   const int n=xv.size(), m=yv.size();
18   vector<DCO_TYPE> x(n), y(m);
19   for (int i=0;i<n;i++) { value(x[i])=xv[i]; derivative(x[i])=xt[i]; }
20   f(x,y);
21   for (int i=0;i<m;i++) { yv[i]=value(y[i]); yt[i]=derivative(y[i]); }
22 }

メインプログラムは、式(2.1)の右辺の変数全てを初期化してから、ドライバーを呼出し、結果を出力する。


 1 #include<iostream>
 2 #include<vector>
 3 using namespace std;
 4
 5 #include "driver.hpp"
 6
 7   int main() {
 8   const int m=2, n=4; cout.precision(15);
 9   vector<double> xv(n), xt(n), yv(m), yt(m);
10   for (int i=0;i<n;i++) { xv[i]=1; xt[i]=1; }
11     driver(xv,xt,yv,yt);
12   for (int i=0;i<m;i++)
13     cout << "y[" << i << "]=" << yv[i] << endl;
14   for (int i=0;i<m;i++)
15     cout << "y^{(1)}[" << i << "]=" << yt[i] << endl;
16   return 0;
17 }

2.2.2 結果

以下の結果が生成される:


y[0]=-2.79401891249195
y[1]=-2.79401891249195
y^{(1)}[0]=14.2435494001203
y^{(1)}[1]=11.4495304876283

2.3 dco/c++の機能

2.3.1 dco.hpp

dco/c++のヘッダーファイル。

2.3.2 DCO_BASE_TYPE

ジェネリック型ADモードの変数インスタンス型。

2.3.3 gt1s<DCO_BASE_TYPE>

ジェネリック型一次タンジェントモード。

2.3.4 DCO_MODE

ジェネリックADモード。例えば、typedef gt1s<DCO_BASE_TYPE> DCO_MODEとする。

2.3.5 gt1s<DCO_BASE_TYPE>::type

ジェネリック型一次タンジェント型。

2.3.6 DCO_TYPE

ジェネリックAD型。例えば、typedef gt1s<DCO_BASE_TYPE> DCO_TYPEとする。

2.3.7 value

引数の値の参照を返す。

2.3.8 derivative

引数の微分値の参照を返す。

3. 一次アジョイントモード

3.1 目的

dco/c++のデータ型ga1s<DCO_BASE_TYPE>::typeにより、一つのグローバル・テープを用いた一次スカラーアジョイントモードが実装される。
多変数関数$f:\mathbb{R}^n \rightarrow \mathbb{R}^m, \mathbf{y}=f(\mathbf{x})$の一次アジョイントモードバージョン:

$$ \left( \begin{array}{c} \mathbf{y} \\ \mathbf{x}_{(1)} \\ \end{array} \right) :=f_{(1)}(\mathbf{x},\mathbf{x}_{(1)},\mathbf{y},\mathbf{y}_{(1)}) $$

これは、多変数関数$f$へのdco/c++の適用により生成され、以下の計算を行う

\begin{eqnarray} \mathbf{y} &:=& f(\mathbf{x}) \\ \mathbf{x}_{(1)} &:=& \mathbf{x}_{(1)}+\langle \mathbf{y}_{(1)},\frac{\partial f}{\partial \mathbf{x}} \rangle \equiv \mathbf{x}_{(1)}+\left(\frac{\partial f}{\partial \mathbf{x}} \right)^T \cdot \mathbf{y}_{(1)}. \tag{3.1} \end{eqnarray}

3.2 例題

3.2.1 例題ソース

ここで、第2章と同じ関数$f:\mathbb{R}^4 \rightarrow \mathbb{R}^2$に対して、以下のドライバーで$\mathbf{x}\hat{=}$xv,$\mathbf{x}_{(1)}\hat{=}$xa,$\mathbf{y}\hat{=}$yv,$\mathbf{y}_{(1)}\hat{=}$yaとして式(3.1)を計算する。


 1 #include <vector>
 2 using namespace std;
 3
 4 #include "dco.hpp"
 5 using namespace dco;
 6
 7 typedef double DCO_BASE_TYPE;
 8 typedef ga1s<DCO_BASE_TYPE> DCO_MODE;
 9 typedef DCO_MODE::type DCO_TYPE;
10 typedef DCO_MODE::tape_t DCO_TAPE_TYPE;
11
12 #include "../gt1s/f.hpp"
13
14 void driver(
15   const vector<double>& xv, const vector<double>& xa,
16   vector<double>& yv, vector<double>& ya
17 ) {
18   DCO_MODE::global_tape=DCO_TAPE_TYPE::create();
19   int n=xv.size(), m=yv.size();
20   vector x(n), y(m);
21   for (int i=0;i<n;i++) {
22     x[i]=xv[i];
23     DCO_MODE::global_tape->register_variable(x[i]);
24   }
25   f(x,y);
26   for (int i=0;i<m;i++) {
27     DCO_MODE::global_tape->register_output_variable(y[i]);
28     yv[i]=value(y[i]); derivative(y[i])=ya[i];
29   }
30   for (int i=0;i<n;i++) derivative(x[i])=xa[i];
31   DCO_MODE::global_tape->interpret_adjoint();
32   for (int i=0;i<n;i++) xa[i]=derivative(x[i]);
33   for (int i=0;i<m;i++) ya[i]=derivative(y[i]);
34   DCO_TAPE_TYPE::remove(DCO_MODE::global_tape);
35 }

メインプログラムは、式(3.1)の右辺の変数全てを初期化してから、ドライバーを呼出し、結果を出力する。


 1 #include<iostream>
 2 #include<vector>
 3 using namespace std;
 4
 5 #include "driver.hpp"
 6
 7   int main() {
 8   const int n=4, m=2; cout.precision(15);
 9   vector<double> xv(n), xa(n), yv(m), ya(m);
10   for (int i=0;i<n;i++) { xv[i]=1; xa[i]=1; }
11   for (int i=0;i<m;i++) ya[i]=1; }
12     driver(xv,xt,yv,yt);
13   for (int i=0;i<m;i++)
14     cout << "y[" << i << "]=" << yv[i] << endl;
15   for (int i=0;i<n;i++)
16     cout << "x_{(1)}[" << i << "]=" << xa[i] << endl;
17   return 0;
18 }

3.2.2 結果

以下の結果が生成される:


y[0]=-2.79401891249195
y[1]=-2.79401891249195
x_{(1)}[0]=-4.5880378249839
x_{(1)}[1]=-11.8190644542335
x_{(1)}[2]=23.050091083483
x_{(1)}[3]=23.050091083483

3.3 dco/c++の機能

3.3.1 ga1s<DCO_BASE_TYPE>

ジェネリック型一次スカラーアジョイントモード。

3.3.2 ga1s<DCO_BASE_TYPE>::type

ジェネリック型一次スカラーアジョイント型。

3.3.3 DCO_TAPE_TYPE

テープ型のDCO_MODE

3.3.4 DCO_MODE::global_tape

グローバルテープポインタのDCO_MODE

3.3.5 DCO_TAPE_TYPE::create

テープクリエイターはテープポインターをDCO_TAPE_TYPE*として返す。

3.3.6 DCO_TAPE_TYPE::register_variable

テープ内の引数(独立変数)エントリーを生成する。

3.3.7 DCO_TAPE_TYPE::register_output_variable

テープ内の従属変数を引数としてマークする。

3.3.8 DCO_TAPE_TYPE::interpret_adjoint

全テープのアジョイント再生。

3.3.9 DCO_TAPE_TYPE::remove

グローバルテープの削除。

4. 二次タンジェントモード

4.1 目的

多変数関数$f:\mathbb{R}^n \rightarrow \mathbb{R}^m, \mathbf{y}=f(\mathbf{x})$の一次タンジェントモードバージョン:

$$ \left( \begin{array}{c} \mathbf{y} \\ \mathbf{y}^{(2)} \\ \mathbf{y}^{(1)} \\ \mathbf{y}^{(1,2)} \\ \end{array} \right) :=f^{(1,2)}(\mathbf{x},\mathbf{x}^{(2)} ,\mathbf{x}^{(1)} ,\mathbf{x}^{(1,2)}) $$

これは、多変数関数$f$へのdco/c++の適用により生成され、以下の計算を行う

\begin{eqnarray} \mathbf{y} &:=& f(\mathbf{x}) \\ \mathbf{y}^{(2)} &:=& \langle \frac{\partial f}{\partial \mathbf{x}},\mathbf{x}^{(2)} \rangle \\ \mathbf{y}^{(1)} &:=& \langle \frac{\partial f}{\partial \mathbf{x}},\mathbf{x}^{(1)} \rangle \\ \mathbf{y}^{(1,2)} &:=& \langle \frac{\partial^2 f}{\partial^2 \mathbf{x}},\mathbf{x}^{(1)},\mathbf{x}^{(2)} \rangle+\langle \frac{\partial f}{\partial \mathbf{x}},\mathbf{x}^{(1,2)} \rangle. \tag{4.1} \end{eqnarray}

4.2 例題

4.2.1 例題ソース

ここで、第2章と同じ関数$f:\mathbb{R}^4 \rightarrow \mathbb{R}^2$に対して、以下のドライバーで$\mathbf{x}\hat{=}$xv,$\mathbf{x}^{(1)}\hat{=}$xt1,$\mathbf{x}^{(2)}\hat{=}$xt2,$\mathbf{x}^{(1,2)}\hat{=}$xt1t2, $\mathbf{y}\hat{=}$yv,$\mathbf{y}^{(1)}\hat{=}$yt1,$\mathbf{y}^{(2)}\hat{=}$yt2,$\mathbf{y}^{(1,2)}\hat{=}$yt1t2 として式(4.1)を計算する。


 1 #include <iostream>
 2 using namespace std;
 3
 4 #include "dco.hpp"
 5 using namespace dco;
 6 typedef gt1s%lt;double> DCO_BASE_MODE;
 7 typedef DCO_BASE_MODE::type DCO_BASE_TYPE;
 8 typedef gt1s<DCO_BASE_TYPE> DCO_MODE;
 9 typedef DCO_MODE::type DCO_TYPE;
10
11 #include "../gt1s/f.hpp"
12
13 void driver(
14   const vector<double>& xv,
15   const vector<double>& xt1,
16   const vector<double>& xt2,
17   const vector<double>& xt1t2,
18   vector<double>& yv, 
19   vector<double>& yt1,
20   vector<double>& yt2,
21   vector<double>& yt1t2,
22 ) {
23   const int n=xv.size(), m=yv.size();
24   vector<DCO_TYPE> x(n), y(m);
25   for (int i=0;i<n;i++) {
26     value(value(x[i]))=xv[i];
27     derivative(value(x[i]))=xt1[i];
28     value(derivative(x[i]))=xt2[i];
29     derivative(derivative(x[i]))=xt1t2[i];
30   }
31   f(x,y);
32   for (int i=0;i<m;i++) {
33     yv[i]=passive_value(y[i]);
34     yt1[i]=derivative(value(y[i]));
35     yt2[i]=value(derivative(y[i]));
36     yt1t2[i]=derivative(derivative(y[i]));
37   }
38 }

メインプログラムは、式(4.1)の右辺の変数全てを初期化してから、ドライバーを呼出し、結果を出力する。


 1 #include<iostream>
 2 #include<vector>
 3 using namespace std;
 4
 5 #include "driver.hpp"
 6
 7   int main() {
 8   const int n=4, m=2; cout.precision(15);
 9   vector<double> xv(n), xt1(n), xt2(n), xt1t2(n);
10   vector<double> yv(m), yt1(m), yt2(m), yt1t2(m);
11   for (int i=0;i<n;i++) { xv[i]=1; xt1[i]=1; xt2[i]=1; xt1t2[i]=1; }
12     driver(xv,xt1,xt2,xt1t2,yv,yt1,yt2,yt1t2);
13   for (int i=0;i<m;i++)
14     cout << "y[" << i << "]=" << yv[i] << endl;
15   for (int i=0;i<m;i++)
16     cout << "y^{(1)}[" << i << "]=" << yt1[i] << endl;
17   for (int i=0;i<m;i++)
18     cout << "y^{(2)}[" << i << "]=" << yt2[i] << endl;
19   for (int i=0;i<m;i++)
20     cout << "y^{(1,2)}[" << i << "]=" << yt1t2[i] << endl;
21   return 0;
22 }

4.2.2 結果

以下の結果が生成される:


y[0]=-2.79401891249195
y[1]=-2.79401891249195
y^{(1)}[0]=14.2435494001203
y^{(1)}[1]=11.4495304876283
y^{(2)}[0]=14.2435494001203
y^{(2)}[1]=11.4495304876283
y^{(1,2)}[0]=-149.94964806237
y^{(1,2)}[1]=-124.256568174621

4.3 dco/c++の機能

4.3.1 passive_value

引数のパッシブ値の参照(二次微分型引数の値)を返す。






補足

以降の説明は文献1に従っています。日本語による自動微分の包括的な解説が文献2に在ります。

文献1:Naumann, U. (2012). The Art of Differentiating Computer Programs. An Introduction to Algorithmic Differentiation. Number 24 in Software, Environments, and Tools. SIAM.
文献2:久保田 光一,伊理 正夫(1998),アルゴリズムの自動微分と応用(現代非線形科学シリーズ 3),コロナ社

タンジェントモード(tangent mode)

セクション1.1.1の式が、関数$f$のtangent-linear modelであるタンジェントモード計算の定義です。これは関数$F$の$\mathbf{x}^{(1)}$方向微分を意味します。
これは、ボトムアップ算法、フォワード算法とも呼ばれます[文献2]。

アジョイントモード(adjoint mode)

セクション1.2.1の式が、関数$f$のadjoint modelであるアジョイントモード計算の定義です。
関数解析においては、一般に共役(アジョイント)作用素(adjoint operator)は、ヒルベルト空間上の適切な内積演算により定義されます。ここでは、$\mathbb{R}^n$上の線形作用素であるヤコビアンのアジョイント作用素がこのモードの計算の基礎となります。
$\nabla f:\mathbb{R}^n \rightarrow \mathbb{R}^m$のアジョイント作用素は、以下の式を満たす$(\nabla f)^{*}:\mathbb{R}^m \rightarrow \mathbb{R}^n$として定義されます。

$$ \langle (\nabla f)^{*} \cdot \mathbf{y}_{(1)},\mathbf{x}^{(1)} \rangle_{\mathbb{R}^n} =\langle \mathbf{y}_{(1)},\nabla f \cdot \mathbf{x}^{(1)} \rangle_{\mathbb{R}^m}, $$

ここで、$\langle.,.\rangle_{\mathbb{R}^n}$および$\langle.,.\rangle_{\mathbb{R}^m}$は、それぞれ$\mathbb{R}^n$および$\mathbb{R}^m$における内積を意味します。また、$(\nabla f)^{*}=(\nabla f)^{T}$であることも証明されます[文献1]。
このヤコビアンの転置は写像$\mathbb{R}^m \rightarrow \mathbb{R}^n$です。アジョイントモードでは、全ヤコビアンは$y_{(1)}$を$\mathbb{R}^m$の基底ベクトル全体に渡らせて計算されます。 これは、高速自動微分法、トップダウン算法、リバース算法とも呼ばれます[文献2]。

チェックポイント

アジョイントモードでは、基本関数(element function)の偏導関数値を"逆向き"に掛け合わせていきます。単純なやり方は、フォワード計算(オーバーローディングされたオリジナルコード計算)の計算過程をテープ"tape"と呼ぶメモリー領域に保存し、リバース計算(アジョイント計算)で必要に応じてそれを用いる事です。コードの規模に依りますが、特に大規模システムでは、一度の計算内では計算に使用するマシン上のメモリーに収まらない場合が生じます。こうした状況に対応するための一般的な処方として、指定された途中計算地点のスナップショットを記憶して、リバース計算ではその各地点から再度計算をするようにしています。これをチェックポイントと呼びます。

基本関数(element function)

プログラムとして与えられた関数$f:\mathbb{R}^n \rightarrow \mathbb{R}^m, \mathbf{y}=f(\mathbf{x})$を、以下のように基本関数$\psi$に分解されると仮定します。 $$ v_{j}=\psi(v_{i}),j=n,...n+p+m-1 $$ ここで、n個の独立な入力変数$x_{i}=v_{i},i=0,...,n-1$は、m個の従属変数$y_{j}=v_{n+p+j},j=0,...,m-1$へ写像差されます。pは、中間変数$v_{k},k=n,...,n+p-1$です。通常は中間変数を四則演算や初頭関数などが想定され、連続で偏微分可能であることが仮定されます。
文献2では「基本演算」と参照しています。

アクティブ(active)とパッシブ(passive)

ヤコビアン行列の全てを計算対象とするのは効率的でありません。一般に、予めゼロであることが明らかな要素や出力が不要な要素が有ります。また偏微分計算に無関係な変数も存在します。アクティブ値とは、少なくとも1つの独立変数の値に依存し、また少なくとも一つの従属変数の値に影響を与える事を言います。アクティブ値を持たないことが明らかでないプログラム変数は全てアクティブです。それ以外はパッシブと呼びます。

テープ(tape)

テープは、オーバーローディングにより静的に割当てられる配列です。アジョイントモードにおけるフォワード計算(オーバーローディングされたオリジナルコード計算)の過程が保存されます。そのエントリーは、演算コード(基本関数マクロ)と、その演算の引数(1か2変数)のアドレス、現在の値とアジョイント値から成ります。後続のリバース計算(アジョイント計算)で逆向きに解釈されます。

$\mathbf{x}^{(1)}$

タンジェントモードの方向微分ベクトル入出力変数については、その識別に上付き括弧を用いて微分次数を添えます。

$\mathbf{y}_{(1)}$

アジョイントモードの方向微分ベクトル入出力変数については、その識別に下付き括弧を用いて微分次数を添えます。

Results matter. Trust NAG.

Privacy Policy | Trademarks