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

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

实际中，我们不知道一些样本对应的真实分布。比如我们手上有$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)

\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)。根据相关概率公式，我们进一步描述其概率意义：
k-th mode:  p ( h n , z n = k ) = p ( z n = k ) p ( h n ∣ z n = k ) all modes:  p ( h n ) = ∑ z p ( h n , z n = k ) = ∑ z p ( z n = k ) p ( h n ∣ z n = k ) posterior of k-th mode:  p ( z n = k ∣ h n ) = p ( z n = k ) p ( h n ∣ z n = k ) p ( h n ) = p ( z n = k ) p ( h n ∣ z n = k ) ∑ z p ( z n = k ) p ( h n ∣ z n = k )

k-th mode: all modes: posterior of k-th mode: ​p(hn​,zn=k)=p(zn=k)p(hn​∣zn=k)p(hn​)=z∑​p(hn​,zn=k)=z∑​p(zn=k)p(hn​∣zn=k)p(zn=k∣hn​)=p(hn​)p(zn=k)p(hn​∣zn=k)​=∑z​p(zn=k)p(hn​∣zn=k)p(zn=k)p(hn​∣zn=k)​​

我们将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，计算后验后验分布，定义为
γ n k t + 1 ≐ p ( z n = k ∣ h n ) = p ( z n = k ) p ( h n ∣ z n = k ) p ( h n ) = α k t C N ( h n ; μ k t , C k t ) ∑ k = 1 K α k t C N ( h n ; μ k t , C k t ) ,    ∀ n , k (22)

\tag{22} γnkt+1​​≐p(zn=k∣hn​)=p(hn​)p(zn=k)p(hn​∣zn=k)​=∑k=1K​αkt​CN(hn​;μkt​,Ckt​)αkt​CN(hn​;μkt​,Ckt​)​,  ∀n,k​(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)

\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)

\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​}的偏导数：
∂ g ∂ α k t + 1 = 0 ∂ g ∂ μ k t + 1 = 0 ∂ g ∂ C k t + 1 = 0 (26)

\tag{26} ∂αkt+1​∂g​∂μkt+1​∂g​∂Ckt+1​∂g​​=0=0=0​(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$求导
∂ g ∂ α k = ∑ n = 1 N p γ n k t + 1 ⋅ 1 α k − λ = 0 ⇒ α k = ∑ n = 1 N p γ n k t + 1 λ

⇒​∂αk​∂g​=n=1∑Np​​γnkt+1​⋅αk​1​−λ=0αk​=λ∑n=1Np​​γnkt+1​​​

又因为${\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}$求导
∂ g ∂ μ k = ∑ n = 1 N p γ n k t + 1 ⋅ ( − 1 σ k 2 ) 2 ( μ k − h n ) = 0 ⇒ ∑ n = 1 N p γ n k t + 1 ⋅ ( μ k − h n ) = 0 ⇒ μ k ⋅ ∑ n = 1 N p γ n k t + 1 = ∑ n = 1 N p γ n k t + 1 h n ⇒ μ k = ∑ n = 1 N p γ n k t + 1 h n ∑ n = 1 N p γ n k t + 1

⇒⇒⇒​∂μk​∂g​=n=1∑Np​​γnkt+1​⋅(−σk2​1​)2(μk​−hn​)=0n=1∑Np​​γnkt+1​⋅(μk​−hn​)=0μk​⋅n=1∑Np​​γnkt+1​=n=1∑Np​​γnkt+1​hn​μk​=∑n=1Np​​γnkt+1​∑n=1Np​​γnkt+1​hn​​​

最后对${\sigma }_k$求导
∂ g ∂ σ k = ∑ n = 1 N p γ n k t + 1 ( 1 σ k 3 ∥ ( h n − μ k t + 1 ) ∥ 2 2 − 2 N U E N B S 1 σ k ) = 0 ⇒ σ k 2 = 1 2 N U E N B S ⋅ ∑ n = 1 N p γ n k t + 1 ∥ ( h n − μ k t + 1 ) ∥ 2 2 ∑ n = 1 N p γ n k t + 1

⇒​∂σk​∂g​=n=1∑Np​​γnkt+1​(σk3​1​∥∥​(hn​−μkt+1​)∥∥​22​−2NUE​NBS​σk​1​)=0σk2​=2NUE​NBS​1​⋅∑n=1Np​​γnkt+1​∑n=1Np​​γnkt+1​∥∥​(hn​−μkt+1​)∥∥​22​​​

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

首先对${\alpha }_k$求导
∂ g ∂ α k = ∑ n = 1 N p γ n k t + 1 ⋅ 1 α k − λ = 0 ⇒ α k = ∑ n = 1 N p γ n k t + 1 λ

⇒​∂αk​∂g​=n=1∑Np​​γnkt+1​⋅αk​1​−λ=0αk​=λ∑n=1Np​​γnkt+1​​​

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

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

​∂x∗∂(bHx)​=0,  ∂x∗∂(xHb)​=b∂x∗∂(xHAx)​=Ax​

进一步，
∂ g ∂ μ k ∗ = ∂ ∂ μ k ∗ ∑ n = 1 N p − γ n k t + 1 ∥ ( C k t + 1 ) − 1 2 ( h n − μ k t + 1 ) ∥ 2 2 = ∑ n = 1 N p − γ n k t + 1 ∂ ∂ μ k ∗ ( h n − μ k t + 1 ) H ( C k t + 1 ) − 1 ( h n − μ k t + 1 ) = ∑ n = 1 N p γ n k t + 1 ( C k t + 1 ) − 1 ( μ k t + 1 − h n ) = ( C k t + 1 ) − 1 ∑ n = 1 N p γ n k t + 1 ( μ k t + 1 − h n ) = 0 ⇒ μ k = ∑ n = 1 N p γ n k t + 1 h n ∑ n = 1 N p γ n k t + 1

∂μk∗​∂g​⇒μk​​=∂μk∗​∂​n=1∑Np​​−γnkt+1​∥∥∥​(Ckt+1​)−21​(hn​−μkt+1​)∥∥∥​22​=n=1∑Np​​−γnkt+1​∂μk∗​∂​(hn​−μkt+1​)H(Ckt+1​)−1(hn​−μkt+1​)=n=1∑Np​​γnkt+1​(Ckt+1​)−1(μkt+1​−hn​)=(Ckt+1​)−1n=1∑Np​​γnkt+1​(μkt+1​−hn​)=0=∑n=1Np​​γnkt+1​∑n=1Np​​γnkt+1​hn​​​

最后对${\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$重写为

g :      ∑ 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 ∣ ) − λ ( ∑ k = 1 K α k − 1 ) = ∑ n = 1 N p ∑ k = 1 K γ n k t + 1 ( − ( h n − μ k t + 1 ) H ( C k t + 1 ) − 1 ( h n − μ k t + 1 ) + log ⁡ α k t + 1 + log ⁡ ∣ ( ( C k t + 1 ) − 1 ) ∗ ∣ ) − λ ( ∑ k = 1 K α k − 1 )

g:​    n=1∑Np​​k=1∑K​γnkt+1​(−∥∥∥​(Ckt+1​)−21​(hn​−μkt+1​)∥∥∥​22​+logαkt+1​−log∣∣​Ckt+1​∣∣​)−λ(k=1∑K​αk​−1)=n=1∑Np​​k=1∑K​γnkt+1​(−(hn​−μkt+1​)H(Ckt+1​)−1(hn​−μkt+1​)+logαkt+1​+log∣∣​((Ckt+1​)−1)∗∣∣​)−λ(k=1∑K​αk​−1)​

根据公式，若$\mathbf{A}\in {\mathbb{C}}^{N\times N}$
∂ y H A H x ∂ A ∗ = x y H ∂ det ( A ) ∂ A = det ( A ) ⋅ ( A T ) − 1

∂A∗∂yHAHx​∂A∂det(A)​​=xyH=det(A)⋅(AT)−1​

对$({\mathbf{C}}_k^{-1}{)}^{\ast }$进行求导:
∂ g ∂ ( C k − 1 ) ∗ = ∑ n = 1 N p γ n k t + 1 ( − ( h n − μ k t + 1 ) ( h n − μ k t + 1 ) H + C k t + 1 ) = C k t + 1 ∑ n = 1 N p γ n k t + 1 − ∑ n = 1 N p γ n k t + 1 ( h n − μ k t + 1 ) ( h n − μ k t + 1 ) H = 0 ⇒ C k t + 1 = ∑ n = 1 N p γ n k t + 1 ( h n − μ k t + 1 ) ( h n − μ k t + 1 ) H ∑ n = 1 N p γ n k t + 1

∂(Ck−1​)∗∂g​⇒Ckt+1​​=n=1∑Np​​γnkt+1​(−(hn​−μkt+1​)(hn​−μkt+1​)H+Ckt+1​)=Ckt+1​n=1∑Np​​γnkt+1​−n=1∑Np​​γnkt+1​(hn​−μkt+1​)(hn​−μkt+1​)H=0=∑n=1Np​​γnkt+1​∑n=1Np​​γnkt+1​(hn​−μkt+1​)(hn​−μkt+1​)H​​

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

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

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

α k = ∑ n = 1 N p γ n k t + 1 N p μ k = ∑ n = 1 N p γ n k t + 1 h n ∑ n = 1 N p γ n k t + 1 σ k 2 = 1 2 N U E N B S ⋅ ∑ n = 1 N p γ n k t + 1 ∥ ( h n − μ k t + 1 ) ∥ 2 2 ∑ n = 1 N p γ n k t + 1

αk​μk​σk2​​=Np​∑n=1Np​​γnkt+1​​=∑n=1Np​​γnkt+1​∑n=1Np​​γnkt+1​hn​​=2NUE​NBS​1​⋅∑n=1Np​​γnkt+1​∑n=1Np​​γnkt+1​∥∥​(hn​−μkt+1​)∥∥​22​​​

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

α k = ∑ n = 1 N p γ n k t + 1 N p μ k = ∑ n = 1 N p γ n k t + 1 h n ∑ n = 1 N p γ n k t + 1 C k = ∑ n = 1 N p γ n k t + 1 ( h n − μ k t + 1 ) ( h n − μ k t + 1 ) H ∑ n = 1 N p γ n k t + 1

αk​μk​Ck​​=Np​∑n=1Np​​γnkt+1​​=∑n=1Np​​γnkt+1​∑n=1Np​​γnkt+1​hn​​=∑n=1Np​​γnkt+1​∑n=1Np​​γnkt+1​(hn​−μkt+1​)(hn​−μkt+1​)H​​