IEEE Trans 2008 Gradient Pursuits论文学习
IEEE Trans 2008 Gradient Pursuits论文学习
闪电gogogo
发布于 2018-01-08 13:27:38
发布于 2018-01-08 13:27:38
之前所学习的论文中求解稀疏解的时候一般采用的都是最小二乘方法进行计算，为了降低计算复杂度和减少内存，这篇论文梯度追踪，属于贪婪算法中一种。主要为三种：梯度（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;
endACGP（代码名称为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
本文参与 腾讯云自媒体分享计划，分享自作者个人站点/博客。
原始发表：2017-12-19 ，如有侵权请联系 cloudcommunity@tencent.com 删除
评论
登录后参与评论
推荐阅读