矩阵求导（下）——矩阵对矩阵的求导

参考：https://zhuanlan.zhihu.com/p/24863977

本篇使用小写字母x表示标量，粗体小写字母$\mathbf{x}$表示列向量，大写字母X表示矩阵。
矩阵对矩阵的求导采用了向量化的思路，常应用于二阶方法求解优化问题。

首先来琢磨一下定义。矩阵对矩阵的导数，需要什么样的定义？
第一，矩阵$F(p\times q)$对矩阵$X(m\times n)$的导数应包含所有mnpq个偏导数$\frac{\partial F_{kl}}{\partial X_{ij}}$，从而不损失信息；
第二，导数与微分有简明的联系，因为在计算导数和应用中需要这个联系；
第三，导数有简明的从整体出发的算法。
我们先定义向量$\mathbf{f}(p\times 1)$对向量$\mathbf{x}(m\times 1$)的导数$\frac{\partial \mathbf{f}}{\partial \mathbf{x}}=\left[\begin{array}{cccc}\frac{\partial f_1}{\partial x_1}& \frac{\partial f_2}{\partial x_1}& \cdots & \frac{\partial f_p}{\partial x_1}\\ \frac{\partial f_1}{\partial x_2}& \frac{\partial f_2}{\partial x_2}& \cdots & \frac{\partial f_p}{\partial x_2}\\ ⋮& ⋮& \ddots & ⋮\\ \frac{\partial f_1}{\partial x_m}& \frac{\partial f_2}{\partial x_m}& \cdots & \frac{\partial f_p}{\partial x_m}\end{array}\right](m\times p)$，有$d\mathbf{f}={\frac{\partial \mathbf{f}}{\partial \mathbf{x}}}^{T}d\mathbf{x}$ ；
再定义矩阵的（按列优先）向量化$\mathrm{v}\mathrm{e}\mathrm{c}(X)=[X_{11},\dots ,X_{m1},X_{12},\dots ,X_{m2},\dots ,X_{1n},\dots ,X_{mn}{]}^T(mn\times 1)$，
并定义矩阵F对矩阵X的导数$\frac{\partial F}{\partial X}=\frac{\partial \mathrm{v}\mathrm{e}\mathrm{c}(F)}{\partial \mathrm{v}\mathrm{e}\mathrm{c}(X)}(mn\times pq)$。

导数与微分有联系$\mathrm{v}\mathrm{e}\mathrm{c}(dF)={\frac{\partial F}{\partial X}}^T\mathrm{v}\mathrm{e}\mathrm{c}(dX)$。几点说明如下：

1. 按此定义，标量f对矩阵X(m×n)的导数$\frac{\partial f}{\partial X}是mn\times 1$向量，与上篇的定义不兼容，不过二者容易相互转换。为避免混淆，用记号${\nabla }_{X}f$表示上篇定义的m×n矩阵，则有$\frac{\partial f}{\partial X}=\mathrm{v}\mathrm{e}\mathrm{c}({\nabla }_{X}f)$。虽然本篇的技术可以用于标量对矩阵求导这种特殊情况，但使用上篇中的技术更方便。读者可以通过上篇中的算例试验两种方法的等价转换。
2. 标量对矩阵的二阶导数，又称Hessian矩阵，定义为${\nabla }_X^{2}f=\frac{{\partial }^{2}f}{\partial X^2}=\frac{\partial {\nabla }_{X}f}{\partial X}(mn\times mn)$，是对称矩阵。对向量$\frac{\partial f}{\partial X}$或矩阵${\nabla }_{X}f$求导都可以得到Hessian矩阵，但从矩阵${\nabla }_{X}f$出发更方便。
3. $\frac{\partial F}{\partial X}=\frac{\partial \mathrm{v}\mathrm{e}\mathrm{c}(F)}{\partial X}=\frac{\partial F}{\partial \mathrm{v}\mathrm{e}\mathrm{c}(X)}=\frac{\partial \mathrm{v}\mathrm{e}\mathrm{c}(F)}{\partial \mathrm{v}\mathrm{e}\mathrm{c}(X)}$，求导时矩阵被向量化，弊端是这在一定程度破坏了矩阵的结构，会导致结果变得形式复杂；好处是多元微积分中关于梯度、Hessian矩阵的结论可以沿用过来，只需将矩阵向量化。例如优化问题中，牛顿法的更新$\Delta X$，满足$\mathrm{v}\mathrm{e}\mathrm{c}(\Delta X)=-({\nabla }_X^{2}f{)}^{-1}\mathrm{v}\mathrm{e}\mathrm{c}({\nabla }_{X}f)$。
4. 在资料中，矩阵对矩阵的导数还有其它定义，比如$\frac{\partial F}{\partial X}=\left[\frac{\partial F_{kl}}{\partial X}\right](mp\times nq)$，它能兼容上篇中的标量对矩阵导数的定义，但微分与导数的联系（$dF等于\frac{\partial F}{\partial X}中每个m\times n子块分别与dX做内积$）不够简明，不便于计算和应用。

然后来建立运算法则。仍然要利用导数与微分的联系$\mathrm{v}\mathrm{e}\mathrm{c}(dF)={\frac{\partial F}{\partial X}}^T\mathrm{v}\mathrm{e}\mathrm{c}(dX)$，求微分的方法与上篇相同，而从微分得到导数需要一些向量化的技巧：

• 线性：$\mathrm{v}\mathrm{e}\mathrm{c}(A+B)=\mathrm{v}\mathrm{e}\mathrm{c}(A)+\mathrm{v}\mathrm{e}\mathrm{c}(B)$。
• 矩阵乘法：$\mathrm{v}\mathrm{e}\mathrm{c}(AXB)=(B^T\otimes A)\mathrm{v}\mathrm{e}\mathrm{c}(X)$，其中$\otimes$表示Kronecker积，$A(m\times n)与B(p\times q)$的Kronecker积是$A\otimes B=[A_{ij}B](mp\times nq)$。此式证明见张贤达《矩阵分析与应用》第107-108页。
• 转置：$\mathrm{v}\mathrm{e}\mathrm{c}(A^T)=K_{mn}\mathrm{v}\mathrm{e}\mathrm{c}(A)$，A是m×n矩阵，其中$K_{mn}(mn\times mn)$是交换矩阵(commutation matrix)。
• 逐元素乘法：$\mathrm{v}\mathrm{e}\mathrm{c}(A\odot X)=\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}(A)\mathrm{v}\mathrm{e}\mathrm{c}(X)$，其中$\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}(A)(mn\times mn)$是用A的元素（按列优先）排成的对角阵。

观察一下可以断言，若矩阵函数F是矩阵X经加减乘法、行列式、逆、逐元素函数等运算构成，则使用相应的运算法则对F求微分，再做向量化并使用技巧将其它项交换至vec(dX)左侧，即能得到导数。

再谈一谈复合：假设已求得$\frac{\partial F}{\partial Y}$，而Y是X的函数，如何求$\frac{\partial F}{\partial X}$呢？从导数与微分的联系入手，$\mathrm{v}\mathrm{e}\mathrm{c}(dF)={\frac{\partial F}{\partial Y}}^T\mathrm{v}\mathrm{e}\mathrm{c}(dY)={\frac{\partial F}{\partial Y}}^T{\frac{\partial Y}{\partial X}}^T\mathrm{v}\mathrm{e}\mathrm{c}(dX)$ ，可以推出链式法则$\frac{\partial F}{\partial X}=\frac{\partial Y}{\partial X}\frac{\partial F}{\partial Y}$。

和标量对矩阵的导数相比，矩阵对矩阵的导数形式更加复杂，从不同角度出发常会得到形式不同的结果。有一些Kronecker积和交换矩阵相关的恒等式，可用来做等价变形：

1. $(A\otimes B{)}^T=A^T\otimes B^T$。
2. $\mathrm{v}\mathrm{e}\mathrm{c}({\mathbf{a}\mathbf{b}}^T)=\mathbf{b}\otimes \mathbf{a}$。
3. $(A\otimes B)(C\otimes D)=(AC)\otimes (BD)$。可以对$F=D^T{B}^{T}XAC$求导来证明，一方面，直接求导得到$\frac{\partial F}{\partial X}=(AC)\otimes (BD)$；另一方面，引入$Y=B^{T}XA，有\frac{\partial F}{\partial Y}=C\otimes D,\frac{\partial Y}{\partial X}=A\otimes B$，用链式法则得到$\frac{\partial F}{\partial X}=(A\otimes B)(C\otimes D)$。
4. $K_{mn}=K_{nm}^T,K_{mn}{K}_{nm}=I$。
5. $K_{pm}(A\otimes B)K_{nq}=B\otimes A$，A是m×n矩阵，B是p×q矩阵。可以对$AX{B}^T$做向量化来证明，一方面，$\mathrm{v}\mathrm{e}\mathrm{c}(AX{B}^T)=(B\otimes A)\mathrm{v}\mathrm{e}\mathrm{c}(X)$；另一方面，$\mathrm{v}\mathrm{e}\mathrm{c}(AX{B}^T)=K_{pm}\mathrm{v}\mathrm{e}\mathrm{c}(B{X}^T{A}^T)=K_{pm}(A\otimes B)\mathrm{v}\mathrm{e}\mathrm{c}(X^T)=K_{pm}(A\otimes B)K_{nq}\mathrm{v}\mathrm{e}\mathrm{c}(X)$。

接下来演示一些算例。

• 例1：$F=AX$，X是m×n矩阵，求$\frac{\partial F}{\partial X}$。

  解：先求微分：$dF=AdX$，再做向量化，使用矩阵乘法的技巧，注意在dX右侧添加单位阵：$\mathrm{v}\mathrm{e}\mathrm{c}(dF)=\mathrm{v}\mathrm{e}\mathrm{c}(AdX)=(I_n\otimes A)\mathrm{v}\mathrm{e}\mathrm{c}(dX)$，对照导数与微分的联系得到$\frac{\partial F}{\partial X}=I_n\otimes A^T$。

• 特例：如果X退化为向量， $\mathbf{f}=A\mathbf{x}$ ，则根据向量的导数与微分的关系 $d\mathbf{f}={\frac{\partial \mathbf{f}}{\partial \mathbf{x}}}^{T}d\mathbf{x}$ ，得到 $\frac{\partial \mathbf{f}}{\partial \mathbf{x}}=A^T$ 。

• 例2：$f=\log |X|$ ，X是n×n矩阵，求${\nabla }_{X}f和{\nabla }_X^{2}f$。

  解：使用上篇中的技术可求得${\nabla }_{X}f=X^{-1T}$ 。
  为求${\nabla }_X^{2}f$，先求微分：$d{\nabla }_{X}f=-(X^{-1}dX{X}^{-1}{)}^T$，再做向量化，使用转置和矩阵乘法的技巧
  $\mathrm{v}\mathrm{e}\mathrm{c}(d{\nabla }_{X}f)=-K_{nn}\mathrm{v}\mathrm{e}\mathrm{c}(X^{-1}dX{X}^{-1})=-K_{nn}(X^{-1T}\otimes X^{-1})\mathrm{v}\mathrm{e}\mathrm{c}(dX)$，
  对照导数与微分的联系，得到${\nabla }_X^{2}f=-K_{nn}(X^{-1T}\otimes X^{-1})$，注意它是对称矩阵。
  在X是对称矩阵时，可简化为${\nabla }_X^{2}f=-X^{-1}\otimes X^{-1}$。

• 例3：$F=A\exp (XB)$，A是l×m，X是m×n，B是n×p矩阵，exp()为逐元素函数，求$\frac{\partial F}{\partial X}$。

  解：先求微分：$dF=A(\exp (XB)\odot (dXB))$，再做向量化，
  使用矩阵乘法的技巧：$\mathrm{v}\mathrm{e}\mathrm{c}(dF)=(I_p\otimes A)\mathrm{v}\mathrm{e}\mathrm{c}(\exp (XB)\odot (dXB))$，
  再用逐元素乘法的技巧：$\mathrm{v}\mathrm{e}\mathrm{c}(dF)=(I_p\otimes A)\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}(\exp (XB))\mathrm{v}\mathrm{e}\mathrm{c}(dXB)$，
  再用矩阵乘法的技巧：$\mathrm{v}\mathrm{e}\mathrm{c}(dF)=(I_p\otimes A)\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}(\exp (XB))(B^T\otimes I_m)\mathrm{v}\mathrm{e}\mathrm{c}(dX)$，
  对照导数与微分的联系得到$\frac{\partial F}{\partial X}=(B\otimes I_m)\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}(\exp (XB))(I_p\otimes A^T)$。

• 例4【一元logistic回归】：$l=-y{\mathbf{x}}^T\mathbf{w}+\log (1+\exp ({\mathbf{x}}^T\mathbf{w}))$，求${\nabla }_{\mathbf{w}}l$和￥${\nabla }_{\mathbf{w}}^{2}l$。其中$y$是取值0或1的标量，$\mathbf{x},\mathbf{w}$是向量。

  解：使用上篇中的技术可求得${\nabla }_{\mathbf{w}}l=\mathbf{x}(\sigma ({\mathbf{x}}^T\mathbf{w})-y)$，其中$\sigma (a)=\frac{\exp (a)}{1+\exp (a)}$ 为$sigmoid$函数。
  为求${\nabla }_{\mathbf{w}}^{2}l$，先求微分：$d{\nabla }_{\mathbf{w}}l$ =$\mathbf{x}{\sigma }^{\prime }({\mathbf{x}}^T\mathbf{w}){\mathbf{x}}^{T}d\mathbf{w}$,
  其中${\sigma }^{\prime }(a)=\frac{\exp (a)}{(1+\exp (a){)}^2}$为sigmoid函数的导数，
  对照导数与微分的联系，得到${\nabla }_w^{2}l=\mathbf{x}{\sigma }^{\prime }({\mathbf{x}}^T\mathbf{w}){\mathbf{x}}^T$。

• 推广：样本$({\mathbf{x}}_1,y_1),\dots ,({\mathbf{x}}_n,y_n)，l={\sum }_{i=1}^N\left(-y_i{\mathbf{x}}_i^T\mathbf{w}+\log (1+\exp ({{\mathbf{x}}_{\mathbf{i}}}^T\mathbf{w}))\right)$，求${\nabla }_{w}l和{\nabla }_w^{2}l$。
  有两种方法，
  方法一：先对每个样本求导，然后相加；
  方法二：定义矩阵$X=\left[\begin{array}{c}{\mathbf{x}}_1^T\\ ⋮\\ {\mathbf{x}}_n^T\end{array}\right]$，向量$\mathbf{y}=\left[\begin{array}{c}{y}_1\\ ⋮\\ y_n\end{array}\right]$，将$l$写成矩阵形式$l=-{\mathbf{y}}^{T}X\mathbf{w}+1^T\log (1+\exp (X\mathbf{w}))$，进而可以求得${\nabla }_{\mathbf{w}}l=X^T(\sigma (X\mathbf{w})-\mathbf{y})，{\nabla }_w^{2}l=X^T\text{diag}({\sigma }^{\prime }(X\mathbf{w}))X$。

• 例5【多元logistic回归】：$l=-{\mathbf{y}}^T\log \text{softmax}(W\mathbf{x})=-{\mathbf{y}}^{T}W\mathbf{x}+\log (1^T\exp (W\mathbf{x}))，求{\nabla }_{W}l和{\nabla }_W^{2}l$ 。

  解：上篇中已求得${\nabla }_{W}l=(\text{softmax}(W\mathbf{x})-\mathbf{y}){\mathbf{x}}^T$。为求${\nabla }_W^{2}l$，
  先求微分：定义$\mathbf{a}=W\mathbf{x}，d\text{softmax}(\mathbf{a})=\frac{\exp (\mathbf{a})\odot d\mathbf{a}}{1^T\exp (\mathbf{a})}-\frac{\exp (\mathbf{a})(1^T(\exp (\mathbf{a})\odot d\mathbf{a}))}{(1^T\exp (\mathbf{a}){)}^2}$，
  这里需要化简去掉逐元素乘法，第一项中$\exp (\mathbf{a})\odot d\mathbf{a}=\text{diag}(\exp (\mathbf{a}))d\mathbf{a}$ ，
  第二项中$1^T(\exp (\mathbf{a})\odot d\mathbf{a})=\exp (\mathbf{a}{)}^{T}d\mathbf{a}，故有d\text{softmax}(\mathbf{a})={\text{softmax}}^{\prime }(\mathbf{a})d\mathbf{a}$，
  其中${\text{softmax}}^{\prime }(\mathbf{a})=\frac{\text{diag}(\exp (\mathbf{a}))}{1^T\exp (\mathbf{a})}-\frac{\exp (\mathbf{a})\exp (\mathbf{a}{)}^T}{(1^T\exp (\mathbf{a}){)}^2}$ ，
  代入有$d{\nabla }_{W}l={\text{softmax}}^{\prime }(\mathbf{a})d\mathbf{a}{\mathbf{x}}^T={\text{softmax}}^{\prime }(W\mathbf{x})dW\mathbf{x}{\mathbf{x}}^T$，
  做向量化并使用矩阵乘法的技巧，得到${\nabla }_W^{2}l=(\mathbf{x}{\mathbf{x}}^T)\otimes {\text{softmax}}^{\prime }(W\mathbf{x})$。

总结

我们发展了从整体出发的矩阵求导的技术，导数与微分的联系是计算的枢纽，标量对矩阵的导数与微分的联系是
$df=\mathrm{t}\mathrm{r}({\nabla }_X^{T}fdX)$，先对f求微分，再使用迹技巧可求得导数，特别地，标量对向量的导数与微分的联系是
$df={\nabla }_{\mathbf{x}}^{T}fd\mathbf{x}；矩阵对矩阵的导数与微分的联系是\mathrm{v}\mathrm{e}\mathrm{c}(dF)={\frac{\partial F}{\partial X}}^T\mathrm{v}\mathrm{e}\mathrm{c}(dX)$，先对F求微分，再使用向量化的技巧可求得导数，特别地，向量对向量的导数与微分的联系是$d\mathbf{f}={\frac{\partial \mathbf{f}}{\partial \mathbf{x}}}^{T}d\mathbf{x}$。