Теория разреженности и сжатые измерения

Представьте: вы хотите записать МРТ-снимок, но у пациента мало времени — нельзя получить все измерения. Или нужно передать сигнал через канал с ограниченной пропускной способностью. Теория compressed sensing (Candès, Romberg, Tao; Donoho, 2006) говорит: если сигнал «разреженный», нескольких линейных измерений достаточно для его точного восстановления.

Основная теорема Compressed Sensing

Постановка задачи: x ∈ ℝⁿ — сигнал (например, изображение из пикселей). A ∈ ℝ^{m×n}, m ≪ n — матрица измерений. Наблюдаем y = Ax (m << n измерений — явно недоопределённая система). Цель: восстановить x.

Условие разреженности: x имеет не более s ненулевых компонент: ||x||₀ = s ≪ n. Большинство сигналов разреженны в каком-то базисе: изображения — в вейвлет-базисе, аудио — в Фурье-базисе.

Теорема (Candès, 2006): Если A удовлетворяет условию RIP(2s, δ):

(1−δ)||x||₂² ≤ ||Ax||₂² ≤ (1+δ)||x||₂² для всех s-разреженных x

то решение x̂ задачи Basis Pursuit Denoising:

x̂ = argmin ||x||₁ subject to ||Ax − y||₂ ≤ ε

удовлетворяет: ||x̂ − x||₂ ≤ C₁ σₛ(x)/√s + C₂ε

Здесь: ||x||₁ = Σᵢ |xᵢ| (L1-норма), σₛ(x) — ошибка наилучшей s-разреженной аппроксимации. Ключевое: задача L1-минимизации (выпуклая!) решает, казалось бы, невозможную задачу восстановления из m ≪ n измерений.

Почему L1? L0-норма (подсчёт ненулевых) — NP-hard оптимизация. L2-норма — находит наименее квадратичное решение, не разреженное. L1 — «ближайшая выпуклая релаксация» L0: сохраняет разреженность, позволяет convex programming.

Случайные матрицы удовлетворяют RIP

Если Aᵢⱼ ~ N(0, 1/m) (гауссова случайная матрица), то A удовлетворяет RIP(s, δ) с вероятностью ≥ 1−ε при:

m ≥ C · s · log(n/s) / δ²

Пример: n = 10000, s = 100, δ = 0.1 → m ≥ C · 100 · log(100) · 100 ≈ 46000/k — достаточно m ≈ 500 измерений вместо n = 10000!

Физический смысл: случайная проекция «видит» все s-разреженные векторы (JL-лемма). Другие матрицы, удовлетворяющие RIP: частичная матрица Фурье (с случайно выбранными строками), матрица Адамара.

Алгоритмы восстановления

LASSO (Tibshirani, 1996): Эквивалент BPDN для регрессии: min_x (1/2)||y − Ax||₂² + λ||x||₁. Аналитического решения нет (кроме ортогональной A), но есть эффективные алгоритмы.

ISTA (Iterative Shrinkage-Thresholding): использует проксимальный оператор для L1: prox_{αλ||·||₁}(v)ᵢ = sign(vᵢ)·max(|vᵢ| − αλ, 0) (soft thresholding)

Шаг: xₜ₊₁ = prox_{αg}(xₜ − α·Aᵀ(Axₜ − y)). Сходимость: O(1/t).

FISTA: Nesterov-ускоренная версия ISTA: O(1/t²).

Matching Pursuit (OMP): жадный алгоритм: на каждом шаге выбираем «наиболее коррелированный» столбец A с остатком. Быстрее, не такая точная теория.

Применения

МРТ: k-пространство (частотные измерения) семплируется случайно по Пуассону → LASSO-реконструкция → в 4−8× меньше времени в аппарате. Первое клиническое применение — 2011 год.

Одно-пиксельная камера: измеряем случайные линейные комбинации пикселей (один фотодиодный датчик) → CS-восстановление. Позволяет работать в диапазонах, где матрица пикселей дорога (инфракрасный, терагерцовый).

Геномика: из m = 500 лабораторных тестов восстанавливаем полный генотип n = 10000 SNPs — каждый человек несёт разреженный набор мутаций.

Структурная разреженность: Group LASSO и Total Variation

Group LASSO: признаки разбиты на группы g. Штраф: λΣ_g ||β_g||₂. Группы выбираются «целиком» — либо все признаки группы ненулевые, либо все нули. Применение: временные ряды (группы лагов), нейронные сети (структурное прунинг — удаляем целые головы внимания).

Total Variation (TV) denoising: x ∈ ℝⁿ — сигнал с шумом. min_u ||y−u||₂² + λ Σᵢ|uᵢ₊₁−uᵢ|. Сохраняет скачки (края), сглаживает шум внутри плоских участков. Стандарт в обработке изображений.

Численный пример

n = 256 (длина сигнала), s = 10 (разреженность в Фурье). m = 60 случайных измерений. Гауссова матрица A (60×256). Сигнал: x с 10 ненулевыми частотами. BPDN восстанавливает x с ошибкой < 0.01 при ε = 0 (без шума). При шуме σ = 0.01: ||x̂ − x||₂ ≈ 0.08. Наивная псевдоинверсия даёт ||x̂_lstsq − x||₂ ≈ 2.1 — в 26 раз хуже.

Задание: Реализуйте CS для звукового сигнала: сгенерируйте сигнал = сумма 5 синусоид, n=1024 точки. Создайте случайную матрицу измерений A (m×1024, m=100). Восстановите через LASSO (sklearn). Постройте: (1) исходный спектр; (2) восстановленный спектр; (3) ошибку восстановления vs m (от 50 до 500). При каком m ошибка опускается ниже 1%?