稀疏矩阵与压缩感知 -- 少即是多的数学奇迹

「少即是多」的奇迹

一张 24 兆像素的原始照片大约 70 MB，JPEG 压到几百 KB——压缩比上百倍——你看不出区别。传统 MRI 扫描要 30 分钟；现在的压缩感知 MRI 只要 5 分钟，图像质量一样。

两个奇迹背后是同一个引擎：稀疏性。绝大多数自然信号，只要换到合适的基底下，真正"有信息量"的系数就那么几个，其余几乎全是零。

压缩感知把这个观察变成了算法。既然信号本身稀疏，那我们用远少于经典采样定理要求的测量，就能把它精确恢复出来。本章讲清楚为什么。

本章你会学到

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

预备知识

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

稀疏性：普遍存在的模式

定义

向量$\vec{x} \in \mathbb{R}^n$$k$-稀疏，是指它至多有$k$个非零元素：$\|\vec{x}\|_0 \;=\; \bigl|\{\,i : x_i \neq 0\,\}\bigr| \;\leq\; k.$$L_0$“范数"只是数非零元素的个数。它其实不是真正的范数——不满足齐次性（$\|2x\|_0 = \|x\|_0$）——但它是衡量稀疏性最直接的量。

实际信号几乎不会严格稀疏，而是可压缩：把系数按大小排序后，序列衰减得很快，截断到前$k$个就能得到几乎完美的近似。JPEG 就是这么干的：算 DCT，留下最大的几百个系数，其余全部丢掉。

现实里的稀疏

字典下的稀疏表示

很多信号在标准基下并不稀疏，但换个基底就稀疏了。挑一个字典$D \in \mathbb{R}^{n \times p}$，它的列称为原子，然后求：$\vec{x} \;=\; D\vec{\alpha}, \qquad \vec{\alpha}\;\text{稀疏}.$常见字典：

当$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$不可解

在字典$D$里找$\vec{x}$的最稀疏表示：$\min_{\vec{\alpha}}\;\|\vec{\alpha}\|_0 \quad \text{s.t.}\quad D\vec{\alpha} = \vec{x}.$这是NP 难的。最朴素的做法要枚举$\binom{p}{k}$种支撑——组合爆炸。一般情形下没有已知的多项式时间算法。

###$L_1$作为凸松弛

把$\|\cdot\|_0$换成$\|\cdot\|_1 = \sum |\alpha_i|$：$\min_{\vec{\alpha}}\;\|\vec{\alpha}\|_1 \quad \text{s.t.}\quad D\vec{\alpha} = \vec{x}.$这是个线性规划，多项式时间内就能解。神奇的地方在于：在对$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$容易在一组相关特征里挑一个、把其余全部踢掉。把两个惩罚加起来：$\min_{\vec{\beta}}\; \tfrac{1}{2}\|X\vec{\beta} - \vec{y}\|^2 + \lambda_1\|\vec{\beta}\|_1 + \lambda_2\|\vec{\beta}\|_2^2$既保留$L_1$的稀疏性，又借$L_2$获得稳定性——这就是Elastic Net，高维回归里更稳妥的默认选择。

LASSO 回归

定义

LASSO——Least Absolute Shrinkage and Selection Operator（Tibshirani, 1996）：$\min_{\vec{\beta}}\;\tfrac{1}{2}\|X\vec{\beta} - \vec{y}\|^2 \;+\; \lambda\|\vec{\beta}\|_1.$同时做两件事：拟合$\vec{y}$，以及把不重要的系数压成精确的零。所以 LASSO 一步完成回归与特征选择。

软阈值

正交设计（$X^\top X = I$）时 LASSO 有闭式解：$\hat\beta_j \;=\; \mathcal{S}_\lambda\!\bigl(\hat\beta_j^{\text{OLS}}\bigr) \;=\; \operatorname{sign}\!\bigl(\hat\beta_j^{\text{OLS}}\bigr)\,\bigl(|\hat\beta_j^{\text{OLS}}| - \lambda\bigr)_+.$这就是软阈值：把每个 OLS 系数往零方向压缩$\lambda$，绝对值小于$\lambda$的直接压成零。对比硬阈值——保留大系数原值、把小的直接归零——软阈值更平滑、更稳定。

LASSO 路径

把$\lambda$从$\infty$扫到$0$：

-$\lambda$很大时：所有系数都是零。 -$\lambda$减小：特征一次进一个，每个特征在自己的"残差相关性超过阈值"时刻进入活动集。 -$\lambda \to 0$：退化为普通最小二乘。

Efron, Hastie, Johnstone, Tibshirani（2004）证明的一个漂亮事实：路径关于$\lambda$分段线性。LARS 算法利用这一点，用一次 OLS 的代价就能算出整条路径。

真正活跃的特征（粗线）很早就进来、值也大；伪特征要么始终不进，要么很晚进来、系数很小。$\lambda$用交叉验证选——通常用"1-SE 法则”：在最小误差$\lambda$附近一个标准误差以内，挑最大的那个$\lambda$，偏向更稀疏的解。

压缩感知

革命性的想法

香农-奈奎斯特告诉你：要恢复一个带限信号，采样率必须不低于带宽的两倍。没有商量。

压缩感知告诉你：如果信号在某个基下是$k$-稀疏的，你只需要$m \;\sim\; k\,\log\!\frac{n}{k}$次线性测量就够了。一段长度 256、小波系数只有 10 个非零的信号，$k\log(n/k) \approx 32$次测量足矣——而不是 256。

这不是免费午餐——它能成立，是因为信号模型更强了。香农-奈奎斯特只假设"带限"，压缩感知假设"$k$-稀疏"。模型越窄，采样越便宜。

测量模型$\vec{y} \;=\; \Phi \vec{x} \;+\; \vec{e}$其中$\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$阶 RIP（常数$\delta_k \in (0,1)$），如果对所有$k$-稀疏向量$\vec{x}$：$(1 - \delta_k)\|\vec{x}\|_2^2 \;\leq\; \|\Phi\vec{x}\|_2^2 \;\leq\; (1 + \delta_k)\|\vec{x}\|_2^2.$直白说：$\Phi$作用在稀疏向量上几乎像等距——它不会把两个不同的稀疏向量挤到（几乎）同一个像。

左：在多个稀疏度下，对大量随机$k$-稀疏向量的$\|\Phi\vec{x}\|/\|\vec{x}\|$作直方图。分布紧密集中在 1 附近——范数被几乎完美保留。右：经验最坏情况畸变$\delta_k$随$k$增长。某个稀疏度之后会越过$\sqrt{2}-1$；在它之下，恢复有保证。

恢复定理

定理（Candes-Tao, 2005）。 若$\delta_{2k}(\Phi) < \sqrt{2} - 1 \approx 0.414$，则任何$k$-稀疏的$\vec{x}$都是下面问题的唯一解：$\min \|\vec{z}\|_1 \quad \text{s.t.}\quad \Phi\vec{z} = \vec{y}, \qquad \vec{y} = \Phi\vec{x}.$有噪声时——用基追踪去噪 BPDN：$\min\|\vec{z}\|_1$s.t.$\|\Phi\vec{z}-\vec{y}\| \leq \epsilon$——恢复误差被一个常数倍$\epsilon$界住，温和地随噪声退化。

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

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

算法

ISTA：迭代收缩-阈值

把近端梯度法套到$\tfrac{1}{2}\|\Phi x - y\|^2 + \lambda\|x\|_1$上：$\boxed{\;\vec{x}^{(t+1)} \;=\; \mathcal{S}_{\lambda/L}\!\left(\vec{x}^{(t)} - \tfrac{1}{L}\Phi^\top\!\bigl(\Phi \vec{x}^{(t)} - \vec{y}\bigr)\right)\;}$其中$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 加速

加一个动量项，收敛率从$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：迭代硬阈值

如果你提前知道稀疏度$k$，把软阈值换成"保留绝对值最大的$k$个"：$\vec{x}^{(t+1)} \;=\; H_k\!\left(\vec{x}^{(t)} - \tfrac{1}{L}\Phi^\top\!\bigl(\Phi \vec{x}^{(t)} - \vec{y}\bigr)\right).$IHT 直接攻打$L_0$问题，在类 RIP 条件下收敛到$k$-稀疏解。每次迭代极便宜：一次矩阵-向量乘加一次部分排序。

四张子图分别是第 0、2、5、30 次迭代的解（淡色背景是真实信号）。两次迭代后支撑集大致正确，30 次后幅值也对上了。残差范数$\|\Phi x_t - y\|_2$下降了好几个数量级。

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$s.t.$\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 对比。

本章要点

系列导航

上一篇： 第十一章 – 矩阵微积分与优化

下一篇： 第十三章 – 张量与多线性代数

本文是「线性代数的本质」系列第 12 章，共 18 章。

参考文献

• 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.