2. 中国运载火箭技术研究院, 北京 100076
2. China Academy of Launch Vehicle Technology, Beijing 100076, China
随着宽带雷达/通信技术的兴起,高精度高采样率模数转换器成为制约宽带系统硬件成本和功耗的重要因素。同时,高精度高采样率模数转换器产生的海量数据也给数据传输、存取和处理带来了挑战。在接收机中采用低精度量化技术是该问题的一种有效的解决方法。此时,目标感知/信道估计问题可以抽象为广义线性模型(generalized linear model,GLM)中的参数估计问题[1-2]。
GLM得到观测值的过程可以描述为:待估计参数
求解广义线性模型的方法可以分为两大类:确定性方法和Bayes方法。确定性方法首先获取模型的似然函数,然后使用信号的l1范数||x||1来施加稀疏性约束,最后利用凸优化或贪婪方法进行求解。这类方法能够保证收敛性,但通常需要手动调整正则化参数以获得适当的稀疏度。Bayes方法则引入信号的稀疏先验分布,利用似然函数和先验分布得到后验均值和后验方差。在迭代过程中,Bayes方法可以通过近似计算进行求解,从而降低计算复杂度。此外,Bayes方法还能通过期望最大化算法自动学习先验分布中的参数和噪声方差参数,并提供估计的不确定度。常见的Bayes方法包括稀疏Bayes学习(sparse Bayesian learning,SBL)[3-4]和广义近似消息传递(generalized approximate message passing,GAMP)[2]算法等。SBL算法性能稳定,但每次迭代过程均涉及矩阵求逆运算,计算复杂度相对较高。GAMP算法计算复杂度较低,但对测量矩阵A相对敏感。当测量矩阵A的元素不满足独立次Gauss分布时,该算法可能会发散。有关GAMP算法收敛性的证明可以参考文[5]。为了提高GAMP的稳健性,阻尼衰减机制被引入GAMP算法中以保证其收敛性[5-6]。然而,阻尼衰减过强也会在一定程度上减缓算法的收敛速度、降低算法的运行效率。为此,学者们提出了自适应GAMP (adaptive GAMP,AD-GAMP)[7],该算法可以通过自适应步长机制调整算法步长。当测量矩阵的元素稍微偏离次Gauss分布时,AD-GAMP仍能正常工作。然而,由于AD-GAMP没有从根本上消除或减小测量矩阵的相关性,当测量矩阵的相关性较强时,AD-GAMP算法的性能仍然较差。广义向量近似消息传递(generalized vector approximate message passing,GVAMP)算法[8]以及正交近似消息传递(orthogonal approximate message passing,OAMP)算法(等价于向量近似消息传递(vector approximate message passing,VAMP)算法[9-10])[11]均通过将因子节点图中的变量节点转换为向量节点,基于期望传播(expectation propagation)方法进行问题求解,可以保证算法在更广泛的一类测量矩阵(右旋转不变矩阵)中收敛。其中:OAMP或VAMP用于求解标准线性模型,而GVAMP用于求解GLM。
酉变换近似消息传递(unitary approximate message passing,UAMP)算法[12-14]则通过对测量矩阵进行奇异值分解(singular value decomposition,SVD)来改善测量矩阵特性,减小其相关性,以达到更加稳健的重构性能。UAMP算法只能用于求解标准线性模型,无法应用于广义线性模型的求解。为此,文[15]将UAMP算法进行推广得到广义酉变换近似消息传递(generalized unitary approximate message passing,GUAMP)算法。然而,GUAMP算法在一些场景下会在平衡点附近来回振荡,即使经过多次迭代也无法收敛到最优值。为此,本文拟提出一种改进算法以提高GUAMP算法的稳健性。
1 AD-GUAMP算法本章给出AD-GUAMP算法的设计原理及设计方法。首先简要介绍GUAMP算法,细节可参考文[15]。
GUAMP算法的核心思想是通过对矩阵A进行奇异值分解,得到A=UΣVT的分解形式。其中:U是由左奇异向量按列排放形成的左奇异矩阵;V是由右奇异向量按列排放形成的右奇异矩阵;Σ是奇异值矩阵,它是一个对角矩阵,对角线上的元素由奇异值构成。在此基础上,引入额外的隐变量
| $ \boldsymbol{x} \sim p(\boldsymbol{x})=\prod\limits_{i=1}^{N} p\left(x_{i}\right), $ | (1a) |
| $ \boldsymbol{b}=\boldsymbol{Q} \boldsymbol{x}, \quad \boldsymbol{Q}=\mathit{\boldsymbol{ \boldsymbol{\varSigma}}} \boldsymbol{V}^{\mathrm{T}}, $ | (1b) |
| $ \mathit{\boldsymbol{z}}= \boldsymbol{Ub}, $ | (1c) |
| $ \boldsymbol{y}=\mathcal{Q}(\boldsymbol{z}+\boldsymbol{w}). $ | (1d) |
其中:x是待估计的信号向量,符合概率分布p(x);b是x通过乘以矩阵Q获得的隐变量。Q的计算公式为Q=ΣVT。进一步,通过矩阵U与隐变量b的乘积得到表示为z的中间变量。最终,通过添加噪声w,将中间变量z进行量化操作
如图 1[15]所示,GUAMP算法可以通过构建因子节点图来说明其推断过程。图 1a中:圆圈表示变量节点,可以推断对应变量边缘分布;方框表示因子节点,包含模型中先验与观测信息;mx→δ(x)表示从节点x传递到节点δ的消息,实际上表征的是变量x的伪线性观测分布,由于该分布为Gauss分布,可以使用分布的均值
首先,初始化模块A传入模块B的外信息(extrinsic information),包括外均值bAext与外方差vA, bext,可以看作变量b的先验信息。此时模块B等价的观测模型为:
| $ \boldsymbol{b} \sim p(\boldsymbol{b})=\mathcal{N}\left(\boldsymbol{b}_{\mathrm{A}}^{\mathrm{ext}}, \operatorname{diag}\left(\boldsymbol{v}_{\mathrm{A}, \boldsymbol{b}}^{\mathrm{ext}}\right)\right), $ | (2a) |
| $ \boldsymbol{z}=\boldsymbol{U} \boldsymbol{b}, $ | (2b) |
| $ p(\boldsymbol{y} \mid \boldsymbol{z})=\prod\limits_{i=1}^{M} p\left(y_{i} \mid z_{i}\right). $ | (2c) |
该模型正好可以通过GAMP算法进行求解。此时等价测量矩阵相关性已去除,因此算法鲁棒性较好。模块B运行GAMP算法会得到变量b的伪观测值bBext和伪方差vB, bext。由于变量b的先验分布(模块A传给模块B的关于变量b的消息)是Gauss分布,根据期望传播方法[18-19],模块B传给模块A的外信息与变量b的伪观测值、伪方差等价,即模块B传给模块A的外均值为bBext、外方差为vB, bext。此时对应模块A的等价观测模型为:
| $ \boldsymbol{x} \sim p(\boldsymbol{x})=\prod\limits_{i=1}^{N} p\left(x_{i}\right), $ | (3a) |
| $ \boldsymbol{b}=\boldsymbol{Q} \boldsymbol{x}, \quad \boldsymbol{Q}=\mathit{\boldsymbol{ \boldsymbol{\varSigma}}} \boldsymbol{V}^{\mathrm{T}}, $ | (3b) |
| $ \boldsymbol{b}_{\mathrm{B}}^{\mathrm{ext}}=\boldsymbol{b}+\widetilde{\boldsymbol{w}}, \quad \widetilde{\boldsymbol{w}} \sim \mathcal{N}\left(\mathbf{0}, \operatorname{diag}\left(\boldsymbol{v}_{\mathrm{B}, \boldsymbol{b}}^{\mathrm{ext}}\right)\right). $ | (3c) |
可以看出,模块B传递给模块A的外部信息可以被用作模块A的伪观测值和伪方差。在模块A中运行AMP算法,并输出变量b的等效先验信息,将其作为模块B的输入,通过这个迭代过程来实现闭环。相较于GAMP和GVAMP算法,GUAMP算法在处理相关矩阵时表现出更好的鲁棒性[11]。然而,本文作者通过数值仿真实验发现,GUAMP算法在平衡点附近会出现振荡现象,并且在低稀疏度情况下性能稍逊于GAMP和GVAMP算法。此外,随着测量矩阵相关性的增加,与AD-GAMP相比,GUAMP算法的性能会迅速恶化。这些仿真实验结果表明,有必要对GUAMP算法进行改进,以提高其稳健性。
由于GUAMP算法可以分解成AMP和GAMP两个模块,因此本文针对这两个模块分别构造相应的目标函数,以寻求特定结构的后验分布,并设计相应的步长。针对模块A的模型参考式(3),可以得到变量x和变量b的联合后验分布为
| $ \begin{gather*} p\left(\boldsymbol{x}, \boldsymbol{b} \mid \boldsymbol{b}_{\mathrm{B}}^{\text {ext }}\right)= \\ p(\boldsymbol{x}) \delta(\boldsymbol{b}-\boldsymbol{Q} \boldsymbol{x}) p\left(\boldsymbol{b}_{\mathrm{B}}^{\text {ext }} \mid \boldsymbol{b}\right) \cdot K^{-1}\left(\boldsymbol{b}_{\mathrm{B}}^{\text {ext }}\right). \end{gather*} $ | (4) |
其中K(bBext)是归一化常数,用于保证后验概率密度函数p(x, b|bBext)积分为1。同时,期望得到x和b的近似后验分布为
| $ \widetilde{p}(\boldsymbol{x}, \boldsymbol{b})=\widetilde{p}(\boldsymbol{x}) \delta(\boldsymbol{b}-\boldsymbol{Q x}). $ | (5) |
因此,为了求得
| $ \min\limits_{\widetilde{p}(\boldsymbol{x}, \boldsymbol{b}), \widetilde{p}(\boldsymbol{x})} D\left(\widetilde{p}(\boldsymbol{x}, \boldsymbol{b}) \| p\left(\boldsymbol{x}, \boldsymbol{b} \mid \boldsymbol{b}_{\mathrm{B}}^{\mathrm{ext}}\right)\right), $ | (6a) |
| $ \text { s.t. } \;\;\;\; \widetilde{p}(\boldsymbol{x}, \boldsymbol{b})=\widetilde{p}(\boldsymbol{x}) \delta(\boldsymbol{b}-\boldsymbol{Q} \boldsymbol{x}) . $ | (6b) |
其中D(p1||p2)表示两个分布p1与p2的Kullback-Leibler(KL)散度。待求解的变量是概率密度函数
| $ \begin{gather*} D\left(\widetilde{p}(\boldsymbol{x}, \boldsymbol{b}) \| p\left(\boldsymbol{x}, \boldsymbol{b} \mid \boldsymbol{b}_{\mathrm{B}}^{\text {ext }}\right)\right) \triangleq \int_{\boldsymbol{x}, \boldsymbol{b}} \widetilde{p}(\boldsymbol{x}, \boldsymbol{b}) \ln \frac{\widetilde{p}(\boldsymbol{x}, \boldsymbol{b})}{p\left(\boldsymbol{x}, \boldsymbol{b} \mid \boldsymbol{b}_{\mathrm{B}}^{\text {ext }}\right)} \mathrm{d} \boldsymbol{x} \mathrm{d} \boldsymbol{b}= \\ \int_{\boldsymbol{x}, \boldsymbol{b}} \widetilde{p}(\boldsymbol{x}, \boldsymbol{b}) \ln \frac{\widetilde{p}(\boldsymbol{x})}{p(\boldsymbol{x}) p\left(\boldsymbol{b}_{\mathrm{B}}^{\text {ext }} \mid \boldsymbol{b}\right) K^{-1}\left(\boldsymbol{b}_{\mathrm{B}}^{\text {ext }}\right)} \mathrm{d} \boldsymbol{x} \mathrm{d}\boldsymbol{b}= \\ \int_{\boldsymbol{x}} \widetilde{p}(\boldsymbol{x}) \ln \frac{\widetilde{p}(\boldsymbol{x})}{p(\boldsymbol{x})} \mathrm{d} \boldsymbol{x}+\int_{\boldsymbol{b}} \widetilde{p}(\boldsymbol{b}) \ln \frac{\widetilde{p}(\boldsymbol{b})}{p\left(\boldsymbol{b}_{\mathrm{B}}^{\text {ext }} \mid \boldsymbol{b}\right) C^{-1}\left(\boldsymbol{b}_{\mathrm{B}}^{\text {ext }}\right)} \mathrm{d} \boldsymbol{b}+\int_{\boldsymbol{b}} \widetilde{p}(\boldsymbol{b}) \ln \frac{1}{\widetilde{p}(\boldsymbol{b})} \mathrm{d} \boldsymbol{b}+\ln \frac{C^{-1}\left(\boldsymbol{b}_{\mathrm{B}}^{\text {ext }}\right)}{K^{-1}\left(\boldsymbol{b}_{\mathrm{B}}^{\text {ext }}\right)}= \\ D\left(\widetilde{p}_{\mathit{\boldsymbol{x}}} \| p_{\mathit{\boldsymbol{x}}}\right)+D\left(\widetilde{p}_{\boldsymbol{b}} \| p_{\boldsymbol{b}_{\mathrm{B}}^{\text {ext }} \mid \boldsymbol{b}} C^{-1}\left(\boldsymbol{b}_{\mathrm{B}}^{\text {ext }}\right)\right)+H_{\mathrm{A}}\left(\widetilde{p}_{\boldsymbol{b}}\right)+\ln \frac{C^{-1}\left(\boldsymbol{b}_{\mathrm{B}}^{\text {ext }}\right)}{K^{-1}\left(\boldsymbol{b}_{\mathrm{B}}^{\text {ext }}\right)}. \end{gather*} $ | (7) |
其中:HA(
| $ \min\limits_{\widetilde{p}(\boldsymbol{b}), \widetilde{p}(\boldsymbol{x})} J_{\mathrm{A}}\left(\widetilde{p}_{\mathit{\boldsymbol{x}}}, \widetilde{p}_{\boldsymbol{b}}\right), $ | (8a) |
| $ \text { s.t. }\;\;\widetilde{p}_{\boldsymbol{b}}=T_{Q} \widetilde{p}_{\mathit{\boldsymbol{x}}} . $ | (8b) |
其中
| $ \begin{gather*} J_{\mathrm{A}}\left(\widetilde{p}_{\boldsymbol{x}}, \widetilde{p}_{\boldsymbol{b}}\right) \triangleq D\left(\widetilde{p}_{\boldsymbol{x}} \| \prod\limits_{j=1}^{N} p_{x}\left(\boldsymbol{x}_{j}\right)\right)+ \\ D\left(\widetilde{p}_{\boldsymbol{b}} \| \prod\limits_{i=1}^{M} p_{b_{\mathrm{B}}^{\mathrm{ext}} \mid b}\left(b_{\mathrm{B}, i}^{\mathrm{ext}} \mid b_{i}\right) C^{-1}\left(b_{\mathrm{B}, i}^{\mathrm{ext}}\right)\right)+H_{\mathrm{A}}\left(\widetilde{p}_{\boldsymbol{b}}\right). \end{gather*} $ | (9) |
TQ为变换函数,即随机向量b的每个元素是通过随机向量x的各元素线性加权得到的,即可以通过一个线性变换矩阵Q和随机向量b的分布p(x)来表征分布p(b)。约束项
| $ \widetilde{p}_{\boldsymbol{x}}(\boldsymbol{x}) =\prod\limits_{j=1}^{N} \widetilde{p}_{x}\left(x_{j}\right), $ | (10) |
| $ \widetilde{p}_{\boldsymbol{b}}(\boldsymbol{b}) =\prod\limits_{i=1}^{M} \widetilde{p}_{b}\left(b_{i}\right). $ | (11) |
此外,当给定分布的一阶矩和二阶矩时,熵函数HA(
| $ \min\limits_{\widetilde{p}(\boldsymbol{b}), \widetilde{p}(\boldsymbol{x})} \widetilde{J}_{\mathrm{A}}\left(\widetilde{p}_{\boldsymbol{x}}, \widetilde{p}_{\boldsymbol{b}}\right), $ | (12a) |
| $ \begin{array}{ll} \text { s. t. } & \mathscr{E}\left[\boldsymbol{b} \mid \widetilde{p}_{\boldsymbol{b}}\right]=\boldsymbol{Q} \mathscr{E}\left[\boldsymbol{x} \mid \widetilde{p}_{\boldsymbol{x}}\right], \\ & \boldsymbol{\tau}_{A}^{p}=\boldsymbol{Q}^{2} \operatorname{var}\left[\boldsymbol{x} \mid \widetilde{p}_{\boldsymbol{x}}\right] . \end{array} $ | (12b) |
其中
| $ \begin{array}{*{20}{c}} \widetilde{J}_{\mathrm{A}}\left(\widetilde{p}_{x}, \widetilde{p}_{b}\right) \triangleq \sum\limits_{j=1}^{N} \underbrace{D\left(\tilde{p}_{x}\left(x_{j}\right) \| p_{x}\left(x_{j}\right)\right)}_{\text {第1项 }}+ \\ \sum\limits_{i=1}^{M} \underbrace{D\left(\widetilde{p}_{b}\left(b_{i}\right) \| p_{b_{\mathrm{B}}^{\mathrm{ext}} \mid b}\left(b_{\mathrm{B}, i}^{\mathrm{ext}} \mid b_{i}\right) C_{i}^{-1}\left(b_{\mathrm{B}, i}^{\mathrm{ext}}\right)\right)}_{\text {第 } 2 \text { 项 }}+ \\ \sum\limits_{i=1}^{M} \underbrace{H_{\mathrm{A}, \mathrm{G}}\left(\widetilde{p}_{b}\left(b_{i}\right), \tau_{\mathrm{A}, i}^{p}\right)}_{\text {第3项 }} . \end{array} $ | (13) |
其中:
| $ \begin{gather*} D\left(\tilde{p}_{x}\left(x_{j}\right) \| p_{x}\left(x_{j}\right)\right)= \\ \int \widetilde{p}_{x}\left(x_{j} ; \hat{r}_{\mathrm{A}, j}, \tau_{\mathrm{A}, j}^{r}\right) \cdot \\ \ln \frac{p_{x}\left(x_{j}\right) \mathcal{N}\left(x_{j} ; \hat{r}_{\mathrm{A}, j}, \tau_{\mathrm{A}, j}^{r}\right)}{p_{x}\left(x_{j}\right) X_{j}\left(\hat{r}_{\mathrm{A}, j}, \tau_{\mathrm{A}, j}^{r}\right)} \mathrm{d} x_{j}= \\ -\ln X_{j}\left(\hat{r}_{\mathrm{A}, j}, \tau_{\mathrm{A}, j}^{r}\right)- \\ \frac{\ln 2 {\mathsf{π}} \tau_{\mathrm{A}, j}^{r}}{2}-\frac{\left|\hat{x}_{j}-\hat{r}_{\mathrm{A}, j}\right|^{2}+\tau_{j}^{x}}{2 \tau_{j}^{x}} . \end{gather*} $ | (14) |
其中:
| $ \hat{x}_{j}=\frac{\lambda}{\lambda+(1-\lambda) \exp \left(-M_{j}\right)} \frac{\tau_{p, j}^{x} \hat{r}_{\mathrm{A}, j}+\tau_{\mathrm{A}, j}^{r} x_{p, j}}{\tau_{\mathrm{A}, j}^{r}+\tau_{p, j}^{x}} . $ | (15) |
| $ \begin{gather*} \tau_{j}^{x}=\frac{\lambda}{\lambda+(1-\lambda) \exp \left(-M_{j}\right)} \cdot \\ \left(\left(\frac{\tau_{p, j}^{x} \hat{r}_{\mathrm{A}, j}+\tau_{\mathrm{A}, j}^{r} x_{p, j}}{\tau_{\mathrm{A}, j}^{r}+\tau_{p, j}^{x}}\right)^{2}+\frac{\tau_{p, j}^{x} \tau_{\mathrm{A}, j}^{r}}{\tau_{p, j}^{x}+\tau_{\mathrm{A}, j}^{r}}\right)- \\ \left(\frac{\lambda}{\lambda+(1-\lambda) \exp \left(-M_{j}\right)} \frac{\tau_{p, j}^{x} \hat{r}_{\mathrm{A}, j}+\tau_{\mathrm{A}, j}^{r} x_{p, j}}{\tau_{\mathrm{A}, j}^{r}+\tau_{p, j}^{x}}\right)^{2} . \end{gather*} $ | (16) |
其中
| $ M_{j}=\frac{1}{2} \ln \frac{\tau_{\mathrm{A}, j}^{r}}{\tau_{\mathrm{A}, j}^{r}+\tau_{p, j}^{x}}+\frac{\left(\hat{r}_{\mathrm{A}, j}\right)^{2}}{2 \tau_{\mathrm{A}, j}^{r}}-\frac{\left(\hat{r}_{\mathrm{A}, j}-x_{p, j}\right)^{2}}{2\left(\tau_{\mathrm{A}, j}^{r}+\tau_{p, j}^{x}\right)} $ | (17) |
对于
| $ \begin{gather*} D\left(\widetilde{p}_{b}\left(b_{i}\right) \| p_{b_{\mathrm{B}}^{\text {ext } \mid b}}\left(\widetilde{b}_{i} \mid b_{i}\right)\right)= \\ -\ln C_{i}\left(\hat{p}_{\mathrm{A}, i}, \tau_{\mathrm{A}, i}^{p}\right)-\frac{1}{2} \ln \left(2 {\mathsf{π}} \tau_{\mathrm{A}, i}^{p}\right)- \\ \frac{\left|\hat{b}_{i}-\hat{p}_{\mathrm{A}, i}\right|^{2}+\tau_{i}^{b}}{2 \tau_{\mathrm{A}, i}^{p}} \end{gather*} $ | (18) |
其中:
| $ \frac{1}{\tau_{i}^{b}}=\frac{1}{\tau_{\mathrm{B}, i}^{r}}+\frac{1}{\tau_{\mathrm{A}, i}^{p}}, $ | (19) |
| $ \frac{\hat{b}_{i}}{\tau_{i}^{b}}=\frac{\hat{r}_{\mathrm{B}, i}}{\tau_{\mathrm{B}, i}^{r}}+\frac{\hat{p}_{\mathrm{A}, i}}{\tau_{\mathrm{A}, i}^{p}} . $ | (20) |
| $ \begin{gather*} H_{\mathrm{A}, \mathrm{G}}\left(\widetilde{p}_{b}\left(b_{i}\right), \tau_{\mathrm{A}, i}^{p}\right) \triangleq \\ \frac{1}{2} \frac{\operatorname{var}\left[b_{i} \mid \widetilde{p}_{b}\left(b_{i}\right)\right]}{\tau_{\mathrm{A}, i}^{p}}+\ln \left(2 {\mathsf{π}} \tau_{\mathrm{A}, i}^{p}\right). \end{gather*} $ | (21) |
其中var[bi|
模块B构造的目标函数与上面的推导过程类似。对于模块B,变量b和变量z的联合后验分布为
| $ \begin{gather*} p\left(\boldsymbol{b}, \boldsymbol{z} \mid \boldsymbol{b}_{\mathrm{A}}^{\text {ext }}, \boldsymbol{y}\right)= \\ p\left(\boldsymbol{b} \mid \boldsymbol{b}_{\mathrm{A}}^{\text {ext }}\right) \delta(\boldsymbol{z}-\boldsymbol{U} \boldsymbol{b}) p(\boldsymbol{y} \mid \boldsymbol{z}) \cdot S^{-1}(\boldsymbol{y}) . \end{gather*} $ | (22) |
其中:S(y)为归一化常数使得式(22)关于变量b、z积分为1。期望得到b和z的后验分布为
| $ \widetilde{p}(\boldsymbol{b}, \boldsymbol{z})=\tilde{p}(\boldsymbol{b}) \delta(\boldsymbol{z}-\boldsymbol{U} \boldsymbol{b}) . $ | (23) |
因此,模块B也可以构造如下的优化问题:
| $ \min\limits_{\widetilde{p}(\boldsymbol{b}, \boldsymbol{z}), \widetilde{p}(\boldsymbol{b})} D\left(\widetilde{p}(\boldsymbol{b}, \boldsymbol{z}) \| p\left(\boldsymbol{b}, \boldsymbol{z} \mid \boldsymbol{y}, \boldsymbol{b}_{\mathrm{A}}^{\mathrm{ext}}\right)\right), $ | (24a) |
| $ \text { s. t. } \quad \widetilde{p}(\boldsymbol{b}, \boldsymbol{z})=\widetilde{p}(\boldsymbol{b}) \delta(\boldsymbol{z}-\boldsymbol{U} \boldsymbol{b}). $ | (24b) |
其中bAext为模块A传给模块B的外信息,包含变量b的先验信息。经过简化后,可以得到与模块B目标函数为:
| $ \min\limits_{\widetilde{p}(\boldsymbol{b}), \widetilde{p}(z)} J_{\mathrm{B}}\left(\widetilde{p}_{\boldsymbol{b}}, \widetilde{p}_{z}\right), $ | (25a) |
| $ \begin{array}{ll} \text { s.t. } & \mathscr{E}\left[\boldsymbol{z} \mid \widetilde{p}_{z}\right]=\boldsymbol{U} \mathscr{E}\left[\boldsymbol{b} \mid \widetilde{p}_{\boldsymbol{b}}\right], \\ & \boldsymbol{\tau}_{\mathrm{B}}^{p}=\boldsymbol{U}^{2} \operatorname{var}\left[\boldsymbol{b} \mid \widetilde{p}_{\boldsymbol{b}}\right] . \end{array} $ | (25b) |
其中
| $ \begin{array}{*{20}{c}} & J_{\mathrm{B}}\left(\widetilde{p}_{b}, \widetilde{p}_{z}\right) \triangleq \sum\limits_{i=1}^{M} \underbrace{D\left(\widetilde{p}_{b}\left(b_{i}\right) \| p_{b \mid b_{\mathrm{A}}^{\mathrm{ext}}}\left(b_{i}\right)\right)}_{\text {第1 项 }}+ \\ & \sum\limits_{i=1}^{M} \underbrace{D\left(\widetilde{p}_{z}\left(z_{i}\right) \| p_{y \mid z}\left(y_{i} \mid z_{i}\right) Z_{i}^{-1}\left(y_{i}\right)\right)}_{\text {第 } 2 \text { 项 }}+ \\ & \sum\limits_{i=1}^{M} \underbrace{H_{\mathrm{B}, \mathrm{G}}\left(\widetilde{p}_{z}\left(z_{i}\right), \tau_{\mathrm{B}, i}^{p}\right)}_{\text {第3项 }} . \end{array} $ | (26) |
Zi(yi)为归一化常数。与模块A的计算方法类似,JB(
| $ \begin{gather*} D\left(\widetilde{p}_{b}\left(b_{i}\right) \| p_{b \mid b_{\mathrm{A}}^{\mathrm{ext}}}\left(b_{i}\right)\right)=-\ln B_{i}\left(\hat{r}_{\mathrm{B}, i}, \tau_{\mathrm{B}, i}^{r}\right)- \\ \frac{\ln \left(2 {\mathsf{π}} \tau_{\mathrm{B}, i}^{r}\right)}{2}-\frac{\left|\hat{b}_{i}-\hat{r}_{\mathrm{B}, i}\right|^{2}+\tau_{i}^{b}}{2 \tau_{\mathrm{B}, i}^{r}} . \end{gather*} $ | (27) |
其中:
| $ \begin{gather*} D\left(\widetilde{p}_{z}\left(z_{i}\right) \| p_{y \mid z}\left(y_{i} \mid z_{i}\right) \cdot Z_{i}^{-1}\right)= \\ -\ln E_{i}\left(\hat{p}_{\mathrm{B}, i}, \tau_{\mathrm{B}, i}^{p}\right)-\frac{1}{2} \ln \left(2 {\mathsf{π}} \tau_{\mathrm{B}, i}^{p}\right)- \\ \frac{\left|\hat{z}_{i}-\hat{p}_{\mathrm{B}, i}\right|^{2}+\tau_{i}^{z}}{2 \tau_{\mathrm{B}, i}^{p}} . \end{gather*} $ | (28) |
其中:
| $ \begin{gather*} H_{\mathrm{B}, \mathrm{G}}\left(\widetilde{p}_{z}\left(z_{i}\right), \tau_{\mathrm{B}, i}^{p}\right) \triangleq \\ \frac{1}{2} \frac{\operatorname{var}\left[z_{i} \mid \widetilde{p}_{z}\right]}{\tau_{\mathrm{B}, i}^{p}}+\ln \left(2 {\mathsf{π}} \tau_{\mathrm{B}, i}^{p}\right) . \end{gather*} $ | (29) |
至此,两个模块的目标函数构造完成。下面介绍如何利用目标函数进行算法步长选择。由于两个模块的步长选择采用相同的策略,仅目标函数与调整的具体变量有所不同。因此,在下面有关步长选择策略的描述中,与步长选择有关的变量不再区分两个模块,用J(t)来统一描述,而不单独区分模块A目标函数JA(t)或模块B目标函数JB(t)。设时刻t目标函数取值为J(t),则时刻t处迭代步长β(t) 的选择规则如下:1) 根据窗长度Win,选择最近的目标函数值共Win个,并求取其中最小值。2) 若当前时刻目标函数J(t)大于上面选取序列的最小值,则认为算法仍在正常迭代,可以增加步长以加快算法收敛速度,步长增量由上一次步长β(t-1)与一个大于1的常数Inc确定,即β(t)=Inc· β(t-1);若不满足上述条件,则认为算法振荡或即将收敛,需要减小步长,步长减少量按照上一次步长β(t-1)与一个小于1的常数Dec确定,即β(t)=Dec·β(t-1)。3) 增大或减小步长均不会超过预设的步长最大值βmax或步长最小值βmin,直至算法收敛到一定精度或者迭代一定次数后结束迭代。
GUAMP算法的详细步骤请参考文[15]。将上述自适应步长选择规则耦合到GUAMP算法中,即可得到AD-GUAMP算法,如图 2所示。图中:gx(
| $ g_{x}\left(\hat{r}, \tau^{r}\right) =\int x \frac{p(x) \mathcal{N}\left(x ; \hat{r}, \tau^{r}\right)}{X\left(\hat{r}, \tau^{r}\right)} \mathrm{d} x, $ | (30) |
| $ g_{z}\left(\hat{p}, \tau^{p}\right)=\int z \frac{p(y \mid z) \mathcal{N}\left(z ; \hat{p}, \tau^{p}\right)}{D\left(\hat{p}, \tau^{p}\right)} \mathrm{d} z . $ | (31) |
|
| 图 2 AD-GUAMP算法流程图 |
其中:X(
本章将通过数值仿真实验来验证所提AD-GUAMP算法的有效性。下面共考虑3种测量矩阵的场景,分别为低秩矩阵、条件数病态矩阵和相关矩阵。
1) 低秩矩阵。
本文使用矩阵秩与其最小维度的比值r,即秩比,来描述矩阵的低秩程度。该值越大表明矩阵越接近满秩,矩阵各行或各列相关性越弱,算法重构性能越好。
矩阵构造方式如下:A=UV。其中:U∈
2) 条件数病态矩阵。
本文使用条件数κ来刻画矩阵的病态特性。条件数越大,表明算法重构问题越困难,信号估计性能越差。矩阵构造方式如下:A=UΣVT。其中:U、V为Gauss矩阵(矩阵各元素从Gauss分布中抽样得到)经奇异值分解产生,U为左奇异矩阵,V为右奇异矩阵。对角矩阵Σ对角线上的元素满足Σi,i/Σi+1,i+1=κ1/min{M,N}。
3) 相关矩阵。
本文使用相关系数ρ来描述一类相关矩阵的相关性。相关系数越大,矩阵各列相关性越强,算法的重构性能越差。相关矩阵构造方式如下[20]:A=RLHRR。其中:RL=R11/2∈
在低秩测量矩阵和条件数病态测量矩阵情况下,理论上可以证明GVAMP算法最优[8]。因此,在这两种情况下,本文所设计的AD-GUAMP算法能够逼近GVAMP算法性能即可表明AD-GUAMP算法的有效性。
仿真条件设置如下:信号x先验服从Bernoulli-Gauss分布,即p(x|q)=
AD-GUAMP算法为两个消息传递模块分别设计了目标函数。下面以相关测量矩阵为例说明目标函数的收敛情况。仿真条件设置如下:信噪比SNR=40 dB,2 bit量化,过采样因子δ=0.8,测量矩阵的相关系数ρ=0.1,其他仿真参数与之前设置的仿真参数保持一致。图 2绘制的函数曲线为-JA(t)与-JB(t),可见目标函数值会随着迭代次数的增加而不断上升。从图 3可以看出,GUMAP算法由于没有进行步长自适应调整,两个模块对应的目标函数值(虚线)会在平衡点附近不断振荡。AD-GUAMP算法则由于增加步长调整机制,两个模块的目标函数值(实线)均稳定上升,经一定迭代次数最后收敛到某个值。此外,从图 2还可以看到,对同样的仿真场景,GUAMP算法两个模块的目标函数值均低于AD-GUAMP算法的两个对应模块。可以预见,采用自适应步长方法后,AD-GUAMP算法性能将优于GUAMP算法性能。
|
| 图 3 模块A与模块B目标函数-JA和-JB随迭代次数的变化曲线 |
2.2 低秩矩阵
图 4展示了随着矩阵秩比的变化算法NMSE(2 bit量化)或dNMSE(1 bit量化)的结果。图 4a和4c中,信噪比为20 dB;图 4b和4d中,信噪比为40 dB。图 4a和4b对应于过采样,过采样因子为δ=4,稀疏度λ=0.15。图 4c和4d对应于欠采样,过采样因子为δ=0.8,稀疏度λ=0.05。从图 4a可以看出,在1 bit量化下,GVAMP算法的NMSE相对于AD-GUAMP算法低约0.4 dB,表明GVAMP算法的性能最佳,AD-GUAMP算法略逊于GVAMP算法,但优于GUAMP和AD-GAMP算法。当秩比r>0.6时,AD-GAMP算法优于GUAMP算法。然而,随着秩比的减小,即矩阵变得更低秩,AD-GAMP算法的发散现象严重,此时GUAMP算法的性能优于AD-GAMP算法。在2 bit量化下,GVAMP算法的NMSE相对于AD-GUAMP算法低约0.25 dB,表明GVAMP算法稍优于AD-GUAMP算法,但GUAMP算法仍优于AD-GAMP算法。图 4b、4c和4d的结果与图 4a基本一致,即矩阵接近满秩时,GUAMP算法的性能仍高于AD-GAMP算法,随着矩阵秩比的减小,GUAMP算法的性能略低于AD-GAMP算法。在欠采样情况下,AD-GUAMP算法和GVAMP算法的性能非常接近。总之,随着秩比的减小,所有算法的性能都会恶化,但在低秩情况下,AD-GAMP算法的恶化现象更为严重,其性能明显低于其他算法,而其他算法在不同秩比下的性能相近。相对于GVAMP和AD-GUAMP算法,GUAMP算法有轻微的性能损失。
|
| 图 4 NMSE、dNMSE随矩阵秩比r的变化曲线 |
2.3 条件数病态矩阵
图 5展示了随着条件数κ的变化算法NMSE或dNMSE的变化曲线。图 5a和5c中,信噪比为20 dB;图 5b和5d中,信噪比为40 dB。图 5a和5b对应于过采样,过采样因子为δ=4,稀疏度λ=0.15。图 5c和5d对应于欠采样,过采样因子为δ=0.8,稀疏度λ=0.05。从图 5a可以看出,在1 bit量化下,GVAMP算法的dNMSE相对于AD-GUAMP算法低约0.6 dB,其重构性能最佳。其次是AD-GUAMP算法,两者的性能相近,并且均优于AD-GAMP和GUAMP算法。对于所有测试条件数,GUAMP算法的性能均优于AD-GAMP算法。在2 bit量化下,结果与1 bit量化场景类似。图 5b(过采样、高信噪比)、5c(欠采样、低信噪比)和5d(欠采样、高信噪比)的结果与图 5a基本一致。在欠采样情况下,AD-GUAMP算法和GVAMP算法的性能非常接近。同时,在高信噪比情况下,AD-GAMP算法相对于其他算法的性能恶化更为严重。总之,随着条件数κ的增加,所有算法的重构性能都变差,其中AD-GAMP算法最不稳健,其次是GUAMP算法。AD-GUAMP算法的重构性能非常接近GVAMP算法,尤其在低过采样因子和低稀疏度的情况下,两者的重构性能几乎完全一致。
|
| 图 5 NMSE、dNMSE随条件数κ的变化曲线 |
2.4 相关矩阵
图 6给出了算法NMSE或dNMSE随相关系数ρ的变化曲线。各分图仿真设置与图 5情况相同。从图 6a可以看出,在过采样和1 bit量化的情况下,AD-GUAMP算法和GUAMP算法的性能相当,均优于AD-GAMP算法和GVAMP算法。此外,当测量矩阵的相关性较低时,GUAMP算法的性能略逊于AD-GAMP算法。随着测量矩阵的相关性增加,AD-GAMP算法的性能急剧恶化,而GUAMP算法的性能优于AD-GAMP算法。当矩阵相关系数ρ>0.1时,GVAMP算法的性能开始急剧恶化。在2 bit量化下,AD-GUAMP算法的性能优于GUAMP算法,且两者均优于GVAMP算法。当测量矩阵相关系数为0<ρ<0.1时,GVAMP算法发散。图 6b(过采样、高信噪比)、6c(欠采样、低信噪比)和6d(欠采样、高信噪比)的结果与图 6a基本一致。总之,当相关系数ρ较小时(ρ<0.05),所有算法的估计性能相近。随着测量矩阵相关性的增加,GUAMP和AD-GUAMP算法的稳健性最高,并且在大多数情况下,AD-GUAMP算法优于GUAMP算法。这表明本文所提出的AD-GUAMP算法的自适应步长衰减机制提升了GUAMP算法的稳健性。
|
| 图 6 NMSE、dNMSE随相关系数ρ的变化曲线 |
2.5 运行时间
表 1给出了在特定条件下各算法的平均运行时间。仿真设置在欠采样场景,其中欠采样因子与之前的仿真设置相同,信噪比为40 dB。不同测量矩阵的参数设置为:低秩矩阵中秩比r=0.8,条件数矩阵中条件数κ=10,相关矩阵中相关系数ρ=0.1。除了上述条件外,其他仿真参数与之前的仿真参数相同。表 1中的时间是通过进行50次Monte Carlo实验并取平均值得到的。
| s | |||||||||||||||||||||||||||||
| 算法 | 相关矩阵 | 低秩矩阵 | 条件数矩阵 | ||||||||||||||||||||||||||
| 1 bit | 2 bit | 1 bit | 2 bit | 1 bit | 2 bit | ||||||||||||||||||||||||
| AD-GAMP | 0.11 | 0.15 | 0.10 | 0.09 | 0.09 | 0.08 | |||||||||||||||||||||||
| GVAMP | 0.09 | 0.15 | 0.08 | 0.07 | 0.07 | 0.07 | |||||||||||||||||||||||
| GUAMP | 0.15 | 0.28 | 0.14 | 0.15 | 0.13 | 0.14 | |||||||||||||||||||||||
| AD-GUAMP | 0.21 | 0.33 | 0.19 | 0.18 | 0.18 | 0.18 | |||||||||||||||||||||||
表 1中的数据表明,AD-GUAMP算法的平均运行时间大约是AD-GAMP算法的2倍。这是由于AD-GUAMP算法需要同时运行GAMP模块和AMP模块,相比之下AD-GAMP算法只需要运行单个GAMP模块,因此AD-GUAMP算法的运行时间会更长。此差异可能是由于AD-GUAMP算法的复杂性和额外的计算步骤所致。在实际应用中,需要综合考虑算法的性能和运行时间之间的权衡,并根据测量矩阵性质选择合适的算法。
3 结论为了解决GUAMP算法容易振荡且在某些情况下性能略低于AD-GAMP算法和GVAMP算法的问题,本文针对GUAMP算法的两个模块分别构建了近似目标函数。通过监测目标函数的下降情况,并按照预设规则对两个模块的参数进行阻尼衰减调整,提出了AD-GUAMP算法,以提升算法的稳健性。通过大量的仿真试验,本文验证了AD-GUAMP算法的有效性。结果显示,在低秩或高条件数测量矩阵的情况下,AD-GUAMP算法的性能几乎与GVAMP相当,并且优于AD-GAMP和GUAMP算法。在相关矩阵场景下,AD-GUAMP算法的性能也优于AD-GAMP、GUAMP和GVAMP算法。本研究结果表明,AD-GUAMP算法在改进GUAMP算法的稳健性方面取得了显著效果。通过引入近似目标函数和参数调整机制,AD-GUAMP算法能够更好地适应不同的场景,并具有出色的性能表现。本文的研究为Bayes稀疏信号重构算法在稳健性和准确性提升方面提供了有益的补充。
| [1] |
XIONG Y Z, WEI N, ZHANG Z P. A low-complexity iterative GAMP-based detection for massive MIMO with low-resolution ADCs[C]//Proceedings of 2017 IEEE Wireless Communications and Networking Conference (WCNC). San Francisco, USA, 2017: 1-6.
|
| [2] |
RANGAN S. Generalized approximate message passing for estimation with random linear mixing[C]//Proceedings of 2011 IEEE International Symposium on Information Theory Proceedings. Saint Petersburg, Russia, 2011: 2168-2172.
|
| [3] |
TIPPING M E. Sparse Bayesian learning and the relevance vector machine[J]. The Journal of Machine Learning Research, 2001, 1: 211-244. |
| [4] |
WIPF D P, RAO B D. Sparse Bayesian learning for basis selection[J]. IEEE Transactions on Signal Processing, 2004, 52(8): 2153-2164. DOI:10.1109/TSP.2004.831016 |
| [5] |
RANGAN S, SCHNITER P, FLETCHER A K, et al. On the convergence of approximate message passing with arbitrary matrices[J]. IEEE Transactions on Information Theory, 2019, 65(9): 5339-5351. DOI:10.1109/TIT.2019.2913109 |
| [6] |
SCHNITER P, RANGAN S. Compressive phase retrieval via generalized approximate message passing[J]. IEEE Transactions on Signal Processing, 2015, 63(4): 1043-1055. DOI:10.1109/TSP.2014.2386294 |
| [7] |
VILA J, SCHNITER P, RANGAN S, et al. Adaptive damping and mean removal for the generalized approximate message passing algorithm[C]//Proceedings of 2015 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP). South Brisbane, Australia, 2015: 2021-2025.
|
| [8] |
SCHNITER P, RANGAN S, FLETCHER A K. Vector approximate message passing for the generalized linear model[C]//Proceedings of the 50th Asilomar Conference on Signals, Systems and Computers. Pacific Grove, USA, 2016: 1525-1529.
|
| [9] |
RANGAN S, SCHNITER P, FLETCHER A K. Vector approximate message passing[J]. IEEE Transactions on Information Theory, 2019, 65(10): 6664-6684. DOI:10.1109/TIT.2019.2916359 |
| [10] |
RUAN C Y, ZHANG Z C, JIANG H, et al. Vector approximate message passing with sparse Bayesian learning for Gaussian mixture prior[J]. China Communications, 2023, 20(5): 57-69. DOI:10.23919/JCC.2023.00.005 |
| [11] |
MA J J, PING L. Orthogonal AMP[J]. IEEE Access, 2017, 5: 2020-2033. DOI:10.1109/ACCESS.2017.2653119 |
| [12] |
GUO Q H, XI J T. Approximate message passing with unitary transformation[Z/OL]. (2015-04-19)[2023-06-27]. https://doi.org/10.48550/arXiv.1504.04799.
|
| [13] |
LUO M, GUO Q H, JIN M, et al. Unitary approximate message passing for sparse Bayesian learning[J]. IEEE Transactions on Signal Processing, 2021, 69: 6023-6039. DOI:10.1109/TSP.2021.3114985 |
| [14] |
YUAN Z D, GUO Q H, LUO M. Approximate message passing with unitary transformation for robust bilinear recovery[J]. IEEE Transactions on Signal Processing, 2021, 69: 617-630. DOI:10.1109/TSP.2020.3044847 |
| [15] |
ZHU J, MENG X M, LEI X P, et al. A unitary transform based generalized approximate message passing[C]//Proceedings of 2023 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP). Rhodes Island, Greece, 2023: 1-5.
|
| [16] |
MENG X M, WU S, ZHU J. A unified Bayesian inference framework for generalized linear models[J]. IEEE Signal Processing Letters, 2018, 25(3): 398-402. DOI:10.1109/LSP.2017.2789163 |
| [17] |
ZHU J. A comment on the "A unified Bayesian inference framework for generalized linear models"[Z/OL]. (2019-04-09)[2023-07-14]. https://doi.org/10.48550/arXiv.1904.04485.
|
| [18] |
MINKA T P. A family of algorithms for approximate Bayesian inference[D]. Cambridge, USA: Massachusetts Institute of Technology, 2001.
|
| [19] |
SEEGER M. Expectation propagation for exponential families[R/OL]. 2005. https://core.ac.uk/download/pdf/147968132.pdf.
|
| [20] |
SHIU D S, FOSCHINI G J, GANS M J, et al. Fading correlation and its effect on the capacity of multielement antenna systems[J]. IEEE Transactions on Communications, 2000, 48(3): 502-513. DOI:10.1109/26.837052 |



