基于高斯混合模型的分布拟合与参数估计：高维复数域的详细推导

基于高斯混合模型的分布拟合

实际中，我们不知道一些样本对应的真实分布。比如我们手上有$N_p$个样本 { h n } n = 1 N p \{ \boldsymbol h_n \}_{n=1}^{N_p} {hn​}n=1Np​​，这些样本都对应到某个特定的label，比如特定的数字或者类。令真实的分布为$f(\mathbf{h})$，我们可以通过经验产生大量样本来近似该分布。定义$f$的近似为$p(\mathbf{h})$，并且限制分布$p$为Gaussian Mixture, i.e.
$$p\left(\mathbf{h}\right)=\sum _{k=1}^K{\alpha }_k\mathcal{C}\mathcal{N}(\mathbf{h};{\mathbf{\mu }}_k,{\mathbf{C}}_k)\text{\hspace{1em}(20)}$$

其中$\mathbf{h}\in {\mathbb{C}}^{L\times 1}$，${\mathbf{\mu }}_{\mathbf{k}}\in {\mathbb{C}}^{L\times 1}$，${\mathbf{C}}_k\in {\mathbb{C}}^{L\times L}$。假设我们固定第$p$个标识label，然后经验生成$N_p$个样本，记为 { h n } n = 1 N p \{ \boldsymbol h_n \}_{n=1}^{N_p} {hn​}n=1Np​​，我们要最大化似然函数：
max ⁡ { α k , μ k , C k } ln ⁡ ( ∏ n = 1 N p ∑ k = 1 K α k C N ( h n ; μ k , C k ) ) = ∑ n = 1 N p ln ⁡ ( ∑ k = 1 K α k C N ( h n ; μ k , C k ) ) (21) \begin{aligned} \max_{\{ \alpha_k, \boldsymbol{\mu}_k, \boldsymbol{C}_k \}} & \ln \left ( \prod_{n=1}^{N_p} \sum_{k=1}^K \alpha_k \mathcal{CN}(\boldsymbol h_n; \boldsymbol \mu_k, \boldsymbol{C}_k) \right ) \\ = & \sum_{n=1}^{N_p} \ln \left( \sum_{k=1}^K \alpha_k \mathcal{CN}(\boldsymbol h_n; \boldsymbol \mu_k, \boldsymbol{C}_k) \right )\\ \end{aligned} \tag{21} {αk​,μk​,Ck​}max​=​ln⎝⎛​n=1∏Np​​k=1∑K​αk​CN(hn​;μk​,Ck​)⎠⎞​n=1∑Np​​ln(k=1∑K​αk​CN(hn​;μk​,Ck​))​(21)

我们采用EM算法进行相关参数的估计。正式计算之前，先引入一个隐变量 z ∈ { 1 , ⋯   , K } z \in \{1,\cdots,K\} z∈{1,⋯,K}表征GMM中的第几个模(mode)。当$K\to \infty$时，我们认为GMM能够很好地逼近原始分布。那么反过来，等效地来看，引入的隐变量$z$可以解释为：先根据隐变量$z$的先验$p(z=k)$确定第$k$个mode，再基于似然$p({\mathbf{h}}_n|z^n=k)$，由第$k$个mode来产生相应的样本生成(sample)。根据相关概率公式，我们进一步描述其概率意义：
$$\begin{array}{rl}\text{k-th mode: }& p({\mathbf{h}}_n,z^n=k)=p(z^n=k)p({\mathbf{h}}_n|z^n=k)\\ \text{all modes: }& p({\mathbf{h}}_n)=\sum _{z}p({\mathbf{h}}_n,z^n=k)=\sum _{z}p(z^n=k)p({\mathbf{h}}_n|z^n=k)\\ \text{posterior of k-th mode: }& p(z^n=k|{\mathbf{h}}_n)=\frac{p(z^n=k)p({\mathbf{h}}_n|z^n=k)}{p({\mathbf{h}}_n)}=\frac{p(z^n=k)p({\mathbf{h}}_n|z^n=k)}{{\sum }_{z}p(z^n=k)p({\mathbf{h}}_n|z^n=k)}\end{array}$$

我们将EM算法的步骤描述为：对于第$t+1$次迭代
(1) E步：固定第$t$次迭代的参数 { α k t , μ k t , C k t } k = 1 K \{\alpha^t_k, \boldsymbol \mu^t_k, \boldsymbol{C}^t_k \}_{k=1}^K {αkt​,μkt​,Ckt​}k=1K​，对$\forall n,k$，即第$n$个sample，第$k$个mode，计算后验后验分布，定义为
$$\begin{array}{rl}{\gamma }_{nk}^{t+1}& \doteq p(z^n=k|{\mathbf{h}}_n)\\ & =\frac{p(z^n=k)p({\mathbf{h}}_n|z^n=k)}{p({\mathbf{h}}_n)}\\ & =\frac{{\alpha }_k^t\mathcal{C}\mathcal{N}({\mathbf{h}}_n;{\mathbf{\mu }}_k^t,{\mathbf{C}}_k^t)}{{\sum }_{k=1}^K{\alpha }_k^t\mathcal{C}\mathcal{N}({\mathbf{h}}_n;{\mathbf{\mu }}_k^t,{\mathbf{C}}_k^t)},\forall n,k\end{array}\text{\hspace{1em}(22)}$$

经过E步，我们便得到了一组后验分布 { γ n k t + 1 } n , k \{\gamma^{t+1}_{nk}\}_{n,k} {γnkt+1​}n,k​

(2) M步：固定后验分布 { γ n k t + 1 } n , k \{\gamma^{t+1}_{nk}\}_{n,k} {γnkt+1​}n,k​，通过最大化证据下界(Evidence Lower BOund)（即对数边际似然$\log p(\mathbf{h};\alpha ,\mathbf{\mu },\mathbf{C})$的下界），来估计第$t+1$次迭代的参数 { α k t + 1 , μ k t + 1 , C k t + 1 } k = 1 K \{ \alpha^{t+1}_k, \boldsymbol \mu^{t+1}_k, \boldsymbol{C}^{t+1}_k \}_{k=1}^K {αkt+1​,μkt+1​,Ckt+1​}k=1K​

     E L B O ( { γ n k t + 1 } n , k , { h n } n ; { α k t + 1 , μ k t + 1 , C k t + 1 } k = 1 K ) = ∑ n = 1 N p ∑ k = 1 K γ n k t + 1 log ⁡ p ( h n , z n = k ) γ n k t + 1 = ∑ n = 1 N p ∑ k = 1 K γ n k t + 1 log ⁡ p ( z n = k ) p ( h n ∣ z n = k ) γ n k t + 1 = ∑ n = 1 N p ∑ k = 1 K γ n k t + 1 log ⁡ α k t + 1 ⋅ C N ( h n ; μ k t + 1 , C k t + 1 ) γ n k t + 1 = ∑ n = 1 N p ∑ k = 1 K γ n k t + 1 ( − ∥ ( C k t + 1 ) − 1 2 ( h n − μ k t + 1 ) ∥ 2 2 + log ⁡ α k t + 1 − log ⁡ ∣ C k t + 1 ∣ ) + C (23) \begin{aligned} & \ \ \ \ ELBO \left (\{\gamma^{t+1}_{nk}\}_{n,k}, \{ \boldsymbol h_n \}_n ; \{ \alpha^{t+1}_k, \boldsymbol \mu^{t+1}_k, \boldsymbol{C}^{t+1}_k \}_{k=1}^K \right ) \\ & = \sum_{n=1}^{N_p} \sum_{k=1}^K \gamma^{t+1}_{nk} \log \frac{p(\boldsymbol h_n,z^n=k)}{ \gamma^{t+1}_{nk}} \\ & = \sum_{n=1}^{N_p} \sum_{k=1}^K \gamma^{t+1}_{nk} \log \frac{p(z^n=k) p(\boldsymbol h_n|z^n=k)}{ \gamma^{t+1}_{nk}} \\ & = \sum_{n=1}^{N_p} \sum_{k=1}^K \gamma^{t+1}_{nk} \log \frac{\alpha^{t+1}_k \cdot \mathcal{CN}(\boldsymbol h_n; \boldsymbol \mu^{t+1}_k, \boldsymbol{C}^{t+1}_k)}{ \gamma^{t+1}_{nk}} \\ &= \sum_{n=1}^{N_p} \sum_{k=1}^K \gamma^{t+1}_{nk} \left ( - \left \Vert (\boldsymbol{C}^{t+1}_k)^{-\frac{1}{2}} \left ( \boldsymbol h_n - \boldsymbol{\mu}^{t+1}_k \right) \right \Vert^2_2 + \log \alpha^{t+1}_k - \log \left \vert \boldsymbol{{C}}^{t+1}_k \right \vert \right) + C \end{aligned} \tag{23} ​    ELBO({γnkt+1​}n,k​,{hn​}n​;{αkt+1​,μkt+1​,Ckt+1​}k=1K​)=n=1∑Np​​k=1∑K​γnkt+1​logγnkt+1​p(hn​,zn=k)​=n=1∑Np​​k=1∑K​γnkt+1​logγnkt+1​p(zn=k)p(hn​∣zn=k)​=n=1∑Np​​k=1∑K​γnkt+1​logγnkt+1​αkt+1​⋅CN(hn​;μkt+1​,Ckt+1​)​=n=1∑Np​​k=1∑K​γnkt+1​(−∥∥∥​(Ckt+1​)−21​(hn​−μkt+1​)∥∥∥​22​+logαkt+1​−log∣∣​Ckt+1​∣∣​)+C​(23)

其中$C$是与 { α k t + 1 , μ k t + 1 , C k t + 1 } k = 1 K \{ \alpha^{t+1}_k, \boldsymbol \mu^{t+1}_k, \boldsymbol{C}^{t+1}_k \}_{k=1}^K {αkt+1​,μkt+1​,Ckt+1​}k=1K​无关的常数。上述最大化ELBO转化为优化问题：
arg max ⁡ { α k t + 1 , μ k t + 1 , C k t + 1 } ∑ n = 1 N p ∑ k = 1 K γ n k t + 1 ( − ∥ ( C k t + 1 ) − 1 2 ( h n − μ k t + 1 ) ∥ 2 2 + log ⁡ α k t + 1 − log ⁡ ∣ C k t + 1 ∣ ) s.t.  ∑ k = 1 K α k = 1 (24) \begin{aligned} & \argmax_{\{ \alpha^{t+1}_k, \boldsymbol \mu^{t+1}_k, \boldsymbol{C}^{t+1}_k \}} \sum_{n=1}^{N_p} \sum_{k=1}^K \gamma^{t+1}_{nk} \left ( - \left \Vert (\boldsymbol{C}^{t+1}_k)^{-\frac{1}{2}} \left (\boldsymbol h_n- \boldsymbol{\mu}^{t+1}_k \right) \right \Vert^2_2 + \log \alpha^{t+1}_k - \log \left \vert \boldsymbol{{C}}^{t+1}_k \right \vert \right) \\ & \text{s.t. } \sum_{k=1}^K \alpha_k = 1 \end{aligned} \tag{24} ​{αkt+1​,μkt+1​,Ckt+1​}argmax​n=1∑Np​​k=1∑K​γnkt+1​(−∥∥∥​(Ckt+1​)−21​(hn​−μkt+1​)∥∥∥​22​+logαkt+1​−log∣∣​Ckt+1​∣∣​)s.t. k=1∑K​αk​=1​(24)

利用拉格朗日函数：
$$g:\sum _{n=1}^{N_p}\sum _{k=1}^K{\gamma }_{nk}^{t+1}\left(-{\Vert ({\mathbf{C}}_k^{t+1}{)}^{-\frac{1}{2}}\left({\mathbf{h}}_n-{\mathbf{\mu }}_k^{t+1}\right)\Vert }_2^2+\log {\alpha }_k^{t+1}-\log \left|{\mathbf{C}}_k^{t+1}\right|\right)-\lambda (\sum _{k=1}^K{\alpha }_k-1)\text{\hspace{1em}(25)}$$

求关于 { α k t + 1 , μ k t + 1 , C k t + 1 } \{ \alpha^{t+1}_k, \boldsymbol \mu^{t+1}_k, \boldsymbol{C}^{t+1}_k \} {αkt+1​,μkt+1​,Ckt+1​}的偏导数：
$$\begin{array}{rl}\frac{\partial g}{\partial {\alpha }_k^{t+1}}& =0\\ \frac{\partial g}{\partial {\mathbf{\mu }}_k^{t+1}}& =0\\ \frac{\partial g}{\partial {\mathbf{C}}_k^{t+1}}& =0\end{array}\text{\hspace{1em}(26)}$$

考虑${\mathbf{C}}_k={\sigma }_k^2\mathbf{I}$

在实际计算时，为了简化，可以取${\mathbf{C}}_k={\sigma }_k^2\mathbf{I}$，即等高线为圆的高斯，这时式(25)可以化简为
$$g:\sum _{n=1}^{N_p}\sum _{k=1}^K{\gamma }_{nk}^{t+1}\left(-\frac{1}{{\sigma }_k^2}{\Vert \left({\mathbf{h}}_n-{\mathbf{\mu }}_k^{t+1}\right)\Vert }_2^2+\log {\alpha }_k^{t+1}-2N_{UE}{N}_{BS}\log {\sigma }_k\right)-\lambda (\sum _{k=1}^K{\alpha }_k-1)$$

将式(26)代入到上式
首先对${\alpha }_k$求导
$$\begin{array}{rl}& \frac{\partial g}{\partial {\alpha }_k^{}}=\sum _{n=1}^{N_p}{\gamma }_{nk}^{t+1}\cdot \frac{1}{{\alpha }_k^{}}-\lambda =0\\ \Rightarrow & {\alpha }_k=\frac{{\sum }_{n=1}^{N_p}{\gamma }_{nk}^{t+1}}{\lambda }\end{array}$$

又因为${\sum }_{k=1}^K{\alpha }_k=1$，因此$\lambda =N_p$，故
$${\alpha }_k=\frac{{\sum }_{n=1}^{N_p}{\gamma }_{nk}^{t+1}}{N_p}$$

然后对${\mathbf{\mu }}_k\in {\mathbb{C}}^{L\times 1}$求导
$$\begin{array}{rl}& \frac{\partial g}{\partial {\mathbf{\mu }}_k^{}}=\sum _{n=1}^{N_p}{\gamma }_{nk}^{t+1}\cdot (-\frac{1}{{\sigma }_k^2})2\left({\mathbf{\mu }}_k-{\mathbf{h}}_n\right)=0\\ \Rightarrow & \sum _{n=1}^{N_p}{\gamma }_{nk}^{t+1}\cdot \left({\mathbf{\mu }}_k-{\mathbf{h}}_n\right)=0\\ \Rightarrow & {\mathbf{\mu }}_k\cdot \sum _{n=1}^{N_p}{\gamma }_{nk}^{t+1}=\sum _{n=1}^{N_p}{\gamma }_{nk}^{t+1}{\mathbf{h}}_n\\ \Rightarrow & {\mathbf{\mu }}_k=\frac{{\sum }_{n=1}^{N_p}{\gamma }_{nk}^{t+1}{\mathbf{h}}_n}{{\sum }_{n=1}^{N_p}{\gamma }_{nk}^{t+1}}\end{array}$$

最后对${\sigma }_k$求导
$$\begin{array}{rl}& \frac{\partial g}{\partial {\sigma }_k}=\sum _{n=1}^{N_p}{\gamma }_{nk}^{t+1}\left(\frac{1}{{\sigma }_k^3}{\Vert \left({\mathbf{h}}_n-{\mathbf{\mu }}_k^{t+1}\right)\Vert }_2^2-2N_{UE}{N}_{BS}\frac{1}{{\sigma }_k}\right)=0\\ \Rightarrow & {\sigma }_k^2=\frac{1}{2N_{UE}{N}_{BS}}\cdot \frac{{\sum }_{n=1}^{N_p}{\gamma }_{nk}^{t+1}{\Vert \left({\mathbf{h}}_n-{\mathbf{\mu }}_k^{t+1}\right)\Vert }_2^2}{{\sum }_{n=1}^{N_p}{\gamma }_{nk}^{t+1}}\end{array}$$

考虑一般的${\mathbf{C}}_k$

首先对${\alpha }_k$求导
$$\begin{array}{rl}& \frac{\partial g}{\partial {\alpha }_k^{}}=\sum _{n=1}^{N_p}{\gamma }_{nk}^{t+1}\cdot \frac{1}{{\alpha }_k^{}}-\lambda =0\\ \Rightarrow & {\alpha }_k=\frac{{\sum }_{n=1}^{N_p}{\gamma }_{nk}^{t+1}}{\lambda }\end{array}$$

又因为${\sum }_{k=1}^K{\alpha }_k=1$，因此$\lambda =N_p$

然后对${\mathbf{\mu }}_k^{\ast }\in {\mathbb{C}}^{L\times 1}$求导，根据公式(复矩阵求导)：
$$\begin{array}{rl}& \frac{\partial ({\mathbf{b}}^H\mathbf{x})}{\partial {\mathbf{x}}^{\ast }}=0,\frac{\partial ({\mathbf{x}}^H\mathbf{b})}{\partial {\mathbf{x}}^{\ast }}=\mathbf{b}\\ & \frac{\partial ({\mathbf{x}}^H\mathbf{A}\mathbf{x})}{\partial {\mathbf{x}}^{\ast }}=\mathbf{A}\mathbf{x}\end{array}$$

进一步，
$$\begin{array}{rl}\frac{\partial g}{\partial {\mathbf{\mu }}_k^{\ast }}& =\frac{\partial }{\partial {\mathbf{\mu }}_k^{\ast }}\sum _{n=1}^{N_p}-{\gamma }_{nk}^{t+1}{\Vert ({\mathbf{C}}_k^{t+1}{)}^{-\frac{1}{2}}\left({\mathbf{h}}_n-{\mathbf{\mu }}_k^{t+1}\right)\Vert }_2^2\\ & =\sum _{n=1}^{N_p}-{\gamma }_{nk}^{t+1}\frac{\partial }{\partial {\mathbf{\mu }}_k^{\ast }}{\left({\mathbf{h}}_n-{\mathbf{\mu }}_k^{t+1}\right)}^H({\mathbf{C}}_k^{t+1}{)}^{-1}\left({\mathbf{h}}_n-{\mathbf{\mu }}_k^{t+1}\right)\\ & =\sum _{n=1}^{N_p}{\gamma }_{nk}^{t+1}({\mathbf{C}}_k^{t+1}{)}^{-1}\left({\mathbf{\mu }}_k^{t+1}-{\mathbf{h}}_n\right)\\ & =({\mathbf{C}}_k^{t+1}{)}^{-1}\sum _{n=1}^{N_p}{\gamma }_{nk}^{t+1}\left({\mathbf{\mu }}_k^{t+1}-{\mathbf{h}}_n\right)=0\\ \Rightarrow {\mathbf{\mu }}_k& =\frac{{\sum }_{n=1}^{N_p}{\gamma }_{nk}^{t+1}{\mathbf{h}}_n}{{\sum }_{n=1}^{N_p}{\gamma }_{nk}^{t+1}}\end{array}$$

最后对${\mathbf{C}}_k$求导，但是实际中，为了方便处理，我们是对$({\mathbf{C}}_k^{-1}{)}^{\ast }$进行求导，在求导之前，我们回顾一些基本概念。

回顾：如果$\mathbf{A}\in {\mathbb{C}}^{N\times N}$是正定矩阵，那么根据$\text{det}(\mathbf{A})={\prod }_n{\lambda }_n$，我们可以得到$\text{det}(\mathbf{A})=\text{det}({\mathbf{A}}^{\ast })=\text{det}({\mathbf{A}}^T)=\text{det}({\mathbf{A}}^H)$。另外，$\text{det}(\mathbf{A})=\frac{1}{\text{det}({\mathbf{A}}^{-1})}$。因此我们可以把拉格朗日函数$g$重写为

$$\begin{array}{rl}g:& \sum _{n=1}^{N_p}\sum _{k=1}^K{\gamma }_{nk}^{t+1}\left(-{\Vert ({\mathbf{C}}_k^{t+1}{)}^{-\frac{1}{2}}\left({\mathbf{h}}_n-{\mathbf{\mu }}_k^{t+1}\right)\Vert }_2^2+\log {\alpha }_k^{t+1}-\log \left|{\mathbf{C}}_k^{t+1}\right|\right)-\lambda (\sum _{k=1}^K{\alpha }_k-1)\\ & =\sum _{n=1}^{N_p}\sum _{k=1}^K{\gamma }_{nk}^{t+1}\left(-{\left({\mathbf{h}}_n-{\mathbf{\mu }}_k^{t+1}\right)}^H({\mathbf{C}}_k^{t+1}{)}^{-1}\left({\mathbf{h}}_n-{\mathbf{\mu }}_k^{t+1}\right)+\log {\alpha }_k^{t+1}+\log \left|(({\mathbf{C}}_k^{t+1}{)}^{-1}{)}^{\ast }\right|\right)-\lambda (\sum _{k=1}^K{\alpha }_k-1)\end{array}$$

根据公式，若$\mathbf{A}\in {\mathbb{C}}^{N\times N}$
$$\begin{array}{rl}\frac{\partial {\mathbf{y}}^H{\mathbf{A}}^H\mathbf{x}}{\partial {\mathbf{A}}^{\ast }}& =\mathbf{x}{\mathbf{y}}^H\\ \frac{\partial \text{det}(\mathbf{A})}{\partial \mathbf{A}}& =\text{det}(\mathbf{A})\cdot ({\mathbf{A}}^T{)}^{-1}\end{array}$$

对$({\mathbf{C}}_k^{-1}{)}^{\ast }$进行求导:
$$\begin{array}{rl}\frac{\partial g}{\partial ({\mathbf{C}}_k^{-1}{)}^{\ast }}& =\sum _{n=1}^{N_p}{\gamma }_{nk}^{t+1}\left(-\left({\mathbf{h}}_n-{\mathbf{\mu }}_k^{t+1}\right){\left({\mathbf{h}}_n-{\mathbf{\mu }}_k^{t+1}\right)}^H+{\mathbf{C}}_k^{t+1}\right)\\ & ={\mathbf{C}}_k^{t+1}\sum _{n=1}^{N_p}{\gamma }_{nk}^{t+1}-\sum _{n=1}^{N_p}{\gamma }_{nk}^{t+1}\left({\mathbf{h}}_n-{\mathbf{\mu }}_k^{t+1}\right){\left({\mathbf{h}}_n-{\mathbf{\mu }}_k^{t+1}\right)}^H=0\\ \Rightarrow {\mathbf{C}}_k^{t+1}& =\frac{{\sum }_{n=1}^{N_p}{\gamma }_{nk}^{t+1}\left({\mathbf{h}}_n-{\mathbf{\mu }}_k^{t+1}\right){\left({\mathbf{h}}_n-{\mathbf{\mu }}_k^{t+1}\right)}^H}{{\sum }_{n=1}^{N_p}{\gamma }_{nk}^{t+1}}\end{array}$$

这与我们的直觉是一致的。

小结：转化为参数估计问题

如果考虑${\mathbf{C}}_k={\sigma }_k^2\mathbf{I}$

$$\begin{array}{rl}{\alpha }_k& =\frac{{\sum }_{n=1}^{N_p}{\gamma }_{nk}^{t+1}}{N_p}\\ {\mathbf{\mu }}_k& =\frac{{\sum }_{n=1}^{N_p}{\gamma }_{nk}^{t+1}{\mathbf{h}}_n}{{\sum }_{n=1}^{N_p}{\gamma }_{nk}^{t+1}}\\ {\sigma }_k^2& =\frac{1}{2N_{UE}{N}_{BS}}\cdot \frac{{\sum }_{n=1}^{N_p}{\gamma }_{nk}^{t+1}{\Vert \left({\mathbf{h}}_n-{\mathbf{\mu }}_k^{t+1}\right)\Vert }_2^2}{{\sum }_{n=1}^{N_p}{\gamma }_{nk}^{t+1}}\end{array}$$

如果考虑一般的${\mathbf{C}}_k$

$$\begin{array}{rl}{\alpha }_k& =\frac{{\sum }_{n=1}^{N_p}{\gamma }_{nk}^{t+1}}{N_p}\\ {\mathbf{\mu }}_k& =\frac{{\sum }_{n=1}^{N_p}{\gamma }_{nk}^{t+1}{\mathbf{h}}_n}{{\sum }_{n=1}^{N_p}{\gamma }_{nk}^{t+1}}\\ {\mathbf{C}}_k^{}& =\frac{{\sum }_{n=1}^{N_p}{\gamma }_{nk}^{t+1}\left({\mathbf{h}}_n-{\mathbf{\mu }}_k^{t+1}\right){\left({\mathbf{h}}_n-{\mathbf{\mu }}_k^{t+1}\right)}^H}{{\sum }_{n=1}^{N_p}{\gamma }_{nk}^{t+1}}\end{array}$$