poly チュートリアル

このチュートリアルでは、luna-poly の多項式をある点で微分します。密な多項式と疎な多項式の値と傾きを求め、より大きな微分対象の式の中で多項式を使い、多項式に対してニュートン法を実行します。数学は poly の設計にあります。

クイックスタート

moon add Luna-Flow/autodiff@0.2.0
moon add Luna-Flow/luna-poly

import {
  "Luna-Flow/autodiff",
  "Luna-Flow/autodiff/poly",
  "Luna-Flow/luna-poly/immut/dense",
}

fn main {
  // p(x) = 1 + 2x + x^2, coefficients from degree 0 upwards
  let p = @dense.DensePolynomial::from_coefficients([1.0, 2.0, 1.0])
  let (value, slope) = @poly.dense_value_and_derivative_at(p, 3.0)
  println("p(3) = \{value}, p'(3) = \{slope}")
}

p(3) = 16, p'(3) = 8

日常的なタスク

導関数の表を作る

fn main {
  // p(x) = x^3 - 3x
  let p = @dense.DensePolynomial::from_coefficients([0.0, -3.0, 0.0, 1.0])
  for x in [-2.0, -1.0, 0.0, 1.0, 2.0] {
    println("p'(\{x}) = \{@poly.dense_derivative_at(p, x)}")
  }
}

p'(-2) = 9
p'(-1) = 0
p'(0) = -3
p'(1) = 0
p'(2) = 9

導関数 $3x^2 - 3$ は臨界点 $\pm 1$ で 0 になります。

整数係数の多項式の厳密な導関数

必要なのは Semiring だけなので、整数係数では厳密な結果が得られます：

fn main {
  // p(x) = 7 + 5x^4
  let p : @dense.DensePolynomial[Int] = @dense.DensePolynomial::from_coefficients([
    7, 0, 0, 0, 5,
  ])
  let (value, slope) = @poly.dense_value_and_derivative_at(p, 3)
  println("p(3) = \{value}, p'(3) = \{slope}")
}

p(3) = 412, p'(3) = 540

次数の間隔が大きい疎な多項式

高次の項が少ない多項式には疎な表現を使います。各項のコストは $O(\log e)$ です：

fn main {
  // p(x) = x^100 + 2x
  let p = @sparse.SparsePolynomial::from_array([([100U], 1.0), ([1U], 2.0)])
  let (value, slope) = @poly.sparse_univariate_value_and_derivative_at(p, 1.0)
  println("p(1) = \{value}, p'(1) = \{slope}")
}

p(1) = 3, p'(1) = 102

より大きな式の中の多項式

eval_dual は双対数の入力を受け取るので、連鎖律はそれを通して続きます。ここでは $x = 0.5$ での $\frac{d}{dx}\,p(\sin x)$ を求めます：

fn main {
  let p = @dense.DensePolynomial::from_coefficients([0.0, 0.0, 1.0]) // s^2
  let d = @autodiff.diff(x => @poly.eval_dual(p, x.sin()), 0.5)
  println("d/dx sin(x)^2 = \{d}")
  println("sin(2x)       = \{@math.sin(1.0)}")
}

d/dx sin(x)^2 = 0.8414709848078965
sin(2x)       = 0.8414709848078965

多項式に対するニュートン法

fn main {
  // p(x) = x^2 - 2
  let p = @dense.DensePolynomial::from_coefficients([-2.0, 0.0, 1.0])
  let mut x = 1.0
  for _ in 0..<6 {
    let (v, d) = @poly.dense_value_and_derivative_at(p, x)
    x = x - v / d
  }
  println("sqrt(2) = \{x}")
}

sqrt(2) = 1.414213562373095

さらに進む

形式的導関数との一致

luna-poly は導関数の多項式を構築することもできます。どちらの方法でも同じ数値が得られます。双対数上での評価なら 2 つ目の多項式が不要です：

fn main {
  let p = @dense.DensePolynomial::from_coefficients([5.0, -1.0, 0.0, 3.0])
  let formal = p.derivative().eval(2.0)
  let dual = @poly.dense_derivative_at(p, 2.0)
  println("formal \{formal}, dual \{dual}")
}

formal 35, dual 35

多変数多項式

疎な多項式のブリッジは 1 変数です。コンテキスト付きの多項式を 1 つの変数について微分するには、まず luna-poly/immut/context の ContextPolynomial::eval_partial で他の変数を評価し、残りの項をその変数の疎な多項式に変換してから、sparse_univariate_value_and_derivative_at を呼び出します。リポジトリの統合テスト src/tests/linalg_poly_test.mbt で変換の全体を示しています。

よくある落とし穴

係数の順序。 DensePolynomial::from_coefficients は定数項から順に $c_0, c_1, \dots$ を受け取ります。
第 2 変数があると中断します。 sparse_univariate_* は 1 つの値で評価します。別の変数の項があるとプログラムが停止します。
古い名前。 derivative_at と sparse_derivative_at は dense_ と sparse_univariate_ の関数のエイリアスであり、多変数版ではありません。
浮動小数点の桁落ち。 重根の近くでは、 $p'(x)$ は大きな項同士の小さな差になります。poly の設計の誤差上界を参照してください。

次のステップ

poly API にはすべての関数とそのエイリアスが記載されています。
poly の設計では双対数上のホーナー法とその誤差上界を導出しています。
luna-poly では多項式型を文書化しています。