贪婪算法中正交匹配追踪算法gOMP的基础学习知识原理及仿真.docx
- 文档编号:24217384
- 上传时间:2023-05-25
- 格式:DOCX
- 页数:15
- 大小:316.11KB
贪婪算法中正交匹配追踪算法gOMP的基础学习知识原理及仿真.docx
《贪婪算法中正交匹配追踪算法gOMP的基础学习知识原理及仿真.docx》由会员分享,可在线阅读,更多相关《贪婪算法中正交匹配追踪算法gOMP的基础学习知识原理及仿真.docx(15页珍藏版)》请在冰豆网上搜索。
贪婪算法中正交匹配追踪算法gOMP的基础学习知识原理及仿真
压缩感知重构算法之广义正交匹配追踪(gOMP)
广义正交匹配追踪(GeneralizedOMP,gOMP)算法可以看作为OMP算法的一种推广,由文献[1]提出,第1作者本硕为哈工大毕业,发表此论文时在KoreaUniversity攻读博士学位。
OMP每次只选择与残差相关最大的一个,而gOMP则是简单地选择最大的S个。
之所以这里表述为“简单地选择”是相比于ROMP之类算法的,不进行任何其它处理,只是选择最大的S个而已。
0、符号说明如下:
压缩观测y=Φx,其中y为观测所得向量M×1,x为原信号N×1(M< x一般不是稀疏的,但在某个变换域Ψ是稀疏的,即x=Ψθ,其中θ为K稀疏的,即θ只有K个非零项。 此时y=ΦΨθ,令A=ΦΨ,则y=Aθ。 (1) y为观测所得向量,大小为M×1 (2)x为原信号,大小为N×1 (3)θ为K稀疏的,是信号在x在某变换域的稀疏表示 (4) Φ称为观测矩阵、测量矩阵、测量基,大小为M×N (5) Ψ称为变换矩阵、变换基、稀疏矩阵、稀疏基、正交基字典矩阵,大小为N×N (6)A称为测度矩阵、传感矩阵、CS信息算子,大小为M×N 上式中,一般有K< 注意: 这里的稀疏表示模型为x=Ψθ,所以传感矩阵A=ΦΨ;而有些文献中稀疏模型为θ=Ψx,而一般Ψ为Hermite矩阵(实矩阵时称为正交矩阵),所以Ψ-1=ΨH (实矩阵时为Ψ-1=ΨT),即x=ΨHθ,所以传感矩阵A=ΦΨH,例如沙威的OMP例程中就是如此。 1、gOMP重构算法流程: 2、广义正交匹配追踪(gOMP)MATLAB代码(CS_gOMP.m) 本代码完全是为了保证和前面的各算法代法格式一致,可以直接使用该实验室网站提供的代码[2]压缩包中的islsp_EstgOMP.m。 [plain] viewplaincopy 1.function [ theta ] = CS_gOMP( y,A,K,S ) 2.%CS_gOMP Summary of this function goes here 3.%Version: 1.0 written by jbb0523 @2015-05-08 4.% Detailed explanation goes here 5.% y = Phi * x 6.% x = Psi * theta 7.% y = Phi*Psi * theta 8.% 令 A = Phi*Psi, 则y=A*theta 9.% 现在已知y和A,求theta 10.% Reference: Jian Wang, Seokbeop Kwon, Byonghyo Shim. Generalized 11.% orthogonal matching pursuit, IEEE Transactions on Signal Processing, 12.% vol. 60, no. 12, pp. 6202-6216, Dec. 2012. 13.% Available at: http: //islab.snu.ac.kr/paper/tsp_gOMP.pdf 14. if nargin < 4 15. S = round(max(K/4, 1)); 16. end 17. [y_rows,y_columns] = size(y); 18. if y_rows 19. y = y';%y should be a column vector 20. end 21. [M,N] = size(A);%传感矩阵A为M*N矩阵 22. theta = zeros(N,1);%用来存储恢复的theta(列向量) 23. Pos_theta = [];%用来迭代过程中存储A被选择的列序号 24. r_n = y;%初始化残差(residual)为y 25. for ii=1: K%迭代K次,K为稀疏度 26. product = A'*r_n;%传感矩阵A各列与残差的内积 27. [val,pos]=sort(abs(product),'descend');%降序排列 28. Sk = union(Pos_theta,pos(1: S));%选出最大的S个 29. if length(Sk)==length(Pos_theta) 30. if ii == 1 31. theta_ls = 0; 32. end 33. break; 34. end 35. if length(Sk)>M 36. if ii == 1 37. theta_ls = 0; 38. end 39. break; 40. end 41. At = A(: Sk);%将A的这几列组成矩阵At 42. %y=At*theta,以下求theta的最小二乘解(Least Square) 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. if norm(r_n)<1e-6 48. break;%quit the iteration 49. end 50. end 51. theta(Pos_theta)=theta_ls;%恢复出的theta 52.end 3、gOMP单次重构测试代码(CS_Reconstuction_Test.m) 以下测试代码基本与OMP单次重构测试代码一样。 也可参考该实验室网站提供的代码[2]压缩包中的Test_gOMP.m。 [plain] viewplaincopy 1.%压缩感知重构算法测试 2.clear all;close all;clc; 3.M = 128;%观测值个数 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*theta 10.Phi = randn(M,N)/sqrt(M);%测量矩阵为高斯矩阵 11.A = Phi * Psi;%传感矩阵 12.y = Phi * x;%得到观测向量y 13.%% 恢复重构信号x 14.tic 15.theta = CS_gOMP( y,A,K); 16.x_r = Psi * theta;% x=Psi * theta 17.toc 18.%% 绘图 19.figure; 20.plot(x_r,'k.-');%绘出x的恢复信号 21.hold on; 22.plot(x,'r');%绘出原信号x 23.hold off; 24.legend('Recovery','Original') 25.fprintf('\n恢复残差: '); 26.norm(x_r-x)%恢复残差 运行结果如下: (信号为随机生成,所以每次结果均不一样) 1) 图: 2)Command windows Elapsedtimeis0.155937seconds. 恢复残差: ans= 2.3426e-014 4、信号稀疏度K与重构成功概率关系曲线绘制例程代码 以下测试代码为了与文献[1]的Fig.1作比较。 由于暂未研究学习LP算法,所以相比于文献[1]的Fig.1)缺少LP算法曲线,加入了SP算法。 以下测试代码与SAMP相应的测试代码基本一致,可以合并在一起运行,只须在主循环内多加几种算法重构就行。 [plain] viewplaincopy 1.%压缩感知重构算法测试CS_Reconstuction_KtoPercentagegOMP.m 2.% 绘制参考文献中的Fig.1 3.% Reference: Jian Wang, Seokbeop Kwon, Byonghyo Shim. Generalized 4.% orthogonal matching pursuit, IEEE Transactions on Signal Processing, 5.% vol. 60, no. 12, pp. 6202-6216, Dec. 2012. 6.% Available at: http: //islab.snu.ac.kr/paper/tsp_gOMP.pdf 7.% Elapsed time is 798.718246 seconds.(@20150509pm) 8.clear all;close all;clc; 9.%% 参数配置初始化 10.CNT = 1000;%对于每组(K,M,N),重复迭代次数 11.N = 256;%信号x的长度 12.Psi = eye(N);%x本身是稀疏的,定义稀疏矩阵为单位阵x=Psi*theta 13.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.tic 19.for mm = 1: length(M_set) 20. M = M_set(mm);%本次测量值个数 21. K_set = 5: 5: 70;%信号x的稀疏度K没必要全部遍历,每隔5测试一个就可以了 22. %存储此测量值M下不同K的恢复成功概率 23. PercentageM = zeros(size(KIND,1),length(K_set)); 24. for kk = 1: length(K_set) 25. K = K_set(kk);%本次信号x的稀疏度K 26. P = zeros(1,size(KIND,1)); 27. fprintf('M=%d,K=%d\n',M,K); 28. for cnt = 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,N)/sqrt(M);%测量矩阵为高斯矩阵 33. A = Phi * Psi;%传感矩阵 34. y = Phi * x;%得到观测向量y 35. % (1)OMP 36. theta = CS_OMP(y,A,K);%恢复重构信号theta 37. x_r = Psi * theta;% x=Psi * theta 38. if norm(x_r-x)<1e-6%如果残差小于1e-6则认为恢复成功 39. P (1) = P (1) + 1; 40. end 41. % (2)ROMP 42. theta = CS_ROMP(y,A,K);%恢复重构信号theta 43. x_r = Psi * theta;% x=Psi * theta 44. if norm(x_r-x)<1e-6%如果残差小于1e-6则认为恢复成功 45. P (2) = P (2) + 1; 46. end 47. %(3)StOMP 48. theta = CS_StOMP(y,A);%恢复重构信号theta 49. x_r = Psi * theta;% x=Psi * theta 50. if norm(x_r-x)<1e-6%如果残差小于1e-6则认为恢复成功 51. P(3) = P(3) + 1; 52. end 53. %(4)SP 54. theta = CS_SP(y,A,K);%恢复重构信号theta 55. x_r = Psi * theta;% x=Psi * theta 56. if norm(x_r-x)<1e-6%如果残差小于1e-6则认为恢复成功 57. P(4) = P(4) + 1; 58. end 59. %(5)CoSaMP 60. theta = CS_CoSaMP(y,A,K);%恢复重构信号theta 61. x_r = Psi * theta;% x=Psi * theta 62. if norm(x_r-x)<1e-6%如果残差小于1e-6则认为恢复成功 63. P(5) = P(5) + 1; 64. end 65. %(6)gOMP,S=3 66. theta = CS_gOMP(y,A,K,3);%恢复重构信号theta 67. x_r = Psi * theta;% x=Psi * theta 68. if norm(x_r-x)<1e-6%如果残差小于1e-6则认为恢复成功 69. P(6) = P(6) + 1; 70. end 71. %(7)gOMP,S=6 72. theta = CS_gOMP(y,A,K,6);%恢复重构信号theta 73. x_r = Psi * theta;% x=Psi * theta 74. if norm(x_r-x)<1e-6%如果残差小于1e-6则认为恢复成功 75. P(7) = P(7) + 1; 76. end 77. %(8)gOMP,S=9 78. theta = CS_gOMP(y,A,K,9);%恢复重构信号theta 79. x_r = Psi * theta;% x=Psi * theta 80. if norm(x_r-x)<1e-6%如果残差小于1e-6则认为恢复成功 81. P(8) = P(8) + 1; 82. end 83. end 84. for iii = 1: size(KIND,1) 85. PercentageM(iii,kk) = P(iii)/CNT*100;%计算恢复概率 86. end 87. end 88. for jjj = 1: size(KIND,1) 89. Percentage(1: length(K_set),mm,jjj) = PercentageM(jjj,: ); 90. end 91.end 92.toc 93.save KtoPercentage1000gOMP %运行一次不容易,把变量全部存储下来 94.%% 绘图 95.S = ['-ks';'-ko';'-yd';'-gv';'-b*';'-r.';'-rx';'-r+']; 96.figure; 97.for mm = 1: length(M_set) 98. M = M_set(mm); 99. K_set = 5: 5: 70; 100. L_Kset = length(K_set); 101. for ii = 1: size(KIND,1) 102. plot(K_set,Percentage(1: L_Kset,mm,ii),S(ii,: ));%绘出x的恢复信号 103. hold on; 104. end 105.end 106.hold off; 107.xlim([5 70]); 108.legend('OMP','ROMP','StOMP','SP','CoSaMP',... 109. 'gOMP(s=3)','gOMP(s=6)','gOMP(s=9)'); 110.xlabel('Sparsity level K'); 111.ylabel('The Probability of Exact Reconstruction'); 112.title('Prob. of exact recovery vs. the signal sparsity K(M=128,N=256)(Gaussian)'); 本程序在联想ThinkPadE430C笔记本(4GBDDR3内存,i5-3210)上运行共耗时798.718246秒,程序中将所有数据均通过“saveKtoPercentage1000gOMP”存储了下来,以后可以再对数据进行分析,只需“loadKtoPercentage1000gOMP”即可。 本程序运行结果: 5、结语 我很好奇: 为什么相比于OMP算法就是简单每次多选几列,重构效果为什么这么好? 居然比复杂的ROMP、CoSaMP、StOMP效果还要好…… 该课题组还提出了MMP算法,可参见文献[3]。 更多关于该课题组的信息可去官方网站查询: http: //islab.snu.ac.kr/,也可直接查看发表的文章: http: //islab.snu.ac.kr/publication.html。 文献[1]最后有两个TABLE,分别是算法的流程和复杂度总结: 谁能告诉我TABLEI表头中的DELETE在这里是什么意思啊? 不是删除的意思么? 6、参考文献 【1】JianWang,SeokbeopKwon,ByonghyoShim. Generalizedorthogonalmatchingpursuit,IEEETransactionsonSignalProcessing,vol.60,no.12,pp.6202-6216,Dec.2012. Availableat: http: //islab.snu.ac.kr/paper/tsp_gOMP.pdf 【2】http: //islab.snu.ac.kr/paper/gOMP.zip 【3】S.Kwon,J.WangandB.Shim,Multipathmatchingpursuit,IEEETransactionsonInformationTheory,vol.60,no.5,pp.2986-3001,May2014. Availableat: http: //islab.snu.ac.kr/paper/TIT_MMP2014.pdf
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 贪婪 算法 中正 匹配 追踪 gOMP 基础 学习 知识 原理 仿真
![提示](https://static.bdocx.com/images/bang_tan.gif)