当前位置：首页 > news >正文

正交匹配追踪（OMP）详解：压缩感知的基石算法

news 2025/9/9 13:16:38

正交匹配追踪（OMP）详解：压缩感知的基石算法

一、引言

在现代信号处理、通信和图像处理等领域，稀疏表示（Sparse Representation） 和 压缩感知（Compressive Sensing, CS） 已成为重要的研究方向。其核心思想是：如果一个信号在某个域内是稀疏的（即绝大部分元素为零或接近零），那么我们可以通过远少于奈奎斯特采样定律要求的观测数据，实现对该信号的准确重构。这一理论大大降低了采样成本，并在雷达成像、无线通信信道估计、图像压缩与恢复等场景中展现出广阔应用前景。

在众多稀疏重构算法中，正交匹配追踪（Orthogonal Matching Pursuit, OMP） 因其 实现简单、直观易懂、计算量相对较低 而被广泛采用。OMP 属于贪婪类算法，其基本思想是：从观测数据中逐步选取与当前残差最相关的基向量，通过迭代方式不断逼近真实稀疏信号。虽然OMP在精度上可能不如某些基于优化的算法（如L1范数最小化方法），但其在工程实践中具有较高性价比，尤其适用于稀疏度较低、实时性要求较高的场景。

本文将围绕 OMP 算法展开详细介绍，包括：

算法原理与数学推导：理解 OMP 如何一步步重构稀疏信号；
数值示例：通过具体例子直观演示迭代过程；
MATLAB 代码实现与解析：展示算法实现细节；
优缺点与应用分析：理解 OMP 的优势与局限。

通过本文，你将能够从 理论 + 实践 + 工程实现 三个维度全面掌握 OMP 算法，并为后续在压缩感知及相关领域的深入研究打下坚实基础。

二、OMP算法原理

2.1 问题背景与数学模型

在压缩感知框架中，我们常遇到如下稀疏重构问题：

$n,\quad h \in \mathbb{R}^K,\; \|h\|_0 = S$

其中：

$\in \mathbb{R}^N$ ：观测向量；
$\in \mathbb{R}^{N \times K}$ ：测量矩阵或字典矩阵（通常 $K > N$ ）；
$\in \mathbb{R}^K$ ：待恢复的稀疏信号；
$n$ ：观测噪声；
$h\|_0 = S$ ：信号稀疏度（非零元素个数）。

OMP（Orthogonal Matching Pursuit）的目标是通过迭代选择与残差最相关的列向量，逐步恢复稀疏信号 $h$ 。

2.2 算法思想

OMP 属于 贪婪类算法（Greedy Algorithm），其基本思路是：

初始化残差为观测信号；
每次迭代选择与残差最相关的列向量；
更新支持集；
使用最小二乘法估计稀疏信号；
更新残差，直至满足停止条件。

2.3 算法步骤与公式推导

设迭代次数为 $k$ ，初始化支持集 $Ω(0)=∅\Omega^{(0)} = \varnothing$ ，稀疏信号初值 $h^{(0)} = 0$ ，残差 $r^{(0)} = y$ 。

每次迭代包含以下步骤：

相关性计算
$c_j = \big| P_j^\top r^{(k)} \big|,\quad j=1,\dots,K$
选择最大相关列
$j^\star = \arg\max_j c_j$
更新支持集
$\Omega^{(k+1)} = \Omega^{(k)} \cup \{j^\star\}$
最小二乘估计
$h_{\Omega^{(k+1)}} = \big(P_{\Omega^{(k+1)}}^\top P_{\Omega^{(k+1)}}\big)^{-1} P_{\Omega^{(k+1)}}^\top y$
其余位置 $j∉Ω(k+1)h_j = 0,\; j\notin \Omega^{(k+1)}$ 。
更新残差
$r^{(k+1)} = y - P_{\Omega^{(k+1)}} h_{\Omega^{(k+1)}}$
停止条件
- 达到最大稀疏度 $S$ ；或残差满足 $∥r(k+1)∥2≤ε⋅∥y∥2\|r^{(k+1)}\|_2 \leq \varepsilon \cdot \|y\|_2$ 。

三、数值例子（详细迭代过程）

下面给出一个小规模、可手算的数值例子，展示 OMP 的每次迭代：相关性计算、索引选择、最小二乘求解、残差更新。为便于手算，我选取了简单的测量矩阵和稀疏向量（均为有理数），并设置稀疏度上限 $S = 3$ 。

3.1 问题设定

设观测维度 $N = 4$ ，字典列数 $K = 6$ 。取测量矩阵 $P=[P1,…,P6]P=[P_1,\dots,P_6]$ （每列为一个长度为 4 的向量）：

$\begin{bmatrix} 1 & 0 & 0 & 0 & 1 & 0\\[4pt] 0 & 1 & 0 & 0 & 1 & 1\\[4pt] 0 & 0 & 1 & 0 & 0 & 1\\[4pt] 0 & 0 & 0 & 1 & 0 & 0 \end{bmatrix},$
即
$\begin{aligned} P_1 &= [1,0,0,0]^\top,\quad P_2=[0,1,0,0]^\top,\quad P_3=[0,0,1,0]^\top,\\ P_4 &= [0,0,0,1]^\top,\quad P_5=[1,1,0,0]^\top,\quad P_6=[0,1,1,0]^\top. \end{aligned}$

真实（稀疏）信号 $htrueh^\text{true}$ 取三个非零元素（索引 1、3、6）：
$h^\text{true} = [\,2,\;0,\;-1,\;0,\;0,\;\tfrac12\,]^\top.$

观测向量（无噪声）：
$P\,h^\text{true}.$

逐项计算得到：
$2P_1 + (-1)P_3 + \tfrac12 P_6 = \begin{bmatrix}2\\[2pt]\tfrac12\\[2pt]-\tfrac12\\[2pt]0\end{bmatrix}.$

我们用 OMP（上限 $S = 3$ 、残差阈值可忽略）从 $y$ 恢复 $h$ ，下面展示每次迭代的详细计算。

3.2 迭代 0（初始化）

初始化残差 $0]⊤r^{(0)}=y=[2,\; \tfrac12,\; -\tfrac12,\;0]^\top$ ；
计算每列与残差的相关性（内积的绝对值）：
$c_j = |P_j^\top r^{(0)}|.$
逐项计算：
- $c1=2P_1^\top r^{(0)} = 2 \quad\Rightarrow\; c_1=2$ .
- $c2=12P_2^\top r^{(0)} = \tfrac12 \quad\Rightarrow\; c_2=\tfrac12$ .
- $c3=12P_3^\top r^{(0)} = -\tfrac12 \quad\Rightarrow\; c_3=\tfrac12$ .
- $c4=0P_4^\top r^{(0)} = 0 \quad\Rightarrow\; c_4=0$ .
- $c5=52P_5^\top r^{(0)} = P_1^\top r^{(0)} + P_2^\top r^{(0)} = 2 + \tfrac12 = \tfrac52 \quad\Rightarrow\; c_5=\tfrac52$ .
- $c6=0P_6^\top r^{(0)} = P_2^\top r^{(0)} + P_3^\top r^{(0)} = \tfrac12 + (-\tfrac12)=0 \quad\Rightarrow\; c_6=0$ .
最大相关性出现在列 5（ $c5=52c_5=\tfrac52$ ），因此选择索引 $j⋆=5j^\star=5$ 。
支持集： $Ω(1)={5}\Omega^{(1)}=\{5\}$ 。

最小二乘解（单列情况）：
$h_{5}^{(1)} = \frac{P_5^\top y}{P_5^\top P_5} = \frac{\tfrac52}{2} = \frac{5}{4}.$

残差更新：
$r^{(1)} = y - P_5 h_5^{(1)} = \begin{bmatrix}2\\[2pt]\tfrac12\\[2pt]-\tfrac12\\[2pt]0\end{bmatrix} -\frac{5}{4}\begin{bmatrix}1\\[2pt]1\\[2pt]0\\[2pt]0\end{bmatrix} = \begin{bmatrix}\tfrac34\\[2pt]-\tfrac34\\[2pt]-\tfrac12\\[2pt]0\end{bmatrix}.$

3.3 迭代 1

当前残差 $0]⊤r^{(1)}=[\tfrac34,\,-\tfrac34,\,-\tfrac12,\,0]^\top$ 。

计算相关性 $cj=∣Pj⊤r(1)∣c_j = |P_j^\top r^{(1)}|$ ：
- $c1=∣P1⊤r(1)∣=34c_1 = |P_1^\top r^{(1)}| = \tfrac34$ .
- $c2=∣P2⊤r(1)∣=34c_2 = |P_2^\top r^{(1)}| = \tfrac34$ .
- $c3=∣P3⊤r(1)∣=12c_3 = |P_3^\top r^{(1)}| = \tfrac12$ .
- $c_4 = 0$ .
- $c5=∣P5⊤r(1)∣=r1+r2=34+(−34)=0c_5 = |P_5^\top r^{(1)}| = r_1 + r_2 = \tfrac34 + (-\tfrac34) = 0$ .
- $∣=54c_6 = |P_6^\top r^{(1)}| = r_2 + r_3 = -\tfrac34 + (-\tfrac12) = -\tfrac54 \Rightarrow |\,\cdot\,|=\tfrac54$ .
最大相关性出现在列 6（ $c6=54c_6=\tfrac54$ ），选择 $j⋆=6j^\star=6$ 。
支持集： $Ω(2)={5,6}\Omega^{(2)}=\{5,6\}$ 。

在支持集 ${5,6\}$ 上做最小二乘 ：

构造矩阵与向量（记 $P_{5,6}=[P_5\; P_6]$ ）：
$P_{5,6}^\top P_{5,6} = \begin{bmatrix}P_5^\top P_5 & P_5^\top P_6\\[4pt] P_6^\top P_5 & P_6^\top P_6\end{bmatrix} =\begin{bmatrix}2 & 1\\[4pt] 1 & 2\end{bmatrix}, \quad P_{5,6}^\top y = \begin{bmatrix} \tfrac52 \\[4pt] 0 \end{bmatrix}.$

对 $2×22\times2$ 矩阵求逆（公式）：
$\begin{bmatrix}2 & 1\\ 1 & 2\end{bmatrix}^{-1} = \frac{1}{(2\cdot2-1\cdot1)} \begin{bmatrix}2 & -1\\ -1 & 2\end{bmatrix} = \frac{1}{3}\begin{bmatrix}2 & -1\\ -1 & 2\end{bmatrix}.$

因此最小二乘解为：
$\begin{bmatrix} h_5^{(2)} \\[4pt] h_6^{(2)} \end{bmatrix} = \frac{1}{3}\begin{bmatrix}2 & -1\\ -1 & 2\end{bmatrix} \begin{bmatrix}\tfrac52\\[4pt] 0\end{bmatrix} =\begin{bmatrix}\tfrac{5}{3}\\[4pt] -\tfrac{5}{6}\end{bmatrix}.$

残差更新：
计算 $r^{(2)} = y - P_5 h_5^{(2)} - P_6 h_6^{(2)}$ 。代入可得（化简为有理数）：
$r^{(2)} = \begin{bmatrix}\tfrac13\\[4pt] -\tfrac13\\[4pt] \tfrac13\\[4pt] 0\end{bmatrix}.$

（可以验证：每分量均为 $±13\pm\frac{1}{3}$ 或 0。）

3.4 迭代 2

当前残差 $0]⊤r^{(2)}=[\tfrac13,\;-\tfrac13,\;\tfrac13,\;0]^\top$ 。

相关性计算：
$c_1 = c_2 = c_3 = \big|P_j^\top r^{(2)}\big| = \tfrac13,\quad c_4=c_5=c_6=0.$
三个列（1、2、3）相关性相同（并列最大）。在实际实现中，若出现并列，常按“第一个最大值”的策略选择；这里为演示我们选择 列 1（也可以选择 2 或 3，结果可能不同，但会亦可达到零残差——见下面讨论）。
选择 $j⋆=1j^\star=1$ ，支持集 $Ω(3)={5,6,1}\Omega^{(3)}=\{5,6,1\}$ 。

在支持集 ${5,6,1\}$ 上做最小二乘：

构造 $3×33\times3$ 正定矩阵 $P_\Omega^\top P_\Omega$ 与 $P_\Omega^\top y$ ：
$\begin{bmatrix} 2 & 1 & 1\\[4pt] 1 & 2 & 0\\[4pt] 1 & 0 & 1 \end{bmatrix}, \qquad b = \begin{bmatrix}\tfrac52\\[4pt]0\\[4pt]2\end{bmatrix}.$

解线性方程 $A x = b$ ，得到（直接给出解）：
$\begin{bmatrix} h_5^{(3)} \\[4pt] h_6^{(3)} \\[4pt] h_1^{(3)} \end{bmatrix} =\begin{bmatrix}1\\[4pt] -\tfrac12\\[4pt]1\end{bmatrix}.$

于是对应的整个估计向量 $h^\hat h$ （按索引 1…6）为：
$h^=[ h1(3), 0, 0, 0, h5(3), h6(3) ]⊤=[ 1, 0, 0, 0, 1, −12 ]⊤. \hat h = \big[\, h_1^{(3)},\;0,\;0,\;0,\; h_5^{(3)},\; h_6^{(3)}\,\big]^\top = \big[\,1,\;0,\;0,\;0,\;1,\; -\tfrac12\,\big]^\top.$

检验残差：
$P_1\cdot1 + P_5\cdot1 + P_6\cdot(-\tfrac12) = \begin{bmatrix}2\\[2pt]\tfrac12\\[2pt]-\tfrac12\\[2pt]0\end{bmatrix} = y,$
因此残差精确为 0，算法可在此停止（或在达到 $S = 3$ 后停止）。

3.5 结果与说明

OMP 在 $3$ 次迭代后得到的估计为
$h^=[ 1, 0, 0, 0, 1, −12 ]⊤, \hat h = [\,1,\;0,\;0,\;0,\;1,\; -\tfrac12\,]^\top,$
满足 $Ph^=yP\hat h = y$ ，残差为零（完美重构）。
注意 原始稀疏向量 与 OMP 得到的向量 不同：
$]⊤而h^=[ 1, 0, 0, 0, 1, −12 ]⊤. h^\text{true} = [\,2,\;0,\;-1,\;0,\;0,\;\tfrac12\,]^\top \quad\text{而}\quad \hat h = [\,1,\;0,\;0,\;0,\;1,\; -\tfrac12\,]^\top.$
这说明在冗余字典（列之间存在线性关系）的情形下，稀疏表示可能并非唯一：不同的支持集与系数组合可以产生相同的观测 $y$ 。在本例中，列之间存在简单线性关系（例如 $P_5=P_1+P_2$ ， $P_6=P_2+P_3$ ），因而出现了这一情况。
另外要注意 并列最大相关性的处理（第 3 次迭代时 $c_1=c_2=c_3$ ）：实现上不同的“破平策略”会影响 OMP 最终选择的支持集，但只要字典满足适当的不相关性条件（例如满足一定的互相关上界或 RIP），OMP 能恢复出原始稀疏支持的唯一解；否则可能得到“等价”的不同稀疏解。

3.6 小结（例子给出的启发）

此例展示了 OMP 的典型流程：相关性选择 → 支持集扩增 → 最小二乘拟合 → 残差更新，并逐步逼近观测向量。
在冗余或高相关的字典下，OMP 可能找到与真实支持不同但等价的稀疏表示；因此在工程实践中要注意字典设计（降低互相关）或限制稀疏度以保证唯一性。
该手算示例适合教学与调试：由于矩阵元素为小有理数，所有中间量都能写成分数，便于逐步验证程序实现是否正确。

四、MATLAB代码实现与解析

本节给出一份 正交匹配追踪（OMP）算法的 MATLAB 实现，该代码为核心功能版本，直接来自于常见的压缩感知重构流程。我们将在代码中添加详细的中文注释，并逐步解释每一部分的作用，使读者能够从数学公式与实现之间建立清晰联系。

4.1 OMP MATLAB 实现

```matlab
function [h_i, Omega] = OMP(y,P,S)
% OMP - 正交匹配追踪算法
% 输入:
%   y : N×1 观测信号向量
%   P : N×K 测量矩阵（字典矩阵）
%   S : 稀疏度上限（最多迭代次数）
% 输出:
%   h_i   : K×1 重构的稀疏信号向量
%   Omega : 被选中的原子（列索引）集合N = size(y,1);        % 观测信号维度
K = size(P,2);        % 字典列数（信号长度）
err = 1e-3;           % 迭代停止阈值（相对残差）h_i = zeros(K, 1);    % Step 0: 初始化重构信号
y_i = y - P * h_i;    % Step 0: 初始化残差 r=yit = 0;               % 迭代计数器
stop = 0;             % 停止标志
P_s = [];             % 支持集对应的子矩阵
pos_array = [];       % 记录选中的列索引while ~stop% Step 1: 计算残差与字典列的相关性for ntl = 1:1:Kproduct(ntl) = abs(P(:,ntl)' * y_i); % 计算每一列与残差的内积并取绝对值endy_i_before = y_i;     % 保存更新前的残差[val,pos] = max(product);  % Step 2: 找到相关性最大的列% Step 3: 更新支持集矩阵及索引记录P_s = [P_s, P(:,pos)];       pos_array = [pos_array, pos]; P(:,pos) = zeros(N,1);       % 防止重复选择同一列% Step 4: 最小二乘估计当前支持集上的系数h_s = (P_s' * P_s)^(-1) * P_s' * y;% h_s = P_s\y; % （可替代方式，更数值稳定）% Step 5: 更新残差y_i = y - P_s * h_s; it = it + 1;  % 迭代次数更新% Step 6: 判断停止条件if (it >= S) || norm(y_i) <= err * norm(y)stop = 1;end
end% Step 7: 构造最终重构信号
h_i(pos_array) = (P_s' * P_s)^(-1) * P_s' * y;    
Omega = pos_array;   % 返回支持集索引
end

4.2 代码解析

初始化
- h_i = zeros(K,1) 对应稀疏信号的初始估计 $h^(0)=0\hat{h}^{(0)} = 0$ 。
- y_i = y - P*h_i 初始化残差 $r^{(0)} = y$ 。
相关性计算
- product(ntl) = abs(P(:,ntl)' * y_i)
  对应数学公式：
  
  $cj=∣Pj⊤r(k)∣c_j = \big| P_j^\top r^{(k)} \big|$
  
  用于衡量每一列与当前残差的相关性。
选择最相关列
- [val,pos] = max(product)
  选择与残差相关性最大的列索引 $j⋆j^\star$ 。
支持集更新
- P_s = [P_s, P(:,pos)]
  累加支持集字典矩阵 $PΩP_{\Omega}$ 。
- pos_array = [pos_array,pos]
  记录支持集索引。
最小二乘估计
- h_s = (P_s' * P_s)^(-1) * P_s' * y
  对应：
  
  $hΩ=(PΩ⊤PΩ)−1PΩ⊤yh_{\Omega} = (P_\Omega^\top P_\Omega)^{-1} P_\Omega^\top y$
残差更新
- y_i = y - P_s * h_s
  对应：
  
  $r(k+1)=y−PΩ(k+1)hΩ(k+1)r^{(k+1)} = y - P_{\Omega^{(k+1)}} h_{\Omega^{(k+1)}}$
停止条件
- 若达到最大稀疏度 S，或残差 $∥r∥2≤10−3∥y∥2\|r\|_2 \leq 10^{-3}\|y\|_2$ 即停止。
信号重构
- h_i(pos_array) = (P_s' * P_s)^(-1) * P_s' * y
  最终在已选支持集上解出稀疏系数，其余位置保持为零。

4.3 算法特点与说明

贪婪选择：每次选取与当前残差最相关的一列；
支持集逐步扩展：保证每一步估计都在当前支持集内进行最优拟合；
最小二乘解法：使用 $(Ps⊤Ps)−1Ps⊤y(P_s^\top P_s)^{-1} P_s^\top y$ 实现线性最优解，但对于大规模问题可用 P_s\y 提升数值稳定性；
残差控制：通过 err 控制精度，可适应噪声场景；
列零化策略：P(:,pos)=zeros(N,1) 用于防止重复选择。

通过该代码，我们可以实现 OMP 算法的基本功能：在已知字典矩阵 P 与观测信号 y 情况下，根据稀疏度约束恢复信号系数向量 h。

五、OMP算法优缺点

正交匹配追踪（Orthogonal Matching Pursuit, OMP）是一种经典的贪婪型稀疏重构算法，在压缩感知、稀疏表示、雷达成像等领域得到广泛应用。尽管其思想简单、实现方便，但在实际工程中，OMP 也存在一定的局限性。以下从优点与缺点两方面进行分析。

5.1 优点

算法实现简单，计算步骤直观
OMP 主要包含相关性计算、列选择、最小二乘估计与残差更新四个步骤，逻辑清晰、易于编程实现。
重构精度较高（相比 MP）
相较于传统 Matching Pursuit (MP) 算法，OMP 在每次迭代中通过最小二乘估计更新已选择的系数，实现“正交化”处理，能够有效避免能量冗余累积，提高重构精度。
计算复杂度较低
对于稀疏度 $\ll K$ 的场景，OMP 的复杂度约为
$O(SKN + S^3)$
其中 $N$ 为观测长度， $K$ 为字典列数。由于 $S$ 较小，OMP 在实际应用中具有较好的计算效率。
较好噪声鲁棒性（结合停止条件）
通过设置残差阈值（如 $∥r∥≤ε∥y∥\|r\| \leq \varepsilon\|y\|$ ），OMP 能在噪声场景下提前终止迭代，避免过拟合。
便于与其他算法结合
OMP 可作为其他压缩感知算法（如 CoSaMP、StOMP）的基础框架，便于扩展与优化。

5.2 缺点

依赖字典列的相关性
当字典矩阵 $P$ 的列高度相关（相干性高）时，OMP 可能在初始迭代阶段选取错误列，从而导致后续重构误差放大。
对稀疏度 $S$ 依赖强
若稀疏度未知，OMP 难以预先确定迭代次数，若 $S$ 设置过大，可能导致过拟合；过小则导致重构不完全。
缺乏全局最优保证
OMP 属于贪婪算法，其每次迭代都是局部最优选择，并不保证最终解为全局最优解。
最小二乘解求解复杂度随迭代增加而上升
随着支持集扩展，每次迭代需重新求解最小二乘问题（复杂度 $O(k^3)$ ），在高稀疏度场景下会造成一定计算负担。
对噪声模型缺乏适应性
当观测噪声较强、且噪声分布未知时，OMP 的停止准则（基于残差阈值）难以精确设定，可能导致误选或漏选。