贪婪算法中ROMP算法的原理介绍及MATLAB仿真

1、压缩感知重构算法之正则化正交匹配追踪(ROMP)正交匹配追踪算法每次迭代均只选择与残差最相关的一列，自然人们会想：“每次迭代是否可以多选几列呢?”，正则化正交匹配追踪(RegularizedOMP)就是其中一种改进方法。本篇将在上一篇压缩感知重构算法之正交匹配追踪(OMP)的基础上给出正则化正交匹配追踪(ROMP)算法的MATLAB函数代码，并且给出单次测试例程代码、测量数M与重构成功概率关系曲线绘制例程代码。0、符号说明如下：        压缩观测y=x，其中y为观测所得向量M×1，x为原信号N×1（M<<N

2、）。x一般不是稀疏的，但在某个变换域是稀疏的，即x=，其中为K稀疏的，即只有K个非零项。此时y=，令A=，则y=A。        (1) y为观测所得向量，大小为M×1        (2)x为原信号，大小为N×1        (3)为K稀疏的，是信号在x在某变换域的稀疏表示        (4)称为观测矩阵、测量矩阵、测量基，大小为M×N 

3、      (5)称为变换矩阵、变换基、稀疏矩阵、稀疏基、正交基字典矩阵，大小为N×N        (6)A称为测度矩阵、传感矩阵、CS信息算子，大小为M×N上式中，一般有K<<M<<N，后面三个矩阵各个文献的叫法不一，以后我将称为测量矩阵、将称为稀疏矩阵、将A称为传感矩阵。1、ROMP重构算法流程：正则化正交匹配追踪（ROMP）算法流程与OMP的最大不同之处：就在于从传感矩阵A中选择列向量的标准，OMP每次只选择与残差内积绝对值最大的那一列，而ROMP则是先

4、选出内积绝对值最大的K列（若所有内积中不够K个非零值则将内积值非零的列全部选出），然后再从这K列中按正则化标准再选择一遍，即为本次迭代选出的列向量（一般并非只有一列）。正则化标准：意思是选择各列向量与残差的内积的绝对值的最大值不能比最小值大两倍以上(comparable coordinates)且能量最大的一组(with the maximal energy)，因为满足条件的子集并非只有一组。似乎用叙述语言描述不清楚，下面给出一种实现第(2)(3)步的算法流程图（此算法并非本人原创，参考网络代码23，本人将代码中的思想进行整理，画出此流程图，方便初学者快速掌握学习ROMP算法）：我将原子选择过

5、程封装成了一个MATLAB函数，代码如下(Regularize.m)：1. function val,pos = Regularize(product,Kin)  2. %Regularize Summary of this function goes here  3. %   Detailed explanation goes here  4. %   p

6、roduct = A'*r_n;%传感矩阵A各列与残差的内积  5. %   K为稀疏度  6. %   pos为选出的各列序号  7. %   val为选出的各列与残差的内积值  8. %   Reference:Needell D，Vershynin R. Uniform uncertainty principle&

7、#160;and  9. %   signal recovery via regularized orthogonal matching pursuit.   10. %   Foundations of Computational Mathematics, 2009,9(3): 317-334.    11.   &

8、#160; productabs = abs(product);%取绝对值  12.     productdes,indexproductdes = sort(productabs,'descend');%降序排列  13.     for ii = length(productdes):-1:1  14.     &#

9、160;   if productdes(ii)>1e-6%判断productdes中非零值个数  15.             break;  16.         end  17.     end  18.  

10、60;  %Identify:Choose a set J of the K biggest coordinates  19.     if ii>=Kin  20.         J = indexproductdes(1:Kin);%集合J  21.   &

11、#160;     Jval = productdes(1:Kin);%集合J对应的序列值  22.         K = Kin;  23.     else%or all of its nonzero coordinates,whichever is smaller

12、0; 24.         J = indexproductdes(1:ii);%集合J  25.         Jval = productdes(1:ii);%集合J对应的序列值  26.         K = ii;  

13、;27.     end  28.     %Regularize:Among all subsets J0J with comparable coordinates  29.     MaxE = -1;%循环过程中存储最大能量值  30.     for kk = 

14、1:K  31.         J0_tmp = zeros(1,K);iJ0 = 1;  32.         J0_tmp(iJ0) = J(kk);%以J(kk)为本次寻找J0的基准(最大值)  33.         Energ

15、y = Jval(kk)2;%本次寻找J0的能量  34.         for mm = kk+1:K  35.             if Jval(kk)<2*Jval(mm)%找到符合|u(i)|<=2|u(j)|的  36.   &#

16、160;             iJ0 = iJ0 + 1;%J0自变量增1  37.                 J0_tmp(iJ0) = J(mm);%更新J0  38.   &

17、#160;             Energy = Energy + Jval(mm)2;%更新能量  39.             else%不符合|u(i)|<=2|u(j)|的  40.      

18、60;          break;%跳出本轮寻找，因为后面更小的值也不会符合要求  41.             end  42.         end  43.      

19、0;  if Energy>MaxE%本次所得J0的能量大于前一组  44.             J0 = J0_tmp(1:iJ0);%更新J0  45.             MaxE = Energy;%更新MaxE，为下次循环做准

20、备  46.         end  47.     end  48.     pos = J0;  49.     val = productabs(J0);  50. end  2、正则化正交匹配追踪(ROMP)MATLAB代码(CS

21、_ROMP.m)这个函数是完全基于上一篇中的CS_OMP.m修改而成的。1. function  theta  = CS_ROMP( y,A,K )  2. %CS_ROMP Summary of this function goes here  3. %Version: 1.0 written by jbb0523 2015-04-24  4. %

22、0;  Detailed explanation goes here  5. %   y = Phi * x  6. %   x = Psi * theta  7. %   y = Phi*Psi * theta  8. %   

23、令 A = Phi*Psi, 则y=A*theta  9. %   现在已知y和A，求theta  10. %   Reference:Needell D，Vershynin RSignal recovery from incomplete and  11. %   inaccurate measurements via

24、0;regularized orthogonal matching pursuitJ  12. %   IEEE Journal on Selected Topics in Signal Processing，2010，4(2)：310316.  13.     y_rows,y_columns = size(y);  14.  

25、0;  if y_rows<y_columns  15.         y = y'%y should be a column vector  16.     end  17.     M,N = size(A);%传感矩阵A为M*N矩阵

26、60; 18.     theta = zeros(N,1);%用来存储恢复的theta(列向量)  19.     At = zeros(M,3*K);%用来迭代过程中存储A被选择的列  20.     Pos_theta = zeros(1,2*K);%用来迭代过程中存储A被选择的列序号  21.     

27、;Index = 0;  22.     r_n = y;%初始化残差(residual)为y  23.     %Repeat the following steps K times(or until |I|>=2K)  24.     for ii=1:K%迭代K次  

28、25.         product = A'*r_n;%传感矩阵A各列与残差的内积  26.         %val,pos = max(abs(product);%找到最大内积绝对值，即与残差最相关的列  27.         val,pos =

29、60;Regularize(product,K);%按正则化规则选择原子  28.         At(:,Index+1:Index+length(pos) = A(:,pos);%存储这几列  29.         Pos_theta(Index+1:Index+length(pos) = pos;%存储这几列的序号  30. &

30、#160;       if Index+length(pos)<=M%At的行数大于列数，此为最小二乘的基础(列线性无关)  31.             Index = Index+length(pos);%更新Index，为下次循环做准备  32.       &#

31、160; else%At的列数大于行数，列必为线性相关的,At(:,1:Index)'*At(:,1:Index)将不可逆  33.             break;%跳出for循环  34.         end  35.        

32、60;A(:,pos) = zeros(M,length(pos);%清零A的这几列(其实此行可以不要因为它们与残差正交)  36.         %y=At(:,1:Index)*theta，以下求theta的最小二乘解(Least Square)  37.         theta_ls = (At(:,1:Index)'*At(:

33、,1:Index)(-1)*At(:,1:Index)'*y;%最小二乘解  38.         %At(:,1:Index)*theta_ls是y在At(:,1:Index)列空间上的正交投影  39.         r_n = y - At(:,1:Index)*theta_ls;%更新残差  40.  

34、60;      if norm(r_n)<1e-6%Repeat the steps until r=0  41.             break;%跳出for循环  42.         end  43. 

35、0;       if Index>=2*K%or until |I|>=2K  44.             break;%跳出for循环  45.         end  46.    &#

36、160;end  47.     theta(Pos_theta(1:Index)=theta_ls;%恢复出的theta  48. end  3、ROMP单次重构测试代码以下测试代码与上一篇中的OMP单次测试代码基本完全一致，除了M和K参数设置及调用CS_ROMP函数三处之外。1. %压缩感知重构算法测试  2. clear all;close all;clc;  3. M = 128;%观测值个数 &#

37、160;4. N = 256;%信号x的长度  5. K = 12;%信号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

38、60; 10. Phi = randn(M,N);%测量矩阵为高斯矩阵  11. A = Phi * Psi;%传感矩阵  12. y = Phi * x;%得到观测向量y  13. % 恢复重构信号x  14. tic  15. theta = CS_ROMP(y,A,K);  16. x_r = Psi 

39、* 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

40、9;,'Original')  25. fprintf('n恢复残差：');  26. norm(x_r-x)%恢复残差  运行结果如下：（信号为随机生成，所以每次结果均不一样）1）图：2）CommandwindowsElapsed time （运行时间）is 0.022132 seconds.恢复残差：ans= 7.8066e-0154、测量数M与重构成功概率关系曲线绘制例程代码以下测试代码与上一篇中的OMP测量数M与重构成功概率关系曲线绘制例程代码基本完全一致。本程序在循环中填加了“kk”一行代码并

41、将“M = M_set(mm)”一行的分号去掉，这是为了在运行过程中可以观察程序运行状态、知道程序到哪一个位置。重构调用CS_ROMP函数并将save命令的文件名改为ROMPMtoPercentage1000，以防止将OMP存储的数据覆盖（因为本人的所有代码都在一个目录下）。1. clear all;close all;clc;  2. % 参数配置初始化  3. CNT = 1000;%对于每组(K,M,N)，重复迭代次数  4. N = 256;%信号x的长度&

42、#160; 5. Psi = eye(N);%x本身是稀疏的，定义稀疏矩阵为单位阵x=Psi*theta  6. K_set = 4,12,20,28,36;%信号x的稀疏度集合  7. Percentage = zeros(length(K_set),N);%存储恢复成功概率  8. % 主循环，遍历每组(K,M,N)  9. tic  10. for kk = 1:length(K_

43、set)  11.     K = K_set(kk);%本次稀疏度  12.     M_set = K:5:N;%M没必要全部遍历，每隔5测试一个就可以了  13.     PercentageK = zeros(1,length(M_set);%存储此稀疏度K下不同M的恢复成功概率  14.    &#

44、160;kk  15.     for mm = 1:length(M_set)  16.        M = M_set(mm)%本次观测值个数  17.        P = 0;  18.       

45、 for cnt = 1:CNT %每个观测值个数均运行CNT次  19.             Index_K = randperm(N);  20.             x = zeros(N,1); 

46、0;21.             x(Index_K(1:K) = 5*randn(K,1);%x为K稀疏的，且位置是随机的                  22.          &

47、#160;  Phi = randn(M,N);%测量矩阵为高斯矩阵  23.             A = Phi * Psi;%传感矩阵  24.             y = Phi * x

48、;%得到观测向量y  25.             theta = CS_ROMP(y,A,K);%恢复重构信号theta  26.             x_r = Psi * theta;% x=Psi * theta&#

49、160; 27.             if norm(x_r-x)<1e-6%如果残差小于1e-6则认为恢复成功  28.                 P = P + 1;  29.   

50、;          end  30.        end  31.        PercentageK(mm) = P/CNT*100;%计算恢复概率  32.     end  33.   &

51、#160; Percentage(kk,1:length(M_set) = PercentageK;  34. end  35. toc  36. save ROMPMtoPercentage1000 %运行一次不容易，把变量全部存储下来  37. % 绘图  38. S = '-ks''-ko''-kd''-kv''-k*' 

52、 39. figure;  40. for kk = 1:length(K_set)  41.     K = K_set(kk);  42.     M_set = K:5:N;  43.     L_Mset = length(M_set);  44.     plot(M_set,Percentage(kk,1:L_Mset),S(kk,:);%绘出x的恢复信号  45.