22.4. 多变量微积分
在 Colab 中打开 Notebook
在 Colab 中打开 Notebook
在 Colab 中打开 Notebook
在 Colab 中打开 Notebook
在 SageMaker Studio Lab 中打开 Notebook

现在我们对单变量函数的导数有了相当深入的理解，让我们回到最初的问题，即我们考虑一个可能包含数十亿权重的损失函数。

22.4.1. 高维微分

第 22.3 节告诉我们，如果我们改变这数十亿个权重中的一个，同时保持其他所有权重不变，我们知道会发生什么！这不过是一个单变量函数，所以我们可以写出

(22.4.1)\[L(w_1+\epsilon_1, w_2, \ldots, w_N) \approx L(w_1, w_2, \ldots, w_N) + \epsilon_1 \frac{d}{dw_1} L(w_1, w_2, \ldots, w_N).\]

我们将固定其他变量时对单个变量的导数称为偏导数，并使用符号\(\frac{\partial}{\partial w_1}\)表示在(22.4.1)中的导数。

现在，让我们把\(w_2\)稍微改变为\(w_2 + \epsilon_2\)

(22.4.2)\[\begin{split}\begin{aligned} L(w_1+\epsilon_1, w_2+\epsilon_2, \ldots, w_N) & \approx L(w_1, w_2+\epsilon_2, \ldots, w_N) + \epsilon_1 \frac{\partial}{\partial w_1} L(w_1, w_2+\epsilon_2, \ldots, w_N+\epsilon_N) \\ & \approx L(w_1, w_2, \ldots, w_N) \\ & \quad + \epsilon_2\frac{\partial}{\partial w_2} L(w_1, w_2, \ldots, w_N) \\ & \quad + \epsilon_1 \frac{\partial}{\partial w_1} L(w_1, w_2, \ldots, w_N) \\ & \quad + \epsilon_1\epsilon_2\frac{\partial}{\partial w_2}\frac{\partial}{\partial w_1} L(w_1, w_2, \ldots, w_N) \\ & \approx L(w_1, w_2, \ldots, w_N) \\ & \quad + \epsilon_2\frac{\partial}{\partial w_2} L(w_1, w_2, \ldots, w_N) \\ & \quad + \epsilon_1 \frac{\partial}{\partial w_1} L(w_1, w_2, \ldots, w_N). \end{aligned}\end{split}\]

我们再次使用了\(\epsilon_1\epsilon_2\)是一个可以舍弃的高阶项的思想，就像我们在上一节可以舍弃\(\epsilon^{2}\)一样，再加上我们在(22.4.1)中看到的。以这种方式继续，我们可以写出

(22.4.3)\[L(w_1+\epsilon_1, w_2+\epsilon_2, \ldots, w_N+\epsilon_N) \approx L(w_1, w_2, \ldots, w_N) + \sum_i \epsilon_i \frac{\partial}{\partial w_i} L(w_1, w_2, \ldots, w_N).\]

这可能看起来很乱，但我们可以通过注意到右边的和完全像一个点积来使其更熟悉，所以如果我们让

(22.4.4)\[\boldsymbol{\epsilon} = [\epsilon_1, \ldots, \epsilon_N]^\top \; \textrm{and} \; \nabla_{\mathbf{x}} L = \left[\frac{\partial L}{\partial x_1}, \ldots, \frac{\partial L}{\partial x_N}\right]^\top,\]

那么

(22.4.5)\[L(\mathbf{w} + \boldsymbol{\epsilon}) \approx L(\mathbf{w}) + \boldsymbol{\epsilon}\cdot \nabla_{\mathbf{w}} L(\mathbf{w}).\]

我们将向量\(\nabla_{\mathbf{w}} L\)称为\(L\)的梯度。

方程(22.4.5)值得我们花点时间思考。它的格式与我们在单维情况下遇到的完全相同，只是我们把所有东西都转换成了向量和点积。它允许我们近似地告诉函数\(L\)在给定任何输入扰动时会如何变化。正如我们将在下一节看到的，这将为我们提供一个重要的工具，从几何上理解我们如何利用梯度中包含的信息进行学习。

但首先，让我们通过一个例子来看看这个近似是如何工作的。假设我们正在处理函数

(22.4.6)\[f(x, y) = \log(e^x + e^y) \textrm{ 梯度为 } \nabla f (x, y) = \left[\frac{e^x}{e^x+e^y}, \frac{e^y}{e^x+e^y}\right].\]

如果我们看一个点，比如\((0, \log(2))\)，我们看到

(22.4.7)\[f(x, y) = \log(3) \textrm{ 梯度为 } \nabla f (x, y) = \left[\frac{1}{3}, \frac{2}{3}\right].\]

因此，如果我们想在\((\epsilon_1, \log(2) + \epsilon_2)\)处近似\(f\)，我们看到我们应该有(22.4.5)的具体实例

(22.4.8)\[f(\epsilon_1, \log(2) + \epsilon_2) \approx \log(3) + \frac{1}{3}\epsilon_1 + \frac{2}{3}\epsilon_2.\]

我们可以在代码中测试一下，看看这个近似有多好。

pytorchmxnettensorflow
%matplotlib inline
import numpy as np
import torch
from IPython import display
from mpl_toolkits import mplot3d
from d2l import torch as d2l


def f(x, y):
    return torch.log(torch.exp(x) + torch.exp(y))
def grad_f(x, y):
    return torch.tensor([torch.exp(x) / (torch.exp(x) + torch.exp(y)),
                     torch.exp(y) / (torch.exp(x) + torch.exp(y))])

epsilon = torch.tensor([0.01, -0.03])
grad_approx = f(torch.tensor([0.]), torch.log(
    torch.tensor([2.]))) + epsilon.dot(
    grad_f(torch.tensor([0.]), torch.log(torch.tensor(2.))))
true_value = f(torch.tensor([0.]) + epsilon[0], torch.log(
    torch.tensor([2.])) + epsilon[1])
f'approximation: {grad_approx}, true Value: {true_value}'

'approximation: tensor([1.0819]), true Value: tensor([1.0821])'

%matplotlib inline
from IPython import display
from mpl_toolkits import mplot3d
from mxnet import autograd, np, npx
from d2l import mxnet as d2l

npx.set_np()

def f(x, y):
    return np.log(np.exp(x) + np.exp(y))
def grad_f(x, y):
    return np.array([np.exp(x) / (np.exp(x) + np.exp(y)),
                     np.exp(y) / (np.exp(x) + np.exp(y))])

epsilon = np.array([0.01, -0.03])
grad_approx = f(0, np.log(2)) + epsilon.dot(grad_f(0, np.log(2)))
true_value = f(0 + epsilon[0], np.log(2) + epsilon[1])
f'approximation: {grad_approx}, true Value: {true_value}'

[21:56:31] ../src/storage/storage.cc:196: Using Pooled (Naive) StorageManager for CPU

'approximation: 1.0819456577301025, true Value: 1.0821242332458496'

%matplotlib inline
import numpy as np
import tensorflow as tf
from IPython import display
from mpl_toolkits import mplot3d
from d2l import tensorflow as d2l


def f(x, y):
    return tf.math.log(tf.exp(x) + tf.exp(y))
def grad_f(x, y):
    return tf.constant([(tf.exp(x) / (tf.exp(x) + tf.exp(y))).numpy(),
                        (tf.exp(y) / (tf.exp(x) + tf.exp(y))).numpy()])

epsilon = tf.constant([0.01, -0.03])
grad_approx = f(tf.constant([0.]), tf.math.log(
    tf.constant([2.]))) + tf.tensordot(
    epsilon, grad_f(tf.constant([0.]), tf.math.log(tf.constant(2.))), axes=1)
true_value = f(tf.constant([0.]) + epsilon[0], tf.math.log(
    tf.constant([2.])) + epsilon[1])
f'approximation: {grad_approx}, true Value: {true_value}'

'approximation: [1.0819457], true Value: [1.0821242]'

22.4.2. 梯度和梯度下降的几何学

再次考虑(22.4.5)中的表达式

(22.4.9)\[L(\mathbf{w} + \boldsymbol{\epsilon}) \approx L(\mathbf{w}) + \boldsymbol{\epsilon}\cdot \nabla_{\mathbf{w}} L(\mathbf{w}).\]

假设我想用这个来帮助最小化我们的损失\(L\)。让我们先从几何上理解在第 2.5 节中首次描述的梯度下降算法。我们将做以下事情

  1. 从随机选择初始参数\(\mathbf{w}\)开始。

  2. 找到使\(L\)在\(\mathbf{w}\)处下降最快的方向\(\mathbf{v}\)。

  3. 朝那个方向迈出一小步：\(\mathbf{w} \rightarrow \mathbf{w} + \epsilon\mathbf{v}\)。

  4. 重复。

我们唯一不确切知道如何做的是在第二步计算向量\(\mathbf{v}\)。我们将这样的方向称为最陡下降方向。使用第 22.1 节中点积的几何理解，我们看到我们可以将(22.4.5)重写为

(22.4.10)\[L(\mathbf{w} + \mathbf{v}) \approx L(\mathbf{w}) + \mathbf{v}\cdot \nabla_{\mathbf{w}} L(\mathbf{w}) = L(\mathbf{w}) + \|\nabla_{\mathbf{w}} L(\mathbf{w})\|\cos(\theta).\]

请注意，为了方便，我们已将我们的方向设为单位长度，并用\(\theta\)表示\(\mathbf{v}\)和\(\nabla_{\mathbf{w}} L(\mathbf{w})\)之间的角度。如果我们想找到使\(L\)下降最快的方向，我们希望使这个表达式尽可能为负。我们选择的方向进入这个方程的唯一方式是通过\(\cos(\theta)\)，因此我们希望使这个余弦尽可能为负。现在，回想一下余弦的形状，我们可以通过使\(\cos(\theta) = -1\)，或等价地使梯度和我们选择的方向之间的角度为\(\pi\)弧度，或等价地\(180\)度，来使其尽可能为负。实现这一点的唯一方法是朝完全相反的方向前进：选择\(\mathbf{v}\)指向与\(\nabla_{\mathbf{w}} L(\mathbf{w})\)完全相反的方向！

这给我们带来了机器学习中最重要的数学概念之一：最陡下降方向指向\(-\nabla_{\mathbf{w}}L(\mathbf{w})\)的方向。因此，我们的非正式算法可以重写如下。

  1. 从随机选择初始参数\(\mathbf{w}\)开始。

  2. 计算\(\nabla_{\mathbf{w}} L(\mathbf{w})\)。

  3. 朝该方向的相反方向迈出一小步：\(\mathbf{w} \leftarrow \mathbf{w} - \epsilon\nabla_{\mathbf{w}} L(\mathbf{w})\)。

  4. 重复。

这个基本算法已经被许多研究者以多种方式修改和调整，但其核心概念在所有这些方法中都保持不变。利用梯度找到使损失下降最快的方向，并更新参数以朝该方向迈出一步。

22.4.3. 关于数学优化的注记

在本书中，我们完全专注于数值优化技术，因为在深度学习设置中我们遇到的所有函数都过于复杂，无法显式地最小化。

然而，考虑一下我们上面获得的几何理解告诉我们如何直接优化函数，这是一个有用的练习。

假设我们希望找到最小化某个函数\(L(\mathbf{x})\)的值\(\mathbf{x}_0\)。再假设有人给了我们一个值，并告诉我们这就是最小化\(L\)的值。我们有什么可以检查的，看看他们的答案是否合理吗？

再次考虑(22.4.5)

(22.4.11)\[L(\mathbf{x}_0 + \boldsymbol{\epsilon}) \approx L(\mathbf{x}_0) + \boldsymbol{\epsilon}\cdot \nabla_{\mathbf{x}} L(\mathbf{x}_0).\]

如果梯度不为零，我们知道我们可以朝着\(-\epsilon \nabla_{\mathbf{x}} L(\mathbf{x}_0)\)的方向迈出一步，以找到一个更小的\(L\)值。因此，如果我们真的处在一个最小值，这种情况就不可能发生！我们可以得出结论，如果\(\mathbf{x}_0\)是一个最小值，那么\(\nabla_{\mathbf{x}} L(\mathbf{x}_0) = 0\)。我们称梯度为\(\nabla_{\mathbf{x}} L(\mathbf{x}_0) = 0\)的点为临界点。

这很好，因为在一些罕见的情况下，我们可以显式地找到所有梯度为零的点，并找到值最小的那个点。

举一个具体的例子，考虑函数

(22.4.12)\[f(x) = 3x^4 - 4x^3 -12x^2.\]

这个函数的导数是

(22.4.13)\[\frac{df}{dx} = 12x^3 - 12x^2 -24x = 12x(x-2)(x+1).\]

最小值的可能位置只有\(x = -1, 0, 2\)，函数在这些点的值分别为\(-5,0, -32\)，因此我们可以得出结论，当\(x = 2\)时，我们的函数达到最小值。一个快速的绘图证实了这一点。

pytorchmxnettensorflow
x = torch.arange(-2, 3, 0.01)
f = (3 * x**4) - (4 * x**3) - (12 * x**2)

d2l.plot(x, f, 'x', 'f(x)')
../_images/output_multivariable-calculus_bf678a_15_0.svg
x = np.arange(-2, 3, 0.01)
f = (3 * x**4) - (4 * x**3) - (12 * x**2)

d2l.plot(x, f, 'x', 'f(x)')
../_images/output_multivariable-calculus_bf678a_18_0.svg
x = tf.range(-2, 3, 0.01)
f = (3 * x**4) - (4 * x**3) - (12 * x**2)

d2l.plot(x, f, 'x', 'f(x)')
../_images/output_multivariable-calculus_bf678a_21_0.svg

这强调了在进行理论或数值工作时需要知道的一个重要事实：我们能够最小化（或最大化）一个函数的唯一可能点是梯度等于零的点，然而，并非每个梯度为零的点都是真正的全局最小值（或最大值）。

22.4.4. 多变量链式法则

假设我们有一个四个变量（\(w, x, y\)和\(z\)）的函数，我们可以通过组合多个项来构成它

(22.4.14)\[\begin{split}\begin{aligned}f(u, v) & = (u+v)^{2} \\u(a, b) & = (a+b)^{2}, \qquad v(a, b) = (a-b)^{2}, \\a(w, x, y, z) & = (w+x+y+z)^{2}, \qquad b(w, x, y, z) = (w+x-y-z)^2.\end{aligned}\end{split}\]

在处理神经网络时，这样的方程链很常见，因此理解如何计算这类函数的梯度是关键。如果我们看一下哪些变量直接相关，我们可以在图 22.4.1中看到这种联系的视觉提示。

../_images/chain-net1.svg

图 22.4.1 上图中的函数关系，其中节点代表值，边显示函数依赖。

没有什么能阻止我们直接将(22.4.14)中的所有内容组合起来，然后写出

(22.4.15)\[f(w, x, y, z) = \left(\left((w+x+y+z)^2+(w+x-y-z)^2\right)^2+\left((w+x+y+z)^2-(w+x-y-z)^2\right)^2\right)^2.\]

然后我们可以只用单变量导数来求导，但如果我们那样做，我们很快就会发现自己被一大堆项淹没，其中许多是重复的！事实上，我们可以看到，例如

(22.4.16)\[\begin{split}\begin{aligned} \frac{\partial f}{\partial w} & = 2 \left(2 \left(2 (w + x + y + z) - 2 (w + x - y - z)\right) \left((w + x + y + z)^{2}- (w + x - y - z)^{2}\right) + \right.\\ & \left. \quad 2 \left(2 (w + x - y - z) + 2 (w + x + y + z)\right) \left((w + x - y - z)^{2}+ (w + x + y + z)^{2}\right)\right) \times \\ & \quad \left(\left((w + x + y + z)^{2}- (w + x - y - z)^2\right)^{2}+ \left((w + x - y - z)^{2}+ (w + x + y + z)^{2}\right)^{2}\right). \end{aligned}\end{split}\]

如果那时我们还想计算\(\frac{\partial f}{\partial x}\)，我们会得到一个类似的方程，再次有许多重复的项，以及两个导数之间许多共享的重复项。这代表了大量的重复工作，如果我们需要用这种方式计算导数，整个深度学习革命可能在开始之前就已经停滞了！

让我们分解这个问题。我们将从尝试理解当我们改变\(a\)时\(f\)如何变化开始，基本上假设\(w, x, y\)和\(z\)都不存在。我们将像我们第一次处理梯度时那样推理。让我们取\(a\)并给它加上一个小的量\(\epsilon\)。

(22.4.17)\[\begin{split}\begin{aligned} & f(u(a+\epsilon, b), v(a+\epsilon, b)) \\ \approx & f\left(u(a, b) + \epsilon\frac{\partial u}{\partial a}(a, b), v(a, b) + \epsilon\frac{\partial v}{\partial a}(a, b)\right) \\ \approx & f(u(a, b), v(a, b)) + \epsilon\left[\frac{\partial f}{\partial u}(u(a, b), v(a, b))\frac{\partial u}{\partial a}(a, b) + \frac{\partial f}{\partial v}(u(a, b), v(a, b))\frac{\partial v}{\partial a}(a, b)\right]. \end{aligned}\end{split}\]

第一行遵循偏导数的定义，第二行遵循梯度的定义。精确追踪每个导数的求值位置在符号上是繁琐的，就像在表达式\(\frac{\partial f}{\partial u}(u(a, b), v(a, b))\)中一样，所以我们经常将其缩写为更易记的

(22.4.18)\[\frac{\partial f}{\partial a} = \frac{\partial f}{\partial u}\frac{\partial u}{\partial a}+\frac{\partial f}{\partial v}\frac{\partial v}{\partial a}.\]

思考这个过程的意义是很有用的。我们试图理解一个形式为\(f(u(a, b), v(a, b))\)的函数如何随着\(a\)的变化而改变其值。这可以通过两条路径发生：一条是\(a \rightarrow u \rightarrow f\)，另一条是\(a \rightarrow v \rightarrow f\)。我们可以通过链式法则计算这两个贡献：分别是\(\frac{\partial w}{\partial u} \cdot \frac{\partial u}{\partial x}\)和\(\frac{\partial w}{\partial v} \cdot \frac{\partial v}{\partial x}\)，然后相加。

想象我们有一个不同的函数网络，其中右边的函数依赖于左边与之相连的函数，如图 22.4.2所示。

../_images/chain-net2.svg

图 22.4.2 另一个更微妙的链式法则例子。

要计算像\(\frac{\partial f}{\partial y}\)这样的东西，我们需要对从\(y\)到\(f\)的所有路径（在这种情况下是\(3\)条）进行求和，得到

(22.4.19)\[\frac{\partial f}{\partial y} = \frac{\partial f}{\partial a} \frac{\partial a}{\partial u} \frac{\partial u}{\partial y} + \frac{\partial f}{\partial u} \frac{\partial u}{\partial y} + \frac{\partial f}{\partial b} \frac{\partial b}{\partial v} \frac{\partial v}{\partial y}.\]

以这种方式理解链式法则，在尝试理解梯度如何通过网络流动，以及为什么各种架构选择（如LSTM（第 10.1 节）或残差层（第 8.6 节））可以通过控制梯度流来帮助塑造学习过程时，将会大有裨益。

22.4.5. 反向传播算法

让我们回到上一节(22.4.14)的例子，其中

(22.4.20)\[\begin{split}\begin{aligned} f(u, v) & = (u+v)^{2} \\ u(a, b) & = (a+b)^{2}, \qquad v(a, b) = (a-b)^{2}, \\ a(w, x, y, z) & = (w+x+y+z)^{2}, \qquad b(w, x, y, z) = (w+x-y-z)^2. \end{aligned}\end{split}\]

如果我们想计算比如\(\frac{\partial f}{\partial w}\)，我们可以应用多变量链式法则看到

(22.4.21)\[\begin{split}\begin{aligned} \frac{\partial f}{\partial w} & = \frac{\partial f}{\partial u}\frac{\partial u}{\partial w} + \frac{\partial f}{\partial v}\frac{\partial v}{\partial w}, \\ \frac{\partial u}{\partial w} & = \frac{\partial u}{\partial a}\frac{\partial a}{\partial w}+\frac{\partial u}{\partial b}\frac{\partial b}{\partial w}, \\ \frac{\partial v}{\partial w} & = \frac{\partial v}{\partial a}\frac{\partial a}{\partial w}+\frac{\partial v}{\partial b}\frac{\partial b}{\partial w}. \end{aligned}\end{split}\]

让我们试着用这种分解来计算\(\frac{\partial f}{\partial w}\)。注意，这里我们只需要各种单步偏导数

(22.4.22)\[\begin{split}\begin{aligned} \frac{\partial f}{\partial u} = 2(u+v), & \quad\frac{\partial f}{\partial v} = 2(u+v), \\ \frac{\partial u}{\partial a} = 2(a+b), & \quad\frac{\partial u}{\partial b} = 2(a+b), \\ \frac{\partial v}{\partial a} = 2(a-b), & \quad\frac{\partial v}{\partial b} = -2(a-b), \\ \frac{\partial a}{\partial w} = 2(w+x+y+z), & \quad\frac{\partial b}{\partial w} = 2(w+x-y-z). \end{aligned}\end{split}\]

如果我们将此写成代码，它会变成一个相当易于管理的表达式。

pytorchmxnettensorflow
# Compute the value of the function from inputs to outputs
w, x, y, z = -1, 0, -2, 1
a, b = (w + x + y + z)**2, (w + x - y - z)**2
u, v = (a + b)**2, (a - b)**2
f = (u + v)**2
print(f'    f at {w}, {x}, {y}, {z} is {f}')

# Compute the single step partials
df_du, df_dv = 2*(u + v), 2*(u + v)
du_da, du_db, dv_da, dv_db = 2*(a + b), 2*(a + b), 2*(a - b), -2*(a - b)
da_dw, db_dw = 2*(w + x + y + z), 2*(w + x - y - z)

# Compute the final result from inputs to outputs
du_dw, dv_dw = du_da*da_dw + du_db*db_dw, dv_da*da_dw + dv_db*db_dw
df_dw = df_du*du_dw + df_dv*dv_dw
print(f'df/dw at {w}, {x}, {y}, {z} is {df_dw}')

    f at -1, 0, -2, 1 is 1024
df/dw at -1, 0, -2, 1 is -4096

# Compute the value of the function from inputs to outputs
w, x, y, z = -1, 0, -2, 1
a, b = (w + x + y + z)**2, (w + x - y - z)**2
u, v = (a + b)**2, (a - b)**2
f = (u + v)**2
print(f'    f at {w}, {x}, {y}, {z} is {f}')

# Compute the single step partials
df_du, df_dv = 2*(u + v), 2*(u + v)
du_da, du_db, dv_da, dv_db = 2*(a + b), 2*(a + b), 2*(a - b), -2*(a - b)
da_dw, db_dw = 2*(w + x + y + z), 2*(w + x - y - z)

# Compute the final result from inputs to outputs
du_dw, dv_dw = du_da*da_dw + du_db*db_dw, dv_da*da_dw + dv_db*db_dw
df_dw = df_du*du_dw + df_dv*dv_dw
print(f'df/dw at {w}, {x}, {y}, {z} is {df_dw}')

    f at -1, 0, -2, 1 is 1024
df/dw at -1, 0, -2, 1 is -4096

# Compute the value of the function from inputs to outputs
w, x, y, z = -1, 0, -2, 1
a, b = (w + x + y + z)**2, (w + x - y - z)**2
u, v = (a + b)**2, (a - b)**2
f = (u + v)**2
print(f'    f at {w}, {x}, {y}, {z} is {f}')

# Compute the single step partials
df_du, df_dv = 2*(u + v), 2*(u + v)
du_da, du_db, dv_da, dv_db = 2*(a + b), 2*(a + b), 2*(a - b), -2*(a - b)
da_dw, db_dw = 2*(w + x + y + z), 2*(w + x - y - z)

# Compute the final result from inputs to outputs
du_dw, dv_dw = du_da*da_dw + du_db*db_dw, dv_da*da_dw + dv_db*db_dw
df_dw = df_du*du_dw + df_dv*dv_dw
print(f'df/dw at {w}, {x}, {y}, {z} is {df_dw}')

    f at -1, 0, -2, 1 is 1024
df/dw at -1, 0, -2, 1 is -4096

然而，请注意，这仍然不便于计算像\(\frac{\partial f}{\partial x}\)这样的东西。原因是我们选择应用链式法则的方式。如果我们看上面我们所做的，我们在可能的情况下总是将\(\partial w\)保留在分母中。通过这种方式，我们选择应用链式法则来看\(w\)如何改变其他所有变量。如果这是我们想要的，这将是一个好主意。然而，回想一下我们来自深度学习的动机：我们想看看每个参数如何改变损失。本质上，我们希望在任何可能的情况下都将\(\partial f\)保留在分子中来应用链式法则！

更明确地说，注意我们可以写出

(22.4.23)\[\begin{split}\begin{aligned} \frac{\partial f}{\partial w} & = \frac{\partial f}{\partial a}\frac{\partial a}{\partial w} + \frac{\partial f}{\partial b}\frac{\partial b}{\partial w}, \\ \frac{\partial f}{\partial a} & = \frac{\partial f}{\partial u}\frac{\partial u}{\partial a}+\frac{\partial f}{\partial v}\frac{\partial v}{\partial a}, \\ \frac{\partial f}{\partial b} & = \frac{\partial f}{\partial u}\frac{\partial u}{\partial b}+\frac{\partial f}{\partial v}\frac{\partial v}{\partial b}. \end{aligned}\end{split}\]

注意，链式法则的这种应用让我们明确地计算了\(\frac{\partial f}{\partial u}, \frac{\partial f}{\partial v}, \frac{\partial f}{\partial a}, \frac{\partial f}{\partial b}, \; \textrm{和} \; \frac{\partial f}{\partial w}\)。没有什么能阻止我们同时包含方程

(22.4.24)\[\begin{split}\begin{aligned} \frac{\partial f}{\partial x} & = \frac{\partial f}{\partial a}\frac{\partial a}{\partial x} + \frac{\partial f}{\partial b}\frac{\partial b}{\partial x}, \\ \frac{\partial f}{\partial y} & = \frac{\partial f}{\partial a}\frac{\partial a}{\partial y}+\frac{\partial f}{\partial b}\frac{\partial b}{\partial y}, \\ \frac{\partial f}{\partial z} & = \frac{\partial f}{\partial a}\frac{\partial a}{\partial z}+\frac{\partial f}{\partial b}\frac{\partial b}{\partial z}. \end{aligned}\end{split}\]

然后跟踪当我们改变整个网络中任何节点时\(f\)的变化情况。让我们来实现它。

pytorchmxnettensorflow
# Compute the value of the function from inputs to outputs
w, x, y, z = -1, 0, -2, 1
a, b = (w + x + y + z)**2, (w + x - y - z)**2
u, v = (a + b)**2, (a - b)**2
f = (u + v)**2
print(f'f at {w}, {x}, {y}, {z} is {f}')

# Compute the derivative using the decomposition above
# First compute the single step partials
df_du, df_dv = 2*(u + v), 2*(u + v)
du_da, du_db, dv_da, dv_db = 2*(a + b), 2*(a + b), 2*(a - b), -2*(a - b)
da_dw, db_dw = 2*(w + x + y + z), 2*(w + x - y - z)
da_dx, db_dx = 2*(w + x + y + z), 2*(w + x - y - z)
da_dy, db_dy = 2*(w + x + y + z), -2*(w + x - y - z)
da_dz, db_dz = 2*(w + x + y + z), -2*(w + x - y - z)

# Now compute how f changes when we change any value from output to input
df_da, df_db = df_du*du_da + df_dv*dv_da, df_du*du_db + df_dv*dv_db
df_dw, df_dx = df_da*da_dw + df_db*db_dw, df_da*da_dx + df_db*db_dx
df_dy, df_dz = df_da*da_dy + df_db*db_dy, df_da*da_dz + df_db*db_dz

print(f'df/dw at {w}, {x}, {y}, {z} is {df_dw}')
print(f'df/dx at {w}, {x}, {y}, {z} is {df_dx}')
print(f'df/dy at {w}, {x}, {y}, {z} is {df_dy}')
print(f'df/dz at {w}, {x}, {y}, {z} is {df_dz}')

f at -1, 0, -2, 1 is 1024
df/dw at -1, 0, -2, 1 is -4096
df/dx at -1, 0, -2, 1 is -4096
df/dy at -1, 0, -2, 1 is -4096
df/dz at -1, 0, -2, 1 is -4096

# Compute the value of the function from inputs to outputs
w, x, y, z = -1, 0, -2, 1
a, b = (w + x + y + z)**2, (w + x - y - z)**2
u, v = (a + b)**2, (a - b)**2
f = (u + v)**2
print(f'f at {w}, {x}, {y}, {z} is {f}')

# Compute the derivative using the decomposition above
# First compute the single step partials
df_du, df_dv = 2*(u + v), 2*(u + v)
du_da, du_db, dv_da, dv_db = 2*(a + b), 2*(a + b), 2*(a - b), -2*(a - b)
da_dw, db_dw = 2*(w + x + y + z), 2*(w + x - y - z)
da_dx, db_dx = 2*(w + x + y + z), 2*(w + x - y - z)
da_dy, db_dy = 2*(w + x + y + z), -2*(w + x - y - z)
da_dz, db_dz = 2*(w + x + y + z), -2*(w + x - y - z)

# Now compute how f changes when we change any value from output to input
df_da, df_db = df_du*du_da + df_dv*dv_da, df_du*du_db + df_dv*dv_db
df_dw, df_dx = df_da*da_dw + df_db*db_dw, df_da*da_dx + df_db*db_dx
df_dy, df_dz = df_da*da_dy + df_db*db_dy, df_da*da_dz + df_db*db_dz

print(f'df/dw at {w}, {x}, {y}, {z} is {df_dw}')
print(f'df/dx at {w}, {x}, {y}, {z} is {df_dx}')
print(f'df/dy at {w}, {x}, {y}, {z} is {df_dy}')
print(f'df/dz at {w}, {x}, {y}, {z} is {df_dz}')

f at -1, 0, -2, 1 is 1024
df/dw at -1, 0, -2, 1 is -4096
df/dx at -1, 0, -2, 1 is -4096
df/dy at -1, 0, -2, 1 is -4096
df/dz at -1, 0, -2, 1 is -4096

# Compute the value of the function from inputs to outputs
w, x, y, z = -1, 0, -2, 1
a, b = (w + x + y + z)**2, (w + x - y - z)**2
u, v = (a + b)**2, (a - b)**2
f = (u + v)**2
print(f'f at {w}, {x}, {y}, {z} is {f}')

# Compute the derivative using the decomposition above
# First compute the single step partials
df_du, df_dv = 2*(u + v), 2*(u + v)
du_da, du_db, dv_da, dv_db = 2*(a + b), 2*(a + b), 2*(a - b), -2*(a - b)
da_dw, db_dw = 2*(w + x + y + z), 2*(w + x - y - z)
da_dx, db_dx = 2*(w + x + y + z), 2*(w + x - y - z)
da_dy, db_dy = 2*(w + x + y + z), -2*(w + x - y - z)
da_dz, db_dz = 2*(w + x + y + z), -2*(w + x - y - z)

# Now compute how f changes when we change any value from output to input
df_da, df_db = df_du*du_da + df_dv*dv_da, df_du*du_db + df_dv*dv_db
df_dw, df_dx = df_da*da_dw + df_db*db_dw, df_da*da_dx + df_db*db_dx
df_dy, df_dz = df_da*da_dy + df_db*db_dy, df_da*da_dz + df_db*db_dz

print(f'df/dw at {w}, {x}, {y}, {z} is {df_dw}')
print(f'df/dx at {w}, {x}, {y}, {z} is {df_dx}')
print(f'df/dy at {w}, {x}, {y}, {z} is {df_dy}')
print(f'df/dz at {w}, {x}, {y}, {z} is {df_dz}')

f at -1, 0, -2, 1 is 1024
df/dw at -1, 0, -2, 1 is -4096
df/dx at -1, 0, -2, 1 is -4096
df/dy at -1, 0, -2, 1 is -4096
df/dz at -1, 0, -2, 1 is -4096

我们从\(f\)向输入方向计算导数，而不是从输入向前到输出（正如我们在上面的第一个代码片段中所做的），这一事实正是这个算法得名的原因：反向传播。注意有两个步骤：1. 计算函数的值，以及从前到后的单步偏导数。虽然上面没有这样做，但这可以组合成一个单一的前向传播。2. 从后到前计算\(f\)的梯度。我们称之为反向传播。

这正是每个深度学习算法所实现的，以允许一次性计算损失相对于网络中每个权重的梯度。我们有这样一种分解是一个惊人的事实。

为了看看如何封装这一点，让我们快速看一下这个例子。

pytorchmxnettensorflow
# Initialize as ndarrays, then attach gradients
w = torch.tensor([-1.], requires_grad=True)
x = torch.tensor([0.], requires_grad=True)
y = torch.tensor([-2.], requires_grad=True)
z = torch.tensor([1.], requires_grad=True)
# Do the computation like usual, tracking gradients
a, b = (w + x + y + z)**2, (w + x - y - z)**2
u, v = (a + b)**2, (a - b)**2
f = (u + v)**2

# Execute backward pass
f.backward()

print(f'df/dw at {w.data.item()}, {x.data.item()}, {y.data.item()}, '
      f'{z.data.item()} is {w.grad.data.item()}')
print(f'df/dx at {w.data.item()}, {x.data.item()}, {y.data.item()}, '
      f'{z.data.item()} is {x.grad.data.item()}')
print(f'df/dy at {w.data.item()}, {x.data.item()}, {y.data.item()}, '
      f'{z.data.item()} is {y.grad.data.item()}')
print(f'df/dz at {w.data.item()}, {x.data.item()}, {y.data.item()}, '
      f'{z.data.item()} is {z.grad.data.item()}')

df/dw at -1.0, 0.0, -2.0, 1.0 is -4096.0
df/dx at -1.0, 0.0, -2.0, 1.0 is -4096.0
df/dy at -1.0, 0.0, -2.0, 1.0 is -4096.0
df/dz at -1.0, 0.0, -2.0, 1.0 is -4096.0

# Initialize as ndarrays, then attach gradients
w, x, y, z = np.array(-1), np.array(0), np.array(-2), np.array(1)

w.attach_grad()
x.attach_grad()
y.attach_grad()
z.attach_grad()

# Do the computation like usual, tracking gradients
with autograd.record():
    a, b = (w + x + y + z)**2, (w + x - y - z)**2
    u, v = (a + b)**2, (a - b)**2
    f = (u + v)**2

# Execute backward pass
f.backward()

print(f'df/dw at {w}, {x}, {y}, {z} is {w.grad}')
print(f'df/dx at {w}, {x}, {y}, {z} is {x.grad}')
print(f'df/dy at {w}, {x}, {y}, {z} is {y.grad}')
print(f'df/dz at {w}, {x}, {y}, {z} is {z.grad}')

df/dw at -1.0, 0.0, -2.0, 1.0 is -4096.0
df/dx at -1.0, 0.0, -2.0, 1.0 is -4096.0
df/dy at -1.0, 0.0, -2.0, 1.0 is -4096.0
df/dz at -1.0, 0.0, -2.0, 1.0 is -4096.0
[21:56:31] ../src/base.cc:48: GPU context requested, but no GPUs found.

# Initialize as ndarrays, then attach gradients
w = tf.Variable(tf.constant([-1.]))
x = tf.Variable(tf.constant([0.]))
y = tf.Variable(tf.constant([-2.]))
z = tf.Variable(tf.constant([1.]))
# Do the computation like usual, tracking gradients
with tf.GradientTape(persistent=True) as t:
    a, b = (w + x + y + z)**2, (w + x - y - z)**2
    u, v = (a + b)**2, (a - b)**2
    f = (u + v)**2

# Execute backward pass
w_grad = t.gradient(f, w).numpy()
x_grad = t.gradient(f, x).numpy()
y_grad = t.gradient(f, y).numpy()
z_grad = t.gradient(f, z).numpy()

print(f'df/dw at {w.numpy()}, {x.numpy()}, {y.numpy()}, '
      f'{z.numpy()} is {w_grad}')
print(f'df/dx at {w.numpy()}, {x.numpy()}, {y.numpy()}, '
      f'{z.numpy()} is {x_grad}')
print(f'df/dy at {w.numpy()}, {x.numpy()}, {y.numpy()}, '
      f'{z.numpy()} is {y_grad}')
print(f'df/dz at {w.numpy()}, {x.numpy()}, {y.numpy()}, '
      f'{z.numpy()} is {z_grad}')

df/dw at [-1.], [0.], [-2.], [1.] is [-4096.]
df/dx at [-1.], [0.], [-2.], [1.] is [-4096.]
df/dy at [-1.], [0.], [-2.], [1.] is [-4096.]
df/dz at [-1.], [0.], [-2.], [1.] is [-4096.]

我们上面所做的一切都可以通过调用f.backwards()自动完成。

22.4.6. Hessian矩阵

与单变量微积分一样，考虑高阶导数是有用的，以便掌握如何获得比单独使用梯度更好的函数近似。

在处理多变量函数的高阶导数时，会立即遇到一个问题，那就是它们的数量很多。如果我们有一个\(n\)个变量的函数\(f(x_1, \ldots, x_n)\)，那么我们可以取\(n^{2}\)个二阶导数，即对于任何\(i\)和\(j\)的选择

(22.4.25)\[\frac{d^2f}{dx_idx_j} = \frac{d}{dx_i}\left(\frac{d}{dx_j}f\right).\]

这通常被组装成一个称为Hessian矩阵的矩阵

(22.4.26)\[\begin{split}\mathbf{H}_f = \begin{bmatrix} \frac{d^2f}{dx_1dx_1} & \cdots & \frac{d^2f}{dx_1dx_n} \\ \vdots & \ddots & \vdots \\ \frac{d^2f}{dx_ndx_1} & \cdots & \frac{d^2f}{dx_ndx_n} \\ \end{bmatrix}.\end{split}\]

这个矩阵的并非所有条目都是独立的。实际上，我们可以证明，只要两个混合偏导数（对多个变量的偏导数）都存在且连续，我们就可以说对于任何\(i\)和\(j\)，

(22.4.27)\[\frac{d^2f}{dx_idx_j} = \frac{d^2f}{dx_jdx_i}.\]

这是通过首先在\(x_i\)方向上扰动函数，然后在\(x_j\)方向上扰动，然后将该结果与如果我们首先扰动\(x_j\)然后扰动\(x_i\)会发生什么进行比较，并且知道这两种顺序都会导致\(f\)的输出发生相同的最终变化。

与单变量一样，我们可以使用这些导数来更好地了解函数在某点附近的行为。特别是，我们可以用它来找到在点\(\mathbf{x}_0\)附近最佳拟合的二次函数，正如我们在单变量中看到的那样。

让我们看一个例子。假设\(f(x_1, x_2) = a + b_1x_1 + b_2x_2 + c_{11}x_1^{2} + c_{12}x_1x_2 + c_{22}x_2^{2}\)。这是两个变量的二次函数的一般形式。如果我们看一下函数的值、它的梯度和它的Hessian矩阵(22.4.26)，都在零点处

(22.4.28)\[\begin{split}\begin{aligned} f(0,0) & = a, \\ \nabla f (0,0) & = \begin{bmatrix}b_1 \\ b_2\end{bmatrix}, \\ \mathbf{H} f (0,0) & = \begin{bmatrix}2 c_{11} & c_{12} \\ c_{12} & 2c_{22}\end{bmatrix}, \end{aligned}\end{split}\]

我们可以通过说

(22.4.29)\[f(\mathbf{x}) = f(0) + \nabla f (0) \cdot \mathbf{x} + \frac{1}{2}\mathbf{x}^\top \mathbf{H} f (0) \mathbf{x}.\]

一般来说，如果我们在任何点\(\mathbf{x}_0\)计算这个展开式，我们看到

(22.4.30)\[f(\mathbf{x}) = f(\mathbf{x}_0) + \nabla f (\mathbf{x}_0) \cdot (\mathbf{x}-\mathbf{x}_0) + \frac{1}{2}(\mathbf{x}-\mathbf{x}_0)^\top \mathbf{H} f (\mathbf{x}_0) (\mathbf{x}-\mathbf{x}_0).\]

这适用于任何维度的输入，并为任何函数在某点提供了最佳的二次近似。举个例子，让我们绘制函数

(22.4.31)\[f(x, y) = xe^{-x^2-y^2}.\]

可以计算出梯度和Hessian矩阵为

(22.4.32)\[\begin{split}\nabla f(x, y) = e^{-x^2-y^2}\begin{pmatrix}1-2x^2 \\ -2xy\end{pmatrix} \; \textrm{和} \; \mathbf{H}f(x, y) = e^{-x^2-y^2}\begin{pmatrix} 4x^3 - 6x & 4x^2y - 2y \\ 4x^2y-2y &4xy^2-2x\end{pmatrix}.\end{split}\]

因此，经过一些代数运算，可以看到在\([-1,0]^\top\)处的近似二次函数是

(22.4.33)\[f(x, y) \approx e^{-1}\left(-1 - (x+1) +(x+1)^2+y^2\right).\]
pytorchmxnettensorflow
# Construct grid and compute function
x, y = torch.meshgrid(torch.linspace(-2, 2, 101),
                   torch.linspace(-2, 2, 101))

z = x*torch.exp(- x**2 - y**2)

# Compute approximating quadratic with gradient and Hessian at (1, 0)
w = torch.exp(torch.tensor([-1.]))*(-1 - (x + 1) + 2 * (x + 1)**2 + 2 * y**2)

# Plot function
ax = d2l.plt.figure().add_subplot(111, projection='3d')
ax.plot_wireframe(x.numpy(), y.numpy(), z.numpy(),
                  **{'rstride': 10, 'cstride': 10})
ax.plot_wireframe(x.numpy(), y.numpy(), w.numpy(),
                  **{'rstride': 10, 'cstride': 10}, color='purple')
d2l.plt.xlabel('x')
d2l.plt.ylabel('y')
d2l.set_figsize()
ax.set_xlim(-2, 2)
ax.set_ylim(-2, 2)
ax.set_zlim(-1, 1)
ax.dist = 12
../_images/output_multivariable-calculus_bf678a_63_0.svg
# Construct grid and compute function
x, y = np.meshgrid(np.linspace(-2, 2, 101),
                   np.linspace(-2, 2, 101), indexing='ij')
z = x*np.exp(- x**2 - y**2)

# Compute approximating quadratic with gradient and Hessian at (1, 0)
w = np.exp(-1)*(-1 - (x + 1) + (x + 1)**2 + y**2)

# Plot function
ax = d2l.plt.figure().add_subplot(111, projection='3d')
ax.plot_wireframe(x.asnumpy(), y.asnumpy(), z.asnumpy(),
                  **{'rstride': 10, 'cstride': 10})
ax.plot_wireframe(x.asnumpy(), y.asnumpy(), w.asnumpy(),
                  **{'rstride': 10, 'cstride': 10}, color='purple')
d2l.plt.xlabel('x')
d2l.plt.ylabel('y')
d2l.set_figsize()
ax.set_xlim(-2, 2)
ax.set_ylim(-2, 2)
ax.set_zlim(-1, 1)
ax.dist = 12
../_images/output_multivariable-calculus_bf678a_66_0.svg
# Construct grid and compute function
x, y = tf.meshgrid(tf.linspace(-2., 2., 101),
                   tf.linspace(-2., 2., 101))

z = x*tf.exp(- x**2 - y**2)

# Compute approximating quadratic with gradient and Hessian at (1, 0)
w = tf.exp(tf.constant([-1.]))*(-1 - (x + 1) + 2 * (x + 1)**2 + 2 * y**2)

# Plot function
ax = d2l.plt.figure().add_subplot(111, projection='3d')
ax.plot_wireframe(x.numpy(), y.numpy(), z.numpy(),
                  **{'rstride': 10, 'cstride': 10})
ax.plot_wireframe(x.numpy(), y.numpy(), w.numpy(),
                  **{'rstride': 10, 'cstride': 10}, color='purple')
d2l.plt.xlabel('x')
d2l.plt.ylabel('y')
d2l.set_figsize()
ax.set_xlim(-2, 2)
ax.set_ylim(-2, 2)
ax.set_zlim(-1, 1)
ax.dist = 12
../_images/output_multivariable-calculus_bf678a_69_0.svg

这构成了第 12.3 节中讨论的牛顿算法的基础，在该算法中，我们迭代地执行数值优化，找到最佳拟合的二次函数，然后精确地最小化该二次函数。

22.4.7. 一点矩阵微积分

涉及矩阵的函数导数结果特别简洁。本节在符号上可能变得繁重，因此在初次阅读时可以跳过，但了解涉及常见矩阵运算的函数导数通常比最初预期的要简洁得多，这很有用，特别是考虑到矩阵运算在深度学习应用中的核心地位。

让我们从一个例子开始。假设我们有一个固定的列向量\(\boldsymbol{\beta}\)，我们想取乘积函数\(f(\mathbf{x}) = \boldsymbol{\beta}^\top\mathbf{x}\)，并理解当我们改变\(\mathbf{x}\)时点积如何变化。

在处理ML中的矩阵导数时，一种有用的符号表示法叫做分母布局矩阵导数，其中我们将偏导数组织成与微分分母中任何向量、矩阵或张量相同的形状。在这种情况下，我们将写成

(22.4.34)\[\begin{split}\frac{df}{d\mathbf{x}} = \begin{bmatrix} \frac{df}{dx_1} \\ \vdots \\ \frac{df}{dx_n} \end{bmatrix},\end{split}\]

这里我们匹配了列向量\(\mathbf{x}\)的形状。

如果我们将我们的函数写成分量形式，即

(22.4.35)\[f(\mathbf{x}) = \sum_{i = 1}^{n} \beta_ix_i = \beta_1x_1 + \cdots + \beta_nx_n.\]

现在如果我们对\(\beta_1\)取偏导数，注意除了第一项外，其他项都为零，第一项只是\(x_1\)乘以\(\beta_1\)，所以我们得到

(22.4.36)\[\frac{df}{dx_1} = \beta_1,\]

或更一般地

(22.4.37)\[\frac{df}{dx_i} = \beta_i.\]

我们现在可以将其重新组合成一个矩阵，看到

(22.4.38)\[\begin{split}\frac{df}{d\mathbf{x}} = \begin{bmatrix} \frac{df}{dx_1} \\ \vdots \\ \frac{df}{dx_n} \end{bmatrix} = \begin{bmatrix} \beta_1 \\ \vdots \\ \beta_n \end{bmatrix} = \boldsymbol{\beta}.\end{split}\]

这说明了我们将在本节中经常遇到的关于矩阵微积分的几个因素

  • 首先，计算会变得相当复杂。

  • 其次，最终结果比中间过程要简洁得多，并且总是看起来与单变量情况相似。在这种情况下，请注意\(\frac{d}{dx}(bx) = b\)和\(\frac{d}{d\mathbf{x}} (\boldsymbol{\beta}^\top\mathbf{x}) = \boldsymbol{\beta}\)都是相似的。

  • 第三，转置常常会看似凭空出现。其核心原因是我们的约定是匹配分母的形状，因此当我们乘以矩阵时，我们需要进行转置以匹配回原始项的形状。

为了继续建立直觉，让我们尝试一个稍微难一点的计算。假设我们有一个列向量\(\mathbf{x}\)和一个方阵\(A\)，我们想要计算

(22.4.39)\[\frac{d}{d\mathbf{x}}(\mathbf{x}^\top A \mathbf{x}).\]

为了使用更易于操作的符号，让我们用爱因斯坦求和约定来考虑这个问题。在这种情况下，我们可以将函数写为

(22.4.40)\[\mathbf{x}^\top A \mathbf{x} = x_ia_{ij}x_j.\]

为了计算我们的导数，我们需要对每个\(k\)，理解

(22.4.41)\[\frac{d}{dx_k}(\mathbf{x}^\top A \mathbf{x}) = \frac{d}{dx_k}x_ia_{ij}x_j.\]

根据乘法法则，这是

(22.4.42)\[\frac{d}{dx_k}x_ia_{ij}x_j = \frac{dx_i}{dx_k}a_{ij}x_j + x_ia_{ij}\frac{dx_j}{dx_k}.\]

对于像\(\frac{dx_i}{dx_k}\)这样的项，不难看出当\(i=k\)时为1，否则为0。这意味着在这个和中，所有\(i\)和\(k\)不同的项都消失了，所以第一个和中剩下的项只有\(i=k\)的那些。同样的推理也适用于第二个项，其中我们需要\(j=k\)。这得到

(22.4.43)\[\frac{d}{dx_k}x_ia_{ij}x_j = a_{kj}x_j + x_ia_{ik}.\]

现在，爱因斯坦求和约定中的索引名称是任意的——\(i\)和\(j\)不同在这一点上对计算无关紧要，所以我们可以重新索引，使它们都使用\(i\)，看到

(22.4.44)\[\frac{d}{dx_k}x_ia_{ij}x_j = a_{ki}x_i + x_ia_{ik} = (a_{ki} + a_{ik})x_i.\]

现在，这里是我们开始需要一些练习才能更进一步的地方。让我们试着从矩阵运算的角度来识别这个结果。\(a_{ki} + a_{ik}\)是\(\mathbf{A} + \mathbf{A}^\top\)的\(k, i\)-th分量。这得到

(22.4.45)\[\frac{d}{dx_k}x_ia_{ij}x_j = [\mathbf{A} + \mathbf{A}^\top]_{ki}x_i.\]

类似地，这个项现在是矩阵\(\mathbf{A} + \mathbf{A}^\top\)与向量\(\mathbf{x}\)的乘积，所以我们看到

(22.4.46)\[\left[\frac{d}{d\mathbf{x}}(\mathbf{x}^\top A \mathbf{x})\right]_k = \frac{d}{dx_k}x_ia_{ij}x_j = [(\mathbf{A} + \mathbf{A}^\top)\mathbf{x}]_k.\]

因此，我们看到(22.4.39)中所需导数的第\(k\)项就是右边向量的第\(k\)项，因此两者是相同的。因此得到

(22.4.47)\[\frac{d}{d\mathbf{x}}(\mathbf{x}^\top A \mathbf{x}) = (\mathbf{A} + \mathbf{A}^\top)\mathbf{x}.\]

这比我们上一个计算需要更多的工作，但最终结果很小。更重要的是，考虑以下传统单变量导数的计算

(22.4.48)\[\frac{d}{dx}(xax) = \frac{dx}{dx}ax + xa\frac{dx}{dx} = (a+a)x.\]

等价地\(\frac{d}{dx}(ax^2) = 2ax = (a+a)x\)。再次，我们得到一个看起来很像单变量结果但加了一个转置的结果。

此时，这个模式应该看起来相当可疑，所以让我们试着找出原因。当我们像这样取矩阵导数时，让我们首先假设我们得到的表达式将是另一个矩阵表达式：一个我们可以用矩阵及其转置的乘积和和来写的表达式。如果存在这样的表达式，它将需要对所有矩阵都成立。特别是，它需要对\(1 \times 1\)矩阵成立，在这种情况下，矩阵乘积就是数的乘积，矩阵和就是和，转置什么也不做！换句话说，我们得到的任何表达式都必须与单变量表达式匹配。这意味着，经过一些练习，人们通常可以通过知道相关的单变量表达式必须是什么样子来猜测矩阵导数！

让我们试试这个。假设\(\mathbf{X}\)是一个\(n \times m\)矩阵，\(\mathbf{U}\)是一个\(n \times r\)矩阵，\(\mathbf{V}\)是一个\(r \times m\)矩阵。让我们试着计算

(22.4.49)\[\frac{d}{d\mathbf{V}} \|\mathbf{X} - \mathbf{U}\mathbf{V}\|_2^{2} = \;?\]

这个计算在一个称为矩阵分解的领域中很重要。然而，对我们来说，它只是一个要计算的导数。让我们试着想象一下对于\(1\times1\)矩阵会是什么样子。在这种情况下，我们得到表达式

(22.4.50)\[\frac{d}{dv} (x-uv)^{2}= -2(x-uv)u,\]

其中，导数是相当标准的。如果我们试着把它转换回一个矩阵表达式，我们得到

(22.4.51)\[\frac{d}{d\mathbf{V}} \|\mathbf{X} - \mathbf{U}\mathbf{V}\|_2^{2}= -2(\mathbf{X} - \mathbf{U}\mathbf{V})\mathbf{U}.\]

然而，如果我们看一下这个，它并不完全正确。回想一下，\(\mathbf{X}\)是\(n \times m\)，\(\mathbf{U}\mathbf{V}\)也是，所以矩阵\(2(\mathbf{X} - \mathbf{U}\mathbf{V})\)是\(n \times m\)。另一方面，\(\mathbf{U}\)是\(n \times r\)，我们不能将一个\(n \times m\)和一个\(n \times r\)矩阵相乘，因为维度不匹配！

我们想得到\(\frac{d}{d\mathbf{V}}\)，它与\(\mathbf{V}\)的形状相同，即\(r \times m\)。所以不知何故，我们需要取一个\(n \times m\)矩阵和一个\(n \times r\)矩阵，将它们相乘（可能带有一些转置）以得到一个\(r \times m\)矩阵。我们可以通过将\(U^\top\)乘以\((\mathbf{X} - \mathbf{U}\mathbf{V})\)来做到这一点。因此，我们可以猜测(22.4.49)的解是

(22.4.52)\[\frac{d}{d\mathbf{V}} \|\mathbf{X} - \mathbf{U}\mathbf{V}\|_2^{2}= -2\mathbf{U}^\top(\mathbf{X} - \mathbf{U}\mathbf{V}).\]

为了证明这是有效的，我们不提供详细的计算是失职的。如果您已经相信这个经验法则有效，可以随时跳过这个推导。要计算

(22.4.53)\[\frac{d}{d\mathbf{V}} \|\mathbf{X} - \mathbf{U}\mathbf{V}\|_2^2,\]

我们必须对每个\(a\)和\(b\)找到

(22.4.54)\[\frac{d}{dv_{ab}} \|\mathbf{X} - \mathbf{U}\mathbf{V}\|_2^{2}= \frac{d}{dv_{ab}} \sum_{i, j}\left(x_{ij} - \sum_k u_{ik}v_{kj}\right)^2.\]

回想一下，就\(\frac{d}{dv_{ab}}\)而言，\(\mathbf{X}\)和\(\mathbf{U}\)的所有条目都是常数，我们可以将导数推入和式中，并对平方应用链式法则得到

(22.4.55)\[\frac{d}{dv_{ab}} \|\mathbf{X} - \mathbf{U}\mathbf{V}\|_2^{2}= \sum_{i, j}2\left(x_{ij} - \sum_k u_{ik}v_{kj}\right)\left(-\sum_k u_{ik}\frac{dv_{kj}}{dv_{ab}} \right).\]

和之前的推导一样，我们可以注意到\(\frac{dv_{kj}}{dv_{ab}}\)仅当\(k=a\)且\(j=b\)时才非零。如果这些条件中的任何一个不成立，和中的项就为零，我们可以自由地丢弃它。我们看到

(22.4.56)\[\frac{d}{dv_{ab}} \|\mathbf{X} - \mathbf{U}\mathbf{V}\|_2^{2}= -2\sum_{i}\left(x_{ib} - \sum_k u_{ik}v_{kb}\right)u_{ia}.\]

这里有一个重要的微妙之处，即\(k=a\)的要求不会出现在内层和中，因为那个\(k\)是一个哑变量，我们在内层项中对其求和。举一个符号更清晰的例子，考虑为什么

(22.4.57)\[\frac{d}{dx_1} \left(\sum_i x_i \right)^{2}= 2\left(\sum_i x_i \right).\]

从这一点开始，我们可以开始识别和的组成部分。首先，

(22.4.58)\[\sum_k u_{ik}v_{kb} = [\mathbf{U}\mathbf{V}]_{ib}.\]

所以和内部的整个表达式是

(22.4.59)\[x_{ib} - \sum_k u_{ik}v_{kb} = [\mathbf{X}-\mathbf{U}\mathbf{V}]_{ib}.\]

这意味着我们现在可以将我们的导数写为

(22.4.60)\[\frac{d}{dv_{ab}} \|\mathbf{X} - \mathbf{U}\mathbf{V}\|_2^{2}= -2\sum_{i}[\mathbf{X}-\mathbf{U}\mathbf{V}]_{ib}u_{ia}.\]

我们希望它看起来像一个矩阵的\(a, b\)元素，这样我们就可以像上一个例子那样使用技术得到一个矩阵表达式，这意味着我们需要交换\(u_{ia}\)上索引的顺序。如果我们注意到\(u_{ia} = [\mathbf{U}^\top]_{ai}\)，我们就可以写

(22.4.61)\[\frac{d}{dv_{ab}} \|\mathbf{X} - \mathbf{U}\mathbf{V}\|_2^{2}= -2\sum_{i} [\mathbf{U}^\top]_{ai}[\mathbf{X}-\mathbf{U}\mathbf{V}]_{ib}.\]

这是一个矩阵乘积，因此我们可以得出结论

(22.4.62)\[\frac{d}{dv_{ab}} \|\mathbf{X} - \mathbf{U}\mathbf{V}\|_2^{2}= -2[\mathbf{U}^\top(\mathbf{X}-\mathbf{U}\mathbf{V})]_{ab}.\]

因此我们可以写出(22.4.49)的解

(22.4.63)\[\frac{d}{d\mathbf{V}} \|\mathbf{X} - \mathbf{U}\mathbf{V}\|_2^{2}= -2\mathbf{U}^\top(\mathbf{X} - \mathbf{U}\mathbf{V}).\]

这与我们上面猜测的解相符！

在这一点上，问“为什么我不能直接写下我学过的所有微积分规则的矩阵版本？很明显这仍然是机械的。为什么我们不一次性解决它！”是合理的。确实有这样的规则，(Petersen and Pedersen, 2008)提供了一个很好的总结。然而，由于矩阵运算可以组合的方式比单值多得多，矩阵导数规则比单变量规则多得多。通常情况下，最好是处理索引，或者在适当的时候留给自动微分。

22.4.8. 小结

  • 在高维空间中，我们可以定义梯度，其作用与一维中的导数相同。这使我们能够看到当对输入进行任意微小改变时，多变量函数如何变化。

  • 反向传播算法可以看作是一种组织多变量链式法则的方法，以允许高效地计算许多偏导数。

  • 矩阵微积分使我们能够以简洁的方式写出矩阵表达式的导数。

22.4.9. 练习

  1. 给定一个列向量\(\boldsymbol{\beta}\)，计算\(f(\mathbf{x}) = \boldsymbol{\beta}^\top\mathbf{x}\)和\(g(\mathbf{x}) = \mathbf{x}^\top\boldsymbol{\beta}\)的导数。为什么你得到相同的答案？

  2. 设\(\mathbf{v}\)是一个\(n\)维向量。那么\(\frac{\partial}{\partial\mathbf{v}}\|\mathbf{v}\|_2\)是什么？

  3. 设\(L(x, y) = \log(e^x + e^y)\)。计算梯度。梯度的分量之和是多少？

  4. 设\(f(x, y) = x^2y + xy^2\)。证明唯一的临界点是\((0,0)\)。通过考虑\(f(x, x)\)，判断\((0,0)\)是最大值、最小值还是都不是。

  5. 假设我们正在最小化一个函数\(f(\mathbf{x}) = g(\mathbf{x}) + h(\mathbf{x})\)。我们如何从几何上解释条件\(\nabla f = 0\)在\(g\)和\(h\)方面的意义？

pytorchmxnettensorflow

讨论

讨论

讨论