广义近似消息传递（GAMP-MMSE-SBL）

1、为了解决什么问题

[1]Rangan, Sundeep. "Generalized approximate message passing for estimation with random linear mixing." 2011 IEEE International Symposium on Information Theory Proceedings. IEEE, 2011.

[2]Andersen, Michael Riis. "Sparse inference using approximate message passing." Technical University of Denmark, Department of Applied Mathematics and Computing (2014).

上一篇介绍的 AMP 算法提供了一种有效解决 BP 或 BPDN 问题的方法。 Sundeep Rangan 表明，该框架可以推广到处理本质上任意的先验分布和任意的噪声分布。对两组分布的唯一要求是它们进行因式分解。论文[1]中介绍了这个通用框架，即广义近似消息传递 (GAMP)。

GAMP 的灵活性使我们能够使用稀疏性进行有效的推理，从而促进先验分布，例如尖峰和平板先验 (其中一个典型的spike and slab 先验就是伯努利高斯（BG）先验)。此外，由于噪声分布也是任意的，因此该框架还可以通过使用二项式噪声分布进行分类。由于通常用到的噪声分布一般是加性高斯白噪声。这里我们也基于该模型。

尽管模型的灵活性大大提高，但计算复杂度保持不变，即 O(nm)。我们在GAMP中用更新的消息只有全链接两侧消息及噪声方差。当然有些算法不需要更新噪声方差，噪声方差作为一个已知项。GAMP 框架既可以配置为基于最大和消息传递执行 MAP 估计，也可以配置为基于和积消息传递执行 MMSE 估计。GAMP 框架的推导比 AMP 的推导更直接一些。推导主要基于泰勒近似和中心极限定理的应用。而在国内的郑大王教授团队，张老师采用因子图上利用BP-MF-SBL及三步近似得到了相同的结果。我认为张老师的推导是出于规则的近似，更像是用利用BP的规则工具导出，而利用泰勒近似和中心极限定理这个推导（数学不好，估计得看很久）。虽然他们都是处于对和积算法的近似。AMP算法可以说是GAMP算法的一个特例，当稀疏先验变成拉普拉斯先验噪声分布选择高斯分布，则由论文[1]中的第二个GAMP算法可以推到出。详细的推到可以在论文[2]中查找。事实上，（更多详细信息请参阅附录 B.2）。

2、算法来源[1,2]

[1]Rangan, Sundeep. "Generalized approximate message passing for estimation with random linear mixing." 2011 IEEE International Symposium on Information Theory Proceedings. IEEE, 2011.

3、先验模型选取与gout和gin的关系

先验模型的选取就有很多选项了，这里我们可以根据系统中的未知参数特性选择合适的先验。比如未知参数是离散的调制符号，此刻系统不再是稀疏问题。就不需要稀疏先验来解决。此时利用GAMP框架和合适的gin函数，可以得到响应的GAMP算法。如果是spike and slab先验或者是sbl先验此刻处理的就是稀疏问题。而gout是

通用框架下可以带这两个参数来更新s和tau_s。

gin函数就是与变量的先验有关了（当然在GAMP-MMSE框架下）

4、具体算法

[1]Luo, Man, et al. "Unitary approximate message passing for sparse Bayesian learning." IEEE transactions on signal processing 69 (2021): 6023-6039.

当然经过酉处理后，可以套用酉近似消息传递算法由论文[1]。不经过酉处理，就是正常的GAMP-SBL。如果用到如下代码，觉得好用还请引用该文献。该文献可以帮助你快速学习酉近似消息传递算法。关于先验的推导文中也有。

5、代码实现（SBL先验）

function [Xhat,Xhatvar,nmseGAMP,gamma_l] = GAMP_vector_SBL(A,y,noise_lamda1,x)
% y = A*x + w=z+w ,where A is a random matrix and w is Gaussian noise.x is a sparse vector
% 未知噪声方差,noise_lamda1为噪声方差初始值
[N, L]=size(A);
noise_lamda=ones(N,1)*noise_lamda1;
%% Define variables  
Xhatvar= ones(L,1); Xhat = zeros(L,1);    % Define variable x (belief) E and Var                
epc = 0;                                  % Initial hyperparameter
eta = 0*ones(L,1);                        % Initial hyperparameter
gamma_l = ones(L,1);                      % Initial hyperparameter
Shat = zeros(N,1);                        % Initial s_n N*1 all zeros
%% Algorithm body
for t = 1:1:30
    %t    
      %% update GAMP left out message phat and phatvar       
         phatvar=(abs(A).^2)*Xhatvar;        
         phat=A*Xhat - phatvar.*Shat;
        %% update intermediate variables Shat 
        Shatvar=1./(noise_lamda+phatvar);
        Shat=(y-phat).*Shatvar;  
         %% update variables z(belief) E and Var
        [Zhat,Zhatvar] =GaussianProduct(y,noise_lamda,phat,phatvar); 
        
        %% update GAMP right out message rhat and  rhatvar     
         rhatvar= 1./(((abs(A).').^2)*Shatvar);
         rhat=rhatvar.*((A')*Shat)+Xhat;%  
         
       %% update x belief: Xhat and Xhatvar
        %[Xhat,Xhatvar] =GaussianProduct(0,gamma_l.^-1,rhat,rhatvar);
        Xhatvar= rhatvar./(1+gamma_l.*rhatvar);
        Xhat  = rhat ./(1+gamma_l.*rhatvar);        
         %% update gamma_l
        gamma_l = (1+epc)./(eta+abs(Xhat).^2+Xhatvar);      

        %% update x belief: Xhat and Xhatvar
         %[Xhat,Xhatvar] =GaussianProduct(0,gamma_l.^-1,rhat,rhatvar);
         Xhatvar= rhatvar./(1+gamma_l.*rhatvar);
         Xhat  = rhat ./(1+gamma_l.*rhatvar);

         nmseGAMP(t) = 20*log10( norm(x-Xhat)/norm(x));
       %% update Noise Accuracy (belief)  
        lamda = N/(sum(Zhatvar)+(y-Zhat)'*(y-Zhat));
        noise_lamda=ones(N,1)/lamda;      
end

function [m3,v3] =GaussianProduct(m1,v1,m2,v2)
v3=(1./(v1)+1./(v2)).^-1;
m3=(m1./v1+m2./v2).*v3;
end
%Msg_RightA_v=(1./AvarRep-1./Msg_LeftA_v).^-1;
%Msg_RightA_m=(AhatRep./AvarRep-Msg_LeftA_m./Msg_LeftA_v).*Msg_RightA_v;

6、实际效果

相对于AMP算法可谓是更进一步，此外AMP算法很不稳定，还需要设置一个正则化项，这个项的设置恰到与否影响算法是否工作

AMP算法（近似消息传递算法）的发散性主要受到与i.i.d.（独立同分布）亚高斯模型的微小偏差的影响。AMP算法是基于一种计算效率很高的迭代方法，对于大I.I.D.亚高斯矩阵A，其每次迭代行为的严格特征是标量状态演化，其不动点在唯一时是Bayes最优的。然而，由于AMP算法的脆弱性，即使是与i.i.d.亚高斯模型的微小偏差也可能导致算法发散。

因此，在使用AMP算法时，需要确保数据满足其假设条件，即数据应接近i.i.d.亚高斯分布。如果数据不满足这些条件，可能需要考虑使用其他算法或对数据进行预处理以满足AMP算法的假设条件。

以上信息仅供参考，如有需要，建议查阅相关文献或咨询专业人士。