分段正交匹配追踪算法StOMP

压缩感知重构算法之分段正交匹配追踪(StOMP)

分段正交匹配追踪(StagewiseOMP)或者翻译为逐步正交匹配追踪，它是OMP另一种改进算法，每次迭代可以选择多个原子。此算法的输入参数中没有信号稀疏度K，因此相比于ROMP及CoSaMP有独到的优势。

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<

注意：这里的稀疏表示模型为x=Ψθ，所以传感矩阵A=ΦΨ；而有些文献中稀疏模型为θ=Ψx，而一般Ψ为Hermite 矩阵(实矩阵时称为正交矩阵)，所以Ψ-1=ΨH(实矩阵时为Ψ-1=ΨT)，即x=ΨHθ，所以传感矩阵A=ΦΨH，例如沙威的OMP例程中就是如此。

1、StOMP重构算法流程：

2、分段正交匹配追踪(StOMP)Matlab代码(CS_StOMP.m)

代码参考了文献[4]中的SolveStOMP.m，也可参考文献[5]中的StOMP.m。其实文献[4]是斯坦福的SparseLab 中的一个函数而已，链接为https://www.doczj.com/doc/362218022.html,/，最新版本为2.1，SolveStOMP.m在目录

SparseLab21-Core\SparseLab2.1-Core\Solvers里面。

1.function[theta]=CS_StOMP(y,A,S,ts)

2.%CS_StOMP Summary ofthisfunctiongoes here

3.%Version:1.0 writtenbyjbb0523@2015-04-29

4.%Detailedexplanationgoeshere

5.%y =Phi*x

6.% x=Psi*theta

7.% y=Phi*Psi*theta

8.% 令 A = Phi*Psi, 则y=A*theta

9.% S is the maximum number of StOMP iterations to perform

10.% ts is the threshold parameter

11.% 现在已知y和A，求theta

12.% Reference:Donoho D L，Tsaig Y，Drori I，Starck J L．Sparse solution of

13.% underdetermined linear equations by stagewise orthogonal matching

14.% pursuit[J]．IEEE Transactions on Information Theory，2012，58(2)：1094—1121

15. if nargin < 4

16. ts = 2.5;%ts范围[2,3],默认值为2.5

17. end

18. if nargin < 3

19. S = 10;%S默认值为10

20. end

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

22. if y_rows

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

24. end

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

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

27. Pos_theta = [];%用来迭代过程中存储A被选择的列序号

28. r_n = y;%初始化残差(residual)为y

29. for ss=1:S%最多迭代S次

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

31. sigma = norm(r_n)/sqrt(M);%参见参考文献第3页Remarks(3)

32. Js = find(abs(product)>ts*sigma);%选出大于阈值的列

33. Is = union(Pos_theta,Js);%Pos_theta与Js并集

34. if length(Pos_theta) == length(Is)

35. if ss==1

36. theta_ls = 0;%防止第1次就跳出导致theta_ls无定义

37. end

38. break;%如果没有新的列被选中则跳出循环

39. end

40. %At的行数要大于列数，此为最小二乘的基础(列线性无关)

41. if length(Is)<=M

42. Pos_theta = Is;%更新列序号集合

43. At = A(:,Pos_theta);%将A的这几列组成矩阵At

44. else%At的列数大于行数，列必为线性相关的,At'*At将不可逆

45. if ss==1

46. theta_ls = 0;%防止第1次就跳出导致theta_ls无定义

47. end

48. break;%跳出for循环

49. end

50. %y=At*theta，以下求theta的最小二乘解(Least Square)

51. theta_ls = (At'*At)^(-1)*At'*y;%最小二乘解

52. %At*theta_ls是y在At列空间上的正交投影

53. r_n = y - At*theta_ls;%更新残差

54. if norm(r_n)<1e-6%Repeat the steps until r=0

55. break;%跳出for循环

56. end

57. end

58. theta(Pos_theta)=theta_ls;%恢复出的theta

59.end

3、StOMP单次重构测试代码

以下测试代码基本与OMP单次重构测试代码一样，除了调用CS_StOMP之外，一定要注意这里的测量矩阵Phi =randn(M,N)/sqrt(M)，一定一定！！！

[plain]view plaincopy

1.%压缩感知重构算法测试

2.clear all;close all;clc;

3.M = 64;%观测值个数

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

10.Phi = randn(M,N)/sqrt(M);%测量矩阵为高斯矩阵

11.A = Phi * Psi;%传感矩阵

12.y = Phi * x;%得到观测向量y

13.%% 恢复重构信号x

14.tic

15.theta = CS_StOMP(y,A);

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

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

1）图：

2）Command windows

Elapsedtime is 0.067904 seconds.

恢复残差：

ans=

6.1267e-015

4、门限参数t s、测量数M与重构成功概率关系曲线绘制例程代码

因为文献[1]中对门限参数ts给出的是一个取值范围，所以有必要仿真ts取不同值时的重构效果，因此以下的代码虽然是基于OMP相应的测试代码修改的，但相对来说改动较大。

[plain]view plaincopy

1.clear all;close all;clc;

2.%% 参数配置初始化

https://www.doczj.com/doc/362218022.html,T = 1000;%对于每组(K,M,N)，重复迭代次数

4.N = 256;%信号x的长度

5.Psi = eye(N);%x本身是稀疏的，定义稀疏矩阵为单位阵x=Psi*theta

6.ts_set = 2:0.2:3;

7.K_set = [4,12,20,28,36];%信号x的稀疏度集合

8.Percentage = zeros(N,length(K_set),length(ts_set));%存储恢复成功概率

9.%% 主循环，遍历每组(ts,K,M,N)

10.tic

11.for tt = 1:length(ts_set)

12. ts = ts_set(tt);

13. for kk = 1:length(K_set)

14. K = K_set(kk);%本次稀疏度

15. %M没必要全部遍历，每隔5测试一个就可以了

16. M_set=2*K:5:N;

17. PercentageK = zeros(1,length(M_set));%存储此稀疏度K下不同M的恢复成功概率

18. for mm = 1:length(M_set)

19. M = M_set(mm);%本次观测值个数

20. fprintf('ts=%f,K=%d,M=%d\n',ts,K,M);

21. P = 0;

22. for cnt = 1:CNT %每个观测值个数均运行CNT次

23. Index_K = randperm(N);

24. x = zeros(N,1);

25. x(Index_K(1:K)) = 5*randn(K,1);%x为K稀疏的，且位置是随机的

26. Phi = randn(M,N)/sqrt(M);%测量矩阵为高斯矩阵

27. A = Phi * Psi;%传感矩阵

28. y = Phi * x;%得到观测向量y

29. theta = CS_StOMP(y,A,10,ts);%恢复重构信号theta

30. x_r = Psi * theta;% x=Psi * theta

31. if norm(x_r-x)<1e-6%如果残差小于1e-6则认为恢复成功

32. P = P + 1;

33. end

34. end

35. PercentageK(mm) = P/CNT*100;%计算恢复概率

36. end

37. Percentage(1:length(M_set),kk,tt) = PercentageK;

38. end

39.end

40.toc

41.save StOMPMtoPercentage1000 %运行一次不容易，把变量全部存储下来

42.%% 绘图

43.for tt = 1:length(ts_set)

44. S = ['-ks';'-ko';'-kd';'-kv';'-k*'];

45. figure;

46. for kk = 1:length(K_set)

47. K = K_set(kk);

48. M_set=2*K:5:N;

49. L_Mset = length(M_set);

50. plot(M_set,Percentage(1:L_Mset,kk,tt),S(kk,:));%绘出x的恢复信号

51. hold on;

52. end

53. hold off;

54. xlim([0 256]);

55. legend('K=4','K=12','K=20','K=28','K=36');

56. xlabel('Number of measurements(M)');

57. ylabel('Percentage recovered');

58. title(['Percentage of input signals recovered correctly(N=256,ts=',...

59. num2str(ts_set(tt)),')(Gaussian)']);

60.end

61.for kk = 1:length(K_set)

62. K = K_set(kk);

63. M_set=2*K:5:N;

64. L_Mset = length(M_set);

65. S = ['-ks';'-ko';'-kd';'-kv';'-k*';'-k+'];

66. figure;

67. for tt = 1:length(ts_set)

68. plot(M_set,Percentage(1:L_Mset,kk,tt),S(tt,:));%绘出x的恢复信号

69. hold on;

70. end

71. hold off;

72. xlim([0 256]);

73. legend('ts=2.0','ts=2.2','ts=2.4','ts=2.6','ts=2.8','ts=3.0');

74. xlabel('Number of measurements(M)');

75. ylabel('Percentage recovered');

76. title(['Percentage of input signals recovered correctly(N=256,K=',...

77. num2str(K),')(Gaussian)']);

78.end

本程序在联想ThinkPadE430C笔记本（4GBDDR3内存，i5-3210）上运行共耗时4707.513276秒，程序中将所有数据均通过“save StOMPMtoPercentage1000”存储了下来，以后可以再对数据进行分析，只需“load

S tOMPMtoPercentage1000”即可。

程序运行结束会出现6+5=11幅图，前6幅图分别是ts分别为2.0、2.2、2.4、2.6、2.8和3.0时的测量数M与重构成功概率关系曲线（类似于OMP此部分，这里只是对每一个不同的ts画出一幅图），后5幅图是分别将稀疏度K为4、12、20、28、32时将六种ts取值的测量数M与重构成功概率关系曲线绘制在一起以比较ts对重构结果的影响。

对于前6幅图这里只给出ts=2.4时的曲线图：

对于后5幅图这里全部给出，为了清楚地看出ts的影响，这里把图的横轴拉伸：

通过对比可以看出，总体上讲ts=2.4或ts=2.6时效果较好，较大和较小重构效果都会降低，这里由于没有ts=2.5的情况，但我们推测ts=2.5应该是一个比较好的值，因此一般默认取为2.5即可。

5、结语

有关StOMP的流程图可参见文献[1]的Fig.1：

有关StOMP门限的选取在文献[1]中也有提及：

关于这个门限的来源文献[1]有也有一个推导，注意推导过程中的N(0,1/n)：

作者在文献[1]中提出StOMP，这篇文章的发表时间是2012年，但看一下这篇文章的左下角会发现一个问题：

注意，文章在2006-04-05就投稿了，直到2011-08-17修回并被接受，然后2012年才发表。也就是说审稿就审了五年多，按说文章第一作者是大牛，虽说IEEE Transactions on InformationTheory是一个顶级期刊，但对Donoho D L来说也应该不算是什么难事，不知道为什么会出现这种现象。当然，英文文献里有个有趣的现象是还未发表就开始被引用，所以你经常会发现参考文献里会有“to be published”或“submittedfor publication”，如果到国内就是参考文献里出现“已录用”或“已投稿”，不知道审稿人看到会是什么心情。不过老外的牛文章似乎都是先在会议上发表再投期刊，如文献[1]首面左下角注明了“Thematerial in thi s paper was presented in part at the Allerton Conference onCommuncation, Control, and Computing, Sept. 2007, Monticello, IL, USA.”，而前面讲过的CoSaMP的提出文章就更夸张了，版本四五个。

尽管StOMP输入参数中不需要信号的稀疏度，但门限设置与测量矩阵有密切的关系，文献[1]中的门限也只适用于随机高斯矩阵而己，因此限制了此算法的应用。

参考文献：

[1]Donoho D L，Tsaig Y，DroriI，Starck J L．Sparsesolution of underdetermined linear equations by stagewise orthogonal matchingpursuit[J]．IEEE Transactions on InformationTheory，2012，58(2)：1094—1121.

[2]杨真真，杨震，孙林慧.信号压缩重构的正交匹配追踪类算法综述[J]. 信号处理，2013，29(4)：486-496.

[3]吴赟.压缩感知测量矩阵的研究[D]. 西安电子科技大学硕士学位论文，2012.

[4]https://www.doczj.com/doc/362218022.html,pared. https://www.doczj.com/doc/362218022.html,/downloads196/sourcecode/graph/detail923222.html

[5]付自杰.cs_matlab. https://www.doczj.com/doc/362218022.html,/downloads641/sourcecode/math/detail2595379.html

一种改进的用于稀疏表示的正交匹配追踪算法

第10卷 第5期 信 息 与 电 子 工 程 Vo1．10，No．5 2012年10月 INFORMATION AND ELECTRONIC ENGINEERING Oct．，2012 文章编号：1672-2892(2012)05-0579-05 一种改进的用于稀疏表示的正交匹配追踪算法 王燕霞，张 弓 (南京航空航天大学 电子信息工程学院，江苏 南京 210016) 摘 要：稀疏表示理论在军事目标识别、雷达目标参数估计等领域应用越来越广，而目标信号的稀疏表示通常不唯一，因此产生了大量的稀疏表示算法。本文基于现有稀疏表示算法的研究，提出一种改进的正交匹配追踪(OMP)算法。首先采用非线性下降的阈值更快速地选择原子，确定备选原子集，提高了算法速度；其次用正则化的二次筛选剔除备选原子集中能量较低的原子，保证了算法精确度；并设置迭代停止条件实现算法的稀疏度自适应。实验结果表明，本文算法可以实现稀疏表示求解精确度和速度上的平衡，求解速度比基追踪(BP)算法快，精确度比OMP 、正则化OMP(ROMP)、基于自适应OMP 回溯(BAOMP)算法高。 关键词：过完备字典；稀疏表示；正交匹配追踪；正则化 中图分类号：TN850.6；TP753 文献标识码：A An improved orthogonal matching pursuit algorithm for sparse representation WANG Yan -xia，ZHANG Gong (College of Electronics and Information Engineering，Nanjing University of Aeronautics and Astronautics，Nanjing Jiangsu 210016，China) Abstract：Usually sparse representation of signal is not unique, which results in a large number of sparse representation algorithms. An improved Orthogonal Matching Pursuit(OMP) algorithm is proposed. The atoms are selected more quickly with nonlinear decline threshold and the set of alternative atoms is determined, which improves the algorithm speed. Regularized secondary screening can remove lower -energy atoms from the alternative atoms set to ensure the accuracy. A stop condition for iteration is preset to realize the adaptive sparsity of new algorithm. Simulation results show that, the improved algorithm can keep a balance between accuracy and speed for sparse solving with a faster speed than Basis Pursuit(BP) algorithm and a higher accuracy than OMP, Regularized OMP(ROMP) and Backtracking -based Adaptive OMP(BAOMP) algorithms. Key words：over -complete dictionary；sparse representation；orthogonal matching pursuit；regularized 稀疏表示理论在军事目标识别[1]、雷达目标参数估计[2]等领域的应用越来越广，2010年Vishal 和Nasser 等人基于科曼奇族前视红外数据库将稀疏表示用于军事目标的识别，稀疏求解的结果表明算法取得了较好的识别效果[1]。对于稀疏表示理论应用于各种目标信号处理来说，信号的表示应该能对自身不同性质，以及各性质间的差异提供明确信息，这就促进了信号在基于大量波形的过完备字典上的分解[1]。然而过完备库中的波形数远远超过信号的维数，因此它对给定信号的表示不是唯一的，由此产生了信号处理工作中大量的稀疏表示算法。 基于过完备字典的稀疏表示起源于20世纪90年代，Mallat 和Zhang 首次提出了匹配追踪(Matching Pursuit ，MP)算法来解决该稀疏分解问题[3?4]，通常这些算法可以分为3类：基于p l 范数正则化的方法、迭代收缩算法和贪婪追踪算法。基于p l 范数正则化的方法通过求解在重构条件下系数的p l 范数最小化来寻找稀疏解，此类比较典型的算法有基追踪(BP)[5]；迭代收缩算法首先要获得系数的初始集，然后迭代修正获得稀疏表示，早期的这类方法有噪声整形(Noise Shaping ，NS)[6]等；贪婪追踪算法力求从字典中选取与余量最匹配的原子作为信号的近似表示，并通过减去与该原子相关的分量来更新余量，此类算法研究的成果较多，典型的有MP [5],OMP [7],ROMP [8]及BAOMP [9]等。 收稿日期：2011-12-14；修回日期：2012-02-09

分段正交匹配追踪算法StOMP

压缩感知重构算法之分段正交匹配追踪(StOMP) 分段正交匹配追踪(StagewiseOMP)或者翻译为逐步正交匹配追踪，它是OMP另一种改进算法，每次迭代可以选择多个原子。此算法的输入参数中没有信号稀疏度K，因此相比于ROMP及CoSaMP有独到的优势。 0、符号说明如下： 压缩观测y=Φx，其中y为观测所得向量M×1，x为原信号N×1（M<

2、分段正交匹配追踪(StOMP)Matlab代码(CS_StOMP.m) 代码参考了文献[4]中的SolveStOMP.m，也可参考文献[5]中的StOMP.m。其实文献[4]是斯坦福的SparseLab 中的一个函数而已，链接为https://www.doczj.com/doc/362218022.html,/，最新版本为2.1，SolveStOMP.m在目录 SparseLab21-Core\SparseLab2.1-Core\Solvers里面。 1.function[theta]=CS_StOMP(y,A,S,ts) 2.%CS_StOMP Summary ofthisfunctiongoes here 3.%Version:1.0 writtenbyjbb0523@2015-04-29 4.%Detailedexplanationgoeshere 5.%y =Phi*x

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

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

我将原子选择过程封装成了一个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);

贪婪算法中正交匹配追踪算法OMP的原理及仿真

压缩感知重构算法之正交匹配追踪(OMP) 前面经过几篇的基础铺垫，本篇给出正交匹配追踪(OMP)算法的MATLAB函数代码，并且给出单次测试例程代码、测量数M与重构成功概率关系曲线绘制例程代码、信号稀疏度K与重构成功概率关系曲线绘制例程代码。 0、符号说明如下： 压缩观测y=Φx，其中y为观测所得向量M×1，x为原信号N×1（M<

2、正交匹配追踪(OMP)MATLAB代码(CS_OMP.m) [plain]view plaincopy 1.function [ theta ] = CS_OMP( y,A,t ) 2.%CS_OMP Summary of this function goes here 3.%Version: 1.0 written by jbb0523 @2015-04-18 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. [y_rows,y_columns] = size(y); 11. if y_rows

mp&omp 匹配追踪 正交匹配追踪

1. 信号的稀疏表示(sparse representation of signals) 给定一个过完备字典矩阵，其中它的每列表示一种原型信号的原子。给定一个信号y，它可以被表示成这些原子的稀疏线性组合。信号y 可以被表达为y = Dx ，或者 。字典矩阵中所谓过完备性，指的是原子的个数远远大于信号y的长度(其长度很显然是n)，即n<