题目：压缩感知重构算法之正则化正交匹配追踪(ROMP)
正交匹配追踪算法每次迭代均只选择与残差最相关的一列，自然人们会想：“每次迭代是否可以多选几列呢?”，正则化正交匹配追踪(RegularizedOMP)就是其中一种改进方法。本篇将在上一篇《压缩感知重构算法之正交匹配追踪(OMP)》的基础上给出正则化正交匹配追踪(ROMP)算法的MATLAB函数代码，并且给出单次测试例程代码、测量数M与重构成功概率关系曲线绘制例程代码。
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称为传感矩阵。
1、ROMP重构算法流程：
正则化正交匹配追踪算法流程与OMP的最大不同之处就在于从传感矩阵A中选择列向量的标准，OMP每次只选择与残差内积绝对值最大的那一列，而ROMP则是先选出内积绝对值最大的K列（若所有内积中不够K个非零值则将内积值非零的列全部选出），然后再从这K列中按正则化标准再选择一遍，，即为本次迭代选出的列向量（一般并非只有一列）。正则化标准意思是选择各列向量与残差内积绝对值的最大值不能比最小值大两倍以上(comparable coordinates)且能量最大的一组(with the maximal energy)，因为满足条件的子集并非只有一组。似乎用叙述语言描述不清楚，下面给出一种实现第(2)(3)步的算法流程图（此算法并非本人原创，参考网络代码[2][3]，本人将代码中的思想进行整理，画出此流程图，方便初学者快速掌握学习ROMP算法）：
我将原子选择过程封装成了一个MATLAB函数，代码如下(Regularize.m)：
function [val,pos] = Regularize(product,Kin)%Regularize Summary of this function goes here% Detailed explanation goes here% product = A'*r_n;%传感矩阵A各列与残差的内积% K为稀疏度% pos为选出的各列序号% val为选出的各列与残差的内积值% Reference:Needell D，Vershynin R. Uniform uncertainty principle and% signal recovery via regularized orthogonal matching pursuit. % Foundations of Computational Mathematics, 2009,9(3): 317-334.productabs = abs(product);%取绝对值[productdes,indexproductdes] = sort(productabs,'descend');%降序排列for ii = length(productdes):-1:1if productdes(ii)>1e-6%判断productdes中非零值个数break;endend%Identify:Choose a set J of the K biggest coordinatesif ii>=KinJ = indexproductdes(1:Kin);%集合JJval = productdes(1:Kin);%集合J对应的序列值K = Kin;else%or all of its nonzero coordinates,whichever is smallerJ = indexproductdes(1:ii);%集合JJval = productdes(1:ii);%集合J对应的序列值K = ii;end%Regularize:Among all subsets J0∈J with comparable coordinatesMaxE = -1;%循环过程中存储最大能量值for kk = 1:KJ0_tmp = zeros(1,K);iJ0 = 1;J0_tmp(iJ0) = J(kk);%以J(kk)为本次寻找J0的基准(最大值)Energy = Jval(kk)^2;%本次寻找J0的能量for mm = kk+1:Kif Jval(kk)<2*Jval(mm)%找到符合|u(i)|<=2|u(j)|的iJ0 = iJ0 + 1;%J0自变量增1J0_tmp(iJ0) = J(mm);%更新J0Energy = Energy + Jval(mm)^2;%更新能量else%不符合|u(i)|<=2|u(j)|的break;%跳出本轮寻找，因为后面更小的值也不会符合要求endendif Energy>MaxE%本次所得J0的能量大于前一组J0 = J0_tmp(1:iJ0);%更新J0MaxE = Energy;%更新MaxE，为下次循环做准备endendpos = J0;val = productabs(J0);end
2、正则化正交匹配追踪(ROMP)MATLAB代码(CS_ROMP.m)
这个函数是完全基于上一篇中的CS_OMP.m修改而成的。
function [ theta ] = CS_ROMP( y,A,K )%CS_ROMP Summary of this function goes here%Version: 1.0 written by jbb0523 @2015-04-24% Detailed explanation goes here% y = Phi * x% x = Psi * theta%y = Phi*Psi * theta% 令 A = Phi*Psi, 则y=A*theta% 现在已知y和A，求theta% Reference:Needell D，Vershynin R．Signal recovery from incomplete and% inaccurate measurements via regularized orthogonal matching pursuit[J]．% IEEE Journal on Selected Topics in Signal Processing，2010，4(2)：310—316.[y_rows,y_columns] = size(y);if y_rows<y_columnsy = y';%y should be a column vectorend[M,N] = size(A);%传感矩阵A为M*N矩阵theta = zeros(N,1);%用来存储恢复的theta(列向量)At = zeros(M,3*K);%用来迭代过程中存储A被选择的列Pos_theta = zeros(1,2*K);%用来迭代过程中存储A被选择的列序号Index = 0;r_n = y;%初始化残差(residual)为y%Repeat the following steps K times(or until |I|>=2K)for ii=1:K%迭代K次product = A'*r_n;%传感矩阵A各列与残差的内积%[val,pos] = max(abs(product));%找到最大内积绝对值，即与残差最相关的列[val,pos] = Regularize(product,K);%按正则化规则选择原子At(:,Index+1:Index+length(pos)) = A(:,pos);%存储这几列Pos_theta(Index+1:Index+length(pos)) = pos;%存储这几列的序号if Index+length(pos)<=M%At的行数大于列数，此为最小二乘的基础(列线性无关)Index = Index+length(pos);%更新Index，为下次循环做准备else%At的列数大于行数，列必为线性相关的,At(:,1:Index)'*At(:,1:Index)将不可逆break;%跳出for循环endA(:,pos) = zeros(M,length(pos));%清零A的这几列(其实此行可以不要因为它们与残差正交)%y=At(:,1:Index)*theta，以下求theta的最小二乘解(Least Square)theta_ls = (At(:,1:Index)'*At(:,1:Index))^(-1)*At(:,1:Index)'*y;%最小二乘解%At(:,1:Index)*theta_ls是y在At(:,1:Index)列空间上的正交投影r_n = y – At(:,1:Index)*theta_ls;%更新残差if norm(r_n)<1e-6%Repeat the steps until r=0break;%跳出for循环endif Index>=2*K%or until |I|>=2Kbreak;%跳出for循环endend
3、ROMP单次重构测试代码
的这一半更多地赢取上帝掌握的那一半。
