7. 稀疏表示之OMP，SOMP算法及openCV實現

一、前言

稀疏表示是自上世紀90年代開始，從人眼的視覺感受野獲得啟示，逐漸被人們所研究。現在已經發展為一種重要的信息表示方法。所謂稀疏表示是指，一個信號在過完備字典中，可以由少數個原子線性表達，

其數學模型可以表達如下：

這個數學模型解算是一個NP-hard問題，也就是說只能通過窮舉去獲得最優解，其時間復雜度很大，幾乎無法獲得其精確的解算。在這種情況下，我們常用貪婪算法去獲得該模型的次最優解。本文介紹一種主流的貪婪算法——

正交匹配追蹤（OMP）。

二、OMP算法

貪婪算法的核心是每次從字典的原子中選擇一個最優原子來表示原始的信號。貪婪算法最大的缺點是，在貪婪算法的思想里，認為全局最優是每個局部最優得到的，這實際上很容易進入局部最優解，無法得到數據的全局最優解。

OMP作為貪婪算法中比較具有代表性的算法，其主要思想在于以下兩點：

1 認為字典原子在信號投影中的越大，對信號的描述越好；

2 每一次選擇的原子都與之前的原子正交。

介于以上兩點，OMP算法的描述如下：

上圖是從網上摘抄下來的。大概就是那樣。但是值得注意的是：這樣的OMP算法在解算的時候其效果往往沒有ORMP算法好，目前好多人說的OMP算法其實質往往是ORMP算法。

比如：開源的工具箱SPAMS上的OMP算法解算部分，其核心就是ORMP算法。ORMP算法相比如OMP算法的不同在于，在計算完殘差后對字典原子進行了另一個的拉伸（具體拉伸見后面代碼部分），如下圖：

三、SOMP算法

SOMP算法又叫同步OMP算法，主要思想為：相似的原子具有相同的稀疏特性。因此在對相似原子進行稀疏表示時，假設稀疏原子位于相同的位置，及其在過完備字典的選擇的原子相同，OMP算法是SOMP算法在原始信號為一個原子

時的特殊情況。OMP算法可以統一到SOMP算法當中，其解算流程幾乎同OMP算法部分。

四、代碼實現

代碼如下：

1 cv::Mat ormpSparseRepresentation::ompSparseL2(const cv::Mat Dict, const cv::Mat Y, const int K) 2 { 3 int r = Dict.rows; 4 int c = Dict.cols; 5 int n = Y.cols; 6 cv::Mat ERR(r,1,CV_32F); 7 ERR = Y; 8 int size[] = {c,n}; 9 cv::Mat A = cv::Mat(2,size,CV_32F,cv::Scalar(0.f)); 10 QVector<int> index; 11 cv::Mat U = cv::Mat::ones(1,c,CV_32F); 12 cv::Mat tmpA; 13 for(int i = 0;i<K;i++) 14 { 15 cv::Mat S = ERR.t()*Dict; 16 cv::pow(S,2,S); 17 if(S.rows != 1) 18 cv::reduce(S,S,0,CV_REDUCE_SUM); 19 cv::sqrt(S,S); 20 S = S/U; 21 if(i!=0) 22 { 23 for(int j = 0;j<index.size();j++) 24 { 25 S.at<float>(0,index[j]) = 0.f; 26 } 27 } 28 29 cv::Point maxLoc; 30 cv::minMaxLoc(S,NULL,NULL,NULL,&maxLoc); 31 int pos = maxLoc.x; 32 index.append(pos); 33 34 cv::Mat subDict; 35 getColDictFormIndex(Dict,index,subDict); 36 37 cv::Mat invSubDict; 38 cv::invert(subDict,invSubDict,cv::DECOMP_SVD); 39 40 tmpA = invSubDict*Y; 41 ERR = Y - subDict*tmpA; 42 43 cv::Mat Dict_T_Dict; 44 cv::mulTransposed(subDict,Dict_T_Dict,1); 45 cv::Mat invDict_T_Dict; 46 cv::invert(Dict_T_Dict,invDict_T_Dict,cv::DECOMP_SVD); 47 48 cv::Mat P = (cv::Mat::eye(r,r,CV_32F) - subDict*invDict_T_Dict*subDict.t())*Dict; 49 cv::pow(P,2,P); 50 cv::reduce(P,P,0,CV_REDUCE_SUM); 51 cv::sqrt(P,U); 52 } 53 for(int i = 0;i<K;i++) 54 { 55 int tmpC = index[i]; 56 const float *inP=tmpA.ptr<float>(i); 57 float *outP=A.ptr<float>(tmpC); 58 for(int j = 0;j<n;j++) 59 { 60 outP[j] = inP[j]; 61 } 62 } 63 return A; 64 } View Code 1 void ormpSparseRepresentation::getColDictFormIndex(const cv::Mat Dict, const QVector<int> index, cv::Mat &res) 2 { 3 if(index.size() == 0) 4 return; 5 if(!Dict.data) 6 return; 7 8 int r = Dict.rows; 9 int c = index.size(); 10 11 cv::Mat Dict_T; 12 cv::transpose(Dict,Dict_T); 13 14 cv::Mat res_T = cv::Mat(c,r,Dict.type()); 15 16 for(int i = 0;i<index.size();i++) 17 { 18 int tmpC = index[i]; 19 const float *inP=Dict_T.ptr<float>(tmpC); 20 float *outP=res_T.ptr<float>(i); 21 for(int j = 0;j<r;j++) 22 { 23 outP[j] = inP[j]; 24 } 25 } 26 cv::transpose(res_T,res); 27 res_T.release(); 28 Dict_T.release(); 29 } View Code

?

posted on 2015-08-28 21:35?機器學習豬 閱讀(...) 評論(...) 編輯 收藏

轉載于:https://www.cnblogs.com/zyore2013/p/4767896.html

總結

以上是生活随笔為你收集整理的7. 稀疏表示之OMP，SOMP算法及openCV实现的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。