Robust PCA：低秩 + 稀疏分解

Art. 4 范数与条件数 建立了 nuclear 范数 $\|A\|_* = \sum_i \sigma_i$ 和 $\ell_1$ 范数 $\|A\|_1 = \sum_{ij} |a_{ij}|$ 两件度量工具，Art. 8 矩阵补全 展示了 nuclear 范数最小化如何从部分观测恢复低秩矩阵。现在我们面对一个不同但同样深刻的问题：如果观测是完整的，但部分元素被大幅度异常值污染了呢？

标准 PCA（Art. 6 PCA）通过 SVD 截断实现降维，隐含假设噪声是小幅度、均匀分布的高斯扰动。但现实数据中的”噪声”往往不均匀——监控视频中行人遮挡了背景，基因表达数据中少数基因发生了大幅度异常表达，传感器矩阵中个别传感器彻底失效。这些异常值虽然稀疏（只影响少数元素），却有大幅度——足以摧毁 SVD 截断的效果。

Robust PCA（鲁棒主成分分析）正是为此设计的：将数据矩阵精确分解为低秩部分 $L$ 和稀疏部分 $S$：

$M = L + S$

低秩 $L$ 捕获数据的全局结构（“背景”），稀疏 $S$ 捕获异常值（“前景”）。

Candès、Li、Ma 和 Wright (2011) 证明了一个优美的定理：在温和条件下，这个分解可以通过一个凸优化问题精确恢复——只需最小化 nuclear 范数与 $\ell_1$ 范数的加权和。 更惊人的是，平衡参数有一个”万能”的选择：$\lambda = 1/\sqrt{\max(m,n)}$——不依赖于秩 $r$ 或稀疏度的先验知识。

问题设定：数据 = 低秩 + 稀疏异常

从标准 PCA 的脆弱性说起

标准 PCA 的数据模型是：

$M = L_0 + N$

其中 $L_0 \in \mathbb{R}^{m \times n}$ 是低秩信号矩阵（$\text{rank}(L_0) = r \ll \min(m,n)$），$N$ 是小幅度噪声（通常假设 i.i.d. 高斯，每个元素 $\sim \mathcal{N}(0, \sigma^2)$）。在这个模型下，SVD 截断到前 $r$ 个奇异值是 $L_0$ 的最佳估计（Eckart-Young 定理，Art. 3 SVD）。

问题：如果少数元素被大幅度污染——哪怕只有一个元素被改为极大值——SVD 截断的效果可能急剧恶化。直觉上，一个极大的异常值会”吸引”主奇异方向偏向它，导致主成分不再反映真实的低维结构。

Robust PCA 的数据模型

Robust PCA 用不同的模型描述数据：

$M = L_0 + S_0 \qquad \text{(1)}$

其中：

• $L_0 \in \mathbb{R}^{m \times n}$：低秩矩阵，$\text{rank}(L_0) = r$，捕获数据的全局低维结构
• $S_0 \in \mathbb{R}^{m \times n}$：稀疏矩阵，大部分元素为零，但非零元素的幅度可以任意大

与标准 PCA 的关键区别：

	标准 PCA	Robust PCA
噪声模型	小幅度、稠密 ($N$ 每个元素都有噪声)	大幅度、稀疏 ($S_0$ 只有少数非零元素)
噪声幅度	$\|N\|_F / \sqrt{mn} \sim \sigma$（小）	$\|S_0\|_\infty$ 可以任意大
方法	SVD 截断	凸优化（nuclear + $\ell_1$）

核心洞察：标准 PCA 假设每个元素都被少量污染（稠密小噪声），Robust PCA 假设少数元素被大量污染（稀疏大异常）。两者是噪声模型的两个极端，对应不同的数学工具。

典型应用：视频背景/前景分离

Robust PCA 最经典的应用是视频监控中的背景/前景分离。将视频的每一帧展平为一个列向量，所有帧排成矩阵 $M \in \mathbb{R}^{p \times t}$（$p$ 是像素数，$t$ 是帧数）：

• 背景（低秩 $L_0$）：静态或缓慢变化的场景。由于背景在各帧间高度相关，背景矩阵的秩很低（几个主成分就能描述）
• 前景（稀疏 $S_0$）：移动的行人、车辆等。在任何一帧中，前景只占画面的小部分（稀疏），但在前景像素上的变化幅度很大

将 $M$ 分解为 $L_0 + S_0$ 就是把静止背景和移动前景分离——不需要任何运动模型或物体检测器。

Principal Component Pursuit

凸优化公式

恢复 $(L_0, S_0)$ 最自然的想法是：

$\min_{L, S} \; \text{rank}(L) + \gamma \|S\|_0 \quad \text{s.t.} \quad M = L + S$

其中 $\|S\|_0$ 是 $S$ 的非零元素个数。但 $\text{rank}(\cdot)$ 和 $\|\cdot\|_0$ 都是非凸、不连续的——这个问题是 NP-hard 的。

Art. 4 范数与条件数 的核心结果提供了解决方案：

• $\text{rank}(L)$ 的凸松弛 → nuclear 范数 $\|L\|_* = \sum_i \sigma_i(L)$
• $\|S\|_0$ 的凸松弛 → $\ell_1$ 范数 $\|S\|_1 = \sum_{ij} |s_{ij}|$

替换后得到 Principal Component Pursuit (PCP)（也称为 Robust PCA 的凸优化形式）：

$\min_{L, S} \; \|L\|_* + \lambda \|S\|_1 \quad \text{s.t.} \quad M = L + S \qquad \text{(2)}$

这是一个凸优化问题：目标函数是两个凸函数之和，约束是仿射等式。

逐项理解：

• $\|L\|_* = \sum_i \sigma_i(L)$：nuclear 范数惩罚 $L$ 的奇异值之和，鼓励 $L$ 低秩——类比 $\ell_1$ 范数鼓励向量稀疏
• $\|S\|_1 = \sum_{ij} |s_{ij}|$：$\ell_1$ 范数惩罚 $S$ 的元素绝对值之和，鼓励 $S$ 稀疏——大部分元素被压到零
• $\lambda > 0$：平衡参数，控制”低秩”和”稀疏”之间的权重
• $M = L + S$：硬约束，要求分解是精确的

为什么 nuclear + $\ell_1$ 是正确的组合

PCP 的凸优化公式同时运用了两个层次的凸松弛：

矩阵层次：$\|L\|_*$ 是 $\text{rank}(L)$ 在 spectral 范数单位球 $\{L : \|L\|_2 \leq 1\}$ 上的凸包（Candès & Recht, 2009）。这意味着在所有凸函数中，nuclear 范数是秩函数的最紧下界——最小化 $\|L\|_*$ 是鼓励低秩的最有效凸方法。

元素层次：$\|S\|_1$ 是 $\|S\|_0$ 在 $\ell_\infty$ 单位球 $\{S : \|S\|_\infty \leq 1\}$ 上的凸包。这是压缩感知中 $\ell_1$ 松弛的经典结果——最小化 $\|S\|_1$ 是鼓励元素稀疏的最有效凸方法。

两个松弛各自独立地是”最紧”的。PCP 将它们组合在同一个目标函数中，同时从奇异值和元素两个角度施加凸正则化。

Robust PCA：低秩 + 稀疏分解

M = L + S — 从混合信号中分离低秩结构和稀疏异常

M（原始 = L + S）

L（低秩恢复，秩 2）

S（稀疏恢复）

-17.3+17.3

左：原始混合矩阵。中：低秩成分（背景）。右：稀疏成分（异常值）。注意 L 的结构平滑，S 只有少数大幅度非零元素。

主定理：精确恢复的理论保证

Incoherence 条件

与矩阵补全类似，精确恢复需要低秩部分 $L_0$ 满足 incoherence 条件。设 $L_0$ 的紧 SVD 为 $L_0 = U \Sigma V^T$，其中 $U \in \mathbb{R}^{m \times r}$，$V \in \mathbb{R}^{n \times r}$。

定义（Incoherence）。矩阵 $L_0$ 满足参数为 $\mu$ 的 incoherence 条件，如果：

$\max_{1 \leq i \leq m} \|U^T \mathbf{e}_i\|_2^2 \leq \frac{\mu r}{m}, \quad \max_{1 \leq j \leq n} \|V^T \mathbf{e}_j\|_2^2 \leq \frac{\mu r}{n} \qquad \text{(3)}$

$\|UV^T\|_\infty \leq \sqrt{\frac{\mu r}{mn}} \qquad \text{(4)}$

直觉与矩阵补全中相同：incoherence 要求低秩矩阵的”能量”不集中在少数行列上，而是均匀分散。如果 $L_0$ 的某一行集中了大部分能量，那么该行上的大值可能被误判为稀疏异常值——分解变得不可能。

稀疏支撑的随机模型

对稀疏矩阵 $S_0$，Candès 等人假设其支撑集（非零元素的位置）$\Omega = \text{supp}(S_0)$ 是随机的：每个元素 $(i,j)$ 独立地以概率 $\rho$ 属于 $\Omega$。非零元素的值可以任意——没有幅度限制。

这个随机支撑假设排除了”对抗性”的稀疏模式。例如，如果 $S_0$ 的非零元素恰好集中在 $L_0$ 的某一列上，那么无法区分 $S_0$ 的贡献和 $L_0$ 中一个额外的秩一成分——分解在信息论上不可能。

Candès, Li, Ma, Wright (2011) 主定理

以下是 Robust PCA 领域最核心的理论结果，严格陈述自 Candès 等人的 Theorem 1.1。

定理（Candès, Li, Ma, Wright, 2011）。设 $M = L_0 + S_0 \in \mathbb{R}^{m \times n}$，其中 $L_0$ 是秩为 $r$ 的矩阵，满足 incoherence 条件 (3)–(4)，参数为 $\mu$。设 $S_0$ 的支撑集中每个元素 $(i,j)$ 独立地以概率 $\rho$ 被选入。设 $n_{(1)} = \max(m, n)$，$n_{(2)} = \min(m, n)$。

取正则化参数

$\lambda = \frac{1}{\sqrt{n_{(1)}}} \qquad \text{(5)}$

则存在数值常数 $c > 0$，使得当

$\text{rank}(L_0) \leq \frac{c \, n_{(2)}}{\mu^2 (\log n_{(1)})^2}, \quad \rho \leq c$

时，Principal Component Pursuit (2) 的唯一解恰好是 $(L_0, S_0)$——即 $\hat{L} = L_0$，$\hat{S} = S_0$——的概率至少为 $1 - c \, n_{(1)}^{-10}$。

定理的逐项解读

$\lambda = 1 / \sqrt{n_{(1)}}$ 的选择。这是定理最令人惊叹的部分：正则化参数只依赖矩阵的维度，不需要知道秩 $r$ 或稀疏度 $\rho$。直觉来源于量纲分析：

• $\|L_0\|_*$ 的量级约为 $\sqrt{n_{(1)} n_{(2)}} \cdot r$（nuclear 范数是奇异值之和，每个奇异值量级 $\sim \sqrt{n_{(1)} n_{(2)}}$，共 $r$ 个——粗略近似）
• $\|S_0\|_1$ 的量级约为 $\rho \cdot m \cdot n \cdot |s|$（$\rho mn$ 个非零元素，每个幅度 $|s|$）
• 要使两项在同一量级上平衡，$\lambda$ 应使 $\lambda \|S_0\|_1 \sim \|L_0\|_*$，粗略地得到 $\lambda \sim 1/\sqrt{n_{(1)}}$

严格证明需要精细的对偶分析，但量纲分析给出了正确的缩放。

秩的上界 $r \leq c \, n_{(2)} / (\mu^2 (\log n_{(1)})^2)$。低秩部分的秩不能太大。当 $\mu = O(1)$（理想 incoherent 情况），秩可以高达 $O(n_{(2)} / \log^2 n_{(1)})$——接近矩阵维度的量级（对数因子除外）。

稀疏度的上界 $\rho \leq c$。常数 $c$ 意味着可以有常数比例的元素是异常值——比如 10% 或 20% 的元素被任意大幅度的异常值污染，仍然可以精确恢复。

精确恢复。$\hat{L} = L_0$，$\hat{S} = S_0$——没有任何误差。这不是渐近结果，是精确等式。

概率 $\geq 1 - c \, n_{(1)}^{-10}$。对于 $n_{(1)} = 1000$，失败概率小于 $10^{-30}$——在任何实际意义上都是确定性的。

与矩阵补全的对比

Robust PCA 和矩阵补全 (Art. 8) 解决的是同一枚硬币的两面：

	矩阵补全	Robust PCA
已知	部分元素 $M_{ij}, (i,j) \in \Omega$	全部元素 $M = L_0 + S_0$
未知	缺失元素	$L_0$ 和 $S_0$ 的分离
核心工具	$\min \|X\|_*$ s.t. 采样约束	$\min \|L\|_* + \lambda\|S\|_1$ s.t. $M = L + S$
随机性	在观测位置	在异常值位置
恢复条件	$m \geq C \mu^2 n r \log^2 n$	$r \leq c \, n / (\mu^2 \log^2 n)$，$\rho \leq c$

两者都基于 nuclear 范数的凸松弛（Art. 4 范数与条件数），都需要 incoherence 条件，都能实现精确恢复。但矩阵补全缺少元素，Robust PCA 多了异常值——一个是”不够”，一个是”太多”。

数值例子

用一个小例子手算体会 PCP 的工作原理。取 $m = n = 4$，$r = 1$：

$L_0 = \begin{bmatrix} 1 \\ 2 \\ 1 \\ 2 \end{bmatrix} \begin{bmatrix} 1 & 2 & 1 & 2 \end{bmatrix} = \begin{bmatrix} 1 & 2 & 1 & 2 \\ 2 & 4 & 2 & 4 \\ 1 & 2 & 1 & 2 \\ 2 & 4 & 2 & 4 \end{bmatrix}$

稀疏矩阵（两个大幅度异常值）：

$S_0 = \begin{bmatrix} 0 & 0 & 0 & 0 \\ 0 & 0 & 10 & 0 \\ 0 & 0 & 0 & 0 \\ -8 & 0 & 0 & 0 \end{bmatrix}$

观测矩阵：

$M = L_0 + S_0 = \begin{bmatrix} 1 & 2 & 1 & 2 \\ 2 & 4 & 12 & 4 \\ 1 & 2 & 1 & 2 \\ -6 & 4 & 2 & 4 \end{bmatrix}$

标准 PCA 的失败：对 $M$ 做 SVD 截断到秩 1，得到的主奇异向量会被 $M_{23} = 12$ 和 $M_{41} = -6$ 严重扭曲——主方向不再是 $(1,2,1,2)^T$，而是偏向异常值的方向。截断后的秩 1 近似既不等于 $L_0$，也无法给出 $S_0$。

PCP 的成功：按照定理取 $\lambda = 1/\sqrt{4} = 0.5$，求解问题 (2)。凸优化器会找到同时使 $\|L\|_*$ 和 $\|S\|_1$ 最小的分解——由于 $L_0$ 是秩 1 的（nuclear 范数很小），$S_0$ 只有 2 个非零元素（$\ell_1$ 范数为 18），任何偏离 $(L_0, S_0)$ 的分解都会增加目标函数值。在满足 incoherence 条件时，PCP 精确恢复 $(L_0, S_0)$。

ALM/ADMM 求解算法

PCP 是凸问题，但目标函数中的 nuclear 范数不可微（需要 SVD 才能计算），直接用通用凸优化器（如内点法）对大规模问题不实用。实践中最常用的是 ALM（Augmented Lagrangian Method） 和 ADMM（Alternating Direction Method of Multipliers）。

增广拉格朗日方法（ALM）

PCP 的约束形式 (2) 的增广拉格朗日函数为：

$\mathcal{L}(L, S, Y, \mu) = \|L\|_* + \lambda\|S\|_1 + \langle Y, M - L - S \rangle + \frac{\mu}{2}\|M - L - S\|_F^2$

其中 $Y \in \mathbb{R}^{m \times n}$ 是拉格朗日乘子（对偶变量），$\mu > 0$ 是惩罚参数。

逐项理解：

• $\|L\|_* + \lambda\|S\|_1$：原始目标函数
• $\langle Y, M - L - S \rangle = \text{tr}(Y^T(M - L - S))$：拉格朗日项，$Y$ “惩罚”约束的违反
• $\frac{\mu}{2}\|M - L - S\|_F^2$：二次惩罚项，进一步”拉紧”约束——这是”增广”的含义

ALM 交替更新 $(L, S)$ 和乘子 $Y$：

1. 更新 $L$（固定 $S$, $Y$）：

$L^{(k+1)} = \arg\min_L \|L\|_* + \frac{\mu}{2}\left\|L - \left(M - S^{(k)} + \frac{Y^{(k)}}{\mu}\right)\right\|_F^2$

这是 nuclear 范数的近端算子（proximal operator），解为奇异值软阈值化（Singular Value Thresholding, SVT）：

$L^{(k+1)} = \text{SVT}_{1/\mu}\left(M - S^{(k)} + \frac{Y^{(k)}}{\mu}\right)$

其中 $\text{SVT}_\tau(X) = U \, \text{diag}((\sigma_i - \tau)_+) \, V^T$，$X = U \, \text{diag}(\sigma_i) \, V^T$ 是 SVD。

2. 更新 $S$（固定 $L$, $Y$）：

$S^{(k+1)} = \arg\min_S \lambda\|S\|_1 + \frac{\mu}{2}\left\|S - \left(M - L^{(k+1)} + \frac{Y^{(k)}}{\mu}\right)\right\|_F^2$

这是 $\ell_1$ 范数的近端算子，解为元素级软阈值化：

$S^{(k+1)} = \mathcal{S}_{\lambda/\mu}\left(M - L^{(k+1)} + \frac{Y^{(k)}}{\mu}\right)$

其中 $\mathcal{S}_\tau(x) = \text{sign}(x) \cdot \max(|x| - \tau, 0)$ 逐元素作用。

3. 更新乘子 $Y$：

$Y^{(k+1)} = Y^{(k)} + \mu\left(M - L^{(k+1)} - S^{(k+1)}\right)$

4. 更新惩罚参数：$\mu \leftarrow \min(\rho \mu, \bar{\mu})$，其中 $\rho > 1$ 是增长因子。

ADMM 视角

上述 ALM 迭代可以理解为 ADMM（Boyd 等, 2011）的一个实例。ADMM 是求解如下分裂凸问题的通用框架：

$\min \; f(x) + g(z) \quad \text{s.t.} \quad Ax + Bz = c$

在 PCP 中，$f(L) = \|L\|_*$，$g(S) = \lambda\|S\|_1$，约束为 $L + S = M$。ADMM 的收敛性由 Boyd 等 (2011) 保证：对任意初始值，迭代收敛到全局最优解。

计算瓶颈与加速

ALM/ADMM 每次迭代的主要计算量在 SVT 步骤的 SVD 计算。对 $m \times n$ 矩阵做完整 SVD 的复杂度是 $O(\min(m,n) \cdot mn)$。

加速策略：

• 部分 SVD：由于 $L$ 的秩 $r \ll \min(m,n)$，只需计算前 $r + \delta$ 个奇异值/奇异向量（使用 Lanczos 算法或随机化 SVD，见 Art. 7 随机化 SVD），复杂度降至 $O(rmn)$
• 热启动：用上一次迭代的秩作为下一次部分 SVD 的目标秩
• 不精确 ALM：Lin 等 (2010) 提出的 IALM（Inexact ALM）放松子问题的精度要求，每步只做一次 SVT 更新（而不是迭代到子问题收敛），实践中收敛速度几乎相同但每步计算量大幅减少

下图直观展示 ALM 每步迭代的两个核心操作：左边是更新 $S$ 时的元素级软阈值化（绝对值小于 $\tau$ 的元素被压到零），右边是更新 $L$ 时的奇异值软阈值化 SVT（小奇异值被压到零，推动 $L$ 低秩）。

算法伪代码

输入：观测矩阵 M ∈ R^{m×n}, 参数 λ = 1/√max(m,n)
初始化：L⁰ = S⁰ = Y⁰ = 0, μ₀ > 0, ρ > 1

重复直到收敛：
  1. L^{k+1} = SVT_{1/μ}(M - Sᵏ + Yᵏ/μ)     // 奇异值软阈值化
  2. S^{k+1} = S_{λ/μ}(M - L^{k+1} + Yᵏ/μ)    // 元素级软阈值化
  3. Y^{k+1} = Yᵏ + μ(M - L^{k+1} - S^{k+1})   // 乘子更新
  4. μ ← min(ρμ, μ̄)                             // 惩罚增长

输出：L^{k+1} (低秩), S^{k+1} (稀疏)

收敛准则通常取 $\|M - L^{(k)} - S^{(k)}\|_F / \|M\|_F < \epsilon$（如 $\epsilon = 10^{-7}$）。

含噪声情形：Stable PCP

实际数据往往同时包含稀疏异常值和稠密小噪声：

$M = L_0 + S_0 + Z$

其中 $Z$ 是稠密噪声（如 $\|Z\|_F \leq \delta$）。Zhou、Li、Wright、Candès 和 Ma (2010) 提出了 Stable Principal Component Pursuit：

$\min_{L, S} \; \|L\|_* + \lambda\|S\|_1 \quad \text{s.t.} \quad \|M - L - S\|_F \leq \delta \qquad \text{(6)}$

或等价的无约束（拉格朗日）形式：

$\min_{L, S} \; \frac{1}{2}\|M - L - S\|_F^2 + \lambda_1 \|L\|_* + \lambda_2 \|S\|_1 \qquad \text{(7)}$

Zhou 等 (2010) 证明了稳定性结果：当 $\delta$ 足够小时，估计误差 $\|\hat{L} - L_0\|_F$ 与 $\delta$ 成比例——算法对小幅度噪声是稳定的。

方法谱系中的位置

Robust PCA 在 Part 1 “拆” 方法谱系中的位置：

$\text{SVD（无约束）} \xrightarrow{+\text{采样约束}} \text{矩阵补全} \xrightarrow{+\text{稀疏异常}} \text{Robust PCA}$

三者共享同一个核心工具——nuclear 范数作为秩的凸松弛（Art. 4 范数与条件数）：

• SVD 截断：给定完整矩阵，找最佳低秩近似
• 矩阵补全：给定部分元素，用 $\min \|X\|_*$ 恢复低秩矩阵
• Robust PCA：给定完整但被污染的矩阵，用 $\min \|L\|_* + \lambda\|S\|_1$ 分离低秩和稀疏

从 SVD → 矩阵补全 → Robust PCA，约束逐步增加，但核心数学工具——nuclear 范数的凸松弛——保持不变。这正是 Art. 4 范数与条件数 中 nuclear 范数理论的威力：一个凸松弛工具，三种不同的应用。

总结与展望

本文深入探讨了 Robust PCA——低秩 + 稀疏分解的理论与算法。回顾关键要点：

• 问题：$M = L_0 + S_0$，从混合信号中精确分离低秩部分和稀疏异常值
• Principal Component Pursuit：$\min_{L,S} \|L\|_* + \lambda\|S\|_1$ s.t. $M = L + S$——nuclear 范数鼓励低秩，$\ell_1$ 范数鼓励稀疏，两者的组合是凸问题
• $\lambda = 1/\sqrt{\max(m,n)}$ 的万能选择：不需要秩或稀疏度的先验知识
• Candès 等 (2011) 主定理：在 incoherence 条件和随机稀疏支撑下，PCP 以压倒性概率精确恢复 $(L_0, S_0)$
• ALM/ADMM 算法：交替执行奇异值软阈值化（更新 $L$）和元素级软阈值化（更新 $S$），收敛到全局最优
• 稳定性：Stable PCP 处理同时存在稀疏异常值和稠密噪声的情形，估计误差与噪声水平成比例

Robust PCA 是 nuclear 范数凸松弛理论的第二个重要应用（第一个是矩阵补全）。它将 Art. 4 范数与条件数 中”nuclear 范数 = 秩的凸松弛”和”$\ell_1$ 范数 = 稀疏的凸松弛”这两个独立工具，首次组合在同一个优化目标中——这种”凸松弛的叠加”思想在后续的结构化信号恢复问题中反复出现。

下一篇我们将矩阵分解推广到更高维度——张量分解：当数据不是二维表格而是三维或更高阶的”数据立方体”时，CP 分解和 Tucker 分解提供了自然的框架。