淺談壓縮感知（二十二）：壓縮感知重構(gòu)算法之正則化正交匹配追蹤（ROMP）
[導(dǎo)讀]淺談壓縮感知（二十二）：壓縮感知重構(gòu)算法之正則化正交匹配追蹤（ROMP）主要內(nèi)容：ROMP的算法流程ROMP的MATLAB實(shí)現(xiàn)一維信號(hào)的實(shí)驗(yàn)與結(jié)果測(cè)量數(shù)M與重構(gòu)成功概率關(guān)系的實(shí)驗(yàn)與結(jié)果一、ROMP的算
淺談壓縮感知（二十二）：壓縮感知重構(gòu)算法之正則化正交匹配追蹤（ROMP）
主要內(nèi)容：
ROMP的算法流程ROMP的MATLAB實(shí)現(xiàn)一維信號(hào)的實(shí)驗(yàn)與結(jié)果測(cè)量數(shù)M與重構(gòu)成功概率關(guān)系的實(shí)驗(yàn)與結(jié)果一、ROMP的算法流程
正則化正交匹配追蹤ROMP算法流程與OMP的最大不同之處就在于從傳感矩陣A中選擇列向量的標(biāo)準(zhǔn)，OMP每次只選擇與殘差內(nèi)積絕對(duì)值最大的那一列，而ROMP則是先選出內(nèi)積絕對(duì)值最大的K列（若所有內(nèi)積中不夠K個(gè)非零值則將內(nèi)積值非零的列全部選出），然后再?gòu)倪@K列中按正則化標(biāo)準(zhǔn)再選擇一遍，即為本次迭代選出的列向量（一般并非只有一列）。正則化標(biāo)準(zhǔn)意思是選擇各列向量與殘差內(nèi)積絕對(duì)值的最大值不能比最小值大兩倍以上(comparable coordinates)且能量最大的一組(with the maximal energy)，因?yàn)闈M足條件的子集并非只有一組。
二、ROMP的MATLAB實(shí)現(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為選出的各列序號(hào) %???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);?%取絕對(duì)值 ????[productdes,indexproductdes]?=?sort(productabs,'descend');?%降序排列 ????for?ii?=?length(productdes):-1:1 ????????if?productdes(ii)>1e-6?%判斷productdes中非零值個(gè)數(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對(duì)應(yīng)的序列值 ????????K?=?Kin; ????else?%?or?all?of?its?nonzero?coordinates,whichever?is?smaller ????????J?=?indexproductdes(1:ii);??%集合J ????????Jval?=?productdes(1:ii);??%集合J對(duì)應(yīng)的序列值 ????????K?=?ii; ????end ????%?Regularize:Among?all?subsets?J0∈J?with?comparable?coordinates ????MaxE?=?-1;??%循環(huán)過(guò)程中存儲(chǔ)最大能量值 ????for?kk?=?1:K ????????J0_tmp?=?zeros(1,K);iJ0?=?1; ????????J0_tmp(iJ0)?=?J(kk);??%以J(kk)為本次尋找J0的基準(zhǔn)(最大值) ????????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)做準(zhǔ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);?%用來(lái)存儲(chǔ)恢復(fù)的theta(列向量) ????At?=?zeros(M,3*K);?%用來(lái)迭代過(guò)程中存儲(chǔ)A被選擇的列 ????pos_num?=?zeros(1,2*K);?%用來(lái)迭代過(guò)程中存儲(chǔ)A被選擇的列序號(hào) ????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)積絕對(duì)值，即與殘差最相關(guān)的列 ????????[val,pos]?=?Regularize(product,K);?%按正則化規(guī)則選擇原子 ????????At(:,Index+1:Index+length(pos))?=?A(:,pos);?%存儲(chǔ)這幾列 ????????pos_num(Index+1:Index+length(pos))?=?pos;?%存儲(chǔ)這幾列的序號(hào) ????????if?Index+length(pos)<=M?%At的行數(shù)大于列數(shù)，此為最小二乘的基礎(chǔ)(列線性無(wú)關(guān)) ????????????Index?=?Index+length(pos);?%更新Index，為下次循環(huán)做準(zhǔn)備 ????????else?%At的列數(shù)大于行數(shù)，列必為線性相關(guān)的,At(:,1:Index)'*At(:,1:Index)將不可逆 ????????????break;?%跳出for循環(huán) ????????end ????????A(:,pos)?=?zeros(M,length(pos));?%清零A的這幾列(其實(shí)此行可以不要,因?yàn)樗鼈兣c殘差正交) ????????%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;%恢復(fù)出的theta end
三、一維信號(hào)的實(shí)驗(yàn)與結(jié)果
%壓縮感知重構(gòu)算法測(cè)試
clear?all;close?all;clc;
M?=?128;%觀測(cè)值個(gè)數(shù)
N?=?256;%信號(hào)x的長(zhǎng)度
K?=?12;%信號(hào)x的稀疏度
Index_K?=?randperm(N);
x?=?zeros(N,1);
x(Index_K(1:K))?=?5*randn(K,1);%x為K稀疏的，且位置是隨機(jī)的
Psi?=?eye(N);%x本身是稀疏的，定義稀疏矩陣為單位陣x=Psi*theta
Phi?=?randn(M,N);%測(cè)量矩陣為高斯矩陣
A?=?Phi?*?Psi;%傳感矩陣
y?=?Phi?*?x;%得到觀測(cè)向量y
%%?恢復(fù)重構(gòu)信號(hào)x
tic
theta?=?CS_ROMP(y,A,K);
x_r?=?Psi?*?theta;%?x=Psi?*?theta
toc
%%?繪圖
figure;
plot(x_r,'k.-');%繪出x的恢復(fù)信號(hào)
hold?on;
plot(x,'r');%繪出原信號(hào)x
hold?off;
legend('Recovery','Original')
fprintf('n恢復(fù)殘差：');
norm(x_r-x)%恢復(fù)殘差四、測(cè)量數(shù)M與重構(gòu)成功概率關(guān)系的實(shí)驗(yàn)與結(jié)果
clear?all;close?all;clc;
%%?參數(shù)配置初始化
CNT?=?1000;?%對(duì)于每組(K,M,N)，重復(fù)迭代次數(shù)
N?=?256;?%信號(hào)x的長(zhǎng)度
Psi?=?eye(N);?%x本身是稀疏的，定義稀疏矩陣為單位陣x=Psi*theta
K_set?=?[4,12,20,28,36];?%信號(hào)x的稀疏度集合
Percentage?=?zeros(length(K_set),N);?%存儲(chǔ)恢復(fù)成功概率
%%?主循環(huán)，遍歷每組(K,M,N)
tic
for?kk?=?1:length(K_set)
????K?=?K_set(kk);%本次稀疏度
????M_set?=?K:5:N;%M沒(méi)必要全部遍歷，每隔5測(cè)試一個(gè)就可以了
????PercentageK?=?zeros(1,length(M_set));%存儲(chǔ)此稀疏度K下不同M的恢復(fù)成功概率
????kk
????for?mm?=?1:length(M_set)
???????M?=?M_set(mm)%本次觀測(cè)值個(gè)數(shù)
???????P?=?0;
???????for?cnt?=?1:CNT?%每個(gè)觀測(cè)值個(gè)數(shù)均運(yùn)行CNT次
????????????Index_K?=?randperm(N);
????????????x?=?zeros(N,1);
????????????x(Index_K(1:K))?=?5*randn(K,1);%x為K稀疏的，且位置是隨機(jī)的????????????????
????????????Phi?=?randn(M,N);%測(cè)量矩陣為高斯矩陣
????????????A?=?Phi?*?Psi;%傳感矩陣
????????????y?=?Phi?*?x;%得到觀測(cè)向量y
????????????theta?=?CS_ROMP(y,A,K);%恢復(fù)重構(gòu)信號(hào)theta
????????????x_r?=?Psi?*?theta;%?x=Psi?*?theta
????????????if?norm(x_r-x)<1e-6%如果殘差小于1e-6則認(rèn)為恢復(fù)成功
????????????????P?=?P?+?1;
????????????end
???????end
???????PercentageK(mm)?=?P/CNT*100;%計(jì)算恢復(fù)概率
????end
????Percentage(kk,1:length(M_set))?=?PercentageK;
end
toc
save?ROMPMtoPercentage?%運(yùn)行一次不容易，把變量全部存儲(chǔ)下來(lái)
%%?繪圖
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的恢復(fù)信號(hào)
????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)');
下載文檔
