参考:https://zhuanlan.zhihu.com/p/24863977
本篇使用小写字母x表示标量,粗体小写字母
x表示列向量,大写字母X表示矩阵。
矩阵对矩阵的求导采用了向量化的思路,常应用于二阶方法求解优化问题。
首先来琢磨一下定义。矩阵对矩阵的导数,需要什么样的定义?
第一,矩阵
F(p×q)对矩阵
X(m×n)的导数应包含所有mnpq个偏导数
∂Xij∂Fkl,从而不损失信息;
第二,导数与微分有简明的联系,因为在计算导数和应用中需要这个联系;
第三,导数有简明的从整体出发的算法。
我们先定义向量
f(p×1)对向量
x(m×1)的导数
∂x∂f=⎣⎢⎢⎢⎢⎡∂x1∂f1∂x2∂f1⋮∂xm∂f1∂x1∂f2∂x2∂f2⋮∂xm∂f2⋯⋯⋱⋯∂x1∂fp∂x2∂fp⋮∂xm∂fp⎦⎥⎥⎥⎥⎤(m×p),有
df=∂x∂fTdx ;
再定义矩阵的(按列优先)向量化
vec(X)=[X11,…,Xm1,X12,…,Xm2,…,X1n,…,Xmn]T(mn×1),
并定义矩阵F对矩阵X的导数
∂X∂F=∂vec(X)∂vec(F)(mn×pq)。
导数与微分有联系
vec(dF)=∂X∂FTvec(dX)。几点说明如下:
- 按此定义,标量f对矩阵X(m×n)的导数
∂X∂f是mn×1向量,与上篇的定义不兼容,不过二者容易相互转换。为避免混淆,用记号
∇Xf表示上篇定义的m×n矩阵,则有
∂X∂f=vec(∇Xf)。虽然本篇的技术可以用于标量对矩阵求导这种特殊情况,但使用上篇中的技术更方便。读者可以通过上篇中的算例试验两种方法的等价转换。
- 标量对矩阵的二阶导数,又称Hessian矩阵,定义为
∇X2f=∂X2∂2f=∂X∂∇Xf(mn×mn),是对称矩阵。对向量
∂X∂f或矩阵
∇Xf求导都可以得到Hessian矩阵,但从矩阵
∇Xf出发更方便。
-
∂X∂F=∂X∂vec(F)=∂vec(X)∂F=∂vec(X)∂vec(F),求导时矩阵被向量化,弊端是这在一定程度破坏了矩阵的结构,会导致结果变得形式复杂;好处是多元微积分中关于梯度、Hessian矩阵的结论可以沿用过来,只需将矩阵向量化。例如优化问题中,牛顿法的更新
ΔX,满足
vec(ΔX)=−(∇X2f)−1vec(∇Xf)。
- 在资料中,矩阵对矩阵的导数还有其它定义,比如
∂X∂F=[∂X∂Fkl](mp×nq),它能兼容上篇中的标量对矩阵导数的定义,但微分与导数的联系(
dF等于∂X∂F中每个m×n子块分别与dX做内积)不够简明,不便于计算和应用。
然后来建立运算法则。仍然要利用导数与微分的联系
vec(dF)=∂X∂FTvec(dX),求微分的方法与上篇相同,而从微分得到导数需要一些向量化的技巧:
- 线性:
vec(A+B)=vec(A)+vec(B)。
- 矩阵乘法:
vec(AXB)=(BT⊗A)vec(X),其中
⊗表示Kronecker积,
A(m×n)与B(p×q)的Kronecker积是
A⊗B=[AijB](mp×nq)。此式证明见张贤达《矩阵分析与应用》第107-108页。
- 转置:
vec(AT)=Kmnvec(A),A是m×n矩阵,其中
Kmn(mn×mn)是交换矩阵(commutation matrix)。
- 逐元素乘法:
vec(A⊙X)=diag(A)vec(X),其中
diag(A)(mn×mn)是用A的元素(按列优先)排成的对角阵。
观察一下可以断言,若矩阵函数F是矩阵X经加减乘法、行列式、逆、逐元素函数等运算构成,则使用相应的运算法则对F求微分,再做向量化并使用技巧将其它项交换至vec(dX)左侧,即能得到导数。
再谈一谈复合:假设已求得
∂Y∂F,而Y是X的函数,如何求
∂X∂F呢?从导数与微分的联系入手,
vec(dF)=∂Y∂FTvec(dY)=∂Y∂FT∂X∂YTvec(dX) ,可以推出链式法则
∂X∂F=∂X∂Y∂Y∂F。
和标量对矩阵的导数相比,矩阵对矩阵的导数形式更加复杂,从不同角度出发常会得到形式不同的结果。有一些Kronecker积和交换矩阵相关的恒等式,可用来做等价变形:
-
(A⊗B)T=AT⊗BT。
-
vec(abT)=b⊗a。
-
(A⊗B)(C⊗D)=(AC)⊗(BD)。可以对
F=DTBTXAC求导来证明,一方面,直接求导得到
∂X∂F=(AC)⊗(BD);另一方面,引入
Y=BTXA,有∂Y∂F=C⊗D,∂X∂Y=A⊗B,用链式法则得到
∂X∂F=(A⊗B)(C⊗D)。
-
Kmn=KnmT,KmnKnm=I。
-
Kpm(A⊗B)Knq=B⊗A,A是m×n矩阵,B是p×q矩阵。可以对
AXBT做向量化来证明,一方面,
vec(AXBT)=(B⊗A)vec(X);另一方面,
vec(AXBT)=Kpmvec(BXTAT)=Kpm(A⊗B)vec(XT)=Kpm(A⊗B)Knqvec(X)。
接下来演示一些算例。
-
例1:
F=AX,X是m×n矩阵,求
∂X∂F。
解:先求微分:
dF=AdX,再做向量化,使用矩阵乘法的技巧,注意在dX右侧添加单位阵:
vec(dF)=vec(AdX)=(In⊗A)vec(dX),对照导数与微分的联系得到
∂X∂F=In⊗AT。
-
特例:如果X退化为向量,
f=Ax ,则根据向量的导数与微分的关系
df=∂x∂fTdx ,得到
∂x∂f=AT 。
-
例2:
f=log∣X∣ ,X是n×n矩阵,求
∇Xf和∇X2f。
解:使用上篇中的技术可求得
∇Xf=X−1T 。
为求
∇X2f,先求微分:
d∇Xf=−(X−1dXX−1)T,再做向量化,使用转置和矩阵乘法的技巧
vec(d∇Xf)=−Knnvec(X−1dXX−1)=−Knn(X−1T⊗X−1)vec(dX),
对照导数与微分的联系,得到
∇X2f=−Knn(X−1T⊗X−1),注意它是对称矩阵。
在X是对称矩阵时,可简化为
∇X2f=−X−1⊗X−1。
-
例3:
F=Aexp(XB),A是l×m,X是m×n,B是n×p矩阵,exp()为逐元素函数,求
∂X∂F。
解:先求微分:
dF=A(exp(XB)⊙(dXB)),再做向量化,
使用矩阵乘法的技巧:
vec(dF)=(Ip⊗A)vec(exp(XB)⊙(dXB)),
再用逐元素乘法的技巧:
vec(dF)=(Ip⊗A)diag(exp(XB))vec(dXB),
再用矩阵乘法的技巧:
vec(dF)=(Ip⊗A)diag(exp(XB))(BT⊗Im)vec(dX),
对照导数与微分的联系得到
∂X∂F=(B⊗Im)diag(exp(XB))(Ip⊗AT)。
-
例4【一元logistic回归】:
l=−yxTw+log(1+exp(xTw)),求
∇wl和¥
∇w2l。其中
y是取值0或1的标量,
x,w是向量。
解:使用上篇中的技术可求得
∇wl=x(σ(xTw)−y),其中
σ(a)=1+exp(a)exp(a) 为
sigmoid函数。
为求
∇w2l,先求微分:
d∇wl =
xσ′(xTw)xTdw,
其中
σ′(a)=(1+exp(a))2exp(a)为sigmoid函数的导数,
对照导数与微分的联系,得到
∇w2l=xσ′(xTw)xT。
-
推广:样本
(x1,y1),…,(xn,yn),l=∑i=1N(−yixiTw+log(1+exp(xiTw))),求
∇wl和∇w2l。
有两种方法,
方法一:先对每个样本求导,然后相加;
方法二:定义矩阵
X=⎣⎢⎡x1T⋮xnT⎦⎥⎤,向量
y=⎣⎢⎡y1⋮yn⎦⎥⎤,将
l写成矩阵形式
l=−yTXw+1Tlog(1+exp(Xw)),进而可以求得
∇wl=XT(σ(Xw)−y),∇w2l=XTdiag(σ′(Xw))X。
-
例5【多元logistic回归】:
l=−yTlogsoftmax(Wx)=−yTWx+log(1Texp(Wx)),求∇Wl和∇W2l 。
解:上篇中已求得
∇Wl=(softmax(Wx)−y)xT。为求
∇W2l,
先求微分:定义
a=Wx,dsoftmax(a)=1Texp(a)exp(a)⊙da−(1Texp(a))2exp(a)(1T(exp(a)⊙da)),
这里需要化简去掉逐元素乘法,第一项中
exp(a)⊙da=diag(exp(a))da ,
第二项中
1T(exp(a)⊙da)=exp(a)Tda,故有dsoftmax(a)=softmax′(a)da,
其中
softmax′(a)=1Texp(a)diag(exp(a))−(1Texp(a))2exp(a)exp(a)T ,
代入有
d∇Wl=softmax′(a)daxT=softmax′(Wx)dWxxT,
做向量化并使用矩阵乘法的技巧,得到
∇W2l=(xxT)⊗softmax′(Wx)。
总结
我们发展了从整体出发的矩阵求导的技术,导数与微分的联系是计算的枢纽,标量对矩阵的导数与微分的联系是
df=tr(∇XTfdX),先对f求微分,再使用迹技巧可求得导数,特别地,标量对向量的导数与微分的联系是
df=∇xTfdx;矩阵对矩阵的导数与微分的联系是vec(dF)=∂X∂FTvec(dX),先对F求微分,再使用向量化的技巧可求得导数,特别地,向量对向量的导数与微分的联系是
df=∂x∂fTdx。