貪婪算法中正交匹配追蹤算法gOMP的原理及仿真.doc

壓縮感知重構(gòu)算法之廣義正交匹配追蹤(gOMP)廣義正交匹配追蹤(Generalized OMP, gOMP)算法可以看作為OMP算法的一種推廣，由文獻(xiàn)1提出，第1作者本碩為哈工大畢業(yè)，發(fā)表此論文時在Korea University攻讀博士學(xué)位。OMP每次只選擇與殘差相關(guān)最大的一個，而gOMP則是簡單地選擇最大的S個。之所以這里表述為“簡單地選擇”是相比于ROMP之類算法的，不進(jìn)行任何其它處理，只是選擇最大的S個而已。0、符號說明如下： 壓縮觀測y=x，其中y為觀測所得向量M1，x為原信號N1（MN）。x一般不是稀疏的，但在某個變換域是稀疏的，即x=，其中為K稀疏的，即只有K個非零項。此時y=，令A(yù)=，則y=A。 (1)y為觀測所得向量，大小為M1 (2)x為原信號，大小為N1 (3)為K稀疏的，是信號在x在某變換域的稀疏表示 (4)稱為觀測矩陣、測量矩陣、測量基，大小為MN (5)稱為變換矩陣、變換基、稀疏矩陣、稀疏基、正交基字典矩陣，大小為NN (6)A稱為測度矩陣、傳感矩陣、CS信息算子，大小為MN上式中，一般有KMN，后面三個矩陣各個文獻(xiàn)的叫法不一，以后我將稱為測量矩陣、將稱為稀疏矩陣、將A稱為傳感矩陣。 注意：這里的稀疏表示模型為x=，所以傳感矩陣A=；而有些文獻(xiàn)中稀疏模型為=x，而一般為Hermite矩陣(實(shí)矩陣時稱為正交矩陣)，所以-1=H(實(shí)矩陣時為-1=T)，即x=H，所以傳感矩陣A=H，例如沙威的OMP例程中就是如此。1、gOMP重構(gòu)算法流程：2、廣義正交匹配追蹤(gOMP)MATLAB代碼(CS_gOMP.m) 本代碼完全是為了保證和前面的各算法代法格式一致，可以直接使用該實(shí)驗(yàn)室網(wǎng)站提供的代碼2壓縮包中的islsp_EstgOMP.m。plainview plaincopy1. functiontheta=CS_gOMP(y,A,K,S)2. %CS_gOMPSummaryofthisfunctiongoeshere3. %Version:1.0writtenbyjbb05232015-05-084. %Detailedexplanationgoeshere5. %y=Phi*x6. %x=Psi*theta7. %y=Phi*Psi*theta8. %令A(yù)=Phi*Psi,則y=A*theta9. %現(xiàn)在已知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/paper/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*theta_ls;%更新殘差46. Pos_theta=Sk;47. ifnorm(r_n)1e-648. break;%quittheiteration49. end50. end51. theta(Pos_theta)=theta_ls;%恢復(fù)出的theta52. end3、gOMP單次重構(gòu)測試代碼(CS_Reconstuction_Test.m) 以下測試代碼基本與OMP單次重構(gòu)測試代碼一樣。也可參考該實(shí)驗(yàn)室網(wǎng)站提供的代碼2壓縮包中的Test_gOMP.m。plainview plaincopy1. %壓縮感知重構(gòu)算法測試2. clearall;closeall;clc;3. M=128;%觀測值個數(shù)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稀疏的，且位置是隨機(jī)的9. Psi=eye(N);%x本身是稀疏的，定義稀疏矩陣為單位陣x=Psi*theta10. Phi=randn(M,N)/sqrt(M);%測量矩陣為高斯矩陣11. A=Phi*Psi;%傳感矩陣12. y=Phi*x;%得到觀測向量y13. %恢復(fù)重構(gòu)信號x14. tic15. theta=CS_gOMP(y,A,K);16. x_r=Psi*theta;%x=Psi*theta17. toc18. %繪圖19. figure;20. plot(x_r,k.-);%繪出x的恢復(fù)信號21. holdon;22. plot(x,r);%繪出原信號x23. holdoff;24. legend(Recovery,Original)25. fprintf(n恢復(fù)殘差：);26. norm(x_r-x)%恢復(fù)殘差 運(yùn)行結(jié)果如下：（信號為隨機(jī)生成，所以每次結(jié)果均不一樣） 1）圖： 2）Command windows Elapsedtime is 0.155937 seconds. 恢復(fù)殘差： ans= 2.3426e-0144、信號稀疏度K與重構(gòu)成功概率關(guān)系曲線繪制例程代碼 以下測試代碼為了與文獻(xiàn)1的Fig.1作比較。由于暫未研究學(xué)習(xí)LP算法，所以相比于文獻(xiàn)1的Fig.1)缺少LP算法曲線，加入了SP算法。以下測試代碼與SAMP相應(yīng)的測試代碼基本一致，可以合并在一起運(yùn)行，只須在主循環(huán)內(nèi)多加幾種算法重構(gòu)就行。plainview plaincopy1. %壓縮感知重構(gòu)算法測試CS_Reconstuction_KtoPercentagegOMP.m2. %繪制參考文獻(xiàn)中的Fig.13. %Reference:JianWang,SeokbeopKwon,ByonghyoShim.Generalized4. %orthogonalmatchingpursuit,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. %參數(shù)配置初始化10. CNT=1000;%對于每組(K,M,N)，重復(fù)迭代次數(shù)11. N=256;%信號x的長度12. Psi=eye(N);%x本身是稀疏的，定義稀疏矩陣為單位陣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);%存儲恢復(fù)成功概率17. %主循環(huán)，遍歷每組(K,M,N)18. tic19. formm=1:length(M_set)20. M=M_set(mm);%本次測量值個數(shù)21. K_set=5:5:70;%信號x的稀疏度K沒必要全部遍歷，每隔5測試一個就可以了22. %存儲此測量值M下不同K的恢復(fù)成功概率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%每個觀測值個數(shù)均運(yùn)行CNT次29. Index_K=randperm(N);30. x=zeros(N,1);31. x(Index_K(1:K)=5*randn(K,1);%x為K稀疏的，且位置是隨機(jī)的32. Phi=randn(M,N)/sqrt(M);%測量矩陣為高斯矩陣33. A=Phi*Psi;%傳感矩陣34. y=Phi*x;%得到觀測向量y35. %(1)OMP36. theta=CS_OMP(y,A,K);%恢復(fù)重構(gòu)信號theta37. x_r=Psi*theta;%x=Psi*theta38. ifnorm(x_r-x)1e-6%如果殘差小于1e-6則認(rèn)為恢復(fù)成功39. P(1)=P(1)+1;40. end41. %(2)ROMP42. theta=CS_ROMP(y,A,K);%恢復(fù)重構(gòu)信號theta43. x_r=Psi*theta;%x=Psi*theta44. ifnorm(x_r-x)1e-6%如果殘差小于1e-6則認(rèn)為恢復(fù)成功45. P(2)=P(2)+1;46. end47. %(3)StOMP48. theta=CS_StOMP(y,A);%恢復(fù)重構(gòu)信號theta49. x_r=Psi*theta;%x=Psi*theta50. ifnorm(x_r-x)1e-6%如果殘差小于1e-6則認(rèn)為恢復(fù)成功51. P(3)=P(3)+1;52. end53. %(4)SP54. theta=CS_SP(y,A,K);%恢復(fù)重構(gòu)信號theta55. x_r=Psi*theta;%x=Psi*theta56. ifnorm(x_r-x)1e-6%如果殘差小于1e-6則認(rèn)為恢復(fù)成功57. P(4)=P(4)+1;58. end59. %(5)CoSaMP60. theta=CS_CoSaMP(y,A,K);%恢復(fù)重構(gòu)信號theta61. x_r=Psi*theta;%x=Psi*theta62. ifnorm(x_r-x)1e-6%如果殘差小于1e-6則認(rèn)為恢復(fù)成功63. P(5)=P(5)+1;64. end65. %(6)gOMP,S=366. theta=CS_gOMP(y,A,K,3);%恢復(fù)重構(gòu)信號theta67. x_r=Psi*theta;%x=Psi*theta68. ifnorm(x_r-x)1e-6%如果殘差小于1e-6則認(rèn)為恢復(fù)成功69. P(6)=P(6)+1;70. end71. %(7)gOMP,S=672. theta=CS_gOMP(y,A,K,6);%恢復(fù)重構(gòu)信號theta73. x_r=Psi*theta;%x=Psi*theta74. ifnorm(x_r-x)1e-6%如果殘差小于1e-6則認(rèn)為恢復(fù)成功75. P(7)=P(7)+1;76. end77. %(8)gOMP,S=978. theta=CS_gOMP(y,A,K,9);%恢復(fù)重構(gòu)信號theta79. x_r=Psi*theta;%x=Psi*theta80. ifnorm(x_r-x)1e-6%如果殘差小于1e-6則認(rèn)為恢復(fù)成功81. P(8)=P(8)+1;82. end83. end84. foriii=1:size(KIND,1)85. PercentageM(iii,kk)=P(iii)/CNT*100;%計算恢復(fù)概率86. end87. end88. forjjj=1:size(KIND,1)89. Percentage(1:length(K_set),mm,jjj)=PercentageM(jjj,:);90. end91. end92. toc93. saveKtoPercentage1000gOMP%運(yùn)行一次不容易，把變量全部存儲下來94. %繪圖95. S=-ks;-ko;-yd;-gv;-b*;-r.;-rx;-r+;96. figure;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的恢復(fù)信號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);110. xlabel(SparsitylevelK);111. ylabel(TheProbabilityofExactReconstruction);112. title(Prob.ofexactrecoveryvs.thesignalsparsityK(M=128,N=256)(Gaussian); 本程序在聯(lián)想ThinkPadE430C筆記本（4GB DDR3內(nèi)存，i5-3210）上運(yùn)行共耗時798.718246秒，程序中將所有數(shù)據(jù)均通過“save KtoPercentage1000gOMP”存儲了下來，以后可以再對數(shù)據(jù)進(jìn)行分析，只需“l(fā)oad KtoPercentage1000gOMP”即可。 本程序運(yùn)行結(jié)果：5、結(jié)語 我很好奇：為什么相比于OMP算法就是簡單每次多選幾列，重構(gòu)效果為什么這么好？居然比復(fù)雜的ROMP、CoSaMP、StOMP效果還要好 該課題