『科学计算_理论』矩阵求导

上篇原文链接

下篇原文链接

标量对矩阵求导

矩阵求导的技术,在统计学、控制论、机器学习等领域有广泛的应用。鉴于我看过的一些资料或言之不详、或繁乱无绪,本文来做个科普,分作两篇,上篇讲标量对矩阵的求导术,下篇讲矩阵对矩阵的求导术。本文使用小写字母x表示标量,粗体小写字母oldsymbol{x} 表示向量,大写字母X表示矩阵。

首先来琢磨一下定义,标量f对矩阵X的导数,定义为frac{partial f}{partial X} := left[frac{partial f }{partial X_{ij}}
ight],即f对X逐元素求导排成与X尺寸相同的矩阵。然而,这个定义在计算中并不好用,实用上的原因是在对较复杂的函数难以逐元素求导;哲理上的原因是逐元素求导破坏了整体性。试想,为何要将f看做矩阵X而不是各元素X_{ij}的函数呢?答案是用矩阵运算更整洁。所以在求导时不宜拆开矩阵,而是要找一个从整体出发的算法。为此,我们来回顾,一元微积分中的导数(标量对标量的导数)与微分有联系:df = f'(x)dx;多元微积分中的梯度(标量对向量的导数)也与微分有联系:df = sum_{i} frac{partial f}{partial x_i}dx_i =  frac{partial f}{partial oldsymbol{x}}^T doldsymbol{x} ,这里第一个等号是全微分公式,第二个等号表达了梯度frac{partial f}{partial oldsymbol{x}}与微分的联系;受此启发,我们将矩阵导数与微分建立联系:df = sum_{i,j} frac{partial f}{partial X_{ij}}dX_{ij} = 	ext{tr}left(frac{partial f}{partial X}^T dX
ight) ,这里tr代表迹(trace)是方阵对角线元素之和,满足性质:对尺寸相同的矩阵A,B,	ext{tr}(A^TB) = sum_{i,j}A_{ij}B_{ij},这用泛函分析的语言来说	ext{tr}(A^TB)是矩阵A,B的内积,因此上式与原定义相容。

然后来建立运算法则。回想遇到较复杂的一元函数如f = log(2+sin x)e^{sqrt{x}},我们是如何求导的呢?通常不是从定义开始求极限,而是先建立了初等函数求导和四则运算、复合等法则,再来运用这些法则。故而,我们来创立常用的矩阵微分的运算法则:

  1. 加减法:d(Xpm Y) = dX pm dY;矩阵乘法:d(XY) = dX Y + X dY ;转置:d(X^T) = (dX)^T;迹:d	ext{tr}(X) = 	ext{tr}(dX)
  2. 逆:dX^{-1} = -X^{-1}dX X^{-1}。此式可在XX^{-1}=I两侧求微分来证明。
  3. 行列式:d|X| = 	ext{tr}(X^{#}dX) ,其中X^{#}表示X的伴随矩阵,在X可逆时又可以写作d|X|= |X|	ext{tr}(X^{-1}dX)。此式可用Laplace展开来证明,详见张贤达《矩阵分析与应用》第279页。
  4. 逐元素乘法:d(Xodot Y) = dXodot Y + Xodot dYodot表示尺寸相同的矩阵X,Y逐元素相乘。
  5. 逐元素函数:dsigma(X) = sigma'(X)odot dX sigma(X) = left[sigma(X_{ij})
ight]是逐元素运算的标量函数。



我们试图利用矩阵导数与微分的联系df = 	ext{tr}left(frac{partial f}{partial X}^T dX
ight) ,在求出左侧的微分df后,该如何写成右侧的形式并得到导数呢?这需要一些迹技巧(trace trick):

  1. 标量套上迹:a = 	ext{tr}(a)
  2. 转置:mathrm{tr}(A^T) = mathrm{tr}(A)
  3. 线性:	ext{tr}(Apm B) = 	ext{tr}(A)pm 	ext{tr}(B)
  4. 矩阵乘法交换:	ext{tr}(AB) = 	ext{tr}(BA)。两侧都等于sum_{i,j}A_{ij}B_{ji}
  5. 矩阵乘法/逐元素乘法交换:	ext{tr}(A^T(Bodot C)) = 	ext{tr}((Aodot B)^TC)。两侧都等于sum_{i,j}A_{ij}B_{ij}C_{ij}

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


在建立法则的最后,来谈一谈复合:假设已求得frac{partial f}{partial Y},而Y是X的函数,如何求frac{partial f}{partial X}呢?在微积分中有标量求导的链式法则frac{partial f}{partial x} = frac{partial f}{partial y} frac{partial y}{partial x},但这里我们不能沿用链式法则,因为矩阵对矩阵的导数frac{partial Y}{partial X}截至目前仍是未定义的。于是我们继续追本溯源,链式法则是从何而来?源头仍然是微分。我们直接从微分入手建立复合法则:先写出df = 	ext{tr}left(frac{partial f}{partial Y}^T dY
ight),再将dY用dX表示出来代入,并使用迹技巧将其他项交换至dX左侧,即可得到frac{partial f}{partial X}

接下来演示一些算例。特别提醒要依据已经建立的运算法则来计算,不能随意套用微积分中标量导数的结论,比如认为AX对X的导数为A,这是没有根据、意义不明的。

例1:f = oldsymbol{a}^T Xoldsymbol{b},求frac{partial f}{partial X}

解:先使用矩阵乘法法则求微分:df =  oldsymbol{a}^T dXoldsymbol{b} ,再套上迹并做交换:df = 	ext{tr}(oldsymbol{a}^TdXoldsymbol{b}) = 	ext{tr}(oldsymbol{b}oldsymbol{a}^TdX),对照导数与微分的联系,得到frac{partial f}{partial X} =  oldsymbol{a}oldsymbol{b}^T

注意:这里不能用frac{partial f}{partial X} =oldsymbol{a}^T frac{partial X}{partial X}oldsymbol{b}=?,导数与乘常数矩阵的交换是不合法则的运算(而微分是合法的)。有些资料在计算矩阵导数时,会略过求微分这一步,这是逻辑上解释不通的。

例2【线性回归】:l = |Xoldsymbol{w}-  oldsymbol{y}|^2,求frac{partial l}{partial oldsymbol{w}}

解:严格来说这是标量对向量的导数,不过可以把向量看做矩阵的特例。将向量范数写成l = (Xoldsymbol{w}-  oldsymbol{y})^T(Xoldsymbol{w}-  oldsymbol{y}),求微分,使用矩阵乘法、转置等法则:dl = (Xdoldsymbol{w})^T(Xoldsymbol{w}-oldsymbol{y})+(Xoldsymbol{w}-oldsymbol{y})^T(Xdoldsymbol{w}) = 2(Xoldsymbol{w}-oldsymbol{y})^TXdoldsymbol{w}。对照导数与微分的联系,得到frac{partial l}{partial oldsymbol{w}}= 2X^T(Xoldsymbol{w}-oldsymbol{y})

例3【多元logistic回归】:l = -oldsymbol{y}^Tlog	ext{softmax}(Woldsymbol{x}),求frac{partial l}{partial W}。其中oldsymbol{y}是除一个元素为1外其它元素为0的向量;	ext{softmax}(oldsymbol{a}) = frac{exp(oldsymbol{a})}{oldsymbol{1}^Texp(oldsymbol{a})},其中exp(oldsymbol{a})表示逐元素求指数,oldsymbol{1}代表全1向量。

解:首先将softmax函数代入并写成l = -oldsymbol{y}^T left(log (exp(Woldsymbol{x}))-oldsymbol{1}log(oldsymbol{1}^Texp(Woldsymbol{x}))
ight) = -oldsymbol{y}^TWoldsymbol{x} + log(oldsymbol{1}^Texp(Woldsymbol{x})),这里要注意向量除标量求逐元素log满足log(oldsymbol{b}/c) = log(oldsymbol{b}) - oldsymbol{1}log(c),以及oldsymbol{y}满足oldsymbol{y}^T oldsymbol{1} = 1。求微分,使用矩阵乘法、逐元素函数等法则:dl =- oldsymbol{y}^TdWoldsymbol{x}+frac{oldsymbol{1}^Tleft(exp(Woldsymbol{x})odot(dWoldsymbol{x})
ight)}{oldsymbol{1}^Texp(Woldsymbol{x})}。再套上迹并做交换,其中第二项的分子是	ext{tr}left(oldsymbol{1}^Tleft(exp(Woldsymbol{x})odot(dWoldsymbol{x})
ight)
ight) = 	ext{tr}left((oldsymbol{1}odotexp(Woldsymbol{x}))^TdWoldsymbol{x}
ight) = 	ext{tr}(exp(Woldsymbol{x})^TdWoldsymbol{x}),故dl = 	ext{tr}left(-oldsymbol{y}^TdWoldsymbol{x}+frac{exp(Woldsymbol{x})^TdWoldsymbol{x}}{oldsymbol{1}^Texp(Woldsymbol{x})}
ight) =	ext{tr}(oldsymbol{x}(	ext{softmax}(Woldsymbol{x})-oldsymbol{y})^TdW)。对照导数与微分的联系,得到frac{partial l}{partial W}= (	ext{softmax}(Woldsymbol{x})-oldsymbol{y})oldsymbol{x}^T

另解:定义oldsymbol{a} = Woldsymbol{x},则l = -oldsymbol{y}^Tlog	ext{softmax}(oldsymbol{a}) ,先如上求出frac{partial l}{partial oldsymbol{a}} = 	ext{softmax}(oldsymbol{a})-oldsymbol{y} ,再利用复合法则:dl = 	ext{tr}left(frac{partial l}{partial oldsymbol{a}}^Tdoldsymbol{a}
ight) = 	ext{tr}left(frac{partial l}{partial oldsymbol{a}}^TdW oldsymbol{x}
ight) = 	ext{tr}left(oldsymbol{x}frac{partial l}{partial oldsymbol{a}}^TdW
ight),得到frac{partial l}{partial W}= frac{partial l}{partialoldsymbol{a}}oldsymbol{x}^T

例4【方差的最大似然估计】:样本oldsymbol{x}_1,dots, oldsymbol{x}_nsim N(oldsymbol{mu}, Sigma),其中Sigma是对称正定矩阵,求方差Sigma的最大似然估计。写成数学式是:l = log|Sigma|+frac{1}{n}sum_{i=1}^n(oldsymbol{x}_i-oldsymbol{ar{x}})^TSigma^{-1}(oldsymbol{x}_i-oldsymbol{ar{x}}),求frac{partial l }{partial Sigma}的零点。

解:首先求微分,使用矩阵乘法、行列式、逆等运算法则,第一项是dlog|Sigma| = |Sigma|^{-1}d|Sigma| = 	ext{tr}(Sigma^{-1}dSigma),第二项是frac{1}{n}sum_{i=1}^n(oldsymbol{x}_i-oldsymbol{ar{x}})^TdSigma^{-1}(oldsymbol{x}_i-oldsymbol{ar{x}}) = -frac{1}{n}sum_{i=1}^n(oldsymbol{x}_i-oldsymbol{ar{x}})^TSigma^{-1}dSigmaSigma^{-1}(oldsymbol{x}_i-oldsymbol{ar{x}})。再给第二项套上迹做交换:dl = 	ext{tr}left(left(Sigma^{-1}-Sigma^{-1}S_nSigma^{-1}
ight)dSigma
ight),其中S_n := frac{1}{n}sum_{i=1}^n(oldsymbol{x}_i-oldsymbol{ar{x}})(oldsymbol{x}_i-oldsymbol{ar{x}})^T定义为样本方差。对照导数与微分的联系,有frac{partial l }{partial Sigma}=(Sigma^{-1}-Sigma^{-1}S_nSigma^{-1})^T,其零点即Sigma的最大似然估计为Sigma = S_n

最后一例留给经典的神经网络。神经网络的求导术是学术史上的重要成果,还有个专门的名字叫做BP算法,我相信如今很多人在初次推导BP算法时也会颇费一番脑筋,事实上使用矩阵求导术来推导并不复杂。为简化起见,我们推导二层神经网络的BP算法。

例5【二层神经网络】:l = -oldsymbol{y}^Tlog	ext{softmax}(W_2sigma(W_1oldsymbol{x})),求frac{partial l}{partial W_1}frac{partial l}{partial W_2}。其中oldsymbol{y}是除一个元素为1外其它元素为0的向量,	ext{softmax}(oldsymbol{a}) = frac{exp(oldsymbol{a})}{oldsymbol{1}^Texp(oldsymbol{a})}同例3,sigma(cdot)是逐元素sigmoid函数sigma(a) = frac{1}{1+exp(-a)}

解:定义oldsymbol{a}_1=W_1oldsymbol{x}oldsymbol{h}_1 = sigma(oldsymbol{a}_1)oldsymbol{a}_2 = W_2 oldsymbol{h}_1,则l =-oldsymbol{y}^Tlog	ext{softmax}(oldsymbol{a}_2)。在例3中已求出frac{partial l}{partial oldsymbol{a}_2} = 	ext{softmax}(oldsymbol{a}_2)-oldsymbol{y} 。使用复合法则,注意此处oldsymbol{h}_1, W_2都是变量:dl = 	ext{tr}left(frac{partial l}{partial oldsymbol{a}_2}^Tdoldsymbol{a}_2
ight) = 	ext{tr}left(frac{partial l}{partial oldsymbol{a}_2}^TdW_2 oldsymbol{h}_1
ight) + 	ext{tr}left(frac{partial l}{partial oldsymbol{a}_2}^TW_2 doldsymbol{h}_1
ight),使用矩阵乘法交换的迹技巧从第一项得到frac{partial l}{partial W_2}= frac{partial l}{partialoldsymbol{a}_2}oldsymbol{h}_1^T,从第二项得到frac{partial l}{partial oldsymbol{h}_1}= W_2^Tfrac{partial l}{partialoldsymbol{a}_2}。接下来求frac{partial l}{partial oldsymbol{a}_1},继续使用复合法则,并利用矩阵乘法和逐元素乘法交换的迹技巧:	ext{tr}left(frac{partial l}{partialoldsymbol{h}_1}^Tdoldsymbol{h}_1
ight) = 	ext{tr}left(frac{partial l}{partialoldsymbol{h}_1}^T(sigma'(oldsymbol{a}_1)odot doldsymbol{a}_1)
ight) = 	ext{tr}left(left(frac{partial l}{partialoldsymbol{h}_1}odot sigma'(oldsymbol{a}_1)
ight)^Tdoldsymbol{a}_1
ight),得到frac{partial l}{partial oldsymbol{a}_1}= frac{partial l}{partialoldsymbol{h}_1}odotsigma'(oldsymbol{a}_1)。为求frac{partial l}{partial W_1},再用一次复合法则:	ext{tr}left(frac{partial l}{partialoldsymbol{a}_1}^Tdoldsymbol{a}_1
ight) = 	ext{tr}left(frac{partial l}{partialoldsymbol{a}_1}^TdW_1oldsymbol{x}
ight) = 	ext{tr}left(oldsymbol{x}frac{partial l}{partialoldsymbol{a}_1}^TdW_1
ight),得到frac{partial l}{partial W_1}= frac{partial l}{partialoldsymbol{a}_1}oldsymbol{x}^T

矩阵对矩阵求导

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

首先来琢磨一下定义。矩阵对矩阵的导数,需要什么样的定义?第一,矩阵F(p×q)对矩阵X(m×n)的导数应包含所有mnpq个偏导数frac{partial F_{kl}}{partial X_{ij}},从而不损失信息;第二,导数与微分有简明的联系,因为在计算导数和应用中需要这个联系;第三,导数有简明的从整体出发的算法。我们先定义向量oldsymbol{f}(p×1)对向量oldsymbol{x}(m×1)的导数frac{partial oldsymbol{f}}{partial oldsymbol{x}} := 
egin{bmatrix}
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}\
vdots & vdots & ddots & vdots\
frac{partial f_1}{partial x_m} & frac{partial f_2}{partial x_m} & cdots & frac{partial f_p}{partial x_m}\
end{bmatrix}(m×p),有doldsymbol{f} = frac{partial oldsymbol{f} }{partial oldsymbol{x} }^T doldsymbol{x} ;再定义矩阵的(按列优先)向量化mathrm{vec}(X) := [X_{11}, ldots, X_{m1}, X_{12}, ldots, X_{m2}, ldots, X_{1n}, ldots, X_{mn}]^T(mn×1),并定义矩阵F对矩阵X的导数frac{partial F}{partial X} := frac{partial mathrm{vec}(F)}{partial mathrm{vec}(X)}(mn×pq)。导数与微分有联系mathrm{vec}(dF) = frac{partial F}{partial X}^T mathrm{vec}(dX)。几点说明如下:

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



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

  1. 线性:mathrm{vec}(A+B) = mathrm{vec}(A) + mathrm{vec}(B)
  2. 矩阵乘法:mathrm{vec}(AXB) = (B^T otimes A) mathrm{vec}(X),其中otimes表示Kronecker积,A(m×n)与B(p×q)的Kronecker积是Aotimes B := [A_{ij}B](mp×nq)。此式证明见张贤达《矩阵分析与应用》第107-108页。
  3. 转置:mathrm{vec}(A^T) = K_{mn}mathrm{vec}(A),A是m×n矩阵,其中K_{mn}(mn×mn)是交换矩阵(commutation matrix)。
  4. 逐元素乘法:mathrm{vec}(Aodot X) = mathrm{diag}(A)mathrm{vec}(X),其中mathrm{diag}(A)(mn×mn)是用A的元素(按列优先)排成的对角阵。

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

再谈一谈复合:假设已求得frac{partial F}{partial Y},而Y是X的函数,如何求frac{partial F}{partial X}呢?从导数与微分的联系入手,mathrm{vec}(dF) = frac{partial F}{partial Y}^Tmathrm{vec}(dY) = frac{partial F}{partial Y}^Tfrac{partial Y}{partial X}^Tmathrm{vec}(dX) ,可以推出链式法则frac{partial F}{partial X} = frac{partial Y}{partial X}frac{partial F}{partial Y}

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

  1. (Aotimes B)^T = A^T otimes B^T
  2. mathrm{vec}(oldsymbol{ab}^T) = oldsymbol{b}otimesoldsymbol{a}
  3. (Aotimes B)(Cotimes D) = (AC)otimes (BD)。可以对F = D^TB^TXAC求导来证明,一方面,直接求导得到frac{partial F}{partial X} = (AC) otimes (BD);另一方面,引入Y = B^T X A,有frac{partial F}{partial Y} = C otimes D, frac{partial Y}{partial X} = A otimes B,用链式法则得到frac{partial F}{partial X} = (Aotimes B)(C otimes D)
  4. K_{mn} = K_{nm}^T, K_{mn}K_{nm} = I
  5. K_{pm}Aotimes B K_{nq} = Botimes A,A是m×n矩阵,B是p×q矩阵。可以对AXB^T做向量化来证明,一方面,mathrm{vec}(AXB^T) = (Botimes A)mathrm{vec}(X);另一方面,mathrm{vec}(AXB^T) = K_{pm}mathrm{vec}(BX^TA^T) = (K_{pm}Aotimes B)mathrm{vec}(X^T) = (K_{pm}Aotimes B K_{nq})mathrm{vec}(X)



接下来演示一些算例。

例1:F = AX,X是m×n矩阵,求frac{partial F}{partial X}

解:先求微分:dF=AdX,再做向量化,使用矩阵乘法的技巧,注意在dX右侧添加单位阵:mathrm{vec}(dF) = mathrm{vec}(AdX) = (I_notimes A)mathrm{vec}(dX),对照导数与微分的关系得到frac{partial F}{partial X} = I_notimes A^T

例2:f = log |X| ,X是n×n矩阵,求
abla_X f
abla^2_X f

解:使用上篇中的技术可求得
abla_X f = X^{-1T} 。为求
abla^2_X f,先求微分:d
abla_X f = -(X^{-1}dXX^{-1})^T,再做向量化,使用转置和矩阵乘法的技巧mathrm{vec}(d
abla_X f)= -K_{nn}mathrm{vec}(X^{-1}dX X^{-1}) = -K_{nn}X^{-1T}otimes X^{-1}mathrm{vec}(dX),对照导数与微分的关系得到
abla^2_X f = -K_{nn}X^{-1T}otimes X^{-1}。注意
abla^2_X f是对称矩阵。在X是对称矩阵时,可简化为
abla^2_X f = -X^{-1}otimes X^{-1}

例3:F = Aexp(XB),A是l×m,X是m×n,B是n×p矩阵,exp()为逐元素函数,求frac{partial F}{partial X}

解:先求微分:dF = A(exp(XB)odot (dXB)),再做向量化,使用矩阵乘法的技巧:mathrm{vec}(dF) = (I_potimes A)mathrm{vec}(exp(XB)odot (dXB)),再用逐元素乘法的技巧:mathrm{vec}(dF) = (I_p otimes A) mathrm{diag}(exp(XB))mathrm{vec}(dXB),再用矩阵乘法的技巧:mathrm{vec}(dF) = (I_potimes A)mathrm{diag}(exp(XB))(B^Totimes I_m)mathrm{vec}(dX),对照导数与微分的关系得到frac{partial F}{partial X} = (Botimes I_m)mathrm{diag}(exp(XB))(I_potimes A^T)



最后做个小结。我们发展了从整体出发的矩阵求导的技术,导数与微分的联系是计算的枢纽,标量对矩阵的导数与微分的联系是df = mathrm{tr}(
abla_X^T f dX),先对f求微分,再使用迹技巧可求得导数;矩阵对矩阵的导数与微分的联系是mathrm{vec}(dF) = frac{partial F}{partial X}^T mathrm{vec}(dX),先对F求微分,再使用向量化的技巧可求得导数。

原文地址:https://www.cnblogs.com/hellcat/p/7216399.html