Orthogonal Matching Pursuit (OMP)
Сжатое зондирование и восстановление сигнала (OMP)
В этом примере мы рассматриваем классическую задачу обработки сигналов: имеется сигнал большой длины n = 256, но известно, что он разрежен — только k = 20 его отсчётов отличны от нуля, остальные равны нулю. Вместо того чтобы измерять все 256 значений, мы делаем лишь m = 128 линейных измерений вида y = A·x, где A — случайная гауссова матрица размера m × n. Задача — по «сжатым» измерениям y восстановить исходный сигнал x.
В чём суть технологии. Классическая теория (теорема Котельникова) требует числа измерений порядка длины сигнала. Однако если сигнал разрежен, то теория сжатого зондирования (Compressed Sensing) гарантирует точное восстановление уже при m ≈ k·log(n), то есть при m ≪ n. Ключевая идея — использовать нелинейный разреженный решатель вместо обычного МНК.
Алгоритм OMP (Orthogonal Matching Pursuit, ортогональное преследование согласованием) — один из самых простых жадных алгоритмов восстановления. На каждой итерации он:
- находит столбец матрицы
A, наиболее коррелирующий с текущим остаткомr = y − A·x; - добавляет его индекс в поддерживаемое множество;
- решает задачу наименьших квадратов только на выбранных столбцах;
- обновляет остаток.
Через k итераций алгоритм возвращает разреженное решение с k ненулевыми компонентами. В нашем примере он восстанавливает сигнал с относительной ошибкой порядка 10⁻¹⁵ — то есть практически идеально, что наглядно демонстрируется графиками исходного и восстановленного сигналов.
Реализация
В примере используются три библиотеки:
- LinearAlgebra — стандартная библиотека линейной алгебры: матричные операции, нормы, решение линейных систем (используется для
norm,A',Asub \ y). - SparseArrays — стандартная библиотека для работы с разреженными векторами и матрицами (используется для
sparsevec,findnz, эффективного хранения сигнала с малым числом ненулевых элементов). - Random — стандартная библиотека генерации случайных чисел (используется для
randn,randpermиRandom.seed!для воспроизводимости).
Pkg.add(["LinearAlgebra", "SparseArrays", "Random"])
using LinearAlgebra, SparseArrays, Random
omp(A, y, k) — реализация алгоритма Orthogonal Matching Pursuit (OMP) для жадного восстановления разреженного вектора 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.
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
Подготовка данных для сжатого зондирования:
- Задаём параметры:
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.
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.
x_rec = omp(A, y, k);
После этого вычисляем относительную ошибку восстановления, выводим её, число ненулевых элементов в исходном и восстановленном сигналах и проверяем совпадение их поддержек.
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)))
Далее строим график исходного сигнала x_true и поверх него пунктиром красного цвета график восстановленного сигнала x_rec для их сравнения.
plot(x_true, label="Исходный", ylabel="Амплитуда",
title="Сравнение сигналов", linewidth=2)
plot!(x_rec, label="Восстановленный",
linewidth=2, linestyle=:dash, color=:red)
Вывод
В этом примере успешно решена задача сжатого зондирования: по неполным измерениям y = A·x (где m = 128 < n = 256) восстановлен разреженный сигнал x с k = 20 ненулевыми компонентами. Алгоритм OMP показал высокую точность:
- относительная ошибка восстановления составила
3.77e-16, что близко к машинной точности; - число ненулевых элементов в исходном и восстановленном сигналах совпадает (20);
- поддержка (позиции ненулевых элементов) восстановлена верно.
Таким образом, пример наглядно подтверждает основные положения теории сжатого зондирования: при достаточной разреженности сигнала его можно точно восстановить по существенно меньшему числу измерений, а жадный алгоритм OMP является эффективным инструментом для решения такой задачи.