Luna Flow
日本語

dual チュートリアル

このチュートリアルでは、Dual[T] を使って手動で導関数を計算します。入力をシードし、通常の算術を実行し、接成分から導関数を読み取ります。最後には、多変数の式を任意の方向に微分し、通常の数と双対数の両方に使える 1 つのジェネリック関数を書き、除算と平方根の失敗をデータとして扱えるようになります。数学は dual の設計 にあります。

クイックスタート

プロジェクトにモジュールを追加します：

moon add Luna-Flow/autodiff@0.2.0

moon.pkg でルートパッケージをインポートします。これは Dual と以下で使うトレイトを再エクスポートしています：

import {
  "Luna-Flow/autodiff",
}

x=1.5x = 1.5 を変数としてシードし、f(x)=x3+2xf(x) = x^3 + 2x を評価します：

fn main {
  let x = @autodiff.Dual::variable(1.5)
  let two = @autodiff.Dual::constant(2.0)
  let y = x * x * x + two * x
  println("f(1.5)  = \{y.value()}")
  println("f'(1.5) = \{y.tangent()}")
}
f(1.5)  = 6.375
f'(1.5) = 8.75

接成分は f′(1.5)=3⋅1.52+2=8.75f'(1.5) = 3 \cdot 1.5^2 + 2 = 8.75 です。

日常的なタスク

計算に定数を混ぜる

微分変数でない値はすべて、接成分 0 の定数として入ります。リテラルは双対数と直接混ぜることができないので、Dual::constant で包みます：

fn main {
  let rate = 0.25
  let x = @autodiff.Dual::variable(2.0)
  let y = @autodiff.Dual::constant(rate) * x.exp()
  println("d/dx 0.25 e^x at 2 = \{y.tangent()}")
}
d/dx 0.25 e^x at 2 = 1.8472640247326626

初等関数を使う

sqrt、exp、exp2、ln、log2、log10、sin、cos、tan は Dual[T] のメソッドであり、連鎖律を自動的に適用します：

fn main {
  let x = @autodiff.Dual::variable(1.0)
  let y = x.sin().exp() // e^(sin x)
  println("value      = \{y.value()}")
  println("derivative = \{y.tangent()}")
  println("cos(1) e^(sin 1) = \{@math.cos(1.0) * @math.exp(@math.sin(1.0))}")
}
value      = 2.319776824715853
derivative = 1.253380767493447
cos(1) e^(sin 1) = 1.253380767493447

2 行目と 3 行目は一致します。接成分は cos⁡(1) esin⁡1\cos(1)\,e^{\sin 1} です。

選んだ方向に微分する

入力が複数ある場合、シードする接成分が方向を決めます。(x,y)=(2,3)(x, y) = (2, 3) を接成分 (1,0)(1, 0) でシードすると ∂f/∂x\partial f/\partial x が得られ、接成分 (v1,v2)(v_1, v_2) では方向微分 ∇f⋅v\nabla f \cdot v が得られます：

fn f(x : @autodiff.Dual[Double], y : @autodiff.Dual[Double]) -> @autodiff.Dual[Double] {
  x * y + x.sin()
}

fn main {
  let dx = f(@autodiff.Dual::new(2.0, 1.0), @autodiff.Dual::new(3.0, 0.0))
  let dy = f(@autodiff.Dual::new(2.0, 0.0), @autodiff.Dual::new(3.0, 1.0))
  let dv = f(@autodiff.Dual::new(2.0, 0.5), @autodiff.Dual::new(3.0, -1.0))
  println("df/dx = \{dx.tangent()}")
  println("df/dy = \{dy.tangent()}")
  println("grad . (0.5, -1) = \{dv.tangent()}")
}
df/dx = 2.5838531634528574
df/dy = 2
grad . (0.5, -1) = -0.7080734182735712

勾配全体やヤコビ行列については、linalg チュートリアル がこのシードを代わりに行います。

有限差分と比較する

中心差分には刻み幅が必要で、桁を失います。双対数の結果は丸めを除いて正確です：

fn g(x : Double) -> Double {
  @math.exp(@math.sin(x))
}

fn main {
  let exact = @math.cos(1.0) * g(1.0)
  let dual = @autodiff.Dual::variable(1.0).sin().exp().tangent()
  for h in [1.0e-2, 1.0e-5, 1.0e-8] {
    let central = (g(1.0 + h) - g(1.0 - h)) / (2.0 * h)
    println("h = \{h}: error \{(central - exact).abs()}")
  }
  println("dual: error \{(dual - exact).abs()}")
}
h = 0.01: error 0.00006752362465478612
h = 0.00001: error 7.4792172455318e-11
h = 1e-8: error 9.450921378828525e-9
dual: error 0

有限差分の誤差は hh とともにまず減少し、その後桁落ちが支配的になると再び増加します。これは dual の設計 で導いたとおりです。

通常の数と双対数のための 1 つの関数を書く

関数を、それが必要とするトレイトに対して書きます。すると同じコードが Double 上で、そして導関数のためには Dual[Double] 上で動作します。整数定数は IntegralHomomorphism::from_integral から得ます：

fn[T : @autodiff.Ring + @autodiff.IntegralHomomorphism + @autodiff.Trigonometric] h(
  x : T,
) -> T {
  let three : T = @autodiff.IntegralHomomorphism::from_integral(3)
  three * x * x + @autodiff.Trigonometric::cos(x)
}

fn main {
  println("h(0.5)  = \{h(0.5)}")
  let d = h(@autodiff.Dual::variable(0.5))
  println("h(0.5)  = \{d.value()} (dual value)")
  println("h'(0.5) = \{d.tangent()}")
}
h(0.5)  = 1.6275825618903728
h(0.5)  = 1.6275825618903728 (dual value)
h'(0.5) = 2.520574461395797

双対数上で計算した値は、Double 上で計算した値とビット単位で一致します。接成分は 0.50.5 における 6x−sin⁡x6x - \sin x です。

除算と平方根の失敗を報告する

div_checked と sqrt_checked は、無限大や NaN の代わりに arithmetic のエラーを持つ Result を返します：

fn main {
  let ctx = @autodiff.ArithmeticContext::new(53)
  for a in [4.0, 0.0, -1.0] {
    match @autodiff.Dual::variable(a).sqrt_checked(ctx) {
      Ok(r) => println("sqrt(\{a}): value \{r.value()}, derivative \{r.tangent()}")
      Err(e) =>
        println(
          "sqrt(\{a}): domain error \{e.is_domain_error()}, division by zero \{e.is_division_by_zero()}",
        )
    }
  }
}
sqrt(4): value 2, derivative 0.25
sqrt(0): domain error false, division by zero true
sqrt(-1): domain error true, division by zero false

00 では平方根は存在しますがその導関数は存在しないので、接成分の除算が失敗します。

さらに進む

入れ子による 2 階導関数

Dual[T] はジェネリックなので、T 自体を双対数にできます。導関数を微分するには内側の関数がジェネリックである必要があり、forward チュートリアル では diff を使ってこれを示しています。手動で行うと次のようになります：

fn[T : @autodiff.Ring] cube(x : T) -> T {
  x * x * x
}

fn main {
  // outer variable: tangent 1 on the outer level
  let outer : @autodiff.Dual[Double] = @autodiff.Dual::variable(2.0)
  // inner variable: the outer number, seeded with tangent 1 on the inner level
  let inner = @autodiff.Dual::new(outer, @autodiff.Dual::constant(1.0))
  let y = cube(inner)
  println("f   = \{y.value().value()}")
  println("f'  = \{y.tangent().value()}")
  println("f'' = \{y.tangent().tangent()}")
}
f   = 8
f'  = 12
f'' = 12

2 つのレベルは異なる型なので、内側と外側の導関数の接成分が混ざることはありません。

独自のスカラー型

呼び出すメソッドの境界を満たす任意の T が使えます。Int のような環だけの型でも、すでに積の規則をサポートしています：

fn main {
  let x : @autodiff.Dual[Int] = @autodiff.Dual::variable(5)
  let y = x * x * x
  println("d/dx x^3 at 5 = \{y.tangent()}")
}
d/dx x^3 at 5 = 75

独自の型で exp や sin を通して微分するには、その型に Luna-Flow/arithmetic の Exponential や Trigonometric を実装してください。

よくある落とし穴

  • リテラルは双対数ではありません。 x が Dual[Double] のとき、x * 2.0 はコンパイルできません。x * @autodiff.Dual::constant(2.0) と書くか、ジェネリックなコードでは IntegralHomomorphism::from_integral(2) を使ってください。
  • 比較は比較したものしか見ません。 Dual[T] には < がありません。x.value() を比較してください。そのとき導関数は選ばれた分岐の導関数になるので、分岐で書かれた abs のような関数は、折れ目で片側導関数を持ちます。
  • 等価性は接成分も含みます。 Dual::new(1.0, 0.0) == Dual::new(1.0, 1.0) は false です。
  • チェックなしの演算は Double に従います。 y.value() == 0.0 での x / y、負の数の ln、0 での sqrt は、接成分に無限大や NaN を生じます。それが許されない場合はチェック付きの形式を使ってください。
  • 非常に小さい除数。 商の規則は c2c^2 で割りますが、これは ∣c∣<1.5×10−162|c| < 1.5 \times 10^{-162} でアンダーフローします。除算の前にスケーリングし直してください。
  • 表示できません。 Dual[T] には Show がありません。value() と tangent() を表示するか、Debug の debug_inspect と @debug.to_string を使ってください。

次のステップ