贪婪算法中,SP算法的原理介绍及MATLAB仿真

1、压缩感知重构算法之子空间追踪(SP)如果掌握了压缩采样匹配追踪(CoSaMP)后，再去学习子空间追踪 (Subspace Pursuit)是一件非常简单的事情,因为它们几乎是完全一样的。SP的提出时间比 CoSaMP提出时间略晚，首个论文版本是参考文献 1,后来更新了两次，最后在 IEEE Transactions on Information Theory 发表2。从算法角度来讲，SP与CoSaMP差别非常小，这一点作者也意识到 了，在文献1首页的左下角就有注释：* AL die time of writing this maiiuript. tl)c audiors became awue

2、 of the related werk by J. Tmpp* D. Needell and R, Vershynin 111 where similar reconstruction alguriilmih are designed. Our results dcvuluped inckpcndenUy, and we believe (hut there arc significant differences in these two proposed rcconstmction ripproaehcs.在文献2第2页提到了 SP与CoSaMP的具体不同：Al the time of w

3、riting of this paper, the authors became aware of the related work by L Tropp. D. Needell, and R, Vershynin 12, describing a similar rcconstniction algorithtiL The main difference between the SP algorithm and the CoSAMP algorithm of 12J is in the manner in which new candidates arc added to the list.

4、 In each iteration, in the SP algorithm, only K new candidates are added, while the CoSAMP algorithm adds 2K vectors. This makes the SP algorithm computationally mare efficient, but the underlying analysis more complex. In addition, (he restricted isomeiry ccrnKtanl for which the SP algorithm is gua

5、r an teed to converge is larger than the one presented in 12. Finally. this paper also contains an analysis of the number of iterations needed for reconstruction of a sparse signal (see Theorem 6 for details), for which there is no counterpart in the CoSAMP study.从上面可以知道，SP 与 CoSaMP 主要区别在于 Ineach it

6、eration, in the SP algorithm, only K new candidates areadded , while theCoSAMP algorithm adds 2K vectors.，即“SP 每次选择 K 个原子，而 CoSaMP 则选择 2K 个原子;这样带来的好处是 “This makes the SP algorithm computationally moreefficient, 。”以下是文献2中的给出的SP算法流程:Algorithm 1 Subspace Pursuit AlgorithmInput: K、典 yInitialization:7* =

7、 indices corresponding to the largest magniude entries in the vector y.2)婢=resid冢).iteration: At the/th iteration, go through the following steps.1)型=T；-i |J K indices corresponding to the largest magnitude entries in the vectnr 中电/1 .Set xp =%V7 ; K indices corresponding to the largest magnitude el

8、ements of 的J.二 resid (If |j|2 “必TJLlet= T$-i and quit theiteration,Ouipui:1) The estimated signal k satisfying 仝1，.少欢-/ = 0 and xT 明妙这个算法流程的初始化（Initialization）其实就是类似于 CoSaMP的第1次迭代，注意第（1）步中选择了 K个原子：“K indices correspo nding to the largest magnitude entries ,在 CoSaMP 里这里要选择 2K 个最大的原子，后面的其 它流程都一样。这里第（5

9、）步增加了一个停止迭代的条件：当残差经过迭代后却变大了的时候就停止迭代。不只是SP作者认识到了自己的算法与CoSaMP的高度相似性，CoSaMP的作者也同样关注到了SP算法，在文献3中就提到：After initially presented this work. Did and Milenkovic developed an Algorithm called Subspace1 Pursuit that is very similar to CoSaMP. They established that their algorithm offers perfamiancc guarantees

10、analogous with those for CoSaMP, Src 10 for details.文献3是CoSaMP原始提出文献的第 2个版本，文献3的早期版本4是没有提及SP算法的。鉴于SP与CoSaMP如此相似，这里不就再单独给出SP的步骤了，参考压缩感知重构算法之压缩采样匹配追踪（CoSaMP） ,只需将第（2）步中的2K改为K即可。引用文献5的3,5节中的几句话：贪婪类算法虽然复杂度低运行速度快，但其重构精度却不如BP类算法,为了寻求复杂度和精度更好地折中，SP算法应运而生”，“SPI法与CoSaMP算法一样其基本思想也是借用回溯的思想，在每步迭代过程中重新估计所有候选者的可信

11、赖性：SFW法与CoSaMP算法有着类似的性质与优缺点”子空间追踪代码可实现如下(CS_SP.m),通过对比可以知道该程序与CoSaMP的代码基本完全一致。本代码未考虑文献2中的给出的SP算法流程的第(5)步。代码可参见参考文献6中的Demo_CS_SP.m 。plain view plaincopyfunction theta = CS_SP( y,A,K )%CS_SP Summary of this function goes here%Version: 1.0 written by jbb0523 2015-05-01% Detailed explanation goes here%

12、y = Phi * x% x = Psi * theta% y = Phi*Psi * theta% 令 A = Phi*Psi,贝U y=A*theta% K is the sparsity level%现在已知y和A,求theta% Reference:Dai W , Milenkovic O . Subspace pursuit for compressive sensing% signal reconstructionJ. IEEE Transactions on Information Theory,% 2009 , 55(5) : 2230-2249.y_rows,y_column

13、s = size(y);if y_rowsy_columnsy = y;%y should be a column vectorendM,N = size(A);%传感矢!阵 A 为 M*N矩阵theta = zeros(N,1);%用来存储恢复的theta( 列向量)Pos_theta = ;%用来迭代过程中存储 A被选择的列序号r_n = y;% 初始化残差(residual) 为 yfor kk=1:K% 最多迭代K次%(1) Identificationproduct = A*r_n;%传感矢1阵A各列与残差的内积val,pos=sort(abs(product),descend);J

14、s = pos(1:K);%选出内积值最大的K列%(2) Support MergerIs = union(Pos_theta,Js);%Pos_theta与Js并集%(3) Estimation%At的行数要大于列数，此为最小二乘的基础 (列线性无关)if length(Is)=MAt = A(:,Is);%将A的这几列组成矩阵 Atelse%At的列数大于行数，列必为线性相关的,At*At 将不可逆break;%跳出 for 循环end%y=At*theta,以下求 theta 的最小二乘解(Least Square)theta_ls = (At*At)A(-1)*At*y;%最小二乘解%

15、(4) Pruningval,pos=sort(abs(theta_ls),descend);%(5) Sample UpdatePos_theta = Is(pos(1:K);theta_ls = theta_ls(pos(1:K);%At(:,pos(1:K)*theta_ls是 y 在 At(:,pos(1:K)列空间上的正交投影r_n = y - At(:,pos(1:K)*theta_ls;%更新残差if norm(r_n)1e-6%Repeat thesteps until r=0break;%跳出 for 循环endendtheta(Pos_theta)=theta_ls;%恢复

16、出的 thetaend鉴于SP与CoSaMP的极其相似性，这里就不再给出单次重构和测量数M与重构成功概率关系曲线绘制例程代码了，只需将 CoSaMP中调用CS_CoSaMP函数的部分改为调用CS_SP即可，无须任何其它改动。这里给出对比两种重构算法所绘制的测量数M与重构成功概率关系曲线的例程代码，只有这样才可以看出两种算法的重构性能优劣，以下是在分别运行完SP与CoSaMP的测量数M与重构成功概率关系曲线绘制例程代码的基础上，即已经存储了数据 CoSaMPMtoPercentage1000.mat 和 SPMtoPercentage1000.mat :plain view plaincopyc

17、lear all;close all;clc;load CoSaMPMtoPercentage1000;PercentageCoSaMP = Percentage;load SPMtoPercentage1000;PercentageSP = Percentage;S1 = -ks;-ko;-kd;-kv;-k*;S2 = -rs;-ro;-rd;-rv;-r*;figure;for kk = 1:length(K_set)K = K_set(kk);M_set = 2*K:5:N;L_Mset = length(M_set);plot(M_set,PercentageCoSaMP(kk,1:

18、L_Mset),S1(kk,:);%绘出 x 的恢复信号hold on;plot(M_set,PercentageSP(kk,1:L_Mset),S2(kk,:);%绘出 x 的恢复信号endhold off;xlim(0 256);legend(CoSaK=4,SPK=4,CoSaK=12,SPK=12,CoSaK=20,.SPK=20,CoSaK=28,SPK=28,CoSaK=36,SPK=36);xlabel(Number of measurements(M);ylabel(Percentage recovered);title(Percentage of input signals

19、recovered correctly(N=256)(Gaussian);运彳谭果如下:Percerrtage of input signals recovered correctly(N=256)(Gajssianj00903070605040302010口f fy /铲现聚a*”售工7呼触摩岫率就值阖tjwte日一CoSaK4 .1f TJ-SPK=4I iJ-e-CQSaK=121TJ-SPK=12JJ8-CoSaK=20-SPK=20,tIj J J7-CoSaK=2611 fI j1V-SPK=2BI1CoS 水 =36-t-SPK=3E1-/I i1L-J1wa150200260Mumber of measuremeritsM可以发现在H较小时SN等好于8s洲R为闻变方寸二者重构性能几乎一样。参考文献：1Dai W , Milenkovic O. Subspace pursuitfor compressive sensing: Closing the gap between performance andcomplexity. HYPERLINK /pdf/0803.0811v1.pdf /pdf/0803.08