函数的线性化形式和反向传播

物理上，如果 f(x) 中的 x 表示时间，即写成 f(t) ，那么导数的意义是瞬时变化率，这个解释下更容易写出数学化后的导数定义：

\[ f'(t) = \lim_{\Delta t\to 0} \frac{f(t+\Delta t)-f(t)}{\Delta t} \]

如果将 y=f(x) 绘制在笛卡尔坐标系下，那么 f'(x) 在几何上表示的是 y=f(x) 图像在给定点的切线的斜率，可以用这种性质去做很多事情，比如线性近似：

在点 x0 附近，函数值 f(x) 可以用该点切线去近似 \( f(x) \approx f(x_0) + f'(x_0)(x - x_0) \)。

注意为了能保持近似 x 需要控制在和 \( x_0 \) 很近的范围，但这个范围是多少需要依据实际问题需求而定。

基于这种切线近似的思路可以发展出函数求根的牛顿法。

值得一提的是，用导数作为切线的斜率实际是一种定义，而不是逻辑推导的结果，因为切线这个概念几乎只是为初等几何里的“圆”量身定制的，它表示的是和圆只有一个交点的直线，但这个定义无法扩展到一般曲线，比如 \( y=x^3 \) 在 x=0 点处的任何直线必然和曲线有一个交点，但我们还是可以用导数来谈论它的切线，导数定义的切线本身是和圆的切线的经典定义相容，且在数学上又是唯一的，可以扩展到所有可微函数上。

现实中我们并不需要这种经典的“和曲线只有一个接触点”的语义来定义或解决现实问题，“作圆的切线”更多是数学练习或者解决和圆有关问题的特殊场景下辅助手段，实际中关注的是函数的变化率和局部线性近似（见后文）。线性化（一种简化问题的手段）才是真正的核心，而不是“切”（恰好接触）的语义概念。

根据几何图像，切线斜率为正说明往 x 大的方向函数值增大，往 x 小的方向则减小，因此它提供了函数极值点相对该点的方向，于是可以发展出迭代优化算法。

举个例子：

\( f(x)=x^2 \) 初始为 x=2 ，用导数信息去迭代寻求最小值，先计算导数 f'(x)=2x, 代入初始点值 f'(2)=4, 它表明如果增大 x 那么 f(x) 也会增大，于是我们应该减小 x ，此时可以用一个很小的步长 \( \eta \) 去乘以 \( -f'(x) \) 并更新 \( x=x-\eta f'(x) \)

可以说，我们用 f'(x) 的符号 sign(f'(x)) 来指示优化的方向，和 f'(x) 的具体大小值无关。

本文中线性化指的是把局部近似变成极限，然后寻找切空间中微分项 (如 df,dx,dy) 之间的线性关系。

切线是对函数的局部近似： \( f(x) \approx f(x_0) + f'(x_0)(x - x_0) \)

这可以图像去直观理解，也可以代数上通过泰勒公式得到 \( f(x) = f(x_0) + f'(x_0)(x - x_0) + o(x-x_0)^2 \)

记 \( \Delta x = x-x_0 \), \( \Delta y = f(x)-f(x_0) \)

理想微积分中，我们令 \( \Delta x \) 无限趋近于零来丢弃高阶项 \( o(\Delta x^2) \) ，因为相比 \( \Delta x \) 来说，高阶项太小了，或者说二者比值极限为 0。

在实践中可以避开无穷小量的争议，比如计算机里浮点数的舍入特性本身就“物理地”实现了丢弃高阶项这一点：

只要选择一个合适的扰动 \( \Delta x \)（例如经典的 \( 10^{-5} \) 或 \( 10^{-8} \)），使得 \( (\Delta x)^2 \) 的量级刚好低于当前浮点数格式（双精度约为 \( 10^{-16} \)）的舍入分辨率。此时，高阶项在硬件乘法运算中会被自然地截断为 0，对结果没有影响；

而线性项 \( f'(x_0) \cdot \Delta x \) 由于量级足够大无法被截断，因此保留在计算结果中。

此时数值差分公式 \( \Delta f \approx f'(x_0) \Delta x \) 在硬件层面退化为一个实际上的等式，因为“约等于”所舍弃的部分（即高于一阶的项），在物理内存中已经彻底消失了，比如 3+1e-16 还是等于 3 。

此时，导数 \( f'(x_0) \) 就是输入扰动 \( \Delta x \) 与输出变化 \( \Delta f \) 之间精确的缩放倍数。

在该精度下，定义自变量的变化量 \( dx = x - x_0 \) ，以及输出变化量 \( df=f(x)-f(x_0) \)，称它们为微分，于是有：

\[ df = f'(x_0) \, dx \]

导数被解释成对自变量微分 dx 的一个变换，它面对的问题是：在输入端施加一个足够小的扰动，会导致输出端多少的线性程度的变化？

而回答是，以导数 \( f'(x_0) \) 的倍数去缩放。

当函数的输入是两个变量比如 z=f(x,y) 时，输入的所有可能组合 (x,y) 构成了一个平面，而导数的“瞬时变化率”解释蕴含着一个特定的时间方向，于是我们只能说在 x,y 平面上某个（被时间参数化）方向函数值 z 的瞬时变化率。

比如在 x 方向（即 \( \vec{u}= 1 \hat{i}+0 \hat{j} \)）的瞬时变化率是偏导 fx

\[ \frac{\partial f}{\partial x} = \lim_{\Delta x\to 0} \frac{f(x+\Delta x, y)-f(x,y)}{\Delta x} \]

类似地可以得到 y 方向的偏导 fy。

而任意方向可以用向量 \( \vec{u} = a\hat{i} + b\hat{j} \) 表示，这里约定 \( a^2+b^2=1 \)，即 u 是单位向量。

函数在 u 方向变化量为 \( \Delta t \) 就对应了 f(x,y) 移动到了 \( f(x+a \Delta t, y+b \Delta t) \), 这样 f(x,y) 在 u 方向的导数就定义为：

\[ D_{u}f = \lim_{\Delta t \to 0} \frac{f(x+a\Delta t, y+b \Delta t) - f(x, y)}{\Delta t} \]

添加了“方向条件”后的导数定义在“瞬时变化率”解释上和单变量函数的导数解释就一致了。

而如果函数可微，可以在分母引入一个中间量：

\[ D_{u}f = \lim_{\Delta t \to 0} \frac{f(x+a\Delta t, y+b \Delta t) - f(x, y+b \Delta t) + f(x, y+b \Delta t) - f(x, y)}{\Delta t} \]

极限部分拆成两个独立的极限，分别对应 x 方向的净变化与 y 方向的净变化：

\[ D_{u}f = \lim_{\Delta t \to 0} \frac{f(x+a\Delta t, y+b\Delta t) - f(x, y+b\Delta t)}{\Delta t} + \lim_{\Delta t \to 0} \frac{f(x, y+b\Delta t) - f(x, y)}{\Delta t} \]

第二个极限完全是固定 x 不动，只让 y 从 y 变化到 \( y+b \Delta t \)。令 \( \Delta y = b \Delta t \)，则分母的 \( \Delta t \) 就等于 \( \Delta y/b \)，于是这一项写成：

\[ b \cdot \lim_{\Delta y \to 0} \frac{f(x, y+\Delta y) - f(x, y)}{\Delta y} \]

而这就是 \( b \cdot f_y \) 。

第一个极限的分子是在 x 方向上的增量，只是此时 y 的坐标被临时固定在了 y+bΔt 位置上。令 Δx = aΔt，则这一项写成：

\[ a \cdot \lim_{\Delta x \to 0} \frac{f(x+\Delta x, y+b\Delta t) - f(x, y+b\Delta t)}{\Delta x} \]

虽然分子里出现的中间坐标是 y+bΔt 而不是 y，但当 Δt 趋于 0 时，bΔt 也必然趋于 0，所以这个临时坐标 y+bΔt 会在极限下回到 y。因此该极限值就是偏导数 \( f_x \)，第一项的结果就是 \( a \cdot f_x \)。

把两部分合并起来就得到了方向导数基于偏导数的定义 \( D_{u}f = a f_x + b f_y \)。

这个式子体现了一种复杂到简单的规约：对于二维输入函数 f(x,y), 给定某个单位向量 u 表示方向，在该方向的导数可以规约成 x 和 y 方向偏导数的线性组合，而偏导数完全是单变量函数求导的结果，因此当建立起了单变量求导函数库和规则后，任何方向导数的计算问题都可以规约到普通单变量求导问题中。

单变量函数导数的几何解释是切线斜率，而双变量函数中与之对应的是特定方向的切线的斜率，如果函数在该点光滑可微，这些切线实际都在一个平面上，这个平面就是切平面。

对此进行简单证明，根据前文方向导数的定义，给定任何一个单位向量 (a,b) 表示方向，可以到该点在该方向上变化率为 \( af_x+bf_y \) ，这意味着该方向的切线的方向可以用 \( (a, b, af_x+bf_y) \) 向量来表示。

所以我们得到了一族直线，它们都经过点 (x,y,f(x,y)) ，方向则被 \( \vec{v} = (a, b, af_x+bf_y) \) 定义。

但如果我们引入一个新的向量 \( \vec{n} = (f_x, f_y, -1) \) 会发现 \( \vec{n} \cdot \vec{v} = 0 \) ，也就是所有经过改点的直线都和向量 \( \vec{n} \) 正交。

这表明各个方向的切线都在同一个以 \( \vec{n} \) 为法向量的平面上（考虑到这些切线都经过同一个点）。

因此我们用 \( (f_x,f_y, -1) \) 法向量以及当前点 (x,y,f(x,y)) 就能刻画所有切线，这个法向量和点组成了切平面的数学表征。

作为对比，用法向量的语言重新表述单变量场景：对于曲线 y = f(x) ，在 x0 点的切线的斜率为 \( f'(x_0) \) ，切线方程就是 \( y - f(x_0) = f'(x_0)(x - x_0) \)。但切线的方向向量可以写为 \( \vec{v} = (1, f'(x_0)) \)，因为当 x 增加 1 个单位时，y 沿着切线增加 \( f'(x_0) \) 个单位。

于是可以找一个向量 \( \vec{n} = (f'(x_0), -1) \)，计算点积：

\[ \vec{n} \cdot \vec{v} = f'(x_0) \cdot 1 + (-1) \cdot f'(x_0) = 0 \]

这说明 \( \vec{n} \) 与切线方向垂直，即 \( \vec{n} \) 是切线的法向量。

注意，双变量的切平面法向量是 \( (f_x, f_y, -1) \)，单变量函数的切线法向量是 \( (f'(x_0), -1) \)，后者是前者的“降维版本”——去掉 y 方向的偏导分量。而由于其中 -1 是常数，因此从切线的属性上看，当在单变量函数中谈论导数 \( f'(x_0) \) 时，在双变量函数中谈论的对应物应该是 \( (f_x, f_y) \) 。（从变化率上看对应的是方向导数）

由于方向导数是 \( D_{u} f = a f_x + b f_y \) ， (a,b) 是单位向量，那么这个结果在几何上解释就是 \( (f_x, f_y) \) 在 (a,b) 方向上的投影的长度，而从几何上就可以看出，只有 (a,b) 方向和 \( (f_x,f_y) \) 一致的时候这个投影最大，即等于 \( (f_x,f_y) \) 模长，这是一个正数，也就是它是函数增长最大的反向。

这意味着在 \( (f_x,f_y) \) 方向上，函数的变化率是最大的，此时我们定义这个方向为梯度，用 \( \nabla f \) 表示。

注意 \( (f_x,f_y) \) 在多个地方出现，首先是方向导数中，它作为基，或者作为一种内积运算的对象，用来把方向转成变化率；在切面中作为法向量的部分分量；而在梯度中，指示的增长最快的方向。虽然三者都涉及偏导，但意义是不同的。

在寻找最值过程中，我们可以找到梯度方向，然后移动一小步向上爬升（寻找极大值），或者反方向移动向下移动寻找极小值。

和单变量函数优化场景类似，我们需要的只是梯度所指示的方向，而不需要它的模长。甚至也不需要严格按照梯度方向去更新，因为它指示的是变化率最大的方向，如果找最大值，只需要保证每次函数值都上升，而更新的方向与梯度方向的夹角为锐角即可满足，因为方向导数 \( \nabla f \cdot \hat{u} = \mid\nabla f\mid \cos(\theta) \) 表明，只要夹角 θ 小于 90 度，方向导数就为正，函数值就会上升。

所以我们甚至可以只朝着 \( (sign(f_x), sign(f_y)) \) 方向小步更新，即只关注各个偏导的符号，这与只关注单变量函数导数的符号进行更新是一致的。

换句话说，梯度方向是上升最快的方向，但绝不是唯一的上升方向——凡是和梯度指向同一半平面内的方向，都能让函数值有所增加，只是增加的快慢不同罢了。

回顾单变量函数的线性化，给定一个极小的扰动 dx（不需要假设无穷小量，而是假设在人类可接受精度边缘的扰动），那么输出端对扰动的响应是 \( df=f'(x)dx \)

这里线性化和导数的区别在于，我们不是关心变化率，而是关心具体的某个微小扰动产生的绝对变化。

因此原本导数可以看作是 \( \frac{\Delta f}{\Delta x} \) 的极限，或者写成 \( \frac{df}{dx} \), 线性化则只关心 \( df \) 和 \( dx \) 的函数关系。

之所以称为线性化，是因为扰动的高阶的部分如 \( \frac{1}{2}f''(x)dx^2 \) 中 \( dx^2 \) 已经比需求精度还要低，不影响实际用途了，应用中完全可以忽略或者浮点截断了，只有关于误差的线性缩放部分精度仍然在人类要求的边缘，因此被保留。

对于二维，目前我们通过“变化率”的解释去定义出了方向导数：

\[ D_{u}f = \lim_{\Delta t \to 0} \frac{f(x+a\Delta t, y+b \Delta t) - f(x, y)}{\Delta t} \]

而如果关注微小扰动的绝对影响，那么应该用变化率乘以时间，对应的数学表达是：

\[ df_u = \lim_{\Delta t \to 0} \Delta t D_{u}f = \lim_{\Delta t \to 0} f(x+a\Delta t, y+b \Delta t) - f(x, y) \]

用和求方向导数类似的技巧拆分：先在 x 方向上单独制造一个变化，再在 y 方向上补上剩余的变化：

\[ \Delta f = f(x+a\Delta t, y+b\Delta t) - f(x, y+b\Delta t) + f(x, y+b\Delta t) - f(x, y) \]

继续沿用和方向导数类似的推导思路，令 \( dx = a\Delta t \)，第一项在极限下是 \( f_x \cdot dx \)，令 \( dy = b\Delta t \)，第二项在极限下等于 \( f_y \cdot dy \)

将两部分加总，得到：

\[ df = f_x dx + f_y dy \]

注意这里的 \( dx = a\Delta t \) 和 \( dy = b\Delta t \) 不再是抽象的单位方向分量 (a,b)，而是自变量空间中施加在 x 和 y 两个坐标轴上的具体扰动量，也就是说 (dx,dy) 里已经蕴含了方向，它是带方向的极小变化。

如果我们把输入 (x,y) 向量化为 \( \vec{x} \), 那么 \( (dx,dy) \) 就可以记为 \( d\vec{x} \), 由于有 dx 和 dy 两个变化方向，所有的微小扰动是在 dx 和 dy 构成的平面里，就像是由标准线性空间的 x 和 y 轴张成的平面一样，这种输入端微小扰动所在的空间称为切空间。

那么双变量函数的线性化数学化表达为 \( df = \nabla f \cdot d\vec{x} \), 这和单变量的 \( df = f'(x)dx \) 形式上是统一的，向量点积对应了标量乘法，梯度对应了导数，梯度的优化意义和导数的优化意义对应了起来，梯度作为切平面方向量的分量和导数作为切线法向量分量的意义也对应了起来。

但关键区分在于， \( \nabla f \) 不是变化率，但 \( f'(x) \) 是变化率。

另外，在二维情形下，输出端的变化量是通过向量内积给出的，因此无法直接写成 \( \frac{df}{d\vec{x}} = \nabla f \) ，因为向量作为除数是没有被定义的数学操作。

若要恢复“变化率”的形式，必须引入一个中间变量 dt 并指定具体的单位方向 \( \hat{u} \)，使得沿该方向的变化率为 \( \frac{df}{dt}\big|_u = \nabla f \cdot \hat{u} \)。

因此在线性化视角下，从单变量升到双变量输入后，导数 f'(x) 的三种解释（物理变化率，切线，对扰动关系的线性化）中的第一种在代数形式上并没有得到自然扩充（必须依赖特定方向），但切平面和线性化是被保留了。

而这种通用性或者一般性是数学所偏爱的，从这个角度看，导数的变化率解释只是一维输入下对函数线性化得到的特例，我们正好能得到一个缩放因子。这也使得在单变量函数求导中把导数看作两个微分的除法 \( \frac{df}{dx} \) 并不会有什么问题。

但在二维中，无法把方向导数写成 \( \frac{df}{d\vec{x}} \), 这个表达式也不能说是梯度，因为标量和向量的除法本身就不是良好定义的，用编程语言的话来说，这是一种类型错误。必须用额外的方向来作为前台谈论导数。

因此最标准的方式就是不去写成变化率的除法形式，而只是展示 df 和 dx 的函数关系： \( df = \nabla f \cdot d\vec{x} \) 。

此时 \( \nabla f \) 被解释为对切空间里扰动的映射，就像 f'(x) 是对单变量扰动 dx 的映射一样。

考虑例子 \( z = x^T y \) ， x 和 y 都是 n 维的向量，如果把输入看作集合对每个 \( x_i,y_i \) 求偏导，会得到对 \( x_i \) 的偏导是 \( y_i \) ，对 \( y_i \) 的偏导则是 \( x_i \) ，因此综合起来看， \( z_x = y \) 而 \( z_y = x \), 即我们可以对整个向量来谈论梯度（这里梯度始终看作列向量）。

考虑直接对 z 线性化，即同时给所有输入一个在需求可接受精度边缘的扰动，这个扰动被函数结构拆分成 dx 和 dy, 于是扰动后结果 \( z'=(x+dx)^T(y+dy) \) ，展开得到：

\[ z' = x^T y + x^T dy + dx^T y + dx^T dy \]

第一项 \( x^T y \) 是原来的输出值 \( z \)，左右两边同时减去 z，左边就是值域的变化量：

\[ dz = x^T dy + dx^T y + dx^T dy \]

最后一项 \( dx^T dy \) 是两个一阶小量的乘积的求和 \( \sum dx_i dy_i \) ，而每个 \( dx_i \) 和 \( dy_i \) 都在需求可接受的精度边缘， \( dx_i dy_i \) 会被 round 到 0 ，不对结果产生任何影响，只留下关于 dx 和 dy 的线性部分：

\[ dz = x^T dy + dx^T y \]

\( dx^Ty \) 是标量，它等于自身转置 \( y^Tdx \) ，于是式子可以写成：

\[ dz = y^T dx + x^T dy \]

如果部分输入是常数，比如 y 向量里所有值是常数，那么 \( dz=y^T dx \) ， 即我们得到了 z 对 x 的梯度是 y （或者说偏导是 \( y^T \)）。

同理，如果 x 向量里所有输入是常数，那么可以看出 z 对 y 的梯度是 x, 偏导则是 \( x^T \)

如果把 x 和 y 看作一个完整的参数列向量 \(\theta = \begin{bmatrix} x \\ y \end{bmatrix}\)，那么对应的扰动就是 \( d\theta = \begin{bmatrix} dx \\ dy \end{bmatrix} \)。此时线性化式子 \( dz = y^T dx + x^T dy \) 可以重新写成：

\[ dz = \begin{bmatrix} y \\ x \end{bmatrix}^T \begin{bmatrix} dx \\ dy \end{bmatrix} \]

也就是说，完整的梯度向量是 \( \nabla_{\theta} z = \begin{bmatrix} y \\ x \end{bmatrix} \)。

这种写法主要是为了用例子说明，无论输入参数有多少、如何分组，只要把它们堆成一个列向量，线性化形式始终保持 \( dz = \nabla z \cdot d\theta \) 这种向量点积形式，而点积的实现是 \( \sum z_i \theta_i \) ，于是它可以在代数上拆解成多个不同组的向量点积，这正是直接用线性化方法得到的 \( dz = y^T dx + x^T dy \) 形式所展现的。

再看一个更复杂的例子 \( f(x,A)=x^TAx \) ，代入扰动后：

\[ f(x+dx, A+dA) = (x+dx)^T (A+dA) (x+dx) \]

展开，保留一阶项：

\[ = x^T A x + x^T A dx + x^T dA x + dx^T A x + \dots \]

第一项 \( x^T A x \) 是原函数值。省略号表示二阶项（包括 \( dx^T A dx \)、\( dx^T dA x \)、\( x^T dA dx \) 等，这些项中一定会包括两个极小量的乘积，而由于两个极小量都已经在精度边缘，因此相乘就被丢弃了），线性部分就是后面三项：

\[ df = x^T A dx + x^T dA x + dx^T A x \]

有两项包含 \( dx \)：\( x^T A dx \) 和 \( dx^T A x \) ，因为两者都是标量，互为转置，因此可以对第二个转置后合并：

\[ x^T A dx + (A x)^T dx = \big( x^T A + (A x)^T \big) dx = \big( A^T x + A x \big)^T dx \]

所以 \( \nabla_x f = (A + A^T)x \)。

再看 dA 相关的项 \( x^T dA x\) 。注意 dA 是一个矩阵，而 \( x^T dA x \) 是标量，为了读出梯度（与 dA 做内积得到这个标量），需要把表达式改写为内积形式。

先利用 tr(ABC)=tr(CAB) 性质有：

\[ x^T dA x = \text{tr}(x^T dA x) = \text{tr}(x x^T dA) \]

还要利用 tr 的这样一种性质：如果 B,A 是矩阵，形状分别为 (m,k) 和 (k,n)，那么 tr(BA) 结果相当于是把 B 和 A 向量化之后求内积。但我们要注意这里 B 和 A 的向量化方式，比如 BA 的对角线上第一个元素是 B 的第一行和 A 的第一列的内积，因此 B 是按行方式展开，A 则是按列方式展开，由于我们是以 A 为准的（它对应了微分向量），因此如果要写成内积，应该是 \( tr(BA) = \langle B^T, A \rangle \)

但由于 \( x x^T \) 的转置还是自身，所以 \( tr(x x^T dA) \) 等价于 \( \langle x^T x, dA \rangle \)

即 \( \nabla_A f = x x^T \)，它是一个 n×n 矩阵，与 dA 同形状。

最终 \( df = (A^Tx + Ax)^Tdx + tr(x x^TdA) \)

如果把输入变量看作一个集合，把 A 看作常数矩阵，那么 \( f(x) = x^T A x \) 本质上就是一个关于 n 个输入分量 \( x_1, x_2, \ldots, x_n \) 的普通多元函数，只不过这些分量通过矩阵 A 耦合在一起而已。如果硬要展开，它就是 \( \sum_{i,j} A_{ij} x_i x_j \)，对每个 \( x_k \) 求偏导，得到 \( \partial f/\partial x_k = \sum_j (A_{kj} + A_{jk}) x_j \)，然后把这些偏导堆成一个向量，就得到了梯度 \( \nabla_x f = (A + A^T)x \)。这和以上线性化结果一致，但在线性化操作中，输入的结构始终被运算结构维持着，不需要先手动去掉结构，然后再组装回去。

如果把 A 看作变量， x 为常向量，那么就是 \( f(A) = x^T A x \) 对 A 求梯度。从集合视角看，它就是一个 n^2 元函数，每个矩阵元素 \( A_{ij} \) 都是一个输入分量。求偏导时，固定其他元素不动，只让 \( A_{ij} \) 变化，得到 \( \frac{\partial f}{\partial A_{ij}} = x_i x_j \)，然后把这些偏导按行和列排成一个矩阵，就得到了梯度矩阵 \( x x^T \)。这和以上线性化结果一致，但线性化操作中始终保持了模块化信息。

由于 \( \nabla f \) 里是由各个变量的偏导组成，因此对于其中元素，其实都是在进行单变量求导，于是我们先理清单变量求导的链式法则

设 \( f: [a, b] \to [c, d] \) 与 \( g: [c, d] \to \mathbb{R} \) 为函数，其中 f 在 x0 处可导，g 在 y0 = f(x0) 处可导。

我们要证明复合函数 \( g \circ f \) 在 x0 处可导，且其导数为： \[ (g \circ f)'(x_0) = g'(y_0) f'(x_0). \]

令 \( y = f(x) \)，\( y_0 = f(x_0) \)。首先假设 \( f'(x_0) \neq 0 \)。在此情形下，当 \( x \neq x_0 \) 但足够接近 x0 时，有 \( y \neq y_0 \)，于是我们可以将差商写成如下形式： \[ \frac{g(f(x)) - g(f(x_0))}{x - x_0} = \frac{g(y) - g(y_0)}{y - y_0} \cdot \frac{f(x) - f(x_0)}{x - x_0}. \] 由于 f 在 x0 处可导，当 \( x \to x_0 \) 时，有 \[ \frac{f(x) - f(x_0)}{x - x_0} \to f'(x_0). \] 当 \( x \to x_0 \) 时，\( f(x) = y \to f(x_0) = y_0 \)。由此推出，当 \( x \to x_0 \) 时， \[ \frac{g(y) - g(y_0)}{y - y_0} \to g'(y_0). \] 综合以上结果，得到： \[ \frac{g(f(x)) - g(f(x_0))}{x - x_0} = \frac{g(y) - g(y_0)}{y - y_0} \cdot \frac{f(x) - f(x_0)}{x - x_0} \to g'(y_0) f'(x_0). \] 这就证明了 \( f'(x_0) \neq 0 \) 情形下的链式法则。

当 \( f'(x_0) = 0 \) 时，我们像上面一样进行论证，但需要更加小心，因为在这种情况下，即使 \( x \neq x_0 \)，也可能出现 \( y = y_0 \)。对于满足 \( y = y_0 \) 的 x，差商为零；而对于满足 \( y \neq y_0 \) 的 x，我们可以像上面那样处理，并将差商重写为两个因子的乘积。在两种情形下，我们都能得到极限为零，从而证明了链式法则的剩余情形。

当要把链式法则扩展到参数模块化的函数时，会遇到这种场景，比如虽然 g(f(x)) 的输出是一个标量，但 f(x) 是一个向量， x 也是向量，那么一般的 \( (g'(f(x_0)) = g'(y_0) f'(x_0). \) 中 \( g'(y_0) \) 是一个向量，而最终梯度也是向量，那么从运算合法性看， f'(x0) 就必须是一个矩阵。

假设 x 是 \( R^n \) 空间里的向量，f(x) 输出两个中间变量，我们先只关注输入中的特定变量 t, 它经过 f(x) 后输出 h(t) 和 k(t), 即对 t 这个输入变量来说，它最终输出是 \( g(t)=f(h(t),k(t)) \)

我们想知道 \( g'(t) \)。按定义写差商，给 t 一个微小增量 \( \Delta t \)：

\[ \frac{g(t+\Delta t) - g(t)}{\Delta t} = \frac{f(h(t+\Delta t), k(t+\Delta t)) - f(h(t), k(t))}{\Delta t} \]

可以用计算方向导数时完全相同的技巧，在分子里插入一个中间项 \( f(h(t), k(t+\Delta t)) \)：

\[ = \frac{ f(h(t+\Delta t), k(t+\Delta t)) - f(h(t), k(t+\Delta t)) + f(h(t), k(t+\Delta t)) - f(h(t), k(t)) }{\Delta t} \]

拆成两个独立的极限：

\[ g'(t) = \lim_{\Delta t \to 0} \frac{ f(h(t+\Delta t), k(t+\Delta t)) - f(h(t), k(t+\Delta t)) }{\Delta t} + \lim_{\Delta t \to 0} \frac{ f(h(t), k(t+\Delta t)) - f(h(t), k(t)) }{\Delta t} \]

这两个极限都是前面证明链式法则时的极限，分别对应 \( \frac{\partial f}{\partial h} \cdot h'(t) \) 和 \( \frac{\partial f}{\partial k} \cdot k'(t) \)

因此最终不同路径上应用链式法则得到的梯度信息用求和的方式汇聚起来：

\[ g'(t) = \frac{\partial f}{\partial h} h'(t) + \frac{\partial f}{\partial k} k'(t) \]

这个模式可以推广到 f(x) 为任意 n 个输出变量的场景。

于是之前线性化的一般公式 \( dg = \nabla_y g \cdot d\vec{y} \) 中如果有 t 对 \( \vec{y} \) 做贡献，那么 g 对 t 的导数就是 \( g_t = \nabla_y g \cdot \frac{d\vec{y}}{dt} \)。

这里 \( \frac{d\vec{y}}{dt} \) 是一个向量，它的第 i 个分量是 \( \frac{dy_i}{dt} \)，即每个中间变量对 t 的变化率。于是上面的求和对所有中间变量求和，本质上就是梯度向量 \( \nabla_y g \) 与 \( \frac{d\vec{y}}{dt} \) 做内积。

而这只是对输入中某个变量 t 的偏导计算，如果输入有更多变量，写成 \( \vec{x} \)，那么对每个 \( \vec{x}_i \) 的偏导都是以上内积形式，而且它们共享同一个 \( \nabla_y g \)，只是每个 \( \frac{d\vec{y}}{d\vec{x}_i} \) 不同。因此这可以看作同一个向量 \( \nabla_y g \) 与多个不同向量并行做内积。

把每个 \( \frac{d\vec{y}}{d\vec{x}_i} \) 向量按行排列成矩阵，再和列向量 \( \nabla_y g \) 做矩阵向量乘法，得到的列向量的各个元素就是对应的各个 \( x_i \) 的梯度：

\[ \nabla_{\vec{x}} g = J_{\vec{y}}^T \nabla_{\vec{y}} g \]

这里 \( J_{\vec{y}} \) 称为雅可比矩阵，我们用它的转置和 \( \nabla_{\vec{y}} g \) 相乘，意味着 \( J_{\vec{y}} \) 的第 j 列是 \( \partial \vec{y}/\partial x_j \)。

如果 \( \vec{x} \in \mathbb{R}^m \), 而 \( \vec{y} \in \mathbb{R}^n \), 那么雅可比矩阵 \( J_{\vec{y}} \) 的形状就是 \( n \times m \)（n 个输出，m 个输入），而 \( \nabla_{\vec{y}} f \) 是 \( n \times 1 \) 的列向量，两者通过转置相乘得到 \( \nabla_{\vec{x}} f \)，其形状为 \( m \times 1 \)，与输入向量 \( \vec{x} \) 的形状完全一致。

为什么 J 的形状约定是 nxm 而不是 mxn 呢？

这是就需要从线性化的一般形式来谈论了。

如果对 \( \vec{y} = f(\vec{x}) \) 进行线性化，根据前文的描述，每个输出 \( y_i \) 的变化量可以写成：

\[ dy_i = \frac{\partial y_i}{\partial x_1} dx_1 + \frac{\partial y_i}{\partial x_2} dx_2 + \cdots + \frac{\partial y_i}{\partial x_m} dx_m \]

封装起来是梯度（列向量）和切空间扰动向量（列向量）的点积： \[ dy_i = \left( \nabla_{\vec{x}} y_i \right) \cdot d\vec{x} \]

但这种写法是两个列向量的点积，而不是一般的矩阵乘法。如果我们要写成矩阵乘法的标准形式 ，那么梯度列向量 \( \nabla_{\vec{x}} y_i \) 就必须转置成行向量，我们一般把转置后的行向量称为导数（或局部线性映射的坐标表示）：

\[ dy_i = \frac{\partial y_i}{\partial \vec{x}} d\vec{x} \]

这里 \( \frac{\partial y_i}{\partial \vec{x}} \) 是一个 \( 1 \times m \) 的行向量，其第 j 个分量正是 \( \partial y_i / \partial x_j \)。在这个形式下，单个输出的线性化表现为“行向量乘以列向量”，结果是一个标量。

然后扩展到 n 个不同的 \( y_i \)，我们只是把多个梯度的转置行向量（导数）上下堆叠，那么和切空间扰动向量的矩阵乘法形式保持不变：

\[ \begin{bmatrix} dy_1 \\ dy_2 \\ \vdots \\ dy_n \end{bmatrix} = \begin{bmatrix} \frac{\partial y_1}{\partial x_1} & \frac{\partial y_1}{\partial x_2} & \cdots & \frac{\partial y_1}{\partial x_m} \\ \frac{\partial y_2}{\partial x_1} & \frac{\partial y_2}{\partial x_2} & \cdots & \frac{\partial y_2}{\partial x_m} \\ \vdots & \vdots & \ddots & \vdots \\ \frac{\partial y_n}{\partial x_1} & \frac{\partial y_n}{\partial x_2} & \cdots & \frac{\partial y_n}{\partial x_m} \end{bmatrix} \begin{bmatrix} dx_1 \\ dx_2 \\ \vdots \\ dx_m \end{bmatrix} \]

写成： \( d\vec{y} = J d\vec{x} \)

因此在这种线性化统一视角下，雅可比矩阵被定义为 \( n \times m \)，即它把输入的 m 维扰动向量线性映射到了输出的 n 维扰动向量，且雅可比矩阵和导数（梯度的行向量）是同一种对象。

对于单层情况，扰动是通过 \( d\vec{y} = J d\vec{x} \) 公式计算的。如果有多层，比如三层，设 \( \vec{y} = f_1(\vec{x}) \)，\( \vec{z} = f_2(\vec{y}) \)，\( \vec{w} = f_3(\vec{z}) \)，那么根据线性映射的复合性质，总的扰动是各层雅可比矩阵的乘积：

\[ d\vec{w} = J_{f_3} d\vec{z} = J_{f_3} J_{f_2} d\vec{y} = J_{f_3} J_{f_2} J_{f_1} d\vec{x} \]

这里雅可比矩阵按前向传播的顺序从左到右乘起来。注意矩阵乘法不满足交换律，所以顺序不能颠倒：第一层雅可比 \( J_{f_1} \) 必须放在最右边（最先作用在 \( d\vec{x} \) 上），最后一层雅可比 \( J_{f_3} \) 放在最左边（最后作用）。

\( J_{f_3},J_{f_2},J_{f_1} \) 只是各个层之间的局部雅克比矩阵，它们只要按照复合顺序进行任意组合得到的都是一个新的雅克比矩阵。

比如乘积 \( J_{f_3} J_{f_2} J_{f_1} \) 就是整个复合映射 \( F = f_3 \circ f_2 \circ f_1 \) 的雅可比矩阵，它把输入端的微小扰动直接映射到输出端的微小扰动。

而 \( J_{f_3} J_{f_2} \) 本身也是一个雅可比矩阵，描述的是从 \( \vec{y} \) 到 \( \vec{w} \) 的复合映射 \( f_3 \circ f_2 \) 的局部线性化。同样，\( J_{f_2} J_{f_1} \) 描述的是从 \( \vec{x} \) 到 \( \vec{z} \) 的复合映射 \( f_2 \circ f_1 \) 的局部线性化。任意连续层的雅可比矩阵按顺序相乘，得到的依然是某个复合映射的雅可比矩阵。

矩阵之间乘法的代价是很高的，对于 nxn 的方阵，复杂度是 \( O(n^3) \)，而输入 \( d\vec{x} \) 是向量，按顺序和一个个局部雅克比矩阵相乘，每次都只是 O(n^2) 的代价，因此跨层的雅克比矩阵基本不会被显式计算出来。

而接下来我们会看到，实际中单层内的雅克比矩阵计算代价也很高，很多时候也不会计算出来。

深度学习中，反向传播目的是计算出损失函数（标量）对特定参数的梯度，并用梯度信息去更新参数。

注意梯度是 \( df = \nabla f \cdot dW \) 中的 \( \nabla f \), 这里 f 必然是标量函数，因为只有标量才能谈论增大和减小方向。

\( \nabla f \) 维度和 dW 一样（这样才能用它去更新 dW），但 dW 更多时候是矩阵，那么 \( \nabla f \) 也是维度和 dW 一样的矩阵，此时可以写成 \( df = tr((\nabla f)^T \quad dW) \), 用迹来表示矩阵的内积，但这里需要对梯度转置，转置的结果称为导数（偏导）或者雅克比矩阵（但雅克比矩阵转置后不一定是梯度，因为雅克比矩阵可以适用于输出为非标量的函数的线性化）。

任意梯度可以通过局部雅克比矩阵复合后和最上游梯度的矩阵乘积形式计算，如三层情况是： \( J_{f_1}^T J_{f_2}^T J_{f_3}^T \nabla_{\vec{y}} L \) 。

这种对各个局部的雅克比矩阵进行转置再相乘的框架看上去非常统一，但具体实现中，追求计算出局部的雅克比矩阵需要很大的代价。

我们以神经网络中最常见的仿射层 y=Wx+b 为例，这里 x 和 y 分别是 n 和 m 维向量， W 是 (m,n) 形状的矩阵。

用线性化的方式去计算 y 对 b 的雅克比矩阵（再次说明，这不是梯度，因为 y 是向量）, 于是对 b 扰动后减去扰动前的值：

\[ dy = Wx+b+db-Wx+b = db \]

这说明 y 对 b 的偏导（雅克比矩阵）是 mxm 形状的单位矩阵 I 。

上游的损失 L 对 y 的梯度假设是列向量 \( \nabla_y L \), 那么 L 对 b 的梯度就是 \( \nabla_b L = I \nabla_y L \)

具体实现中，我们并不需要写出单位矩阵和上游梯度的乘积，只需要返回 \( \nabla_y L \) 即可。

同样的方式计算 y 对 x 的雅克比矩阵：

dy = W(x+dx)+b-(Wx+b) = Wdx 。

因此有 \( \frac{\partial{y}}{\partial{x}} = W \) ， 通过梯度传播误差后 \( \nabla_xL = W^T \nabla_y L \)

由于雅克比矩阵 W 本身就是参数，不需要额外的计算，所以具体实现中会用上游误差梯度乘以这个雅克比矩阵。

但我们看 y 对 W 的雅克比矩阵：

dy = (W+dW)x+b-(Wx+b) = dW x 。

注意，如果要提取出雅克比矩阵 J, 我们要将它写成 dy=J dW' 的形式，这里 W' 是 W 拉成的 mn 维的列向量，而 J 会是一个 (m,mn) 形状的非常大但有非常稀疏的矩阵，这种计算代价完全是没有必要的。

这时候需要用迹技巧，先写出 L 对 y 的线性化： \( dL = \nabla_y L^T dy \), 而因为有 \( dy=dW x \), 同时为了方便阅读，重命名 L 对 y 的导数为 g ，即 \( g = \nabla_y L \), 那么

\( dL = g^T dW x \)

由于 dL 是一个标量，于是可以对它取矩阵的迹：

\[ dL = tr(g^T dW x) \]

然后这时候利用迹的循环置换性质：

\( dL = tr(xg^TdW) \)

根据 tr 的内积形式，这意味着 \( g x^T \) 就是 L 对 W 的梯度，而其转置就是它的导数（雅克比矩阵），即 \( xg^T \) 。

因此我们绕过了 y 对 W 的局部雅克比矩阵，直接得到更为精简的 L 对 W 的全局雅克比矩阵，它就是上游梯度向量和输入 x 的外积。

可以看到，具体实现中，我们不需要把反向传播都规约到统一的局部雅克比矩阵和上游梯度进行向量乘法的形式，许多时候其他等价方式可以更为直接高效地计算出最终结果。

也不需要遵循这种线性化的一般方式，对 W 导数计算反而用单个变量的思路去理解更为简单。

再看前文用线性化方法计算过的 \( f(x,A)=x^TAx \) 的例子，如果用反向传播，可以把它看作 \( x^Tg(x,A) \) 这里 \( g(x,A)=Ax \)

要计算 f 对 A 的偏导，先计算 \( \nabla_g f = x \), 然后根据刚才讨论的 y=Wx 中 y 对 W 的梯度信息可以知道，f 对 A 偏导是 \( x x^T \)

要计算 f 对 x 的偏导，先计算第一层中对 \( x^T \) 的偏导为 Ax, 然后加上内部一层 \( A^Tx \), 最终是 \( (A+A^T)x \)

以上我们讨论的都是数学约定下梯度和雅克比矩阵的计算方式，比如我们总是假设上游误差梯度是列向量，而且矩阵和向量乘积总是标准的 Wx 形式。

但现实实现中，输入往往是 (N,d) 格式或者 (N,d1,d2…) 格式，这里 N 是 batch size 。

我们先关注 X 是 (N,d) 形状下的仿射变换，且输出维度是 (N, h)。

此时一般的实现是（以 numpy 语法为例）： X @ W + b;

或者数学上写成 \( Y=XW+b \)

这种情况下本质上是多行并行的 \( y^T = x^T W+b^T \) 计算。

假设 L 对 Y 的梯度是 G; 由于 Y 是 (N, h) 形状，那么其梯度也是 (N,h) 形状矩阵。

我们用线性化的思路计算损失 L 对 W 梯度，先对仿射层线性化：

\( dY = X(W+dW)+b - (XW+b) = XdW \)

和前文中类似，dW 是矩阵，要写成雅克比矩阵会很复杂，于是直接线性化仿射输入到 L 的所有层：

\( dL = tr(G^T dY) = tr(G^T X dW)\)

利用迹技巧得到 \( dL = tr((X^TG)^T dW) \)

因此梯度 \( \nabla_W L= X^TG \)

要对 b 线性化，先要把 b 扩充为 (N,h) 形状，即每行都复制 N 次，得到 B 矩阵，这相当于 b 经过一个函数输出了 N 个变量，每个输出都是 b:

对 Y 和 B 的关系线性化：

dY = I dB, 直接线性化仿射输入到 L 的所有层：

\( dL = tr(G^T dY)= tr(G^T I dB) \)

此时 \( G^T \) 就是 L 对 B 的梯度，而 L 对 b 的梯度根据多变量输出梯度的累加原则，就是 \( G^T \mathbf{1} \), 即对 G 每行累加得到 (1,h) 形状。

Y 对 X 线性化：

\( dY = (X+dX)W - XW = dX W \)

仍然线性化仿射输入到 L 的所有层：

\( dL = tr(G^TdY) = tr(G^T dX W) \)

利用旋转性质得到： \( dL = tr(WG^T dX) \)

因此 \( \nabla_X L = GW^T \)

如果实现起来，前向传播是：

def affine_forward(x, w, b):
    out = x @ w + b
    cache = (x, w, b)
    return out, cache

这里返回了输出和反向传播所需要的中间变量，反向传播则接受一个 (N,h) 形状的上游梯度 dout, 以及获取中间变量进行计算：

def affine_backward(dout, cache):
    x, w, b = cache
    dw = x.T @ dout
    db = dout.sum(axis=0)
    dx = dout @ w.T
    return dx, dw, db