Frontiers in New-Quality Communication Technology

Multiuser detection algorithm based on an efficient Laplacian scale mixture prior

  • Pingping CHEN 1, 2 ,
  • Wei LIN 1, 2 ,
  • Changwei SHI 2 ,
  • Yukai FENG 1, 2 ,
  • Zhijian LIN , 1, * ,
  • Yi FANG 3
Expand
  • 1. School of Advanced Manufacturing, Fuzhou University, Quanzhou 362251, China
  • 2. Fujian Provincial Key Laboratory of Intelligent Processing and Wireless Transmission of Media Information, Fuzhou University, Fuzhou 350108, China
  • 3. College of Information Engineering, Guangdong University of Technology, Guangzhou 510006, China

Received date: 2025-02-26

  Online published: 2025-11-07

Copyright

All rights reserved. Unauthorized reproduction is prohibited.

Abstract

Objective: With the rapid evolution of fifth-generation (5G) mobile communication technologies, massive machine-type communication (mMTC) has become a pivotal application scenario in modern networks. This paradigm shift presents significant challenges in multiuser detection, particularly due to the exponential growth in user connections and heightened signal activity. Traditional orthogonal multiple access schemes, while ensuring minimal interuser interference, inherently limit the number of supported users by relying on orthogonal resource allocation, thereby failing to meet the scalability demands of mMTC. Consequently, grant-free nonorthogonal multiple access has emerged as a key enabler for Internet of Things communications, allowing nonorthogonal data superposition on limited resource blocks to enhance access capacity. However, existing sparse Bayesian learning (SBL) algorithms—although capable of achieving optimal sparse solutions—suffer from high computational complexity, primarily due to matrix inversion operations during expectation-maximization iterations. This complexity impedes real-time deployment in large-scale mMTC systems. To address this gap, this work proposes a novel SBL framework leveraging an efficient Laplace scale mixture (ELSM) prior, aiming to simultaneously enhance detection performance, reduce computational overhead, and adapt to dynamic multimeasurement scenarios. Methods: This paper proposes an ELSM-SBL algorithm to overcome the limitations of conventional SBL methods. First, a hierarchical Bayesian model is constructed using a Laplace scale mixture prior, which leverages the sharp peaks and heavy-tailed properties of Laplace distributions to promote sparsity and robustness against outliers. To avoid computationally expensive matrix inversions, a surrogate function is introduced to approximate the Gaussian likelihood function. This approximation is optimized within a majorization-minimization (MM) framework, where a block coordinate descent (BCD) algorithm solves the resulting nonconvex optimization problem. For single measurement vector (SMV) scenarios, the ELSM-SBL-SMV algorithm optimizes hyperparameters via evidence maximization, while an MM framework with BCD resolves nonconvexity in the joint cost function. For multiple measurement vector (MMV) scenarios, the ELSM-SBL-MMV scheme exploits temporal correlations among active user sets across consecutive time slots by sharing sparsity-controlling hyperparameters, thereby enhancing reconstruction performance. Results: Extensive simulations were conducted under mMTC settings with a total user count of K=108, subcarriers N=72, and active users M=12 using BPSK modulation and repeated over 1, 000 trials. For MMV scenarios, the number of measurement vectors was set to T=7. Compared with state-of-the-art SBL algorithms (e.g., GIG-SBL, BGIG-SBL, and LSM-SBL), the proposed ELSM-SBL-SMV algorithm can achieve a performance gain of about 2 dB, while the ELSM-SBL-MMV algorithm can achieve a gain of 3 dB. Meanwhile, the computational complexity of ELSM-SBL-SMV is $ {\cal{O}}$(K2), which is superior to existing SBL algorithms, demonstrating the superiority of the proposed algorithm in terms of performance and efficiency. Conclusions: The ELSM-SBL algorithm approximates the Gaussian likelihood function of the model by introducing a surrogate function, avoiding matrix inversion and reducing algorithm complexity, thereby considerably improving multiuser detection in mMTC systems. The SMV and MMV extensions demonstrate robust performance gains of 2 and 3 dB, respectively, while achieving $ {\cal{O}}$(K2) complexity, which is optimal for large-scale deployments. The experimental results confirm the algorithm's superiority in terms of BER, AER, and runtime, making it a viable solution for 5G and beyond-5G networks. Future work will explore real-time hardware implementations and extensions to massive Multiple-Input Multiple-Output scenarios to further enhance applicability in dynamic wireless environments.

Cite this article

Pingping CHEN , Wei LIN , Changwei SHI , Yukai FENG , Zhijian LIN , Yi FANG . Multiuser detection algorithm based on an efficient Laplacian scale mixture prior[J]. Journal of Tsinghua University(Science and Technology), 2025 , 65(11) : 2067 -2079 . DOI: 10.16511/j.cnki.qhdxxb.2025.27.048

随着第五代移动通信技术的发展,大规模机器类通信(massive machine type communication,mMTC)已成为5G网络的关键应用场景之一[1]。随之而来的挑战在于:如何设计高效的多用户检测算法以应对mMTC场景下用户数量的激增和信号活动性高的问题[2-3]。传统的正交多址接入(orthogonal multiple access,OMA)方案,通过对资源的正交复用来减少用户之间的相互干扰,提高通信的稳定性。
然而,系统所支持的用户连接数量严格受到正交资源数量的限制,无法满足mMTC通信场景的需求[4-5]。因此,免调度非正交多址接入(grant-free non-orthogonal multiple access,GF-NOMA)技术应运而生,被认为是物联网通信关键技术之一[6-7],将用户的传输数据以非正交的形式叠加到有限的资源块上,提高接入数量,而基站则通过多用户检测(multi-user detection,MUD)和编译码算法恢复多用户信息[4, 8]
在mMTC通信场景中,活跃用户数在繁忙的时段通常不会超过总用户数的10%[9],这一稀疏特性为多用户检测提供了优化空间。通过利用压缩感知(compressive sensing,CS)算法,可有效提高用户活跃性的检测性能[10]。Shim等[11]将正交匹配追踪(orthogonal matching pursuit,OMP) 算法与线性最小均方误差检测算法的优势相结合,在低复杂度下能够实现良好的检测性能。Monsees等[12]提出将基于压缩感知的多用户检测方法与结合CDMA的正交多载波系统相结合,提高了时频资源的灵活性和频谱效率。Wang等[13]利用压缩采样匹配追踪(compressive sampling matching pursuit,CoSaMP)算法先进行活跃度检测,再利用匹配追踪(matching pursuit,MP)进行数据检测,所提CS-MPA算法的性能优于CoSaMP算法[14]和MP算法。上述贪婪类算法的计算效率较高,但在信号重构的鲁棒性和精度方面却相对较低。
稀疏Bayes学习(sparse Bayesian learning,SBL)算法能够获得最优稀疏解,引起了广泛的关注[15]。SBL通过对模型参数施加先验分布的方式来实现稀疏性,由一组超参数进行控制[16]。Wipf等[17]通过期望最大化(expectation maximization,EM)算法更新超参数,取得良好的信号重建效果,但矩阵求逆的操作大大降低了运算效率。Tipping等[18]利用贪婪算法思想,逐步添加和删除基函数的方法来减少计算时间,但是仅适用于稳态噪声场景。Xiao等[19]通过Stackelberg博弈模型揭示了攻击者-防御者的策略交互,为对抗环境下检测算法的设计提供了理论支撑。Zhang等[20]和Abebe等[21]利用活跃用户的相邻时隙传输特性将单测量向量(single measurement vector,SMV)通信推广到更通用的多测量向量(multiple measurement vector,MMV)场景。Zhang等[22]将块稀疏结构模型与正交匹配追踪算法相结合提出了块正交匹配追踪(block OMP,BOMP)算法,应用在GF-NOMA场景,算法复杂度较高。Al-Shoukairi等[23]提出基于消息传递的块稀疏Bayes学习(message passing based block SBL,MP-BSBL)算法,将置信传播(belief propagation,BP)和平均场(mean field,MF)理论结合,对概率分布进行Gauss近似以降低算法复杂度。当前的SBL算法大多是基于Gauss逆Gamma(Gaussian inverse Gamma,GIG)先验模型[24-25],模型中的多随机变量可能导致尾部较厚的分布,难以用Gauss混合函数的先验分布来准确拟合。
相比之下,Laplace先验分布具有尖锐的峰值和重尾特性,更容易提取远离均值的样本,促进信号的稀疏性并对异常值具有更高的鲁棒性[26]。因此,本文提出了基于高效Laplace尺度混合(efficient Laplacian scale mixture,ELSM)先验的ELSM-SBL算法,引入代理函数来近似模型的Gauss似然函数,从而规避矩阵求逆。接着,通过证据最大化方法优化决定模型稀疏度的超参数。引入数据依赖项构造惩罚项可分离的联合成本函数。最后,在主要—最小化(majorization-minimization,MM)框架中,使用块坐标下降(block coordinate descent,BCD)方法解决联合成本函数产生的非凸优化问题。此外,本文利用活跃用户集在相邻时隙相关性,共享一个控制稀疏度的超参数来增强该重构性质,推广到更符合实际通信的MMV场景。实验结果表明,本文提出的ELSM-SBL算法在SMV场景和MMV场景下的检测性能,相较SBL算法分别有2 dB和3 dB的性能增益。同时所提ELSM-SBL-SMV算法能够实现$ {\cal{O}}$(K2)的计算复杂度,优于目前主流的SBL算法。

1 系统模型

考虑到通信场景的一般性,本文考虑一个上行免调度NOMA通信系统。该系统由一个基站和K个用户组成,每个用户与基站均配备一根天线用于发送和接收信号,系统模型如图 1所示。鉴于mMTC多用户通信场景中活跃用户的稀疏特性[16],假设某一个时刻处于活跃状态的用户有m个,则满足关系mK
图 1 上行免调度NOMA通信系统模型
在帧稀疏的通信场景下,用户以帧为基本时间单位传输信号,一帧包括T个时隙。相应地,用户的活跃状态也以帧为基本单位进行描述,即活跃用户在一帧的时间内连续保持活跃状态,不活跃的用户则一直处于不活跃状态。第k个用户在t时刻的活跃状态可以表示为
$u_{k}(t)= \begin{cases}0, & k \notin \mathcal{K}_{\mathrm{a}} \text { 且 } t=1, \cdots, T ; \\ 1, & k \in \mathcal{K}_{\mathrm{a}} \text { 且 } t=1, \cdots, T .\end{cases}$
其中$ \mathcal{K}_{\mathrm{a}}$表示活跃用户的索引集合。根据用户的活跃状态,其发送数据符号xk(t)可以表示为
$\begin{cases}x_{k}(t)=0, & u_{k, t}=0 ; \\ x_{k}(t) \in \mathcal{X}, & u_{k, t}=1 .\end{cases}$
其中$ \mathcal{X}$表示归一化星座图。
考虑到大规模用户随机接入的场景[7, 27],该系统采用过载的方式进行数据传输,即子载波个数NK。基站在第n个子载波上接收到的t时刻信号表示为
$\begin{gather*}y_{n}(t)=\sum\limits_{k=1}^{K} g_{n k}(t) s_{n k} x_{k}(t)+w_{n}(t) ,\\n=1, \cdots, N .\end{gather*}$
其中:snk表示第k个用户在第n个子载波上的扩频分量。gnk(t)表示t时刻第k个用户在第n个子载波的信道增益,服从独立同分布gnk(t)~N(0, 1)。wn(t)表示在t时刻第n个子载波上的加性白噪声,服从均值为零、方差为σ2的Gauss分布。
为体现系统零星通信的特点,定义$ \boldsymbol{H}_{t} \in \mathbb{C}^{N \times K}$t时刻的等效信道矩阵,即t时刻的信道增益与扩频序列融合成等效信道Ht。在t时刻的等效信道Ht中第n行,第k列元素hnk(t)可以表示为
$\begin{equation*}h_{n k}(t)=g_{n k}(t) \cdot s_{n k} \cdot \end{equation*}$
因此,基站接收端合并N个子载波的接收信号,t时刻的接收信号向量可以表示为
$\boldsymbol{y}_{t}=\boldsymbol{H}_{t} \boldsymbol{x}_{t}+\boldsymbol{w}_{t} . $
其中:yt=(y1(t), y2(t), …, yN(t))T表示t时刻基站接收到的信号向量;xt=(x1(t), x2(t), …, xK(t))T表示t时隙发送的信号向量;wt=(w1(t), w2(t), …, wN(t))T表示t时隙的Gauss噪声向量。此外,假设多用户在多时隙传输的活跃用户索引集合保持不变,这有利于增强信号的联合稀疏特性,并提升稀疏信号的重构性能[27]

2 LSM-SBL算法

由于Laplace分布在零点附近具有尖锐的峰值和显著重尾特性,相较于Gauss分布具有更高的鲁棒性,在稀疏信号的模型建立中更为适用。考虑到系统通信通常具有零星通信的特点,仅有少量的用户处于活跃状态,符合稀疏信号的特点。因此,本文采用基于Laplace尺度混合先验(Laplacian scale mixture,LSM)的稀疏Bayes学习(LSM-SBL)模型来提升模型的准确性[26]

2.1 LSM-SBL的分层模型

对于接收信号y,本文采用分层Bayes模型进行建模。由于信号模型的观测噪声为加性Gauss白噪声,即服从均值为零,方差为σ2的Gauss分布。因此,似然函数p(y|x, λ)~$ {\cal{N}}$(Hx, λ-1I),其中λ=σ-2,表示噪声信号方差的倒数;I表示N×N的单位矩阵,即
$p(\boldsymbol{y} \mid \boldsymbol{x}, \lambda)=\left(\frac{\lambda}{2 {\rm{ \mathsf{ π}}}}\right)^{\frac{N}{2}} \exp \left(-\frac{\lambda}{2}\|\boldsymbol{y}-\boldsymbol{H} \boldsymbol{x}\|_{2}^{2}\right) . $
由于Gauss分布方差倒数的共轭先验分布为Gamma分布,故假设λ的超先验服从给定参数的Gamma分布[15],即λ的概率密度函数p(λ)可表示为
$p(\lambda)=\frac{b_{0}^{a_{0}}}{\boldsymbol{\varLambda}\left(a_{0}\right)} \lambda^{a_{0}-1} \exp \left(-b_{0} \lambda\right) . $
其中:$ \boldsymbol{\varLambda}(a)=\int_{0}^{\infty} x^{a-1} \exp (-x) \mathrm{d} x, a_{0}$表示尺度参数,b0表示形状参数,且a0>0,b0>0。使用l1范数求解信号x正则化的过程相当于对x作Laplace先验[26]
$p(\boldsymbol{x} \mid \beta)=\frac{\beta}{2} \exp \left(-\frac{\beta}{2}\|\boldsymbol{x}\|_{1}\right) . $
其中β是Laplace分布的尺度参数。令τ=β/λ,可利用最大后验(maximum a posteriori,MAP)估计即可求解出式(6)和式(7)。然而,式(8)的Laplace先验与式(6)的Gauss分布并非共轭先验分布,导致后验分布的参数无法直接通过更新先验参数获得,使得Bayes迭代计算过程非常复杂。
为了简化Bayes计算过程,下文引入分层先验进行建模[23]。在第一层的Bayes模型中,假设p(x|γ)服从Gauss分布,即
$p(\boldsymbol{x} \mid \boldsymbol{\gamma})=\prod\limits_{i=1}^{K} \mathcal{N}\left(x_{i} \mid 0, \gamma_{i}\right) . $
其中γ =(γ1, γ2, …, γK),γ是决定x先验分布的超参数。
在第二层中,通过引入超参数β控制γi,其中i=1, 2, …, K。假设p(γi)服从参数β的指数分布,则有
$p\left(\gamma_{i}\right)=\frac{\beta}{2} \exp \left(-\frac{\beta \gamma_{i}}{2}\right) . $
其中: γi≥0, β≥0。
根据Bayes联合概率分布公式,则
$\begin{gather*}p(\boldsymbol{x} \mid \beta)=\int p(\boldsymbol{x} \mid \boldsymbol{\gamma}) p(\boldsymbol{\gamma}) \mathrm{d} \boldsymbol{\gamma}= \\\prod\limits_{i=1}^{K} \int p\left(x_{i} \mid \gamma_{i}\right) p\left(\gamma_{i}\right) \mathrm{d} \gamma_{i}= \\\frac{\beta^{K / 2}}{2^{K}} \exp \left(-\sqrt{\beta} \sum\limits_{i=1}^{K}\left|x_{i}\right|\right) . \end{gather*}$
最后,假设超参数β服从Gamma先验,则
$p(\beta \mid \theta)=\varGamma(\beta \mid \theta / 2, \theta / 2) . $
上述的分层Laplace模型如图 2所示。前两层的式(9)和式(10)是为了得到Laplace分布,第三层用于计算参数β
图 2 分层Laplace框图

2.2 LSM-SBL推断

Bayes推断是基于后验分布进行的,即
$p(\boldsymbol{x}, \boldsymbol{\gamma}, \lambda, \beta \mid \boldsymbol{y})=\frac{p(\boldsymbol{x}, \boldsymbol{\gamma}, \lambda, \beta, \boldsymbol{y})}{p(\boldsymbol{y})}, $
$p(\boldsymbol{y})=\iiint \int p(\boldsymbol{x}, \boldsymbol{\gamma}, \lambda, \beta, \boldsymbol{y}) \mathrm{d} \boldsymbol{x} \mathrm{~d} \boldsymbol{\gamma} \mathrm{~d} \lambda \mathrm{~d} \beta .$
但在实际应用中,由于边际似然函数p(y)无法通过积分直接获得,使得后验分布p(x, γ, λ, β|y)难以计算。因此,通常将后验分布分解为
$\begin{gather*}p(\boldsymbol{x}, \boldsymbol{\gamma}, \lambda, \beta \mid \boldsymbol{y})= \\p(\boldsymbol{x} \mid \boldsymbol{y}, \boldsymbol{\gamma}, \beta, \lambda) p(\boldsymbol{\gamma}, \beta, \lambda \mid \boldsymbol{y}) . \end{gather*}$
其中:变量x的后验分布p(x|y, γ, β, λ)服从均值为μ、方差为Σ的多元Gauss分布$ {\cal{N}}$(x|μ, Σ),有
$\hat{\boldsymbol{\mu}}=\boldsymbol{\varSigma} \lambda \boldsymbol{H}^{\mathrm{T}} \boldsymbol{y}, $
$\boldsymbol{\varSigma}=\left[\lambda \boldsymbol{H}^{\mathrm{T}} \boldsymbol{H}+\boldsymbol{\varLambda}\right]^{-1} . $
其中Λ=diag(1/γi)。
根据Bayes定理,式(15)的p(γ, β, λ|y)满足
$\begin{gather*}p(\boldsymbol{\gamma}, \beta, \lambda \mid \boldsymbol{y})= \\\frac{p(\boldsymbol{y}, \boldsymbol{\gamma}, \beta, \lambda)}{p(\boldsymbol{y})} \propto p(\boldsymbol{y}, \boldsymbol{\gamma}, \beta, \lambda) . \end{gather*}$
因此,将最大化后验概率分布p(γ, β, λ|y)估计超参数的方式转化为对联合概率分布p(y, γ, β, λ)求最大值。根据Bayes定理,则有
$\begin{gather*}p(\boldsymbol{y}, \boldsymbol{\gamma}, \beta, \lambda)= \\\int p(\boldsymbol{y} \mid \boldsymbol{x}, \beta) p(\boldsymbol{x} \mid \boldsymbol{\gamma}) p(\boldsymbol{\gamma} \mid \lambda) \cdot p(\lambda) p(\beta) \mathrm{d} \boldsymbol{x}= \\\left(\frac{1}{2 {\rm{ \mathsf{ π}}}}\right)^{K / 2}|\boldsymbol{C}|^{-1 / 2} \exp \left[-\frac{1}{2} \boldsymbol{y}^{\mathrm{T}} \boldsymbol{C}^{-1} \boldsymbol{y}\right] \cdot \\p(\boldsymbol{\gamma} \mid \lambda) p(\lambda) p(\beta) . \end{gather*}$
其中C =(λ-1I + -1HT)。
对式(19) 中的联合概率分布p(y, γ, β, λ)左右两边取对数,则有
$\begin{gather*}\mathcal{L}=-\frac{1}{2} \log|\boldsymbol{C}|-\frac{1}{2} \mathrm{y}^{\mathrm{T}} \boldsymbol{C}^{-1} \boldsymbol{y}+K \log \frac{\beta}{2}- \\\frac{\beta}{2} \sum\limits_{i=1}^{K} \gamma_{i}+\frac{\theta}{2} \log \frac{\theta}{2}-\log \varGamma(\theta / 2)+\left(\frac{\theta}{2}-1\right) \log \beta- \\\frac{\theta}{2} \beta+\left(a_{0}-1\right) \log \lambda-b_{0} \lambda .\end{gather*}$
因此,将对联合概率分布求最大值的形式转换为对数函数求最大值。通过对C取绝对值并利用行列式恒等式,可得
$\begin{gather*}|\boldsymbol{C}|=|\boldsymbol{\varLambda}|^{-1}\left|\lambda^{-1} \boldsymbol{I} \|\boldsymbol{\varLambda}+\lambda \boldsymbol{H}^{\mathrm{T}} \boldsymbol{H}\right|= \\|\boldsymbol{\varLambda}|^{-1}\left|\lambda^{-1} \boldsymbol{I} \|\boldsymbol{\varSigma}^{-1}\right|, \end{gather*}$
$\begin{equation*}\log|\boldsymbol{C}|=-\log|\boldsymbol{\varLambda}|-K \log \lambda-\log|\boldsymbol{\varSigma}|. \end{equation*}$
此外,利用Woodbury公式[28],则有
$\begin{gather*}\boldsymbol{C}^{-1}=\left(\lambda^{-1} \boldsymbol{I}+\boldsymbol{H} \boldsymbol{\varLambda}^{-1} \boldsymbol{H}^{\mathrm{T}}\right)^{-1}= \\\lambda \boldsymbol{I}-\lambda \boldsymbol{H}\left(\boldsymbol{\varLambda}+\lambda \boldsymbol{H} \boldsymbol{H}^{\mathrm{T}}\right)^{-1} \boldsymbol{H}^{\mathrm{T}} \lambda= \\\lambda \boldsymbol{I}-\lambda \boldsymbol{H} \boldsymbol{\varSigma} \boldsymbol{H}^{\mathrm{T}} \lambda . \end{gather*}$
因此,
$\begin{gather*}\boldsymbol{y}^{\mathrm{T}} \boldsymbol{C}^{-1} \boldsymbol{y}=\lambda \boldsymbol{y}^{\mathrm{T}} \boldsymbol{y}-\lambda \boldsymbol{y}^{\mathrm{T}} \boldsymbol{H} \boldsymbol{\varSigma} \boldsymbol{H}^{\mathrm{T}} \lambda \boldsymbol{y}= \\\lambda \boldsymbol{y}^{\mathrm{T}}(\boldsymbol{y}-\boldsymbol{H} \boldsymbol{\mu})= \\\lambda\|\boldsymbol{y}-\boldsymbol{H} \boldsymbol{\mu}\|_{2}^{2}+\lambda \boldsymbol{\mu}^{\mathrm{T}} \boldsymbol{H}^{\mathrm{T}}(\boldsymbol{y}-\boldsymbol{H} \boldsymbol{\mu})= \\\lambda\|\boldsymbol{y}-\boldsymbol{H} \boldsymbol{\mu}\|_{2}^{2}+\boldsymbol{\mu}^{\mathrm{T}} \boldsymbol{\varLambda} \boldsymbol{\mu} . \end{gather*}$
将式(21)—(24)代入式(20),并对$ {\cal{L}}$γi偏导,可得
$\begin{equation*}\frac{\mathrm{d} \mathcal{L}}{\mathrm{~d} \gamma_{i}}=\frac{1}{2}\left[-\frac{1}{\gamma_{i}}+\frac{\left\langle w_{i}^{2}\right\rangle}{\gamma_{i}^{2}}-\beta\right] . \end{equation*}$
其中〈wi2〉=μi2+ΣiiΣii表示Σ的第i个对角元素。令等式两边等于0,可得
$\hat{\gamma}_{i}=-\frac{1}{2 \beta}+\sqrt{\frac{1}{4 \beta^{2}}+\frac{\left\langle\omega_{i}^{2}\right\rangle}{\beta}} . $
同样,分别对$ {\cal{L}}$βλ偏导并令等式等于0,可得
$\hat{\beta}=\frac{K-1+\theta / 2}{\sum\limits_{i=1}^{K} \gamma_{i} / 2+\theta / 2}, $
$\hat{\lambda}=\frac{K / 2+a_{0}}{\left\langle\|\boldsymbol{y}-\boldsymbol{H} \boldsymbol{x}\|_{2}^{2}\right\rangle / 2+b_{0}} . $
图 2可以看出,变量xi由每个γi控制,而发送信号向量x是稀疏的,意味着向量x中许多元素值为0。为了促进向量的稀疏性并降低计算的复杂度,每次迭代通过判别的方法更新超参数γi,而不用式(26)对每个参数γi进行计算。
利用稀疏Bayes模型快速边际似然最大化算法[18],将式(19)中的C进行基分解,可得
$\begin{gather*}\boldsymbol{C}=\lambda^{-1} \boldsymbol{I}+\sum\limits_{i} \gamma_{i} \boldsymbol{\eta}_{i} \boldsymbol{\eta}_{i}^{\mathrm{T}}= \\\lambda^{-1} \boldsymbol{I}+\sum\limits_{j \neq i} \gamma_{j} \boldsymbol{\eta}_{j} \boldsymbol{\eta}_{j}^{\mathrm{T}}+\gamma_{i} \boldsymbol{\eta}_{i} \boldsymbol{\eta}_{i}^{\mathrm{T}}= \\\boldsymbol{C}_{-i}+\gamma_{i} \boldsymbol{\eta}_{i} \boldsymbol{\eta}_{i}^{\mathrm{T}} . \end{gather*}$
其中:ηi表示包括矩阵H中第i列所有元素的基向量,Ci表示删除第i个基函数后的C矩阵。利用Woodbury恒等式,可得
$\boldsymbol{C}^{-1}=\boldsymbol{C}_{-i}^{-1}-\frac{\boldsymbol{C}_{-i}^{-1} \boldsymbol{\eta}_{i} \boldsymbol{\eta}_{i}^{\mathrm{T}} \boldsymbol{C}_{-i}^{-1}}{1 / \gamma_{i}+\boldsymbol{\eta}_{i}^{\mathrm{T}} \boldsymbol{C}_{-i}^{-1} \boldsymbol{\eta}_{i}}, $
$|\boldsymbol{C}|=\left|\boldsymbol{C}_{-i}\right|\left|1+\gamma_i \boldsymbol{\eta}_i^{\mathrm{T}} \boldsymbol{C}_{-i}^{-1} \boldsymbol{\eta}_i\right|.$
将式(30)和式(31)代入式(20),可得
$\begin{gather*}\mathcal{L}(\boldsymbol{\gamma})=-\frac{1}{2}\left(\log \left|\boldsymbol{C}_{-i}\right|+\boldsymbol{y}^{\mathrm{T}} \boldsymbol{C}_{-i}^{-1} \boldsymbol{y}+\frac{\lambda}{2} \sum\limits_{j \neq i} \gamma_{j}\right)+ \\\frac{1}{2}\left(\log \frac{1}{1+\gamma_{i} s_{i}}+\frac{q_{i}^{2} \gamma_{i}}{1+\gamma_{i} s_{i}}-\lambda \gamma_{i}\right)= \\\mathcal{L}\left(\boldsymbol{\gamma}_{-i}\right)+l\left(\gamma_{i}\right) . \end{gather*}$
其中l(γi)=(1/2)[log(1/(1+γisi))+(qi2γi/(1+γisi))-λγi],并且qi, si分别定义为
$s_{i} =\boldsymbol{\eta}_{i}^{\mathrm{T}} \boldsymbol{C}_{-i}^{-1} \boldsymbol{\eta}_{i}, $
$q_{i} =\boldsymbol{\eta}_{i}^{\mathrm{T}} \boldsymbol{C}_{-i}^{-1} \boldsymbol{y} . $
同样对$ {\cal{L}}$γi的偏导,并令等式为0,可得到超参数γi的更新式为
$\hat{\gamma}_{i}= \begin{cases}\frac{-s_{i}\left(s_{i}+2 \beta\right)+s_{i} \sqrt{\Delta}}{2 \beta s_{i}^{2}}, & q_{i}^{2}-s_{i}>\lambda ; \\ 0, & \text { else. }\end{cases}$
其中Δ=(si+2β)2-4β(siqi2+β),参数qi, si同样可以利用添加或删除基函数的方式快速计算,则有
$S_{i}=\lambda \boldsymbol{\eta}_{i}^{\mathrm{T}} \boldsymbol{\eta}_{i}-\lambda^{2} \boldsymbol{\eta}_{i}^{\mathrm{T}} \boldsymbol{H} \boldsymbol{\varSigma} \boldsymbol{H}^{\mathrm{T}} \boldsymbol{\eta}_{i}, $
$Q_{i}=\lambda \boldsymbol{\eta}_{i}^{\mathrm{T}} \boldsymbol{y}-\lambda^{2} \boldsymbol{\eta}_{i}^{\mathrm{T}} \boldsymbol{H} \boldsymbol{\varSigma} \boldsymbol{H}^{\mathrm{T}} \boldsymbol{y}, $
$s_{i}=\frac{S_{i}}{1-\gamma_{i} S_{i}} q_{i}, $
$q_{i}=\frac{Q_{i}}{1-\gamma_{i} S_{i}} . $
LSM-SBL算法总结如图 3所示。

3 ELSM-SBL算法

现有的SBL算法在迭代过程中通常需要进行矩阵的求逆计算,如式(17)所示。这种K×K的矩阵求逆计算会导致$ {\cal{O}}$(K3)的计算复杂度,限制该算法在大规模数据集上的广泛应用。为了解决这一问题,LSM-SBL算法利用稀疏Bayes模型快速边际似然最大化算法进行基分解,减少计算所需要的时间。然而,这种简化方式需要对基向量做出强稀疏性的假设,这限制了该算法的实际应用。
本文在LSM-SBL框架上,提出了ELSM-SBL算法,通过引入代理函数近似模型的Gauss似然函数,从而规避矩阵的求逆运算。然后,通过证据最大化方法确定模型的超参数并构造一个可分离的联合成本函数。最后在MM框架中使用BCD算法解决由此产生的非凸优化问题。此外,本文利用发送信号在相邻连续时隙内具有一定的时间相关性这一特点,相邻时隙通过共享一个控制稀疏度的超参数来提升接收信号的重构性能,使得ELSM-SBL算法能够更好地适应MMV场景。

3.1 ELSM-SBL-SMV算法

由于基于单观测向量y进行信号恢复的SBL算法能够直接适用于SMV场景,因此在本文未特别说明的情况下,各种SBL算法同样指SBL-SMV算法。ELSM-SBL-SMV算法的概率图模型如图 4所示。
图 4 ELSM-SBL-SMV算法的概率图模型
根据Zhou等[29]的凸优化引理公式,可推导出Gauss似然函数具有严格的下界。因此,引入代理函数$ \hat{p}(\boldsymbol{y} \mid \boldsymbol{x}, \lambda ; \boldsymbol{\alpha})$来近似替代Gauss似然p(y|x, λ),式(7)可表示为
$p(\boldsymbol{y} \mid \boldsymbol{x}, \lambda)=\mathop {\max }\limits_{\boldsymbol{\alpha}} \hat{p}(\boldsymbol{y} \mid \boldsymbol{x}, \lambda ; \boldsymbol{\alpha}) . $
其中α为用于构建Gauss似然的上界辅助变量,α$ \mathbb{R}^n$
变量x的近似后验密度通过Bayes推理,可得
$p(\boldsymbol{x} \mid \boldsymbol{y}, \lambda, \boldsymbol{\gamma}) \approx \frac{\hat{p}(\boldsymbol{y} \mid \boldsymbol{\gamma}, \lambda ; \hat{\boldsymbol{\alpha}}) p(\boldsymbol{x} \mid \boldsymbol{\gamma})}{\int \hat{p}(\boldsymbol{y} \mid \boldsymbol{\gamma}, \lambda ; \hat{\boldsymbol{\alpha}}) p(\boldsymbol{x} \mid \boldsymbol{\gamma}) \mathrm{d} \boldsymbol{x}} . $
此时,有
$\begin{gather*}\hat{p}(\boldsymbol{y} \mid \boldsymbol{x}, \lambda ; \hat{\boldsymbol{\alpha}}) p(\boldsymbol{x} \mid \boldsymbol{\gamma})= \\\left(\frac{\lambda}{2 {\rm{ \mathsf{ π}}}}\right)^{\frac{N}{2}}|2 {\rm{ \mathsf{ π}}} \boldsymbol{\varLambda}|^{-\frac{1}{2}} \exp \left\{-\frac{1}{2}\left[\lambda R(\boldsymbol{x}, \hat{\boldsymbol{\alpha}})+\boldsymbol{x}^{\mathrm{T}} \boldsymbol{\varLambda}^{-1} \boldsymbol{x}\right]\right\} . \end{gather*}$
公式(42)中的$ \lambda R(\boldsymbol{x}, \hat{\boldsymbol{\alpha}})+\boldsymbol{x}^{\mathrm{T}} \boldsymbol{\varLambda}^{-1} \boldsymbol{x}$可具体表示为
$\begin{gather*}\lambda R(\boldsymbol{x}, \hat{\boldsymbol{\alpha}})+\boldsymbol{x}^{\mathrm{T}} \boldsymbol{\varLambda}^{-1} \boldsymbol{x}= \\\lambda\|\boldsymbol{y}-\boldsymbol{H} \hat{\boldsymbol{\alpha}}\|_{2}^{2}+2 \lambda(\boldsymbol{x}-\hat{\boldsymbol{\alpha}})^{\mathrm{T}} \boldsymbol{H}^{\mathrm{T}}(\boldsymbol{H} \hat{\boldsymbol{\alpha}}-\boldsymbol{y})+ \\\varepsilon \lambda\|\boldsymbol{x}-\hat{\boldsymbol{\alpha}}\|_{2}^{2}+\boldsymbol{x}^{\mathrm{T}} \boldsymbol{\varLambda}^{-1} \boldsymbol{x}= \\\boldsymbol{x}^{\mathrm{H}}\left(\boldsymbol{\varLambda}^{-1}+\varepsilon \lambda \boldsymbol{I}\right) \boldsymbol{x}-2 \lambda \boldsymbol{x}^{\mathrm{T}}\left(\varepsilon \hat{\boldsymbol{\alpha}}-\boldsymbol{H}^{\mathrm{T}} \boldsymbol{H} \hat{\boldsymbol{\alpha}}+\boldsymbol{H}^{\mathrm{T}} \boldsymbol{y}\right)+ \\\varepsilon \lambda \hat{\boldsymbol{\alpha}}^{\mathrm{T}} \hat{\boldsymbol{\alpha}}-\lambda \hat{\boldsymbol{\alpha}}^{\mathrm{T}} \boldsymbol{H}^{\mathrm{T}} \boldsymbol{H} \hat{\boldsymbol{\alpha}}+\lambda \boldsymbol{y}^{\mathrm{T}} \boldsymbol{y}= \\\left(\boldsymbol{x}-\lambda \boldsymbol{A}^{-1}\epsilon\right)^{\mathrm{T}} \boldsymbol{A}\left(\boldsymbol{x}-\lambda \boldsymbol{A}^{-1} \epsilon\right)+G(\lambda, \boldsymbol{\gamma}, \hat{\boldsymbol{\alpha}}), \end{gather*}$
其中:$ \boldsymbol{A}=\boldsymbol{\varLambda}^{-1}+\varepsilon \lambda \boldsymbol{I}, { }{\epsilon}=\varepsilon \hat{\boldsymbol{\alpha}}-\boldsymbol{H}^{\mathrm{T}} \boldsymbol{H} \hat{\boldsymbol{\alpha}}+\boldsymbol{H}^{\mathrm{T}} \boldsymbol{y}, G(\lambda$$ \boldsymbol{\gamma}, \hat{\boldsymbol{\alpha}}):=\varepsilon \lambda \hat{\boldsymbol{\alpha}}^{\mathrm{T}} \hat{\boldsymbol{\alpha}}-\lambda \hat{\boldsymbol{\alpha}}^{\mathrm{T}} \boldsymbol{H}^{\mathrm{T}} \boldsymbol{H} \hat{\boldsymbol{\alpha}}+\lambda \boldsymbol{y}^{\mathrm{T}} \boldsymbol{y}-\lambda^{2} \epsilon^{\mathrm{T}} \boldsymbol{A}^{-1} \epsilon $
由于x$ G(\lambda, \boldsymbol{\gamma}, \hat{\boldsymbol{\alpha}})$无关,则有
$\begin{gather*}\int \hat{p}(\boldsymbol{y} \mid \boldsymbol{x}, \lambda ; \hat{\boldsymbol{\alpha}}) p(\boldsymbol{x} \mid \boldsymbol{\gamma}) \mathrm{d} \boldsymbol{x}= \\\left(\frac{\lambda}{2 {\rm{ \mathsf{ π}}}}\right)^{\frac{N}{2}}|\boldsymbol{\varLambda} \boldsymbol{A}|^{-\frac{1}{2}} \exp \left\{-\frac{1}{2} G(\lambda, \boldsymbol{\gamma}, \hat{\boldsymbol{\alpha}})\right\} . \end{gather*}$
最后,将式(42)至式(44)代入式(41),可得
$\begin{gather*}\frac{\hat{p}(\boldsymbol{y} \mid \boldsymbol{x}, \lambda ; \hat{\boldsymbol{\alpha}}) p(\boldsymbol{x} \mid \boldsymbol{\gamma})}{\int \hat{p}(\boldsymbol{y} \mid \boldsymbol{x}, \lambda ; \hat{\boldsymbol{\alpha}}) p(\boldsymbol{x} \mid \boldsymbol{\gamma}) \mathrm{d} \boldsymbol{x}}= \\(2 {\rm{ \mathsf{ π}}})^{-\frac{{\rm{ \mathsf{ π}}}}{2}}|\boldsymbol{A}|^{\frac{1}{2}} \exp \left\{-\frac{1}{2}\left(\boldsymbol{x}-\lambda \boldsymbol{A}^{-1} \epsilon\right)^{\mathrm{T}} \boldsymbol{A}\left(\boldsymbol{x}-\lambda \boldsymbol{A}^{-1} \epsilon\right)\right\} \end{gather*}$
从式(45)的具体表达式中,可得到此时的方差Σ和均值μ分别为
$\boldsymbol{\varSigma}=(\boldsymbol{\varLambda}+\varepsilon \lambda \boldsymbol{I})^{-1}, $
$\hat{\boldsymbol{\mu}}=\lambda \boldsymbol{\varSigma}\left(\varepsilon \hat{\boldsymbol{\alpha}}-\boldsymbol{H}^{\mathrm{T}} \boldsymbol{H} \hat{\boldsymbol{\alpha}}+\boldsymbol{H}^{\mathrm{T}} \boldsymbol{y}\right) .$
其中ε=eig(HΤH)+ε0, ε0为某一大于零的常数。
接下来通过最大化联合概率分布估计超参数,可表示为
$(\lambda, \boldsymbol{\gamma}, \boldsymbol{\alpha})={\rm{arg}}\underset{\lambda, \gamma, \beta}{\operatorname{max}} \hat{p}(\boldsymbol{y}, \lambda, \boldsymbol{\gamma} ; \boldsymbol{\alpha}), $
$\begin{gather*}\hat{p}(\boldsymbol{y}, \lambda, \boldsymbol{\gamma} ; \boldsymbol{\alpha})= \\\int \hat{p}(\boldsymbol{y} \mid \boldsymbol{x}, \lambda ; \boldsymbol{\alpha}) p(\boldsymbol{x} \mid \boldsymbol{\gamma}) p(\lambda) p(\boldsymbol{\gamma}) \mathrm{d} \boldsymbol{x}=\\p(\lambda) p(\boldsymbol{\gamma})\left(\frac{\lambda}{2 {\rm{ \mathsf{ π}}}}\right)^{\frac{N}{2}}\left|\boldsymbol{\varLambda} \boldsymbol{\varSigma}^{-1}\right|^{-\frac{1}{2}} \exp \left\{-\frac{1}{2} S(\lambda, \boldsymbol{\gamma}, \boldsymbol{\alpha})\right\}.\end{gather*}$
其中$ S(\lambda, \gamma, \boldsymbol{\alpha})$可表示为
$\begin{gather*}S(\lambda, \boldsymbol{\gamma}, \boldsymbol{\alpha})= \\\lambda\left(\|\boldsymbol{y}-\boldsymbol{H} \boldsymbol{\alpha}\|_{2}^{2}-2 \boldsymbol{\alpha}^{\mathrm{T}} \boldsymbol{H}^{\mathrm{T}}(\boldsymbol{H} \boldsymbol{\alpha}-\boldsymbol{y})+\right. \\\left.\varepsilon \boldsymbol{\alpha}^{\mathrm{T}} \boldsymbol{\alpha}\right)-\boldsymbol{\mu}^{\mathrm{T}} \boldsymbol{\varSigma}^{-1} \boldsymbol{\mu}= \\\min _{{\boldsymbol{x}}} \lambda\left(\|\boldsymbol{y}-\boldsymbol{H} \boldsymbol{\alpha}\|_{2}^{2}+2(\boldsymbol{\mu}-\boldsymbol{\alpha})^{\mathrm{T}} \boldsymbol{H}^{\mathrm{T}}(\boldsymbol{H} \boldsymbol{\alpha}-\boldsymbol{y})+\right. \\\left.\varepsilon\|\boldsymbol{\mu}-\boldsymbol{\alpha}\|_{2}^{2}\right)+\boldsymbol{\mu}^{\mathrm{T}} \boldsymbol{\varLambda}^{-1} \boldsymbol{\mu} . \end{gather*}$
值得注意的是,数据相关项S(λ, γ, α)与x无关。对式(49)取对数,得到最小化成本函数,则有
$\begin{array}{c} \mathcal{L}(\boldsymbol{\alpha}, \lambda, \boldsymbol{\gamma})=\psi(\boldsymbol{\gamma})+2 b_{0} \lambda-\left(N+2 a_{0}-2\right) \ln \lambda+\\\sum\limits_{i=1}^{K} \ln \left(1+\varepsilon \lambda \gamma_{i}\right)+S(\lambda, \boldsymbol{\gamma}, \boldsymbol{\alpha}) . \end{array}$
其中ψ(γ)=-2lnp(γ)。将式(50)代入式(51),即得到可分离的联合成本函数,则有
$\begin{gather*}\mathcal{L}\left(\boldsymbol{x}, \boldsymbol{\alpha}, \sigma^{2}, \boldsymbol{\gamma}\right)= \\\sum\limits_{i=1}^{K} \ln \left(\sigma^{2}+\varepsilon \gamma_{i}\right)+\psi(\boldsymbol{\gamma})-n_{1} \ln \sigma^{2}+ \\\frac{R(\boldsymbol{x}, \boldsymbol{\alpha})+2 b_{0}}{\sigma^{2}}+\boldsymbol{x}^{\mathrm{T}} \boldsymbol{\varLambda} \boldsymbol{x} . \end{gather*}$
其中:$ n_{1}=K+2-N-2 a_{0}, R(\boldsymbol{x}, \boldsymbol{\alpha})=\|\boldsymbol{y}-\boldsymbol{H} \boldsymbol{\alpha}\|_{2}^{2} $+$ 2(\boldsymbol{x}-\boldsymbol{\alpha})^{\mathrm{T}} \boldsymbol{H}^{\mathrm{T}}(\boldsymbol{H} \boldsymbol{\alpha}-\boldsymbol{y})+\varepsilon\|\boldsymbol{x}-\boldsymbol{\alpha}\|_{2}^{2}$
对式(10)服从指数分布的p(γ)取对数得到ψ(γ),再将ψ(γ)代入式(51)得到Laplace先验的联合成本函数
$\begin{gather*}\mathcal{L}\left(\boldsymbol{x}, \boldsymbol{\alpha}, \sigma^{2}, \boldsymbol{\gamma}\right)=\sum\limits_{i=1}^{K}\left(\ln \left(\sigma^{2}+\varepsilon \gamma_{i}\right)+\varepsilon \gamma_{i}+\frac{x_{i}^{2}}{\gamma_{i}}\right)+ \\\frac{R(\boldsymbol{x}, \boldsymbol{\alpha})+2 b_{0}}{\sigma^{2}}-n_{1} \ln \sigma^{2}. \end{gather*}$
由于函数ψ(γ)具有下界,可以分解为凹函数项f1(γ)=0和凸函数项f2(γ)=$ \sum\limits_{i=1}^{K} \varepsilon \gamma_{i}$之和的形式。联合成本函数同样可分解为凹函数项f和凸函数项g之和的形式,具体式子为
$f(\boldsymbol{u})=f_{1}(\boldsymbol{\gamma})+\sum\limits_{i=1}^{K}\left(\ln \left(\sigma^{2}+\varepsilon \gamma_{i}\right)+\varepsilon \gamma_{i}+\frac{x_{i}^{2}}{\gamma_{i}}\right), $
$g(\boldsymbol{x}, \boldsymbol{\alpha}, \boldsymbol{u})=f_{2}(\boldsymbol{\gamma})+\frac{R(\boldsymbol{x}, \boldsymbol{\alpha})+2 b_{0}}{\sigma^{2}}-n_{1} \ln \sigma^{2}.$
其中u =(σ2, γ)。
为了解决联合成本函数产生的凹凸函数项问题,本文采用MM算法[29],并得到代理函数$ \hat{\mathcal{L}}_{k}$,具体为
$\hat{\mathcal{L}}_{k}(\boldsymbol{x}, \boldsymbol{\alpha}, \boldsymbol{u})=h\left(\boldsymbol{u}, \boldsymbol{u}^{(k)}\right)+g(\boldsymbol{x}, \boldsymbol{\alpha}, \boldsymbol{u}), $
$\begin{equation*}h\left(\boldsymbol{u}, \boldsymbol{u}^{(k)}\right)=\left\langle\boldsymbol{u}-\boldsymbol{u}^{(k)}, \nabla f\left(\boldsymbol{u}^{(k)}\right)\right\rangle+f\left(\boldsymbol{u}^{(k)}\right) . \end{equation*}$
其中h(u, u(k))表示在点u(k)处线性优化f(u)得到的控制函数。
由于代理函数具有部分可分解性,因此使用BCD算法进行分步求解超参数,具体步骤为
$\begin{equation*}\boldsymbol{x}^{(k+1)} \in \operatorname{argmin}_{{\boldsymbol{x}}} \hat{\mathcal{L}}_{k}\left(\boldsymbol{x}, \boldsymbol{\alpha}^{(k)}, \boldsymbol{u}^{(k)}\right), \end{equation*}$
$ \boldsymbol{\alpha}^{(k+1)} \in \operatorname{argmin}_{\boldsymbol{\alpha}} \hat{\mathcal{L}}_{k}\left(\boldsymbol{x}^{(k+1)}, \boldsymbol{\alpha}, \boldsymbol{u}^{(k)}\right), $
$\boldsymbol{u}^{(k+1)} \in \operatorname{argmin}_{u} \hat{\mathcal{L}}_{k}\left(\boldsymbol{x}^{(k+1)}, \boldsymbol{\alpha}{ }^{(k+1)}, \boldsymbol{u}\right) . $
上述的目标函数在零梯度点处达到最小值。因此,通过将梯度调节为零,得到更新规则为
$\hat{\boldsymbol{\varSigma}}=\left(\boldsymbol{\varLambda}^{-1}+\frac{\varepsilon}{\left(\sigma^{2}\right)} \boldsymbol{I}\right)^{-1} , $
$\hat{\boldsymbol{x}}=\frac{1}{\left(\sigma^{2}\right)} \boldsymbol{\varSigma}\left(\varepsilon \boldsymbol{\alpha}-\boldsymbol{H}^{T} \boldsymbol{H} \boldsymbol{\alpha}+\boldsymbol{H}^{T} \boldsymbol{y}\right) , $
$\hat{\boldsymbol{\alpha}}=\hat{\boldsymbol{x}} , $
$\hat{\gamma}_{i}=\sqrt{\frac{\left(\sigma^{2}\right)+\varepsilon \gamma_{i}}{\varepsilon+\beta\left(\sigma^{2}\right)+\beta \varepsilon \gamma_{i}}}\left|x_{i}\right|, $
$\rho=\sum\limits_{i=1}^{n} \frac{1}{\left(\sigma^{2}\right)+\varepsilon \gamma_{i}} , $
$\left(\hat{\sigma}^{2}\right)=\frac{n_{1}+\sqrt{n_{1}^{2}+4 \rho\left(\|\boldsymbol{y}-\boldsymbol{H} \boldsymbol{\alpha}\|_{2}^{2}+2 b_{0}\right)}}{2 \rho} $
ELSM-SBL-SMV算法总结如图 5所示。
值得注意的是,ELSM-SBL-SMV算法迭代式(61)中Σ的更新仅需要对角矩阵求逆的操作,因此其复杂度为$ {\cal{O}}$(K)。而式(62)中的HT是该算法的主要计算步骤,使得ELSM-SBL-SMV算法每次迭代的复杂度为$ {\cal{O}}$(K2)。而文[30]中GIG-SBL算法每次迭代需要对K×K维矩阵进行求逆,所需的复杂度为$ {\cal{O}}$(K3)。陈平平等[30]提出的基于Bernoulli Gauss逆Gamma(BGIG)先验的BGIG-SBL算法虽然通过引入Bernoulli先验的二元向量提高了算法的重构性能,但也将每次迭代的复杂度提升至$ {\cal{O}}$(K6)。算法1的LSM-SBL通过添加和删除基向量的方式来减少计算时间,其计算复杂度随着基函数数目的增多而增加,因此每次迭代的复杂度为$ {\cal{O}}$(Nv2K2),其中0<v<1。对于大规模数据集,所提出的ELSM-SBL-SMV算法在每次迭代中具有最低的计算复杂度。上述4种SBL算法每次迭代所需的计算复杂度如表 1所示。
表 1 4种SBL算法每次迭代的计算复杂度
SBL算法 每次迭代的计算复杂度
GIG-SBL[30] $ {\cal{O}}$(K3)
BGIG-SBL[30] $ {\cal{O}}$(K6)
LSM-SBL $ {\cal{O}}$(Nv2K2), 0<v<1
ELSM-SBL-SMV $ {\cal{O}}$(K2)

3.2 ELSM-SBL-MMV算法

在实际通信场景下,用户往往会在连续的几个相邻时隙内传输信息。因此,本节考虑能够利用时隙之间强相关性的多测量向量MMV场景进行建模。在MMV场景下,SBL算法通过将观测向量y扩展到T个接收时隙,能够同时对多个接收向量[y1, y2, …, yT]进行处理。
另外,在之前的mMTC传输场景中,基于SBL的多用户检测算法大多是假设信道增益g和扩频矩阵φ在多个时隙内保持不变,因此等效信道矩阵H在多个时隙内也是不变的。然而,实际mMTC场景的信道状态在多个时隙内是变化的。因此,本文考虑信道增益g在多时隙发送随时间变化,尽管扩频矩阵S依旧保持不变,但由于扩频矩阵S与不同的信道增益gnk(t)相乘,等效信道矩阵H被划分成了多个不同时隙[H1, H2, …, HT]。
Wipf等[31]通过共享Gamma分布中的公共超参数τ并利用Bayes推理生成了一种MSBL(multiple SBL,MSBL)算法。本文通过共享一组控制变量x稀疏性的超参数γ,将ELSM-SBL算法推广到上述的MMV模型中,提出了ELSM-SBL-MMV算法,进一步提高了稀疏信号的重构性能。在该模型中,ELSM-SBL-MMV算法仍假设各参数具有与上一节3.1相同的先验分布,即
$p\left(\lambda_{t}\right) \sim \operatorname{Gamma}\left(a_{0}, b_{0}\right), $
$p\left(\gamma_{i}\right)=\frac{\beta}{2} \exp \left(-\frac{\beta \gamma_{i}}{2}\right) . $
定义对于任意向量xt,有标量xi, txt, ∀i=1, 2, …, Kxi, t为第t个向量xi中的第i个元素,t=1, 2, …, T。基于MMV的ELSM-SBL-MMV算法的先验概率模型如图 6所示。
图 6 ELSM-SBL-MMV算法的概率图模型
按照3.1节的Bayes推理过程,同样利用构造Laplace先验的可分离联合成本函数$ {\cal{L}}_t$来估计MMV场景下的超参数,$ {\cal{L}}_i$可具体表示为
$\begin{gather*}\mathcal{L}_{t}\left(\boldsymbol{x}_{t}, \boldsymbol{\alpha}_{t}, \sigma_{t}^{2}, \boldsymbol{\gamma}_{t}\right)= \\\sum\limits_{i=1}^{K}\left(\ln \left(\sigma_{t}^{2}+\varepsilon_{t} \gamma_{i, t}\right)+\varepsilon_{t} \gamma_{i, t}+\frac{\boldsymbol{x}_{i, t}^{2}}{\gamma_{i, t}}\right)+ \\\frac{R\left(\boldsymbol{x}_{t}, \boldsymbol{\alpha}_{t}\right)+2 b_{0}}{\sigma_{t}^{2}}-n_{1} \ln \sigma_{t}^{2} .\end{gather*}$
其中εt=eig(HtTHt)+ε0, ε0为某一大于零的常数。
紧接着,利用MM算法得到代理函数并使用BCD算法分步求解超参数,得到变量和参数的更新规则为
$\hat{\boldsymbol{\varSigma}}_{t}=\left(\boldsymbol{\varLambda}_{t}^{-1}+\frac{\varepsilon_{t}}{\left(\sigma_{t}^{2}\right)} \boldsymbol{I}\right)^{-1}, $
$\hat{\boldsymbol{x}}_{t}=\frac{1}{\left(\sigma_{t}^{2}\right)} \boldsymbol{\varSigma}_{t}\left(\varepsilon_{t} \boldsymbol{\alpha}_{t}-\boldsymbol{H}_{t}^{\mathrm{T}} \boldsymbol{H}_{t} \boldsymbol{\alpha}_{t}+\boldsymbol{H}_{t}^{\mathrm{T}} \boldsymbol{y}_{t}\right), $
$\hat{\boldsymbol{\alpha}}_{t}=\hat{\boldsymbol{x}}_{t}, $
$\hat{\gamma}_{i}=\frac{\sum\limits_{t=1}^{T}\left(\left|\boldsymbol{x}_{i, t}\right|\right)}{\sum\limits_{t=1}^{T}\left(\sqrt{\beta+\frac{\varepsilon_{t}}{\sigma_{t}^{2}+\varepsilon_{t} \gamma_{i}}}\right)}, $
$\rho_{t}=\sum\limits_{i=1}^{n} \frac{1}{\left(\sigma_{t}^{2}\right)+\varepsilon_{t} \gamma_{i}}, $
$\left(\hat{\sigma}_{t}^{2}\right)=\frac{n_{1}+\sqrt{n_{1}^{2}+4 \rho_{t}\left(\left\|\boldsymbol{y}_{t}-\boldsymbol{H}_{t} \boldsymbol{\alpha}_{t}\right\|_{2}^{2}+2 b_{0}\right)}}{2 \rho_{t}} . $
从ELSM-SBL-MMV算法更新规则(见式(73))可以看出,算法的关键在于通过共享一组控制信号稀疏性的超参数γ =(γ1, γ2, …, γK)实现多时隙联合优化。具体而言,超参数γi的更新融合了T个时隙的信号xi, t、噪声方差σt2以及信道参数εt,从而增强多时隙信号的联合稀疏性。
ELSM-SBL-MMV算法总结如图 7所示。
图 7 ELSM-SBL-MMV算法
值得注意的是,ELSM-SBL-MMV算法通过共享超参数的方式提升信号的重构性能,主要体现在式(73)上,其复杂度为$ {\cal{O}}$(TK)。与SMV场景下类似,式(71)才是该算法的主要计算步骤,故ELSM-SBL-MMV算法的复杂度为$ {\cal{O}}$(TK2)。在MMV场景下,几种SBL算法每次迭代所需的计算复杂度如下表 2所示,可以看到所提算法具有最低的复杂度。
表 2 MMV场景下不同SBL算法每次迭代的计算复杂度
SBL算法 每次迭代的计算复杂度
GIG-SBL-MMV[30] $ {\cal{O}}$(TK3)
BGIG-SBL_MMV[30] $ {\cal{O}}$(TK6)
LSM-SBL_MMV $ {\cal{O}}$(TK2v2K2), 0<v<1
ELSM-SBL-MMV $ {\cal{O}}$(TK2)

4 算法仿真结果及分析

在仿真实验参数设置方面,本文设置多用户mMTC通信中的总活跃用户数K=108,子载波数N=72,活跃用户数M=12,采用BPSK调制,信道增益系数gnk(t)~N(0, 1),重复试验次数为1 000次。对于MMV模型,本文选择测量向量的数量T =7。采用GIG-SBL、BGIG-SBL以及LSM-SBL算法作为本实验的主要对比算法[30]。在所有实验中,算法的最大迭代次数设置为500次,初始化参数设置为a0=b0=1×10-6β=0及所有γi=0,迭代停止阈值设置为1×10-5[26]
本文通过比较不同算法分别在SMV和MMV场景中的检测性能,所采用的检测性能指标包括误比特率(bit error rate,BER)和活跃度差错率(activity error rate,AER)。
BER和AER可以分别表示为
$\mathrm{BER}=\frac{1}{N_{s}} \sum\limits_{n=1}^{N_{s}}\frac{{{\rm{be}}{{\rm{r}}_n}}}{{KT}} , $
$\mathrm{AER}=\frac{1}{N_{s}} \sum\limits_{n=1}^{N_{s}} \frac{{{\rm{aer}}{_n}}}{{MT}}.$
其中:Ns表示仿真的实验次数,bern表示第n次重构信号后的错误比特数,aern表示第n次估计的活跃用户索引出错个数。

4.1 SMV传输场景下的性能仿真分析

在SMV传输场景下,随着信噪比SNR的增加, 不同SBL算法的BER表现如图 8a所示。图中可以得知,随着SNR的增加,各个算法的重构性能都会有所提升,本文所提出的ELSM-SBL-SMV相较于改进前的LSM-SBL-SMV约有2 dB的提升。同时,与传统的GIG-SBL-SMV和改进BGIG-SBL-SMV算法相比,分别实现了约4 dB和2 dB的性能提升。这是因为,相比于BGIG-SBL-SMV,本文所提的ELSM-SBL-SMV利用了Laplace分布具有重尾的特性,其概率密度函数在零点附近有一个尖锐的峰值,这使得它能够更好地促进信号的稀疏性,从而提高稀疏信号的重构性能。此外,本文提出的ELSM-SBL-SMV算法通过构造惩罚项可分离的联合成本函数,并结合BCD方法和MM策略解决了其产生的非凸优化问题,相较于其他SBL算法,表现出了更优的信号重构性能。从图 8b可以进一步看出,在SMV传输模型下,本文提出的ELSM-SBL-SMV在检测用户活跃索引方面,均明显优于其他算法,验证了所提算法的信号估计和重构性能优势。
图 8 SMV场景下不同SBL算法的BER和AER对比

4.2 MMV传输场景下的性能仿真分析

在MMV传输场景下,不同SBL算法在不同信噪比SNR条件下的误比特率BER如图 9a所示。从图中可以观察到,在MMV构架下,当BER达到10-5时,本文提出的ELSM-SBL-MMV算法相较于最优的BGIG-SBL-MMV约有3 dB的提升,与传统的GIG-SBL-MMV和改进前的LSM-SBL-MMV相比,分别实现了约5 dB和6 dB的性能增益。这主要是因为ELSM-SBL-MMV在处理多时隙信号时,通过共享控制模型稀疏性的超参数,利用多时隙信号向量之间的时间相关性,提高了信号重构的性能。
图 9 MMV场景下不同SBL算法的BER和AER对比
图 9b展示了不同SBL算法在MMV场景下不同SNR的AER。可以看出,在相同AER水平时,本文提出ELSM-SBL-MMV在检测用户方面相较于传统的GIG-SBL-MMV和LSM-SBL-MMV法均实现了约2 dB的性能提升。这是因为ELSM-SBL-MMV算法通过共享控制模型稀疏性的超参数将多时隙联合优化,提高了稀疏信号的时间相关性,提升了多时隙检测性能。此外,在MMV场景下,用户活跃状态在时隙间高度相关,SBL算法能够通过能量累积排序检测,提升了AER性能,由于BGIG-SBL-MMV与ELSM-SBL-MMV采用了类似的共享超参数方案,因此两者算法提升后的性能差距较小。

4.3 活跃用户度对SBL算法的性能影响

该小节探究活跃用户占比对重构信号的性能影响,假定活跃用户度Pa=M/K。在SMV传输场景下,SNR=3 dB时,不同SBL算法在不同Pa下的BER性能如图 10a所示。由图可知,本文所提的ELSM-SBL-SMV算法在Pa较小(Pa<0.15)时,展现出了良好的检测性能。例如,在Pa=0.05时,本文提出的ELSM-SBL-SMV算法相较于LSM-SBL-SMV算法和BGIG-SBL-SMV算法具有更优的重构性能。然而,随着Pa的增加,当Pa≥0.15时,几种SBL算法的重构性能均缓慢降低并趋于一致。这表明,SBL算法恢复效果依赖于信号的稀疏性质,当用户总数保持不变,而活跃用户数量增加时,信号的稀疏性降低,从而导致信号重构性能下降。本文提出的ELSM-SBL-SMV算法对于相同的Pa,相较于其他算法仍能保持较高的重构性能水平。从图 10b可以看出,在MMV传输场景下,SNR=8 dB时,在不同Pa条件下,本文提出的ELSM-SBL-MMV算法的BER性能均明显优于其他算法,验证了所提算法的优势和稳定性。
图 10 不同活跃用户度下的BER仿真结果

4.4 算法复杂度的仿真分析

为分析ELSM-SBL算法的计算复杂度,假定Pa保持0.1不变,过载率ρ=K/N=100%不变,用户总数K从10逐渐增加到100,观察重构信号每次需要的平均运行时间[32]。实验在16GB RAM与Intel(R) Core(TM) i7-10700 CPU的MATLAB R2021b上运行实现。从图 11a中可以看出,在SMV场景下,随着用户总数的增加,SBL算法每次重构信号所需的时间均有所增加,但本文所提ELSM-SBL-SMV算法在用户总数相同的情况下,相较于其他算法的运行时间有明显优势。这是因为ELSM-SBL-SMV算法通过引入代理函数近似模型参数的后验分布,从而规避矩阵的求逆操作,使得算法的复杂度降低至$ {\cal{O}}$(K2),验证了所提算法在计算复杂度上的优势。图 11b展示了MMV场景下的不同算法运行时间,可以看出本文所提的ELSM-SBL-MMV相较于其他算法的计算时间优势。
图 11 不同用户总数下不同SBL算法的运行时间

5 结论

本文针对传统稀疏Bayes学习算法复杂度高的问题,提出一种基于高效Laplace尺度混合先验的SBL算法。ELSM-SBL算法通过引入代理函数近似模型参数的后验分布,设计构造一个惩罚项可分离的联合成本函数,并利用BCD方法解决由此产生的非凸优化问题,从而规避了矩阵求逆的操作。主要实验结果如下:
1) 在mMTC多用户检测场景,ELSM-SBL算法在SMV场景下,在AER和BER核心指标上都显著优于现有的LSM-SBL算法、BGIG-SBL算法和GIG-SBL算法。
2) ELSM-SBL算法通过共享控制稀疏解的超参数,利用相邻向量之间的联合稀疏性,提升了算法在MMV传输场景下的信号重构性能。
3) 在算法复杂度方面,ELSM-SBL能够实现降低至$ {\cal{O}}$(K2)级别,重构信号每次所需要的平均运行时间低于其他算法。
1
MAHMOOD N H , BÖCKER S , MOERMAN I , et al. Machine type communications: Key drivers and enablers towards the 6G era[J]. EURASIP Journal on Wireless Communications and Networking, 2021, 2021 (1): 134.

DOI

2
EJAZ W , ANPALAGAN A , IMRAN M A , et al. Internet of Things (IoT) in 5G wireless communications[J]. IEEE Access, 2016, 4, 10310- 10314.

DOI

3
石昌伟, 郭里婷, 康芃, 等. 可学习阈值优化的大规模动态多用户接入检测[J]. 电子学报, 2025, 53 (5): 1436- 1444.

SHI C W , GUO L T , KANG P , et al. Learnable threshold optimization for massive dynamic multi-user access detection[J]. Acta Electronica Sinica, 2025, 53 (5): 1436- 1444.

4
王平, 孙臻, 殷柳国, 等. 用于量子安全直接通信的空间耦合LDPC-BCH码[J]. 清华大学学报(自然科学版), 2019, 59 (9): 737- 743.

WANG P , SUN Z , YIN L G , et al. Spatially coupled LDPC-BCH codes in quantum secure direct communications[J]. Journal of Tsinghua University (Science and Technology), 2019, 59 (9): 737- 743.

5
崔兆阳, 黄容兰, 万德焕. 一种有效改善系统性能的非正交传输方案[J]. 电子学报, 2020, 48 (10): 1915- 1922.

CUI Z Y , HUANG R L , WAN D H . A promising non-orthogonal multiple access scheme: Devoting to system performance improvement[J]. Acta Electronica Sinica, 2020, 48 (10): 1915- 1922.

6
DONOHO D L , TSAIG Y , DRORI I , et al. Sparse solution of underdetermined systems of linear equations by stagewise orthogonal matching pursuit[J]. IEEE Transactions on Information Theory, 2012, 58 (2): 1094- 1121.

DOI

7
丁青锋, 李怡浩, 徐梦引. 去蜂窝大规模MIMO-NOMA系统能效优化算法[J]. 电子学报, 2023, 51 (8): 2020- 2029.

DING Q F , LI Y H , XU M Y . Energy efficiency optimization algorithm for cell-free massive MIMO-NOMA Systems[J]. Acta Electronica Sinica, 2023, 51 (8): 2020- 2029.

8
陈平平, 张旭, 谢肇鹏, 等. 基于多智能体近端策略优化的多信道动态频谱接入[J]. 电子学报, 2024, 52 (6): 1824- 1831.

CHEN P P , ZHANG X , XIE Z P , et al. Multi-channel dynamic spectrum access based on multi-agent proximal policy optimization[J]. Acta Electronica Sinica, 2024, 52 (6): 1824- 1831.

9
HONG J P , CHOI W , RAO B D . Sparsity controlled random multiple access with compressed sensing[J]. IEEE Transactions on Wireless Communications, 2015, 14 (2): 998- 1010.

DOI

10
石光明, 刘丹华, 高大化, 等. 压缩感知理论及其研究进展[J]. 电子学报, 2009, 37 (5): 1070- 1081.

SHI G M , LIU D H , GAO D H , et al. Advances in theory and application of compressed sensing[J]. Acta Electronica Sinica, 2009, 37 (5): 1070- 1081.

11
SHIM B , SONG B . Multiuser detection via compressive sensing[J]. IEEE Communications Letters, 2012, 16 (7): 972- 974.

DOI

12
MONSEES F, WOLTERING M, BOCKELMANN C, et al. Compressive sensing multi-user detection for multicarrier systems in sporadic machine type communication[C]//2015 IEEE 81st Vehicular Technology Conference (VTC Spring). Glasgow, UK: IEEE, 2015: 1-5.

13
WANG B C, DAI L L, YUAN Y F, et al. Compressive sensing based multi-user detection for uplink grant-free non-orthogonal multiple access[C]//2015 IEEE 82nd Vehicular Technology Conference (VTC2015-Fall). Boston, USA: IEEE, 2015: 1-5.

14
DAVENPORT M A , NEEDELL D , WAKIN M B . Signal space CoSaMP for sparse recovery with redundant dictionaries[J]. IEEE Transactions on Information Theory, 2013, 59 (10): 6820- 6829.

DOI

15
TIPPING M E . Sparse Bayesian learning and the relevance vector machine[J]. The Journal of Machine Learning Research, 2001, 1, 211- 244.

16
ZHANG X X , FAN P Z , LIU J Q , et al. Bayesian learning-based multiuser detection for grant-free NOMA systems[J]. IEEE Transactions on Wireless Communications, 2022, 21 (8): 6317- 6328.

DOI

17
WIPF D P , RAO B D . Sparse Bayesian learning for basis selection[J]. IEEE Transactions on Signal Processing, 2004, 52 (8): 2153- 2164.

DOI

18
TIPPING M E, FAUL A C. Fast marginal likelihood maximisation for sparse Bayesian models[C]// Proceedings of the Ninth International Workshop on Artificial Intelligence and Statistics. PMLR, 2003: 276-283.

19
XIAO L , XU D J , MANDAYAM N B , et al. Attacker-centric view of a detection game against advanced persistent threats[J]. IEEE Transactions on Mobile Computing, 2018, 17 (11): 2512- 2523.

DOI

20
ZHANG Z L , RAO B D . Sparse signal recovery with temporally correlated source vectors using sparse Bayesian learning[J]. IEEE Journal of Selected Topics in Signal Processing, 2011, 5 (5): 912- 926.

DOI

21
ABEBE A T , KANG C G . Iterative order recursive least square estimation for exploiting frame-wise sparsity in compressive sensing-based MTC[J]. IEEE Communications Letters, 2016, 20 (5): 1018- 1021.

DOI

22
ZHANG Y Y , GUO Q H , WANG Z Y , et al. Block sparse Bayesian learning based joint user activity detection and channel estimation for grant-free NOMA systems[J]. IEEE Transactions on Vehicular Technology, 2018, 67 (10): 9631- 9640.

DOI

23
AL-SHOUKAIRI M , SCHNITER P , RAO B D . A GAMP-based low complexity sparse Bayesian learning algorithm[J]. IEEE Transactions on Signal Processing, 2018, 66 (2): 294- 308.

DOI

24
TAN X , LI J . Computationally efficient sparse Bayesian learning via belief propagation[J]. IEEE Transactions on Signal Processing, 2010, 58 (4): 2010- 2021.

DOI

25
SHEKARAMIZ M, MOON T K. Compressive sensing via variational Bayesian inference[C]//2020 Intermountain Engineering, Technology and Computing (IETC). Orem, USA: IEEE, 2020: 1-6.

26
BABACAN S D , MOLINA R , KATSAGGELOS A K . Bayesian compressive sensing using Laplace priors[J]. IEEE Transactions on Image Processing, 2010, 19 (1): 53- 63.

DOI

27
高鹏宇. 压缩感知辅助的非正交多址接入检测技术研究[D]. 成都: 电子科技大学, 2019.

GAO P Y. Research on non-orthogonal multiple access detection technology assisted by compressed sensing[D]. Chengdu: University of Electronic Science and Technology of China, 2019. (in Chinese)

28
张贤达. 矩阵分析与应用[M]. 北京: 清华大学出版社, 2004.

ZHANG X D . Matrix analysis and applications[M]. Beijing: Tsinghua University Press, 2004.

29
ZHOU W , ZHANG H T , WANG J . An efficient sparse Bayesian learning algorithm based on Gaussian-scale mixtures[J]. IEEE Transactions on Neural Networks and Learning Systems, 2022, 33 (7): 3065- 3078.

DOI

30
陈平平, 王宣达, 谢肇鹏, 等. 基于稀疏贝叶斯学习的大规模多用户检测算法[J]. 通信学报, 2023, 44 (10): 186- 197.

CHEN P P , WANG X D , XIE Z P , et al. Sparse Bayesian learning-based massive multi-user detection algorithm[J]. Journal on Communications, 2023, 44 (10): 186- 197.

31
WIPF D P , RAO B D . An empirical Bayesian strategy for solving the simultaneous sparse approximation problem[J]. IEEE Transactions on Signal Processing, 2007, 55 (7): 3704- 3716.

DOI

32
HASAN S M , MAHATA K , HYDER M M . Uplink grant-free NOMA with sinusoidal spreading sequences[J]. IEEE Transactions on Communications, 2021, 69 (6): 3757- 3770.

DOI

Outlines

/