机器学习——回归算法之正交匹配追踪算法

前言

学习正交匹配追踪算法之前需要看懂匹配追踪算法(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$
在每一次迭代时，都要计算：
$$f=f_k+R_{k}f$$
其中$f_k表示当前近似，R_{k}f表示当前的残差（下面的算法流程中有计算公式）$
算法具体过程：

1. 初始化参数：$k=1,f_0=0,R_{0}f=f=x_i,其中x_i表示随机在原经归一化后的数据集中选择一个样本$
2. 计算内积 { ⟨ R k f , x n ⟩ } n \{\left \langle R_{k}f,x_{n} \right \rangle\}_{n} {⟨Rk​f,xn​⟩}n​
3. 按照下式选择下一个样本(即从剩下的样本中选择大于上一次所选样本的所有内积的样本)：
4. $$\Vert ⟨R_{k}f,x_{n_{k+1}}⟩\Vert \ge \alpha \underset{j}{\sup }\Vert ⟨R_{k}f,x_j⟩\Vert ,0<\alpha \le 1$$
5. 更新参数：$$f_{k+1}=f_k+⟨R_{k}f,x_{n_{k+1}}⟩x_{n_{k+1}}$$
   $$R_{k+1}f=R_{k}f-⟨R_{k}f,x_{n_{k+1}}⟩x_{n_{k+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={\varphi }_t({\varphi }_t^T{\varphi }_t{)}^{-1}{\varphi }_t^T$$
   通过从e0减去其在Φt所张成空间上的正交投影得到残差e1:$$e_1=e_0-P{e}_0=(I-P)e_0$$
4. 对残差迭代执行(2)、(3)步$$e_{m+1}=e_m-P{e}_m=(I-P)e_m$$其中I为单位阵。需要注意的是在迭代过程中Φt为所有被选择过的原子组成的矩阵，因此每次都是不同的，所以由它生成的正交投影算子矩阵P每次都是不同的。
5. 直到达到某个指定的停止准则后停止算法

OMP减去的$P{e}_m$是$e_m$在所有被选择过的原子组成的矩阵${\varphi }_t$所张成空间上的正交投影，而MP减去的$P{e}_m$是$e_m$在本次被选择的原子${\varphi }_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