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

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

2、的，即。只有K个非零项。此时 丫二令人二3，则y=A0o(1) y为观测所得向量，大小为Mx 1(2)x为原信号，大小为 NX 1(3)。为K稀疏的，是信号在 x在某变换域的稀疏表示(4)称为观测矩阵、 测量矩阵、测量基，大小为 MX N(5) W称为变换矩阵、变换基、稀疏矩阵、稀疏基、正交基字典矩阵，大小为NXN(6)A称为测度矩阵、 传感矩阵、CS信息算子，大小为 MX N上式中，一般有KMN ,后面三个矩阵各个文献的叫法不一，以后我将 称为测量矩阵、将W称为稀疏矩阵、将A称为传感矩阵。注意：这里的稀疏表示模型为 x=WQ所以传感矩阵 A二照 而有些文献中稀疏模型为 0 =3*而一般 W为

3、 Hermite矩阵(实矩阵时称为正交矩阵),所以W-1=WH (实矩阵时为 W-1=WT),即x=H 0,所以传感矩阵 A二,例 如沙威的OMP例程中就是如此。1、gOMP重构算法流程:输入士(IJAfxJV的传感矩阵A =中死(2) N xl维观颖恂量yG)信号的稀疏度K(4)每次选择的原子个数弘默认取值为K4,若KdCW默认取1输出口)信号稀疏表示系数估计8(4 万幻维度差m=f-4a以下流程中；6表示残差，表示迭代次数，。表示空集，4表示f次迭代的索引(列 序号)集合，4表示第八欠迭代找到的索弓1(列序号)，与表示矩阵的第/列，4表示搜索 引4选出的矩阵川的列集合(大小为财的比阵斗为发

4、1的列向量，符号u表示集合 并运算，丫)表示求向量内积.(1泡始化弓=，4三0，4二0工二1 ；2)计算值=1加j。5一(即计算：11田；)/M J E M),选辉观中最大的个值，将这些值对应A的列序号构成集合/ (列序号集合)；! ?)告/1t = 4d4 Aj = Udty(for ally 17)；(4)求,三金的最小二乘解：g三argn由司|三父人父,； a(5)更新麻片=尸_4=1_4|j(6) = +1,如果 M匕则退回第步I否则停止迭代进入第17)步； 重构所得自在几处有非零项，其值分别为最后一次迭代厮得& 0注I：得到启后，利用稀疏矩阵可得重构信号至二汨12、广义正交匹配追踪

5、(gOMP)MATLAB 代码(CS_gOMP.m)本代码完全是为了保证和前面的各算法代法格式一致，可以直接使用该实验室网站提供的代码2压缩包中的islsp_EstgOMP.m 。plain view plaincopyfunction theta = CS_gOMP( y,A,K,S )%CS_gOMP Summary of this function goes here%Version: 1.0 written by jbb0523 2015-05-08% Detailed explanation goes here% y = Phi * x% x = Psi * theta% y = P

6、hi*Psi * theta% 令 A = Phi*Psi, 贝U y=A*theta%现在已知y和A,求theta% Reference: Jian Wang, Seokbeop Kwon, Byonghyo Shim. Generalized% orthogonal matching pursuit, IEEE Transactions on Signal Processing,% vol. 60, no. 12, pp. 6202-6216, Dec. 2012.% Available at: HYPERLINK http:/islab.snu.ac.kr/paper/tsp_gOMP.

7、pdf http:/islab.snu.ac.kr/paper/tsp_gOMP.pdfif nargin 4S = round(max(K/4, 1);endy_rows,y_columns = size(y);if y_rowsMif ii = 1theta_ls = 0;endbreak;endAt = A(:,Sk);%将A的这几列组成矩阵 At%y=At*theta,以下求 theta 的最小二乘解(Least Square)theta_ls = (At*At)A(-1)*At*y;%最小二乘解%At*theta_ls是y在At)列空间上的正交投影r_n = y - At*theta

8、_ls;%更新残差Pos_theta = Sk;if norm(r_n)1e-6break;%quit the iterationendendtheta(Pos_theta)=theta_ls;%恢复出的 thetaend3、gOMP 单次重构测试代码(CS_Reconstuction_Test.m)以下测试代码基本与 OMP单次重构测试代码一样。也可参考该实验室网站提供的代码2压缩包中的Test_gOMP.m 。plain view plaincopy%压缩感知重构算法测试clear all;close all;clc;M = 128;%观测值个数N = 256;% 信号x的长度K = 30

9、;% 信号x的稀疏度Index_K = randperm(N);x = zeros(N,1);x(Index_K(1:K) = 5*randn(K,1);%x为 K稀疏的，且位置是随机的Psi = eye(N);%x本身是稀疏的，定义稀疏矩阵为单位阵x=Psi*thetaPhi = randn(M,N)/sqrt(M);%测量矩阵为高斯矩阵A = Phi * Psi;% 传感矩阵y = Phi * x;% 得到观测向量 y%恢复重构信号xtictheta = CS_gOMP( y,A,K);x_r = Psi * theta;% x=Psi * thetatoc%绘图figure;plot(x

10、_r,k.-);% 绘出x的恢复信号hold on;plot(x,r);%绘出原信号 xhold off;legend(Recovery,Original)fprintf(n恢复残差：);norm(x_r-x)% 恢复残差运行结果如下：(信号为随机生成，所以每次结果均不一样)1 )图:1图:Command windowsElapsedtime is 0.155937 seconds.恢复残差：ans=2.3426e-0144、信号稀疏度K与重构成功概率关系曲线绘制例程代码以下测试代码为了与文献1的Fig.1作比较。由于暂未研究学习LP算法，所以相比于文献1的Fig.1)缺少LP算法曲线，加入了

11、 SP算法。以下测试代码与 SAMP相应的测试代码基本一致，可以合并在一起运行，只须在 主循环内多加几种算法重构就行。plain view plaincopy c%E 缩感知重构算法测试 CS_Reconstuction_KtoPercentagegOMP.m% 绘制参考文献中的 Fig.1% Reference: Jian Wang, Seokbeop Kwon, Byonghyo Shim. Generalized% orthogonal matching pursuit, IEEE Transactions on Signal Processing,% vol. 60, no. 12,

12、pp. 6202-6216, Dec. 2012.% Available at: HYPERLINK http:/islab.snu.ac.kr/paper/tsp_gOMP.pdf http:/islab.snu.ac.kr/paper/tsp_gOMP.pdf.% Elapsed time is 798.718246 seconds.(20150509pm)clear all;close all;clc;%参数配置初始化CNT = 1000;% 对于每组(K,M,N),重复迭代次数N = 256;% 信号x的长度Psi = eye(N);%x本身是稀疏的，定义稀疏矩阵为单位阵x=Psi*t

13、hetaM_set = 128;%测量值集合KIND = OMP ;ROMP ;StOMP ;SP ;CoSaMP ;.gOMP(s=3);gOMP(s=6);gOMP(s=9);Percentage = zeros(N,length(M_set),size(KIND,1);%存储恢复成功概率%主循环，遍历每组(K,M,N)ticfor mm = 1:length(M_set)M = M_set(mm);%本次测量值个数K_set = 5:5:70;% 信号x的稀疏度K没必要全部遍历，每隔5测试一个就可以了%存储此测量值M下不同K的恢复成功概率PercentageM = zeros(size(

14、KIND,1),length(K_set);for kk = 1:length(K_set)K = K_set(kk);%本次信号x的稀疏度KP = zeros(1,size(KIND,1);fprintf(M=%d,K=%dn,M,K);for cnt = 1:CNT %每个观测值个数均运行CNT次Index_K = randperm(N);x = zeros(N,1);x(Index_K(1:K) = 5*randn(K,1);%x为K稀疏的，且位置是随机的Phi = randn(M,N)/sqrt(M);%测量矩阵为高斯矩阵A = Phi * Psi;%传感矩阵y = Phi * x;%

15、得到观测向量 y%(1)OMPtheta = CS_OMP(y,A,K);%x_r = Psi * theta;% x=Psi * thetaif norm(x_r-x)1e-6%P(1) = P(1) + 1;end%(2)ROMPtheta = CS_ROMP(y,A,K);%x_r = Psi * theta;% x=Psi * thetaif norm(x_r-x)1e-6%P(2) = P(2) + 1;end%(3)StOMPtheta = CS_StOMP(y,A);%x_r = Psi * theta;% x=Psi * thetaif norm(x_r-x)1e-6%P(3)

16、 = P(3) + 1;end%(4)SPtheta = CS_SP(y,A,K);%x_r = Psi * theta;% x=Psi * theta恢复重构信号theta如果残差小于1e-6则认为恢复成功恢复重构信号theta如果残差小于1e-6则认为恢复成功恢复重构信号theta如果残差小于1e-6则认为恢复成功恢复重构信号theta如果残差小于1e-6则认为恢复成功恢复重构信号theta如果残差小于1e-6则认为恢复成功恢复重构信号theta如果残差小于1e-6则认为恢复成功恢复重构信号theta如果残差小于1e-6则认为恢复成功恢复重构信号theta如果残差小于1e-6则认为恢复成功

17、PercentageM(iii,kk) = P(iii)/CNT*100;%计算恢复概率end end for jjj = 1:size(KIND,1)Percentage(1:length(K_set),mm,jjj) = PercentageM(jjj,:);endend toc save KtoPercentage1000gOMP %运行一次不容易，把变量全部存储下来if norm(x_r-x)1e-6%P(4) = P(4) + 1;end%(5)CoSaMPtheta = CS_CoSaMP(y,A,K);%x_r = Psi * theta;% x=Psi* thetaif nor

18、m(x_r-x)1e-6%P(5) = P(5) + 1;end%(6)gOMP,S=3theta = CS_gOMP(y,A,K,3);%x_r = Psi * theta;% x=Psi* thetaif norm(x_r-x)1e-6%P(6) = P(6) + 1;end%(7)gOMP,S=6theta = CS_gOMP(y,A,K,6);%x_r = Psi * theta;% x=Psi* thetaif norm(x_r-x)1e-6%P(7) = P(7) + 1;end%(8)gOMP,S=9theta = CS_gOMP(y,A,K,9);%x_r = Psi * th

19、eta;% x=Psi* thetaif norm(x_r-x)1e-6%P(8) = P(8) + 1;endendfor iii = 1:size(KIND,1)%绘图S = -ks;-ko;-yd;-gv;-b*;-r.;-rx;-r+;figure;for mm = 1:length(M_set)M = M_set(mm);K_set = 5:5:70;L_Kset = length(K_set);for ii = 1:size(KIND,1)plot(K_set,Percentage(1:L_Kset,mm,ii),S(ii,:);%绘出 x 的恢复信号hold on;endend1

20、009.hold off;xlim(5 70);legend(OMP,ROMP,StOMP,SP,CoSaMP,.gOMP(s=3),gOMP(s=6),gOMP(s=9);xlabel(Sparsity level K);ylabel(The Probability of Exact Reconstruction);title(Prob. of exact recovery vs. the signal sparsity K(M=128,N=256)(Gaussian);本程序在联想 ThinkPadE430C 笔记本（4GB DDR3内存，i5-3210 ）上运行 共

21、耗时798.718246秒，程序中 将所有数据均通过save KtoPercentage1000gOMP”存储了下来，以后可以再对数据进行分析，只需 “loadKtoPercentage1000gOMP ” 即可。本程序运行结果：ProL Mod耻t recovery vs. the signal sparsity K(M=128=256(Gaussian!Sparsity level K文献口中的Fig 1：二OUF1N 6u口口clfr0EK3 一0 AgunbsJLSparsityFig . 1. RectumLmcticn perfurmuiwe for L -sparse Gauss

22、ian signal vector as a Hinciion of sparsity A.5、结语我很好奇：为什么相比于 OMP算法就是简单每次多选几列，重构效果为什么这么好？居然比复杂的ROMP、CoSaMP、StOMP效果还要好该课题组还提出了MMP算法，可参见文献3。更多关于该课题组的信息可去官方网站查询： HYPERLINK http:/islab.snu.ac.kr/ http:/islab.snu.ac.kr/ ,也可直接查看发表的文章： HYPERLINK http:/islab.snu.ac.kr/publication.html http:/islab.snu.ac.kr/

23、publication.html 。文献1最后有两个TABLE ,分别是算法的流程和复杂度总结:TABLE I-DELETLtrOMP AlgorithmInput; mcusurcmcnu y G R:*nxing matrix 小 亮皿 乂 spars-ity K“number of indices for each selection N (N K and N .InilialLz4:iteration wuiil k 0*residual, vector r = yT EstinuikiJ -ppon ssl A = *While l|N g f and * minK, mN do f

24、c = Ar 4 1.ifdendtication) Select indices (i)h=i/2一一，,v corresjxjnding to jV largest crtTries (in mag力inide in.(AumentatioD) A* = A* U 风1，！雄(). (Estimadon of xAt xAt = arg nijn y 中a依】丁 Ikaidual Update) r1 = y EndOutput x = arg min |y - u|2. U：Kupp U)二h菇TABLE 11COMflEXJTY OFTHt CjOM * ALliORITHM (A-TH S