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

1、壓縮感知重構(gòu)算法之子空間追蹤(SP)如果掌握了壓縮采樣匹配追蹤(CoSaMP)后，再去學(xué)習(xí)子空間追蹤 (Subspace Pursuit)是一件非常簡(jiǎn)單的事情,因?yàn)樗鼈儙缀跏峭耆粯拥摹P的提出時(shí)間比 CoSaMP提出時(shí)間略晚，首個(gè)論文版本是參考文獻(xiàn) 1,后來(lái)更新了兩次，最后在 IEEE Transactions on Information Theory 發(fā)表2。從算法角度來(lái)講，SP與CoSaMP差別非常小，這一點(diǎn)作者也意識(shí)到 了，在文獻(xiàn)1首頁(yè)的左下角就有注釋?zhuān)? AL die time of writing this maiiuript. tl)c audiors became awue

2、 of the related werk by J. Tmpp* D. Needell and R, Vershynin 111 where similar reconstruction alguriilmih are designed. Our results dcvuluped inckpcndenUy, and we believe (hut there arc significant differences in these two proposed rcconstmction ripproaehcs.在文獻(xiàn)2第2頁(yè)提到了 SP與CoSaMP的具體不同：Al the time of w

3、riting of this paper, the authors became aware of the related work by L Tropp. D. Needell, and R, Vershynin 12, describing a similar rcconstniction algorithtiL The main difference between the SP algorithm and the CoSAMP algorithm of 12J is in the manner in which new candidates arc added to the list.

4、 In each iteration, in the SP algorithm, only K new candidates are added, while the CoSAMP algorithm adds 2K vectors. This makes the SP algorithm computationally mare efficient, but the underlying analysis more complex. In addition, (he restricted isomeiry ccrnKtanl for which the SP algorithm is gua

5、r an teed to converge is larger than the one presented in 12. Finally. this paper also contains an analysis of the number of iterations needed for reconstruction of a sparse signal (see Theorem 6 for details), for which there is no counterpart in the CoSAMP study.從上面可以知道，SP 與 CoSaMP 主要區(qū)別在于 Ineach it

6、eration, in the SP algorithm, only K new candidates areadded , while theCoSAMP algorithm adds 2K vectors.，即“SP 每次選擇 K 個(gè)原子，而 CoSaMP 則選擇 2K 個(gè)原子;這樣帶來(lái)的好處是 “This makes the SP algorithm computationally moreefficient, ?！币韵率俏墨I(xiàn)2中的給出的SP算法流程:Algorithm 1 Subspace Pursuit AlgorithmInput: K、典 yInitialization:7* =

7、 indices corresponding to the largest magniude entries in the vector y.2)婢=resid冢).iteration: At the/th iteration, go through the following steps.1)型=T；-i |J K indices corresponding to the largest magnitude entries in the vectnr 中電/1 .Set xp =%V7 ; K indices corresponding to the largest magnitude el

8、ements of 的J.二 resid (If |j|2 “必TJLlet= T$-i and quit theiteration,Ouipui:1) The estimated signal k satisfying 仝1，.少歡-/ = 0 and xT 明妙這個(gè)算法流程的初始化（Initialization）其實(shí)就是類(lèi)似于 CoSaMP的第1次迭代，注意第（1）步中選擇了 K個(gè)原子：“K indices correspo nding to the largest magnitude entries ,在 CoSaMP 里這里要選擇 2K 個(gè)最大的原子，后面的其 它流程都一樣。這里第（5

9、）步增加了一個(gè)停止迭代的條件：當(dāng)殘差經(jīng)過(guò)迭代后卻變大了的時(shí)候就停止迭代。不只是SP作者認(rèn)識(shí)到了自己的算法與CoSaMP的高度相似性，CoSaMP的作者也同樣關(guān)注到了SP算法，在文獻(xiàn)3中就提到：After initially presented this work. Did and Milenkovic developed an Algorithm called Subspace1 Pursuit that is very similar to CoSaMP. They established that their algorithm offers perfamiancc guarantees

10、analogous with those for CoSaMP, Src 10 for details.文獻(xiàn)3是CoSaMP原始提出文獻(xiàn)的第 2個(gè)版本，文獻(xiàn)3的早期版本4是沒(méi)有提及SP算法的。鑒于SP與CoSaMP如此相似，這里不就再單獨(dú)給出SP的步驟了，參考?jí)嚎s感知重構(gòu)算法之壓縮采樣匹配追蹤（CoSaMP） ,只需將第（2）步中的2K改為K即可。引用文獻(xiàn)5的3,5節(jié)中的幾句話(huà)：貪婪類(lèi)算法雖然復(fù)雜度低運(yùn)行速度快，但其重構(gòu)精度卻不如BP類(lèi)算法,為了尋求復(fù)雜度和精度更好地折中，SP算法應(yīng)運(yùn)而生”，“SPI法與CoSaMP算法一樣其基本思想也是借用回溯的思想，在每步迭代過(guò)程中重新估計(jì)所有候選者的可信

11、賴(lài)性：SFW法與CoSaMP算法有著類(lèi)似的性質(zhì)與優(yōu)缺點(diǎn)”子空間追蹤代碼可實(shí)現(xiàn)如下(CS_SP.m),通過(guò)對(duì)比可以知道該程序與CoSaMP的代碼基本完全一致。本代碼未考慮文獻(xiàn)2中的給出的SP算法流程的第(5)步。代碼可參見(jiàn)參考文獻(xiàn)6中的Demo_CS_SP.m 。plain view plaincopyfunction theta = CS_SP( y,A,K )%CS_SP Summary of this function goes here%Version: 1.0 written by jbb0523 2015-05-01% Detailed explanation goes here%

12、y = Phi * x% x = Psi * theta% y = Phi*Psi * theta% 令 A = Phi*Psi,貝U y=A*theta% K is the sparsity level%現(xiàn)在已知y和A,求theta% Reference:Dai W , Milenkovic O . Subspace pursuit for compressive sensing% signal reconstructionJ. IEEE Transactions on Information Theory,% 2009 , 55(5) : 2230-2249.y_rows,y_column

13、s = size(y);if y_rowsy_columnsy = y;%y should be a column vectorendM,N = size(A);%傳感矢!陣 A 為 M*N矩陣theta = zeros(N,1);%用來(lái)存儲(chǔ)恢復(fù)的theta( 列向量)Pos_theta = ;%用來(lái)迭代過(guò)程中存儲(chǔ) A被選擇的列序號(hào)r_n = y;% 初始化殘差(residual) 為 yfor kk=1:K% 最多迭代K次%(1) Identificationproduct = A*r_n;%傳感矢1陣A各列與殘差的內(nèi)積val,pos=sort(abs(product),descend);J

14、s = pos(1:K);%選出內(nèi)積值最大的K列%(2) Support MergerIs = union(Pos_theta,Js);%Pos_theta與Js并集%(3) Estimation%At的行數(shù)要大于列數(shù)，此為最小二乘的基礎(chǔ) (列線(xiàn)性無(wú)關(guān))if length(Is)=MAt = A(:,Is);%將A的這幾列組成矩陣 Atelse%At的列數(shù)大于行數(shù)，列必為線(xiàn)性相關(guān)的,At*At 將不可逆break;%跳出 for 循環(huán)end%y=At*theta,以下求 theta 的最小二乘解(Least Square)theta_ls = (At*At)A(-1)*At*y;%最小二乘解%

15、(4) Pruningval,pos=sort(abs(theta_ls),descend);%(5) Sample UpdatePos_theta = Is(pos(1:K);theta_ls = theta_ls(pos(1:K);%At(:,pos(1:K)*theta_ls是 y 在 At(:,pos(1:K)列空間上的正交投影r_n = y - At(:,pos(1:K)*theta_ls;%更新殘差if norm(r_n)1e-6%Repeat thesteps until r=0break;%跳出 for 循環(huán)endendtheta(Pos_theta)=theta_ls;%恢復(fù)

16、出的 thetaend鑒于SP與CoSaMP的極其相似性，這里就不再給出單次重構(gòu)和測(cè)量數(shù)M與重構(gòu)成功概率關(guān)系曲線(xiàn)繪制例程代碼了，只需將 CoSaMP中調(diào)用CS_CoSaMP函數(shù)的部分改為調(diào)用CS_SP即可，無(wú)須任何其它改動(dòng)。這里給出對(duì)比兩種重構(gòu)算法所繪制的測(cè)量數(shù)M與重構(gòu)成功概率關(guān)系曲線(xiàn)的例程代碼，只有這樣才可以看出兩種算法的重構(gòu)性能優(yōu)劣，以下是在分別運(yùn)行完SP與CoSaMP的測(cè)量數(shù)M與重構(gòu)成功概率關(guān)系曲線(xiàn)繪制例程代碼的基礎(chǔ)上，即已經(jīng)存儲(chǔ)了數(shù)據(jù) CoSaMPMtoPercentage1000.mat 和 SPMtoPercentage1000.mat :plain view plaincopyc

17、lear all;close all;clc;load CoSaMPMtoPercentage1000;PercentageCoSaMP = Percentage;load SPMtoPercentage1000;PercentageSP = Percentage;S1 = -ks;-ko;-kd;-kv;-k*;S2 = -rs;-ro;-rd;-rv;-r*;figure;for kk = 1:length(K_set)K = K_set(kk);M_set = 2*K:5:N;L_Mset = length(M_set);plot(M_set,PercentageCoSaMP(kk,1:

18、L_Mset),S1(kk,:);%繪出 x 的恢復(fù)信號(hào)hold on;plot(M_set,PercentageSP(kk,1:L_Mset),S2(kk,:);%繪出 x 的恢復(fù)信號(hào)endhold off;xlim(0 256);legend(CoSaK=4,SPK=4,CoSaK=12,SPK=12,CoSaK=20,.SPK=20,CoSaK=28,SPK=28,CoSaK=36,SPK=36);xlabel(Number of measurements(M);ylabel(Percentage recovered);title(Percentage of input signals

19、recovered correctly(N=256)(Gaussian);運(yùn)彳譚果如下:Percerrtage of input signals recovered correctly(N=256)(Gajssianj00903070605040302010口f fy /鏟現(xiàn)聚a*”售工7呼觸摩岫率就值闔tjwte日一CoSaK4 .1f TJ-SPK=4I iJ-e-CQSaK=121TJ-SPK=12JJ8-CoSaK=20-SPK=20,tIj J J7-CoSaK=2611 fI j1V-SPK=2BI1CoS 水 =36-t-SPK=3E1-/I i1L-J1wa150200260Mumber of measuremeritsM可以發(fā)現(xiàn)在H較小時(shí)SN等好于8s洲R為聞變方寸二者重構(gòu)性能幾乎一樣。參考文獻(xiàn)：1Dai W , Milenkovic O. Subspace pursuitfor compressive sensing: Closing the gap between performance andcomplexity. HYPERLINK /pdf/0803.0811v1.pdf /pdf/0803.08