iSMART〜The Hyper-Dual Numbers in Fortran
４ヶ月に亘る冬籠りでPython，Juliaの自動微分の実装を終え、Fortran Simulator “iSMART”をJulia言語に変換中にある。スタンフォード大の博士論文を渉猟している折、Hyper-dual numbersなる超実数を用いて厳密に一階偏導関数、加えて二階偏導関数も導出できる論文に出会う。急遽変換を中断し、FortranによるHyper-dual numbersの実装を試みることにする。あの冬籠りは一体何だったんだろうか？After four months of dedication during the winter, I finished the implementation in Python and Julia for the automatic differentiation. I am currently working on converting the in-house Fortran program “iSMART” to Julia. While browsing papers at Stanford University, I came across a Ph.D. thesis on hyper-dual numbers, which is an extension of real numbers that can derive the first-order partial derivatives as well as the second-order partial derivatives in a strict manner. Suddenly I stopped converting and started implementing hyper-dual numbers in Fortran. What were those four months in winter?
Golden Week
ゴールデンウィークの五月晴れの中日、コーディングに疲れて気分転換に隣県へ大好物のおろしそばを食べにツーリングに繰り出す。行き先は天空の城で勝手知った越前大野へ。途中観光名所に立ち寄る。
先ず勝山の恐竜博物館へ向かう。何とリニューアル中で閉鎖中でした。楽しみにしていただけに、情けない。広場では催し物が開催中で家族連れて混雑している。人気のある化石掘りはやってみたいと思う。
次は名水百選に選ばれた御清水へ向かう。
イケる。丸みの有る口当たりでスッキリしている。嫌味は全ない。雑多な山の水よりも格段に上質。殿様に献上されただけは有る。湧いてくる水を杓子ですくいペットボトルに入れて、水割り用に持ち帰る。
昼ご飯の時間になったので久ちゃん食堂へ向かう。
良心的な値段設定もあって混んでいる。蕎麦が茹で上がるのに15分程度待つ。
大盛を注文するが食べ足りない。２杯は注文すべきだ。
１０割蕎麦故か、茹で時間が足らないか故か、蕎麦が固い。
おろしが一般的な味で辛さに欠ける。汁も深み・甘みに欠ける。好みの味ではない。
育った三国の新保屋か石勝食堂の方がよっぽど好みの味。
野宿する予定でテントを持参して来たものの、おろしそばがイマイチだったこともあって野宿する意気込みが失せる。宿泊地予定だった奥越ふれあい公園で二時間ほど昼寝をして気力を回復後、三国湊座、新設の道の駅蓮如に立ち寄ってから帰宅する。
連休最終日、以前から行きたいと思っていた航空プラザへ出向く。
館内見学中、一斉に地震アラートが鳴り出す。十五秒ちょっと横揺れを感じた程度で何事もなくほっとする。暫くして館内放送で能登で震度６強の揺れとのアナウンスを聞き驚く。自衛隊機が被害状況の確認のため離陸したとの放送も受ける。その後も数回揺れを感じる。
館内入口にあるジャイロプレーン。これが一番興味をもった。
ローターに直接動力が繋がっておらず空気の反作用でローターが回転し揚力を確保している。推進力はプロペラで得る。従って、ヘリコプターとは違い離陸するにはかなりの滑走が必要となる。考案した人は偉い。
圧巻は一階の実機展示場。ブルーインパルス(T-2)や第２世代超音速Jet戦闘機F-104戦闘機といったさまざまな機体が展示されている。俺はそれほど実機には興味はないが、好きな人に撮っては堪らないであろう。
二階は飛行機の発展の歴史、模型が展示されている。航空力学を学ぶ実験設備、前政府専用機の貴賓室、シミュレーターを使った飛行体験もある。特に、精巧に作られた航空機の模型は必見。夥しい数の模型が展示されている。
帰宅して甚大な被害が発生した能登地震の実情を知る。
孫が里帰りの折、皆で見学に行こうと考えているが？
Derivation of Partial Derivatives
Forward Difference
PrimitiveなInhouse Simulator “iSMART”ではヤコビアンを構築するために、実装しやすい前進差分法を用いていた。
$$\frac{\partial f(\mathbf{x})}{\partial x_{j}} = \frac{f(\mathbf{x}+h\mathbf{x}{\epsilon_{j}})-f(\mathbf{x})}{h}+O(h)\tag{1}$$
この算法は、演算負荷が高いばかりでなく最適な微小値\(h\)を求める事が難しいデメリットがある。具体的には、SPE#9の課題に於いて真の解からかけ離れた数値解に収束する現象の発生に悩まされた経験が有る。このため、已む無く解析微分に基づいた高精度シミュレータSMARTの開発に至ることとなる。
Automatic Differentiation
自動微分法は単純な原理に基づく。詳細はPart 1をご参照のこと。
\(x\)の偏微分が\(“1″\)であることから、変数値\(x\)と変数の偏微分値のペアを(\(x,1)\)と定義する。定数\(a\)に対しては\((a,0)\)となる。このペアを使って忍耐のいる偏微分導出の計算をコンピュータの演算させるだけだ。
所詮、本法は算法の一種であるがAIの分野でよく使われる。実装方法は大きく分けてForward型とReverse型がある。俺が志向する実装方法はForward型である。
Dual Numbers -> Hyper-dual Numbers
一方、Dual numbers(双対数、乃至は二重数)、其の拡張型のHyper-dual numbers(超双対数)は数学的に秀麗な導出手法に基づき厳密に一階及び二階偏導関数を求めることができる優位性がある。
双対数とは複素数と同じような考えに立ち超実数として定義する。複素数では虚部を表すために\(i\)を用いるが、双対数では\(\epsilon\)なる非実数部を示す表現を用いて、超実数\(x+h\epsilon\)を定義する。双対数が複素数と異なる点は、複素数が\(i^2=−1\)を満たすのに対して\(\epsilon^2=0\)と定めることである。これによって、厳密に一階、二階偏導関数が求めることができる。自乗で零になる定義付とは感動すべき手法である。
\(x,h\in\mathbb{R}\)に対して、双対数\(x+h\epsilon\)は行列を用いて以下のように表せる。
$$x+h\epsilon=\left(\begin{array}{cc}x & h \\ 0 & x \end{array}\right)\tag{2-1}$$
$$\epsilon=\left(\begin{array}{cc}1 & 0 \\ 0 & 1 \end{array}\right)
,~~~~\epsilon^2=0~~(\epsilon\neq 0)\tag{2-2}$$
実関数\(f:\mathbb {R} \to \mathbb{R}\)に対して任意の初期値\(\mathbf{x}_k\in \mathbb{R}\)にて\(f(\mathbf{x})\)をテイラー展開して得られる二次多項近似式はPart 2の式(3-3)で与えられたが、\(x,h\in\mathbb{R}\)に対しては次式が成立する。
$$\begin{align}
f(x+h\epsilon)&=f(x)+f'(x) h\epsilon+\frac{1}{2!}f^{\prime\prime}(x){h}^2\epsilon^2+\cdots \\
&=f(x)+f'(x) h\epsilon \\
∴f'(x)&=\frac{1}{h}f(x+h\epsilon)\end{align}\tag{2-3}$$
即ち、\(\epsilon\)項の係数が元関数の偏微分となる事が分かる。詰り、任意の関数に対して双対数を扱うと自動的にその関数の偏微分値を得ることができることになる。
Hyper-Dual Numbers
これ迄双対数の概念に触れて超実数を定義したが、双対数の特性から\(\epsilon^2=0\)は零になるので、二階以上の高階偏導関数を求めることはできない。そこで、二種類の双対数を扱うことで二階偏微分を求める。
具体的には、
$$x=x_0+x_1\epsilon_1+x_2\epsilon_2+x_3\epsilon_1\epsilon_2 \\
\epsilon_1^2=\epsilon_2^2=0,~~\epsilon_1 \ne \epsilon_2 \ne 0, ~~~~\epsilon_1\epsilon_2=\epsilon_2\epsilon_1\ne 0\tag{3-1}$$
と実数を拡張する。この実数\(x\)のことを超双対数と云う。
超双対数を用いたテイラー展開式は次式で与えられるが、\(\epsilon1\epsilon2=0\)は零になるので、二階迄偏導関数を求めることが可能である。
$$f(x+h_1\epsilon_1+h_2\epsilon_2+0\epsilon_1\epsilon_2) \\
=f(x)+h_1f'(x)\epsilon_1+h_2f'(x)\epsilon_2+h_1h_2f”\epsilon_1\epsilon_2\tag{3-2}$$
次式の関数に於いて、
$$f(x+h_1\epsilon_1e_1+h_2\epsilon_2e_2+0\epsilon_1\epsilon_2)\tag{3-3}$$
偏導関数を評価すると次式らで与えられる。
$$\begin{align}
\frac{\partial f(x)}{\partial x_i} &=\frac{\epsilon_1part[f(x+h_1\epsilon_1e_i+h_2\epsilon_2e_j+0\epsilon_1\epsilon_2)]}{h_1}\\
\frac{\partial f(x)}{\partial x_j} &=\frac{\epsilon_2part[f(x+h_1\epsilon_1e_i+h_2\epsilon_2e_j+0\epsilon_1\epsilon_2)]}{h_2} \\
\frac{\partial^2 f(x)}{\partial x_i \partial x_j} &=\frac{\epsilon_1\epsilon_2part[f(x+h_1\epsilon_1e_i+h_2\epsilon_2e_j+0\epsilon_1\epsilon_2)]}{h_1h_2}
\end{align}\tag{3-4}$$
実際のコーディングに於いて関数値、偏導関数を其々J.Fike氏の博士論文に準じ次のように定める。ここで、\(f0\)は関数値、\(f1\)が一階偏導関数(Gradient)、\(f12\)が二階偏導関数(Hessian)となる。因みに、単変数からなる関数に於いては\(f1\)と\(f2\)は等しくなる。
$$\begin{align}
f0&=f(x+h_1\epsilon_1e_1+h_2\epsilon_2e_2+0\epsilon_1\epsilon_2) \\
f1&= \frac{\partial f(x)}{\partial x_i} \\
f2&=\frac{\partial f(x)}{\partial x_j} \\
f12&=\frac{\partial^2 f(x)}{\partial x_i \partial x_j}
\end{align}\tag{3-5}$$
ここで、次式で与える二種類の超双対数から、
$$\begin{align}
a&=a_0+a_1\epsilon_1+a_2\epsilon_2+a_3\epsilon_1\epsilon_2 \\
b&=b_0+b_1\epsilon_1+b_2\epsilon_2+b_3\epsilon_1\epsilon_2
\end{align}\tag{3-6}$$
加減算、乗算、逆数は次式で与えられる。
$$\begin{align}
a±b&=(a_0±b_0) + (a_1±b_1)\epsilon_1 + (a_2±b_2)\epsilon_2 + (a_3±b_3)\epsilon_1\epsilon_2 \\ \\
ab&=(a_0b_0) + (a_0b_1+a_1b_0)\epsilon_1) + (a_0b_2+a_2b_0)\epsilon_2 \\
&+(a_0b_3+a_1b_2+a_2b_1+a3b_0)\epsilon_1\epsilon_2) \\ \\
\frac{b}{a}&=\frac{b_0}{a_0} ~ – ~b_1 \frac{a_1}{a_0^{2}}\epsilon_1 ~ – ~b_2 \frac{a_2}{a_0^{2}}\epsilon_2 + b_3( \frac{2a_1a_2}{a_0^3}-\frac{a_3}{a_0^2})\epsilon_1\epsilon_2
\end{align}\tag{3-7}$$
ここで、\(\epsilon1\epsilon2\)の係数として二階偏導関数が求まる。特に、除算に於いて\(a\)の場合、除算が実行されないので実数部\(a_0\)の除算は影響されない。
Implementation in Fortran
Fortranへの実装法についてはDual Numbers用がGitHubで公開されている。同FortranコードはYu教授が開発した著名なコードで確かなものであると思慮。これを参考に組込む。
Fortran2003の機能強化を盛り込んだポイントは以下の通り。
- Juliaの構造体”struct”構文に該当する”type”構文にて超双対数のメンバ(成分)を宣言する
- “operator”(四則演算子,+,-,*,/)は他の言語同様型形式に応じて作成してOverloadさせる。型形式(hyper, double precision, integer, 配列等)毎に其々約150種にもなり、コードの総行数は2500行にものぼる。実際全ては必要としない。
- “interface”構文を用いる。interface構文は四則演算子を含む一般的な演算子(+,-,*,/,=等)に加えてユーザー定義した演算子、関数、サブルーチンをOverloadできる機能を持つ。Overloadは適切な引数の関数を自動で見つけて呼び出す総称名関数で定義する。Fortran90で導入された構文との事であるが俺の記憶には残って居らず。
- typeの成分宣言に初期値を設定することができ、又デフォルトの初期化を持つメンバに対しては値を省略できる。
主なコードを下記に紹介する。
局所的にヤコビ行列を構築して演算負荷の低減を図ることを目指したJulia版自動微分モジュールで実装した”LocalAD”と同じ仕組みで非常にシンプルなものである。
module HyperDN_mod
implicit none
type, public:: hyper
! x = f0 + f1 e1 + f2 e2 + f12 e1e2
sequence
double precision:: f0
double precision:: f1 = 1.0d0
double precision:: f2 = 1.0d0
double precision:: f12 = 0.0d0
end type
! Assignment
interface assignment(=)
module procedure assign_hh
module procedure assign_hd
end interface
...
contains
...
end module因みに、type (hyper)に参照(アクセス)するには、
- public指定したメンバを直に参照する
- (built-in)constructorにて直に参照する
の二通りの方法がある。メンバを直にアクセスする方法は、例えば、メンバ変数xhに対してxh%f0, xh%f1, xh%f2, xh%f12等を参照する為見た目は良くない。俺は主にbuilt-in constructorで設定している。
実は当初variableHDNなる関数を組込み、同関数にてメンバの設定を行っていたが、Fortran2003を読み込んでtypeの成分宣言に初期値を設定することができることを知る。
Verification of Fortran code
Python･Julia版同様の課題にてコードの妥当性を検証した。特記事項は以下の通り。
１．以前Part 2で述べたように、ニュートン法は、一階偏導関数\(\nabla f(\mathbf{x})\)勾配ベクトル(Gradient)、二階偏導関数\(\nabla^{2} f(\mathbf{x})\)ヘッセ行列(Hessian)を用いて、
$$\nabla^{2} f(\mathbf{x})= \begin{pmatrix} f_{xx} & f_{xy} \\ f_{yx} & f_{yy} \end{pmatrix}, ~~\nabla f(\mathbf{x})= \begin{pmatrix} f_x \\ f_y \end{pmatrix}\tag{4-1}$$
次式の探索方向にて、
$$\mathbf{x}_{n+1} = \mathbf{x}_{n}-(\nabla^{2} f(\mathbf{x}_{n}) )^{-1} \nabla f(\mathbf{x}_{n})\tag{4-2}$$
逐次的に更新して真の最適解に近づいていくアルゴリズムである。
$$\begin{align}
\begin{pmatrix} x_{n+1} \\ y_{n+1} \end{pmatrix} &=\begin{pmatrix} x_{n} \\ y_{n} \end{pmatrix}
-{\begin{pmatrix} f_{xx} & f_{xy} \\ f_{yx} & f_{yy} \end{pmatrix}}^{-1}
\begin{pmatrix} f_x \\ f_y \end{pmatrix}
\end{align}\tag{4-3}$$
特に、ヘッセ行列の導出にはひと工夫がいる。具体的には、局所的にヤコビ行列を構築している為、”LocalAD”同様、Type(hyper)メンバを適宜”On”、”Off”することでヘッセ行列を求めている。
subroutine newtonRosenbrock(N, xv, f, grad, hess, itr)
use hyperDN_mod
implicit none
integer,intent(in):: N
double precision,intent(inout):: xv(1:N)
double precision,intent(out):: f
integer,intent(out):: itr
type(hyper):: xh(1:N), fh
integer:: itrMax=10
double precision:: resi, upd(1:N), grad(1:N), hess(1:N,1:N)
double precision:: tol=1.0d-6
do itr=1,itrMax
! calculate gradient and diagonal on Hessian
xh(1)%f0 = xv(1); xh(2)%f0 = xv(2)
xh(1) = Hyper(f0=xv(1), f1=1.0d0, f2=1.0d0)
xh(2) = Hyper(f0=xv(2), f1=0.0d0, f2=0.0d0)
fh = funcRosenbrock(N, xh)
grad(1) = fh%f1; hess(1,1) = fh%f12
xh(1) = Hyper(f0=xv(1), f1=0.0d0, f2=0.0d0)
xh(2) = Hyper(f0=xv(2), f1=1.0d0, f2=1.0d0)
fh = funcRosenbrock(N, xh)
grad(2) = fh%f1; hess(2,2) = fh%f12
! calculate off-diagonal elements on Hessian
! Hessian is a symmetric matrix likely under continous derivative play case
xh(1) = Hyper(f0=xv(1), f1=1.0d0, f2=0.0d0)
xh(2) = Hyper(f0=xv(2), f1=0.0d0, f2=1.0d0)
fh = funcRosenbrock(N, xh)
hess(1,2) = fh%f12; hess(2,1) = hess(1,2)
! Solve A*x = B
call linsolve(hess, grad, upd)
f = fh%f0; resi = norm2(grad)
! Update solution vector: x_k+1 = x_k - f'(x_k) / f''(x_k)
xv(1) = xv(1) - upd(1); xv(2) = xv(2) - upd(2)
if (abs(resi) < tol) return
enddo
end subroutine
function funcRosenbrock(N, xh) result(res)
use hyperDN_mod
implicit none
integer, intent(in ) :: N
type(hyper), intent(in ):: xh(1:N)
type(hyper) :: res
double precision:: a=1.0d0, b=100.0d0
res = (a - xh(1))**2 + b * (xh(2) - xh(1)**2)**2
end function２．J.Fike氏の博士論文に有る例題を用いて、一階、二階偏導関数値の検証を試みる。偏導関数は論文に記載された式をそのまま流用している。
$$\begin{align}
f(x)&=\frac{e^x}{\sqrt{\sin^3{x}+\sin^3{x}}} \tag{4-4} \\ \\
f'(x)&=\frac{e^{x}(3\cos x+5\cos 3x+9\sin x+\sin 3x)}{8(\sin ^{3}x+\cos ^{3}x)^{\frac{3}{2}}} \tag{4-5} \\ \\
f^{\prime\prime}(x)&=\left(\frac{e^{x}}{64(\sin ^{3}x+\cos ^{3}x)^{\frac{5}{2}}}\right) \\ \\
&(130-12\cos 2x+30\cos 4x+12\cos 6x – 111\sin 2x+48\sin 4x+5\sin 6x) \\ \\
&= 4.4978 + 4.0534\epsilon1 + 4.0534\epsilon2 + 9.4631\epsilon1\epsilon2~~@ x= 1.5
\end{align}\tag{4-6}$$
比較結果は下表。ヘッセ行列の計算誤差は大きい気がするが、こんなもん何だろうか？
Hyperと自動微分との結果に差異はないことから、計算上の間違いはないものと思料。
| Item | Hyper-dual Numbers | Automatic Different’n | Analytical Different’n |
| Function | 4.49778 | <- | <- |
| Gradient | 4.05343 | <- | 4.03893 |
| Hessian | 9.4631 | <- | 22.0677 |
尚、同関数の最小化問題は制約なしニュートン法では計算途中で平方根内の逐次解が負になってしまい解を求めることはできない。
３．Python･Julia版では自動微分変数の結果を印字する際、型形式を含めて表記されウンザリしていたが、Fortranでは超双対数を用いた為か其の点スッキリしていて可読性は優れている。
D-2.Rosenbrock Minimization by Hyper-dual numbers example
Newton Method by Hyper-dual numbers
i = 1, x,y = -0.500000D+00 0.400000D+01, dx,y = 0.200267D-02 0.374800D+01
norm = 0.105854D+04
value= 0.140850D+04
grad = 0.747000D+03 0.750000D+03
hess = -0.129800D+04 0.200000D+03; 0.200000D+03 0.200000D+03
i = 2, x,y = -0.502003D+00 0.252003D+00, dx,y = -0.150080D+01 0.150681D+01
norm = 0.300481D+01
value= 0.225601D+01
grad = -0.300481D+01 -0.802138D-03
hess = 0.203607D+03 0.200801D+03; 0.200801D+03 0.200000D+03
i = 3, x,y = 0.998796D+00 -0.125480D+01, dx,y = -0.266645D-05 -0.225240D+01
norm = 0.100633D+04
value= 0.507329D+03
grad = 0.899872D+03 -0.450479D+03
hess = 0.170103D+04 -0.399518D+03; -0.399518D+03 0.200000D+03
i = 4, x,y = 0.998799D+00 0.997599D+00, dx,y = -0.120118D-02 -0.239948D-02
norm = 0.240236D-02
value= 0.144284D-05
grad = -0.240236D-02 -0.142182D-08
hess = 0.800079D+03 -0.399520D+03; -0.399520D+03 0.200000D+03
i = 5, x,y = 0.100000D+01 0.999999D+00, dx,y = -0.170747D-11 -0.144284D-05
norm = 0.645255D-03
value= 0.208177D-09
grad = 0.577134D-03 -0.288567D-03
hess = 0.802001D+03 -0.400000D+03; -0.400000D+03 0.200000D+03
i = 6, x,y = 0.100000D+01 0.100000D+01, dx,y = -0.444089D-15 -0.999201D-15
norm = 0.488579D-13
value= 0.142981D-29
grad = 0.435207D-13 -0.222045D-13
hess = 0.802000D+03 -0.400000D+03; -0.400000D+03 0.200000D+03
a = 1, b = 100, x0,y0 = -0.500000D+00 0.400000D+01-> x*,y* = 0.100000D+01 0.100000D+01 @itr = 6
@Rosenbrock value (x*) = 0.142981D-29
@Rosenbrock gradient (x*) = 0.435207D-13 -0.222045D-13
@Rosenbrock hessian (x*) = 0.802000D+03 -0.400000D+03 -0.400000D+03 0.200000D+03OpenAI
Python、Juliaでは文法書を精力的に読み込んでからコーディングを始めたものだが、Fortranはかなり知った言語であるため、今更Fortran2003文法書を読み込む気にもならない。そこで、超双対数Fortranコードを作成するに当たって積極的にopenAIを活用してみた。
Julia Partで触れたが、個々のHDN変数を合体して一元化する関数 ”assembleVariableHDN” の実装方法に置ける複数引数(multiple arguments)への対応に難儀する。chatGPTのアドバイスはWorkableではなく、逆の意味で考えさせるSuggestionであった。
具体的には、旧形式？のAsterisk(*)を用いるassumed-size array(サイズ引継配列)を組込むことで複数引数の入力を受付を可能とする提案であった。これ迄assumed-shape array(形状引継配列)を用いて配列入力を受ける関数を実装していた。尤も、”interface”構文にて総称名の機能を用いて関数、サブルーチンを多重定義することも可能だが、スマートでないので実装を諦めていた。
ある晩酌中に一日を振り替える中急遽思いついた手法でトライした処一発で解決する。俺は古い形式の配列構成子”/…/”しか使ってこなかったが(“/…/”はCompile Errorになる)、Fortran2003で導入された構成子”[…]”を思い出してトライした次第。関数に渡す複数引数を配列構成子にて渡す手法であった。因みに、関数 ”assembleVariableHDN”は配列入力仕様で構成されている。
加えてFortran2003機能強化で適切なサイズに自動的に再配列される機能が盛り込まれたことを踏まえ、”append”構文をモジュールに組込む。
又、呼出時に省略することを可能とする引数option属性付き引数を組込関数”present”と共に用いるという助言は有効であった。
!multiple arguments vs. vector argument
assembleVariableHDN([m1h, m2h, m3h]) !multiple x/m1h,m2h,m3h/
assembleVariableHDN(vh) !vector
!append
function append_iiv(lhs, rhs) result(res)
integer,intent(in) :: lhs
integer,dimension(:),allocatable,intent(in):: rhs
integer,dimension(size(rhs)+1) :: res
res = [lhs, rhs]
end function
...
!optional keyword
function variableHDN(f0, f1, f2, f12) result(res) !deprecated now
! public constructor
double precision,intent(in) :: f0
double precision,intent(in),optional:: f1, f2, f12
type(hyper) :: res
res%f0 = f0
if (present(f1)) then; res%f1 = f1
else; res%f1 = 1.0d0; endif
if (present(f2)) then; res%f2 = f2
else; res%f2 = 1.0d0; endif
if (present(f12)) then; res%f12 = f12
else; res%f12 = 1.0d0; endif
end function
!Please note that there is one fatal error in "variableHDN" function in term of application of the function.Closing
Python/Juliaに実装した自動微分法は算法の一種であるが、超双対数は数学的にも秀悦で美しく、Modern Fortranで対応できる実にしっくりくる手法である。特に、Fortranは配列処理がシンプル故逆に優れていると思う。初動を見誤り、Python、Juliaに自動微分を組込むことになったが、何故最初に超双対数に出わなかったのか残念である。
正に「井の中の蛙」
Python･Julia版ADモジュールを作成したものの、言語自体の使い勝手が悪くこれ以上追求する気になれない。PythonしてもJuliaにしてもRunさせてみなければ分からない。其の点、Fortranは型宣言が徹底しているので安心感がある。
俺的は言語使用に関して次のように位置づけている。
| 言語 | 位置づけ |
| Fortran | 科学技術計算に特化して一生使い続ける |
| Python | 日産Consult用として本格的に活用。ライブラリ・可視化機能も豊富で汎用性の高い言語故必要あれば使う。現状使用するネタはない |
| C++ | 若気の見栄で使い始める。何にでも対応できるが可読性が最悪。ボケが始まったらコードを見直すことは略不可能。今後使用することはなかろう |
| Arduino | 機器に接続しIOTとして用いている。必要とあれば今後も利用する |
| Julia | Python･C言語を意識した作りだが、中途半端な位置づけ。俺好みではない |
Modern Fortranは実に進化していることを改めて実感。俺にはこれ一つで足りうる環境に有る。
Juliaで実装した自動微分系”LocalAD”をFortranで作成することもありうる選択肢であるが、先ずiSMARTに超双対数を実装することに専念したい。ライフワークである自家製シミュレータ”iSMART”はヘッセ行列を必要としないが、演算速度にそれほど負荷が掛かるとも思わないので、comment outしないで”as is”で使用するつもりである。
一体あの4ヶ月の冬籠りは何だったんだろうか？
「冬籠り 無駄骨だった 微分かな」
以上
on the move!
References
- Yu, Wenbin and Blair, Maxwell, “DNAD, a Simple Tool for Automatic Differentiation of Fortran Codes Using Dual Numbers” (2013)
- J. A. Fike, Multi-Objective Optimization using Hyper-Dual Numbers: A Dissertation Submitted to the Department of Aeronaitics and Astronautics and the Committee on Graduate Studies of Stanford University in Partial Fulfillment of the Requirements for the Degree of Doctor of Philosophy, 2013
- Jeffrey A. Fike and Juan J. Alonso,~”The Development of Hyper-Dual Numbers for Exact Second-Derivative Calculations”, 49th AIAA Aerospace Sciences Meeting(2011)
- J. A. Fike and J. J. Alonso. The Development of Hyper-Dual Numbers for Exact Second Derivative Calculations. AIAA paper 2011-886, 49th AIAA Aerospace Sciences Meeting, January 4-7, 2011
- J. A. Fike and J. J. Alonso. The development of hyper-dual numbers for exact second-derivative calculations. In AIAA paper 2011-886, 49th AIAA Aerospace Sciences Meeting, 2011.
- J. A. Fike and J. J. Alonso. Automatic differentiation through the use of hyper-dual numbers for second derivatives. In S. Forth, P. Hovland, E. Phipps, J. Utke, A. Walther, T. J.Barth, M. Griebel, D. E. Keyes, R. M. Nieminen, D. Roose, and T. Schlick, editors, Recent Advances in Algorithmic Differentiation, volume 87 of Lecture Notes in Computational Science and Engineering, pages 163–173. Springer Berlin Heidelberg, 2012.
- J. A. Fike, S. Jongsma, J. J. Alonso, and E. van der Weide. Optimization with gradient and Hessian information calculated using hyper-dual numbers. In AIAA paper 2011-3807, 29th AIAA Applied Aerodynamics Conference, 2011.
〆