分段正交匹配追蹤算法StOMP

壓縮感知重構(gòu)算法之分段正交匹配追蹤(StOMP)分段正交匹配追蹤(StagewiseOMP)或者翻譯為逐步正交匹配追蹤，它是OMP另一種改進(jìn)算法，每次迭代可以選擇多個(gè)原子。此算法的輸入?yún)?shù)中沒(méi)有信號(hào)稀疏度K，因此相比于ROMP及CoSaMP有獨(dú)到的優(yōu)勢(shì)。0、符號(hào)說(shuō)明如下：壓縮觀測(cè)y=Φx，其中y為觀測(cè)所得向量M×1，x為原信號(hào)N×1（M<<N）。x一般不是稀疏的，但在某個(gè)變換域Ψ是稀疏的，即x=Ψθ，其中θ為K稀疏的，即θ只有K個(gè)非零項(xiàng)。此時(shí)y=ΦΨθ，令A(yù)=ΦΨ，則y=Aθ。(1)y為觀測(cè)所得向量，大小為M×1(2)x為原信號(hào)，大小為N×1(3)θ為K稀疏的，是信號(hào)在x在某變換域的稀疏表示(4)Φ稱為觀測(cè)矩陣、測(cè)量矩陣、測(cè)量基，大小為M×N(5)Ψ稱為變換矩陣、變換基、稀疏矩陣、稀疏基、正交基字典矩陣，大小為N×N(6)A稱為測(cè)度矩陣、傳感矩陣、CS信息算子，大小為M×N上式中，一般有K<<M<<N，后面三個(gè)矩陣各個(gè)文獻(xiàn)的叫法不一，以后我將Φ稱為測(cè)量矩陣、將Ψ稱為稀疏矩陣、將A稱為傳感矩陣。注意：這里的稀疏表示模型為x=Ψθ，所以傳感矩陣A=ΦΨ；而有些文獻(xiàn)中稀疏模型為θ=Ψx，而一般Ψ為Hermite矩陣(實(shí)矩陣時(shí)稱為正交矩陣)，所以Ψ-1=ΨH(實(shí)矩陣時(shí)為Ψ-1=ΨT)，即x=ΨHθ，所以傳感矩陣A=ΦΨH，例如沙威的OMP例程中就是如此。1、StOMP重構(gòu)算法流程：2、分段正交匹配追蹤(StOMP)Matlab代碼(CS_StOMP.m)代碼參考了文獻(xiàn)[4]中的SolveStOMP.m，也可參考文獻(xiàn)[5]中的StOMP.m。其實(shí)文獻(xiàn)[4]是斯坦福的SparseLab中的一個(gè)函數(shù)而已，鏈接為/，最新版本為2.1，SolveStOMP.m在目錄SparseLab21-Core\SparseLab2.1-Core\Solvers里面。function[theta]=CS_StOMP(y,A,S,ts)%CS_StOMPSummaryofthisfunctiongoeshere%Version:1.0writtenbyjbb0523@2015-04-29%Detailedexplanationgoeshere%y=Phi*x%x=Psi*theta%y=Phi*Psi*theta%

令

A

=

Phi*Psi,

則y=A*theta

%

S

is

the

maximum

number

of

StOMP

iterations

to

perform

%

ts

is

the

threshold

parameter

%

現(xiàn)在已知y和A，求theta

%

Reference:Donoho

D

L，Tsaig

Y，Drori

I，Starck

J

L．Sparse

solution

of

%

underdetermined

linear

equations

by

stagewise

orthogonal

matching

%

pursuit[J]．IEEE

Transactions

on

Information

Theory，2012，58(2)：1094—1121

if

nargin

<

4

ts

=

2.5;%ts范圍[2,3],默認(rèn)值為2.5

end

if

nargin

<

3

S

=

10;%S默認(rèn)值為10

end

[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);%用來(lái)存儲(chǔ)恢復(fù)的theta(列向量)

Pos_theta

=

[];%用來(lái)迭代過(guò)程中存儲(chǔ)A被選擇的列序號(hào)

r_n

=

y;%初始化殘差(residual)為y

for

ss=1:S%最多迭代S次

product

=

A'*r_n;%傳感矩陣A各列與殘差的內(nèi)積

sigma

=

norm(r_n)/sqrt(M);%參見參考文獻(xiàn)第3頁(yè)Remarks(3)

Js

=

find(abs(product)>ts*sigma);%選出大于閾值的列

Is

=

union(Pos_theta,Js);%Pos_theta與Js并集

if

length(Pos_theta)

==

length(Is)

if

ss==1

theta_ls

=

0;%防止第1次就跳出導(dǎo)致theta_ls無(wú)定義

end

break;%如果沒(méi)有新的列被選中則跳出循環(huán)

end

%At的行數(shù)要大于列數(shù)，此為最小二乘的基礎(chǔ)(列線性無(wú)關(guān))

if

length(Is)<=M

Pos_theta

=

Is;%更新列序號(hào)集合

At

=

A(:,Pos_theta);%將A的這幾列組成矩陣At

else%At的列數(shù)大于行數(shù)，列必為線性相關(guān)的,At'*At將不可逆

if

ss==1

theta_ls

=

0;%防止第1次就跳出導(dǎo)致theta_ls無(wú)定義

end

break;%跳出for循環(huán)

end

%y=At*theta，以下求theta的最小二乘解(Least

Square)

theta_ls

=

(At'*At)^(-1)*At'*y;%最小二乘解

%At*theta_ls是y在At列空間上的正交投影

r_n

=

y

-

At*theta_ls;%更新殘差

if

norm(r_n)<1e-6%Repeat

the

steps

until

r=0

break;%跳出for循環(huán)

end

end

theta(Pos_theta)=theta_ls;%恢復(fù)出的theta

end

3、StOMP單次重構(gòu)測(cè)試代碼

以下測(cè)試代碼基本與OMP單次重構(gòu)測(cè)試代碼一樣，除了調(diào)用CS_StOMP之外，一定要注意這里的測(cè)量矩陣Phi=randn(M,N)/sqrt(M)，一定一定?。?！[plain]

\o"viewplain"viewplain\o"copy"copy%壓縮感知重構(gòu)算法測(cè)試

clear

all;close

all;clc;

M

=

64;%觀測(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)/sqrt(M);%測(cè)量矩陣為高斯矩陣

A

=

Phi

*

Psi;%傳感矩陣

y

=

Phi

*

x;%得到觀測(cè)向量y

%%

恢復(fù)重構(gòu)信號(hào)x

tic

theta

=

CS_StOMP(y,A);

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ù)殘差

運(yùn)行結(jié)果如下：（信號(hào)為隨機(jī)生成，所以每次結(jié)果均不一樣）

1）圖：

2）Command

windows

Elapsedtimeis0.067904seconds.

恢復(fù)殘差：

ans=

6.1267e-0154、門限參數(shù)ts、測(cè)量數(shù)M與重構(gòu)成功概率關(guān)系曲線繪制例程代碼

因?yàn)槲墨I(xiàn)[1]中對(duì)門限參數(shù)ts給出的是一個(gè)取值范圍，所以有必要仿真ts取不同值時(shí)的重構(gòu)效果，因此以下的代碼雖然是基于OMP相應(yīng)的測(cè)試代碼修改的，但相對(duì)來(lái)說(shuō)改動(dòng)較大。[plain]

\o"viewplain"viewplain\o"copy"copyclear

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

ts_set

=

2:0.2:3;

K_set

=

[4,12,20,28,36];%信號(hào)x的稀疏度集合

Percentage

=

zeros(N,length(K_set),length(ts_set));%存儲(chǔ)恢復(fù)成功概率

%%

主循環(huán)，遍歷每組(ts,K,M,N)

tic

for

tt

=

1:length(ts_set)

ts

=

ts_set(tt);

for

kk

=

1:length(K_set)

K

=

K_set(kk);%本次稀疏度

%M沒(méi)必要全部遍歷，每隔5測(cè)試一個(gè)就可以了

M_set=2*K:5:N;

PercentageK

=

zeros(1,length(M_set));%存儲(chǔ)此稀疏度K下不同M的恢復(fù)成功概率

for

mm

=

1:length(M_set)

M

=

M_set(mm);%本次觀測(cè)值個(gè)數(shù)

fprintf('ts=%f,K=%d,M=%d\n',ts,K,M);

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)/sqrt(M);%測(cè)量矩陣為高斯矩陣

A

=

Phi

*

Psi;%傳感矩陣

y

=

Phi

*

x;%得到觀測(cè)向量y

theta

=

CS_StOMP(y,A,10,ts);%恢復(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(1:length(M_set),kk,tt)

=

PercentageK;

end

end

toc

save

StOMPMtoPercentage1000

%運(yùn)行一次不容易，把變量全部存儲(chǔ)下來(lái)

%%

繪圖

for

tt

=

1:length(ts_set)

S

=

['-ks';'-ko';'-kd';'-kv';'-k*'];

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,Percentage(1:L_Mset,kk,tt),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,ts=',...

num2str(ts_set(tt)),')(Gaussian)']);

end

for

kk

=

1:length(K_set)

K

=

K_set(kk);

M_set=2*K:5:N;

L_Mset

=

length(M_set);

S

=

['-ks';'-ko';'-kd';'-kv';'-k*';'-k+'];

figure;

for

tt

=

1:length(ts_set)

plot(M_set,Percentage(1:L_Mset,kk,tt),S(tt,:));%繪出x的恢復(fù)信號(hào)

hold

on;

end

hold

off;

xlim([0

256]);

legend('ts=2.0','ts=2.2','ts=2.4','ts=2.6','ts=2.8','ts=3.0');

xlabel('Number

of

measurements(M)');

ylabel('Percentage

recovered');

title(['Percentage

of

input

signals

recovered

correctly(N=256,K=',...

num2str(K),')(Gaussian)']);

end

本程序在聯(lián)想ThinkPadE430C筆記本（4GBDDR3內(nèi)存，i5-3210）上運(yùn)行共耗時(shí)4707.513276秒，程序中將所有數(shù)據(jù)均通過(guò)“saveStOMPMtoPercentage1000”存儲(chǔ)了下來(lái)，以后可以再對(duì)數(shù)據(jù)進(jìn)行分析，只需“l(fā)oadStOMPMtoPercentage1000”即可。

程序運(yùn)行結(jié)束會(huì)出現(xiàn)6+5=11幅圖，前6幅圖分別是ts分別為2.0、2.2、2.4、2.6、2.8和3.0時(shí)的測(cè)量數(shù)M與重構(gòu)成功概率關(guān)系曲線（類似于OMP此部分，這里只是對(duì)每一個(gè)不同的ts畫出一幅圖），后5幅圖是分別將稀疏度K為4、12、20、28、32時(shí)將六種ts取值的測(cè)量數(shù)M與重構(gòu)成功概率關(guān)系曲線繪制在一起以比較ts對(duì)重構(gòu)結(jié)果的影響。

對(duì)于前6幅圖這里只給出ts=2.4時(shí)的曲線圖：

對(duì)于后5幅圖這里全部給出，為了清楚地看出ts的影響，這里把圖的橫軸拉伸：

通過(guò)對(duì)比可以看出，總體上講ts=2.4或ts=2.6時(shí)效果較好，較大和較小重構(gòu)效果都會(huì)降低，這里由于沒(méi)有ts=2.5的情況，但我們推測(cè)ts=2.5應(yīng)該是一個(gè)比較好的值，因此一般默認(rèn)取為2.5即可。5、結(jié)語(yǔ)

有關(guān)StOMP的流程圖可參見文獻(xiàn)[1]的Fig.1：

有關(guān)StOMP門限的選取在文獻(xiàn)[1]中也有提及：關(guān)于這個(gè)門限的來(lái)源文獻(xiàn)[1]有也有一個(gè)推導(dǎo)，注意推導(dǎo)過(guò)程中的N(0,1/n)：

作者在文獻(xiàn)[1]中提出StOMP，這篇文章的發(fā)表時(shí)間是2012年，但看一下這篇文章的左下角會(huì)發(fā)現(xiàn)一個(gè)問(wèn)題：注意，文章在2006-04-05就投稿了，直到2011-08-17修回并被接受，然后2012年才發(fā)表。也就是說(shuō)審稿就審了五年多，按說(shuō)文章第一作者是大牛，雖說(shuō)IEEETransactionsonInformationTheory是一個(gè)頂級(jí)期刊，但對(duì)DonohoDL來(lái)說(shuō)也應(yīng)該不算是什么難事，不知道為什么會(huì)出現(xiàn)這種現(xiàn)象。當(dāng)然，英文文獻(xiàn)里有個(gè)有趣的現(xiàn)