Hyper-dual numbersによる二階偏微分の計算

Hyper-dual numbersと呼ばれる、実数を拡張した考えを導入すると導関数が計算できます。

あらかじめ、引数がHyper-dual numbersである時の数々の関数の定義を実装しておけば、その関数の組み合わせで作られる関数の一階導関数、二階導関数のほぼ厳密な答えを得ることが出来ます。
この考えに従う関数の微分方法は、Forward型の自動微分と呼ばれます。

  1. Dual number
    1. Dual numberと導関数
  2. Hyper-dual numbers
    1. Hyper-dual numbersの演算
    2. Hyper-dual numbersのプログラム
    3. ヘッセ行列
  3. 参考文献

Dual number


Dual numberという考えがあります[2,3,5]。

これは実数を拡張する、という考えで複素数に似た考え方です。

良く知られている実数の拡張方法の一つは、複素数です。
通常の実数に\(i^2=-1\)を満たす\(i\)という数を付加するのが複素空間です。

拡張の方法は何も\(i\)だけではありません。
例えば、\(\epsilon^2=0, (\epsilon\ne 0)\)を満たす\(\epsilon\)という数を付加することもできます。

この\(\epsilon\)は複素数ではありません
複素平面上でこの性質を満たす数は無いことからも、この\(\epsilon\)を追加するということは新しい方向への数の拡張です。
実際、複素平面上で\(\epsilon^2=0, (\epsilon\ne 0)\)を満たす数があるのかを探しても
\(
\begin{align}
(a+ib)^2=a^2-b^2+i2ab=0
\end{align}
\)

であるので、これを満たすのは\(a=b=0\)しか存在せず、複素平面上の数ではないことが分かります。

そんな新しい方向\(\epsilon\)の平面で定義された数をDual number(二重数、または双対数)と呼びます。

二重数\(a\)は実数部と非実数部から構成されており、
\(
a=a_0+a_1\epsilon
\)

と書くことが出来ます。ここで、\(a_0, a_1\)は実数であり、 \(\epsilon\)は、虚数単位\(i\)に倣って二重数単位とでも名付けておきましょう。

Dual numberと導関数


二重数は関数の導関数と大きな関係がある事を示しましょう。
テーラー展開を行います。
関数\(f\)の\(x\)周りの展開は刻み幅\(\Delta\)を用いて
\(\displaystyle
f(x+\Delta)=f(x)+f'(x) \Delta+\frac{1}{2!}f^{\prime\prime}(x)\Delta^2+\cdots
\)

と記述することが出来ます。もし仮に\(\Delta\)が純非実数だとしましょう。すなわち、\(\Delta=h\epsilon\)とします。ここで、\(h\)は実数です。テーラー展開の式に代入すれば、
\(
\begin{eqnarray}
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
\end{eqnarray}
\)

となるわけです。注記しますが、上の式は\(h\)の1次で打ち切っているのではなく、厳密にイコールが成り立っています。
この式が言っているのは、関数\(f(a), (a\text{は二重数})\)を計算し、その非実数部を\(h\)で割ると関数の導関数になっているということです。

二重数の非実数部を取り出す関数\(\text{Dual}\)を定義します。
すると、導関数は
\(\displaystyle
f'(x)=\frac{1}{h}\text{Dual}(f(x+h\epsilon))
\)

として得られます。
言葉で書けば、二重数空間に拡張した\(f(x+h\epsilon)\)を計算すると、その非実数部に導関数が現れる、ということです。

二重数のプログラムはJeffrey Fikeさんによる[4]にありますので、そちらをご参考ください。

Hyper-dual numbers


さて、ここまでで二重数の概念を簡単に説明し、新しい実数の拡張を行いました。
Dual Numberのままでは高階導関数は得られません。なぜなら、導関数の二次以降は二重数の性質\(\epsilon^2=0\)によってゼロになるからです。

そこで、若干工夫します。
Dual Numberを二種類用意することで二階微分を得ることが出来ます[2,3]。
すなわち、
\(\begin{gather}
a=a_0+a_1\epsilon_1+a_2\epsilon_2+a_3\epsilon_1\epsilon_2\\
\epsilon_1^2=\epsilon_2^2=(\epsilon_1\epsilon_2)^2=0,~\epsilon_1\ne 0,~~\epsilon_2\ne 0,~~\epsilon_1\epsilon_2 \ne 0
\end{gather}
\)

と実数を拡張します。この様に拡張した数\(a\)をHyper-dual numbersと呼びます[2,3]。

Hyper-dual numbersの演算


Hyper-dual numbersである
\(
\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}
\)

を用意します。和、積、商はそれぞれ
\(
\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+a_3b_0)\epsilon_1\epsilon_2\\
\frac{1}{a}&=\frac{1}{a_0}-\frac{a_1}{a_0^2}\epsilon_1-\frac{a_2}{a_0^2}\epsilon_2+\left(-\frac{a_3}{a_0^2}+\frac{2a_1a_2}{a_0}\right)\epsilon_1\epsilon_2
\end{align}
\)
と定義されます[2]。任意の関数は
\(
f(x)=f(x_0)+x_1f'(x_0)\epsilon_1+x_2f'(x_0)\epsilon_2
+\left(x_3f'(x_0)+x_1x_2f^{\prime\prime}(x_0)\right)\epsilon_1\epsilon_2
\)

として計算することが出来ます。なので、この結果から\(\epsilon_1\epsilon_2\)の係数として関数の二階微分が得られます。

Hyper-dual numbersのプログラム


実際にプログラムを組みましょう。
Hyper-dual numbersの型を持つ変数はFortran90では定義できません。
なので、構造体を利用して自分で型と、演算を定義します。

Fortran90のプログラムは以下の通りになるかと思います。
基本的な四則演算、基本的な初等関数のHyper-dual numbersの演算をモジュールとして書いています。

上のモジュールと下のメインプログラムを一緒にコンパイルすることにより、関数の二階微分が得られます。

例として
2変数関数
\(\displaystyle
f(x,y)=\frac{\text{ln}(xy^2)e^x}{\sqrt{\sin^3{x}+\sin^3{y}}}
\)

の偏微分
\(\begin{align}
\frac{\partial f}{\partial x},~~\frac{\partial f}{\partial y},~~\frac{\partial^2 f}{\partial x\partial y}
\end{align}
\)

を得ることを考えます。
\(
f(x+1\epsilon_1,y+1\epsilon_2)
\)

を計算すると
\(
f(x+1\epsilon_1,y+1\epsilon_2)=a_0+a_1\epsilon_1+a_2\epsilon_2+a_3\epsilon_1\epsilon_2
\)

のように実数係数\(a_0, a_1, a_2, a_3\)が得られます。
すると、
\(
\begin{align}
\frac{\partial f}{\partial x}=a_1 \\
\frac{\partial f}{\partial y}=a_2 \\
\frac{\partial^2 f}{\partial x\partial y} =a_3
\end{align}
\)

として偏微分が得られます。
ちなみに、二階微分が欲しい場合は
\(
f(x+1\epsilon_1+1\epsilon_2, y)
\)

を考えると
\(
\begin{align}
\frac{\partial f}{\partial x}=a_1 \\
\frac{\partial f}{\partial x}=a_2 \\
\frac{\partial^2 f}{\partial x^2} =a_3
\end{align}
\)

として得られます。\(a_1,~a_2\)はどちらを採用しても構いません。

プログラムでは変数の型Hyperdualを持つ入力変数を\(\text{xH,yH}\),出力を\(\text{rH}\)と置いています。

program main
  use Hyperdualmod
  implicit none
  type(Hyperdual)::xH,yH,rH
 
  xH%x0 = 0.3d0 ! real part
  xH%x1 = 1d0   ! unreal part \epsilon_1
  xH%x2 = 0d0   ! unreal part \epsilon_2
  xH%x3 = 0d0   ! unreal part \epsilon_1\epsilon_2

  yH%x0 = 0.4d0 ! real part
  yH%x1 = 0d0   ! unreal part \epsilon_1
  yH%x2 = 1d0   ! unreal part \epsilon_2
  yH%x3 = 0d0   ! unreal part \epsilon_1\epsilon_2

  write(6,'(4f23.16)')xH%x0,xH%x1,xH%x2,xH%x3
  write(6,'(4f23.16)')yH%x0,yH%x1,yH%x2,yH%x3  
 
  rH = log(xH*yH**2)*exp(xH)/sqrt(sin(xH)**3+cos(yH)**3)
  !rH = asin(2d0*xH)*acos(yH)/atan(xH*yH)
  !rH = xH**yH
 
  write(6,'(4f23.16)')rH%x0,rH%x1,rH%x2,rH%x3  

  stop
end program main

ヘッセ行列


二階偏微分の計算が出来たので、ヘッセ行列が簡単に計算できます。
ルーチンを作れば、以下の通りになるかと思います。
下のプログラムは3変数関数
\(\displaystyle
f(x,y,z)=\exp(xy)\tan(z)
\)

の\(x=-2,~y=3,~z=1\)におけるヘッセ行列を計算します。

program main
  use Hyperdualmod
  implicit none
  integer::N
  double precision::x,y,z,f
  double precision,allocatable::nabla(:),Hesse(:,:)
  double precision,allocatable::w(:)
  external::func

  N=3

  allocate(nabla(1:N),Hesse(1:N,1:N))
  nabla=0d0
  Hesse=0d0
  allocate(w(1:N))
  w=0d0

  x = -2d0
  y =  3d0
  z =  1d0
  w(1) = x
  w(2) = y
  w(3) = z
  call Hessian(N,w,func,nabla,Hesse)

  f=exp(x*y)*tan(z)
  write(6,'(1e24.16)')nabla(1)
  write(6,'(1e24.16)')y*f
  write(6,'(1e24.16)')nabla(2)
  write(6,'(1e24.16)')x*f
  write(6,'(1e24.16)')nabla(3)
  write(6,'(1e24.16)')exp(x*y)/(cos(z)**2)
  write(6,'(3e24.16)')Hesse(1,1:3)
  write(6,'(3e24.16)')f*y**2, f*(1+x*y), f*y/(cos(z)*sin(z))
  write(6,'(3e24.16)')Hesse(2,1:3)
  write(6,'(3e24.16)')f*(1+x*y), x**2*f, f*x/(cos(z)*sin(z))
  write(6,'(3e24.16)')Hesse(3,1:3)
  write(6,'(3e24.16)')f*y/(cos(z)*sin(z)), f*x/(cos(z)*sin(z)), f*2d0/(cos(z)**2)

  stop
end program main

subroutine func(N,x,f)
  use Hyperdualmod
  implicit none
  integer::N
  type(Hyperdual),intent(in)::x(1:N)
  type(Hyperdual),intent(out)::f
 
  f = exp(x(1)*x(2))*tan(x(3))
 
  return
end subroutine func

subroutine Hessian(N,x,func,nabla,Hesse)
  use Hyperdualmod
  implicit none
  integer,intent(in)::N
  double precision,intent(in)::x(1:N)
  double precision,intent(out)::nabla(1:N)
  double precision,intent(out)::Hesse(1:N,1:N)
  external::func

  integer::i,j,k
  type(Hyperdual),allocatable::xH(:)
  type(Hyperdual)::rH
 
  allocate(xH(1:N))
  do i=1,N
     xH(i)%x0 = 0d0
     xH(i)%x1 = 0d0
     xH(i)%x2 = 0d0
     xH(i)%x3 = 0d0
  enddo
 
  do i=1,N
     xH(i)%x0 = x(i)
     xH(i)%x3 = 0d0
  enddo

  do i=1,N
     do k=1,N
        xH(k)%x1 = 0d0
        xH(k)%x2 = 0d0
     enddo
     xH(i)%x1 = 1d0
     xH(i)%x2 = 1d0
     
     call func(N,xH,rH)
     nabla(i) = rH%x1
     Hesse(i,i) = rH%x3
  enddo

  do i=1,N
     do j=i+1,N
        do k=1,N
           xH(k)%x1 = 0d0
           xH(k)%x2 = 0d0
        enddo
        xH(i)%x1 = 1d0
        xH(i)%x2 = 0d0
        xH(j)%x1 = 0d0
        xH(j)%x2 = 1d0
        call func(N,xH,rH)
        Hesse(i,j) = rH%x3
        Hesse(j,i) = Hesse(i,j)
     enddo
  enddo  
 
  return
end subroutine Hessian

Hyperdual.f90にモジュールを入れ、メインプログラムをmain.f90に入れたとすると、以下の結果を得ます。

> gfortran Hyperdual.f90 main.f90  
> ./a.out
  0.1158128336233602E-01
  0.1158128336233602E-01
 -0.7720855574890680E-02
 -0.7720855574890680E-02
  0.8491012233306163E-02
  0.8491012233306162E-02
  0.3474385008700806E-01 -0.1930213893722670E-01  0.2547303669991849E-01
  0.3474385008700806E-01 -0.1930213893722670E-01  0.2547303669991849E-01
 -0.1930213893722670E-01  0.1544171114978136E-01 -0.1698202446661233E-01
 -0.1930213893722670E-01  0.1544171114978136E-01 -0.1698202446661233E-01
  0.2547303669991849E-01 -0.1698202446661233E-01  0.2644793608458059E-01
  0.2547303669991849E-01 -0.1698202446661233E-01  0.2644793608458058E-01

実行結果の奇数行目はHyper-dual Numberによる計算結果、偶数行目は解析解を表します。
また、6行目までは一階微分、7行目以降はヘッセ行列を表します。

参考文献


[1]関数
\(\displaystyle
f(x,y)=\frac{\text{ln}(xy^2)e^x}{\sqrt{\sin^3{x}+\sin^3{y}}}
\)

の偏微分\(\frac{\partial^2 f}{\partial x\partial y}\)
の計算
https://www.wolframalpha.com/input/?i=D%5Bln(x*y%5E2)*e%5E(x)%2Fsqrt(sin%5E3(x)%2Bcos%5E3(y)),+y,x%5D

[2]Jeffrey A. Fike and Juan J. Alonso,~”The Development of Hyper-Dual Numbers for Exact Second-Derivative Calculations”, 49th AIAA Aerospace Sciences Meeting(2011)
http://adl.stanford.edu/hyperdual/fike_aiaa-2011-886_slides.pdf,
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.http://adl.stanford.edu/hyperdual/Fike_AIAA-2011-886.pdf

[3]JeffreyA.Fike,~”Derivative Calculations Using Hyper-Dual Numbers”, Sandia National Laboratories (2016)
https://www.osti.gov/servlets/purl/1368722

[4]Jeffrey Fike, Aerospace Design Lab, http://adl.stanford.edu/hyperdual/

[5]松本佳彦, 新しい数をつくる, (2018) http://ymatz.net/assets/docs/20180629-jtpa-slide-mod


コメントを残す

メールアドレスが公開されることはありません。 * が付いている欄は必須項目です