做数值计算的人大概都有过这么一次经历手头有个函数 $f(x,y)$ 在某点附近变化很小想用低阶多项式替它做快速估算于是顺手把一元泰勒公式的套路搬过来写成 $f \approx f_0 f_x h f_y k$觉得很合理。可一旦你把展开阶数提高到二阶问题就来了——交叉项前面那个系数到底是 1 还是 2符号表里查不到教材上一晃而过自己推又容易在两个变量的求导顺序里绕晕。更麻烦的是这个系数错了不一定报错程序照跑结果照出只是误差估计悄悄偏了一倍等到某个模型开始偶尔不准的时候你已经很难回溯到这一步了。二元函数的泰勒展开就是干这件事的把两个变量的函数在某个点的邻域内用一个多项式去逼近它逼近的好坏由阶数和扰动幅度决定。它是误差传播分析、Hessian 极值判定、牛顿迭代、灵敏度分析、数值线性化这些工作的共同底座。下面我按自己教学和工程中反复用到的顺序把这条公式从为什么长这样一直讲到怎么算才不会错中间穿插两道手算例题和一次数值复核。不管你是正在啃高数的学生还是需要写误差分析代码的工程师看完应该都能直接上手。1. 先别急着背公式把二元展开压回一元的那条直线上1.1 一个把误差估计算错一倍的现场我最早踩的坑发生在做传感器标定的误差传递分析时。当时有两个输入量 $x$ 和 $y$各自带有小幅随机扰动我想估计输出量的偏差就写了一阶近似 $f_x \Delta x f_y \Delta y$这一层没问题。但当我想进一步估计二阶偏差时脑子里默认写了 $\frac{1}{2}(f_{xx}\Delta x^2 f_{xy}\Delta x\Delta y f_{yy}\Delta y^2)$把交叉项系数漏掉了。算出来的偏差比实测小了一截我一度怀疑是测量噪声的问题换了几组数据都对不上。后来把公式老实展开一遍才发现交叉项前面确实是 2不是 1。这不是排版笔误而是两个求导路径叠加的结果。这个错误的隐蔽性在于它不影响一阶项只影响二阶及以上如果你的模型里一阶项占主导你根本发现不了但只要扰动稍微大一点二阶项开始起作用误差就浮出来了。所以我现在的习惯是凡是涉及二元的二阶近似交叉项系数必须现场推一遍不靠记忆。1.2 沿直线取样把二元的难题降成一元理解二元展开最省脑子的办法不是去背那个长长的求和式而是把它降维回一元。具体做法是固定一个方向从展开点 $(x_0, y_0)$ 出发沿着向量 $(h, k)$ 走一条直线把点写成参数形式$$x x_0 th, \quad y y_0 tk, \quad t \in [0, 1]$$然后把 $f$ 沿着这条直线限制成一个单变量函数$$g(t) f(x_0 th,\ y_0 tk)$$这个操作的关键在于$g(t)$ 是一个货真价实的一元函数$t0$ 对应展开点$t1$ 对应我们要估算的目标点。于是二元展开的全部困难就转化成了对 $g(t)$ 在 $t0$ 处做一元泰勒展开然后取 $t1$这么一件事。用一句大白话说二元泰勒展开本质上就是沿各个方向做一元展开之后的结果被统一写成了一个式子。这个视角还有一个额外好处——它顺便解释了为什么二元展开的余项形式和一元几乎一模一样。因为 $g(t)$ 的余项就是一元余项只是里面的导数换成了沿方向 $(h,k)$ 的方向导数而已。你不需要为二元的余项单独发明一套理论。1.3 g(t) 求导交叉项的系数 2 是这么冒出来的现在把 $g(t)$ 求导让系数 2 自己现身。根据链式法则$$g(t) f_x(x_0th,\ y_0tk)\cdot h f_y(x_0th,\ y_0tk)\cdot k$$这一步是干净的一阶项就是 $f_x h f_y k$没有争议。麻烦在二阶。再对 $t$ 求一次导注意 $f_x$ 和 $f_y$ 本身也是二元函数还要继续用链式法则$$\frac{d}{dt}f_x(\cdot) f_{xx}h f_{xy}k$$$$\frac{d}{dt}f_y(\cdot) f_{yx}h f_{yy}k$$把这两个结果分别乘上外面的 $h$ 和 $k$$$g(t) (f_{xx}h f_{xy}k)h (f_{yx}h f_{yy}k)k f_{xx}h^2 f_{xy}hk f_{yx}hk f_{yy}k^2$$到这里答案已经很清楚了$hk$ 这一项出现了两次一次来自对 $f_x$ 求导、一次来自对 $f_y$ 求导。如果函数是二阶连续可微的那么 $f_{xy} f_{yx}$两项合并就是 $2f_{xy}hk$。这就是那个 2 的全部来历——它不是人为规定的系数而是两条求导路径的合并结果。我一直觉得只要你在纸上把这一段亲手写一遍这个 2 就再也忘不掉了。反过来如果只背结论很容易在写代码或者列公式的时候顺手把它省掉因为它看起来应该像一元那样是 1/2 系的整齐结构。数学的整齐有时候是假象这里是真有两份。2. 完整公式长什么样从一阶切平面到 n 阶算子写法2.1 记号约定h、k 是增量不是微分先把记号钉死很多混乱其实来自记号。设 $h x - x_0$$k y - y_0$它们是把目标点 $(x,y)$ 拉回到展开点 $(x_0,y_0)$ 的有限增量不是无穷小量。这一点很关键泰勒展开之所以有实用价值恰恰是因为它处理的是有限的、不太小的扰动而不是微分那种极限意义下的无穷小。展开点通常记成下标 0比如 $f(x_0,y_0)$ 写成 $f_0$$f_x(x_0,y_0)$ 写成 $f_{x0}$依此类推。不带下标的偏导数符号表示在展开点取值。写代码的时候我建议把展开点坐标和增量分成两组变量管理比如x0, y0和dx, dy避免在推导过程中把 $h$ 误当成 $x$ 的一部分去求导那是另一个高频错误。另外提醒一句$h$ 和 $k$ 是有量纲的它们的量纲和被展开的函数自变量一致。这一点在第 4 章会展开讲因为当两个自变量的物理量纲差别很大时比如一个是长度、一个是角度展开阶数的选择会变得很微妙。2.2 一阶展开与切平面近似一阶展开最直观公式是$$f(x_0h,\ y_0k) \approx f_0 f_{x0}h f_{y0}k$$几何意义特别干净右边是关于 $h,k$ 的线性函数在三维空间里画出来就是一个平面而且这个平面在点 $(x_0,y_0,f_0)$ 处与曲面相切所以叫切平面近似。曲面在这一点附近有多像这个平面取决于曲面弯曲得有多厉害。如果曲面在该点附近几乎是平的一阶就够如果弯得明显一阶会系统性偏低或偏高因为有二阶项在单边贡献。工程上用一阶展开的场景非常多因为它天然就是线性化。控制系统里把非线性环节在工作点附近线性化本质就是一阶泰勒展开电路的小信号模型也是一阶展开测量里的误差传递一阶公式还是它。一阶展开的好处不只是简单而是它让叠加原理成立——多个小扰动的作用可以直接相加。但一阶展开有个系统性的短板它给出的近似值与真实值之间的偏差符号是有规律的不随机。也就是说一阶近似在展开点的一侧偏高、另一侧偏低如果你做的是一组重复测量然后取平均这个偏差不会自动抵消会作为系统偏差留在结果里。要处理它就得上二阶。2.3 二阶展开的显式形式和余项把 $g(t)$ 和 $g(t)$ 在 $t0$ 处取值代入一元泰勒并令 $t1$得到二元二阶展开$$f(x_0h,\ y_0k) f_0 f_{x0}h f_{y0}k \frac{1}{2}\left(f_{xx0}h^2 2f_{xy0}hk f_{yy0}k^2\right) R_2$$把这三项的角色分清楚很有必要。$f_{xx0}h^2/2$ 是沿 $x$ 方向的纯曲率贡献$f_{yy0}k^2/2$ 是沿 $y$ 方向的纯曲率贡献而 $f_{xy0}hk$ 是交叉贡献它描述的是$x$ 方向变化会不会改变 $y$ 方向的斜率也就是两个变量之间的耦合程度。如果 $f_{xy0}0$说明在该点附近两个方向可以解耦处理如果 $f_{xy0}$ 很大说明耦合强烈你不能分别调 $x$ 和 $y$。余项 $R_2$ 用拉格朗日形式写是$$R_2 \frac{1}{6}\Big(h\frac{\partial}{\partial x} k\frac{\partial}{\partial y}\Big)^3 f(x_0\theta h,\ y_0\theta k), \quad 0\theta1$$注意展开成六个三阶偏导的线性组合后前系数是 $1/6$这个 $1/6$ 来自 $3!$。余项的意义在于给出误差的上界实际做误差估计时我们通常不知道 $\theta$但可以用区域内三阶偏导的最大值 $M_3$ 去卡一个保守界$$|R_2| \le \frac{M_3}{6}\left(|h| |k|\right)^3$$这个不等式在实际工作里非常好用因为它只需要估一个三阶导数的上界不需要知道具体的 $\theta$ 位置。2.4 n 阶通式算子幂写法要把结论推广到任意阶最省事的写法是引入方向导数算子$$D h\frac{\partial}{\partial x} k\frac{\partial}{\partial y}$$那么 $n$ 阶展开就是$$f(x_0h,\ y_0k) \sum_{m0}^{n}\frac{1}{m!}D^m f(x_0,y_0) R_n$$算子写法看着抽象但它的好处是把公式长什么样的问题一次性解决了。你把 $D^m$ 当作二项式展开处理就行$D^2$ 展开出 $h^2\partial_{xx} 2hk\partial_{xy} k^2\partial_{yy}$那两倍的系数自然来自二项式系数 $\binom{2}{1}2$$D^3$ 展开出 $h^3\partial_{xxx} 3h^2k\partial_{xxy} 3hk^2\partial_{xyy} k^3\partial_{yyy}$系数 3 来自 $\binom{3}{1}\binom{3}{2}3$。所以之前那个神秘的 2在算子视角下就是二项式系数。这个统一性很有价值三元、四元函数的展开公式不需要另记你只要把算子写成 $D h_1\partial_1 h_2\partial_2 h_3\partial_3$然后做多项式展开系数由多项式定理给出。我平时写多变量近似代码时就是按多项式系数批量生成各项而不是手工逐项列。顺带说一个实际取舍工程计算里用到四阶以上展开的情形很少因为高维情况下项数增长太快。二元函数到二阶有 6 个系数含常数、两个一阶、三个二阶到三阶再加 4 个三阶偏导一共 10 个到四阶再加 5 个是 15 个。三元的话二阶就有 10 项。所以实践中二元通常停在二阶或三阶需要更高精度时改换别的数值手段而不是硬堆阶数。3. 两道手算例题从原点展开到非原点展开3.1 例一$e^x \ln(1y)$ 在原点展到二阶选这个函数是因为它足够温和——在原点附近可以展开且两个方向的耦合在一阶就出现了很适合看清楚交叉项的作用。取 $f(x,y) e^x\ln(1y)$展开点 $(0,0)$。先算函数值和一阶偏导$$f_0 e^0 \ln 1 0$$$$f_x e^x\ln(1y) \Rightarrow f_{x0} 1 \cdot 0 0$$$$f_y \frac{e^x}{1y} \Rightarrow f_{y0} \frac{1}{1} 1$$一阶项只有一个 $k$$h$ 前面的系数是 0。这在几何上说明沿 $x$ 方向的一阶变化率为零因为 $\ln(1y)$ 在 $y0$ 处是零$x$ 变大也乘不出来东西。接着算二阶偏导$$f_{xx} e^x\ln(1y) \Rightarrow f_{xx0} 0$$$$f_{xy} \frac{e^x}{1y} \Rightarrow f_{xy0} 1$$$$f_{yy} -\frac{e^x}{(1y)^2} \Rightarrow f_{yy0} -1$$代入二阶展开$$f \approx 0 0\cdot h 1\cdot k \frac{1}{2}\left(0\cdot h^2 2\cdot 1\cdot hk (-1)k^2\right) k hk - \frac{k^2}{2}$$写成原始变量就是 $y xy - y^2/2$。注意这里交叉项 $xy$ 的系数是 1因为 $2 \times 1/2 1$而 $y^2$ 项是 $-1/2$。这两个系数符号相反恰好反映了 $\ln(1y)$ 在正向是凹的、在负向增长更快的特性。3.2 例二$\sin x\cos y$ 在 $(\pi/3, \pi/4)$ 展到二阶第二题换个非原点的情况因为实际工作中很少正好在原点展开。取 $f(x,y)\sin x\cos y$展开点 $(x_0,y_0)(\pi/3,\pi/4)$。先把三角函数值备好$\sin\frac{\pi}{3}\frac{\sqrt{3}}{2}$$\cos\frac{\pi}{3}\frac{1}{2}$$\sin\frac{\pi}{4}\cos\frac{\pi}{4}\frac{\sqrt{2}}{2}$。函数值$f_0 \frac{\sqrt{3}}{2}\cdot\frac{\sqrt{2}}{2} \frac{\sqrt{6}}{4}$。一阶偏导$$f_x \cos x\cos y \Rightarrow f_{x0} \frac{1}{2}\cdot\frac{\sqrt{2}}{2} \frac{\sqrt{2}}{4}$$$$f_y -\sin x\sin y \Rightarrow f_{y0} -\frac{\sqrt{3}}{2}\cdot\frac{\sqrt{2}}{2} -\frac{\sqrt{6}}{4}$$二阶偏导$$f_{xx} -\sin x\cos y \Rightarrow f_{xx0} -\frac{\sqrt{3}}{2}\cdot\frac{\sqrt{2}}{2} -\frac{\sqrt{6}}{4}$$$$f_{xy} -\cos x\sin y \Rightarrow f_{xy0} -\frac{1}{2}\cdot\frac{\sqrt{2}}{2} -\frac{\sqrt{2}}{4}$$$$f_{yy} -\sin x\cos y \Rightarrow f_{yy0} -\frac{\sqrt{6}}{4}$$代入$$f \approx \frac{\sqrt{6}}{4} \frac{\sqrt{2}}{4}h - \frac{\sqrt{6}}{4}k - \frac{\sqrt{6}}{8}h^2 - \frac{\sqrt{2}}{4}hk - \frac{\sqrt{6}}{8}k^2$$可以看到这里 $f_{xx0}$ 和 $f_{yy0}$ 相等这是巧合还是必然对 $\sin x\cos y$ 来说这两个偏导在任意点都恒等都等于 $-\sin x\cos y$所以并不是展开点的特殊性质。相比之下交叉项系数是 $-\sqrt{2}/4$量级上比 $-\sqrt{6}/8 \approx -0.306$ 小一些$-\sqrt{2}/4\approx-0.354$两者接近说明在这个点附近两个方向的耦合与自身曲率同等重要不能忽略交叉项。3.3 用 sympy 和 numpy 复核系数手推的系数我一般不放心直接拿去用尤其是符号容易看串的时候。下面这两段代码是我常用来复核的模板。第一段用 sympy 做符号展开逐次对单个变量做级数展开效果等同于二元泰勒import sympy as sp x, y sp.symbols(x y) f sp.exp(x) * sp.log(1 y) # 先对 x 展开到 3 阶再对 y 展开到 3 阶 fx sp.series(f, x, 0, 3).removeO() fxy sp.series(fx, y, 0, 3).removeO() print(sp.expand(fxy))跑出来的多项式里常数项和一次项之后就是 $xy - y^2/2$与我们手算完全一致还会带上一些三阶及以上的残留项取决于阶数设置需要自己按阶截断。这里有个小坑sympy 的series对某个变量展开时会把另一个变量当作常数处理所以必须逐个变量依次展开而且两次的阶数都要给够否则高阶项会被误截断。我见过有人只对一个变量展开就直接用结果漏掉了 $x$ 方向的二阶项。第二段用 numpy 做数值验证直接比较近似多项式和真实值import numpy as np def f(x, y): return np.exp(x) * np.log(1 y) def p2(x, y): # 在原点展开的二阶近似y xy - y^2/2 return y x * y - y * y / 2 for h, k in [(0.01, 0.01), (0.05, 0.08), (0.1, 0.15), (0.2, 0.2)]: real f(h, k) approx p2(h, k) print(fh{h}, k{k}, real{real:.10f}, approx{approx:.10f}, err{real-approx:.3e})这段代码能直接告诉你近似在什么扰动幅度下还站得住。三类扰动下的误差量级分析、以及一个提升精度的技巧放在第 6 章细说因为那部分需要先有误差界的概念。4. 五个反复踩到的细节交叉项、余项、可交换性与展开点4.1 交叉项前的 2 什么时候会消失交叉项的系数在标准写法里是 $2f_{xy}$但如果你换一种写法把它写成 $\frac{1}{2}f_{xy}hk$ 加 $\frac{1}{2}f_{yx}hk$那个 2 就消失了——因为被拆成了两份。这是很多人抄公式时出错的根源不同教材的写法不同有的写 $\frac{1}{2}(f_{xx}h^2 2f_{xy}hk f_{yy}k^2)$有的写 $\frac{1}{2}f_{xx}h^2 f_{xy}hk \frac{1}{2}f_{yy}k^2$两种完全等价但如果你把第一种的括号拆开又忘了把 $\frac{1}{2}$ 继续乘到 $2f_{xy}hk$ 上就会得到 $f_{xy}hk$ 的两倍或一半。我自己的防错习惯是永远把二阶项写成括号形式括号外统一乘 $\frac{1}{2}$括号内交叉项写 $2f_{xy}hk$。这样三项对称检查时一眼就能看出是否漏了系数。写代码时也是同理用一个统一的二次型表达$$Q(\mathbf{v}) \frac{1}{2}\mathbf{v}^T H \mathbf{v}, \quad \mathbf{v} \begin{pmatrix} h \ k \end{pmatrix}, \quad H \begin{pmatrix} f_{xx} f_{xy} \ f_{xy} f_{yy} \end{pmatrix}$$用矩阵形式写交叉项系数 2 的问题被自动处理掉——因为 $H$ 是对称矩阵$\mathbf{v}^T H \mathbf{v}$ 展开时 $f_{xy}hk$ 会自动出现两次。这是我在代码里最推荐的写法几乎不可能出错。4.2 展开点不是原点时的写法与常见移项错误展开点不在原点时最常见的错误是把增量写成 $x$ 而不是 $h x - x_0$或者反过来把偏导数算在了错误的点上。举个具体的错法有人对 $f(x,y)x^2y$ 在 $(1,1)$ 处展开直接写 $f \approx 1 3h k \cdots$但他算 $f_x 2xy$ 时代入的是 $(0,0)$ 而不是 $(1,1)$得到的是 0 而不是 $2$整条式子就废了。我的做法是分三步走每一步都在纸上或代码里显式写清楚先算展开点处的所有偏导值列成一张小表表头是 $f, f_x, f_y, f_{xx}, f_{xy}, f_{yy}$值全部是常数。再定义增量 $h x-x_0$、$k y-y_0$把目标点用增量的形式写出来。最后把表里的常数代进公式得到一个关于 $h,k$ 的多项式。这三步分开做比一边推导一边代入要可靠得多。尤其第二步如果你要用这个近似去解方程或做迭代一定要在心里清楚公式是关于 $h,k$ 的多项式回代成 $x,y$ 的时候要把增量补回来别把两套坐标系混着用。4.3 阶数与余项要匹配二阶展开配三阶余项余项的阶数总比展开阶数高一阶这不是可选项。二阶展开配的是三阶余项 $R_2$三阶展开配的是四阶余项 $R_3$。这个规则来自一元泰勒二元情况下不变。为什么这件事值得单独拎出来讲因为做误差估计时如果你在二阶展开后面随手写了个误差是 $O(h^2)$那大概率是错的——除非展开点到目标点的路径上有三阶偏导为零的特殊结构。正确的说法是二阶展开的截断误差是 $O(|(h,k)|^3)$也就是三次量级。很多人被二阶近似这个词误导以为误差也是二阶结果误差界给宽了。量级验证很简单如果 $h,k$ 都是 $10^{-2}$ 量级那么 $|(h,k)|^3 \sim 10^{-6}$绝对误差在这个量级。如果你实测出来的误差是 $10^{-4}$说明有别的因素在起作用比如展开点选得不好、或者二阶项本身被算错了值得回头查一查。这个量级对照表我在调试时用得最多它能快速判断一个近似是否正常扰动幅度 $|(h,k)|$二阶展开预期误差量级一阶展开预期误差量级$10^{-3}$$10^{-9}$$10^{-6}$$10^{-2}$$10^{-6}$$10^{-4}$$10^{-1}$$10^{-3}$$10^{-2}$$1$很大不可用不可用4.4 混合偏导相等不是白给的公式里能把 $f_{xy}hk$ 和 $f_{yx}hk$ 合并成 $2f_{xy}hk$依赖的是 $f_{xy} f_{yx}$。这个结论克莱罗定理不是无条件成立的它要求两个混合偏导在展开点的一个邻域内存在且连续。绝大多数工程中遇到的函数初等函数、有理函数、常见的物理场都满足所以实践中默认成立没问题。但有两类情形需要留神。第一类是分段定义或者带绝对值的函数在分段边界上混合偏导可能不连续甚至不存在这时泰勒展开在边界附近本身就不适用谈系数合并没意义。第二类是数值计算中由离散数据插值出来的函数它的偏导数是我们自己用差分估的$f_{xy}$ 和 $f_{yx}$ 的数值估计会因为差分方向和步长不同而有差异这时强行合并会引入额外误差。我的建议是如果能解析求导就解析求导插值得到的函数就老老实实接受两个方向导数的估计值有差异别硬凑成 2 倍。还有一点实操提醒在用有限差分估算混合偏导时常用的九点公式是$$f_{xy} \approx \frac{f(xa, yb) - f(xa, y-b) - f(x-a, yb) f(x-a, y-b)}{4ab}$$步长 $a,b$ 的选取需要权衡。太大截断误差按 $O(a^2b^2)$ 增长太小浮点相减的舍入误差被放大。经验上取 $a \approx \epsilon^{1/3}\cdot \max(|x|,1)$ 比较稳其中 $\epsilon$ 是双精度机器精度约 $2.2\times10^{-16}$算下来 $a$ 大概在 $10^{-5}$ 量级。这个数字比很多人直觉上的 $10^{-8}$ 要大不少原因就是舍入误差的压制。4.5 变量的尺度差异会毁掉近似精度这条是我在真实项目里吃过亏之后才重视起来的。设想一个函数$x$ 是长度典型量级 $10^{-3}$ 米$y$ 是角度典型量级 $10^{-1}$ 弧度。两者的扰动幅度差了两个数量级。这时候二阶展开里$f_{xx}h^2/2$ 的量级是 $10^{-6}$而 $f_{yy}k^2/2$ 的量级是 $10^{-2}$差了四个数量级。如果你按统一的阶数去截断会出现一个问题为了让 $k$ 方向达到可接受的精度你可能需要展开到四阶但在 $x$ 方向上二阶早就绰绰有余了。处理办法有两个看你的场景选。如果要的是公式的简洁那就对两个方向分别截断$x$ 方向取二阶、$y$ 方向取四阶写出来的式子不对称但每一项都有用。如果要的是实现的统一那就先做无量纲化——把两个变量各自除以其典型尺度变成量级相当的归一化变量再统一展开到该阶数。后者在数值代码里更常用因为它避免了项数结构的分裂代价是需要维护一套归一化因子。这个细节在很多教材里不会强调因为纯数学讨论里变量没有量纲。但只要是做实际的工程计算量纲差异几乎是必然的稍微注意一下能省掉不少为什么我的近似在某些工况下突然不准的困惑。5. 工程里的实际用途误差传播、极值判定与牛顿迭代5.1 误差传播公式本质上就是二阶展开测量误差传递的经典公式是$$\sigma_f^2 \approx f_x^2\sigma_x^2 2f_xf_y\rho\sigma_x\sigma_y f_y^2\sigma_y^2$$很多人把它当作一条独立的统计公式来记其实它就是 $f$ 的一阶泰勒展开加上方差运算的直接结果。设 $f \approx f_0 f_x\Delta x f_y\Delta y$求方差不就得到上面这个式子$\rho$ 是两个输入量的相关系数。理解这层关系的好处是你能立刻知道这个公式的适用范围——它要求扰动足够小使得一阶项主导如果扰动较大、二阶项不可忽略这个公式就会系统性低估或高估方差。需要更精确的时候就要把二阶展开拿进来得到所谓的二阶误差传播。此时 $f$ 的表达式多出二次项方差的推导里会出现三阶和四阶矩偏度、峰度式子变得很长。我在实际项目中的折中做法是先算一阶传播再用一阶传播的残差去估二阶修正量如果修正量小于一阶结果的 5%就直接用一阶超过 5% 才去老老实实推二阶。这样避免了在不必要的地方把公式复杂化。顺带一提牛顿迭代法判断收敛二阶、梯度下降收敛一阶这些收敛阶的定义本身就来自泰勒展开的余项分析。梯度下降的每步残差按 $O(\alpha)$ 递减$\alpha$ 是步长而牛顿法按 $O(\alpha^2)$ 递减原因就在于牛顿法用到了二阶信息把一次项的残差消掉了。这些看似独立的数值算法结论源头都在同一个公式上。5.2 Hessian 正定性与极值判定二阶展开在优化里最重要的一处应用是极值判定。把展开点取在驻点$f_xf_y0$一阶项消失于是$$f(x_0h,\ y_0k) - f_0 \approx \frac{1}{2}\left(f_{xx}h^2 2f_{xy}hk f_{yy}k^2\right) \frac{1}{2}\mathbf{v}^T H\mathbf{v}$$右端是一个二次型。二次型的符号完全由矩阵 $H$ 的定性决定若 $H$ 正定所有特征值为正二次型对任意非零 $\mathbf{v}$ 都为正说明函数在这一点取极小值。若 $H$ 负定所有特征值为负取极大值。若 $H$ 不定特征值有正有负该点是鞍点既不是极大也不是极小。对二元情况正定性可以用行列式判据快速判断$f_{xx} 0$ 且 $\det H f_{xx}f_{yy} - f_{xy}^2 0$ 时正定$f_{xx} 0$ 且 $\det H 0$ 时负定$\det H 0$ 时不定。只有 $\det H 0$ 的退化情形才需要看更高阶这种情况在实际数据里基本遇不到除非函数有特殊结构。我把这套判据用在曲面参数拟合的稳定性检查上有过一次很有价值的收获。当时拟合出一个驻点只看一阶导确实为零但算完 Hessian 发现 $\det H 0$是个鞍点。如果按极小值去处理参数会沿着负曲率方向跑飞。这类问题在做最小二乘时尤其要留神——迭代算法给出的驻点不一定是我们要的极小值点用 Hessian 验一下只要几行代码。5.3 二维牛顿法二阶展开的直接产物牛顿法求方程根本质上就是拿一阶展开去近似函数然后令其为零。多元情形同理。设要求解方程组 $\mathbf{F}(\mathbf{x}) \mathbf{0}$在当前点 $\mathbf{x}_n$ 处做一阶展开$$\mathbf{F}(\mathbf{x}_n \Delta\mathbf{x}) \approx \mathbf{F}(\mathbf{x}_n) J(\mathbf{x}_n)\Delta\mathbf{x}$$令右端为零解出增量$$\Delta\mathbf{x} -J^{-1}\mathbf{F}(\mathbf{x}n), \quad \mathbf{x}{n1} \mathbf{x}_n \Delta\mathbf{x}$$这就是多元牛顿迭代。对优化问题目标是求 $\nabla f 0$把它代入上式$J$ 就成了 Hessian $H$于是得到$$\mathbf{x}_{n1} \mathbf{x}_n - H^{-1}\nabla f(\mathbf{x}_n)$$需要注意的实操细节是牛顿法在用泰勒展开时要求 Hessian 正定才能保证沿下降方向走否则更新方向可能朝上去。工程上常见的是加上线搜索或者阻尼因子其实就是信赖域方法的简化版。另外 Hessian 的解析表达式不一定好求这时候要么用第 4.4 节的差分公式估要么用拟牛顿法BFGS 之类迭代近似 —— 后者本质上是在用一系列一阶信息拼出一个越来越准的二阶模型思路仍然是泰勒展开只是把二阶导的获取方式换了。6. 数值实验把近似多项式和真实值摆在一起看6.1 实验代码理论说了一堆最后还是要看数字。用第 3 章那个 $f(x,y)e^x\ln(1y)$ 做一次完整的数值对照把一阶、二阶近似和真实值三列摆在一起同时给出误差比值看它和理论量级是否吻合。import numpy as np def f(x, y): return np.exp(x) * np.log(1 y) def p1(x, y): # 一阶近似y return y def p2(x, y): # 二阶近似y xy - y^2/2 return y x * y - y * y / 2 def p3(x, y): # 三阶近似用于看误差是否按预期下降 return y x * y - y * y / 2 (x * x * y) / 2 - (x * y * y) / 2 (y ** 3) / 3 cases [(0.01, 0.01), (0.05, 0.08), (0.1, 0.15), (0.2, 0.2)] print(f{h:6} {k:6} {real:14} {p1:14} {p2:14} {err1:10} {err2:10}) for h, k in cases: r f(h, k) a1 p1(h, k) a2 p2(h, k) print(f{h:6.2f} {k:6.2f} {r:14.10f} {a1:14.10f} {a2:14.10f} f{abs(r-a1):10.2e} {abs(r-a2):10.2e})三阶近似里我手工补了展开到三次的项用于验证误差下降的规律$x^2y/2$ 来自 $e^x$ 的二次项乘 $\ln(1y)$ 的一次项$-xy^2/2$ 来自两个一次项相乘$y^3/3$ 来自 $\ln(1y)$ 的三次项。6.2 结果与截断误差量级分析按上面的参数跑会看到几个很有说明性的现象。以 $(h,k)(0.05,0.08)$ 这组为例真实值大约是 $0.08090692$一阶近似是 $0.08$误差约 $9.07\times10^{-4}$二阶近似是 $0.0808$误差约 $1.07\times10^{-4}$。误差从 $10^{-3}$ 掉到 $10^{-4}$正好差一个量级和一阶误差 $O(|v|^2)$、二阶误差 $O(|v|^3)$的预测一致这里 $|v|\sim0.1$平方是 $10^{-2}$立方是 $10^{-3}$再乘上偏导系数的量级落到 $10^{-3}$ 和 $10^{-4}$ 很合理。当扰动扩大到 $(0.2,0.2)$ 时变化就明显了。$|v|\sim 0.28$二阶误差理论上应该落在 $10^{-2}$ 量级实测也确实在这个范围。这时候一阶近似已经开始不能用了相对误差接近 15%。这个数字给我的直观感受是二阶展开的舒适区大约在扰动小于展开点尺度 10% 的时候超过 20%就该考虑换展开点或者提高阶数了。另外有一个容易被忽略的现象误差的符号。在上述几组数据里误差始终是正的真实值大于近似值说明二阶近似在这个区域内是系统性低估。这不是随机误差不会因为多测几次而抵消。如果拿这个近似去做批量计算的中间步骤这个系统偏差会累积起来最终结果的偏差方向是可预测的但幅度可能被放大。修正办法很简单就是再补上三阶项或者把展开点挪到更靠近实际工作区间的位置。6.3 一个改善精度的实用技巧最后分享一个我在做参数化模型时用得比较顺手的技巧不要在固定点上展开而是在工作区间的中点展开。原因不难理解。截断误差和展开点到目标点的距离的三次方成正比。如果你把展开点选在区间一端那么另一端的误差就最大整个区间的最大误差由这个最远点决定。移到中点后最远距离减半误差按三次方算就降到原来的八分之一。这是零成本的改进——偏导数的表达式不用改只需要换一下取值点。实测中这条经验的效果比想象中明显。同样是一段 $[0, 0.2]$ 的参数区间端点展开时最大误差是 $8\times10^{-3}$ 量级中点展开降到 $1\times10^{-3}$ 量级正好八倍关系。所以现在我做任何局部近似第一件事就是把展开点定在区间的几何中心然后把整个区间两端都验一遍误差取最大值作为这个近似模型的精度指标写进注释里。这样后来接手的人一眼就能知道这个近似在什么范围内可信比写一句误差较小有用得多。如果要追求更高的精度但不想增加阶数还有一个办法是把大区间切成若干小区间每段各自在中点展开然后在边界处平滑衔接。代价是系数表变长好处是每一段里都是三次量级的小误差。这套思路其实就是典型的查表法配合局部多项式在嵌入式设备上做超越函数计算时非常常见我后来在做实时估算模块时也是这么处理的。
二元函数泰勒展开全解析：交叉项系数2与工程应用
做数值计算的人大概都有过这么一次经历手头有个函数 $f(x,y)$ 在某点附近变化很小想用低阶多项式替它做快速估算于是顺手把一元泰勒公式的套路搬过来写成 $f \approx f_0 f_x h f_y k$觉得很合理。可一旦你把展开阶数提高到二阶问题就来了——交叉项前面那个系数到底是 1 还是 2符号表里查不到教材上一晃而过自己推又容易在两个变量的求导顺序里绕晕。更麻烦的是这个系数错了不一定报错程序照跑结果照出只是误差估计悄悄偏了一倍等到某个模型开始偶尔不准的时候你已经很难回溯到这一步了。二元函数的泰勒展开就是干这件事的把两个变量的函数在某个点的邻域内用一个多项式去逼近它逼近的好坏由阶数和扰动幅度决定。它是误差传播分析、Hessian 极值判定、牛顿迭代、灵敏度分析、数值线性化这些工作的共同底座。下面我按自己教学和工程中反复用到的顺序把这条公式从为什么长这样一直讲到怎么算才不会错中间穿插两道手算例题和一次数值复核。不管你是正在啃高数的学生还是需要写误差分析代码的工程师看完应该都能直接上手。1. 先别急着背公式把二元展开压回一元的那条直线上1.1 一个把误差估计算错一倍的现场我最早踩的坑发生在做传感器标定的误差传递分析时。当时有两个输入量 $x$ 和 $y$各自带有小幅随机扰动我想估计输出量的偏差就写了一阶近似 $f_x \Delta x f_y \Delta y$这一层没问题。但当我想进一步估计二阶偏差时脑子里默认写了 $\frac{1}{2}(f_{xx}\Delta x^2 f_{xy}\Delta x\Delta y f_{yy}\Delta y^2)$把交叉项系数漏掉了。算出来的偏差比实测小了一截我一度怀疑是测量噪声的问题换了几组数据都对不上。后来把公式老实展开一遍才发现交叉项前面确实是 2不是 1。这不是排版笔误而是两个求导路径叠加的结果。这个错误的隐蔽性在于它不影响一阶项只影响二阶及以上如果你的模型里一阶项占主导你根本发现不了但只要扰动稍微大一点二阶项开始起作用误差就浮出来了。所以我现在的习惯是凡是涉及二元的二阶近似交叉项系数必须现场推一遍不靠记忆。1.2 沿直线取样把二元的难题降成一元理解二元展开最省脑子的办法不是去背那个长长的求和式而是把它降维回一元。具体做法是固定一个方向从展开点 $(x_0, y_0)$ 出发沿着向量 $(h, k)$ 走一条直线把点写成参数形式$$x x_0 th, \quad y y_0 tk, \quad t \in [0, 1]$$然后把 $f$ 沿着这条直线限制成一个单变量函数$$g(t) f(x_0 th,\ y_0 tk)$$这个操作的关键在于$g(t)$ 是一个货真价实的一元函数$t0$ 对应展开点$t1$ 对应我们要估算的目标点。于是二元展开的全部困难就转化成了对 $g(t)$ 在 $t0$ 处做一元泰勒展开然后取 $t1$这么一件事。用一句大白话说二元泰勒展开本质上就是沿各个方向做一元展开之后的结果被统一写成了一个式子。这个视角还有一个额外好处——它顺便解释了为什么二元展开的余项形式和一元几乎一模一样。因为 $g(t)$ 的余项就是一元余项只是里面的导数换成了沿方向 $(h,k)$ 的方向导数而已。你不需要为二元的余项单独发明一套理论。1.3 g(t) 求导交叉项的系数 2 是这么冒出来的现在把 $g(t)$ 求导让系数 2 自己现身。根据链式法则$$g(t) f_x(x_0th,\ y_0tk)\cdot h f_y(x_0th,\ y_0tk)\cdot k$$这一步是干净的一阶项就是 $f_x h f_y k$没有争议。麻烦在二阶。再对 $t$ 求一次导注意 $f_x$ 和 $f_y$ 本身也是二元函数还要继续用链式法则$$\frac{d}{dt}f_x(\cdot) f_{xx}h f_{xy}k$$$$\frac{d}{dt}f_y(\cdot) f_{yx}h f_{yy}k$$把这两个结果分别乘上外面的 $h$ 和 $k$$$g(t) (f_{xx}h f_{xy}k)h (f_{yx}h f_{yy}k)k f_{xx}h^2 f_{xy}hk f_{yx}hk f_{yy}k^2$$到这里答案已经很清楚了$hk$ 这一项出现了两次一次来自对 $f_x$ 求导、一次来自对 $f_y$ 求导。如果函数是二阶连续可微的那么 $f_{xy} f_{yx}$两项合并就是 $2f_{xy}hk$。这就是那个 2 的全部来历——它不是人为规定的系数而是两条求导路径的合并结果。我一直觉得只要你在纸上把这一段亲手写一遍这个 2 就再也忘不掉了。反过来如果只背结论很容易在写代码或者列公式的时候顺手把它省掉因为它看起来应该像一元那样是 1/2 系的整齐结构。数学的整齐有时候是假象这里是真有两份。2. 完整公式长什么样从一阶切平面到 n 阶算子写法2.1 记号约定h、k 是增量不是微分先把记号钉死很多混乱其实来自记号。设 $h x - x_0$$k y - y_0$它们是把目标点 $(x,y)$ 拉回到展开点 $(x_0,y_0)$ 的有限增量不是无穷小量。这一点很关键泰勒展开之所以有实用价值恰恰是因为它处理的是有限的、不太小的扰动而不是微分那种极限意义下的无穷小。展开点通常记成下标 0比如 $f(x_0,y_0)$ 写成 $f_0$$f_x(x_0,y_0)$ 写成 $f_{x0}$依此类推。不带下标的偏导数符号表示在展开点取值。写代码的时候我建议把展开点坐标和增量分成两组变量管理比如x0, y0和dx, dy避免在推导过程中把 $h$ 误当成 $x$ 的一部分去求导那是另一个高频错误。另外提醒一句$h$ 和 $k$ 是有量纲的它们的量纲和被展开的函数自变量一致。这一点在第 4 章会展开讲因为当两个自变量的物理量纲差别很大时比如一个是长度、一个是角度展开阶数的选择会变得很微妙。2.2 一阶展开与切平面近似一阶展开最直观公式是$$f(x_0h,\ y_0k) \approx f_0 f_{x0}h f_{y0}k$$几何意义特别干净右边是关于 $h,k$ 的线性函数在三维空间里画出来就是一个平面而且这个平面在点 $(x_0,y_0,f_0)$ 处与曲面相切所以叫切平面近似。曲面在这一点附近有多像这个平面取决于曲面弯曲得有多厉害。如果曲面在该点附近几乎是平的一阶就够如果弯得明显一阶会系统性偏低或偏高因为有二阶项在单边贡献。工程上用一阶展开的场景非常多因为它天然就是线性化。控制系统里把非线性环节在工作点附近线性化本质就是一阶泰勒展开电路的小信号模型也是一阶展开测量里的误差传递一阶公式还是它。一阶展开的好处不只是简单而是它让叠加原理成立——多个小扰动的作用可以直接相加。但一阶展开有个系统性的短板它给出的近似值与真实值之间的偏差符号是有规律的不随机。也就是说一阶近似在展开点的一侧偏高、另一侧偏低如果你做的是一组重复测量然后取平均这个偏差不会自动抵消会作为系统偏差留在结果里。要处理它就得上二阶。2.3 二阶展开的显式形式和余项把 $g(t)$ 和 $g(t)$ 在 $t0$ 处取值代入一元泰勒并令 $t1$得到二元二阶展开$$f(x_0h,\ y_0k) f_0 f_{x0}h f_{y0}k \frac{1}{2}\left(f_{xx0}h^2 2f_{xy0}hk f_{yy0}k^2\right) R_2$$把这三项的角色分清楚很有必要。$f_{xx0}h^2/2$ 是沿 $x$ 方向的纯曲率贡献$f_{yy0}k^2/2$ 是沿 $y$ 方向的纯曲率贡献而 $f_{xy0}hk$ 是交叉贡献它描述的是$x$ 方向变化会不会改变 $y$ 方向的斜率也就是两个变量之间的耦合程度。如果 $f_{xy0}0$说明在该点附近两个方向可以解耦处理如果 $f_{xy0}$ 很大说明耦合强烈你不能分别调 $x$ 和 $y$。余项 $R_2$ 用拉格朗日形式写是$$R_2 \frac{1}{6}\Big(h\frac{\partial}{\partial x} k\frac{\partial}{\partial y}\Big)^3 f(x_0\theta h,\ y_0\theta k), \quad 0\theta1$$注意展开成六个三阶偏导的线性组合后前系数是 $1/6$这个 $1/6$ 来自 $3!$。余项的意义在于给出误差的上界实际做误差估计时我们通常不知道 $\theta$但可以用区域内三阶偏导的最大值 $M_3$ 去卡一个保守界$$|R_2| \le \frac{M_3}{6}\left(|h| |k|\right)^3$$这个不等式在实际工作里非常好用因为它只需要估一个三阶导数的上界不需要知道具体的 $\theta$ 位置。2.4 n 阶通式算子幂写法要把结论推广到任意阶最省事的写法是引入方向导数算子$$D h\frac{\partial}{\partial x} k\frac{\partial}{\partial y}$$那么 $n$ 阶展开就是$$f(x_0h,\ y_0k) \sum_{m0}^{n}\frac{1}{m!}D^m f(x_0,y_0) R_n$$算子写法看着抽象但它的好处是把公式长什么样的问题一次性解决了。你把 $D^m$ 当作二项式展开处理就行$D^2$ 展开出 $h^2\partial_{xx} 2hk\partial_{xy} k^2\partial_{yy}$那两倍的系数自然来自二项式系数 $\binom{2}{1}2$$D^3$ 展开出 $h^3\partial_{xxx} 3h^2k\partial_{xxy} 3hk^2\partial_{xyy} k^3\partial_{yyy}$系数 3 来自 $\binom{3}{1}\binom{3}{2}3$。所以之前那个神秘的 2在算子视角下就是二项式系数。这个统一性很有价值三元、四元函数的展开公式不需要另记你只要把算子写成 $D h_1\partial_1 h_2\partial_2 h_3\partial_3$然后做多项式展开系数由多项式定理给出。我平时写多变量近似代码时就是按多项式系数批量生成各项而不是手工逐项列。顺带说一个实际取舍工程计算里用到四阶以上展开的情形很少因为高维情况下项数增长太快。二元函数到二阶有 6 个系数含常数、两个一阶、三个二阶到三阶再加 4 个三阶偏导一共 10 个到四阶再加 5 个是 15 个。三元的话二阶就有 10 项。所以实践中二元通常停在二阶或三阶需要更高精度时改换别的数值手段而不是硬堆阶数。3. 两道手算例题从原点展开到非原点展开3.1 例一$e^x \ln(1y)$ 在原点展到二阶选这个函数是因为它足够温和——在原点附近可以展开且两个方向的耦合在一阶就出现了很适合看清楚交叉项的作用。取 $f(x,y) e^x\ln(1y)$展开点 $(0,0)$。先算函数值和一阶偏导$$f_0 e^0 \ln 1 0$$$$f_x e^x\ln(1y) \Rightarrow f_{x0} 1 \cdot 0 0$$$$f_y \frac{e^x}{1y} \Rightarrow f_{y0} \frac{1}{1} 1$$一阶项只有一个 $k$$h$ 前面的系数是 0。这在几何上说明沿 $x$ 方向的一阶变化率为零因为 $\ln(1y)$ 在 $y0$ 处是零$x$ 变大也乘不出来东西。接着算二阶偏导$$f_{xx} e^x\ln(1y) \Rightarrow f_{xx0} 0$$$$f_{xy} \frac{e^x}{1y} \Rightarrow f_{xy0} 1$$$$f_{yy} -\frac{e^x}{(1y)^2} \Rightarrow f_{yy0} -1$$代入二阶展开$$f \approx 0 0\cdot h 1\cdot k \frac{1}{2}\left(0\cdot h^2 2\cdot 1\cdot hk (-1)k^2\right) k hk - \frac{k^2}{2}$$写成原始变量就是 $y xy - y^2/2$。注意这里交叉项 $xy$ 的系数是 1因为 $2 \times 1/2 1$而 $y^2$ 项是 $-1/2$。这两个系数符号相反恰好反映了 $\ln(1y)$ 在正向是凹的、在负向增长更快的特性。3.2 例二$\sin x\cos y$ 在 $(\pi/3, \pi/4)$ 展到二阶第二题换个非原点的情况因为实际工作中很少正好在原点展开。取 $f(x,y)\sin x\cos y$展开点 $(x_0,y_0)(\pi/3,\pi/4)$。先把三角函数值备好$\sin\frac{\pi}{3}\frac{\sqrt{3}}{2}$$\cos\frac{\pi}{3}\frac{1}{2}$$\sin\frac{\pi}{4}\cos\frac{\pi}{4}\frac{\sqrt{2}}{2}$。函数值$f_0 \frac{\sqrt{3}}{2}\cdot\frac{\sqrt{2}}{2} \frac{\sqrt{6}}{4}$。一阶偏导$$f_x \cos x\cos y \Rightarrow f_{x0} \frac{1}{2}\cdot\frac{\sqrt{2}}{2} \frac{\sqrt{2}}{4}$$$$f_y -\sin x\sin y \Rightarrow f_{y0} -\frac{\sqrt{3}}{2}\cdot\frac{\sqrt{2}}{2} -\frac{\sqrt{6}}{4}$$二阶偏导$$f_{xx} -\sin x\cos y \Rightarrow f_{xx0} -\frac{\sqrt{3}}{2}\cdot\frac{\sqrt{2}}{2} -\frac{\sqrt{6}}{4}$$$$f_{xy} -\cos x\sin y \Rightarrow f_{xy0} -\frac{1}{2}\cdot\frac{\sqrt{2}}{2} -\frac{\sqrt{2}}{4}$$$$f_{yy} -\sin x\cos y \Rightarrow f_{yy0} -\frac{\sqrt{6}}{4}$$代入$$f \approx \frac{\sqrt{6}}{4} \frac{\sqrt{2}}{4}h - \frac{\sqrt{6}}{4}k - \frac{\sqrt{6}}{8}h^2 - \frac{\sqrt{2}}{4}hk - \frac{\sqrt{6}}{8}k^2$$可以看到这里 $f_{xx0}$ 和 $f_{yy0}$ 相等这是巧合还是必然对 $\sin x\cos y$ 来说这两个偏导在任意点都恒等都等于 $-\sin x\cos y$所以并不是展开点的特殊性质。相比之下交叉项系数是 $-\sqrt{2}/4$量级上比 $-\sqrt{6}/8 \approx -0.306$ 小一些$-\sqrt{2}/4\approx-0.354$两者接近说明在这个点附近两个方向的耦合与自身曲率同等重要不能忽略交叉项。3.3 用 sympy 和 numpy 复核系数手推的系数我一般不放心直接拿去用尤其是符号容易看串的时候。下面这两段代码是我常用来复核的模板。第一段用 sympy 做符号展开逐次对单个变量做级数展开效果等同于二元泰勒import sympy as sp x, y sp.symbols(x y) f sp.exp(x) * sp.log(1 y) # 先对 x 展开到 3 阶再对 y 展开到 3 阶 fx sp.series(f, x, 0, 3).removeO() fxy sp.series(fx, y, 0, 3).removeO() print(sp.expand(fxy))跑出来的多项式里常数项和一次项之后就是 $xy - y^2/2$与我们手算完全一致还会带上一些三阶及以上的残留项取决于阶数设置需要自己按阶截断。这里有个小坑sympy 的series对某个变量展开时会把另一个变量当作常数处理所以必须逐个变量依次展开而且两次的阶数都要给够否则高阶项会被误截断。我见过有人只对一个变量展开就直接用结果漏掉了 $x$ 方向的二阶项。第二段用 numpy 做数值验证直接比较近似多项式和真实值import numpy as np def f(x, y): return np.exp(x) * np.log(1 y) def p2(x, y): # 在原点展开的二阶近似y xy - y^2/2 return y x * y - y * y / 2 for h, k in [(0.01, 0.01), (0.05, 0.08), (0.1, 0.15), (0.2, 0.2)]: real f(h, k) approx p2(h, k) print(fh{h}, k{k}, real{real:.10f}, approx{approx:.10f}, err{real-approx:.3e})这段代码能直接告诉你近似在什么扰动幅度下还站得住。三类扰动下的误差量级分析、以及一个提升精度的技巧放在第 6 章细说因为那部分需要先有误差界的概念。4. 五个反复踩到的细节交叉项、余项、可交换性与展开点4.1 交叉项前的 2 什么时候会消失交叉项的系数在标准写法里是 $2f_{xy}$但如果你换一种写法把它写成 $\frac{1}{2}f_{xy}hk$ 加 $\frac{1}{2}f_{yx}hk$那个 2 就消失了——因为被拆成了两份。这是很多人抄公式时出错的根源不同教材的写法不同有的写 $\frac{1}{2}(f_{xx}h^2 2f_{xy}hk f_{yy}k^2)$有的写 $\frac{1}{2}f_{xx}h^2 f_{xy}hk \frac{1}{2}f_{yy}k^2$两种完全等价但如果你把第一种的括号拆开又忘了把 $\frac{1}{2}$ 继续乘到 $2f_{xy}hk$ 上就会得到 $f_{xy}hk$ 的两倍或一半。我自己的防错习惯是永远把二阶项写成括号形式括号外统一乘 $\frac{1}{2}$括号内交叉项写 $2f_{xy}hk$。这样三项对称检查时一眼就能看出是否漏了系数。写代码时也是同理用一个统一的二次型表达$$Q(\mathbf{v}) \frac{1}{2}\mathbf{v}^T H \mathbf{v}, \quad \mathbf{v} \begin{pmatrix} h \ k \end{pmatrix}, \quad H \begin{pmatrix} f_{xx} f_{xy} \ f_{xy} f_{yy} \end{pmatrix}$$用矩阵形式写交叉项系数 2 的问题被自动处理掉——因为 $H$ 是对称矩阵$\mathbf{v}^T H \mathbf{v}$ 展开时 $f_{xy}hk$ 会自动出现两次。这是我在代码里最推荐的写法几乎不可能出错。4.2 展开点不是原点时的写法与常见移项错误展开点不在原点时最常见的错误是把增量写成 $x$ 而不是 $h x - x_0$或者反过来把偏导数算在了错误的点上。举个具体的错法有人对 $f(x,y)x^2y$ 在 $(1,1)$ 处展开直接写 $f \approx 1 3h k \cdots$但他算 $f_x 2xy$ 时代入的是 $(0,0)$ 而不是 $(1,1)$得到的是 0 而不是 $2$整条式子就废了。我的做法是分三步走每一步都在纸上或代码里显式写清楚先算展开点处的所有偏导值列成一张小表表头是 $f, f_x, f_y, f_{xx}, f_{xy}, f_{yy}$值全部是常数。再定义增量 $h x-x_0$、$k y-y_0$把目标点用增量的形式写出来。最后把表里的常数代进公式得到一个关于 $h,k$ 的多项式。这三步分开做比一边推导一边代入要可靠得多。尤其第二步如果你要用这个近似去解方程或做迭代一定要在心里清楚公式是关于 $h,k$ 的多项式回代成 $x,y$ 的时候要把增量补回来别把两套坐标系混着用。4.3 阶数与余项要匹配二阶展开配三阶余项余项的阶数总比展开阶数高一阶这不是可选项。二阶展开配的是三阶余项 $R_2$三阶展开配的是四阶余项 $R_3$。这个规则来自一元泰勒二元情况下不变。为什么这件事值得单独拎出来讲因为做误差估计时如果你在二阶展开后面随手写了个误差是 $O(h^2)$那大概率是错的——除非展开点到目标点的路径上有三阶偏导为零的特殊结构。正确的说法是二阶展开的截断误差是 $O(|(h,k)|^3)$也就是三次量级。很多人被二阶近似这个词误导以为误差也是二阶结果误差界给宽了。量级验证很简单如果 $h,k$ 都是 $10^{-2}$ 量级那么 $|(h,k)|^3 \sim 10^{-6}$绝对误差在这个量级。如果你实测出来的误差是 $10^{-4}$说明有别的因素在起作用比如展开点选得不好、或者二阶项本身被算错了值得回头查一查。这个量级对照表我在调试时用得最多它能快速判断一个近似是否正常扰动幅度 $|(h,k)|$二阶展开预期误差量级一阶展开预期误差量级$10^{-3}$$10^{-9}$$10^{-6}$$10^{-2}$$10^{-6}$$10^{-4}$$10^{-1}$$10^{-3}$$10^{-2}$$1$很大不可用不可用4.4 混合偏导相等不是白给的公式里能把 $f_{xy}hk$ 和 $f_{yx}hk$ 合并成 $2f_{xy}hk$依赖的是 $f_{xy} f_{yx}$。这个结论克莱罗定理不是无条件成立的它要求两个混合偏导在展开点的一个邻域内存在且连续。绝大多数工程中遇到的函数初等函数、有理函数、常见的物理场都满足所以实践中默认成立没问题。但有两类情形需要留神。第一类是分段定义或者带绝对值的函数在分段边界上混合偏导可能不连续甚至不存在这时泰勒展开在边界附近本身就不适用谈系数合并没意义。第二类是数值计算中由离散数据插值出来的函数它的偏导数是我们自己用差分估的$f_{xy}$ 和 $f_{yx}$ 的数值估计会因为差分方向和步长不同而有差异这时强行合并会引入额外误差。我的建议是如果能解析求导就解析求导插值得到的函数就老老实实接受两个方向导数的估计值有差异别硬凑成 2 倍。还有一点实操提醒在用有限差分估算混合偏导时常用的九点公式是$$f_{xy} \approx \frac{f(xa, yb) - f(xa, y-b) - f(x-a, yb) f(x-a, y-b)}{4ab}$$步长 $a,b$ 的选取需要权衡。太大截断误差按 $O(a^2b^2)$ 增长太小浮点相减的舍入误差被放大。经验上取 $a \approx \epsilon^{1/3}\cdot \max(|x|,1)$ 比较稳其中 $\epsilon$ 是双精度机器精度约 $2.2\times10^{-16}$算下来 $a$ 大概在 $10^{-5}$ 量级。这个数字比很多人直觉上的 $10^{-8}$ 要大不少原因就是舍入误差的压制。4.5 变量的尺度差异会毁掉近似精度这条是我在真实项目里吃过亏之后才重视起来的。设想一个函数$x$ 是长度典型量级 $10^{-3}$ 米$y$ 是角度典型量级 $10^{-1}$ 弧度。两者的扰动幅度差了两个数量级。这时候二阶展开里$f_{xx}h^2/2$ 的量级是 $10^{-6}$而 $f_{yy}k^2/2$ 的量级是 $10^{-2}$差了四个数量级。如果你按统一的阶数去截断会出现一个问题为了让 $k$ 方向达到可接受的精度你可能需要展开到四阶但在 $x$ 方向上二阶早就绰绰有余了。处理办法有两个看你的场景选。如果要的是公式的简洁那就对两个方向分别截断$x$ 方向取二阶、$y$ 方向取四阶写出来的式子不对称但每一项都有用。如果要的是实现的统一那就先做无量纲化——把两个变量各自除以其典型尺度变成量级相当的归一化变量再统一展开到该阶数。后者在数值代码里更常用因为它避免了项数结构的分裂代价是需要维护一套归一化因子。这个细节在很多教材里不会强调因为纯数学讨论里变量没有量纲。但只要是做实际的工程计算量纲差异几乎是必然的稍微注意一下能省掉不少为什么我的近似在某些工况下突然不准的困惑。5. 工程里的实际用途误差传播、极值判定与牛顿迭代5.1 误差传播公式本质上就是二阶展开测量误差传递的经典公式是$$\sigma_f^2 \approx f_x^2\sigma_x^2 2f_xf_y\rho\sigma_x\sigma_y f_y^2\sigma_y^2$$很多人把它当作一条独立的统计公式来记其实它就是 $f$ 的一阶泰勒展开加上方差运算的直接结果。设 $f \approx f_0 f_x\Delta x f_y\Delta y$求方差不就得到上面这个式子$\rho$ 是两个输入量的相关系数。理解这层关系的好处是你能立刻知道这个公式的适用范围——它要求扰动足够小使得一阶项主导如果扰动较大、二阶项不可忽略这个公式就会系统性低估或高估方差。需要更精确的时候就要把二阶展开拿进来得到所谓的二阶误差传播。此时 $f$ 的表达式多出二次项方差的推导里会出现三阶和四阶矩偏度、峰度式子变得很长。我在实际项目中的折中做法是先算一阶传播再用一阶传播的残差去估二阶修正量如果修正量小于一阶结果的 5%就直接用一阶超过 5% 才去老老实实推二阶。这样避免了在不必要的地方把公式复杂化。顺带一提牛顿迭代法判断收敛二阶、梯度下降收敛一阶这些收敛阶的定义本身就来自泰勒展开的余项分析。梯度下降的每步残差按 $O(\alpha)$ 递减$\alpha$ 是步长而牛顿法按 $O(\alpha^2)$ 递减原因就在于牛顿法用到了二阶信息把一次项的残差消掉了。这些看似独立的数值算法结论源头都在同一个公式上。5.2 Hessian 正定性与极值判定二阶展开在优化里最重要的一处应用是极值判定。把展开点取在驻点$f_xf_y0$一阶项消失于是$$f(x_0h,\ y_0k) - f_0 \approx \frac{1}{2}\left(f_{xx}h^2 2f_{xy}hk f_{yy}k^2\right) \frac{1}{2}\mathbf{v}^T H\mathbf{v}$$右端是一个二次型。二次型的符号完全由矩阵 $H$ 的定性决定若 $H$ 正定所有特征值为正二次型对任意非零 $\mathbf{v}$ 都为正说明函数在这一点取极小值。若 $H$ 负定所有特征值为负取极大值。若 $H$ 不定特征值有正有负该点是鞍点既不是极大也不是极小。对二元情况正定性可以用行列式判据快速判断$f_{xx} 0$ 且 $\det H f_{xx}f_{yy} - f_{xy}^2 0$ 时正定$f_{xx} 0$ 且 $\det H 0$ 时负定$\det H 0$ 时不定。只有 $\det H 0$ 的退化情形才需要看更高阶这种情况在实际数据里基本遇不到除非函数有特殊结构。我把这套判据用在曲面参数拟合的稳定性检查上有过一次很有价值的收获。当时拟合出一个驻点只看一阶导确实为零但算完 Hessian 发现 $\det H 0$是个鞍点。如果按极小值去处理参数会沿着负曲率方向跑飞。这类问题在做最小二乘时尤其要留神——迭代算法给出的驻点不一定是我们要的极小值点用 Hessian 验一下只要几行代码。5.3 二维牛顿法二阶展开的直接产物牛顿法求方程根本质上就是拿一阶展开去近似函数然后令其为零。多元情形同理。设要求解方程组 $\mathbf{F}(\mathbf{x}) \mathbf{0}$在当前点 $\mathbf{x}_n$ 处做一阶展开$$\mathbf{F}(\mathbf{x}_n \Delta\mathbf{x}) \approx \mathbf{F}(\mathbf{x}_n) J(\mathbf{x}_n)\Delta\mathbf{x}$$令右端为零解出增量$$\Delta\mathbf{x} -J^{-1}\mathbf{F}(\mathbf{x}n), \quad \mathbf{x}{n1} \mathbf{x}_n \Delta\mathbf{x}$$这就是多元牛顿迭代。对优化问题目标是求 $\nabla f 0$把它代入上式$J$ 就成了 Hessian $H$于是得到$$\mathbf{x}_{n1} \mathbf{x}_n - H^{-1}\nabla f(\mathbf{x}_n)$$需要注意的实操细节是牛顿法在用泰勒展开时要求 Hessian 正定才能保证沿下降方向走否则更新方向可能朝上去。工程上常见的是加上线搜索或者阻尼因子其实就是信赖域方法的简化版。另外 Hessian 的解析表达式不一定好求这时候要么用第 4.4 节的差分公式估要么用拟牛顿法BFGS 之类迭代近似 —— 后者本质上是在用一系列一阶信息拼出一个越来越准的二阶模型思路仍然是泰勒展开只是把二阶导的获取方式换了。6. 数值实验把近似多项式和真实值摆在一起看6.1 实验代码理论说了一堆最后还是要看数字。用第 3 章那个 $f(x,y)e^x\ln(1y)$ 做一次完整的数值对照把一阶、二阶近似和真实值三列摆在一起同时给出误差比值看它和理论量级是否吻合。import numpy as np def f(x, y): return np.exp(x) * np.log(1 y) def p1(x, y): # 一阶近似y return y def p2(x, y): # 二阶近似y xy - y^2/2 return y x * y - y * y / 2 def p3(x, y): # 三阶近似用于看误差是否按预期下降 return y x * y - y * y / 2 (x * x * y) / 2 - (x * y * y) / 2 (y ** 3) / 3 cases [(0.01, 0.01), (0.05, 0.08), (0.1, 0.15), (0.2, 0.2)] print(f{h:6} {k:6} {real:14} {p1:14} {p2:14} {err1:10} {err2:10}) for h, k in cases: r f(h, k) a1 p1(h, k) a2 p2(h, k) print(f{h:6.2f} {k:6.2f} {r:14.10f} {a1:14.10f} {a2:14.10f} f{abs(r-a1):10.2e} {abs(r-a2):10.2e})三阶近似里我手工补了展开到三次的项用于验证误差下降的规律$x^2y/2$ 来自 $e^x$ 的二次项乘 $\ln(1y)$ 的一次项$-xy^2/2$ 来自两个一次项相乘$y^3/3$ 来自 $\ln(1y)$ 的三次项。6.2 结果与截断误差量级分析按上面的参数跑会看到几个很有说明性的现象。以 $(h,k)(0.05,0.08)$ 这组为例真实值大约是 $0.08090692$一阶近似是 $0.08$误差约 $9.07\times10^{-4}$二阶近似是 $0.0808$误差约 $1.07\times10^{-4}$。误差从 $10^{-3}$ 掉到 $10^{-4}$正好差一个量级和一阶误差 $O(|v|^2)$、二阶误差 $O(|v|^3)$的预测一致这里 $|v|\sim0.1$平方是 $10^{-2}$立方是 $10^{-3}$再乘上偏导系数的量级落到 $10^{-3}$ 和 $10^{-4}$ 很合理。当扰动扩大到 $(0.2,0.2)$ 时变化就明显了。$|v|\sim 0.28$二阶误差理论上应该落在 $10^{-2}$ 量级实测也确实在这个范围。这时候一阶近似已经开始不能用了相对误差接近 15%。这个数字给我的直观感受是二阶展开的舒适区大约在扰动小于展开点尺度 10% 的时候超过 20%就该考虑换展开点或者提高阶数了。另外有一个容易被忽略的现象误差的符号。在上述几组数据里误差始终是正的真实值大于近似值说明二阶近似在这个区域内是系统性低估。这不是随机误差不会因为多测几次而抵消。如果拿这个近似去做批量计算的中间步骤这个系统偏差会累积起来最终结果的偏差方向是可预测的但幅度可能被放大。修正办法很简单就是再补上三阶项或者把展开点挪到更靠近实际工作区间的位置。6.3 一个改善精度的实用技巧最后分享一个我在做参数化模型时用得比较顺手的技巧不要在固定点上展开而是在工作区间的中点展开。原因不难理解。截断误差和展开点到目标点的距离的三次方成正比。如果你把展开点选在区间一端那么另一端的误差就最大整个区间的最大误差由这个最远点决定。移到中点后最远距离减半误差按三次方算就降到原来的八分之一。这是零成本的改进——偏导数的表达式不用改只需要换一下取值点。实测中这条经验的效果比想象中明显。同样是一段 $[0, 0.2]$ 的参数区间端点展开时最大误差是 $8\times10^{-3}$ 量级中点展开降到 $1\times10^{-3}$ 量级正好八倍关系。所以现在我做任何局部近似第一件事就是把展开点定在区间的几何中心然后把整个区间两端都验一遍误差取最大值作为这个近似模型的精度指标写进注释里。这样后来接手的人一眼就能知道这个近似在什么范围内可信比写一句误差较小有用得多。如果要追求更高的精度但不想增加阶数还有一个办法是把大区间切成若干小区间每段各自在中点展开然后在边界处平滑衔接。代价是系数表变长好处是每一段里都是三次量级的小误差。这套思路其实就是典型的查表法配合局部多项式在嵌入式设备上做超越函数计算时非常常见我后来在做实时估算模块时也是这么处理的。
