线性代数（十二）稀疏矩阵与压缩感知

线性代数的本质（十二）：稀疏矩阵与压缩感知

稀疏性是自然界和数据中普遍存在的特征，压缩感知利用这种特性实现了"少即是多"的奇迹。

引言：从 JPEG 到 MRI ——稀疏性无处不在

你是否想过，为什么一张 10MB 的原始照片可以压缩成几百 KB 的 JPEG 文件而几乎看不出差别？为什么 MRI 扫描时间可以从 30 分钟缩短到 5 分钟？

答案就是稀疏性——信号在某个合适的表示下，大部分系数为零或接近零。

生活中的稀疏性： - 图像：自然图像在小波基下是稀疏的（大部分小波系数很小） - 音频：语音信号在傅里叶基下是稀疏的（只有少数频率有能量） - 文本：一篇文章只使用词汇表中很小一部分词汇 - 社交网络：一个人只认识所有人中很小一部分

压缩感知的革命性思想：既然信号是稀疏的，我们能否用比传统方法更少的测量来恢复它？

答案是肯定的！这就是本章要探讨的核心内容。

本章学习目标

1. 稀疏表示：理解信号在字典下的稀疏分解
2. L0 与 L1 范数：为什么 L1 能促进稀疏？
3. 压缩感知理论： RIP 条件、测量矩阵设计
4. 算法实现： LASSO 、 OMP 、 ISTA
5. 应用： MRI 成像、单像素相机、图像压缩

一、稀疏表示的数学基础

1.1 什么是稀疏性？

定义：向量 是k-稀疏的（ k-sparse），如果它最多有 个非零元素。

数学表示：

这里 称为"L0 范数"（虽然严格来说它不是真正的范数，因为不满足齐次性）。

例子：

• 是 3-稀疏的
• 严格来说是 2-稀疏的

近似稀疏：实际中，信号通常是近似稀疏的——大部分系数很小但不严格为零。

比喻：想象一个班级的成绩分布。如果只有少数学生得了高分，大部分人是零分或低分，这个成绩向量就是稀疏的。

1.2 稀疏表示与字典

大多数信号本身并不稀疏，但在某个变换域或字典下可能是稀疏的。

数学表达：给定信号 和字典 ，找稀疏系数 ：

字典的列称为原子（ atoms）。

常见字典：

字典类型	适用信号	特点
傅里叶基	周期信号、音频	频域稀疏
小波基	图像、非平稳信号	时频局部化
DCT 基	图像（ JPEG）	类傅里叶，实数
学习字典	特定领域数据	数据驱动

过完备字典：当 时，字典是过完备的，意味着有更多的原子可以选择，表示更灵活，但也带来了唯一性问题。

1.3 稀疏编码问题

问题：给定信号 和字典 ，找最稀疏的表示：

困难：这是一个NP-hard问题！

为什么 NP-hard？ - 需要搜索所有可能的非零元素位置组合 - 对于 维向量，有 种可能的稀疏模式 - 没有已知的多项式时间算法

比喻：这就像在图书馆找一本书，但你只知道书的内容大意，而且书可能放在任何位置。你需要检查每一个位置组合，这在书很多时是不可行的。

1.4 稀疏表示的唯一性

问题：即使找到了稀疏解，它是唯一的吗？

定理（唯一性条件）：如果 且，则 是唯一的最稀疏表示。

其中 是字典的互相干性（ mutual coherence）：

直觉：互相干性衡量字典中原子之间的最大相似度。如果原子之间非常不同（低互相干性），稀疏解更容易唯一。

二、 L1 正则化：稀疏性的凸松弛

2.1 从 L0 到 L1

既然 L0 优化是 NP-hard 的，我们需要一个可计算的替代。

关键观察： L1 范数是 L0 范数的最佳凸松弛。

凸松弛问题：

这是一个线性规划问题，可以在多项式时间内求解！

2.2 为什么 L1 促进稀疏？——几何解释

这是本章最重要的直觉之一。

L1 球的形状：在 2D 中， L1 单位球是菱形（正方形旋转 45 度）：

L2 球的形状：在 2D 中， L2 单位球是圆：

关键区别： L1 球在坐标轴上有"尖角"（ corner），而 L2 球是光滑的。

几何直觉： 1. 约束条件 定义了一个仿射子空间（直线、平面等） 2. 最小化 等价于找这个子空间与最小 球的交点 3. L1 球的"尖角"在坐标轴上，所以交点更可能在角点（稀疏解） 4. L2 球光滑，交点一般不在坐标轴上（非稀疏解）

比喻：想象你在一个房间里，需要触碰墙壁（约束）。如果你手里拿着一个菱形气球（ L1 球）慢慢膨胀，气球最先碰到墙的地方很可能是尖角。但如果拿着圆形气球（ L2 球），碰到墙的点可能在任何地方。

2.3 L1 与 L2 正则化的比较

特性	L1 正则化	L2 正则化
单位球形状	菱形（有尖角）	球形（光滑）
解的特点	稀疏（很多精确零）	小但非零
几何直觉	倾向于在坐标轴上	倾向于均匀缩小
优化	非光滑，需要特殊算法	光滑，梯度下降可用
应用	特征选择、压缩感知	防过拟合、岭回归

2.4 Elastic Net：结合 L1 和 L2

弹性网络结合了两者的优点：

优点： - 继承 L1 的稀疏性 - L2 部分提供稳定性，特别是当特征高度相关时 - 可以选择特征组（ grouping effect）

三、 LASSO 回归

3.1 LASSO 的定义

LASSO（ Least Absolute Shrinkage and Selection Operator）是统计学中最重要的稀疏方法：

等价形式（约束形式）：

3.2 LASSO 的解

子梯度条件： LASSO 的最优解 满足： $$

X^T(X - ) + = 0 $$

其中 是 的子梯度： $$

s_j =

$$

3.3 软阈值操作

对于简单情况（正交设计）， LASSO 有闭式解：

这称为软阈值（ soft thresholding）操作。

对比硬阈值： - 软阈值：平滑地将小系数缩小到零 - 硬阈值：直接将小系数设为零，大系数不变

可视化：软阈值函数像一个"挤压"操作，将所有系数向零推，小系数完全变成零。

3.4 LASSO 的性质

稀疏性： LASSO 自动将不重要的特征系数设为精确的零。

特征选择：非零系数对应被选中的特征，这是自动的特征选择。

偏差-方差权衡： - 大 → 更稀疏，偏差大，方差小 - 小 → 更稠密，偏差小，方差大

选择：通常使用交叉验证。

3.5 LASSO 路径

LASSO 路径是随 变化的轨迹。

特点： - 当 很大时，所有系数为零 - 随着 减小，系数逐渐变为非零 - 路径是分段线性的（可以高效计算）

LARS 算法（ Least Angle Regression）可以高效计算整个 LASSO 路径。

四、压缩感知理论

4.1 核心思想

传统采样：香农-奈奎斯特定理要求采样率至少是信号最高频率的两倍。

压缩感知：如果信号是稀疏的，可以用远少于奈奎斯特率的测量来精确恢复信号！

革命性含义： - MRI 扫描时间可以大幅缩短 - 传感器可以更简单、更便宜 - 数据存储和传输更高效

4.2 测量模型

其中： - ：原始信号（假设 k-稀疏或在某基下 k-稀疏） - ：测量矩阵（） - ：测量值 - ：噪声

问题：从 个测量恢复 维信号（欠定系统，无穷多解！）

关键洞察：稀疏性提供了额外约束，使恢复成为可能。

4.3 限制等距性质（ RIP）

定义：矩阵 满足k 阶限制等距性质（ RIP），如果存在 使得对所有 k-稀疏向量：

直觉： 作用在稀疏向量上时，近似保持向量的长度。

几何意义： 将稀疏向量从 映射到 时，不会过度压缩或拉伸。

比喻：想象把一个高维物体投影到低维。 RIP 保证投影不会把不同的稀疏向量映射到同一点（否则无法区分）。

4.4 恢复定理

定理（ Cand è s-Tao）：如果 满足，则从 可以通过 L1 最小化精确恢复任何 k-稀疏信号：

测量数要求：满足 RIP 的随机矩阵存在，只需要 个测量。

意义： - 传统方法需要 个测量 - 压缩感知只需要 个测量 - 当 时，这是巨大的节省！

4.5 测量矩阵的设计

随机矩阵（最常用）： - 高斯随机矩阵： - 伯努利随机矩阵： - 部分傅里叶矩阵：随机选择傅里叶矩阵的行

随机矩阵的优点： - 以高概率满足 RIP - 与信号无关（通用性） - 容易生成

确定性矩阵：存在但构造复杂，常用于特定应用。

4.6 噪声情况

Basis Pursuit Denoising（ BPDN）：

恢复保证：如果 满足 RIP 且噪声有界， BPDN 的解与真实信号的误差是量级。

五、算法

5.1 迭代收缩阈值算法（ ISTA）

ISTA 是求解 LASSO 问题的经典算法。

更新规则：

其中： - 是 的最大特征值（ Lipschitz 常数） - 是软阈值函数

算法解释：

1. 梯度下降步：
2. 软阈值步：促进稀疏性

收敛速度：

5.2 快速迭代收缩阈值算法（ FISTA）

FISTA 使用 Nesterov 加速技术：

收敛速度：——平方级加速！

5.3 正交匹配追踪（ OMP）

OMP 是一种贪心算法：

算法步骤：

1. 初始化：残差，支撑集
2. 重复直到满足停止条件：
   • 选择：
   • 更新支撑集：
   • 最小二乘更新：
   • 更新残差：

优点：简单、快速

缺点：贪心可能找不到全局最优

5.4 坐标下降法

思想：每次只优化一个坐标。

对于 LASSO：

优点：对于大规模稀疏问题非常高效。

六、应用

6.1 医学成像：压缩感知 MRI

背景： MRI 扫描需要采集大量 k 空间数据，耗时长。

压缩感知 MRI： - 只采集部分 k 空间数据（随机欠采样） - 利用图像在小波域的稀疏性 - 通过 L1 优化重建图像

效果：扫描时间缩短 2-8 倍，对运动敏感的成像特别有用。

实际系统： 2017 年起， FDA 批准了多款压缩感知 MRI 系统。

6.2 单像素相机

原理：用单个光电探测器拍摄图像。

工作方式： 1. 用数字微镜阵列（ DMD）生成随机图案 2. 每个图案产生一个测量值（所有像素的加权和） 3. 个测量足以重建 像素图像（）

应用： - 红外成像（探测器昂贵） - 太赫兹成像 - 高速成像

6.3 图像压缩与恢复

图像修复： - 已知部分像素，恢复完整图像 - 利用图像的稀疏先验（小波、 DCT 、学习字典）

超分辨率： - 从低分辨率图像恢复高分辨率 - 假设高分辨率图像是稀疏可表示的

6.4 基因组学

问题：识别与疾病相关的基因。

数据： 个基因， 个样本（）

方法： LASSO 回归 - 稀疏性假设：只有少数基因与疾病相关 - 自动选择相关基因

6.5 金融

稀疏投资组合： - 只投资少数资产（降低交易成本） - L1 正则化促进稀疏持仓

高频交易信号： - 从大量特征中选择关键指标

七、稀疏矩阵存储与计算

7.1 稀疏矩阵格式

COO 格式（ Coordinate）： 存储三元组 表示非零元素。

CSR 格式（ Compressed Sparse Row）： - data：非零元素值 - indices：列索引 - indptr：行指针

CSC 格式（ Compressed Sparse Column）： 类似 CSR，但按列存储。

7.2 稀疏矩阵运算

矩阵-向量乘法： - 稠密： - 稀疏（ nnz 个非零元）： 稀疏矩阵求解： - 直接方法：稀疏 LU 分解 - 迭代方法：共轭梯度、 GMRES

八、深入理解：为什么压缩感知有效？

8.1 高维几何的反直觉

在高维空间中，我们的直觉经常失效：

高维球面的体积： - 大部分体积集中在球面附近（不在中心） - 高维球的体积趋向于零

稀疏向量的结构： - k-稀疏向量生活在 个 k 维子空间的并集中 - 这个并集的"维度"远小于

8.2 约翰逊-林登斯特劳斯引理

JL 引理： 个点可以被随机投影到 维空间，同时保持两两距离（相对误差）。

与压缩感知的联系： RIP 是 JL 引理在稀疏向量上的加强版。

8.3 信息论视角

稀疏信号的信息量： k-稀疏向量的自由度约为。

测量数的下界： 是必要的。

结论：压缩感知的测量数是信息论最优的（ up to 常数）！

九、 Python 实现

9.1 LASSO 回归

import numpy as np
from sklearn.linear_model import Lasso

# 生成稀疏数据
np.random.seed(42)
n, p, k = 100, 500, 10  # 样本数、特征数、稀疏度

# 真实稀疏系数
beta_true = np.zeros(p)
beta_true[:k] = np.random.randn(k) * 5

# 设计矩阵和响应
X = np.random.randn(n, p) / np.sqrt(n)
y = X @ beta_true + np.random.randn(n) * 0.5

# LASSO 回归
lasso = Lasso(alpha=0.1)
lasso.fit(X, y)

# 比较
print(f"真实非零系数数: {np.sum(beta_true != 0)}")
print(f"估计非零系数数: {np.sum(np.abs(lasso.coef_) > 0.01)}")
print(f"恢复误差: {np.linalg.norm(lasso.coef_ - beta_true):.4f}")

9.2 OMP 算法

def omp(Phi, y, k):
    """正交匹配追踪算法"""
    n = Phi.shape[1]
    r = y.copy()  # 残差
    S = []  # 支撑集
    x = np.zeros(n)
    
    for _ in range(k):
        # 选择最相关的原子
        correlations = np.abs(Phi.T @ r)
        j = np.argmax(correlations)
        S.append(j)
        
        # 最小二乘更新
        Phi_S = Phi[:, S]
        x_S = np.linalg.lstsq(Phi_S, y, rcond=None)[0]
        
        # 更新残差
        r = y - Phi_S @ x_S
    
    # 构造稀疏解
    x[S] = x_S
    return x

9.3 ISTA 算法

def soft_threshold(x, tau):
    """软阈值函数"""
    return np.sign(x) * np.maximum(np.abs(x) - tau, 0)

def ista(Phi, y, lambda_param, max_iter=1000, tol=1e-6):
    """迭代收缩阈值算法"""
    m, n = Phi.shape
    x = np.zeros(n)
    
    # Lipschitz 常数
    L = np.linalg.norm(Phi.T @ Phi, 2)
    
    for i in range(max_iter):
        x_old = x.copy()
        
        # 梯度下降步
        grad = Phi.T @ (Phi @ x - y)
        x = x - (1/L) * grad
        
        # 软阈值步
        x = soft_threshold(x, lambda_param / L)
        
        # 收敛检查
        if np.linalg.norm(x - x_old) < tol:
            break
    
    return x

十、总结与展望

本章关键要点

1. 稀疏性是自然信号的普遍特征

2. L1 范数是 L0 的凸松弛，几何上促进稀疏

3. LASSO同时实现稀疏估计和特征选择

4. 压缩感知：从 测量恢复 k-稀疏信号

5. RIP 条件保证恢复的正确性

6. 算法： ISTA 、 FISTA 、 OMP 、坐标下降

与其他章节的联系

• 第 9 章（ SVD）：低秩是另一种"稀疏性"
• 第 10 章（范数）： L1/L2 范数的性质
• 第 11 章（优化）： LASSO 是凸优化问题
• 第 15 章（机器学习）： LASSO 用于特征选择

进一步学习

1. 结构稀疏：群稀疏、树稀疏
2. 字典学习：数据驱动的字典设计
3. 矩阵补全：低秩矩阵的恢复
4. 深度压缩感知：神经网络+稀疏恢复

下一章预告

《张量与多线性代数》

• 张量是向量和矩阵的推广
• CP 分解、 Tucker 分解
• 应用：推荐系统、化学计量学

练习题

基础题

1. 证明 L1 球在坐标轴上有角点。

2. 计算向量 的 L0 、 L1 、 L2 范数。

3. 解释为什么 L0 优化是 NP-hard 的。

进阶题

4. 推导软阈值操作的公式。

5. 证明：如果 满足 RIP，则不同的稀疏向量被映射到不同的测量。

6. 实现 FISTA 算法并与 ISTA 比较收敛速度。

编程题

7. 在 MNIST 数据集上用 LASSO 进行特征选择。

8. 实现压缩感知图像重建。

9. 比较 OMP 和 LASSO 的恢复性能。

参考资料

1. Cand è s, E. & Wakin, M. (2008). "An Introduction to Compressive Sampling". IEEE Signal Processing Magazine.

2. Tibshirani, R. (1996). "Regression Shrinkage and Selection via the Lasso". JRSS-B.

3. Elad, M. (2010). Sparse and Redundant Representations. Springer.

4. Hastie, T., Tibshirani, R., & Wainwright, M. (2015). Statistical Learning with Sparsity. CRC Press.

本文是《线性代数的本质与应用》系列的第 12 章，共 18 章。 作者：[Your Name] | 发布日期： 2019-05-20