综合微分法 [synthetic-differential]

  • May 9, 2024 — May 14, 2024
  • kokic

Exegesis 1. 数值微分的误差 [numerical]

$\gdef\spaces#1{~ #1 ~}$ $\gdef\quads#1{\quad #1 \quad}$ $\gdef\eqq{\quads=}$

我们首先观察朴素的数值微分的计算过程. 待微分函数为 $(x+3)^2$, 此处不妨取 $x=-2$, $h=10^{-5}$, 这一位置的导数 $A$ 近似为

$$ A \quads\approx \frac{f(x+h)-f(x)}h \eqq \frac{(1 + 10^{-5})^2 - 1^2}{10^{-5}} \eqq 2.00001 $$

并且如果试图得到与预期的 $A=2$ 更接近的结果, 则须将无穷小近似 $h$ 选取为更靠近 $0$ 的数值, 这对于浮点精度有限的计算机来说极为不利.

而常规符号微分程序视表达式为树结构, 通过遍历树并对符合模式的节点应用求导规则 $(uv)' = u'v+uv'$. 即

$$ \begin{aligned} (x+3)^2 & \spaces\to (x+3)'(x+3) + (x+3)(x+3)' \\ & \spaces\to 1 \cdot (x+3) + (x+3) \cdot 1 \\ & \spaces\to 2x+6 \end{aligned} $$

此时再应用 $x=-2$ 得到 $A=2$. 这个做法在结果上可行, 但处理复杂的表达式时, 复合求导的中间过程本身会产生大量新的结点, 使递归遍历的实际时间达到指数级 $\mathcal{O}(2^N)$.

有没有方法可以既避免递归又减少中间结点呢? 我们先看一个涉及复数域的观察.

Exegesis 2. 复步微分法 [complex-step]

$\gdef\quads#1{\quad #1 \quad}$ $\gdef\spaces#1{~ #1 ~}$

Mathworks 的创始人 Cleve Moler 约 60 年 1 前提出了借用虚数单位 $i$ 的数值微分, 相关文献称为 complex step 微分法 2. 考虑光滑函数 $f(x)$, 其在 $x=a$ 处可表为关于 $X$ 的 Taylor 级数 3 $f(a) + f'(a)(X-a) + \text{etc.}$, 注意

$$ f(a+i h) \spaces= f(a) + i h f'(a) - \frac{h^2}{2!}f''(a) - \frac{ih^3}{3!}f'''(a) + \cdots $$

这个其实就是将 Taylor 级数的每一项完全展开, 在 后续 的计算中也会继续用到. 如此便有 $$ \frac{\partial f}{\partial x} \spaces= \frac{\Im f(x+ih)}{h} + \mathcal{O}(h^2) \quads\implies \frac{\partial f}{\partial x} \quads\approx \frac{\Im f(x+ih)}{h} $$

这个方法最初被设计用于处理数值微分问题, 但稍加思考就能发现, 该过程也适用于符号微分. 与前一个问题所改进的结果的精度不同, 用于 符号微分 时, 所取得的优势是更精简的中间表达式和非递归的计算过程.

1

即 1967 年. 不过这个名字要等到 1998 年, William Squire 和 George Trapp 才正式提出.

2

按方法的流程看, complex step 微分 可以翻译为 “复步微分”.

3

习惯上会混淆 $x$ 与 $X$.

Exegesis 3. 误差修正 [moler-correction]

$\gdef\C{\mathbf{C}}$ $\gdef\quads#1{\quad #1 \quad}$ $\gdef\eqq{\quads=}$

选取同样的 $h=10^{-5}$, 以下计算表明 Moler 的方法 相比 此处 基于有限差分的数值微分方法能够很好的避免精度损失.

$$ \frac{\Im(1 + 10^{-5}i)^{2}}{10^{-5}} \eqq \frac{\Im(\square + 2 \times 10^{- 5}i)}{10^{- 5}} \eqq 2 $$

实际上此处 $h$ 的选取, 只要其非零, 从而带有 $i$ 的项不消失, 就不会影响到最终的结果. 换言之, 这个方法从根本上来说, 是无关于 $h \in \R\smallsetminus\{0\}$ 和 $i$ 的, 只是在数值计算时借由复数算术和复函数值可以免去一些幂零结构的讨论. 一旦意识到这一点, 我们便可以将 复步微分法 的

$$ \frac{\partial f}{\partial x} \quads\approx \frac{\Im f(x+ih)}{h} $$

中的近似 “$\approx$” 修正 为严格等于 “$=$”, 这只需要将 $\{z \in \C : z^2 = -1\}$ 替换为 $D = \{x \in R : x^2 = 0\}$, 问题就来到了如何构造这样的 $R$ 使得 $D$ 中有非零的元素, 这样的结构实际上会 动摇 经典逻辑中对于排中律的看法.

Definition 4. 对偶数 [dual-number]

$\gdef\R{\mathbf{R}}$ $\gdef\C{\mathbf{C}}$

一种通俗的讲法是, 认为对偶数 $R[x]/(x^2)$ 可以作为复数 $\R[x]/(x^2+1)$ 的类比. 如果有读者能轻率地暂时忽略 $R[x]/(x^2)$ 这种结构的存在性, 或者说接受了下面的 $d$ 的定义

$$ \exists ~ d \in R\smallsetminus\{0\}, ~ d^2 = 0 $$

不能接受也很正常, 因为这轻易就会导致某个经典逻辑中的 矛盾. 相比之下满足 $x^2 = -1$ 的 $x$ 对于大众而言要更容易接受的多, 甚至对于三次方程是 必须品. 为此我们需要下面的准备工作.

Axiom 4.1. Kock–Lawvere [kock-lawvere]

$\gdef\spaces#1{~ #1 ~}$ $\gdef\R{\mathbf{R}}$

首先需要强调的是, 这个公理所涉及到的结构的有效性全都依赖于 意象理论, 而相关基础的严格性保证基本上在 1970 年之前就已经由 William Lawvere 完工. 记 $\mathcal{E}$ 是光滑空间之间的光滑映射构成的范畴, 同时假定 $\mathcal{E}$ 还是笛卡尔闭范畴, 也就是 $\mathcal{E}$ 中的箭头对笛卡尔积封闭.

我们可以从 $\mathcal{E}$ 中选出一条几何直线 $R$, 通过指定 $R$ 上两点 $0$ 和 $1$ 之间的距离作为单位长度来确定其他线段的长度. 发挥一些古希腊精神, 线段的移动可以给出 $R$ 上的加法, 尺规作图所构造的 相似三角形 能给出 $R$ 上的乘法 $\gamma = \alpha\beta$.

Methods in Algebra - Volume 1, p. 369

因此 $R$ 带有交换环结构, 并且允许经典数学中的对象实数环 $\R$ 成为 $R$ 的模型, 注意这里我们只能考虑 $\R$ 作为环的部分, 因为 $R$ 中存在着幂零元.

Kock–Lawvere 公理说的是, 对任意映射 $f: D \to R$, 存在唯一的 $a,b \in R$, 使得

$$ f(\epsilon) \spaces= a + b \epsilon, \quad \forall \epsilon \in D $$

将这里的 $a$ 换成 $f(0)$, 并完全使用量词叙述, 则是

$$ \forall ~ f \in R^D, ~ \exists! ~ b ~ \text{s.t.} ~ \forall ~ d \in D \quad (f(d) = f(0) + b d) $$

如果说这个公理的形式还不足以暗示它的目的, 那么下面这个推论就完全能做到了.

Corollary. 所有函数都光滑 [smooths]

$\gdef\quads#1{\quad #1 \quad}$ $\gdef\spaces#1{~ #1 ~}$

对任意函数 $f:R \rightarrow R$, 存在唯一的函数 $f':R \rightarrow R$ 满足

$$ f(x + \varepsilon) \spaces= f(x) + f'(x)\varepsilon,\quad \forall ~ x \in R, ~ \forall ~ \varepsilon \in D $$

这个通常作为 Kock–Lawvere 公理 的推论出现, 不过从 $\mathcal{E}$ 的定义来看, 这才是整套框架真正的目的之一. 即, 为全体函数恢复牛顿时期 “幂零无穷小量” 计算上的直观, 而不导致矛盾.

Axiom is incompatible with the law of excluded middle. Either the one or the other has to leave the scene. In Part I of this book, the law of excluded middle has to leave, being incompatible with the natural synthetic reasoning on smooth geometry to be presented here. In the terms which the logicians use, this means that the logic employed is ‘constructive’ or ‘intuitionistic’. We prefer to think of it just as ‘that reasoning which can be carried out in all sufficiently good cartesian closed categories’.

— Anders Kock, Synthetic Differential Geometry

无论是单纯接受这个 公理 还是认可该公理存在的舞台 $\mathcal{E}$, 其实都会导致完全相同的结果, 那就是 $\mathcal{E}$ 当中一般函数的性质发生了变化, 让所有的函数都变得光滑. 与之相对的, 这样的好性质所要求的代价是, $\mathcal{E}$ 当中不能使用经典逻辑中的选择公理、排中律、反证法等命题.

Exegesis 4.2. 矩阵表示 [matrix-represent]

$\gdef\quads#1{\quad #1 \quad}$ $\gdef\str#1{{\footnotesize #1}}$ $\gdef\C{\mathbf{C}}$ $\gdef\spaces#1{~ #1 ~}$

这样一来, 我们所说的 “将近似修正为严格等于” 就可以精确表示为

$$ \frac{\partial f}{\partial x} \quads= f(x+\varepsilon) ~ \str{中} ~ \varepsilon ~ \str{项的系数} $$

复步微分法 的来源是完全解析的, 但其得到的微分方法却能代数地刻画. 正如我们可以考虑复数的 $z = a+b\sqrt{-1} \in \C$ 矩阵表示

$$ A \spaces= \begin{pmatrix} a & -b \\ b & ~~~a \\ \end{pmatrix} ,\quad \det A \spaces= (a+b\sqrt{-1})(a-b\sqrt{-1}) $$

我们也可以构造 对偶数 $a + b \varepsilon \in R$ 的矩阵表示, 并自动得到微分运算的加法和乘法规则, 这可以被简单地视为 Moler 方法的矩阵版本.

$$ A \spaces= \begin{pmatrix} a & a' \\ 0 & a \\ \end{pmatrix} ,\quad \det A \spaces= (a+b\varepsilon)(a-b\varepsilon) $$

Example 4.3. 导出微分运算 [derived-differential]

$\gdef\Mat{\operatorname{Mat}}$ $\gdef\d{\operatorname{d}}$ $\gdef\Z{\mathbf{Z}}$ $\gdef\spaces#1{~ #1 ~}$

设 $u, v \in R^R$, $n \in \Z$. 记 $A(u, u') = (\begin{smallmatrix} u & u' \\ 0 & u \end{smallmatrix})$. 从矩阵运算中立刻得到

  1. $A(u, u') + A(v, v') \spaces= A(u+v, u'+v')$
  2. $A(u, u') - A(v, v') \spaces= A(u-v, u'-v')$
  3. $A(u, u') \cdot A(v, v') \spaces= A(uv, u'v + uv')$
  4. $A(u, u') \cdot (A(v, v'))^{-1} \spaces= A(uv^{-1}, (u'v - uv')v^{-2})$
  5. $(A(u, u'))^n = A(u^n, nu^{n-1}u')$

对于 $e^{A(u,u')}$, 注意 $A(u, u') = A(u, 0) + A(0, u')$, 这里 $(A(0, u))^2 = 0$. 同样可以计算得到 $e^{A(u,u')} = A(e^u, 0) \cdot A(1, u') = A(e^u, e^u u')$. 而对于 $\log A(u,u')$ 我们直接展开为 $\sum_{k \ge 1}\frac{(-1)^{k+1}}{k} \cdot (A-1)^k$ 即可验证:

$$ \begin{pmatrix} \sum_{k \ge 1}\frac{(-1)^{k+1}}{k} \cdot (u-1)^k & \sum_{k \ge 0} (-1)^k (u-1)^k \cdot u' \\ 0 & \sum_{k \ge 1}\frac{(-1)^{k+1}}{k} \cdot (u-1)^k \end{pmatrix} \spaces= \begin{pmatrix} \log u & u'u^{-1} \\ 0 & \log u \end{pmatrix} $$

Exegesis 4.4. 链式法则 [chain-rule]

$\gdef\Mat{\operatorname{Mat}}$ $\gdef\spaces#1{~ #1 ~}$ $\gdef\d{\operatorname{d}}$ $\gdef\eqdef{\overset{\scriptscriptstyle\text{def}}{=}}$

与 $\Mat_{2 \times 2}(R^R)$ 上的 运算 不同, 复合运算是只位于 $R^R$ 上的, 或者说 $\Mat_{2 \times 2}(R^R)$ 上并没有周知的复合运算. 为了自然地给出 $A: J^1(R, R) \to \Mat_{2 \times 2}(R^R)$ 上的复合, 这里 $J^1$ 是一阶 射流丛. 记 $\square(x)$ 为 $\square_x $, $A(u)$ 为 $A(u,u')$, 考虑显式应用的 $A(u)$, 即

$$ A(u)_x \spaces\eqdef \begin{pmatrix} u_x & u'_x \\ 0 & u_x \end{pmatrix} \spaces\in \Mat_{2 \times 2}(R), \quad x \in R $$

对于 $u: Y \to Z$, $v: X \to Y$. $A$ 上的复合 $A(u) \circ A(v)$ 可定义为 $A(u)_{v_x}$, $(v_x, v'_x) \in A(v)_x$. 我们来验证 $A(u) \circ A(v) = A(u \circ v)$.

$$ \begin{aligned} A(u)_{v_x} &\spaces= \Big(u \cdot \frac{\partial A}{\partial u} + \frac{\d u}{\d x} \cdot \frac{\partial A}{\partial u'}\Big)\Big|_{v_x} \\ &\spaces= u(v_x) \cdot \frac{\partial A}{\partial u} + \frac{\d u(v_x)}{\d v_x} \frac{\d v_x}{\d x} \cdot \frac{\partial A}{\partial u'} \\ &\spaces= A(u(v(x)), u'(v(x))v'(x)) \\ \end{aligned} $$