压缩感知初入门小白，如有错误，欢迎指正交流~

原理

在上篇博文中，指出当前压缩感知中常用的重构算法，其中贪婪算法中有使用广泛的两种算法，即匹配追踪算法（Matching Pursuit，MP）和正交匹配追踪算法（Orthogonal Matching Pursuit，OMP）。本文对以上两种算法进行总结，用于存档~

提出

匹配追踪算法（MP）起初提出是基于稀疏分解，在压缩感知概念还未被提出之前，稀疏分解要解决的问题是在冗余字典A中选出k列，用这k列的线性组合近似表达待稀疏分解信号y，可以表示为y=Aθ，求θ；
而压缩感知所要解决的问题是事先存在一个θ和矩阵A，然后得到y=Aθ（压缩观测），即在已知y和A的情况下要重构θ。二者在解决的问题上基本一致，都是已知y和A求θ。

即：已知传感矩阵A={ai}(i=1,2…N)，其中ai代表的是一个列向量，即每一列称作一个原子。我们认为观测向量y是由A中的原子构成，而θ即可看作是每个原子对y的贡献，即系数。而已知y和A求θ，就是求每个原子对y的贡献程度。

核心思想：二者的核心思想基本一致，即遍历字典A中的每一个原子ai，根据内积最大化的原则找到与贡献最大的原子，作为当前的匹配原子。之后从信号中减去与该原子相近的成分，对剩余部分（即残差）与其他原子求贡献，找到贡献次大的原子，依次迭代，直至满足迭代停止条件（最大次数/迭代阈值等），停止迭代。

区别：二者的区别在于同一个原子的成分是否会被重复选择。OMP保证了残差与所有选择过的原子正交，即同一个原子不会被二次选择。而MP中不能保证正交，所以会有已选择过的原子的相近成分被二次选择，增加了收敛速度。

步骤

匹配追踪（MP）

Input：传感矩阵A（MxN），观测向量y（Mx1）
Output：稀疏系数向量s（Mx1）
Step1：初始化：残差r=y，系数x[N]={0}；
Step2：计算A中每个原子ai与残差r的内积<ai，r>。求出内积最大的原子ai，并记录这个原子对应的索引 i；
Step3：计算θ：系数向量θ的第i个元素θi=<ai，r>。（可以理解为残差在di这个方向上的投影）；
Step4：更新残差r：更新后残差r(i+1)=当前残差ri-di*θi;
Step5：判断迭代终止条件你，不满足则转至Step2。

正交匹配追踪（OMP）

Input：传感矩阵A（MxN），观测向量y（Mx1）
Output：稀疏系数向量s（Mx1）
Step1：初始化：残差r=y，系数x[N]={0}，被选择原子所构成的字典子集S=Ø（每次迭代被选中的原子添加进这个子集作为新的一列）；
Step2：计算A中每个原子ai与残差r的内积<ai，r>。求出内积最大的原子ai，并记录这个原子对应的索引 i。将ai添加到字典子集S中（每迭代一次，S中就会增加一个新的原子，这个加入的新原子就是此次迭代中发现的与残差内积最大的那个原子）；
Step3：计算θ：即y=Sθ的最小二乘解（最小二乘解保证了残差与选择过的原子正交）。计算得到的最小二乘解维度和S一样，是小于稀疏系数θ的维度的，这个只需要把最小二乘解中的各个元素赋值到θ的对应元素就行了（第一步中记录了被选中原子的索引，就是为了实现这个功能）;

注：区别：MP稀疏系数的计算是每次迭代只计算系数向量θ中的一个元素θi；而OMP中每次都要根据新加入的原子所更新的子集字典S重新计算整个系数向量θ(θ为向量)
Step4：更新残差r：更新后残差r(i+1)=当前残差ri-Sθ;
Step5：判断迭代终止条件你，不满足则转至Step2。

代码实现

MP较为简单，这里只放上了OMP的C语言实现代码：

#include<stdio.h>
#include<stdlib.h>
#include<malloc.h>
#include <windows.h>
#include <math.h>
#include <assert.h>//定义矩阵数据类型
typedef struct
{double **mat;int m, n;
}matrix;matrix OMP(matrix y, matrix A, int t);//为矩阵申请存储空间
void initial_mat(matrix *T, int m, int n)
{int i;(*T).mat = (double**)malloc(m * sizeof(double*));for (i = 0; i < m; i++){(*T).mat[i] = (double*)malloc(n * sizeof(double));}(*T).m = m;(*T).n = n;
}
//初始化矩阵
void initzero(matrix *T, int m, int n)
{int i, j;initial_mat(T, m, n);for (i = 0; i < m; i++){for (j = 0; j < n; j++){(*T).mat[i][j] = 0;}}
}
//释放存储空间
void destroy(matrix *T)
{int i;for (i = 0; i < (*T).m; i++){free((*T).mat[i]);}free((*T).mat);
}
//变换为单位矩阵
void set_identity_matrix(matrix m) {int i;int j;assert(m.m == m.n);for (i = 0; i < m.m; ++i) {for (j = 0; j < m.n; ++j) {if (i == j) {m.mat[i][j] = 1.0;}else {m.mat[i][j] = 0.0;}}}
}
//矩阵转置
void transpose_matrix(matrix input, matrix output)
{int i, j;assert(input.m == output.n);assert(input.n == output.m);for (i = 0; i < input.m; i++){for (j = 0; j < input.n; j++){output.mat[j][i] = input.mat[i][j];}}
}
//矩阵相乘
void multiply_matrix(matrix a, matrix b, matrix output)
{int i, j, k;assert(a.n == b.m);assert(output.m == a.m);assert(output.n == b.n);//printf("\n");for (i = 0; i < output.m; i++){for (j = 0; j < output.n; j++){output.mat[i][j] = 0.0;for (k = 0; k < a.n; k++){//printf("a%lf b%lf", a.mat[i][k], b.mat[k][j]);output.mat[i][j] += a.mat[i][k] * b.mat[k][j];}//printf("%lf ", output.mat[i][j]);}//printf("\n");}
}
/* 交换矩阵的两行 */
void swap_rows(matrix m, int r1, int r2) {double *tmp;assert(r1 != r2);tmp = m.mat[r1];m.mat[r1] = m.mat[r2];m.mat[r2] = tmp;
}
/*矩阵某行乘以一个系数  */
void scale_row(matrix m, int r, double scalar) {int i;assert(scalar != 0.0);for (i = 0; i < m.n; ++i) {m.mat[r][i] *= scalar;}
}/* Add scalar * row r2 to row r1. */
void shear_row(matrix m, int r1, int r2, double scalar) {int i;assert(r1 != r2);for (i = 0; i < m.n; ++i) {m.mat[r1][i] += scalar * m.mat[r2][i];}
}//矩阵求逆
int matrix_inversion(matrix input, matrix output)
{int i, j, r;double scalar, shear_needed;assert(input.m == input.n);assert(input.m == output.m);assert(input.m == output.n);set_identity_matrix(output);/* Convert input to the identity matrix via elementary row operations.The ith pass through this loop turns the element at i,i to a 1and turns all other elements in column i to a 0. */for (i = 0; i < input.m; ++i) {if (input.mat[i][i] == 0.0) {/* We must swap m to get a nonzero diagonal element. */for (r = i + 1; r < input.m; ++r) {if (input.mat[r][i] != 0.0) {break;}}if (r == input.m) {/* Every remaining element in this column is zero, so thismatrix cannot be inverted. */return 0;}swap_rows(input, i, r);swap_rows(output, i, r);}/* Scale this row to ensure a 1 along the diagonal.We might need to worry about overflow from a huge scalar here. */scalar = 1.0 / input.mat[i][i];scale_row(input, i, scalar);scale_row(output, i, scalar);/* Zero out the other elements in this column. */for (j = 0; j < input.m; ++j) {if (i == j) {continue;}shear_needed = -input.mat[j][i];shear_row(input, j, i, shear_needed);shear_row(output, j, i, shear_needed);}}return 1;
}
matrix OMP(matrix y, matrix A, int t)
{int M = A.m = y.m;int N = A.n;matrix s;initzero(&s, N, 1);matrix At;initzero(&At, M, t);matrix Pos_s;initzero(&Pos_s, 1, t);matrix r_n;initzero(&r_n, M, 1);//printf("\nr_n列向量：\n");for (int i = 0; i < M; i++){r_n.mat[i][0] = y.mat[i][0];//printf("%lf ", r_n.mat[i][0]);}matrix s_ls;initzero(&s_ls, t, 1);for (int d = 0; d < t; d++){matrix A_T;initzero(&A_T, N, M);transpose_matrix(A, A_T);matrix product;initzero(&product, N, 1);multiply_matrix(A_T, r_n, product);/*printf("\n product列向量：\n");for (int i = 0; i < N; i++){printf("%lf ", product.mat[i][0]);}*/int pos = 0;double max = fabs(product.mat[0][0]);for (int i = 1; i < N; i++){if (max < fabs(product.mat[i][0])){max = fabs(product.mat[i][0]);pos = i;}}//printf("\n pos：%d\n",pos);matrix Atd;initzero(&Atd, M, d+1);for (int i = 0; i < M; i++){Atd.mat[i][d] = A.mat[i][pos];}Pos_s.mat[0][d] = pos;for (int i = 0; i < M; i++){A.mat[i][pos] = 0;}matrix Atd_T;initzero(&Atd_T, d+1, M);transpose_matrix(Atd, Atd_T);matrix temp1;initzero(&temp1, d+1, d+1);multiply_matrix(Atd_T, Atd, temp1);/*printf("\n乘积：\n");for (int i = 0; i < d+1; i++){for (int j = 0; j < d+1; j++){printf("%lf ", temp1.mat[i][j]);}}*/matrix temp2;initzero(&temp2, d+1, d+1);matrix_inversion(temp1, temp2);/*printf("\n求逆：\n");for (int i = 0; i < d+1; i++){for (int j = 0; j < d+1; j++){printf("%lf ", temp2.mat[i][j]);}}*/matrix temp3;initzero(&temp3, d+1, M);multiply_matrix(temp2, Atd_T, temp3);/*printf("\n乘ATD_T：\n");for (int i = 0; i < d + 1; i++){for (int j = 0; j < M; j++){printf("%lf ", temp3.mat[i][j]);}}*/matrix s_ls_d;initzero(&s_ls_d, d + 1, 1);multiply_matrix(temp3, y, s_ls_d);/*printf("\ns：\n");for (int i = 0; i < d + 1; i++){for (int j = 0; j < 1; j++){printf("%lf ", s_ls_d.mat[i][j]);}}*/for (int i = 0; i < d + 1; i++){s_ls.mat[i][0] = s_ls_d.mat[i][0];}matrix temp4;initzero(&temp4, M, 1);multiply_matrix(Atd, s_ls_d, temp4);for (int i = 0; i < M; i++){r_n.mat[i][0] = y.mat[i][0] - temp4.mat[i][0];}}/*printf("\ns_ls:\n");for (int i = 0; i < t; i++){printf("%lf ", s_ls.mat[i][0]);}*/for (int i = 0; i < t; i++){int index = Pos_s.mat[0][i];//printf("[%d]%lf ", index, Pos_s.mat[0][i]);//printf("\n");s.mat[index][0] = s_ls.mat[i][0];//printf("[%d]%lf ",index, s_ls.mat[i][0]);}return s;
}
void main()
{matrix A;initzero(&A, 2, 5);A.mat[0][0] = 0.0591; A.mat[0][1] = -1.6258; A.mat[0][2] = 2.6052; A.mat[0][3] = 0.2570; A.mat[0][4] = -1.1464; A.mat[1][0] = -1.4669; A.mat[1][1] = -1.9648; A.mat[1][2] = 0.9724; A.mat[1][3] = -0.9742; A.mat[1][4] = 0.5476; printf("传感矩阵A:\n");for (int i = 0; i < 2; i++){for (int j = 0; j < 5; j++){printf("%lf ", A.mat[i][j]);}printf("\n");}matrix y;initzero(&y, 2, 1);y.mat[0][0] = 7.9498;y.mat[1][0] = 2.9672;printf("\n观测值y:\n");for (int i = 0; i < 2; i++){for (int j = 0; j < 1; j++){printf("%lf ", y.mat[i][j]);}printf("\n");}int t = 1;matrix s;initzero(&s, 5, 1);s = OMP(y, A, t);matrix PSi;initzero(&PSi, 5, 5);for (int i = 0; i < 5; i++){PSi.mat[i][i] = 1;}printf("\n稀疏基PSi:\n");for (int i = 0; i < 5; i++){for (int j = 0; j <5; j++){printf("%lf ", PSi.mat[i][j]);}printf("\n");}matrix x_r;initzero(&x_r, 5, 1);multiply_matrix(PSi, s, x_r);printf("\ns:\n");for (int i = 0; i < 5; i++){printf("%lf ", s.mat[i][0]);}matrix x;initzero(&x, 5, 1);x.mat[2][0] = 3.0515;printf("\n原始信号x:\n");for (int i = 0; i < 5; i++){printf("%lf ", x_r.mat[i][0]);}printf("\n恢复信号x_r:\n");for (int i = 0; i < 5; i++){printf("%lf ", x_r.mat[i][0]);}getchar();
}

结果

参考

【1】一分钟了解”匹配追踪算法（Matching Pursuit，MP） 和 正交匹配追踪算法（Orthogonal Matching Pursuit，OMP）“
【2】匹配追踪MP和正交匹配追踪OMP算法
【3】知乎：正交匹配追踪算法OMP
【4】知乎：匹配追踪算法（match pursuit, MP）实现信号的稀疏表示