压缩感知重构算法之广义正交匹配追踪(gOMP)

转载自彬彬有礼的专栏

题目：压缩感知重构算法之广义正交匹配追踪(gOMP)

        广义正交匹配追踪(Generalized OMP, gOMP)算法可以看作为OMP算法的一种推广，由文献[1]提出，第1作者本硕为哈工大毕业，发表此论文时在Korea University攻读博士学位。OMP每次只选择与残差相关最大的一个，而gOMP则是简单地选择最大的S个。之所以这里表述为“简单地选择”是相比于ROMP之类算法的，不进行任何其它处理，只是选择最大的S个而已。

0、符号说明如下：

        压缩观测y=Φx，其中y为观测所得向量M×1，x为原信号N×1（M<<N）。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<<M<<N，后面三个矩阵各个文献的叫法不一，以后我将Φ称为测量矩阵、将Ψ称为稀疏矩阵、将A称为传感矩阵。

        注意：这里的稀疏表示模型为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] view plain copy

print?

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<y_columns  
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  

function [ 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 = Phi*Psi * theta%   令 A = Phi*Psi, 则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: http://islab.snu.ac.kr/paper/tsp_gOMP.pdfif nargin < 4S = round(max(K/4, 1));end    [y_rows,y_columns] = size(y);    if y_rows<y_columns        y = y';%y should be a column vector    end    [M,N] = size(A);%传感矩阵A为M*N矩阵    theta = zeros(N,1);%用来存储恢复的theta(列向量)    Pos_theta = [];%用来迭代过程中存储A被选择的列序号    r_n = y;%初始化残差(residual)为y    for ii=1:K%迭代K次，K为稀疏度        product = A'*r_n;%传感矩阵A各列与残差的内积        [val,pos]=sort(abs(product),'descend');%降序排列        Sk = union(Pos_theta,pos(1:S));%选出最大的S个        if length(Sk)==length(Pos_theta)            if ii == 1                theta_ls = 0;            end            break;        end        if length(Sk)>M            if ii == 1                theta_ls = 0;            end            break;        end        At = A(:,Sk);%将A的这几列组成矩阵At        %y=At*theta，以下求theta的最小二乘解(Least Square)        theta_ls = (At'*At)^(-1)*At'*y;%最小二乘解        %At*theta_ls是y在At)列空间上的正交投影        r_n = y - At*theta_ls;%更新残差        Pos_theta = Sk;        if norm(r_n)<1e-6            break;%quit the iteration        end    end    theta(Pos_theta)=theta_ls;%恢复出的thetaend

3、gOMP单次重构测试代码(CS_Reconstuction_Test.m)

        以下测试代码基本与OMP单次重构测试代码一样。也可参考该实验室网站提供的代码^[2]压缩包中的Test_gOMP.m。

[plain] view plain copy

print?

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)%恢复残差  

%压缩感知重构算法测试clear all;close all;clc;M = 128;%观测值个数N = 256;%信号x的长度K = 30;%信号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_r,'k.-');%绘出x的恢复信号hold on;plot(x,'r');%绘出原信号xhold off;legend('Recovery','Original')fprintf('\n恢复残差：');norm(x_r-x)%恢复残差

        运行结果如下：（信号为随机生成，所以每次结果均不一样）

        1） 图：

        2）Command  windows

        Elapsedtime is 0.155937 seconds.

        恢复残差：

        ans=

          2.3426e-014

4、信号稀疏度K与重构成功概率关系曲线绘制例程代码

        以下测试代码为了与文献[1]的Fig.1作比较。由于暂未研究学习LP算法，所以相比于文献[1]的Fig.1)缺少LP算法曲线，加入了SP算法。以下测试代码与SAMP相应的测试代码基本一致，可以合并在一起运行，只须在主循环内多加几种算法重构就行。

[plain] view plain copy

print?

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)');  

%压缩感知重构算法测试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, pp. 6202-6216, Dec. 2012. %   Available at: 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*thetaM_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(KIND,1),length(K_set));    for kk = 1:length(K_set)       K = K_set(kk);%本次信号x的稀疏度K       P = zeros(1,size(KIND,1));       fprintf('M=%d,K=%d\n',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;%得到观测向量y            %(1)OMP            theta = CS_OMP(y,A,K);%恢复重构信号theta            x_r = Psi * theta;% x=Psi * theta            if norm(x_r-x)<1e-6%如果残差小于1e-6则认为恢复成功                P(1) = P(1) + 1;            end            %(2)ROMP            theta = CS_ROMP(y,A,K);%恢复重构信号theta            x_r = Psi * theta;% x=Psi * theta            if norm(x_r-x)<1e-6%如果残差小于1e-6则认为恢复成功                P(2) = P(2) + 1;            end            %(3)StOMP            theta = CS_StOMP(y,A);%恢复重构信号theta            x_r = Psi * theta;% x=Psi * theta            if norm(x_r-x)<1e-6%如果残差小于1e-6则认为恢复成功                P(3) = P(3) + 1;            end            %(4)SP            theta = CS_SP(y,A,K);%恢复重构信号theta            x_r = Psi * theta;% x=Psi * theta            if norm(x_r-x)<1e-6%如果残差小于1e-6则认为恢复成功                P(4) = P(4) + 1;            end            %(5)CoSaMP            theta = CS_CoSaMP(y,A,K);%恢复重构信号theta            x_r = Psi * theta;% x=Psi * theta            if norm(x_r-x)<1e-6%如果残差小于1e-6则认为恢复成功                P(5) = P(5) + 1;            end            %(6)gOMP,S=3            theta = CS_gOMP(y,A,K,3);%恢复重构信号theta            x_r = Psi * theta;% x=Psi * theta            if norm(x_r-x)<1e-6%如果残差小于1e-6则认为恢复成功                P(6) = P(6) + 1;            end            %(7)gOMP,S=6            theta = CS_gOMP(y,A,K,6);%恢复重构信号theta            x_r = Psi * theta;% x=Psi * theta            if norm(x_r-x)<1e-6%如果残差小于1e-6则认为恢复成功                P(7) = P(7) + 1;            end            %(8)gOMP,S=9            theta = CS_gOMP(y,A,K,9);%恢复重构信号theta            x_r = Psi * theta;% x=Psi * theta            if norm(x_r-x)<1e-6%如果残差小于1e-6则认为恢复成功                P(8) = P(8) + 1;            end       end       for iii = 1:size(KIND,1)           PercentageM(iii,kk) = P(iii)/CNT*100;%计算恢复概率       end    end    for jjj = 1:size(KIND,1)        Percentage(1:length(K_set),mm,jjj) = PercentageM(jjj,:);    endendtocsave KtoPercentage1000gOMP %运行一次不容易，把变量全部存储下来%% 绘图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;    endendhold 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）上运行共耗时798.718246秒，程序中将所有数据均通过“save KtoPercentage1000gOMP”存储了下来，以后可以再对数据进行分析，只需“load KtoPercentage1000gOMP”即可。

        本程序运行结果：

        文献[1]中的Fig.1：

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, Seokbeop 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

  相关文章推荐

• • 压缩感知重构算法之压缩采样匹配追踪(CoSaMP)
• • 深度学习部署系统构建--刘文志
• • 压缩感知重构算法之子空间追踪(SP)
• • 搜狗机器翻译技术分享--陈伟
• • matlab绘制多曲线对比图
• • Hadoop生态系统零基础入门
• • 使用l1-magic工具箱求解基追踪(BP)和基追踪降噪(BPDN)
• • 最懂程序员的学习方式 TensorFlow入门

• • 压缩感知稀疏基之离散小波变换
• • Retrofit 从入门封装到源码解析
• • 压缩感知重构算法之Gradient Descent with Sparsification(GraDeS)
• • 程序员如何转型AI工程师--蒋涛
• • 压缩感知重构算法之迭代硬阈值(IHT)
• • 压缩感知重构算法之稀疏度自适应匹配追踪(SAMP)
• • 为什么正交匹配追踪(OMP)一定能恢复信号？
• • 压缩感知正交匹配追踪重构算法