AnyMath 文档
Notebook

压缩传感和信号恢复(OMP)

在这个例子中，我们正在考虑一个经典的信号处理问题：存在一个长长度的信号。 n = 256 但已知它是稀疏的-只 k = 20 它的计数非零，其余为零。 我们只测量所有256个值，而不是测量所有256个值 m = 128 表格的线性测量 y = A·x，在哪里 A -大小的随机高斯矩阵 m × n. 该任务基于"压缩"测量 y 恢复原始信号 x.

技术的本质是什么。 经典理论（Kotelnikov定理）需要对信号长度的量级进行多次测量。 然而，如果信号是稀疏的，压缩传感的理论保证了准确的恢复已经在 m ≈ k·log(n) 也就是说，当 m ≪ n. 关键思想是使用非线性稀疏求解器而不是传统的OLS。

Omp算法（正交匹配算法）是最简单的贪婪恢复算法之一。 在每次迭代中，他:

  1. 查找矩阵的列 A，与当前余额最相关 r = y − A·x;
  2. 将其索引添加到支持的集;
  3. 仅在选定列上求解最小二乘问题;
  4. 更新余额。

通过 k 算法的迭代返回一个稀疏解 k 非零分量。 在我们的例子中，它以相对顺序误差恢复信号。 10⁻¹⁵ -也就是说，它几乎是完美的，这是由原始和恢复信号的图表清楚地表明。

实现；实现

该示例使用三个库:

-LinearAlgebra-线性代数的标准库：矩阵运算，规范，求解线性系统（用于 norm, A', Asub \ y).
-SparseArrays是用于处理稀疏向量和矩阵的标准库（用于 sparsevec, findnz，有效地存储具有少量非零元素的信号）。
-Random-用于生成随机数的标准库（用于 randn, randperm 和 Random.seed! 于再现性）。

In [ ]: 
Pkg.add(["LinearAlgebra", "SparseArrays", "Random"])
using LinearAlgebra, SparseArrays, Random

omp(A, y, k) -贪婪稀疏向量重建的正交匹配追求算法的实现 x 从方程 y = A·x，在哪里 k -非零元素的预期数量。

算法:

-初始化: r = y, support = [], x = zeros(n).
-重复 k 一次:
-计算相关性 A' * r 并采取他们的模块;
-重置已选择的索引;
-选择索引 j 具有最大相关性并将其添加到 support;
-建立子矩阵 Asub = A[:, support];
-解决了跨国公司的问题 Asub \ y，获取系数 coef;
-更新 x：重置和写入 coef 的位置上 support;
-重新计算余额 r = y - Asub * coef;
-如果 norm(r) < 1e-12，中断循环。
-回报 x.

In [ ]: 
function omp(A, y, k)
    m, n = size(A)
    r = copy(y)
    support = Int[]
    x = zeros(n)
    for _ in 1:k
        c = abs.(A' * r)
        c[support] .= 0
        j = argmax(c)
        push!(support, j)
        Asub = A[:, support]
        coef = Asub \ y
        x .= 0
        x[support] = coef
        r = y - Asub * coef
        if norm(r) < 1e-12
            break
        end
    end
    return x
end
Out[0]:
omp (generic function with 1 method)

为压缩传感准备数据:

-设置参数: n = 256, m = 128, k = 20，在哪里 m < n 提供压缩。
-固定发电机: Random.seed!(42) 于再现性。
-创建稀疏信号: x_true = zeros(n)，选择 k 通过随机索引 randperm(n)[1:k] 我们把它写下来 randn(k).
-我们正在建立一个测量矩阵 A = randn(m, n) 我们规范化它的列: A = A ./ sqrt.(sum(A.^2, dims=1)) -使OMP不会给具有更多能量的柱带来优势。
-计算压缩测量: y = A * x_true.

In [ ]: 
n = 256                # 信号长度
m = 128                # 测量次数(m<n)
k = 20                 # 稀疏性（非零元素数）

Random.seed!(42)
x_true = zeros(n)
support = randperm(n)[1:k]
x_true[support] = randn(k)
A = randn(m, n)
A = A ./ sqrt.(sum(A.^2, dims=1))
y = A * x_true;

计算重建信号: x_rec = omp(A, y, k)，参数:

  • A -测量矩阵;
  • y -测量向量;
  • k -稀疏。

函数返回 x_rec 长度 n 不超过 k 非零分量。 接下来，我们比较 x_rec 与 x_true.

In [ ]: 
x_rec = omp(A, y, k);

之后，我们计算相对恢复误差，输出它，原始和恢复信号中非零元素的数量，并检查它们的支持是否匹配。

In [ ]: 
err = norm(x_true - x_rec) / norm(x_true)
println("相对恢复错误: ", err)
println("原始信号中非零元素的数量: ", k)
println("重建的非零元素的数目: ",
        count(!iszero, x_rec))
println("支持已正确恢复: ",
        Set(findall(!iszero, x_true)) == Set(findall(!iszero, x_rec)))

Относительная ошибка восстановления: 3.772586544601344e-16
Число ненулевых элементов в исходном сигнале: 20
Число ненулевых элементов в восстановленном: 20
Поддержка восстановлена верно: true

接下来，我们绘制原始信号 x_true 在它的顶部，恢复信号的图形用红色点缀。 x_rec 来比较它们。

In [ ]: 
plot(x_true, label="原版", ylabel="的振幅",
          title="信号比较", linewidth=2)
plot!(x_rec, label="修复后",
      linewidth=2, linestyle=:dash, color=:red)
Out[0]:

结论

在这个例子中，压缩传感的问题已经成功地解决了：通过不完整的测量 y = A·x （哪里 m = 128 < n = 256）已恢复稀疏信号 x 与 k = 20 非零分量。 OMP算法具有很高的准确性:

-相对恢复错误是 3.77e-16 接近机器精度;
-原始和重建信号中非零元素的数量相同（20）;
-支持（非零元素的位置）已正确恢复。

因此，该示例清楚地证实了压缩感测理论的主要规定：在信号足够稀疏的情况下，可以从显着较少数量的测量中精确重建，贪婪OMP算法是解决此类问题的有