在看代碼之前，先拜讀了ROMP的經典文章：Needell D，VershyninR．Signal recovery from incompleteand inaccurate measurements via regularized orthogonal matching pursuit[J]. IEEEJournal on Selected Topics in Signal Processing, 2010, 4(2): 310—316.
看完一臉懵逼，真的沒看懂啥，雖然頁數不多，在下文中就單純的借鑒文章中的算法流程。
正交匹配追蹤算法每次迭代均只選擇與殘差最相關的一列，自然人們會想：“每次迭代是否可以多選幾列呢?”，正則化正交匹配追蹤(RegularizedOMP)就是其中一種改進方法。本篇將在上一篇《壓縮感知重構算法之正交匹配追蹤(OMP)》的基礎上給出正則化正交匹配追蹤(ROMP)算法的MATLAB函數代碼，並且給出單次測試例程代碼、測量數M與重構成功概率關系曲線繪制例程代碼。
0、符號說明如下
壓縮觀測y=Φx，其中y為觀測所得向量M×1，x為原信號N×1（M<<N）。x一般不是稀疏的，但在某個變換域Ψ是稀疏的，即x=Ψθ，其中θ為K稀疏的，即θ只有K個非零項。此時y=ΦΨθ，令A=ΦΨ，則y=Aθ。
(1) y為觀測所得向量，大小為M×1
(2) x為原信號，大小為N×1
(3) θ為K稀疏的，是信號在x在某變換域的稀疏表示
(4) Φ稱為觀測矩陣、測量矩陣、測量基，大小為M×N
(5) Ψ稱為變換矩陣、變換基、稀疏矩陣、稀疏基、正交基字典矩陣，大小為N×N
(6) A稱為測度矩陣、傳感矩陣、CS信息算子，大小為M×N
上式中，一般有K<<M<<N，后面三個矩陣各個文獻的叫法不一，以后我將Φ稱為測量矩陣、將Ψ稱為稀疏矩陣、將A稱為傳感矩陣。
1、ROMP重構算法流程
正則化正交匹配追蹤算法流程與OMP的最大不同之處就在於從傳感矩陣A中選擇列向量的標准，OMP每次只選擇與殘差內積絕對值最大的那一列，而ROMP則是先選出內積絕對值最大的K列（若所有內積中不夠K個非零值則將內積值非零的列全部選出），然后再從這K列中按正則化標准再選擇一遍，即為本次迭代選出的列向量（一般並非只有一列）。正則化標准意思是選擇各列向量與殘差內積絕對值的最大值不能比最小值大兩倍以上(comparable coordinates)且能量最大的一組(with the maximal energy)，因為滿足條件的子集並非只有一組。似乎用敘述語言描述不清楚，下面給出一種實現第(2)(3)步的算法流程圖：
貼出文獻[1]中的算法流程：
看完論文后對算法的理解並不是很深入，下面結合博客中的算法流程來對ROMP算法流程進行解釋。上述流程圖講的是正則化的過程，最多經過K次迭代可選出全部所需的原子。在Identify中首先將所得到的內積值按降序排列，然計算內積中非零元素的個數，然后選取前K個內積值或者所有非零值（也就是論文中提到的選擇集合比較小的那個），記錄選取的內積值所對應的列序號，構成集合J，相應的內積值構成集合Jval。
正則化的過程結合代碼來解釋。
function[val,pos]=Regularize(product,Kin)
% Regularize Summary of this function goes here
% 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
第12行代碼中用到了函數sort進行排序，采用的是降序方式，indexproductdes索引中保存的是排序后的內積值productdes在原來集合productabs中原來所在的位置。
第13-17行判斷大於0的內積值的個數，並在第19到27行中進行選擇，將內積值所對應的列序號形成集合J，並將所選擇的內積值組成集合Jval。
第29行，首先初始化 MaxE為-1.
第30行，接下來是在第某次選擇出的J中選擇子集J0 ，總共迭代K次，K為原始信號非零元素的個數。
接着聊聊如何選擇J0 ，首先選擇Jval（kk）（為與K區分，選用與代碼中一樣的kk形式）為基准，初始化m=kk，然后遍歷m+1即（k+1，也就是此次k的下一個內積值）到K，判斷Jval(kk)<=2*Jval(mm)，也就是論文中所說的|u(i)|<=2|u(j)|，這個地方想了很久，然后去運行代碼調試才搞明白。這里先假設我已經找到了能量最大的一組。然后我選擇出來的J0 所包含的列向量的序號有此次的k，還有滿足Jval(kk)<=2*Jval(mm)的mm，在代碼中開始已經將J（kk）的值賦給了J0_tmp(iJ0)（初始iJ0=1），也就是代碼的第32行，后續滿足條件的J（mm）也分別賦值給了J0_tmp（iJ0）（iJ0=iJ0+1），所以最后的J0 =J0_tmp(1:iJ0)（也就是初始的基准Jval（kk）和后面滿足條件的m），在流程圖中J0=J（1：m），我覺得這是不夠嚴謹的，一開始想了很久想不明白，難道變換的只是m而已嗎，所以我覺得在這里是不太正確的。
另外仔細觀察論文中的這一項
，也是困擾了我很久的地方，注意這里的u（i）相當於此次我們的Jval（k），u（j）是我們要遍歷判斷的，最后i，j都是屬於J0的，所以也說明了上一段J0所應包含的列向量序號的集合。
，也是困擾了我很久的地方，注意這里的u（i）相當於此次我們的Jval（k），u（j）是我們要遍歷判斷的，最后i，j都是屬於J0的，所以也說明了上一段J0所應包含的列向量序號的集合。
接着說明J0的選擇，應該是在所有滿足條件的J的子集中能量最大的一組，第43到46行進行了能量的比較，如果能量比上一次的能量大才會進行J0的賦值，否則進入下一次循環直至結束。到這里一次的J0選擇完畢。
2、正則化正交匹配追蹤(ROMP)MATLAB代碼(CS_ROMP.m)
這個函數是完全基於上一篇中的CS_OMP.m修改而成的。
function [ theta ] = CS_ROMP( y,A,K )
%CS_ROMP Summary of this function goes here
%Version: 1.0 written by jbb0523 @2015-04-24
% 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.
[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(列向量)
At = zeros(M,3*K);%用來迭代過程中存儲A被選擇的列
Pos_theta = zeros(1,2*K);%用來迭代過程中存儲A被選擇的列序號
Index = 0;
r_n = y;%初始化殘差(residual)為y
%Repeat the following steps K times(or until |I|>=2K)
for ii=1:K%迭代K次
product = A'*r_n;%傳感矩陣A各列與殘差的內積
%[val,pos] = max(abs(product));%找到最大內積絕對值，即與殘差最相關的列
[val,pos] = Regularize(product,K);%按正則化規則選擇原子
At(:,Index+1:Index+length(pos)) = A(:,pos);%存儲這幾列
Pos_theta(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)列空間上的正交投影
r_n = y - At(:,1:Index)*theta_ls;%更新殘差
if norm(r_n)<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_theta(1:Index))=theta_ls;%恢復出的theta
end
相比CS_OMP.m，ROMP所做的修改已用顏色標出。
首先解釋下第19行和20行，博客中的解釋是：
然后我還是沒有太明白，但是傳感矩陣滿足2K階RIP，滿足2K階RIP的矩陣任意2K列線性無關。可能跟這個有關系，以后再看看。
接着是第21行，為什么索引值Index不直接設置為1呢，每次選擇的原子有可能為幾列，則這次所選擇出來的原子存放的位置，應該從上次存放的最后一列的位置+1到這次所選擇的原子長度加上上次存放的最后一列的位置。比如第一次存放的位置應該為0+1：0+length（pos），下一次存放的位置表示為Index+1：Index+length(pos)。
繼續解釋第30到33行，這里是判斷我們所選擇出的原子構成的矩陣At行數與列數比較的關系。At選擇的列向量都是非零的，也就是說At是列滿秩的矩陣。（列滿秩就是列秩等於列數,就是初等變換以后沒有一列全為0. 滿秩矩陣是一個很重要的概念, 它是判斷一個矩陣是否可逆的充分必要條件）看了下線性代數，還沒有看懂。。。
第40行到第44行是對循環結束條件的判斷，或者殘差小於一定范圍，或者是索引集合Index>=2K。
3、ROMP單次重構測試代碼
以下測試代碼與上一篇中的OMP單次測試代碼基本完全一致，除了M和K參數設置及調用CS_ROMP函數三處之外。
%壓縮感知重構算法測試
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)%恢復殘差
運行結果如下：（信號為隨機生成，所以每次結果均不一樣）
1）圖：
2）Commandwindows
Elapsedtime is 0.022132 seconds.
恢復殘差：
ans=
7.8066e-015
4、測量數M與重構成功概率關系曲線繪制例程代碼
以下測試代碼與上一篇中的OMP測量數M與重構成功概率關系曲線繪制例程代碼基本完全一致。本程序在循環中填加了“kk”一行代碼並將“M = M_set(mm)”一行的分號去掉，這是為了在運行過程中可以觀察程序運行狀態、知道程序到哪一個位置。
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
forkk=1:length(K_set)
K=K_set(kk);%本次稀疏度
M_set=K:5:N;%M沒必要全部遍歷，每隔5測試一個就可以了
PercentageK=zeros(1,length(M_set));%存儲此稀疏度K下不同M的恢復成功概率
kk
formm=1:length(M_set)
M=M_set(mm)%本次觀測值個數
P=0;
forcnt=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
ifnorm(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 ROMPMtoPercentage1000%運行一次不容易，把變量全部存儲下來
%% 繪圖
S=['-ks';'-ko';'-kd';'-kv';'-k*'];
figure;
forkk=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([0256]);
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)');
本程序運行結果：
參考文獻：
[1] Needell D，VershyninR．Signal recovery from incompleteand inaccurate measurements via regularized orthogonal matching pursuit[J]. IEEEJournal on Selected Topics in Signal Processing, 2010, 4(2): 310—316.