贪婪算法中正交匹配追踪算法gOMP的基础学习知识原理及仿真.docx

资源描述

贪婪算法中正交匹配追踪算法gOMP的基础学习知识原理及仿真.docx

《贪婪算法中正交匹配追踪算法gOMP的基础学习知识原理及仿真.docx》由会员分享，可在线阅读，更多相关《贪婪算法中正交匹配追踪算法gOMP的基础学习知识原理及仿真.docx（15页珍藏版）》请在冰豆网上搜索。

贪婪算法中正交匹配追踪算法gOMP的基础学习知识原理及仿真.docx

贪婪算法中正交匹配追踪算法gOMP的基础学习知识原理及仿真

压缩感知重构算法之广义正交匹配追踪（gOMP）

广义正交匹配追踪（GeneralizedOMP,gOMP）算法可以看作为OMP算法的一种推广，由文献[1]提出，第1作者本硕为哈工大毕业，发表此论文时在KoreaUniversity攻读博士学位。

OMP每次只选择与残差相关最大的一个，而gOMP则是简单地选择最大的S个。

之所以这里表述为“简单地选择”是相比于ROMP之类算法的，不进行任何其它处理，只是选择最大的S个而已。

0、符号说明如下：

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

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<

注意：

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

1、gOMP重构算法流程：

2、广义正交匹配追踪（gOMP）MATLAB代码（CS_gOMP.m）

本代码完全是为了保证和前面的各算法代法格式一致，可以直接使用该实验室网站提供的代码[2]压缩包中的islsp_EstgOMP.m。

[plain] viewplaincopy

1.function [ theta ] = CS_gOMP（ y,A,K,S ）

2.%CS_gOMP Summary of this function goes here

3.%Version:

1.0 written by jbb0523 @2015-05-08

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:

Jian Wang, Seokbeop Kwon, Byonghyo Shim. Generalized

11.% orthogonal matching pursuit, IEEE Transactions on Signal Processing,

12.% vol. 60, no. 12, pp. 6202-6216, Dec. 2012.

13.% Available at:

http:

//islab.snu.ac.kr/paper/tsp_gOMP.pdf

14. if nargin < 4

15. S = round（max（K/4, 1））;

16. end

17. [y_rows,y_columns] = size（y）;

18. if y_rows

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

20. end

21. [M,N] = size（A）;%传感矩阵A为M*N矩阵

22. theta = zeros（N,1）;%用来存储恢复的theta（列向量）

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

24. r_n = y;%初始化残差（residual）为y

25. for ii=1:

K%迭代K次，K为稀疏度

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

27. [val,pos]=sort（abs（product）,'descend'）;%降序排列

28. Sk = union（Pos_theta,pos（1:

S））;%选出最大的S个

29. if length（Sk）==length（Pos_theta）

30. if ii == 1

31. theta_ls = 0;

32. end

33. break;

34. end

35. if length（Sk）>M

36. if ii == 1

37. theta_ls = 0;

38. end

39. break;

40. end

41. At = A（:

Sk）;%将A的这几列组成矩阵At

42. %y=At*theta，以下求theta的最小二乘解（Least Square）

43. theta_ls = （At'*At）^（-1）*At'*y;%最小二乘解

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

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

46. Pos_theta = Sk;

47. if norm（r_n）<1e-6

48. break;%quit the iteration

49. end

50. end

51. theta（Pos_theta）=theta_ls;%恢复出的theta

52.end

3、gOMP单次重构测试代码（CS_Reconstuction_Test.m）

以下测试代码基本与OMP单次重构测试代码一样。

也可参考该实验室网站提供的代码[2]压缩包中的Test_gOMP.m。

[plain] viewplaincopy

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

2.clear all;close all;clc;

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

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

5.K = 30;%信号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_gOMP（ y,A,K）;

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

Elapsedtimeis0.155937seconds.

恢复残差：

ans=

2.3426e-014

4、信号稀疏度K与重构成功概率关系曲线绘制例程代码

以下测试代码为了与文献[1]的Fig.1作比较。

由于暂未研究学习LP算法，所以相比于文献[1]的Fig.1）缺少LP算法曲线，加入了SP算法。

以下测试代码与SAMP相应的测试代码基本一致，可以合并在一起运行，只须在主循环内多加几种算法重构就行。

[plain] viewplaincopy

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

2.% 绘制参考文献中的Fig.1

3.% Reference:

Jian Wang, Seokbeop Kwon, Byonghyo Shim. Generalized

4.% orthogonal matching pursuit, IEEE Transactions on Signal Processing,

5.% vol. 60, no. 12, pp. 6202-6216, Dec. 2012.

6.% Available at:

http:

//islab.snu.ac.kr/paper/tsp_gOMP.pdf

7.% Elapsed time is 798.718246 seconds.（@20150509pm）

8.clear all;close all;clc;

9.%% 参数配置初始化

10.CNT = 1000;%对于每组（K,M,N），重复迭代次数

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

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

13.M_set = [128];%测量值集合

14.KIND = ['OMP ';'ROMP ';'StOMP ';'SP ';'CoSaMP ';...

15. 'gOMP（s=3）';'gOMP（s=6）';'gOMP（s=9）'];

16.Percentage = zeros（N,length（M_set）,size（KIND,1））;%存储恢复成功概率

17.%% 主循环，遍历每组（K,M,N）

18.tic

19.for mm = 1:

length（M_set）

20. M = M_set（mm）;%本次测量值个数

21. K_set = 5:

70;%信号x的稀疏度K没必要全部遍历，每隔5测试一个就可以了

22. %存储此测量值M下不同K的恢复成功概率

23. PercentageM = zeros（size（KIND,1）,length（K_set））;

24. for kk = 1:

length（K_set）

25. K = K_set（kk）;%本次信号x的稀疏度K

26. P = zeros（1,size（KIND,1））;

27. fprintf（'M=%d,K=%d\n',M,K）;

28. for cnt = 1:

CNT %每个观测值个数均运行CNT次

29. Index_K = randperm（N）;

30. x = zeros（N,1）;

31. x（Index_K（1:

K）） = 5*randn（K,1）;%x为K稀疏的，且位置是随机的

32. Phi = randn（M,N）/sqrt（M）;%测量矩阵为高斯矩阵

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

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

35. %

（1）OMP

36. theta = CS_OMP（y,A,K）;%恢复重构信号theta

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

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

39. P

（1） = P

（1） + 1;

40. end

41. %

（2）ROMP

42. theta = CS_ROMP（y,A,K）;%恢复重构信号theta

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

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

45. P

（2） = P

（2） + 1;

46. end

47. %（3）StOMP

48. theta = CS_StOMP（y,A）;%恢复重构信号theta

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

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

51. P（3） = P（3） + 1;

52. end

53. %（4）SP

54. theta = CS_SP（y,A,K）;%恢复重构信号theta

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

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

57. P（4） = P（4） + 1;

58. end

59. %（5）CoSaMP

60. theta = CS_CoSaMP（y,A,K）;%恢复重构信号theta

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

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

63. P（5） = P（5） + 1;

64. end

65. %（6）gOMP,S=3

66. theta = CS_gOMP（y,A,K,3）;%恢复重构信号theta

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

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

69. P（6） = P（6） + 1;

70. end

71. %（7）gOMP,S=6

72. theta = CS_gOMP（y,A,K,6）;%恢复重构信号theta

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

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

75. P（7） = P（7） + 1;

76. end

77. %（8）gOMP,S=9

78. theta = CS_gOMP（y,A,K,9）;%恢复重构信号theta

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

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

81. P（8） = P（8） + 1;

82. end

83. end

84. for iii = 1:

size（KIND,1）

85. PercentageM（iii,kk） = P（iii）/CNT*100;%计算恢复概率

86. end

87. end

88. for jjj = 1:

size（KIND,1）

89. Percentage（1:

length（K_set）,mm,jjj） = PercentageM（jjj,:

）;

90. end

91.end

92.toc

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

94.%% 绘图

95.S = ['-ks';'-ko';'-yd';'-gv';'-b*';'-r.';'-rx';'-r+'];

96.figure;

97.for mm = 1:

length（M_set）

98. M = M_set（mm）;

99. K_set = 5:

70;

100. L_Kset = length（K_set）;

101. for ii = 1:

size（KIND,1）

102. plot（K_set,Percentage（1:

L_Kset,mm,ii）,S（ii,:

））;%绘出x的恢复信号

103. hold on;

104. end

105.end

106.hold off;

107.xlim（[5 70]）;

108.legend（'OMP','ROMP','StOMP','SP','CoSaMP',...

109. 'gOMP（s=3）','gOMP（s=6）','gOMP（s=9）'）;

110.xlabel（'Sparsity level K'）;

111.ylabel（'The Probability of Exact Reconstruction'）;

112.title（'Prob. of exact recovery vs. the signal sparsity K（M=128,N=256）（Gaussian）'）;

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

本程序运行结果：

5、结语

我很好奇：

为什么相比于OMP算法就是简单每次多选几列，重构效果为什么这么好？

居然比复杂的ROMP、CoSaMP、StOMP效果还要好……

该课题组还提出了MMP算法，可参见文献[3]。

更多关于该课题组的信息可去官方网站查询：

http:

//islab.snu.ac.kr/，也可直接查看发表的文章：

http:

//islab.snu.ac.kr/publication.html。

文献[1]最后有两个TABLE，分别是算法的流程和复杂度总结：

谁能告诉我TABLEI表头中的DELETE在这里是什么意思啊？

不是删除的意思么？

6、参考文献

【1】JianWang,SeokbeopKwon,ByonghyoShim. Generalizedorthogonalmatchingpursuit,IEEETransactionsonSignalProcessing,vol.60,no.12,pp.6202-6216,Dec.2012.

Availableat:

http:

//islab.snu.ac.kr/paper/tsp_gOMP.pdf

【2】http:

//islab.snu.ac.kr/paper/gOMP.zip

【3】S.Kwon,J.WangandB.Shim,Multipathmatchingpursuit,IEEETransactionsonInformationTheory,vol.60,no.5,pp.2986-3001,May2014.

Availableat:

http:

//islab.snu.ac.kr/paper/TIT_MMP2014.pdf

展开阅读全文