矩阵求导笔记

最近因为课程需要，学习了一下矩阵求导的方法，做个记录。

之前一直感觉矩阵求导很神奇又高大上，毕竟对一个矩阵求导很难像对一个数求导那样有一个直观的理解。后来发现对矩阵求导不过就是多元函数的求导，只是为了记法上的统一性 ，需要有一套对矩阵求导自己的方式，因为在对矩阵求导时有一个重要的前提，即独立性基本假设：假定实值函数的向量变元 $x=[x _ i]^{m} _ {i=1} \in \mathbb{R}^m$ 或者矩阵变元 $ X=[x _ {ij}]^{m,n} _ {i=1,j=1} \in \mathbb{R}^{m \times n} $ 本身无任何特殊结构，即向量或矩阵变元的元素之间是各自独立的。

上述独立性基本假设可以用数学公式表示成

$$ \begin{equation} \frac{\partial x _ i}{\partial x _ j} = \delta _ {ij} = \begin{cases} 1 , & \ i=j \\
0 , & \ \text{others} \end{cases} \end{equation} \tag{1} $$

以及

$$ \begin{equation} \frac{\partial x _ {kl}}{\partial x _ {ij} } = \delta _ {ki}\delta _ {lj} = \begin{cases} 1 , &\ k=i\ \text{and} \ l=j \\
0 , &\ \text{others} \end{cases} \end{equation} \tag{2} $$

根据以上的基本公式可以发现，矩阵求导和对多个多元函数求导并没有本质的区别。对于以下的多元函数组：

$$ \begin{equation} \begin{cases} y _ 1 &= f _ 1(x _ 1, x _ 2, \cdots,x _ n) \\
y _ 2 &= f _ 2(x _ 1, x _ 2, \cdots,x _ n) \\
&\vdots \\
y _ m &= f _ m(x _ 1, x _ 2, \cdots,x _ n) \end{cases} \end{equation} \tag{3} $$

我们求导需要让每一个 $y$ 对每个 $x$ 进行求偏导，得到 $$ \frac{\partial y_1}{\partial x_1}, \frac{\partial y_1}{\partial x_2}, \cdots \frac{\partial y_1}{\partial x_n}, \frac{\partial y_2}{\partial x_1}, \cdots \frac{\partial y_m}{\partial x_n} \tag{4} $$ 对矩阵求导也是同样的道理，只不过是因为我们上面的写法破坏了矩阵的统一性，所以需要使用矩阵的记法来表示相同的含义。

为了方便理解，对变元和函数作统一的符号规定：

$\boldsymbol{x}=[x_1,\cdots,x_m]^{\rm T}\in \mathbb{R}^m$ 为实向量变元；

$\boldsymbol{X}=[\boldsymbol{x_1},\cdots,\boldsymbol{x_n}]\in \mathbb{R}^{m \times n}$ 为实矩阵变元；

$f(\boldsymbol{x}) \in \mathbb{R}$ 为实值标量函数，其变元为 $m \times 1$ 实值向量 $\boldsymbol{x}$，记作 $f:\ \mathbb{R}^m \rightarrow \mathbb{R}$；

$f(\boldsymbol{X}) \in \mathbb{R}$ 为实值标量函数，其变元为 $m \times n$ 实值矩阵 $\boldsymbol{X}$，记作 $f:\ \mathbb{R}^{m \times n} \rightarrow \mathbb{R}$；

$\boldsymbol{f(\boldsymbol{x})} \in \mathbb{R}^p$ 为 $p$ 维实列向量函数，其变元为 $m \times 1$ 实值向量 $\boldsymbol{x}$，记作 $\boldsymbol{f}:\ \mathbb{R}^{m} \rightarrow \mathbb{R}^{p}$；

$\boldsymbol{f(\boldsymbol{X})} \in \mathbb{R}^p$ 为 $p$ 维实列向量函数，其变元为 $m \times n$ 实值矩阵 $\boldsymbol{X}$，记作 $\boldsymbol{f}:\ \mathbb{R}^{m} \rightarrow \mathbb{R}^{p}$；

$\boldsymbol{F(\boldsymbol{x})} \in \mathbb{R}^{p \times q}$ 为 $p \times q$ 实矩阵函数，其变元为 $m \times 1$ 实值向量 $\boldsymbol{x}$，记作 $\boldsymbol{F}:\ \mathbb{R}^{m} \rightarrow \mathbb{R}^{p \times q}$；

$\boldsymbol{F(\boldsymbol{X})} \in \mathbb{R}^{p \times q}$ 为 $p \times q$ 实矩阵函数，其变元为 $m \times n$ 实矩阵 $\boldsymbol{X}$，记作 $\boldsymbol{F}:\ \mathbb{R}^{m \times n} \rightarrow \mathbb{R}^{p \times q}$。

向量的求导

偏导向量

$1 \times m$ 行向量偏导算子记为 $$ {\rm D} _ \boldsymbol{x} \overset{{\rm def}}{=} \frac{\partial}{\partial \boldsymbol{x}^{\rm T}} = [\frac{\partial}{\partial x _ 1}, \cdots ,\frac{\partial}{\partial x _ m}] \tag{5} $$ 于是，实值标量函数 $f(\boldsymbol{x})$ 在 $x$ 的偏导向量由 $1 \times m$ 行向量 $$ {\rm D} _ \boldsymbol{x}f(\boldsymbol{x}) = \frac{\partial f(\boldsymbol{x})}{\partial \boldsymbol{x}^{\rm T}} = [\frac{\partial f(\boldsymbol{x})}{\partial x _ 1}, \cdots ,\frac{\partial f(\boldsymbol{x})}{\partial x _ m}] \tag{6} $$ 给出。行偏导向量即向量的 $Jacobian$ 矩阵。

梯度向量

采用列向量形式定义的偏导算子称为列向量偏导算子，习惯称为梯度算子。

$m \times 1$ 列向量偏导算子即梯度算子记作 $\nabla _\boldsymbol{x}$，定义为 $$ \nabla _\boldsymbol{x} \overset{{\rm def}}{=} \frac{\partial}{\partial \boldsymbol{x}} = [\frac{\partial}{\partial x_1}, \cdots ,\frac{\partial}{\partial x_m}]^{\rm T} \tag{7} $$ 因此，实值标量函数 $f(\boldsymbol{x})$ 的梯度向量 $\nabla _\boldsymbol{x} f(\boldsymbol{x})$ 为 $m \times 1$ 列向量，定义为 $$ \nabla _\boldsymbol{x}f(\boldsymbol{x}) \overset{{\rm def}}{=} [\frac{\partial f(\boldsymbol{x})}{\partial x_1}, \cdots ,\frac{\partial f(\boldsymbol{x})}{\partial x_m}]^{\rm T} = \frac{\partial f(\boldsymbol{x})}{\partial \boldsymbol{x}} \tag{8} $$ 容易得到， $$ \nabla _ \boldsymbol{x}f(\boldsymbol{x}) = {\rm D} _ \boldsymbol{x}^{\rm T}f(\boldsymbol{x}) \tag{9} $$ 即，实值标量函数 $f(\boldsymbol{x})$ 的梯度矩阵等于 $Jacobian$ 矩阵的转置。

实值标量函数求导

标量函数 $f(\boldsymbol{x})$ 的 $Jacobian$ 矩阵与微分矩阵之间存在等价关系 $$ {\rm d}f(\boldsymbol{x}) = {\rm tr}(\boldsymbol{A}{\rm d}\boldsymbol{x}) \Leftrightarrow {\rm D}_\boldsymbol{x}f(\boldsymbol{x}) = \frac{\partial f(\boldsymbol{x})}{\partial \boldsymbol{x}^{\rm T}} = \boldsymbol{A} \tag{10} $$ 因此，只要将函数 $f(\boldsymbol{x})$ 的微分写作 ${\rm d}f(\boldsymbol{x}) = {\rm tr}(\boldsymbol{A}{\rm d}\boldsymbol{x}) $ 的形式，则矩阵 $\boldsymbol{A}$ 就是函数 $f(\boldsymbol{x})$ 关于其变元向量 $\boldsymbol{x}$ 的 $Jacobian$ 矩阵，$\boldsymbol{A}^{\rm T}$ 即微分矩阵，也即梯度矩阵。

在对 $f(\boldsymbol{x})$ 套上迹后，可以使用一些迹技巧进行运算：

• ${\rm tr}(c) = c$，即标量的迹等于标量本身
• ${\rm tr}(\boldsymbol{A} + \boldsymbol{B}) = {\rm tr}(\boldsymbol{A}) + {\rm tr}(\boldsymbol{B})$
• ${\rm tr}(c\boldsymbol{A}) = c {\rm tr}(\boldsymbol{A})$
• ${\rm tr}(\boldsymbol{A}\boldsymbol{B}) = {\rm tr}(\boldsymbol{B}\boldsymbol{A})$
• ${\rm tr}(\boldsymbol{A}_1\boldsymbol{A}_2\cdots\boldsymbol{A}_n) = {\rm tr}(\boldsymbol{A}_n \boldsymbol{A} _ 1 \cdots \boldsymbol{A} _ {n-1})$
• ${\rm tr}(\boldsymbol{A}) = {\rm tr}(\boldsymbol{A}^{\rm T})$
• ${\rm tr}(\boldsymbol{A}^T(\boldsymbol{B}\odot\boldsymbol{C})) = {\rm tr}((\boldsymbol{A}\odot\boldsymbol{B})^T\boldsymbol{C})$，$\odot$ 表示矩阵的 Hadamard 积，即尺寸相同的矩阵的逐元素乘法，有的教材也用 $\circ$ 表示

此外，我们还需要一些求导运算法则：

• ${\rm d}(\boldsymbol{X} \pm \boldsymbol{Y}) = {\rm d}\boldsymbol{X} \pm {\rm d}\boldsymbol{Y}$
• ${\rm d}(\boldsymbol{X}\boldsymbol{Y}) = ({\rm d}\boldsymbol{X})\boldsymbol{Y} + \boldsymbol{X}{\rm d}\boldsymbol{Y}$
• ${\rm d}(\boldsymbol{X}^{\rm T})=({\rm d}\boldsymbol{X})^{\rm T}$
• ${\rm dtr}(\boldsymbol{X}) = {\rm tr}({\rm d}\boldsymbol{X})$
• ${\rm d}(\boldsymbol{X}\odot\boldsymbol{Y}) = {\rm d}\boldsymbol{X}\odot\boldsymbol{Y} + \boldsymbol{X}\odot{\rm d}\boldsymbol{Y}$
• ${\rm d}\sigma(\boldsymbol{X})=\sigma'(\boldsymbol{X})\odot{\rm d}\boldsymbol{X}$，$\sigma(\boldsymbol{X})=[\sigma(X _ {ij})]$ 是逐元素标量函数运算，$\sigma'(\boldsymbol{X})=[\sigma'(X _ {ij})]$ 是逐元素求导运算

有了以上的法则，我们可以通过迹技巧和求导运算，将函数 $f(\boldsymbol{x})$ 的微分写作 ${\rm d}f(\boldsymbol{x}) = {\rm tr}(\boldsymbol{A}{\rm d}\boldsymbol{x}) $ 的形式，$\boldsymbol{A}$ 即为 $Jacobian$ 矩阵。

例题

例1 求实值函数 $f(\boldsymbol{x})=\boldsymbol{x}^{\rm T}\boldsymbol{A}\boldsymbol{x}$ 的 $Jacobian$ 矩阵。 $$ \begin{aligned} f(\boldsymbol{x})&=\boldsymbol{x}^{\rm T}\boldsymbol{A}\boldsymbol{x} \\
{\rm d}f(\boldsymbol{x})&={\rm dtr}(\boldsymbol{x}^{\rm T}\boldsymbol{A}\boldsymbol{x}) \\
&={\rm tr}(({\rm d}\boldsymbol{x}^{\rm T})\boldsymbol{A}\boldsymbol{x} + \boldsymbol{x}^{\rm T}\boldsymbol{A}{\rm d}\boldsymbol{x}) \\
&={\rm tr}(\boldsymbol{x}^{\rm T}\boldsymbol{A}^{\rm T}{\rm d}\boldsymbol{x} + \boldsymbol{x}^{\rm T}\boldsymbol{A}{\rm d}\boldsymbol{x}) \\
&= {\rm tr}(\boldsymbol{x}^{\rm T}(\boldsymbol{A}^{\rm T} + \boldsymbol{A}){\rm d}\boldsymbol{x}) \end{aligned} \tag{11} $$ 因此，$f(\boldsymbol{x})$ 的 $Jacobian$ 矩阵为 $\boldsymbol{x}^{\rm T}(\boldsymbol{A}^{\rm T} + \boldsymbol{A})$，梯度矩阵为 $(\boldsymbol{A} + \boldsymbol{A}^{\rm T})\boldsymbol{x}$。

例2 推导 MLP 中反向传播的随机梯度下降公式。

在如图所示的 MLP 中，展示了包括输出层、最后一个隐藏层以及倒数第二个隐藏层的三层神经网络，其余部分未画出。其中，$\boldsymbol{W}^{(l)}$ 和 $\boldsymbol{B}^{(l)}$ 表示第 $l$ 层神经元的权重和偏置项参数，$\boldsymbol{x} _ {out} ^ {(l)}$ 表示第 $l$ 层神经元的输出向量，$\boldsymbol{v} ^ {(l)}$ 表示第 $l$ 层神经元在经过激活函数之前的输出向量，$\phi(\cdot)$ 表示激活函数，则有 $$ \boldsymbol{x} _ {out} ^ {(l)} = \phi (\boldsymbol{v} ^ {(l)}) \tag{12} $$ 输出层神经元的输出为 $\boldsymbol{y}$，即 $\boldsymbol{x} _ {out} ^ {(L)}$，预期输出即真值为 $\boldsymbol{d}$，设损失函数为 $$ E=\frac{1}{2}(\boldsymbol{y} - \boldsymbol{d})^{\rm T}(\boldsymbol{y} - \boldsymbol{d}) \tag{13} $$ 在这里，我们引入一个新的中间变量 $$ \boldsymbol{\delta}^{(l)} = \frac{\partial E}{\partial \boldsymbol{v} ^ {(l)}} \tag{14} $$ 它表示损失函数对第 $l$ 层神经元经过激活函数前的输出值的偏导，之所以引入这个值是为了更加方便地计算 $\boldsymbol{W}^{(l)}$ 和 $\boldsymbol{B}^{(l)}$ 的梯度，我们也将很快看到，$\boldsymbol{\delta}^{(l)}$ 可以很方便地通过神经网络进行反向传播，使我们得到非常简洁的表示形式。

1. 首先我们从输出层开始，求 $\boldsymbol{\delta}^{(L)}$。

$$ \begin{aligned} E&=\frac{1}{2}(\boldsymbol{x} _ {out} ^ {(L)} - \boldsymbol{d})^{\rm T}(\boldsymbol{x} _ {out} ^ {(L)} - \boldsymbol{d}) \\
{\rm d}E &= {\rm dtr}(\frac{1}{2}(\boldsymbol{x} _ {out} ^ {(L)} - \boldsymbol{d})^{\rm T}(\boldsymbol{x} _ {out} ^ {(L)} - \boldsymbol{d})) \\
&= {\rm tr} (\frac{1}{2}{\rm d}\boldsymbol{x} _ {out} ^ {(L)\rm T}(\boldsymbol{x} _ {out} ^ {(L)} - \boldsymbol{d})+\frac{1}{2}(\boldsymbol{x} _ {out} ^ {(L)} - \boldsymbol{d})^{\rm T}{\rm d}\boldsymbol{x} _ {out} ^ {(L)}) \\
&= {\rm tr}((\boldsymbol{x} _ {out} ^ {(L)} - \boldsymbol{d})^{\rm T}{\rm d}\boldsymbol{x} _ {out} ^ {(L)}) \end{aligned} \tag{15} $$

又根据 $$ \boldsymbol{x} _ {out} ^ {(l)} = \phi (\boldsymbol{v} ^ {(l)}) \tag{16} $$ 可以得到 $$ {\rm d}\boldsymbol{x} _ {out} ^ {(L)} = \phi^{'} (\boldsymbol{v} ^ {(L)})\odot{\rm d}\boldsymbol{v} ^ {(L)} \tag{17} $$ 将 $(17)$ 式中的 ${\rm d}\boldsymbol{x} _ {out} ^ {(L)}$ 代入 $(15)$ 式，可以得到 $$ \begin{aligned} {\rm d}E &= {\rm tr}((\boldsymbol{x} _ {out} ^ {(L)} - \boldsymbol{d})^{\rm T}{\rm d}\boldsymbol{x} _ {out} ^ {(L)}) \\
&= {\rm tr}((\boldsymbol{x} _ {out} ^ {(L)} - \boldsymbol{d})^{\rm T} [\phi^{'} (\boldsymbol{v} ^ {(L)})\odot{\rm d}\boldsymbol{v} ^ {(L)}]) \\
&= {\rm tr}([(\boldsymbol{x} _ {out} ^ {(L)} - \boldsymbol{d})\odot\phi^{'} (\boldsymbol{v} ^ {(L)})]^{\rm T}{\rm d}\boldsymbol{v} ^ {(L)}) \end{aligned} \tag{18} $$ 则有 $$ \begin{aligned} \boldsymbol{\delta}^{(L)} &= \frac{\partial E}{\partial \boldsymbol{v} ^ {(L)}} \\
&= (\boldsymbol{x} _ {out} ^ {(L)} - \boldsymbol{d})\odot\phi^{'} (\boldsymbol{v} ^ {(L)}) \\
&= (\boldsymbol{y} - \boldsymbol{d})\odot\phi^{'} (\boldsymbol{v} ^ {(L)}) \end{aligned} \tag{19} $$ 至此，我们计算出了最后一层的 $\boldsymbol{\delta}^{(L)}$，对于中间层 $l$，我们使用相同的方式求出 $\boldsymbol{\delta}^{(l+1)}$ 和 $\boldsymbol{\delta}^{(l)}$ 的递推关系。有了最后一层 $\boldsymbol{\delta}^{(L)}$ 和层之间的递推关系，我们就能从后向前计算出每一层的 $\boldsymbol{\delta}^{(l)}$，这也正是反向传播的灵魂所在。在得到了 $\boldsymbol{\delta}^{(l)}$ 之后，我们很快将看到如何用 $\boldsymbol{\delta}^{(l)}$ 方便快速地计算 $\boldsymbol{W}^{(l)}$ 和 $\boldsymbol{B}^{(l)}$ 的梯度，这也是我们的最终目标。

2. 求 $\boldsymbol{\delta}^{(l+1)}$ 和 $\boldsymbol{\delta}^{(l)}$ 的递推关系。

$$ \begin{aligned} {\rm d}E &= {\rm tr}( (\frac{\partial E}{\partial \boldsymbol{v} ^ {(l+1)}})^{\rm T} {\rm d}\boldsymbol{v} ^ {(l+1)}) \\
&= {\rm tr}( \boldsymbol{\delta}^{(l+1)\rm T}{\rm d}\boldsymbol{v} ^ {(l+1)}) \end{aligned} \tag{20} $$

对于第 $l+1$ 层，有 $$ \begin{aligned} \boldsymbol{v} ^ {(l+1)} &= \boldsymbol{W}^{(l+1)}\phi(\boldsymbol{v} ^ {(l)}) + \boldsymbol{B}^{(l+1)} \\
{\rm d}\boldsymbol{v} ^ {(l+1)} &= \boldsymbol{W}^{(l+1)}{\rm d}\phi(\boldsymbol{v} ^ {(l)}) \\
&= \boldsymbol{W}^{(l+1)} (\phi^{'} (\boldsymbol{v} ^ {(l)})\odot{\rm d}\boldsymbol{v} ^ {(l)}) \end{aligned} \tag{21} $$ 将 $(21)$ 代入 $(20)$，有 $$ \begin{aligned} {\rm d}E &= {\rm tr}( \boldsymbol{\delta}^{(l+1){\rm T}}{\rm d}\boldsymbol{v} ^ {(l+1)}) \\
&= {\rm tr}(\boldsymbol{\delta}^{(l+1){\rm T}}\boldsymbol{W}^{(l+1)} [\phi^{'} (\boldsymbol{v} ^ {(l)})\odot{\rm d}\boldsymbol{v} ^ {(l)}]) \\
&= {\rm tr}([\boldsymbol{W}^{(l+1){\rm T}}\boldsymbol{\delta}^{(l+1)} \odot \phi^{'} (\boldsymbol{v} ^ {(l)})]^{\rm T}{\rm d}\boldsymbol{v} ^ {(l)}) \\
\boldsymbol{\delta}^{(l)} &= \boldsymbol{W}^{(l+1){\rm T}}\boldsymbol{\delta}^{(l+1)} \odot \phi^{'} (\boldsymbol{v} ^ {(l)}) \end{aligned} \tag{22} $$ 至此，我们得到了从 $\boldsymbol{\delta}^{(l+1)}$ 递推到 $\boldsymbol{\delta}^{(l)}$ 的表达式，这正是反向传播的最核心的过程。

3. 计算 $\boldsymbol{W}^{(l)}$ 和 $\boldsymbol{B}^{(l)}$ 的梯度。

对于第 $l$ 层，计算 $E$ 对 $\boldsymbol{W}^{(l)}$ 的梯度

$$ \begin{aligned} \boldsymbol{v} ^ {(l)} &= \boldsymbol{W}^{(l)} \boldsymbol{x} _ {out} ^ {(l-1)} + \boldsymbol{B}^{(l)} \\
{\rm d}\boldsymbol{v} ^ {(l)} &= {\rm d}\boldsymbol{W}^{(l)} \boldsymbol{x} _ {out} ^ {(l-1)} \\
{\rm d}E &= {\rm tr}( \boldsymbol{\delta}^{(l){\rm T}}{\rm d}\boldsymbol{v} ^ {(l)}) \\
&= {\rm tr}( \boldsymbol{\delta}^{(l){\rm T}}{\rm d}\boldsymbol{W}^{(l)} \boldsymbol{x} _ {out} ^ {(l-1)}) \\
&= {\rm tr}( \boldsymbol{x} _ {out} ^ {(l-1)} \boldsymbol{\delta}^{(l){\rm T}}{\rm d}\boldsymbol{W}^{(l)}) \\
&= {\rm tr}( [\boldsymbol{\delta}^{(l)}\boldsymbol{x} _ {out} ^ {(l-1){\rm T}}]^{\rm T} {\rm d}\boldsymbol{W}^{(l)}) \\
\frac{\partial E}{\partial \boldsymbol{W}^{(l)}} &= \boldsymbol{\delta}^{(l)}\boldsymbol{x} _ {out} ^ {(l-1){\rm T}} \end{aligned} \tag{23} $$

计算 $E$ 对 $\boldsymbol{B}^{(l)}$ 的梯度

$$ \begin{aligned} \boldsymbol{v} ^ {(l)} &= \boldsymbol{W}^{(l)} \boldsymbol{x} _ {out} ^ {(l-1)} + \boldsymbol{B}^{(l)} \\
{\rm d}\boldsymbol{v} ^ {(l)} &= {\rm d}\boldsymbol{B}^{(l)} \\
{\rm d}E &= {\rm tr}( \boldsymbol{\delta}^{(l){\rm T}}{\rm d}\boldsymbol{v} ^ {(l)}) \\
&= {\rm tr}( \boldsymbol{\delta}^{(l){\rm T}}{\rm d}\boldsymbol{B}^{(l)}) \\
\frac{\partial E}{\partial \boldsymbol{B}^{(l)}} &= \boldsymbol{\delta}^{(l)} \end{aligned} \tag{24} $$

4. 梯度更新。

$$ \begin{aligned} \boldsymbol{W}^{(l)} _ {t+1} &= \boldsymbol{W}^{(l)} _ t- \eta \frac{\partial E}{\partial \boldsymbol{W}^{(l)} _ t} = \boldsymbol{W}^{(l)} _ t - \eta \boldsymbol{\delta}^{(l)}\boldsymbol{x} _ {out} ^ {(l-1){\rm T}} \\
\boldsymbol{B}^{(l)} _ {t+1} &= \boldsymbol{B}^{(l)} _ t - \eta \frac{\partial E}{\partial \boldsymbol{B}^{(l)} _ t} = \boldsymbol{B}^{(l)} _ t - \eta \boldsymbol{\delta}^{(l)} \end{aligned} \tag{25} $$

其中 $\eta$ 为学习率。

基于以上推导的公式，我用 numpy 实现了一个简单的只含有一层隐藏层的 MLP，输入为一个标量，输出也为一个标量，用它拟合了函数 $f(x) = 1.2sin(\pi x)-cos(2.4\pi x),\ \text{for} \ x\in[-1.6,1.6]$ 的图像。训练前和训练后的图像如下如示，可以看出来 MLP 通过反向传播与梯度下降成功学习到了合适的参数，python 代码也附在了后面。在实现时，有以下几点需要注意：

• 在实际中，一般会用 mini-batch 而不是一条数据作为一个输入，最终的 loss 为这个 mini-batch 所有 loss 的平均值。在使用 $(19)$ 计算最后一层的误差向量时，$(\boldsymbol{y} - \boldsymbol{d})$ 也应该再对整个 mini-batch 计算一个平均值，再反向传播给 mini-batch 中每一条数据，计算得到各自的梯度，将梯度相加得到最终的梯度。也就是说，在一个 mini-batch 中，只有一次梯度的更新，mini-batch 内部中数据的顺序对这次梯度的更新没有任何影响。
• 在实际中，输入 $\boldsymbol{X}$ 一般会写成 [batch_size * input_vectors] 的形式，因此在进行前向传播时需要右乘参数矩阵 $\boldsymbol{W}$，而不是上面推导过程中的左乘。因此在反向传播时，矩阵的乘法也都调了个位置，核心过程是一样的，但写的时候要特别注意矩阵的顺序。
• 在下面的实现中，我在 loss 函数中添加了正则项 $\frac{1}{2}\boldsymbol{w}_1\boldsymbol{w}_2$，因此在计算梯度时还有一部分是正则项对权重参数的梯度，这部分的值通过参数 lam 进行调节，当 lam 为 0 时，等于 loss 函数中没有正则项。
• 拟合函数属于回归问题，回归问题应当直接在输出层进行输出而不能再经过一次激活函数，因此输出层的 $\boldsymbol{v} ^ {(L)}$ 就是 $\boldsymbol{y}$，在用 $(19)$ 计算最后一层的误差向量时，注意到 $\phi^{'} (\boldsymbol{v} ^ {(L)})$ 分量全是 1，因此可以省去这一 Hadamard 乘积的步骤。

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103

import numpy as np
from math import sin, cos, pi
import matplotlib.pyplot as plt

def f(x):
    return 1.2 * sin(pi * x) - cos(2.4 * pi * x)

def generate_data(shuffle=True):
    data = []
    for x in np.arange(-1.6, 1.6, 0.01):
        data.append([x, f(x)])
    data = np.array(data)
    if shuffle:
        np.random.shuffle(data)
    return data[:, 0], data[:, 1]

class MLP(object):
    def __init__(self, n_hidden, batch_size, lr=1e-3, lam=0):
        self.n_hidden   = n_hidden
        self.batch_size = batch_size
        self.lr         = lr
        self.lam        = lam
        self.loss       = 0
        self.x      = np.zeros([batch_size, 1])                                 # [batch_size * 1]
        self.w1     = 2 * np.random.random([1, n_hidden]) - 1                   # [1 * n]
        self.b1     = 2 * np.random.random([1, n_hidden]) - 1                   # [1 * n]
        self.dif_v1 = np.zeros([batch_size, n_hidden])                          # [batch_size * n]
        self.o1     = np.zeros([batch_size, n_hidden])                          # [batch_size * n]
        self.delta1 = np.zeros([batch_size, n_hidden])                          # [batch_size * n]
        self.w2     = 2 * np.random.random([n_hidden, 1])                       # [n * 1]
        self.b2     = 2 * np.random.random([1, 1]) - 1                          # [1 * 1]
        self.delta2 = np.zeros([batch_size, 1])                                 # [batch_size * 1]
        self.y      = np.zeros([batch_size, 1])                                 # [batch_size * 1]
        self.d      = np.zeros([batch_size, 1])                                 # [batch_size * 1]

    def _tanh(self, x):
        return (np.exp(x) - np.exp(-x)) / (np.exp(x) + np.exp(-x))
    
    def _tanh_dif(self, x):
        return 1 - self._tanh(x) ** 2
    
    def _forward(self):
        # self.x: [batch_size * 1]
        v1          = np.dot(self.x, self.w1) + np.tile(self.b1, (self.batch_size, 1))     # [batch_size * n]
        self.o1     = self._tanh(v1)                                                       # [batch_size * n]
        self.dif_v1 = self._tanh_dif(v1)                                                   # [batch_size * n]
        v2          = np.dot(self.o1, self.w2) + np.tile(self.b2, (self.batch_size, 1))    # [batch_size * 1]
        self.y      = v2                                                                   # [batch_size * 1]

        return np.mean(1 / 2 * (self.y - self.d) ** 2) # MSE loss

    def _bp(self):
        self.delta2 = np.zeros([self.batch_size, 1]) + np.mean(self.y - self.d)                    # [batch_size * 1]
        self.delta1 = np.dot(self.delta2, self.w2.T) * self.dif_v1                                 # [batch_size * n]
        w1          = self.w1
        w2          = self.w2
        # self.lam / 2 * w2.T and self.lam / 2 * w1.T are the regularization terms
        self.w1    -= self.lr * (np.dot(self.x.T, self.delta1) + self.lam / 2 * w2.T)
        self.b1    -= self.lr * self.delta1.sum(0)
        self.w2    -= self.lr * (np.dot(self.o1.T, self.delta2) + self.lam / 2 * w1.T)
        self.b2    -= self.lr * self.delta2.sum(0)
    
    def train(self, epochs=200):
        steps = 0
        for epoch in range(epochs):
            x, d = generate_data()
            num = x.shape[0]
            for i in range(0, num, self.batch_size):
                self.x = x[i:(i+self.batch_size)].reshape(-1, 1)
                self.d = d[i:(i+self.batch_size)].reshape(-1, 1)
                loss = self._forward()
                self.loss = self.loss + 0.2 * (loss - self.loss) if self.loss else loss
                if steps % 10000 == 0:
                    print("Steps: {}\tLoss:{:.6f}".format(steps, self.loss))
                self._bp()
                steps += 1
    
    def predict(self, x):
        x = x.reshape(-1, 1)
        batch_size = x.shape[0]
        v1 = np.dot(x, self.w1) + np.tile(self.b1, (batch_size, 1))
        v2 = np.dot(self._tanh(v1), self.w2) + np.tile(self.b2, (batch_size, 1))
        
        return v2.reshape(batch_size)

    def plot(self):
        labels = ['x', 'y']
        x, d = generate_data(False)
        y = self.predict(x)
        fig, ax = plt.subplots()
        ax.plot(x, y, label='MLP output')
        ax.plot(x, d, label='Desired line')
        ax.set_xlabel('x')
        ax.set_ylabel('y')
        ax.set_title("Fitting result of MLP")
        ax.legend()
        plt.show()

if __name__ == '__main__':
    mlp = MLP(10, 5, lr=0.005, lam=0)
    mlp.plot()
    mlp.train(10000)
    mlp.plot()

向量对向量求导

向量函数 $\boldsymbol{f}(\boldsymbol{x}):\ \mathbb{R}^{m} \rightarrow \mathbb{R}^{p}$ 的 $Jacobian$ 矩阵可以通过下式辨识 $$ \begin{aligned} {\rm d}\boldsymbol{f}(\boldsymbol{x}) &= \boldsymbol{A}{\rm d}\boldsymbol{x} \\
\Leftrightarrow\ {\rm D}_\boldsymbol{x}\boldsymbol{f}(\boldsymbol{x}) &= \frac{\partial \boldsymbol{f}(\boldsymbol{x})}{\partial \boldsymbol{x}^{\rm T}} = \boldsymbol{A} \\
\Leftrightarrow\ \nabla _ \boldsymbol{x}\boldsymbol{f}(\boldsymbol{x}) &= \frac{\partial \boldsymbol{f}(\boldsymbol{x})^{\rm T}}{\partial \boldsymbol{x}} = \boldsymbol{A}^{\rm T} \end{aligned} \tag{26} $$ 其中，$\boldsymbol{A}$ 是一个 $p \times m$ 的矩阵。 $$ \boldsymbol{A} = \left[ \begin{matrix} \frac{\partial f_1(\boldsymbol{x})}{\partial x_1} & \frac{\partial f_1(\boldsymbol{x})}{\partial x_2} & \cdots & \frac{\partial f_1(\boldsymbol{x})}{\partial x_m} \\
\frac{\partial f_2(\boldsymbol{x})}{\partial x_1} & \frac{\partial f_2(\boldsymbol{x})}{\partial x_2} & \cdots & \frac{\partial f_2(\boldsymbol{x})}{\partial x_m} \\
\vdots & \vdots & & \vdots \\
\frac{\partial f_p(\boldsymbol{x})}{\partial x_1} & \frac{\partial f_p(\boldsymbol{x})}{\partial x_2} & \cdots & \frac{\partial f_p(\boldsymbol{x})}{\partial x_m} \end{matrix} \right] \tag{27} $$ 因此，只要将函数 $\boldsymbol{f}(\boldsymbol{x})$ 的微分写作 ${\rm d}\boldsymbol{f}(\boldsymbol{x}) = \boldsymbol{A}{\rm d}\boldsymbol{x} $ 的形式，则矩阵 $\boldsymbol{A}$ 就是函数 $\boldsymbol{f}(\boldsymbol{x})$ 关于其变元向量 $\boldsymbol{x}$ 的 $Jacobian$ 矩阵，$\boldsymbol{A}^{\rm T}$ 即微分矩阵，也即梯度矩阵。

例题

例3 设 MLP 中损失函数为 $softmax$ 后的交叉熵，推导最后一层的误差向量 $\boldsymbol{\delta}^{(L)}$。

从 例2 中我们已经看到，对于最小二乘法形式的损失函数，最后一层的误差向量为 $$ \boldsymbol{\delta}^{(L)} = (\boldsymbol{y} - \boldsymbol{d})\odot\phi^{'} (\boldsymbol{v} ^ {(L)}) \tag{28} $$ 其中 $\boldsymbol{y}$ 为网络的输出向量。实际上，通过之前的推导，我们可以很容易地得到以下结论：设损失函数 $E=f(\boldsymbol{y})$，最后一层的误差向量即为 $$ \boldsymbol{\delta}^{(L)} = \frac{\partial E}{\partial \boldsymbol{y}} \odot\phi^{'} (\boldsymbol{v} ^ {(L)}) \tag{29} $$ 因此，对于不同的损失函数，只需要计算 $\frac{\partial E}{\partial \boldsymbol{y}}$，即可快速得到最后一层的误差向量，并通过式 $(22)$ 向前进行反向传播。

我们使用 $s(\cdot)$ 来表示 $softmax$ 函数，用 $\boldsymbol{a}$ 表示网络原始输出 $\boldsymbol{y}$ 经过 $softmax$ 后的概率输出，$\boldsymbol{d}$ 是预期输出，为 one-hot 表示的向量，那么我们可以得到交叉熵的损失函数为 $$ \begin{aligned} \boldsymbol{a} &= s(\boldsymbol{y}) \\
E &= - \boldsymbol{d}^{\rm T} log(\boldsymbol{a}) \end{aligned} \tag{30} $$ 其中，$log(\cdot)$ 为逐元素的 $log$ 函数。我们的目标是求 $\frac{\partial E}{\partial \boldsymbol{y}}$，可以先将 $\boldsymbol{a}$ 作为中间变量，求 $\frac{\partial E}{\partial \boldsymbol{a}}$，这属于标量对向量的求导，参考上部分实值标量函数求导的过程，容易得到 $$ \begin{aligned} E &= - \boldsymbol{d}^{\rm T} log(\boldsymbol{a}) \\
{\rm d}E &= {\rm dtr}(- \boldsymbol{d}^{\rm T} {\rm d}log(\boldsymbol{a})) \\
&= {\rm dtr}(- \boldsymbol{d}^{\rm T} [log'(\boldsymbol{a}) \odot {\rm d}\boldsymbol{a}]) \\
&= {\rm dtr}([-\boldsymbol{d} \odot log'(\boldsymbol{a})]^{\rm T} {\rm d}\boldsymbol{a}) \end{aligned} \tag{31} $$ 其中，$log'(\cdot)$ 为 $log(\cdot)$ 的导函数，即逐元素求倒数。

现在我们来求 ${\rm d}\boldsymbol{a}$。由 $(26)$ 和 $(30)$ 可知，$\boldsymbol{a}$ 对 $\boldsymbol{y}$ 属于向量对向量的求导，得到的应该是一个 $n \times n$ 的矩阵，$n$ 为网络输出向量的维度。我们设 $\boldsymbol{a}$ 的 $Jacobian$ 矩阵为 $\boldsymbol{A}$，那么根据 $(26)$，有以下关系 $$ {\rm d}\boldsymbol{a} = \boldsymbol{A} {\rm d}\boldsymbol{y} \tag{32} $$ 我们将 $\boldsymbol{A}$ 的计算放一放，先将 $(32)$ 代入 $(31)$ 可以得到 $$ \begin{aligned} {\rm d}E &= {\rm dtr}([-\boldsymbol{d} \odot log'(\boldsymbol{a})]^{\rm T} {\rm d}\boldsymbol{a}) \\
&= {\rm dtr}([-\boldsymbol{d} \odot log'(\boldsymbol{a})]^{\rm T} \boldsymbol{A} {\rm d}\boldsymbol{y}) \\
&= {\rm dtr}([-\boldsymbol{A}^{\rm T}(\boldsymbol{d} \odot log'(\boldsymbol{a}))]^{\rm T} {\rm d}\boldsymbol{y}) \\
\frac{\partial E}{\partial \boldsymbol{y}} &= -\boldsymbol{A}^{\rm T}(\boldsymbol{d} \odot log'(\boldsymbol{a})) \end{aligned} \tag{33} $$ 此时，我们已经得到了 $\frac{\partial E}{\partial \boldsymbol{y}}$ 的表示形式。但显然，我们对这个表示是不满意的，因为我们还没有研究 $\boldsymbol{A}$ 的内部结构，$(33)$ 现在还无法用做实际的计算。很快我们将看到，利用 $\boldsymbol{d}$ 的 one-hot 特点，我们会得到非常简洁的最终结果。

$\boldsymbol{A}^{\rm T}$ 是一个 $n \times n$ 的矩阵，$\boldsymbol{d} \odot log'(\boldsymbol{a})$ 是一个 $n \times 1$ 的向量，因此最终的结果是一个 $n \times 1$ 的向量，这也与 $\boldsymbol{y}$ 的维度相同。因为 $\boldsymbol{d}$ 是 one-hot 编码的向量，这就意味着 $\boldsymbol{d} \odot log'(\boldsymbol{a})$ 只有一个非零值。我们设当前的正确分类为 $i$，那么 $\boldsymbol{d} \odot log'(\boldsymbol{a})$ 除了第 $i$ 维的分量为 $\frac{1}{a_i}$，其余全是 0。$-\boldsymbol{A}^{\rm T}(\boldsymbol{d} \odot log'(\boldsymbol{a}))$ 实际上就等于 $-\frac{1}{a_i}\boldsymbol{v}_i$，$\boldsymbol{v}_i$ 为 $\boldsymbol{A}^{\rm T}$ 第 $i$ 个列向量。

由 $(27)$ 和 $(32)$ 可知，$\boldsymbol{A}^{\rm T}$ 第 $i$ 个列向量即 $\boldsymbol{A}$ 第 $i$ 个行向量的转置，为

$$ \boldsymbol{v}_i = [\frac{\partial a_i}{\partial y_1}, \frac{\partial a_i}{\partial y_2}, \cdots, \frac{\partial a_i}{\partial y_n}]^{\rm T} \tag{34} $$

由 $softmax$ 定义可知 $$ a_i = \frac{e^{y_i}}{\sum_{k=1}^{n}e^{y_k}} \tag{35} $$ 对于 $\frac{\partial a_i}{\partial y_j}$，我们分两种情况讨论。

1. $j = i$，即 $\frac{\partial a_i}{\partial y_i}$

$$ \begin{aligned} \frac{\partial a_i}{\partial y_i} &= (\frac{e^{y_i}}{\sum_{k=1}^{n}e^{y_k}})' \\
&= \frac{e^{y_i}(\sum_{k=1}^{n}e^{y_k})-(e^{y_i})^2}{(\sum_{k=1}^{n}e^{y_k})^2} \\
&= \frac{e^{y_i}}{\sum_{k=1}^{n}e^{y_k}}-(\frac{e^{y_i}}{\sum_{k=1}^{n}e^{y_k}})^2 \\
&= a_i-a_i^2 \end{aligned} \tag{36} $$

2. $j \neq i$

$$ \begin{aligned} \frac{\partial a_i}{\partial y_j} &= (\frac{e^{y_i}}{\sum_{k=1}^{n}e^{y_k}})' \\
&= \frac{-e^{y_i}e^{y_j}}{(\sum_{k=1}^{n}e^{y_k})^2} \\
&= -\frac{e^{y_i}}{\sum_{k=1}^{n}e^{y_k}}\frac{e^{y_j}}{\sum_{k=1}^{n}e^{y_k}} \\
&= -a_ia_j \end{aligned} \tag{37} $$

因此对于 $-\frac{1}{a_i}\boldsymbol{v}_i$，当 $j=i$ 时，分量为 $a_j-1$；当 $j \neq i$ 时，分量为 $a_j$。又因为 $\boldsymbol{d}$ 是 one-hot 编码的向量，有 $$ d_j = \begin{cases} \begin{aligned} 1 &,& \text{j} = \text{i} \\
0 &,& \text{j} \neq \text{i} \end{aligned} \end{cases} \tag{38} $$ 我们可以得到最终的表示 $$ \frac{\partial E}{\partial \boldsymbol{y}} = -\frac{1}{a_i}\boldsymbol{v}_i = \boldsymbol{a} - \boldsymbol{d} \tag{39} $$ 可以看到，结果异常地简洁优雅。最后一层的误差向量可以根据 $(29)$ 很容易地计算得到 $$ \begin{aligned} \boldsymbol{\delta}^{(L)} &= \frac{\partial E}{\partial \boldsymbol{y}} \odot\phi^{'} (\boldsymbol{v} ^ {(L)}) \\
&= (\boldsymbol{a} - \boldsymbol{d}) \odot\phi^{'} (\boldsymbol{v} ^ {(L)}) \end{aligned} \tag{40} $$ 极具数学之美，听懂掌声。