压缩感知模型FOCUSS算法与Python实现：从理论到实践

一、压缩感知理论概述：突破奈奎斯特的革命

压缩感知（Compressive Sensing, CS）理论由Donoho、Candes等学者于2006年提出，其核心思想在于：若信号在某个变换域下具有稀疏性，则可通过远低于奈奎斯特采样率的非自适应线性投影实现信号的准确重建。这一理论颠覆了传统信号处理中”先采样后压缩”的范式，直接在采样阶段完成压缩，为高维数据采集（如医学影像、遥感监测）提供了革命性解决方案。

1.1 数学模型构建

压缩感知的数学框架可表示为：
[ y = \Phi x ]
其中，( x \in \mathbb{R}^N )为原始信号，( \Phi \in \mathbb{R}^{M \times N} )（( M \ll N )）为测量矩阵，( y \in \mathbb{R}^M )为观测向量。当( x )在某个基( \Psi )下稀疏（即( x = \Psi \theta )，( \theta )中非零元素极少）时，可通过求解以下优化问题重建信号：
[ \min{\theta} |\theta|_0 \quad \text{s.t.} \quad y = \Phi \Psi \theta ]
由于( \ell_0 )范数优化为NP难问题，实际中常用( \ell_1 )范数松弛：
[ \min{\theta} |\theta|_1 \quad \text{s.t.} \quad y = \Phi \Psi \theta ]

1.2 测量矩阵设计准则

测量矩阵( \Phi )需满足有限等距性质（RIP），即对任意( k )-稀疏信号( x )，存在常数( \delta_k \in (0,1) )使得：
[ (1-\delta_k)|x|_2^2 \leq |\Phi x|_2^2 \leq (1+\delta_k)|x|_2^2 ]
常用随机矩阵（如高斯矩阵、伯努利矩阵）和结构化矩阵（如部分傅里叶矩阵）均可满足RIP条件。

二、FOCUSS算法：加权最小二乘的稀疏解

FOCUSS（FOCal Underdetermined System Solver）算法由Gorodnitsky等人提出，其核心思想是通过迭代加权最小二乘法逐步逼近稀疏解。与传统( \ell_1 )优化不同，FOCUSS直接优化( \ell_p )（( 0 < p < 1 )）准范数，在每次迭代中根据当前解的幅度调整权重，强制解向稀疏方向收敛。

2.1 算法推导

给定观测方程( y = A\theta )（其中( A = \Phi \Psi )），FOCUSS在第( k )次迭代中求解：
[ \theta^{(k)} = \arg\min{\theta} \sum{i=1}^N w_i^{(k-1)} |\theta_i|^2 \quad \text{s.t.} \quad y = A\theta ]
权重更新规则为：
[ w_i^{(k-1)} = \left( |\theta_i^{(k-1)}|^2 + \epsilon \right)^{p/2 - 1} ]
其中，( \epsilon )为防止除零的小常数，( p )控制稀疏性强度（通常取( p \in [0.1, 0.5] )）。

2.2 算法步骤

1. 初始化：( \theta^{(0)} = A^T y )（最小二乘初始解）
2. 迭代更新：
   • 计算权重矩阵( W^{(k-1)} = \text{diag}\left( \left( |\theta_i^{(k-1)}|^2 + \epsilon \right)^{p/2 - 1} \right) )
   • 求解加权最小二乘问题：
     [ \theta^{(k)} = \left( A^T W^{(k-1)} A \right)^{-1} A^T W^{(k-1)} y ]
3. 终止条件：当( |\theta^{(k)} - \theta^{(k-1)}|_2 / |\theta^{(k-1)}|_2 < \tau )时停止。

三、Python实现与代码解析

以下代码使用NumPy和SciPy实现FOCUSS算法，并对比OMP（正交匹配追踪）的重建效果。

3.1 环境准备

import numpy as np
from scipy.linalg import lstsq
from sklearn.linear_model import OrthogonalMatchingPursuit
import matplotlib.pyplot as plt

3.2 FOCUSS算法实现

def focuss(A, y, p=0.5, epsilon=1e-6, max_iter=100, tol=1e-4):
    """
    FOCUSS算法实现
    参数:
        A: 测量矩阵 (M x N)
        y: 观测向量 (M,)
        p: 稀疏性参数 (0 < p < 1)
        epsilon: 防止除零的小常数
        max_iter: 最大迭代次数
        tol: 收敛阈值
    返回:
        theta: 重建的稀疏系数 (N,)
    """
    M, N = A.shape
    theta = lstsq(A, y)[0]  # 初始解（最小二乘）
    for _ in range(max_iter):
        # 计算权重
        abs_theta = np.abs(theta)
        weights = (abs_theta**2 + epsilon)**(p/2 - 1)
        W = np.diag(weights)
        # 加权最小二乘更新
        A_weighted = A @ W @ A.T
        b_weighted = A @ W @ theta
        theta_new = lstsq(A_weighted, b_weighted)[0]
        # 检查收敛
        if np.linalg.norm(theta_new - theta) / np.linalg.norm(theta) < tol:
            break
        theta = theta_new
    return theta

3.3 实验对比：FOCUSS vs OMP

# 生成稀疏信号
N = 100  # 信号维度
k = 5    # 稀疏度
x_true = np.zeros(N)
x_true[np.random.choice(N, k, replace=False)] = np.random.randn(k)
# 生成测量矩阵（高斯随机矩阵）
M = 30  # 测量数
Phi = np.random.randn(M, N)
Phi = Phi / np.linalg.norm(Phi, axis=0)  # 列归一化
y = Phi @ x_true
# FOCUSS重建
theta_focuss = focuss(Phi, y, p=0.3)
x_focuss = theta_focuss  # 假设Psi为单位矩阵
# OMP重建
omp = OrthogonalMatchingPursuit(n_nonzero_coefs=k)
omp.fit(Phi, y)
x_omp = omp.coef_
# 计算重建误差
error_focuss = np.linalg.norm(x_true - x_focuss) / np.linalg.norm(x_true)
error_omp = np.linalg.norm(x_true - x_omp) / np.linalg.norm(x_true)
print(f"FOCUSS重建误差: {error_focuss:.4f}")
print(f"OMP重建误差: {error_omp:.4f}")
# 可视化
plt.figure(figsize=(10, 6))
plt.stem(x_true, label='原始信号', markerfmt='bo')
plt.stem(np.arange(N), x_focuss, label='FOCUSS重建', markerfmt='rx')
plt.stem(np.arange(N), x_omp, label='OMP重建', markerfmt='g^')
plt.legend()
plt.title('压缩感知重建对比')
plt.show()

3.4 结果分析

实验表明，当测量数( M )较小时（如( M=30 )），FOCUSS的重建误差通常低于OMP。这是因为FOCUSS通过迭代加权更有效地利用了信号的稀疏性，而OMP作为贪心算法可能陷入局部最优。但FOCUSS的计算复杂度高于OMP，适合对重建质量要求较高的场景。

四、优化策略与应用建议

1. 参数选择：

   • ( p )值：通常取( p \in [0.1, 0.5] )，( p )越小稀疏性越强，但收敛速度越慢。
   • ( \epsilon )：建议取( 10^{-6} \sim 10^{-4} )，避免数值不稳定。

2. 加速技巧：

   • 使用Cholesky分解加速矩阵求逆：

     L = np.linalg.cholesky(A_weighted).T
     theta_new = lstsq(L, lstsq(L.T, b_weighted)[0])[0]

   • 对大规模问题，可采用随机化FOCUSS（如Stochastic FOCUSS）。

3. 应用场景：

   • 医学影像：MRI加速成像（通过部分傅里叶采样减少扫描时间）。
   • 无线传感：低功耗数据采集（减少传感器采样率）。
   • 图像处理：超分辨率重建（利用图像在小波域的稀疏性）。

五、总结与展望

FOCUSS算法通过迭代加权最小二乘法，在压缩感知框架下实现了高效的稀疏信号重建。其优势在于对稀疏性的直接优化，但需注意参数调优和计算复杂度。未来研究可结合深度学习（如用神经网络替代迭代步骤）进一步提升重建速度和质量。对于开发者，建议从简单问题（如一维信号）入手，逐步掌握FOCUSS的核心思想后再扩展至复杂应用。