稀疏贝叶斯谱估计及EM算法求解

稀疏贝叶斯

稀疏贝叶斯学习(sparse bayes learning，SBL)最早被提出是作为一种机器学习算法[1]。但是在这里我们主要用它来做谱估计，作为求解稀疏重构问题的方法[2]。稀疏重构还有个更好听的名字叫压缩感知，但我既不知道他哪里压缩了也不明白他怎么个感知法，也有人说这是两回事，在此咱们不纠结，就叫他稀疏重构了。

稀疏重构问题

对于如下问题：

$$\begin{array}{c}\mathbf{t}=\mathbf{\Phi }\mathbf{w}+\mathbf{ϵ}\end{array}\text{\hspace{1em}(1)}$$

其中$\mathbf{\Phi }\in {\mathbb{C}}^{N\times M}$ 为过完备字典，即$rank(\mathbf{\Phi })=NandM>N$ 。$\mathbf{t}\in {\mathbb{C}}^{N\times 1}$ 为观测信号，$\mathbf{w}\in {\mathbb{C}}^{M\times 1}$为权重向量，$\mathbf{ϵ}\in {\mathbb{C}}^{N\times 1}$为观测噪声，希望求得$\mathbf{w}$满足：

$$\begin{array}{c}\mathbf{w}=\arg \underset{\mathbf{w}}{\min }||\mathbf{t}-\mathbf{\Phi }\mathbf{w}|{|}_2^2+\lambda ||\mathbf{w}|{|}_0\end{array}\text{\hspace{1em}(2)}$$

其中$||\mathbf{w}|{|}_0$ 为0范数，即$\mathbf{w}$中非零元的个数。通俗来讲就是用字典矩阵 $\mathbf{\Phi }$ 中最少的向量对$\mathbf{t}$进行表示，所求的$\mathbf{w}$就是求得的权重系数。该问题属于$NP-hard$问题，无法直接求解，但有许多近似解法，包括LASSO，OMP等，本篇所介绍的SBL也是求解算法之一。

稀疏贝叶斯

式(1)中的噪声 $\mathbf{ϵ}$ 通常被认为服从0均值高斯分布，即$ϵ\sim \mathcal{N}(0,{\sigma }^2\mathbf{I})$，则$\mathbf{t}$的条件概率密度函数可以写作:

$$\begin{array}{c}p(\mathbf{t}|\mathbf{w};{\sigma }^2)=(2\pi {\sigma }^2{)}^{-\frac{N}{2}}\exp (-\frac{1}{2{\sigma }^2}||\mathbf{t}-\mathbf{\Phi }\mathbf{w}|{|}_2^2)\end{array}\text{\hspace{1em}(3)}$$

同时假定$\mathbf{w}$的先验服从高斯分布：

$$\begin{array}{c}p(\mathbf{w};\mathbf{\Gamma })=(2\pi {)}^{-\frac{M}{2}}|\mathbf{\Gamma }{|}^{-\frac{1}{2}}\exp (-\frac{1}{2}{\mathbf{w}}^H{\mathbf{\Gamma }}^{-1}\mathbf{w})\end{array}\text{\hspace{1em}(4)}$$

其中$\mathbf{\Gamma }=diag([{\gamma }_1,{\gamma }_2,...,{\gamma }_M{]}^T)$为对角阵，是权重的方差。这里通过式(3)和式(4)可以求得边际概率：

$$\begin{array}{rl}p(\mathbf{t};\mathbf{\Gamma },{\mathbf{\sigma }}^2)=& \int p(\mathbf{t}|\mathbf{w};{\sigma }^2)p(\mathbf{w};\mathbf{\Gamma })d\mathbf{w}\\ =& \int (2\pi {\sigma }^2{)}^{-\frac{N}{2}}(2\pi {)}^{-\frac{M}{2}}|\mathbf{\Gamma }{|}^{-\frac{1}{2}}\exp (-\frac{1}{2{\sigma }^2}||\mathbf{t}-\mathbf{\Phi }\mathbf{w}|{|}_2^2-\frac{1}{2}{\mathbf{w}}^H{\mathbf{\Gamma }}^{-1}\mathbf{w})d\mathbf{w}\end{array}\text{\hspace{1em}(5)}$$

对于上式，先重点看指数里面部分：

$$\begin{array}{rl}& -\frac{1}{2{\sigma }^2}||\mathbf{t}-\mathbf{\Phi }\mathbf{w}|{|}_2^2-\frac{1}{2}{\mathbf{w}}^H{\mathbf{\Gamma }}^{-1}\mathbf{w}\\ =& -\frac{1}{2}[{\sigma }^{-2}{\mathbf{t}}^H\mathbf{t}-2{\sigma }^{-2}{\mathbf{t}}^H\mathbf{\Phi }\mathbf{w}+{\mathbf{w}}^H({\sigma }^{-2}{\mathbf{\Phi }}^H\mathbf{\Phi }+{\mathbf{\Gamma }}^{-1})\mathbf{w}]\\ =& -\frac{1}{2}[{\sigma }^{-2}{\mathbf{t}}^H\mathbf{t}-2{\sigma }^{-2}{\mathbf{t}}^H\mathbf{\Phi }\mathbf{w}+{\mathbf{w}}^H{{\mathbf{\Sigma }}_{\mathbf{w}}}^{-1}\mathbf{w}]\\ =& -\frac{1}{2}({\sigma }^{-2}{{\mathbf{\Sigma }}_{\mathbf{w}}\mathbf{\Phi }}^H\mathbf{t}-\mathbf{w}{)}^H{{\mathbf{\Sigma }}_{\mathbf{w}}}^{-1}({\sigma }^{-2}{{\mathbf{\Sigma }}_{\mathbf{w}}\mathbf{\Phi }}^H\mathbf{t}-\mathbf{w})-\frac{1}{2}[{\mathbf{t}}^H(\mathbf{I}-\mathbf{\Phi }{\mathbf{\Sigma }}_{\mathbf{w}}{\mathbf{\Phi }}^{\mathbf{H}})\mathbf{t}]\end{array}\text{\hspace{1em}(6)}$$

其中${\mathbf{\Sigma }}_{\mathbf{w}}=({\sigma }^{-2}{\mathbf{\Phi }}^H\mathbf{\Phi }+{\mathbf{\Gamma }}^{-1}{)}^{-1}$，并令$(\mathbf{I}-\mathbf{\Phi }{\mathbf{\Sigma }}_{\mathbf{w}}{\mathbf{\Phi }}^{\mathbf{H}}{)}^{-1}={\mathbf{\Sigma }}_{\mathbf{t}}$，根据矩阵求逆引理可以推导出${\mathbf{\Sigma }}_{\mathbf{t}}={\sigma }^2\mathbf{I}+{\mathbf{\Phi \Gamma \Phi }}^H$。
然后将上述结果带入(5):

$$\begin{array}{rl}& p(\mathbf{t};\mathbf{\Gamma },{\sigma }^2)\\ =& (2\pi {\sigma }^2{)}^{-\frac{N}{2}}(2\pi {)}^{-\frac{M}{2}}|\mathbf{\Gamma }{|}^{-\frac{1}{2}}\exp (-\frac{1}{2}{\mathbf{t}}^H{{\mathbf{\Sigma }}_{\mathbf{t}}}^{-1}\mathbf{t})\\ & \int \exp [-\frac{1}{2}({\sigma }^{-2}{{\mathbf{\Sigma }}_{\mathbf{w}}\mathbf{\Phi }}^H\mathbf{t}-\mathbf{w}{)}^H{{\mathbf{\Sigma }}_{\mathbf{w}}}^{-1}({\sigma }^{-2}{{\mathbf{\Sigma }}_{\mathbf{w}}\mathbf{\Phi }}^H\mathbf{t}-\mathbf{w})]d\mathbf{w}\end{array}\text{\hspace{1em}(7)}$$

可以看出指数上凑出了高斯分布的形式，但还差指数外的系数。要凑成完整的高斯分布，注意到${\sigma }^{-\frac{N}{2}}|{\mathbf{\Sigma }}_{\mathbf{w}}{|}^{\frac{1}{2}}|\mathbf{\Gamma }{|}^{-\frac{1}{2}}=|{\mathbf{\Sigma }}_{\mathbf{t}}{|}^{-\frac{1}{2}}$，可得：

$$\begin{array}{rl}& p(\mathbf{t};\mathbf{\Gamma },{\sigma }^2)\\ =& (2\pi {)}^{-\frac{N}{2}}|{\mathbf{\Sigma }}_{\mathbf{t}}{|}^{-\frac{1}{2}}\exp (-\frac{1}{2}{\mathbf{t}}^H{{\mathbf{\Sigma }}_{\mathbf{t}}}^{-1}\mathbf{t})\\ & \int (2\pi {)}^{-\frac{M}{2}}|{\mathbf{\Sigma }}_{\mathbf{w}}{|}^{-\frac{1}{2}}\exp [-\frac{1}{2}({\sigma }^{-2}{{\mathbf{\Sigma }}_{\mathbf{w}}\mathbf{\Phi }}^H\mathbf{t}-\mathbf{w}{)}^H{{\mathbf{\Sigma }}_{\mathbf{w}}}^{-1}({\sigma }^{-2}{{\mathbf{\Sigma }}_{\mathbf{w}}\mathbf{\Phi }}^H\mathbf{t}-\mathbf{w})]d\mathbf{w}\end{array}\text{\hspace{1em}(8)}$$

积分内积完结果为$1$，所以：

$$\begin{array}{rl}& p(\mathbf{t};\mathbf{\Gamma },{\sigma }^2)=(2\pi {)}^{-\frac{N}{2}}|{\mathbf{\Sigma }}_{\mathbf{t}}{|}^{-\frac{1}{2}}\exp (-\frac{1}{2}{\mathbf{t}}^H{{\mathbf{\Sigma }}_{\mathbf{t}}}^{-1}\mathbf{t})\end{array}\text{\hspace{1em}(9)}$$

另外，积分内积掉的为$\mathbf{w}$的后验分布：

$$\begin{array}{rl}& p(\mathbf{w}|\mathbf{t};\mathbf{\Gamma },{\sigma }^2)\\ =& (2\pi {)}^{-\frac{M}{2}}|{\mathbf{\Sigma }}_{\mathbf{w}}{|}^{-\frac{1}{2}}\exp [-\frac{1}{2}({\sigma }^{-2}{{\mathbf{\Sigma }}_{\mathbf{w}}\mathbf{\Phi }}^H\mathbf{t}-\mathbf{w}{)}^H{{\mathbf{\Sigma }}_{\mathbf{w}}}^{-1}({\sigma }^{-2}{{\mathbf{\Sigma }}_{\mathbf{w}}\mathbf{\Phi }}^H\mathbf{t}-\mathbf{w})]\\ =& \mathcal{N}({\mathbf{\mu }}_{\mathbf{w}},{\mathbf{\Sigma }}_{\mathbf{w}})\end{array}\text{\hspace{1em}(10)}$$

其中${\mathbf{\mu }}_{\mathbf{w}}={\sigma }^{-2}{{\mathbf{\Sigma }}_{\mathbf{w}}\mathbf{\Phi }}^H\mathbf{t}$，${\mathbf{\Sigma }}_{\mathbf{w}}=({\sigma }^{-2}{\mathbf{\Phi }}^H\mathbf{\Phi }+{\mathbf{\Gamma }}^{-1}{)}^{-1}$。我们可以选择$\mathbf{\mu }$作为$\mathbf{w}$的求解结果（均值嘛，合情合理），但是算$\mathbf{\mu }$就要知道${\sigma }^2$和$\mathbf{\Gamma }$，这俩可以通过对$p(\mathbf{t};\mathbf{\Gamma },{\sigma }^2)$进行最大似然估计求得，但是很不幸，直接求解似然函数是求不出来的（你可以试试，我没算出来，当然参考文献里也没算出来），而且这大概率不是个凸函数，最优解也算不出来，算个次优解意思意思得了。下面介绍下用EM算法求解这个问题。

EM算法

EM算法是怎么回事，什么思想，什么原理网上其他帖子已经讲的很好了[3]，这里直接介绍如何求解上述问题。

理论推导

我们的的核心问题还是要对$p(\mathbf{t};\mathbf{\Gamma },{\sigma }^2)$最大似然求解参数，即：

$$\begin{array}{rl}{\sigma }^2,\mathbf{\Gamma }=& \arg \underset{{\sigma }^2,\mathbf{\Gamma }}{\max }\mathcal{L}({\sigma }^2,\mathbf{\Gamma })\\ \mathcal{L}({\sigma }^2,\mathbf{\Gamma })& =\log [p(\mathbf{t};\mathbf{\Gamma },{\sigma }^2)]\end{array}\text{\hspace{1em}(11)}$$

然后引入$\mathbf{w}$，将似然函数写作：

$$\begin{array}{rl}\log [p(\mathbf{t};\mathbf{\Gamma },{\sigma }^2)]& =\log [\int p(\mathbf{t},\mathbf{w};\mathbf{\Gamma },{\sigma }^2)d\mathbf{w}]\\ =& \log [\int \frac{Q(\mathbf{w})}{Q(\mathbf{w})}p(\mathbf{t},\mathbf{w};\mathbf{\Gamma },{\sigma }^2)d\mathbf{w}]\end{array}\text{\hspace{1em}(12)}$$

上式中$Q(\mathbf{w})$是$\mathbf{w}$的一个函数，显然它可以是任意值不为零的函数，此时我们认为它是$\mathbf{w}$的某个分布函数，因此上式的积分可以写作期望形式：

$$\begin{array}{rl}\log [p(\mathbf{t};\mathbf{\Gamma },{\sigma }^2)]=& \log [{\mathbf{E}}_{\mathbf{w}\sim Q(\mathbf{w})}\frac{p(\mathbf{t},\mathbf{w};\mathbf{\Gamma },{\sigma }^2)}{Q(\mathbf{w})}]\end{array}\text{\hspace{1em}(13)}$$

对上式进一步操作需要用到Jensen不等式，Jensen不等式网上有详细讲解[4]，大概就是对于一个随机变量$X$ 和一个函数$f(X)$ ，当$f$的二阶导小于0（上凸）时$f(\mathbf{E}(X))\ge \mathbf{E}f(X)$，等号成立的条件是$X=\mathbf{E}(X)$。则式(13)可写为：

$$\begin{array}{rl}\log [p(\mathbf{t};\mathbf{\Gamma },{\sigma }^2)]=& \log [{\mathbf{E}}_{\mathbf{w}\sim Q(\mathbf{w})}\frac{p(\mathbf{t},\mathbf{w};\mathbf{\Gamma },{\sigma }^2)}{Q(\mathbf{w})}]\\ \ge & {\mathbf{E}}_{\mathbf{w}\sim Q(\mathbf{w})}\log [p(\mathbf{t},\mathbf{w};\mathbf{\Gamma },{\sigma }^2)]-{\mathbf{E}}_{\mathbf{w}\sim Q(\mathbf{w})}\log [Q(\mathbf{w})]\\ =& \mathcal{J}({\sigma }^2,\mathbf{\Gamma })\end{array}\text{\hspace{1em}(14)}$$

这里我们得到了$\mathcal{L}({\sigma }^2,\mathbf{\Gamma })$的一个下界$\mathcal{J}({\sigma }^2,\mathbf{\Gamma })$，理论上可以通过最大化$\mathcal{J}({\sigma }^2,\mathbf{\Gamma })$来最大化$\mathcal{L}({\sigma }^2,\mathbf{\Gamma })$，但是一眼望去$\mathcal{J}({\sigma }^2,\mathbf{\Gamma })$更布嚎算，至少$Q(\mathbf{w})$是啥咱还不知道呢！所以得确定下$Q(\mathbf{w})$。

首先可以确定的是，在固定${\sigma }^2$和$\mathbf{\Gamma }$时，根据Jensen不等式取得等号的条件，$\mathcal{J}({\sigma }^2,\mathbf{\Gamma })$能取到的最大值就是$\mathcal{L}({\sigma }^2,\mathbf{\Gamma })$，取得最大值的条件是期望内的数是个常数（与所求期望的随机变量无关）即$\frac{p(\mathbf{t},\mathbf{w};\mathbf{\Gamma },{\sigma }^2)}{Q(\mathbf{w})}$与$\mathbf{w}$无关。不难看出，满足此条件的$Q(\mathbf{w})$为：

$$\begin{array}{r}Q(\mathbf{w})=p(\mathbf{w}|\mathbf{t};\mathbf{\Gamma },{\sigma }^2)\end{array}\text{\hspace{1em}(15)}$$

但是这样直接把$Q(\mathbf{w})$代回$\mathcal{J}({\sigma }^2,\mathbf{\Gamma })$让他等于$\mathcal{L}({\sigma }^2,\mathbf{\Gamma })$那不是就又绕回去了，所以不能这么干。聪明的人们想到了两步走的方法：
第一步，用当前现有的${\sigma }_{(k)}^2$和${\mathbf{\Gamma }}_{(k)}$，得到$Q_{(k)}(\mathbf{w})=p(\mathbf{w}|\mathbf{t};{\mathbf{\Gamma }}_{(k)},{\sigma }_{(k)}^2)$，并用$Q_{(k)}(\mathbf{w})$计算出${\mathcal{J}}_{(k)}({\sigma }^2,\mathbf{\Gamma })={\mathbf{E}}_{\mathbf{w}\sim Q_{(k)}(\mathbf{w})}\log [p(\mathbf{t},\mathbf{w};\mathbf{\Gamma },{\sigma }^2)]-{\mathbf{E}}_{\mathbf{w}\sim Q_{(k)}(\mathbf{w})}\log [Q_{(k)}(\mathbf{w})]$ ，这就是EM算法中的E步，这一步是从Jensen不等式的层面对$\mathcal{J}({\sigma }^2,\mathbf{\Gamma })$进行最大化，使${\mathcal{J}}_{(k)}({\sigma }^2,\mathbf{\Gamma })$逼近$\mathcal{L}({\sigma }^2,\mathbf{\Gamma })$。
第二步，固定$Q_{(k)}(\mathbf{w})$不动，通过最大化${\mathcal{J}}_{(k)}({\sigma }^2,\mathbf{\Gamma })$求${\sigma }_{(k+1)}^2$和${\mathbf{\Gamma }}_{(k+1)}$。由于固定$Q_{(k)}(\mathbf{w})$以后${\mathcal{J}}_{(k)}({\sigma }^2,\mathbf{\Gamma })$中的${\mathbf{E}}_{\mathbf{w}\sim Q_{(k)}(\mathbf{w})}\log [Q_{(k)}(\mathbf{w})]$始终为常数，所以在第一步其实就用不到，直接丢掉。于是乎有：

$$\begin{array}{rl}{\sigma }_{(k+1)}^2,{\mathbf{\Gamma }}_{(k+1)}=& \arg \underset{{\sigma }^2,\mathbf{\Gamma }}{\max }{\mathbf{E}}_{\mathbf{w}\sim Q_{(k)}(\mathbf{w})}\log [p(\mathbf{t},\mathbf{w};\mathbf{\Gamma },{\sigma }^2)]\end{array}\text{\hspace{1em}(16)}$$

这就是EM算法中的M步，这一步是对逼近$\mathcal{L}({\sigma }^2,\mathbf{\Gamma })$的${\mathcal{J}}_{(k)}({\sigma }^2,\mathbf{\Gamma })$最大化，让得到的${\sigma }_{(k+1)}^2,{\mathbf{\Gamma }}_{(k+1)}$进一步趋近最优解。把M步中算出的${\sigma }_{(k+1)}^2,{\mathbf{\Gamma }}_{(k+1)}$再拿回E步，就完成了一次迭代。

目前为止还是理论层面，我们得回到我们的问题，看看${\mathbf{E}}_{\mathbf{w}\sim Q_{(k)}(\mathbf{w})}\log [p(\mathbf{t},\mathbf{w};\mathbf{\Gamma },{\sigma }^2)]$以及${\sigma }_{(k+1)}^2,{\mathbf{\Gamma }}_{(k+1)}$到底怎么算。

EM算法求解SBL

先看${\mathbf{E}}_{\mathbf{w}\sim Q_{(k)}(\mathbf{w})}\log [p(\mathbf{t},\mathbf{w};\mathbf{\Gamma },{\sigma }^2)]$，其中:

$$p(\mathbf{t},\mathbf{w};\mathbf{\Gamma },{\sigma }^2)=p(\mathbf{t};\mathbf{\Gamma },{\sigma }^2)p(\mathbf{w}|\mathbf{t};\mathbf{\Gamma },{\sigma }^2)=p(\mathbf{t}|\mathbf{w};{\sigma }^2)p(\mathbf{w};\mathbf{\Gamma })\text{\hspace{1em}(17)}$$

$p(\mathbf{t};\mathbf{\Gamma },{\sigma }^2)$在式(7)，$p(\mathbf{w}|\mathbf{t};\mathbf{\Gamma },{\sigma }^2)$在式(10)，$p(\mathbf{t}|\mathbf{w};{\sigma }^2)$在式(3)，$p(\mathbf{w};\mathbf{\Gamma })$在式(4)，随便选一对就能算出$p(\mathbf{t},\mathbf{w};\mathbf{\Gamma },{\sigma }^2)$，我这里就不算了。
对计算的结果取对数后得到：

$$\begin{array}{rl}& \log [p(\mathbf{t},\mathbf{w};\mathbf{\Gamma },{\sigma }^2)]\\ =& -\frac{N}{2}log({\sigma }^2)-\frac{1}{2}log|\Gamma |-\frac{1}{2{\sigma }^2}({\mathbf{t}}^H\mathbf{t}+2{\mathbf{t}}^H\mathbf{\Phi }\mathbf{w}+{\mathbf{w}}^H{\mathbf{\Phi }}^H\mathbf{\Phi }\mathbf{w})-\frac{1}{2}{\mathbf{w}}^H{\mathbf{\Gamma }}^{-1}\mathbf{w}+C\end{array}\text{\hspace{1em}(18)}$$

然后求期望，为了避免式子冗长，后面把没用的常数$C$扔了，并将${\mathbf{E}}_{\mathbf{w}\sim Q_{(k)}(\mathbf{w})}$简写为${\mathbf{E}}_{(k)}$：

$$\begin{array}{rl}& {\mathbf{E}}_{(k)}\log [p(\mathbf{t},\mathbf{w};\mathbf{\Gamma },{\sigma }^2)]\\ =& -\frac{N}{2}log({\sigma }^2)-\frac{1}{2}log|\mathbf{\Gamma }|\\ & -\frac{1}{2{\sigma }^2}({\mathbf{t}}^H\mathbf{t}+2{\mathbf{E}}_{(k)}[{\mathbf{t}}^H\mathbf{\Phi }w]+{\mathbf{E}}_{(k)}[{\mathbf{w}}^H{\mathbf{\Phi }}^H\mathbf{\Phi }\mathbf{w}])-\frac{1}{2}{\mathbf{E}}_{(k)}[{\mathbf{w}}^H{\mathbf{\Gamma }}^{-1}\mathbf{w}]\\ =& {\mathcal{J}}_{(k)}({\sigma }^2,\mathbf{\Gamma })\end{array}\text{\hspace{1em}(19)}$$

由于此刻$\mathbf{w}\sim Q_{(k)}(\mathbf{w}),Q_{(k)}(\mathbf{w})=p(\mathbf{w}|\mathbf{t};{\mathbf{\Gamma }}_{(k)},{\sigma }_{(k)}^2)$，并结合式(10)，可得出${\mathbf{E}}_{(k)}[{\mathbf{t}}^H\mathbf{\Phi }w]={\mathbf{t}}^H\mathbf{\Phi }{\mu }_{(k)}$，${\mathbf{E}}_{(k)}[{\mathbf{w}}^H{\mathbf{\Phi }}^H\mathbf{\Phi }\mathbf{w}]=tr[{\mathbf{\Phi }}^H{\mathbf{\Phi }{\mathbf{\Sigma }}_{\mathbf{w}}}^{(k)}]+{\mathbf{\mu }}_{(k)}^H{\mathbf{\Phi }}^H\mathbf{\Phi }{\mathbf{\mu }}_{(k)}$，${\mathbf{E}}_{(k)}[{\mathbf{w}}^H{\mathbf{\Gamma }}^{-1}\mathbf{w}]=tr[{\mathbf{\Gamma }}_{(k)}^{-1}{\mathbf{\Sigma }}_{\mathbf{w}}^{(\mathbf{k})}]+{\mathbf{\mu }}_{(k)}^H{\mathbf{\Gamma }}^{-1}{\mathbf{\mu }}_{(k)}$，其中${{\mathbf{\Sigma }}_{\mathbf{w}}}^{(k)}=({\sigma }_{(k)}^{-2}{\mathbf{\Phi }}^H\mathbf{\Phi }+{\mathbf{\Gamma }}_{(k)}^{-1}{)}^{-1}$，${\mathbf{\mu }}_{(k)}={\sigma }_{(k)}^{-2}{\mathbf{\Sigma }}_{\mathbf{w}}^{(\mathbf{k})}{\mathbf{\Phi }}^H\mathbf{t}$，$tr(\cdot )$为求矩阵的迹，即主对角线元素的和。对于二次型求期望可以参考matrix cookbook[5]。然后有：

$$\begin{array}{rl}& {\mathcal{J}}_{(k)}({\sigma }^2,\mathbf{\Gamma })\\ =& -\frac{N}{2}log({\sigma }^2)-\frac{1}{2}log|\mathbf{\Gamma }|\\ & -\frac{1}{2{\sigma }^2}({\mathbf{t}}^H\mathbf{t}+{\mathbf{t}}^H\mathbf{\Phi }{\mu }_{(k)}+tr[{\mathbf{\Phi }}^H{\mathbf{\Phi }{\mathbf{\Sigma }}_{\mathbf{w}}}^{(k)}]+{\mathbf{\mu }}_{(k)}^H{\mathbf{\Phi }}^H\mathbf{\Phi }{\mathbf{\mu }}_{(k)})\\ & -\frac{1}{2}(tr[{\mathbf{\Gamma }}_{(k)}^{-1}{{\mathbf{\Sigma }}_{\mathbf{w}}}^{(k)}]+{\mathbf{\mu }}_{(k)}^H{\mathbf{\Gamma }}^{-1}{\mathbf{\mu }}_{(k)})\end{array}\text{\hspace{1em}(20)}$$

上述部分算是EM算法中的E步，下面看M步，即最大化${\mathcal{J}}_{(k)}({\sigma }^2,\mathbf{\Gamma })$来求${\sigma }_{(k+1)}^2$和${\mathbf{\Gamma }}_{(k+1)}$。最大化的方式就是求导，导函数为0的点即为最大值。当然，严谨的来讲导函数为0的点是否为最大值还需要经过一些验证，这里就不验证了。求导过程比较简单，也不详细讲了，其中涉及到的矩阵求导可以参考matrix cookbook，这里直接给出结果：

$$\begin{array}{r}\frac{\partial {\mathcal{J}}_{(k)}({\sigma }^2,\mathbf{\Gamma })}{\partial {\gamma }_i}=-\frac{1}{2{\gamma }_i}+\frac{1}{2{\gamma }_i^2}{{\mathbf{\Sigma }}_{\mathbf{w}}}^{(k)}(i,i)+\frac{1}{2{\gamma }_i^2}[{\mathbf{\mu }}_{(k)}(i){]}^2\end{array}\text{\hspace{1em}(21)}$$

其中${\gamma }_i$表示对角阵$\mathbf{\Gamma }$的第$i$个元素，${{\mathbf{\Sigma }}_{\mathbf{w}}}^{(k)}(i,i)$为${{\mathbf{\Sigma }}_{\mathbf{w}}}^{(k)}$的主对角线上第$i$个元素，${\mathbf{\mu }}_{(k)}(i)$为${\mathbf{\mu }}_{(k)}$的第$i$个元素。令式(21)为0可以得到：

$${\gamma }_i^{(k+1)}={{\mathbf{\Sigma }}_{\mathbf{w}}}^{(k)}(i,i)+[{\mathbf{\mu }}_{(k)}(i){]}^2\text{\hspace{1em}(22)}$$

其中${\gamma }_i^{(k+1)}$为${\mathbf{\Gamma }}_{(k+1)}$的第$i$元素。${\sigma }_{(k+1)}^2$同理：

$$\begin{array}{rl}& \frac{\partial {\mathcal{J}}_{(k)}({\sigma }^2,\mathbf{\Gamma })}{\partial {\sigma }^2}=-\frac{N}{2{\sigma }^2}+\frac{1}{2({\sigma }^2{)}^2}(||t-\Phi {\mathbf{\mu }}_{(k)}|{|}_2^2+tr[{\mathbf{\Phi \Phi }}^h{{\mathbf{\Sigma }}_{\mathbf{w}}}^{(k)}])=0\\ & {\sigma }_{(k+1)}^2=\frac{||t-\Phi {\mathbf{\mu }}_{(k)}|{|}_2^2+tr[{\mathbf{\Phi \Phi }}^h{{\mathbf{\Sigma }}_{\mathbf{w}}}^{(k)}]}{N}\end{array}\text{\hspace{1em}(23)}$$

其中的$tr[{\mathbf{\Phi \Phi }}^h{{\mathbf{\Sigma }}_{\mathbf{w}}}^{(k)}]$还可以写作${\sigma }_{(k)}^{2}tr[\mathbf{I}-{\mathbf{\Gamma }}_{(k)}^{-1}{{\mathbf{\Sigma }}_{\mathbf{w}}}^{(k)}]$，读者自证不难）。
至此我们已经完成了全部推导。

回顾式(17)到式(23)不难发现，实际进行EM算法求解的过程中并不需要真正的去求期望，和最大化损失函数，这些步骤已经蕴含在推导中了，而实际要做的只是根据(10)利用${\sigma }_{(k)}^2$和${\mathbf{\Gamma }}_{(k)}$算出${\mathbf{\mu }}_{(k)}$和${{\mathbf{\Sigma }}_{\mathbf{w}}}^{(k)}$，再根据(22)和(23)算出${\sigma }_{(k+1)}^2$和${\mathbf{\Gamma }}_{(k+1)}$就行了。鉴于这仍然是两步走，所以第一步是实际中的E步，第二步为M步。当迭代一步新算出来的结果和上一步相差很小的时候，就可以认为结果收敛了。

代码及结果

matlab 代码如下

function [w,ii] = sbl(t,Phi,gpu,sigma,Gamma)
% t is the received signal
% Phi is the dictionary
% gpu means if you want to accelerate your code by gpu
% sigma is the initial value of sigma, usually set 1
% Gamma is the initial value of Gamma, usually set eye(M)

    % initial
    [N,M] = size(Phi);
    if(strcmp(gpu,'gpu'))
        t = gpuArray(t);
        Phi = gpuArray(Phi);
        Gamma = gpuArray(Gamma);
        matEyeN = gpuArray(eye(N));
        vecOne = gpuArray(ones(M,1));
    else
        matEyeN = eye(N);
        vecOne = ones(M,1);
    end
    sigmaS = sigma; 
    sigmaP = sigmaS;
    
    GammaP = Gamma;

    ii = 0;
    iterErrNow = 100;

    while((iterErrNow>2e-3)&&(ii<500))
        ii = ii+1;
        % e-step
        Sigmat = sigmaS.*matEyeN + Phi*Gamma*Phi';
        Sigmaw = Gamma-Gamma*Phi'*Sigmat^(-1)*Phi*Gamma;
        mu = sigmaS^(-1).*Sigmaw*(Phi')*t;

        % m-step
        Gamma = diag(abs(diag(Sigmaw)+abs(mu).^2));
        sigmaS = abs((t-Phi*mu)'*(t-Phi*mu)+sigmaS*sum(vecOne-diag(Gamma).^(-1).*diag(Sigmaw)))/N;

        iterErrNow = norm(sigmaP-sigmaS)+norm(diag(GammaP)-diag(Gamma)); % stop when variable's change become small enough

        sigmaP = sigmaS;
        GammaP = Gamma;
        w = mu;
    end
    
end

实验代码

close all

%% parameter
N = 64;
M = 16*N;
fs = 1;
f = 0:1/M:(M-1)/M;
t = 0:N-1;
window = rectwin(N);
Phi = exp(2*pi*1i*t'*f);



%% signal module
% fx = 0.3*rand(1,3);
fx = [0.1,0.2,0.21];
a = 0.3+0.7*rand(3,1);
x = exp(2*pi*1i*t'*fx)*a;
x = x./sqrt(var(x));
SNR = 5;
noise = 2^(-0.5)*(randn(N,1)+1i*randn(N,1));
x = x + 10^(-SNR/20)*noise;
% x = x.*window;
xf = fft(x,M)/N;

fig = figure;
tile = tiledlayout(1,1,"TileSpacing","tight");
nexttile
hold on
plot(0:1/M:(M-1)/M,20*log10(abs(xf)))


%% sbl
tic
[wEst,iterNum] = sbl(x,Phi,'gpu',5,eye(M));
toc

plot(0:1/M:(M-1)/M,20*log10(abs(wEst)))
scatter(fx,20*log10(a))
xlim([0,0.5])
ylim([-60,10])
xlabel(tile,'rad')
ylabel(tile,'dB')
legend('fft','sbl','truth')

figNormalize(fig,'save','sblresult');

结果如下：

参考文献

[1]TIPPING M E. Sparse Bayesian Learning and the Relevance Vector Machine[J]. Journal of Machine Learning Research, 2001, 1(Jun): 211-244
[2]WIPF D P, RAO B D. Sparse Bayesian learning for basis selection[J/OL]. IEEE Transactions on Signal Processing, 2004, 52(8): 2153-2164. DOI:10.1109/TSP.2004.831016
[3]https://zhuanlan.zhihu.com/p/78311644
[4]https://zhuanlan.zhihu.com/p/39315786
[5]https://www.math.uwaterloo.ca/~hwolkowi/matrixcookbook.pdf

原创作者: IrregularRhythm 转载于: https://www.cnblogs.com/IrregularRhythm/p/18882385