淺談壓縮感知(九):正交匹配追蹤算法OMP


主要內容:

  1. OMP算法介紹
  2. OMP的MATLAB實現
  3. OMP中的數學知識

一、OMP算法介紹

來源:http://blog.csdn.net/scucj/article/details/7467955

1、信號的稀疏表示(sparse representation of signals)

給定一個過完備字典矩陣,其中它的每列表示一種原型信號的原子。給定一個信號y,它可以被表示成這些原子的稀疏線性組合。信號 y 可以被表達為 y = Dx ,或者字典矩陣中所謂過完備性,指的是原子的個數遠遠大於信號y的長度(其長度很顯然是n),即n<<k

2、MP算法(匹配追蹤算法)

2.1 算法描述

作為對信號進行稀疏分解的方法之一,將信號在完備字典庫上進行分解。

假定被表示的信號為y,其長度為n。假定H表示Hilbert空間,在這個空間H里,由一組向量構成字典矩陣D,其中每個向量可以稱為原子(atom),其長度與被表示信號 y 的長度n相同,而且這些向量已作為歸一化處理,即|,也就是單位向量長度為1MP算法的基本思想:從字典矩陣D(也稱為過完備原子庫中),選擇一個與信號 y 最匹配的原子(也就是某列),構建一個稀疏逼近,並求出信號殘差,然后繼續選擇與信號殘差最匹配的原子,反復迭代,信號y可以由這些原子來線性和,再加上最后的殘差值來表示。很顯然,如果殘差值在可以忽略的范圍內,則信號y就是這些原子的線性組合。如果選擇與信號y最匹配的原子?如何構建稀疏逼近並求殘差?如何進行迭代?我們來詳細介紹使用MP進行信號分解的步驟:[1] 計算信號 y 與字典矩陣中每列(原子)的內積,選擇絕對值最大的一個原子,它就是與信號 y 在本次迭代運算中最匹配的。用專業術語來描述:令信號,從字典矩陣中選擇一個最為匹配的原子,滿足r0 表示一個字典矩陣的列索引。這樣,信號 y 就被分解為在最匹配原子的垂直投影分量和殘值兩部分,即:[2]對殘值R1f進行步驟[1]同樣的分解,那么第K步可以得到:

 其中 滿足。可見,經過K步分解后,信號 y 被分解為:,其中

2.2 繼續討論

(1)為什么要假定在Hilbert空間中?Hilbert空間就是定義了完備的內積空。很顯然,MP中的計算使用向量的內積運算,所以在在Hilbert空間中進行信號分解理所當然了。什么是完備的內積空間?篇幅有限就請自己搜索一下吧。

(2)為什么原子要事先被歸一化處理了,即上面的描述。內積常用於計算一個矢量在一個方向上的投影長度,這時方向的矢量必須是單位矢量。MP中選擇最匹配的原子是,是選擇內積最大的一個,也就是信號(或是殘值)在原子(單位的)垂直投影長度最長的一個,比如第一次分解過程中,投影長度就是,三個向量,構成一個三角形,且正交(不能說垂直,但是可以想象二維空間這兩個矢量是垂直的)。

(3)MP算法是收斂的,因為正交,由這兩個可以得出,得出每一個殘值比上一次的小,故而收斂。

2.3 MP算法的缺點

如上所述,如果信號(殘值)在已選擇的原子進行垂直投影是非正交性的,這會使得每次迭代的結果並不少最優的而是次最優的,收斂需要很多次迭代。舉個例子說明一下:在二維空間上,有一個信號 y D=[x1, x2]來表達,MP算法迭代會發現總是在x1x2上反復迭代,即,這個就是信號(殘值)在已選擇的原子進行垂直投影的非正交性導致的。再用嚴謹的方式描述[1]可能容易理解:Hilbert空間H中,,定義,就是它是這些向量的張成中的一個,MP構造一種表達形式:;這里的Pvf表示 fV上的一個正交投影操作,那么MP算法的第 k 次迭代的結果可以表示如下(前面描述時信號為y,這里變成f了,請注意)

如果  是最優的k項近似值,當且僅當。由於MP僅能保證,所以一般情況下是次優的。這是什么意思呢?k個項的線性表示,這個組合的值作為近似值,只有在第k個殘差和正交,才是最優的。如果第k個殘值與正交,意味這個殘值與fk的任意一項都線性無關,那么第k個殘值在后面的分解過程中,不可能出現fk中已經出現的項,這才是最優的。而一般情況下,不能滿足這個條件,MP一般只能滿足第k個殘差和xk正交,這也就是前面為什么提到"信號(殘值)在已選擇的原子進行垂直投影是非正交性的"的原因。如果第k個殘差和fk不正交,那么后面的迭代還會出現fk中已經出現的項,很顯然fk就不是最優的,這也就是為什么說MP收斂就需要更多次迭代的原因。不是說MP一定得到不到最優解,而且其前面描述的特性導致一般得到不到最優解而是次優解。那么,有沒有辦法讓第k個殘差與正交,方法是有的,這就是下面要談到的OMP算法。

3、OMP算法

3.1 算法描述

OMP算法的改進之處在於:在分解的每一步對所選擇的全部原子進行正交化處理,這使得在精度要求相同的情況下,OMP算法的收斂速度更快。

那么在每一步中如何對所選擇的全部原子進行正交化處理呢?在正式描述OMP算法前,先看一點基礎思想。

先看一個 k  階模型,表示信號 f 經過 k 步分解后的情況,似乎很眼熟,但要注意它與MP算法不同之處,它的殘值與前面每個分量正交,這就是為什么這個算法多了一個正交的原因,MP中僅與最近選出的的那一項正交。

(1)

 k + 1 階模型如下:

(2)

應用 k + 1階模型減去k 階模型,得到如下:

(3)

我們知道,字典矩陣D的原子是非正交的,引入一個輔助模型,它是表示對前k個項的依賴,描述如下:

(4)

和前面描述類似,span(x1, ...xk)之一上的正交投影操作,后面的項是殘值。這個關系用數學符號描述:

請注意,這里的 a b 的上標表示第 k 步時的取值。

(4)帶入(3)中,有:

5

如果一下兩個式子成立,(5)必然成立。

(6)

(7)

,有

其中

ak的值是由求法很簡單,通過對(7)左右兩邊添加作內積消減得到:

后邊的第二項因為它們正交,所以為0,所以可以得出ak的第一部分。對於,在(4)左右兩邊中與作內積,可以得到ak的第二部分。

對於(4),可以求出,求的步驟請參見參考文件的計算細節部分。為什么這里不提,因為后面會介紹更簡單的方法來計算。

3.2 收斂性證明

通過(7),由於正交,將兩個殘值移到右邊后求二范的平方,並將ak的值代入可以得到:

可見每一次殘差比上一次殘差小,可見是收斂的。

3.3 算法步驟

整個OMP算法的步驟如下:

由於有了上面的來龍去脈,這個算法就相當好理解了。

到這里還不算完,后來OMP的迭代運算用另外一種方法可以計算得知,有位同學的論文[2]描述就非常好,我就直接引用進來:

對比中英文描述,本質都是一樣,只是有細微的差別。這里順便貼出網一哥們寫的OMP算法的代碼,源出處不得而知,共享給大家。

二、OMP的MATLAB實現

 1、一維信號重建

代碼:

%  1-D信號壓縮傳感的實現(正交匹配追蹤法Orthogonal Matching Pursuit)
%  測量數M>=K*log(N/K),K是稀疏度,N信號長度,可以近乎完全重構
%  編程人--香港大學電子工程系 沙威  Email: wsha@eee.hku.hk
%  編程時間:2008年11月18日
%  文檔下載: http://www.eee.hku.hk/~wsha/Freecode/freecode.htm 
%  參考文獻:Joel A. Tropp and Anna C. Gilbert 
%  Signal Recovery From Random Measurements Via Orthogonal Matching
%  Pursuit,IEEE TRANSACTIONS ON INFORMATION THEORY, VOL. 53, NO. 12,
%  DECEMBER 2007.

clc;clear

%%  1. 時域測試信號生成
K=7;      %  稀疏度(做FFT可以看出來)
N=256;    %  信號長度
M=64;     %  測量數(M>=K*log(N/K),至少40,但有出錯的概率)
f1=50;    %  信號頻率1
f2=100;   %  信號頻率2
f3=200;   %  信號頻率3
f4=400;   %  信號頻率4
fs=800;   %  采樣頻率
ts=1/fs;  %  采樣間隔
Ts=1:N;   %  采樣序列
x=0.3*cos(2*pi*f1*Ts*ts)+0.6*cos(2*pi*f2*Ts*ts)+0.1*cos(2*pi*f3*Ts*ts)+0.9*cos(2*pi*f4*Ts*ts);  %  完整信號,由4個信號疊加而來

%%  2.  時域信號壓縮傳感
Phi=randn(M,N);                                   %  測量矩陣(高斯分布白噪聲)64*256的扁矩陣,Phi也就是文中說的D矩陣
s=Phi*x.';                                        %  獲得線性測量 ,s相當於文中的y矩陣

%%  3.  正交匹配追蹤法重構信號(本質上是L_1范數最優化問題)
%匹配追蹤:找到一個其標記看上去與收集到的數據相關的小波;在數據中去除這個標記的所有印跡;不斷重復直到我們能用小波標記“解釋”收集到的所有數據。

m=2*K;                                            %  算法迭代次數(m>=K),設x是K-sparse的
Psi=fft(eye(N,N))/sqrt(N);                        %  傅里葉正變換矩陣
T=Phi*Psi';                                       %  恢復矩陣(測量矩陣*正交反變換矩陣)

hat_y=zeros(1,N);                                 %  待重構的譜域(變換域)向量                     
Aug_t=[];                                         %  增量矩陣(初始值為空矩陣)
r_n=s;                                            %  殘差值

for times=1:m;                                    %  迭代次數(有噪聲的情況下,該迭代次數為K)
    for col=1:N;                                  %  恢復矩陣的所有列向量
        product(col)=abs(T(:,col)'*r_n);          %  恢復矩陣的列向量和殘差的投影系數(內積值) 
    end
    [val,pos]=max(product);                       %  最大投影系數對應的位置,即找到一個其標記看上去與收集到的數據相關的小波
    Aug_t=[Aug_t,T(:,pos)];                       %  矩陣擴充    
    
    T(:,pos)=zeros(M,1);                          %  選中的列置零(實質上應該去掉,為了簡單我把它置零),在數據中去除這個標記的所有印跡
    aug_y=(Aug_t'*Aug_t)^(-1)*Aug_t'*s;           %  最小二乘,使殘差最小
    r_n=s-Aug_t*aug_y;                            %  殘差
    pos_array(times)=pos;                         %  紀錄最大投影系數的位置
end
hat_y(pos_array)=aug_y;                           %  重構的譜域向量
hat_x=real(Psi'*hat_y.');                         %  做逆傅里葉變換重構得到時域信號

%%  4.  恢復信號和原始信號對比
figure(1);
hold on;
plot(hat_x,'k.-')                                 %  重建信號
plot(x,'r')                                       %  原始信號
legend('Recovery','Original')
norm(hat_x.'-x)/norm(x)                           %  重構誤差

結果:

2、二維信號重建

代碼:

%  本程序實現圖像LENA的壓縮傳感
%  程序作者:沙威,香港大學電氣電子工程學系,wsha@eee.hku.hk
%  算法采用正交匹配法,參考文獻 Joel A. Tropp and Anna C. Gilbert 
%  Signal Recovery From Random Measurements Via Orthogonal Matching
%  Pursuit,IEEE TRANSACTIONS ON INFORMATION THEORY, VOL. 53, NO. 12,
%  DECEMBER 2007.
%  該程序沒有經過任何優化

function Wavelet_OMP

clc;clear

%  讀文件
X=imread('lena256.bmp');
X=double(X);
[a,b]=size(X);

%  小波變換矩陣生成
ww=DWT(a);

%  小波變換讓圖像稀疏化(注意該步驟會耗費時間,但是會增大稀疏度)
X1=ww*sparse(X)*ww';
% X1=X;
X1=full(X1);

%  隨機矩陣生成
M=190;
R=randn(M,a);
% R=mapminmax(R,0,255);
% R=round(R);

%  測量值
Y=R*X1;

%  OMP算法
%  恢復矩陣
X2=zeros(a,b); 
%  按列循環
for i=1:b 
    %  通過OMP,返回每一列信號對應的恢復值(小波域)
    rec=omp(Y(:,i),R,a);
    %  恢復值矩陣,用於反變換
    X2(:,i)=rec;
end


%  原始圖像
figure(1);
imshow(uint8(X));
title('原始圖像');

%  變換圖像
figure(2);
imshow(uint8(X1));
title('小波變換后的圖像');

%  壓縮傳感恢復的圖像
figure(3);
%  小波反變換
X3=ww'*sparse(X2)*ww; 
% X3=X2;
X3=full(X3);
imshow(uint8(X3));
title('恢復的圖像');

%  誤差(PSNR)
%  MSE誤差
errorx=sum(sum(abs(X3-X).^2));        
%  PSNR
psnr=10*log10(255*255/(errorx/a/b))   


%  OMP的函數
%  s-測量;T-觀測矩陣;N-向量大小
function hat_y=omp(s,T,N)
Size=size(T);                                     %  觀測矩陣大小
M=Size(1);                                        %  測量
hat_y=zeros(1,N);                                 %  待重構的譜域(變換域)向量                     
Aug_t=[];                                         %  增量矩陣(初始值為空矩陣)
r_n=s;                                            %  殘差值

for times=1:M;                                  %  迭代次數(稀疏度是測量的1/4)
    for col=1:N;                                  %  恢復矩陣的所有列向量
        product(col)=abs(T(:,col)'*r_n);          %  恢復矩陣的列向量和殘差的投影系數(內積值) 
    end
    [val,pos]=max(product);                       %  最大投影系數對應的位置
    Aug_t=[Aug_t,T(:,pos)];                       %  矩陣擴充
    T(:,pos)=zeros(M,1);                          %  選中的列置零(實質上應該去掉,為了簡單我把它置零)
    aug_y=(Aug_t'*Aug_t)^(-1)*Aug_t'*s;           %  最小二乘,使殘差最小
    r_n=s-Aug_t*aug_y;                            %  殘差
    pos_array(times)=pos;                         %  紀錄最大投影系數的位置
    
    if (norm(r_n)<0.9)                            %  殘差足夠小
        break;
    end
end
hat_y(pos_array)=aug_y;                           %  重構的向量

結果:

三、OMP算法中的數學知識

算法會有一些線性代數中的概念和知識,主要是關於子空間、最小二乘、投影矩陣等,詳情請參考

最小二乘的幾何意義及投影矩陣


免責聲明!

本站轉載的文章為個人學習借鑒使用,本站對版權不負任何法律責任。如果侵犯了您的隱私權益,請聯系本站郵箱yoyou2525@163.com刪除。



 
粵ICP備18138465號   © 2018-2025 CODEPRJ.COM