贪婪算法中正交匹配追踪算法gOMP的原理及仿真

《贪婪算法中正交匹配追踪算法gOMP的原理及仿真》由会员分享，可在线阅读，更多相关《贪婪算法中正交匹配追踪算法gOMP的原理及仿真（10页珍藏版）》请在装配图网上搜索。

1、压缩感知重构算法之广义正交匹配追踪(gOMP)广义正交匹配追踪(Generalized OMP, gOMP)算法可以看作为OMP算法的一种推广，由文献1提出，第1作者本硕为哈工大毕业，发表此论文时在Korea University攻读博士学位。OMP每次只选择与残差相关最大的一个，而gOMP则是简单地选择最大的S个。之所以这里表述为“简单地选择”是相比于ROMP之类算法的，不进行任何其它处理，只是选择最大的S个而已。0、符号说明如下： 压缩观测y=x，其中y为观测所得向量M1，x为原信号N1（MN）。x一般不是稀疏的，但在某个变换域是稀疏的，即x=，其中为K稀疏的，即只有K个非零项。此时y=，

2、令A=，则y=A。 (1)y为观测所得向量，大小为M1 (2)x为原信号，大小为N1 (3)为K稀疏的，是信号在x在某变换域的稀疏表示 (4)称为观测矩阵、测量矩阵、测量基，大小为MN (5)称为变换矩阵、变换基、稀疏矩阵、稀疏基、正交基字典矩阵，大小为NN (6)A称为测度矩阵、传感矩阵、CS信息算子，大小为MN上式中，一般有KMN，后面三个矩阵各个文献的叫法不一，以后我将称为测量矩阵、将称为稀疏矩阵、将A称为传感矩阵。 注意：这里的稀疏表示模型为x=，所以传感矩阵A=；而有些文献中稀疏模型为=x，而一般为Hermite矩阵(实矩阵时称为正交矩阵)，所以-1=H(实矩阵时为-1=T)，即x=

3、H，所以传感矩阵A=H，例如沙威的OMP例程中就是如此。1、gOMP重构算法流程：2、广义正交匹配追踪(gOMP)MATLAB代码(CS_gOMP.m) 本代码完全是为了保证和前面的各算法代法格式一致，可以直接使用该实验室网站提供的代码2压缩包中的islsp_EstgOMP.m。plainview plaincopy1. functiontheta=CS_gOMP(y,A,K,S)2. %CS_gOMPSummaryofthisfunctiongoeshere3. %Version:1.0writtenbyjbb05232015-05-084. %Detailedexplanationgoes

4、here5. %y=Phi*x6. %x=Psi*theta7. %y=Phi*Psi*theta8. %令A=Phi*Psi,则y=A*theta9. %现在已知y和A，求theta10. %Reference:JianWang,SeokbeopKwon,ByonghyoShim.Generalized11. %orthogonalmatchingpursuit,IEEETransactionsonSignalProcessing,12. %vol.60,no.12,pp.6202-6216,Dec.2012.13. %Availableat:http:/islab.snu.ac.kr/pa

5、per/tsp_gOMP.pdf14. ifnargin415. S=round(max(K/4,1);16. end17. y_rows,y_columns=size(y);18. ify_rowsM36. ifii=137. theta_ls=0;38. end39. break;40. end41. At=A(:,Sk);%将A的这几列组成矩阵At42. %y=At*theta，以下求theta的最小二乘解(LeastSquare)43. theta_ls=(At*At)(-1)*At*y;%最小二乘解44. %At*theta_ls是y在At)列空间上的正交投影45. r_n=y-At

6、*theta_ls;%更新残差46. Pos_theta=Sk;47. ifnorm(r_n)1e-648. break;%quittheiteration49. end50. end51. theta(Pos_theta)=theta_ls;%恢复出的theta52. end3、gOMP单次重构测试代码(CS_Reconstuction_Test.m) 以下测试代码基本与OMP单次重构测试代码一样。也可参考该实验室网站提供的代码2压缩包中的Test_gOMP.m。plainview plaincopy1. %压缩感知重构算法测试2. clearall;closeall;clc;3. M=12

7、8;%观测值个数4. N=256;%信号x的长度5. K=30;%信号x的稀疏度6. Index_K=randperm(N);7. x=zeros(N,1);8. x(Index_K(1:K)=5*randn(K,1);%x为K稀疏的，且位置是随机的9. Psi=eye(N);%x本身是稀疏的，定义稀疏矩阵为单位阵x=Psi*theta10. Phi=randn(M,N)/sqrt(M);%测量矩阵为高斯矩阵11. A=Phi*Psi;%传感矩阵12. y=Phi*x;%得到观测向量y13. %恢复重构信号x14. tic15. theta=CS_gOMP(y,A,K);16. x_r=Psi

8、*theta;%x=Psi*theta17. toc18. %绘图19. figure;20. plot(x_r,k.-);%绘出x的恢复信号21. holdon;22. plot(x,r);%绘出原信号x23. holdoff;24. legend(Recovery,Original)25. fprintf(n恢复残差：);26. norm(x_r-x)%恢复残差 运行结果如下：（信号为随机生成，所以每次结果均不一样） 1）图： 2）Command windows Elapsedtime is 0.155937 seconds. 恢复残差： ans= 2.3426e-0144、信号稀疏度K与

9、重构成功概率关系曲线绘制例程代码 以下测试代码为了与文献1的Fig.1作比较。由于暂未研究学习LP算法，所以相比于文献1的Fig.1)缺少LP算法曲线，加入了SP算法。以下测试代码与SAMP相应的测试代码基本一致，可以合并在一起运行，只须在主循环内多加几种算法重构就行。plainview plaincopy1. %压缩感知重构算法测试CS_Reconstuction_KtoPercentagegOMP.m2. %绘制参考文献中的Fig.13. %Reference:JianWang,SeokbeopKwon,ByonghyoShim.Generalized4. %orthogonalmatch

10、ingpursuit,IEEETransactionsonSignalProcessing,5. %vol.60,no.12,pp.6202-6216,Dec.2012.6. %Availableat:http:/islab.snu.ac.kr/paper/tsp_gOMP.pdf7. %Elapsedtimeis798.718246seconds.(20150509pm)8. clearall;closeall;clc;9. %参数配置初始化10. CNT=1000;%对于每组(K,M,N)，重复迭代次数11. N=256;%信号x的长度12. Psi=eye(N);%x本身是稀疏的，定义稀

11、疏矩阵为单位阵x=Psi*theta13. M_set=128;%测量值集合14. KIND=OMP;ROMP;StOMP;SP;CoSaMP;.15. gOMP(s=3);gOMP(s=6);gOMP(s=9);16. Percentage=zeros(N,length(M_set),size(KIND,1);%存储恢复成功概率17. %主循环，遍历每组(K,M,N)18. tic19. formm=1:length(M_set)20. M=M_set(mm);%本次测量值个数21. K_set=5:5:70;%信号x的稀疏度K没必要全部遍历，每隔5测试一个就可以了22. %存储此测量值M下

12、不同K的恢复成功概率23. PercentageM=zeros(size(KIND,1),length(K_set);24. forkk=1:length(K_set)25. K=K_set(kk);%本次信号x的稀疏度K26. P=zeros(1,size(KIND,1);27. fprintf(M=%d,K=%dn,M,K);28. forcnt=1:CNT%每个观测值个数均运行CNT次29. Index_K=randperm(N);30. x=zeros(N,1);31. x(Index_K(1:K)=5*randn(K,1);%x为K稀疏的，且位置是随机的32. Phi=randn(M

13、,N)/sqrt(M);%测量矩阵为高斯矩阵33. A=Phi*Psi;%传感矩阵34. y=Phi*x;%得到观测向量y35. %(1)OMP36. theta=CS_OMP(y,A,K);%恢复重构信号theta37. x_r=Psi*theta;%x=Psi*theta38. ifnorm(x_r-x)1e-6%如果残差小于1e-6则认为恢复成功39. P(1)=P(1)+1;40. end41. %(2)ROMP42. theta=CS_ROMP(y,A,K);%恢复重构信号theta43. x_r=Psi*theta;%x=Psi*theta44. ifnorm(x_r-x)1e-6

14、%如果残差小于1e-6则认为恢复成功45. P(2)=P(2)+1;46. end47. %(3)StOMP48. theta=CS_StOMP(y,A);%恢复重构信号theta49. x_r=Psi*theta;%x=Psi*theta50. ifnorm(x_r-x)1e-6%如果残差小于1e-6则认为恢复成功51. P(3)=P(3)+1;52. end53. %(4)SP54. theta=CS_SP(y,A,K);%恢复重构信号theta55. x_r=Psi*theta;%x=Psi*theta56. ifnorm(x_r-x)1e-6%如果残差小于1e-6则认为恢复成功57.

15、P(4)=P(4)+1;58. end59. %(5)CoSaMP60. theta=CS_CoSaMP(y,A,K);%恢复重构信号theta61. x_r=Psi*theta;%x=Psi*theta62. ifnorm(x_r-x)1e-6%如果残差小于1e-6则认为恢复成功63. P(5)=P(5)+1;64. end65. %(6)gOMP,S=366. theta=CS_gOMP(y,A,K,3);%恢复重构信号theta67. x_r=Psi*theta;%x=Psi*theta68. ifnorm(x_r-x)1e-6%如果残差小于1e-6则认为恢复成功69. P(6)=P(6

16、)+1;70. end71. %(7)gOMP,S=672. theta=CS_gOMP(y,A,K,6);%恢复重构信号theta73. x_r=Psi*theta;%x=Psi*theta74. ifnorm(x_r-x)1e-6%如果残差小于1e-6则认为恢复成功75. P(7)=P(7)+1;76. end77. %(8)gOMP,S=978. theta=CS_gOMP(y,A,K,9);%恢复重构信号theta79. x_r=Psi*theta;%x=Psi*theta80. ifnorm(x_r-x)1e-6%如果残差小于1e-6则认为恢复成功81. P(8)=P(8)+1;82

17、. end83. end84. foriii=1:size(KIND,1)85. PercentageM(iii,kk)=P(iii)/CNT*100;%计算恢复概率86. end87. end88. forjjj=1:size(KIND,1)89. Percentage(1:length(K_set),mm,jjj)=PercentageM(jjj,:);90. end91. end92. toc93. saveKtoPercentage1000gOMP%运行一次不容易，把变量全部存储下来94. %绘图95. S=-ks;-ko;-yd;-gv;-b*;-r.;-rx;-r+;96. fig

18、ure;97. formm=1:length(M_set)98. M=M_set(mm);99. K_set=5:5:70;100. L_Kset=length(K_set);101. forii=1:size(KIND,1)102. plot(K_set,Percentage(1:L_Kset,mm,ii),S(ii,:);%绘出x的恢复信号103. holdon;104. end105. end106. holdoff;107. xlim(570);108. legend(OMP,ROMP,StOMP,SP,CoSaMP,.109. gOMP(s=3),gOMP(s=6),gOMP(s=9

19、);110. xlabel(SparsitylevelK);111. ylabel(TheProbabilityofExactReconstruction);112. title(Prob.ofexactrecoveryvs.thesignalsparsityK(M=128,N=256)(Gaussian); 本程序在联想ThinkPadE430C笔记本（4GB DDR3内存，i5-3210）上运行共耗时798.718246秒，程序中将所有数据均通过“save KtoPercentage1000gOMP”存储了下来，以后可以再对数据进行分析，只需“load KtoPercentage1000g

20、OMP”即可。 本程序运行结果：5、结语 我很好奇：为什么相比于OMP算法就是简单每次多选几列，重构效果为什么这么好？居然比复杂的ROMP、CoSaMP、StOMP效果还要好 该课题组还提出了MMP算法，可参见文献3。更多关于该课题组的信息可去官方网站查询：http:/islab.snu.ac.kr/，也可直接查看发表的文章：http:/islab.snu.ac.kr/publication.html。 文献1最后有两个TABLE，分别是算法的流程和复杂度总结：谁能告诉我TABLE I表头中的DELETE在这里是什么意思啊？不是删除的意思么？6、参考文献【1】Jian Wang, Seokbe

21、op Kwon,Byonghyo Shim. Generalized orthogonalmatching pursuit, IEEE Transactions on Signal Processing, vol. 60, no. 12, pp.6202-6216, Dec. 2012.Available at: http:/islab.snu.ac.kr/paper/tsp_gOMP.pdf【2】http:/islab.snu.ac.kr/paper/gOMP.zip【3】S. Kwon, J. Wang andB. Shim, Multipath matching pursuit, IEEE Transactions on Information Theory,vol. 60, no. 5, pp. 2986-3001, May 2014.Available at: http:/islab.snu.ac.kr/paper/TIT_MMP2014.pdf