贪婪算法中正交匹配追踪算法gOMP的原理及仿真文档格式.docx

贪婪算法中正交匹配追踪算法gOMP的原理及仿真文档格式.docx

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

3.%Version:

1.0 

written 

by 

jbb0523 

@2015-05-08 

4.% 

Detailed 

explanation 

5.% 

y 

Phi 

* 

x 

6.% 

Psi 

7.% 

Phi*Psi 

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. 

y_rows<

y_columns 

19. 

y'

;

%y 

should 

be 

a 

column 

vector 

20. 

21. 

[M,N] 

size（A）;

%传感矩阵A为M*N矩阵 

22. 

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. 

length（Sk）==length（Pos_theta） 

30. 

ii 

== 

1 

31. 

theta_ls 

0;

32. 

33. 

break;

34. 

35. 

length（Sk）>

M 

36. 

37. 

38. 

39. 

40. 

41. 

At 

A（:

Sk）;

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

42. 

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

Square） 

43. 

（At'

*At）^（-1）*At'

*y;

%最小二乘解 

44. 

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

45. 

- 

At*theta_ls;

%更新残差 

46. 

Sk;

47. 

norm（r_n）<

1e-6 

48. 

%quit 

the 

iteration 

49. 

50. 

51. 

theta（Pos_theta）=theta_ls;

%恢复出的theta 

52.end 

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

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

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

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

2.clear 

all;

close 

clc;

3.M 

128;

%观测值个数 

4.N 

256;

%信号x的长度 

5.K 

30;

%信号x的稀疏度 

6.Index_K 

randperm（N）;

7.x 

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 

Psi;

%传感矩阵 

12.y 

x;

%得到观测向量y 

13.%% 

恢复重构信号x 

14.tic 

15.theta 

y,A,K）;

16.x_r 

theta;

% 

x=Psi 

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相应的测试代码基本一致，可以合并在一起运行，只须在主循环内多加几种算法重构就行。

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

2.% 

绘制参考文献中的Fig.1 

3.% 

Elapsed 

time 

is 

798.718246 

seconds.（@20150509pm） 

8.clear 

9.%% 

参数配置初始化 

10.CNT 

1000;

%对于每组（K,M,N），重复迭代次数 

11.N 

12.Psi 

13.M_set 

[128];

%测量值集合 

14.KIND 

['

OMP 

ROMP 

StOMP 

SP 

CoSaMP 

... 

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） 

M_set（mm）;

%本次测量值个数 

K_set 

5:

70;

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

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

PercentageM 

zeros（size（KIND,1）,length（K_set））;

kk 

length（K_set） 

K 

K_set（kk）;

%本次信号x的稀疏度K 

P 

zeros（1,size（KIND,1））;

fprintf（'

M=%d,K=%d\n'

M,K）;

cnt 

CNT 

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

Index_K 

x（Index_K（1:

%

（1）OMP 

CS_OMP（y,A,K）;

%恢复重构信号theta 

x_r 

norm（x_r-x）<

1e-6%如果残差小于1e-6则认为恢复成功 

P

（1） 

+ 

1;

%

（2）ROMP 

CS_ROMP（y,A,K）;

P

（2） 

%（3）StOMP 

CS_StOMP（y,A）;

P（3） 

52. 

53. 

%（4）SP 

54. 

CS_SP（y,A,K）;

55. 

56. 

57. 

P（4） 

58. 

59. 

%（5）CoSaMP 

60. 

CS_CoSaMP（y,A,K）;

61. 

62. 

63. 

P（5） 

64. 

65. 

%（6）gOMP,S=3 

66. 

CS_gOMP（y,A,K,3）;

67. 

68. 

69. 

P（6） 

70. 

71. 

%（7）gOMP,S=6 

72. 

CS_gOMP（y,A,K,6）;

73. 

74. 

75. 

P（7） 

76. 

77. 

%（8）gOMP,S=9 

78. 

CS_gOMP（y,A,K,9）;

79. 

80. 

81. 

P（8） 

82. 

83. 

84. 

iii 

size（KIND,1） 

85. 

PercentageM（iii,kk） 

P（iii）/CNT*100;

%计算恢复概率 

86. 

87. 

88. 

jjj 

89. 

Percentage（1:

length（K_set）,mm,jjj） 

PercentageM（jjj,:

90. 

91.end 

92.toc 

93.save 

KtoPercentage1000gOMP 

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

94.%% 

95.S 

-ks'

-ko'

-yd'

-gv'

-b*'

-r.'

-rx'

-r+'

96.figure;

97.for 

98. 

99. 

100. 

L_Kset 

length（K_set）;

101. 

102. 

plot（K_set,Percentage（1:

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

））;

103. 

hold 

104. 

105.end 

106.hold 

107.xlim（[5 

70]）;

108.legend（'

OMP'

ROMP'

StOMP'

SP'

CoSaMP'

... 

109. 

110.xlabel（'

Sparsity 

level 

K'

111.ylabel（'

The 

Probability 

Exact 

Reconstruction'

112.title（'

Prob. 

exact 

recovery 

vs. 

signal 

sparsity 

K（M=128,N=256）（Gaussian）'

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

本程序运行结果：

5、结语

我很好奇：

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

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

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

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

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

//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.

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