poly 設計

このページでは、luna-poly の多項式を双対数上で評価するとその導関数が得られる理由、密な評価順序と疎な評価順序が何を計算するか、結果の精度、そしてブリッジが 1 変数に限定されている理由を説明します。

設計目標

luna-poly が既に持っている評価アルゴリズムを再利用することで、新しい多項式表現も記号的な処理もなしに、luna-poly のユーザーにある点での多項式の値と導関数を提供します。

数学的背景

双対数上での評価が導関数を与える

可換半環 $R$ 上の $p(x) = \sum_{k=0}^{n} c_k x^k$ について、dual の設計は次を示しています

p(a + b\varepsilon) = p(a) + p'(a)\,b\,\varepsilon, \qquad p'(x) = \sum_{k=1}^{n} k\,c_k x^{k-1},

ここで $k\,c_k$ は $c_k$ を $k$ 回足したものを意味します。この恒等式は代数的なので、整数係数の多項式では厳密に、浮動小数点の多項式では丸めの範囲内で成り立ちます。極限も除算も必要ありません。

双対数上のホーナー法

DensePolynomial::eval は $r_{n+1} = 0$ 、 $k = n, \dots, 0$ について $r_k = r_{k+1}\,x + c_k$ を計算し、 $r_0 = p(x)$ を返します。 $x = a + b\varepsilon$ 、持ち上げた係数 $c_k + 0\varepsilon$ 、 $r_k = p_k + q_k\varepsilon$ とすると、双対数の積と和から次が得られます

\begin{aligned} p_k &= p_{k+1}\,a + c_k, \\ q_k &= p_{k+1}\,b + a\,q_{k+1}, \qquad p_{n+1} = q_{n+1} = 0 . \end{aligned}

$b = 1$ の場合、これは多項式とその導関数を同時に評価する古典的な方式です。¹1 D. E. Knuth, The Art of Computer Programming, vol. 2, 3rd ed., section 4.6.4. 帰納法により $p_k = \sum_{i \ge k} c_i a^{i-k}$ かつ $q_k = b\sum_{i > k} (i - k)\,c_i a^{i-k-1}$ なので、 $p_0 = p(a)$ かつ $q_0 = p'(a)\,b$ です。

累乗による疎な項の評価

SparsePolynomial::eval は格納された項について $c\,x^e$ を合計し、 $x^e$ を二分累乗法で計算します。双対数上ではすべての積が積の規則を適用するので、 $x^e$ は $a^e + e\,a^{e-1} b\,\varepsilon$ になり（dual の設計の二項恒等式）、項の和は再び $p(a) + p'(a)\,b\,\varepsilon$ を与えます。

設計上の判断

記号的に微分せずに評価する

問題。 ユーザーは各点での $p'(x)$ を必要とします。

選択肢。 DensePolynomial::derivative で導関数の多項式を構築して評価する、または Dual[T] 上で $p$ を評価する。

選択。 Dual[T] 上で評価します。1 回のパスで $p(x)$ と $p'(x)$ を返し、Semiring だけを必要とし（luna-poly の形式的導関数は $k\,c_k$ を作るために NatHomomorphism も必要とします）、導関数の階数ごとに 2 つ目の多項式を確保することもなく、eval_dual を通じて他の双対数計算と合成できます。導関数の多項式そのものが欲しい場合は形式的導関数が引き続き適した手法であり、テストスイートは両者が一致することを確認しています。

`luna-poly` の評価を再利用する

ブリッジは係数を Dual::constant で持ち上げ、通常のコンストラクタで DensePolynomial[Dual[T]] または SparsePolynomial[Dual[T]] を構築し、eval を呼び出します。ホーナー法を再実装しないので、luna-poly の評価（および 0 係数の正規化。これが T : Eq が必要な理由です）への変更はここにもそのまま反映されます。

疎な多項式は 1 変数

問題。 SparsePolynomial は多変数ですが、偏導関数には変数の選択が必要であり、変数の同一性は luna-poly のコンテキスト（VariableContext、ContextPolynomial）の中にあります。

選択。 疎な多項式のブリッジは単一の代入 [x] で評価し、sparse_univariate_* と命名されています。他の変数を使う多項式では eval が中断します。多変数の API は変数コンテキストを受け取って勾配を返す必要があり、今後の課題です。コンテキスト付きの多項式を持つ呼び出し側は、リポジトリの統合テストのように、まず他の変数を部分評価しておくことができます。

互換用の名前

最初のリリースでは関数名を derivative_at、value_and_derivative_at、eval_sparse_dual、sparse_derivative_at、sparse_value_and_derivative_at としていました。dense_ と sparse_univariate_ の名前は、どの表現とどの種類の導関数を意味するかを示します。既存のコードが動作し続けるよう、古い名前は単純なエイリアスとして残しています。

正しさと不変条件

厳密な環での厳密性

厳密な算術を持つ T（Int、BigInt、厳密な有理数）では、上の恒等式により結果はちょうど $p(x)$ と $p'(x)$ になります。

密な評価の丸め誤差

$\mathrm{fl}(x \circ y) = (x \circ y)(1 + \delta)$ 、 $|\delta| \le u$ と、積の補題 $\prod (1 + \delta_i)^{\pm1} = 1 + \theta_k$ 、 $|\theta_k| \le \gamma_k = ku/(1 - ku)$ を用います。²2 N. J. Higham, Accuracy and Stability of Numerical Algorithms, 2nd ed., SIAM, 2002, Lemma 3.1 and chapter 5.

値。射影準同型により、 $\hat p_0$ は通常のホーナー法の結果です。各ステップは累積和に $(1 + \delta_{\times})(1 + \delta_{+})$ を掛け、1 回丸めて $c_k$ を加えるので、 $k < n$ では $c_k a^k$ は $2k + 1$ 個の因子を持ちます。最初のステップ $0 \cdot a + c_n$ は厳密なので、最高次の係数は $2n$ 個の因子を持ちます：

\hat p_0 = \sum_{k=0}^{n} c_k a^k (1 + \theta^{(k)}), \quad |\theta^{(k)}| \le \gamma_{2n}, \qquad |\hat p_0 - p(a)| \le \gamma_{2n} \sum_{k=0}^{n} |c_k|\,|a|^k .

導関数（ $b = 1$ ）。 接成分のステップは $\mathrm{fl}\big(\hat p_{k+1} \cdot 1 + \mathrm{fl}(a\,\hat q_{k+1})\big)$ を計算し、定数 $c_k$ の加算は接成分に厳密な 0 を加えます。導関数には $c_k a^{k-1}$ が $k$ 個含まれ、それぞれは $c_k$ が値の連鎖から接成分の連鎖に移るステップ $j < k$ に対応します。そのような経路に沿って、 $c_k$ は入るときに 1 回、値のステップ $k - 1, \dots, j + 1$ ごとに 2 回、接成分に移るときに 1 回、接成分のステップ $j - 1, \dots, 0$ ごとに 2 回丸められます。すなわち $1 + 2(k - 1 - j) + 1 + 2j = 2k$ 個の因子です。したがって

|\hat q_0 - p'(a)| \le \gamma_{2n} \sum_{k=1}^{n} k\,|c_k|\,|a|^{k-1} ,

となり、これは形式的導関数をホーナー法で評価する場合と同じ上界です。項が打ち消し合わない限り、すなわち $p'$ が $a$ で悪条件でない限り、相対誤差は小さくなります。

コスト

密：双対数の積和ステップが $n + 1$ 回、すなわち T の乗算約 $3n$ 回と加算約 $3n$ 回に加え、持ち上げた $n + 1$ 個の係数の配列 1 つ。疎：項ごとに二分累乗 1 回（ $O(\log e)$ 回の双対数の積）に加え、双対数の積 1 回と加算 1 回。

却下した代替案

記号的に微分してから評価する。 パスが 2 回、余分な多項式が 1 つ必要で、T に対してより強い境界が必要です。多項式そのものが必要なユーザー向けには、luna-poly 自身の DensePolynomial::derivative として残っています。
独自のホーナー法のループ。 定数倍だけ速くなりますが、luna-poly の意味論と重複し、やがて食い違っていきます。
多変数の疎な多項式で変数を推測する。 「第 1 変数」に関する偏導関数を黙って返すと誤用されやすいため、ブリッジは代わりに中断します。

境界

1 変数のみです。多変数多項式の偏導関数や勾配はなく、ContextPolynomial のサポートもありません。
各点での 1 階導関数のみです。導関数の多項式（luna-poly を使ってください）や高階導関数はありません。
チェック付きの版はありません。密な評価と疎な評価は、上で説明した疎な多項式の中断を除いて失敗しません。
密と疎の immut 表現のみです。

D. E. Knuth, The Art of Computer Programming, vol. 2, 3rd ed., section 4.6.4. ↩
N. J. Higham, Accuracy and Stability of Numerical Algorithms, 2nd ed., SIAM, 2002, Lemma 3.1 and chapter 5. ↩