linalg 設計

このページでは、linalg のドライバーがスカラーの双対数から勾配とヤコビ行列をどのように得るか、それぞれのコスト、そしてパッケージが linear-algebra への薄いブリッジとして構成されている理由を説明します。

設計目標

linear-algebra のイミュータブルなベクトルで書かれた関数 $f : T^n \to T$ と $f : T^n \to T^m$ を、Dual[T] をそのまま再利用し、結果が通常の行列の慣習に従う小さな API で微分します。

数学的背景

1 回のパスによる方向微分

$f : \mathbb R^n \to \mathbb R^m$ が $x$ で微分可能であるとし、 $v \in \mathbb R^n$ とします。各入力を $x_j + v_j\varepsilon$ としてシードします。dual の設計により、各出力は次のようになります

f_i(x + v\varepsilon) = f_i(x) + \sum_{j} \frac{\partial f_i}{\partial x_j}(x)\,v_j\,\varepsilon = f_i(x) + \big(J_f(x)\,v\big)_i\,\varepsilon ,

したがって双対数上での 1 回の評価でヤコビ行列ベクトル積 $J_f(x) v$ が得られます。これは多変数の連鎖律です。 $g(t) = f(x + tv)$ に対して $g'(0) = J_f(x) v$ となります。

単位ベクトルのシードによる列

単位ベクトル $e_j$ をシードすると、ヤコビ行列の第 $j$ 列 $J_f(x) e_j = \partial f / \partial x_j$ が返ります。ドライバーはこのようなパスを $n$ 回行います：

J_f(x) = \big[\, J_f(x) e_0 \;\big|\; J_f(x) e_1 \;\big|\; \cdots \;\big|\; J_f(x) e_{n-1} \,\big] .

$m = 1$ のとき、唯一の行は $\nabla f(x)^{\mathsf T}$ であり、gradient はそれをベクトルとして返します。

リバースモードとのコスト比較

$T$ 上で $f$ を評価するコストを $C(f)$ とします。双対数のパス 1 回のコストは高々 $C(f)$ の小さな定数倍です（dual の設計を参照）。したがって

C(\texttt{gradient}) \approx n \cdot c \cdot C(f), \qquad C(\texttt{jacobian}) \approx (n + 1) \cdot c \cdot C(f), \qquad c \approx 3 .

リバースモードはベクトル・ヤコビ行列積 $u^{\mathsf T} J_f(x)$ をパスごとに 1 行ずつ計算し、計算を記録するという代償と引き換えに、 $n$ によらない $C(f)$ の定数倍で勾配を得ます。¹1 これは「安価な勾配の原理」です。A. Griewank and A. Walther, Evaluating Derivatives, 2nd ed., SIAM, 2008, section 4.6 を参照してください。したがってフォワードモードは $n$ が小さいか $n \lesssim m$ の場合に適した手法であり、多変数関数の勾配には遅いほうの手法です。

設計上の判断

スカラーの接成分と $n$ 回のパス

問題。 勾配には $n$ 個の方向微分が必要です。

選択肢。 $n$ 個の接成分のベクトルを持つ双対数型（1 回のパス、演算ごとに $O(n)$ の処理）、またはスカラーの Dual[T] による $n$ 回のパス。

選択。 $n$ 回のパスです。算術の総量は同じオーダー $O(n \cdot C(f))$ であり、スカラー型では演算ごとのメモリ確保も新しい数値型も不要です。ベクトルの接成分を持つ型は今後の課題です。

出力×入力の慣習

行列は $m \times n$ で、成分は $(i, j) = \partial f_i / \partial x_j$ です。これは連鎖律 $J_{f \circ g} = J_f\,J_g$ の慣習であり、結果を通常の行列積で合成できます。第 $j$ 列はパス $j$ に対応します。実装は $n$ 個の出力ベクトルを集めてから、Matrix::make(m, n, (i, j) => columns[j][i].tangent()) で行列を構築します。

値と $m$ のための独立したパス

value_and_gradient と value_and_jacobian は、シードしたパスのいずれかから値を読み取るのではなく、すべての接成分を 0 にして f をもう 1 回評価します。射影準同型によりどのパスも同じ値を持つので、これは評価 1 回分のコストがかかりますが、ドライバーを互いに独立に保ち、 $n = 0$ でも正しく動作させます。jacobian も、行列を確保する前に $m$ を知るために同じ接成分 0 の呼び出しを使います。

形状の検証はまだない

問題。 関数が入力の範囲外を読んだり、出力の長さを変えたりすることがあります。

選択肢。 形状エラーを持つ Result を返す、中断する、事前条件として文書化する。

選択。 文書化された事前条件です。linear-algebra はまだこのリポジトリと形状エラーの型を共有しておらず、独自のエラー型は後で置き換えられることになります。ソースにはチェック付きの版のための TODO があります。

`Dual[T]` 上の環演算のみ

ドライバーが必要とするのはシードのための One + Zero だけです。f 内部のベクトル演算と行列演算は Dual[T] の Add、Mul、Neg インスタンスを使うので、環構造だけを必要とする linear-algebra の演算はすべて双対数のベクトル上で動作します。スカラーに Field、Inverse や順序を必要とするアルゴリズム（たとえばピボット選択）は Dual[T] でインスタンス化できません。これは意図的なものです。dual の設計を参照してください。

正しさと不変条件

事前条件を満たすプログラムについて、gradient(f, x)[j] $= \partial f / \partial x_j (x)$ かつ jacobian(f, x)[i][j] $= \partial f_i / \partial x_j (x)$ であり、各成分には dual の設計の丸めの上界が適用されます。
value_and_gradient(f, x) $= (f(x), \nabla f(x))$ かつ value_and_jacobian(f, x) $= (f(x), J_f(x))$ です。値は f が T 上で計算するものとちょうど一致します。
評価回数：gradient は $n$ 、value_and_gradient は $n + 1$ 、jacobian は $n + 1$ 、value_and_jacobian は $n + 2$ 。
メモリ：jacobian は $m \times n$ 行列を構築する前に、長さ $m$ の出力ベクトルを $n$ 個保持します。

却下した代替案

勾配のためのリバースモード。 $n$ が大きい場合は漸近的に安価ですが、計算の記録が必要であり、実装されていません。
公開のヤコビ行列ベクトル積。 Dual::new(x[j], v[j]) による 1 回のパスで既に $J_f(x) v$ が得られるので（linalg チュートリアルを参照）、専用の関数を追加してもほとんど得るものはありません。
ミュータブルな行列。 ドライバーは immut の値を返し、これはリポジトリの他の部分の値の意味論と一致します。

境界

フォワードモードのみです。リバースモードも、ヘッセ行列や高階のドライバーもありません。
形状チェックもチェック付きの版もありません。
密な immut のベクトルと行列のみです。疎なコンテナやミュータブルなコンテナはありません。
Dual[T] 上に体や順序を必要とするアルゴリズムはありません。

これは「安価な勾配の原理」です。A. Griewank and A. Walther, Evaluating Derivatives, 2nd ed., SIAM, 2008, section 4.6 を参照してください。 ↩