入门：The Augmented Lagrange Multiplier Method for Exact Recovery of Corrupted Low-Rank Matrices辅助阅读+总结

精确恢复损坏低阶矩阵的增强拉格朗日乘子法

摘要

本文提出了可扩展且快速的算法来解决鲁棒 PCA 问题，即恢复一个低秩矩阵，其条目的未知部分被任意损坏。这个问题在许多应用中都会出现，例如图像处理、网络数据排名和生物信息数据分析。

最近的研究表明，在令人惊讶的广泛条件下，鲁棒 PCA 问题可以通过最小化核范数和 l1-范数的组合的凸优化来精确解决。

在本文中，我们应用增强拉格朗日乘子（ALM）的方法来求解这个凸规划。由于目标函数是非光滑的，我们展示了如何将 ALM 的经典分析扩展到此类新的目标函数，并证明所提出算法的最优性并表征其收敛速度。根据经验，所提出的新算法比之前最先进的鲁棒 PCA 算法（例如加速近端梯度（APG）算法）快五倍以上。此外，新算法实现了更高的精度，但对存储/内存的要求较低。我们还表明，ALM 技术可用于解决（相关但稍微简单的）矩阵完成问题，并获得相当有希望的结果。所讨论的所有算法的 Matlab 代码均可在 http://perception.csl.illinois.edu/matrix-rank/home.html 上找到。

1 引言

主成分分析（PCA）作为一种流行的高维数据处理、分析、压缩和可视化工具，在科学和工程领域有着广泛的应用[13]。它假设给定的高维数据位于低维线性子空间附近。在很大程度上，PCA的目标是高效、准确地估计这个低维子空间。

假设给定的数据被排列为大矩阵 D ∈ Rm×n 的列。估计低维子空间的数学模型是找到一个低秩矩阵A，使得A和D之间的差异最小化，从而导致以下约束优化：

其中 r ≪ min(m, n) 是子空间的目标维度， ‖ · ‖F 是 Frobenius 范数，其对应于假设数据被 i.i.d的 高斯噪声损坏。这个问题可以通过先计算 D 的奇异值分解（SVD），然后将 D 的列向量投影到 D 的 r 个主左奇异向量所组成的子空间上来得到方便地解决[13]。

当损失是由附加的i.i.d 高斯噪声引起时，PCA 会给出最佳估计。只要噪声的幅度很小，它在实践中就可以很好地发挥作用。然而，它在严重损失的情况下就会崩溃，即使这种损失只影响极少数的观测结果。事实上，即使只有 A 的一项被任意破坏，经典 PCA 得到的估计 ^ A 也可能与真实 A 相差任意远。因此，有必要研究低秩矩阵 A 是否仍然可以有效且从损坏的数据矩阵 D = A + E 中准确恢复，其中附加误差 E 的某些条目可能是任意大。

最近，赖特等人。 [22]表明，在相当广泛的条件下，答案是肯定的：只要误差矩阵 E 足够稀疏（相对于 A 的秩），通过求解以下凸优化问题就可以从 D = A + E 精确地恢复低秩矩阵 A：

其中 ‖ · ‖* 表示矩阵的核范数（即其奇异值之和）， ‖ · ‖1 表示矩阵条目的绝对值之和， λ 是正权重参数。由于即使存在较大错误或异常值，也能够准确恢复数据中的底层低秩结构，因此这种优化在[22]中被称为鲁棒 PCA (RPCA)（一个流行术语，已被一系列旨在使 PCA 对异常值和严重损失具有鲁棒性的工作广泛使用）。 RPCA 的几个应用，例如背景建模以及从面部图像中去除阴影和镜面反射，已在[23]中进行了演示，以显示 RPCA 的优势。

优化 (2) 可以被视为一般的凸优化问题，并在被重新表述为半定规划 [10] 后，可以由任何现成的内点求解器（“interior point solver”）（例如 CVX [12]）求解。然而，虽然内点方法通常需要很少的迭代来收敛，但它们在处理大型矩阵时遇到困难，因为计算步长方向的复杂度为 O(m6)，其中 m 是矩阵的维度。因此，在典型的个人计算机 (PC) 上，通用内点求解器无法处理维度大于 m = 102 的矩阵。相比之下，图像和视频处理中的应用通常涉及维度 m = 104 到 105 的矩阵；网络搜索和生物信息学中的应用很容易涉及维度 m = 106 及以上的矩阵。因此，通用内点求解器对于鲁棒 PCA 来说过于有限，无法在许多实际应用中实用。

内点求解器对于大型矩阵不能很好地扩展，因为它们依赖于目标函数的二阶信息。为了克服可扩展性问题，我们应该仅使用一阶信息并充分利用此类凸优化问题的特殊性质。例如，最近的研究表明，（一阶）迭代阈值（IT）算法对于压缩感知中出现的 l1 范数最小化问题非常有效 [24,4,25,8]。 [7] 中还表明，相同的技术可用于最小化矩阵完成 (MC) 问题的核范数，即从不完整但干净的条目子集恢复低秩矩阵 [18, 9 ]。

由于矩阵恢复（鲁棒PCA）问题（2）涉及最小化l1范数和核范数的组合，在初始论文[22]中，作者还采用了迭代阈值技术来解决（2）和获得了类似的收敛性和可扩展性。然而，[22]中提出的迭代阈值方案收敛速度非常慢。通常，它需要大约 104 次迭代才能收敛，每次迭代的成本与一次 SVD 相同。因此，即使对于维度小至 m = 800 的矩阵，该算法也必须在典型 PC 上运行 8 小时。为了缓解迭代阈值方法[22]收敛缓慢的问题，Lin 等人 [15]提出了两种解决问题（2）的新算法，它们在某种意义上是相互补充的：第一个是应用于原始的加速近端梯度（APG）算法，它是[4] 引入的FISTA框架的直接应用，加上快速延续技术1；第二个是应用于问题 (2) 的对偶问题的梯度上升算法。从维度高达 m = 1, 000 的矩阵的模拟来看，这两种方法都比迭代阈值方法至少快 50 倍（更多详细信息请参阅 [15]）。

在本文中，我们提出了利用增强拉格朗日乘子（ALM）技术的矩阵恢复新算法。这里提出的精确ALM（EALM）方法被证明具有令人满意的Q线性收敛速度，而APG在理论上只是次线性的。略微改进精确 ALM 导致了不精确 ALM (IALM) 方法，其收敛速度实际上与精确 ALM 一样快，但所需的部分 SVD 数量明显较少。实验结果表明，IALM比APG至少快五倍，而且精度也更高。特别是，由 IALM 计算的 E 中非零项的数量比由 APG 计算的准确（实际上，通常是精确）得多，后者通常会在 E 中留下许多小的非零项。

在本文的其余部分中，为了完整起见，我们将首先在第 2 节中概述之前的工作。然后我们在第 3 节中介绍基于 ALM 的新算法并分析其收敛特性（同时将所有技术证明留给附录 A）。我们还将快速说明如何轻松地采用相同的 ALM 方法来解决（相关但稍微简单的）矩阵完成 (MC) 问题。然后，我们将在第 4 节中讨论算法的一些实现细节。接下来在第 5 节中，我们通过对随机生成的矩阵进行广泛的模拟，比较矩阵恢复和矩阵补全的新算法和其他现有算法。最后我们在第 6 节给出一些结论。

2 先前的矩阵恢复算法

在本节中，为了完整性和比较的目的，我们简要介绍和总结其他现有的解决矩阵恢复问题（2）的算法。

2.1 迭代阈值方法The Iterative Thresholding Approach

[22] 中提出的 IT 方法解决了 (2) 的松弛凸问题：

其中 τ 是一个较大的正标量，因此目标函数仅受到轻微扰动。通过引入拉格朗日乘数 Y 去除等式约束，得到式（3）的拉格朗日函数：

然后 IT 方法迭代更新 A、E 和 Y。它通过最小化 L(A, E, Y ) 更新 A 和 E，对于 A 和 E， Y 固定。然后违反约束 A + E = D 的量用于更新 Y 。

为了方便起见，我们引入以下软阈值（收缩）运算符：

其中 x ∈ R 且 ε > 0。通过逐元素应用，该运算符可以将其扩展到向量和矩阵。然后 IT 方法的工作原理如算法 1 中所述，其中阈值直接来自众所周知的分析 [7, 24]：

其中 USV T 是 W 的 SVD。 IT算法虽然极其简单且可证明正确，但需要大量的迭代才能收敛，并且难以选择步长δk来加速，因此其适用性受到限制。

IT算法原理：

2.2 加速近端梯度法The Accelerated Proximal Gradient Approach

加速近端梯度法的一般理论可以在[21,4,17]中找到。为了解决以下无约束凸问题：

其中 H 是一个实希尔伯特空间，具有内积 〈·,·〉 和相应的范数 ‖ · ‖，g 和 f 都是凸的，并且 f 进一步是 Lipschitz 连续的： ‖∇f (X1) − ∇f (X2) ‖ ≤ Lf ‖X1 − X2‖，可以将 f (X) 局部近似为二次函数并求解

它被假设为很容易更新解 X 。该迭代的收敛行为很大程度上依赖于形成近似值 Q(X, Yk) 的点 Yk。自然选择 Yk = Xk（例如，由[11]提出）可以解释为梯度算法，并且收敛速度不低于 O(k−1) [4]。然而，对于平滑 g, Nesterov 表明，对于满足 t2k+1 − tk+1 ≤ t2k 的序列 {tk}，设置 Yk = Xk + tk−1−1 *(Xk − Xk−1) / tk 可以将收敛速度提高到 O (k−2) [17]。最近，Beck 和 Teboulle 将这个方案推广到非光滑 g，再次证明了在 F (Xk) − F (X*) ≤ Ck−2的意义下，能有O(k−2) 的收敛速度，[4]。

上述加速近端梯度APG方法可以在满足

直接应用于 RPCA 问题的宽松版本，其中 μ 是一个小的正标量。连续技术[19]（可以改变μ，从一个大的初始值μ0开始，并在每次迭代中以几何方式减小它，直到达到下限μ_）可以大大加快收敛速度​​。 RPCA 的 APG 方法在算法 2 中描述（有关详细信息，请参阅 [15, 23]）。

2.3 The Dual Approach 双重方法

我们早期的工作  中提出的对偶方法[15]通过其对偶解决了问题 (2)。也就是说，首先为了最优拉格朗日乘子 Y，求解对偶问题

其中

‖ · ‖∞ 是最大的矩阵元素的绝对值。求解式(9)可以采用约束在曲面{Y |J(Y ) = 1}上的最速上升算法，其中约束最速上升方向是通过将D投影到凸的域{Y |J (Y)≤1}的切锥上得到的。事实证明，在寻找约束最速上升方向的过程中可以得到原问题(2)的最优解。最终算法的细节可以参考[15]。

对偶方法的优点是仅需要与最大奇异值 1 相关的主奇异空间。理论上，计算这个特殊的主奇异空间应该比计算与未知的前导奇异值相关的主奇异空间更容易。因此，如果可以获得一种计算与已知最大奇异值相关的主奇异空间的有效方法，那么对偶方法是有前途的。

3 增广拉格朗日乘子方法 The Methods of Augmented Lagrange Multipliers

在[5]中，引入了增广拉格朗日乘子的通用方法来解决此类约束优化问题：

其中 f : Rn → R 和 h : Rn → Rm。可以定义增广拉格朗日函数：

其中μ是正标量，然后可以通过增强拉格朗日乘子的方法来解决优化问题，概述为算法3（更多详细信息请参见[6]）。

在一些相当一般的条件下，当{μk}是递增序列并且f和h都是连续可微函数时，文献[5]已经证明，当{μk}有界时，算法3产生的拉格朗日乘子Yk Q-linearly 收敛到最优解；当{μk}无界时，算法3产生的拉格朗日乘子Yk Super-Q-linearly 收敛到最优解。 ALM 的这种卓越的收敛特性使其非常有吸引力。 ALM的另一个优点是更新Yk的最佳步长被证明是所选择的惩罚参数μk，这使得参数调整比迭代阈值算法容易得多。 ALM 的第三个优点是，即使不需要 μk 接近无穷大，算法也能收敛到精确的最优解 [5]。相反，严格来说，前面提到的迭代阈值和 APG 方法只能找到问题的近似解。最后，ALM 算法的（收敛）分析和实现相对简单，我们将在矩阵恢复和矩阵补全问题上进行演示。

注：Super-Q-linearly是数学中用来描述一个数列收敛到极限的速度的一个术语。它意味着数列收敛得比Q-linearly快，而Q-linearly是指相邻两项之间的误差的比值有一个小于一的常数上界。Super-Q-linear收敛发生在误差的比值有一个前一项误差的幂次的上界的时候。

例如，如果x_k是一个收敛到L的数列，那么它是super-Q-linear收敛的，如果

lim k → inf | x k + 1 − L | / | x k − L | ^ q = μ

对于某些q > 1 和 μ > 0。

3.1 两种用于鲁棒PCA（矩阵恢复）的ALM算法

对于 RPCA 问题（2），我们可以通过满足条件：

来应用增强拉格朗日乘子法。那么拉格朗日函数为：

解决 RPCA 问题的 ALM 方法可以在算法 4 中描述，我们将其称为精确 ALM (EALM) 方法，原因很快就会清楚。

算法中的初始值 Y0* = sgn(D)/J(sgn(D)) 受到对偶问题 (9) 的启发，因为它很可能使目标函数值<D,  Y0* >相当大。

尽管RPCA问题（2）的目标函数是非光滑的，因此[5]中的结果不能直接适用于此，但我们仍然可以证明算法4具有同样出色的收敛性。更准确地说，我们做出了以下声明。

定理1.

对于算法4，在

的意义下，其中f*是RPCA问题的最优值。(Ak*,Ek*)的任意累加点(A*,E*)都是RPCA问题的最优解，且收敛速度至少为O(μk−1)。

证明见附录A.3。

由定理1可知，如果μk呈几何增长，则EALM方法将Q线性收敛；如果 μk 增长得更快，那么 EALM 方法也会收敛得更快。然而，数值测试表明，对于较大的 μk，用于解决解决子问题 (A* k+1, E* k+1) = arg min A,E L(A, E, Y k*  , μk)的迭代阈值IT方法会收敛得更慢。由于 SVD 占据了大部分计算负载，因此 {μk} 的选择应该合适，以使 SVD 的总数最小。

幸运的是，事实证明，我们不必求出子问题

的精确解。相反，在解决这个子问题时更新一次 Ak 和 Ek 就足以让 Ak 和 Ek 收敛到 RPCA 问题的最优解。这导致了不精确ALM (IALM) 方法，如算法 5 中所述。

算法5的有效性和最优性由以下定理保证。

定理2

对于算法 5，如果 {μk} 是非递减的且

，则 (Ak, Ek) 收敛到 RPCA 问题的最优解 (A*, E*)。

证明见附录A.5。

我们可以进一步证明，条件   也是保证收敛所必需的，如下定理所述。

定理3

如果  

则算法5 产生的序列{(Ak, Ek)} 可能不会收敛到RPCA 问题的最优解。

证明见附录A.6。

请注意，与精确 ALM 方法的定理 1 不同，定理 2 只保证收敛，但没有指定非精确 ALM 方法的收敛速度。条件   意味着 μk 不能增长得太快。 μk 的选择将在第 4 节中详细讨论。

3.2 矩阵补全的ALM算法 An ALM Algorithm for Matrix Completion

矩阵补全（MC）问题可以被视为矩阵恢复问题的一种特殊情况，其中在给定有限数量的已知条目的情况下，必须恢复矩阵中丢失的条目。这样的问题无处不在，例如在机器学习[1,2,3]、控制[16]和计算机视觉[20]中。在许多应用中，假设【要恢复的矩阵是低秩的】是合理的。在最近的一篇论文 [9] 中，Cand`es 和 Recht 证明了大多数秩为 r 的矩阵 A 可以通过解决以下优化问题来完美恢复：

假设对于某些正常数 C，样本数 p 满足 p ≥ Crn6/5 ln n，其中 Ω 是样本索引集。此后，其他几个人的工作改进了这个界限。解决 MC 问题 (14) 的最先进算法包括 APG 方法 [19] 和奇异值阈值 (SVT) 方法 [7]。由于 RPCA 问题与 MC 问题密切相关，因此很自然地相信 ALM 方法在 MC 问题上也能同样有效。

我们可以将 MC 问题表述如下

其中 πΩ : Rm×n → Rm×n 是一个线性运算符，它保持 Ω 中的条目不变，并将 Ω 之外的条目（即 ̄ Ω 中的条目）设置为零。由于E将补充D的未知条目，因此D的未知条目被简单地设置为零。

那么(15)的部分增广拉格朗日函数（[5]的2.4节）是

然后类似地，我们可以针对 MC 问题使用精确和不精确的 ALM 方法，其中为了更新 E，在最小化 L(A, E, Y, μ) 时应强制执行约束 πΩ(E) = 0。算法 6 描述了不精确的 ALM 方法。

请注意，由于 Ek 的选择，π ̄ Ω(Yk) = 0 在整个迭代过程中保持不变，即未知条目处的 Yk 值始终为零。定理 1 和 2 对于矩阵补全问题也成立。由于证明类似，这里省略。定理 3 对于矩阵补全问题也成立，因为很容易验证 Yk = πΩ( ˆ Yk)。由于 { ˆ Yk} 是有界的（参见引理 1），因此 {Yk} 也是有界的。所以定理3的证明对于矩阵补全问题仍然有效。

4 实施细节 Implementation Details

预测主奇异空间的维数。显然，计算 RPCA 和 MC 问题的完整 SVD 是不必要的：我们只需要那些大于特定阈值的奇异值及其相应的奇异向量。因此一个软件包PROPACK[14]在社区得到了广泛推荐。要使用 PROPACK，必须预测奇异值大于给定阈值的主奇异空间的维数。对于算法 5，预测相对容易，因为观察到 Ak 的秩是单调递增的，并且在真实秩上变得稳定。所以预测规则是：

其中 d = min(m, n)，svk 为预测维度，svpk 为 svk 奇异值中大于 μk−1 的奇异值个数，sv0 = 10。算法 4 也采用上述预测策略对于求解 (A* k+1, E* k+1) 的内循环。对于外层循环，预测规则就是 svk+1 = min(svpk + round(0.1d), d)。对于算法6，预测要困难得多，因为Ak的秩经常振荡。通常情况下，对于较小的 k，Ak 的秩接近 d，然后逐渐下降到真实的秩，使得部分 SVD 效率低下数值测试表明，当我们想要计算超过 0.2d 个主奇异向量/值时，使用 PROPACK 通常比计算完整的 SVD 慢。。为了解决这个问题，我们将 Y 和 A 都初始化为零矩阵，并采用下面的（与[19]中策略类似的）截断策略：

sv0 = 5且

其中maxgapk和maxidk分别是连续奇异值（将计算出的svk奇异值按降序排列）与相应索引之间的最大比率。我们利用间隙信息是因为我们观察到奇异值很快被分为两组，它们之间的间隙很大，使得排名显示快速且可靠。通过上述预测方案，Ak的秩变得单调递增并且稳定在真实的秩上。

Order of Updating A and E.更新 A 和 E 的顺序。虽然理论上先更新 A 和 E 中的任何一个不会影响收敛速度，但数值测试表明，这确实会导致达到相同精度的迭代次数略有不同。考虑到 SVD 对于大维矩阵的巨大复杂性，也应该考虑这种微小的差异。通过广泛的数值测试，我们建议在算法 4 和 5 中首先更新 E。同样重要的是，首先更新 E 也使得 Ak 的秩更有可能单调递增，这对于部分 SVD 的有效至关重要，因为已在上一段中详细阐述。

Memory Saving for Algorithm 6.算法 6 的内存节省。在算法 6 的实际实现中，稀疏矩阵用于存储 D 和 Yk，如[19]中所做的那样，A 表示为 A = LRT ，其中 L 和 R 都是大小为 m × svpk的矩阵。Ek 没有通过观察

明确地存储。这样，计算Yk和D − Ek + μk−1  Yk只有πΩ (Ak)是必要的。由于样本比例较小，因此可以节省大量内存。

Stopping Criteria停止标准。对于RPCA问题，KKT条件为：

当且仅当∂‖A*‖* ∩ ∂(‖λE*‖1) ≠ ∅ 时，后两个条件成立。因此我们可以将以下条件作为算法4和5的停止标准：

其中 dist(X, Y ) = min{‖x − y‖F |x ∈X, y ∈Y }。对于算法 4，第二个条件始终由内循环保证。所以我们只需要检查第一个条件。不幸的是，对于算法 5，计算 dist(∂‖Ak‖*, ∂(‖λEk‖1)) 的成本很高，因为到 ∂‖Ak‖* 上的投影成本很高。因此，我们可以通过 ‖ ˆ Yk − Yk‖F = μk−1‖Ek − Ek−1‖F 来估计 dist(∂‖Ak‖*, ∂(‖λEk‖1))，因为 ˆ Yk ∈ ∂‖Ak‖* 并且Yk ∈ ∂(‖λEk‖1)。

类似地，对于MC问题我们可以将以下条件作为算法6的停止标准：

其中 S = {Y |Yij = 0 if (i, j) 不属于 Ω}。同样，由于计算 dist(∂‖Ak‖*, S) 的成本很高，我们可以通过 ‖ ˆ Yk−Yk‖F = μk−1‖Ek−Ek−1‖F 来估计 dist(∂‖Ak‖*, S) 。由于μk−1‖Ek−Ek−1‖F  实际上显着高估了 dist(∂‖Ak‖*, S)，在实际计算中我们可以对算法 6 使用以下停止标准：

注意，通过 (20) ‖Ek−Ek−1‖F可以方便地计算为

Updating μk. 更新μk。对于RPCA问题，更新规则为：

其中ρ > 1。不难看出，这个更新规则与定理2是一致的。对于MC问题，更新规则为：

Choosing Parameters.选择参数。对于算法 4，我们设置 μ0 = 0.5/‖ sgn(D)‖2 和 ρ = 6。内循环的停止标准为 ‖Ak j+1 − Ak j‖F /‖D‖F < 10−6 且‖Ek j +1 −Ek j ‖F /‖D‖F < 10−6。外层迭代的停止标准是 ‖D−Ak*  −Ek*‖F /‖D‖F < 10−7。对于算法 5，我们设置 μ0 = 1.25/‖D‖2 且 ρ = 1.6。停止准则中的参数为 ε1 = 10−7 和 ε2 = 10−5。对于算法6，我们设置μ0 = 1/‖D‖2和ρ = 1.2172 + 1.8588ρs，其中ρs = |Ω|/(mn)是采样密度，通过回归得到ρ和ρs之间的关系。停止准则中的参数为 ε1 = 10−7 和 ε2 = 10−6。

5 Simulations 模拟

在本节中，使用数值模拟，对于 RPCA 问题，我们将所提出的 ALM 算法与[15]中提出的 APG 算法进行比较；对于MC问题，我们将不精确ALM算法与SVT算法[7]和APG算法[19]进行比较。所有模拟均在同一工作站上进行和计时，该工作站配备 Intel Xeon E5540 2.53GHz CPU，具有 4 核和 24GB 内存3，运行 Windows 7 和 Matlab（版本 7.7）。

I. Comparison on the Robust PCA Problem. 鲁棒PCA问题的比较。对于 RPCA 问题，我们使用随机生成的方阵进行模拟。我们用有序对 (A*, E*) ∈ Rm×m × Rm×m 表示真实解。我们生成 R 阶矩阵 A* 作为乘积 LRT ，其中 L 和 R 是独立的 m × r 矩阵，其元素是独立同分布的。具有零均值和单位方差的高斯随机变量。 我们生成 E* 作为稀疏矩阵，其支持度是随机均匀选择的，其非零项是 i.i.d 于均匀分布在区间 [−500, 500] 内。矩阵 D. = A* + E* 是算法的输入，( ˆ A, ˆ E) 表示输出。对于给定的问题，我们选择固定的权重参数 λ = m−1/2。

我们使用[15]作者提供的算法 2的最新代码，并应用预测规则 (17)，其中 sv0 = 5，以便可以利用部分 SVD。使用部分 SVD，APG 比 2.3 节中的双重方法更快。所以我们不需要涉及双重方法来进行比较。

表 1 和表 2 给出了这三种算法的简要比较。我们可以看到，APG 和 IALM 算法都停止在相对恒定的迭代次数上，并且 IALM 至少比 APG 快五倍。而且，EALM和IALM的精度高于APG。特别是，APG 经常过高估计||E*||0，即 E* 中非零的数量。而 EALM 和 IALM 估计的 ‖E*‖0 总是非常接近真实情况。

II. Comparison on the Matrix Completion Problem.矩阵补全问题的比较。对于 MC 问题，首先生成真正的低秩矩阵 A*，就像 RPCA 问题一样。然后我们从 A* 中均匀采样 p 个元素，形成 D 中的已知样本。一个有用的参考量是 dr = r(2m − r)，它是秩为 r 的 m × m 矩阵中的自由度数 [ 19]。

SVT和APGL（APG with line search7）代码分别由[7]和[19]的作者提供。表 3 对这三种算法进行了简要比较。可以看出，IALM 始终比 SVT 更快。当采样密度 p/m2 相对较高时，例如 p/m2 > 10%，它也比 APGL 更有优势。这种现象实际上与 RPCA 问题的结果一致，其中 D 的大多数样本被假设是准确的，尽管准确样本的位置是先验未知的。

6 Conclusions结论

在本文中，我们提出了两种基于增强拉格朗日乘子的算法，即 EALM 和 IALM，用于解决鲁棒 PCA 问题 (2)。这两种算法都比之前最先进的 APG 算法 [15] 更快。特别是，在所有模拟中，IALM 始终比 APG 快五倍以上。

我们还将增广拉格朗日乘子的方法应用于矩阵补全问题。相应的 IALM 算法比著名的 SVT 算法 [7] 快得多。当可用条目的百分比不太低（例如 > 10%）时，它也比最先进的 APGL 算法 [19] 更快。

与基于加速近端梯度的方法相比，基于增强拉格朗日乘子的算法更易于分析且更易于实现。此外，它们的准确性也高得多，因为即使惩罚参数不接近无穷大，迭代也被证明可以收敛到问题的精确解[5]。相比之下，APG 方法通常通过解决松弛问题来找到近似解。最后，对于 RPCA 和 MC 问题，ALM 算法比 APG 需要更少的存储/内存8。对于大规模应用程序，例如 Web 数据分析，这可能被证明是 ALM 类型算法的一大优势。

为了帮助读者比较和使用所有算法，我们在网站上发布了所有算法的 Matlab 代码：“http://perception.csl.illinois.edu/matrix-rank/home.html” (Lin 等, 2013, p. 14)

总结

1.PCA

假设给定的数据被排列为大矩阵 D ∈ Rm×n 的列。估计低维子空间的数学模型是找到一个低秩矩阵A，使得A和D之间的差异最小化，从而导致以下约束优化：

其中 r ≪ min(m, n) 是子空间的目标维度， ‖ · ‖F 是 Frobenius 范数，其对应于假设数据被 i.i.d的 高斯噪声损坏。这个问题可以通过先计算 D 的奇异值分解（SVD），然后将 D 的列向量投影到 D 的 r 个主左奇异向量所组成的子空间上来得到方便地解决[13]。

当损失是由附加的i.i.d 高斯噪声引起时，PCA 会给出最佳估计。只要噪声的幅度很小，它在实践中就可以很好地发挥作用。

2.RPCA

PCA在严重损失的情况下就会崩溃，即使这种损失只影响极少数的观测结果。事实上，即使只有 A 的一项被任意破坏，经典 PCA 得到的估计 ^ A 也可能与真实 A 相差任意远。所以我们研究秩矩阵 A 是否仍然可以有效且从损坏的数据矩阵 D = A + E 中准确恢复，其中附加误差 E 的某些条目可能是任意大。

最近，赖特等人。 [22]表明，在相当广泛的条件下，答案是肯定的：只要误差矩阵 E 足够稀疏（相对于 A 的秩），通过求解以下凸优化问题就可以从 D = A + E 精确地恢复低秩矩阵 A：

其中 ‖ · ‖* 表示矩阵的核范数（即其奇异值之和）， ‖ · ‖1 表示矩阵条目的绝对值之和， λ 是正权重参数。由于即使存在较大错误或异常值，也能够准确恢复数据中的底层低秩结构，因此这种优化在[22]中被称为鲁棒 PCA (RPCA)

3.增广拉格朗日乘子方法 The Methods of Augmented Lagrange Multipliers

在[5]中，引入了增广拉格朗日乘子的通用方法来解决此类约束优化问题：

其中 f : Rn → R 和 h : Rn → Rm。可以定义增广拉格朗日函数：

其中μ是正标量，然后可以通过增强拉格朗日乘子的方法来解决优化问题，概述为算法3（更多详细信息请参见[6]）。

3.1 两种用于鲁棒PCA（矩阵恢复）的ALM算法

对于 RPCA 问题（2），我们可以通过满足条件：

来应用增强拉格朗日乘子法。那么拉格朗日函数为：

解决 RPCA 问题的 ALM 方法可以在算法 4 中描述，我们将其称为精确 ALM (EALM) 方法，原因很快就会清楚。

算法中的初始值 Y0* = sgn(D)/J(sgn(D)) 受到对偶问题 (9) 的启发，因为它很可能使目标函数值<D,  Y0* >相当大。

尽管RPCA问题（2）的目标函数是非光滑的，因此[5]中的结果不能直接适用于此，但我们仍然可以证明算法4具有同样出色的收敛性。更准确地说，我们做出了以下声明。

定理1.

对于算法4，在

的意义下，其中f*是RPCA问题的最优值。(Ak*,Ek*)的任意累加点(A*,E*)都是RPCA问题的最优解，且收敛速度至少为O(μk−1)。

证明见附录A.3。

由定理1可知，如果μk呈几何增长，则EALM方法将Q线性收敛；如果 μk 增长得更快，那么 EALM 方法也会收敛得更快。然而，数值测试表明，对于较大的 μk，用于解决解决子问题 (A* k+1, E* k+1) = arg min A,E L(A, E, Y k*  , μk)的迭代阈值IT方法会收敛得更慢。由于 SVD 占据了大部分计算负载，因此 {μk} 的选择应该合适，以使 SVD 的总数最小。

幸运的是，事实证明，我们不必求出子问题

的精确解。相反，在解决这个子问题时更新一次 Ak 和 Ek 就足以让 Ak 和 Ek 收敛到 RPCA 问题的最优解。这导致了不精确ALM (IALM) 方法，如算法 5 中所述。

算法5的有效性和最优性由以下定理保证。

定理2

对于算法 5，如果 {μk} 是非递减的且

，则 (Ak, Ek) 收敛到 RPCA 问题的最优解 (A*, E*)。

证明见附录A.5。

我们可以进一步证明，条件   也是保证收敛所必需的。

实施细节：

对于RPCA和MC问题，我们只需要那些大于特定阈值的奇异值及其相应的奇异向量。因此一个软件包PROPACK[14]在社区得到了广泛推荐。要使用 PROPACK，必须预测奇异值大于给定阈值的主奇异空间的维数。对于算法 5，预测相对容易，因为观察到 Ak 的秩是单调递增的，并且在真实秩上变得稳定。所以预测规则是：

其中 d = min(m, n)，svk 为预测维度，svpk 为 svk 奇异值中大于 μk−1 的奇异值个数，sv0 = 10。算法 4 也采用上述预测策略对于求解 (A* k+1, E* k+1) 的内循环。对于外层循环，预测规则就是 svk+1 = min(svpk + round(0.1d), d)。

a.建议在算法 4 和 5 中首先更新 E。

b.Stopping Criteria停止标准。对于RPCA问题，KKT条件为：

当且仅当∂‖A*‖* ∩ ∂(‖λE*‖1) ≠ ∅ 时，后两个条件成立。因此我们可以将以下条件作为算法4和5的停止标准：

其中 dist(X, Y ) = min{‖x − y‖F |x ∈X, y ∈Y }。对于算法 4，第二个条件始终由内循环保证。所以我们只需要检查第一个条件。不幸的是，对于算法 5，计算 dist(∂‖Ak‖*, ∂(‖λEk‖1)) 的成本很高，因为到 ∂‖Ak‖* 上的投影成本很高。因此，我们可以通过 ‖ ˆ Yk − Yk‖F = μk−1‖Ek − Ek−1‖F 来估计 dist(∂‖Ak‖*, ∂(‖λEk‖1))，因为 ˆ Yk ∈ ∂‖Ak‖* 并且Yk ∈ ∂(‖λEk‖1)。

c.Updating μk. 更新μk。对于RPCA问题，更新规则为：

其中ρ > 1。不难看出，这个更新规则与定理2是一致的。

d.Choosing Parameters.选择参数。

对于算法 4，我们设置 μ0 = 0.5/‖ sgn(D)‖2 和 ρ = 6。内循环的停止标准为 ‖Ak j+1 − Ak j‖F /‖D‖F < 10−6 且‖Ek j +1 −Ek j ‖F /‖D‖F < 10−6。外层迭代的停止标准是 ‖D−Ak*  −Ek*‖F /‖D‖F < 10−7。

对于算法 5，我们设置 μ0 = 1.25/‖D‖2 且 ρ = 1.6。停止准则中的参数为 ε1 = 10−7 和 ε2 = 10−5。

3.2 矩阵补全的ALM算法

大多数秩为 r 的矩阵 A 可以通过解决以下优化问题来完美恢复：

假设对于某些正常数 C，样本数 p 满足 p ≥ Crn6/5 ln n，其中 Ω 是样本索引集。此后，其他几个人的工作改进了这个界限。解决 MC 问题 (14) 的最先进算法包括 APG 方法 [19] 和奇异值阈值 (SVT) 方法 [7]。由于 RPCA 问题与 MC 问题密切相关，因此很自然地相信 ALM 方法在 MC 问题上也能同样有效。

我们可以将 MC 问题表述如下

其中 πΩ : Rm×n → Rm×n 是一个线性运算符，它保持 Ω 中的条目不变，并将 Ω 之外的条目（即 ̄ Ω 中的条目）设置为零。由于E将补充D的未知条目，因此D的未知条目被简单地设置为零。

那么(15)的部分增广拉格朗日函数（[5]的2.4节）是

然后类似地，我们可以针对 MC 问题使用精确和不精确的 ALM 方法，其中为了更新 E，在最小化 L(A, E, Y, μ) 时应强制执行约束 πΩ(E) = 0。算法 6 描述了不精确的 ALM 方法。

实施细节：

对于算法6，预测要困难得多，因为Ak的秩经常振荡。通常情况下，对于较小的 k，Ak 的秩接近 d，然后逐渐下降到真实的秩，使得部分 SVD 效率低下数值测试表明，当我们想要计算超过 0.2d 个主奇异向量/值时，使用 PROPACK 通常比计算完整的 SVD 慢。。为了解决这个问题，我们将 Y 和 A 都初始化为零矩阵，并采用下面的（与[19]中策略类似的）截断策略：

sv0 = 5且

其中maxgapk和maxidk分别是连续奇异值（将计算出的svk奇异值按降序排列）与相应索引之间的最大比率。我们利用间隙信息是因为我们观察到奇异值很快被分为两组，它们之间的间隙很大，使得排名显示快速且可靠。通过上述预测方案，Ak的秩变得单调递增并且稳定在真实的秩上。

a.在实际计算中我们可以对算法 6 使用以下停止标准：

注意，通过 (20) ‖Ek−Ek−1‖F可以方便地计算为

b.对于MC问题，更新规则为：

c.对于算法6，我们设置μ0 = 1/‖D‖2和ρ = 1.2172 + 1.8588ρs，其中ρs = |Ω|/(mn)是采样密度，通过回归得到ρ和ρs之间的关系。停止准则中的参数为 ε1 = 10−7 和 ε2 = 10−6。