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算法部分。

四、代码实现

代码如下：

cv::Mat ormpSparseRepresentation::ompSparseL2(const cv::Mat Dict, const cv::Mat Y, const int K)
{
    int r = Dict.rows;
    int c = Dict.cols;
    int n = Y.cols;
    cv::Mat ERR(r,1,CV_32F);
    ERR = Y;
    int size[] = {c,n};
    cv::Mat A = cv::Mat(2,size,CV_32F,cv::Scalar(0.f));
    QVector<int> index;
    cv::Mat U = cv::Mat::ones(1,c,CV_32F);
    cv::Mat tmpA;
    for(int i = 0;i<K;i++)
    {
        cv::Mat S = ERR.t()*Dict;
        cv::pow(S,2,S);
        if(S.rows != 1)
            cv::reduce(S,S,0,CV_REDUCE_SUM);
        cv::sqrt(S,S);
        S = S/U;
        if(i!=0)
        {
            for(int j = 0;j<index.size();j++)
            {
                S.at<float>(0,index[j]) = 0.f;
            }
        }

        cv::Point maxLoc;
        cv::minMaxLoc(S,NULL,NULL,NULL,&maxLoc);
        int pos = maxLoc.x;
        index.append(pos);

        cv::Mat subDict;
        getColDictFormIndex(Dict,index,subDict);

        cv::Mat invSubDict;
        cv::invert(subDict,invSubDict,cv::DECOMP_SVD);

        tmpA = invSubDict*Y;
        ERR = Y - subDict*tmpA;

        cv::Mat Dict_T_Dict;
        cv::mulTransposed(subDict,Dict_T_Dict,1);
        cv::Mat invDict_T_Dict;
        cv::invert(Dict_T_Dict,invDict_T_Dict,cv::DECOMP_SVD);

        cv::Mat P = (cv::Mat::eye(r,r,CV_32F) - subDict*invDict_T_Dict*subDict.t())*Dict;
        cv::pow(P,2,P);
        cv::reduce(P,P,0,CV_REDUCE_SUM);
        cv::sqrt(P,U);
    }
    for(int i = 0;i<K;i++)
    {
        int tmpC = index[i];
        const float *inP=tmpA.ptr<float>(i);
        float *outP=A.ptr<float>(tmpC);
        for(int j = 0;j<n;j++)
        {
            outP[j] = inP[j];
        }
    }
    return A;
}

void ormpSparseRepresentation::getColDictFormIndex(const cv::Mat Dict, const QVector<int> index, cv::Mat &res)
{
    if(index.size() == 0)
        return;
    if(!Dict.data)
        return;

    int r = Dict.rows;
    int c = index.size();

    cv::Mat Dict_T;
    cv::transpose(Dict,Dict_T);

    cv::Mat res_T = cv::Mat(c,r,Dict.type());

    for(int i = 0;i<index.size();i++)
    {
        int tmpC = index[i];
        const float *inP=Dict_T.ptr<float>(tmpC);
        float *outP=res_T.ptr<float>(i);
        for(int j = 0;j<r;j++)
        {
            outP[j] = inP[j];
        }
    }
    cv::transpose(res_T,res);
    res_T.release();
    Dict_T.release();
}