トップ サイトマップ
関連情報

自動微分の事例

数理ファイナンスにおける典型的な数値計算パターンへのアジョイント自動微分ツールの適用(前半部)

Uwe Naumann, Jacques du Toit(*)
RWTH Aachen University, Germany,
(*)Numerical Algorithms Group(NAG) Ltd., UK


[2016年1月掲載]

概要

本レポートでは、数理ファイナンスに現れる数値計算パターン(カーネル)に対する、オーバーローディング技術を基礎にした柔軟で使いやすいC++自動微分(AD)ツールの適用事例を紹介します。金融分野においてはアジョイントあるいは通常のAD手法が知られていますが、製品としてC++コードを統合するツールは僅かしかありません。アジョイントADは極めて強力な手法として知られていますが、これは巨大なメモリーを必要とします。本レポートではこの問題に対処する幾つかの手法を紹介し、金融分野において頻繁に用いられる数値計算カーネルにこれを適用した事例を紹介します。本レポートは、任意のC++コードを高度に柔軟に扱うことが可能な我々が開発したADツール'dco/c++'に焦点を当てています。このツールのコンセプトはまた、内製のツールを含めた他のADソリューションへも転用可能なものです。ここで議論されるADソリューションや数値計算カーネルのソースコード群は、弊社サイトからダウンロード可能です。これには適用事例のコードとdco/c++ツールが含まれており、dco/c++のトライアルライセンスもNAGが提供しています。


1. イントロダクション

ブラック・ショールズおよびマートンの業績以来、条件付請求権のプライシングとリスクマネジメントは、数学的な微分計算と密接に関係してきました。完全市場モデルにおいて、債権をヘッジするにはその価格の微分値を知れば十分であり、結果的に定量ファイナンス実務者の間ではリスクという言葉が微分と同義になりつつあります。

残念ながら多くのファイナンス問題は閉じた形の解になりません。ほとんどの場合、得られるのはコンピュータプログラムによる近似的な解です。数学的な微分を得ることは困難であり、もっとも単純には有限差分といった数学的なテクニックを用いなければなりません。しかし有限差分は (グラジェントの量によっては)コストが掛かり、高次の微分に対しては誤差を生じます。

自動微分(AD) [Griewank and Walter, 2008, Naumann, 2012]は、与えられたコンピュータプログラムの正確な数学的微分を与える数値計算技術です。これは、与えられたコンピュータプログラムを任意の次数かつマシン精度で正確に微分し、いかなる近似も用いません。さらにアジョイントAD(ADD)はグラジェントのサイズに依らず、典型的には元のプログラムの計算コストの定数倍のコストでグラジェントを計算します。このためAD(特にADD)は、計算機科学や計算工学の多くの領域で、例えばADに関する最近の学会[Bischof et al., 2008, Bücker et al., 2005, Forth et al., 2012]においても高く注目されています。ADDにおける主要な問題は、極めて大きなメモリーが必要になる事です。このメモリーを小さくする手法の議論に本レポートの大部分を割きます。

ADはファイナンス領域でしばしば議論されていますが、その適用は比較的にゆっくりとされてきました。[Giles and Glasserman, 2006, Capriotti, 2011]。その理由として以下のことが挙げられます。


@この手法をよく知らない

A大きなプログラムコードに'手作業で'ADを実装することが困難

BアジョイントADがメモリー上限を超えることに対する不安

C現実的なファイナンス問題に如何にADを適用すればよいか分からない


本レポートはこれらの課題に焦点を当てています。最大の難点はAからのものでしょう。実務コードに手動でADを適用することは、単純に現実的でなく、現代的なソフトウェア工学的観点からも適合しません。実務コードに対してはソフトウェアツールを用いるべきです。これらツールは、如何に使い易いか、どのような言語とその特徴をサポートするか、柔軟に対応できるか、そして如何に効率的な微分計算コードを出力するかにより評価されるべきです。我々の提案するツール、dco/c++(オーバーローディングによる微分計算)はNAGとRWTHアーヘン大学の共同作業により開発されました。dco/c++は、数百万行の工業分野におけるC/C++コードへ成功裏に適用されており、使い易く効果的で、かつ高度な柔軟性を持ち合わせています。後述のセクションでdco/c++の完全な記述が示されます。

Cに対するアプローチとして、ファイナンス領域のコードに頻繁に現れる幾つかの共通の'数値計算パターン(カーネル)'を特定しました。それを'エボリューション'(シリアルループ)、'アンサンブル'(並列ループ)、線形ソルバー、求根、および最適化といった基本的な構成ブロックに分けました。本レポートは最初の2つに焦点を当てます。これらのカーネルを利用する小さなファンナンス・アプリケーションを、テスト用ですが実用に値する処方で作成し、これにADを適用しました。そしてメモリー要求および/又は計算コストの削減のために、このカーネルの特性を利用する様々な戦略について紹介します。

後の議論で明確になりますが、ADは本質的にはソースコードの変換技術です。よって全アプリケーションのソースコードに対して適用可能です。しかしこれは、ソースが公開されていないサードパーティのライブラリーを用いる場合には困難となります。後述のセクション6では、これらソースの無いライブラリーのdco/c++による取り扱いを含めた議論から生じる様々な話題について簡単に紹介します。

このテストアプリケーションとADを適用した全てのソースコードは以下のアドレスから入手可能です。

http://www.nag.com/adtoolsforfinance/sourcematerials

このコードにはドキュメント(dco/c++に関する文書)も含まれ、コードの構造は本レポートの数学的表現に準じています。これにより本レポートでのAD適用の詳細と共にコードを勉強でき、さらに自身のコードに適用可能なdco/c++トライアル版も含まれます。dco/c++のトライアルライセンスは上記リンクを通じてNAGから入手できます。


2. ADの要約

ADはセマンティック(意味を変える)変換です。与えられたコードはプライマルコードあるいはプライマル関数と呼ばれます。プライマル関数値の計算に加えて、その変換されたコードはプライマル関数のあるパラメータ群に関する微分値も計算します。$\mathbb{R}^{n+\bar{n}}$から$\mathbb{R}^{m+\bar{m}}$への写像の関数$F$の実装を考えましょう。ここで、$n$個の入力 (アクティブ入力)に関する$F$の$m$個の出力 (アクティブ出力)の微分を知りたいとします。$F$の第二の$\bar{m}$個の出力と$\bar{n}$個の入力は、それぞれパッシブ出力、パッシブ入力と呼ばれます。例えば、アクティブ出力があるオプションのモンテカルロ価格の時は、パッシブ出力は対応する信頼区間となります(セクション3をご覧ください)。アクティブ入力が初期資産価格$S_0$の場合、パッシブ入力はモンテカルロ・シミュレーションで用いる乱数のセットです。

一般性を保ち、記法の単純化のために、議論をスカラーのアクティブ出力、即ち$m=1$とします($m>1$への拡張は容易です)。そこで、以下の多変数関数を考えます。

\begin{equation*} F:\mathbb{R}^{n} \times \mathbb{R}^{\bar{n}} \rightarrow \mathbb{R} \times \mathbb{R}^{\bar{m}}, \quad (y,\mathbf{\tilde{y}})=F(\mathbf{x,\tilde{x}}),  \tag{1} \end{equation*}

ここで、$F$とその実装は必要なオーダーまで微分可能であると仮定します。(補足:微分可能性は(自動)微分に必須の前提条件です。ある点において微分可能でない場合には、ADに局所的な有限差分近似を組み合わせるなどの平滑化手法により扱うことが出来ます。この問題に対する詳細な議論は本レポートの範囲を超えますが、その結果は、一般化された微分情報を扱うことが必要であるにもかかわらず、微分不可能な点の存在下でも有効なものです。最近の、微分不可能な数値計算シミュレーションコードへの拡張に関しては[Griewank, 2013]をご覧ください)

ここで以下の、アクティブ入力$\mathbf{x}$に関するアクティブ出力$y$の全偏微分ベクトルを生成する、(半)自動的な方法を求めることとします。

\begin{equation*} \nabla F=\nabla F(\mathbf{x},\mathbf{\tilde{x}}) \equiv \left( \frac{\partial y}{\partial x_i} \right)_{i=0,…,n-1} \in \mathbb{R^n}, \tag{2} \end{equation*}

ここで、このグラジェントはアクティブ入力とパッシブ入力の関数としての全てのアクティブ出力とパッシブ出力に呼応するものです。同様に、$\mathbf{x}$に関する$y$の全ての二次微分であるヘッシアンは以下のようになります。

\begin{equation*} \nabla^2 F=\nabla^2 F(\mathbf{x},\mathbf{\tilde{x}}) \equiv \left( \frac{\partial^2 y}{\partial x_j \partial x_i} \right)_{j,i=0,…,n-1} \in \mathbb{R^{n \times n}}. \tag{3} \end{equation*}

2.1 TangentモードAD

以下では文献[Naumann, 2012]の記法を用います。'フォワード(tangent)モードAD'はプライマル関数$F$のtangent版$F^{(1)}:\mathbb{R^n} \times \mathbb{R^n} \times \mathbb{R}^{\bar{n}} \rightarrow \mathbb{R} \times \mathbb{R} \times \mathbb{R}^{\bar{m}}$を生成します。Tangentコードは関数$(y,y^{(1)},\mathbf{\tilde{y}}):=F^{(1)}(\mathbf{x},\mathbf{x}^{(1)},\mathbf{\tilde x})$であり、ここで、

\begin{equation*} y^{(1)}:=\nabla F(\mathbf{x},\mathbf{\tilde{x}})^T \cdot \mathbf{x}^{(1)} \tag{4} \end{equation*}

は、$\mathbf{x}^{(1)}$方向への$F$の方向微分です。ここで、上肩記号(1)は一次のtangentを示すこととします。演算子:=は命令型プログラミング言語のアサインを示し、数学的な等号と混同しないよう用いる事としました。全グラジェントは、'シード'と'ハーベスト'として知られる処方でn回のtangentコード実行で一つづつ計算することが出来ます:(4)式のベクトル$\mathbf{x^{(1)}}$は$\mathbb{R^n}$のカルテシアン基底ベクトルとして次々に種まき(シード)されて、tangentコードが実行され、対応するグラジェントが$y^{(1)}$から収穫(ハーベスト)されます。このグラジェントはマシン精度で、有限差分近似での場合と同様な計算コストO(n)・Cost(F)で計算されます。

2.2 アジョイントモードAD

リバース(アジョイント)モードADはプライマル関数$F$のアジョイント版を生成します。

\begin{equation*} F_{(1)}:\mathbb{R^n} \times \mathbb{R^n} \times \mathbb{R}^{\bar{n}} \times \mathbb{R} \rightarrow \mathbb{R} \times \mathbb{R}^{\bar{m}} \times \mathbb{R^n} \times \mathbb{R} \end{equation*}

アジョイントコードは関数$(y,\mathbf{\tilde{y}},\mathbf{x_{(1)}},y_{(1)}):=F_{(1)}(\mathbf{x},\mathbf{x_{(1)}},\mathbf{\tilde x},y_{(1)})$であり、各$\mathbf{x}_{(1)},y_{(1)}$は次のようになります。

\begin{eqnarray} x_{(1)} &:=& x_{(1)}+\nabla F(\mathbf{x},\mathbf{\tilde{x}})^T \cdot \mathbf{y}_{(1)} \\ y_{(1)} &:=& 0 \tag{5} \end{eqnarray}

アジョイントコードは、グラジェントと与えられたアクティブ出力のアジョイント$y_{(1)}$との内積により、与えられたアクティブ入力のアジョイント$\mathbf{x}_{(1)}$をインクリメントします。アジョイントのアウトプットは0にリセットされます。これはプログラム変数$y$で表される以前にプライマルコードで用いた値のアジョイントを正しくインクリメントするためです。詳しくは文献[Naumann, 2012]をご覧ください。ここで、下付き添え字の(1)は、一次のアジョイントを表します。初期化(シード)$\mathbf{x}_{(1)}=0$および$y_{(1)}=1$により、1回のアジョイントコード実行から$\mathbf{x}_{(1)}$のグラジェントが生成されます。またこのグラジェントはマシン精度で計算されます。計算コストはグラジェントのサイズ$n$に依存しません。

2.3 二次微分

二次微分は一次微分の一次微分です。よって二次微分コードはフォワードモードとリバースモードの如何なる組み合わせでも得ることが可能です。フォワード-フォワードモードでは、セクション2.1のtangentコードへフォワードモードADを適用します。このモデルは'シード'される3つのベクトルがあり、全ヘッシアンは二次のtangentコードの$n^2$回の計算から'ハーベスト'されます。計算コストは二次の有限差分近似と同等です。有限差分の場合は通常満足できる精度でないことが多く、単精度の場合は特にそうです。フォワード-フォワードモードADの詳細については文献[Naumann, 2012]をご覧ください。

残りの3つの二次アジョイントモードは数学的に同じものです:詳細は文献[Naumann, 2012]をご覧ください。よってこれらの内一つのみを議論することとします。フォワード-リバースモードADにおいては、セクション2のアジョイントコードへフォワードモードを適用し、プライマルコードの二次のアジョイント版を生成します。このモデルは'シード'される4つ入力を持ち、$n$回の二次のアジョイントコードにより全ヘッシアンを生成します。プライマルコードに対してまずアジョイントモードを適用することにより、計算の複雑さが減少し、それは最終的なヘッシアンにも引き継がれます。ヘッシアンのスパース性がその要因です[Gebremedhin et al., 2005]。最後になりますが、ヘッシアン・ベクトル積は1回の二次アジョイントコードの実行で計算することが出来ます。これはヘッシアンの対角成分がドミナントかつ、対角部のみが要求される場合に役立つでしょう。

2.4 高次微分

高次のADコードは、プライマルコードへtangentあるいはアジョイントモードを適用することにより得られます。詳細については文献[Naumann, 2012]をご覧ください。

2.5 効率とメモリー

ADコードの効率(あるいはコスト)は通常以下の比で評価されます。

\begin{equation*} \mathcal{R}=\frac{アジョイントコード実行時間}{プライマルコードの実行時間} \tag{6} \end{equation*}

一次のアジョイントコードでは、比$\mathcal{R}$は1から100の間の範囲になりますが、これはプライマルコードの数学的かつ構造的な性質と、それらがADの中で如何にうまく利用されたかに依存します。フォワードモードADあるいは有限差分との比較では、比$\mathcal{R}$は、アジョイントを用いたほうが好ましくなるには、どのくらい大きなグラジェントが必要かという指標になります。ツールdco/c++では、この比は、小さなグラジェントの場合でも、通常15以下で10以下になるケースも多くあります。

アジョイントAD手法の複雑さの特徴は、プライマルコードを実質的に逆向き実行しなくてはならないことです。この詳細の議論は、本レポートのスコープから外れるため、文献[Naumann, 2012]を参照してください。プライマルコードの逆向き実行は、プライマルコードで生成されるある種の中間的な値を、逆向き(逆向きのデータフロー)に保存しなければならないことを意味します。これは極めて単純なシミュレーションにおいてさえ、巨大なメモリー要求を生じます。

さらに、演算子オーバーロードに基づく(dco/c++のような)ADの実装は、計算の正確性を期するために実行時に全計算のイメージ(通常'テープ'と呼ばれます)を保存かつ再生することが必要となるため、この傾向をさらに増大させます。このテープは、アジョイントコードがフォワード実行されたときに生成・増加し、データフローを逆向きにするために解釈(プレイバック)され、アジョイントを実行します。結果的に、アジョイントコード(特にテープを用いる場合)は、実際にADを適用する際に非現実的なメモリー要求をもたらすことがあります。セクション4.1において、dco/c++ツールで如何にこの問題に対処するかを説明します。

3 数値計算パターンとテストアプリケーション

本レポートで扱う数値計算パターンは、モンテカルロ・シミュレーションでの'アンサンブル'(並列ループ)と、放物型PDEソルバーでの'エボリューション'(シリアルループ)です。ツールdco/c++を用いたAD実装における、線形システム、反復法による求根および最適化に対する直接解法のシンボリックな微分については、別のレポートで扱います。

ADの観点からこれらの数値計算パターンへのアプローチを説明するために、簡単なプライシング問題に対する2つの解法を実装しました。これに関連する議論については文献[W. Xu, 2014]を参照してください。本レポートは、すぐに利用可能なソフトウェアツールと関連する実装のドキュメントをサポートしていることを特徴としています。

3.1 モンテカルロ・プライシング

モンテカルロ法は説明を要しないでしょう。ADの観点からは、特にアジョイント計算においてモンテカルロ法は、サンプルパス数に必要メモリーが比例するという問題を持ちます。性能に影響を与えずにメモリー利用を制御することが重要です。

テストアプリケーションとして、ローカルボラティリティ・プロセスでの単純なヨーロピアン・コールオプションを考えましょう。$S=(S_t)_{t\geq 0}$を次のSDEの解とします。

\begin{equation*} dS_t=rS_t dt+\sigma(log(S_t),t)S_t dW_t \tag{7} \end{equation*}

ここで、$W=(W_t)_{t\geq 0}$は標準ブラウン運動、$r\gt0$はリスクフリー金利、$\sigma$はローカルボラティリティ関数です。コールオプションの価格はこのとき、満期$T\gt 0$かつストライク$K\gt 0$に対して、以下で与えられます。

\begin{equation*} V=e^{-rT} \mathbb{E}(S_T-K)^+. \tag{8} \end{equation*}

実際には$\sigma$は、通常は市場観測予測ボラティリティ・サーフェースから計算され、バイキュービック・スプラインか一次元スプライン列として表現されることがあります。そのようなケースを扱ったとすれば、必要以上に複雑なテストアプリケーションになる可能性がありました。その代わりに、コードを単純にするために、$\sigma$を下記のように表現してみます。

\begin{equation*} \sigma(x,t)=g(x) \cdot t \tag{9} \end{equation*}

ここで、$g:\mathbb{R} \mapsto \mathbb{R}_+$は、m次の多項式$p_m$とn次の多項式$q_n$を用いて、以下のパデ近似で与えられます。

\begin{equation*} g(x)=\frac{p_m(x)}{q_n(x)}= \frac{a_0+a_1 x+ \cdots a_m x^m}{b_0+b_1 x+ \cdots b_n x^n } \tag{10} \end{equation*}

この近似を選んだのは、$g$が正値で、変数$x$の全体に渡り良好な近似形を持つためです。なお、$\sigma$をこのように仮定した理由は、コードを簡単にしたいためです:ADの手法は$\sigma$の表現方法に依存しません。

式(8)からコールプライス$V$を計算するために、以下のSDEを満たすlogプロセス$X_t=log(S_t)$へオイラー-丸山スキームを適用します。

\begin{equation*} dX_t=(r-\frac{1}{2}\sigma^2(X_t,t))dt+\sigma(X_t,t)dW_t. \tag{11} \end{equation*}

次に、整数$M$を用いて$\Delta=T/M$として、モンテカルロ・時間ステップ列$t_i=i\Delta,i=1,…,M$を定義して、以下のようにします。

\begin{equation*} X_{t_{i+1}}= X_{t_i} + (r-\frac{1}{2} \sigma^2 (X_{t_i},t_i)) \Delta + \sigma (X_{t_i},t_i)\sqrt{\Delta} Z_i \tag{12} \end{equation*}

ここで、各$Z_i$は正規乱数で、$X_{t_0}=log(S_0)$が初期値です。生成されたlogプロセスの$N$サンプルパスは、価格$V$の推定と信頼区間を得る際の(8)式のモンテカルロ積分で用いられます。ここで、入力パラメータ$K,T,r,S_0,a_0,…,a_m,b_0,…,b_n$に関する$V$の感度を求めるために、ADを利用します。パデ近似の次数を変えれば、このモデルにさらに多くの入力パラメータを簡単に追加することが可能です。

3.2 PDE法

PDE法は扱いやすさに限界があるモンテカルロ法よりも好まれる場合が多いでしょう。放物型PDEソルバーのコアは、時間の一点における解の値を示す、状態の時間発展法です。ある時間から次のステップへの移動には、通常例えば陰解法や混合法(クランク・ニコルソン等)等の線形システム解法が含まれます。

アジョイントADの観点からは、この状態の時間発展形式は、時間ステップに比例したメモリーが必要になるという問題が生じます。よってこのメモリー利用について制御が必要です。さらに、普通に整形ソルバーを用いる事(詳しくは後程明らかにします)は、O($n^3$)オーダーのメモリーと計算コスト要求を生じます。ソルバーの構造をうまく考えれば、これらのコストはO($n^2$)まで減らせます。大きく稠密な系では大幅な改善が得られます。我々のスティフネス行列は三重対角に纏められますが、離散過程での偏微分積分方程式(PIDE)では大規模密行列が生じます(例えば、[Cont and Tankov, 2004, Matache et al., 2004]を参照してください)。密な線形システムは、最近傍相関行列問題(例えば、[Borsdorf and Higham, 2010, Qi and Sun, 2006]を参照してください)あるいは非線形最適化のような他の領域でも見られます。線形ソルバーを如何に正確に扱うかは各方面で役に立ちます。

テストアプリケーションは、セクション3.1で述べた同じ問題をクランク・ニコルソン・スキームで解くものです。特に、式(9)と(10)で与えられる$\sigma$を用いたSDE(12)を採用し、式(8)のプライシング問題を考察します。プライシング問題をPDEで考察するために、値$V$を以下の関数$V(x,t)$へ拡張します。

\begin{equation*} V(x,t)=e^{-r(T-t)} \mathbb{E}_{x,t} (e^{X_{T}}-K)^+ \tag{13} \end{equation*}

ここで$\mathbb{P}_{x,t}$を、値$x\in \mathbb{R}$の時間$t\in [0,T]$で開始するマルコフ過程$X$の測度とすれば、$\mathbb{E}_{x,t}$は$\mathbb{P}_{x,t}$に関する期待値を示します。するとマルコフ過程の基本的な考え方により、$V$は次の放物型PDEを満たすことが示されます。

\begin{eqnarray} & & 0=\frac{\partial}{\partial t}V(x,t) +(r-\frac{1}{2}\sigma^2(x,t))\frac{\partial}{\partial x}V(x,t) \\ & & \quad +\frac{1}{2}\sigma^2(x,t)\frac{\partial^2}{\partial x^2}V(x,t)-rV(x,t) \quad for \ (x,t) \in \mathbb{R} \times [0,T), \tag{14} \\ & & (e^x-K)^+=V(x,T) \quad for \ all \ x \in \mathbb{R}. \tag{15} \end{eqnarray}

この式に対して、以下の漸近的な境界条件を課します。

\begin{eqnarray} & & \displaystyle \lim_{ n \to -\infty } V(x,t)=0 \quad for \ all \ t \in [0,T], \tag{16} \\ & & \displaystyle \lim_{ n \to \infty } V(x,t)=e^{-r(T-t)}(e^x-K) \quad for \ all \ t \in [0,T]. \tag{17} \end{eqnarray}

このシステムは文献[Andersen and Brotherton-Ratcliffe, 2000]で示されているようなクランク-ニコルソン・スキームにより解くことが出来ます。ここで以前と同様に、感度計算にADを用います。

4 C++におけるオーバーローディングを用いたAAD

理想的にはADはソース変換により実装され、静的プログラム構造解析とコンパイラの最適化技術を全面的に利用します。手動によるソース変換は、単純なコード以外では面倒でエラーを生じやすい作業である事は明らかです。残念ながら、C++に対するソース変換ツールの現状はまだ初等的な段階です。現存するツール[Hascoët and Pascual, 2013, Schmitz et al., 2011, Voßbeck et al., 2008]は、C言語で書かれたコードに対しては満足できるレベルにありますが、C++に対してはそうではありません。メタプログラミング技法と合わせた演算子オーバーローディングは、いかなるC++コードへもADを適用できる有用な別の方法です。C++に対するオーバーローディングADツールは最先端の技術である以下の特徴を持ちます:

  1. ・一般的なプライマルコード(テンプレート)の微分タイプのインスタンス化により、一次および高次のtangentとアジョイントモデルを提供します。プライマルコードとその任意次数の微分コードの最適化は、一般的により高次の微分コードでも再利用可能です。例えば、(最適化されたチェックポイントおよび/又はシンボリックな微分カーネルを伴う場合でもおそらく)最適化された一次微分アジョイントコードの'基本型'(セクション4.2をご覧ください)を、浮動小数点型から一次微分タイプへ変更することにより、効率的な二次微分アジョイントコードを得ることが出来ます。

  2. ・一次のアジョイントコードの式(6)のコスト比$\mathcal{R}$は最小になります。テープがすべてメモリーに入りきれば、理想的には10以下です。プライマルコードに対し極限まで最適化を進めればそこに到達できます。このことをアジョイントコードの実行時間を確認するときはいつも頭に入れておくべきです。如何なる場合でも、$\mathcal{R}$がないレポートは価値が半減します。

  3. ・一次のアジョイントコードで生成されるテープのメモリーは最小化されます。低コストの圧縮技術とキャッシュ最適化されたデータ構造を用いて、理想的には最小限の項目のみ記録されます。

これらの課題は常に、様々な学術的なADツール開発の事柄として取り扱われてきたものです。ですがオーバーローディング・ツールの性能(アジョイントモードでのメモリーや実行時間)は様々です。他の研究者による多くの比較や実験により、dco/c++は同等というよりもむしろ最も実行時間が少なく、メモリーが最も小さいことが示されてします。しかしながらここで、ファイナンス実務者にとってより重大と考えられる2つの追加の課題を見つけることが出来ました:

  1. ・ 実際のシミュレーションコードでは、完全に自動でアジョイントモードを得ることは稀です。チェックポイント、ユーザー定義アジョイントコード、および他の進んだ技術(続く本節に記しています)は、ADツールにより生成される内部データに対する、柔軟かつ直観的な操作が必要となります。このような操作が不要なのはプログラムが小さい場合のみです。定量的ファイナンス・ライブラリーは多分、その実際の対象にはならないでしょう。

  2. ・ その微分コードの保守の簡易さも重要です。このため、既存の使用中の(複雑な)ソフトウェア・プラットフォームへADツールを完全に組み込むことが強く望まれます。そうすれば開発や、保守およびそれらの配置に関して多大な簡略化がもたらされます。この際には、テストや、システマチックなバグ追跡と修正や、バージョン管理および専門的なサポートを通して、ADツールが最新のソフトウェア工学的手法を保有することが必要となります。ソフトウェア開発工程内にアジョイントADを中心的なコンポーネントとすることは、コードの複雑性を増やすことになりますが、出来る限りこれを小さくすべきです。

ツールdco/c++はこれらの課題に対して開発されました。分散あるいはシェアードメモリーの並列計算において効果的に用いられ[Schanen et al., 2012]、アクセラレータ(GPU)と併用され[Du Toit et al., 2014]、多くのTier1投資銀行や他の工業界(自動車、航空宇宙産業など)において現在利用されています。

4.1 dco/c++とは

dco/c++は演算子オーバーローディング・ADツールです。セクション2.5では、このようなツールが直面する(主にメモリー要求)課題を議論しましたが、これからdco/c++がそれらを如何に扱うかを見ていきましょう。

アジョイントモードで用いる場合は、dco/c++は実行時に主要な計算結果をテープに保存しなくてはなりません。一次のアジョイントdco/c++データ型を持つジェネリック型(テンプレート) プライマルコードの実行時間は、テープへ記録しない場合でさえ2倍まで増加します。主メモリーへテープを保存する場合は、プライマルコードの構造に依存して、その他に2倍から6倍の実行時間が追加されます(dco/c++の代入レベルの圧縮アルゴリズムでは、記録コストはその代入の右辺の長さに依存します)。デフォルトでは、テープは動的に必要に応じてメモリーを増大させますが、静的なメモリーサイズ指定オプションもあります(テープがどのくらい大きくなるかを知る必要があります)。

テープの解釈計算には、キャッシュデータの最適化を行えば、1.5倍以下の追加で通常は十分です。現代的なCPUベースのワークステーション上で、dco/c++が生成した一次のアジョイントコードの実行時間オーバーヘッド比$\mathcal{R}$は、テープがメモリーに収まれば、通常5から15の間になると考えられます。もしテープのサイズが主メモリーサイズを超えてしまう場合、これは実際に最も良く起こるケースですが、さらなる$\mathcal{R}$の増加や減少時においても結果を出すために、さらに追加の手段が必要になります。

可能な手段として、'エボリューション'(シリアルループ)や'アンサンブル'(並列ループ)、またシンボリックな数値微分カーネル(例えば線形システムソルバー)に対するチェックポイント機能があります。dco/c++のテープベースのアジョイント計算にこれらの技術を、如何に組み込むことが出来るかを以下で示します。これら全てを可能にする鍵となるのはユーザー定義アジョイント関数で、これはセクション4.3で議論します。dco/c++が生成するアジョイントコードの効率性とスケーラビリティをさらに改善する方法に、テープへの並列記録や、アクセラレータ(GPU)の利用、および事前評価のテクニックがあります。これらのトピックの幾つかを、本レポートのスコープを超えない範囲でセクション6で簡単に紹介します。

4.2 dco/c++の利用

dco/c++をセクション3.1で扱ったプライシング問題に当てはめてみましょう。このレポートで用いるプライシング問題を解くコードは、以下に示すコードを扱うこととします。

以下のようなアクティブ入出力のジェネリック(テンプレート)データ型とパッシブ入出力のレギュラー(ノン・テンプレート)データ型を考えましょう。


 1 template<typename ATYPE>
 2 struct ACTIVE_INPUTS {
 3  ATYPE S0 , r ,K,T;
 4  LocalVolSurface<ATYPE> sigmaSq ;
 5 } ;
 6 template<typename ATYPE>
 7 struct ACTIVE_OUTPUTS {
 8  ATYPE V;
 9 } ;
10
11 struct PASSIVE_INPUTS {
12  int N,M;
13  double rngseed [ 6 ] ;
14 } ;
15 struct PASSIVE_OUTPUTS {
16  double ci ;
17 } ;

LocalVolSurfaceは、ローカルボラティリティ・サーフェースのジェネリック型を示し、パデ近似(10)のパラメータを含みます。Vはモンテカルロ・オプション価格、ciは信頼区間の半値です。ここでジェネリック(テンプレート) プライマル関数を考えましょう。


1 template<typename ATYPE>
2 void price (
3  const ACTIVE_INPUTS<ATYPE> &X,
4  const PASSIVE_INPUTS &XP,
5  ACTIVE_OUTPUTS<ATYPE> &Y,
6  PASSIVE_OUTPUTS &YP
7 ) ;

これは入力XPとXを出力YPとYへマッピングします。

4.2.1 一次のtanegentモード

一次のtangent ADモードでは、アクティブ入出力は、dco/c++の一次のtangentデータ型である、dco::gt1s<double>::type over base type doubleでインスタント化されます。そして、以下の$\partial V/\partial S_0$と$\partial V/\partial r$を計算するコードのように、一次のtangentコードをシードおよびハーベストします。


 1 #include "dco.hpp"
 2 typedef dco::gt1s<double> ADMODE;
 3 typedef ADMODE::type AD_TYPE;
 4
 5  ACTIVE_INPUTS<AD_TYPE> X;
 6  PASSIVE_INPUTS XP;
 7  ACTIVE_OUTPUTS<AD_TYPE> Y;
 8  PASSIVE_OUTPUTS YP;
 9
10  dco::derivative(X.S0)=1;
11  price(X,XP,Y,YP) ;
12  cout << "Y=" << dco::value(Y.V) << endl ;
13  cout << "dY/dX. S0=" << dco::derivative(Y.V) << endl ;
14
15  dco::derivative(X.S0)=0;
16  dco::derivative(X.r)=1;
17  price(X,XP,Y,YP) ;
18  cout << "dY/dX.r=" << dco::derivative(Y.V) << endl ;

10行目の呼出しdco::derivativeは、(4)式の$x^{(1)}$を$\partial V/\partial S_0$に相当するカルテシアン基底ベクトルへセットすることと同等のものです。オプション値は12行目で得られ、微分値は13行目で得られます。15-16行目は、$\partial V/\partial r$にカルテシアン基底ベクトル$x^{(1)}$をシードします。微分が18行目でハーベストされるように、その前に再度プライマルコードを実行するようになっています。

4.2.2 一次のアジョイントモード

一次のアジョイントモードは、上述のものと同じジェネリックデータ型とプライマル関数を用います。アクティブデータ型と同様にテープは、dco/c++一次アジョイント型 dco::ga1s::type でインスタント化されます。全てのアクティブ変数はテープへ保存され(13-15行)、プライマルコードはテープを埋めつつ実行します(17行)。


 1 #include "dco.hpp"
 2 typedef dco::ga1s<double> AD_MODE;
 3 typedef ADMODE::type AD_TYPE;
 4 typedef ADMODE::tape_t AD_TAPE_TYPE;
 5 AD_TAPE_TYPE* & AD_TAPE_POINTER=ADMODE::globaltape ;
 6
 7  ACTIVE INPUTS<AD_TYPE> X;
 8  PASSIVE INPUTS XP;
 9  ACTIVE OUTPUTS<AD_TYPE> Y;
10  PASSIVE OUTPUTS YP;
11
12  AD_TAPE_POINTER = AD_TAPE_TYPE::create() ;
13  AD_TAPE_POINTER->registervariable(X.S0) ;
14  AD_TAPE_POINTER->registervariable(X.r) ;
15 . . .
16
17  price(X,XP,Y,YP) ;
18
19  AD_TAPE_POINTER->registeroutputvariable(Y.V) ;
20  dco::derivative(Y.V)=1;
21  AD_TAPE_POINTER->interpretadjoint() ;
22
23  cout << "Y=" << dco::value(Y.V) << endl ;
24  cout << "dY/dX.S0=" << dco::derivative(X.S0) << endl ;
25  cout << "dY/dX.r=" << dco::derivative(X.r) << endl ;
26 . . .
27
28  AD_TAPE_TYPE::remove(AD_TAPE_POINTER) ;

Y.Vは19行目でアクティブ出力としてマークされます。20行目で式(5)の入力アジョイント$y_{(1)}$は、21行目のテープのプレイバック前に1にセットされます。この値は、24-26行で全ての微分が出力される前に、23行目で出力されます。全てのグラジェントは一回のアジョイントコードで得ることが出来ます。テープのメモリーは最終的に28行目に解放されます。

4.3 ユーザー定義アジョイント: テープのギャップに注意

アジョイントモードのデフォルトにおいてdco/c++は、プライマルコードが順方向に実行される間、中間値がテープに記録され、その後テープが解釈される時、実行は逆向きになります。アジョイント計算で利用したいと思う如何なる付加的な情報や構造も、必然的にこのプロセスの妨げになります: dco/c++が通常行うことを止め、その代わりにより洗練された方法でやりたいことを明確にしなくてはなりません。結果としてdco/c++のテープにギャップが生じます。

概念的には次のようなことが発生します。プライマルコード中のセクションgを考えましょう。

\begin{eqnarray} F:(\mathbf{x},\tilde{\mathbf{x}}) \buildrel f_1 \over \longrightarrow (\mathbf{u},\tilde{\mathbf{u}}) \buildrel g \over \longrightarrow (\mathbf{v},\tilde{\mathbf{v}}) \buildrel f_2 \over \longrightarrow (y,\tilde{\mathbf{y}}). \tag{18} \\ \end{eqnarray}

このgは、ユーザーが特別に扱いたいと思うセクション(例えば、モンテカルロ計算ループ)を示しています。順方向の計算は$f_1$を計算し、テープに記録します。実行がgに達したときに、dco/c++をいったん止めて、テープに保存せずにgを計算する、即ち'ギャップ'を作らねばなりません。このセクションの終わりでは、テープを展開して通常のdco/c++の$f_2$計算へ制御が戻されます。このプロセスを図1に示されています。gを計算してテープへの記録を止める時にdco/c++に割込みをかけるのは、dco/c++データ型を用い無いようにするためです。図1では$\mathbf{u}'$と$\mathbf{v}'$で示していますが、gは標準的な浮動小数点データ型で計算する必要があるからです。

テープが解釈される時には、dco/c++は$f_2$のアジョイント計算から始まります。そしてテープのギャップに達したら、微分計算におけるチェーンルールに従ってgのアジョイントを計算する「ユーザー定義アジョイント関数」を呼び出します。こうすれば、次の$f_1$のアジョイントを計算するdco/c++へ制御を渡すことが出来て、テープの解釈はうまく完了します。

ユーザー定義アジョイント関数は、dco/c++の「外部関数インターフェイス」により登録できます。dco/c++はテープに関連付けられる外部関数オブジェクト(EFO)を扱う「ファクトリー」機能([Gamma et al., 1995]を参照してください)を提供し、これはC++の例外回避とメモリーリーク回避を保証します。上述の関数$(\mathbf{v},\tilde{\mathbf{v}})=g(\mathbf{u},\tilde{\mathbf{u}})$のアクティブ入力$\mathbf{u}$とアクティブ出力$\mathbf{v}$は、EFOに登録されます。さらにギャップのヤコビアン$\nabla g \equiv \partial \mathbf{v}/\partial \mathbf{u}$の計算に必要とされる他のデータと共に、パッシブ入力$\tilde{\mathbf{u}}$とパッシブ出力$\tilde{\mathbf{v}}$もEFOに登録されます。

関数gのユーザーによる実装(例えば、g_make_gapと呼びましょう)により順方向で実行される時に、EFOは生成され配置されます。ギャップを埋めるユーザー定義アジョイント関数(g_fill_gapと呼びましょう)はEFOに登録される必要があります。テープが解釈される時、dco/c++は$f_2$のアジョイント$\mathbf{v}_{(1)}$の計算から開始します。この解釈がギャップに達したとき、gのアクティブ入力のアジョイント$\mathbf{u}_{(1)}:= \mathbf{u}_{(1)}+(\nabla \mathbf{v})^T \cdot \mathbf{v}_{(1)}$が計算されます。そしてdco/c++に制御が戻され、$\mathbf{x}_{(1)}:= \mathbf{x}_{(1)}+(\nabla F)^T \cdot y_{(1)}$を得るためのテープの解釈を完了します。

外部関数インターフェイスは、ユーザーの情報をdco/c++アジョイントADソリューションへ組込むためのオプションです。セクション5において我々のテストアプリケーションにおける技法として、モンテカルロ・シミュレーションでの独立性の利用やチェックポイントを用いた浮動小数点演算のためのメモリー流用について議論します。外部関数インターフェイスはdco/c++を用いたAADのアプローチの様々な他の実装を簡便にします。その内の幾つかに関してセクション6で簡単に議論します。

ここではまだ、ユーザーが如何にしてギャップを埋めて、g_fill_gapにおいてgのアジョイントを計算するかを説明していません。マニュアル作成アジョイントやソース変換技術等の解析的な結果が選べ(セクション6をご覧ください)、最終的にdco/c++それ自身を用いてギャップを埋めることが出来ることが明確になります。これは以下のセクション5で議論します。


図1: dco/c++の外部関数インターフェイス: 図(a)はセクション4.3でのテープ生成を表しています。その計算のアクティブ部分(実線の円)がテープです。破線は(ゼロでない)入力に関する出力の偏微分を表します。点線はテープから外部関数、あるいはその逆へのアサインで、偏微分の単位を表します。図(b)はセクション4.3で議論したテープの解釈を表します。アクティブな計算は、g_fill_gapがテープのギャップを埋める間にdco/c++により行われます。


セクション5以降の後半の議論はリンク:「数理ファイナンスにおける典型的な数値計算パターンへのアジョイント自動微分ツールの適用(後半部)」をご覧ください。

Results matter. Trust NAG.

Privacy Policy | Trademarks