线性代数（十二）：稀疏矩阵与压缩感知——少即是多的数学奇迹

「少即是多」的奇迹#

一张 2400 万像素的原始照片，大小约为 70 MB；用 JPEG 压缩后，只有几百 KB，压缩比高达 100 倍，但你看不出区别。传统 MRI 扫描需要 30 分钟，而现代压缩感知 MRI 只需 5 分钟，图像质量完全一样。

这两个奇迹的背后，靠的是同一个核心：稀疏性。大多数自然信号，只要用合适的基表示，真正有意义的系数就那么几个，其余几乎全是零。

压缩感知把这一特性变成了算法。由于信号本身稀疏，测量次数可以远低于经典采样定理的要求，但仍能精确恢复信号。本章将讲清楚背后的原理。

本章你会学到#

• 稀疏性的数学定义，以及它为什么无处不在
• $L_0$ /$L_1$ /$L_2$ 范数，以及 $L_1$ 产生稀疏解的几何原因
• LASSO：一次优化完成回归和特征选择
• 压缩感知理论：RIP 条件、测量矩阵、恢复保证
• 算法：ISTA、FISTA、迭代硬阈值（IHT）、OMP
• 实际应用：MRI、单像素相机、基因组学、金融

预备知识#

• 范数与条件数（第 10 章 ）
• 凸优化（第 11 章 ）
• SVD 与低秩逼近（第 9 章 ）

稀疏性：无处不在的模式#

定义#

这里的 $L_0$ “范数”其实就是统计非零元素的数量。严格来说，它并不是真正的范数——因为它不满足齐次性（比如 $\|2x\|_0 = \|x\|_0$ ）——但它是最直观的稀疏性度量。

实际中的信号很少是完全稀疏的，更多时候是可压缩的。把信号的系数按大小排序后，序列会迅速衰减。只保留前 $k$ 个最大的系数，就能得到一个几乎完美的近似。JPEG 就是这么工作的：计算 DCT，保留几百个最大的系数，其余的直接扔掉。

现实中的稀疏性#

字典下的稀疏表示#

常见字典如下：

当 $p > n$ 时，字典是过完备的：信号可以有多种表示方式，而我们需要找到其中最稀疏的那个。

稀疏矩阵的存储#

在讨论恢复之前，先想想一个大部分元素为零的矩阵该怎么存。三种主流格式各有优劣，分别在构建、读取和矩阵-向量乘法上做了不同的权衡。

稀疏矩阵与向量相乘的复杂度是 $O(\text{nnz})$ ，而不是稠密矩阵的 $O(mn)$ 。对于一个 $10^6 \times 10^6$ 且包含 100 万个非零元素的矩阵，这一步就能提速百万倍。

图中展示了一个 $8 \times 8$ 的矩阵，包含 13 个非零元素。COO 直接列出每个非零元素的行、列和值。CSR 将行索引压缩成指针：indptr[i]:indptr[i+1] 表示第 $i$ 行在 cols 和 vals 中的切片范围。CSC 对列做类似处理。

什么时候真能省内存？CSR 每个非零元素需要 8 字节存值，4 字节存列索引，总共 12 字节；而稠密存储每个元素固定占 8 字节，无论是否为零。

CSR 在约 67% 密度（8 / (8+4)）时与稠密存储持平。密度低于 10% 时，大多数实际稀疏矩阵都处于这个区间，内存和计算节省能达到数量级。

$L_1$ 正则化：让稀疏可计算的关键#

$L_0$ 不可解#

这是个 NP 难 问题。最直接的方法是枚举所有 $\binom{p}{k}$ 种可能的支撑组合，结果就是组合爆炸。目前没有已知的多项式时间算法能解决一般情况。

$L_1$ 作为凸松弛#

这是个 线性规划 问题，能在多项式时间内求解。神奇的是，在对 $D$ 施加一些温和假设后，$L_1$ 的解与 $L_0$ 的解完全一致。凸松弛在这里是 精确 的。

为什么 $L_1$ 能促进稀疏？几何视角#

把约束 $\{\vec{\alpha} : D\vec{\alpha} = \vec{x}\}$ 看作 $\mathbb{R}^p$ 中的一个仿射子空间。最小化 $\|\vec{\alpha}\|_1$ 相当于从原点向外膨胀一个 $L_1$ 球，直到它刚好碰到约束面。

• $L_1$ 球 是个菱形，尖角正好落在坐标轴上。
• $L_2$ 球 是个圆球，表面处处光滑。

菱形膨胀时，通常会先用尖角碰到一条普通直线——而尖角在坐标轴上，意味着某些坐标会被强制为零。圆球碰到的则是表面上的一个普通点，没有任何坐标被强制为零。

这张图就是全部直觉所在——$L_1$ 、LASSO、基追踪、压缩感知背后的原理都在这里。高维版本只是同一个故事的延续：$\mathbb{R}^n$ 中的 $L_1$ 球在坐标轴上有 $2n$ 个尖角，“先碰到尖角”的几何特性依然成立。

$L_1$ 与 $L_2$ 对比#

Elastic Net#

既保留了 $L_1$ 的稀疏性，又借 $L_2$ 获得了稳定性。这就是 Elastic Net，在高维回归中更稳妥的默认选择。

LASSO 回归#

定义#

这个公式同时做了两件事：用数据项拟合 $\vec{y}$ ，用 $L_1$ 惩罚把不重要的系数直接压到 零。所以 LASSO 一次搞定回归和特征选择。

软阈值#

这就是 软阈值 操作：每个 OLS 系数都往零方向压缩 $\lambda$ ，绝对值小于 $\lambda$ 的直接变成零。对比 硬阈值，硬阈值会保留大系数不变，小系数直接归零，显得突兀且不稳定。

LASSO 路径#

让 $\lambda$ 从 $\infty$ 逐渐减小到 $0$ ：

• $\lambda$ 很大时，所有系数都是零。
• $\lambda$ 减小时，特征一个接一个进入活动集，每个特征在残差相关性超过阈值时加入。
• $\lambda \to 0$ 时，退化为普通最小二乘法。

Efron、Hastie、Johnstone 和 Tibshirani（2004）证明了一个漂亮的结论：路径关于 $\lambda$ 是 分段线性的。LARS 算法利用这一点，用一次普通最小二乘法的计算成本就能算出整条路径。

真正活跃的特征（粗线）很早就进入模型，系数也迅速增大；伪特征要么始终不出现，要么很晚才进来，系数也很小。我们通常用交叉验证选 $\lambda$ ，并采用“1-SE 法则”：在最小误差 $\lambda$ 附近一个标准误差范围内，挑最大的 $\lambda$ ，这样更倾向于稀疏解。

压缩感知#

革命性的想法#

香农-奈奎斯特 告诉我们：要恢复一个带限信号，采样率必须达到其最高频率的两倍。这是铁律。

压缩感知 则说：如果信号在某个基下是 $k$ -稀疏的，只需要 $m \sim k\,\log\!\frac{n}{k}$ 次线性测量就够了。比如，一个 256 像素的信号，只有 10 个非零小波系数，那么 $k\log(n/k) \approx 32$ 次测量就足够了——不需要 256 次。

这可不是天上掉馅饼。它之所以成立，是因为信号模型更强了。香农-奈奎斯特假设“任何带限信号”，而压缩感知假设“任何 $k$ -稀疏信号”。更窄的模型让采样变得更便宜。

测量模型#

其中 $\vec{x} \in \mathbb{R}^n$ 是 $k$ -稀疏的，$\Phi \in \mathbb{R}^{m \times n}$ 且 $m \ll n$ ，$\vec{e}$ 是噪声。系统 $\Phi \vec{x} = \vec{y}$ 是欠定的，解有无穷多个。稀疏性就是那个挑选正确解的约束条件。

图中用 ISTA 算法测试了随机高斯矩阵 $\Phi$ ，参数为 $m=64$ 、$n=256$ 、$k=10$ 。恢复信号与真实信号的相对误差约为 $10^{-2}$ ，支撑集和幅值完全一致。回想一下，$k\log(n/k) \approx 32$ ，64 次测量远超阈值，因此恢复效果几乎完美。

限制等距性质（RIP）#

什么时候欠定的 $\Phi$ 能保留足够的信息，唯一确定一个稀疏的 $\vec{x}$ ？Candes 和 Tao 的答案是 RIP。

简单来说，$\Phi$ 在稀疏向量上几乎像等距映射一样工作——不会把两个不同的稀疏向量挤到同一个像。

左图：对大量随机 $k$ -稀疏向量计算 $\|\Phi\vec{x}\|/\|\vec{x}\|$ 的直方图，不同稀疏度下的分布都紧密集中在 1 附近——范数几乎被完美保留。右图：经验最坏情况畸变 $\delta_k$ 随 $k$ 增长的变化。超过某个稀疏度后，$\sqrt{2}-1$ 的界限会被打破；低于这个值时，恢复有保证。

恢复定理#

恢复误差被控制在常数倍的 $\epsilon$ 内，随噪声温和退化。

如何设计 $\Phi$ #

构造满足 RIP 的确定性矩阵很难。诀窍是抛硬币：随机矩阵以压倒性概率满足 RIP。

通用性——同一个 $\Phi$ 对任何 $k$ -稀疏信号都适用——是随机测量矩阵在工程中无可替代的原因。

算法#

ISTA：迭代收缩-阈值算法#

其中 $L = \|\Phi^\top\Phi\|_2$ 是光滑部分的 Lipschitz 常数。每次迭代先做一步梯度下降，再进行软阈值操作。收敛速度为 $O(1/t)$ 。

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15

import numpy as np

def soft_threshold(x, tau):
    return np.sign(x) * np.maximum(np.abs(x) - tau, 0.0)

def ista(Phi, y, lam, max_iter=2000, tol=1e-6):
    n = Phi.shape[1]
    x = np.zeros(n)
    L = float(np.linalg.norm(Phi.T @ Phi, 2))
    for _ in range(max_iter):
        x_old = x.copy()
        x = soft_threshold(x - (Phi.T @ (Phi @ x - y)) / L, lam / L)
        if np.linalg.norm(x - x_old) < tol:
            break
    return x

FISTA：Nesterov 加速版 ISTA#

在 ISTA 的基础上加入动量项，收敛速度从 $O(1/t)$ 提升到 $O(1/t^2)$ ：

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15

def fista(Phi, y, lam, max_iter=2000, tol=1e-6):
    n = Phi.shape[1]
    x = np.zeros(n)
    z = x.copy()
    t = 1.0
    L = float(np.linalg.norm(Phi.T @ Phi, 2))
    for _ in range(max_iter):
        x_old = x.copy()
        x = soft_threshold(z - (Phi.T @ (Phi @ z - y)) / L, lam / L)
        t_new = (1.0 + np.sqrt(1.0 + 4.0 * t * t)) / 2.0
        z = x + ((t - 1.0) / t_new) * (x - x_old)
        t = t_new
        if np.linalg.norm(x - x_old) < tol:
            break
    return x

IHT：迭代硬阈值算法#

IHT 直接求解 $L_0$ 问题，在类 RIP 条件下能收敛到 $k$ -稀疏解。每次迭代只需一次矩阵-向量乘法和一次部分排序，计算成本极低。

四张子图分别展示了第 0、2、5、30 次迭代的结果（淡色背景是真实信号）。两次迭代后支撑集基本正确，30 次迭代后幅值也接近真实值。残差范数 $\|\Phi x_t - y\|_2$ 下降了多个数量级。

OMP：正交匹配追踪算法#

OMP 是一种纯贪心算法：每一步选择与当前残差最相关的原子，加入支撑集，然后在支撑集上重新求解最小二乘问题。

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12

def omp(Phi, y, k):
    n = Phi.shape[1]
    r = y.copy()
    S = []
    x = np.zeros(n)
    for _ in range(k):
        S.append(int(np.argmax(np.abs(Phi.T @ r))))
        Phi_S = Phi[:, S]
        x_S, *_ = np.linalg.lstsq(Phi_S, y, rcond=None)
        r = y - Phi_S @ x_S
    x[S] = x_S
    return x

OMP 简单易懂、速度快，但贪心策略容易导致早期错误无法纠正。相比之下，LASSO 和 FISTA 虽然较慢，但能保证全局最优解。

应用#

压缩感知 MRI#

MRI 在 k 空间 采样，也就是图像的傅里叶变换，每次采一行。全采 k 空间是 MRI 慢的主要原因。

CS-MRI 的做法：

1. 只随机采 k 空间的一部分行（通常 $m \approx 0.25 n$ ）。
2. 利用 MR 图像在小波域稀疏的特点。
3. 解优化问题 $\min \|\Psi \vec{x}\|_1$ ，约束条件是 $\Phi \vec{x} \approx \vec{y}$ ，其中 $\Psi$ 是小波变换。

实际效果：扫描时间缩短 2 到 8 倍，诊断质量不变。FDA 从 2017 年开始批准了多款 CS-MRI 产品，比如 Siemens 的 “Compressed Sensing Cardiac Cine” 和 GE 的 “HyperSense”。

单像素相机#

直接省掉百万像素传感器。数字微镜阵列把一系列随机 $\pm 1$ 图案投到场景上，单个光电二极管将每张图案的光强积分成一个数值。用 $m \sim k \log(n/k)$ 个图案就能重建出 $n$ 像素的图像。在红外和太赫兹波段，像素阵列贵得离谱，这种技术不可或缺。

基因组学#

基因研究通常涉及 $n \sim 20{,}000$ 个候选基因和 $m \sim 100$ 个患者——典型的 $n \gg m$ 场景，OLS 根本没法用。LASSO 假设只有少数基因对性状起作用，并自动筛选出这些基因。同样的方法还能用于数量性状定位、GWAS 后处理以及药物响应预测。

稀疏投资组合#

经典的 Markowitz 投资组合是稠密的——你得持有每种资产的一小部分。在均值-方差目标函数中加入 $\lambda \|\vec{w}\|_1$ ，得到的投资组合通常只有几十个有效仓位。交易成本大幅降低，风险调整后的收益几乎不受影响。

为什么压缩感知能成立？再深一层#

高维几何帮了大忙#

一个 $k$ -稀疏向量位于 $\binom{n}{k}$ 个 $k$ 维坐标子空间的并集中——这在 $\mathbb{R}^n$ 中占比极小，几乎可以忽略不计。稀疏向量集合非常稀疏，随机低维投影几乎总能把任意两个稀疏向量分开。这就是 RIP 的几何核心。

与 Johnson-Lindenstrauss 引理的联系#

JL 引理 指出：可以把 $N$ 个点投影到 $O(\log N / \epsilon^2)$ 维空间，同时将所有 $\binom{N}{2}$ 对点之间的距离误差控制在 $1\pm\epsilon$ 范围内。RIP 是 JL 引理针对（不可数的）$k$ -稀疏向量集合的特例。两者的证明都依赖于相同的测度集中工具。

信息论的最优性#

$\mathbb{R}^n$ 中的 $k$ -稀疏向量包含大约 $k \log(n/k)$ 比特的结构信息（哪 $k$ 个位置非零），再加上 $k$ 个实数值。压缩感知用 $O(k \log(n/k))$ 次测量就能完成恢复——在常数倍范围内达到了信息论下界。本质上，不可能做得更好。

练习题#

基础题#

1. 给定 $\vec{x} = (0, 3, -1, 0, 0, 2, 0)$ ，计算 $\|\vec{x}\|_0$ 、$\|\vec{x}\|_1$ 和 $\|\vec{x}\|_2$ 。
2. 用一段话说明为什么 $L_0$ 最小化是 NP 难问题，以及为什么凸松弛能起作用。
3. 最小化 $\tfrac{1}{2}(z-a)^2 + \lambda|z|$ ，推导出软阈值公式。

进阶题#

4. 证明 $\mathbb{R}^n$ 中的 $L_1$ 单位球恰好有 $2n$ 个顶点，并且每个顶点都位于某条坐标轴上。
5. 证明：如果 $\Phi$ 满足 $\delta_{2k} < 1$ ，那么任意两个不同的 $k$ -稀疏向量会被映射到不同的测量结果。
6. 推导 LASSO 的坐标下降更新公式，并解释为什么单坐标的优化问题会归结为软阈值操作。
7. 证明 LASSO 路径关于 $\lambda$ 是分段线性的。

编程题#

8. 在同一个压缩感知问题上实现 ISTA 和 FISTA。在对数坐标下绘制目标函数随迭代次数的变化曲线，观察 $O(1/t)$ 和 $O(1/t^2)$ 的收敛速度差异。
9. 对比 OMP、IHT 和 LASSO 在合成数据上的恢复成功率，$k$ 从 1 到 $m/2$ 变化，绘制“成功率 vs $k$ ”的曲线。
10. 取一张 $64 \times 64$ 的灰度图像，采集 $m = 1500$ 个随机傅里叶样本，基于小波系数进行 $L_1$ 重建，并与零填充逆 FFT 的结果对比。

总结#

参考文献#

• Candes, E. & Wakin, M. (2008). “An Introduction to Compressive Sampling.” IEEE Signal Processing Magazine 25(2), 21–30.
• Candes, E., Romberg, J. & Tao, T. (2006). “Robust Uncertainty Principles.” IEEE Trans. Information Theory 52(2), 489–509.
• Tibshirani, R. (1996). “Regression Shrinkage and Selection via the Lasso.” JRSS-B 58, 267–288.
• Efron, B., Hastie, T., Johnstone, I. & Tibshirani, R. (2004). “Least Angle Regression.” Annals of Statistics 32(2), 407–499.
• Beck, A. & Teboulle, M. (2009). “A Fast Iterative Shrinkage-Thresholding Algorithm for Linear Inverse Problems.” SIAM J. Imaging Sciences 2(1), 183–202.
• Blumensath, T. & Davies, M. (2009). “Iterative Hard Thresholding for Compressed Sensing.” Applied and Computational Harmonic Analysis 27(3), 265–274.
• Foucart, S. & Rauhut, H. (2013). A Mathematical Introduction to Compressive Sensing. Birkhauser.
• Hastie, T., Tibshirani, R. & Wainwright, M. (2015). Statistical Learning with Sparsity: The Lasso and Generalizations. CRC Press.