稀疏信号重构（SSR）详解

稀疏信号重构（SSR）详解

目录

1. 引言
2. 稀疏信号及变换域表示
3. 压缩感知与测量模型
4. 稀疏信号重构的优化模型
   1. $L_0$ 范数最小化
   2. $L_1$ 范数最小化：Basis Pursuit (BP)
   3. 带噪声情形与BPDN
   4. Lasso (Least Absolute Shrinkage and Selection Operator)
5. 典型求解算法
   1. 内点法、坐标下降法等凸优化方法
   2. 贪婪算法：OMP为例
   3. 迭代阈值算法 (ISTA/FISTA)
6. 稀疏重构的可行性：RIP与相干性
7. 总结
8. 示例代码

---

1. 引言

在很多应用中（图像处理、医学成像、雷达信号处理、机器学习等），我们经常会遇到以下需求：利用有限观测数据来重构或恢复原始高维信号。
传统采样理论（例如香农采样定理）会要求较高采样率，然而有些情况下，采样成本很高或测量不完整，这就需要一种新的理论与方法——利用信号的稀疏性来进行信号重构，也就是我们所说的稀疏信号重构 (Sparse Signal Recovery, SSR)。

---

2. 稀疏信号及变换域表示

2.1 稀疏性定义

假设有一个$N$维向量
$$\mathbf{x}=\left[\begin{array}{c}{x}_1\\ x_2\\ ⋮\\ x_N\end{array}\right]\in {\mathbb{R}}^N.$$

• 如果仅有 $K\ll N$ 个分量为非零，则称$\mathbf{x}$是$K$-稀疏 (K-sparse)。

在实际应用中，很多信号本身并不一定在时域/空域就稀疏，但往往在某个变换域（如傅里叶变换、小波变换、离散余弦变换等）下可以逼近稀疏。

• 例如，令$\mathbf{\Psi }$为某个$N\times N$的可逆(或正交)变换矩阵，令$\mathbf{\alpha }$为$\mathbf{x}$在该变换域下的系数向量，则
  $$\mathbf{x}=\mathbf{\Psi }\mathbf{\alpha }.$$
• 如果$\mathbf{\alpha }$是$K$-稀疏的，即$\Vert \mathbf{\alpha }{\Vert }_0=K$，就意味着$\mathbf{x}$在$\mathbf{\Psi }$域上是稀疏可表示的。

---

3. 压缩感知与测量模型

3.1 测量模型

一般地，我们进行线性测量，得到一个观测向量 $\mathbf{y}$:
$$\mathbf{y}=A\mathbf{x}+\mathbf{e},$$
其中：

• $A\in {\mathbb{R}}^{M\times N}$ 为观测矩阵，行数 $M$ 通常远小于列数 $N$，即“欠定系统”；
• $\mathbf{y}\in {\mathbb{R}}^M$ 为测量值；
• $\mathbf{e}$ 为噪声或测量误差（在无噪声时可视为 $\mathbf{e}=0$）。

3.2 欠定系统的困难

当 $M<N$ 时，单从$\mathbf{y}=A\mathbf{x}$来解$\mathbf{x}$是不可能唯一的；理论上有无穷多解。但如果事先知道$\mathbf{x}$是稀疏的，我们就可利用这一先验信息构造出一个“稀疏性约束”的优化问题，来找出唯一解或近似唯一解。

3.3 压缩感知的核心

• 信号具有稀疏表示：即$\mathbf{x}$在某个域（或字典）中可以用少量非零系数来表示；
• 随机或良性构造的测量矩阵 $A$：“捕捉”信号主体信息，同时保证可以从有限测量中稳定重构。

---

4. 稀疏信号重构的优化模型

4.1 $L_0$ 范数最小化

为了获得最稀疏解，最直接的想法是最小化 $L_0$ 范数（非零元素个数）：
∥ x ∥ 0 = ∑ i = 1 N 1 { x i ≠ 0 } , \|\mathbf{x}\|_0 = \sum_{i=1}^{N} \mathbf{1}_{\{x_i \neq 0\}}, ∥x∥0​=i=1∑N​1{xi​=0}​,
其中 1 { ⋅ } \mathbf{1}_{\{ \cdot \}} 1{⋅}​ 是指示函数。
于是我们想解以下问题：
$$\underset{\mathbf{x}}{\min }\Vert \mathbf{x}{\Vert }_0\text{subject to}\mathbf{y}=A\mathbf{x}.$$
但是，$\Vert \mathbf{x}{\Vert }_0$ 不仅是非凸函数，而且是NP-困难问题——在大规模应用中几乎不可能直接求解。

---

4.2 $L_1$ 范数最小化：Basis Pursuit (BP)

一个著名的替代策略是：用$L_1$范数来近似$L_0$范数。

• $L_1$ 范数定义：
  $$\Vert \mathbf{x}{\Vert }_1=\sum _{i=1}^N|x_i|.$$
• 由几何直观可知，最小化 $\Vert \mathbf{x}{\Vert }_1$ 也会倾向产生稀疏解，同时它是凸的。

这样，便得到Basis Pursuit (BP)：
$$\underset{\mathbf{x}}{\min }\Vert \mathbf{x}{\Vert }_1\text{subject to}\mathbf{y}=A\mathbf{x}.$$
由于这成为凸优化问题，可以运用很多成熟的凸优化算法（如内点法、线性规划等）来求解。

---

4.3 带噪声情形与BPDN

如果测量包含噪声项 $\mathbf{e}$，或者我们仅仅希望“近似”地满足$\mathbf{y}=A\mathbf{x}$，常见做法是将约束放宽为：
$$\Vert \mathbf{y}-A\mathbf{x}{\Vert }_2\le ϵ,$$
这里 $\Vert \cdot {\Vert }_2$ 是欧几里得 (L2) 范数。这个问题的形式常被称为Basis Pursuit De-Noising (BPDN) 或者噪声下的稀疏重构：
$$\underset{\mathbf{x}}{\min }\Vert \mathbf{x}{\Vert }_1\text{subject to}\Vert \mathbf{y}-A\mathbf{x}{\Vert }_2\le ϵ.$$

---

4.4 Lasso (Least Absolute Shrinkage and Selection Operator)

在统计或机器学习中，经常使用另一种形式：
$$\underset{\mathbf{x}}{\min }\frac{1}{2}\Vert \mathbf{y}-A\mathbf{x}{\Vert }_2^2+\lambda \Vert \mathbf{x}{\Vert }_1,$$

• 其中第一项 $\frac{1}{2}\Vert \mathbf{y}-A\mathbf{x}{\Vert }_2^2$ 用来衡量拟合误差，
• 第二项 $\lambda \Vert \mathbf{x}{\Vert }_1$ 相当于**$L_1$ 正则化项**，鼓励解的稀疏性，
• $\lambda$ 是一个正则化系数，在拟合度与稀疏度之间作平衡。

当 $\lambda$ 越大时，算法更侧重“稀疏”，会更多地将$x_i$压缩为零；当$\lambda$ 越小，则更注重让 $\Vert \mathbf{y}-A\mathbf{x}{\Vert }_2^2$ 变小。

---

5. 典型求解算法

5.1 内点法、坐标下降法等凸优化方法

对 $L_1$ 范数最小化、Lasso 等凸问题，常见的通用算法包括：

• 内点法 (Interior-point method)：可以看作在约束与目标函数的可行域里“内部”迭代寻优；
• 坐标下降法 (Coordinate Descent)：每次在一个坐标方向上寻找最优解，再轮流迭代更新；
• 近端梯度法 (Proximal Gradient Method)：针对带$\Vert \mathbf{x}{\Vert }_1$正则项的优化，有专门的“软阈值”操作来实现稀疏解更新。

这些方法在中等规模上有很好的性能，但对于极大规模问题，仍可能面临较高的计算成本。

---

5.2 贪婪算法：OMP为例

为在实际问题中快速求解，很多贪婪(Greedy)算法被提出。其中最经典的一个是正交匹配追踪(OMP, Orthogonal Matching Pursuit)。

以无噪声简单场景为例，给出OMP算法核心流程：

1. 初始化：
   • 残差向量 ${\mathbf{r}}^{(0)}=\mathbf{y}$。
   • 支持集(已选列集合) $S^{(0)}=\varnothing$。
   • 迭代计数 $t=0$。
2. 迭代：在第$t$次迭代中：
   1. 找到与当前残差相关性最大的列索引：
      $$j^{\star }=\arg \underset{j}{\max }|⟨{\mathbf{a}}_j,{\mathbf{r}}^{(t)}⟩|,$$
      其中${\mathbf{a}}_j$为观测矩阵 $A$ 的第$j$列。
   2. 更新支持集：
      S ( t + 1 ) = S ( t ) ∪ { j ⋆ } . S^{(t+1)} = S^{(t)} \cup \{j^\star\}. S(t+1)=S(t)∪{j⋆}.
   3. 在当前支持集上做最小二乘估计：
      $${\mathbf{x}}^{(t+1)}=\arg \underset{\mathbf{z}}{\min }\Vert \mathbf{y}-A_{S^{(t+1)}}\mathbf{z}{\Vert }_2^2,$$
      其中 $A_{S^{(t+1)}}$ 表示从 $A$ 中选出列索引属于 $S^{(t+1)}$ 的子矩阵。
      该解可显式写为：
      $${\mathbf{z}}^{\ast }=(A_{S^{(t+1)}}^{\mathsf{T}}{A}_{S^{(t+1)}}{)}^{-1}{A}_{S^{(t+1)}}^{\mathsf{T}}\mathbf{y}.$$
      然后把 ${\mathbf{x}}^{(t+1)}$ 的对应分量置为 ${\mathbf{z}}^{\ast }$，其他位置为0。
   4. 更新残差：
      $${\mathbf{r}}^{(t+1)}=\mathbf{y}-A{\mathbf{x}}^{(t+1)}.$$
   5. 迭代次数 $t\leftarrow t+1$，直到达到设定终止条件（如迭代次数 $=K$ 或$\Vert {\mathbf{r}}^{(t)}{\Vert }_2$ 小于一定阈值）为止。

OMP实现速度快、实现简单的优点。但它是贪婪策略，**有时在测量矩阵较“差”**时恢复性能不如$L_1$凸优化方法。

---

5.3 迭代阈值算法 (ISTA/FISTA)

对于Lasso或带$\Vert \mathbf{x}{\Vert }_1$正则项的最小化问题，还可以使用迭代软阈值算法 (ISTA)。
例如，考虑
$$\underset{\mathbf{x}}{\min }\frac{1}{2}\Vert \mathbf{y}-A\mathbf{x}{\Vert }_2^2+\lambda \Vert \mathbf{x}{\Vert }_1.$$
每次迭代大致做以下更新：
$${\mathbf{x}}^{(k+1)}={\mathrm{soft}}_{\lambda \cdot \eta }({\mathbf{x}}^{(k)}+\eta A^{\mathsf{T}}(\mathbf{y}-A{\mathbf{x}}^{(k)})),$$

• $\eta$ 是一个步长（需要合适选取，一般与$\Vert A{\Vert }_2^2$有关），
• ${\mathrm{soft}}_{\theta }(\cdot )$ 是软阈值操作：
  s o f t θ ( z ) = s i g n ( z )   max ⁡ { ∣ z ∣ − θ ,   0 } . \mathrm{soft}_{\theta}(z) = \mathrm{sign}(z)\,\max\{|z| - \theta,\,0\}. softθ​(z)=sign(z)max{∣z∣−θ,0}.

若使用Nesterov加速技术，则是FISTA (Fast ISTA)，在很多问题中收敛更快。

---

6. 稀疏重构的可行性：RIP与相干性

6.1 RIP (Restricted Isometry Property)

为什么用$L_1$范数就能重构出原信号？这离不开对观测矩阵 $A$ 的一定要求。

• 限制等距性质 (RIP)：若对所有$K$-稀疏向量 $\mathbf{x}$，有
  $$(1-{\delta }_K)\Vert \mathbf{x}{\Vert }_2^2\le \Vert A\mathbf{x}{\Vert }_2^2\le (1+{\delta }_K)\Vert \mathbf{x}{\Vert }_2^2,$$
  则称 $A$ 满足RIP常数${\delta }_K$。当${\delta }_K$足够小（小于某个特定阈值），可保证**$L_1$最小化**精确重构$K$-稀疏向量。

6.2 相干性 (Coherence)

另一种衡量矩阵“好坏”的简单指标是相干性(coherence)：
$$\mu (A)=\underset{i\ne j}{\max }\frac{|⟨{\mathbf{a}}_i,{\mathbf{a}}_j⟩|}{\Vert {\mathbf{a}}_i{\Vert }_2\Vert {\mathbf{a}}_j{\Vert }_2},$$

• ${\mathbf{a}}_i$ 表示矩阵 $A$ 的第$i$列。
• 如果相干性 $\mu (A)$ 很小，则列之间不太“相似”，有助于稀疏重构的成功。

理论上，随机高斯矩阵或随机部分傅里叶矩阵往往满足良好的RIP或相干性小，从而保证重构效果。

---

7. 总结

• 稀疏信号重构 (SSR) 依赖于信号稀疏性这一重要先验；
• 通过 $L_0$ 范数最小化 可以直接获得最稀疏解，但因其NP-困难不可直接求解；
• $L_1$ 范数最小化（Basis Pursuit / Lasso）是业界最常用的替代方案，具有良好的理论保证（RIP、相干性等）；
• 实际应用中，既有凸优化求解 (内点法、坐标下降、ISTA/FISTA) 也有贪婪算法 (OMP、CoSaMP等)；
• SSR极大拓展了传统采样理论的边界，已经成为压缩感知、高维数据分析、机器学习稀疏建模等领域的核心技术。

---

8. 示例代码

下面给出一个简化的OMP示例伪代码（仅作演示）：

def OMP(y, A, K):
    """
    y: 测量向量, shape = (M,)
    A: 观测矩阵, shape = (M, N)
    K: 稀疏度（最多选K个原子/列）
    """
    import numpy as np
    
    M, N = A.shape
    
    # 初始化
    r = y.copy()              # 残差
    S = []                    # 当前支持集，初始为空
    x_hat = np.zeros(N)       # 最终重构向量
    
    for k in range(K):
        # 1. 找到与当前残差最相关的列
        corr = np.abs(A.T @ r)   # 相关性
        j_star = np.argmax(corr)
        
        # 2. 将该列放入支持集
        if j_star not in S:
            S.append(j_star)
        
        # 3. 在当前S上做最小二乘求解
        A_S = A[:, S]  # 仅取支持集对应列
        x_S = np.linalg.lstsq(A_S, y, rcond=None)[0]  # 求解最小二乘
        
        # 4. 更新残差
        r = y - A_S @ x_S
        
        # 可加终止条件：若残差足够小则提前停止
        
    # 将计算得到的x_S写回x_hat
    x_hat[S] = x_S
    return x_hat