分段正交匹配追踪算法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 中的一个函数而已，链接为/，最新版本为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