笔记：Online robust principal component analysis via truncated nuclear norm regularization

Hong, B., Wei, L., Hu, Y., Cai, D., & He, X. (2016). Online robust principal component analysis via truncated nuclear norm regularization. Neurocomputing, 175, 216-222.

本文是这篇 Neurocomputing 期刊论文的笔记，主要是对文中的理论方法进行展开详解。本人学术水平有限，文中如有错误之处，敬请指正。

摘要： Robust principal component analyssi (RPCA) 已经被广泛用于处理高维的噪声数据，在许多应用功中。传统的 RPCA 方法考虑所有的样本恢复低维的子空间用批量的方式，导致了昂贵的存储代价并且不能有效地更新流数据的低维子空间。所以有必要设计一种在线的 RPCA 方法。此文中，提出了一种新颖的 online RPCA 算法，采取最近提出的 truncated nuclear norm 作为低秩约束的更好的近似。这里将目标函数按样本分解代价的和，并设计了 online 有效的交替优化方法。

1 简介

在许多机器学习和数据挖掘问题中，经常遇到高维的样本，包含一些噪声（损坏或奇异点）。为了恢复内部的低维子空间，从全部的样本集中， RPCA 被大量地研究，应用于视频监控 1，图像配准 2，文本语料建模 3 和 音频处理 4 。

原理上，典型的 RPCA 方法假设样本可以被分为低秩的部分和稀疏的部分。正式地，给定一个样本 Z∈Rm×n， RPCA 尝试将 Z 分解为一个低秩的矩阵 X 和一个稀疏的矩阵 E 的和

minX,Es.t. rank(X)+λ||E||0 Z=X+E,(1)

其中 λ 是一个约束的参数。

已经被证明低维的子空间可以在合适的条件下，被精确地、有效地恢复。然而，该问题是高度非凸的，不易处理的，因为秩函数和 ℓ0 范数。大多数研究在寻找合适的秩函数和 ℓ0 范数的替代，将原问题转化为一个凸的优化问题。其中，Lin et al. 应用増广 Lagrange 乘子来得到凸问题 5。Shang et al. 6 和 Tao et al. 7 考虑更一般的情况，观测的数据是缺失的并被严重破坏的，提出了一种统一的框架，结合了 RPCA 和矩阵补全方法。

以上所有的方法都是处理批量数据的。也就是每一次迭代中，所有的样本都是需要使用的，这造成了两种限制。首先，存储代价是昂贵的，需要内存中有所有的样本在优化过程中，尤其是对于大规模数据是不可接受的。另一方面，如果数据是以流的形式获得，这些方法不能有效地处理低维子空间当一个新样本到来时。

为了解决这个问题，online RPCA 方法出现了。内存消耗与样本的规模是无关的，并发现到的低维子空间可以快速更新。另一个重要的 online RPCA 的优势是它可以跟踪动态的低维子空间，当其会随着时间变化时。所以 online RPCA 可以被用于移动摄像头的视频跟踪 8 。Goes et al. 扩展了批量版本的 PRCA 到随机，并提供了一个子线性收敛保证 9，明显地减少了存储的要求和时间复杂度。He et al. 提出了在线自适应子空间跟踪算法基于 Grassmannian 流形 10，其结合了増广 Lagrangian 和经典的随机梯度框架。Mairal 提出了更一般的 在线字典学习机制为了稀疏编码基于随机近似 11 。受到次启发之后，Feng et al. 12 和 Shen et al. 13 尝试用在线方式解决 RPCA 问题。他们分别采用了 核范数和最大范数，作为秩函数的代替，两者都可以被表示为顺序数据的矩阵分解的形式。尽管核范数和最大范数是矩阵的秩函数的凸包络，但是也导致了不能忽视的近似误差，在真实的应用中 14 。所以，一些研究者尝试设计非凸的代替，来实现更精确的近似 15 。

此文的目标是解决解决 RPCA 问题，通过一个在线非凸的优化框架。特别地，此文用最小化一个最近提出的 truncated nuclear norm 16 来代替目标函数，最小化矩阵的秩。此范数也可以被表示为矩阵分解的形式，其提供了思路来估计每一个样本对于 truncated 范数的增量的贡献。基于此，此文提出了一种用新样本更新低维空间的 online 机制。接着设计了一种有效的、迭代优化方法的实现。通过 truncated 范数，此算法的优化可以更接近矩阵的秩，子空间恢复也可以更精确。此文的主要贡献是两方面：

此文提出了一个 online 机制来解决 RPCA 问题，通过采用矩阵的非凸的近似，相比于凸的代替更为精确。

此文设计了一个高效的优化算法解决提出的目标函数。

2 预定义

大写加粗字母表示矩阵，小写加粗字母表示向量。||X||1，||X||∗ 和 ||X||F 分别表示 ℓ1，核范数和 Frobenius 范数。tr(⋅) 表示方阵的迹函数。||v||1 和 ||v||p 表示向量的 ℓp 范数。⟨⋅,⋅⟩ 表示内积。I 表示单位矩阵。

给定一个矩阵 X∈Rm×n 和一个非负的整数 s<min(m,n)，truncated 范数 ||X||s 定义为最小的 min(m,n)−s 个奇异值之和，也就是 ||X||s=∑min(m,n)i=s+1σi(X) ，其中 σ1(X)≥⋯≥σmin(m,n)(X) 。换句话说，||X||s 不关心最大的 s 个奇异值，两者的关系阐述为如下

||X||s=||X||∗−maxUUT=I, VVT=Itr(UXVT),(2)

其中 U∈Rs×m， V∈Rs×n 。公式中并不能明显看出范数和每一个样本的关系，很难估计每一个样本对范数的单独的贡献。幸运的是，核范数可以被分解为

||X||∗=minX=LRT 12(||L||2F+||R||2F),(3)

其中 L∈Rm×d，R∈Rn×d 对任意的 d≥rank(X)。

Lemma 2.1 truncated 范数可以分解为

||X||ss.t.=minL,R,U,V 12||L||2F+12||R||2F−tr(ULRTVT), X=LRT, UUT=I, VVT=I,(4)

其中 U∈Rs×m，V∈Rs×n，L∈Rm×d，R∈Rn×d，d≥rank(X) 。

Proof 对于任意的 U,V,L,R 满足 X=LRT，UUT=I，VVT=I，

||X||s=||X||∗−maxU,V tr(UXVT)≤12||L||2F+12||R||2F−maxU,V tr(UXVT)≤12||L||2F+12||R||2F−tr(UXVT).(5)

另一方面，假设矩阵的奇异值分解 X=PΣQT，其中 P=(p1,⋯,pm)∈Rm×m，Q=(q1,⋯,qn)∈Rn×n 和 Σ∈Rm×n。令 U^=(p1,⋯,ps)T 和 V^=(q1,⋯,qs)T，然后

tr(U^XV^T)=∑i=1s(X).

令 L^=PΣ1/2 和 R^=QΣ1/2，可以直接得到 X=L^R^T 和 ||X||∗=12||L^||2F+12||R^||2F 。

||X||s=||X||∗−∑i=1s(X)=12||L^||2F+12||R^||2F−tr(U^XV^T).(6)

该分解将 ||X||s 基于维度削减。如此，L 可以被看做字典，而 R 的每一列可以看做系数。

3 此文提出的算法

此算法的目标是：给定一个 m 维的数据集 Z=(z1,⋯,zn)∈Rm×n，要将其分解为低秩矩阵 X 和稀疏矩阵 E，对每一个样本 zi=xi+ei 。不同于传统的方法，采用核范数作为秩函数的近似，truncated 核范数在本文中采用。所以此文的目标函数可以写为

minX,E 12||Z−X−E||2F+λ1||X||s+λ2||E||1,(7)

其中 λ1,λ2 是约束系数。注意使用 ℓ1 范数而不是 ℓ0 范数来约束稀疏项 E，因为 ℓ1 范数计算更易处理，通常在实际方法中被采用，获得稀疏解。||X||s 是一个整体的形式。为了获得更多关于低维空间 X 的结构信息，将其分解 X=LRT, L∈Rm×d, R∈Rn×d, d≥rank(X) 。在 online RPCA 方法中，L 视为字典，X 的每一列都当成 L 的元素关于 R 的每一行的系数的线性组合。结合了矩阵的分解，原目标函数可以转化为如下的形式

minL,R,U,V,Es.t. 12||Z−LRT−E|<
24000
span class="msubsup" id="MathJax-Span-1963">|2F+λ1(12||L||2F+12||R||2F−tr(ULRTVT))+λ2||E||1 UUT=I, VVT=I.(8)

该形式提供了一种解释：每一个样本 zi 近似 Lri+ei，其中 rTi 是 R 的第 i 行。根据 ||⋅||F 和 ||⋅||1 的加法性质，以上的问题可以分解为每一个样本的形式

minL,R,U,V,Es.t. 12||zi−Lri−ei||22+λ1(12||L||2F+12∑i=1n||ri||22−∑i=1nwTiri)+λ2∑i=1n||ei||1 UUT=I, VVT=I,(9)

其中 wi 是矩阵 W=VTUL∈Rn×d 的第 i 行。这里使用了如下的迹函数的交换性质：tr(ABC)=tr(CAB) 。为了简化形式，定义 f(L,zi,ri,ei)≜12||zi−Lri−ei||22+λ1(12||ri||22−wTiri)+λ2||ei||1 来整合一个样本 zi 对目标函数的贡献。可以将以上的目标函数化简为

minL,R,U,V,Es.t. ∑i=1nf(L,zi,ri,ei)+λ12||L||2F UUT=I, VVT=I.(10)

从中，可以看出目标函数是样本逐渐累加起来的，给定字典 L，就等价于最小化平均代价

J(L,n)=≜1n∑i=1nf~(L,zi)+λ12n||L||2F,(11)

其中 f~ 是每一个样本的损失函数，在最优的字典表示下

f~(L,z)s.t. =minr,e,U,V f(L,z,r,e) UUT=I, VVT=I.(12)

至此，已经将原优化问题转化为平均代价的最小化问题。其中每一个样本是在已知字典 L 的情况下获得。

4 优化

此文采用在线的方式交替地更新变量 L,R,U,V,E 假设样本是以流的形式到来，并且当前的样本是 zt，优化步骤可以分为两个连续的部分。第一，首先优化向量 rt,et 在已知 Lt−1,Ut−1,Vt−1 的情况下，通过求解如下的优化问题

{rt,et}=argminr,e 12||zt−Lt−1r−e||22+λ1(12||r||22−wTtr)+λ2||e||1,(13)

其中 wt 是 矩阵 Wt−1=VTt−1Ut−1Lt−1 的第 t 行。第二步，优化变量 Lt,Vt,Ut，使用之前已知的 {ri}ti=1,{ei}ti=1，通过求解以下的优化问题（无关项已删除）

{Lt,Vt,Ut}s.t. =argminL,V,U 12∑i=1t||zi−Lri−ei||22+λ1(12||L||2F−tr(ULRTtVT)), UUT=I, VVT=I,(14)

其中 RTt=(r1,⋯,rt,0,⋯,0)∈Rd×n 。值得注意的是，对于每一个新的样本 zt，Lt,Vt,Ut 是完全更新的（其中所有的元素都改变），而最优的 rt 只是增加到 R 的第 t 行之中。类似地，et 增加到 E 的第 t 列。

更新 rt ：

f(r)=12||zt−Lt−1r−ekt||22+λ1(12||r||22−wTtrt).(15)

令 ∂f/∂r=0, 可以得到如下的闭式解

rk+1t=(LTt−1Lt−1+λ1I)−1(LTt−1(zt−ekt)+λ1wt).(16)

更新 et：

g(e)=12||e||22−(zt−Lt−1rk+1t)Te+λ2||e||1.(17)

求解 e 可以使用标准的内点法因为 g(e) 是凸的。然而此方法是很费时的。注意到 g(e) 是两个凸函数的和，涉及 ℓ1 范数约束，可以使用分离固定点算法。定义 shrinkage 操作

Sλ(x)=⎧⎩⎨x−λ,x+λ,0,if x>λ,if x<−λ,otherwise.

此函数是 element-wise 的。获得如下的闭式解

ek+1t=Sλ2(zt−Lt−1rk+1t).(18)

更新 Lt：

h(L)=12∑i=1t||zi−Lri−ei||22+λ1(12||L||2F−tr(Ut−1LRTtVTt−1)).(19)

使用块坐标下降法更新字典的每一列，令 A=λ1I+∑ti=1rirTi=(a1,⋯,ad)，B=∑ti=1(zi−ei)rTi=(b1,⋯,bd)，C=UTVRt=(c1,⋯,cd)，Lt=(lt,1,⋯,lt,d)，那么字典 Lt 的每一列都可以更新

lt,j=1Ajj(bj+λ1cj−Lt−1aj)+lt−1,j, j=1,⋯,d.(20)

更新 Ut：

Ut=s.t. argmaxU tr(ULtRTtVTt−1) UUT=I.(21)

这是一个正交约束问题，通常是很困难的因为其非凸性质，保证的代价太昂贵在迭代中。这里提出了一个简单、但是有效的算法求解该问题，基于以下法则：

Lemma 4.1 假设 X∈Rm×n (m<n) 满足 XXT=I 。则可以获得其中一个最优解

maxXs.t. tr(XM) XXT=I.(22)

是 X∗=(Q,0)PT，其中 P,Q 由 M 奇异值分解得到：M=PΣQT，P∈Rn×n，Q∈Rm×m，Σ∈Rn×m，PTP=I，QTQ=I 。

Proof 通过假设，tr(XM)=tr(XPΣQT)=tr(QTXPΣ) 。令 X~=QTXP，接着有 X~X~T=QTXPPTXTQ=I 。所以，有 tr(XM)=tr(X~Σ)=∑mi=1X~iiσi，其中 σi 是矩阵 M 的奇异值分解。由于 |X~ij|<1 和 σi≥0，∀i,j，tr(XM) 取得最大值在集合 {X~∣X~X~T=I, Xii~=1, if σi>0}。一种特殊情况就是 X~∗=(I,0) 。这样的话，可以获得最优解之一

X~∗=QQTX∗PPT=QX~∗PT=Q(I,0)PT=(Q,0)PT.(23)

由于 Ut∈Rs×m, UtUTt=I, s<m 。原问题的形式与该定理一致，可以直接求解 LtRtVTt−1 的奇异值分解。Ut 可以可以直接由该 Lemma 得到。

更新 Vt：

Vt=s.t. argmaxV tr(VRtLTtUTt) VVT=I.(24)

该问题与求解 Ut 形式一致，同样可以使用 Lemma 求解。总的算法流程总结于 Algorithm 1 中。

Algorithm 1 Online RPCA 通过 Truncated nuclear norm

Input: 数据 Z=(z1,⋯,zn)∈Rm×n, 约束系数 λ1,λ2，矩阵 L0∈Rm×d,U0∈Rs×m,V0∈Rs×n；

Initialize: 随机初始化 L0，随机单位化 U0,V0 。

for t=1,⋯,n do

   Step 1: 计算 rt,et；

   初始化 rt=0,et=0；

   令 wt 取自 VTt−1Ut−1Lt−1 的第 t 行；

   repeat

     计算 rt←(LTt−1Lt−1+λ1I)−1(LTt−1(zt−et−1)+λ1wt)；

     计算 et←Sλ2(zt−Lt−1rt)；

   until 收敛

   Step 2 更新 Lt,Ut,Vt；

   repeat

   令 RTt=(r1,⋯,rt,0,⋯,0)∈Rd×n 更新 Lt 的列；

     [PU,ΣU,QU]=svd(LtRTtVTt−1)，

     Ut←(QU,0)PTU，

     [PV,ΣV,QV]=svd(RtLTtUTt)，

     Vt←(QV,0)PTV 。

   until 收敛 。

end for

Output: Ln,Rn 。

4 实验

此文的实验过于简单。略

J. Wright, A. Ganesh, S. Rao, Y. Peng, Y. Ma, Robust principal component analysis: exact recovery of corrupted low-rank matrices via convex optimization, in: Advances in Neural Information Processing Systems, 2009, pp. 2080–2088. ↩
Peng, Y., Ganesh, A., Wright, J., Xu, W., & Ma, Y. (2012). RASL: Robust alignment by sparse and low-rank decomposition for linearly correlated images. IEEE Transactions on Pattern Analysis and Machine Intelligence, 34(11), 2233-2246. ↩
Min, K., Zhang, Z., Wright, J., & Ma, Y. (2010, October). Decomposing background topics from keywords by principal component pursuit. In Proceedings of the 19th ACM international conference on Information and knowledge management (pp. 269-278). ↩
Huang, P. S., Chen, S. D., Smaragdis, P., & Hasegawa-Johnson, M. (2012, March). Singing-voice separation from monaural recordings using robust principal component analysis. In Acoustics, Speech and Signal Processing (ICASSP), IEEE International Conference on (pp. 57-60). ↩
Lin, Z., Chen, M., & Ma, Y. (2010). The augmented lagrange multiplier method for exact recovery of corrupted low-rank matrices. arXiv preprint arXiv:1009.5055. ↩
Shang, F., Liu, Y., Cheng, J., & Cheng, H. (2014, November). Robust principal component analysis with missing data. In Proceedings of the 23rd ACM International Conference on Conference on Information and Knowledge Management (pp. 1149-1158). ↩
Tao, M., & Yuan, X. (2011). Recovering low-rank and sparse components of matrices from incomplete and noisy observations. SIAM Journal on Optimization, 21(1), 57-81. ↩
Feng, J., Xu, H., & Yan, S. (2013). Online robust pca via stochastic optimization. In Advances in Neural Information Processing Systems (pp. 404-412). ↩
Goes, J., Zhang, T., Arora, R., & Lerman, G. (2014). Robust Stochastic Principal Component Analysis. In AISTATS (pp. 266-274). ↩
He, J., Balzano, L., & Lui, J. (2011). Online robust subspace tracking from partial information. arXiv preprint arXiv:1109.3827. ↩
Min, K., Zhang, Z., Wright, J., & Ma, Y. (2010, October). Decomposing background topics from keywords by principal component pursuit. In Proceedings of the 19th ACM international conference on Information and knowledge management (pp. 269-278). ACM. ↩
Feng, J., Xu, H., & Yan, S. (2013). Online robust pca via stochastic optimization. In Advances in Neural Information Processing Systems (pp. 404-412). ↩
Shen, J., Xu, H., & Li, P. (2014). Online optimization for max-norm regularization. In Advances in Neural Information Processing Systems (pp. 1718-1726). ↩
Recht, B., Fazel, M., & Parrilo, P. A. (2010). Guaranteed minimum-rank solutions of linear matrix equations via nuclear norm minimization. SIAM review, 52(3), 471-501. ↩
Srebro, N., Rennie, J. D., & Jaakkola, T. S. (2004, December). Maximum-Margin Matrix Factorization. In NIPS (Vol. 17, pp. 1329-1336). ↩
Zhang, D., Hu, Y., Ye, J., Li, X., & He, X. (2012, June). Matrix completion by truncated nuclear norm regularization. In Computer Vision and Pattern Recognition (CVPR), IEEE Conference on (pp. 2192-2199). ↩