淺談壓縮感知（二十三）：壓縮感知重構(gòu)算法之壓縮采樣匹配追蹤（CoSaMP）

淺談壓縮感知（二十三）：壓縮感知重構(gòu)算法之壓縮采樣匹配追蹤（CoSaMP）

主要內(nèi)容：

CoSaMP的算法流程CoSaMP的MATLAB實現(xiàn)一維信號的實驗與結(jié)果測量數(shù)M與重構(gòu)成功概率關(guān)系的實驗與結(jié)果一、CoSaMP的算法流程

壓縮采樣匹配追蹤(CompressiveSampling MP)是D. Needell繼ROMP之后提出的又一個具有較大影響力的重構(gòu)算法。CoSaMP也是對OMP的一種改進，每次迭代選擇多個原子，除了原子的選擇標(biāo)準之外，它有一點不同于ROMP：ROMP每次迭代已經(jīng)選擇的原子會一直保留，而CoSaMP每次迭代選擇的原子在下次迭代中可能會被拋棄。

二、CS_CoSaMP的MATLAB實現(xiàn)(CS_CoSaMP.m)

function?[?theta?]?=?CS_CoSaMP(?y,A,K?)
%???CS_CoSaOMP
%???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
%???現(xiàn)在已知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.
????[m,n]?=?size(y);
????if?m<n
????????y?=?y';?%y?should?be?a?column?vector
????end
????[M,N]?=?size(A);?%傳感矩陣A為M*N矩陣
????theta?=?zeros(N,1);?%用來存儲恢復(fù)的theta(列向量)
????pos_num?=?[];?%用來迭代過程中存儲A被選擇的列序號
????res?=?y;?%初始化殘差(residual)為y
????for?kk=1:K?%最多迭代K次
????????%(1)?Identification
????????product?=?A'*res;?%傳感矩陣A各列與殘差的內(nèi)積
????????[val,pos]=sort(abs(product),'descend');
????????Js?=?pos(1:2*K);?%選出內(nèi)積值最大的2K列
????????%(2)?Support?Merger
????????Is?=?union(pos_num,Js);?%Pos_theta與Js并集
????????%(3)?Estimation
????????%At的行數(shù)要大于列數(shù)，此為最小二乘的基礎(chǔ)(列線性無關(guān))
????????if?length(Is)<=M
????????????At?=?A(:,Is);?%將A的這幾列組成矩陣At
????????else?%At的列數(shù)大于行數(shù)，列必為線性相關(guān)的,At'*At將不可逆
????????????if?kk?==?1
????????????????theta_ls?=?0;
????????????end
????????????break;?%跳出for循環(huán)
????????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_num?=?Is(pos(1:K));
????????theta_ls?=?theta_ls(pos(1:K));
????????%At(:,pos(1:K))*theta_ls是y在At(:,pos(1:K))列空間上的正交投影
????????res?=?y?-?At(:,pos(1:K))*theta_ls;?%更新殘差?
????????if?norm(res)<1e-6?%Repeat?the?steps?until?r=0
????????????break;?%跳出for循環(huán)
????????end
????end
????theta(pos_num)=theta_ls;?%恢復(fù)出的theta
end

三、一維信號的實驗與結(jié)果

%壓縮感知重構(gòu)算法測試
clear?all;close?all;clc;
M?=?64;?%觀測值個數(shù)
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

%%?恢復(fù)重構(gòu)信號x
tic
theta?=?CS_CoSaMP(?y,A,K?);
x_r?=?Psi?*?theta;?%?x=Psi?*?theta
toc

%%?繪圖
figure;
plot(x_r,'k.-');?%繪出x的恢復(fù)信號
hold?on;
plot(x,'r');?%繪出原信號x
hold?off;
legend('Recovery','Original')
fprintf('n恢復(fù)殘差：');
norm(x_r-x)?%恢復(fù)殘差

四、測量數(shù)M與重構(gòu)成功概率關(guān)系的實驗與結(jié)果

clear?all;close?all;clc;

%%?參數(shù)配置初始化
CNT?=?1000;?%對于每組(K,M,N)，重復(fù)迭代次數(shù)
N?=?256;?%信號x的長度
Psi?=?eye(N);?%x本身是稀疏的，定義稀疏矩陣為單位陣x=Psi*theta
K_set?=?[4,12,20,28,36];?%信號x的稀疏度集合
Percentage?=?zeros(length(K_set),N);?%存儲恢復(fù)成功概率

%%?主循環(huá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的恢復(fù)成功概率
????for?mm?=?1:length(M_set)
???????M?=?M_set(mm);?%本次觀測值個數(shù)
???????fprintf('K=%d,M=%dn',K,M);
???????P?=?0;
???????for?cnt?=?1:CNT?%每個觀測值個數(shù)均運行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);?%恢復(fù)重構(gòu)信號theta
????????????x_r?=?Psi?*?theta;?%?x=Psi?*?theta
????????????if?norm(x_r-x)<1e-6?%如果殘差小于1e-6則認為恢復(fù)成功
????????????????P?=?P?+?1;
????????????end
???????end
???????PercentageK(mm)?=?P/CNT*100;?%計算恢復(fù)概率
????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的恢復(fù)信號
????hold?on;
end
hold?off;
xlim([0?256]);
legend('K=4','K=12','K=20','K=28','K=36');
xlabel('Number?of?measurements(M)');
ylabel('Percentage?recovered');
title('Percentage?of?input?signals?recovered?correctly(N=256)(Gaussian)');