Bao, B.K., et al., Inductive robust principal component analysis. IEEE Transactions on Image Processing, 2012. 21(8): p. 3794-3800.
本文针对经典的 Inductive Robust Principal Component Analysis 的理论方法进行展开详解。本人学术水平有限，文中如有错误之处，敬请指正。

摘要：RPCA 是一种直推式方法，并不能很好地处理新的样本。如果一个新的数据加入，RPCA需要重新计算所有的数据，导致非常高的计算代价。所以，RPCA 是一个不适用于在线计算（数据会分批按顺序加入）的应用。为了克服这个难题，此文提出了一个 inductive robust principal component analysis （IRPCA）方法，给定一组训练数据，和 RPCA 不同的是，（RPCA 旨在恢复原始数据矩阵），IRPCA 目标是学习一个潜在的投影矩阵，可以有效地去除数据中的损坏。

1 简介

假设有一个误差修正的问题：数据组成部分如下

x=y+e,x∈Rd,(1)

\begin{equation} \tag{1} x = y + e, \quad x \in \mathbb{R}^d, \end{equation}
其中 yy 是低秩的子空间部分，ee 是误差项。给定任意的数据向量 xx，求解的目标是从原始数据中分离出向量 yy 和误差 ee。

上述问题等价于找出数据中的主成分，因为低秩的子空间可以被退化的高斯分布很好地建模。前提是误差服从小方差的高斯分布，广泛使用的主成分分析（PCA）可以有效地处理上述问题。有 nn 个训练数据， X=[x1,x2,⋯,xn]X = [x_1,x_2,\cdots,x_n]，PCA 通过最小化重建误差，来学习低秩投影，

minU ||X−UUTX||2F , s.t. UTU=Ir ,(2)

\begin{equation} \tag{2}\min_U \ || X-UU^TX ||_F^2 \ , \ \ \mathrm{s.t.} \ U^T U = I_r \ , \end{equation}
其中 IrI_r 是一个 r×rr \times r 的单位矩阵， ||⋅||2F||\cdot||_F^2 表示 Frobenious 范数， rr 是 UU 的秩（ UU 的列数）。这个优化问题可以有效地用奇异值分解（SVD）求出。假设 U∗U^* 是从训练数据 XX 求出的解，有了一个数据 xx，它的主成分可以估计出

y=U∗(U∗)Tx.(3)

\begin{equation} \tag{3}y = U^*(U^*)^T x. \end{equation}
这个方法计算很有效，而且稳定，使 PCA 被广泛用于误差修正。然而，在实际应用中，PCA 对严重的损坏效果很差，通常是偏离真实的子空间。

为了克服 PCA 的缺陷，Wright et al. 1 2 提出了 robust principal componen analysis (RPCA) ，成功应用于各个方面。
同样给定一组数据 X=[x1,x2,⋯,xn]X = [x_1, x_2, \cdots, x_n]，RPCA 找到其主成分 Y=[y1,y2,⋯,yn]Y = [y_1, y_2, \cdots, y_n]，通过如下的凸优化问题

minY,E ||Y||∗+λ||E||1, s.t. X=Y+E,(4)

\begin{equation} \tag{4}\min_{Y,E} \ ||Y||_* + \lambda ||E||_1, \ \ \mathrm{s.t.} \ X = Y + E, \end{equation}
其中 ||⋅||∗||\cdot||_* 表示核范数（矩阵的奇异值之和）， ||⋅||1||\cdot||_1 表示 ℓ1\ell_1 范数， λ\lambda 是权衡参数。假设损坏是足够稀疏的。这个稀疏假设是合理的、现实中普遍存在的。然而 RPCA 很难扩展已学习到的模型到新的数据上。可能处理测试数据的步骤： Y∗,E∗Y^*,E^* 是从训练数据学习到的解，计算 Y∗Y^* 的 SVD， Y∗=U∗Σ∗(V∗)TY^* = U^* \Sigma^* (V^*)^T；然后用 y=U∗(U∗)Txy = U^*(U^*)^T x 处理新数据。这个方法的实际效果很差。更具体来说， (Y∗,E∗)↔(U∗(U∗)TX,X−U∗(U∗)TX)(Y^*,E^*) \leftrightarrow (U^*(U^*)^T X, X - U^*(U^*)^T X) 其实不接近，重构的解是不正确的。并且 E∗E^* 不是稀疏的。RPCA 不适用于在线计算的应用。

于是，此文提出了 inductive principal componen analysis (IRPCA)，不仅可以处理严重的损坏（相比于 RPCA），并且有良好的泛化能力。关键在于：从训练数据中学习出一个低秩投影，它能有效地移除误差，并把数据投影到其潜在的子空间中。

2 相关工作

Liu et al. 3 （基于 RPCA）提出了 low-rank representation (LRR) 模型

minZ,E ||Z||∗+λ||E||ℓ, s.t. X=AZ+E,(5)

\begin{equation} \tag{5}\min_{Z,E} \ ||Z||_* + \lambda ||E||_\ell, \ \ \mathrm{s.t.} \ X = AZ + E, \end{equation}
其中 AA 是数据空间的字典，||⋅||ℓ||\cdot||_\ell 根据具体情况选择合适的范数。如果 A=IA = I，那么 LRR 就变成了 RPCA 。子空间分割可以用 LRR 解决

minZ,E ||Z|∗+λ||E||2,1, s.t. X=XZ+E,(6)

\begin{equation} \tag{6}\min_{Z,E} \ ||Z|_* + \lambda ||E||_{2,1}, \ \ \mathrm{s.t.} \ X = XZ + E, \end{equation}
其中区别是 A=XA = X， ||E||2,1=∑nj=1∑di=1([E]ij)2−−−−−−−−−−√||E||_{2,1} = \sum_{j=1}^{n} \sqrt{\sum_{i=1}^{d} ([E]_{ij})^2}。这个方法也有 RPCA 同样的缺点：不能很好地处理新数据，需要重新计算，代价很高。

3 目标优化求解

3.1 模型建立

如果损坏没有任何的限制，那么 IRPCA 也是不实用的。比如，损坏是无序的，一般情况下没有简单的模型可以拟合它。幸运的是，IRPCA 是可行的。首先，即使损坏是乱序的，也能存在一个线性投影 P0P_0 把数据投影到子空间中，能正确地恢复出数据（即使不是完全精确的恢复）；

其次，两个高维的向量，通常是独立的，近似是相互正交 4，也就是说，损坏通常不在正确的子空间中。这样的情况下， P0P_0 严格地从数据中去除损坏。只要有数据 xx，其主成分就可以通过 y=P0xy = P_0 x 获得。

开始建立优化模型，训练数据 X=[x1,x2,⋯,xn]X = [x_1,x_2,\cdots,x_n]。目的是学习低秩投影 P0P_0，给出如下优化函数

minP,E rank(P)+λ||E||0, s.t. X=PX+E ,(7)

\begin{equation} \tag{7}\min_{P,E} \ \mathrm{rank}(P) + \lambda ||E||_0, \ \ \mathrm{s.t.} \ X = PX + E \ , \end{equation}
其中 λ\lambda 是系数， ||⋅||0||\cdot||_0 是 ℓ0\ell_0 范数。该目标函数是不连续性的，因为 rank(⋅)\mathrm{rank}(\cdot) 和 ℓ0\ell_0 存在。根据通常的方法，采用核范数代替。同时， ℓ0\ell_0 也用 ℓ1\ell_1 代替。重新构建的凸的优化问题

minP,E ||P||∗+λ||E||1, s.t. X=PX+E .(8)

\begin{equation} \tag{8}\min_{P,E} \ ||P||_* + \lambda ||E||_1, \ \ \mathrm{s.t.} \ X = PX + E \ . \end{equation}
假设最优解 P∗P^* 已经得到。对于一个新的数据 xx，我们估计其主成分 y=P∗xy = P^* x 和误差 e=x−P∗xe = x - P^* x，计算很快速、简便。

3.2 优化求解

首先将上述问题转化为其转置的形式

minP,E ||PT||∗+λ||ET||1, s.t. XT=PTXT+ET .(9)

\begin{equation} \tag{9}\min_{P,E} \ ||P^T||_* + \lambda ||E^T||_1, \ \ \mathrm{s.t.} \ X^T = P^T X^T + E^T \ . \end{equation}
根据论文 5 中提到，其计算复杂度是 O(d3)O(d^3)，对高维数据代价更大。考虑到计算的效率，此文不直接求解该问题，而是到一个更简单的形式（参考 Theorem 1 6）： P∗P^* 总是在 XX 的各列分布的字空间中。P∗P^* 可以被分解为 P∗=L∗(Q∗)TP^* = L^* (Q^*)^T，其中 Q∗Q^* 可以由正交化 XX 的各列得到。于是，问题被转化为如下的形式

minL,E ||J||∗+λ||E||1, s.t. X=LA+E, L=J,(10)

\begin{equation} \tag{10}\min_{L,E} \ ||J||_* + \lambda ||E||_1, \ \ \mathrm{s.t.} \ X = LA + E , \ L = J, \end{equation}
其中 A=(Q∗)TXA = (Q^*)^T X， ||⋅||∗||\cdot||_* 表示核范数， ||⋅||1||\cdot||_1 表示 ℓ1\ell_1 范数， λ\lambda 是权衡参数。

构建 Lagrangian 函数

L(L,E,J)= ||J||∗+λ||E||1+tr(YT1(X−LA−E))+tr(YT2(L−J))+μ2(||X−LA−E||2F+||L−J||2F),(11)

\begin{align} \mathcal{L}(L,E,J) =&\ ||J||_*+\lambda ||E||_1 + \mathrm{tr}\left( Y_1^T (X-LA-E) \right) + \mathrm{tr} \left( Y_2^T (L-J) \right) \\ & + \frac{\mu}{2} \left( ||X-LA-E||_F^2 + ||L-J||_F^2 \right), \tag{11} \end{align}
其中 Y1,Y2∈Rm×nY_1,Y_2 \in \mathbb{R}^{m \times n} 分别是 Lagrange 乘子矩阵， tr(⋅)\mathrm{tr}(\cdot) 是迹函数， μ\mu 是惩罚项系数（非变量）， ||⋅||F||\cdot||_F 表示 Frobenius 范数。显然，优化目标函数可以给出

(Lk+1,Ek+1,Jk+1)Y1k+1Y2k+1μ=argminL,E,J L(L,E,J,Y1k,Y2k),=Y1k+μ(X−LA−E),=Y2k+μ(L−J),=min(ρμ,μmax).

\begin{align} (L_{k+1}, E_{k+1},J_{k+1}) &= \arg \min_{L,E,J} \ \mathcal{L}(L,E,J, Y_{1\,k},Y_{2\,k}), \\ Y_{1\,k+1} &= Y_{1\,k} + \mu (X-LA-E), \\ Y_{2\,k+1} &= Y_{2\,k} + \mu (L-J), \\ \mu &= \min(\rho \mu, \mu_{\max}). \end{align}
其中 ρ>1\rho > 1 是一个常数，用于不断增加 μ\mu 的值。
由于优化目标函数中含有多个变量，通常的做法是每次最小化求解一个变量，而固定其他变量，之后更新 Lagrange 乘子，反复迭代直至收敛。论文使用经典的 inexact Augmented Lagrange Multiplier Method，也叫交替方向乘子法（Alternating Direction Multiplier Method），

Lk+1Ek+1Jk+1=argminL L(L,Ek,Jk,Y1k,Y2k),=argminE L(Lk+1,E,Jk,Y1k,Y2k),=argminJ L(Lk+1,Ek+1,J,Y1k,Y2k).(12)(13)(14)

\begin{align} L_{k+1} &= \arg \min_{L} \ \mathcal{L} (L, E_k, J_k, Y_{1\,k}, Y_{2\,k}), \tag{12} \\ E_{k+1} &= \arg \min_{E} \ \mathcal{L} (L_{k+1}, E, J_k, Y_{1\,k}, Y_{2\,k}), \tag{13} \\ J_{k+1} &= \arg \min_{J} \ \mathcal{L} (L_{k+1}, E_{k+1}, J, Y_{1\,k}, Y_{2\,k}). \tag{14} \end{align}
迭代更新参数中还加入了 singular value thresholding (SVT) 和 shrinkage 操作，

Sα(x)=sign(x)⋅max{|x|−α, 0},(15)

\begin{equation} \tag{15} S_\alpha (x) = \mathrm{sign}(x) \cdot \max \{ |x| - \alpha,\ 0 \}, \end{equation}
其中 α>0\alpha > 0 是一个设定的阈值。这是一个标量函数，对于矩阵或向量的操作都是 elementwise 的。
给出具体的迭代公式，推导过程见 Appendix，

Lk+1=((X−Ek)AT+Jk+Y1kAT−Y2kμ)(I+AAT)−1;Ek+1=Sλμ(X−Lk+1A+Y1kμ);(U,Σ,V)=svd(Lk+1+Y2kμ),Jk+1=US1μ(Σ)VT;(16)(17)(18)(19)

\begin{align} & L_{k+1} = \left((X-E_k) A^T + J_k + \frac{Y_{1\,k} A^T - Y_{2\,k}}{\mu}\right) \left(I + AA^T \right)^{-1}; \tag{16} \\ &E_{k+1} = S_{\frac{\lambda}{\mu}} \left( X - L_{k+1}A + \frac{Y_{1\,k}}{\mu} \right) ; \tag{17} \\ & (U, \Sigma, V) = \mathrm{svd} \left( L_{k+1} + \frac{Y_{2\,k}}{\mu} \right), \tag{18} \\ & J_{k+1} = U S_{\frac{1}{\mu}}(\Sigma)V^T; \tag{19} \end{align}
其中 Sα(⋅)S_\alpha(\cdot) 就是为了近似代替优化求解过程中的 ||J||∗+λ||E||1||J||_* + \lambda ||E||_1 而加入的；也就是说，在梯度求解 Lk+1,Ek+1,Jk+1L_{k+1},E_{k+1},J_{k+1} 的过程中，并没有考虑这两项，而是用减去 EkE_k 的较小的元素值 Sλ/μ(Ek)S_{\lambda/\mu}(E_k)，和减去 JkJ_k 的部分（较小的）奇异值 S1/μ(Σ)S_{{1}/{\mu}}(Σ) 来代替，直接求闭式解表达式。
优化步骤的迭代停止条件：

||X−Lk+1A−Ek+1||∞<ε and ||Z−Jk+1||∞<ε.(20)

\begin{equation} \tag{20} ||X-L_{k+1}A-E_{k+1}||_\infty
其中 ||⋅||∞||\cdot||_\infty 表示最大范数，定义为矩阵中最大的元素绝对值 ||A||∞=maxi,j{|aij|}||A||_\infty = \max_{i,j} \left\{ |a_{ij}| \right\}。

4 实验

略

Appendix

原论文中的公式有明显的符号错误。

求解 Lk+1L_{k+1}, 当 ∂L∂L=0\frac{\partial \mathcal{L}}{\partial L} = 0 时，（ ||J||∗+λ||E||1||J||_*+\lambda ||E||_1 与 LL 无关）

∂L∂L∂L∂LLk+1=∂∂L{tr(YT1k(X−LA−Ek))+tr(YT2k(L−Jk))+μ2(||X−LA−Ek||2F+||L−Jk||2F)}=−Y1kAT+Y2k−μ(X−LA−Ek)AT+μ(L−Jk)=μL(I+AAT)+(Y2k−Y1kAT)+μ((Ek−X)AT−Jk),=0,=((X−Ek)AT+Jk+Y1kAT−Y2kμ)(I+AAT)−1.

\begin{align} \frac{\partial \mathcal{L}}{\partial L} &= \frac{\partial}{\partial L} \left\{ \mathrm{tr}\left( Y_{1\,k}^T (X-LA-E_k) \right) + \mathrm{tr} \left( Y_{2\,k}^T (L-J_k) \right) + \frac{\mu}{2} \left( ||X-LA-E_k||_F^2 + ||L-J_k||_F^2 \right) \right\} \\ &= -Y_{1\,k} A^T + Y_{2\,k} - \mu (X - LA - E_k) A^T + \mu (L - J_k) \\ &= \mu L (I + AA^T) + (Y_{2\,k} - Y_{1\,k} A^T) + \mu \left((E_k - X) A^T - J_k \right), \\ \frac{\partial \mathcal{L}}{\partial L} &= 0, \\ L_{k+1} &= \left( (X - E_k) A^T + J_k + \frac{Y_{1\,k} A^T - Y_{2\,k}}{\mu}\right) \left(I + AA^T \right)^{-1}. \end{align}
求解 Ek+1E_{k+1}, 当 ∂L∂E=0\frac{\partial \mathcal{L}}{\partial E} = 0 时，（ ||J||∗||J||_* 与 EE 无关，λ||E||1\lambda ||E||_1 用 Sλ/μ(⋅)S_{\lambda / \mu} (\cdot) 近似）

∂L∂E∂L∂EEk+1=∂L∂E{tr(YT1k(X−Lk+1A−E))+μ2||X−Lk+1A−E||2F}=−Y1k+μ(E+Lk+1A−X),=0,=Sλμ(X−Lk+1A+Y1kμ).

\begin{align} \frac{\partial \mathcal{L}}{\partial E} &= \frac{\partial \mathcal{L}}{\partial E} \left\{ \mathrm{tr}\left( Y_{1\,k}^T (X-L_{k+1}A-E) \right) + \frac{\mu}{2} ||X-L_{k+1}A-E||_F^2 \right\} \\ &= -Y_{1\,k} + \mu (E + L_{k+1}A - X), \\ \frac{\partial \mathcal{L}}{\partial E} &= 0, \\ E_{k+1} &= S_{\frac{\lambda}{\mu}} \left(X - L_{k+1}A + \frac{Y_{1\,k}}{\mu} \right). \end{align}
求解 Jk+1J_{k+1}, 当 ∂L∂J=0\frac{\partial \mathcal{L}}{\partial J} = 0 时，（ λ||E||1\lambda ||E||_1 与 JJ 无关，||J||∗||J||_* 用 S1/μ(Σ), Σ=svd(⋅)S_{1/ \mu} (\Sigma),\ \Sigma = \mathrm{svd}(\cdot) 近似）

∂L∂J∂L∂J(U,Jk+1=∂L∂J{tr(YT2k(Lk+1−J))+μ2||Lk+1−J||2F}=−Y2k+μ(J−Lk+1),=0,Σ,V)=svd(Lk+1+Y2kμ),=US1μ(Σ)VT.

\begin{align} \frac{\partial \mathcal{L}}{\partial J} &= \frac{\partial \mathcal{L}}{\partial J} \left\{ \mathrm{tr} \left( Y_{2\,k}^T (L_{k+1}-J) \right) + \frac{\mu}{2} ||L_{k+1}-J||_F^2 \right\} \\ &= -Y_{2\,k} + \mu (J - L_{k+1}), \\ \frac{\partial \mathcal{L}}{\partial J} &= 0, \\(U, &\, \Sigma, V) = \mathrm{svd} \left( L_{k+1} +\frac{Y_{2\,k}}{\mu} \right), \\J_{k+1} &= U S_{\frac{1}{\mu}}(\Sigma)V^T. \end{align}

E. J. Candes, X. Li, Y. Ma, and J. Wright. (2009, Dec.). Robust Principal Component Analysis? [Online]. Available: http://wwwstat.stanford.edu/∼candes/papers/RobustPCA.pdf ↩
J. Wright, A. Ganesh, S. Rao, Y. Peng, and Y. Ma, “Robust principal component analysis: Exact recovery of corrupted low-rank matrices via convex optimization,” in Proc. Neural Inf. Process. Syst., 2009, pp. 1–9. ↩
G. Liu, Z. Lin, and Y. Yu, “Robust subspace segmentation by low-rank representation,” in Proc. Int. Conf. Mach. Learn., 2010, pp. 1–8. ↩
R. R. Hamming, Art of Doing Science and Engineering: Learning to Learn. Boca Raton, FL: CRC Press, 1997, p. 364. ↩
G. Liu, Z. Lin, and Y. Yu, “Robust subspace segmentation by low-rank representation,” in Proc. Int. Conf. Mach. Learn., 2010, pp. 1–8. ↩
G. Liu, Z. Lin, S. Yan, J. Sun, Y. Yu, and Y. Ma. (2010). Robust Recovery of Subspace Structures by Low-Rank Representation [Online]. Available: http://arxiv.org/pdf/1010.2955.pdf ↩

笔记：Inductive Robust Principal Component Analysis相关推荐

Inductive Robust Principal Component Analysis
本来想转载原文,不过复制好像有问题,附上链接: https://blog.csdn.net/xueshengke/article/details/53206726?utm_source=blogxgw ...
笔记：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 ...
Robust principal component analysis?(RPCA简单理解)
参考文献:Candès, E.J., Li, X., Ma, Y., and Wright, J.: 'Robust principal component analysis?', J. ACM, 2 ...
Robust Principal Component Analysis?（PCP）
目录引一些微弱的假设: 问题的解决理论去随机 Dual Certificates(对偶保证?) Golfing Scheme 数值实验代码 Candes E J, Li X, Ma Y, e ...
Principal Component Analysis 主成分分析【学习笔记】
Principal Component Analysis 主成分分析笔记文章目录 Principal Component Analysis 主成分分析笔记 PCA是什么? 数据的线性变换拉伸操作: ...
PCA（Principal Component Analysis）的原理、算法步骤和实现。
PCA的原理介绍: PCA(Principal Component Analysis)是一种常用的数据分析方法.PCA通过线性变换将原始数据变换为一组各维度线性无关的表示,可用于提取数据的主要特征分 ...
机器学习与高维信息检索 - Note 7 - 核主成分分析（Kernel Principal Component Analysis，K-PCA）
Note 7 - 核主成分分析(Kernel Principal Component Analysis) 核主成分分析 Note 7 - 核主成分分析(Kernel Principal Compone ...
机器学习与高维信息检索 - Note 4 - 主成分分析及其现代解释（Principal Component Analysis, PCA）及相关实例
主成分分析及其现代解释 4. 主成分分析及其现代解释 Principal Component Analysis and Its Modern Interpretations 4.1 几何学解释 The ...
PCA（principal component analysis）主成分分析降维和KPCA（kernel principal component analysis）核
PCA(principal component analysis)主成分分析降维和KPCA(kernel principal component analysis)核主成分分析降维方法详解及实战 PC ...

笔记：Inductive Robust Principal Component Analysis