5.8 多元微分学与微分形式初步

前四节讨论的函数都只有一个自变量。最优化、机器学习与统计中的函数几乎都是多变量的：目标函数依赖成千上万个参数，损失函数是参数向量的函数。把微分学推广到多维，需要处理新的问题：多变量的函数没有单一的「变化率」，每个方向都有变化率；局部形状由一个矩阵（Hessian）而不是一个数描述；方程组定义的隐式关系需要用矩阵的逆来表达。

本节的组织沿四条线索展开。微分线索定义偏导数、方向导数、全微分与梯度，其中可微是一个比「偏导数存在」更强的条件，实验一用中心差分核对 $f(x,y)=x^{2}y+\sin (xy)$ 在 $(1.2,0.7)$ 处的梯度（偏差 $3.56\times {10}^{-10}$）并验证「梯度方向是方向导数的最大值方向」。二阶线索引入 Hessian 矩阵与二阶 Taylor 公式：误差项从一元情形的 $O(h^3)$ 变成 $O(\parallel h{\parallel }^3)$（实验二中误差除以 $\parallel h{\parallel }^3$ 稳定在 $0.52$），Hessian 的特征值符号给出极值判别（正定是极小、不定是鞍点、半定需要更高阶信息）。隐函数线索给出反函数定理与隐函数定理，它们把「方程组有没有解」这个存在性问题化为「Jacobian 矩阵是否可逆」这个线性代数问题，是 Lagrange 乘子法（第三章 3.6 节）的理论来源。微分形式线索把梯度、旋度、散度统一成外微分算子，Green、Gauss、Stokes 三大公式在下一节统一为一条公式。这条线索对计算不直接有用，但提供了理解向量分析的结构性视角。

5.8.1 多元函数的极限与连续

多元函数的极限沿用数列的语言，但「接近」的含义变了：点 $x$ 接近 $x_0$ 指欧氏距离 $\parallel x-x_0\parallel \to 0$，而空间中的接近路径有无穷多条。定义：对任意 $\epsilon >0$ 存在 $\delta >0$，使得 $0<\parallel x-x_0\parallel <\delta$ 时 $|f(x)-f(x_0)|<\epsilon$。

路径的多样性带来一元情形没有的现象：沿所有路径的极限都存在且相等才可能有极限，而检查所有路径是不可能的。判定极限不存在的方法比较直接：找两条路径给出不同的极限。例：$f(x,y)=\frac{xy}{x^2+y^2}$ 在原点。沿 $y=kx$ 接近时 $f=\frac{k}{1+k^2}$，取值随 $k$ 变化（$k=1$ 时是 $\frac{1}{2}$，$k=0$ 时是 $0$），极限不存在。例：$f(x,y)=\frac{x^{2}y}{x^4+y^2}$ 沿 $y=kx$ 得 $\frac{kx}{1+k^2{x}^2}\to 0$，沿抛物线 $y=x^2$ 得 $\frac{1}{2}$，极限仍不存在，说明只检查直线路径不够。

连续的定义与一元相同：$f$ 在 $x_0$ 连续指 $\underset{x\to x_0}{lim}f(x)=f(x_0)$。多元连续函数在有界闭集（紧集）上仍然取到最大值与最小值（5.3 节结论的推广），这条性质在最优化中保证「最小值存在」，是分析算法收敛性的前提。

5.8.2 全微分、偏导数与梯度

偏导数与方向导数

偏导数是把其他变量固定后的一元导数：

$$\frac{\partial f}{\partial x_i}(x_0)=\underset{h\to 0}{lim}\frac{f(x_0+h{e}_i)-f(x_0)}{h}$$

其中 $e_i$ 是第 $i$ 个坐标方向的单位向量。偏导数度量沿坐标轴的变化率，几何含义是曲面被坐标平面截出的曲线在该点的切线斜率。计算上，求 $\frac{\partial f}{\partial x_1}$ 时把其余变量当作常数，一元求导法则全部适用。

方向导数把方向从坐标轴推广到任意单位向量 $v$：

$$D_{v}f(x_0)=\underset{t\to 0}{lim}\frac{f(x_0+tv)-f(x_0)}{t}$$

偏导数是方向导数在 $v=e_i$ 时的特例。方向导数对每个方向独立定义，方向有无穷多个，因此「所有方向导数存在」是一条很强的要求，而「所有偏导数存在」是很弱的要求。

可微与全微分

可微的定义是一元情形（5.4 节）的直接推广：存在向量 $A$ 使得

$$f(x_0+h)=f(x_0)+A\cdot h+o(\parallel h\parallel ),h\to 0$$

右端的线性部分 $A\cdot h$ 称全微分。可微蕴含所有方向导数存在，且 $D_{v}f(x_0)=A\cdot v$；取 $v=e_i$ 得 $A_i=\frac{\partial f}{\partial x_i}(x_0)$。因此可微时

$$D_{v}f(x_0)=\sum _{i=1}^n\frac{\partial f}{\partial x_i}(x_0)v_i=\mathrm{\nabla }f(x_0)\cdot v$$

其中 $\mathrm{\nabla }f=(\frac{\partial f}{\partial x_1},\dots ,\frac{\partial f}{\partial x_n})$ 是梯度。这条公式说明方向导数关于方向是线性的，只要知道 $n$ 个偏导数就能算出所有方向的导数，这是梯度作为完整一阶信息的含义。

偏导数存在不蕴含可微。反例是

$$f(x,y)=\{\begin{array}{ll}{\displaystyle \frac{xy}{\sqrt{x^2+y^2}}},& (x,y)\ne (0,0)\\ 0,& (x,y)=(0,0)\end{array}$$

两个偏导数在原点都存在且等于零（沿坐标轴函数恒为零），但沿 $y=x$ 时 $f=\frac{|x|}{\sqrt{2}}$，$\frac{f(x,x)}{\parallel (x,x)\parallel }\to \frac{1}{2}$，不满足可微定义中的 $o(\parallel h\parallel )$ 条件。充分条件：若偏导数在 $x_0$ 的邻域内存在且在 $x_0$ 处连续，则 $f$ 在 $x_0$ 可微。证明用一元中值定理把差 $f(x_0+h)-f(x_0)$ 沿坐标方向逐项拆开，每一项等于对应偏导数在某中间点的值乘以该分量，再利用偏导数的连续性把中间点的值换成 $x_0$ 处的值。

可微的两个推论在计算中常用：可微蕴含连续；可微函数沿水平集移动时函数值的一阶变化为零，因此梯度与水平集（等高线、等值面）正交。

梯度与最速上升方向

由 $D_{v}f=\mathrm{\nabla }f\cdot v=\parallel \mathrm{\nabla }f\parallel \cos \theta$（$\theta$ 是 $v$ 与梯度的夹角），方向导数在 $v$ 与梯度同向时取最大值 $\parallel \mathrm{\nabla }f\parallel$，在反向时取最小值 $-\parallel \mathrm{\nabla }f\parallel$，在与梯度垂直时为零。这三条结论给出梯度的完整几何含义：

1. 梯度的方向是函数上升最快的方向（最速上升方向），负梯度方向是下降最快方向。
2. 梯度垂直于水平集，因为沿水平集移动时方向导数为零。
3. 梯度的范数是最大方向导数的值，度量该点最陡的程度。

梯度下降法（第三章 3.3 节）的全部依据就是第 1 条：沿着 $-\mathrm{\nabla }f$ 走使函数值下降，步长的选择则依赖第 3 条与二阶信息。实验一的数据给出这三条结论的具体数值：梯度 $(2.147223978089,2.240955391010)$、范数 $3.103619157786$，方向导数在角度 $0^{\circ }$（沿 $x$ 轴）时为 $2.1472$，${90}^{\circ }$ 时为 $2.2410$，${120}^{\circ }$ 时降到 $0.8671$，${180}^{\circ }$ 时取 $-2.1472$（反向下降最快），而沿梯度方向的值恰好等于范数 $3.1036$。

5.8.3 链式法则与 Jacobian 矩阵

链式法则

设 $f:{\mathbb{R}}^n\to \mathbb{R}$ 可微，$g:{\mathbb{R}}^m\to {\mathbb{R}}^n$ 可微，复合函数 $h=f\circ g$ 的梯度是

$$\mathrm{\nabla }h(x)=J_g(x{)}^{\mathsf{T}}\mathrm{\nabla }f(g(x))$$

其中 $J_g$ 是 $g$ 的 Jacobian 矩阵，$(J_g{)}_{ij}=\frac{\partial g_i}{\partial x_j}$。分量形式是

$$\frac{\partial h}{\partial x_j}=\sum _{i=1}^n\frac{\partial f}{\partial y_i}\frac{\partial g_i}{\partial x_j}$$

这条公式是一元链式法则在多变量下的形式，求和出现在中间变量上。神经网络的反向传播就是这条公式的递归应用：把网络看作多层复合函数，从输出层往输入层逐层用链式法则传递梯度，每层的 Jacobian 由该层的运算决定。

变换与 Jacobian 行列式

坐标变换 $x=\mathrm{\Phi }(u)$ 的局部行为由 Jacobian 矩阵描述：小体积元的变换是

$$d{x}_{1}d{x}_2\cdots d{x}_n=|det{J}_{\mathrm{\Phi }}|d{u}_{1}d{u}_2\cdots d{u}_n$$

$|det{J}_{\mathrm{\Phi }}|$ 是体积膨胀率，这个因子在重积分的换元公式（5.9 节）中起核心作用。例：极坐标 $x=r\cos \theta$、$y=r\sin \theta$ 的 Jacobian 行列式是 $r$，因此 $dxdy=rdrd\theta$；球坐标的 Jacobian 行列式是 $r^2\sin \phi$，对应体积元 $r^2\sin \phi drd\phi d\theta$。5.8.7 节实验一用的数值梯度核对方式（中心差分）也是链式法则的离散版本。

5.8.4 高阶导数、Hessian 与 Taylor 公式

混合偏导数与 Hessian 矩阵

二阶偏导数记作 $\frac{{\partial }^{2}f}{\partial x_i\partial x_j}$。Clairaut-Schwarz 定理：若两个混合偏导数在 $x_0$ 的邻域内存在且连续，则它们相等：

$$\frac{{\partial }^{2}f}{\partial x_i\partial x_j}(x_0)=\frac{{\partial }^{2}f}{\partial x_j\partial x_i}(x_0)$$

因此二阶偏导数构成一个对称矩阵

$$H={\left(\frac{{\partial }^{2}f}{\partial x_i\partial x_j}\right)}_{n\times n}$$

称 Hessian 矩阵。对称性使 $H$ 的特征值都是实数，可以按符号分类，这正是一元情形中 $f^{″}$ 的符号在多维下的推广。

二阶 Taylor 公式与极值判别

$f$ 在 $x_0$ 附近二阶可微时

$$f(x_0+h)=f(x_0)+\mathrm{\nabla }f(x_0)\cdot h+\frac{1}{2}{h}^{\mathsf{T}}H(x_0)h+o(\parallel h{\parallel }^2)$$

极值判别（$x_0$ 是临界点，即 $\mathrm{\nabla }f(x_0)=0$）：

1. $H$ 正定（所有特征值为正）时 $x_0$ 是局部极小点。
2. $H$ 负定（所有特征值为负）时 $x_0$ 是局部极大点。
3. $H$ 不定（特征值有正有负）时 $x_0$ 是鞍点。
4. $H$ 半定（有零特征值）时二阶信息不足，需要更高阶项。

实验二的三行数据给出三种情形：$x^2+3xy+5y^2$ 的特征值是 $0.5$ 与 $5.5$（正定），$x^2-y^2$ 的是 $-1$ 与 $1$（不定，鞍点），$x^2-2xy+y^2$ 的是 $0$ 与 $2$（半定，二次型退化为 $(x-y{)}^2$，沿直线 $x=y$ 方向上二阶项为零，需要四阶项判断：沿该直线函数恒为零，整条直线都是极小点，属于退化极小）。

有一点需要区分：Hessian 的特征值符号只在临界点给出极值判据。实验二在 $(1.2,0.7)$ 处算出的 Hessian 特征值是 $-2.67818145$ 与 $2.64102023$（不定），但该点不是临界点（梯度范数 $3.1036$），特征值的含义是该点附近沿某些方向的曲率为负，这类信息在信赖域方法与鞍点逃离（第三章 3.8 节）中有直接用途：负曲率方向指向函数值可以继续下降的方向。

二阶 Taylor 的误差阶

二阶 Taylor 的误差是 $O(\parallel h{\parallel }^3)$，与一元情形相同。实验二的五组数据验证这一点：把 $|h|$ 从 ${10}^{-1}$ 减到 $5\times {10}^{-3}$，误差从 $5.34\times {10}^{-4}$ 降到 $6.50\times {10}^{-8}$，而误差除以 $\parallel h{\parallel }^3$ 稳定在 $0.52$ 附近（$0.534$、$0.527$、$0.522$、$0.520$、$0.520$），说明误差确实按三次方衰减，系数趋于一个常数。这条数据的用法与一元情形相同：用「误差除以步长的高次幂」是否稳定来判断收敛阶。

5.8.5 反函数定理与隐函数定理

反函数定理

设 $F:{\mathbb{R}}^n\to {\mathbb{R}}^n$ 在 $x_0$ 处连续可微，且 Jacobian 矩阵 $J_F(x_0)$ 可逆，则 $F$ 在 $x_0$ 附近存在连续可微的局部反函数 $F^{-1}$，且

$$J_{F^{-1}}(y_0)={[J_F(x_0)]}^{-1},y_0=F(x_0)$$

定理的内容可以概括为「可逆的线性化给出局部的可逆性」。证明的思路是用压缩映射原理（5.1 节）把求解 $F(x)=y$ 化为不动点问题：由 $J_F(x_0)$ 可逆，把方程改写成

$$x=x_0+J_F(x_0{)}^{-1}[y-F(x_0)-(F(x)-F(x_0)-J_F(x_0)(x-x_0))]$$

右端关于 $x$ 是压缩映射（可微性保证方括号内的项是 $o(\parallel x-x_0\parallel )$，系数可以做到小于 $1$），由压缩映射原理得到不动点存在唯一。这条路线把非线性方程组的局部可解性归约到线性代数的可逆性，是分析隐式关系的基本手法。

隐函数定理

隐函数定理：设 $F:{\mathbb{R}}^n\times {\mathbb{R}}^m\to {\mathbb{R}}^m$ 连续可微，$F(x_0,y_0)=0$（$x\in {\mathbb{R}}^n$ 是自变量，$y\in {\mathbb{R}}^m$ 是待解变量），且关于 $y$ 的 Jacobian $\frac{\partial F}{\partial y}(x_0,y_0)$ 可逆，则存在 $x_0$ 的邻域与唯一的连续可微函数 $y=\phi (x)$ 使得 $F(x,\phi (x))=0$，且

$$\frac{\partial \phi }{\partial x}(x_0)=-{[\frac{\partial F}{\partial y}(x_0,y_0)]}^{-1}\frac{\partial F}{\partial x}(x_0,y_0)$$

公式的来源是把 $F(x,\phi (x))=0$ 对 $x$ 求导（链式法则）：

$$\frac{\partial F}{\partial x}+\frac{\partial F}{\partial y}\frac{\partial \phi }{\partial x}=0⟹\frac{\partial \phi }{\partial x}=-{\left(\frac{\partial F}{\partial y}\right)}^{-1}\frac{\partial F}{\partial x}$$

它把隐函数的导数化成了偏导数矩阵的逆乘偏导数矩阵，全部是线性代数运算。实验三用三个例子核对这条公式。

一元隐函数：$F(x,y)=x^2+y^2-1=0$（上半圆）。$\frac{\partial F}{\partial y}=2y$，$\frac{\partial F}{\partial x}=2x$，故

$$\frac{dy}{dx}=-\frac{2x}{2y}=-\frac{x}{y}$$

在 $x=0.6$、$y=0.8$ 处是 $-0.75$，数值差分给出 $-0.75000000$，偏差 $4.93\times {10}^{-11}$。另一条信息藏在公式里：$\frac{\partial F}{\partial y}=2y$ 在 $y=0$ 处为零，公式失效，这一点对应 $x=\pm 1$ 处（圆的左右端点）$y$ 不再是 $x$ 的函数（同一个 $x$ 附近有两个 $y$ 分支），与定理的条件一致。

二元隐函数：$F(x,y,z)=x^2+y^2+z^2-1=0$，$\frac{\partial F}{\partial z}=2z$，因此

$$\frac{\partial z}{\partial x}=-\frac{x}{z},\frac{\partial z}{\partial y}=-\frac{y}{z}$$

在 $(x,y)=(0.3,0.4)$ 处 $z=0.8660254038$，解析值 $(-0.34641016,-0.46188022)$ 与数值差分逐位一致（偏差 $6.9\times {10}^{-11}$ 与 $5.1\times {10}^{-13}$）。

方程组：两个方程决定两个未知数

$$\{\begin{array}{l}{u}^2+v^2-1=0\\ u-v-x=0\end{array}$$

关于 $(u,v)$ 的 Jacobian 是 $J=\left(\begin{array}{cc}2u& 2v\\ 1& -1\end{array}\right)$，$\frac{\partial F}{\partial x}=\left(\begin{array}{c}0\\ -1\end{array}\right)$，由隐函数定理

$$\left(\begin{array}{c}du/dx\\ dv/dx\end{array}\right)=-J^{-1}\left(\begin{array}{c}0\\ -1\end{array}\right)$$

在 $x=0$ 处 $u=v=\frac{1}{\sqrt{2}}=0.7071067812$，$J=\left(\begin{array}{cc}\sqrt{2}& \sqrt{2}\\ 1& -1\end{array}\right)$，$detJ=-2\sqrt{2}$，解得 $\frac{du}{dx}=0.5$、$\frac{dv}{dx}=-0.5$，与解析解 $u=\frac{x+\sqrt{2-x^2}}{2}$ 的导数一致。数据表还给出 $x=0.5$ 处的值 $(\frac{du}{dx},\frac{dv}{dx})=(0.31101776,-0.68898224)$。

与 Lagrange 乘子的连接

隐函数定理是 Lagrange 乘子法的证明工具。在约束 $g(x)=0$ 下求 $f$ 的极值，若 $\mathrm{\nabla }g(x_0)\ne 0$，则由隐函数定理可以把约束曲面局部参数化（把 $x_n$ 表成其余 $n-1$ 个变量的函数），代入 $f$ 后化为无约束问题，求导得 $\mathrm{\nabla }f=\lambda \mathrm{\nabla }g$。第三章 3.6 节的 KKT 条件、对偶理论与约束优化的收敛性分析都建立在这条定理之上。约束资格条件（如线性无关约束资格）的实质就是保证相关的 Jacobian 行向量线性无关，从而隐函数定理可用。

5.8.6 微分形式初步

外微分与微分形式

微分形式是可以在曲面上积分、并且对坐标变换有确定规律的对象。0 形式是函数 $f$；1 形式写成 $\omega =\sum _i{a}_i(x)d{x}_i$（$d{x}_i$ 是坐标微分）；2 形式写成 $\sum _{i<j}{a}_{ij}(x)d{x}_i\wedge d{x}_j$，其中 $\wedge$ 是外积，满足反交换律

$$d{x}_i\wedge d{x}_j=-d{x}_j\wedge d{x}_i,d{x}_i\wedge d{x}_i=0$$

反交换律是微分形式的核心代数性质，它自动处理了定向与符号。外微分算子 $d$ 把 $k$ 形式映为 $k+1$ 形式，规则是

$$d(\sum a_{I}d{x}_I)=\sum _i\sum _I\frac{\partial a_I}{\partial x_i}d{x}_i\wedge d{x}_I$$

三条性质定义了 $d$ 的行为：线性、Leibniz 律 $d(\omega \wedge \eta )=d\omega \wedge \eta +(-1{)}^k\omega \wedge d\eta$（$\omega$ 是 $k$ 形式）、以及

$$d(d\omega )=0$$

最后一条由混合偏导数的对称性与外积的反交换性相消得到，是最重要的结构性结论。

梯度、旋度、散度的统一

在 ${\mathbb{R}}^3$ 中，外微分把三个经典算子统一起来。对函数 $f$（0 形式），$df=\frac{\partial f}{\partial x}dx+\frac{\partial f}{\partial y}dy+\frac{\partial f}{\partial z}dz$ 的系数是梯度。对 1 形式 $\omega =Pdx+Qdy+Rdz$，

$$d\omega =(\frac{\partial R}{\partial y}-\frac{\partial Q}{\partial z})dy\wedge dz+(\frac{\partial P}{\partial z}-\frac{\partial R}{\partial x})dz\wedge dx+(\frac{\partial Q}{\partial x}-\frac{\partial P}{\partial y})dx\wedge dy$$

三个系数是旋度的分量。对 2 形式 $\eta =Pdy\wedge dz+Qdz\wedge dx+Rdx\wedge dy$，

$$d\eta =(\frac{\partial P}{\partial x}+\frac{\partial Q}{\partial y}+\frac{\partial R}{\partial z})dx\wedge dy\wedge dz$$

系数是散度。三个算子由同一个算子 $d$ 在不同阶的形式上统一起来，各自独立的定义只是坐标表达的选择。

$d(d\omega )=0$ 在 ${\mathbb{R}}^3$ 中给出两条熟悉的关系：梯度场无旋（$\mathrm{\nabla }\times (\mathrm{\nabla }f)=0$），旋度场无散（$\mathrm{\nabla }\cdot (\mathrm{\nabla }\times F)=0$）。

闭合形式与恰当形式

若 $d\omega =0$，称 $\omega$ 闭合；若存在 $\eta$ 使 $\omega =d\eta$，称 $\omega$ 恰当。由 $d(d\eta )=0$，恰当蕴含闭合。Poincaré 引理给出反向的结论（局部）：在单连通区域内，闭合形式一定恰当。

在 ${\mathbb{R}}^3$ 中，这条引理对应两个熟悉的结论：无旋场（旋度为零）局部是某个标量函数的梯度（存在势函数），无散场局部是某个向量场的旋度（存在向量势）。全局与局部的差别由一个拓扑量描述，典型例子是

$$\omega =\frac{-ydx+xdy}{x^2+y^2}$$

在 ${\mathbb{R}}^2\setminus \{0\}$ 上闭合（直接计算得 $d\omega =0$），但不恰当：沿单位圆周积分得 $2\pi \ne 0$，而若 $\omega =d\eta$，由 Stokes 公式（下一节）该积分应当为零。这个例子说明 Poincaré 引理的「单连通」条件不能去掉，${\mathbb{R}}^2\setminus \{0\}$ 上的这个形式测度的正是绕原点的圈数。同一现象在复分析中表现为 $\frac{1}{z}$ 的积分，在物理中表现为磁通量与相位的关系。

与三大积分公式的关系

广义 Stokes 公式把微积分基本定理、Green 公式、Gauss 公式、Stokes 公式统一为一条：

$${\int }_{\partial M}\omega ={\int }_{M}d\omega$$

右端是 $k$ 形式在 $k$ 维区域上的积分，左端是它在区域边界（$k-1$ 维）上的积分。$k=1$、区域为区间时得到 ${\int }_a^{b}df=f(b)-f(a)$（微积分基本定理）；$k=2$、区域为平面区域时得到 Green 公式；$k=3$、区域为曲面时得到 Stokes 公式；区域为三维体时得到 Gauss 公式。下一节会给出这些公式的具体形式与数值验证。

这条线索的用途

微分形式的语言在日常计算中不是必需的（大多数工程计算用梯度、旋度、散度的坐标公式就够了），它的价值有两处。一是结构性：把四五个看似无关的积分公式归为一条，把三个微分算子归为一个，理解成本大幅下降。二是几何与拓扑：闭合但不恰当的形式度量空间的洞，这条思路在电磁学（磁通量量子化）、规范场论（曲率形式与陈类）与机器学习的信息几何（第四章 4.3 节的 Fisher 信息矩阵作为度量）中都出现。

5.8.7 数值实验

实验一：梯度、方向导数与最速上升方向

核对解析梯度与中心差分，并验证方向导数关于方向的线性与极值性质。

import math

def f(x, y):
    return x * x * y + math.sin(x * y)

def grad_f(x, y):
    return (2 * x * y + y * math.cos(x * y), x * x + x * math.cos(x * y))

x0, y0 = 1.2, 0.7
ga = grad_f(x0, y0)
h = 1e-6
gn = ((f(x0 + h, y0) - f(x0 - h, y0)) / (2 * h),
      (f(x0, y0 + h) - f(x0, y0 - h)) / (2 * h))
print(f"f(x, y) = x^2*y + sin(x*y)，在 ({x0}, {y0}) 处")
print(f"  解析梯度     = ({ga[0]:.12f}, {ga[1]:.12f})")
print(f"  中心差分数值 = ({gn[0]:.12f}, {gn[1]:.12f})")
print(f"  最大偏差     = {max(abs(ga[0] - gn[0]), abs(ga[1] - gn[1])):.3e}")

nrm = math.hypot(*ga)
print(f"梯度范数 |grad f| = {nrm:.12f}")
print(f"{'角度(度)':>10} {'方向向量':>22} {'解析方向导数':>16} {'数值方向导数':>16}")
for deg in [0, 30, 60, 90, 120, 180, 270]:
    th = math.radians(deg)
    vx, vy = math.cos(th), math.sin(th)
    dan = ga[0] * vx + ga[1] * vy
    t = 1e-6
    dnu = (f(x0 + t * vx, y0 + t * vy) - f(x0 - t * vx, y0 - t * vy)) / (2 * t)
    print(f"{deg:>10} ({vx:>8.5f}, {vy:>8.5f}) {dan:>16.10f} {dnu:>16.10f}")

输出：

f(x, y) = x^2*y + sin(x*y)，在 (1.2, 0.7) 处
  解析梯度     = (2.147223978089, 2.240955391010)
  中心差分数值 = (2.147223977733, 2.240955391031)
  最大偏差     = 3.558e-10
梯度范数 |grad f| = 3.103619157786
     角度(度)                   方向向量           解析方向导数           数值方向导数
         0 ( 1.00000,  0.00000)     2.1472239781     2.1472239777
        30 ( 0.86603,  0.50000)     2.9800282081     2.9800282082
        60 ( 0.50000,  0.86603)     3.0143362864     3.0143362867
        90 ( 0.00000,  1.00000)     2.2409553910     2.2409553910
       120 (-0.50000,  0.86603)     0.8671123083     0.8671123082
       180 (-1.00000,  0.00000)    -2.1472239781    -2.1472239777
       270 (-0.00000, -1.00000)    -2.2409553910    -2.2409553910

前三行给出梯度核对的结果：解析值 $(2.147223978089,2.240955391010)$ 与中心差分 $(2.147223977733,2.240955391031)$ 的最大偏差 $3.56\times {10}^{-10}$。这个量级与 5.4 节中心差分的最优精度（约 ${10}^{-11}$）相比略差，因为这里取固定步长 ${10}^{-6}$ 而该函数的最优步长略有不同，量级上一致。

方向导数的表格验证了三条结论。解析值与数值逐位相符（偏差在 ${10}^{-10}$ 以内），说明 $D_{v}f=\mathrm{\nabla }f\cdot v$ 这条线性关系在实际计算中成立。数值的分布验证了极值性质：角度 $0^{\circ }$（沿 $x$ 轴）时 $2.1472$ 就等于 $\frac{\partial f}{\partial x}$，${90}^{\circ }$ 时 $2.2410$ 等于 $\frac{\partial f}{\partial y}$，两者都小于范数 $3.1036$；角度 ${60}^{\circ }$ 时 $3.0143$ 已经很接近最大值（梯度方向的角度是 $\arctan \frac{2.2410}{2.1472}={46.2}^{\circ }$，${60}^{\circ }$ 与它相差不大）；${120}^{\circ }$ 时降到 $0.8671$；${180}^{\circ }$ 与 ${270}^{\circ }$ 分别是负的偏导数值，对应下降方向。沿梯度方向的方向导数等于范数这一点在表里没有直接列出（表中取的是整数角度），但表中的最大值 $3.0143$ 与范数 $3.1036$ 的接近程度说明极值就在附近。

实验二：Hessian 矩阵与二阶 Taylor

核对解析 Hessian 与数值差分，测量二阶 Taylor 的误差阶，并用二次型的特征值判断极值类型。

import math
import numpy as np

def g(x, y):
    return x * x * y + math.sin(x * y)

def hess_analytic(x, y):
    return np.array([
        [2 * y - y * y * math.sin(x * y), 2 * x + math.cos(x * y) - x * y * math.sin(x * y)],
        [2 * x + math.cos(x * y) - x * y * math.sin(x * y), -x * x * math.sin(x * y)],
    ])

def hess_numeric(x, y, hh=1e-4):
    def gx(a, b):
        return 2 * a * b + b * math.cos(a * b)
    def gy(a, b):
        return a * a + a * math.cos(a * b)
    H = np.zeros((2, 2))
    H[0, 0] = (gx(x + hh, y) - gx(x - hh, y)) / (2 * hh)
    H[0, 1] = (gx(x, y + hh) - gx(x, y - hh)) / (2 * hh)
    H[1, 0] = (gy(x + hh, y) - gy(x - hh, y)) / (2 * hh)
    H[1, 1] = (gy(x, y + hh) - gy(x, y - hh)) / (2 * hh)
    return H

x0, y0 = 1.2, 0.7
ga = (2 * x0 * y0 + y0 * math.cos(x0 * y0), x0 * x0 + x0 * math.cos(x0 * y0))
Ha = hess_analytic(x0, y0)
Hn = hess_numeric(x0, y0)
print(f"  数值差分的最大偏差 = {float(np.max(np.abs(Ha - Hn))):.3e}")
print(f"  对称性偏差 = {abs(Ha[0, 1] - Ha[1, 0]):.3e}")
print(f"  特征值 = {np.linalg.eigvalsh(Ha)[0]:.8f}, {np.linalg.eigvalsh(Ha)[1]:.8f}")

base = g(x0, y0)
rng = np.random.default_rng(1)
dirv = rng.normal(size=2)
dirv /= np.linalg.norm(dirv)
print(f"{'|h|':>10} {'误差':>14} {'误差/|h|^3':>14}")
for sc in [1e-1, 5e-2, 2e-2, 1e-2, 5e-3]:
    hx, hy = sc * dirv[0], sc * dirv[1]
    approx = (base + ga[0] * hx + ga[1] * hy
              + 0.5 * (Ha[0, 0] * hx * hx + 2 * Ha[0, 1] * hx * hy + Ha[1, 1] * hy * hy))
    err = abs(g(x0 + hx, y0 + hy) - approx)
    print(f"{sc:>10.0e} {err:>14.4e} {err / sc ** 3:>14.6f}")

for name, Q in [("x^2 + 3xy + 5y^2", np.array([[1.0, 1.5], [1.5, 5.0]])),
                ("x^2 - y^2", np.array([[1.0, 0.0], [0.0, -1.0]])),
                ("x^2 - 2xy + y^2", np.array([[1.0, -1.0], [-1.0, 1.0]]))]:
    ev = np.linalg.eigvalsh(Q)
    if np.all(ev > 0):
        kind = "正定，局部极小"
    elif np.all(ev < 0):
        kind = "负定，局部极大"
    elif np.all(ev >= 0) or np.all(ev <= 0):
        kind = "半定，需要高阶信息"
    else:
        kind = "不定，鞍点"
    print(f"  {name:>16}: 特征值 ({ev[0]:>8.4f}, {ev[1]:>8.4f})  {kind}")

输出：

  数值差分的最大偏差 = 3.306e-09
  对称性偏差 = 0.000e+00
  特征值 = -2.67818145, 2.64102023
       |h|             误差       误差/|h|^3
     1e-01     5.3412e-04       0.534124
     5e-02     6.5839e-05       0.526714
     2e-02     4.1763e-06       0.522035
     1e-02     5.2044e-07       0.520440
     5e-03     6.4954e-08       0.519635
  x^2 + 3xy + 5y^2: 特征值 (  0.5000,   5.5000)  正定，局部极小
         x^2 - y^2: 特征值 ( -1.0000,   1.0000)  不定，鞍点
     x^2 - 2xy + y^2: 特征值 (  0.0000,   2.0000)  半定，需要高阶信息

解析 Hessian 在 $(1.2,0.7)$ 处是 $\left(\begin{array}{cc}1.035125& 2.441963\\ 2.441963& -1.072286\end{array}\right)$，与数值差分一致（最大偏差 $3.31\times {10}^{-9}$）。对称性偏差为零是结构性的：混合偏导数的对称性来自函数本身的二阶可微性，两个数值差分公式（一个先对 $x$ 再对 $y$，另一个反过来）用的是同一个网格与同一组 $f$ 的偏导数值，因此等号在浮点意义上精确成立。这条观察给出了检查 Hessian 实现的一种方式：若数值 Hessian 的对称性偏差明显大于截断误差量级，说明实现中有错误。

特征值 $-2.67818145$ 与 $2.64102023$ 一正一负，Hessian 不定。需要强调该点不是临界点（梯度范数 $3.1036$），所以这条信息不用于极值判别，而是说明该点附近存在曲率为负的方向；沿负曲率方向移动可以比梯度下降更快地降低函数值，这正是共轭梯度法、信赖域方法与鞍点逃离算法的出发点（第三章 3.5 节与 3.8 节）。

二阶 Taylor 的五组数据验证误差阶：$|h|$ 每次减半，误差降到约八分之一（$5.34\times {10}^{-4}\to 6.58\times {10}^{-5}$ 是 $\frac{1}{8.1}$，$2\times {10}^{-2}\to {10}^{-2}$ 是 $\frac{1}{8.0}$，${10}^{-2}\to 5\times {10}^{-3}$ 是 $\frac{1}{8.0}$），第三列稳定在 $0.52$（从 $0.534$ 缓慢下降到 $0.520$）。$O(\parallel h{\parallel }^3)$ 的结论成立，系数趋于常数。

二次型的三行数据把极值判别整理成表。$x^2+3xy+5y^2$ 的两个特征值为正，原点是严格局部极小（二次型正定意味着全局极小）；$x^2-y^2$ 的特征值一正一负，原点是鞍点（沿 $x$ 轴上升、沿 $y$ 轴下降）；$x^2-2xy+y^2=(x-y{)}^2$ 的一个特征值为零，二阶信息不足，直接观察函数形式可以看出沿直线 $x=y$ 上函数恒为零，因此整条直线都是极小点，属于退化极小（Hessian 在该方向为零，需要四阶项确认）。三种情形对应优化算法中的三类临界点，梯度下降在鞍点附近会变慢（第三章 3.3 节讨论了这一点）。

实验三：隐函数定理的数值验证

用解析公式与数值差分分别计算隐函数的导数，观察两者的一致性与公式失效的位置。

import math

print("(a) F(x, y) = x^2 + y^2 - 1 = 0，y(x) = sqrt(1 - x^2)")
for xt in [0.0, 0.3, 0.6, 0.8]:
    yt = math.sqrt(1 - xt * xt)
    dydx_analytic = -xt / yt
    eps = 1e-6
    yp = math.sqrt(1 - (xt + eps) ** 2)
    ym = math.sqrt(1 - (xt - eps) ** 2)
    dydx_num = (yp - ym) / (2 * eps)
    print(f"    x = {xt:>4.1f}: y = {yt:.10f}, dy/dx 解析 = {dydx_analytic:>12.8f}, "
          f"数值 = {dydx_num:>12.8f}, 偏差 = {abs(dydx_analytic - dydx_num):.2e}")

xa, ya = 0.3, 0.4
za = math.sqrt(1 - xa * xa - ya * ya)
zx_a, zy_a = -xa / za, -ya / za
eps = 1e-6
zx_n = (math.sqrt(1 - (xa + eps) ** 2 - ya ** 2) - math.sqrt(1 - (xa - eps) ** 2 - ya ** 2)) / (2 * eps)
zy_n = (math.sqrt(1 - xa ** 2 - (ya + eps) ** 2) - math.sqrt(1 - xa ** 2 - (ya - eps) ** 2)) / (2 * eps)
print(f"    在 ({xa}, {ya}) 处：z = {za:.10f}")
print(f"    (dz/dx, dz/dy) 解析 = ({zx_a:>12.8f}, {zy_a:>12.8f})")
print(f"    (dz/dx, dz/dy) 数值 = ({zx_n:>12.8f}, {zy_n:>12.8f})")

print("(c) 方程组 u^2 + v^2 - 1 = 0, u - v - x = 0")
print(f"{'x':>8} {'u':>14} {'v':>14} {'du/dx':>12} {'dv/dx':>12}")
for xt in [0.0, 0.2, 0.4, 0.5]:
    ut = (xt + math.sqrt(2 - xt * xt)) / 2.0
    vt = ut - xt
    dudx = (1 - xt / math.sqrt(2 - xt * xt)) / 2.0
    print(f"{xt:>8.1f} {ut:>14.10f} {vt:>14.10f} {dudx:>12.8f} {dudx - 1.0:>12.8f}")

输出：

    x =  0.0: y = 1.0000000000, dy/dx 解析 =  -0.00000000, 数值 =   0.00000000, 偏差 = 0.00e+00
    x =  0.3: y = 0.9539392014, dy/dx 解析 =  -0.31448545, 数值 =  -0.31448545, 偏差 = 2.60e-12
    x =  0.6: y = 0.8000000000, dy/dx 解析 =  -0.75000000, 数值 =  -0.75000000, 偏差 = 4.93e-11
    x =  0.8: y = 0.6000000000, dy/dx 解析 =  -1.33333333, 数值 =  -1.33333333, 偏差 = 1.98e-11
    在 (0.3, 0.4) 处：z = 0.8660254038
    (dz/dx, dz/dy) 解析 = ( -0.34641016,  -0.46188022)
    (dz/dx, dz/dy) 数值 = ( -0.34641016,  -0.46188022)
       x              u              v        du/dx        dv/dx
     0.0   0.7071067812   0.7071067812   0.50000000  -0.50000000
     0.2   0.8000000000   0.6000000000   0.42857143  -0.57142857
     0.4   0.8782329983   0.4782329983   0.35255804  -0.64744196
     0.5   0.9114378278   0.4114378278   0.31101776  -0.68898224

三组数据都验证了 $\frac{\partial \phi }{\partial x}=-{\left(\frac{\partial F}{\partial y}\right)}^{-1}\frac{\partial F}{\partial x}$ 这条公式。第一组（圆）的偏差在 ${10}^{-11}$ 到 ${10}^{-12}$ 量级（$x=0$ 处精确为零，因为公式与差分都给出零）；导数值从 $0$ 增到 $-1.3333$，绝对值随 $x$ 增大而增大，原因在公式 $-\frac{x}{y}$ 里：$x$ 增大时分子增大、分母 $y=\sqrt{1-x^2}$ 减小，两者共同作用。曲线在 $x\to 1$ 处趋于竖直（导数趋于负无穷），而 $\frac{\partial F}{\partial y}=2y\to 0$，正是定理条件失效的边界。

第二组（球面）在 $(0.3,0.4)$ 处的两个偏导数与数值差分逐位一致，偏差 $6.9\times {10}^{-11}$ 与 $5.1\times {10}^{-13}$（后者更小是因为 $y$ 方向的差分涉及的函数值变化更平缓）。$(\frac{\partial z}{\partial x},\frac{\partial z}{\partial y})=(-\frac{x}{z},-\frac{y}{z})$ 的比值 $\frac{0.3}{0.4}=0.75$ 与实测的 $\frac{0.34641016}{0.46188022}=0.75$ 相符，这条关系在约束优化中表现为「约束曲面切空间内的下降方向」。

第三组（方程组）的数据展示了矩阵形式的公式。$x=0$ 处的 $\frac{du}{dx}=0.5$ 与手工用 $J^{-1}$ 计算的结果一致；数据表的第四列与第五列之差恒为 $1$（$\frac{du}{dx}-\frac{dv}{dx}=1$），这是第二个方程 $u-v=x$ 求导的直接结果，可以作为独立的一致性检查。随着 $x$ 从 $0$ 增到 $0.5$，$\frac{du}{dx}$ 从 $0.5$ 降到 $0.31101776$，$\frac{dv}{dx}$ 从 $-0.5$ 降到 $-0.68898224$，两者的变化方向相反，符合几何直观：$x$ 增大时解沿着约束圆向右下方移动，$v$ 的变化幅度更大。

实验四：反向模式自动微分

用计算图实现反向模式自动微分，与解析梯度、数值梯度对照，并观察梯度计算代价随维数的变化。

import math
import time

class Var:
    __slots__ = ("v", "parents", "local", "g")
    def __init__(self, v, parents=(), local=()):
        self.v = v
        self.parents = parents
        self.local = local
        self.g = 0.0

def vadd(a, b):
    return Var(a.v + b.v, (a, b), (1.0, 1.0))

def vmul(a, b):
    return Var(a.v * b.v, (a, b), (b.v, a.v))

def vexp(a):
    return Var(math.exp(a.v), (a,), (math.exp(a.v),))

def vsin(a):
    return Var(math.sin(a.v), (a,), (math.cos(a.v),))

def backward(out):
    order, seen = [], set()
    def build(node):
        if id(node) in seen:
            return
        seen.add(id(node))
        for p in node.parents:
            build(p)
        order.append(node)
    build(out)
    out.g = 1.0
    for node in reversed(order):
        for p, l in zip(node.parents, node.local):
            p.g += node.g * l

def make_tape(xs):
    """f(x) = sum_i exp(x_i) * sin(x_i)"""
    x = [Var(v) for v in xs]
    terms = [vmul(vexp(xi), vsin(xi)) for xi in x]
    out = terms[0]
    for t in terms[1:]:
        out = vadd(out, t)
    return x, out

xs = [0.3, -0.7, 1.1]
xv, out = make_tape(xs)
backward(out)
g_ad = [xi.g for xi in xv]
g_exact = [math.exp(v) * (math.sin(v) + math.cos(v)) for v in xs]
print(f"  f = {out.v:.12f}")
print(f"  反向模式梯度 = ({g_ad[0]:.12f}, {g_ad[1]:.12f}, {g_ad[2]:.12f})")
print(f"  解析梯度     = ({g_exact[0]:.12f}, {g_exact[1]:.12f}, {g_exact[2]:.12f})")
print(f"  最大偏差     = {max(abs(a - b) for a, b in zip(g_ad, g_exact)):.3e}")

def f_plain(xs):
    s = 0.0
    for v in xs:
        s += math.exp(v) * math.sin(v)
    return s

def grad_fd(xs, hh=1e-6):
    g = []
    for i in range(len(xs)):
        xp = list(xs); xm = list(xs)
        xp[i] += hh; xm[i] -= hh
        g.append((f_plain(xp) - f_plain(xm)) / (2 * hh))
    return g

xs_big = [0.01 * i - 3.0 for i in range(200)]
xv_big, out_big = make_tape(xs_big)
backward(out_big)
g_ad_big = [xi.g for xi in xv_big]
g_fd_big = grad_fd(xs_big)
print(f"\n维数 200 时两种方法的一致性：最大偏差 = "
      f"{max(abs(a - b) for a, b in zip(g_ad_big, g_fd_big)):.3e}")

print("梯度计算的代价随维数的变化（同一实现下的计时）")
print(f"{'n':>6} {'反向模式 (ms)':>16} {'有限差分 (ms)':>16} {'代价比':>10}")
for n in [10, 50, 200, 800]:
    xt = [0.01 * i - 2.0 for i in range(n)]
    t0 = time.perf_counter()
    for _ in range(20):
        xv2, out2 = make_tape(xt)
        backward(out2)
    t1 = time.perf_counter()
    for _ in range(20):
        grad_fd(xt)
    t2 = time.perf_counter()
    tr = (t1 - t0) / 20 * 1000
    tf = (t2 - t1) / 20 * 1000
    print(f"{n:>6} {tr:>16.3f} {tf:>16.3f} {tf / tr:>10.1f}")

输出：

  f = 2.756336389241
  反向模式梯度 = (1.688479927823, 0.059900354000, 4.040012927968)
  解析梯度     = (1.688479927823, 0.059900354000, 4.040012927968)
  最大偏差     = 0.000e+00

维数 200 时两种方法的一致性：最大偏差 = 1.741e-09
     n        反向模式 (ms)        有限差分 (ms)        代价比
    10            0.032            0.013        0.4
    50            0.138            0.296        2.1
   200            0.535            4.000        7.5
   800            2.383           68.789       28.9

第一段数据验证反向模式的正确性：三个分量的梯度与解析值逐位一致（最大偏差精确为零），与 5.4 节前向模式双数的结果性质相同。反向模式的实现只有二十余行，核心是两步：前向构造计算图（每个运算记录父节点与局部偏导），反向按逆拓扑序累加伴随量。累加的规则是链式法则：节点把它的伴随量乘以局部偏导后加到每个父节点上。拓扑序的作用是保证一个节点被处理时它所有的下游贡献都已经累加完毕（同一个节点可能被多个下游节点使用，例如本例中每个 $x_i$ 同时进入 $\exp$ 与 $\sin$ 两个分支）。

第二段数据把代价的结构写清楚了。维数 $n=200$ 时反向模式与中心差分的一致性到 $1.74\times {10}^{-9}$（与中心差分本身的截断误差量级一致）。耗时的变化是本节最有实用价值的一组数据：$n=10$ 时反向模式 $0.032$ 毫秒、有限差分 $0.013$ 毫秒，代价比 $0.4$（反向模式反而更慢，因为构造计算图的对象开销在 $n$ 小时占主导）；$n=50$ 时比值 $2.1$；$n=200$ 时 $7.5$；$n=800$ 时 $28.9$。比值随 $n$ 近似线性增长，原因可以直接从复杂度看出：反向模式的代价与计算图中的运算个数成正比，本例中每个变量对应 4 个运算（$\exp$、$\sin$、乘、加），总量 $O(n)$；有限差分需要 $2n$ 次前向计算，每次前向的代价是 $O(n)$，总量 $O(n^2)$。两者的比值因此按 $n$ 增长。

这组数据的实用含义是一个明确的判据：变量少时用有限差分（实现简单、常数小），变量多时用反向模式自动微分（代价与变量个数几乎无关）。深度学习中损失函数有 ${10}^6$ 到 ${10}^9$ 个参数而输出是标量，反向模式的代价与参数个数无关（只与前向计算的代价同阶），这是它在框架中成为标准做法的全部原因。5.4 节的前向模式在「输入少、输出多」时占优，两种模式的选择依据是输入个数与输出个数哪个更大。

5.8.8 本节小结

多元微分的四层内容

本节把多元微分学组织为四层。微分层给出偏导数、方向导数、全微分与梯度：可微是比「偏导数存在」更强的条件（反例 $\frac{xy}{\sqrt{x^2+y^2}}$ 的所有偏导数存在但不可微），可微时 $D_{v}f=\mathrm{\nabla }f\cdot v$，方向导数关于方向线性，这使 $n$ 个偏导数成为完整的一阶信息；梯度的三条几何含义分别是上升最快的方向、与水平集正交、范数等于最大方向导数（实验一：梯度 $(2.1472,2.2410)$、范数 $3.1036$，方向导数从 $0^{\circ }$ 的 $2.1472$ 到 ${120}^{\circ }$ 的 $0.8671$，与解析公式的偏差在 ${10}^{-10}$ 以内）。二阶层给出 Hessian 矩阵与二阶 Taylor：混合偏导数的对称性来自 Clairaut-Schwarz 定理（实验二的数值 Hessian 对称性偏差精确为零），二阶 Taylor 的误差是 $O(\parallel h{\parallel }^3)$（实测误差除以 $\parallel h{\parallel }^3$ 稳定在 $0.52$），极值判别由特征值符号给出（正定极小、负定极大、不定鞍点、半定需要高阶信息，实验二的三行数据分别对应正定、不定与半定）。隐函数层给出反函数定理与隐函数定理，核心公式是 $\frac{\partial \phi }{\partial x}=-{\left(\frac{\partial F}{\partial y}\right)}^{-1}\frac{\partial F}{\partial x}$，它把隐式关系的求导化为矩阵求逆（实验三验证了一元隐函数、二元隐函数与方程组三种情形，偏差在 ${10}^{-11}$ 量级；$y=0$ 与 $x=\pm 1$ 处 $\frac{\partial F}{\partial y}=0$ 对应公式失效，与定理条件一致）；这条定理是 Lagrange 乘子法与 KKT 条件的理论来源。微分形式层把梯度、旋度、散度统一为外微分算子 $d$ 在不同阶形式上的作用，$d(d\omega )=0$ 给出梯度场无旋与旋度场无散，Poincaré 引理说明单连通区域上闭合蕴含恰当（$\omega =\frac{-ydx+xdy}{x^2+y^2}$ 是拓扑障碍的标准例子），广义 Stokes 公式 ${\int }_{\partial M}\omega ={\int }_{M}d\omega$ 把微积分基本定理与三大积分公式统一为一条。计算层的四组实验给出完整的数值对照：解析与数值的梯度、Hessian、隐函数导数都一致到截断误差量级；反向模式自动微分的梯度与解析值逐位一致（偏差为零），代价与变量个数基本无关，$n=800$ 时比有限差分快 $28.9$ 倍，而 $n=10$ 时反而更慢。

多元的微分说清楚了，下一节处理多元的积分。重积分的换元与 Jacobian 行列式、曲线曲面积分、Green 与 Gauss 与 Stokes 公式、以及高维积分的数值方法（网格法与 Monte Carlo）构成了下一节的内容，其间会用到本节的全部结论。

练习题

第 1 题 概念推导

(a) 说明为什么「所有偏导数存在」不蕴含可微，并举出反例。(b) 证明梯度与水平集正交。(c) 说明隐函数定理中条件 $\frac{\partial F}{\partial y}\ne 0$ 的必要性，用圆的例子说明条件失效时会发生什么。

参考答案

(a) 偏导数只沿 $n$ 个坐标方向考察函数的行为，可微要求函数在所有方向上都与同一个线性函数相差 $o(\parallel h\parallel )$。反例是

$$f(x,y)=\{\begin{array}{ll}{\displaystyle \frac{xy}{\sqrt{x^2+y^2}}},& (x,y)\ne (0,0)\\ 0,& (x,y)=(0,0)\end{array}$$

沿 $x$ 轴 $f(x,0)=0$，故 $\frac{\partial f}{\partial x}(0,0)=0$；同理 $\frac{\partial f}{\partial y}(0,0)=0$。若可微，全微分必须是零线性映射，要求 $\frac{f(h_1,h_2)}{\parallel h\parallel }\to 0$。取 $h_1=h_2=t$：$\frac{f(t,t)}{\sqrt{2}t}=\frac{t^2/(\sqrt{2}t)}{\sqrt{2}t}=\frac{1}{2}$，不趋于零，故不可微。这说明偏导数只提供 $n$ 个方向的信息，不足以确定函数的局部线性行为。

(b) 设水平集是 $\{x:f(x)=c\}$，取其中一条可微曲线 $\gamma (t)$ 满足 $f(\gamma (t))=c$ 与 $\gamma (0)=x_0$。对 $t$ 求导（链式法则）：

$$\mathrm{\nabla }f(\gamma (0))\cdot {\gamma }^{\prime }(0)=0$$

${\gamma }^{\prime }(0)$ 是水平集在 $x_0$ 处的任意切向量，因此梯度与水平集的所有切向量正交，即梯度垂直于水平集。推论是 $\mathrm{\nabla }f(x_0)=0$ 时水平集在该点可能退化（不再是光滑曲面），这是临界点的几何特征。等值线的疏密程度也由梯度的范数决定：$\parallel \mathrm{\nabla }f\parallel$ 大处等值线密（函数变化快），小处疏。

(c) 公式 $\frac{dy}{dx}=-\frac{F_x}{F_y}$ 在 $F_y=0$ 处失效，原因不只是分母为零：此时隐函数可能根本不存在（无法在 $x_0$ 附近把 $y$ 表示为 $x$ 的函数）。以 $F(x,y)=x^2+y^2-1$ 为例，$F_y=2y$ 在 $y=0$ 处为零，对应点 $(\pm 1,0)$。在 $x_0=1$ 附近，$F(x,y)=0$ 的解是 $y=\pm \sqrt{1-x^2}$：对每个接近 $1$ 的 $x<1$，有两个 $y$ 值，且这两个值在 $x\to 1$ 时都趋于零并在一起形成竖直切线，无法在 $x=1$ 附近选出唯一的连续分支。这正是定理的结论：$F_y$ 可逆是「唯一分支存在」的充分条件，条件失效处的解集可能出现分支合并或消失。

第 2 题 计算推理

(a) 求 $f(x,y,z)=x^{2}yz+\ln (1+x^2)$ 在 $(1,2,3)$ 处的梯度、Hessian 与沿方向 $v=\frac{1}{3}(2,2,1)$ 的方向导数。(b) 对隐式曲面 $F(x,y,z)=x^2+2y^2+3z^2-6=0$，求在点 $(1,1,1)$ 处的切平面方程与法向量。(c) 说明为什么 Jacobian 行列式在换元积分中取绝对值。

参考答案

(a) 一阶偏导数：

$$\frac{\partial f}{\partial x}=2xyz+\frac{2x}{1+x^2},\frac{\partial f}{\partial y}=x^{2}z,\frac{\partial f}{\partial z}=x^{2}y$$

在 $(1,2,3)$ 处：$\frac{\partial f}{\partial x}=12+1=13$，$\frac{\partial f}{\partial y}=3$，$\frac{\partial f}{\partial z}=2$，梯度为 $(13,3,2)$。二阶偏导数：$\frac{{\partial }^{2}f}{\partial x^2}=2yz+\frac{2(1-x^2)}{(1+x^2{)}^2}$（在 $(1,2,3)$ 处为 $12$），$\frac{{\partial }^{2}f}{\partial y^2}=0$，$\frac{{\partial }^{2}f}{\partial z^2}=0$，$\frac{{\partial }^{2}f}{\partial x\partial y}=2xz=6$，$\frac{{\partial }^{2}f}{\partial x\partial z}=2xy=4$，$\frac{{\partial }^{2}f}{\partial y\partial z}=x^2=1$。Hessian 是

$$H=\left(\begin{array}{ccc}12& 6& 4\\ 6& 0& 1\\ 4& 1& 0\end{array}\right)$$

方向导数：$D_{v}f=\mathrm{\nabla }f\cdot v=\frac{1}{3}(13\times 2+3\times 2+2\times 1)=\frac{34}{3}=11.3333$。这个值小于 $\parallel \mathrm{\nabla }f\parallel =\sqrt{169+9+4}=\sqrt{182}=13.4907$，符合方向导数的极值性质。注意该 Hessian 不是半正定（对角元 $(2,2)$ 与 $(3,3)$ 为零），说明该点不是极值点的候选。

(b) 由隐函数定理，曲面 $F=0$ 的法向量是 $\mathrm{\nabla }F=(2x,4y,6z)$，在 $(1,1,1)$ 处为 $(2,4,6)$。单位法向量是 $\frac{1}{\sqrt{56}}(2,4,6)=\frac{1}{\sqrt{14}}(1,2,3)$。切平面方程：$2(x-1)+4(y-1)+6(z-1)=0$，化简得

$$x+2y+3z=6$$

这条结论的用法是把「切平面」化为「梯度与切向量正交」这个线性条件。切向量可以取为 $\frac{\partial }{\partial x}(x,y(x,z),z)$ 的形式，通过隐函数求导得到，结果给出同一个切平面。

(c) Jacobian 行列式的绝对值 $|detJ|$ 是局部体积的膨胀率，体积本身非负，因此换元公式中取绝对值：

$${\int }_{\mathrm{\Omega }}f(x)dx={\int }_{{\mathrm{\Phi }}^{-1}(\mathrm{\Omega })}f(\mathrm{\Phi }(u))|det{J}_{\mathrm{\Phi }}(u)|du$$

行列式的符号记录方向：$detJ>0$ 时变换保持定向，$detJ<0$ 时反转定向（例如把 $x$ 换成 $-x$，或交换两个坐标轴）。不带绝对值的版本在微分形式的积分中自然出现（$\int \omega$ 对定向敏感，$k$ 形式在定向反转时变号），这正是微分形式语言的便利之处：它把 Jacobian 行列式的符号自动处理在 $\wedge$ 的外积规则里。数值计算中取绝对值，避免把定向信息误当作体积信息。

第 3 题 代码验证

用本节代码做三组实验并解释：(a) 在实验一里把函数换成 $f(x,y)=e^{xy}\ln (x^2+y^2+1)$，检查梯度与数值差分的一致性，并找出方向导数取最大值的角度；(b) 在实验二里求 $g(x,y)=x^{2}y+\sin (xy)$ 的临界点，用 Hessian 特征值判断该点的类型；(c) 在实验四里把 $\exp$ 与 $\sin$ 换成需要更多运算的函数（例如 $\tanh$ 与 $\log$），观察反向模式的时间如何随计算图规模增长。

参考答案

(a) 解析梯度为

$$\frac{\partial f}{\partial x}=y{e}^{xy}\ln (x^2+y^2+1)+\frac{2x{e}^{xy}}{x^2+y^2+1},\frac{\partial f}{\partial y}=x{e}^{xy}\ln (x^2+y^2+1)+\frac{2y{e}^{xy}}{x^2+y^2+1}$$

在 $(1.2,0.7)$ 处，中心差分（步长 ${10}^{-6}$）与解析值的偏差应在 ${10}^{-10}$ 量级（与实验一相同，两条公式都用中心差分格式，误差量级只依赖步长与函数的光滑程度）。最大方向导数的角度是 $\arctan \frac{\partial f/\partial y}{\partial f/\partial x}$ 对应的方向（按象限调整），在该角度处的方向导数等于梯度的范数。验证方式是在 $[0^{\circ },{360}^{\circ })$ 上以 $1^{\circ }$ 的步长扫描方向导数并找最大值，最大值应当出现在上述角度附近且数值等于 $\parallel \mathrm{\nabla }f\parallel$。这条练习说明最速上升方向不限于坐标轴方向，对一般函数它可能落在任意角度。

(b) 求临界点需要解 $\mathrm{\nabla }g=0$：

$$2xy+y\cos (xy)=0,x^2+x\cos (xy)=0$$

$y\ne 0$ 时第一个方程给出 $2x+\cos (xy)=0$，$x\ne 0$ 时第二个方程给出 $x+\cos (xy)=0$，两式相减得 $x=0$，与 $x\ne 0$ 矛盾。因此临界点落在坐标轴上：$x=0$ 时两个方程都化为 $0=0$（第二个方程 $x^2+x\cos (xy)=0$ 自动满足），故整条直线 $x=0$ 上都是临界点。这类退化临界点集合（整条直线而非孤立点）是 Hessian 半定的典型情形，二阶判别失效，需要更高阶的展开或直接考察函数结构。这条练习的结论是：二阶判别只适用于孤立的非退化临界点，遇到连续多解的临界点集合时应当先做结构分析。

(c) 反向模式的时间与计算图的运算个数成正比，与变量的个数无关（变量对图规模的贡献是它们参与的运算数）。把 $\exp$ 与 $\sin$ 换成 $\tanh$ 与 $\log$ 后，每个变量对应的运算数增加（$\tanh$ 的实现需要两次 $\exp$ 加一次除法，$\log$ 的导数 $\frac{1}{x}$ 需要一次除法），因此时间按运算个数同比例增长。变量数固定时，时间随每个变量的运算个数线性增长，比例系数就是单次运算的平均代价（对象创建、属性访问、浮点运算）。这条结论对框架实现有直接指导：加快反向传播主要是减少计算图中的节点数与降低单节点开销，与变量个数没有直接关系。

第 4 题 综合应用

某数值程序需要求解一个由 $n=500$ 个变量、$m=500$ 个方程构成的非线性方程组 $F(x)=0$（$F$ 的 Jacobian 矩阵稠密，可解析求导）。请比较三种求解方案并给出推荐：Newton 法、拟 Newton 法（BFGS）、固定点迭代。说明各自的每步代价、收敛阶与适用条件。

参考答案

方案一：Newton 法。每步求解线性方程组 $J(x_k)\mathrm{\Delta }x=-F(x_k)$，$J$ 是 $500\times 500$ 的稠密矩阵。代价分三部分：构造 Jacobian（${500}^2=2.5\times {10}^5$ 个元素，解析求导的代价约 $O(n^2)$）；LU 分解（$O(n^3)\approx 1.25\times {10}^8$ 次运算，稠密矩阵的常数约 $\frac{2}{3}$，实际约 $8\times {10}^7$）；回代（$O(n^2)$）。总代价由 LU 分解主导，单步约 ${10}^8$ 次浮点运算，毫秒级到数十毫秒级。收敛是二次的（解附近误差满足 $e_{k+1}\approx C{e}_k^2$，5.4 节的实验三给出了量化的例子），从合理初值出发通常不超过 $20$ 步。条件：$J(x_0)$ 可逆（由隐函数定理保证局部存在解），初值足够接近解。

方案二：拟 Newton 法（BFGS）。不计算 Jacobian，用函数值之差构造近似的 Jacobian（逆）的更新。每步代价 $O(n^2)$（约 $2.5\times {10}^5$ 次运算，比 Newton 法低三个数量级），收敛是超线性（比线性快、比二次慢，典型阶数在 $1.5$ 到 $1.8$ 之间），步数比 Newton 法多（通常 $30$ 到 $100$ 步）。总代价通常低于 Newton 法（单步便宜三个数量级，步数只多几倍）。条件是：不需要解析 Jacobian（当 Jacobian 难以计算时这是决定性优势），但需要存储 $n\times n$ 的近似矩阵（内存 $O(n^2)$）。

方案三：固定点迭代（例如把方程整理成 $x=G(x)$）。每步代价 $O(n)$，是三者中最低的，但收敛是线性的，且收敛常数可能接近 $1$（甚至在 $1$ 以上发散）。收敛条件由压缩映射原理（5.1 节）给出：$\parallel J_G\parallel <1$ 在迭代区域内成立。对一般的非线性方程组这个条件通常不满足，需要把方程重新整理（预条件）才能使用。

推荐：Jacobian 可解析求得、$n=500$ 时用 Newton 法。依据是稠密 LU 分解的代价 $O(n^3)) 在这个规模下可接受（通常 $n$ 上万时才需要改用迭代法求解线性方程组），二次收敛把步数压到几十步以内，总代价与精度都最优。若 Jacobian 只能数值近似或 $n$ 增大到 ${10}^4$ 以上，改用拟 Newton 法（单步 $O(n^2)$ 且不需要 Jacobian），并考虑有限内存版本把存储降到 $O(n)$。固定点迭代只适用于已知压缩性质的场合。实现细节上要注意两点：Newton 法需要阻尼或信赖域（防止从远处初值出发时步长过大而发散）；判断收敛时同时检查残差 $\parallel F(x_k)\parallel$ 与步长 $\parallel \mathrm{\Delta }{x}_k\parallel$，两者都小才停止。

常见错误

错误 1 · 用偏导数存在作为可微的判据

原因：偏导数只沿坐标方向定义，$n$ 个偏导数存在不足以确定函数在所有方向上的局部线性行为。反例 $\frac{xy}{\sqrt{x^2+y^2}}$ 的所有偏导数在原点存在（都为零），但沿对角线方向 $\frac{f}{\parallel h\parallel }\to \frac{1}{2}$，不可微。若在算法中直接用偏导数构造梯度并进行梯度下降，函数可能在某些点处不满足下降方向的估计，导致步长分析失效。

解决：区分三个条件（偏导数存在、可微、偏导数连续），后者蕴含前者、前者不蕴含后者。可微性的验证方法是检查 $f(x_0+h)-f(x_0)-\mathrm{\nabla }f(x_0)\cdot h=o(\parallel h\parallel )$，实际做法是取几个方向（坐标轴、对角线、随机方向）比较差商与方向导数。工程实用的充分条件是偏导数连续，函数由初等函数复合而成时这条一般成立。涉及分段定义的函数（ReLU、绝对值、最大值）要单独检查分段边界处的可微性。

错误 2 · 在非临界点用 Hessian 的特征值判断极值

原因：二阶判别的前提是 $\mathrm{\nabla }f(x_0)=0$。在梯度不为零的点，Hessian 的特征值只描述局部曲率（沿各方向的二阶变化率），与「该点是不是极值点」无关。本节实验二在 $(1.2,0.7)$ 处算出 Hessian 的特征值 $-2.67818145$ 与 $2.64102023$，若据此判断「该点是鞍点」就错了：该点的梯度范数是 $3.1036$，函数在这里沿负梯度方向仍在下降，它不满足极值点的必要条件。这类错误在调试优化程序时会造成误判（把迭代中间的某个点当作收敛点）。

解决：判断极值前先检查梯度范数（是否小于容差），确认是临界点后再看 Hessian 的特征值。梯度不为零时 Hessian 的信息用于其他目的：负特征值给出可以加速下降的方向（信赖域与鞍点逃离），特征值的条件数给出问题的病态程度（第三章 3.5 节的二阶方法）。

错误 3 · 忽略隐函数定理的条件直接套公式

原因：$\frac{dy}{dx}=-\frac{F_x}{F_y}$ 只在 $F_y\ne 0$（更一般地，关于待解变量的 Jacobian 可逆）时成立。条件失效处的表现形式有三类：解不存在、解不唯一（分支合并）、导数无穷大。圆的例子中 $F_y=2y$ 在 $y=0$ 处为零，对应圆的左右端点，那里 $y$ 不是 $x$ 的单值函数。数值计算中忽略这条条件会在奇异点附近得到巨大的导数（$F_y$ 很小时公式给出很大的值），进而使迭代步长失控。

解决：套用公式前检查 Jacobian 是否可逆（计算行列式或条件数）。条件接近失效时用正则化或改用参数化方式（圆的情形改用角度参数，避免在奇点处求显式函数）。在方程组情形，若 Jacobian 奇异，说明约束之间存在线性相关（第三章 3.6 节的约束资格条件讨论的正是这种情形），需要检查约束的线性无关性。

错误 4 · 把反向模式与前向模式的适用场景搞混

原因：两种模式的代价结构不同，选错会在高维问题上付出巨大代价。前向模式的代价与输入个数成正比（每个输入需要一次传播），适合输入少、输出多；反向模式的代价与输出个数成正比（每个输出需要一次反向传播），适合输入多、输出少。本节实验四的数据给出量化的对照：$n=800$ 时反向模式 $2.383$ 毫秒、有限差分（可视为代价更高的前向型方法）$68.789$ 毫秒，比值 $28.9$；而 $n=10$ 时反向模式 $0.032$ 毫秒、有限差分 $0.013$ 毫秒，反向模式反而慢。若把梯度下降这样的多输入单输出问题交给前向模式，代价会随参数个数线性增长（${10}^6$ 个参数需要 ${10}^6$ 次传播），不可行。

解决：按输入与输出的个数选择。损失函数（参数极多、输出一个标量）用反向模式；计算 Jacobian 的全部元素（输入少、输出多）用前向模式；$n$ 很小时（不超过几十）用有限差分，实现最简单且常数因子小。报告中写明所用方法的代价量级与它依赖的规模参数。