壓縮采樣匹配追蹤(CompressiveSampling MP)是D. Needell繼ROMP之后提出的又一個具有較大影響力的重構算法。CoSaMP也是對OMP的一種改進，每次迭代選擇多個原子，除了原子的選擇標准之外，它有一點不同於ROMP：ROMP每次迭代已經選擇的原子會一直保留，而CoSaMP每次迭代選擇的原子在下次迭代中可能會被拋棄。
在這之前先讀了下參考論文[1]，論文前面還是看得懂一點的，講了一些壓縮感知的基礎知識，還聊到了壓縮重構方法主要分為三類，但是到了第2部分介紹算法的時候又看不懂了，感覺符號都還沒聊清楚就開始講流程了。佩服看得懂的博主，還說很容易就看懂了。。。希望自己好好努力也能看懂這些外文文獻，fighting啦！
那么我先把論文中的流程貼出來，並沒有對符號作過多的解釋。。。 supp(x)表示x的支撐集
然后在網上找到了符合論文中符號的代碼。
function Sest = cosaomp(Phi,u,K,tol,maxiterations) Sest = zeros(size(Phi,2),1); v = u; t = 1; numericalprecision = 1e-12; T = []; while (t <= maxiterations) && (norm(v)/norm(u) > tol) y = abs(Phi'*v); [vals,z] = sort(y,'descend'); Omega = find(y >= vals(2*K) & y > numericalprecision); T = union(Omega,T); b = pinv(Phi(:,T))*u; [vals,z] = sort(abs(b),'descend'); Kgoodindices = (abs(b) >= vals(K) & abs(b) > numericalprecision); T = T(Kgoodindices); Sest = zeros(size(Phi,2),1); phit = Phi(:,T); b = pinv(phit)*u; Sest(T) = b; v = u - phit*b; t = t+1; end
接下來綜合代碼我准備強行解釋一波論文算法的偽代碼流程，哎呀半懂半懂希望以后要全懂全懂。
1.Identification（識別）
大意是說要構造一個signal proxy，在偽代碼中構造signal proxy是y=Phi*v，下圖是從論文中摘出來的，突然明白了這段話的意思，首先翻譯一下。信號重構的最大難點在於找到目標信號中這些最大項所在的位置。CoSaOMP受到RIP的啟發，假設字典矩陣的RIP常數為遠遠小於1的一個值，對s稀疏的信號x，y=Phi*Phi x可以作為信號的一個代理。因為y的每一個s向量的結合的能量與信號x中s個向量的能量相對應。（我覺得這里的Phi應該是理解為字典矩陣的，因為計算內積的時候我們是選擇將字典矩陣與殘差相乘，殘差初始化為觀測向量也就是Phi*x）。這一步對應着代碼的第8行。
接着是偽代碼中所說的Identify large components，也就是找到內積值中最大的2K項，復制給Ω，對應上述代碼的第10行。
2.Support Merger（合並支撐集）
代碼第11行。
3.Estimation
這里是求解一個最小二乘問題，pinv是求偽逆矩陣。
“b|Tc←0”中的“Tc”應該是T的補集(complementary set)，向量b的元素序號為全集，子集T對應的元素等於最小二乘解，補集對應的元素為零。
4.Pruning（修剪）
代碼第13行，選出b中K個最大項。
5.Sample Update（更新）
強行解釋結束了，接下來貼出博主的解釋。
1、CoSaMP重構算法流程
步驟（5）稍微有點繞，綜合代碼理解一下還是不難的。
2、壓縮采樣匹配追蹤(CoSaOMP)Matlab代碼(CS_CoSaMP.m)
function [ theta ] = CS_CoSaMP( y,A,K )
%CS_CoSaOMP Summary of this function goes here
%Created by jbb0523@@2015-04-29
%Version: 1.1 modified by jbb0523 @2015-05-09
% Detailed explanation goes here
% y = Phi * x
% x = Psi * theta
% y = Phi*Psi * theta
% 令 A = Phi*Psi, 則y=A*theta
% K is the sparsity level
% 現在已知y和A，求theta
% Reference:Needell D，Tropp J A．CoSaMP：Iterative signal recovery from
% incomplete and inaccurate samples[J]．Applied and Computation Harmonic
% Analysis，2009，26：301-321.
[y_rows,y_columns] = size(y);
if y_rows<y_columns
y = y';%y should be a column vector
end
[M,N] = size(A);%傳感矩陣A為M*N矩陣
theta = zeros(N,1);%用來存儲恢復的theta(列向量)
Pos_theta = [];%用來迭代過程中存儲A被選擇的列序號
r_n = y;%初始化殘差(residual)為y
for kk=1:K%最多迭代K次
%(1) Identification
product = A'*r_n;%傳感矩陣A各列與殘差的內積
[val,pos]=sort(abs(product),'descend');
Js = pos(1:2*K);%選出內積值最大的2K列
%(2) Support Merger
Is = union(Pos_theta,Js);%Pos_theta與Js並集
%(3) Estimation
%At的行數要大於列數，此為最小二乘的基礎(列線性無關)
if length(Is)<=M
At = A(:,Is);%將A的這幾列組成矩陣At
else%At的列數大於行數，列必為線性相關的,At'*At將不可逆
if kk == 1
theta_ls = 0;
end
break;%跳出for循環
end
%y=At*theta，以下求theta的最小二乘解(Least Square)
theta_ls = (At'*At)^(-1)*At'*y;%最小二乘解
%(4) Pruning
[val,pos]=sort(abs(theta_ls),'descend');
%(5) Sample Update
Pos_theta = Is(pos(1:K));
theta_ls = theta_ls(pos(1:K));
%At(:,pos(1:K))*theta_ls是y在At(:,pos(1:K))列空間上的正交投影
r_n = y - At(:,pos(1:K))*theta_ls;%更新殘差
if norm(r_n)<1e-6%Repeat the steps until r=0
break;%跳出for循環
end
end
theta(Pos_theta)=theta_ls;%恢復出的theta
end
3、CoSaMP單次重構測試代碼
以下測試代碼基本與OMP單次重構測試代碼一樣。
%壓縮感知重構算法測試
clear all;close all;clc;
M = 64;%觀測值個數
N = 256;%信號x的長度
K = 12;%信號x的稀疏度
Index_K = randperm(N);
x = zeros(N,1);
x(Index_K(1:K)) = 5*randn(K,1);%x為K稀疏的，且位置是隨機的
Psi = eye(N);%x本身是稀疏的，定義稀疏矩陣為單位陣x=Psi*theta
Phi = randn(M,N);%測量矩陣為高斯矩陣
A = Phi * Psi;%傳感矩陣
y = Phi * x;%得到觀測向量y
%% 恢復重構信號x
tic
theta = CS_CoSaMP( y,A,K );
x_r = Psi * theta;% x=Psi * theta
toc
%% 繪圖
figure;
plot(x_r,'k.-');%繪出x的恢復信號
hold on;
plot(x,'r');%繪出原信號x
hold off;
legend('Recovery','Original')
fprintf('\n恢復殘差：');
norm(x_r-x)%恢復殘差
運行結果如下：（信號為隨機生成，所以每次結果均不一樣）
1）圖：
2）Command windows
Elapsedtime is 0.073375 seconds.
恢復殘差：
ans=
7.3248e-015
4、測量數M與重構成功概率關系曲線繪制例程代碼
以下測試代碼基本與OMP測量數M與重構成功概率關系曲線繪制代碼一樣。增加了“fprintf('K=%d,M=%d\n',K,M);”，可以觀察程序運行進度。
clear all;close all;clc;
%% 參數配置初始化
CNT = 1000;%對於每組(K,M,N)，重復迭代次數
N = 256;%信號x的長度
Psi = eye(N);%x本身是稀疏的，定義稀疏矩陣為單位陣x=Psi*theta
K_set = [4,12,20,28,36];%信號x的稀疏度集合
Percentage = zeros(length(K_set),N);%存儲恢復成功概率
%% 主循環，遍歷每組(K,M,N)
tic
for kk = 1:length(K_set)
K = K_set(kk);%本次稀疏度
M_set = 2*K:5:N;%M沒必要全部遍歷，每隔5測試一個就可以了
PercentageK = zeros(1,length(M_set));%存儲此稀疏度K下不同M的恢復成功概率
for mm = 1:length(M_set)
M = M_set(mm);%本次觀測值個數
fprintf('K=%d,M=%d\n',K,M);
P = 0;
for cnt = 1:CNT %每個觀測值個數均運行CNT次
Index_K = randperm(N);
x = zeros(N,1);
x(Index_K(1:K)) = 5*randn(K,1);%x為K稀疏的，且位置是隨機的
Phi = randn(M,N)/sqrt(M);%測量矩陣為高斯矩陣
A = Phi * Psi;%傳感矩陣
y = Phi * x;%得到觀測向量y
theta = CS_CoSaMP(y,A,K);%恢復重構信號theta
x_r = Psi * theta;% x=Psi * theta
if norm(x_r-x)<1e-6%如果殘差小於1e-6則認為恢復成功
P = P + 1;
end
end
PercentageK(mm) = P/CNT*100;%計算恢復概率
end
Percentage(kk,1:length(M_set)) = PercentageK;
end
toc
save CoSaMPMtoPercentage1000 %運行一次不容易，把變量全部存儲下來
%% 繪圖
S = ['-ks';'-ko';'-kd';'-kv';'-k*'];
figure;
for kk = 1:length(K_set)
K = K_set(kk);
M_set = 2*K:5:N;
L_Mset = length(M_set);
plot(M_set,Percentage(kk,1:L_Mset),S(kk,:));%繪出x的恢復信號
hold on;
end
本程序運行結果：
最后摘出文獻[4]中關於ROMP優缺點的分析：
ROMP 算法雖具有貪婪算法的速度以及凸優化算法的強有力的理論保證，但其與StOMP 算法一樣，對稀疏度K 的依賴性太大，稀疏度估計的准確與否，將會影響到算法的收斂性、收斂速度、魯棒性以及重構信號的性能。針對此缺點，提出了正則化自適應匹配追蹤( Regularized adaptive matching pursuit，RAMP) 算法，該算法將ROMP 算法進行改進，引入回溯的思想，自適應調節候選集原子的數目。
CoSaMP 算法為了提高算法的收斂速度和算法效率，通過回溯的思想從原子庫里選擇多個相關原子同時剔除部分不相關原子。設算法的迭代步長為K，候選集中最多有3K 個原子，每次最多剔除K個原子，以保證支撐集中有2K 個原子。CoSaMP 算法雖結合了組合算法的思想保證了算法的收斂速度，而且提供了嚴格的誤差界，但其與ROMP 算法一樣需要已知信號的稀疏度K。然而實際應用中的信號往往是近似稀疏的，其稀疏度K需要去估計。如果低估了K 的值，算法精確重構的能力會下降甚至可能會消除，導致算法不再收斂;若高估了K 的值，算法的魯棒性和重構精度都將會下降，使得重構的信號的誤差增大，導致最終重構得到的信號失真。此外，CoSaMP 算法每步迭代時增加與剔除原子所依據的原則不一樣，增加原子時依據的是原子與余量的相關性原則，剔除原子時依據的是重構信號的大小，標准的不同可能會導致對支撐集的估計不准確。針對第一個問題，仍可以采取改變信號的稀疏度，多次運行，使得重構精度最高的稀疏度即作為信號的稀疏度的思想，也可以引入回溯思想自適應選擇原子。針對第二個問題，可以改變剔除原子的標准，采用依據原子與余量的相關性原則，剔除相關性小的原子，使得每步迭代時增加與剔除原子所依據的標准一樣。
參考文獻：
[1]D. Needell, J.A. Tropp．CoSaMP: Iterative signal recoveryfrom incomplete and inaccurate samples.http://arxiv.org/pdf/0803.2392v2.pdf
[2]D.Needell, J.A. Tropp．CoSaMP: Iterative signal recoveryfrom incomplete and inaccurate samples[J]. Communications of theACM，2010，53(12)：93-100.
(http://dl.acm.org/citation.cfm?id=1859229)
[4] 楊真真，楊震，孫林慧.信號壓縮重構的正交匹配追蹤類算法綜述[J]. 信號處理，2013，29(4)：486-496.