主要內容：
- ROMP的算法流程
- ROMP的MATLAB實現
- 一維信號的實驗與結果
- 測量數M與重構成功概率關系的實驗與結果
一、ROMP的算法流程
正則化正交匹配追蹤ROMP算法流程與OMP的最大不同之處就在於從傳感矩陣A中選擇列向量的標准，OMP每次只選擇與殘差內積絕對值最大的那一列，而ROMP則是先選出內積絕對值最大的K列（若所有內積中不夠K個非零值則將內積值非零的列全部選出），然后再從這K列中按正則化標准再選擇一遍，即為本次迭代選出的列向量（一般並非只有一列）。正則化標准意思是選擇各列向量與殘差內積絕對值的最大值不能比最小值大兩倍以上(comparable coordinates)且能量最大的一組(with the maximal energy)，因為滿足條件的子集並非只有一組。
二、ROMP的MATLAB實現(CS_ROMP.m)
1、正則化代碼Regularize.m
function [val,pos] = Regularize(product,Kin) % Regularize % Detailed explanation goes here % product = A'*r_n;%傳感矩陣A各列與殘差的內積 % K為稀疏度 % pos為選出的各列序號 % val為選出的各列與殘差的內積值 % 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中非零值個數 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; %循環過程中存儲最大能量值 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)|<=2|u(j)|的 break; %跳出本輪尋找，因為后面更小的值也不會符合要求 end end if Energy>MaxE %本次所得J0的能量大於前一組 J0 = J0_tmp(1:iJ0); %更新J0 MaxE = Energy; %更新MaxE，為下次循環做准備 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 % 現在已知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各列與殘差的內積 %[val,pos] = max(abs(product)); %找到最大內積絕對值，即與殘差最相關的列 [val,pos] = Regularize(product,K); %按正則化規則選擇原子 At(:,Index+1:Index+length(pos)) = A(:,pos); %存儲這幾列 pos_num(Index+1:Index+length(pos)) = pos; %存儲這幾列的序號 if Index+length(pos)<=M %At的行數大於列數，此為最小二乘的基礎(列線性無關) Index = Index+length(pos); %更新Index，為下次循環做准備 else %At的列數大於行數，列必為線性相關的,At(:,1:Index)'*At(:,1:Index)將不可逆 break; %跳出for循環 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)<1e-6 %Repeat the steps until r=0 break; %跳出for循環 end if Index>=2*K %or until |I|>=2K break; %跳出for循環 end end theta(pos_num(1:Index))=theta_ls;%恢復出的theta end
三、一維信號的實驗與結果
%壓縮感知重構算法測試 clear all;close all;clc; M = 128;%觀測值個數 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)%恢復殘差
四、測量數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 = K:5:N;%M沒必要全部遍歷，每隔5測試一個就可以了 PercentageK = zeros(1,length(M_set));%存儲此稀疏度K下不同M的恢復成功概率 kk for mm = 1:length(M_set) M = M_set(mm)%本次觀測值個數 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);%測量矩陣為高斯矩陣 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)');
