IEEE Trans 2008 Gradient Pursuits论文学习

之前所学习的论文中求解稀疏解的时候一般采用的都是最小二乘方法进行计算,为了降低计算复杂度和减少内存,这篇论文梯度追踪,属于贪婪算法中一种。主要为三种:梯度(gradient)、共轭梯度(conjugate gradient)、近似共轭梯度(an approximation to the conjugate gradient),看师兄之前做压缩感知的更新点就是使用近似共轭梯度方法代替了StOMP中的最小二乘的步骤。

首先说明一下论文中的符号表示:

则GP算法流程如下:

在作者的主页上下载了论文中的代码,(附下载地址为:http://www.personal.soton.ac.uk/tb1m08/sparsify/sparsify.html

对其中的GP算法进行了自行仿写,但是不知道为什么跟作者的代码比较起来效果很差,代码所在文件夹为sparsify_0_5,将自己写的代码贴出如下,用了两种不同的写法,其实原理都是一样的,仿真结果很差,残差很大。

function[theta]=CS_GP(y,A,t)
    %CS_GP Summary of this function goes here
    %Version: 1.0 written by wwf @2017-04-28
    %   Detailed explanation goes here
    %   y = Phi * x
    %   x = Psi * theta
    %   y = Phi*Psi * theta
    %   令 A = Phi*Psi, 则y=A*theta
    %   现在已知y和A,求theta
    [y_rows,y_columns]= size(y);
    if y_rows<y_columns
        y = y';%y should be a column vector
    end
    [M,N]= size(A);                            %传感矩阵A为M*N矩阵
    theta = zeros(N,1);                         %用来存储恢复的theta(列向量)
    aug_y=[];
    r_n = y;                                    %初始化残差(residual)为y
    Aug_t=[];
    for ii =1:t                                  %迭代t次,t为输入参数
        for col=1:N;
            product(col)=abs(A(:,col)'*r_n);             %传感矩阵A各列与残差的内积
        end
        [val,pos]= max(abs(product));        %找到最大内积绝对值,即与残差最相关的列
        Aug_t=[Aug_t,A(:,pos)]; 
        pos_array(ii)=pos;                             
        g_n=Aug_t'*r_n;                               %  梯度方向
        c_n=Aug_t*g_n;
        w_n=(r_n'*c_n)/(c_n'*c_n);                   %  最速下降步长
        d_n=w_n*g_n;
        [x1,x2]=size(d_n);
        [y1,y2]=size(aug_y);
        D=aug_y;
        aug_y=zeros(x1,x2);
        aug_y(1:y1,1:y2)=D;
        aug_y=aug_y+d_n ;                                %  最小二乘,使残差最小
        r_n=r_n-(w_n)*(c_n);                               %  残差
    end
   theta(pos_array)=aug_y;
end

第二种:

function[theta]=GP_test(y,A,t)
%CS_GP Summary of this function goes here
%Version: 1.0 written by wwf @2017-04-28
%   Detailed explanation goes here
%   y = Phi * x
%   x = Psi * theta
%   y = Phi*Psi * theta
%   令 A = Phi*Psi, 则y=A*theta
%   现在已知y和A,求theta
    [y_rows,y_columns]= size(y);
    if y_rows<y_columns
        y = y';%y should be a column vector
    end
    [M,N]= size(A);                            %传感矩阵A为M*N矩阵
    theta = zeros(N,1);                         %用来存储恢复的theta(列向量)
    r_n = y;                                    %初始化残差(residual)为y
    d_n=zeros(N,1);
    P =@(z) A*z;
    Pt =@(z) A'*z;
    IN=[];
   
    for ii =1:50                                  %迭代t次,t为输入参数
        product=Pt(r_n);
        [v I]=max(abs(product));
     
        if isempty(find (IN==I))
            IN=[IN I];
        else
            break;
        end
        d_n(IN)=product(IN);
        c_n=P(d_n);
        a_n=r_n'*c_n/(c_n'*c_n);
        theta=theta+a_n*d_n;
        r_n=r_n-a_n*c_n;
    end
end

测试代码和前几次的一样,实验结果如下所示:

输出残差为:

调用作者写的代码的时候发现,有的时候恢复效果比较好,残差很小,但有的时候也会出现残差比较大的情况,猜测可能和生成的信号有关系,因为每次信号是随机生成的。结果如下图所示:

在参考文献[2]中对论文进行了翻译,读过一遍之后对文献[1]也更加理解了,这里就不再单独给出,摘取文献[2]中的内容进行解释,在以后的论文编写中若需要文献[1],建议自己再翻译比较好。

接下来解释一下CGP和ACGP的主要代码,CGP的代码没有看太懂。

CGP(代码名称为greed_omp_cgp)

while~done
     DR(IN)=0;
    [v I]=max(abs(DR));
     IN=[IN I];
     k=k+1;
        ifk==1
             d(IN)=1;
             PG(1)=1;
        else
   %%%% Calculate P'*G, but only need new column and new row %%%%%
             mask=zeros(m,1);
             mask(IN(k))=1;%将mask中对应的第k次迭代所选出的内积所在的列序号的项置为1
             new_element=P(mask);%选出此次迭代所选择出的原子
             gnew=Pt(new_element);%gnew相当于G,即Psi'*Psi
             PG(k-1,k-1)=D(1:k,k-1)'*[g;1];
             g=gnew(IN);
            %PG计算的是D’*G
             PG(:,k)=D'*[g;zeros(maxM-k,1)]; % 1 general mult.
   %%%% Calculate conjugate directions %%%
             b=(PG(1:k-1,1:k)*DR(IN))./(dPPd(1:k-1));
             d(IN)=DR(IN)-D(1:k,1:k-1)*b;%d should be orthogonal to the first k-1 columns of G.
        end
         D(1:k,k)=d(IN);%D是由n-1次的更新方向组成的矩阵
         Pd=P(d);
         dPPd(k)=Pd'*Pd;
         a=(DR'*d)/dPPd(k);
         s=s+a*d;
         Residual=Residual-a*Pd;
         DR=Pt(Residual);
    
    ERR=Residual'*Residual/n;
    ifcomp_err
         err_mse(iter)=ERR;
    end
    
    ifcomp_time
         iter_time(iter)=toc;
    end

ACGP(代码名称为greed_nomp)

tic
t=0;
p=zeros(m,1);
DR=Pt(Residual);
[v I]=max(abs(DR));
ifweakness~=1
   [vals inds]=sort(abs(DR),'descend');
    I=inds(find(vals>=alpha*v));
end
   
IN=union(IN,I);
ifstrcmp(STOPCRIT,'M')&length(IN)>=STOPTOL
    IN=IN(1:STOPTOL);
end
MASK=zeros(size(DR));
pDDp=1;
done=0;
iter=1;
while~done
   
   % Select new element
   ifisa(GradSteps,'char')
       ifstrcmp(GradSteps,'auto')
            
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  Iteration to automatic selection of the number of gradient steps
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%             finished=0;   
%             while ~finished
           % Update direction   
                ifiter==1
                     p(IN)=DR(IN);   %p相当于论文中的d_n,当迭代次数为1时,d_n等于内积
                     Dp=P(p);        %Dp相当于论文中的c_n,即Psi与d_n的乘积
                else
                     MASK(IN)=1;     %IN为此次迭代选出的内积值最大的列序号,将MASK的该项置为1
                     PDR=P(DR.*MASK);%取出最大的内积值,与字典矩阵Psi相乘
                     b=-Dp'*PDR/pDDp;%计算系数b1
                     p(IN)=DR(IN)+b*p(IN);%计算更新的方向d_n
                     Dp=PDR+b*Dp;  %c_n是Psi与d_n-1的乘积,即P(d_n-1),将d_n-1展开带入即得
                end
            % Step size
%                  Dp=P(p); % =P(DR(IN)) +b P(p(IN));
                 pDDp=Dp'*Dp; 
                 a=Residual'*Dp/(pDDp);
            % Update coefficients  
                 s=s+a*p;
            % New Residual and inner products
                 Residual=Residual-a*Dp;
                 DR=Pt(Residual);
                % select new element
                    [v I]=max(abs(DR));
                    ifweakness~=1
                        [vals inds]=sort(abs(DR),'descend');
                         I=inds(find(vals>=alpha*v));
                    end
                     IN=union(IN,I);
                    ifstrcmp(STOPCRIT,'M')&length(IN)>=STOPTOL
                        IN=IN(1:STOPTOL);
                    end
%                  % Only if we select new element do we leave the loop   
%                      if isempty(find (IN==I, 1))
%                         IN=[IN I];
%                         finished=1;
%                      end
%             end

参考文献:

[1] Blumensath  T,  Davies  M.  Gradient  pursuits[J].  IEEE  Transactions  on  Signal  Processing,  2008, 56(6):2370-2382. 

[2] 刘盼盼.压缩感知中梯度追踪算法的研究[D].南京:南京邮电大学,2015:7-21

本文参与腾讯云自媒体分享计划,欢迎正在阅读的你也加入,一起分享。

发表于

我来说两句

0 条评论
登录 后参与评论

相关文章

来自专栏AI研习社

教你几招搞定 LSTMs 的独门绝技(附代码)

如果你用过 PyTorch 进行深度学习研究和实验的话,你可能经历过欣喜愉悦、能量爆棚的体验,甚至有点像是走在阳光下,感觉生活竟然如此美好 。但是直到你试着用 ...

6561
来自专栏数据派THU

为何RNN能够在众多机器学习方法中脱颖而出?(附指南)

来源:机器人圈 作者:BaymaxZ 本文长度为5000字,建议阅读20分钟 本文介绍RNN的重要性和先进性,并详细阐释几种用于深度学习中的RNN模型。 近年来...

2705
来自专栏深度学习自然语言处理

实例讲解决策树分类器

决策树是一种简单高效并且具有强解释性的模型,广泛应用于数据分析领域。其本质是一颗由多个判断节点组成的树,如:

1144
来自专栏UAI人工智能

连载 | 深度学习入门第六讲

1356
来自专栏机器之心

谷歌开放GNMT教程:如何使用TensorFlow构建自己的神经机器翻译系统

选自谷歌 机器之心编译 参与:机器之心编辑部 近日,谷歌官方在 Github 开放了一份神经机器翻译教程,该教程从基本概念实现开始,首先搭建了一个简单的NMT模...

5786
来自专栏技术随笔

[福利] 深入理解 RNNs & LSTM 网络学习资料图解

7379
来自专栏专注研发

扩展欧几里得算法

    有两个数 a b,现在,我们要求 a b 的最大公约数,怎么求?枚举他们的因子?不现实,当 a b 很大的时候,枚举显得那么的naïve ,那怎么做?

4813
来自专栏AI科技大本营的专栏

如何快速优化机器学习的模型参数

【导读】一般来说机器学习模型的优化没什么捷径可循。用什么架构,选择什么优化算法和参数既取决于我们对数据集的理解,也要不断地试错和修正。所以快速构建和测试模型的能...

1212
来自专栏谭正中的专栏

TensorFlow入门(3):使用神经网络拟合N元一次方程

现实中大部分情况是不能简单使用 N 元一次方程这样的公式表达的,神经网络的出现,给这类问题提供了一个很好的解决方法。本文继续给出一个简单的例子,使用 Tenso...

2.8K1
来自专栏闪电gogogo的专栏

压缩感知重构算法之压缩采样匹配追踪(CoSaMP)

压缩采样匹配追踪(CompressiveSampling MP)是D. Needell继ROMP之后提出的又一个具有较大影响力的重构算法。CoSaMP也是对OM...

38910

扫码关注云+社区

领取腾讯云代金券