Fortran で自動微分

前置き　

暑いので肌色多めの涼しげな画像を眺めつつ BGM がわりに Youtube 動画を垂れ流していたところ、MIT 公開講座で Julia 言語で数値微分の話になって、面白くないので全く聞いていなかったのですが、数値の組に適当な代数則を与えてニュートン法すると、ついでに微分値が求まるという話になって、「マジか！？」となったので Fortran で再現してみます。

1. Alan Edelman and Julia Language
   https://youtu.be/rZS2LGiurKY?t=972

無限小解析による数値微分？

以前、実関数の定義域・値域を複素数に拡張して正則な関数ならばどの方向から微分しても微係数が一致することを利用して虚軸側から近づくことで桁落ちを防ぐ微分法を紹介しました。
数値微分で定義域を複素数に拡張して

ここでの方法はそれより先を行くようで、ざっとみたところ無限小解析を応用して、数値微分なれども桁落ちによる誤差のみならず、無限小極限を取るところを有限で打ち切る誤差もなくすことができるようです。計算機の方では dual number と呼んでいるようです。

形式的には、分数や複素数は、二つの実数の組にそれぞれ加減乗除に相応の演算則を与えることで計算できますが、ここでも二つの実数の組に加減乗除の演算則を与えています。そうして、二数のうちの普通の実数の演算則を満たす方で x を関数値 f(x) にする計算を行うと、もう一方の数が f'(x) の値になっているというからくりになっています。

無限小解析

無限小 $\epsilon$ 但し $\epsilon^2=0$ となる数を加えて代数系を作ると、一般の数は、複素数みたいな普通の実数 $a, b$ 二数の組で $z=a+b\epsilon$ のように書き表されます。

加減乗除を考えると
$$
(x+x'\epsilon) + (y+y'\epsilon) = (x+y) + (x'+y')\epsilon
$$
$$
(x+x'\epsilon) - (y+y'\epsilon) = (x-y) + (x'-y')\epsilon
$$
$$
(x+x'\epsilon) * (y+y'\epsilon) = (xy) + (x'y+xy')\epsilon
$$
$$
(x+x'\epsilon) / (y+y'\epsilon) = (x/y) + \frac{(x'y-xy')}{y^2}\epsilon
$$

Taylor 展開を一次で打ち切ったようなものと思えば、規則から見て無限小にかかる係数の方を微分値だと考えることが出来ます。

Fortran プログラム

Fortran90 で導入された演算子オーバーロードを用いて、二数の組の派生型に上の演算則を定義し、Youtube 例題でやっていたニュートン法に適用してみることにします。

ソース

    module deriv_m
        implicit none
        integer, parameter :: kd = kind(0.0d0)
        
        type :: pair_t
            real(kd) :: x, dx
        end type pair_t
        
        interface operator(+)
            module procedure :: add_p
        end interface
        
        interface operator(-)
            module procedure :: sub_p
        end interface
        
        interface operator(*)
            module procedure :: mul_p
        end interface

        interface operator(/)
            module procedure :: div_p
        end interface
    contains
    
        pure elemental type (pair_t) function add_p(a, b) result(res)
             type(pair_t), intent(in) :: a, b
             res%x  = a%x  + b%x
             res%dx = a%dx + b%dx
        end function add_p
    
        pure elemental type (pair_t) function sub_p(a, b) result(res)
             type(pair_t), intent(in) :: a, b
             res%x  = a%x  - b%x
             res%dx = a%dx - b%dx
        end function sub_p
    
        pure elemental type (pair_t) function mul_p(a, b) result(res)
             type(pair_t), intent(in) :: a, b
             res%x  = a%x  * b%x
             res%dx = a%dx * b%x + a%x * b%dx
        end function mul_p
    
        pure elemental type (pair_t) function div_p(a, b) result(res)
             type(pair_t), intent(in) :: a, b
             res%x  = a%x  / b%x
             res%dx = (a%dx * b%x - a%x * b%dx) / b%x**2
        end function div_p
    
    end module deriv_m

    
    
    
    program auto_derivative
        use :: deriv_m
        implicit none
        
        type(pair_t) :: a, x 
        integer :: i
        
! Newton's method for SQRT
        a = pair_t(20.0_kd, 1.0_kd) 
        print *, 'Square Root of ', a%x

        x = pair_t(1.0_kd, 0.0_kd)
        do i = 1, 10
            x = pair_t(0.5_kd, 0.0_kd) * (x + a / x)
            print *, x
        end do    
        print *, 'Theoretical values'
        print *, sqrt(a%x), 0.5_kd / sqrt(a%x)
        print *
        
! Newton's method for Qubic Root
        a = pair_t(20.0_kd, 1.0_kd) 
        print *, 'Qubic Root of ', a%x
        
        x = pair_t(1.0_kd, 0.0_kd)
        do i = 1, 10
            x = (pair_t(2.0_kd, 0.0_kd) * x + a / (x * x)) / pair_t(3.0_kd, 0.0_kd) 
            print *, x
        end do    
        print *, 'Theoretical values'
        print *, a%x**(1.0_kd / 3.0_kd), a%x**(-2.0_kd / 3.0_kd) / 3.0_kd
        
    end program auto_derivative

実行結果

20 のルートと三乗根を求めています。副産物として √20 と 20^(1/3) の地点での導関数の値がそれぞれ求まっています。　

a^n を求める時、微係数の方はつねに 1 を与えるだけでいいようで、結果はこの値に比例しています。解の反復の初期値は何でもいいようです。

最後に理論値も書き出して比較をしています。

 Square Root of    20.0000000000000
   10.5000000000000       0.500000000000000
   6.20238095238095       0.252267573696145
   4.71347454528837       0.141172045622539
   4.47831444547438       0.113122077400955
   4.47214021706570       0.111805110525440
   4.47213595500161       0.111803398876570
   4.47213595499958       0.111803398874989
   4.47213595499958       0.111803398874989
   4.47213595499958       0.111803398874989
   4.47213595499958       0.111803398874989
 Theoretical values
   4.47213595499958       0.111803398874989

 Qubic Root of    20.0000000000000
   7.33333333333333       0.333333333333333
   5.01285583103765       0.217150847316137
   3.60720453602482       0.135047257987521
   2.91715357089440       7.728606445168938E-002
   2.72818093146380       4.918374657634904E-002
   2.71448693395070       4.527873410110104E-002
   2.71441761836499       4.524029554360360E-002
   2.71441761659491       4.524029360991512E-002
   2.71441761659491       4.524029360991511E-002
   2.71441761659491       4.524029360991511E-002
 Theoretical values
   2.71441761659491       4.524029360991511E-002

まとめ

正直、まだスッキリ分かった気がしないですがｗ ゆっくり考えてゆきたいと思います。

複素数や分数は、派生型と加減乗除の演算子オーバーロードの教材向きですが、自明すぎてやや面白みに欠けます。無限小解析への適用での自動微分はとても興味深い題材だと感じられました。歴史の遺物か色物というイメージで無限小解析をみていましたが、数値計算にこんな形で応用できるとは考えもよりませんでした。