Matlab实现压缩感知重建算法：ROMP正则化的正交匹配追踪算法

ROMP正则化的正交匹配追踪算法

$ROMP$算法是一种贪婪的压缩感知重建算法，在$OMP$的基础上加入了正则化过程。此外，相比于$OMP$算法，$ROMP$在每次迭代过程中选择与残差最相关的$m$个列向量。下面是$DeannaNeedell$和$RomanVershynin$的文章中给出的$ROMP$的算法。

在正则化过程，需要寻找索引集合的所有子集，在这里可以使用$Matlab$中的$ff2n$函数实现。对于正则化条件，只用满足$u$中的最大元素绝对值小于最小元素绝对值的$2$倍即可。
本文中用高斯随机矩阵作为测量矩阵，待观测目标信号长度$N$为$512$，且元素值只有$-1,0，1$，观测值个数$M$为$120$，目标信号稀疏度$K$为$10$。$Matlab$实验结果如下：

clear;
clc;
N=512;%待观测目标信号长度
K=10;%目标信号稀疏度
M=120;%观测值个数
maxnum=100;
x=zeros(N,1);
q=randperm(N);%随机打乱排序
t=1:N;
x(q(1:K))=sign(randn(K,1));%生成幅值为1，-1的K个非0元素的目标信号y
A=randn(M,N);%生成M*N的服从高斯分布的测量矩阵
A=orth(A')';%把测量矩阵按行正交化
y=A*x;%获得观测向量
x2=A'*y;%基于L2范数最小化估算的初值
S=ff2n(K);%20个元素所有可能的组合,(2^K)*K的矩阵
I=[];
A1=[];
r=y;%初始化时，残差r=y
for i=1:K
    u=A'*r;
    u1=abs(u);
    [temp,index]=sort(u1);%对向量u从小到大排序,index储存索引
    len1=length(u);
    J=index(len1-K+1:end);%储存K个的最大系数索引,此处存在问题，应该是非0
   if length(I)>=2*K
       break;
   else
       %正则化标准意思是选择各列向量与残差内积绝对值的最大值不能比最小值大两倍以上
       %(comparable coordinates)且能量最大的一组(with the maximal energy)
       [row,col]=size(S);%记录S的行数
       uu=[0;u];
       index1=S.*J';%存储所有子集在u中对应的索引
       temp2=uu(index1+1);%存储所有子集中的元素值，(2^K)个子集，一行看成一个向量
       abstemp=abs(temp2);
       max1=max(abstemp');
       abstemp1=abstemp';
       %abstemp1(find(abstemp1==0))=maxnum;
       min1=min(abstemp1(find(abstemp1~=0)));
       min1(1)=0;
       diff=abs(max1)-2*abs(min1);%最大值大于最小值两倍以上，则>0，
       %2倍关系可以换成其他的条件吗？
       index2=find(diff>0);%寻找满足最大值小于最小值两倍以上的索引
       if isempty(find(diff>0,1))
          disp('未发现满足正则条件的集合');
          break;
       end
       temp3=temp2(index2,:);%保留满足条件的子集
       temp4=temp3.*temp3;%中间变量，存储元素平方
       temp5=sqrt(sum(temp4,2));%按行求和，再开方，即求二范数
       [maxnorm,index3]=max(temp5);%会不会存在多个最大值？
       temp6=temp3(index3(1),:);%存储第一个最大范数的坐标向量
       temp7=temp2-temp6;%若temp6是temp2中某一行，则那一行全为0
       temp8=temp7;
       temp8(1,:)=[];
       J0=find(all(temp8== 0,2));%all(temp7== 0,2)寻找全是0的行 
       J0=index1(J0+1,:);
       if isempty(intersect(I,J0))~=0
           T=intersect(I,J0);%查找重复元素,并删除
           for jj=1:length(T)
               T1=find(J0==T(jj));
               J0(T1)=[];
           end
       end
       I=[I J0];%更新索引集合
       A1=A(:,I);
       y2=inv(A1'*A1)*A1'*y;%最小二乘法求y2
       r=y-A1*y2;
   end
end
y3=zeros(N,1);
y3(I)=y2;%储存恢复结果
plot(t,x,'r',t,y3,'g  *');
legend('red 代表原始信号','green 代表恢复的信号');
length(find(y3~=0));%在能正确恢复信号时，y3中非0系数是2K
% A1(:,1)'*A1(:,10)
% A(:,1)'*A(:,200)

Needell D , Vershynin R . Uniform Uncertainty Principle and signal recovery via Regularized Orthogonal Matching Pursuit[J]. 2007.