精华内容
下载资源
问答
  • 压缩感知的稀疏重构中广泛应用的正交匹配追踪(OMP)算法matlab程序,该算法由香港大学电子工程系 沙威老师开发,代码注释详细,便于读者理解。已测试,可以正常运行。读者通过代码可以加深对该算法以及压缩感知、...
  • 正交匹配追踪法,可用于一维信号的重构,
  • matlab正交匹配追踪算法
  • 受压缩感知理论的启发, 基于参数向量所具有的稀疏特性, 提出一种新的阈值正交匹配追踪算法辨识系统的参数和时滞. 仿真实验表明, 所提出的算法能在少量采样时有效地辨识系统参数、估计未知时滞, 同时验证了算法的...
  • 现有的正交匹配追踪(OMP)算法都是在给定迭代次数(待重建图像的稀疏度)的条件下重建,这使其需要通过非常多的线性测量来保证精确重建.为此,文中提出一种改进的后退型最优OMP方法:首先利用最优正交匹配追踪(OOMP)算法在...
  • omp(正交匹配追踪算法、MATLAB编写)
  • 匹配追踪MP、正交匹配追踪算法OMP,稀疏表示里的基本算法
  • 2.MP算法(匹配追踪算法) 2.1 算法描述 作为对信号进行稀疏分解的方法之一,将信号在完备字典库上进行分解。 假定被表示的信号为y,其长度为n。假定H表示Hilbert空间,在这个空间H里,由一组向量 构成字典...

    主要介绍MP(Matching Pursuits)算法和OMP(Orthogonal Matching Pursuit)算法[1],这两个算法虽然在90年代初就提出来了,但作为经典的算法,国内文献(可能有我没有搜索到)都仅描述了算法步骤和简单的应用,并未对其进行详尽的分析,国外的文献还是分析的很透彻,所以我结合自己的理解,来分析一下写到博客里,算作笔记。

    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算法的代码,源出处不得而知,共享给大家。

    再贴另外一个洋牛paper[3]中关于OMP的描述,之所以引入,是因为它描述的非常严谨,但是也有点苦涩难懂,不过有了上面的基础,就容易多了。


    它的描述中的Sweep步骤就是寻找与当前残差最大的内积时列在字典矩阵D中的索引,它的这个步骤描述说明为什么要选择内积最大的以及如何选择。见下图,说的非常清晰。


    它的算法步骤Update Provisional Solution中求很简单,就是在 b = Ax 已知 A和b求x, 在x的最小二范就是A的伪逆与b相乘,即:


    看上去头疼,其实用matlab非常简单,看看上面的matlab的代码就明白了。

    我们可以看得出来,算法流程清晰明了,还是很好理解的。这正是OMP算法的魅力,作为工具使用简单,背后却隐藏着很有趣的思想。

    写这篇博客的目的,是因为搜索了一下,MP和OMP没有人很详细的介绍。文献[1]讲的很清楚的,大家有兴趣可以找来看看。不要被老板发现我居然在搜中文文献还写中文博客。



    参考文献:

    [1] Orthogonal Matching Pursuit:Recursive Function Approximat ion with Applications to Wavelet Decomposition
    [2]http://wenku.baidu.com/view/22f3171614791711cc7917e4.html

    [3] From Sparse Solutions of Systems of Equations to Sparse Modeling of Signals and Images

    展开全文
  • omp正交匹配追踪算法

    2013-11-19 16:54:22
    该资源详细描述了OMP代码的matlab程序和c语言程序(矩阵的求逆采用LU分解法),并且对两者结果进行了比较,恢复的信号可以精确到小数点5位,误差非常小,测量矩阵采用随机高斯矩阵,程序里面还有matlab和c语言版对...
  • 在分析分段正交匹配追踪StOMP算法迭代过程中多原子匹配方法的基础上,为了进一步减少算法迭代次数,提高重构精度,提出了基于互相关向量的自适应极限因子选取和迭代结束条件重设的优化方法。通过MATLAB编程实验,在...
  • 机器学习——回归算法之正交匹配追踪算法(Orthogonal Matching Pursuit,OMP)前言匹配追踪算法(Matching Pursuit,MP)1、算法思想2、算法过程3、MP算法的问题正交匹配追踪算法(Orthogonal Matching Pursuit,OMP)1、...

    前言

    学习正交匹配追踪算法之前需要看懂匹配追踪算法(MP),显然OMP是在MP的基础上添加一个正交约束,具体是怎样的请往下看。一定要明白稀疏分解和压缩感知的区别和联系

    匹配追踪算法(Matching Pursuit,MP)

    1、算法思想

    MP算法的基本思想:从原始数据集X(也称为过完备原子库中),选择一个与信号 y(标签) 最匹配的原子(也就是某些列),构建一个稀疏逼近,并求出信号残差,然后继续选择与信号残差最匹配的原子,反复迭代,信号y可以由这些原子的线性和,再加上最后的残差值来表示。很显然,如果残差值在可以忽略的范围内(迭代停止条件),则信号y就是这些原子的线性组合。

    2、算法过程

    假设原始数据集有n个样本 X = { x 1 , x 2 , . . , x n } X= \{x_{1},x_{2},..,x_{n}\} X={x1,x2,..,xn},并且每个样本属性已经归一化 ∣ ∣ x i ∣ ∣ = 1 ||x_{i}||=1 xi=1
    在每一次迭代时,都要计算:
    f = f k + R k f f=f_{k}+R_{k}f f=fk+Rkf
    其中 f k 表 示 当 前 近 似 , R k f 表 示 当 前 的 残 差 ( 下 面 的 算 法 流 程 中 有 计 算 公 式 ) f_{k}表示当前近似,R_{k}f表示当前的残差(下面的算法流程中有计算公式) fkRkf
    算法具体过程:

    1. 初始化参数: k = 1 , f 0 = 0 , R 0 f = f = x i , 其 中 x i 表 示 随 机 在 原 经 归 一 化 后 的 数 据 集 中 选 择 一 个 样 本 k =1,f_{0}=0,R_{0}f=f=x_{i},其中x_{i}表示随机在原经归一化后的数据集中选择一个样本 k=1,f0=0,R0f=f=xi,xi
    2. 计算内积 { ⟨ R k f , x n ⟩ } n \{\left \langle R_{k}f,x_{n} \right \rangle\}_{n} {Rkf,xn}n
    3. 按照下式选择下一个样本(即从剩下的样本中选择大于上一次所选样本的所有内积的样本):
    4. ∥ ⟨ R k f , x n k + 1 ⟩   ∥ ≥ α s u p j ∥ ⟨ R k f , x j ⟩   ∥ , 0 < α ≤ 1 \left\|\left \langle R_{k}f,x_{n_{k+1}} \right \rangle\ \right\|\geq \alpha \underset{j}{sup}\left\|\left \langle R_{k}f,x_{j} \right \rangle\ \right\|,0<\alpha \leq 1 Rkf,xnk+1 αjsupRkf,xj ,0<α1
    5. 更新参数: f k + 1 = f k + ⟨ R k f , x n k + 1 ⟩ x n k + 1 f_{k+1} =f_{k}+ \left \langle R_{k}f,x_{n_{k+1}} \right \rangle x_{n_{k+1}} fk+1=fk+Rkf,xnk+1xnk+1
      R k + 1 f = R k f − ⟨ R k f , x n k + 1 ⟩ x n k + 1 R_{k+1}f = R_{k}f-\left \langle R_{k}f,x_{n_{k+1}} \right \rangle x_{n_{k+1}} Rk+1f=RkfRkf,xnk+1xnk+1
    6. k=k+1,循环迭代2~5,直到残差满足一定阈值为止。

    3、MP算法的问题

    在描述MP算法时,有类似这样的话:在匹配追踪(MP)中,字典原子不是相互正交的向量。因此上面减去投影计算残差的过程中会再次引入与前面使用的原子不正交的成分。或者是:信号(残值)在已选择的原子进行垂直投影是非正交性的,这会使得每次迭代的结果并不是最优的而是次最优的,收敛需要很多次迭代。
    j解决办法:简要概括就是我们是在找字典中的一组基进行线性组合后来作为f的最接近表示,每次的残差也就是Rkf会和f在当前所选择的基xk上的正交投影垂直,上式中的 fk是多个选择后的基的线性组合,不和残差项垂直。也就是下面OMP算法使用的正交。

    正交匹配追踪算法(Orthogonal Matching Pursuit,OMP)

    1、算法思想

    在正交匹配追踪OMP中,残差是总与已经选择过的原子正交的。这意味着一个原子不会被选择两次,结果会在有限的几步收敛。

    2、算法流程

    1. 用x表示初始样本,初始化残差e0=x;
    2. 在剩余的样本中选择与e0内积绝对值最大的原子,表示为φ1;
    3. 将选择的原子作为列组成矩阵Φt,定义Φt列空间的正交投影算子为: P = ϕ t ( ϕ t T ϕ t ) − 1 ϕ t T P =\phi_{t}(\phi_{t}^{T}\phi_{t})^{-1}\phi_{t}^{T} P=ϕt(ϕtTϕt)1ϕtT
      通过从e0减去其在Φt所张成空间上的正交投影得到残差e1: e 1 = e 0 − P e 0 = ( I − P ) e 0 e_{1}=e_{0}-Pe_{0}=(I-P)e_{0} e1=e0Pe0=(IP)e0
    4. 对残差迭代执行(2)、(3)步 e m + 1 = e m − P e m = ( I − P ) e m e_{m+1}=e_{m}-Pe_{m}=(I-P)e_{m} em+1=emPem=(IP)em其中I为单位阵。需要注意的是在迭代过程中Φt为所有被选择过的原子组成的矩阵,因此每次都是不同的,所以由它生成的正交投影算子矩阵P每次都是不同的。
    5. 直到达到某个指定的停止准则后停止算法

    OMP减去的 P e m Pe_{m} Pem e m e_{m} em在所有被选择过的原子组成的矩阵 ϕ t \phi_{t} ϕt所张成空间上的正交投影,而MP减去的 P e m Pe_{m} Pem e m e_{m} em在本次被选择的原子 ϕ m \phi_{m} ϕm所张成空间上的正交投影。

    3、代码

    代码好像不对,望高人指点!!!!

    def OMP(X,y,K):#k表示稀疏度
        #初始化phi
        tmp = []
        #归一化
        _range = np.max(X) - np.min(X)
        X = (X - np.min(X)) / _range
        #初始化残差e0
        np.random.seed(0)
        p = np.random.permutation(range(len(X)))#一同打乱x和y
        X, Y = X[p], y[p]
        residual_e = Y.reshape(X.shape[0],-1)
        t = 0
        while (t<K or np.any(X)):
            t += 1
            inner_products = np.abs(np.dot(residual_e.T, X))
            max_index = inner_products.argmax()
            tmp.append(X[:, max_index])
            phi  = np.array(tmp).reshape(X.shape[0],t)
            X = np.delete(X, [max_index], axis=1)
            P = np.dot(np.dot(phi, np.linalg.inv(np.dot(phi.T, phi))), phi.T)
            residual_e = np.dot((np.eye(P.shape[0], P.shape[1])-P),residual_e)
        return residual_e
    
    展开全文
  • omp算法matlab代码稀疏重建算法 此Repo由适用于以下稀疏信号重构算法的MATLAB代码组成: 正交匹配追踪-OMP 展望OMP-LAOMP 批处理OMP-BLAOMP 协同软件 子空间追踪 流式细胞仪
  • 压缩感知的omp算法的源代码,绝对简单,初学者一看就懂,绝对不骗人
  • 正交匹配追踪

    千次阅读 2018-08-05 09:01:07
    OrthogonalMatchingPursuit (正交匹配追踪法)和 orthogonal_mp 使用了 OMP 算法近似拟合了一个带限制的线性模型,该限制影响于模型的非 0 系数(例:L0 范数)。 就像最小角回归一样,作为一个前向特征选择方法,...

    OrthogonalMatchingPursuit (正交匹配追踪法)和 orthogonal_mp
    使用了 OMP 算法近似拟合了一个带限制的线性模型,该限制影响于模型的非 0 系数(例:L0 范数)。

    就像最小角回归一样,作为一个前向特征选择方法,正交匹配追踪法可以近似一个固定非 0 元素的最优向量解:

    \text{arg,min,} ||y - X\gamma||_2^2 \text{ subject to } \
    ||\gamma||0 \leq n{nonzero_coefs}

    正交匹配追踪法也可以针对一个特殊的误差而不是一个特殊的非零系数的个数。可以表示为:

    \text{arg,min,} ||\gamma||_0 \text{ subject to } ||y-X\gamma||_2^2 \
    \leq \text{tol}

    OMP 是基于每一步的贪心算法,其每一步元素都是与当前残差高度相关的。它跟较为简单的匹配追踪(MP)很相似,但是相比 MP 更好,在每一次迭代中,可以利用正交投影到之前选择的字典元素重新计算残差。

    示例:

    Orthogonal Matching Pursuit

    Using orthogonal matching pursuit for recovering a sparse signal from a noisy measurement encoded with a dictionary

    A

    print(__doc__)
    
    import matplotlib.pyplot as plt
    import numpy as np
    from sklearn.linear_model import OrthogonalMatchingPursuit
    from sklearn.linear_model import OrthogonalMatchingPursuitCV
    from sklearn.datasets import make_sparse_coded_signal
    
    n_components, n_features = 512, 100
    n_nonzero_coefs = 17
    
    # generate the data
    ###################
    
    # y = Xw
    # |x|_0 = n_nonzero_coefs
    
    y, X, w = make_sparse_coded_signal(n_samples=1,
                                       n_components=n_components,
                                       n_features=n_features,
                                       n_nonzero_coefs=n_nonzero_coefs,
                                       random_state=0)
    
    idx, = w.nonzero()
    
    # distort the clean signal
    y_noisy = y + 0.05 * np.random.randn(len(y))
    
    # plot the sparse signal
    plt.figure(figsize=(7, 7))
    plt.subplot(4, 1, 1)
    plt.xlim(0, 512)
    plt.title("Sparse signal")
    plt.stem(idx, w[idx])
    
    # plot the noise-free reconstruction
    展开全文
  • 稀疏正交匹配追踪

    2018-05-02 10:48:04
    稀疏正交匹配追踪(OMP算法)正交匹配追踪(OMP)算法属于贪婪算法。而贪婪算法是一种不追求最优解,只希望得到较为满意解的方法。贪婪法一般可以快速得到满意的解,因为它省去了为找最优解要穷尽所有可能而必须耗费...
  • 匹配追踪MP和正交匹配追踪OMP算法

    万次阅读 多人点赞 2018-03-01 15:25:38
    匹配追踪MP和正交匹配追踪OMP算法...

    匹配追踪MP和正交匹配追踪OMP算法

    http://blog.csdn.net/wwf_lightning/article/details/70142985

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

    http://blog.csdn.net/jbb0523/article/details/45130793

    MP:matching pursuit匹配追踪
    OMP:正交匹配追踪
    主要介绍MP与OMP算法的思想与流程,解释为什么需要引入正交?


        !!今天发现一个重大问题,是在读了博主的正交匹配追踪(OMP)在稀疏分解与压缩感知重构中的异同http://blog.csdn.net/jbb0523/article/details/45100659之后一脸懵逼,CS中的稀疏表示不就是把信号转换到另一个变换域中吗?怎么跑出来一个稀疏分解里面又有MP和OMP算法!!后来发现原来稀疏分解先于压缩感知提出,信号稀疏表示的目的就是在给定的超完备字典中用尽可能少的原子来表示信号,可以获得信号更为简洁的表示方式,从而使我们更容易地获取信号中所蕴含的信息,更方便进一步对信号进行加工处理,如压缩、编码等。后面的学者用稀疏分解的思想应用于压缩感知重构中。其实两者解决的问题是一样的。

        从数学模型来入手分析这个问题:

         1)稀疏分解要解决的问题是在冗余字典A中选出k列,用这k列的线性组合近似表达待稀疏分解信号y,可以用表示为y=Aθ,求θ。
             
        2)压缩感知重构要解决的问题是事先存在一个θ和矩阵A,然后得到y=Aθ(压缩观测),现在是在已知y和A的情况下要重构θ。

         看到了没?实际上它们要解决的问题都是对已知y和A的情况下求y=Aθ中的θ。
         上面各式中,AM×N矩阵(M>>N,稀疏分解中为冗余字典,压缩感知中为传感矩阵A=ΦΨ,即测量矩阵Φ乘以稀疏矩阵Ψ),yM×1的列向量(稀疏分解中为待稀疏分解信号,压缩感知中为观测向量),θN×1的列向量(稀疏分解中为待求分解系数,压缩感知中为信号x在变换域Ψ的系数x=Ψθ)。
         所不同的是,在稀疏分解中θ是事先不存在的,我们要去求一个θ近似表示y,求出的θ并不能说对与错;在压缩感知中,θ是事先存在的,只是现在不知道,我们要通过某种方法如OMP去把θ求出来,求出的θ应该等于原先的θ的,然后可求原信号x=Ψθ
        其实MP也好,改进后的OMP也罢,最初提出都是面向稀疏分解的,当时还没有压缩感知的概念,只是后来压缩感知提出后将其引入到了压缩感知重构中,因为前面也说了,其实他们的本质是一样的,都是已知y和A的情况下求y=Aθ中的θ。


    1.冗余字典与稀疏表示   
         作为对信号进行稀疏分解的方法之一,将信号在完备字典库上进行分解。即在字典中找到一组基来表示信号,而用一组特定基表达一个信号其实就是找到相应的一组展开系数一组基表达信号的能力取决于信号的特性是否与基向量的特性相吻合。例如,光滑连续信号可以被傅里叶基稀疏的表达,但脉冲信号不行。再如,带有孤立不连续点的平滑信号可用小波基稀疏表达,但小波基在表达傅里叶频谱中有窄带高频支撑的信号时却是无效的。
         现实世界中的信号经常包含有用单一基所不能表达的特征。对于这些信号,你或许希望可以选择来自不同基的向量(如用小波基和傅里叶基来联合表达一个信号)。因为你想保证你可以表达一个信号空间的所有信号向量,所以由所有可选向量组成的字典应该能够张成这个信号空间。然而由于这组字典中的向量来自不同的基,它们可能不是线性独立的,会造成用这组字典做信号表达时系数不唯一。然而如果创建一组冗余字典,你就可以把你的信号展开在一组可以适应各种时频或时间-尺度特性的向量上。你可以自由的创建包含多个基的字典。例如,你可以构造一组表达平方可积空间的基,这组基包含小波包基和局部余弦基。这样构造的字典可以极大地增加你稀疏表达各种特性信号的能力。

    2.字典非线性近似
         定义表达你的信号空间的归一化基本模块作为字典。这些归一化向量叫做原子。如果字典的原子张成了整个信号空间,那么字典就是完全的。如果有原子之间线性相关,那么字典就是冗余的在大多数匹配追踪的应用中,字典都是完全冗余的。

         假设{φk}表示字典的原子。假设字典是完全且冗余的。使用字典的线性组合表达信号将是不唯一的。


        一个重要的问题是是否存在一种最好的表达方式,一种直观的最好方式是选择φk使得近似信号和原始信号有最大的内积,如最好的φk满足


         即对于正交原子,为投影到由φk张成的子空间上的幅值。

         匹配追踪的中心问题是你如何选择信号在字典中最优的M个展开项。

    3.MP算法

    》基本思想
         MP算法的基本思想:从字典矩阵D(也称为过完备原子库中),选择一个与信号 y 最匹配的原子(也就是某列),构建一个稀疏逼近,并求出信号残差,然后继续选择与信号残差最匹配的原子,反复迭代,信号y可以由这些原子的线性和,再加上最后的残差值来表示。很显然,如果残差值在可以忽略的范围内,则信号y就是这些原子的线性组合。

    》算法流程
        
         用Φ={φk}表示一个原子归一化的字典,x表示信号。
        (1) 首先初始化残差 e 0 = x
        (2)匹配追踪的第一步是从字典中找到与e0内积绝对值最大的原子,表示为φ1
        (3)通过从e0减去其在φ1所张成空间上的正交投影得到残差e1

        其中<e0, φ1>表示e0φ1的内积。由于已经说明Φ为原子归一化的字典,即φ1为单位列向量,所以e0φ1所张成空间上的正交投影可以表示为<e0, φ1>φ1 (由于为一个向量,所以e0φ1所张成空间上的正交投影即为e0φ1方向上的正交投影分量),若φ1不是单位列向量,则e0φ1方向上的正交投影分量为<e0, φ1>φ1/(φ1Tφ1),其中上标T表示转置。

         (4) 对残差迭代执行 (2) (3) 步;

         其中φm+1是从字典中找到与em的内积绝对值最大的原子。
         (5)直到达到某个指定的停止准则后停止算法。
        这里要从数学上说明一点:由于内积 < e m ,   φ m+1 > 实际上为一个数 ( 标量 ) ,所以

    分母中乘上φTm+1φm+1 ,一开始没想明白,注意这里的φm+1  是单位向量,则转置乘上自己是一个数值,也就是1
        若令 P = φ m+1  ( φ T m+1 φ m+1 ) -1 φ T m+1   ,(这里的P是用每次选择出的与残差内积最大的原子来计算的,而下面OMP中的P是用每次选择过的原子组成的矩阵来计算的!!!)则第 (4) 步的迭代公式为

         注意矩阵P每次迭代都是不同的。
    》提出一个问题        
        在描述MP算法时,有类似这样的话: 在匹配追踪 (MP) 中,字典原子不是相互正交的向量。因此上面减去投影计算残差的过程中会再次引入与前面使用的原子不正交的成分。或者是:信号(残值)在已选择的原子进行垂直投影是非正交性的,这会使得每次迭代的结果并不是最优的而是次最优的,收敛需要很多次迭代。
         那么如何理解什么是最优的,什么是次最优的?为什么不是正交的呢?
      
        首先回顾下正交投影, 一个向量 (b) 在另一个向量 (a) 上的投影:


        p称为b在a上的正交投影,正交投影就是法线垂直于a的投影(不知道这种说法对不对,有待考究) 实际上就是寻找在 a 上离 b 最近的点。如果我们把 p 看作是 a 的估计值,那么我们定义 e = b - p ,称 e 为误差 (error)
        现在,我们的任务是找到这样的p,我们可以规定p = xax是某个数),那么e = b - xa,因为ep也就是a垂直,所以有aT(b - xa) = 0,展开化简得到:
           x代入到p中,得到:我们发现,如果改变b,那么p相对应改变,然而改变ap无变化。
       
        有了上面的背景知识,我们可以正式进入主题了,投影矩阵(projection matrix)  p = Pb 显然这里有:
        这个投影到底有什么用呢,从线性代数的角度来说,Ax = b并不一定总有解,这在实际情况中会经常遇到(m > n)。所以我们就把b投影到向量p上,求解Ax = p
        接下来,我们可以考虑更高维度的投影,三维空间的投影是怎么样的呢,我们可以想象一个三维空间内的向量在该空间内的一个平面上的投影:

       我们假设这个平面的基 (basis)a1, a2.那么矩阵A = [a1 a2]的列空间就是该平面。假设一个不在该平面上的向量b在该平面上的投影是p,那么我们就有下面这个表达式p = x1a1+ x2a2= Ax我们的任务就是找到这样的x。这里有一个关键的地方:e = b - Ax与该平面垂直,所以a1T(b - Ax) = 0a2T(b- Ax) = 0.用矩阵的形式表达就是:AT(b - Ax) = 0.

          文献[3]中与文献[1]中所表示的算法流程符号略有不同,但意思是一样的,由于要借鉴这篇文章来解释上述提出的问题,所以先简要介绍一下该文献中所说明的MP的算法流程。
         假定被表示的信号为y,其长度为n。假定H表示Hilbert空间(希尔伯特空间即完备的内积空间),在这个空间H里,由一组向量{x1,x2,...xn}成字典矩阵D,其中每个向量可以称为原子(atom),其长度与被表示信号 y 的长度n相同,而且这些向量已归一化处理,即 ‖xi‖=1,也就是单位向量长度为1。
         (1)计算信号y与字典矩阵中每列(原子)的内积,选择绝对值最大的一个原子,它就是与信号y在本次迭代运算中最匹配的。用专业术语来描述:令信号y∈H,从字典矩阵中选择一个最为匹配的原子,满足,(sup为上确界)r0表示一个字典矩阵的列索引。这样,信号y就被分解为在最匹配原子x的垂直投影分量和残值两部分,即:
         之前想不明白为什么还要乘上xr0,联系一下傅里叶系数和傅里叶级数,此处y和xr0的内积是系数,而y是要由一组基的线性组合来表示的。
         (2)对残值R1f进行步骤[1]同样的分解,那么第K步可以得到:   其中 满足
    在别的参考文献中有形如这样的形式 之前没有搞明白,其实也 就是移项了而已 可见,经过 K 步分解后,信号 y 被分解为 ,  其中
         我们要对信号进行稀疏表示,即在字典中找到一组最适合描述信号的基,把信号表示成这组基的线性组合。那为什么会造成不正交呢?举个例子说明一下:在二维空间上,有一个信号yD=[x1, x2]来表达,MP算法迭代会发现总是在x1x2上反复迭代,即,这个就是信号(残值)在已选择的原子进行垂直投影的非正交性导致的。再用严谨的方式描述可能容易理解在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算法。
         简要概括就是我们是在找字典中的一组基进行线性组合后来作为f的最接近表示,每次的残差也就是Rkf会和f在当前所选择的基xk上的正交投影垂直,上式中的 fk是多个选择后的基的线性组合,不和残差项垂直。

    4.OMP算法

    》算法流程

         在正交匹配追踪OMP中,残差是总与已经选择过的原子正交的。这意味着一个原子不会被选择两次,结果会在有限的几步收敛。OMP的算法如下
        
        (1) x表示你的信号,初始化残差e0=x
        (2)选择与e0内积绝对值最大的原子,表示为φ1
        (3)将选择的原子作为列组成矩阵Φt,定义Φt列空间的正交投影算子为
         通过从e0减去其在Φt所张成空间上的正交投影得到残差e1
        (4)对残差迭代执行(2)(3)步;

         其中I为单位阵。需要注意的是在迭代过程中Φt为所有被选择过的原子组成的矩阵,因此每次都是不同的,所以由它生成的正交投影算子矩阵P每次都是不同的。
        (5)直到达到某个指定的停止准则后停止算法。

        OMP减去的 Pemem所有被选择过的原子组成的矩阵Φt所张成空间上的正交投影,而MP减去的Pemem本次被选择的原子φm所张成空间上的正交投影
    》提出一个问题  

        OMP是怎么实现与所有选择过的原子正交的?
         →施密特正交化
         在现代数学引论中有学习过,但是和线性代数中的表达式不太一样,对两者进行了比较,发现其实本质是一样的。所选择的一组基是线性无关的,我们可以通过施密特正交化来将这组选择的基转换为正交基。
     


        那么具体在OMP算法中是如何体现的?

        文献[4]中给出了施密特(Schimidt)正交化的过程:

        上面的的 [x,y]表示向量内积,[x,y]=xTy=yTx=[x,y]。施密特正交化公式中的br实际上可写为:
        分子之所以可以这么变化是由于[x,y]实际上为一个数,因此[x,y]x=x[x,y]= xxTy
        首先给出一个结论:

        设OMP共从冗余字典中选择了r个原子,分别是a1a2……ar,根据正交匹配追踪的流程可以知道待分解信号x最后剩余的残差eromp
             (1)
        该残差也可以表示为
             (2)

        其中矩阵A为选择的r个原子组成的矩阵,e(r-1)omp为选择(r-1)个原子时的残差。
        将选择的r个原子a1a2……ar进行施密特正交化:
         则残差eromp还可以写为
             (3)
         (1)一般出现在稀疏分解算法中,(2)一般出现在重构算法中,(3)是自己琢磨出来的(受到沙威的文档中提到的施密特正交化的启发,但沙威只限于向量情况下,详情可参见[6],此处相当于一个推广)
        OMP分解过程,实际上是将所选原子依次进行Schimidt正交化,然后将待分解信号减去在正交化后的原子上各自的分量即可得残差。其实(式3)求残差的过程也是在进行施密特正交化

        有个关键问题还是要说的,分解后在所选择各原子上的系数是多少呢?答案其实也很简单,各个系数是(ATA)-1ATx,即最小二乘解,这个解是一个列向量,每一个元素分别是组成矩阵A的各原子的线性组合系数,这个在下一篇《正交匹配追踪(OMP)在稀疏分解与压缩感知重构中的异同》也会明确再次说明。

        同理,若设MP共从冗余字典中选择了r个原子,分别是a1a2,……,ar,根据匹配追踪的流程可以知道待分解信号x每次迭代后剩余的残差ermp

        
        ·比较式(3)的第2个等号表示的 e romp与此处的ermp也可以体会出OMP与MP的区别吧。


         
    参考文献:
    [4] 同济大学数学系. 线性代数(第五版)[M].高等教育出版社,2007:114.
    [6] 沙威. “压缩传感”引论.http://www.eee.hku.hk/~wsha/Freecode/Files/Compressive_Sensing.pdf


    展开全文

空空如也

空空如也

1 2 3 4 5 ... 20
收藏数 1,870
精华内容 748
关键字:

匹配追踪与正交匹配追踪