压缩感知与稀疏建模
本教程定位: Steve Brunton 和 张颢 网课课程的学习笔记。
阅读方式：第 1–5 章建立模型、采样理论与恢复条件；第 6–12 章系统学习凸优化与贪婪恢复；第 13–18 章开展实验并了解工程应用；第 19–22 章用于复习、练习及延伸阅读；第 23 章附录提供补充的数学证明与完整实验代码。
0. 阅读指南
0.1 教程导航
| 内容 | 章节 | 核心 |
|---|---|---|
| 第一部分：为什么可能 | 1–5 | 稀疏表示、压缩测量、欠定模型、相干性与 RIP |
| 第二部分：凸优化恢复 | 6–7 | BP、BPDN、LASSO、LP 与 ISTA |
| 第三部分：贪婪恢复 | 8–12 | OST、OMP、CoSaMP 的步骤、证明与数值细节 |
| 第四部分：可复现实验 | 13–16 | 七个实验、随机采样、图像恢复、稳健回归 |
| 第五部分：应用扩展 | 17–18 | SRC、RPCA、传感器优化与 SSPOC |
| 第六部分：复习与参考 | 19–22 | 工程适用性、8 周安排与练习 |
| 附录：证明与代码 | 23 | OST、OMP、NSP、RIP 稳定性等证明及七个实验源码 |
0.2 课程来源与内容对照
| 来源 | 核心内容 | 本教程位置 |
|---|---|---|
| Steve Brunton 第 1–4 讲 | 信号空间、稀疏表示、压缩感知模型、欠定逆问题 | 第 1–4 章 |
| Brunton 第 5–6 讲 | L1 与 L2 欠定系统求解、Python 和 MATLAB | 第 6、14 章 |
| Brunton 第 7–10 讲 | 奈奎斯特采样、随机与均匀采样、双频恢复 | 第 3、15 章 |
| Brunton 第 11–12 讲 | L1 直觉、非相干性、采样数、RIP、何时可用 | 第 5–9 章 |
| 张颢《稀疏恢复》课堂 | OST、OMP、CoSaMP；相干性、支撑集、投影及恢复条件 | 第 5、8–12 章 |
| Brunton 第 13–16 讲 | L1 鲁棒回归、LASSO、可解释模型 | 第 16 章 |
| Brunton 第 17–20 讲 | SRC、鲁棒 PCA、SVD+QR 传感器布局 | 第 17–18 章 |
| Brunton 第 21–22 讲 | 简约性、压缩感知与工程系统 | 第 1、19–22 章 |
0.3 学习成果
完成后应能够：（1）解释普通数据压缩与压缩感知的区别；（2）从采样系统写出 $y=C\Psi s+e$；（3）理解欠定问题为何需要结构先验；（4）区分 $\ell_0,\ell_1,\ell_2$ 与残差损失；（5）理解非相干性、spark、RIP 及其相互关系；（6）推导 OST/OMP 的相干性充分条件；（7）实现 BP、ISTA、OST、OMP、CoSaMP；（8）重现随机采样双音恢复；（9）辨析 LASSO、SRC、RPCA 和定制传感与狭义压缩感知的联系。
课程模块｜第一部分：压缩感知为什么可能
1. 从自然信号的可压缩性说起
1.1 一个 20×20 的思想实验
Brunton 的课程首先问：为什么图像、音频、视频往往可以压缩？ 考虑一张 $20\times20$ 的二值图像，400 个像素各有黑白两种选择，因此可能的图像数量为
\[2^{400}\approx2.58\times10^{120}.\]随机独立地设置每个像素通常得到类似雪花的噪声图案，而具有物体、边缘、纹理和相互关系的自然图像包含大量统计结构。自然图像在这个极大的表示空间中不是均匀分布的。因此，我们希望找到一种能够利用结构的紧凑描述。
直觉不能代替数学定理。 “自然图像占据小区域”是一种启发，并不意味着每张自然图像在任意傅里叶基上都恰好稀疏。某一变换是否有效，必须由具体数据与近似误差验证。
1.2 稀疏与可压缩并不相同
设系数向量 $s\in\mathbb R^P$，零“范数”定义为非零分量个数：
\[\lVert s\rVert_0=\#\{i:s_i\neq0\}.\]若 $\lVert s\rVert_0\le K$，称为 $K$-稀疏；非零指标集为 $S=\operatorname{supp}(s)$。例如 $s=(0,4,0,0,-2,0)^T$ 是 2-稀疏。
实际信号不一定严格只有 $K$ 个非零系数，而可能有一串逐渐变小的系数。令 $s_K$ 为只保留 $K$ 个绝对值最大分量的向量，则
\[\sigma_K(s)_p=\min_{\lVert z\rVert_0\le K}\lVert s-z\rVert_p=\lVert s-s_K\rVert_p\]表示最佳 $K$ 项近似误差。若 $\sigma_K(s)_p$ 随 $K$ 增加迅速下降，我们称信号可压缩。压缩感知的稳定恢复理论可以同时处理严格稀疏和近似稀疏信号。
1.3 稀疏性与简约建模
Brunton 第 21–22 讲把稀疏性放在更大的背景下：一个描述模型如果仅需少数自由度、少量特征或少量有效参数，往往更易解释、存储和估计。但模型不能无限简化，否则会出现欠拟合。稀疏性是一种有用的先验约束，而非“参数越少越好”的普遍法则。
思考：一段含 1000 个样本、只有 3 个明显频率分量的音频，在时域稀疏吗？在傅里叶域呢？两者为何不同？
2. 稀疏表示、变换基与超完备字典
2.1 通用表示
把原始信号记为 $x\in\mathbb R^N$，选取字典 $\Psi\in\mathbb R^{N\times P}$，使
\[\boxed{x=\Psi s,\qquad \lVert s\rVert_0\le K.}\]这里 $s\in\mathbb R^P$ 是稀疏系数。当 $\Psi$ 是正交基（$P=N$，$\Psi^T\Psi=I$）时，$s=\Psi^Tx$，变换与逆变换都很直接。对复数傅里叶基，需把转置改为共轭转置 $\Psi^*$。
典型字典：
| 变换 | 常见适用信号 | 需要注意的限制 |
|---|---|---|
| DFT / FFT | 少量正弦波、周期成分 | 边缘、突变可能使系数分散 |
| DCT | 具有局部平滑性的图像块 | 纹理、块边界可能导致伪影 |
| 小波 | 边缘、突变、分段平滑信号 | 小波类型、层级与边界处理影响精度 |
| PCA / SVD 模态 | 已知类别的图像、流场、温度场 | 依赖训练集；低维表示不等于系数稀疏 |
| 学习字典 | 特定领域数据的局部结构 | 训练成本与泛化能力需要验证 |
勘误： Brunton 口误将 JPEG 与傅里叶变换混在一起。通常的基线 JPEG 采用分块 DCT，而 JPEG 2000 采用小波变换；FFT 是计算 DFT 的快速算法，不等同于 DCT。
2.2 低维与稀疏的区别
如果通过 SVD 得到 $r$ 维基 $U_r\in\mathbb R^{N\times r}$，可以写成
\[x\approx U_r a,\qquad r\ll N.\]这是一种低维表示，但 $a\in\mathbb R^r$ 的全部分量都可能非零。它与 $N$ 维变换中的 $K$-稀疏向量并不完全相同。Brunton 关于定制基和传感器布局的讲授将利用这一点。
2.3 超完备字典与稀疏编码
张颢老师的课程特别重视超完备字典。当 $P>N$ 时，$\Psi$ 拥有比信号维数更多的列向量（原子），因此不能由两两正交的非零列组成。丰富的原子有助于描述各种信号，但稀疏编码问题更难：
\[\min_s\lVert s\rVert_0\quad\text{subject to}\quad x=\Psi s.\]较大的字典不自动保证更好的恢复：高度相关的原子可能难以区分，且系数可能不唯一。后面将用互相干性 $\mu$ 和 RIP 定量分析这一风险。
3. 普通压缩、奈奎斯特采样与压缩感知
3.1 “先测量再压缩”与“直接少量测量”
传统变换编码通常遵循：完整测量 $x$ → 变换得到 $s$ → 保留重要系数、量化和编码 → 解码并逆变换。它利用的是测量后的可压缩性。
压缩感知则希望：直接采集少量线性测量 $y$，然后从这些测量重建稀疏表示 $s$，最终恢复 $x$。 两者对测量硬件、数据采集成本和可恢复性提出不同要求。
3.2 奈奎斯特—香农采样定理
对严格带限至 $B$ Hz 的连续时间信号，经典均匀采样定理要求采样率高于 $2B$ Hz（在边界处需注意数学假设）。若均匀采样率不足，不同频率可能在样本上产生完全相同的观测，这称为混叠。
例如以 $f_s=64$ Hz 均匀采样时，对整数 $n$，
\[\cos\!\left(2\pi\frac{37}{64}n\right)=\cos\!\left(2\pi\frac{27}{64}n\right).\]仅凭这些均匀样本，无法区分 37 Hz 和 27 Hz 的分量。
3.3 突破奈奎斯特
Brunton 第 7–10 讲用随机时间点观测一个仅含 97 Hz 和 777 Hz 的双音信号，尽管平均测量速率明显低于通常按最高频率计算的均匀采样速率，仍可以在知道精确采样时刻、给定高分辨率候选频率网格的条件下，恢复少数频域系数。
这不是推翻奈奎斯特定理。 原定理对给定带宽内的通用带限信号提供采样保证；压缩感知增加了“频域高度稀疏”等先验，并采用特殊测量设计和非线性重建算法，处理的是更小的信号类别。低平均采样率并不等于硬件可以忽略前端带宽、时间分辨率或噪声。
4. 统一测量模型：从数据到欠定逆问题
4.1 统一符号约定
全文采用
\[\boxed{y=Cx+e=C\Psi s+e=As+e},\qquad A=C\Psi.\]| 量 | 维度 | 物理或数学意义 |
|---|---|---|
| $x$ | $N\times1$ | 原始信号或高维状态 |
| $\Psi$ | $N\times P$ | 变换基或过完备字典 |
| $s$ | $P\times1$ | 待恢复的稀疏系数 |
| $C$ | $M\times N$ | 测量或采样算子 |
| $A=C\Psi$ | $M\times P$ | 有效恢复矩阵 |
| $y$ | $M\times1$ | 已知观测数据 |
| $e$ | $M\times1$ | 噪声、量化误差或模型误差 |
通常 $M<P$，所以 $As=y$ 是欠定方程组。张颢老师直接使用 $A$ 与稀疏未知向量；Brunton 常显式分开采样算子 $C$、表示字典 $\Psi$ 和系数向量 $s$。两种写法是一回事。
4.2 为什么少量方程可以恢复较多未知数
单独看 $As=y$，若 $A$ 行满秩且 $M<P$，零空间维数至少 $P-M$，因此无穷多组 $s$ 可以产生相同的 $y$。压缩感知引入稀疏先验，在无穷多解中寻找满足 $\lVert s\rVert_0\le K$ 的向量。
例子：一个 $M=80$、$P=400$ 的线性系统，单从方程数看无法唯一识别 400 个变量；但如果只有 $K=5$ 个变量非零，而且矩阵满足足够好的可区分性，稀疏恢复可能成功。
4.3 三种常见观测
- 随机像素/时间点采样：$C$ 由单位矩阵的某些行构成；它在原始信号坐标中选取位置。
- 高斯随机投影：$C_{ij}$ 或有效矩阵 $A_{ij}$ 服从适当缩放的独立高斯分布；数学分析方便，但硬件实现可能昂贵。
- 部分傅里叶测量：MRI 等应用可能在傅里叶域采集部分样本，测量算子与信号稀疏字典不能混为一谈。
“随机”不是每一种问题的必要条件；核心是采样算子与相关的信号结构之间是否具备足够可辨识性。
5. 唯一性、互相干性、非相干性与 RIP
5.1 无法恢复的简单例子
若 $A$ 存在两列完全相同的原子 $a_1=a_2$，那么 $s=(1,0,\ldots)$ 和 $s’=(0,1,\ldots)$ 的观测相同。即使两个向量都是 1-稀疏，也无法从 $y$ 判断真实系数的位置。
所以我们需要的不是“单纯足够多的测量”，而是“不同稀疏表示在测量后仍可区分”。
5.2 spark 给出的唯一性判据
定义 $\operatorname{spark}(A)$ 为 $A$ 中线性相关的最少列数。若
\[\boxed{\operatorname{spark}(A)>2K},\]则每一个至多 $K$-稀疏的解（若存在）都是唯一的。证明：如果 $s\ne s’$ 均为 $K$-稀疏且 $As=As’$，则 $h=s-s’$ 至多 $2K$-稀疏，满足 $Ah=0$，这与任意少于 $\operatorname{spark}(A)$ 列都线性无关相矛盾。
该条件容易理解，但实际计算 spark 一般很困难。其与相干性的关系 $\operatorname{spark}(A)\ge 1+1/\mu(A)$ 的证明见附录第 23.2 节。
5.3 列互相干性：张颢课程的重要工具
归一化 $A=[a_1,\ldots,a_P]$ 的列，使 $\lVert a_i\rVert_2=1$，定义
\[\boxed{\mu(A)=\max_{i\ne j}\lvert a_i^*a_j\rvert.}\]$\mu$ 越小，不同原子越容易区分。经典充分条件
\[\boxed{(2K-1)\mu(A)<1}\]可以保证 OMP 在无噪声、已知稀疏度等条件下精确恢复。注意这是充分而非必要条件。
当 $P>M$ 时，单位范数列不可能全部两两正交。Welch 下界给出：
\[\mu(A)\ge\sqrt{\frac{P-M}{M(P-1)}}.\]这表明相干性既受测量维数影响，也受字典列数影响；不能把它统一近似成 $1/\sqrt P$。
5.4 两种“非相干”不要混淆
张颢分析的 $\mu(A)$ 测量的是有效恢复矩阵内部不同列的相似性。
Brunton 第 12 讲强调的是测量基与表示基之间的非相干性。例如位置采样向量 $e_i$ 与傅里叶基 $\psi_j$ 满足
\[\lvert e_i^*\psi_j\rvert=\frac1{\sqrt N},\]因此单点测量不会只“看见”某一个傅里叶频率。更一般地，可定义交叉相干性 $\max_{i,j}\lvert\langle\phi_i,\psi_j\rangle\rvert$（对两组单位范数基）。二者关联很深，但不是同一个数值指标。
5.5 受限等距性质（RIP）
当所有 $K$-稀疏向量 $h$ 都满足
\[\boxed{(1-\delta_K)\lVert h\rVert_2^2\le\lVert Ah\rVert_2^2\le(1+\delta_K)\lVert h\rVert_2^2}\]时，称 $A$ 满足 $K$ 阶 RIP，$\delta_K$ 为其受限等距常数。它说的不是 $A$ 在整个高维空间上近似正交，而是对稀疏向量近似保持长度。
若 $\delta_{2K}<1$，不同的 $K$-稀疏向量不可能映射为同一观测，故得到唯一性；然而唯一性还不意味着具体算法一定能找出该解。
5.6 测量次数：为什么出现 $K\log(P/K)$
对适当归一化的高斯或某些次高斯随机矩阵，存在保证 RIP/稳定稀疏恢复的经典结果，其测量规模约为
\[\boxed{M\gtrsim c\,K\log(P/K)}\]这里 $c$ 取决于目标成功概率、所需的 RIP 常数和矩阵分布，$\gtrsim$ 是渐近量级，不是可直接代入得到工程采样率的精确公式。某些结构化矩阵（包括随机部分傅里叶）可能需要额外的对数因子。
例题。 $P=10000$、$K=20$ 时，$K\log(P/K)\approx20\log(500)\approx124.3$。这表示一个测量规模的基准量级，不能据此断言“124 个测量一定能恢复”。对高斯随机矩阵，$K\log(P/K)$ 量级的证明概要见附录第 23.10 节。
课程模块｜第二部分：稀疏恢复的优化视角
6. 为什么选 $\ell_1$ 而不是 $\ell_2$？
6.1 最理想却难求的 $\ell_0$ 优化
无噪声模型中的最稀疏恢复为
\[\boxed{\min_s\lVert s\rVert_0\quad\mathrm{s.t.}\quad As=y.}\]若寻找大小为 $K$ 的支撑集，潜在组合数为 $\binom PK$。这一非凸组合问题对一般大规模矩阵难以穷举。
6.2 为什么伪逆不是稀疏恢复的终点
若 $A$ 行满秩，最小 $\ell_2$ 范数解为
\[s_{\ell_2}=A^\dagger y=A^T(AA^T)^{-1}y.\]它在满足测量的所有向量中欧氏长度最短，但通常会分散到许多坐标，不必稀疏。Brunton 第 5–6 讲的数值实验对比了这个密集解与 $\ell_1$ 稀疏解。
6.3 基追踪 BP
用凸函数 $\lVert s\rVert_1=\sum_i\lvert s_i\rvert$ 替代 $\lVert s\rVert_0$：
\[\boxed{\min_s\lVert s\rVert_1\quad\mathrm{s.t.}\quad As=y.}\]称为 Basis Pursuit（BP）。它并非对所有矩阵都与 $\ell_0$ 等价，但在满足适当 RIP、零空间性质等条件时，可以精确恢复真正的稀疏向量。BP 的充要恢复条件之一——零空间性质（NSP）——将在附录第 23.5 节证明。
二维几何直觉：$\ell_2$ 范数球的边界光滑，最优接触点未必位于坐标轴；$\ell_1$ 范数球在坐标轴方向存在尖角，线性约束更容易在尖角处相切，从而得到零系数。这个直觉说明“促进稀疏”，但不是完整证明。
6.4 有噪声时：BPDN、LASSO 与约束形式
若 $y=As+e$ 且已知 $\lVert e\rVert_2\le\varepsilon$，可用
\[\boxed{\min_s\lVert s\rVert_1\quad\mathrm{s.t.}\quad\lVert As-y\rVert_2\le\varepsilon.}\]称为 Basis Pursuit Denoising（BPDN）。另一常见表述为
\[\boxed{\min_s\frac12\lVert As-y\rVert_2^2+\lambda\lVert s\rVert_1.}\]即 LASSO 的惩罚形式。$\lambda$ 控制数据拟合与系数稀疏度的权衡：更大 $\lambda$ 通常带来更多系数收缩；过大则欠拟合。
这里不能混淆两种 $\ell_1$ 用法：
- 惩罚 系数 $\lVert s\rVert_1$：主要目的是稀疏表示、变量选择。
- 惩罚 残差 $\lVert As-y\rVert_1$：主要目的是降低大异常值对拟合的影响。
同一种范数出现的位置不同，统计含义就不同。这是理解 Brunton 第 13–16 讲的关键。
6.5 稳定恢复的定量表述
在适当 RIP 等假设下，BPDN 的典型误差界为
\[\boxed{\lVert\hat s-s\rVert_2\le C_0\varepsilon+C_1\frac{\sigma_K(s)_1}{\sqrt K}.}\]第一项来自测量噪声，第二项来自原信号不完全稀疏。若无噪声且严格 $K$-稀疏，两项均为零，就可得到精确恢复。经典的一个保守充分条件是 $\delta_{2K}<\sqrt2-1$；文献中还有更精细的界。附录第 23.6 节提供这一充分条件下的完整“锥约束—管约束—RIP”推导及可取的误差常数。
7. 如何实际求解凸问题：LP、软阈值、ISTA
7.1 把 BP 改写成线性规划
对于实数向量，引入辅助变量 $t_i\ge0$，使 $-t_i\le s_i\le t_i$。则
\[\begin{aligned}\min_{s,t}&\quad\mathbf1^Tt\\ \mathrm{s.t.}&\quad As=y,\\&\quad -t\le s\le t,\quad t\ge0. \end{aligned}\]这就是标准线性规划，可用 SciPy 的 scipy.optimize.linprog 求解。复数系数需要单独处理实虚部或采用支持复数的凸优化工具，不能直接套用上述实数 LP。
7.2 软阈值算子
定义
\[\operatorname{soft}(z,\tau)=\operatorname{sign}(z)\max(\lvert z\rvert-\tau,0).\]对向量按坐标应用。它正是 $\ell_1$ 范数的近端映射：
\[\operatorname{prox}_{\tau\lVert\cdot\rVert_1}(v)=\arg\min_s\left(\frac12\lVert s-v\rVert_2^2+\tau\lVert s\rVert_1\right).\]证明可以逐坐标考虑 $s>0$、$s<0$ 和 $s=0$ 的次梯度最优性条件。对正值区域，一阶条件为 $s-v+\tau=0$；对负值区域为 $s-v-\tau=0$；中间区域的最优解为零，从而拼成软阈值公式。
7.3 ISTA：迭代阈值收缩
设 $f(s)=\tfrac12\lVert As-y\rVert_2^2$，则 $\nabla f(s)=A^T(As-y)$；其梯度 Lipschitz 常数可取 $L=\lVert A\rVert_2^2$。因此 ISTA 的一步为
\[\boxed{s^{t+1}=\operatorname{soft}\!\left(s^t-\frac1L A^T(As^t-y),\,\frac\lambda L\right).}\]这是一阶凸优化算法，区别于后面直接操作支撑集的贪婪算法。FISTA 在 ISTA 基础上引入加速步骤，通常能更快逼近最优解。
重要区分： OMP 或 CoSaMP 并不是“快速求解 BP 的另一种实现”；它们分别对应不同的迭代恢复策略和理论保证。
课程模块｜第三部分：张颢老师重点讲授的贪婪恢复算法
8. 一步阈值法 OST：从相关性开始
8.1 目标：找出支撑集
令 $A=[a_1,a_2,\ldots,a_P]$，列已经归一化，真实信号为
\[y=As=\sum_{j\in S}s_j a_j.\]如果 $a_i$ 与测量向量 $y$ 内积较大，可以把它看作重要原子的一个候选。因此 OST（One-Step Thresholding）的核心步骤非常简单：
\[z=A^*y,\qquad \hat S=\operatorname{TopK}(\lvert z\rvert).\]只用相关性挑选支撑集，随后可用最小二乘在 $\hat S$ 上重新估计系数：
\[\hat s_{\hat S}=A_{\hat S}^{\dagger}y,\qquad\hat s_{\hat S^c}=0.\]算法伪代码
1
2
3
4
5
6
输入：A, y, 稀疏度 K
z ← A* y # 共轭转置；实数时为 Aᵀy
S ← |z| 最大的 K 个指标
s[S] ← least_squares(A[:, S], y) # 可选的系数重估
s[其余] ← 0
输出：s, S
8.2 OST 的支撑恢复充分条件
设 $\mu=\mu(A)$，$s_{\min}=\min_{i\in S}\lvert s_i\rvert$。如果真实支撑集 $S$ 上的所有相关性都严格超过虚假原子的相关性，即
\[\min_{i\in S}\lvert a_i^*y\rvert>\max_{j\notin S}\lvert a_j^*y\rvert,\]则 OST 选择前 $K$ 个相关性时必然得到正确支撑集。利用三角不等式与列互相干性，可以证明：在无噪声、列归一化条件下，只要
\[\boxed{\frac{s_{\min}}{\lVert s\rVert_1}>\frac{2\mu(A)}{1+\mu(A)},} \tag{OST}\]OST 就能精确识别支撑集。这个条件同时限制矩阵列相关性和非零系数的动态范围：若存在极弱的真实分量，一次相关性排序仍可能错选。
当所有非零系数的幅值相等时，上式化为 $(2K-1)\mu<1$。逐项下界、错误原子上界和含噪声推广的完整证明放在附录第 23.3 节。
8.3 OST 的失败机制
一次排序没有扣除已解释的信号。强原子的相关性会“污染”与它相近的其他原子，从而掩盖真实但幅度较弱的原子。一旦选错，OST 没有第二轮纠正机会。其优点是速度快、结构简单，适合初始化和比较基线。
9. OMP：正交匹配追踪与精确恢复条件
9.1 为什么要从 $y$ 改成残差 $r_t$
OST 始终只看 $A^*y$。OMP（Orthogonal Matching Pursuit）希望先找到最重要的一个原子，扣除它已经解释的部分，再研究剩余信号。
初始化
\[S_0=\varnothing,\qquad r_0=y.\]在第 $t$ 轮：
\[i_t=\arg\max_{i\notin S_t}\lvert a_i^*r_t\rvert,\qquad S_{t+1}=S_t\cup\{i_t\},\]然后求最小二乘解
\[\hat s_{S_{t+1}}=\arg\min_z\lVert y-A_{S_{t+1}}z\rVert_2^2,\]再更新
\[\boxed{r_{t+1}=y-A_{S_{t+1}}\hat s_{S_{t+1}}.}\]这里“匹配”指按残差相关性寻找原子；“追踪”指逐轮更新；“正交”指最小二乘使残差垂直于已选原子的张成空间。
9.2 为什么残差会与所选原子正交？
若 $A_S$ 列满秩，最小二乘的正规方程为
\[A_S^*(y-A_S\hat s_S)=0.\]因此
\[\boxed{A_S^*r=0.}\]又令 $P_S=A_SA_S^\dagger$ 为到 $\operatorname{span}(A_S)$ 的正交投影，便有
\[r=(I-P_S)y.\]它可以帮助我们避开已解释的成分。在实际代码里还应显式排除已选指标，以防浮点误差导致异常重复选择。
9.3 OMP 的精确恢复条件：ERC
设真实支撑集为 $S$、$\lvert S\rvert=K$，$A_S$ 列满秩。若
\[\boxed{\max_{j\notin S}\lVert A_S^\dagger a_j\rVert_1<1,}\tag{ERC}\]则称该支撑集满足 Exact Recovery Condition（ERC）。它意味着：对任何属于真实列空间且尚未消失的残差，真实支撑集上的最大相关性严格大于支撑集外的最大相关性。由此可以归纳证明，OMP 逐次选择的原子均来自 $S$，最终恢复正确系数。
ERC 是依赖真实支撑集 $S$ 的条件；分析实际问题时，由于 $S$ 往往未知，很难直接检查。完整归纳证明和关键矩阵范数不等式见附录第 23.4 节。
9.4 从 ERC 推出容易理解的相干性条件
对列归一化字典，令 $G=A_S^*A_S=I+E$。由于不同列内积的模最多为 $\mu(A)$，可知 $\lVert E\rVert_1\le (K-1)\mu(A)$。利用 Neumann 级数控制 $G^{-1}$，可以进一步证明
\[\max_{j\notin S}\lVert A_S^\dagger a_j\rVert_1\le \frac{K\mu(A)}{1-(K-1)\mu(A)}.\]若右侧小于 1，即
\[\boxed{(2K-1)\mu(A)<1,}\tag{OMP}\]便可保证无噪声 OMP 的精确恢复。它是充分条件而非必要条件：不满足条件不代表 OMP 一定失败。Neumann 级数的合法性、逆矩阵范数界及其对 ERC 的推导均放在附录第 23.4 节。
9.5 OMP 的致命限制：不撤销早期误选
标准 OMP 的支撑集满足
\[S_t\subseteq S_{t+1}.\]如果第 2 次迭代错误选中一个与真实原子接近但不正确的列，它将一直留在解中。面对高相干、低信噪比或误差严重的场景，这种不可逆的选择会降低恢复性能。
10. CoSaMP：从不可撤销选择到可修正支撑集
10.1 思路
CoSaMP（Compressive Sampling Matching Pursuit）把 OST 的“批量寻找候选”与 OMP 的“最小二乘分析”结合起来，再增加剪枝，让旧指标有机会退出。张颢课程用“候选人面试”的类比解释了这一机制。
10.2 标准流程
给定目标稀疏度 $K$，初始化 $s^{(0)}=0,r^{(0)}=y$。每一轮：
- 识别候选：$\Omega_t=\operatorname{Top}_{2K}(\lvert A^*r^{(t)}\rvert)$；
- 合并支持集：$T_t=\Omega_t\cup\operatorname{supp}(s^{(t)})$，其大小至多 $3K$；
- 联合估计：$b_{T_t}=\arg\min_z\lVert y-A_{T_t}z\rVert_2^2$，其他坐标为零；
- 剪枝保留：$s^{(t+1)}=H_K(b)$，$H_K$ 只保留最大的 $K$ 个绝对值分量；
- 更新残差：$r^{(t+1)}=y-As^{(t+1)}$。
实际实现中也常在第 4 步完成后对保留的支撑集再做一次最小二乘来减小系数偏差。配套实验代码使用了这个消偏变体；它不应与原论文的逐步表述混为一谈。
10.3 一组 3-稀疏支撑集的演化
设原子指标有 100 个。第一次选出 6 个候选，最小二乘后保留 $S_1={1,7,18}$。第二次由残差识别新候选，合并后最小二乘，留下的可能变成 $S_2={4,7,25}$。这样，原本选中的 1、18 被移除，4、25 被加入。CoSaMP 的关键不是仅仅“多选”，而是“候选扩展+竞争+淘汰”。
10.4 性能保证与失败场景
在足够强的 RIP 条件下，CoSaMP 对噪声和非严格稀疏信号具有收敛及稳定误差界；具体常数与所引用的定理版本有关。经典来源是 Needell 与 Tropp 的论文 CoSaMP: Iterative Signal Recovery from Incomplete and Inaccurate Samples。
但当 $A_T$ 高度病态时，最小二乘系数可能非常不稳定；有严重离群值时，平方误差会放大污染；剪枝还可能导致某轮残差增大。它具有修正能力，但不能无条件修正所有错误，也不保证所有任务中都优于 OMP。收敛证明中为何需要“候选识别—联合估计—剪枝”三个步骤，见附录第 23.8 节。
11. 三种贪婪算法与凸优化如何选择
| 方法 | 每轮候选 | 需要估计 $K$ 吗 | 能否撤销选中指标 | 主要目标 | 典型特性 |
|---|---|---|---|---|---|
| OST | 一次选 $K$ 个 | 通常需要 | 无迭代 | 快速支撑估计 | 成本低、可能受强系数干扰 |
| OMP | 每次选 1 个 | 通常需要 | 否 | 逐步稀疏逼近 | 直观、可解释、早期误选不可逆 |
| CoSaMP | 每轮选约 $2K$ 个候选 | 需要 | 是 | 迭代修正 $K$-稀疏估计 | 可剪枝、需要较大的 LS 子问题 |
| BP | 全部坐标参与凸规划 | 不要求先知 $K$ | 不适用 | 约束下最小 $\ell_1$ | 理论成熟，内存/求解开销需考虑 |
| BPDN / LASSO | 全部坐标，软阈值等 | 不必先知 $K$ | 不适用 | 噪声下稀疏估计 | 可正则化、需要调整 $\varepsilon$ 或 $\lambda$ |
工程建议。 目标 $K$ 已知、维度适中时，OMP 是清晰的入门基线；候选较大、希望允许纠错时尝试 CoSaMP；需要明确噪声模型或支持稳定凸优化时考虑 BPDN/LASSO；OST 适合作为成本很低的首轮筛选。最终应以同一数据集上的验证误差、恢复率、时间和内存占用比较，而不是仅凭理论充分条件作排名。
12. 实现细节：最小二乘、停止条件与复杂度
12.1 不要在代码里显式求逆
数学上可写 $\hat s_S=(A_S^A_S)^{-1}A_S^y$，但数值上推荐 QR、SVD 或 numpy.linalg.lstsq，因为直接形成 $A_S^*A_S$ 会平方条件数，使病态问题更严重。
12.2 规范化
对使用相关性选择原子的算法，应该归一化 $A$ 的列。否则大范数的列天然获得更大内积，相关性分数不可比较。需要记录归一化缩放，并在最后把系数还原至原始物理单位。
12.3 停止条件
已知 $K$ 时可以迭代至 $K$ 个指标；有噪声时可根据已知噪声尺度停止，例如 $\lVert r_t\rVert_2\le\varepsilon$；也可结合相邻迭代变化、交叉验证和最大迭代次数。不能在含噪数据上强制要求残差为零。
12.4 恢复效果不能只看残差
若 $A$ 有相似或重复的列，错误支撑集也可能使 $\lVert y-A\hat s\rVert_2$ 很小。应区分：
\[E_{\rm meas}=\frac{\lVert y-A\hat s\rVert_2}{\lVert y\rVert_2},\qquad E_{\rm coef}=\frac{\lVert s-\hat s\rVert_2}{\lVert s\rVert_2}.\]如果知道真实 $S$，再计算支撑集精确恢复率、Precision、Recall。物理状态重建还应检查 $\lVert x-\hat x\rVert_2$。
12.5 信号与硬件的现实因素
理论上的随机高斯矩阵可能要求每个测量都综合大量原始位置，物理成本高。非均匀时间采样可能要求较高的时间戳精度。MRI 的欠采样发生于 $k$ 空间，测量设计还受扫描轨迹、运动、噪声与硬件限制。采样数量少不必然意味着总体系统成本低；重建计算、稳定性和设备成本需要一起衡量。
课程模块｜第四部分：Python 动手实验（七个相互衔接的实验）
13. 环境配置与复现约定
本教程的七个实验只依赖 NumPy 与 SciPy：
1
2
python -m pip install numpy scipy
python compressed_sensing_labs.py
配套脚本为 compressed_sensing_labs.py，提供函数 basis_pursuit、ista、ost、omp、cosamp、least_absolute_deviations 和 rpca_ialm，并包含七个实验的入口函数。所用随机数种子已固定，可在常见 Python 环境直接运行。
实现规定：数组从 0 开始索引；OMP 与 CoSaMP 的相关性计算采用共轭转置，以兼容复数傅里叶字典；最小二乘使用 numpy.linalg.lstsq。七个实验的数字是这组代码与参数的可复现输出，不是算法优劣的普遍定理。
14. 实验一与实验二：L2、BP、OST、OMP、CoSaMP
14.1 实验一：所有方程都拟合，却未必恢复正确系数
取 $P=128$、$M=48$、$K=6$。先构造已知的 6-稀疏真实系数 $s_0$，再生成高斯测量矩阵 $A$，最后设置 $y=As_0$。
对比两种求解：
\[s_2=A^\dagger y,\qquad s_1=\arg\min_s\lVert s\rVert_1\quad\text{s.t.}\quad As=y.\]在配套脚本 experiment_1() 的固定实验中：
| 指标 | 最小 L2 范数解 | BP（最小 L1 范数解） |
|---|---|---|
| 系数相对误差 | 0.731562 | 接近 0 |
| 非零个数（阈值 $10^{-6}$） | 128 | 6 |
结论： 两种方法都可以拟合观测，但 L2 的最小范数解通常是密集的。BP 在这个满足良好测量条件的构造实例中恢复了真实稀疏系数。
请注意一个与原 Brunton 字幕相关的细节：如果随便生成随机 $y$，而没有先定义某个真实稀疏 $s_0$，只能够比较“L1 解更稀疏”，不能声称它恢复了某个真实信号。这里我们显式构造了地面真值，避免这一逻辑混淆。
14.2 实验二：三种贪婪策略的对照
取 $M=64$、$P=192$、$K=8$，独立重复 25 组高斯测量与稀疏系数。每组分别使用 OST、OMP、CoSaMP 恢复。设“成功”为系数相对误差小于 $10^{-4}$。
一次固定种子的运行结果：
| 算法 | 25 组中的成功率 | 平均相对误差 |
|---|---|---|
| OST | 0% | 0.5234 |
| OMP | 96% | 0.0105 |
| CoSaMP | 100% | 接近 0 |
这个场景对 OST 的一次性筛选比较困难。OMP 通常显著改善性能，CoSaMP 在这 25 组中表现最好。但是绝不能据此认定 CoSaMP 总是比 OMP 好：改变矩阵、信噪比、幅值分布和计算资源后，结果可以变化。
练习建议： 把 $K$ 从 2 增加至 20，计算不同 $K$ 下的成功率曲线；把 $M$ 从 24 增加至 100，查看不同算法的相变趋势；再用相关性较强的字典替换高斯矩阵，检验理论分析与经验性能的区别。
15. 实验三：随机时间采样为什么能识别高频双音？
15.1 对 Brunton 实验的缩小复现
Brunton 原例使用 4096 个高分辨率点和 128 个随机采样点，合成 97 Hz 与 777 Hz 的双音信号。为了使实验运行更迅速，本教程使用 $N=1024$、$M=64$，频率换为 37 Hz 和 211 Hz：
\[x(t)=\cos(2\pi\cdot37t)+0.8\cos(2\pi\cdot211t),\qquad t_n=n/1024.\]在完整 1024 点 DFT 上，非零系数位于频率索引 ${37,211,1024-211,1024-37}$。需要四个复数 DFT 原子，是因为实信号的正负频率成对出现。
分别采用：
- 随机采样：从 1024 个候选时刻随机抽取 64 个；
- 低速均匀采样：每隔 16 个原始网格点取一个，总计 64 个。
二者都记录每个测量点的确切时刻，使用对应的部分 DFT 字典与 OMP 恢复 4-稀疏频域系数。采样位置不同，矩阵 $A$ 不同，这是能否区分频率的关键。
15.2 实验结果
配套脚本 experiment_3() 的固定种子输出：
| 采样 | 找到的频域索引 | 全信号相对误差 | 测量残差 |
|---|---|---|---|
| 64 个随机点 | 37、211、813、987 | 接近 0 | 约 $6.2\times10^{-13}$ |
| 64 个均匀点 | 83、539、805、877 | 约 1.41421 | 约 $5.5\times10^{-13}$ |
均匀采样同样可以把观测点拟合得非常好，但恢复出来的完整信号是错的。这一现象正是“测量拟合好”不等于“原始信号恢复好”的生动例证。
15.3 理解何为真正的稀疏先验
如果你早已知道信号只能有 37 Hz 和 211 Hz 两个频率，那么只估计两个幅度，问题比未知支撑集容易得多。压缩感知真正要解决的是：在可能存在大量候选频率的情况下，由少量精心设计的观测自行识别活跃频率及其系数。
这里的“低于奈奎斯特”是对总观测数量或平均采样速率而言，不意味着高频信号在低通前端或粗糙时间戳下仍可无损识别。若时钟不精确、信号不在离散频率网格上、噪声很强，恢复将变困难。
15.4 实验七：二维图像只采集部分像素
为把音频的思想迁移到二维空间，配套脚本还构造了一张 $24\times24$ 的合成 DCT 稀疏图像。在二维 DCT 域仅设置 9 个非零系数；从全部 576 个像素中随机观测 150 个（约 26.0%），再用部分 DCT 测量矩阵与 OMP 重建。
固定种子运行得到相对图像误差接近 0。它演示了从“少量像素”反推“全部变换系数”和“完整图像”的过程；但这是严格稀疏的人工示例，并不代表真实照片只采样 26% 像素就一定能够无损恢复。
进阶操作： 改用真实灰度图像，将 DCT 系数截断为最大的 $K$ 项，然后比较正交小波、DCT 和不同采样掩码；分别记录测量误差、完整图像 PSNR 与 SSIM，而不是只看采样位置是否拟合。
16. Brunton 的 L1 鲁棒回归与 LASSO
16.1 同样叫 L1，却解决两种不同问题
考虑回归模型 $b=X\beta+\varepsilon$。
最小绝对偏差（LAD）回归：
\[\boxed{\min_\beta\lVert X\beta-b\rVert_1}\]惩罚的是残差，相较于平方误差，对少量大型离群值通常更稳健。
LASSO：
\[\boxed{\min_\beta\frac{1}{2n}\lVert X\beta-b\rVert_2^2+\lambda\lVert\beta\rVert_1}\]惩罚的是系数，促进变量选择与稀疏解释。由于拟合损失仍为平方损失，LASSO 并不自动具有 LAD 的抗异常值性质。
16.2 一个极简反例：均值与中位数
如果只有常数模型 $b_i\approx\beta$，L2 最小二乘的解是样本均值，L1 绝对偏差的解是样本中位数。对样本 $[1,1,1,1,100]$：
\[\hat\beta_{L2}=20.8,\qquad\hat\beta_{L1}=1.\]一个异常值即可显著拉偏均值，而中位数在这个例子中不受其影响。这是理解 Brunton 第 13–15 讲的最好入门例子。
16.3 带异常值的直线拟合
配套脚本 experiment_4() 对 $b\approx2t+1$ 生成 81 个含小高斯噪声的观测，并向其中 4 个点加入大异常值。结果为：
| 方法 | 估计斜率 | 估计截距 |
|---|---|---|
| 真实参数 | 2.0000 | 1.0000 |
| L2 最小二乘 | 1.8497 | 1.0009 |
| LAD（L1 残差） | 2.0077 | 0.9967 |
该例展示了平方损失对异常值敏感，而绝对偏差更接近多数正常数据点。
16.4 LASSO 与解释性建模
当 $X$ 的列代表年龄、温度、压力、不同多项式项或候选物理特征时，LASSO 可促使某些 $\beta_j=0$，从大量特征中选出较少的一组。它与压缩感知共享 $\ell_1$ 稀疏化，但任务通常是特征选择或参数估计，不一定是“少量采样恢复高维信号”。
实际使用时必须在训练数据上拟合标准化参数，再同样处理验证和测试数据；选择 $\lambda$ 时使用交叉验证或模型选择；面对高度共线的解释变量，系数选择可能不稳定。Brunton 第 16 讲还讨论了动力系统建模中的稀疏模型发现与 SR3 之类的扩展方法，属于这一建模分支。
16.5 同时处理稀疏信号与稀疏异常值
设 $y=As+o+\xi$，其中 $s$ 是稀疏信号，$o$ 是少量受污染的测量位置，$\xi$ 是小噪声。一个可行模型是
\[\boxed{\min_{s,o}\frac12\lVert As+o-y\rVert_2^2+\lambda_s\lVert s\rVert_1+\lambda_o\lVert o\rVert_1.}\]两个 L1 惩罚分别对应信号稀疏性和异常值稀疏性。它提供了连接压缩感知、稳健回归与稀疏错误修正的统一视角。
课程模块｜第五部分：把稀疏思想推广到工程与机器学习
17. 稀疏表示分类（SRC）与鲁棒 PCA
17.1 SRC：分类不必先恢复“完整清晰的人脸”
Brunton 第 17 讲讨论了稀疏表示分类（Sparse Representation-based Classification，SRC）。把不同类别的训练向量按列拼成字典
\[D=[D_1,D_2,\ldots,D_C],\]对一个新样本 $y$，寻找一个主要集中于少数训练原子或某一类别的稀疏系数 $\alpha$：
\[\hat\alpha=\arg\min_\alpha\lVert\alpha\rVert_1\quad\mathrm{s.t.}\quad D\alpha\approx y.\]得到系数后，按类别保留对应部分 $\delta_c(\hat\alpha)$，计算
\[r_c=\lVert y-D\delta_c(\hat\alpha)\rVert_2,\qquad\hat c=\arg\min_c r_c.\]分类依据是“哪一类原子最能解释新样本”。对于稀疏遮挡或损坏，可改为 $y=D\alpha+o$，联合估计稀疏表示系数与异常像素。
该方法依赖训练字典能否代表测试数据，以及图像对齐、光照变化等条件；演示中的遮挡鲁棒性不能无条件推广到任意图像或任意分类任务。
17.2 RPCA：低秩结构加稀疏错误
Brunton 第 18 讲讨论鲁棒主成分分析（Robust PCA）。设观测矩阵 $D\in\mathbb R^{n\times m}$ 可分解为
\[\boxed{D=L+S},\]其中 $L$ 是低秩的主要结构，$S$ 是元素稀疏的污染或异常值。直接最小化 $\operatorname{rank}(L)+\lambda\lVert S\rVert_0$ 不易求解，因此使用凸替代：
\[\boxed{\min_{L,S}\lVert L\rVert_*+\lambda\lVert S\rVert_1\quad\mathrm{s.t.}\quad D=L+S.}\]$\lVert L\rVert_*=\sum_i\sigma_i(L)$ 是核范数（奇异值之和），鼓励低秩；$\lVert S\rVert_1=\sum_{ij}\lvert S_{ij}\rvert$ 鼓励元素稀疏。这个优化问题称为 Principal Component Pursuit。
若真实信号被少量坏像素污染，普通 PCA 可能被异常值拖偏，而 RPCA 可以在适当结构假设下分离低秩背景与稀疏误差。若异常呈整行、整列、成块遮挡或低秩结构本身并不符合假设，简单的元素稀疏模型可能不合适，需改为结构化稀疏模型。
17.3 Python 实验五：拆分背景与污染
配套脚本 experiment_5() 生成 $50\times55$ 的秩 3 矩阵 $L_0$，随机选择约 5% 元素加入大幅异常值，得到 $D=L_0+S_0$。采用奇异值软阈值 + 元素软阈值的增广拉格朗日迭代分离 $L,S$。
固定随机种子的计算结果，两部分相对误差均小于 $10^{-6}$。这是人为构造、符合模型假设的实例，表明在适当条件下可做到近乎精确的分离。真实数据应同时验证低秩性、污染结构与分离结果的可解释性。
18. 稀疏传感器优化：从重建到分类
18.1 一定要随机采样吗？
Brunton 第 19–20 讲提出另一个非常有用的问题：如果工程师事先知道测量对象，例如某类流场、海面温度或翼面应变，是否能找到比随机部署更有效的固定传感器位置？
这不是简单重复传统压缩感知，而是数据驱动的定制稀疏传感：先从训练数据学习一个低维子空间，再优化在哪些位置测量。
18.2 用 SVD 建立状态模态
把训练时刻的全状态排列为
\[X=[x_1,\ldots,x_T]\in\mathbb R^{N\times T}.\]对 $X$ 做 SVD：$X=U\Sigma V^T$。截取前 $r$ 个模态 $U_r\in\mathbb R^{N\times r}$，则新状态近似为
\[x\approx U_r a.\]选择 $m$ 个位置得到采样矩阵 $C$，观测为
\[y=Cx+e\approx CU_r a+e.\]若 $m\ge r$，且 $CU_r$ 列满秩，可通过
\[\hat a=(CU_r)^\dagger y,\qquad\hat x=U_r\hat a\]重构全状态。传感器设计应让 $CU_r$ 条件良好，从而降低噪声放大。
18.3 为什么对 $U_r^T$ 做带主元 QR 分解
带列主元的 QR 分解可写成
\[U_r^T P=QR.\]置换矩阵 $P$ 对应的主元顺序，提供一个近似有效的传感器位置排序。原因是 $U_r^T$ 的每一列对应一个可测量的空间位置；主元选择会优先挑选能提供新增线性独立信息的行。
实际应用中还可以考虑噪声协方差、部署成本、可观测性约束和传感器间的最小距离。QR 主元是一种高效启发式近似，不等于对每个问题求出了全局最优部署。
18.4 Python 实验六：QR 与随机传感器比较
配套脚本 experiment_6() 采用 $N=150$ 的空间系统、$r=6$ 个模态和 $m=10$ 个传感器，在 200 个不同状态上加入测量噪声。固定的结果为：
| 部署方式 | 200 次试验的平均重建相对误差 |
|---|---|
| QR 主元选址 | 0.0357 |
| 随机位置 | 0.0663 |
在这组设置中，QR 选址在有噪声时更好地重构了低维状态。改变模态、噪声和候选位置后，结果会有所不同。
18.5 如果目标只是分类，传感器位置要重新设计
完整重建与分类不是同一个任务。假如只需要区分两种状态，那么传感器应关注类别之间差异最大的位置，而不一定重建所有变化模式。
Brunton 第 20 讲介绍的 SSPOC（Sparse Sensor Placement Optimization for Classification）利用低维特征及判别方向 $w$，通过类似
\[\min_q\lVert q\rVert_1\quad\mathrm{s.t.}\quad U_r^Tq=w\]的稀疏设计模型，寻求少量可测位置来实现分类判别。$q$ 中非零条目对应被选中的传感器。这个问题依赖已学习的判别模型与训练数据；它不同于随机测量恢复全部 $x$ 的压缩感知。
18.6 应用地图
| 应用 | 稀疏结构 | 测量设计 | 目标 |
|---|---|---|---|
| 双音音频 | 频域少数非零系数 | 随机时间点、精确时间戳 | 重建全时域信号 |
| 图像子采样 | DCT/小波/傅里叶近似稀疏 | 适当设计的像素或线性测量 | 重建图像 |
| MRI 欠采样 | 变换域稀疏或梯度稀疏 | $k$ 空间欠采样 | 加速高质量成像 |
| 稳健回归 | 大误差位置少 | 普通回归观测 | 抗离群值参数估计 |
| LASSO | 有效参数少 | 特征矩阵 | 变量选择与预测 |
| SRC | 测试样本由少量训练原子表示 | 训练字典 | 分类与稀疏遮挡处理 |
| RPCA | 低秩背景 + 稀疏异常 | 已观测数据矩阵 | 背景和异常分离 |
| 稀疏传感 | 低维状态子空间 | QR 等优化位置 | 低成本状态重建 |
| SSPOC | 低维判别方向 | 稀疏选择测量位置 | 少量传感器完成分类 |
课程模块｜第六部分：学习路径、练习与知识边界
19. 压缩感知在真实工程中何时值得用？
压缩感知的重要场景通常不是“计算机压缩文件时少占磁盘”——那里常规变换编码和熵编码已非常有效。更有价值的场景是：原始测量非常昂贵、缓慢、危险或不可重复，而计算资源相对可获得。
例如，在 MRI 中欠采样 $k$ 空间可以减少采集负担，并以稀疏约束或总变分重建弥补未测量的信息。在无线频谱监测中，若活跃频段稀疏，可用适当的压缩采样方案捕捉宽频带稀疏结构。在流场和环境监测中，若完整空间传感器部署昂贵，可以利用低维模态与合理位置选择重建特定场景。
然而，工程可行性取决于模拟前端是否能产生所需测量、采样时刻是否准确、噪声分布、模型与真实信号是否匹配、计算成本以及错误恢复的后果。与医学和安全相关的系统尤其需要真实数据验证和完整的风险评估，不能仅凭合成实验声称临床可用。
19.1 事先必须问的五个问题
- 到底稀疏在哪里？ 原始坐标、傅里叶、小波、学习字典，还是仅仅近似低秩？
- 谁决定测量算子？ 如果不能设计采样，仅能使用已有 $C$，就必须检验其可辨识性，而不能假设随机。
- 真正的误差分布是什么？ 小高斯噪声、脉冲异常值、缺失数据、量化误差、系统偏差，对应的损失函数不同。
- 需要恢复什么？ 完整信号、活跃频率、某个参数、分类结果，目标决定最优传感器和优化模型。
- 如何验证？ 仅拟合观测 $y$ 不够，还需要留出测量、独立标注、模拟真值或下游性能检验。
19.2 哪些情形不宜直接使用基础压缩感知
信号没有可靠稀疏先验；测量矩阵有重复列或过高相干性；采样存在未知偏差；字典与真实物理过程严重失配；信号变化快于测量周期；计算资源无法满足实时约束；对误恢复代价极高而缺乏外部验证——都需要格外谨慎。
20. 8 周学习安排
| 周次 | 学习主题 | 讲座对应 | 动手任务 |
|---|---|---|---|
| 1 | 结构、稀疏与变换基 | Brunton 1–2、21–22 | 对一维信号做 DFT/DCT，绘制系数衰减 |
| 2 | 采样、欠定逆问题、奈奎斯特 | Brunton 3–4、7–10 | 复现随机与均匀采样的双音实验 |
| 3 | $\ell_0/\ell_1/\ell_2$、LP、BPDN | Brunton 5–6、11 | 复现实验一，手写 BP 线性规划 |
| 4 | 相干性、spark、RIP、测量复杂度 | Brunton 12，张颢理论讨论 | 实测不同随机矩阵的相干性和恢复率 |
| 5 | OST、OMP、CoSaMP | 张颢课堂 | 实现与比较三种贪婪算法 |
| 6 | ISTA、LASSO、LAD | Brunton 13–16 | 复现稳健回归并考察异常值比例 |
| 7 | SRC 与 RPCA | Brunton 17–18 | 复现低秩+稀疏分离并解释失败案例 |
| 8 | QR 传感器与分类任务 | Brunton 19–20 | 复现 QR 布点，撰写小型实验报告 |
每周不仅要“能够跑通程序”，还应写下相应假设、矩阵维度、误差指标和一个使算法失效的反例。
21. 练习题与参考提示
21.1 基础题
题 1： 对 $s=(0,2,0,-3,0,1)^T$，求 $\lVert s\rVert_0,\lVert s\rVert_1,\lVert s\rVert_2$ 和支撑集。
提示： 分别得到 $3,6,\sqrt{14}$；支撑集为 ${2,4,6}$（从 1 开始编号）。
题 2： 设一个 $K$-稀疏信号的测量矩阵满足 $\operatorname{spark}(A)>2K$。证明在无噪声情况下该 $K$-稀疏解唯一。
提示： 假设两个不同 $K$-稀疏解，考察其差在 $A$ 的零空间中且至多有 $2K$ 个非零分量。
题 3： 某列归一化矩阵 $\mu(A)=0.05$。根据 $(2K-1)\mu<1$，最多能保证多大整数 $K$？
答案： $K\le10$。边界严格不等号不可随意改成等号。
题 4： 说明为什么 $A_S^*r=0$ 可以避免 OMP 重新选择已选原子。
提示： 已选列与当前残差的相关性理论上为零；但数值实现仍应显式排除它们。
21.2 进阶题
题 5： 重新独立推导 OST 条件 $s_{\min}/\lVert s\rVert_1>2\mu/(1+\mu)$，并解释它为何取决于真实系数幅度分布。
提示： 分别用三角不等式对真实列给出下界、对虚假列给出上界。
题 6： 对 $M=64,P=192,K=8$，使用 100 个随机实验测量 OMP 和 CoSaMP 的成功率；然后构造两列近乎相同的矩阵，观察结果变化。
提示： 计算列互相干性 $\max_{i\neq j}\lvert a_i^*a_j\rvert$，比较不同矩阵的相关性。
题 7： 对双音示例，仅将 64 个采样点改为等距，比较完整重建误差与测量残差；解释为什么它们差异很大。
提示： 研究部分 DFT 矩阵中频率别名列的关系。
题 8： 假设 $D=L+S$，但是低秩成分 $L$ 也恰好只有一个位置非零。此时能否唯一拆分 $L$ 与 $S$？
提示： 同一矩阵既低秩又稀疏时，分解可能不可辨识；RPCA 需要额外假设。
题 9： 为什么“只要分类”有时需要的传感器更少？
提示： 只要保留与决策边界有关的信息，就不必重建所有背景变化。
21.3 综合课程项目（推荐）
任选双音检测、稀疏图像重建、稀疏异常修正或传感器布点之一。提交一个可运行脚本、一张结构示意图、一组误差随测量次数/噪声变化的实验图，以及 1000–1500 字说明：测量矩阵如何形成、稀疏先验是什么、算法为什么可能成功、哪里会失败、什么情况下应使用不同方法。含失败案例的报告比只展示成功图更能证明理解。
22. 推荐阅读
Brunton 课程回答了：为什么自然数据有结构、怎样借助先验设计少量测量、何时在科学与工程中值得付出重建代价？
张颢课程回答了：一旦写成稀疏欠定方程，我们怎么从字典中找到原子、怎样使用投影和最小二乘、哪些相干性条件保证选对？
将两条路线放在一起，完整的压缩感知链条就是：
\[\text{结构先验}\longrightarrow\text{测量设计}\longrightarrow \text{欠定逆问题}\longrightarrow\text{稀疏求解}\longrightarrow \text{恢复保证与应用验证}.\]以下为两位老师相关课程之外的补充阅读。
- E. Candès, J. Romberg, T. Tao. Stable Signal Recovery from Incomplete and Inaccurate Measurements. arXiv:math/0503066. 稳定稀疏恢复与噪声误差界。
- D. Needell, J. A. Tropp. CoSaMP: Iterative Signal Recovery from Incomplete and Inaccurate Samples. 期刊 DOI. CoSaMP 原始理论。
- R. Tibshirani. Regression Shrinkage and Selection via the Lasso. DOI. LASSO 与变量选择。
- J. Wright et al. Robust Face Recognition via Sparse Representation. DOI. SRC 与遮挡处理。
- E. J. Candès et al. Robust Principal Component Analysis? DOI. 低秩加稀疏的主成分追踪。
- K. Manohar et al. Data-Driven Sparse Sensor Placement for Reconstruction. arXiv:1701.07569. SVD、QR 与定制传感器。
- B. W. Brunton et al. Sparse Sensor Placement Optimization for Classification. DOI. SSPOC 与分类导向测量。
23. 附录：数学证明补充与 Python 实验代码
阅读说明。 本附录把正文需要反复使用、但证明较长的结论集中起来。前七节提供可逐行复核的证明或推导；CoSaMP 的完整收缩定理需要一组技术引理，因此附录第 23.8 节明确区分已证明的工具与证明框架，避免将直觉误作完整定理。推导统一假设列向量已归一化；复数情形用共轭转置 $A^*$，实数情形即 $A^T$。
23.1 预备：支撑集、矩阵范数与 Gram 矩阵
设 $S\subseteq{1,\ldots,P}$，$A_S$ 是取 $A$ 中指标为 $S$ 的列组成的子矩阵。若 $A_S$ 列满秩，记
\[A_S^\dagger=(A_S^*A_S)^{-1}A_S^*.\]对矩阵 $B$，诱导 $\ell_1$ 范数和诱导 $\ell_\infty$ 范数分别是
\[\lVert B\rVert_1=\max_j\sum_i\lvert B_{ij}\rvert,\qquad\lVert B\rVert_\infty=\max_i\sum_j\lvert B_{ij}\rvert.\]并有 $\lVert Bz\rVert_1\le\lVert B\rVert_1\lVert z\rVert_1$。当 $A$ 的列归一化且 $\lvert S\rvert=K$，$G=A_S^*A_S=I+E$ 中 $E$ 的对角元为零，非对角元的模至多为 $\mu(A)$，因此
\[\lVert E\rVert_1\le(K-1)\mu(A).\]若此值小于 1，Neumann 级数 $\sum_{n=0}^\infty(-E)^n$ 在矩阵范数意义下收敛，并等于 $(I+E)^{-1}$；同时
\[\boxed{\lVert G^{-1}\rVert_1\le\frac1{1-(K-1)\mu(A)}.}\tag{A1}\]证明：有限部分和 $Q_n=\sum_{j=0}^n(-E)^j$ 满足 $(I+E)Q_n=I-(-E)^{n+1}$。令 $n\to\infty$，余项范数至多 $\lVert E\rVert_1^{n+1}\to0$，故 $Q_n\to G^{-1}$。三角不等式及几何级数直接给出 (A1)。
23.2 证明：spark、互相干性与稀疏解唯一性
命题 1： 如果 $\operatorname{spark}(A)>2K$，则对任意观测 $y$，至多存在一个满足 $As=y$ 的 $K$-稀疏解。
证明。 反设 $s,s’$ 是不同的两个 $K$-稀疏解，令 $h=s-s’\ne0$。则 $Ah=0$ 且 $\lVert h\rVert_0\le2K$，表示不超过 $2K$ 个列存在非平凡线性关系，违背 spark 的定义。证毕。
命题 2： 对单位范数列矩阵，若 $\mu(A)>0$，则
\[\boxed{\operatorname{spark}(A)\ge1+\frac1{\mu(A)}.}\tag{A2}\]证明。 设最少 $q=\operatorname{spark}(A)$ 列线性相关。取非零 $h\in\ker(A_T)$，其中 $\lvert T\rvert=q$，并选取 $i\in T$ 使 $\lvert h_i\rvert=\lVert h\rVert_\infty$。由 $A_T h=0$，左乘 $a_i^*$，得
\[0=h_i+\sum_{j\in T\setminus\{i\}}(a_i^*a_j)h_j.\]所以
\[\lvert h_i\rvert\le\sum_{j\ne i}\lvert a_i^*a_j\rvert\,\lvert h_j\rvert\le(q-1)\mu(A)\lvert h_i\rvert.\]由于 $\lvert h_i\rvert>0$，有 $1\le(q-1)\mu(A)$，即 $q\ge1+1/\mu(A)$。证毕。
推论。 若 $(2K-1)\mu(A)<1$，则 $2K<1+1/\mu(A)\le\operatorname{spark}(A)$，所以所有 $K$-稀疏解唯一。注意这一结论保证“解可区分”，并不自动保证任何具体求解方法都会找到它。
若 $\mu(A)=0$，所有列两两正交，在有限列数情形下直接依据列独立性判断，无须除以 $\mu$。
命题 3（Welch 下界）。 对 $A\in\mathbb C^{M\times P}$，$P>M$，若所有列单位归一化，则
\[\boxed{\mu(A)^2\ge\frac{P-M}{M(P-1)}.}\tag{A2b}\]证明。 记 $G=A^A$。由单位范数列，$\operatorname{tr}(G)=P$。设 $AA^$ 的 $M$ 个非负特征值为 $\lambda_1,\ldots,\lambda_M$（不足 $M$ 个正特征值时用零补齐）。由柯西—施瓦茨不等式，
\[\lVert G\rVert_F^2=\lVert AA^*\rVert_F^2=\sum_{i=1}^M\lambda_i^2\ge\frac{(\sum_i\lambda_i)^2}{M}=\frac{P^2}{M}.\]另一方面，$G_{ii}=1$，且每个非对角元素的模至多为 $\mu(A)$，故
\[\lVert G\rVert_F^2=\sum_{i,j}\lvert a_i^*a_j\rvert^2\le P+P(P-1)\mu(A)^2.\]合并两式，整理得到 (A2b)。证毕。这说明在固定测量维数和字典大小时，互相干性存在不可突破的几何下界。
23.3 证明：OST 的恢复充分条件及含噪声版本
设 $y=As+e$、$S=\operatorname{supp}(s)$、$\lvert S\rvert=K$、$s_{\min}=\min_{i\in S}\lvert s_i\rvert$，定义噪声引起的最大相关性扰动
\[\tau=\lVert A^*e\rVert_\infty.\]真实列的相关性下界。 对 $i\in S$：
\[\begin{aligned} \lvert a_i^*y\rvert&=\left\lvert s_i+\sum_{j\in S\setminus\{i\}}s_j a_i^*a_j+a_i^*e\right\rvert\\ &\ge \lvert s_i\rvert-\mu\sum_{j\ne i}\lvert s_j\rvert-\tau\\ &\ge (1+\mu)s_{\min}-\mu\lVert s\rVert_1-\tau. \end{aligned}\]虚假列的相关性上界。 对 $j\notin S$：
\[\lvert a_j^*y\rvert=\left\lvert\sum_{i\in S}s_i a_j^*a_i+a_j^*e\right\rvert\le\mu\lVert s\rVert_1+\tau.\]若前者严格大于后者，则 Top-$K$ 相关性必然来自真实支撑集，故只需
\[(1+\mu)s_{\min}-\mu\lVert s\rVert_1-\tau>\mu\lVert s\rVert_1+\tau.\]整理得到含噪声 OST 的充分条件：
\[\boxed{(1+\mu)s_{\min}>2\mu\lVert s\rVert_1+2\tau.}\tag{A3}\]无噪声时 $\tau=0$，即正文第 8.2 节的公式。若仅知道 $\lVert e\rVert_2\le\varepsilon$，因为 $\lVert a_i\rVert_2=1$，柯西—施瓦茨不等式给出 $\tau\le\varepsilon$，因此把 (A3) 右端的 $2\tau$ 换成 $2\varepsilon$ 是一个更保守、但便于检查的条件。
若 $K$ 个非零系数等幅，$s_{\min}=\lVert s\rVert_1/K$，无噪声条件即 $(2K-1)\mu<1$。证明的关键是相关性间隔，不是“相关性大的列必然正确”这一未经量化的直觉。
23.4 证明：OMP 的 ERC、Neumann 级数与逐轮选对
设 $S$ 为真实支撑集，$\lvert S\rvert=K$，$A_S$ 列满秩且 $y\in\operatorname{span}(A_S)$。若当前选中集合 $T\subseteq S$，OMP 通过最小二乘得到
\[r=(I-P_T)y,\qquad P_T=A_TA_T^\dagger.\]从而 $A_T^*r=0$；又因为 $y\in\operatorname{span}(A_S)$ 且 $P_Ty\in\operatorname{span}(A_S)$，所以 $r\in\operatorname{span}(A_S)$。
引理（错误相关性由真实相关性控制）。 对任意 $r\in\operatorname{span}(A_S)$，都有
\[\boxed{\lVert A_{S^c}^*r\rVert_\infty\le \gamma_S\lVert A_S^*r\rVert_\infty,\quad \gamma_S=\max_{j\notin S}\lVert A_S^\dagger a_j\rVert_1.}\tag{A4}\]证明。 由于 $A_S$ 满列秩，$r=A_S(A_S^\dagger r)$；对任意 $j\notin S$，
\[\begin{aligned} \lvert a_j^*r\rvert &=\lvert a_j^*A_S(A_S^*A_S)^{-1}A_S^*r\rvert\\ &=\lvert(A_S^\dagger a_j)^*A_S^*r\rvert\\ &\le\lVert A_S^\dagger a_j\rVert_1\lVert A_S^*r\rVert_\infty. \end{aligned}\]对 $j\notin S$ 取最大值即得 (A4)。
归纳证明 OMP 选对。 设 $\gamma_S<1$。初始 $T=\varnothing$，$r=y$。若 $r\ne0$，由于 $A_S$ 满列秩，$A_S^*r$ 不会全为零；结合 (A4) 得到
\[\max_{j\notin S}\lvert a_j^*r\rvert<\max_{i\in S}\lvert a_i^*r\rvert.\]而已经选中的指标 $i\in T$ 满足 $a_i^*r=0$，所以下一轮最大相关性必出现在 $S\setminus T$。于是每次都增加一个真实指标。由于 $A_S$ 列满秩、真实 $s_i$ 均非零，当 $T\subsetneq S$ 时 $y$ 不属于 $\operatorname{span}(A_T)$，故 $r\ne0$；迭代 $K$ 轮后 $T=S$，最小二乘解就是 $s_S$。证毕。
从相干性得到 ERC。 根据附录第 23.1 节的 (A1) 和 $\lVert A_S^*a_j\rVert_1\le K\mu$，
\[\begin{aligned} \lVert A_S^\dagger a_j\rVert_1 &=\lVert(A_S^*A_S)^{-1}A_S^*a_j\rVert_1\\ &\le\frac{K\mu}{1-(K-1)\mu}. \end{aligned}\]当 $(2K-1)\mu<1$ 时，分母为正且该上界严格小于 1，因此 $\gamma_S<1$，ERC 成立，OMP 精确恢复。这里同时解释了正文里“Neumann 级数可用”的条件并非无条件成立。
23.5 证明：BP 精确恢复的零空间性质（NSP）
假设 $A\in\mathbb R^{M\times P}$ 或 $\mathbb C^{M\times P}$。对给定整数 $K$，称 $A$ 满足 $K$ 阶零空间性质，若对每个非零 $h\in\ker A$ 以及任意 $\lvert S\rvert\le K$ 都有
\[\boxed{\lVert h_S\rVert_1<\lVert h_{S^c}\rVert_1.}\tag{NSP}\]定理。 BP 对每一个至多 $K$-稀疏向量 $s_0$ 都能唯一精确恢复，当且仅当 $A$ 满足 (NSP)。
充分性。 设 $s_0$ 是任意 $K$-稀疏向量，$S=\operatorname{supp}(s_0)$。任意其他可行解可写成 $s_0+h$，其中 $0\ne h\in\ker A$。按支撑分块并利用反三角不等式：
\[\begin{aligned} \lVert s_0+h\rVert_1 &=\lVert(s_0+h)_S\rVert_1+\lVert h_{S^c}\rVert_1\\ &\ge\lVert s_0\rVert_1-\lVert h_S\rVert_1+\lVert h_{S^c}\rVert_1\\ &>\lVert s_0\rVert_1. \end{aligned}\]因此 $s_0$ 是严格唯一的最小 $\ell_1$ 范数可行解。
必要性。 反设 NSP 不成立，则存在某非零 $h\in\ker A$ 与 $\lvert S\rvert\le K$，使 $\lVert h_S\rVert_1\ge\lVert h_{S^c}\rVert_1$。取 $s_0=-h_S$，它是 $K$-稀疏的；而 $s_1=s_0+h=h_{S^c}$ 是另一个可行解，且
\[As_1=As_0,\qquad\lVert s_1\rVert_1=\lVert h_{S^c}\rVert_1\le\lVert h_S\rVert_1=\lVert s_0\rVert_1.\]故 BP 不能唯一恢复 $s_0$，矛盾。证毕。
理解。 这一命题把“BP 何时真的等价于寻找最稀疏表示”转化为零空间的几何不等式。它并不要求知道真实支撑集，但直接检查 NSP 仍往往很困难，因此实际理论常进一步使用 RIP 等可分析的矩阵性质。
23.6 证明：BPDN 在 RIP 下的稳定恢复误差界
这里证明一个经典的、常数不追求最优的版本。设
\[y=As+e,\qquad\lVert e\rVert_2\le\varepsilon,\]$\hat s$ 是 BPDN 问题 $\min_z\lVert z\rVert_1$、满足 $\lVert Az-y\rVert_2\le\varepsilon$ 的任一最优解。设 $\delta=\delta_{2K}<\sqrt2-1$；令 $S$ 是 $s$ 的 $K$ 个最大绝对值坐标的集合，$\sigma=\lVert s_{S^c}\rVert_1=\sigma_K(s)_1$，$h=\hat s-s$。
23.6.1 管约束（tube constraint）
因为真实 $s$ 与求得的 $\hat s$ 均满足 BPDN 约束，
\[\boxed{\lVert Ah\rVert_2\le\lVert A\hat s-y\rVert_2+\lVert As-y\rVert_2\le2\varepsilon.}\tag{A5}\]23.6.2 锥约束（cone constraint）
由 $\lVert\hat s\rVert_1\le\lVert s\rVert_1$ 及分块三角不等式，
\[\begin{aligned} \lVert s\rVert_1&\ge\lVert s+h\rVert_1\\ &\ge\lVert s_S\rVert_1-\lVert h_S\rVert_1+\lVert h_{S^c}\rVert_1-\lVert s_{S^c}\rVert_1, \end{aligned}\]整理得到
\[\boxed{\lVert h_{S^c}\rVert_1\le\lVert h_S\rVert_1+2\sigma.}\tag{A6}\]23.6.3 依幅度对误差坐标分块
把 $h_{S^c}$ 按绝对值递减顺序分成每组最多 $K$ 个坐标的 $S_1,S_2,\ldots$，令 $T=S\cup S_1$，则 $\lvert T\rvert\le2K$。排序保证对 $j\ge2$，
\[\lVert h_{S_j}\rVert_2\le\sqrt K\,\lVert h_{S_j}\rVert_\infty\le\frac{\lVert h_{S_{j-1}}\rVert_1}{\sqrt K}.\]对所有 $j\ge2$ 求和可得
\[\boxed{\sum_{j\ge2}\lVert h_{S_j}\rVert_2\le\frac{\lVert h_{S^c}\rVert_1}{\sqrt K}.}\tag{A7}\]23.6.4 RIP 的受限正交性
若 $u,v$ 支撑不相交且各自至多有 $K$ 个非零分量，由 RIP 和极化恒等式可得
\[\lvert\langle Au,Av\rangle\rvert\le\delta_{2K}\lVert u\rVert_2\lVert v\rVert_2.\]这个结论可以通过把 $u/\lVert u\rVert_2$ 与 $v/\lVert v\rVert_2$ 分别加减、应用两次 RIP 再相减证明。把 $h_T=h_S+h_{S_1}$ 拆成两块，与每个 $h_{S_j}$（$j\ge2$）配对，再用柯西—施瓦茨不等式，得到
\[\lvert\langle Ah_T,Ah_{S_j}\rangle\rvert\le\sqrt2\,\delta\,\lVert h_T\rVert_2\lVert h_{S_j}\rVert_2.\]23.6.5 主要误差块的上界
注意 $Ah_T=Ah-\sum_{j\ge2}Ah_{S_j}$。结合 RIP 下界、(A5)、(A7) 和上式：
\[\begin{aligned} (1-\delta)\lVert h_T\rVert_2^2 &\le\lVert Ah_T\rVert_2^2\\ &=\left\lvert\langle Ah_T,Ah\rangle-\sum_{j\ge2}\langle Ah_T,Ah_{S_j}\rangle\right\rvert\\ &\le2\varepsilon\sqrt{1+\delta}\,\lVert h_T\rVert_2 +\sqrt2\,\delta\,\lVert h_T\rVert_2\frac{\lVert h_{S^c}\rVert_1}{\sqrt K}. \end{aligned}\]若 $h_T=0$，后面的上界直接成立。否则除以 $\lVert h_T\rVert_2$，再由 (A6) 与 $\lVert h_S\rVert_1\le\sqrt K\lVert h_T\rVert_2$ 得
\[\bigl[1-(1+\sqrt2)\delta\bigr]\lVert h_T\rVert_2 \le2\sqrt{1+\delta}\,\varepsilon +2\sqrt2\,\delta\frac{\sigma}{\sqrt K}.\]因为 $\delta<\sqrt2-1=1/(1+\sqrt2)$，系数
\[D_\delta=1-(1+\sqrt2)\delta>0,\]所以
\[\boxed{\lVert h_T\rVert_2\le\frac{2\sqrt{1+\delta}}{D_\delta}\varepsilon+\frac{2\sqrt2\,\delta}{D_\delta}\frac{\sigma}{\sqrt K}.}\tag{A8}\]23.6.6 总体误差界
由三角不等式、(A7) 和 (A6)，
\[\begin{aligned} \lVert h\rVert_2 &\le\lVert h_T\rVert_2+\sum_{j\ge2}\lVert h_{S_j}\rVert_2\\ &\le\lVert h_T\rVert_2+\frac{\lVert h_S\rVert_1+2\sigma}{\sqrt K}\\ &\le2\lVert h_T\rVert_2+\frac{2\sigma}{\sqrt K}. \end{aligned}\]代入 (A8)，得到所需稳定恢复界：
\[\boxed{\lVert\hat s-s\rVert_2\le \underbrace{\frac{4\sqrt{1+\delta}}{D_\delta}}_{C_0}\varepsilon +\underbrace{\left(2+\frac{4\sqrt2\,\delta}{D_\delta}\right)}_{C_1}\frac{\sigma_K(s)_1}{\sqrt K}.}\tag{A9}\]因此，严格 $K$-稀疏且无噪声时误差为零；近似稀疏或含噪声时，误差由噪声能量和最佳 $K$ 项截断误差控制。这就是正文第 6.5 节结论的具体证明。
23.7 证明：软阈值与 ISTA 更新式的来历
先考虑标量问题
\[\min_z\left\{\frac12(z-v)^2+\tau\lvert z\rvert\right\},\quad\tau>0.\]当 $z>0$ 时，一阶条件是 $z-v+\tau=0$，故 $z=v-\tau$ 且必须满足 $v>\tau$；当 $z<0$ 时一阶条件是 $z-v-\tau=0$，故 $z=v+\tau$ 且必须满足 $v<-\tau$；在 $z=0$ 时次梯度条件为 $0\in-v+\tau[-1,1]$，等价于 $\lvert v\rvert\le\tau$。三种情况合并，得
\[\operatorname{prox}_{\tau\lvert\cdot\rvert}(v)=\operatorname{sign}(v)\max(\lvert v\rvert-\tau,0).\]向量 $\ell_1$ 惩罚是坐标可分的，因此向量近端映射就是逐坐标软阈值。
对 $F(s)=f(s)+\lambda\lVert s\rVert_1$，$f(s)=\tfrac12\lVert As-y\rVert_2^2$，有 $\nabla f(s)=A^*(As-y)$，且梯度的 Lipschitz 常数为 $L=\lVert A\rVert_2^2$。对任何 $u$，下降引理给出
\[f(u)\le f(s^t)+\langle\nabla f(s^t),u-s^t\rangle+\frac L2\lVert u-s^t\rVert_2^2.\]以右端加上 $\lambda\lVert u\rVert_1$ 为上界模型求最小值，并配方：
\[\begin{aligned} s^{t+1}&=\arg\min_u\left\{\frac L2\left\lVert u-\left(s^t-\frac1L\nabla f(s^t)\right)\right\rVert_2^2+\lambda\lVert u\rVert_1\right\}\\ &=\operatorname{soft}\left(s^t-\frac1L A^*(As^t-y),\frac\lambda L\right). \end{aligned}\]这证明了 ISTA 更新式。该算法求解的是带 $\ell_1$ 惩罚的复合凸目标，而不是通过一步选择支撑集来求解；两类算法应在理论与实验中分别讨论。
23.8 CoSaMP 的证明框架：为什么扩展、最小二乘和剪枝缺一不可
CoSaMP 的完整收缩定理涉及多轮不等式和明确的 RIP 常数。这里列出能独立验证的基础工具，再解释定理证明的逻辑，不将下列步骤冒充为完整的收缩定理证明。
23.8.1 RIP 控制残差代理量
若 $z$ 是支撑不超过 $K$ 的误差向量，代理量 $u=A^*Az$ 在相关支撑上的分量与 $z$ 接近：对任意集合 $T$，且 $\lvert T\cup\operatorname{supp}(z)\rvert\le q$，
\[\boxed{\lVert(A^*Az-z)_T\rVert_2\le\delta_q\lVert z\rVert_2.}\]证明：将差写为 $A_T^A_Uz_U-z_T$，其中 $U=\operatorname{supp}(z)$。它正是子 Gram 矩阵 $A_{T\cup U}^A_{T\cup U}-I$ 的一个坐标块，其谱范数不超过 $\delta_q$。
23.8.2 最小二乘提供最佳局部拟合
若 $b_T=\arg\min_v\lVert y-A_Tv\rVert_2$，则对任何支撑包含于 $T$ 的候选 $w$，有 $\lVert y-A_Tb_T\rVert_2\le\lVert y-A_Tw_T\rVert_2$。正规方程还保证残差与 $\operatorname{span}(A_T)$ 正交。
23.8.3 Top-$K$ 剪枝的误差界
对真实 $K$-稀疏 $s$、任意近似 $b$，令 $H_K(b)$ 保留 $b$ 中最大 $K$ 个幅值，则
\[\boxed{\lVert s-H_K(b)\rVert_2\le2\lVert s-b\rVert_2.}\]证明：由于 $H_K(b)$ 是 $b$ 的最佳 $K$ 项欧氏近似，$\lVert b-H_K(b)\rVert_2\le\lVert b-s\rVert_2$；再使用三角不等式即可。
23.8.4 收缩证明的结构
在足够强的 RIP 下，首先证明残差相关性排序能捕捉真实支撑集中尚未解释的大部分能量；然后利用最小二乘和 RIP 控制候选集合上的系数误差；接着用上面的剪枝界，保证回到 $K$-稀疏模型不会使误差无限放大；最后合成形如
\[\lVert s-s^{t+1}\rVert_2\le\rho\lVert s-s^t\rVert_2+c\lVert e\rVert_2,\qquad0<\rho<1\]的迭代误差关系。仅靠工具 1–3 还不足以证明特定 $\rho$ 和 RIP 阈值；要得到论文中的严格常数，必须加入候选识别与局部拟合误差的完整技术引理。建议在掌握本教程附录第 23.6 节后，对照 Needell–Tropp 原论文逐条研究这些引理。
23.9 证明：低速均匀采样为何产生频率混叠
以正文第 15 章实验为例，原始离散网格为 $N=1024$，只选择 $M=64$ 个等间隔时刻 $n=16m$（$m=0,\ldots,63$）。第 $k$ 个离散傅里叶原子在这些位置上的复数相位为
\[\psi_k(16m)=\frac1{\sqrt{1024}}\exp\left(2\pi i\frac{k\,16m}{1024}\right)=\frac1{\sqrt{1024}}\exp\left(2\pi i\frac{km}{64}\right).\]对于任何整数 $q$，
\[\psi_{k+64q}(16m)=\frac1{\sqrt{1024}}\exp\left(2\pi i\frac{(k+64q)m}{64}\right)=\psi_k(16m),\]因为 $e^{2\pi iqm}=1$。因此有效恢复矩阵 $A=C\Psi$ 的第 $k$ 和第 $k+64q$ 列相同，它们的相干性达到 1、spark 为 2，甚至单频支撑都可能不可辨识。
随机选取时刻会打破这种固定周期采样的列重复结构，因此对稀疏频率识别通常有利；但随机性只是建立较好矩阵性质的一种机制，不意味着任意随机采样都必然成功。
23.10 证明概要：高斯随机矩阵为何只需 $M=O(K\log(P/K))$ 个测量
这里证明一个固定 RIP 精度 $\delta\in(0,1)$ 和失败概率 $\eta\in(0,1)$ 下的基本量级结论，重点说明概率集中、有限网与并集界如何结合。完整尾界的最佳常数不是本教程的目标。
定理（高斯 RIP 的数量级）。 令 $A\in\mathbb R^{M\times P}$ 的元素独立服从 $\mathcal N(0,1/M)$。存在普适常数 $C>0$，使得当
\[\boxed{M\ge C\delta^{-2}\left(K\log\frac{eP}{K}+\log\frac2\eta\right),}\tag{A10}\]时，$A$ 以至少 $1-\eta$ 的概率满足 $K$ 阶 RIP，常数不超过 $\delta$。这里 $1\le K\le P$；要求 $2K$ 阶 RIP 时，原则上把上式中的 $K$ 换成 $2K$。
23.10.1 固定向量的长度集中
给定任意固定的单位向量 $z\in\mathbb R^P$，由于高斯分布的旋转不变性，$Az$ 的每个坐标独立服从 $\mathcal N(0,1/M)$，故
\[\lVert Az\rVert_2^2\overset{d}=\frac1M\sum_{i=1}^M Z_i^2,\quad Z_i\overset{\mathrm{i.i.d.}}\sim\mathcal N(0,1).\]对 $t\in(0,1)$，由卡方分布的矩母函数和 Chernoff 上界，存在普适常数 $c>0$，使
\[\boxed{\Pr\bigl(\lvert\lVert Az\rVert_2^2-1\rvert>t\bigr)\le2e^{-cMt^2}.}\tag{A11}\]其中矩母函数可直接由 $\mathbb E e^{u(Z_i^2-1)}=e^{-u}(1-2u)^{-1/2}$（$u<1/2$）推导，因此 (A11) 是标准概率集中结论，而不是随意作出的“随机矩阵通常很好”的假设。
23.10.2 在固定支撑集上构造有限网
固定一个大小为 $K$ 的支撑集 $S$。其单位球面上的一个 $1/4$-网 $\mathcal N_S$ 可选择为不超过 $9^K$ 个点，即对任意支撑位于 $S$ 的单位向量 $u$，都存在 $v\in\mathcal N_S$ 使 $\lVert u-v\rVert_2\le1/4$。这是球体体积比较给出的标准覆盖数上界。
把 (A11) 中的阈值选为 $t=\delta/2$，对全部网点使用并集界，得到
\[\Pr\left(\exists v\in\mathcal N_S:\ \lvert\lVert Av\rVert_2^2-1\rvert>\delta/2\right)\le2\cdot9^K e^{-cM\delta^2/4}.\]记 $B_S=A_S^A_S-I$，它是实对称矩阵。若所有网点满足 $\lvert v^B_Sv\rvert\le\delta/2$，取达到 $\lVert B_S\rVert_2$ 的单位向量 $u$，选择距离不超过 $1/4$ 的网点 $v$，可得
\[\lvert u^*B_Su\rvert\le\lvert v^*B_Sv\rvert+2\lVert u-v\rVert_2\lVert B_S\rVert_2\le\delta/2+\frac12\lVert B_S\rVert_2.\]因此 $\lVert B_S\rVert_2\le\delta$。这就把有限网上的测量长度保证推广到了固定支撑集上的全部向量。
23.10.3 对所有可能支撑集取并集界
大小为 $K$ 的支撑集数量为 $\binom PK\le(eP/K)^K$。再次应用并集界，出现至少一个不符合 RIP 的支撑集的概率至多
\[2\binom PK9^K e^{-cM\delta^2/4}\le2\exp\left(K\log\frac{9eP}{K}-\frac{cM\delta^2}{4}\right).\]要求该概率不超过 $\eta$，只需 $M$ 大于某个普适常数倍的 $\delta^{-2}[K\log(eP/K)+\log(2/\eta)]$。这证明了 (A10) 的量级。最后，RIP 必须作用于所有至多 $K$-稀疏向量，因此“固定向量集中”本身远远不够；从固定向量到全部稀疏向量，需要上述两次并集控制。
适用边界。 这个证明针对独立高斯测量矩阵；结构化部分傅里叶、随机卷积、确定性矩阵、带依赖的物理传感器，并不能直接套用同一个独立性和覆盖数论证，额外对数因子与硬件约束也可能出现。
23.11 Python 实验代码
以下代码就是配套文件 compressed_sensing_labs.py 的完整内容。复制后可直接运行；需要安装 NumPy 和 SciPy。使用实数 BP 与 LAD 时采用线性规划，含复数频谱的双音实验使用支持复数共轭内积的 OMP。
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""压缩感知与稀疏建模：配套可复现实验。
依赖：numpy, scipy。运行：python compressed_sensing_labs.py
所有随机实验均使用固定随机种子。数值结果取决于算法参数与具体样本，
不构成算法普遍优劣的证明。
"""
import numpy as np
from scipy.optimize import linprog
from scipy.linalg import qr
def topk(z, k):
"""返回绝对值最大的 k 个索引。"""
k = min(max(int(k), 0), z.size)
return np.argsort(np.abs(z), kind="stable")[-k:] if k else np.array([], dtype=int)
def refit(A, y, support):
"""在给定支撑集上执行最小二乘，其他系数为零。"""
s = np.zeros(A.shape[1], dtype=np.result_type(A.dtype, y.dtype))
if len(support):
s[support] = np.linalg.lstsq(A[:, support], y, rcond=None)[0]
return s
def ost(A, y, k):
"""One-step thresholding + 支撑集内最小二乘重估。"""
return refit(A, y, topk(A.conj().T @ y, k))
def omp(A, y, k, tol=1e-10):
"""Orthogonal matching pursuit，假设各列具有相同范数。"""
residual = y.copy()
support = []
estimate = np.zeros(A.shape[1], dtype=np.result_type(A.dtype, y.dtype))
for _ in range(min(k, A.shape[1])):
if np.linalg.norm(residual) < tol:
break
scores = np.abs(A.conj().T @ residual)
scores[support] = -np.inf
j = int(np.argmax(scores))
if not np.isfinite(scores[j]):
break
support.append(j)
estimate = refit(A, y, np.array(support))
residual = y - A @ estimate
return estimate
def cosamp(A, y, k, iters=15, tol=1e-10):
"""CoSaMP 变体：每次剪枝后附加一次最小二乘消偏。"""
n = A.shape[1]
estimate = np.zeros(n, dtype=np.result_type(A.dtype, y.dtype))
residual = y.copy()
for _ in range(iters):
if np.linalg.norm(residual) < tol:
break
candidates = topk(A.conj().T @ residual, min(2 * k, n))
union = np.union1d(np.flatnonzero(estimate), candidates)
b = refit(A, y, union)
support = topk(b, k)
new_estimate = refit(A, y, support)
new_residual = y - A @ new_estimate
if (np.array_equal(np.flatnonzero(estimate), np.flatnonzero(new_estimate))
and np.linalg.norm(new_residual-residual) < tol):
estimate = new_estimate
break
estimate, residual = new_estimate, new_residual
return estimate
def basis_pursuit(A, y):
"""实数 BP: min ||s||_1 subject to A s=y，线性规划。"""
if np.iscomplexobj(A) or np.iscomplexobj(y):
raise ValueError("本示例只实现了实数 BP")
m, n = A.shape
eye = np.eye(n)
c = np.r_[np.zeros(n), np.ones(n)]
A_ub = np.vstack((np.hstack((eye, -eye)), np.hstack((-eye, -eye))))
b_ub = np.zeros(2*n)
A_eq = np.hstack((A, np.zeros((m, n))))
bounds = [(None, None)]*n + [(0, None)]*n
result = linprog(c, A_ub=A_ub, b_ub=b_ub,
A_eq=A_eq, b_eq=y, bounds=bounds, method="highs")
if not result.success:
raise RuntimeError("BP 求解失败: " + result.message)
return result.x[:n]
def soft_threshold(x, lam):
return np.sign(x) * np.maximum(np.abs(x)-lam, 0)
def ista(A, y, lam, iterations=800):
"""实值 LASSO：min .5||As-y||_2^2 + lam||s||_1。"""
lipschitz = np.linalg.norm(A, 2)**2
s = np.zeros(A.shape[1], dtype=float)
for _ in range(iterations):
z = s - A.T @ (A @ s-y)/lipschitz
s = soft_threshold(z, lam/lipschitz)
return s
def least_absolute_deviations(D, y):
"""LAD: min_beta ||D beta-y||_1；beta 不受符号约束。"""
m, p = D.shape
c = np.r_[np.zeros(p), np.ones(2*m)]
A_eq = np.hstack((D, np.eye(m), -np.eye(m)))
bounds = [(None, None)]*p + [(0, None)]*(2*m)
result = linprog(c, A_eq=A_eq, b_eq=y, bounds=bounds, method="highs")
if not result.success:
raise RuntimeError("LAD 求解失败: " + result.message)
return result.x[:p]
def rpca_ialm(M, max_iter=350, tol=1e-7):
"""教学版主成分追踪（Principal Component Pursuit）的 IALM 迭代。"""
lam = 1/np.sqrt(max(M.shape))
sigma_max = np.linalg.norm(M, 2)
if sigma_max == 0:
return np.zeros_like(M), np.zeros_like(M)
dual_norm = max(sigma_max, np.max(np.abs(M))/lam)
Y = M / dual_norm
mu = 1.25 / sigma_max
mu_max = 1e7 * mu
L = np.zeros_like(M)
S = np.zeros_like(M)
for _ in range(max_iter):
U, singular, Vh = np.linalg.svd(M-S+Y/mu, full_matrices=False)
rank = np.sum(singular > 1/mu)
if rank:
L = (U[:, :rank] * (singular[:rank]-1/mu)) @ Vh[:rank, :]
else:
L = np.zeros_like(M)
S = soft_threshold(M-L+Y/mu, lam/mu)
Z = M-L-S
Y += mu*Z
if np.linalg.norm(Z, "fro") / np.linalg.norm(M, "fro") < tol:
break
mu = min(mu*1.5, mu_max)
return L, S
def relative_error(a, b):
return float(np.linalg.norm(a-b)/np.linalg.norm(a))
def experiment_1():
print("\n[实验1] 最小 L2 解与 BP 的差异")
rng = np.random.default_rng(13)
m, n, k = 48, 128, 6
A = rng.normal(size=(m, n))
A /= np.linalg.norm(A, axis=0)
truth = np.zeros(n)
idx = rng.choice(n, k, replace=False)
truth[idx] = rng.choice([-1, 1], k)*rng.uniform(0.7, 1.5, k)
y = A@truth
l2 = np.linalg.pinv(A)@y
bp = basis_pursuit(A,y)
print(" 真实支撑:", sorted(idx.tolist()))
print(" L2 误差:", round(relative_error(truth, l2),6),
"L2 非零分量数(>1e-6):",int(np.sum(np.abs(l2)>1e-6)))
print(" BP 误差:", round(relative_error(truth, bp),6),
"BP 非零分量数(>1e-6):",int(np.sum(np.abs(bp)>1e-6)))
def experiment_2():
print("\n[实验2] OST / OMP / CoSaMP: 多次随机试验")
rng = np.random.default_rng(21)
m, n, k = 64, 192, 8
methods = {"OST":ost,"OMP":omp,"CoSaMP":cosamp}
scores = {name: [] for name in methods}
n_trials = 25
for _ in range(n_trials):
A = rng.standard_normal((m, n))
A /= np.linalg.norm(A, axis=0)
truth = np.zeros(n)
inds = rng.choice(n,k,replace=False)
truth[inds] = rng.choice([-1,1], k)*rng.uniform(.6,1.4,k)
y = A@truth
for name,algo in methods.items():
recovered = algo(A,y,k)
scores[name].append(relative_error(truth,recovered))
for name,vals in scores.items():
vals=np.asarray(vals)
print(f" {name:7s}: 恢复率(error<1e-4)={np.mean(vals<1e-4):.0%}, "
f"平均相对误差={np.mean(vals):.4f}")
def experiment_3():
print("\n[实验3] 双频信号：随机采样 vs 低速均匀采样")
# 参考 Brunton 的 97 Hz / 777 Hz 例子；此处使用更小规模的 37 Hz / 211 Hz
N, M = 1024, 64
t=np.arange(N)/N
truth=np.cos(2*np.pi*37*t)+0.8*np.cos(2*np.pi*211*t)
rng=np.random.default_rng(2)
groups={"随机":np.sort(rng.choice(N,M,replace=False)),
"均匀":np.arange(0,N,N//M)}
kvec=np.arange(N)
for label,indices in groups.items():
# 已知不均匀采样位置时，使用部分 DFT 测量矩阵
Phi=np.exp(2j*np.pi*np.outer(indices,kvec)/N)/np.sqrt(N)
estimate=omp(Phi,truth[indices].astype(complex),k=4,tol=1e-8)
reconstructed=np.fft.ifft(estimate)*np.sqrt(N)
dominant=sorted(topk(estimate,4).tolist())
print(f" {label}: 频域索引={dominant}, "
f"全信号误差={relative_error(truth,reconstructed):.5f}, "
f"测量残差={np.linalg.norm(truth[indices]-reconstructed[indices]):.3g}")
def experiment_4():
print("\n[实验4] 鲁棒回归：L2 残差 vs L1 残差")
rng=np.random.default_rng(8)
t=np.linspace(-2,2,81)
D=np.column_stack((t,np.ones(len(t))))
y=2*t+1+rng.normal(scale=.08,size=len(t))
y[[4,17,55,73]] += np.array([16,-18,15,-13])
beta_l2=np.linalg.lstsq(D,y,rcond=None)[0]
beta_lad=least_absolute_deviations(D,y)
print(" 真参数: [2, 1]")
print(" L2参数:",np.round(beta_l2,4))
print(" LAD参数:",np.round(beta_lad,4))
def experiment_5():
print("\n[实验5] 鲁棒 PCA：低秩背景与稀疏污染")
rng=np.random.default_rng(10)
n, p, r = 50, 55, 3
U=rng.standard_normal((n,r))/np.sqrt(r)
V=rng.standard_normal((r,p))/np.sqrt(r)
L0=U@V
S0=np.zeros((n,p))
mask=rng.random((n,p))<.05
S0[mask]=rng.uniform(-8,8,np.count_nonzero(mask))
M=L0+S0
L,S=rpca_ialm(M)
print(" 低秩部分相对误差:",round(relative_error(L0,L),6))
print(" 稀疏异常值部分相对误差:",round(relative_error(S0,S),6))
def experiment_6():
print("\n[实验6] 稀疏传感器：QR 主元选址与随机选址")
rng=np.random.default_rng(15)
n, r, m = 150, 6, 10
modes,_ = np.linalg.qr(rng.standard_normal((n,r)))
q,rmat,pivots=qr(modes.T,pivoting=True,mode="economic")
idx_qr=pivots[:m]
idx_random=rng.choice(n,m,replace=False)
errors={"QR优化":[],"随机":[]}
for _ in range(200):
alpha=rng.standard_normal(r)
truth=modes@alpha
for key,idx in [("QR优化",idx_qr),("随机",idx_random)]:
observed=truth[idx]+rng.normal(scale=.01,size=m)
estimate=modes@np.linalg.lstsq(modes[idx,:],observed,rcond=None)[0]
errors[key].append(relative_error(truth,estimate))
for key,vals in errors.items():
print(f" {key}: 平均误差={np.mean(vals):.4f}")
def experiment_7():
print("\n[实验7] 二维图像：仅凭部分像素恢复 DCT 稀疏图像")
from scipy.fft import idctn
rng = np.random.default_rng(31)
side, m = 24, 150
n = side*side
coefficients = np.zeros((side,side))
locations = [(0,0),(0,1),(1,0),(1,1),(2,0),
(0,2),(2,2),(3,0),(1,3)]
amplitudes = [5.,-2.,1.5,1.2,-1.1,.9,.8,-.7,.6]
for ij, amp in zip(locations, amplitudes):
coefficients[ij] = amp
truth = idctn(coefficients, norm="ortho")
# 每一个 DCT 单位系数对应一张基图像；以列排列成合成字典
dictionary = idctn(
np.eye(n).reshape(n,side,side), axes=(1,2), norm="ortho"
).reshape(n,n).T
idx = np.sort(rng.choice(n,m,replace=False))
A = dictionary[idx,:]
y = truth.ravel()[idx]
column_norms = np.linalg.norm(A,axis=0)
normalized_A = A/column_norms
s_scaled = omp(normalized_A,y,k=len(locations))
recovered = dictionary @ (s_scaled/column_norms)
print(" 像素采样率:",f"{m/n:.1%}")
print(" 图像相对误差:",round(relative_error(truth.ravel(),recovered),6))
if __name__ == "__main__":
experiment_1()
experiment_2()
experiment_3()
experiment_4()
experiment_5()
experiment_6()
experiment_7()
阅读结束后的检验标准：能够解释公式里的每一个矩阵从哪里来；能够写出 OMP 与 CoSaMP 的伪代码；知道“L1 惩罚系数”和“L1 惩罚残差”的区别；能够指出在高相干、噪声、模型失配和均匀混叠场景下恢复为什么会失败；能够通过可复现实验验证或证伪自己的判断。
