通俗理解正交匹配追踪（OMP）算法及MATLAB代码实现_omp算法

2023.12.16 在阅读OTFS相关论文的过程中，总是被数学知识绊住，因此在这里从通信小白的视角尝试理解一下基础的相关算法//部分内容有参考ChatGPT。

1 背景：压缩感知问题模型

1.1 定义

压缩感知（Compressed Sensing，CS）是一种用于信号处理的技术，它利用信号的稀疏性来从远少于Nyquist采样定理所需的样本数中重构信号。压缩感知的核心是解决一个优化问题，通常表述为以下形式：

$$\underset{\mathbf{x}}{\min }\Vert \mathbf{x}{\Vert }_0\text{subject to}\mathbf{y}=\mathbf{Ax}\text{(1)}$$

符号含义：

• $\mathbf{x}$: $N\times 1$ 维稀疏信号向量。在大多数情况下，$\mathbf{x}$ 表示原始信号，其大部分元素都是0或接近0（即稀疏的）。
• $\Vert \mathbf{x}{\Vert }_0$: $\mathbf{x}$的$0$-范数，表示$\mathbf{x}$中非零元素的数量。压缩感知的目的是最小化这个数量，寻找最稀疏的解。
• $\mathbf{y}$: $M\times 1$ 维观测向量。这是通过对原始信号进行少量测量获得的，其中$M\ll N$。
• $\mathbf{A}$: $M\times N$ 维测量矩阵。这个矩阵用于从高维信号空间$N$转换到低维测量空间$M$。理想的测量矩阵应满足一定的条件（如限制等距性质），以确保信号能够从较少的测量中准确重建。

1.2 解决方法

直接解决上述优化问题是NP难的，因此通常采用以下方法近似解决：

1. $l_1$-范数松弛:

   • 替换$0$-范数为$1$-范数（即使用$\Vert \mathbf{x}{\Vert }_1$），将问题转化为可通过各种有效算法求解的凸优化问题。
   • 优化问题变为:
     $$\underset{\mathbf{x}}{\min }\Vert \mathbf{x}{\Vert }_1\text{subject to}\mathbf{y}=\mathbf{Ax}$$

2. 贪婪算法:

   • 如匹配追踪（Matching Pursuit, MP）和正交匹配追踪（Orthogonal Matching Pursuit, OMP）等，这些算法通过迭代选择信号的稀疏表示来逼近原始信号。

2 简介

正交匹配追踪（OMP）算法是一种用于信号处理和稀疏信号重建的强大工具。它通过以下方式工作：首先，从字典中选择最相关的原子，然后使用它们来逼近目标信号。这种方法在多个领域都有广泛的应用，包括图像处理、压缩感知、和机器学习。

2.1 通俗解释：什么是OMP算法？

  OMP算法就像是一位画家，试图用尽可能少的颜色来绘制一幅画。它首先选择最重要的颜色（原子），然后根据这些颜色逼近整幅画（目标信号）。这种方法可以在保持画面质量的同时（重构出的信号尽可能逼近原始信号），减少颜色的使用量。

2.2 基本概念

过完备字典：字典中的原子数量多于信号的维度。在正交匹配追踪（OMP）算法中，过完备字典通常用于处理高维信号，其中信号的维度小于字典中的原子数量。这种情况下，信号无法由字典中的一个原子线性表示，需要通过组合多个原子来逼近信号。过完备字典的使用可以提高信号重建的灵活性，允许更精确地表示和重建信号，特别是在信号稀疏性较高的情况下。
稀疏表示：稀疏表示意味着，虽然你有非常多的原子可用，但是你只需要很少的一部分就能准确地表示一个特定的信号。这就好比用最少的词汇来准确表达一个复杂的想法。

2.3 如何工作？

1. 初始化：首先，我们有一幅待绘制的画（目标信号）和一盘各种颜色的调色板（原子库）。

2. 原子选择：然后，画家会挑选出最适合绘制当前部分画面的颜色（原子），通常是最能突出特征的颜色。

3. 绘制：画家使用选中的颜色来绘制一部分画面，并观察结果。

4. 更新：如果需要，画家会再次挑选颜色，绘制下一部分画面，然后不断重复这个过程，直到整幅画完成。

3 数学公式推导

3.1 问题设定

我们可以用数学公式来描述OMP算法。首先，我们有以下关系：

$$\mathbf{b}=\mathbf{Dx}\text{(1)}$$

其中，$\mathbf{b}$是目标信号，$\mathbf{D}$是原子库/过完备字典矩阵（每个列向量称为原子），$\mathbf{x}$是稀疏表示。

3.2 目标函数

我们的目标是找到一个稀疏表示$\mathbf{x}$，使得下面的目标函数最小化：

$$\underset{\mathbf{x}}{\min }\Vert \mathbf{b}-\mathbf{Dx}{\Vert }_2^2\text{subject to}\Vert \mathbf{x}{\Vert }_0\le K\text{(2)}$$

这个目标函数的意思是，我们要找到尽量少的原子${\mathbf{d}}_i$的线性组合，使得与目标信号$\mathbf{y}$的误差最小。

3.3 算法步骤

详细解释：
1.初始化：我们开始时什么都不知道，所以初始残差${\mathbf{r}}^0$等于目标信号$\mathbf{b}$，选中原子集合${\Lambda }_0$为空，迭代次数$k$为0。
5.原子选择：我们选择一个与当前残差最相关的原子，即内积$|⟨{\mathbf{a}}_j,{\mathbf{r}}_{k-1}⟩|$最大的原子${\lambda }_k$。
6. 更新选中原子集合：将选中的原子${\lambda }_k$加入选中原子集合${\Lambda }^k$中。
7. 求解最小二乘问题：通过最小二乘法求解稀疏表示${\mathbf{x}}^k$。
8/9. 更新残差：更新残差${\mathbf{r}}^{k+1}$，计算$\mathbf{y}$与${\mathbf{D}}_{{\Lambda }^{k+1}}{\mathbf{x}}^{k+1}$之间的差异。
迭代条件：如果未达到稀疏度$K$或残差足够小，则增加迭代次数$k$，重复步骤2-5。

4 对比MP与OMP算法

4.1 MP算法的基本思想：

该算法是一种迭代的贪婪策略，旨在从一个预先定义的过完备字典中选择原子来逐步构建信号的稀疏近似表示。在每次迭代中，算法选择与当前残差最相关的原子（即，使得该原子与残差的内积最大），构建一个稀疏逼近，并求出信号残差，然后继续选择与信号残差最匹配的原子，反复迭代。通过这种方式，MP算法能够逐步减少残差的能量，直到满足特定的停止条件，从而实现信号的有效稀疏分解。

4.2 MP算法的缺点：

1. 非正交选择：MP算法在每一步只是简单地选择与残差相关性最高的字典元素，但不从残差中移除该元素的全部贡献，导致后续迭代中原子之间可能会存在重叠或共享的成分。这意味着算法效率低下，因为它可能需要更多的迭代才能达到相同的表示准确性。
2. 冗余表示：由于字典元素间的非正交性，MP算法在表示信号时可能会引入冗余，使得最终的稀疏表示并不是最优的。
3. 收敛速度：MP算法的收敛速度可能比较慢，因为每次只减去了所选元素在当前残差上的投影，而不是整个信号。
4. 算法稳定性：MP算法可能受到噪声的影响更大，因为它没有考虑到字典元素间的相互关系，容易在噪声的影响下选择错误的元素。

4.3 OMP算法的优势：

6. 正交化过程：OMP算法在每次迭代中，选择与残差最相关的字典元素后，会进行一个正交化过程，即更新所有已选择元素的系数，以确保它们与新的残差正交。这避免了冗余，并提高了算法的效率。
7. 更优的稀疏表示：由于正交化过程，OMP能够提供更加准确且稀疏的信号表示。
8. 更快的收敛：相比MP，OMP算法通常需要更少的迭代步骤就能收敛到一个稳定的解，因为它在每一步都减少了更多的信号能量。
9. 稳定性和鲁棒性：正交化使得OMP算法对噪声具有更强的鲁棒性，因为它能更好地区分信号与噪声，减少误选。

5 应用举例

5.1 图像压缩

OMP算法可以用于图像压缩。它选择一小部分重要的像素点（原子），并使用它们来重建整幅图像，从而减小图像文件的大小，同时保持图像质量。

5.2 图像去噪

通过选择最相关的图像特征（原子），重建图像以去除噪声，从而获得清晰的图像。

5.3 信号处理

在信号处理中，OMP算法可以用于从嘈杂的测量中恢复原始信号。它帮助我们找到最相关的信号分量，从而更准确地还原原始信号。

5.4 机器学习

OMP算法在稀疏表示学习中也有广泛应用。它可以用于特征选择，帮助机器学习模型识别最重要的特征。

6 代码实现（matlab）

%% 该代码来自chatgpt4.0，Cuby!在该基础上进行修改
% 初始化参数
N = 100; % 信号维度
M = 50;  % 观测维度(测量矩阵A的行数)
K = 10;  % 稀疏度：生成稀疏信号x，这里假设只有10个非零元素
% 生成稀疏信号 x
x = zeros(N, 1); % 初始化 x 为 N 维零向量
nonzero_indices = randperm(N, K); % 随机选择 K 个不同的索引
x(nonzero_indices) = randn(K, 1); % 将这些索引处的 x 值设为随机数

% 生成测量矩阵A
A = randn(M, N); % 测量矩阵用于从原始信号 x 中产生观测信号 y

% 生成测量信号y
y = A * x;

% 使用OMP算法恢复稀疏信号
x_hat = omp(A, y, M);

% 显示原始信号和估计信号
figure;
subplot(2, 1, 1);
stem(x, 'r', 'LineWidth', 1.5, 'MarkerSize', 5);
title('原始信号 x');

subplot(2, 1, 2);
stem(x_hat, 'b', 'LineWidth', 1.5, 'MarkerSize', 5);
title('估计信号 x_hat');

% 计算估计误差
mse = mean((x - x_hat).^2);
snr = 10 * log10(var(x) / mse);
relative_error = norm(x - x_hat) / norm(x);

% 打印误差结果
disp(['MSE: ', num2str(mse)]);
disp(['SNR: ', num2str(snr), ' dB']);
disp(['相对误差: ', num2str(relative_error)]);

% OMP算法实现
function x_hat = omp(A, y, M)
    N = size(A, 2);
    x_hat = zeros(N, 1);
    residual = y;
    selected_indices = [];
    tolerance = 1e-6; % 设置残差阈值

    for i = 1:M  % 为避免过拟合，迭代次数不能超过观测数量
        % 计算测量残差的投影
        projections = abs(A' * residual);

        % 选择最大投影对应的索引
        [~, index] = max(projections);

        % 添加选定的原子索引
        selected_indices = [selected_indices, index];

        % 更新估计的稀疏信号
        x_hat(selected_indices) = A(:, selected_indices) \ y;

        % 更新残差
        residual = y - A * x_hat;

        % 检查残差是否足够小//停止条件
        if norm(residual) < tolerance
            break;
        end
    end
end
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69

仿真结果如下：

参考文献

[1] J. Wang, S. Kwon and B. Shim, “Generalized Orthogonal Matching Pursuit,” in IEEE Transactions on Signal Processing, vol. 60, no. 12, pp. 6202-6216, Dec. 2012, doi: 10.1109/TSP.2012.2218810.
[2] https://angms.science/doc/RM/OMP.pdf
[3] Chatgpt 4.0