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

压缩感知重构算法之正则化正交匹配追踪(ROMP)

正交匹配追踪算法每次迭代均只选择与残差最相关的一列，自然人们会想：“每次迭代是否可以多选几列呢”，正则化正交匹配追踪(RegularizedOMP)就是其中一种改进方法。本篇将在上一篇《压缩感知重构算法之正交匹配追踪(OMP)》的基础上给出正则化正交匹配追踪(ROMP)算法的MATLAB函数代码，并且给出单次测试例程代码、测量数M与重构成功概率关系曲线绘制例程代码。

0、符号说明如下：

压缩观测y=Φx，其中y为观测所得向量M×1，x为原信号N×1（M<

(1) y为观测所得向量，大小为M×1

(2)x为原信号，大小为N×1

(3)θ为K稀疏的，是信号在x在某变换域的稀疏表示

(4)Φ称为观测矩阵、测量矩阵、测量基，大小为M×N

(5)Ψ称为变换矩阵、变换基、稀疏矩阵、稀疏基、正交基字典矩阵，大小为N×N

(6)A称为测度矩阵、传感矩阵、CS信息算子，大小为M×N

上式中，一般有K<

1、ROMP重构算法流程：

我将原子选择过程封装成了一个MATLAB函数，代码如下(Regularize.m)：

1.function [val,pos] = Regularize(product,Kin)

2.%Regularize Summary of this function goes here

3.% Detailed explanation goes here

4.% product = A'*r_n;%传感矩阵A各列与残差的内积

5.% K为稀疏度

6.% pos为选出的各列序号

7.% val为选出的各列与残差的内积值

8.% Reference:Needell D，Vershynin R. Uniform uncertainty principle and

9.% signal recovery via regularized orthogonal matching pursuit.

10.% Foundations of Computational Mathematics, 2009,9(3): 317-334.

11. productabs = abs(product);%取绝对值

12. [productdes,indexproductdes] = sort(productabs,'descend');%降序排列

13. for ii = length(productdes):-1:1

14. if productdes(ii)>1e-6%判断productdes中非零值个数

15. break;

16. end

17. end

18. %Identify:Choose a set J of the K biggest coordinates

19. if ii>=Kin

20. J = indexproductdes(1:Kin);%集合J

21. Jval = productdes(1:Kin);%集合J对应的序列值

22. K = Kin;

23. else%or all of its nonzero coordinates,whichever is smaller

24. J = indexproductdes(1:ii);%集合J

25. Jval = productdes(1:ii);%集合J对应的序列值

26. K = ii;

27. end

28. %Regularize:Among all subsets J0∈J with comparable coordinates

29. MaxE = -1;%循环过程中存储最大能量值

30. for kk = 1:K

31. J0_tmp = zeros(1,K);iJ0 = 1;

32. J0_tmp(iJ0) = J(kk);%以J(kk)为本次寻找J0的基准(最大值)

33. Energy = Jval(kk)^2;%本次寻找J0的能量

34. for mm = kk+1:K

35. if Jval(kk)<2*Jval(mm)%找到符合|u(i)|<=2|u(j)|的

36. iJ0 = iJ0 + 1;%J0自变量增1

37. J0_tmp(iJ0) = J(mm);%更新J0

38. Energy = Energy + Jval(mm)^2;%更新能量

39. else%不符合|u(i)|<=2|u(j)|的

40. break;%跳出本轮寻找，因为后面更小的值也不会符合要求

41. end

42. end

43. if Energy>MaxE%本次所得J0的能量大于前一组

44. J0 = J0_tmp(1:iJ0);%更新J0

45. MaxE = Energy;%更新MaxE，为下次循环做准备

46. end

47. end

48. pos = J0;

49. val = productabs(J0);

50.end

2、正则化正交匹配追踪(ROMP)MATLAB代码(CS_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.% 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:Needell D，Vershynin R．Signal recovery from incomplete and

11.% inaccurate measurements via regularized orthogonal matching pursuit[J]．

12.% IEEE Journal on Selected Topics in Signal Processing，2010，4(2)：310—316.

13. [y_rows,y_columns] = size(y);

14. if y_rows

15. y = y';%y should be a column vector

16. end

17. [M,N] = size(A);%传感矩阵A为M*N矩阵

18. theta = zeros(N,1);%用来存储恢复的theta(列向量)

19. At = zeros(M,3*K);%用来迭代过程中存储A被选择的列

20. Pos_theta = zeros(1,2*K);%用来迭代过程中存储A被选择的列序号

21. 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次

25. product = A'*r_n;%传感矩阵A各列与残差的内积

26. %[val,pos] = max(abs(product));%找到最大内积绝对值，即与残差最相关的列

27. [val,pos] = Regularize(product,K);%按正则化规则选择原子

28. At(:,Index+1:Index+length(pos)) = A(:,pos);%存储这几列

29. Pos_theta(Index+1:Index+length(pos)) = pos;%存储这几列的序号

30. if Index+length(pos)<=M%At的行数大于列数，此为最小二乘的基础(列线性无关)

31. Index = Index+length(pos);%更新Index，为下次循环做准备

32. else%At的列数大于行数，列必为线性相关的,At(:,1:Index)'*At(:,1:Index)将不可逆

33. break;%跳出for循环

34. end

35. A(:,pos) = zeros(M,length(pos));%清零A的这几列(其实此行可以不要因为它们与残差正交)

36. %y=At(:,1:Index)*theta，以下求theta的最小二乘解(Least Square)

37. theta_ls = (At(:,1:Index)'*At(:,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. if norm(r_n)<1e-6%Repeat the steps until r=0

41. break;%跳出for循环

42. end

43. if Index>=2*K%or until |I|>=2K