linalg 设计

本页解释 linalg 驱动函数如何从标量对偶数得到梯度和雅可比矩阵、各自的代价，以及为什么这个包被设计成通往 linear-algebra 的一层薄桥接。

设计目标

对使用 linear-algebra 不可变向量编写的函数 $f : T^n \to T$ 和 $f : T^n \to T^m$ 求导，原样复用 Dual[T]，并提供一个结果遵循常规矩阵约定的小型 API。

数学背景

一遍计算得到方向导数

设 $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 ,

因此在对偶数上求值一次即可得到雅可比-向量积 $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 将其作为向量返回。

与反向模式的代价比较

设 $C(f)$ 为在 $T$ 上对 $f$ 求值的代价。一遍对偶计算的代价至多为 $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)$ ，每遍得到一行，以与 $n$ 无关的 $C(f)$ 常数倍代价得到梯度，代价是需要存储计算过程。¹1 这就是”廉价梯度原理”；见 A. Griewank 与 A. Walther 所著 Evaluating Derivatives（第 2 版，SIAM，2008）第 4.6 节。因此当 $n$ 较小或 $n \lesssim m$ 时，前向模式是合适的工具；而对于多变量函数的梯度，它是较慢的那一种。

设计决策

标量切向分量与 $n$ 遍计算

问题。 一个梯度需要 $n$ 个方向导数。

选项。 带有 $n$ 个切向分量向量的对偶类型（一遍计算，每次运算 $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 会以全零切向分量再对 f 求值一次，而不是从某一遍带种子的计算中读取值。根据投影同态，每一遍的值都相同，因此这会多花一次求值，但能让各驱动函数彼此独立，并在 $n = 0$ 时依然正确。jacobian 用同样的零切向分量调用在分配矩阵之前得知 $m$ 。

暂不校验形状

问题。 函数可能越界读取输入，或改变输出长度。

选项。 返回带形状错误的 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$ 矩阵之前要保存 $n$ 个长度为 $m$ 的输出向量。

被否决的方案

用反向模式计算梯度。 当 $n$ 很大时渐近代价更低，但它需要记录计算过程，目前尚未实现。
公开的雅可比-向量积函数。 用 Dual::new(x[j], v[j]) 做一遍计算就已经能得到 $J_f(x) v$ （见 linalg 教程）；专门的函数增益甚微。
可变矩阵。 驱动函数返回 immut 值，与本仓库其余部分的值语义保持一致。

边界

仅支持前向模式：没有反向模式，也没有黑塞矩阵或高阶驱动函数。
没有形状检查，也没有带检查的变体。
仅支持稠密的 immut 向量和矩阵；不支持稀疏或可变容器。
不支持需要 Dual[T] 上的域结构或序的算法。

这就是”廉价梯度原理”；见 A. Griewank 与 A. Walther 所著 Evaluating Derivatives（第 2 版，SIAM，2008）第 4.6 节。 ↩