貪婪算法中ROMP算法的原理介紹及MATLAB仿真

1、壓縮感知重構(gòu)算法之正則化正交匹配追蹤(ROMP)正交匹配追蹤算法每次迭代均只選擇與殘差最相關(guān)的一列，自然人們會想：“每次迭代是否可以多選幾列呢”，正則化正交匹配追蹤(RegularizedOMP)就是其中一種改進(jìn)方法。本篇將在上一篇壓縮感知重構(gòu)算法之正交 匹配追蹤(OMP)的基礎(chǔ)上給出正則化正交匹配追蹤(ROMP)算法的MATLAB函數(shù)代碼，并且給出單次測試?yán)?代碼、測量數(shù)M與重構(gòu)成功概率關(guān)系曲線繪制例程代碼。0、符號說明如下：壓縮觀測丫=血，其中j為觀測所得向量Mx1，*為原信號Nx1(M1e-6%判斷 productdes 中非零值個數(shù)break;endend%Identify:Choo

2、sea set J of the K biggest coordinatesif ii=KinJ = indexproductdes(1:Kin);%集合 JJval = productdes(1:Kin);%集合 J 對應(yīng)的序列值K=Kin;else%or all of its nonzero coordinates,whichever is smaller TOC o 1-5 h z J=indexproductdes(1:ii);%集合JJval = productdes(1:ii);%集合J 對應(yīng)的序列值K=ii;end%Regularize:AmongallsubsetsJ0 e J

3、 with comparable coordinatesMaxE =-1;%循環(huán)過程中存儲最大能量值for kk= 1:KJ0_tmp = zeros(1,K);iJ0 = 1;J0_tmp(iJ0) = J(kk);%以J(kk)為本次尋找J0的基準(zhǔn)(最大值)Energy = Jval(kk)A2;%本次尋找 J0 的能量for mm = kk+1:Kif Jval(kk)2*Jval(mm)%找到符合|u(i)|=2|u(j)| 的iJ0 = iJ0 + 1;%J0 自變量增 1J0_tmp(iJ0) = J(mm);%更新 J0Energy = Energy + Jval(mm)A2;%

4、 更新能量else%不符合|u(i)|MaxE%本次所得J0的能量大于前一組J0 =J0_tmp(1:iJ0);%更新J0MaxE= Energy;%更新MaxE,為下次循環(huán)做準(zhǔn)備endendpos=J0;val=productabs(J0);2、正則化正交匹配追蹤(ROMP)MATLAB代碼(CS_ROMP.m)這個函數(shù)是完全基于上一篇中的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-

5、04-24% 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 pursuitJ.% IEEE Journal on Selected To

6、pics in Signal Processing, 2010, 4(2)： 310316.y_rows,y_columns = size(y);if y_rows=2K)for ii=1:K%迭代 K 次product = A*r_n;%傳感矩陣A各列與殘差的內(nèi)積%val,pos = max(abs(product);%找到最大內(nèi)積絕對值，即與殘差最相關(guān)的列val,pos = Regularize(product,K);% 按正則化規(guī)則選擇原子At(:,Index+1:Index+length(pos) = A(:,pos);% 存儲這幾列Pos_theta(Index+1:Index+le

7、ngth(pos) = pos;% 存儲這幾列的序號if Index+length(pos)=M%At的行數(shù)大于列數(shù)，此為最小二乘的基礎(chǔ)(列線性無關(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)endA(:,pos) = zeros(M,length(pos);%清零A的這幾列(其實此行可以不要因為它們與殘差正交)%y=At(:,1:Index)*theta，以下求 theta 的最小二乘解(Least Squ

8、are)theta_ls = (At(:,1:Index)*At(:,1:Index)A(-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)=2*K%or until |I|=2K44.break;%44.break;%跳出for循環(huán)end45.endendtheta(Pos_theta(1:Index)=theta_ls;%恢復(fù)出的 thetaend3、ROMP單次重構(gòu)測試代碼以下測試

9、代碼與上一篇中的OMP單次測試代碼基本完全一致，除了 M和K參數(shù)設(shè)置及調(diào)用CS_ROMP函數(shù)三處 之外。%壓縮感知重構(gòu)算法測試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*thetaPhi = randn(M,N);%測量矩陣為高斯矩陣A = Phi * Psi;

10、%傳感矩陣y = Phi * x;%得到觀測向量y%恢復(fù)重構(gòu)信號xtictheta = CS_ROMP(y,A,K);x_r = Psi * theta;% x=Psi * thetatoc%繪圖figure;plot(x_r,k.-);%繪出x的恢復(fù)信號hold on;plot(x,r);%繪出原信號 xhold off;legend(Recovery,Original)fprintf(n 恢復(fù)殘差：);norm(x_r-x)%恢復(fù)殘差運行結(jié)果如下：(信號為隨機生成，所以每次結(jié)果均不一樣)1)圖:1510-550-101510-550-102) CommandwindowsElapsed t

11、ime (運行時間)is 0.022132 seconds.恢復(fù)殘差：ans=7.8066e-0154、測量數(shù)M與重構(gòu)成功概率關(guān)系曲線繪制例程代碼以下測試代碼與上一篇中的OMP測量數(shù)M與重構(gòu)成功概率關(guān)系曲線繪制例程代碼基本完全一致。本程序在 循環(huán)中填加了 “kk ”一行代碼并將“M = M_set(mm) ”一行的分號去掉，這是為了在運行過程中可以觀察程序運行狀 態(tài)、知道程序到哪一個位置。重構(gòu)調(diào)用CS_ROMP函數(shù)并將save命令的文件名改為 ROMPMtoPercentage1000，以防止將OMP存儲的數(shù)據(jù)覆蓋(因為本人的所有代碼都在一個目錄下)。clear all;close all;c

12、lc;%參數(shù)配置初始化CNT = 1000;%對于每組(K,M,N)，重復(fù)迭代次數(shù)N = 256;%信號x的長度Psi = eye(N);%x本身是稀疏的，定義稀疏矩陣為單位陣x=Psi*thetaK_set = 4,12,20,28,36;%信號 x 的稀疏度集合Percentage = zeros(length(K_set),N);% 存儲恢復(fù)成功概率%主循環(huán)，遍歷每組(K,M,N)ticfor kk = 1:length(K_set)K = K_set(kk);% 本次稀疏度M_set = K:5:N;%M沒必要全部遍歷，每隔5測試一個就可以了PercentageK = zeros(1,

13、length(M_set);%存儲此稀疏度K下不同M的恢復(fù)成功概率kkfor 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;%得到觀測向量ytheta = CS_ROMP(y,A,K);%恢復(fù)重構(gòu)信號 thetax_

14、r = Psi * theta;% x=Psi * thetaif norm(x_r-x)1e-6%如果殘差小于1e-6則認(rèn)為恢復(fù)成功P = P + 1;endendPercentageK(mm) = P/CNT*100;% 計算恢復(fù)概率endPercentage(kk,1:length(M_set) = PercentageK;endtocsave ROMPMtoPercentage1000 %運行一次不容易，把變量全部存儲下來%繪圖S = -ks;-ko;-kd;-kv;-k*;figure;for kk = 1:length(K_set)K = K_set(kk);M_set = K:5

15、:N;L_Mset = length(M_set);plot(M_set,Percentage(kk,1:L_Mset),S(kk,:);%繪出 x 的恢復(fù)信號hold on;endhold 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);本程序在聯(lián)想ThinkPadE430C筆記本(

16、4GBDDR3內(nèi)存，i5-3210)上運行共耗時871.829395秒(上篇中OMP運行共耗時1171.606254秒)，程序中將所有數(shù)據(jù)均通“save ROMPMtoPercentage1000存儲了下來， 以后可以再對數(shù)據(jù)進(jìn)行分析，只需“l(fā)oad ROMPMtoPercentage1000”即可。程序運行結(jié)果比文獻(xiàn)4的Fig.1要差（OMP仿真結(jié)果比文獻(xiàn)中的要好），原因不詳。我調(diào)用了文獻(xiàn)2中的ROMP函數(shù)，效果和這里的CS ROMP差不多。本程序運行結(jié)果;HEaxooaCDEElLICDECDdHEaxooaCDEElLICDECDd文獻(xiàn)4中的Fig.1：Percentage 寸 站胡si

17、gnals rqver-eciNurnber ofN)Percentage 寸 站胡signals rqver-eciNurnber ofN)Fig. 1 The percentage cf spar泥 fiat signals exactly rectjvtjred by ROMP as a fiinctiDn of the number of meaurmtnts N in dimensLon d = 256 for various levels of s-pusi ty np$#tl 必。Elc5、結(jié)語國內(nèi)引用ROMP文獻(xiàn)時一般為4和5,文獻(xiàn)4更早一些，其實它們有關(guān)ROMP的敘述基本是一

18、樣子的， 下面分別是文獻(xiàn)4和5中的ROMP步驟：文獻(xiàn)4：Regularized Orthogonal Matching Pursuit (ROMP)Input: Measurement vector x E and sparsity level nOutput: Index set I C L ，dInitialized Let the index set 7 = 0 and the Tesidual r jc.Repeat the following steps until r = 0;Identify: Choose a set J of the n biggest coordijiate

19、s in magiutude of the observation vector u =如*廣，or all of its nonzero coordinates, whichever set is smaller.Regularize: Among all subsets J()C J with comparable coordinates2w(0| 2|m(j)| for all i, j 如choose 島 with the maximal energy l.Update: Add the set Jo to the index set: 1 JU 如 and update the re

20、sidual:y = arg min |x z|c； =工一y.zgRj文獻(xiàn)5：其實它們基本是一樣的，除了綜上迭代的條件略有不同，本文的CS_ROMP中都進(jìn)行了考慮。其實若要看ROMP的原理，建議看文獻(xiàn)4的1.4節(jié)：When we are trying to recover the signal u from its measurements x = we can use the observation vector u = *1： as a good local approximation to the signal v. Namely, the observation vector u e

21、ncodes correlations of the measurement vector x with the columns of .Note that 0 is a dictionary, and so since the signal v sparse, x has a sparse representation with respect to the dictionary. By the restricted isometry condition, every n columns form approximately an orthonormn system. Therefore,

22、every n coordinates of th巳 observation vector u look like correlations of the measurement vector x with the orthonormal basis audq therefore, are close in the Euclidean norm to the corresponding n coefficients of u. This is documented in Proposition 3.2 below.The local approximation property sugge&t

23、s to make us-e af the n biggest coordinates of the observation vector w, rather than one biggest coordinate as OMP did. We thus force the selected coordinates to be more regular (i.巳、closer to uniform) by selecting only the coordinates with comparable sizes. To this end, a new regufarizatim step will be needed lo ens.ure that each of these coordinates gets an even share of information. This leads to the following algorithm for sparse recovery:另外，針對ROMP的改進(jìn)算法也很多，后面看些文獻(xiàn)再敘吧。推薦看看文獻(xiàn)67兩篇重構(gòu)算法的綜述性文 獻(xiàn)，本文第一句話就來自文獻(xiàn)6。比較本文的ROMP仿真結(jié)果與上一篇中OMP的仿真結(jié)果可以發(fā)現(xiàn)效果遠(yuǎn)沒有OMP好，當(dāng)然本文文獻(xiàn)4 中的Fig.1比上一篇文獻(xiàn)1中的Fig.1,其中原因不詳。參考文獻(xiàn)