淺談壓縮感知（二十二）：壓縮感知重構算法之正則化正交匹配追蹤（ROMP）

主要內(nèi)容：

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

正則化正交匹配追蹤ROMP算法流程與OMP的最大不同之處就在于從傳感矩陣A中選擇列向量的標準，OMP每次只選擇與殘差內(nèi)積絕對值最大的那一列，而ROMP則是先選出內(nèi)積絕對值最大的K列（若所有內(nèi)積中不夠K個非零值則將內(nèi)積值非零的列全部選出），然后再從這K列中按正則化標準再選擇一遍，即為本次迭代選出的列向量（一般并非只有一列）。正則化標準意思是選擇各列向量與殘差內(nèi)積絕對值的最大值不能比最小值大兩倍以上(comparable coordinates)且能量最大的一組(with the maximal energy)，因為滿足條件的子集并非只有一組。

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

1、正則化代碼Regularize.m

function?[val,pos]?=?Regularize(product,Kin)
%???Regularize
%???Detailed?explanation?goes?here
%???product?=?A'*r_n;%傳感矩陣A各列與殘差的內(nèi)積
%???K為稀疏度
%???pos為選出的各列序號
%???val為選出的各列與殘差的內(nèi)積值
%???Reference:Needell?D，Vershynin?R.?Uniform?uncertainty?principle?and
%???signal?recovery?via?regularized?orthogonal?matching?pursuit.?
%???Foundations?of?Computational?Mathematics,?2009,9(3):?317-334.??
????productabs?=?abs(product);?%取絕對值
????[productdes,indexproductdes]?=?sort(productabs,'descend');?%降序排列
????for?ii?=?length(productdes):-1:1
????????if?productdes(ii)>1e-6?%判斷productdes中非零值個數(shù)
????????????break;
????????end
????end
????%?Identify:Choose?a?set?J?of?the?K?biggest?coordinates
????if?ii>=Kin
????????J?=?indexproductdes(1:Kin);?%集合J
????????Jval?=?productdes(1:Kin);?%集合J對應的序列值
????????K?=?Kin;
????else?%?or?all?of?its?nonzero?coordinates,whichever?is?smaller
????????J?=?indexproductdes(1:ii);??%集合J
????????Jval?=?productdes(1:ii);??%集合J對應的序列值
????????K?=?ii;
????end
????%?Regularize:Among?all?subsets?J0∈J?with?comparable?coordinates
????MaxE?=?-1;??%循環(huán)過程中存儲最大能量值
????for?kk?=?1:K
????????J0_tmp?=?zeros(1,K);iJ0?=?1;
????????J0_tmp(iJ0)?=?J(kk);??%以J(kk)為本次尋找J0的基準(最大值)
????????Energy?=?Jval(kk)^2;??%本次尋找J0的能量
????????for?mm?=?kk+1:K
????????????if?Jval(kk)<2*Jval(mm)?%找到符合|u(i)|<=2|u(j)|的
????????????????iJ0?=?iJ0?+?1;?%J0自變量增1
????????????????J0_tmp(iJ0)?=?J(mm);?%更新J0
????????????????Energy?=?Energy?+?Jval(mm)^2;?%更新能量
????????????else?%不符合|u(i)|MaxE?%本次所得J0的能量大于前一組
????????????J0?=?J0_tmp(1:iJ0);?%更新J0
????????????MaxE?=?Energy;?%更新MaxE，為下次循環(huán)做準備
????????end
????end
????pos?=?J0;
????val?=?productabs(J0);
end

2、ROMP代碼CS_ROMP.m

function?[?theta?]?=?CS_ROMP(?y,A,K?)
%???CS_ROMP
%???Detailed?explanation?goes?here
%???y?=?Phi?*?x
%???x?=?Psi?*?theta
%????y?=?Phi*Psi?*?theta
%???令?A?=?Phi*Psi,?則y=A*theta
%???現(xiàn)在已知y和A，求theta
%???Reference:Needell?D，Vershynin?R．Signal?recovery?from?incomplete?and
%???inaccurate?measurements?via?regularized?orthogonal?matching?pursuit[J]．
%???IEEE?Journal?on?Selected?Topics?in?Signal?Processing，2010，4(2)：310—316.
????[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);?%用來存儲恢復的theta(列向量)
????At?=?zeros(M,3*K);?%用來迭代過程中存儲A被選擇的列
????pos_num?=?zeros(1,2*K);?%用來迭代過程中存儲A被選擇的列序號
????Index?=?0;
????res?=?y;?%初始化殘差(residual)為y
????%Repeat?the?following?steps?K?times(or?until?|I|>=2K)
????for?ii=1:K?%迭代K次
????????product?=?A'*res;?%傳感矩陣A各列與殘差的內(nèi)積
????????%[val,pos]?=?max(abs(product));?%找到最大內(nèi)積絕對值，即與殘差最相關的列
????????[val,pos]?=?Regularize(product,K);?%按正則化規(guī)則選擇原子
????????At(:,Index+1:Index+length(pos))?=?A(:,pos);?%存儲這幾列
????????pos_num(Index+1:Index+length(pos))?=?pos;?%存儲這幾列的序號
????????if?Index+length(pos)<=M?%At的行數(shù)大于列數(shù)，此為最小二乘的基礎(列線性無關)
????????????Index?=?Index+length(pos);?%更新Index，為下次循環(huán)做準備
????????else?%At的列數(shù)大于行數(shù)，列必為線性相關的,At(:,1:Index)'*At(:,1:Index)將不可逆
????????????break;?%跳出for循環(huán)
????????end
????????A(:,pos)?=?zeros(M,length(pos));?%清零A的這幾列(其實此行可以不要,因為它們與殘差正交)
????????%y=At(:,1:Index)*theta，以下求theta的最小二乘解(Least?Square)
????????theta_ls?=?(At(:,1:Index)'*At(:,1:Index))^(-1)*At(:,1:Index)'*y;?%最小二乘解
????????%At(:,1:Index)*theta_ls是y在At(:,1:Index)列空間上的正交投影
????????res?=?y?-?At(:,1:Index)*theta_ls;?%更新殘差
????????if?norm(res)=2*K?%or?until?|I|>=2K
????????????break;?%跳出for循環(huán)
????????end
????end
????theta(pos_num(1:Index))=theta_ls;%恢復出的theta
end

三、一維信號的實驗與結果

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

%%?恢復重構信號x
tic
theta?=?CS_ROMP(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)%恢復殘差

四、測量數(shù)M與重構成功概率關系的實驗與結果

clear?all;close?all;clc;
%%?參數(shù)配置初始化
CNT?=?1000;?%對于每組(K,M,N)，重復迭代次數(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);?%存儲恢復成功概率

%%?主循環(huán)，遍歷每組(K,M,N)
tic
for?kk?=?1:length(K_set)
????K?=?K_set(kk);%本次稀疏度
????M_set?=?K:5:N;%M沒必要全部遍歷，每隔5測試一個就可以了
????PercentageK?=?zeros(1,length(M_set));%存儲此稀疏度K下不同M的恢復成功概率
????kk
????for?mm?=?1:length(M_set)
???????M?=?M_set(mm)%本次觀測值個數(shù)
???????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);%測量矩陣為高斯矩陣
????????????A?=?Phi?*?Psi;%傳感矩陣
????????????y?=?Phi?*?x;%得到觀測向量y
????????????theta?=?CS_ROMP(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?ROMPMtoPercentage?%運行一次不容易，把變量全部存儲下來

%%?繪圖
S?=?['-ks';'-ko';'-kd';'-kv';'-k*'];
figure;
for?kk?=?1:length(K_set)
????K?=?K_set(kk);
????M_set?=?K:5:N;
????L_Mset?=?length(M_set);
????plot(M_set,Percentage(kk,1:L_Mset),S(kk,:));%繪出x的恢復信號
????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)');