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]。その理由として以下のことが挙げられます。
① この手法をよく知らない
② 大きなプログラムコードに '手作業で' AD を実装することが困難
③ アジョイント AD がメモリー上限を超えることに対する不安
④ 現実的なファイナンス問題に如何に AD を適用すればよいか分からない
本レポートはこれらの課題に焦点を当てています。最大の難点は ② からのものでしょう。実務コードに手動で AD を適用することは、単純に現実的でなく、現代的なソフトウェア工学的観点からも適合しません。実務コードに対してはソフトウェアツールを用いるべきです。これらツールは、如何に使い易いか、どのような言語とその特徴をサポートするか、柔軟に対応できるか、そして如何に効率的な微分計算コードを出力するかにより評価されるべきです。我々の提案するツール、dco/c++(オーバーローディングによる微分計算)は nAG と RWTH アーヘン大学の共同作業により開発されました。dco/c++ は、数百万行の工業分野における C/C++ コードへ成功裏に適用されており、使い易く効果的で、かつ高度な柔軟性を持ち合わせています。後述のセクションで dco/c++ の完全な記述が示されます。
④ に対するアプローチとして、ファイナンス領域のコードに頻繁に現れる幾つかの共通の '数値計算パターン'(カーネル)を特定しました。それを 'エボリューション'(シリアルループ)、'アンサンブル'(並列ループ)、線形ソルバー、求根、および最適化といった基本的な構成ブロックに分けました。本レポートは最初の2つに焦点を当てます。これらのカーネルを利用する小さなファンナンス・アプリケーションを、テスト用ですが実用に値する処方で作成し、これに AD を適用しました。そしてメモリー要求および/又は計算コストの削減のために、このカーネルの特性を利用する様々な戦略について紹介します。
後の議論で明確になりますが、AD は本質的にはソースコードの変換技術です。よって全アプリケーションのソースコードに対して適用可能です。しかしこれは、ソースが公開されていないサードパーティのライブラリーを用いる場合には困難となります。後述のセクション 6 では、これらソースの無いライブラリーの dco/c++ による取り扱いを含めた議論から生じる様々な話題について簡単に紹介します。
このテストアプリケーションと AD を適用した全てのソースコードは以下のアドレスから入手可能です。
このコードにはドキュメント(dco/c++ に関する文書)も含まれ、コードの構造は本レポートの数学的表現に準じています。これにより本レポートでの AD 適用の詳細と共にコードを勉強でき、さらに自身のコードに適用可能な dco/c++ トライアル版も含まれます。dco/c++ のトライアルライセンスは上記リンクを通じて nAG から入手できます。
2 AD の要約
AD はセマンティック(意味を変える)変換です。与えられたコードはプライマルコードあるいはプライマル関数と呼ばれます。プライマル関数値の計算に加えて、その変換されたコードはプライマル関数のあるパラメータ群に関する微分値も計算します。Rn+ˉn から Rm+ˉm への写像の関数Fの実装を考えましょう。ここで、n 個の入力 (アクティブ入力)に関する F の m 個の出力 (アクティブ出力)の微分を知りたいとします。F の第二の ˉm 個の出力と ˉn 個の入力は、それぞれパッシブ出力、パッシブ入力と呼ばれます。例えば、アクティブ出力があるオプションのモンテカルロ価格の時は、パッシブ出力は対応する信頼区間となります(セクション 3 をご覧ください)。アクティブ入力が初期資産価格 S0 の場合、パッシブ入力はモンテカルロ・シミュレーションで用いる乱数のセットです。
一般性を保ち、記法の単純化のために、議論をスカラーのアクティブ出力、即ち m=1 とします(m>1 への拡張は容易です)。そこで、以下の多変数関数を考えます。
F:Rn×Rˉn→R×Rˉm,(y,˜y)=F(x,˜x),
ここで、F とその実装は必要なオーダーまで微分可能であると仮定します。(補足:微分可能性は(自動)微分に必須の前提条件です。ある点において微分可能でない場合には、AD に局所的な有限差分近似を組み合わせるなどの平滑化手法により扱うことが出来ます。この問題に対する詳細な議論は本レポートの範囲を超えますが、その結果は、一般化された微分情報を扱うことが必要であるにもかかわらず、微分不可能な点の存在下でも有効なものです。最近の、微分不可能な数値計算シミュレーションコードへの拡張に関しては [Griewank, 2013] をご覧ください。)
ここで以下の、アクティブ入力 x に関するアクティブ出力 y の全偏微分ベクトルを生成する、(半)自動的な方法を求めることとします。
∇F=∇F(x,˜x)≡(∂y∂xi)i=0,…,n−1∈Rn,
ここで、このグラジェントはアクティブ入力とパッシブ入力の関数としての全てのアクティブ出力とパッシブ出力に呼応するものです。同様に、x に関する y の全ての二次微分であるヘッシアンは以下のようになります。
∇2F=∇2F(x,˜x)≡(∂2y∂xj∂xi)j,i=0,…,n−1∈Rn×n.
2.1 Tangent モード AD
以下では文献 [Naumann, 2012] の記法を用います。'フォワード(tangent)モード AD' はプライマル関数 F の tangent 版 F(1):Rn×Rn×Rˉn→R×R×Rˉm を生成します。Tangent コードは関数 (y,y(1),˜y):=F(1)(x,x(1),˜x) であり、ここで、
y(1):=∇F(x,˜x)T⋅x(1)
は、x(1) 方向への F の方向微分です。ここで、上肩記号 (1) は一次の tangent を示すこととします。演算子 := は命令型プログラミング言語のアサインを示し、数学的な等号と混同しないよう用いる事としました。全グラジェントは、'シード' と 'ハーベスト' として知られる処方で n 回の tangent コード実行で一つづつ計算することが出来ます:(4) 式のベクトル x(1) は Rn のカルテシアン基底ベクトルとして次々に種まき(シード)されて、tangent コードが実行され、対応するグラジェントが y(1) から収穫(ハーベスト)されます。このグラジェントはマシン精度で、有限差分近似での場合と同様な計算コスト O(n)⋅Cost(F) で計算されます。
2.2 アジョイントモード AD
リバース(アジョイント)モード AD はプライマル関数 F のアジョイント版を生成します。
F(1):Rn×Rn×Rˉn×R →R×Rˉm×Rn×R
アジョイントコードは関数 (y,˜y,x(1),y(1)):=F(1)(x,x(1),˜x,y(1)) であり、各 x(1),y(1) は次のようになります。
x(1):=x(1)+∇F(x,˜x)T⋅y(1)y(1):=0
アジョイントコードは、グラジェントと与えられたアクティブ出力のアジョイント y(1) との内積により、与えられたアクティブ入力のアジョイント x(1) をインクリメントします。アジョイントのアウトプットは 0 にリセットされます。これはプログラム変数 y で表される以前にプライマルコードで用いた値のアジョイントを正しくインクリメントするためです。詳しくは文献 [Naumann, 2012] をご覧ください。ここで、下付き添え字の (1) は、一次のアジョイントを表します。初期化(シード)x(1)=0 および y(1)=1 により、1 回のアジョイントコード実行から x(1) のグラジェントが生成されます。またこのグラジェントはマシン精度で計算されます。計算コストはグラジェントのサイズ n に依存しません。
2.3 二次微分
二次微分は一次微分の一次微分です。よって二次微分コードはフォワードモードとリバースモードの如何なる組み合わせでも得ることが可能です。フォワード-フォワードモードでは、セクション 2.1 の tangent コードへフォワードモード AD を適用します。このモデルは 'シード' される3つのベクトルがあり、全ヘッシアンは二次の tangent コードの n2 回の計算から 'ハーベスト' されます。計算コストは二次の有限差分近似と同等です。有限差分の場合は通常満足できる精度でないことが多く、単精度の場合は特にそうです。フォワード-フォワードモード AD の詳細については文献 [Naumann, 2012] をご覧ください。
残りの3つの二次アジョイントモードは数学的に同じものです:詳細は文献 [Naumann, 2012] をご覧ください。よってこれらの内一つのみを議論することとします。フォワード-リバースモード AD においては、セクション 2 のアジョイントコードへフォワードモードを適用し、プライマルコードの二次のアジョイント版を生成します。このモデルは 'シード' される4つ入力を持ち、n 回の二次のアジョイントコードにより全ヘッシアンを生成します。プライマルコードに対してまずアジョイントモードを適用することにより、計算の複雑さが減少し、それは最終的なヘッシアンにも引き継がれます。ヘッシアンのスパース性がその要因です [Gebremedhin et al., 2005] 。最後になりますが、ヘッシアン・ベクトル積は 1 回の二次アジョイントコードの実行で計算することが出来ます。これはヘッシアンの対角成分がドミナントかつ、対角部のみが要求される場合に役立つでしょう。
2.4 高次微分
高次の AD コードは、プライマルコードへ tangent あるいはアジョイントモードを適用することにより得られます。詳細については文献 [Naumann, 2012] をご覧ください。
2.5 効率とメモリー
AD コードの効率(あるいはコスト)は通常以下の比で評価されます。
R=アジョイントコード実行時間プライマルコードの実行時間
一次のアジョイントコードでは、比 R は 1 から 100 の間の範囲になりますが、これはプライマルコードの数学的かつ構造的な性質と、それらが AD の中で如何にうまく利用されたかに依存します。フォワードモード AD あるいは有限差分との比較では、比 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=(St)t≥0 を次の SDE の解とします。
dSt=rStdt+σ(log(St),t)StdWt
ここで、W=(Wt)t≥0 は標準ブラウン運動、r>0 はリスクフリー金利、σ はローカルボラティリティ関数です。コールオプションの価格はこのとき、満期 T>0 かつストライク K>0 に対して、以下で与えられます。
V=e−rTE(ST−K)+.
実際には σ は、通常は市場観測予測ボラティリティ・サーフェースから計算され、バイキュービック・スプラインか一次元スプライン列として表現されることがあります。そのようなケースを扱ったとすれば、必要以上に複雑なテストアプリケーションになる可能性がありました。その代わりに、コードを単純にするために、σ を下記のように表現してみます。
σ(x,t)=g(x)⋅t
ここで、g:R↦R+ は、m次の多項式 pm とn次の多項式 qn を用いて、以下のパデ近似で与えられます。
g(x)=pm(x)qn(x)=a0+a1x+⋯amxmb0+b1x+⋯bnxn
この近似を選んだのは、g が正値で、変数 x の全体に渡り良好な近似形を持つためです。なお、σ をこのように仮定した理由は、コードを簡単にしたいためです:AD の手法は σ の表現方法に依存しません。
式 (8) からコールプライス V を計算するために、以下の SDE を満たす log プロセス Xt=log(St) へオイラー-丸山スキームを適用します。
dXt=(r−12σ2(Xt,t))dt+σ(Xt,t)dWt.
次に、整数Mを用いて Δ=T/M として、モンテカルロ・時間ステップ列 ti=iΔ,i=1,…,M を定義して、以下のようにします。
Xti+1=Xti+(r−12σ2(Xti,ti))Δ+σ(Xti,ti)√ΔZi
ここで、各 Zi は正規乱数で、Xt0=log(S0) が初期値です。生成された log プロセスの N サンプルパスは、価格 V の推定と信頼区間を得る際の (8) 式のモンテカルロ積分で用いられます。ここで、入力パラメータ K,T,r,S0,a0,…,am,b0,…,bn に関する V の感度を求めるために、AD を利用します。パデ近似の次数を変えれば、このモデルにさらに多くの入力パラメータを簡単に追加することが可能です。
3.2 PDE 法
PDE 法は扱いやすさに限界があるモンテカルロ法よりも好まれる場合が多いでしょう。放物型 PDE ソルバーのコアは、時間の一点における解の値を示す、状態の時間発展法です。ある時間から次のステップへの移動には、通常例えば陰解法や混合法(クランク・ニコルソン等)等の線形システム解法が含まれます。
アジョイント AD の観点からは、この状態の時間発展形式は、時間ステップに比例したメモリーが必要になるという問題が生じます。よってこのメモリー利用について制御が必要です。さらに、普通に整形ソルバーを用いる事(詳しくは後程明らかにします)は、O(n3) オーダーのメモリーと計算コスト要求を生じます。ソルバーの構造をうまく考えれば、これらのコストは O(n2) まで減らせます。大きく稠密な系では大幅な改善が得られます。我々のスティフネス行列は三重対角に纏められますが、離散過程での偏微分積分方程式(PIDE)では大規模密行列が生じます(例えば、[Cont and Tankov, 2004, Matache et al., 2004] を参照してください)。密な線形システムは、最近傍相関行列問題(例えば、[Borsdorf and Higham, 2010, Qi and Sun, 2006] を参照してください)あるいは非線形最適化のような他の領域でも見られます。線形ソルバーを如何に正確に扱うかは各方面で役に立ちます。
テストアプリケーションは、セクション 3.1 で述べた同じ問題をクランク・ニコルソン・スキームで解くものです。特に、式 (9) と (10) で与えられる σ を用いた SDE (12) を採用し、式 (8) のプライシング問題を考察します。プライシング問題を PDE で考察するために、値 V を以下の関数 V(x,t) へ拡張します。
V(x,t)=e−r(T−t)Ex,t(eXT−K)+
ここで Px,t を、値 x∈R の時間 t∈[0,T] で開始するマルコフ過程 X の測度とすれば、Ex,t は Px,t に関する期待値を示します。するとマルコフ過程の基本的な考え方により、V は次の放物型 PDE を満たすことが示されます。
0=∂∂tV(x,t)+(r−12σ2(x,t))∂∂xV(x,t)+12σ2(x,t)∂2∂x2V(x,t)−rV(x,t)for (x,t)∈R×[0,T),(ex−K)+=V(x,T)for all x∈R.
この式に対して、以下の漸近的な境界条件を課します。
limn→−∞V(x,t)=0for all t∈[0,T],limn→∞V(x,t)=e−r(T−t)(ex−K)for all t∈[0,T].
このシステムは文献 [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 ツールは最先端の技術である以下の特徴を持ちます:
- ・一般的なプライマルコード(テンプレート)の微分タイプのインスタンス化により、一次および高次の tangent とアジョイントモデルを提供します。プライマルコードとその任意次数の微分コードの最適化は、一般的により高次の微分コードでも再利用可能です。例えば、(最適化されたチェックポイントおよび/又はシンボリックな微分カーネルを伴う場合でもおそらく)最適化された一次微分アジョイントコードの '基本型'(セクション 4.2 をご覧ください)を、浮動小数点型から一次微分タイプへ変更することにより、効率的な二次微分アジョイントコードを得ることが出来ます。
- ・一次のアジョイントコードの式 (6) のコスト比 R は最小になります。テープがすべてメモリーに入りきれば、理想的には 10 以下です。プライマルコードに対し極限まで最適化を進めればそこに到達できます。このことをアジョイントコードの実行時間を確認するときはいつも頭に入れておくべきです。如何なる場合でも、 R がないレポートは価値が半減します。
- ・一次のアジョイントコードで生成されるテープのメモリーは最小化されます。低コストの圧縮技術とキャッシュ最適化されたデータ構造を用いて、理想的には最小限の項目のみ記録されます。
これらの課題は常に、様々な学術的な AD ツール開発の事柄として取り扱われてきたものです。ですがオーバーローディング・ツールの性能(アジョイントモードでのメモリーや実行時間)は様々です。他の研究者による多くの比較や実験により、dco/c++ は同等というよりもむしろ最も実行時間が少なく、メモリーが最も小さいことが示されてします。しかしながらここで、ファイナンス実務者にとってより重大と考えられる2つの追加の課題を見つけることが出来ました:
- ・実際のシミュレーションコードでは、完全に自動でアジョイントモードを得ることは稀です。チェックポイント、ユーザー定義アジョイントコード、および他の進んだ技術(続く本節に記しています)は、AD ツールにより生成される内部データに対する、柔軟かつ直観的な操作が必要となります。このような操作が不要なのはプログラムが小さい場合のみです。定量的ファイナンス・ライブラリーは多分、その実際の対象にはならないでしょう。
- ・その微分コードの保守の簡易さも重要です。このため、既存の使用中の(複雑な)ソフトウェア・プラットフォームへ 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++ が生成した一次のアジョイントコードの実行時間オーバーヘッド比 R は、テープがメモリーに収まれば、通常 5 から 15 の間になると考えられます。もしテープのサイズが主メモリーサイズを超えてしまう場合、これは実際に最も良く起こるケースですが、さらなる 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 でインスタント化されます。そして、以下の ∂V/∂S0 と ∂V/∂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::derivativeX.S0=1;
11 priceX,XP,Y,YP ;
12 cout << "Y=" << dco::valueY.V << endl ;
13 cout << "dY/dX. S0=" << dco::derivativeY.V << endl ;
14
15 dco::derivativeX.S0=0;
16 dco::derivativeX.r=1;
17 priceX,XP,Y,YP ;
18 cout << "dY/dX.r=" << dco::derivativeY.V << endl ;
10 行目の呼出し dco::derivative は、(4) 式の x(1) を ∂V/∂S0 に相当するカルテシアン基底ベクトルへセットすることと同等のものです。オプション値は 12 行目で得られ、微分値は 13 行目で得られます。15-16 行目は、∂V/∂r にカルテシアン基底ベクトル x(1) をシードします。微分が 18 行目でハーベストされるように、その前に再度プライマルコードを実行するようになっています。
4.2.2 一次のアジョイントモード
一次のアジョイントモードは、上述のものと同じジェネリックデータ型とプライマル関数を用います。アクティブデータ型と同様にテープは、dco/c++ 一次アジョイント型 dco::ga1s
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->registervariableX.S0 ;
14 AD_TAPE_POINTER->registervariableX.r ;
15 . . .
16
17 priceX,XP,Y,YP ;
18
19 AD_TAPE_POINTER->registeroutputvariableY.V ;
20 dco::derivativeY.V=1;
21 AD_TAPE_POINTER->interpretadjoint ;
22
23 cout << "Y=" << dco::valueY.V << endl ;
24 cout << "dY/dX.S0=" << dco::derivativeX.S0 << endl ;
25 cout << "dY/dX.r=" << dco::derivativeX.r << endl ;
26 . . .
27
28 AD_TAPE_TYPE::removeADTAPEPOINTER ;
Y.V は 19 行目でアクティブ出力としてマークされます。20 行目で式 (5) の入力アジョイント y(1) は、21 行目のテープのプレイバック前に 1 にセットされます。この値は、24-26 行で全ての微分が出力される前に、23 行目で出力されます。全てのグラジェントは一回のアジョイントコードで得ることが出来ます。テープのメモリーは最終的に 28 行目に解放されます。
4.3 ユーザー定義アジョイント:テープのギャップに注意
アジョイントモードのデフォルトにおいて dco/c++ は、プライマルコードが順方向に実行される間、中間値がテープに記録され、その後テープが解釈される時、実行は逆向きになります。アジョイント計算で利用したいと思う如何なる付加的な情報や構造も、必然的にこのプロセスの妨げになります: dco/c++ が通常行うことを止め、その代わりにより洗練された方法でやりたいことを明確にしなくてはなりません。結果として dco/c++ のテープにギャップが生じます。
概念的には次のようなことが発生します。プライマルコード中のセクション g を考えましょう。
F:(x,˜x)f1⟶(u,˜u)g⟶(v,˜v)f2⟶(y,˜y).
この g は、ユーザーが特別に扱いたいと思うセクション(例えば、モンテカルロ計算ループ)を示しています。順方向の計算は f1 を計算し、テープに記録します。実行が g に達したときに、dco/c++ をいったん止めて、テープに保存せずに g を計算する、即ち 'ギャップ' を作らねばなりません。このセクションの終わりでは、テープを展開して通常の dco/c++ の f2 計算へ制御が戻されます。このプロセスを図 1 に示されています。g を計算してテープへの記録を止める時に dco/c++ に割込みをかけるのは、dco/c++ データ型を用い無いようにするためです。図 1 では u′ と v′ で示していますが、g は標準的な浮動小数点データ型で計算する必要があるからです。
テープが解釈される時には、dco/c++ は f2 のアジョイント計算から始まります。そしてテープのギャップに達したら、微分計算におけるチェーンルールに従って g のアジョイントを計算する「ユーザー定義アジョイント関数」を呼び出します。こうすれば、次の f1 のアジョイントを計算する dco/c++ へ制御を渡すことが出来て、テープの解釈はうまく完了します。
ユーザー定義アジョイント関数は、dco/c++ の「外部関数インターフェイス」により登録できます。dco/c++ はテープに関連付けられる外部関数オブジェクト(EFO)を扱う「ファクトリー」機能([Gamma et al., 1995] を参照してください)を提供し、これは C++ の例外回避とメモリーリーク回避を保証します。上述の関数 (v,˜v)=g(u,˜u) のアクティブ入力 u とアクティブ出力 v は、EFO に登録されます。さらにギャップのヤコビアン ∇g≡∂v/∂u の計算に必要とされる他のデータと共に、パッシブ入力 ˜u とパッシブ出力 ˜v も EFO に登録されます。
関数 g のユーザーによる実装(例えば、g_make_gap と呼びましょう)により順方向で実行される時に、EFO は生成され配置されます。ギャップを埋めるユーザー定義アジョイント関数(g_fill_gap と呼びましょう)は EFO に登録される必要があります。テープが解釈される時、dco/c++ は f2 のアジョイント v(1) の計算から開始します。この解釈がギャップに達したとき、g のアクティブ入力のアジョイント u(1):=u(1)+(∇v)T⋅v(1) が計算されます。そして dco/c++ に制御が戻され、x(1):=x(1)+(∇F)T⋅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++ により行われます。
後半の議論は 数理ファイナンスにおける典型的な数値計算パターンへのアジョイント自動微分ツールの適用(後半部) をご覧ください。