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

压缩感知重构算法之子空间追踪(SP)

如果掌握了压缩采样匹配追踪(CoSaMP)后，再去学习子空间追踪(Subspace Pursuit)是一件非常简单的事情，因为它们几乎是完全一样的。

SP的提出时间比CoSaMP提出时间略晚，首个论文版本是参考文献[1]，后来更新了两次，最后在IEEE Transactions on Information Theory发表[2]。从算法角度来讲，SP与CoSaMP差别非常小，这一点作者也意识到了，在文献[1]首页的左下角就有注释：

在文献[2]第2页提到了SP与CoSaMP的具体不同：

从上面可以知道，SP与CoSaMP主要区别在于“Ineach iteration, in the SP algorithm, only K new candidates are added, while theCoSAMP algorithm adds 2K vectors.”，即SP每次选择K个原子，而CoSaMP则选择2K个原子；这样带来的好处是“This makes the SP algorithm computationally moreefficient,”。

以下是文献[2]中的给出的SP算法流程：

这个算法流程的初始化(Initialization)其实就是类似于CoSaMP的第1次迭代，注意第(1)步中选择了K个原子：“K indices correspo nding to the largest magnitude entries”，在CoSaMP里这里要选择2K个最大的原子，后面的其它流程都一样。这里第(5)步增加了一个停止迭代的条件：当残差经过迭代后却变大了的时候就停止迭代。

不只是SP作者认识到了自己的算法与CoSaMP的高度相似性，CoSaMP的作者也同样关注到了SP算法，在文献[3]中就提到：

文献[3]是CoSaMP原始提出文献的第2个版本，文献[3]的早期版本[4]是没有提及SP算法的。

鉴于SP与CoSaMP如此相似，这里不就再单独给出SP的步骤了，参考《压缩感知重构算法之压缩采样匹配追踪(CoSaMP)》，只需将第(2)步中的2K改为K即可。

引用文献[5]的3.5节中的几句话：“贪婪类算法虽然复杂度低运行速度快，但其重构精度却不如BP类算法，为了寻求复杂度和精度更好地折中，SP算法应运而生”，“SP算法与CoSaMP算法一样其基本思想也是借用回溯的思想，在每步迭代过程中重新估计所有候选者的可信赖性”，“SP算法与CoSaMP算法有着类似的性质与优缺点”。

子空间追踪代码可实现如下(CS_SP.m)，通过对比可以知道该程序与CoSaMP的代码基本完全一致。本代码未考虑文献[2]中的给出的SP算法流程的第(5)步。代码可参见参考文献[6]中的Demo_CS_SP.m。

[plain]view plaincopy

1.function [ theta ] = CS_SP( y,A,K )

2.%CS_SP Summary of this function goes here

3.%Version: 1.0 written by jbb0523 @2015-05-01

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.% K is the sparsity level

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

11.% Reference:Dai W，Milenkovic O．Subspace pursuit for compressive sensing

12.% signal reconstruction[J]．IEEE Transactions on Information Theory，

13.% 2009，55(5)：2230-2249.

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

15. if y_rows

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

17. end

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

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

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

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

22. for kk=1:K%最多迭代K次

23. %(1) Identification

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

25. [val,pos]=sort(abs(product),'descend');

26. Js = pos(1:K);%选出内积值最大的K列

27. %(2) Support Merger

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

29. %(3) Estimation

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

31. if length(Is)<=M

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

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

34. break;%跳出for循环

35. end

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

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

38. %(4) Pruning

39. [val,pos]=sort(abs(theta_ls),'descend');

40. %(5) Sample Update