• 问答
  • 技术
  • 实践
  • 资源
深入浅出多维高斯分布的最大似然估计矩阵推导
技术讨论

作者丨善财童子@知乎(已授权)
来源丨https://zhuanlan.zhihu.com/p/301741267
编辑丨极市平台

之前搞特征分析时,对一个样本矩阵 $X \in R^{m \times n}$ 求协方差矩阵,在 $X$ 中行表示样本数,列表示特征数。直接套协方差的公式: $\Sigma = X^T X$ 。对于这个公式什么来的,我当时没有具体去推导过。今天就从矩阵的角度推导一下(好像这样说不太严谨,欢迎批评指正)。

首先给出多元高斯分布的公式:

$$
G(x_t|\mu, \Sigma) = \frac{1}{(2\pi)^{d/2} det(\Sigma)^{1/2}}exp\left( -\frac{1}{2}(x_t - \mu)^T \Sigma^{-1} (x_t - \mu) \right) \tag{1}
$$

其中 $x_t \in R^n$ 。

定义关于参数 $\theta = { \mu, \Sigma }$ 的似然函数为:

$$
p(X|\theta) = \prod_{t=1}^{m}G(x_t|\theta) \tag{2}
$$

为了方便计算,我们对公式(2)两边取自然对数,得到:

$$
L = ln(p(X|\theta)) = \sum_{t=1}^{m}ln(G(x_t|\theta)) \tag{3}
$$

将公式(1)代入公式(3),并利用对数的性质,将乘法展开成加法,得到:

$$
L = -\frac{md}{2} ln(2\pi) - \frac{m}{2}ln(|\Sigma|) - \frac{1}{2} \sum_{t=1}^{m}(x_t - \mu)^T \Sigma^{-1}(x_t - \mu) \tag{4}
$$

OK,到此为止,我们的准备工作完成了,下面分别对 $\mu$ 和 $\Sigma$ 求导。

这里贴出三条求导法则:

(1) $\frac{\partial (x^T A x)}{\partial} = (A + A^T)x$ ,若A是对称矩阵,则 $\frac{\partial (x^T A x)}{\partial} = 2Ax$ ;

(2) $\frac{\partial ln(|X|)}{\partial X} = (X^{-1})^T$ ,来自微分形式 $d ln(|X|) = tr(X^{-1}dX)$ ;

(3) $dX^{-1} = -X^{-1}dXX^{-1}$ ,可由 $0=dI=d(XX^{-1})=dXX^{-1}+XdX^{-1}$ 推导得到。

1.对 $\mu$ 求偏导数

这里套用法则(1),可得:

$$
\frac{\partial L}{\partial \mu} = \sum_{t=1}^{m}\Sigma^{-1}(x_t - \mu) \tag{5}
$$

注意,公式(5)套用了链式求导法则,即:

$$
\frac{\partial (- \frac{1}{2} \sum_{t=1}^{m}(x_t - \mu)^T \Sigma^{-1}(x_t - \mu))}{\partial \mu} =\frac{\partial (xt - \mu)}{\partial \mu} \frac{- \frac{1}{2} \sum{t=1}^{m} \partial((x_t - \mu)^T \Sigma^{-1}(x_t - \mu))}{\partial(x_t-\mu)}
$$

$\Rightarrow$ $=-1 * (-\frac{1}{2}) \sum_{t=1}^{m} 2 \Sigma^{-1}(xt - \mu)=\sum{t=1}^{m}\Sigma^{-1}(x_t-\mu)$

注意,按照协方差矩阵的定义,协方差矩阵是一个对称矩阵,它的逆矩阵也是对称的。

令 $\frac{\partial L}{\partial \mu} = 0$ ,可得到:

$$
\mu = \frac{1}{m}\sum_{t=1}^{m}x_t \tag{6}
$$

2.对 $\Sigma$ 求偏导数

这里运用矩阵微分来做,可得:

$$
\begin{aligned} \mathrm{d} L &= tr(\mathrm{d}L) \ &=-\frac{m}{2}tr(\Sigma^{-1}\mathrm{d}\Sigma) - \frac{1}{2} \sum_{t=1}^{m}tr((x_t-\mu)^T\mathrm{d}\Sigma^{-1}(xt-\mu)) \ &=-\frac{m}{2}tr(\Sigma^{-1}\mathrm{d}\Sigma) - \frac{1}{2}\sum{t=1}^{m}tr((x_t-\mu)^T(-\Sigma^{-1}\mathrm{d}\Sigma\Sigma^{-1})(xt-\mu))\ &=-\frac{m}{2}tr(\Sigma^{-1}\mathrm{d}\Sigma) + \frac{1}{2}\sum{t=1}^{m}tr\left( \Sigma^{-1}(x_t-\mu)(x_t-\mu)^T\Sigma^{-1}\mathrm{d}\Sigma \right) \end{aligned} \ \tag{7}
$$

根据矩阵微分与导数的关系: $df = tr\left(\left(\frac{df}{\partial x}\right)^Tdx\right)$ ,可得:

$$
\frac{\partial L}{\partial \Sigma} = -\frac{m}{2}\Sigma^{-1} + \frac{1}{2}\sum_{t=1}^{m}\Sigma^{-1}(x_t-\mu)(x_t -\mu)^T\Sigma^{-1} \tag{8}
$$

注意,公式(8)的推导用到了 $\Sigma^{-1} $ 是对称矩阵的事实。

令 $\frac{\partial L}{\partial\Sigma} = 0$ ,整理可得到:

$$
\Sigma = \frac{1}{m} \sum_{t=1}^{m}(x_t-\mu)(x_t-\mu)^T \tag{9}
$$

OK,对比文章开始的协方差矩阵公式 $\Sigma = X^T X$ ,好像和公式(9)不太一样?其实公式(9)首先对每个样本做了中心化处理,即: $x_t - \mu$ 。我们仔细观察一下公式(9)的求和项 $(x_t-\mu)(x_t-\mu)^T $ 。为方便理解,假如 $\mu$ 刚好是0,那么不就是 $xx_t^T$ 了吗,这个是一个 $n \times n$ 矩阵,每个样本都对应这样一个矩阵,将它们相加起来再除以样本数就是协方差矩阵了。

所以可以看到原文开始的公式 $\Sigma = X^T X$ 只是少了中心化和除以样本数而已。

3.验证

我们在Python中验证上面的结论是否正确

import numpy as np
# 100*10的矩阵,标准正态分布,即均值为0,方差为1,所以下面没有中心化操作
X = np.random.randn(100, 10) 
S = X.T @ X # X^T * X

#按照公式(9)计算协方差矩阵,
s = np.zeros((10, 10))
# 注意,下面的for循环计算完后,我没有除以样本数100
for it in X:
    s += it[:,None] @ it[None,:]

# 验证结果是否一致,即S是否等于s
np.allclose(S,s)

运行结果截图:

可见,输出为True,所以得到验证。

  • 0
  • 0
  • 595
收藏
暂无评论