正交匹配追踪算法OMP
淺談壓縮感知（九）：正交匹配追蹤算法OMP
</h1><div class="clear"></div><div class="postBody">主要內容：
一、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相同，而且這些向量已作為歸一化處理，即|，也就是單位向量長度為1。MP算法的基本思想：從字典矩陣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算法迭代會發現總是在x1和x2上反復迭代，即，這個就是信號(殘值)在已選擇的原子進行垂直投影的非正交性導致的。再用嚴謹的方式描述[1]可能容易理解:在Hilbert空間H中，，，定義，就是它是這些向量的張成中的一個，MP構造一種表達形式：;這里的Pvf表示 f在V上的一個正交投影操作，那么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>=Klog(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.3cos(2pif1Tsts)+0.6cos(2pif2Tsts)+0.1cos(2pif3Tsts)+0.9cos(2pif4Ts*ts); % 完整信號,由4個信號疊加而來
%% 2. 時域信號壓縮傳感
Phi=randn(M,N); % 測量矩陣(高斯分布白噪聲)64256的扁矩陣，Phi也就是文中說的D矩陣
s=Phix.’; % 獲得線性測量 ，s相當于文中的y矩陣
%% 3. 正交匹配追蹤法重構信號(本質上是L_1范數最優化問題)
%匹配追蹤：找到一個其標記看上去與收集到的數據相關的小波；在數據中去除這個標記的所有印跡；不斷重復直到我們能用小波標記“解釋”收集到的所有數據。
m=2K; % 算法迭代次數(m>=K)，設x是K-sparse的
Psi=fft(eye(N,N))/sqrt(N); % 傅里葉正變換矩陣
T=PhiPsi’; % 恢復矩陣(測量矩陣*正交反變換矩陣)
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)]; % 矩陣擴充
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?;
% 測量值
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=10log10(255255/(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_taug_y; % 殘差
pos_array(times)=pos; % 紀錄最大投影系數的位置
end
hat_y(pos_array)=aug_y; % 重構的向量
結果：
三、OMP算法中的數學知識
算法會有一些線性代數中的概念和知識，主要是關于子空間、最小二乘、投影矩陣等，詳情請參考
最小二乘的幾何意義及投影矩陣
分類: 壓縮感知compressive sensing, 數學Mathematics <div id="blog_post_info"> 好文要頂 關注我 收藏該文 AndyJee關注 - 0
粉絲 - 324 +加關注 4 0 <div class="clear"></div> <div id="post_next_prev"><a href="https://www.cnblogs.com/AndyJee/p/5045668.html" class="p_n_p_prefix">? </a> 上一篇： <a href="https://www.cnblogs.com/AndyJee/p/5045668.html" title="發布于 2015-12-14 16:29">淺談壓縮感知（八）：兩篇科普文章</a> <br> <a href="https://www.cnblogs.com/AndyJee/p/5048235.html" class="p_n_p_prefix">? </a> 下一篇： <a href="https://www.cnblogs.com/AndyJee/p/5048235.html" title="發布于 2015-12-15 14:30">淺談壓縮感知（十）：范數與稀疏性</a> posted @ 2015-12-15 09:29?AndyJee 閱讀(28779) 評論(1) 編輯 收藏 </div><!--end: topics 文章、評論容器-->
評論列表 </div>
#1樓
2017-06-20 12:59
<a id="a_comment_author_3718412" href="https://home.cnblogs.com/u/1178954/" target="_blank">數理小先生</a></div><div class="feedbackCon"> 樓主在嗎？急需要求助 支持(0) 反對(0) </div></div> 刷新評論刷新頁面返回頂部 注冊用戶登錄后才能發表評論，請 登錄 或 注冊， 訪問 網站首頁。 【推薦】超50萬C++/C#源碼: 大型實時仿真組態圖形源碼【推薦】零基礎輕松玩轉華為云產品，獲壕禮加返百元大禮
【推薦】新手上天翼云，數十款云產品、新一代主機0元體驗
【推薦】華為IoT平臺開發者套餐9.9元起，購買即送免費課程
相關博文：
· MP和OMP算法
· MP算法和OMP算法及其思想
· MP算法和OMP算法及其思想
· MP算法和OMP算法及其思想
· MP算法和OMP算法及其思想
<div id="google_ads_iframe_/1090369/C2_0__container__" style="border: 0pt none;"><iframe id="google_ads_iframe_/1090369/C2_0" title="3rd party ad content" name="google_ads_iframe_/1090369/C2_0" width="468" height="60" scrolling="no" marginwidth="0" marginheight="0" frameborder="0" srcdoc="" style="border: 0px; vertical-align: bottom;" data-google-container-id="2" data-load-complete="true"></iframe></div></div> </div> <div id="under_post_kb"> 最新 IT 新聞:
· Windows 10終于拿下操作系統市場半壁江山
· 藏在1.85億人體內的隱形致癌病毒，有人確診即是晚期
· 嘀嗒出行：用戶突破1.3億，車主突破1500萬
· 無人機泡沫破裂：創企紛紛關門 風投遭重創
· 基因編輯明星動物遇挫：無角牛體內發現細菌基因污染 | 專訪 ?
? 更多新聞...
總結
以上是生活随笔為你收集整理的正交匹配追踪算法OMP的全部內容，希望文章能夠幫你解決所遇到的問題。
- 上一篇： 邮箱服务器未运行,smtp服务器,smt
- 下一篇： VBA操作WORD（二）：替换字符（含空