壓縮感知初入門小白，如有錯誤，歡迎指正交流~

原理

在上篇博文中，指出當前壓縮感知中常用的重構演算法，其中貪婪演算法中有使用廣泛的兩種演算法，即匹配追蹤演算法（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 1
	   and 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 this
				   matrix 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）實現信號的稀疏表示